text:
module PyCall

export pyinitialize, pyfinalize, pycall, pyimport, pybuiltin, PyObject,
       pysym, PyPtr, pyincref, pydecref, pyversion, PyArray, PyArray_Info,
       pyerr_check, pyerr_clear, pytype_query, PyAny, @pyimport, PyDict,
       pyisinstance, pywrap, pytypeof, pyeval, pyhassym, PyVector, pystring,
       pyraise, pytype_mapping, pygui, pygui_start, pygui_stop,
       pygui_stop_all, @pylab, set!, PyTextIO

import Base: size, ndims, similar, copy, getindex, setindex!, stride,
       convert, pointer, summary, convert, show, haskey, keys, values,
       eltype, get, delete!, empty!, length, isempty, start, done,
       next, filter!, hash, splice!, pop!, isequal, help, push!,
       unshift!, shift!, append!, insert!, prepend!

# Python C API is not interrupt-save.  In principle, we should
# use sigatomic for every ccall to the Python library, but this
# should really be fixed in Julia (#2622).  However, we will
# use the sigatomic_begin/end functions to protect pycall and
# similar long-running (or potentially long-running) code.
import Base: sigatomic_begin, sigatomic_end

#########################################################################

# Mirror of C PyObject struct (for non-debugging Python builds).
# We won't actually access these fields directly; we'll use the Python
# C API for everything.  However, we need to define a unique Ptr type
# for PyObject*, and we might as well define the actual struct layout
# while we're at it.
immutable PyObject_struct
    ob_refcnt::Int
    ob_type::Ptr{Void}
end

typealias PyPtr Ptr{PyObject_struct} # type for PythonObject* in ccall

#########################################################################

# Global configuration variables.  Note that, since Julia does not allow us
# to attach type annotations to globals, we need to annotate these explicitly
# as initialized::Bool and libpython::Ptr{Void} when we use them.
initialized = false # whether Python is initialized
finalized = false # whether Python has been finalized
libpython = C_NULL # Python shared library (from dlopen)

pysym(func::Symbol) = dlsym(libpython::Ptr{Void}, func)
pysym_e(func::Symbol) = dlsym_e(libpython::Ptr{Void}, func)
pyhassym(func::Symbol) = pysym_e(func) != C_NULL

# call pysym_e on the arguments and return the first non-NULL result
function pysym_e(funcs...)
    for func in funcs
        p = pysym_e(func)
        if p != C_NULL
            return p
        end
    end
    return C_NULL
end

# Macro version of pysym to cache dlsym lookup (thanks to vtjnash)
macro pysym(func)
    z = gensym(string(func))
    @eval global $z = C_NULL
    quote begin
        global $z
        if $z::Ptr{Void} == C_NULL
            $z::Ptr{Void} = dlsym(libpython::Ptr{Void}, $(esc(func)))
        end
        $z::Ptr{Void}
    end end
end

# Macro version of pyinitialize() to inline initialized? check
macro pyinitialize()
    :(initialized::Bool ? nothing : pyinitialize())
end

#########################################################################
# Wrapper around Python's C PyObject* type, with hooks to Python reference
# counting and conversion routines to/from C and Julia types.

type PyObject
    o::PyPtr # the actual PyObject*
    function PyObject(o::PyPtr)
        po = new(o)
        finalizer(po, pydecref)
        return po
    end
    PyObject() = PyObject(convert(PyPtr, C_NULL))
end

function pydecref(o::PyObject)
    if initialized::Bool # don't decref after pyfinalize!
        ccall((@pysym :Py_DecRef), Void, (PyPtr,), o.o)
    end
    o.o = convert(PyPtr, C_NULL)
    o
end

function pyincref(o::PyObject)
    ccall((@pysym :Py_IncRef), Void, (PyPtr,), o)
    o
end

# doing an incref *before* creating a PyObject may safer in the
# case of borrowed references, to ensure that no exception or interrupt
# induces a double decref.
function pyincref(o::PyPtr)
    ccall((@pysym :Py_IncRef), Void, (PyPtr,), o)
    PyObject(o)
end

pyisinstance(o::PyObject, t::PyObject) =
  t.o != C_NULL && ccall((@pysym :PyObject_IsInstance), Cint, (PyPtr,PyPtr), o, t.o) == 1

pyisinstance(o::PyObject, t::Ptr{Void}) =
  t != C_NULL && ccall((@pysym :PyObject_IsInstance), Cint, (PyPtr,PyPtr), o, t) == 1

pyquery(q::Ptr{Void}, o::PyObject) =
  ccall(q, Cint, (PyPtr,), o) == 1

pytypeof(o::PyObject) = o.o == C_NULL ? throw(ArgumentError("NULL PyObjects have no Python type")) : pycall(TypeType, PyObject, o)

# conversion to pass PyObject as ccall arguments:
convert(::Type{PyPtr}, po::PyObject) = po.o

# use constructor for generic conversions to PyObject
convert(::Type{PyObject}, o) = PyObject(o)
PyObject(o::PyObject) = o

#########################################################################

inspect = PyObject() # inspect module, needed for module introspection
builtin = PyObject() # __builtin__ module, needed for pybuiltin

# Python has zillions of types that a function be, in addition to the FunctionType
# in the C API.  We have to obtain these at runtime and cache them in globals
BuiltinFunctionType = PyObject()
TypeType = PyObject() # "type" function
MethodType = PyObject()
MethodWrapperType = PyObject()
# also WrapperDescriptorType = type(list.__add__) and
#      MethodDescriptorType = type(list.append) ... is it worth detecting these?

# special function type used in NumPy and SciPy (if available)
ufuncType = PyObject()

# cache Python None
pynothing = convert(PyPtr,C_NULL) # PyPtr, not PyObject, to prevent it from
                                  # being finalized prematurely on exit

# Python 2/3 compatibility: cache dlsym for renamed functions
pystring_fromstring = C_NULL
pystring_asstring = C_NULL
pystring_size = C_NULL
pystring_type = C_NULL
pyint_type = C_NULL
pyint_from_size_t = C_NULL
pyint_from_ssize_t = C_NULL
pyint_as_ssize_t = C_NULL

pyxrange = PyObject() # xrange type (or range in Python 3)

# cache ctypes.c_void_p type and function if available
c_void_p_Type = PyObject()
py_void_p = p::Ptr -> PyObject(uint(p))

# PyCObject_Check and PyCapsule_CheckExact are actually macros
# that check against PyCObject_Type and PyCapsule_Type global variables,
# which we cache if they are available:
PyCObject_Type = C_NULL
PyCapsule_Type = C_NULL

# Py_SetProgramName needs its argument to persist as long as Python does
pyprogramname = bytestring("")

# PyUnicode_* may actually be a #define for another symbol,
# so we cache the correct dlsym
PyUnicode_AsUTF8String = C_NULL
PyUnicode_DecodeUTF8 = C_NULL

# mpmath mpf/mpc types
mpmath = PyObject()
mpf = PyObject()
mpc = PyObject()

# whether to use unicode for strings by default, ala Python 3
pyunicode_literals = false

# traceback.format_tb function, for show(PyError)
format_traceback = PyObject()

# cache whether Python hash values are Clong (Python < 3.2) or Int (>= 3.2)
pyhashlong = false

# Return the Python version as a Julia VersionNumber
pyversion = v"0"

# low-level initialization, given a pointer to dlopen result on libpython,
# or C_NULL if python symbols are in the global namespace:
# initialize the Python interpreter (no-op on subsequent calls)
function pyinitialize(libpy::Ptr{Void})
    global initialized
    global finalized
    global libpython
    global pyprogramname
    global inspect
    global BuiltinFunctionType
    global TypeType
    global MethodType
    global MethodWrapperType
    global ufuncType
    global PyUnicode_AsUTF8String
    global PyUnicode_DecodeUTF8
    global pyunicode_literals
    global pynothing
    global pyxrange
    global pyhashlong
    global pystring_fromstring
    global pystring_asstring
    global pystring_size
    global pystring_type
    global pyint_type
    global pyint_from_size_t
    global pyint_from_ssize_t
    global pyint_as_ssize_t
    global pyversion
    global PyCObject_Type
    global PyCapsule_Type
    global c_void_p_Type
    global py_void_p
    global format_traceback
    if !initialized::Bool
        if finalized::Bool
            # From the Py_Finalize documentation:
            #    "Some extensions may not work properly if their
            #     initialization routine is called more than once; this
            #     can happen if an application calls Py_Initialize()
            #     and Py_Finalize() more than once."
            # For example, numpy and scipy seem to crash if this is done.
            error("Calling pyinitialize after pyfinalize is not supported")
        end
        libpython::Ptr{Void} = libpy == C_NULL ? ccall(:jl_load_dynamic_library, Ptr{Void}, (Ptr{Uint8},Cuint), C_NULL, 0) : libpy
        already_inited = 0 != ccall((@pysym :Py_IsInitialized), Cint, ())
        if !already_inited
            if !isempty(pyprogramname)
                if pyversion.major < 3
                    ccall((@pysym :Py_SetProgramName), Void, (Ptr{Uint8},),
                          pyprogramname::ASCIIString)
                else
                    ccall((@pysym :Py_SetProgramName), Void, (Ptr{Cwchar_t},),
                          pyprogramname::Array{Cwchar_t})
                end
            end
            ccall((@pysym :Py_InitializeEx), Void, (Cint,), 0)
        end
        initialized::Bool = true
        inspect::PyObject = pyimport("inspect")
        types = pyimport("types")
        BuiltinFunctionType::PyObject = types["BuiltinFunctionType"]
        TypeType::PyObject = pybuiltin("type")
        MethodType::PyObject = types["MethodType"]
        MethodWrapperType::PyObject = pytypeof(PyObject(PyObject[])["__add__"])
        try
            ufuncType = pyimport("numpy")["ufunc"]
        catch
            ufuncType = PyObject() # NumPy not available
        end
        PyUnicode_AsUTF8String::Ptr{Void} =
          pysym_e(:PyUnicode_AsUTF8String,
                  :PyUnicodeUCS4_AsUTF8String,
                  :PyUnicodeUCS2_AsUTF8String)
        PyUnicode_DecodeUTF8::Ptr{Void} =
          pysym_e(:PyUnicode_DecodeUTF8,
                  :PyUnicodeUCS4_DecodeUTF8,
                  :PyUnicodeUCS2_DecodeUTF8)
        pynothing::PyPtr = convert(PyPtr, pysym(:_Py_NoneStruct))
        pyxrange::PyObject = pyincref(convert(PyPtr, pysym(:PyRange_Type)))
        if pyhassym(:PyString_FromString)
            pystring_fromstring::Ptr{Void} = pysym(:PyString_FromString)
            pystring_asstring::Ptr{Void} = pysym(:PyString_AsString)
            pystring_size::Ptr{Void} = pysym(:PyString_Size)
            pystring_type::Ptr{Void} = pysym(:PyString_Type)
        else
            pystring_fromstring::Ptr{Void} = pysym(:PyBytes_FromString)
            pystring_asstring::Ptr{Void} = pysym(:PyBytes_AsString)
            pystring_size::Ptr{Void} = pysym(:PyBytes_Size)
            pystring_type::Ptr{Void} = pysym(:PyBytes_Type)
        end
        if pyhassym(:PyInt_Type)
            pyint_type::Ptr{Void} = pysym(:PyInt_Type)
            pyint_from_size_t::Ptr{Void} = pysym(:PyInt_FromSize_t)
            pyint_from_ssize_t::Ptr{Void} = pysym(:PyInt_FromSsize_t)
            pyint_as_ssize_t::Ptr{Void} = pysym(:PyInt_AsSsize_t)
        else
            pyint_type::Ptr{Void} = pysym(:PyLong_Type)
            pyint_from_size_t::Ptr{Void} = pysym(:PyLong_FromSize_t)
            pyint_from_ssize_t::Ptr{Void} = pysym(:PyLong_FromSsize_t)
            pyint_as_ssize_t::Ptr{Void} = pysym(:PyLong_AsSsize_t)
        end
        PyCObject_Type::Ptr{Void} = pysym_e(:PyCObject_Type)
        PyCapsule_Type::Ptr{Void} = pysym_e(:PyCapsule_Type)
        try
            c_void_p_Type::PyObject = pyimport("ctypes")["c_void_p"]
            py_void_p::Function = p::Ptr -> pycall(c_void_p_Type::PyObject, PyObject, uint(p))
        catch # fallback to CObject
            pycobject_new = pysym(:PyCObject_FromVoidPtr)
            py_void_p::Function = p::Ptr -> PyObject(ccall(pycobject_new, PyPtr, (Ptr{Void}, Ptr{Void}), p, C_NULL))
        end
        pyversion::VersionNumber =
          VersionNumber(convert((Int,Int,Int,String,Int),
                                pyimport("sys")["version_info"])[1:3]...)
        pyhashlong::Bool = pyversion::VersionNumber < v"3.2"
        pyunicode_literals::Bool = pyversion::VersionNumber >= v"3.0"
        format_traceback::PyObject = pyimport("traceback")["format_tb"]
        pyexc_initialize()
        if !already_inited
            # some modules (e.g. IPython) expect sys.argv to be set
            if pyversion.major < 3
                argv_s = bytestring("")
                argv = convert(Ptr{Uint8}, argv_s)
                ccall(pysym(:PySys_SetArgvEx), Void, (Cint,Ptr{Ptr{Uint8}},Cint), 1, &argv, 0)
            else
                argv_s = Cwchar_t[0]
                argv   = convert(Ptr{Cwchar_t}, argv_s)
                ccall(pysym(:PySys_SetArgvEx), Void, (Cint, Ptr{Ptr{Cwchar_t}}, Cint), 1, &argv, 0)
            end
        end
    end
    return
