text:
module BinaryPlatforms

export platform_key_abi, platform_dlext, valid_dl_path, arch, libc, compiler_abi,
       libgfortran_version, libstdcxx_version, cxxstring_abi, parse_dl_name_version,
       detect_libgfortran_version, detect_libstdcxx_version, detect_cxxstring_abi,
       call_abi, wordsize, triplet, select_platform, platforms_match,
       CompilerABI, Platform, UnknownPlatform, Linux, MacOS, Windows, FreeBSD
import Base: show
import Libdl

abstract type Platform end

"""
    UnknownPlatform

A placeholder `Platform` that signifies an unknown platform.
"""
struct UnknownPlatform <: Platform
    # Just swallow up whatever arguments get passed to you
    UnknownPlatform(args...; kwargs...) = new()
end

# We need to track our compiler ABI compatibility.
struct CompilerABI
    # libgfortran SOVERSION we're linked against (if any)
    libgfortran_version::Union{Nothing,VersionNumber}

    # libstdc++ SOVERSION we're linked against (if any)
    libstdcxx_version::Union{Nothing,VersionNumber}

    # Whether we're using cxx11abi strings, not using them, or don't care
    # This is only relevant when linked against `libstdc++`, when linked against
    # libc++ or none (because it's not C++ code) we don't care.
    # Valid Symbol values are `:cxx03` and `:cxx11`
    cxxstring_abi::Union{Nothing,Symbol}

    function CompilerABI(;libgfortran_version::Union{Nothing, VersionNumber} = nothing,
                         libstdcxx_version::Union{Nothing, VersionNumber} = nothing,
                         cxxstring_abi::Union{Nothing, Symbol} = nothing)
        if libgfortran_version !== nothing && (libgfortran_version < v"3" ||
                                              libgfortran_version >= v"6")
            throw(ArgumentError("Unsupported libgfortran '$libgfortran_version'"))
        end

        if libstdcxx_version !== nothing && (libstdcxx_version < v"3.4.0" ||
                                            libstdcxx_version >= v"3.5")
            throw(ArgumentError("Unsupported libstdc++ '$libstdcxx_version'"))
        end

        if cxxstring_abi !== nothing && !in(cxxstring_abi, [:cxx03, :cxx11])
            throw(ArgumentError("Unsupported string ABI '$cxxstring_abi'"))
        end

        return new(libgfortran_version, libstdcxx_version, cxxstring_abi)
    end
end

# Easy replacement constructor
function CompilerABI(cabi::CompilerABI; libgfortran_version=nothing,
                                        libstdcxx_version=nothing,
                                        cxxstring_abi=nothing)
    lgv = something(libgfortran_version, Some(cabi.libgfortran_version))
    lsv = something(libstdcxx_version, Some(cabi.libstdcxx_version))
    ca = something(cxxstring_abi, Some(cabi.cxxstring_abi))
    return CompilerABI(;libgfortran_version=lgv, libstdcxx_version=lsv, cxxstring_abi=ca)
end

libgfortran_version(cabi::CompilerABI) = cabi.libgfortran_version
libstdcxx_version(cabi::CompilerABI) = cabi.libstdcxx_version
cxxstring_abi(cabi::CompilerABI) = cabi.cxxstring_abi

function show(io::IO, cabi::CompilerABI)
    args = String[]
    if cabi.libgfortran_version !== nothing
        push!(args, "libgfortran_version=$(repr(cabi.libgfortran_version))")
    end
    if cabi.libstdcxx_version !== nothing
        push!(args, "libstdcxx_version=$(repr(cabi.libstdcxx_version))")
    end
    if cabi.cxxstring_abi !== nothing
        push!(args, "cxxstring_abi=$(repr(cabi.cxxstring_abi))")
    end
    write(io, "CompilerABI($(join(args, ", ")))")
end

struct Linux <: Platform
    arch::Symbol
    libc::Union{Nothing,Symbol}
    call_abi::Union{Nothing,Symbol}
    compiler_abi::CompilerABI

    function Linux(arch::Symbol;
                   libc::Union{Nothing,Symbol}=nothing,
                   call_abi::Union{Nothing,Symbol}=nothing,
                   compiler_abi::CompilerABI=CompilerABI())
        if !in(arch, [:i686, :x86_64, :aarch64, :powerpc64le, :armv6l, :armv7l])
            throw(ArgumentError("Unsupported architecture '$arch' for Linux"))
        end

        # The default libc on Linux is glibc
        if libc === nothing
            libc = :glibc
        end

        if !in(libc, [:glibc, :musl])
            throw(ArgumentError("Unsupported libc '$libc' for Linux"))
        end

        # Auto-map the `call_abi` to be `eabihf` on armv6l/armv7l
        if call_abi === nothing && arch in (:armv6l, :armv7l)
            call_abi = :eabihf
        end

        if !in(call_abi, [:eabihf, nothing])
            throw(ArgumentError("Unsupported calling abi '$call_abi' for Linux"))
        end

        # If we're constructing for armv7l/armv6l, we MUST have the eabihf abi
        if arch in (:armv6l, :armv7l) && call_abi !== :eabihf
            throw(ArgumentError("armv6l/armv7l Linux must use eabihf, not '$call_abi'"))
        end
        # ...and vice-versa
        if !(arch in (:armv6l, :armv7l)) && call_abi === :eabihf
            throw(ArgumentError("eabihf Linux is only supported on armv6l/armv7l, not '$arch'!"))
        end

        return new(arch, libc, call_abi, compiler_abi)
    end
end

struct MacOS <: Platform
    arch::Symbol
    libc::Nothing
    call_abi::Nothing
    compiler_abi::CompilerABI

    # Provide defaults for everything because there's really only one MacOS
    # target right now.  Maybe someday iOS.  :fingers_crossed:
    function MacOS(arch::Symbol=:x86_64;
                   libc::Union{Nothing,Symbol}=nothing,
                   call_abi::Union{Nothing,Symbol}=nothing,
                   compiler_abi::CompilerABI=CompilerABI())
        if arch !== :x86_64
            throw(ArgumentError("Unsupported architecture '$arch' for macOS"))
        end
        if libc !== nothing
            throw(ArgumentError("Unsupported libc '$libc' for macOS"))
        end
        if call_abi !== nothing
            throw(ArgumentError("Unsupported abi '$call_abi' for macOS"))
        end

        return new(arch, libc, call_abi, compiler_abi)
    end
end

struct Windows <: Platform
    arch::Symbol
    libc::Nothing
    call_abi::Nothing
    compiler_abi::CompilerABI

    function Windows(arch::Symbol;
                     libc::Union{Nothing,Symbol}=nothing,
                     call_abi::Union{Nothing,Symbol}=nothing,
                     compiler_abi::CompilerABI=CompilerABI())
        if !in(arch, [:i686, :x86_64])
            throw(ArgumentError("Unsupported architecture '$arch' for Windows"))
        end
        # We only support the one libc/abi on Windows, so no need to play
        # around with "default" values.
        if libc !== nothing
            throw(ArgumentError("Unsupported libc '$libc' for Windows"))
        end
        if call_abi !== nothing
            throw(ArgumentError("Unsupported abi '$call_abi' for Windows"))
        end

        return new(arch, libc, call_abi, compiler_abi)
    end
end

struct FreeBSD <: Platform
    arch::Symbol
    libc::Nothing
    call_abi::Union{Nothing,Symbol}
    compiler_abi::CompilerABI

    function FreeBSD(arch::Symbol=:x86_64;
                     libc::Union{Nothing,Symbol}=nothing,
                     call_abi::Union{Nothing,Symbol}=nothing,
                     compiler_abi::CompilerABI=CompilerABI())
        # `uname` on FreeBSD reports its architecture as amd64 and i386 instead of x86_64
        # and i686, respectively. In the off chance that Julia hasn't done the mapping for
        # us, we'll do it here just in case.
        if arch === :amd64
            arch = :x86_64
        elseif arch === :i386
            arch = :i686
        elseif !in(arch, [:i686, :x86_64, :aarch64, :powerpc64le, :armv6l, :armv7l])
            throw(ArgumentError("Unsupported architecture '$arch' for FreeBSD"))
        end

        # The only libc we support on FreeBSD is the blank libc, which corresponds to
        # FreeBSD's default libc
        if libc !== nothing
            throw(ArgumentError("Unsupported libc '$libc' for FreeBSD"))
        end

        # Auto-map the `call_abi` to be `eabihf` on armv6l/armv7l
        if call_abi === nothing && arch in (:armv6l, :armv7l)
            call_abi = :eabihf
        end

        if !in(call_abi, [:eabihf, nothing])
            throw(ArgumentError("Unsupported calling abi '$call_abi' for FreeBSD"))
        end

        # If we're constructing for armv7l, we MUST have the eabihf abi
        if arch in (:armv6l, :armv7l) && call_abi !== :eabihf
            throw(ArgumentError("armv6l/armv7l FreeBSD must use eabihf, not '$call_abi'"))
        end
        # ...and vice-versa
        if !(arch in (:armv6l, :armv7l)) && call_abi === :eabihf
            throw(ArgumentError("eabihf FreeBSD is supported only on armv6l/armv7l, not '$arch'!"))
        end

        return new(arch, libc, call_abi, compiler_abi)
    end
end

"""
    platform_name(p::Platform)

Get the "platform name" of the given platform.  E.g. returns "Linux" for a
`Linux` object, or "Windows" for a `Windows` object.
"""
platform_name(p::Linux) = "Linux"
platform_name(p::MacOS) = "MacOS"
platform_name(p::Windows) = "Windows"
platform_name(p::FreeBSD) = "FreeBSD"
platform_name(p::UnknownPlatform) = "UnknownPlatform"

"""
    arch(p::Platform)

Get the architecture for the given `Platform` object as a `Symbol`.

# Examples
```jldoctest
julia> arch(Linux(:aarch64))
:aarch64

julia> arch(MacOS())
:x86_64
```
"""
arch(p::Platform) = p.arch
arch(u::UnknownPlatform) = nothing

"""
    libc(p::Platform)

Get the libc for the given `Platform` object as a `Symbol`.

# Examples
```jldoctest
julia> libc(Linux(:aarch64))
:glibc

julia> libc(FreeBSD(:x86_64))
```
"""
libc(p::Platform) = p.libc
libc(u::UnknownPlatform) = nothing

"""
   call_abi(p::Platform)

Get the calling ABI for the given `Platform` object, returns either `nothing` (which
signifies a "default choice") or a `Symbol`.

# Examples
```jldoctest
julia> call_abi(Linux(:x86_64))

julia> call_abi(FreeBSD(:armv7l))
:eabihf
```
"""
call_abi(p::Platform) = p.call_abi
call_abi(u::UnknownPlatform) = nothing

"""
    compiler_abi(p::Platform)

Get the compiler ABI object for the given `Platform`
# Examples
```jldoctest
julia> compiler_abi(Linux(:x86_64))
CompilerABI()
```
"""
compiler_abi(p::Platform) = p.compiler_abi
compiler_abi(p::UnknownPlatform) = CompilerABI()

# Also break out CompilerABI getters for our platforms
libgfortran_version(p::Platform) = libgfortran_version(compiler_abi(p))
libstdcxx_version(p::Platform) = libstdcxx_version(compiler_abi(p))
cxxstring_abi(p::Platform) = cxxstring_abi(compiler_abi(p))

"""
    wordsize(platform)

Get the word size for the given `Platform` object.

# Examples
```jldoctest
julia> wordsize(Linux(:armv7l))
32

julia> wordsize(MacOS())
64
```
"""
wordsize(p::Platform) = (arch(p) === :i686 || arch(p) === :armv7l) ? 32 : 64
wordsize(u::UnknownPlatform) = 0

"""
    triplet(platform)

Get the target triplet for the given `Platform` object as a `String`.

# Examples
```jldoctest
julia> triplet(MacOS())
"x86_64-apple-darwin14"

julia> triplet(Windows(:i686))
"i686-w64-mingw32"

julia> triplet(Linux(:armv7l; compiler_abi=CompilerABI(;libgfortran_version=v"3")))
"armv7l-linux-gnueabihf-libgfortran3"
```
"""
triplet(p::Platform) = string(
    arch_str(p),
    vendor_str(p),
    libc_str(p),
    call_abi_str(p),
    compiler_abi_str(p),
)
vendor_str(p::Windows) = "-w64-mingw32"
vendor_str(p::MacOS) = "-apple-darwin14"
vendor_str(p::Linux) = "-linux"
vendor_str(p::FreeBSD) = "-unknown-freebsd11.1"

# Special-case UnknownPlatform
triplet(p::UnknownPlatform) = "unknown-unknown-unknown"

# Helper functions for Linux and FreeBSD libc/abi mishmashes
arch_str(p::Platform) = string(arch(p))
function libc_str(p::Platform)
    if libc(p) === nothing
        return ""
    elseif libc(p) === :glibc
        return "-gnu"
    else
        return "-$(libc(p))"
    end
end
call_abi_str(p::Platform) = (call_abi(p) === nothing) ? "" : string(call_abi(p))
function compiler_abi_str(cabi::CompilerABI)
    str = ""
    if cabi.libgfortran_version !== nothing
        str *= "-libgfortran$(cabi.libgfortran_version.major)"
    end

    if cabi.libstdcxx_version !== nothing
        str *= "-libstdcxx$(libstdcxx_version(cabi).patch)"
    end

    if cabi.cxxstring_abi !== nothing
        str *= "-$(cabi.cxxstring_abi)"
    end
    return str
end
compiler_abi_str(p::Platform) = compiler_abi_str(compiler_abi(p))

Sys.isapple(p::Platform) = p isa MacOS
Sys.islinux(p::Platform) = p isa Linux
Sys.iswindows(p::Platform) = p isa Windows
Sys.isbsd(p::Platform) = (p isa FreeBSD) || (p isa MacOS)