end

pyconfigvar(python::String, var::String) = chomp(readall(`$python -c "import distutils.sysconfig; print(distutils.sysconfig.get_config_var('$var'))"`))
pysys(python::String, var::String) = chomp(readall(`$python -c "import sys; print(sys.$var)"`))
pyconfigvar(python, var, default) = let v = pyconfigvar(python, var)
    v == "None" ? default : v
end

function dlopen_libpython(python::String)
    lib = pyconfigvar(python, "LDLIBRARY")
    if lib == "None"
        lib = @windows ? "python" * pyconfigvar(python,"VERSION","") * ".dll" : (@osx ? "libpython.dylib" : "libpython.so")
    end
    libpaths = [pyconfigvar(python, "LIBDIR"),
                (@windows ? dirname(pysys(python, "executable")) : joinpath(dirname(dirname(pysys(python, "executable"))), "lib"))]
    @osx_only push!(libpaths, pyconfigvar(python, "PYTHONFRAMEWORKPREFIX"))

    # `prefix` and `exec_prefix` are the path prefixes where python should look for python only and compiled libraries, respectively.
    # These are also changed when run in a virtualenv.
    exec_prefix = pyconfigvar(python, "exec_prefix")
    # Since we only use `libpaths` to find the python dynamic library, we should only add `exec_prefix` to it.
    push!(libpaths, exec_prefix)
    if !haskey(ENV, "PYTHONHOME")
        prefix = pyconfigvar(python, "prefix")
        # PYTHONHOME tells python where to look for both pure python and binary modules.
        # When it is set, it replaces both `prefix` and `exec_prefix` and we thus need to set it to both
        # in case they differ. This is also what the documentation recommends.
        ENV["PYTHONHOME"] = prefix * ":" * exec_prefix
    end
    # TODO: look in python-config output? pyconfigvar("LDFLAGS")?
    for libpath in libpaths
        if isfile(joinpath(libpath, lib))
            try
                return dlopen(joinpath(libpath, lib),
                              RTLD_LAZY|RTLD_DEEPBIND|RTLD_GLOBAL)
            end
        end
    end
    lib = splitext(lib)[1]
    try
        # We do this *last* because the libpython in the system
        # library path might be the wrong one if multiple python
        # versions are installed; we prefer the one in LIBDIR.
        return dlopen(lib, RTLD_LAZY|RTLD_DEEPBIND|RTLD_GLOBAL)
    catch
        error("Couldn't find libpython ($lib); try pyinitialize(\"/path/to/$lib\")")
    end
end

# Python 3.x uses wchar_t arrays for some string arguments
function wbytestring(s::String)
    if pyversion.major < 3 || sizeof(Cwchar_t) == 1 # ASCII (or UTF8)
        bytestring(s)
    else # UTF16 or UTF32, presumably
        n = length(s)
        w = Array(Cwchar_t, n + 1)
        if sizeof(Cwchar_t) == 4 # UTF32
            for i = 1:n
                w[i] = s[i]
            end
        else # UTF16, presumably
            @assert sizeof(Cwchar_t) == 2
            if defined(UTF16String) # Julia 0.3
                s16 = utf16(s)
                for i = 1:n
                    w[i] = s16.data[i]
                end
            else
                for i = 1:n
                    # punt on multi-byte encodings
                    s[i] > 0xffff && error("unsupported char in \"$s\"")
                    w[i] = s[i]
                end
            end
        end
        w[n+1] = 0
        w
    end
end

# initialize the Python interpreter (no-op on subsequent calls)
function pyinitialize(python::String)
    global initialized
    if !initialized::Bool
        libpy = try
            lib = dlopen_libpython(python)
            global pyprogramname = wbytestring(pysys(python, "executable"))
            lib
        catch
            # perhaps we were passed library name and not executable?
            lib = dlopen(python, RTLD_LAZY|RTLD_DEEPBIND|RTLD_GLOBAL)
            global pyprogramname = wbytestring(python)
            lib
        end
        pyinitialize(libpy)
    end
    return
end

pyinitialize() = pyinitialize("python")
dlopen_libpython() = dlopen_libpython("python")

# end the Python interpreter and free associated memory
function pyfinalize()
    global initialized
    global finalized
    global libpython
    global inspect
    global builtin
    global BuiltinFunctionType
    global TypeType
    global MethodType
    global MethodWrapperType
    global ufuncType
    global pyxrange
    global c_void_p_Type
    global format_traceback
    global mpmath
    global mpf
    global mpc
    if initialized::Bool
        pygui_stop_all()
        pydecref(mpc::PyObject)
        pydecref(mpf::PyObject)
        pydecref(mpmath::PyObject)
        pydecref(ufuncType::PyObject)
        npyfinalize()
        pydecref(format_traceback::PyObject)
        pydecref(c_void_p_Type::PyObject)
        pydecref(pyxrange::PyObject)
        pydecref(BuiltinFunctionType::PyObject)
        pydecref(TypeType::PyObject)
        pydecref(MethodType::PyObject)
        pydecref(MethodWrapperType::PyObject)
        pydecref(builtin::PyObject)
        pydecref(inspect::PyObject)
        pyexc_finalize()
        pycallback_finalize()
        pyio_finalize()
        pygc_finalize()
        gc() # collect/decref any remaining PyObjects
        ccall((@pysym :Py_Finalize), Void, ())
        dlclose(libpython::Ptr{Void})
        libpython::Ptr{Void} = C_NULL
        initialized::Bool = false
        finalized::Bool = true
    end
    return
end

#########################################################################

include("exception.jl")
include("gui.jl")

#########################################################################

include("gc.jl")

# make a PyObject that embeds a reference to keep, to prevent Julia
# from garbage-collecting keep until o is finalized.
PyObject(o::PyPtr, keep::Any) = pyembed(PyObject(o), keep)

#########################################################################

include("conversions.jl")

include("pytype.jl")

include("callback.jl")

include("io.jl")

#########################################################################
# Pretty-printing PyObject

function pystring(o::PyObject)
    if o.o == C_NULL
        return "NULL"
    else
        s = ccall((@pysym :PyObject_Repr), PyPtr, (PyPtr,), o)
        if (s == C_NULL)
            pyerr_clear()
            s = ccall((@pysym :PyObject_Str), PyPtr, (PyPtr,), o)
            if (s == C_NULL)
                pyerr_clear()
                return string(o.o)
            end
        end
        return convert(String, PyObject(s))
    end
end

function show(io::IO, o::PyObject)
    print(io, "PyObject $(pystring(o))")
end

function help(o::PyObject)
    pycall(pybuiltin("help"), PyAny, o)
end

#########################################################################
# computing hashes of PyObjects

pysalt = hash("PyCall.PyObject") # "salt" to mix in to PyObject hashes

function hash(o::PyObject)
    if o.o == C_NULL
        bitmix(pysalt::Uint, hash(C_NULL))
    elseif is_pyjlwrap(o)
        # call native Julia hash directly on wrapped Julia objects,
        # since on 64-bit Windows the Python 2.x hash is only 32 bits
        bitmix(pysalt::Uint, hash(unsafe_pyjlwrap_to_objref(o.o)))
    else
        h = pyhashlong::Bool ? # changed to Py_hash_t in Python 3.2
               ccall((@pysym :PyObject_Hash), Clong, (PyPtr,), o) :
               ccall((@pysym :PyObject_Hash), Int, (PyPtr,), o)
        if h == -1 # error
            pyerr_clear()
            return bitmix(pysalt::Uint, hash(o.o))
        end
        bitmix(pysalt::Uint, uint(h))
    end
end

#########################################################################
# PyObject equality

function isequal(o1::PyObject, o2::PyObject)
    # TODO: use PyObject_RichCompareBool?  Should
    #       we also add a __cmp__ method for pyjlwrap objects?
    return o1.o == o2.o
end

#########################################################################
# For o::PyObject, make o["foo"] and o[:foo] equivalent to o.foo in Python,
# with the former returning an raw PyObject and the latter giving the PyAny
# conversion.

function getindex(o::PyObject, s::String)
    if (o.o == C_NULL)
        throw(ArgumentError("ref of NULL PyObject"))
    end
    p = ccall((@pysym :PyObject_GetAttrString), PyPtr,
              (PyPtr, Ptr{Uint8}), o, bytestring(s))
    if p == C_NULL
        pyerr_clear()
        throw(KeyError(s))
    end
    return PyObject(p)
end

getindex(o::PyObject, s::Symbol) = convert(PyAny, getindex(o, string(s)))

function setindex!(o::PyObject, v, s::String)
    if (o.o == C_NULL)
        throw(ArgumentError("assign of NULL PyObject"))
    end
    if -1 == ccall((@pysym :PyObject_SetAttrString), Cint,
                   (PyPtr, Ptr{Uint8}, PyPtr), o, bytestring(s), PyObject(v))
        pyerr_clear()
        throw(KeyError(s))
    end
    o
end

setindex!(o::PyObject, v, s::Symbol) = setindex!(o, v, string(s))

function haskey(o::PyObject, s::String)
    if (o.o == C_NULL)
        throw(ArgumentError("haskey of NULL PyObject"))
    end
    return 1 == ccall((@pysym :PyObject_HasAttrString), Cint,
                      (PyPtr, Ptr{Uint8}), o, bytestring(s))
end

haskey(o::PyObject, s::Symbol) = haskey(o, string(s))

#########################################################################
# keys(o) should return an iterator over the keys (members) of o

type PyObjectMembers
    members::PyVector{(Symbol,PyObject)}
end

function show(io::IO, m::PyObjectMembers)
    print(io, "PyObjectMembers iterator (", length(m), " members)")
end

keys(o::PyObject) = PyObjectMembers(pycall(inspect["getmembers"],
                                           PyVector{(Symbol,PyObject)}, o))

start(m::PyObjectMembers) = 1
done(m::PyObjectMembers, i::Int) = i > length(m.members)
next(m::PyObjectMembers, i::Int) = (m.members[i][1], i+1)
length(m::PyObjectMembers) = length(m.members)

#########################################################################
# Create anonymous composite w = pywrap(o) wrapping the object o
# and providing access to o's members (converted to PyAny) as w.member.

# we skip wrapping Julia reserved words (which cannot be type members)
const reserved = Set("while", "if", "for", "try", "return", "break",
                     "continue", "function", "macro", "quote", "let", "local",
                     "global", "const", "abstract", "typealias", "type",
                     "bitstype", "immutable", "ccall", "do", "module",
                     "baremodule", "using", "import", "export", "importall",
                     "pymember", "false", "true")

function pywrap(o::PyObject, mname::Symbol=:__anon__)
    @pyinitialize
    members = convert(Vector{(String,PyObject)},
                      pycall(inspect["getmembers"], PyObject, o))
    filter!(m -> !(m[1] in reserved), members)
    m = Module(mname)
    consts = [Expr(:const, Expr(:(=), symbol(x[1]), convert(PyAny, x[2]))) for x in members]
    exports = try
                  convert(Vector{Symbol}, o["__all__"])
              catch
                  [symbol(x[1]) for x in filter(x -> x[1][1] != '_', members)]
              end
    eval(m, Expr(:toplevel, consts..., :(pymember(s) = getindex($(o), s)),
                 Expr(:export, exports...)))
    m
end

#########################################################################

pyimport(name::String) =
    PyObject(@pycheckn ccall((@pysym :PyImport_ImportModule), PyPtr,
                             (Ptr{Uint8},), bytestring(name)))

pyimport(name::Symbol) = pyimport(string(name))

# convert expressions like :math or :(scipy.special) into module name strings
modulename(s::Symbol) = string(s)
function modulename(e::Expr)
    if e.head == :.
        string(modulename(e.args[1]), :., modulename(e.args[2]))
    elseif e.head == :quote
        modulename(e.args...)
    else
        throw(ArgumentError("invalid module"))
    end
end

# separate this function in order to make it easier to write more
# pyimport-like functions
pyimport_name(name, optional_varname) =
    let len = length(optional_varname)
        len > 0 && (len != 2 || optional_varname[1] != :as) ?
        throw(ArgumentError("usage @pyimport module [as name]")) :
        (len == 2 ? optional_varname[2] :
         typeof(name) == Symbol ? name :
         throw(ArgumentError("$mname is not a valid module variable name, use @pyimport $mname as <name>")))
    end

macro pyimport(name, optional_varname...)
    mname = modulename(name)
    Name = pyimport_name(name, optional_varname)
    quote
        if !isdefined($(Expr(:quote, Name)))
            const $(esc(Name)) = pywrap(pyimport($mname))
        elseif !isa($(esc(Name)), Module)
            error("@pyimport: ", $(Expr(:quote, Name)), " already defined")
        end
        nothing
    end
end

#########################################################################

# look up a global builtin
function pybuiltin(name)
    global builtin
    if (builtin::PyObject).o == C_NULL
        builtin::PyObject = try
            pyimport("__builtin__")
        catch
            pyimport("builtins") # renamed in Python 3
        end
    end
    (builtin::PyObject)[name]
end

#########################################################################

typealias TypeTuple Union(Type,NTuple{Type})

function pycall(o::PyObject, returntype::TypeTuple, args...; kwargs...)
    oargs = map(PyObject, args)
    nargs = length(args)
    sigatomic_begin()
    try
        arg = PyObject(@pycheckn ccall((@pysym :PyTuple_New), PyPtr, (Int,),
                                       nargs))
        for i = 1:nargs
            @pycheckzi ccall((@pysym :PyTuple_SetItem), Cint,
                             (PyPtr,Int,PyPtr), arg, i-1, oargs[i])
            pyincref(oargs[i]) # PyTuple_SetItem steals the reference
        end
        if isempty(kwargs)
            ret = PyObject(@pycheckni ccall((@pysym :PyObject_Call), PyPtr,
                                          (PyPtr,PyPtr,PyPtr), o, arg, C_NULL))
        else
            kw = PyObject((String=>Any)[string(k) => v for (k, v) in kwargs])
            ret = PyObject(@pycheckni ccall((@pysym :PyObject_Call), PyPtr,
                                            (PyPtr,PyPtr,PyPtr), o, arg, kw))
        end
        jret = convert(returntype, ret)
        return jret
    finally
        sigatomic_end()
    end
end

#########################################################################
# Once Julia lets us overload ".", we will use [] to access items, but
# for now we can define "get"

function get(o::PyObject, returntype::TypeTuple, k, default)
    r = ccall((@pysym :PyObject_GetItem), PyPtr, (PyPtr,PyPtr), o,PyObject(k))
    if r == C_NULL
        pyerr_clear()
        default
    else
        convert(returntype, PyObject(r))
    end
end

get(o::PyObject, returntype::TypeTuple, k) =
    convert(returntype, PyObject(@pycheckni ccall((@pysym :PyObject_GetItem),
                                 PyPtr, (PyPtr,PyPtr), o, PyObject(k))))

get(o::PyObject, k, default) = get(o, PyAny, k, default)
get(o::PyObject, k) = get(o, PyAny, k)

function delete!(o::PyObject, k)
    e = ccall((@pysym :PyObject_DelItem), Cint, (PyPtr, PyPtr), o, PyObject(k))
    if e == -1
        pyerr_clear() # delete! ignores errors in Julia
    end
    return o
end

function set!(o::PyObject, k, v)
    @pycheckzi ccall((@pysym :PyObject_SetItem), Cint, (PyPtr, PyPtr, PyPtr),
                     o, PyObject(k), PyObject(v))
    v
end

#########################################################################

const Py_eval_input = 258 # from Python.h
const pyeval_fname = bytestring("PyCall.jl") # filename for pyeval

# evaluate a python string, returning PyObject, given a dictionary
# (string/symbol => value) of local variables to use in the expression
function pyeval_(s::String, locals::PyDict)
    sb = bytestring(s) # use temp var to prevent gc before we are done with o
    sigatomic_begin()
    try
        o = PyObject(@pycheckn ccall((@pysym :Py_CompileString), PyPtr,
                                     (Ptr{Uint8}, Ptr{Uint8}, Cint),
                                     sb, pyeval_fname, Py_eval_input))
        main = @pycheckni ccall((@pysym :PyImport_AddModule),
                                PyPtr, (Ptr{Uint8},),
                                bytestring("__main__"))
        maindict = @pycheckni ccall((@pysym :PyModule_GetDict), PyPtr,
                                    (PyPtr,), main)
        return PyObject(@pycheckni ccall((@pysym :PyEval_EvalCode),
                                         PyPtr, (PyPtr, PyPtr, PyPtr),
                                         o, maindict, locals))
    finally
        sigatomic_end()
    end
end

function pyeval(s::String, returntype::TypeTuple=PyAny; kwargs...)
    locals = PyDict{String,PyObject}()
    for (k, v) in kwargs
        locals[string(k)] = v
    end
    return convert(returntype, pyeval_(s, locals))
end

#########################################################################

end # module PyCall


text:
using .VegaLite

export DensityShader, CurrentShader

#######################################################################
# Shaders
#######################################################################
struct CurrentShader{K,T}
    axis::Int
    kernel::K
    transform::T
end

"""
    CurrentShader(axis = 0; kernel = 1, transform = identity)

Construct a `CurrentShader` object that can be used in `vlplot(h, psi; ...)` to visualize
the current of a state `psi` along a given `axis` under Hamiltonian `h`, optionally
transformed by `transform`. An `axis = 0` represent the norm of the current.

The current along a link `i=>k` is defined as

    j_ki = (r[k][axis]-r[i][axis]) * imag(psi[k]' * h[k,i] * kernel * psi[j])

The current at a site `i` is the sum of links incident on that site

    j_i = ∑_k j_ki

The visualized current is `transform(j_ki)` and `transform(j_i)`, respectively.
"""
CurrentShader(axis = 0; kernel = 1, transform = identity) =
    CurrentShader(axis, kernel, transform)

struct DensityShader{K,T}
    kernel::K
    transform::T
end

"""
    DensityShader(kernel = 1, transform = identity)

Construct a `DensityShader` object that can be used in `vlplot(h, psi; ...)` to visualize
the density of a state `psi`, optionally transformed by `transform`.

The visualized density at a site `i` is defined as

    ρ_i = transform(psi[i]' * kernel * psi[i])

The density at a half-link is equal to the density at the originating site.
"""
DensityShader(; kernel = 1, transform = identity) = DensityShader(kernel, transform)

site_shader(shader, h, psi) = i -> applyshader(shader, h, psi, i)
link_shader(shader, h, psi) = (i, j, dn) -> applyshader(shader, h, psi, i, j, dn)

applyshader(shader, h, psi, i...) = applyshader_vector(shader, h, psi, i...)
applyshader(shader, h, psi::AbstractMatrix, i...) = sum(col -> applyshader_vector(shader, h, view(psi, :, col), i...), axes(psi, 2))