"""
    platform_key_abi(machine::AbstractString)

Returns the platform key for the current platform, or any other though the
the use of the `machine` parameter.
"""
function platform_key_abi(machine::AbstractString)
    # We're going to build a mondo regex here to parse everything:
    arch_mapping = Dict(
        :x86_64 => "(x86_|amd)64",
        :i686 => "i\\d86",
        :aarch64 => "aarch64",
        :armv7l => "arm(v7l)?", # if we just see `arm-linux-gnueabihf`, we assume it's `armv7l`
        :armv6l => "armv6l",
        :powerpc64le => "p(ower)?pc64le",
    )
    platform_mapping = Dict(
        :darwin => "-apple-darwin[\\d\\.]*",
        :freebsd => "-(.*-)?freebsd[\\d\\.]*",
        :mingw32 => "-w64-mingw32",
        :linux => "-(.*-)?linux",
    )
    libc_mapping = Dict(
        :libc_nothing => "",
        :glibc => "-gnu",
        :musl => "-musl",
    )
    call_abi_mapping = Dict(
        :call_abi_nothing => "",
        :eabihf => "eabihf",
    )
    libgfortran_version_mapping = Dict(
        :libgfortran_nothing => "",
        :libgfortran3 => "(-libgfortran3)|(-gcc4)", # support old-style `gccX` versioning
        :libgfortran4 => "(-libgfortran4)|(-gcc7)",
        :libgfortran5 => "(-libgfortran5)|(-gcc8)",
    )
    libstdcxx_version_mapping = Dict(
        :libstdcxx_nothing => "",
        # This is sadly easier than parsing out the digit directly
        (Symbol("libstdcxx$(idx)") => "-libstdcxx$(idx)" for idx in 18:26)...,
    )
    cxxstring_abi_mapping = Dict(
        :cxxstring_nothing => "",
        :cxx03 => "-cxx03",
        :cxx11 => "-cxx11",
    )

    # Helper function to collapse dictionary of mappings down into a regex of
    # named capture groups joined by "|" operators
    c(mapping) = string("(",join(["(?<$k>$v)" for (k, v) in mapping], "|"), ")")

    triplet_regex = Regex(string(
        "^",
        c(arch_mapping),
        c(platform_mapping),
        c(libc_mapping),
        c(call_abi_mapping),
        c(libgfortran_version_mapping),
        c(libstdcxx_version_mapping),
        c(cxxstring_abi_mapping),
        "\$",
    ))

    m = match(triplet_regex, machine)
    if m !== nothing
        # Helper function to find the single named field within the giant regex
        # that is not `nothing` for each mapping we give it.
        get_field(m, mapping) = begin
            for k in keys(mapping)
                if m[k] !== nothing
                    strk = string(k)
                    # Convert our sentinel `nothing` values to actual `nothing`
                    if endswith(strk, "_nothing")
                        return nothing
                    end
                    # Convert libgfortran/libstdcxx version numbers
                    if startswith(strk, "libgfortran")
                        return VersionNumber(parse(Int,strk[12:end]))
                    elseif startswith(strk, "libstdcxx")
                        return VersionNumber(3, 4, parse(Int,strk[10:end]))
                    else
                        return k
                    end
                end
            end
        end

        # Extract the information we're interested in:
        arch = get_field(m, arch_mapping)
        platform = get_field(m, platform_mapping)
        libc = get_field(m, libc_mapping)
        call_abi = get_field(m, call_abi_mapping)
        libgfortran_version = get_field(m, libgfortran_version_mapping)
        libstdcxx_version = get_field(m, libstdcxx_version_mapping)
        cxxstring_abi = get_field(m, cxxstring_abi_mapping)

        # First, figure out what platform we're dealing with, then sub that off
        # to the appropriate constructor.  If a constructor runs into trouble,
        # catch the error and return `UnknownPlatform()` here to be nicer to client code.
        ctors = Dict(:darwin => MacOS, :mingw32 => Windows, :freebsd => FreeBSD, :linux => Linux)
        try
            T = ctors[platform]
            compiler_abi = CompilerABI(;
                libgfortran_version=libgfortran_version,
                libstdcxx_version=libstdcxx_version,
                cxxstring_abi=cxxstring_abi
            )
            return T(arch, libc=libc, call_abi=call_abi, compiler_abi=compiler_abi)
        catch
        end
    end

    @warn("Platform `$(machine)` is not an officially supported platform")
    return UnknownPlatform()
end

# Define show() for these Platform objects for two reasons:
#  - I don't like the `BinaryProvider.` at the beginning of the types;
#    it's unnecessary as these are exported
#  - I like to auto-expand non-`nothing` arguments
function show(io::IO, p::Platform)
    write(io, "$(platform_name(p))($(repr(arch(p)))")

    if libc(p) !== nothing
        write(io, ", libc=$(repr(libc(p)))")
    end
    if call_abi(p) !== nothing
        write(io, ", call_abi=$(repr(call_abi(p)))")
    end
    if compiler_abi(p) != CompilerABI()
        write(io, ", compiler_abi=$(repr(compiler_abi(p)))")
    end
    write(io, ")")
end

"""
    platform_dlext(platform::Platform = platform_key_abi())

Return the dynamic library extension for the given platform, defaulting to the
currently running platform.  E.g. returns "so" for a Linux-based platform,
"dll" for a Windows-based platform, etc...
"""
platform_dlext(::Linux) = "so"
platform_dlext(::FreeBSD) = "so"
platform_dlext(::MacOS) = "dylib"
platform_dlext(::Windows) = "dll"
platform_dlext(::UnknownPlatform) = "unknown"
platform_dlext() = platform_dlext(platform_key_abi())

"""
    parse_dl_name_version(path::AbstractString, platform::Platform)

Given a path to a dynamic library, parse out what information we can
from the filename.  E.g. given something like "lib/libfoo.so.3.2",
this function returns `"libfoo", v"3.2"`.  If the path name is not a
valid dynamic library, this method throws an error.  If no soversion
can be extracted from the filename, as in "libbar.so" this method
returns `"libbar", nothing`.
"""
function parse_dl_name_version(path::AbstractString, platform::Platform)
    dlext_regexes = Dict(
        # On Linux, libraries look like `libnettle.so.6.3.0`
        "so" => r"^(.*?).so((?:\.[\d]+)*)$",
        # On OSX, libraries look like `libnettle.6.3.dylib`
        "dylib" => r"^(.*?)((?:\.[\d]+)*).dylib$",
        # On Windows, libraries look like `libnettle-6.dll`
        "dll" => r"^(.*?)(?:-((?:[\.\d]+)*))?.dll$"
    )

    # Use the regex that matches this platform
    dlregex = dlext_regexes[platform_dlext(platform)]
    m = match(dlregex, basename(path))
    if m === nothing
        throw(ArgumentError("Invalid dynamic library path '$path'"))
    end

    # Extract name and version
    name = m.captures[1]
    version = m.captures[2]
    if version === nothing || isempty(version)
        version = nothing
    else
        version = VersionNumber(strip(version, '.'))
    end
    return name, version
end

"""
    valid_dl_path(path::AbstractString, platform::Platform)

Return `true` if the given `path` ends in a valid dynamic library filename.
E.g. returns `true` for a path like `"usr/lib/libfoo.so.3.5"`, but returns
`false` for a path like `"libbar.so.f.a"`.
"""
function valid_dl_path(path::AbstractString, platform::Platform)
    try
        parse_dl_name_version(path, platform)
        return true
    catch
        return false
    end
end

"""
    detect_libgfortran_version(libgfortran_name::AbstractString)

Examines the given libgfortran SONAME to see what version of GCC corresponds
to the given libgfortran version.
"""
function detect_libgfortran_version(libgfortran_name::AbstractString, platform::Platform = default_platkey)
    name, version = parse_dl_name_version(libgfortran_name, platform)
    if version === nothing
        # Even though we complain about this; we allow it to continue, in the hopes
        # that we shall march on to a BRIGHTER TOMORROW, one in which we are not shackled
        # by the constraints of libgfortran compiler ABIs on our precious programming
        # languages; one where the mistakes of yesterday are mere memories and not
        # continual maintenance burdens upon the children of tomorrow; one where numeric
        # code can be cleanly implemented in a modern language and not bestowed onto the
        # next generation by grizzled ancients, documented only with a faded yellow
        # sticky note that says simply "good luck".
        @warn("Unable to determine libgfortran version from '$(libgfortran_name)'")
    end
    return version
end

"""
    detect_libgfortran_version()

Inspects the current Julia process to determine the libgfortran version this Julia is
linked against (if any).
"""
function detect_libgfortran_version(;platform::Platform = default_platkey)
    libgfortran_paths = filter(x -> occursin("libgfortran", x), Libdl.dllist())
    if isempty(libgfortran_paths)
        # One day, I hope to not be linking against libgfortran in base Julia
        return nothing
    end
    return detect_libgfortran_version(first(libgfortran_paths), platform)
end

"""
    detect_libstdcxx_version()

Inspects the currently running Julia process to find out what version of libstdc++
it is linked against (if any).
"""
function detect_libstdcxx_version()
    libstdcxx_paths = filter(x -> occursin("libstdc++", x), Libdl.dllist())
    if isempty(libstdcxx_paths)
        # This can happen if we were built by clang, so we don't link against
        # libstdc++ at all.
        return nothing
    end

    # Brute-force our way through GLIBCXX_* symbols to discover which version we're linked against
    hdl = Libdl.dlopen(first(libstdcxx_paths))
    for minor_version in 26:-1:18
        if Libdl.dlsym(hdl, "GLIBCXX_3.4.$(minor_version)"; throw_error=false) !== nothing
            Libdl.dlclose(hdl)
            return VersionNumber("3.4.$(minor_version)")
        end
    end
    Libdl.dlclose(hdl)
    return nothing
end

"""
    detect_cxxstring_abi()

Inspects the currently running Julia process to see what version of the C++11 string ABI
it was compiled with (this is only relevant if compiled with `g++`; `clang` has no
incompatibilities yet, bless its heart).  In reality, this actually checks for symbols
within LLVM, but that is close enough for our purposes, as you can't mix configurations
between Julia and LLVM; they must match.
"""
function detect_cxxstring_abi()
    # First, if we're not linked against libstdc++, then early-exit because this doesn't matter.
    libstdcxx_paths = filter(x -> occursin("libstdc++", x), Libdl.dllist())
    if isempty(libstdcxx_paths)
        # We were probably built by `clang`; we don't link against `libstdc++`` at all.
        return nothing
    end

    function open_libllvm(f::Function)
        for lib_name in ("libLLVM", "LLVM", "libLLVMSupport")
            hdl = Libdl.dlopen_e(lib_name)
            if hdl != C_NULL
                try
                    return f(hdl)
                finally
                    Libdl.dlclose(hdl)
                end
            end
        end
        error("Unable to open libLLVM!")
    end

    return open_libllvm() do hdl
        # Check for llvm::sys::getProcessTriple(), first without cxx11 tag:
        if Libdl.dlsym_e(hdl, "_ZN4llvm3sys16getProcessTripleEv") != C_NULL
            return :cxx03
        elseif Libdl.dlsym_e(hdl, "_ZN4llvm3sys16getProcessTripleB5cxx11Ev") != C_NULL
            return :cxx11
        else
            @warn("Unable to find llvm::sys::getProcessTriple() in libLLVM!")
            return nothing
        end
    end
end

function detect_compiler_abi(platform::Platform=default_platkey)
    return CompilerABI(;
        libgfortran_version=detect_libgfortran_version(;platform=platform),
        libstdcxx_version=detect_libstdcxx_version(),
        cxxstring_abi=detect_cxxstring_abi(),
    )
end

# Cache the default platform_key_abi() since that's by far the most common way
# we call platform_key_abi(), and we don't want to parse the same thing over
# and over and over again.  Note that we manually slap on a compiler abi
# string onto the end of Sys.MACHINE, like we expect our triplets to be encoded.
# Note futher that manually pass in an incomplete platform_key_abi() to the `detect_*()`
# calls, because we need to know things about dynamic library naming rules and whatnot.
default_platkey = platform_key_abi(string(
    Sys.MACHINE,
    compiler_abi_str(detect_compiler_abi(platform_key_abi(Sys.MACHINE))),
))
function platform_key_abi()
    global default_platkey
    return default_platkey
end

function platforms_match(a::Platform, b::Platform)
    # Check to see if a and b  satisfy the rigid constraints first, these are
    # things that are simple equality checks:
    function rigid_constraints(a, b)
        return (typeof(a) <: typeof(b) || typeof(b) <: typeof(a)) &&
               (arch(a) == arch(b)) && (libc(a) == libc(b)) &&
               (call_abi(a) == call_abi(b))
    end

    # The flexible constraints are ones that can do equals, but also have things
    # like "any" values, etc....
    function flexible_constraints(a, b)
        ac = compiler_abi(a)
        bc = compiler_abi(b)
        gcc_match = (ac.libgfortran_version === nothing
                  || bc.libgfortran_version === nothing
                  || ac.libgfortran_version == bc.libgfortran_version)
        cxx_match = (ac.cxxstring_abi === nothing
                  || bc.cxxstring_abi === nothing
                  || ac.cxxstring_abi == bc.cxxstring_abi)
        return gcc_match && cxx_match
    end

    return rigid_constraints(a, b) && flexible_constraints(a, b)
end

function platforms_match(a::AbstractString, b::Platform)
    @nospecialize b
    return platforms_match(platform_key_abi(a), b)
end
function platforms_match(a::Platform, b::AbstractString)
    @nospecialize a
    return platforms_match(a, platform_key_abi(b))
end
platforms_match(a::AbstractString, b::AbstractString) = platforms_match(platform_key_abi(a), platform_key_abi(b))