applyshader_vector(shader::Number, h, x...) = shader
applyshader_vector(shader::Missing, h, x...) = 0.0
applyshader_vector(shader::Function, h, psi, i) = shader(psi[i])
applyshader_vector(shader::Function, h, psi, i, j, dn) = shader(psi[i], psi[j])
applyshader_vector(shader::DensityShader, h, psi, i) = real(shader.transform(psi[i]' * shader.kernel * psi[i]))
applyshader_vector(shader::DensityShader, h, psi, i, j, dn) = real(shader.transform(psi[j]' * shader.kernel * psi[j]))

function applyshader_vector(shader::CurrentShader, h, psi::AbstractVector{T}, src) where {T}
    pos = allsitepositions(h.lattice)
    c = zero(real(eltype(T)))
    for har in h.harmonics
        rows = rowvals(har.h)
        for ptr in nzrange(har.h, src)
            dst = rows[ptr]
            c += applyshader_vector(shader, h, psi, dst, src, har.dn)
        end
    end
    return c
end

function applyshader_vector(shader::CurrentShader, h, psi, i, j, dn)
    pos = allsitepositions(h.lattice)
    # Assuming that only links from base unit cell are plotted
    dr = pos[i] - pos[j] + bravais(h) * dn
    dc = imag(psi[i]' * shader.kernel * h[dn, i, j] * psi[j])
    c = iszero(shader.axis) ? shader.transform(norm(dr * dc)) : shader.transform(dr[shader.axis] * dc)
    return real(c)
end

#######################################################################
# vlplot
#######################################################################
"""
    vlplot(b::Bandstructure{1}; kw...)

Plot the 1D bandstructure `b` using VegaLite in 2D.

    vlplot(h::Hamiltonian; kw...)

Plot the the Hamiltonian lattice projected along `axes` using VegaLite.

    vlplot(h::Hamiltonian, psi::AbstractVector; kw...)

Plot an eigenstate `psi` on the lattice, using one of various possible visualization
shaders: `sitesize`, `siteopacity`, `sitecolor`, `linksize`, `linkopacity` or `linkcolor`
(see keywords below). If `psi` is obtained in a flattened form, it will be automatically
"unflattened" to restore the orbital structure of `h`. The vector `psi` must have length
equal to the number of sites (after unflattening if required).

If `h` is defined on an unbounded `L`-dimensional lattice, and `psi` is of the form

    psi = [cell₁ => psi₁, cell₂ => psi₂,...]

where `cellᵢ::NTuple{L,Int}` are indices of unit cells, said unit cells will all be plotted
together with their corresponding `psiᵢ`. Otherwise, a single cell at the origin will be
plotted.

    vlplot(h::Hamiltonian, psi::AbstractMatrix; kw...)

Same as above, but columns of psi will be summed over after applying shaders to each.

    vlplot(h::Hamiltonian, psi::Subspace; kw...)

Equivalent to `vlplot(h, psi.basis; kw...)`.

# Keyword arguments and defaults

## Common
    - `size = 800`: the `(width, height)` of the plot (or `max(width, height)` if a single number)
    - `labels`: labels for the x and y plot axes. Defaults `("φ/2π", "ε")` and `("x", "y")` respectively
    - `scaling = (1/2π, 1)`: `(scalex, scaley)` scalings for the x (Bloch phase) and y (energy) variables
    - `xlims = missing` and `ylims = missing`: `(xmin, xmax)` and `(ymin, ymax)` to constrain plot range

## Bandstructure-specifc
    - `points = false`: whether to plot points on line plots
    - `bands = missing`: bands to plot (or all if `missing`)
    - `thickness = 2`: thickness of band lines

## Hamiltonian-specific
    - `axes = (1,2)`:lattice axes to project onto the plot x-y plane
    - `digits = 4`: number of significant digits to show in onsite energy and hopping tooltips
    - `plotsites = true`: whether to plot sites
    - `plotlinks = true`: whether to plot links
    - `onlyonecell = false`: whether to omit plotting links to neighboring unit cells in unbounded lattices
    - `sitestroke = :white`: the color of site outlines. If `nothing`, no outline will be plotted.
    - `colorscheme = "redyellowblue"`: Color scheme from `https://vega.github.io/vega/docs/schemes/`
    - `discretecolorscheme = "category10"`: Color scheme from `https://vega.github.io/vega/docs/schemes/`
    - `maxdiameter = 15`: maximum diameter of sites, useful when `sitesize` is not a constant.
    - `mindiameter = 0`: site diameters below this value will be excluded (may be used with `plotlinks = false` to increase performance)
    - `maxthickness = 0.25`: maximum thickness of links, as a fraction of `maxdiameter`, useful when `linksize` is not a constant.
    - `sitesize = maxdiameter`: diameter of sites.
    - `siteopacity = 0.9`: opacity of sites.
    - `sitecolor = missing`: color of sites. If `missing`, colors will encode the site sublattice following `discretecolorscheme`
    - `linksize = maxthickness`: thickness of hopping links.
    - `linkopacity = 1.0`: opacity of hopping links.
    - `linkcolor = missing`: color of links. If `missing`, colors will encode source site sublattice, following `discretecolorscheme`
    - `colorrange = missing`: the range of values to encode using `sitecolor` and `linkcolor`. If `missing` the minimum and maximum values will be used.

Apart from a number, all site shaders (`sitesize`, `siteopacity` and `sitecolor`) can also
be a real-valued function of the wavefunction amplitude `psi` at each site, `psi -> value(psi)`.
Similarly, link shaders (`linksize`, `linkopacity` and `linkcolor`) can be a function of the
wavefunction amplitudes at end and origin of the link, `(psi´, psi) -> value(psi´, psi)`.

Any shader can also be a `CurrentShader` or a `DensityShader` object, that computes the
current or density at each site or link (see `CurrentShader` and `DensityShader`).

Note that values of sizes are relative, so there is no need to normalize shader values. The
maximum visual value of site diameter and link thickness is given by `maxdiameter` (in
pixels) and `maxthickness` (as a fraction of `maxdiameter`)
"""
function VegaLite.vlplot(b::Bandstructure;
    labels = ("φ", "ε"), scaling = (1, 1), size = 640, points = false, xlims = missing, ylims = missing, bands = missing, thickness = 2)
    labelx, labely = labels
    table = bandtable(b, make_it_two(scaling), bands)
    sizes = make_it_two(size)
    corners = _corners(table)
    plotrange = (xlims, ylims)
    (domainx, domainy), _ = domain_size(corners, size, plotrange)
    p = table |> vltheme(sizes, points) + @vlplot(
        mark = {:line, strokeWidth = thickness},
        x = {:x, scale = {domain = domainx, nice = false}, title = labelx, sort = nothing},
        y = {:y, scale = {domain = domainy, nice = false}, title = labely, sort = nothing},
        color = "band:n",
        selection = {grid = {type = :interval, bind = :scales}})
    return p
end

# We optimize 1D band plots by collecting connected simplices into line segments (because :rule is considerabley slower than :line)
function bandtable(b::Bandstructure{1,C,T}, (scalingx, scalingy), bandsiter) where {C,T}
    bandsiter´ = bandsiter === missing ? eachindex(bands(b)) : bandsiter
    NT = typeof((;x = zero(T), y = zero(T), band = 1, tooltip = 1))
    table = NT[]
    for nb in bandsiter´
        band = bands(b, nb)
        verts = vertices(band)
        simps = band.simps
        isempty(simps) && continue
        s0 = (0, 0)
        for s in simps
            if first(s) == last(s0) || s0 == (0, 0)
                push!(table, (; x = verts[first(s)][1] * scalingx, y = verts[first(s)][2] * scalingy, band = nb, tooltip = degeneracy(band, first(s))))
            else
                push!(table, (; x = verts[last(s0)][1] * scalingx, y = verts[last(s0)][2] * scalingy, band = nb, tooltip = degeneracy(band, last(s0))))
                push!(table, (; x = T(NaN), y = T(NaN), band = nb, tooltip = 0))  # cut
                push!(table, (; x = verts[first(s)][1] * scalingx, y = verts[first(s)][2] * scalingy, band = nb, tooltip = degeneracy(band, first(s))))
            end
            s0 = s
        end
        push!(table, (; x = verts[last(s0)][1] * scalingx, y = verts[last(s0)][2] * scalingy, band = nb, tooltip = degeneracy(band, last(s0))))
    end
    return table
end

function VegaLite.vlplot(h::Hamiltonian{LA}, psi = missing;
                         labels = ("x","y"), size = 800, axes::Tuple{Int,Int} = (1,2), xlims = missing, ylims = missing, digits = 4,
                         onlyonecell = false, plotsites = true, plotlinks = true,
                         maxdiameter = 20, mindiameter = 0, maxthickness = 0.25,
                         sitestroke = :white, sitesize = maxdiameter, siteopacity = 0.9, sitecolor = missing,
                         linksize = maxthickness, linkopacity = 1.0, linkcolor = missing,
                         colorrange = missing,
                         colorscheme = "redyellowblue", discretecolorscheme = "category10") where {E,LA<:Lattice{E}}

    table = linkstable(h, psi; axes, digits, onlyonecell, plotsites, plotlinks,
                               sitesize, siteopacity, sitecolor, linksize, linkopacity, linkcolor)
    maxsitesize    = plotsites ? maximum(s.scale for s in table if !s.islink)   : 0.0
    maxlinksize    = plotlinks ? maximum(s.scale for s in table if s.islink)    : 0.0

    mindiameter > 0 && filter!(s -> s.islink || s.scale >= mindiameter*maxsitesize/maxdiameter, table)

    corners    = _corners(table)
    plotrange  = (xlims, ylims)
    (domainx, domainy), sizes = domain_size(corners, size, plotrange)

    if (plotsites && sitecolor !== missing) || (plotlinks && linkcolor !== missing)
        colorfield = :color
        colorscheme´ = colorscheme
        colorrange´ = sanitize_colorrange(colorrange, table)
    else
        colorfield = :sublat
        colorscheme´ = discretecolorscheme
        colorrange´ = nothing
    end

    p = vltheme(sizes)
    if plotlinks
        p += @vlplot(
            mark = {:rule},
            strokeWidth = {:scale,
                scale = {range = (0, maxthickness * maxdiameter), domain = (0, maxlinksize)},
                legend = needslegend(linksize)},
            color = {colorfield,
                scale = {domain = colorrange´, scheme = colorscheme´},
                legend = needslegend(linkcolor)},
            strokeOpacity = {:opacity,
                scale = {range = (0, 1), domain = (0, 1)},
                legend = needslegend(linkopacity)},
            transform = [{filter = "datum.islink"}],
            selection = {grid2 = {type = :interval, bind = :scales}},
            encoding = {
                x = {:x, scale = {domain = domainx, nice = false}, axis = {grid = false}, title = labels[1]},
                y = {:y, scale = {domain = domainy, nice = false}, axis = {grid = false}, title = labels[2]},
                x2 = {:x2, scale = {domain = domainx, nice = false}, axis = {grid = false}},
                y2 = {:y2, scale = {domain = domainy, nice = false}, axis = {grid = false}}})
    end
    if plotsites
        p += @vlplot(
            mark = {:circle, stroke = sitestroke},
            size = {:scale,
                scale = {range = (0, maxdiameter^2), domain = (0, maxsitesize), clamp = false},
                legend = needslegend(sitesize)},
            color = {colorfield,
                scale = {domain = colorrange´, scheme = colorscheme´},
                legend = needslegend(sitecolor)},
            opacity = {:opacity,
                scale = {range = (0, 1), domain = (0, 1)},
                legend = needslegend(siteopacity)},
            selection = {grid1 = {type = :interval, bind = :scales}},
            transform = [{filter = "!datum.islink"}],
            encoding = {
                x = {:x, scale = {domain = domainx, nice = false}, axis = {grid = false}, title = labels[1]},
                y = {:y, scale = {domain = domainy, nice = false}, axis = {grid = false}, title = labels[2]}
            })
    end
    return table |> p
end

function vltheme((sizex, sizey), points = false)
    p = @vlplot(
        tooltip = :tooltip,
        width = sizex, height = sizey,
        config = {
            circle = {strokeWidth = 1, size = 200},
            line = {point = points},
            scale = {minOpacity = 0, maxOpacity = 1.0}})
    return p
end

sanitize_colorrange(::Missing, table) = extrema(s -> s.color, table)
sanitize_colorrange(r::AbstractVector, table) = convert(Vector, r)
sanitize_colorrange(r::Tuple, table) = [first(r), last(r)]

needslegend(x::Number) = nothing
needslegend(x) = true

unflatten_orbitals_or_reinterpret_or_missing(psi::Missing, h) = missing
unflatten_orbitals_or_reinterpret_or_missing(psi::AbstractArray, h) = unflatten_orbitals_or_reinterpret(psi, orbitalstructure(h))
unflatten_orbitals_or_reinterpret_or_missing(s::Subspace, h) = unflatten_orbitals_or_reinterpret(s.basis, orbitalstructure(h))

checkdims_psi(h, psi) = size(h, 2) == size(psi, 1) ||
    throw(ArgumentError("The eigenstate length $(size(psi,1)) must match the Hamiltonian dimension $(size(h, 2))"))
checkdims_psi(h, ::Missing) = nothing

function linkstable(h, psi; kw...)
    T = numbertype(h.lattice)
    table = NamedTuple{(:x, :y, :x2, :y2, :sublat, :tooltip, :scale, :color, :opacity, :islink),
                       Tuple{T,T,T,T,NameType,String,Float64,Float64,Float64,Bool}}[]
    linkstable!(table, h, psi; kw...)
    return table
end

function linkstable!(table, h, cellpsi::AbstractVector{<:Pair}; kw...)
    foreach(cellpsi) do (cell, psi)
        linkstable!(table, h, psi, cell; kw..., onlyonecell = true)
    end
    return table
end

function linkstable!(table, h, psi, cell = missing;
                     axes, digits, onlyonecell, plotsites, plotlinks,
                     sitesize, siteopacity, sitecolor, linksize, linkopacity, linkcolor) where {LA,L}
    psi´ = unflatten_orbitals_or_reinterpret_or_missing(psi, h)
    checkdims_psi(h, psi´)
    d = (; axes = axes, digits = digits,
        sitesize_shader    = site_shader(sitesize, h, psi´),
        siteopacity_shader = site_shader(siteopacity, h, psi´),
        sitecolor_shader   = site_shader(sitecolor, h, psi´),
        linksize_shader    = link_shader(linksize, h, psi´),
        linkopacity_shader = link_shader(linkopacity, h, psi´),
        linkcolor_shader   = link_shader(linkcolor, h, psi´))
    (a1, a2) = d.axes
    lat = h.lattice
    br = bravais(h)
    T = numbertype(lat)
    r0 = cell === missing ? br * zero(SVector{latdim(h),T}) : br * toSVector(T,cell)
    slats = sublats(lat)
    rs = allsitepositions(lat)
    h0 = h.harmonics[1].h
    rows = Int[] # list of plotted destination sites for dn != 0

    for har in h.harmonics
        !iszero(har.dn) && onlyonecell && continue
        resize!(rows, 0)
        ridx = 0
        for ssrc in slats
            if iszero(har.dn)
                for (i, r) in enumeratesites(lat, ssrc)
                    x = get(r + r0, a1, zero(T)); y = get(r + r0, a2, zero(T))
                    ridx += 1
                    plotsites && push!(table, (x = x, y = y, x2 = x, y2 = y,
                                sublat = sublatname(lat, ssrc), tooltip = matrixstring_inline(i, h0[i, i], d.digits),
                                scale = d.sitesize_shader(ridx), color = d.sitecolor_shader(ridx),
                                opacity = d.siteopacity_shader(ridx), islink = false))
                end
            end
            for sdst in slats
                itr = nonzero_indices(har, siterange(lat, sdst), siterange(lat, ssrc))
                for (row, col) in itr
                    iszero(har.dn) && col == row && continue
                    rsrc = rs[col] + r0
                    rdst = (rs[row] + bravais(lat) * har.dn + r0)
                    # sites in neighboring cells connected to dn=0 cell. But add only once.
                    if !iszero(har.dn) && !in(row, rows)
                        x = get(rdst, a1, zero(T)); y = get(rdst, a2, zero(T))
                        plotsites && push!(table,
                            (x = x, y = y, x2 = x, y2 = y,
                            sublat = sublatname(lat, sdst), tooltip = matrixstring_inline(row, h0[row, row], d.digits),
                            scale =  d.sitesize_shader(row), color = d.sitecolor_shader(row),
                            opacity = 0.5 * d.siteopacity_shader(row), islink = false))
                        push!(rows, row)
                    end
                    plotlinks || continue
                    # draw half-links but only intracell
                    rdst = ifelse(cell !== missing || iszero(har.dn), (rdst + rsrc) / 2, rdst)
                    x  = get(rsrc, a1, zero(T)); y  = get(rsrc, a2, zero(T))
                    x´ = get(rdst, a1, zero(T)); y´ = get(rdst, a2, zero(T))
                    # Exclude links perpendicular to the screen
                    rdst ≈ rsrc || push!(table,
                        (x = x, y = y, x2 = x´, y2 = y´,
                        sublat = sublatname(lat, ssrc), tooltip = matrixstring_inline(row, col, har.h[row, col], d.digits),
                        scale = d.linksize_shader(row, col, har.dn), color = d.linkcolor_shader(row, col, har.dn),
                        opacity = ifelse(iszero(har.dn), 1.0, 0.5) * d.linkopacity_shader(row, col, har.dn), islink = true))
                end
            end
        end
    end
    return table
end

function _corners(table)
    min´ = max´ = (first(table).x, first(table).y)
    for row in table
        x, y = row.x, row.y
        (isnan(x) || isnan(y)) && continue
        min´ = min.(min´, (x, y))
        max´ = max.(max´, (x, y))
        if isdefined(row, :x2)
            x2, y2 = row.x2, row.y2
            min´ = min.(min´, (x2, y2))
            max´ = max.(max´, (x2, y2))
        end
    end
    return min´, max´
end

function domain_size(corners, size, (rangex, rangey))
    domainx = iclamp((corners[1][1], corners[2][1]), rangex)
    domainy = iclamp((corners[1][2], corners[2][2]), rangey)
    dx = domainx[2]-domainx[1]
    dy = domainy[2]-domainy[1]
    sizex, sizey = compute_sizes(size, (dx, dy))
    return (domainx, domainy), (sizex, sizey)
end

make_it_two(x::Number) = (x, x)
make_it_two(x::Tuple{Number,Number}) = x

function compute_sizes(size::Number, (dx, dy))
    if dx > dy
        sizex = size
        sizey = size * dy/dx
    else
        sizex = size * dx/dy
        sizey = size
    end
    return sizex, sizey
end

compute_sizes(ss::Tuple{Number,Number}, d) = ss

ext: jl

max_stars_repo_name: pablosanjose/Quantica

max_issues_repo_name: pablosanjose/Quantica

max_forks_repo_name: pablosanjose/Quantica


text:
using DSGE, FileIO, JLD2, ModelConstructors, Test, Random, Dates
path = dirname(@__FILE__)

# Set up arguments
m = AnSchorfheide(testing = true)
m <= Setting(:date_forecast_start, quartertodate("2015-Q4"))
m <= Setting(:forecast_horizons, 1)

system, histshocks = JLD2.jldopen("$path/../reference/forecast_args.jld2","r") do file
    read(file, "system"), read(file, "histshocks")
end

# Read expected output
exp_states, exp_obs, exp_pseudo =
    JLD2.jldopen("$path/../reference/shock_decompositions_out.jld2", "r") do file
        read(file, "exp_states"), read(file, "exp_obs"), read(file, "exp_pseudo")
    end

# With shockdec_startdate not null
states, obs, pseudo = shock_decompositions(m, system, histshocks)

@testset "Test shockdec with non-null startdate" begin
    @test @test_matrix_approx_eq exp_states[:startdate] states
    @test @test_matrix_approx_eq exp_obs[:startdate]    obs
    @test @test_matrix_approx_eq exp_pseudo[:startdate] pseudo
end

# With shockdec_startdate null
#m <= Setting(:shockdec_startdate, Nullable{Date}())
m <= Setting(:shockdec_startdate, nothing)
states, obs, pseudo = shock_decompositions(m, system, histshocks)

@testset "Test shockdec with null startdate" begin
    @test @test_matrix_approx_eq exp_states[:no_startdate] states
    @test @test_matrix_approx_eq exp_obs[:no_startdate]    obs
    @test @test_matrix_approx_eq exp_pseudo[:no_startdate] pseudo
end

@testset "Deterministic trends" begin
    dettrend = DSGE.deterministic_trends(m, system, zeros(8))
    @test @test_matrix_approx_eq dettrend[1] zeros(size(dettrend[1]))
    @test @test_matrix_approx_eq dettrend[2] zeros(size(dettrend[2]))
    @test @test_matrix_approx_eq dettrend[3] zeros(size(dettrend[3]))
end

@testset "Trends" begin
    out = DSGE.trends(system)
    @test @test_matrix_approx_eq out[1] system[:CCC]
    @test @test_matrix_approx_eq out[2] system[:ZZ] * system[:CCC] + system[:DD]
    @test @test_matrix_approx_eq out[3] system[:ZZ_pseudo] * system[:CCC] + system[:DD_pseudo]
end

# Check trends for regime-switching system
reg_sys = RegimeSwitchingSystem(System{Float64}[system, system])
@testset "Regime-switching trends" begin
    out = DSGE.trends(reg_sys)
    for i in 1:2
        @test @test_matrix_approx_eq out[1][:, i] reg_sys[i, :CCC]
        @test @test_matrix_approx_eq out[2][:, i] reg_sys[i, :ZZ] * reg_sys[i, :CCC] + reg_sys[i, :DD]
        @test @test_matrix_approx_eq out[3][:, i] reg_sys[i, :ZZ_pseudo] * reg_sys[i, :CCC] + reg_sys[i, :DD_pseudo]
    end
end

m.settings[:regime_dates] = Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2010, 3, 31)))
m.test_settings[:regime_dates] = Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2010, 3, 31)))
@testset "Regime-switching deterministic trends" begin
    dettrend = DSGE.deterministic_trends(m, reg_sys, zeros(8))
    @test @test_matrix_approx_eq dettrend[1] zeros(size(dettrend[1]))
    @test @test_matrix_approx_eq dettrend[2] zeros(size(dettrend[2]))
    @test @test_matrix_approx_eq dettrend[3] zeros(size(dettrend[3]))
end

# Check shock decompositions for regime-switching system
out_shockdec1 = shock_decompositions(m, reg_sys, histshocks,
                                    date_presample_start(m), date_forecast_start(m)) # check we only return the shockdecs for requested period
out_shockdec2 = shock_decompositions(m, reg_sys, histshocks[:, 1:end - 1]) # check default

@testset "Shock decompositions with regime-switching" begin
    @test @test_matrix_approx_eq out_shockdec1[1] states
    @test @test_matrix_approx_eq out_shockdec1[2] obs
    @test @test_matrix_approx_eq out_shockdec1[3] pseudo
    @test @test_matrix_approx_eq out_shockdec2[1] states
    @test @test_matrix_approx_eq out_shockdec2[2] obs
    @test @test_matrix_approx_eq out_shockdec2[3] pseudo
end

nothing


text:
function x_crosser(r::Robot)
    side=Nord
    for a ∈ 1:4
        while isborder(r,side)==false && isborder(r,next(side))==false #Проверка двух сторон
            move!(r,side)
            move!(r,next(side))
            putmarker!(r)
        end
        side=inverse(side)
        while ismarker(r)
            move!(r,side)
            move!(r,next(side))
        end
        side=inverse(next(side))
    end
    putmarker!(r)
end

function inverse(side::HorizonSide) #Взято из библиотечного файла
    HorizonSide(mod(Int(side)+2,4))
end
function next(side::HorizonSide)
    HorizonSide(mod(Int(side)+1,4))
end
#= Возможна реализация с помощью кортежа #=

max_stars_repo_name: alexgruzberg/mirea-progs

max_issues_repo_name: alexgruzberg/mirea-progs

max_forks_repo_name: alexgruzberg/mirea-progs


text:
using GaussQuadrature
using SpecialFunctions: besseli
using Test

include("test_eigen.jl")
include("test_legendre.jl")
include("test_chebyshev.jl")
include("test_jacobi.jl")
include("test_laguerre.jl")
include("test_hermite.jl")
include("test_logweight.jl")


text:
module Dec8Part2
function get_new_index_compute_accumulator(op, value, current_index, accumulator)
    if op == "jmp"
        current_index += value
    elseif op == "acc"
        accumulator += value
        current_index += 1
    else
        current_index += 1
    end
    return current_index, accumulator
end

function parse_operation(x::String)::Array{Union{String,Int64},1}
    op, value = split(x, ' ')
    return [String(op), parse(Int64, value)]
end

function get_accumulator_run_operations(operations::Array, index::Int64, accumulator::Int64)::Union{Nothing,Int64}
    visited_index = Int64[]
    push!(visited_index, index)
    n_ops = length(operations)
    while index <= n_ops
        op, value = operations[index]
        index, accumulator = get_new_index_compute_accumulator(op, value, index, accumulator)
        if index ∉ visited_index
            push!(visited_index, index)
        else
            return nothing
        end
    end
    if n_ops ∈ visited_index
        return accumulator
    else
        return nothing
    end
end

function correct_get_accumulator(operations::Array, index::Int64, accumulator::Int64)
    for (index_op, operation) in enumerate(operations)
        ops = deepcopy(operations)
        op = operation[1]
        value = operation[2]
        if op == "jmp"
            ops[index_op][1] = "nop"
            acc = get_accumulator_run_operations(ops, index, accumulator)
            if acc === nothing
                continue
            else
                return acc
            end
        elseif op == "nop"
            ops[index_op][1] = "jmp"
            acc = get_accumulator_run_operations(ops, index, accumulator)
            if acc === nothing
                continue
            else
                return acc
            end
        end
    end
end

function main(file_path::String)::Int64
    operations = readlines(file_path)
    current_index = 1
    accumulator = 0
    operations = parse_operation.(operations)
    return correct_get_accumulator(operations, current_index, accumulator)
end
end


text:
## ComponentArrays
# Similar
Base.similar(x::ComponentArray) = ComponentArray(similar(getdata(x)), getaxes(x)...)
Base.similar(x::ComponentArray, ::Type{T}) where T = ComponentArray(similar(getdata(x), T), getaxes(x)...)
Base.similar(x::ComponentArray, ::Type{T}, ax::Tuple{Vararg{Int64,N}}) where {T,N} = similar(x, T, ax...)
Base.similar(x::ComponentArray, ::Type{T}, ax::Union{Integer, Base.OneTo}...) where T =
    similar(getdata(x), T, ax...)
## TODO: write length method for AbstractAxis so we can do this?
    # function Base.similar(::Type{CA}) where CA<:ComponentArray{T,N,A,Axes} where {T,N,A,Axes}
#     axs = getaxes(CA)
#     return ComponentArray(similar(A, length.(axs)...), axs...)
# end

Base.zero(x::ComponentArray) = zero.(x)

## FIXME: waiting on similar(::Type{<:ComponentArray})
# Base.zeros(CA::Type{<:ComponentArray}) = (similar(CA) .= 0)

# Base.ones(CA::Type{<:ComponentArray}) = (similar(CA) .= 1)

# Copy
Base.copy(x::ComponentArray) = ComponentArray(copy(getdata(x)), getaxes(x))

Base.copyto!(dest::AbstractArray, src::ComponentArray) = copyto!(dest, getdata(src))
function Base.copyto!(dest::ComponentArray, src::AbstractArray)
    copyto!(getdata(dest), src)
    return dest
end
function Base.copyto!(dest::ComponentArray, src::ComponentArray)
    copyto!(getdata(dest), getdata(src))
    return dest
end

Base.deepcopy(x::ComponentArray) = ComponentArray(deepcopy(getdata(x)), getaxes(x))