"""
    select_platform(download_info::Dict, platform::Platform = platform_key_abi())

Given a `download_info` dictionary mapping platforms to some value, choose
the value whose key best matches `platform`, returning `nothing` if no matches
can be found.

Platform attributes such as architecture, libc, calling ABI, etc... must all
match exactly, however attributes such as compiler ABI can have wildcards
within them such as `nothing` which matches any version of GCC.
"""
function select_platform(download_info::Dict, platform::Platform = platform_key_abi())
    @nospecialize platform
    ps = collect(filter(p -> platforms_match(p, platform), keys(download_info)))

    if isempty(ps)
        return nothing
    end

    # At this point, we may have multiple possibilities.  E.g. if, in the future,
    # Julia can be built without a direct dependency on libgfortran, we may match
    # multiple tarballs that vary only within their libgfortran ABI.  To narrow it
    # down, we just sort by triplet, then pick the last one.  This has the effect
    # of generally choosing the latest release (e.g. a `libgfortran5` tarball
    # rather than a `libgfortran3` tarball)
    p = last(sort(ps, by = p -> triplet(p)))
    return download_info[p]
end

end # module


text:
@testset "Boson measurements" begin
    copyto!(mc.p.hsfield, randconf)

    @test isapprox(measure_chi_dynamic(mc.p.hsfield), load("data/O3.jld", "chi_dyn"))
    @test isapprox(measure_chi_static(mc.p.hsfield), 12.420575691388407)
    @test all(isapprox.(measure_binder_factors(mc.p.hsfield), (0.7762859807117749, 0.6026199238496421)))
    a,b = measure_op(mc.p.hsfield)
    @test isapprox(a, 0.5080018386458331)
    @test isapprox(b, [0.514548, 0.536541, 0.472916], atol=1e-4)

    qy, qx, ws = get_momenta_and_frequencies(mc)
    @test qy == qx
    @test isapprox(qy, [0.0, 1.5708, 3.14159], atol=1e-4)
    @test isapprox(ws, [0.0, 6.28319, 12.5664, 18.8496, 25.1327, 31.4159], atol=1e-4)

    qy, qx, ws = get_momenta_and_frequencies(mc.p.hsfield)
    @test qy == qx
    @test isapprox(qy, [0.0, 1.5708, 3.14159], atol=1e-4)
    @test isapprox(ws, [0.0, 6.28319, 12.5664, 18.8496, 25.1327, 31.4159], atol=1e-4)
end

@testset "Fermion measurements" begin
    init!(mc)

    @testset "Helpers" begin
        @testset "Accessing Greens" begin
            sql = reshape(1:16, (4,4))
            N = mc.l.sites
            greens = mc.s.greens[1:2*N,1:2*N] # fake O(2) greens

            f = () -> begin
                for x in 1:4
                    for y in 1:4
                        siteidx(mc, sql, x, y) == sql[y,x] || return false
                    end
                end
                return true
            end
            @test f()

            # PBC
            @test siteidx(mc, sql, 5, 1) == 1
            @test siteidx(mc, sql, 1, 5) == 1
            @test siteidx(mc, sql, 5, 5) == 1
            @test siteidx(mc, sql, 3, 7) == 11
            @test siteidx(mc, sql, 128, 32) == 16

            # greensidx
            f = () -> begin
                for i in 1:N
                    greensidx(N, 1, i) == i || return false
                    greensidx(N, 2, i) == N+i || return false
                    greensidx(N, 3, i) == 2*N+i || return false
                    greensidx(N, 4, i) == 3*N+i || return false
                end
                return true
            end
            @test f()

            # G for O(3) equal to indexing greens
            f = () -> begin
                for i in 1:4*N, j in 1:4*N
                    G(mc, i, j) == mc.s.greens[i, j] || return false
                end
                return true
            end
            @test f()

            # Gtilde for O(3) equal to indexing I - greens
            f = () -> begin
                for i in 1:4*N, j in 1:4*N
                    Gtilde(mc, i, j) == kd(i,j) - mc.s.greens[i, j] || return false
                end
                return true
            end
            @test f()

            @test isapprox(fullG(mc, mc.s.greens), mc.s.greens)
            @test isapprox(fullGtilde(mc, mc.s.greens), I - mc.s.greens)

            mc.p.opdim = 2

            # G for O(2)/O(1)
            f = () -> begin
                for i in 1:2*N, j in 1:2*N
                    isapprox(G(mc, 2*N+i, 2*N+j, greens), conj(greens[i, j])) || return false
                    isapprox(Gtilde(mc, 2*N+i, 2*N+j, greens), kd(2*N+i, 2*N+j) - conj(greens[i, j])) || return false
                end
                return true
            end
            @test f()

            @test isapprox(fullG(mc, greens), I - fullGtilde(mc, greens))

            mc.p.opdim = 3
        end

        @testset "Permute and FFT Greens" begin
            N = mc.l.sites
            L = mc.p.L
            xu, yd, xd, yu = 1, 2, 3, 4 # definition = old order

            flvtrafo = Dict{Int64, Int64}(xu => xu, yd => yu, xd => xd, yu => yd)

            @test isapprox(permute_greens(permute_greens(mc.s.greens)), mc.s.greens)

            gperm = PseudoBlockArray(permute_greens(mc.s.greens), [N,N,N,N], [N,N,N,N])
            g = PseudoBlockArray(mc.s.greens, [N,N,N,N], [N,N,N,N])

            f = () -> begin
                for j in (xu, yd, xd, yu)
                    for i in (xu, yd, xd, yu)
                        isapprox(getblock(gperm, i, j), getblock(g, flvtrafo[i], flvtrafo[j])) || return false
                    end
                end
                return true
            end
            @test f()

            # FFT
            @test isapprox(ifft_greens(mc, fft_greens(mc, mc.s.greens; fftshift=false); ifftshift=false), mc.s.greens)
            @test isapprox(ifft_greens(mc, fft_greens(mc, mc.s.greens; fftshift=true); ifftshift=true), mc.s.greens)

            gk = fft_greens(mc, mc.s.greens)
            @test isapprox(gk, load("data/O3.jld", "greens_fft"))

            qs = collect(range(-pi, pi, length=L+1))[1:end-1]
            qs7 = [-2.6927937030769655, -1.7951958020513104, -0.8975979010256552, 0.0, 0.8975979010256552, 1.7951958020513104, 2.6927937030769655]
            @test all(isapprox.(fftmomenta(7; fftshift=true), qs7))
            @test all(isapprox.(fftmomenta(mc; fftshift=true), qs))
            @test all(isapprox.(fftmomenta(7; fftshift=false), ifftshift(qs7)))
            @test all(isapprox.(fftmomenta(mc; fftshift=false), ifftshift(qs)))
        end
    end

    # @testset "Effective -> Actual GF" begin
    #     # TODO
    # end

    # Calculate ETGF and logdet
    g = effective_greens2greens(mc, mc.s.greens)
    lgdet = calculate_logdet(mc)

    # Calculate TDGFs
    allocate_tdgfs!(mc)
    measure_tdgfs!(mc)
    Gt0 = mc.s.meas.Gt0
    G0t = mc.s.meas.G0t

    @testset "ETGF(s)" begin
        gmeas = measure_greens(mc, mc.p.safe_mult, mc.s.current_slice);
        gmeas2, lgdetmeas = measure_greens_and_logdet(mc, mc.p.safe_mult, mc.s.current_slice);

        @test isapprox(g, gmeas)
        @test isapprox(gmeas, gmeas2)
        @test isapprox(lgdet, lgdetmeas)

        # ETGFs at all time slices
        a = calc_all_greens_explicit(mc);
        b = calc_all_greens(mc);

        @test maximum(maximum.(real.(a .- b))) < 1e-12
        @test maximum(maximum.(imag.(a .- b))) < 1e-12
    end

    @testset "Occupation" begin
        @test isapprox(occupation(mc, g), 0.4216, atol=1e-4)
        @test all(isapprox.(occupations_flv(mc, g), 0.4216, atol=1e-4))
    end

    @testset "TDGFs" begin
        @test isapprox(calc_tdgf_direct(mc, 1), Gt0[1])
        # @test_broken isapprox(calc_tdgf_direct(mc, 4), Gt0[4])

        check_unitarity = (u_stack) -> begin
          for i in 1:length(u_stack)
            U = u_stack[i]
            !isapprox(U * adjoint(U), I) && (return false)
          end
          return true
        end

        @test check_unitarity(mc.s.meas.BT0Inv_u_stack)
        @test check_unitarity(mc.s.meas.BBetaT_u_stack)
        @test check_unitarity(mc.s.meas.BT0_u_stack)
        @test check_unitarity(mc.s.meas.BBetaTInv_u_stack)

        # TODO: Check stacks values?

        @test isapprox(estimate_memory_usage_tdgfs(mc), 1.8, atol=1e-1)
        @test isapprox(memory_usage_tdgfs(mc), 1.8, atol=1e-1)
    end

    @testset "ETPC" begin
        allocate_etpc!(mc)
        @test mc.s.meas.etpc_minus == zero(mc.s.meas.etpc_minus)
        @test mc.s.meas.etpc_plus == zero(mc.s.meas.etpc_plus)

        measure_etpc!(mc, g)
        Pm = mc.s.meas.etpc_minus
        Pp = mc.s.meas.etpc_plus

        @test isreal(Pm)
        @test isreal(Pp)
        @test is_reflection_symmetric(Pm, tol=1e-3)
        @test is_reflection_symmetric(Pp, tol=1e-3)
    end

    @testset "ETCDC" begin
        allocate_etcdc!(mc)
        @test mc.s.meas.etcdc_minus == zero(mc.s.meas.etcdc_minus)
        @test mc.s.meas.etcdc_plus == zero(mc.s.meas.etcdc_plus)

        measure_etcdc!(mc, g)
        Cm = mc.s.meas.etcdc_minus
        Cp = mc.s.meas.etcdc_plus

        @test maximum(imag(Cm)) < 1e-12
        @test maximum(imag(Cp)) < 1e-12
        @test is_reflection_symmetric(Cm, tol=1e-3)
        @test is_reflection_symmetric(Cp, tol=1e-3)
    end

    @testset "ZFPC" begin
        allocate_zfpc!(mc)
        @test mc.s.meas.zfpc_minus == zero(mc.s.meas.zfpc_minus)
        @test mc.s.meas.zfpc_plus == zero(mc.s.meas.zfpc_plus)

        measure_zfpc!(mc, Gt0)
        Pm = mc.s.meas.zfpc_minus
        Pp = mc.s.meas.zfpc_plus

        @test isreal(Pm)
        @test isreal(Pp)
        @test is_reflection_symmetric(Pm, tol=1e-3)
        @test is_reflection_symmetric(Pp, tol=1e-3)

        M = mc.p.slices
        measure_zfpc!(mc, [g for _ in 1:M]) # ~ τ = 0
        @test Pm/M ≈ mc.s.meas.etpc_minus
        @test Pp/M ≈ mc.s.meas.etpc_plus
    end

    @testset "ZFCDC" begin
        allocate_zfcdc!(mc)
        @test mc.s.meas.zfcdc_minus == zero(mc.s.meas.zfcdc_minus)
        @test mc.s.meas.zfcdc_plus == zero(mc.s.meas.zfcdc_plus)

        measure_zfcdc!(mc, g, Gt0, G0t)
        Cm = mc.s.meas.zfcdc_minus
        Cp = mc.s.meas.zfcdc_plus

        @test maximum(imag(Cm)) < 1e-12
        @test maximum(imag(Cp)) < 1e-12
        @test is_reflection_symmetric(Cm, tol=1e-3)
        @test is_reflection_symmetric(Cp, tol=2e-3) # is 2e-3 a problem?
    end

    @testset "ZFCCC + Superfluid density" begin
        allocate_zfccc!(mc)
        @test mc.s.meas.zfccc_xx == zero(mc.s.meas.zfccc_xx)
        @test mc.s.meas.zfccc_yy == zero(mc.s.meas.zfccc_yy)

        measure_zfccc!(mc, g, Gt0, G0t)
        Λxx = mc.s.meas.zfccc_xx
        Λyy = mc.s.meas.zfccc_yy

        @test maximum(imag(Λxx)) < 1e-12
        @test is_reflection_symmetric(Λxx, tol=2e-3) # is 2e-3 a problem?
        @test maximum(imag(Λyy)) < 1e-12
        # TODO: For our model Λxx ≈ Λyy'. perhaps we should test that?
        @test is_reflection_symmetric(Λyy, tol=1e-2) # is 1e-2 a problem?

        @test isapprox(measure_sfdensity(mc, Λxx, Λyy), 0.05157097614282691) # This value is only self-consistent.
    end

    @test isnothing(deallocate_tdgfs_stacks!(mc))
    @test length(mc.s.meas.BT0Inv_u_stack) == 0
    @test length(mc.s.meas.BT0Inv_d_stack) == 0
    @test length(mc.s.meas.BT0Inv_t_stack) == 0
    @test length(mc.s.meas.BBetaT_u_stack) == 0
    @test length(mc.s.meas.BBetaT_d_stack) == 0
    @test length(mc.s.meas.BBetaT_t_stack) == 0
    @test length(mc.s.meas.BT0_u_stack) == 0
    @test length(mc.s.meas.BT0_d_stack) == 0
    @test length(mc.s.meas.BT0_t_stack) == 0
    @test length(mc.s.meas.BBetaTInv_u_stack) == 0
    @test length(mc.s.meas.BBetaTInv_d_stack) == 0
    @test length(mc.s.meas.BBetaTInv_t_stack) == 0

end

nothing

ext: jl

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
@testset "Database IO" begin

    @testset "SQLite" begin
        fname = joinpath(testpath, "test.db")
        db = SQLite.DB(fname)

        write_table!(db, "test1", df)
        @test filesize(fname) > 0
        df_recovered = DataFrame(read_table(fname, "test1"); copycols=false)
        @test df == df_recovered

        df_sql = DataFrame(read_sql(db, "select * from test1 where a < 5"); copycols=false)
        @test df[df.a .< 5, :] == df_sql

        write_table!(fname, "test2", nt)
        nt_recovered = read_table(db, "test2")
        @test DataFrame(nt) == DataFrame(nt_recovered)
    end

    @testset "PostgreSQL" begin
        # the following tests require a running PostgreSQL database.
        # `docker run --rm --detach --name test-libpqjl -e POSTGRES_HOST_AUTH_METHOD=trust -p 5432:5432 postgres`
        conn = LibPQ.Connection("dbname=postgres user=postgres")

        execute(conn, """CREATE TEMPORARY TABLE test1 (
            a integer PRIMARY KEY,
            b numeric,
            c character varying,
            d boolean,
            e date,
            f character varying
            );""")
        write_table!(conn, "test1", df)
        df_recovered = DataFrame(read_table(conn, "test1"); copycols=false)
        @test df == df_recovered

        df_sql = DataFrame(read_sql(conn, "select * from test1 where a < 5"); copycols=false)
        @test df[df.a .< 5, :] == df_sql

        execute(conn, """CREATE TEMPORARY TABLE test2 (
            a integer PRIMARY KEY,
            b numeric,
            c character varying
            );""")
        write_table!(conn, "test2", nt)
        nt_recovered = read_table(conn, "test2")
        @test DataFrame(nt) == DataFrame(nt_recovered)

        close(conn)
    end
end


text:
"""
Write-only Unix File Descriptor.
"""

raw_transfer(fd::FD{Out}, buf, count) = C.write(fd, buf, count)

shutdown(fd::FD{Out}) = shutdown(fd, C.SHUT_WR)

@static if isdefined(Base, :shutdown)
    Base.shutdown(fd::FD{Out}) = UnixIO.shutdown(fd, C.SHUT_WR)
end

Base.iswritable(fd::FD{Out}) = Base.isopen(fd)

Base.isreadable(::FD{Out}) = false

@db 1 function Base.unsafe_write(fd::FD{Out}, buf::Ptr{UInt8}, nbytes::UInt)
    @require !fd.isclosed
    nwritten = 0                   ;if nbytes < 100
                                        @db 1 repr(unsafe_string(buf, nbytes))
                                    end
    while nwritten < nbytes
        n = transfer(fd, buf + nwritten, nbytes - nwritten)
        nwritten += n
    end
    @ensure nwritten == nbytes
    @db 1 return Int(nwritten)
end

Base.write(fd::FD{Out}, x::UInt8) = transfer(fd, [x], 1, 1)

# End of file: WriteFD.jl


text:
struct PgramData <: JAXTAMData
    mission    :: Mission
    instrument :: Symbol
    obsid      :: String
    e_range    :: Tuple{Float64,Float64}
    bin_time   :: Real
    pg_type    :: Symbol
    power      :: Array
    freq       :: Array
    group      :: Int
end

function _pgram(counts, times, bin_time, pg_type=:standard; maximum_frequency=:auto)
    if maximum_frequency == :auto
        maximum_frequency = 0.5/bin_time
    elseif maximum_frequency > 0.5/bin_time
        @warn "Maximum frequency cannot be greater than 0.5/bin_time ($(0.5/bin_time))"
        maximum_frequency = 0.5/bin_time
    end

    pg_plan = LombScargle.plan(
        times, float(counts),
        normalization=pg_type,
        minimum_frequency=1/(length(counts)*bin_time), # TODO: Check if setting min freq to 1/exposure is correct
        maximum_frequency=maximum_frequency,
        # samples_per_peak=1,
    )

    pg = lombscargle(pg_plan)

    freq, power = freqpower(pg)

    return freq, power
end

function _pgram(lc::BinnedData, pg_type, group; maximum_frequency=:auto)
    freq, power = (missing, missing)
    if group == 0 # Using whole lightcurve
        lc_groups = _group_return(lc; min_length_sec=16)
        lc_group_counts = vcat([lc[2].counts for lc in lc_groups]...)
        lc_group_times  = vcat([lc[2].times for lc in lc_groups]...)
        freq, power = _pgram(lc_group_counts, lc_group_times, lc.bin_time, pg_type; maximum_frequency=maximum_frequency)
    else
        freq, power = _pgram(lc.counts, lc.times, lc.bin_time, pg_type; maximum_frequency=maximum_frequency)
    end

    return PgramData(lc.mission, lc.instrument, lc.obsid, lc.e_range, lc.bin_time,
        pg_type, power, freq, group)
end

function _pgram_lc_group_pad(group_lc::Dict{Int64,JAXTAM.BinnedData})
    groups = keys(group_lc)
    longest_group = maximum([length(x.counts) for x in values(group_lc)])

    padded_groups = Dict{Int64,JAXTAM.BinnedData}()

    for group in groups
        mission    = group_lc[group].mission
        instrument = group_lc[group].instrument
        obsid      = group_lc[group].obsid
        bin_time   = group_lc[group].bin_time
        counts     = group_lc[group].counts
        times      = group_lc[group].times
        gtis       = group_lc[group].gtis

        if length(counts) < longest_group
            counts = [counts; zeros(longest_group-length(counts))]
        end

        padded_groups[group] = BinnedData(mission, instrument, obsid, bin_time, counts, times, gtis)
    end

    return padded_groups
end

function pgram(instrument_lc::Dict{Symbol,BinnedData}; pg_type=:standard, per_group=false, maximum_frequency=:auto)
    instruments = keys(instrument_lc)

    instrument_pgram = per_group ? Dict{Symbol,Dict{Int,PgramData}}() : Dict{Symbol,PgramData}()
    for instrument in instruments
        if per_group
            group_pgram = Dict{Int,PgramData}()

            lc_groups = _group_return(instrument_lc[instrument])
            #lc_groups = _pgram_lc_group_pad(lc_groups)
            groups = keys(lc_groups)

            for group in groups
                group_pgram[group] = _pgram(lc_groups[group], pg_type, group, maximum_frequency=maximum_frequency)
            end

            #group_pgram[-1] = PgramData(group_pgram[1].mission, group_pgram[1].instrument, group_pgram[1].obsid,
            #    group_pgram[1].bin_time, pg_type, mean([pgram.power for pgram in values(group_pgram)]), group_pgram[1].freq, -1)

            instrument_pgram[instrument] = group_pgram
        else
            instrument_pgram[instrument] = _pgram(instrument_lc[instrument], pg_type, 0, maximum_frequency=maximum_frequency)
        end
    end

    return instrument_pgram
end


text:
"""
Tools for the exploration of chaos and nonlinear dynamics
"""
module ChaosTools

using Reexport
@reexport using DynamicalSystemsBase

using DynamicalSystemsBase: DS, DDS, CDS
using DynamicalSystemsBase: MDI, TDI
using DynamicalSystemsBase: stateeltype
using DynamicalSystemsBase.DiffEqBase: AbstractODEIntegrator, u_modified!, DEIntegrator
DEI = DEIntegrator

include("orbitdiagrams/discrete_diagram.jl")
include("orbitdiagrams/poincare.jl")

include("dimensions/entropies.jl")
include("dimensions/dims.jl")
include("dimensions/correlationdim.jl")

include("nlts.jl")

include("period_return/periodic_points.jl")
include("period_return/period.jl")
include("period_return/transit_times_statistics.jl")

include("chaosdetection/lyapunovs.jl")
include("chaosdetection/gali.jl")
include("chaosdetection/expansionentropy.jl")
include("chaosdetection/partially_predictable.jl")
include("chaosdetection/testchaos01.jl")

# Ugly methods that shouldn't exist:
include("ugliness.jl")

display_update = false
version = "1.10.0"
update_name = "update_v$version"

if display_update
if !isfile(joinpath(@__DIR__, update_name))
printstyled(stdout,
"""
\nUpdate message: ChaosTools v$version

A method to calculate the expansion entropy for discrete and
continuous systems is now included as `expansionentropy`!

See B. Hunt & E. Ott, ‘Defining Chaos’, Chaos 25.9 (2015).
\n
"""; color = :light_magenta)
touch(joinpath(@__DIR__, update_name))
end
end

end # module


text:
### Test error exits
mf = ModelFrame(Yield ~ 1 + (1|Batch), ds)
Λ = [LowerTriangular(eye(1)) for i = 1]
y = convert(Vector, ds[:Yield])
Rem = push!([], remat(:(1|Batch), ds))
X = ones(30,1)

@test_throws ArgumentError LinearMixedModel(mf, ones(1), Λ, X, y, Float64[])
@test_throws DimensionMismatch LinearMixedModel(mf, Rem, Λ, X, y[1:20], Float64[])
@test_throws DimensionMismatch LinearMixedModel(mf, Rem, Λ, X[1:20, :], y, Float64[])
@test_throws DimensionMismatch LinearMixedModel(mf, append!(copy(Rem), Rem), Λ, X, y, Float64[])
@test_throws DimensionMismatch LinearMixedModel(mf, Rem, append!(copy(Λ), Λ), X, y, Float64[])
@test_throws DimensionMismatch LinearMixedModel(mf,Rem, [LowerTriangular(eye(2)) for i in 1], X, y, Float64[])
@test_throws DomainError LinearMixedModel(mf, Rem, Λ, X, y, fill(-1., 30))

@test_throws ArgumentError lmm(Yield ~ 1, ds)

modl = lmm(Yield ~ 1 + (1|Batch), ds);

@test_throws ArgumentError modl[:ϕ] = [1.]
@test_throws DimensionMismatch modl[:θ] = [0., 1.]

@test_throws DimensionMismatch Base.LinAlg.Ac_ldiv_B!(UpperTriangular(fm3.R[1, 1]), ones(30, 1))

const tri = LowerTriangular(eye(3))
const hblk = MixedModels.HBlkDiag(ones(2,3,2))
const dd3 = Diagonal(ones(3))
@test_throws KeyError tri[:foo]
@test_throws KeyError tri[:foo] = ones(3)
@test_throws DimensionMismatch tri[:θ] = ones(2)
@test_throws DimensionMismatch MixedModels.tscale!(tri, ones(2))
@test_throws DimensionMismatch MixedModels.tscale!(ones(2,2), tri)
@test_throws DimensionMismatch MixedModels.tscale!(tri, hblk)
@test_throws ArgumentError Base.cholfact!(hblk)
@test_throws DimensionMismatch MixedModels.downdate!(dd3, Diagonal(ones(2)))
@test_throws DimensionMismatch MixedModels.downdate!(dd3, sparse(ones(4,1)))
@test_throws DimensionMismatch MixedModels.downdate!(ones(2,2), dd3, ones(3,3))
@test_throws DimensionMismatch MixedModels.downdate!(ones(2,2), dd3, ones(2,2))
@test_throws DimensionMismatch MixedModels.inject!(ones(2,2), dd3)
@test_throws DimensionMismatch MixedModels.inject!(ones(3,2), dd3)
@test_throws DimensionMismatch Base.LinAlg.A_ldiv_B!(dd3, Diagonal(ones(2)))

const speye4 = speye(4)
@test_throws ArgumentError MixedModels.inject!(speye(4), sparse(ones(4,4)))
@test_throws DimensionMismatch Base.LinAlg.A_ldiv_B!(dd3, speye4)
@test_throws DimensionMismatch Base.LinAlg.Ac_ldiv_B!(UpperTriangular(hblk), ones(1,1))
@test_throws DimensionMismatch MixedModels.unscaledre!(ones(1), modl.wttrms[1], ones(1, 6))

@test_throws DimensionMismatch MixedModels.downdate!(ones(1,1), speye4, ones(4,1))
@test_throws DimensionMismatch MixedModels.downdate!(ones(1,1), speye4, speye4)
@test_throws DimensionMismatch MixedModels.downdate!(ones(1,1), speye(3), speye4)
@test_throws DimensionMismatch MixedModels.downdate!(ones(1,1), speye4)

@test_throws DimensionMismatch MixedModels.inject!(speye4, speye(3))

@test_throws DimensionMismatch MixedModels.unscaledre!(zeros(size(slp, 1)),
  remat(:(1 + Days | Subject), slp), ones(1,30))

@test_throws DimensionMismatch MixedModels.tscale!(LowerTriangular(eye(2)), Diagonal(ones(30)))
@test_throws DimensionMismatch MixedModels.tscale!(Diagonal(ones(30)), LowerTriangular(eye(2)))

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
#
# Potential energy function using r2
#

function upair3(r2, p1, p2, data)
    dx = p1[1] - p2[1]
    dy = p1[2] - p2[2]
    r6 = r2^3
    r8 = r6*r2
    r12 = r6^2
    r14 = r12*r2
    u = data.ε4*(data.σ12/r12 - data.σ6/r6)
    f = data.ε24*(2*data.σ12/r14 - data.σ6/r8) * [dx, dy]
    return u, f
end


text:
using IonBase

const test_dir = joinpath(pkgdir(IonBase), "test")
const test_ion_dir = joinpath(test_dir, ".ion")
const test_template = joinpath(test_dir, "acdemic.toml")

function with_test_ion(f)
    cd(test_dir) do
        withenv(f, "DOT_ION_PATH"=>test_ion_dir, "COMONICON_DEBUG"=>"ON")
    end
end

with_test_ion() do
    IonBase.copy_assets()
end


text:
if conf["is_skeletal"]
    dir_master = conf["master"]
    @load string(ckpt_path, "/../../", dir_master, "/conds.bson") conds
    n_exp = size(conds)[1]
else
    n_exp = Int64(conf["n_exp"])
    Random.seed!(1);
    plan = randomLHC(n_exp, 4) ./ n_exp

    conds = scaleLHC(plan,[(Float64(conf["temperature"]["lb"]),
                            Float64(conf["temperature"]["ub"])),
                           (Float64(conf["pressure"]["lb"]),
                            Float64(conf["pressure"]["ub"])),
                           (Float64(conf["phi"]["lb"]),
                            Float64(conf["phi"]["ub"])), (0.0, 1.0)])

    for i_exp in 1:n_exp
        local T0 = conds[i_exp, 1]
        global P = conds[i_exp, 2] * one_atm
        local phi = conds[i_exp, 3]
        conds[i_exp, 4] = f_idt(T0, P, phi, p; dT=dT, tfinal=10.0)
        if i_exp < 10
            check_sol(T0, P, phi, p; i_exp=i_exp)
        end
        @printf("(%d) T %.1f [K] P %.1f [atm] phi %.1f IDT0 = %.2e [s] \n",
                i_exp, conds[i_exp, 1], conds[i_exp, 2],
                conds[i_exp, 3], conds[i_exp, 4])
    end
    @save string(ckpt_path, "/../conds.bson") conds
end