Base.convert(::Type{CA}, A::AbstractArray) where CA<:ComponentArray = ComponentArray(A, getaxes(CA))
Base.convert(::Type{CA}, x::ComponentArray) where CA<:ComponentArray = ComponentArray(getdata(x), getaxes(CA))
Base.convert(::Type{<:Array}, x::ComponentArray) = convert(Array, getdata(x))

# Conversion to from ComponentArray to NamedTuple (note, does not preserve numeric types of
# original NamedTuple)
function _namedtuple(x::ComponentVector)
    data = []
    idxmap = indexmap(getaxes(x)[1])
    for key in valkeys(x)
        idx = idxmap[getval(key)]
        # if idx isa AliasAxis
        #     val = idx.f
        # else
            val = getproperty(x, key) |> _namedtuple
        # end
        push!(data, getval(key) => val)
    end
    return (; data...)
end
_namedtuple(v::AbstractVector) = _namedtuple.(v)
_namedtuple(x) = x

Base.convert(::Type{NamedTuple}, x::ComponentVector) = _namedtuple(x)
Base.NamedTuple(x::ComponentVector) = _namedtuple(x)

## AbstractAxis conversion and promotion
# Base.convert(TypeAx::Type{<:Ax1}, ::Ax2) where {Ax1<:AbstractAxis,Ax2<:AbstractAxis} = promote_type(TypeAx,Ax2)()
# Base.convert(::Type{Ax1}, ax::Ax2) where {Ax1<:AbstractAxis,Ax2<:AbstractAxis} = Ax1()

# Could not figure out promote_rule for these.
Base.promote(ax::AbstractAxis, ::NullAxis) = ax
Base.promote(ax::AbstractAxis, ::FlatAxis) = ax
Base.promote(::NullAxis, ax::AbstractAxis) = ax
Base.promote(::FlatAxis, ax::AbstractAxis) = ax
Base.promote(::FlatAxis, ::NullAxis) = FlatAxis()
Base.promote(::NullAxis, ::FlatAxis) = FlatAxis()
Base.promote(::FlatAxis, ::FlatAxis) = FlatAxis()
Base.promote(::NullAxis, ::NullAxis) = NullAxis()
Base.promote(ax::Ax, ::Ax) where {Ax<:AbstractAxis} = ax
function Base.promote(ax1::AbstractAxis, ax2::AbstractAxis)
    @warn "Tried to promote dissimilar axes of types $(typeof(ax1)) and $(typeof(ax2)). Falling back to FlatAxis."
    return FlatAxis()
end


text:
export supnorm

"""
    refine(tt, n)

Refine range by decreasing stepsize by a factor `n`.
"""
refine(tt, n) =  first(tt):(Base.step(tt)/n):last(tt)

"""
    rescale(x, a=>b, u=>v)

Linearly map the interval [a,b] to [u,v].

Example:

    rescale.(x, Ref(extrema(x)))
"""
rescale(x, (a,b)::Union{Pair,NTuple{2}}, (u,v)::Union{Pair,NTuple{2}}=(zero(x)=>one(x))) = (x-a)/(b-a)*(v-u) + u

"""
    runmean(x)

Running mean of the vector `x`.
"""
runmean(x, cx = cumsum(x)) = [cx[n]/n for n in 1:length(x)]

function runmean(xx::Matrix)
    yy = copy(xx) / 1
    m = 0 * (copy(xx[1,:])/1)
    for i in 1:size(yy, 1)
        m[:] = m + (xx[i, :] - m)/i
        yy[i, :] = m
    end
    yy
end

"""
    cumsum0(x)

Cumulative sum starting at 0 such that `cumsum0(diff(x)) ≈ x`.
"""
function cumsum0(dx::Vector)
        n = length(dx) + 1
        x = similar(dx, n)
        x[1] = 0.0
        for i in 2:n
                x[i] = x[i-1] + dx[i-1]
        end
        x
end

"""
    Supremum norm
"""
supnorm(x) = maximum(abs.(x))

"""
    outer(x[, y])

Short-hand for quadratic form xx' (or xy').
"""
outer(x) = x*x'
outer(x,y) = x*y'

"""
    inner(x[, y])

Short-hand for quadratic form x'x (or x'y).
"""
inner(x) = dot(x,x)
inner(x,y) = dot(x,y)

"""
    mat(yy::Vector{SVector})

Reinterpret `X` or `yy` to an array without change in memory.
"""
mat(yy::Vector{SVector{d,T}}) where {d,T} = reshape(reinterpret(T, yy), d, length(yy))

unmat(A::Matrix{T}) where {T} = reinterpret(SVector{size(A, 1), T}, A, (size(A, 2),))
unmat(::Type{SVector{d,T}}, A::Matrix{T}) where {d,T} = reinterpret(SVector{d,T}, A, (div(length(A), d),))

"""
    quaternion(m::SMatrix{3,3})

Compute the (rotation-) quarternion of a 3x3 rotation matrix. Useful to create
isodensity ellipses from spheres in GL visualizations.
"""
function quaternion(m)
    qw2 = 1 + m[1,1] + m[2,2] + m[3,3]
    if qw2 < 0
        #if abs(det(m) - 1) < sqrt(eps()) # might fail with svd() for any reasonable values of `eps`
            qw2 = -qw2
        #else
        #    throw(DomainError())
        #end
    end

    qw = sqrt(qw2)/2
    qx = (m[3,2] - m[2,3])/(4qw)
    qy = (m[1,3] - m[3,1])/(4qw)
    qz = (m[2,1] - m[1,2])/(4qw)
    SVector{4}(qx,qy,qz,qw)
end

"""
    piecewise(X::SamplePath, [endtime]) -> tt, xx

If X is a jump process with piecewise constant paths and jumps in `X.tt`,
piecewise returns coordinates path for plotting purposes. The second argument
allows to choose the right endtime of the last interval.
"""
function piecewise(Y::SamplePath, tend = Y.tt[end])
    tt = [Y.tt[1]]
    n = length(Y.yy)
    append!(tt, repeat(Y.tt[2:n], inner=2))
    push!(tt, tend)
    tt, repeat(Y.yy, inner=2)
end

"""
    upsample(x, td, t)

If `x` is piecewise constant with jumps at `td`, return values of `x`
at times `t`.
"""
upsample(x, td, t) = [x[s == td[end] ? end : searchsortedlast(td, s)] for s in t]

"""
    viridis

Map `s` onto the first `maxviri` viridis colors
"""
function viridis(s, alpha = 0.9f0, maxviri = 200, RGBA=RGBA)
    l, u = extrema(s)
    map(x->RGBA(Float32.(Bridge._viridis[1+floor(Int,x*maxviri)])..., alpha), (s .- l)/(u .- l))
end

function iterateall(it)
    u = iterate(it)
    u === nothing && error("Empty iterator")
    x, state = u
    while true
        u = iterate(it, state)
        u === nothing && return x
        x, state = u
    end
end