n_train = Int64(n_exp * 0.7)

max_stars_repo_name: DENG-MIT/DeepReduction

max_issues_repo_name: DENG-MIT/DeepReduction

max_forks_repo_name: DENG-MIT/DeepReduction


text:
__precompile__()

module MUMPSjInv

using LinearAlgebra
using SparseArrays
using Distributed
using Printf
using Libdl
import jInv.LinearSolvers

mutable struct MUMPSfactorization{T}
	ptr::Int64     # pointer to factorization
	worker::Int64  # id of worker that holds factorization
	n::Int64       # matrix size
	a11::T         # first element (HACK for generating the right parametric type)
	time::Float64  # factorization time
end
	const MUMPSlibPath  = abspath(joinpath(splitdir(Base.source_path())[1],"..","lib","MUMPS"))
	const hasMUMPS = find_library([MUMPSlibPath])!=""
	println("MUMPSlibPath = ",MUMPSlibPath)
	println("hasMUMPS = ",hasMUMPS)

	include("MUMPSfuncs.jl")
	include("mumpsWrapper.jl")

	export solveMUMPS,solveMUMPS!, factorMUMPS, applyMUMPS,applyMUMPS!,destroyMUMPS, MUMPSfactorization

end


text:
#=
Project:      MIRTdev
File:         mri_geom_compress.jl
Description:  MRI Coil Compression
Author:       Rajas Gupta
Date Created: October 8th 2019
Date Update:  October 8th 2019
=#

export ir_mri_coil_compress

using LinearAlgebra: svd, norm

"""
    (odata, σ, Vr) = ir_mri_coil_compress(idata ; ncoil)

MRI coil compression via PCA.
Given multiple MRI surface coil images (idata), use SVD/PCA
to find a smaller number of virtual coil images (odata).

In:
* `idata` `[(N) n_in]`: noisy complex images (2D or 3D) for each coil

Option:
* `ncoil` Desired # of virtual coils (default: 1)

Out:
* `odata` `[(N) ncoil]`: virtual coil images
* `σ`     `[n_in]`: singular values.
* `Vr`    `[n_in, ncoil]`: compression matrix for reducing other data.

todo: currently ignores noise correlations

"""
function ir_mri_coil_compress(idata::AbstractArray{<:Number} ; ncoil::Int = 1)

    idim = size(idata)
    n_in = idim[end]
    idata = reshape(idata, :, n_in)

    (_, σ, V) = svd(idata)

    Vr = V[:, 1:ncoil]
    odata = idata * Vr
    odata = reshape(odata, idim[1:end-1]..., ncoil)
    return (odata, σ, Vr)
end

"""
    snr2sigma(db, yb)
convert SNR in dB to noise σ for complex gaussian noise
"""
function snr2sigma(db, yb::AbstractArray{<:Complex})
    10^(-db/20) * norm(yb) / sqrt(length(yb)) / sqrt(2)
end


text:
##
##  l a . j l  Linear algebra routines
##

function trisolve(d0::Vector{T},
                             d1::Vector{T},d2::Vector{T},rhs::Vector{T}) where {T<:Number}
    local n = length(d0)
    if n < 3
        error("Length of 'a' must be greater or equal to 3.")
    end
    if length(d1) != n-1 || length(d2) != n-1
        error("Vectors 'dup' and 'dlo' must be of a length of length(a)-1.")
    end
    if length(rhs) != n
        error("Vector 'rhs' must be of the same length as 'a'.")
    end

    local x = copy(rhs); local a = copy(d0); local b = [d1, 0]; local d = copy(d2)
    for i = 1:(n-1)
        if d[i] != 0
            t = a[i]/d[i]
            si = 1/sqrt(1+t*t)
            co = t*si
            a[i] = a[i]*co + d[i]*si

            h = b[i]
            b[i] = h*co + a[i+1]*si
            a[i+1] = -h*si + a[i+1]*co
            d[i] = b[i+1]*si
            b[i+1] = b[i+1]*co

            h = x[i]
            x[i] = h*co + x[i+1]*si
            x[i+1] = -h*si + x[i+1]*co
        end
    end

    if any(a .== 0.0)
        error("Triangular matrix is singular---system not solvable.")
    end

    x[n] = x[n]/a[n]
    x[n-1] = ( x[n-1] - b[n-1]*x[n] ) / a[n-1]
    for i = (n-2):-1:1
        x[i] = (x[i] - b[i]*x[i+1] - d[i]*x[i+2]) / a[i]
    end

    return x
end


text:
using ITensors, Test, Random

@testset "Basic DMRG" begin

  @testset "Spin-one Heisenberg" begin
    N = 10
    sites = siteinds("S=1",N)

    ampo = AutoMPO()
    for j=1:N-1
      add!(ampo,"Sz",j,"Sz",j+1)
      add!(ampo,0.5,"S+",j,"S-",j+1)
      add!(ampo,0.5,"S-",j,"S+",j+1)
    end
    H = MPO(ampo,sites)

    psi = randomMPS(sites)

    sweeps = Sweeps(3)
    @test length(sweeps) == 3
    maxdim!(sweeps,10,20,40)
    mindim!(sweeps,1,10)
    cutoff!(sweeps,1E-11)
    noise!(sweeps,1E-10)
    str = split(sprint(show, sweeps), '\n')
    @test length(str) > 1
    energy,psi = dmrg(H, psi, sweeps; outputlevel=0)
    @test energy < -12.0
  end

  @testset "QN-conserving Spin-one Heisenberg" begin
    N = 10
    sites = siteinds("S=1",N; conserve_qns=true)

    ampo = AutoMPO()
    for j=1:N-1
      ampo += "Sz",j,"Sz",j+1
      ampo += 0.5,"S+",j,"S-",j+1
      ampo += 0.5,"S-",j,"S+",j+1
    end
    H = MPO(ampo,sites)

    state = [isodd(n) ? "Up" : "Dn" for n in 1:N]
    psi = randomMPS(sites,state,4)

    sweeps = Sweeps(3)
    @test length(sweeps) == 3
    maxdim!(sweeps,10,20,40)
    mindim!(sweeps,1,10)
    cutoff!(sweeps,1E-11)
    noise!(sweeps,1E-10)
    str = split(sprint(show, sweeps), '\n')
    @test length(str) > 1
    energy,psi = dmrg(H, psi, sweeps; outputlevel=0)
    @test energy < -12.0
  end

  @testset "Transverse field Ising" begin
    N = 32
    sites = siteinds("S=1/2",N)
    Random.seed!(432)
    psi0 = randomMPS(sites)

    ampo = AutoMPO()
    for j = 1:N
      j < N && add!(ampo,-1.0,"Sz",j,"Sz",j+1)
      add!(ampo,-0.5,"Sx",j)
    end
    H = MPO(ampo,sites)

    sweeps = Sweeps(5)
    maxdim!(sweeps,10,20)
    cutoff!(sweeps,1E-12)
    noise!(sweeps,1E-10)
    energy,psi = dmrg(H,psi0,sweeps,outputlevel=0)

    # Exact energy for transverse field Ising model
    # with open boundary conditions at criticality
    energy_exact = 0.25 - 0.25/sin(π/(2*(2*N+1)))
    @test abs((energy-energy_exact)/energy_exact) < 1e-4
  end

  @testset "Transverse field Ising, conserve Sz parity" begin
    N = 32
    sites = siteinds("S=1/2", N, conserve_szparity = true)
    Random.seed!(432)

    state = [isodd(j) ? "↑" : "↓" for j in 1:N]
    psi0 = randomMPS(sites, state)

    ampo = AutoMPO()
    for j = 1:N
      j < N && add!(ampo,-1.0,"Sx",j,"Sx",j+1)
      add!(ampo,-0.5,"Sz",j)
    end
    H = MPO(ampo,sites)

    sweeps = Sweeps(5)
    maxdim!(sweeps,10,20)
    cutoff!(sweeps,1E-12)
    noise!(sweeps,1E-10)
    energy,psi = dmrg(H,psi0,sweeps,outputlevel=0)

    # Exact energy for transverse field Ising model
    # with open boundary conditions at criticality
    energy_exact = 0.25 - 0.25/sin(π/(2*(2*N+1)))
    @test abs((energy-energy_exact)/energy_exact) < 1e-4
  end

  @testset "DMRGObserver" begin
    N = 10
    sites = siteinds("S=1/2",N)
    Random.seed!(42)
    psi0 = randomMPS(sites)

    ampo = AutoMPO()
    for j = 1:N-1
      ampo += -1,"Sz",j,"Sz",j+1
    end
    for j = 1:N
      ampo += -0.2,"Sx",j
    end
    H = MPO(ampo,sites)

    sweeps = Sweeps(3)
    maxdim!(sweeps,10)
    cutoff!(sweeps,1E-12)

    observer = DMRGObserver(["Sz","Sx"], sites)

    E,psi = dmrg(H,psi0,sweeps,observer=observer,outputlevel=0)
    @test length(measurements(observer)["Sz"])==3
    @test length(measurements(observer)["Sx"])==3
    @test all(length.(measurements(observer)["Sz"]) .== N)
    @test all(length.(measurements(observer)["Sx"]) .== N)
    @test length(energies(observer))==3
    @test length(truncerrors(observer))==3
    @test energies(observer)[end]==E
    @test all(truncerrors(observer) .< 1E-9)

    orthogonalize!(psi,1)
    m = scalar(dag(psi[1])*noprime(op(sites, "Sz", 1)*psi[1]))
    @test measurements(observer)["Sz"][end][1] ≈ m
  end

  @testset "Sum of MPOs (ProjMPOSum)" begin
    N = 10
    sites = siteinds("S=1",N)

    ampoZ = AutoMPO()
    for j=1:N-1
      ampoZ += "Sz",j,"Sz",j+1
    end
    HZ = MPO(ampoZ,sites)

    ampoXY = AutoMPO()
    for j=1:N-1
      ampoXY += 0.5,"S+",j,"S-",j+1
      ampoXY += 0.5,"S-",j,"S+",j+1
    end
    HXY = MPO(ampoXY,sites)

    psi = randomMPS(sites)

    sweeps = Sweeps(3)
    maxdim!(sweeps,10,20,40)
    mindim!(sweeps,1,10,10)
    cutoff!(sweeps,1E-11)
    noise!(sweeps,1E-10)
    energy,psi = dmrg([HZ,HXY], psi, sweeps; outputlevel=0)
    @test energy < -12.0
  end

  @testset "Excited-state DMRG" begin
    N = 10
    weight = 15.0

    sites = siteinds("S=1",N)
    sites[1] = Index(2,"S=1/2,n=1,Site")
    sites[N] = Index(2,"S=1/2,n=$N,Site")

    ampo = AutoMPO()
    for j=1:N-1
      ampo += "Sz",j,"Sz",j+1
      ampo += 0.5,"S+",j,"S-",j+1
      ampo += 0.5,"S-",j,"S+",j+1
    end
    H = MPO(ampo,sites)

    psi0i = randomMPS(sites,10)

    sweeps = Sweeps(4)
    maxdim!(sweeps, 10,20,100,100)
    cutoff!(sweeps, 1E-11)
    noise!(sweeps,1E-10)

    energy0, psi0 = dmrg(H,psi0i, sweeps; outputlevel=0)
    @test energy0 < -11.5

    psi1i = randomMPS(sites,10)
    energy1,psi1 = dmrg(H,[psi0],psi1i,sweeps;outputlevel=0,weight=weight)

    @test energy1 > energy0
    @test energy1 < -11.1

    @test inner(psi1,psi0) < 1E-5
  end

  @testset "Fermionic Hamiltonian" begin
    N = 10
    t1 = 1.0
    t2 = 0.5
    V = 0.2
    s = siteinds("Fermion", N; conserve_qns = true)

    state = fill(1,N)
    state[1] = 2
    state[3] = 2
    state[5] = 2
    state[7] = 2
    psi0 = productMPS(s,state)

    ampo = AutoMPO()
    for j=1:N-1
      ampo += -t1,"Cdag",j,"C",j+1
      ampo += -t1,"Cdag",j+1,"C",j
      ampo += V,"N",j,"N",j+1
    end
    for j=1:N-2
      ampo += -t2,"Cdag",j,"C",j+2
      ampo += -t2,"Cdag",j+2,"C",j
    end
    H = MPO(ampo, s)

    sweeps = Sweeps(5)
    maxdim!(sweeps, 10, 20, 100, 100, 200)
    cutoff!(sweeps, 1E-8)
    noise!(sweeps, 1E-10)

    energy, psi = dmrg(H, psi0, sweeps; outputlevel = 0)
    @test (-6.5 < energy < -6.4)
  end

  @testset "Hubbard model" begin
    N = 10
    Npart = 8
    t1 = 1.0
    U  = 1.0
    V1 = 0.5
    sites = siteinds("Electron",N; conserve_qns=true)
    ampo = AutoMPO()
    for i=1:N
      ampo += (U,"Nupdn",i)
    end
    for b=1:N-1
      ampo += -t1,"Cdagup",b,"Cup",b+1
      ampo += -t1,"Cdagup",b+1,"Cup",b
      ampo += -t1,"Cdagdn",b,"Cdn",b+1
      ampo += -t1,"Cdagdn",b+1,"Cdn",b
      ampo += V1,"Ntot",b,"Ntot",b+1
    end
    H = MPO(ampo,sites)
    sweeps = Sweeps(6)
    maxdim!(sweeps,50,100,200,400,800,800)
    cutoff!(sweeps,1E-10)
    state = ["Up",
             "Dn",
             "Dn",
             "Up",
             "Emp",
             "Up",
             "Up",
             "Emp",
             "Dn",
             "Dn"
            ]
    psi0 = randomMPS(sites,state,10)
    energy,psi = dmrg(H,psi0,sweeps; outputlevel = 0)
    @test (-8.02 < energy < -8.01)
  end

end

nothing


text:
##########################################################################################
### Logo
using JLD2, Plots
function createLogo()

    # Load stored data
    untarget = load("./src/data/calib1_untarget_small.jld2");
    target = load("./src/data/calib1_target_small.jld2");
    isolate = load("./src/data/calib1_isolate_small.jld2");

    # Create variables
    TT = untarget["TT"];
    IN_3 = 1 .- untarget["ST"] .- untarget["RT"];
    IN_4 = 1 .- target["ST"] .- target["RT"];
    IN_5 = 1 .- isolate["ST"] .- isolate["RT"];
    Yss = 1.01;
    Cdev_3 = ((untarget["YT"]-untarget["C1T"]).-Yss)./Yss.*100;
    Cdev_4 = ((target["YT"]-target["C1T"]).-Yss)./Yss.*100;
    Cdev_5 = ((isolate["YT"]-isolate["C1T"]).-Yss)./Yss.*100;

    # Figure 7
    args_first = (titlefontsize=16, guidefontsize=10, linecolor=["red" "blue" "purple"], legendfontsize=12, label=["untargeted" "targeted" "isolate"]);
    args_rest  = (titlefontsize=16, guidefontsize=10, linecolor=["red" "blue" "purple"], legend=false);
    p1 = plot(1:TT, [IN_3.*100 IN_4.*100 IN_5.*100], title="Infected, I_t", ylims=(-2,30); args_first...)
    plot!(size=(800,400))
    p6 = plot(1:TT, [Cdev_3 Cdev_4 Cdev_5], title="Agg. Consumption, C_t", ylabel="% dev. from steady state"; args_rest...)
    plot!(size=(800,400))
    plot(p1, p6, layout=(1,2))
    savefig("docs/src/assets/logo.png")

end

#createLogo()

function createHeader()
    # Load stored data
    untarget = load("./src/data/calib1_untarget_small.jld2");
    target = load("./src/data/calib1_target_small.jld2");
    isolate = load("./src/data/calib1_isolate_small.jld2");

    # Create variables
    TT = untarget["TT"];
    IN_3 = 1 .- untarget["ST"] .- untarget["RT"];
    IN_4 = 1 .- target["ST"] .- target["RT"];
    IN_5 = 1 .- isolate["ST"] .- isolate["RT"];
    death_3 = cumsum(untarget["Par"].gamma .* untarget["Par"].delta .* IN_3);
    death_4 = cumsum(target["Par"].gamma .* target["Par"].delta .* IN_4);
    death_5 = cumsum(isolate["Par"].gamma .* isolate["Par"].delta .* IN_5);
    Yss = 1.01;
    Cdev_3 = ((untarget["YT"]-untarget["C1T"]).-Yss)./Yss.*100;
    Cdev_4 = ((target["YT"]-target["C1T"]).-Yss)./Yss.*100;
    Cdev_5 = ((isolate["YT"]-isolate["C1T"]).-Yss)./Yss.*100;

    # Figure 7
    args_first = (titlefontsize=8, guidefontsize=6, linecolor=["red" "blue" "purple"], labelfontsize=1, label=["untargeted" "targeted" "isolate"]);
    args_rest  = (titlefontsize=8, guidefontsize=6, linecolor=["red" "blue" "purple"], legend=false);
    p1 = plot(1:TT, [IN_3.*100 IN_4.*100 IN_5.*100], title="Infected, I_t", ylabel=""; args_first...)
    plot!(size=(800,200))
    p2 = plot(1:TT, [untarget["ST"].*100 target["ST"].*100 isolate["ST"].*100], title="Susceptible, S_t", ylabel="%"; args_rest...)
    plot!(size=(800,200))
    p3 = plot(1:TT, [death_3.*100 death_4.*100 death_5.*100], title="Cumulative Deaths", ylabel="%"; args_rest...)
    plot!(size=(800,200))
    p6 = plot(1:TT, [Cdev_3 Cdev_4 Cdev_5], title="Agg. Consumption, C_t", ylabel="% dev. from steady state"; args_rest...)
    plot!(size=(800,200))
    p = plot(p1, p2, p3, p6, layout = (1, 4));
    savefig("docs/src/assets/header.png")

end

createHeader()


text:
using ControlSystems
using ControlSystems: AbstractStateSpace, ssdata, ninputs, noutputs

function basis_factor(p::Number)
    (A=p, B=sqrt(1-abs2(p)), C=sqrt(1-abs2(p)), D=-p)
end

# eq (3) in https://people.kth.se/~andersb/research/BF03.pdf
# IDENTIFICATION OF RATIONAL SPECTRAL DENSITIES USING ORTHONORMAL BASIS FUNCTIONS
function expand_basis(l, s) # large, small
    A = [l.A zeros(size(l.A,1), 1); s.B*l.C s.A]
    B = [l.B; s.B*l.D]
    C = [s.D*l.C s.C]
    D = s.D*l.D
    (; A, B, C, D)
end

function kautz(a::AbstractVector)
    h = 1/(pi*maximum(a))
    b = kautz(a, h)
    d2c(b)
end

function kautz(a::AbstractVector, h)
    if maximum(abs, a) > 1
        if all(>(0) ∘ real, a)
            a = -a
        end
        a = exp.(h .* a)
        @assert all(<(1) ∘ abs, a) "Ustable poles does not make sense"
    end
    factors = basis_factor.(a)
    Q = factors[1]
    for f in factors[2:end]
        Q = expand_basis(Q, f)
    end
    T = (Q...,)
    # A,B,C,D = real.(T)
    A,B,C,D = T
    C2 = I(size(A,1))
    ss(A,B,C2,0, h)
end

function laguerre_oo(a::Number, Nq)
    A = diagm(fill(-a, Nq))
    for i = 2:Nq
        A[diagind(A, i-1)] .+= 2a*(-1)^i
    end
    C = sqrt(2a)*(-1).^(0:Nq-1)'
    ss(A,I(Nq),C,0)
end

function laguerre(a, Nq)
    A = diagm(fill(-a, Nq))
    A[diagind(A, 1)] .= 1
    B = diagm(fill(1, Nq))
    for i = 1:Nq-1
        for (k,j) in enumerate(diagind(A, i))
            B[j] = (-1)^i * binomial(i+k-1, k-1)
        end
    end
    Ma = sqrt(2a)*diagm((2a).^(0:Nq-1))
    C = [1 zeros(1, Nq-1)]
    ss(A,Ma*B,C,0)
end

function adhocbasis(a)
    b = map(1:length(a)-1) do i
        ss(leadlinkat(a[i], 4, 1/2)*laglink(a[i+1], 2))
    end
    [b..., leadlinkat(a[end], 4, 1/2)]
end

function simplified_laguerre(a)
    map(enumerate(a)) do (i,ai)
        tf(1, [1.0, ai])^i
    end
end

# Utils ===============================================================================
const SomeBasis = Union{AbstractVector{<:LTISystem}, LTISystem}

basislength(v::AbstractVector) = length(v)
function basislength(v::LTISystem)
    n = ninputs(v)
    n == 1 ? noutputs(v) : n
end

function add_poles(basis::AbstractStateSpace, sys::LTISystem)
    add_poles(basis, ss(sys).A)
end

function add_poles(basis::AbstractStateSpace, Ad)
    A,B,C,D = ssdata(basis)
    T = eltype(A)
    nx,nu,ny = basis.nx,basis.nu,basis.ny
    nn = size(Ad, 1)
    @assert size(A) == size(B) "only supports laguerre_oo basis"
    B == I || @warn("B was not identity, the B matrix will be overwritten with I")

    Ae = [A zeros(size(A,1), nn); zeros(T, nn, nx) Ad]
    # Be = [B; zeros(T, nn, nu)]
    Ce = [C ones(T, ny, nn)]
    ss(Ae,I(nn+nx),Ce,0)
end

function ControlSystems.freqresp(P::LTISystem, basis::SomeBasis, ω::AbstractVector, p)
    @assert ninputs(P) == noutputs(P) == 1 "Only supports SISO"
    resp_P = freqresp(P, ω) |> vec
    basis_resp = basis_responses(basis, ω, inverse=false)
    @assert length(basis_resp) == length(p) "p does not have the same length as the basis"

    FPs = map(basis_resp) do Fi
        Fi .* resp_P
    end

    Gmat = reduce(hcat, FPs)
    Gmat*p
end

function sum_basis(basis::AbstractVector, p::AbstractVector)
    out = ss(p[1]*basis[1])
    for i in 2:length(p)
        out = out + balreal(p[i] * basis[i])[1]
    end
    out
end

function sum_basis(basis::AbstractStateSpace, p::AbstractVector)
    A,B,C,D = ssdata(basis)
    @assert all(iszero, D)
    if ninputs(basis) > 1
        isdiscrete(basis) ? ss(A,B*p,C,0,basis.Ts) : ss(A,B*p,C,0)
    else
        isdiscrete(basis) ? ss(A,B,p'C,0,basis.Ts) : ss(A,B,p'C,0)
    end
end

function basis_responses(basis::AbstractVector, ω; inverse=false)
    inverse && (basis = inv.(basis))
    Gs = freqresp.(basis, Ref(ω))
    vec.(Gs)
end

function basis_responses(basis::AbstractStateSpace, ω; inverse=false)
    inverse && (basis = inv(basis))
    Gs = freqresp(basis, ω)
    if ninputs(basis) > 1
        vec.(eachslice(Gs, dims=3))
    else
        vec.(eachslice(Gs, dims=2))
    end
end

function filter_bank(basis::AbstractStateSpace{<:Discrete}, signal::AbstractVector)
    no = ninputs(basis)
    t = range(0, step=basis.Ts, length=length(signal))
    if no > 1
        y,_,_ = lsim(basis, repeat(signal, 1, no), t)
    else
        y,_,_ = lsim(basis, signal, t)
    end
    y
end

function ωζ2complex(ω, ζ)
    p = ω*cis(pi-acos(ζ))
    p
end

Base.one(::TransferFunction{Continuous, ControlSystems.SisoRational{Float64}}) = tf(1)

reflect(x) = complex(-abs(real(x)), imag(x))
function minimum_phase(G::TransferFunction)
    z,p,k = zpkdata(G) .|> first
    z = reflect.(z)
    p = reflect.(p)
    Gmp = tf(zpk(z,p,k))
    if sign(dcgain(G)[]) != sign(dcgain(Gmp)[])
        Gmp = -Gmp
    end
    Gmp
end

function minimum_phase(G::StateSpace)
    H = tf(G)
    G = ss(minimum_phase(H))
    G = ControlSystems.balance_statespace(G)[1]
    try
        G = balreal(G)[1]
    catch
    end
    G
end


text:
mutable struct CoreVisualizer
    scope::WebIO.Scope
    tree::SceneNode
    command_channel::Observable{Vector{UInt8}}
    request_channel::Observable{String}
    controls_channel::Observable{Vector{Any}}
    controls::Dict{String, Tuple{Observable, AbstractControl}}

    function CoreVisualizer()
        scope = WebIO.Scope(imports=ASSET_KEYS)
        command_channel = Observable(scope, "meshcat-command", UInt8[])
        request_channel = Observable(scope, "meshcat-request", "")
        controls_channel = Observable(scope, "meshcat-controls", [])

        onimport(scope, @js function(mc)
            @var element = this.dom.children[0]
            this.viewer = @new mc.Viewer(element)
            $request_channel[] = String(Date.now())
            window.document.body.style.margin = "0"
            window.meshcat_viewer = this.viewer
        end)

        onjs(command_channel, @js function(val)
            this.viewer.handle_command_message(Dict(:data => val))
        end)
        scope = scope(dom"div.meshcat-viewer"(
            style=Dict(
                :width => "100vw",
                :height => "100vh",
                :position => "absolute",
                :left => 0,
                :right => 0,
                :margin => 0,
            )
        ))
        scope.dom.props[:style][:overflow] = "hidden"

        tree = SceneNode()
        controls = Dict{String, Observable}()
        vis = new(scope, tree, command_channel, request_channel, controls_channel, controls)
        on(request_channel) do x
            send_scene(vis)
        end

        on(controls_channel) do msg
            name::String, value = msg
            if haskey(vis.controls, name)
                @async vis.controls[name] = value
                # Base.invokelatest(setindex!, vis.controls[name], value)
            end
        end
        vis
    end
end

WebIO.render(core::CoreVisualizer) = core.scope

function WebIO.iframe(core::CoreVisualizer; height="100%", width="100%", minHeight="400px")
    ifr = WebIO.iframe(core.scope)
    onimport(ifr, @js function()
        this.dom.style.height = "100%"
        window.foo = this
        this.dom.children[0].children[0].style.flexGrow = "1"
    end)
    style = get!(Dict, ifr.dom.props, :style)
    style[:height] = height
    style[:minHeight] = minHeight
    style[:width] = width
    style[:display] = "flex"
    style[:flexDirection] = "column"
    ifr
end

function update_tree!(core::CoreVisualizer, cmd::SetObject, data)
    core.tree[cmd.path].object = data
end

function update_tree!(core::CoreVisualizer, cmd::SetTransform, data)
    core.tree[cmd.path].transform = data
end

function update_tree!(core::CoreVisualizer, cmd::Delete, data)
    if length(cmd.path) == 0
        core.tree = SceneNode()
    else
        delete!(core.tree, cmd.path)
    end
end

update_tree!(core::CoreVisualizer, cmd::SetControl, data) = nothing
update_tree!(core::CoreVisualizer, cmd::SetAnimation, data) = nothing
update_tree!(core::CoreVisualizer, cmd::SetProperty, data) = nothing

function send_scene(core::CoreVisualizer)
    foreach(core.tree) do node
        if node.object !== nothing
            core.command_channel[] = node.object
        end
        if node.transform !== nothing
            core.command_channel[] = node.transform
        end
    end
    for (name, (obs, control)) in core.controls
        send(core, SetControl(control))
    end
end

function send(c::CoreVisualizer, cmd::AbstractCommand)
    data = pack(lower(cmd))
    update_tree!(c, cmd, data)
    c.command_channel[] = data
    nothing
end

function Base.wait(c::CoreVisualizer)
    pool = c.scope.pool
    while true
        while isready(pool.new_connections)
            push!(pool.connections, take!(pool.new_connections))
        end
        if !isempty(pool.connections)
            break
        else
            sleep(0.25)
        end
    end