"""
    _viridis

Color data of the Viridis map by Nathaniel J. Smith, Stefan van Der Walt, Eric Firing
from https://github.com/BIDS/colormap/blob/master/colormaps.py .
"""
_viridis = [[0.267004, 0.004874, 0.329415],
[0.268510, 0.009605, 0.335427],
[0.269944, 0.014625, 0.341379],
[0.271305, 0.019942, 0.347269],
[0.272594, 0.025563, 0.353093],
[0.273809, 0.031497, 0.358853],
[0.274952, 0.037752, 0.364543],
[0.276022, 0.044167, 0.370164],
[0.277018, 0.050344, 0.375715],
[0.277941, 0.056324, 0.381191],
[0.278791, 0.062145, 0.386592],
[0.279566, 0.067836, 0.391917],
[0.280267, 0.073417, 0.397163],
[0.280894, 0.078907, 0.402329],
[0.281446, 0.084320, 0.407414],
[0.281924, 0.089666, 0.412415],
[0.282327, 0.094955, 0.417331],
[0.282656, 0.100196, 0.422160],
[0.282910, 0.105393, 0.426902],
[0.283091, 0.110553, 0.431554],
[0.283197, 0.115680, 0.436115],
[0.283229, 0.120777, 0.440584],
[0.283187, 0.125848, 0.444960],
[0.283072, 0.130895, 0.449241],
[0.282884, 0.135920, 0.453427],
[0.282623, 0.140926, 0.457517],
[0.282290, 0.145912, 0.461510],
[0.281887, 0.150881, 0.465405],
[0.281412, 0.155834, 0.469201],
[0.280868, 0.160771, 0.472899],
[0.280255, 0.165693, 0.476498],
[0.279574, 0.170599, 0.479997],
[0.278826, 0.175490, 0.483397],
[0.278012, 0.180367, 0.486697],
[0.277134, 0.185228, 0.489898],
[0.276194, 0.190074, 0.493001],
[0.275191, 0.194905, 0.496005],
[0.274128, 0.199721, 0.498911],
[0.273006, 0.204520, 0.501721],
[0.271828, 0.209303, 0.504434],
[0.270595, 0.214069, 0.507052],
[0.269308, 0.218818, 0.509577],
[0.267968, 0.223549, 0.512008],
[0.266580, 0.228262, 0.514349],
[0.265145, 0.232956, 0.516599],
[0.263663, 0.237631, 0.518762],
[0.262138, 0.242286, 0.520837],
[0.260571, 0.246922, 0.522828],
[0.258965, 0.251537, 0.524736],
[0.257322, 0.256130, 0.526563],
[0.255645, 0.260703, 0.528312],
[0.253935, 0.265254, 0.529983],
[0.252194, 0.269783, 0.531579],
[0.250425, 0.274290, 0.533103],
[0.248629, 0.278775, 0.534556],
[0.246811, 0.283237, 0.535941],
[0.244972, 0.287675, 0.537260],
[0.243113, 0.292092, 0.538516],
[0.241237, 0.296485, 0.539709],
[0.239346, 0.300855, 0.540844],
[0.237441, 0.305202, 0.541921],
[0.235526, 0.309527, 0.542944],
[0.233603, 0.313828, 0.543914],
[0.231674, 0.318106, 0.544834],
[0.229739, 0.322361, 0.545706],
[0.227802, 0.326594, 0.546532],
[0.225863, 0.330805, 0.547314],
[0.223925, 0.334994, 0.548053],
[0.221989, 0.339161, 0.548752],
[0.220057, 0.343307, 0.549413],
[0.218130, 0.347432, 0.550038],
[0.216210, 0.351535, 0.550627],
[0.214298, 0.355619, 0.551184],
[0.212395, 0.359683, 0.551710],
[0.210503, 0.363727, 0.552206],
[0.208623, 0.367752, 0.552675],
[0.206756, 0.371758, 0.553117],
[0.204903, 0.375746, 0.553533],
[0.203063, 0.379716, 0.553925],
[0.201239, 0.383670, 0.554294],
[0.199430, 0.387607, 0.554642],
[0.197636, 0.391528, 0.554969],
[0.195860, 0.395433, 0.555276],
[0.194100, 0.399323, 0.555565],
[0.192357, 0.403199, 0.555836],
[0.190631, 0.407061, 0.556089],
[0.188923, 0.410910, 0.556326],
[0.187231, 0.414746, 0.556547],
[0.185556, 0.418570, 0.556753],
[0.183898, 0.422383, 0.556944],
[0.182256, 0.426184, 0.557120],
[0.180629, 0.429975, 0.557282],
[0.179019, 0.433756, 0.557430],
[0.177423, 0.437527, 0.557565],
[0.175841, 0.441290, 0.557685],
[0.174274, 0.445044, 0.557792],
[0.172719, 0.448791, 0.557885],
[0.171176, 0.452530, 0.557965],
[0.169646, 0.456262, 0.558030],
[0.168126, 0.459988, 0.558082],
[0.166617, 0.463708, 0.558119],
[0.165117, 0.467423, 0.558141],
[0.163625, 0.471133, 0.558148],
[0.162142, 0.474838, 0.558140],
[0.160665, 0.478540, 0.558115],
[0.159194, 0.482237, 0.558073],
[0.157729, 0.485932, 0.558013],
[0.156270, 0.489624, 0.557936],
[0.154815, 0.493313, 0.557840],
[0.153364, 0.497000, 0.557724],
[0.151918, 0.500685, 0.557587],
[0.150476, 0.504369, 0.557430],
[0.149039, 0.508051, 0.557250],
[0.147607, 0.511733, 0.557049],
[0.146180, 0.515413, 0.556823],
[0.144759, 0.519093, 0.556572],
[0.143343, 0.522773, 0.556295],
[0.141935, 0.526453, 0.555991],
[0.140536, 0.530132, 0.555659],
[0.139147, 0.533812, 0.555298],
[0.137770, 0.537492, 0.554906],
[0.136408, 0.541173, 0.554483],
[0.135066, 0.544853, 0.554029],
[0.133743, 0.548535, 0.553541],
[0.132444, 0.552216, 0.553018],
[0.131172, 0.555899, 0.552459],
[0.129933, 0.559582, 0.551864],
[0.128729, 0.563265, 0.551229],
[0.127568, 0.566949, 0.550556],
[0.126453, 0.570633, 0.549841],
[0.125394, 0.574318, 0.549086],
[0.124395, 0.578002, 0.548287],
[0.123463, 0.581687, 0.547445],
[0.122606, 0.585371, 0.546557],
[0.121831, 0.589055, 0.545623],
[0.121148, 0.592739, 0.544641],
[0.120565, 0.596422, 0.543611],
[0.120092, 0.600104, 0.542530],
[0.119738, 0.603785, 0.541400],
[0.119512, 0.607464, 0.540218],
[0.119423, 0.611141, 0.538982],
[0.119483, 0.614817, 0.537692],
[0.119699, 0.618490, 0.536347],
[0.120081, 0.622161, 0.534946],
[0.120638, 0.625828, 0.533488],
[0.121380, 0.629492, 0.531973],
[0.122312, 0.633153, 0.530398],
[0.123444, 0.636809, 0.528763],
[0.124780, 0.640461, 0.527068],
[0.126326, 0.644107, 0.525311],
[0.128087, 0.647749, 0.523491],
[0.130067, 0.651384, 0.521608],
[0.132268, 0.655014, 0.519661],
[0.134692, 0.658636, 0.517649],
[0.137339, 0.662252, 0.515571],
[0.140210, 0.665859, 0.513427],
[0.143303, 0.669459, 0.511215],
[0.146616, 0.673050, 0.508936],
[0.150148, 0.676631, 0.506589],
[0.153894, 0.680203, 0.504172],
[0.157851, 0.683765, 0.501686],
[0.162016, 0.687316, 0.499129],
[0.166383, 0.690856, 0.496502],
[0.170948, 0.694384, 0.493803],
[0.175707, 0.697900, 0.491033],
[0.180653, 0.701402, 0.488189],
[0.185783, 0.704891, 0.485273],
[0.191090, 0.708366, 0.482284],
[0.196571, 0.711827, 0.479221],
[0.202219, 0.715272, 0.476084],
[0.208030, 0.718701, 0.472873],
[0.214000, 0.722114, 0.469588],
[0.220124, 0.725509, 0.466226],
[0.226397, 0.728888, 0.462789],
[0.232815, 0.732247, 0.459277],
[0.239374, 0.735588, 0.455688],
[0.246070, 0.738910, 0.452024],
[0.252899, 0.742211, 0.448284],
[0.259857, 0.745492, 0.444467],
[0.266941, 0.748751, 0.440573],
[0.274149, 0.751988, 0.436601],
[0.281477, 0.755203, 0.432552],
[0.288921, 0.758394, 0.428426],
[0.296479, 0.761561, 0.424223],
[0.304148, 0.764704, 0.419943],
[0.311925, 0.767822, 0.415586],
[0.319809, 0.770914, 0.411152],
[0.327796, 0.773980, 0.406640],
[0.335885, 0.777018, 0.402049],
[0.344074, 0.780029, 0.397381],
[0.352360, 0.783011, 0.392636],
[0.360741, 0.785964, 0.387814],
[0.369214, 0.788888, 0.382914],
[0.377779, 0.791781, 0.377939],
[0.386433, 0.794644, 0.372886],
[0.395174, 0.797475, 0.367757],
[0.404001, 0.800275, 0.362552],
[0.412913, 0.803041, 0.357269],
[0.421908, 0.805774, 0.351910],
[0.430983, 0.808473, 0.346476],
[0.440137, 0.811138, 0.340967],
[0.449368, 0.813768, 0.335384],
[0.458674, 0.816363, 0.329727],
[0.468053, 0.818921, 0.323998],
[0.477504, 0.821444, 0.318195],
[0.487026, 0.823929, 0.312321],
[0.496615, 0.826376, 0.306377],
[0.506271, 0.828786, 0.300362],
[0.515992, 0.831158, 0.294279],
[0.525776, 0.833491, 0.288127],
[0.535621, 0.835785, 0.281908],
[0.545524, 0.838039, 0.275626],
[0.555484, 0.840254, 0.269281],
[0.565498, 0.842430, 0.262877],
[0.575563, 0.844566, 0.256415],
[0.585678, 0.846661, 0.249897],
[0.595839, 0.848717, 0.243329],
[0.606045, 0.850733, 0.236712],
[0.616293, 0.852709, 0.230052],
[0.626579, 0.854645, 0.223353],
[0.636902, 0.856542, 0.216620],
[0.647257, 0.858400, 0.209861],
[0.657642, 0.860219, 0.203082],
[0.668054, 0.861999, 0.196293],
[0.678489, 0.863742, 0.189503],
[0.688944, 0.865448, 0.182725],
[0.699415, 0.867117, 0.175971],
[0.709898, 0.868751, 0.169257],
[0.720391, 0.870350, 0.162603],
[0.730889, 0.871916, 0.156029],
[0.741388, 0.873449, 0.149561],
[0.751884, 0.874951, 0.143228],
[0.762373, 0.876424, 0.137064],
[0.772852, 0.877868, 0.131109],
[0.783315, 0.879285, 0.125405],
[0.793760, 0.880678, 0.120005],
[0.804182, 0.882046, 0.114965],
[0.814576, 0.883393, 0.110347],
[0.824940, 0.884720, 0.106217],
[0.835270, 0.886029, 0.102646],
[0.845561, 0.887322, 0.099702],
[0.855810, 0.888601, 0.097452],
[0.866013, 0.889868, 0.095953],
[0.876168, 0.891125, 0.095250],
[0.886271, 0.892374, 0.095374],
[0.896320, 0.893616, 0.096335],
[0.906311, 0.894855, 0.098125],
[0.916242, 0.896091, 0.100717],
[0.926106, 0.897330, 0.104071],
[0.935904, 0.898570, 0.108131],
[0.945636, 0.899815, 0.112838],
[0.955300, 0.901065, 0.118128],
[0.964894, 0.902323, 0.123941],
[0.974417, 0.903590, 0.130215],
[0.983868, 0.904867, 0.136897],
[0.993248, 0.906157, 0.143936]]


text:
module AtomicDFT
using DocStringExtensions
using ArgCheck
using Unitful, UnitfulHartree
using LibXC: DFTUnits, Units, XCFunctional
using AxisArrays
using Iterators: chain
using SimpleTraits

export zeros_axisarray, density_array, is_spin_polarized, radial_hartree

const HartreeUnits = Units;

export OrbitalIndex, @nl_str

macro lintpragma(s) end

include("OrbitalIndexing.jl")
include("Traits.jl")
include("ArrayCreation.jl")
# include("Potentials.jl")
include("Radial.jl")
end # module


text:
mean(d::LogitNormal; kwargs...) = estimateMean(d,kwargs...)

function mode(d::LogitNormal)
    (μ, σ) = params(d)
    # if mu<0 then maximum is left of median, if mu>0 right of median
    # if mu=0 the maximum is either at mu or there are two maxima of the same height
    # here, report the left mode
    interval = μ <= 0.0 ? (0.0,logistic(μ)) : (logistic(μ),1.0)
    resOpt = optimize(x -> -pdf(d, x), interval[1], interval[2])
    resOpt.minimizer
end

function var(d::LogitNormal; mean=missing, kwargs...)
    # estimateVariance does not pass kwargs to mean, need to do it here
    m = ismissing(mean) ? estimateMean(d,kwargs...) : mean
    estimateVariance(d; kwargs..., mean=m)
end

# function fit(::Type{LogitNormal}, m::AbstractMoments)
#     error("Fitting LogitNormal to moments not implemented.")
# end

function fit(::Type{LogitNormal}, lower::QuantilePoint, upper::QuantilePoint)
    lower_logit = QuantilePoint(lower, q = logit(lower.q))
    upper_logit = QuantilePoint(upper, q = logit(upper.q))
    DN = fit(Normal, lower_logit, upper_logit)
    LogitNormal(params(DN)...)
end

function fit_mode_quantile(::Type{LogitNormal}, mode::Real, qp::QuantilePoint)
    matchModeUpper(mode, qp.q, Val(40); perc = qp.p)
end

function matchModeUpper(
    mode::T, upper::T
    , ::Val{nTry}
    ; perc::Real = 0.99
) where {nTry, T<:Real}
    mode == 0.5 && return matchMedianUpper(LogitNormal, 0.5, upper; perc=perc)
    # for given mu we can compute sigma by mode and upper quantile
    # hence univariate search for mu
    # we now that mu is in (\code{logit(mode)},0) for \code{mode < 0.5} and in
    # \code{(0,logit(mode))} for \code{mode > 0.5} for unimodal distribution
    # there might be a maximum in the middle and optimized misses the low part
    # hence, first get near the global minimum by a evaluating the cost at a
    # grid that is spaced narrower at the edge
    #
    logitMode = logit(mode)
    logitUpper = logit(upper)
    upperMu = abs(logitMode) - eps()
    muTry = SVector{nTry}(
        sign(logitMode) .* log.(range(1,stop=exp(upperMu),length=nTry)))
    oF(mu) = ofLogitNormalModeUpper(mu, mode, logitMode, logitUpper, perc)
    ofMuTry = oF.(muTry)
    iMin = argmin(ofMuTry)
    # on postive side muTry are increasing, on negative side muTry decreasing
    # neet to have the lower value at the beginning of the interval
    interval = (logitMode >= 0) ?
        (muTry[max(1,iMin-1)], muTry[min(nTry,iMin+1)]) :
        (muTry[max(1,iMin+1)], muTry[min(nTry,max(1,iMin-1))])
    resOpt = optimize(oF, interval...)
    resOpt.converged || error("could not find minimum")
    μ = resOpt.minimizer
    σ = sqrt((logitMode - μ)/(2*mode - 1))
    LogitNormal{T}(μ,σ)
end

function matchModeUpper(
    mode::S, upper::T, ::Val{nTry}; kwargs...
) where {nTry, T<:Real, S<:Real}
    p = promote(mode, upper)
    matchModeUpper(p[1], p[2], Val(nTry); kwargs...)
end

"objective function used by `matchModeUpper(LogitNormal,...)`"
function ofLogitNormalModeUpper(mu, mode, logitMode, logitUpper, perc)
  # given mu and mode, we can calculate sigma,
  # predict the percentile of logitUpper
  # and return the squared difference as cost to be minimized
  sigma2 = (logitMode - mu)/(2.0*mode - 1.0)
  sigma2 > 0.0 || return prevfloat(Inf)
  predp = normcdf(mu, sqrt(sigma2), logitUpper)
  diff = predp - perc
  diff*diff
end

"""
    fit_mode_flat(::Type{LogitNormal}, mode::T; peakedness = 1.0)

Find the maximum-spread logitnormal distribution that has a single mode at given location.

More peaked distributions with given single mode can be optained by increasing
argument peakedness. They will have a spread by originally inferred σ² devidied
by peakedness.

# Examples
```jldoctest; output=false
m = 0.6
d = fit_mode_flat(LogitNormal, m; peakedness = 1.5)
mode(d) ≈ m
# output
true
```
"""
function fit_mode_flat(::Type{LogitNormal}, mode::T, ::Val{nTry} = Val(40); peakedness = 1.0) where {T<:Real,nTry}
    mode == 0.5 && return(LogitNormal{T}(0.0, sqrt(2)/peakedness))
    is_right = mode > 0.5
    mode_r = is_right ? mode : 1.0 - mode
    res_opt = optimize(x -> of_mode_flat(x, mode_r, logit(mode_r)), 0.0, 0.5)
    res_opt.converged || error("could not find minimum")
    xt = res_opt.minimizer
    σ2 = (1/xt + 1/(1-xt))/2/peakedness^2
    μr = logit(mode_r) - σ2*(2.0*mode_r - 1.0)
    μ = is_right ? μr : -μr
    LogitNormal{T}(μ,sqrt(σ2))