end

"""
    vis = Visualizer()

Construct a new MeshCat visualizer instance.

Useful methods:

    vis[:group1] # get a new visualizer representing a sub-tree of the scene
    setobject!(vis, geometry) # set the object shown by this visualizer's sub-tree of the scene
    settransform!(vis], tform) # set the transformation of this visualizer's sub-tree of the scene
"""
struct Visualizer <: AbstractVisualizer
    core::CoreVisualizer
    path::Path
end

Visualizer() = Visualizer(CoreVisualizer(), ["meshcat"])

"""
$(SIGNATURES)

Wait until at least one browser has connected to the
visualizer's server. This is useful in scripts to delay
execution until the browser window has opened.
"""
Base.wait(v::Visualizer) = wait(v.core)

IJuliaCell(vis::Visualizer; kw...) = iframe(vis.core; kw...)

Base.show(io::IO, v::Visualizer) = print(io, "MeshCat Visualizer with path $(v.path)")

"""
$(SIGNATURES)

Set the object at this visualizer's path. This replaces whatever
geometry was presently at that path. To draw multiple geometries,
place them at different paths by using the slicing notation:

    setobject!(vis[:group1][:box1], geometry1)
    setobject!(vis[:group1][:box2], geometry2)
"""
function setobject!(vis::Visualizer, obj::AbstractObject)
    send(vis.core, SetObject(obj, vis.path))
    vis
end

"""
$(SIGNATURES)

Set the transform of this visualizer's path. This can be done
before or after adding an object at that path. The overall transform
of an object is the composition of the transforms of all of its parents,
so setting the transform of `vis[:group1]` affects the poses of the objects
at `vis[:group1][:box1]` and `vis[:group1][:box2]`.
"""
function settransform!(vis::Visualizer, tform::Transformation)
    send(vis.core, SetTransform(tform, vis.path))
    vis
end

"""
$(SIGNATURES)

Delete the geometry at this visualizer's path and all of its descendants.
"""
function delete!(vis::Visualizer)
    send(vis.core, Delete(vis.path))
    vis
end

"""
$(SIGNATURES)

Set a single property for the object at the given path.

(this is named setprop! instead of setproperty! to avoid confusion
with the Base.setproperty! function introduced in Julia v0.7)
"""
function setprop!(vis::Visualizer, property::AbstractString, value)
    send(vis.core, SetProperty(vis.path, property, value))
    vis
end

function setcontrol!(vis::Visualizer, name::AbstractString, obs::Observable)
    control = Button(vis.core.controls_channel, name)
    vis.core.controls[name] = (obs, control)
    send(vis.core, SetControl(control))
    vis
end

function setcontrol!(vis::Visualizer, name::AbstractString, obs::Observable, value, min=zero(value), max=one(value))
    control = NumericControl(vis.core.controls_channel, name, value, min, max)
    vis.core.controls[name] = (obs, control)
    send(vis.core, SetControl(control))
    vis
end

function setanimation!(vis::Visualizer, anim::Animation; play::Bool=true, repetitions::Integer=1)
    cmd = SetAnimation(anim, play, repetitions)
    send(vis.core, cmd)
end

Base.getindex(vis::Visualizer, path...) =
    Visualizer(vis.core, joinpath(vis.path, path...))


text:
#Function to calculate triangle coefficients for angular momenta
#Based on Matlab code written by Amita B. Deb, Clarendon Lab. 2007.

function triangle_coeff(a,b,c)
  return tri = factorial(a+b-c)*factorial(a-b+c)*factorial(-a+b+c)/factorial(a+b+c+1)
end

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
"""
    principal_value(R::Rotation{3})

**Background:** All non `RotMatrix` rotation types can represent the same `RotMatrix` in two or more ways. Sometimes a
particular set of numbers is better conditioned (e.g. `MRP`) or obeys a particular convention (e.g. `AngleAxis` has
non-negative rotation). In order to preserve differentiability it is necessary to allow rotation representations to
travel slightly away from the nominal domain; this is critical for applications such as optimization or dynamics.

This function takes a rotation type (e.g. `QuatRotation`, `RotXY`) and outputs a new rotation of the same type that corresponds
to the same `RotMatrix`, but that obeys certain conventions or is better conditioned. The outputs of the function have
the following properties:

- all angles are between between `-pi` to `pi` (except for `AngleAxis` which is between `0` and `pi`).
- all `QuatRotation have non-negative real part
- the components of all `MRP` have a norm that is at most 1.
- the `RotationVec` rotation is at most `pi`

"""
principal_value(r::Rotation) = r
principal_value(q::Q) where Q <: QuatRotation = real(q.q) < zero(eltype(q)) ? Q(-q.q) : q
function principal_value(spq::MRP{T}) where {T}
    # A quat with positive real part: QuatRotation( qw,  qx,  qy,  qz)
    #
    # A mrp corresponding to the QuatRotation with a positive real part:
    # MRP( qx / (1 + qw),  qy / (1 + qw),  qz / (1 + qw)) ≡ MRP(spx, spy, spz)
    #
    # A mrp corresponding to the QuatRotation with a negative real part:
    # MRP(-qx / (1 - qw), -qy / (1 - qw), -qz / (1 - qw)) ≡ MRP(snx, sny, snz)
    #
    # Claim:         spx / snx = -1 / (spx^2 + spy^2 + spz^2)
    #     -(1 + qw) / (1 - qw) = -1 / (spx^2 + spy^2 + spz^2)
    #      (1 - qw) / (1 + qw) = (spx^2 + spy^2 + spz^2)
    #      (1 - qw) / (1 + qw) = (qx^2 + qy^2 + qz^2) / (1 + qw)^2
    #      (1 - qw) * (1 + qw) =  qx^2 + qy^2 + qz^2
    #                 1 - qw^2 =  qx^2 + qy^2 + qz^2 (Q.E.D.)
    alpha_2 = spq.x^2 + spq.y^2 + spq.z^2
    if one(T) < alpha_2
        scale = -one(T) / alpha_2
        return MRP(scale * spq.x, scale * spq.y, scale * spq.z)
    else
        return spq
    end
end

function principal_value(aa::AngleAxis{T}) where {T}
    theta = rem2pi(aa.theta, RoundNearest)
    if theta < zero(T)
        return AngleAxis(-theta, -aa.axis_x, -aa.axis_y, -aa.axis_z, false)
    else
        return AngleAxis( theta,  aa.axis_x,  aa.axis_y,  aa.axis_z, false)
    end
end

function principal_value(rv::RotationVec{T}) where {T}
    theta = rotation_angle(rv)
    if pi < theta
        re_s = rem2pi(theta, RoundNearest) / theta
        return RotationVec(re_s * rv.sx, re_s * rv.sy, re_s * rv.sz)
    else
        return rv
    end
end

for rot_type in [:RotX, :RotY, :RotZ, :Angle2d]
    @eval begin
        function principal_value(r::$rot_type{T}) where {T}
            return $(rot_type){T}(rem2pi(r.theta, RoundNearest))
        end
    end
end

for rot_type in [:RotXY, :RotYX, :RotZX, :RotXZ, :RotYZ, :RotZY]
    @eval begin
        function principal_value(r::$rot_type{T}) where {T}
            theta1 = rem2pi(r.theta1, RoundNearest)
            theta2 = rem2pi(r.theta2, RoundNearest)
            return $(rot_type){T}(theta1, theta2)
        end
    end
end

for rot_type in [:RotXYX, :RotYXY, :RotZXZ, :RotXZX, :RotYZY, :RotZYZ, :RotXYZ, :RotYXZ, :RotZXY, :RotXZY, :RotYZX, :RotZYX]
    @eval begin
        function principal_value(r::$rot_type{T}) where {T}
            theta1 = rem2pi(r.theta1, RoundNearest)
            theta2 = rem2pi(r.theta2, RoundNearest)
            theta3 = rem2pi(r.theta3, RoundNearest)
            return $(rot_type){T}(theta1, theta2, theta3)
        end
    end
end


text:
module SnowflakeWrapper

using DataFramesMeta, DataFrames, ODBC, CSV, IterTools, Arrow, Base.Iterators

INSERT_CHUNK_LIMIT = 16000
REGEX = r"({ \w+ })"

struct dbClass
    conn::ODBC.Connection
end

function initWrapper()
    env_keys = keys(ENV)
    if !("ACCOUNT" in env_keys) || !("DBUSER" in env_keys) || !("PASSWD" in env_keys)
        @warn "Missing value in env Either ACCOUNT , USER or PASSWD"
        return nothing
    end
    @debug "Opening connection from ENV settings"
    dbClass(ODBC.Connection("Driver=$(ENV["DRIVER_PATH"]);
    Server=$(ENV["ACCOUNT"]).snowflakecomputing.com;Warehouse=$(get(ENV, "WAREHOUSE", "USER_WH"));Role=$(get(ENV, "ROLE", "DATA_ANALYST"));Database=$(get(ENV, "DATABASE", "MODEL_DB"));Schema=$(get(ENV, "SCHEMA", "PUBLIC"))", ENV["DBUSER"], ENV["PASSWD"]))
end

function DbConstructor(databaseUrl::String)
    @debug "Opening connection from databaseUrl : $databaseUrl"
    dbClass(ODBC.Connection(databaseUrl, ENV["DBUSER"], ENV["PASSWD"]))
end

function execQuery(db, query::String)
    @debug "Exec Query $query"
    DBInterface.execute(db.conn, query)
end

function execQueryFromFile(db, filePath::String, replaceFromEnv=true)
    query = open(f -> read(f, String), filePath)
    if replaceFromEnv
        for m in eachmatch(REGEX, query)
            val = get(ENV, chop(m.match, head=2, tail=2), missing)
            query = replace(query, m.match => val)
        end
    end
    @debug "Exec QueryFromFile $query"
    DBInterface.execute(db.conn, query)
end

toChunk(df::DataFrame, n) = [df[i:min(i + n - 1, end),:] for i in 1:n:nrow(df)]

function write(db, df, tableName)
    _prepare_field(x::Any) = x
    _prepare_field(x::Missing) = ""
    _prepare_field(x::Nothing) = ""
    _prepare_field(x::AbstractString) = "'" * replace(x, "'" => '"') * "'"
    _prepare_field(x::AbstractFloat) = isnan(x) ? 0.0 : round(x, digits=2)

    params = join(names(df), ',')
    chunks = toChunk(df, INSERT_CHUNK_LIMIT)
    for chunk in chunks
        rowsStrings = replace(join(map(x -> '(' * join(x, ',') * ')', eachrow(_prepare_field.(chunk))), ','), ",," => ",'',")
        query = "INSERT INTO $tableName($params) VALUES $rowsStrings"
        if nrow(chunk) > 0
            DBInterface.execute(db.conn, query)
        end
    end
    return nrow(df)
end

function close!(db::dbClass)
    DBInterface.close!(db.conn)
end

end


text:
function parse(obj::EOPHistoryLoader_Parser, arg0::InputStream, arg1::JString)
    return jcall(obj, "parse", Collection, (InputStream, JString), arg0, arg1)
end


text:
int main () {
printInt(fact(7)) ;
printInt(factr(7)) ;
return 0 ;
}
// iterative factorial
int fact (int n) {
int i,r ;
i = 1 ;
r = 1 ;
while (i <= n) {
r = r * i ;
i++ ;
}
return r ;
}
// recursive factorial
int factr (int n) {
if (n < 2)
return 1 ;
else
return n * factr(n-1) ;
}

max_stars_repo_name: fredefox/javalette

max_issues_repo_name: fredefox/javalette

max_forks_repo_name: fredefox/javalette


text:
module Fitting

using Compat
using Distributions
using Optim
using GLM
using Formatting

export fit_weibul, simple_fit, r2, min_sqerror_with, ff

function ff(f::Real, n::Int=1)
    abs(f) < 1/10^n || abs(f) > 10^n ? fmt(FormatSpec("$(n)e"), f) : fmt(FormatSpec("$(n)f"), f)
end

function fit_weibul(data)

    y = log(sort(data))
    rank = 1:length(data)  # ranks = {1, 2, 3, ... 10}
    median_rank = (rank - 0.3)/(length(rank) + 0.4)
    x = log(-log(1 - median_rank))

    intercept, slope = linreg(x,y)

    shape = 1/slope
    x_intercept = - intercept / shape
    scale = exp(-x_intercept/slope)

#     @show(shape, scale)
    return Weibull(shape, scale)
end

function simple_fit(x, y)

    fit = GLM.lm([ones(x) x], y)
    b,a = fit.pp.beta0;

    fit_r2 = r2(y, fit);

    return @compat Dict( :fit => fit, :r2 => fit_r2, :a => a, :b => b, :f => (x -> a*x + b),
        :fstr => "\$$( ff(a,1)) \\,x + $(ff(b,1))\$"  )
end

function r2(ydata, fit)
    1-sumabs2(GLM.residuals(fit))/sumabs2(ydata - mean(ydata)) # or 1-var(residuals(mod))/var(y)
end

function min_sqerror_with(func::Function, x, y)
    function sqerror(C)
        err = 0.0
        for i in 1:length(x)
            ŷ = func(x[i], C)
            err += (y[i] - ŷ)^2
        end
        return err
    end

    return sqerror
end

function optimize(func::Function, x::Array, y::Array, params)
    return optimize(min_sqerror_with(func, x, y), Float64[params...])
end

end


text:
# p15.jl - solve eigenvalue BVP u_xx = lambda*u, u(-1)=u(1)=0

N = 36; (D,x) = cheb(N); D2 = D^2; D2 = D2[2:N,2:N];
(lam,V) = eig(D2);
ii = sortperm(-lam);          # sort eigenvalues and -vectors
lam = lam[ii]; V = V[:,ii]; clf();
for j = 5:5:30                  # plot 6 eigenvectors
    u = [0;V[:,j];0]; subplot(7,1,j/5)
    plot(x,u,".",markersize=6); grid(true);
    xx = -1:.01:1; uu = polyval(polyfit(x,u),xx);
    plot(xx,uu); axis("off");
    text(-.4,.1,"eig $j = $(lam[j]*4/pi^2) π^2/4");
    text(.7,.1,"$(signif(4*N/(pi*j),2))  ppw");
end


text:
"""
#### Read the leaf-cells of the hydro-data:
- select variables
- limit to a maximum level
- limit to a spatial range
- set a minimum density or sound speed
- check for negative values in density and thermal pressure
- print the name of each data-file before reading it
- toggle verbose mode

```julia
gethydro(   dataobject::InfoType;
            lmax::Real=dataobject.levelmax,
            vars::Array{Symbol,1}=[:all],
            xrange::Array{<:Any,1}=[missing, missing],
            yrange::Array{<:Any,1}=[missing, missing],
            zrange::Array{<:Any,1}=[missing, missing],
            center::Array{<:Any,1}=[0., 0., 0.],
            range_unit::Symbol=:standard,
            smallr::Real=0.,
            smallc::Real=0.,
            check_negvalues::Bool=false,
            print_filenames::Bool=false,
            verbose::Bool=verbose_mode )
```
#### Returns an object of type HydroDataType, containing the hydro-data table, the selected options and the simulation ScaleType and summary of the InfoType
```julia
return HydroDataType()

# get an overview of the returned fields:
# e.g.:
julia> info = getinfo(100)
julia> gas  = gethydro(info)
julia> viewfields(gas)
#or:
julia> fieldnames(gas)
```

#### Arguments
##### Required:
- **`dataobject`:** needs to be of type: "InfoType", created by the function *getinfo*
##### Predefined/Optional Keywords:
- **`lmax`:** the maximum level to be read from the data
- **`var(s)`:** the selected hydro variables in arbitrary order: :all (default), :cpu, :rho, :vx, :vy, :vz, :p, :var6, :var7...
- **`xrange`:** the range between [xmin, xmax] in units given by argument `range_unit` and relative to the given `center`; zero length for xmin=xmax=0. is converted to maximum possible length
- **`yrange`:** the range between [ymin, ymax] in units given by argument `range_unit` and relative to the given `center`; zero length for ymin=ymax=0. is converted to maximum possible length
- **`zrange`:** the range between [zmin, zmax] in units given by argument `range_unit` and relative to the given `center`; zero length for zmin=zmax=0. is converted to maximum possible length
- **`range_unit`:** the units of the given ranges: :standard (code units), :Mpc, :kpc, :pc, :mpc, :ly, :au , :km, :cm (of typye Symbol) ..etc. ; see for defined length-scales viewfields(info.scale)
- **`center`:** in units given by argument `range_unit`; by default [0., 0., 0.]; the box-center can be selected by e.g. [:bc], [:boxcenter], [value, :bc, :bc], etc..
- **`smallr`:** set lower limit for density; zero means inactive
- **`smallc`:** set lower limit for thermal pressure; zero means inactive
- **`check_negvalues`:** check loaded data of "rho" and "p" on negative values; false by default
- **`print_filenames`:** print on screen the current processed hydro file of each CPU
- **`verbose`:** print timestamp, selected vars and ranges on screen; default: set by the variable `verbose_mode`

### Defined Methods - function defined for different arguments
- gethydro( dataobject::InfoType; ...) # no given variables -> all variables loaded
- gethydro( dataobject::InfoType, var::Symbol; ...) # one given variable -> no array needed
- gethydro( dataobject::InfoType, vars::Array{Symbol,1}; ...)  # several given variables -> array needed

#### Examples
```julia
# read simulation information
julia> info = getinfo(420)

# Example 1:
# read hydro data of all variables, full-box, all levels
julia> gas = gethydro(info)

# Example 2:
# read hydro data of all variables up to level 8
# data range 20x20x4 kpc; ranges are given in kpc relative to the box (here: 48 kpc) center at 24 kpc
julia> gas = gethydro(    info,
                          lmax=8,
                          xrange=[-10.,10.],
                          yrange=[-10.,10.],
                          zrange=[-2.,2.],
                          center=[24., 24., 24.],
                          range_unit=:kpc )

# Example 3:
# give the center of the box by simply passing: center = [:bc] or center = [:boxcenter]
# this is equivalent to center=[24.,24.,24.] in Example 2
# the following combination is also possible: e.g. center=[:bc, 12., 34.], etc.
julia> gas = gethydro(    info,
                          lmax=8,
                          xrange=[-10.,10.],
                          yrange=[-10.,10.],
                          zrange=[-2.,2.],
                          center=[33., bc:, 10.],
                          range_unit=:kpc )

# Example 4:
# read hydro data of the variables density and the thermal pressure, full-box, all levels
julia> gas = gethydro( info, [:rho, :p] ) # use array for the variables

# Example 5:
# read hydro data of the single variable density, full-box, all levels
julia> gas = gethydro( info, :rho ) # no array for a single variable needed
...
```

"""
function gethydro( dataobject::InfoType, var::Symbol;
                    lmax::Real=dataobject.levelmax,
                    xrange::Array{<:Any,1}=[missing, missing],
                    yrange::Array{<:Any,1}=[missing, missing],
                    zrange::Array{<:Any,1}=[missing, missing],
                    center::Array{<:Any,1}=[0., 0., 0.],
                    range_unit::Symbol=:standard,
                    smallr::Real=0.,
                    smallc::Real=0.,
                    check_negvalues::Bool=false,
                    print_filenames::Bool=false,
                    verbose::Bool=verbose_mode )
                    #, progressbar::Bool=show_progressbar)

    return gethydro(dataobject, vars=[var],
                    lmax=lmax,
                    xrange=xrange, yrange=yrange, zrange=zrange, center=center,
                    range_unit=range_unit,
                    smallr=smallr,
                    smallc=smallc,
                    check_negvalues=check_negvalues,
                    print_filenames=print_filenames,
                    verbose=verbose)
end

function gethydro( dataobject::InfoType, vars::Array{Symbol,1};
                    lmax::Real=dataobject.levelmax,
                    xrange::Array{<:Any,1}=[missing, missing],
                    yrange::Array{<:Any,1}=[missing, missing],
                    zrange::Array{<:Any,1}=[missing, missing],
                    center::Array{<:Any,1}=[0., 0., 0.],
                    range_unit::Symbol=:standard,
                    smallr::Real=0.,
                    smallc::Real=0.,
                    check_negvalues::Bool=false,
                    print_filenames::Bool=false,
                    verbose::Bool=verbose_mode )
                    #, progressbar::Bool=show_progressbar)

    return gethydro(dataobject,
                    vars=vars,
                    lmax=lmax,
                    xrange=xrange, yrange=yrange, zrange=zrange, center=center,
                    range_unit=range_unit,
                    smallr=smallr,
                    smallc=smallc,
                    check_negvalues=check_negvalues,
                    print_filenames=print_filenames,
                    verbose=verbose)
end

function gethydro( dataobject::InfoType;
                    lmax::Real=dataobject.levelmax,
                    vars::Array{Symbol,1}=[:all],
                    xrange::Array{<:Any,1}=[missing, missing],
                    yrange::Array{<:Any,1}=[missing, missing],
                    zrange::Array{<:Any,1}=[missing, missing],
                    center::Array{<:Any,1}=[0., 0., 0.],
                    range_unit::Symbol=:standard,
                    smallr::Real=0.,
                    smallc::Real=0.,
                    check_negvalues::Bool=false,
                    print_filenames::Bool=false,
                    verbose::Bool=verbose_mode )
                    #, progressbar::Bool=show_progressbar)

    printtime("Get hydro data: ", verbose)
    checkfortype(dataobject, :hydro)
    checklevelmax(dataobject, lmax)
    isamr = checkuniformgrid(dataobject, lmax)

    # create variabe-list and vector-mask (nvarh_corr) for gethydrodata-function
    # print selected variables on screen
    nvarh_list, nvarh_i_list, nvarh_corr, read_cpu, used_descriptors = prepvariablelist(dataobject, :hydro, vars, lmax, verbose)

    # convert given ranges and print overview on screen
    ranges = prepranges(dataobject, range_unit, verbose, xrange, yrange, zrange, center)

    # read hydro-data of the selected variables
    if read_cpu
        vars_1D, pos_1D, cpus_1D = gethydrodata( dataobject, length(nvarh_list),
                                         nvarh_corr, lmax, ranges,
                                         print_filenames, read_cpu, isamr )
    else
        vars_1D, pos_1D          = gethydrodata( dataobject, length(nvarh_list),
                                         nvarh_corr, lmax, ranges,
                                         print_filenames, read_cpu, isamr )
    end

    # set minimum density in cells and check vor negative values
    vars_1D = manageminvalues(vars_1D, check_negvalues, smallr, smallc, nvarh_list, nvarh_corr)

    # prepare column names for the data table
    names_constr = preptablenames(dataobject.nvarh, nvarh_list, used_descriptors, read_cpu, isamr)

    # create data table
    # decouple pos_1D/vars_1D from ElasticArray with ElasticArray.data
    if read_cpu # load also cpu number related to cell
        if isamr
            data = @inbounds table(pos_1D[4,:].data, cpus_1D[:], pos_1D[1,:].data, pos_1D[2,:].data, pos_1D[3,:].data,
                        [vars_1D[nvarh_corr[i],: ] for i in nvarh_i_list]...,
                        names=collect(names_constr), pkey=[:level,:cx, :cy, :cz], presorted = false ) #[names_constr...]
        else # if uniform grid
            data = @inbounds table(cpus_1D[:], pos_1D[1,:].data, pos_1D[2,:].data, pos_1D[3,:].data,
                        [vars_1D[nvarh_corr[i],: ].data for i in nvarh_i_list]...,
                        names=collect(names_constr), pkey=[:cx, :cy, :cz], presorted = false ) #[names_constr...]
        end
    else
        if isamr
            data = @inbounds table(pos_1D[4,:].data, pos_1D[1,:].data, pos_1D[2,:].data, pos_1D[3,:].data,
                        [vars_1D[ nvarh_corr[i],: ].data for i in nvarh_i_list]...,
                        names=collect(names_constr), pkey=[:level,:cx, :cy, :cz], presorted = false ) #[names_constr...]
        else # if uniform grid
            data = @inbounds table(pos_1D[1,:].data, pos_1D[2,:].data, pos_1D[3,:].data,
                        [vars_1D[ nvarh_corr[i],: ].data for i in nvarh_i_list]...,
                        names=collect(names_constr), pkey=[:cx, :cy, :cz], presorted = false ) #[names_constr...]
        end
    end

    printtablememory(data, verbose)

    # Return data
    hydrodata = HydroDataType()
    hydrodata.data = data
    hydrodata.info = dataobject
    hydrodata.lmin = dataobject.levelmin
    hydrodata.lmax = lmax
    hydrodata.boxlen = dataobject.boxlen
    hydrodata.ranges = ranges
    if read_cpu
        hydrodata.selected_hydrovars = [-1, nvarh_list...]
    else
        hydrodata.selected_hydrovars = nvarh_list
    end
    hydrodata.used_descriptors = used_descriptors
    hydrodata.smallr = smallr
    hydrodata.smallc = smallc
    hydrodata.scale = dataobject.scale
    return hydrodata
end

function manageminvalues(vars_1D::ElasticArray{Float64,2,1}, check_negvalues::Bool, smallr::Real, smallc::Real, nvarh_list::Array{Int,1}, nvarh_corr::Array{Int,1})

    # set minimum density in cells
    if smallr != 0. && in(1, nvarh_list)

        vars_1D[1,:] =clamp.(vars_1D[nvarh_corr[1],:], smallr, maximum(vars_1D[nvarh_corr[1],:]) + 1 )

    else
        # check for negative values in density
        if check_negvalues == true
            if in(1, nvarh_list)
                count_nv = count(x->x<0., vars_1D[nvarh_corr[1],:])
                if count_nv > 0
                    println()
                    println("[Mera]: Found $count_nv negative value(s) in density data.")
                end
            end
        end
    end

    # set minimum thermal pressure in cells
    if smallc != 0.  && in(5, nvarh_list)
        vars_1D[5,:] =clamp.(vars_1D[nvarh_corr[5],:], smallr, maximum(vars_1D[nvarh_corr[5],:]) + 1 )

    else
        # check for negative values in thermal pressure
        if check_negvalues == true
            if in(5, nvarh_list)
                count_nv = count(x->x<0., vars_1D[nvarh_corr[5],:])
                if count_nv > 0
                    println()
                    println("[Mera]: Found $count_nv negative value(s) in thermal pressure data.")
                end
            end
        end
    end

    return vars_1D
end

function preptablenames(nvarh::Int, nvarh_list::Array{Int, 1}, used_descriptors::Dict{Any,Any}, read_cpu::Bool, isamr::Bool)

    if read_cpu
        if isamr
            names_constr = [Symbol("level") ,Symbol("cpu"), Symbol("cx"), Symbol("cy"), Symbol("cz")]
        else    #if uniform grid
            names_constr = [Symbol("cpu"), Symbol("cx"), Symbol("cy"), Symbol("cz")]
        end
                    #, Symbol("x"), Symbol("y"), Symbol("z")
    else
        if isamr
            names_constr = [Symbol("level") , Symbol("cx"), Symbol("cy"), Symbol("cz")]
        else    #if uniform grid
            names_constr = [Symbol("cx"), Symbol("cy"), Symbol("cz")]
        end
    end

    for i=1:nvarh
        if in(i, nvarh_list)
            if length(used_descriptors) == 0 || !haskey(used_descriptors, i)
                if i == 1
                    append!(names_constr, [Symbol("rho")] )
                elseif i == 2
                    append!(names_constr, [Symbol("vx")] )
                elseif i == 3
                    append!(names_constr, [Symbol("vy")] )
                elseif i == 4
                    append!(names_constr, [Symbol("vz")] )
                elseif i == 5
                    append!(names_constr, [Symbol("p")] )
                elseif i > 5
                    append!(names_constr, [Symbol("var$i")] )
                end
            else append!(names_constr, [used_descriptors[i]] )

            end
        end
    end

    return names_constr
end