end

function of_mode_flat(x, m, logitm = logit(m))
  # lhs = 1/x + 1/(1-x)
  # rhs = (logitm - logit(x))/(m-x)
  # lhs = (m-x)/x + (m-x)/(1-x)
  # rhs = (logitm - logit(x))
  lhs = m/x + (m-1)/(1-x)
  rhs = logitm - logit(x)
  d = lhs - rhs
  d*d
end

"""
    shifloNormal(lower,upper)

Get a Shifted Flat LogitNormal distribution that is most spread
with an extent between lower and upper.
This is a more smooth alternative to the bounded uniform distribution.
"""
function shifloNormal(lower,upper)
    lower,upper = promote(lower,upper/1)
    dln = LogitNormal{typeof(lower)}(0.0, sqrt(2))
    LocationScale(lower, (upper -lower), dln)
end


text:
module HonlyHelperFlagSoundSource

using ..Ahorn, Maple

@mapdef Entity "HonlyHelper/FlagSoundSource" FlagSoundSource(x::Integer, y::Integer, sound::String = "sound", flag::String = "flag", AllowFade::Bool=false)

const placements = Ahorn.PlacementDict(
    "Flag Sound Source (HonlyHelper)" => Ahorn.EntityPlacement(
        FlagSoundSource
    )
)

function Ahorn.selection(entity::FlagSoundSource)
    x, y = Ahorn.position(entity::FlagSoundSource)
    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

Ahorn.editingOptions(entity::FlagSoundSource) = Dict{String, Any}(
    "sound" => EnvironmentSounds.sounds
)

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::FlagSoundSource, room::Maple.Room) = Ahorn.drawImage(ctx, Ahorn.Assets.speaker, -12, -12)

end

max_stars_repo_name: HollyMagala/HonlyHelper

max_issues_repo_name: HollyMagala/HonlyHelper

max_forks_repo_name: HollyMagala/HonlyHelper


text:
using SententialDecisionDiagrams
const SDD = SententialDecisionDiagrams

# set up vtree and manager
var_count = 4
vtree_type = "right"
vtree = SDD.vtree(var_count, vtree_type)
manager = SDD.sdd_manager(vtree)

x = [SDD.literal(i,manager) for i in 1:5]

# construct the term X_1 ^ X_2 ^ X_3 ^ X_4
α =  x[1] &  x[2] & x[3] & x[4]

# construct the term ~X_1 ^ X_2 ^ X_3 ^ X_4
β  = ~x[1] &  x[2] & x[3] & x[4]

# construct the term ~X_1 ^ ~X_2 ^ X_3 ^ X_4
γ = ~x[1] & ~x[2] & x[3] & x[4]

println("before referencing:")
println("live sdd size = $(SDD.live_size(manager))")
println("dead sdd size = $(SDD.dead_size(manager))")

# ref SDDs so that they are not garbage collected
SDD.ref(α)
SDD.ref(β)
SDD.ref(γ)
println("after referencing:")
println("live sdd size = $(SDD.live_size(manager))")
println("dead sdd size = $(SDD.dead_size(manager))")

# garbage collect
SDD.garbage_collect(manager)
println("after garbage collection:");
println("live sdd size = $(SDD.live_size(manager))")
println("dead sdd size = $(SDD.dead_size(manager))")

SDD.deref(α)
SDD.deref(β)
SDD.deref(γ)

println("saving vtree & shared sdd ...")

SDD.dot("$(@__DIR__)/output/shared-vtree.dot", vtree)
SDD.shared_save_as_dot("$(@__DIR__)/output/shared.dot", manager)


text:
@doc raw"""
    typical(ohlc::TimeArray; h = :High, l = :Low, c = :Close)

Typical Price

```math
\text{Typical Price} = \frac{P^\text{High} + P^\text{Low} + P^\text{Close}}{3}
```
"""
function typical(ohlc::TimeArray; h = :High, l = :Low, c = :Close)
  val = (ohlc[h] .+ ohlc[l] .+ ohlc[c]) ./ 3
  rename(val, [:typical])
end

function mean_abs_dev(a::AbstractArray{T,1}, scale::Bool=false) where T
  c = scale ? (1 / (-sqrt(2) * erfcinv(3 * 1 / 2))) : 1

  res = ones(length(a))
  for i in 1:length(a)
    @inbounds res[i] = abs(a[i] - nanmean(a))
  end

  nanmean(res) * c
end

"""
Parameters:
    * orig: the original colnames
    * suffix: the suffix applied to `orig`

```jldoctest
julia> MarketTechnicals.gen_colnames([:Open, :Close], [:macd, :dif, :sig])
6-element Array{Symbol,1}:
 :Open_macd
 :Close_macd
 :Open_dif
 :Close_dif
 :Open_sig
 :Close_sig
```
"""
gen_colnames(orig::Vector{Symbol}, suffix::Vector{Symbol}) =
  vec([Symbol(o, "_", s) for o ∈ orig, s ∈ suffix])

relu(x) = max(x, 0)

relu(ta::TimeArray) = TimeArray(timestamp(ta), relu.(values(ta)), colnames(ta), meta(ta))

"""
The smooth method used by ADX
"""
function wilder_smooth(ta::TimeArray, n::Integer; padding::Bool = false, dx::Bool = false)
  val = similar(Array{Float64}, axes(values(ta)))

  first_cal = (dx ? nanmean : nansum)

  ncols = size(val, 2)
  _vals_ta = values(ta)

  for i ∈ 1:size(val, 1)
    if i < n
      val[i, :] .= NaN
    elseif i == n
      _fv = first_cal(values(ta[1:n]))
      val[i, :] = isa(_fv, Number) ? reshape([_fv], (1,ncols)) : _fv
    elseif dx
      _nv = (val[i-1, :] .* (n - 1) .+ _vals_ta[i, :]) ./ n
      _nv[isnan.(_nv)] .= val[i-1, :][isnan.(_nv)]

      val[i, :] = reshape(_nv, (1, ncols))
    else
      _nv = (val[i-1, :] .* (n - 1) ./ n) .+ _vals_ta[i, :]
      _nv[isnan.(_nv)] .= val[i-1, :][isnan.(_nv)]

      val[i, :] = reshape(_nv, (1, ncols))
    end
  end

  ret = TimeArray(timestamp(ta), val, colnames(ta), meta(ta))

  if padding
    ret
  else
    dropnan(ret)
  end
end

safediv(x, y) = ifelse(iszero(x) && iszero(y), x, x / y)

function lagfill(ta::TimeArray, r1::Integer, fill::AbstractFloat)
  _lta = lag(ta, r1, padding = true)
  _lta_values = values(_lta)
  _lta_values[1:r1] .= fill

  TimeArray(timestamp(ta), _lta_values, colnames(ta))
end

_nanmean(x) = mean(filter(!isnan, x))
nanmean(x; dims = 1) = ndims(x) > 1 ? mapslices(_nanmean, x, dims = dims) : _nanmean(x)

_nansum(x) = sum(filter(!isnan,x))
nansum(x; dims = 1) = ndims(x) > 1 ? mapslices(_nansum, x, dims = dims) : _nansum(x)

_nanstd(x) = std(filter(!isnan,x))
nanstd(x; dims = 1) = ndims(x) > 1 ? mapslices(_nanstd, x, dims = dims) : _nanstd(x)

function nancumsum(x)
	x[isnan.(x)] .= 0
	cumsum(x)
end

function nanargmax(x)
	x[isnan.(x)] .= -Inf
	argmax(x)
end

function nanargmin(x)
	x[isnan.(x)] .= Inf
	argmin(x)
end

function nanmax(x)
	x[isnan.(x)] .= -Inf
	maximum(x)
end

function nanmin(x)
	x[isnan.(x)] .= Inf
	minimum(x)
end


text:
function ConstantProcessNoise(arg0::RealMatrix)
    return ConstantProcessNoise((RealMatrix,), arg0)
end

function ConstantProcessNoise(arg0::RealMatrix, arg1::RealMatrix)
    return ConstantProcessNoise((RealMatrix, RealMatrix), arg0, arg1)
end

function get_process_noise_matrix(obj::ConstantProcessNoise, arg0::SpacecraftState, arg1::SpacecraftState)
    return jcall(obj, "getProcessNoiseMatrix", RealMatrix, (SpacecraftState, SpacecraftState), arg0, arg1)
end


text:
# Some julia inline functions to replace C Macros

"Mark TeX code group"
macro group(name, code)
    code |> esc
end

"x--"
macro decr!(x)
    :(x -= 1) |> esc
end

"x = -x"
macro negate!(x)
    :(x = -x) |> esc
end


text:
using Test
using QCMaterial

@testset "util.topylist" begin
    org_array = [1, 2, 3.0]
    pylist = topylist(org_array)
    @test all(org_array .== pylist)
end

@testset "util.doublefunc" begin
    org_array = [1, 2, 3]
    @test all(doublefunc(org_array) .== 2 .* org_array)
end

@testset "util.numerical_grad" begin
    # Some test for numerical_grad
    f(x) = x[1] + 2*x[2]
    deriv = numerical_grad(f, zeros(2))
    @test deriv ≈ [1.0, 2.0]
end

@testset "util.generate_numerical_grad" begin
    # Some test for numerical_grad
    N = 20
    f(x) = sum(collect(1:N) .* x)
    x = 1. .* collect(1:N)
    deriv = numerical_grad(f, x)
    grad = QCMaterial.generate_numerical_grad(f)
    @test isapprox(grad(x), collect(1:N), rtol=1e-5)
end

@testset "util.fit_svd" begin
    #test for truncated svd
    #y = Ax
    #A = U∑V^{dagger}
    #x = V∑^{-1}(U^{dagger}y)
    N = 10
    A = randn(Float64, (N, N))
    y = ones(Float64, N)
    x = QCMaterial.fit_svd(y, A, 1e-10)
    @test y ≈ A * x
end

@testset "util.tikhonov" begin
    N = 10
    A = rand(Float64, (N,N))
    y = ones(Float64, N)
    x = QCMaterial.tikhonov(y, A, 1e-15)
    @test y ≈ A * x
end

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
abstract type AbstractCoriolis <: AbstractPhysicsComponent end

Base.@kwdef struct DeepShellCoriolis{FT} <: AbstractCoriolis
    Ω::FT # s⁻¹
end

Base.@kwdef struct ThinShellCoriolis{FT} <: AbstractCoriolis
    Ω::FT # s⁻¹
end

Base.@kwdef struct BetaPlaneCoriolis{FT} <: AbstractCoriolis
    f₀::FT # s⁻¹
    β::FT # s⁻¹m⁻¹
end

@inline calc_force!(source, ::Nothing, state, _...) = nothing

@inline function calc_force!(source, coriolis::DeepShellCoriolis, state, _...)
    Ω⃗ = @SVector [-0, -0, coriolis.Ω]
    ρu⃗ = state.ρu

    source.ρu -= 2Ω⃗ × ρu⃗

    return nothing
end

@inline function calc_force!(source, coriolis::ThinShellCoriolis, state, aux, orientation, _...)
    Ω⃗ = @SVector [-0, -0, coriolis.Ω]
    ρu⃗ = state.ρu

    k̂ = vertical_unit_vector(orientation, aux)

    source.ρu -= (2Ω⃗ ⋅ k̂) * (k̂ × ρu⃗)

    return nothing
end

@inline function calc_force!(source, coriolis::BetaPlaneCoriolis, state, aux, orientation, _...)
    f₀ = coriolis.f₀
    β = coriolis.β
    ρu⃗ = state.ρu
    y = aux.y

    f = f₀ + β * y
    k̂ = vertical_unit_vector(orientation, aux)

    source.ρu -= f * (k̂ × ρu⃗)

    return nothing
end

max_stars_repo_name: CliMA/ClimaCoupler