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

# editing files

"""
    editor()

Determines the editor to use when running functions like `edit`. Returns an Array compatible
for use within backticks. You can change the editor by setting `JULIA_EDITOR`, `VISUAL` or
`EDITOR` as an environment variable.
"""
function editor()
    if Sys.iswindows() || Sys.isapple()
        default_editor = "open"
    elseif isfile("/etc/alternatives/editor")
        default_editor = "/etc/alternatives/editor"
    else
        default_editor = "emacs"
    end
    # Note: the editor path can include spaces (if escaped) and flags.
    args = shell_split(get(ENV,"JULIA_EDITOR", get(ENV,"VISUAL", get(ENV,"EDITOR", default_editor))))
    isempty(args) && error("editor is empty")
    return args
end

"""
    edit(path::AbstractString, line::Integer=0)

Edit a file or directory optionally providing a line number to edit the file at.
Returns to the `julia` prompt when you quit the editor. The editor can be changed
by setting `JULIA_EDITOR`, `VISUAL` or `EDITOR` as an environment variable.
"""
function edit(path::AbstractString, line::Integer=0)
    command = editor()
    name = basename(first(command))
    issrc = length(path)>2 && path[end-2:end] == ".jl"
    if issrc
        f = find_source_file(path)
        f !== nothing && (path = f)
    end
    background = true
    line_unsupported = false
    if startswith(name, "vim.") || name == "vi" || name == "vim" || name == "nvim" ||
            name == "mvim" || name == "nano" ||
            name == "emacs" && contains(in, command, ["-nw", "--no-window-system" ]) ||
            name == "emacsclient" && contains(in, command, ["-nw", "-t", "-tty"])
        cmd = line != 0 ? `$command +$line $path` : `$command $path`
        background = false
    elseif startswith(name, "emacs") || name == "gedit" || startswith(name, "gvim")
        cmd = line != 0 ? `$command +$line $path` : `$command $path`
    elseif name == "textmate" || name == "mate" || name == "kate"
        cmd = line != 0 ? `$command $path -l $line` : `$command $path`
    elseif startswith(name, "subl") || startswith(name, "atom")
        cmd = line != 0 ? `$command $path:$line` : `$command $path`
    elseif name == "code" || (Sys.iswindows() && uppercase(name) == "CODE.EXE")
        cmd = line != 0 ? `$command -g $path:$line` : `$command -g $path`
    elseif startswith(name, "notepad++")
        cmd = line != 0 ? `$command $path -n$line` : `$command $path`
    elseif Sys.isapple() && name == "open"
        cmd = `open -t $path`
        line_unsupported = true
    else
        cmd = `$command $path`
        background = false
        line_unsupported = true
    end

    if Sys.iswindows() && name == "open"
        @static Sys.iswindows() && # don't emit this ccall on other platforms
            systemerror(:edit, ccall((:ShellExecuteW, "shell32"), stdcall, Int,
                                     (Ptr{Void}, Cwstring, Cwstring, Ptr{Void}, Ptr{Void}, Cint),
                                     C_NULL, "open", path, C_NULL, C_NULL, 10) ≤ 32)
    elseif background
        spawn(pipeline(cmd, stderr=STDERR))
    else
        run(cmd)
    end
    line != 0 && line_unsupported && println("Unknown editor: no line number information passed.\nThe method is defined at line $line.")

    nothing
end

"""
    edit(function, [types])

Edit the definition of a function, optionally specifying a tuple of types to
indicate which method to edit. The editor can be changed by setting `JULIA_EDITOR`,
`VISUAL` or `EDITOR` as an environment variable.
"""
edit(f)                   = edit(functionloc(f)...)
edit(f, @nospecialize t)  = edit(functionloc(f,t)...)
edit(file, line::Integer) = error("could not find source file for function")

# terminal pager

if Sys.iswindows()
    function less(file::AbstractString, line::Integer)
        pager = shell_split(get(ENV, "PAGER", "more"))
        g = pager[1] == "more" ? "" : "g"
        run(Cmd(`$pager +$(line)$(g) \"$file\"`, windows_verbatim = true))
    end
else
    function less(file::AbstractString, line::Integer)
        pager = shell_split(get(ENV, "PAGER", "less"))
        run(`$pager +$(line)g $file`)
    end
end

"""
    less(file::AbstractString, [line::Integer])

Show a file using the default pager, optionally providing a starting line number. Returns to
the `julia` prompt when you quit the pager.
"""
less(file::AbstractString) = less(file, 1)

"""
    less(function, [types])

Show the definition of a function using the default pager, optionally specifying a tuple of
types to indicate which method to see.
"""
less(f)                   = less(functionloc(f)...)
less(f, @nospecialize t)  = less(functionloc(f,t)...)
less(file, line::Integer) = error("could not find source file for function")

# clipboard copy and paste

if Sys.isapple()
    function clipboard(x)
        open(pipeline(`pbcopy`, stderr=STDERR), "w") do io
            print(io, x)
        end
    end
    clipboard() = read(`pbpaste`, String)

elseif Sys.islinux()
    _clipboardcmd = nothing
    function clipboardcmd()
        global _clipboardcmd
        _clipboardcmd !== nothing && return _clipboardcmd
        for cmd in (:xclip, :xsel)
            success(pipeline(`which $cmd`, DevNull)) && return _clipboardcmd = cmd
        end
        error("no clipboard command found, please install xsel or xclip")
    end
    function clipboard(x)
        c = clipboardcmd()
        cmd = c == :xsel  ? `xsel --nodetach --input --clipboard` :
              c == :xclip ? `xclip -silent -in -selection clipboard` :
            error("unexpected clipboard command: $c")
        open(pipeline(cmd, stderr=STDERR), "w") do io
            print(io, x)
        end
    end
    function clipboard()
        c = clipboardcmd()
        cmd = c == :xsel  ? `xsel --nodetach --output --clipboard` :
              c == :xclip ? `xclip -quiet -out -selection clipboard` :
            error("unexpected clipboard command: $c")
        read(pipeline(cmd, stderr=STDERR), String)
    end

elseif Sys.iswindows()
    # TODO: these functions leak memory and memory locks if they throw an error
    function clipboard(x::AbstractString)
        if containsnul(x)
            throw(ArgumentError("Windows clipboard strings cannot contain NUL character"))
        end
        systemerror(:OpenClipboard, 0==ccall((:OpenClipboard, "user32"), stdcall, Cint, (Ptr{Void},), C_NULL))
        systemerror(:EmptyClipboard, 0==ccall((:EmptyClipboard, "user32"), stdcall, Cint, ()))
        x_u16 = cwstring(x)
        # copy data to locked, allocated space
        p = ccall((:GlobalAlloc, "kernel32"), stdcall, Ptr{UInt16}, (UInt16, Int32), 2, sizeof(x_u16))
        systemerror(:GlobalAlloc, p==C_NULL)
        plock = ccall((:GlobalLock, "kernel32"), stdcall, Ptr{UInt16}, (Ptr{UInt16},), p)
        systemerror(:GlobalLock, plock==C_NULL)
        ccall(:memcpy, Ptr{UInt16}, (Ptr{UInt16},Ptr{UInt16},Int), plock, x_u16, sizeof(x_u16))
        systemerror(:GlobalUnlock, 0==ccall((:GlobalUnlock, "kernel32"), stdcall, Cint, (Ptr{Void},), plock))
        pdata = ccall((:SetClipboardData, "user32"), stdcall, Ptr{UInt16}, (UInt32, Ptr{UInt16}), 13, p)
        systemerror(:SetClipboardData, pdata!=p)
        ccall((:CloseClipboard, "user32"), stdcall, Void, ())
    end
    clipboard(x) = clipboard(sprint(print, x)::String)
    function clipboard()
        systemerror(:OpenClipboard, 0==ccall((:OpenClipboard, "user32"), stdcall, Cint, (Ptr{Void},), C_NULL))
        pdata = ccall((:GetClipboardData, "user32"), stdcall, Ptr{UInt16}, (UInt32,), 13)
        systemerror(:SetClipboardData, pdata==C_NULL)
        systemerror(:CloseClipboard, 0==ccall((:CloseClipboard, "user32"), stdcall, Cint, ()))
        plock = ccall((:GlobalLock, "kernel32"), stdcall, Ptr{UInt16}, (Ptr{UInt16},), pdata)
        systemerror(:GlobalLock, plock==C_NULL)
        # find NUL terminator (0x0000 16-bit code unit)
        len = 0
        while unsafe_load(plock, len+1) != 0; len += 1; end
        # get Vector{UInt16}, transcode data to UTF-8, make a String of it
        s = transcode(String, unsafe_wrap(Array, plock, len))
        systemerror(:GlobalUnlock, 0==ccall((:GlobalUnlock, "kernel32"), stdcall, Cint, (Ptr{UInt16},), plock))
        return s
    end

else
    clipboard(x="") = error("`clipboard` function not implemented for $(Sys.KERNEL)")
end

"""
    clipboard(x)

Send a printed form of `x` to the operating system clipboard ("copy").
"""
clipboard(x)

"""
    clipboard() -> AbstractString

Return a string with the contents of the operating system clipboard ("paste").
"""
clipboard()

# system information

# used by sysinfo.jl
function _show_cpuinfo(io::IO, info::Sys.CPUinfo, header::Bool=true, prefix::AbstractString="    ")
    tck = Sys.SC_CLK_TCK
    if header
        println(io, info.model, ": ")
        print(io, " "^length(prefix))
        if tck > 0
            @printf(io, "    %5s    %9s    %9s    %9s    %9s    %9s\n",
                    "speed", "user", "nice", "sys", "idle", "irq")
        else
            @printf(io, "    %5s    %9s  %9s  %9s  %9s  %9s ticks\n",
                    "speed", "user", "nice", "sys", "idle", "irq")
        end
    end
    print(io, prefix)
    if tck > 0
        @printf(io, "%5d MHz  %9d s  %9d s  %9d s  %9d s  %9d s",
                info.speed, info.cpu_times!user / tck, info.cpu_times!nice / tck,
                info.cpu_times!sys / tck, info.cpu_times!idle / tck, info.cpu_times!irq / tck)
    else
        @printf(io, "%5d MHz  %9d  %9d  %9d  %9d  %9d ticks",
                info.speed, info.cpu_times!user, info.cpu_times!nice,
                info.cpu_times!sys, info.cpu_times!idle, info.cpu_times!irq)
    end
end

"""
    versioninfo(io::IO=STDOUT; verbose::Bool=false, packages::Bool=false)

Print information about the version of Julia in use. The output is
controlled with boolean keyword arguments:

- `packages`: print information about installed packages
- `verbose`: print all additional information
"""
function versioninfo(io::IO=STDOUT; verbose::Bool=false, packages::Bool=false)
    println(io, "Julia Version $VERSION")
    if !isempty(GIT_VERSION_INFO.commit_short)
        println(io, "Commit $(GIT_VERSION_INFO.commit_short) ($(GIT_VERSION_INFO.date_string))")
    end
    if ccall(:jl_is_debugbuild, Cint, ())!=0
        println(io, "DEBUG build")
    end
    println(io, "Platform Info:")
    println(io, "  OS: ", Sys.iswindows() ? "Windows" : Sys.isapple() ?
        "macOS" : Sys.KERNEL, " (", Sys.MACHINE, ")")

    if verbose
        lsb = ""
        if Sys.islinux()
            try lsb = readchomp(pipeline(`lsb_release -ds`, stderr=DevNull)) end
        end
        if Sys.iswindows()
            try lsb = strip(read(`$(ENV["COMSPEC"]) /c ver`, String)) end
        end
        if !isempty(lsb)
            println(io, "      ", lsb)
        end
        if Sys.isunix()
            println(io, "  uname: ", readchomp(`uname -mprsv`))
        end
    end

    if verbose
        cpuio = IOBuffer() # print cpu_summary with correct alignment
        Sys.cpu_summary(cpuio)
        for (i, line) in enumerate(split(String(take!(cpuio)), "\n"))
            prefix = i == 1 ? "  CPU: " : "       "
            println(io, prefix, line)
        end
    else
        cpu = Sys.cpu_info()
        println(io, "  CPU: ", cpu[1].model)
    end

    if verbose
        println(io, "  Memory: $(Sys.total_memory()/2^30) GB ($(Sys.free_memory()/2^20) MB free)")
        try println(io, "  Uptime: $(Sys.uptime()) sec") end
        print(io, "  Load Avg: ")
        print_matrix(io, Sys.loadavg()')
        println(io)
    end
    println(io, "  WORD_SIZE: ", Sys.WORD_SIZE)
    if Base.libblas_name == "libopenblas" || BLAS.vendor() == :openblas || BLAS.vendor() == :openblas64
        openblas_config = BLAS.openblas_get_config()
        println(io, "  BLAS: libopenblas (", openblas_config, ")")
    else
        println(io, "  BLAS: ",libblas_name)
    end
    println(io, "  LAPACK: ",liblapack_name)
    println(io, "  LIBM: ",libm_name)
    println(io, "  LLVM: libLLVM-",libllvm_version," (", Sys.JIT, ", ", Sys.cpu_name, ")")

    println(io, "Environment:")
    for (k,v) in ENV
        if ismatch(r"JULIA", String(k))
            println(io, "  $(k) = $(v)")
        end
    end
    if verbose
        for (k,v) in ENV
            if ismatch(r"PATH|FLAG|^TERM$|HOME", String(k))
                println(io, "  $(k) = $(v)")
            end
        end
    end
    if packages || verbose
        println(io, "Packages:")
        println(io, "  Package Directory: ", Pkg.dir())
        print(io, "  Package Status:")
        if isdir(Pkg.dir())
            println(io, "")
            Pkg.status(io)
        else
            println(io, " no packages installed")
        end
    end
end

# displaying type-ambiguity warnings
"""
    code_warntype([io::IO], f, types)

Prints lowered and type-inferred ASTs for the methods matching the given generic function
and type signature to `io` which defaults to `STDOUT`. The ASTs are annotated in such a way
as to cause "non-leaf" types to be emphasized (if color is available, displayed in red).
This serves as a warning of potential type instability. Not all non-leaf types are particularly
problematic for performance, so the results need to be used judiciously.
See [`@code_warntype`](@ref man-code-warntype) for more information.
"""
function code_warntype(io::IO, f, @nospecialize(t))
    function slots_used(ci, slotnames)
        used = falses(length(slotnames))
        scan_exprs!(used, ci.code)
        return used
    end

    function scan_exprs!(used, exprs)
        for ex in exprs
            if isa(ex, Slot)
                used[ex.id] = true
            elseif isa(ex, Expr)
                scan_exprs!(used, ex.args)
            end
        end
    end

    emph_io = IOContext(io, :TYPEEMPHASIZE => true)
    for (src, rettype) in code_typed(f, t)
        println(emph_io, "Variables:")
        slotnames = sourceinfo_slotnames(src)
        used_slotids = slots_used(src, slotnames)
        for i = 1:length(slotnames)
            print(emph_io, "  ", slotnames[i])
            if used_slotids[i]
                if isa(src.slottypes, Array)
                    show_expr_type(emph_io, src.slottypes[i], true)
                end
            else
                print(emph_io, " <optimized out>")
            end
            print(emph_io, '\n')
        end
        print(emph_io, "\nBody:\n  ")
        body = Expr(:body)
        body.args = src.code
        body.typ = rettype
        # Fix slot names and types in function body
        show_unquoted(IOContext(emph_io, :SOURCEINFO => src, :SOURCE_SLOTNAMES => slotnames),
                      body, 2)
        print(emph_io, '\n')
    end
    nothing
end
code_warntype(f, @nospecialize(t)) = code_warntype(STDOUT, f, t)

typesof(args...) = Tuple{Any[ Core.Typeof(a) for a in args ]...}

function gen_call_with_extracted_types(__module__, fcn, ex0)
    if isa(ex0, Expr)
        if any(a->(Meta.isexpr(a, :kw) || Meta.isexpr(a, :parameters)), ex0.args)
            # remove keyword args, but call the kwfunc
            args = filter(a->!(Meta.isexpr(a, :kw) || Meta.isexpr(a, :parameters)), ex0.args)
            return quote
                local arg1 = $(esc(args[1]))
                $(fcn)(Core.kwfunc(arg1),
                       Tuple{Vector{Any}, Core.Typeof(arg1),
                             $(typesof)($(map(esc, args[2:end])...)).parameters...})
            end
        elseif ex0.head == :call
            return Expr(:call, fcn, esc(ex0.args[1]),
                        Expr(:call, typesof, map(esc, ex0.args[2:end])...))
        end
    end
    if isa(ex0, Expr) && ex0.head == :macrocall # Make @edit @time 1+2 edit the macro by using the types of the *expressions*
        return Expr(:call, fcn, esc(ex0.args[1]), Tuple{#=__source__=#LineNumberNode, #=__module__=#Module, Any[ Core.Typeof(a) for a in ex0.args[3:end] ]...})
    end
    ex = expand(__module__, ex0)
    exret = Expr(:none)
    if !isa(ex, Expr)
        exret = Expr(:call, :error, "expression is not a function call or symbol")
    elseif ex.head == :call
        if any(e->(isa(e, Expr) && e.head==:(...)), ex0.args) &&
            (ex.args[1] === GlobalRef(Core,:_apply) ||
             ex.args[1] === GlobalRef(Base,:_apply))
            # check for splatting
            exret = Expr(:call, ex.args[1], fcn,
                        Expr(:tuple, esc(ex.args[2]),
                            Expr(:call, typesof, map(esc, ex.args[3:end])...)))
        else
            exret = Expr(:call, fcn, esc(ex.args[1]),
                         Expr(:call, typesof, map(esc, ex.args[2:end])...))
        end
    elseif ex.head == :body
        a1 = ex.args[1]
        if isa(a1, Expr) && a1.head == :call
            a11 = a1.args[1]
            if a11 == :setindex!
                exret = Expr(:call, fcn, a11,
                             Expr(:call, typesof, map(esc, a1.args[2:end])...))
            end
        end
    end
    if ex.head == :thunk || exret.head == :none
        exret = Expr(:call, :error, "expression is not a function call, "
                                  * "or is too complex for @$fcn to analyze; "
                                  * "break it down to simpler parts if possible")
    end
    return exret
end

for fname in [:which, :less, :edit, :functionloc, :code_warntype,
              :code_llvm, :code_llvm_raw, :code_native]
    @eval begin
        macro ($fname)(ex0)
            gen_call_with_extracted_types(__module__, $(Expr(:quote, fname)), ex0)
        end
    end
end
macro which(ex0::Symbol)
    ex0 = QuoteNode(ex0)
    return :(which_module($__module__, $ex0))
end

for fname in [:code_typed, :code_lowered]
    @eval begin
        macro ($fname)(ex0)
            thecall = gen_call_with_extracted_types(__module__, $(Expr(:quote, fname)), ex0)
            quote
                results = $thecall
                length(results) == 1 ? results[1] : results
            end
        end
    end
end

"""
    @which

Applied to a function or macro call, it evaluates the arguments to the specified call, and
returns the `Method` object for the method that would be called for those arguments. Applied
to a variable, it returns the module in which the variable was bound. It calls out to the
`which` function.
"""
:@which

"""
    @less

Evaluates the arguments to the function or macro call, determines their types, and calls the `less`
function on the resulting expression.
"""
:@less

"""
    @edit

Evaluates the arguments to the function or macro call, determines their types, and calls the `edit`
function on the resulting expression.
"""
:@edit

"""
    @functionloc

Applied to a function or macro call, it evaluates the arguments to the specified call, and
returns a tuple `(filename,line)` giving the location for the method that would be called for those arguments.
It calls out to the `functionloc` function.
"""
:@functionloc

"""
    @code_typed

Evaluates the arguments to the function or macro call, determines their types, and calls
[`code_typed`](@ref) on the resulting expression.
"""
:@code_typed

"""
    @code_warntype

Evaluates the arguments to the function or macro call, determines their types, and calls
[`code_warntype`](@ref) on the resulting expression.
"""
:@code_warntype

"""
    @code_lowered

Evaluates the arguments to the function or macro call, determines their types, and calls
[`code_lowered`](@ref) on the resulting expression.
"""
:@code_lowered

"""
    @code_llvm

Evaluates the arguments to the function or macro call, determines their types, and calls
[`code_llvm`](@ref) on the resulting expression.
"""
:@code_llvm

"""
    @code_native

Evaluates the arguments to the function or macro call, determines their types, and calls
[`code_native`](@ref) on the resulting expression.
"""
:@code_native

function type_close_enough(@nospecialize(x), @nospecialize(t))
    x == t && return true
    return (isa(x,DataType) && isa(t,DataType) && x.name === t.name &&
            !_isleaftype(t) && x <: t) ||
           (isa(x,Union) && isa(t,DataType) && (type_close_enough(x.a, t) || type_close_enough(x.b, t)))
end

# `methodswith` -- shows a list of methods using the type given
"""
    methodswith(typ[, module or function][, showparents::Bool=false])

Return an array of methods with an argument of type `typ`.

The optional second argument restricts the search to a particular module or function
(the default is all modules, starting from Main).

If optional `showparents` is `true`, also return arguments with a parent type of `typ`,
excluding type `Any`.
"""
function methodswith(t::Type, f::Function, showparents::Bool=false, meths = Method[])
    for d in methods(f)
        if any(function (x)
                   let x = rewrap_unionall(x, d.sig)
                       (type_close_enough(x, t) ||
                        (showparents ? (t <: x && (!isa(x,TypeVar) || x.ub != Any)) :
                         (isa(x,TypeVar) && x.ub != Any && t == x.ub)) &&
                        x != Any)
                   end
               end,
               unwrap_unionall(d.sig).parameters)
            push!(meths, d)
        end
    end
    return meths
end

function methodswith(t::Type, m::Module, showparents::Bool=false)
    meths = Method[]
    for nm in names(m)
        if isdefined(m, nm)
            f = getfield(m, nm)
            if isa(f, Function)
                methodswith(t, f, showparents, meths)
            end
        end
    end
    return unique(meths)
end

function methodswith(t::Type, showparents::Bool=false)
    meths = Method[]
    mainmod = Main
    # find modules in Main
    for nm in names(mainmod)
        if isdefined(mainmod, nm)
            mod = getfield(mainmod, nm)
            if isa(mod, Module)
                append!(meths, methodswith(t, mod, showparents))
            end
        end
    end
    return unique(meths)
end

# file downloading

downloadcmd = nothing
if Sys.iswindows()
    function download(url::AbstractString, filename::AbstractString)
        res = ccall((:URLDownloadToFileW,:urlmon),stdcall,Cuint,
                    (Ptr{Void},Cwstring,Cwstring,Cuint,Ptr{Void}),C_NULL,url,filename,0,C_NULL)
        if res != 0
            error("automatic download failed (error: $res): $url")
        end
        filename
    end
else
    function download(url::AbstractString, filename::AbstractString)
        global downloadcmd
        if downloadcmd === nothing
            for checkcmd in (:curl, :wget, :fetch)
                if success(pipeline(`which $checkcmd`, DevNull))
                    downloadcmd = checkcmd
                    break
                end
            end
        end
        if downloadcmd == :wget
            try
                run(`wget -O $filename $url`)
            catch
                rm(filename)  # wget always creates a file
                rethrow()
            end
        elseif downloadcmd == :curl
            run(`curl -g -L -f -o $filename $url`)
        elseif downloadcmd == :fetch
            run(`fetch -f $filename $url`)
        else
            error("no download agent available; install curl, wget, or fetch")
        end
        filename
    end
end
function download(url::AbstractString)
    filename = tempname()
    download(url, filename)
end

"""
    download(url::AbstractString, [localfile::AbstractString])

Download a file from the given url, optionally renaming it to the given local file name.
Note that this function relies on the availability of external tools such as `curl`, `wget`
or `fetch` to download the file and is provided for convenience. For production use or
situations in which more options are needed, please use a package that provides the desired
functionality instead.
"""
download(url, filename)

# workspace management

"""
    workspace()

Replace the top-level module (`Main`) with a new one, providing a clean workspace. The
previous `Main` module is made available as `LastMain`. A previously-loaded package can be
accessed using a statement such as `using LastMain.Package`.

This function should only be used interactively.
"""
function workspace()
    last = Core.Main # ensure to reference the current Main module
    b = Base # this module
    ccall(:jl_new_main_module, Any, ()) # make Core.Main a new baremodule
    m = Core.Main # now grab a handle to the new Main module
    ccall(:jl_add_standard_imports, Void, (Any,), m)
    eval(m, Expr(:toplevel,
                 :(const Base = $b),
                 :(const LastMain = $last),
                 :(using Base.MainInclude)))
    empty!(package_locks)
    return m
end

# testing

"""
    runtests([tests=["all"] [, numcores=ceil(Int, Sys.CPU_CORES / 2) ]])

Run the Julia unit tests listed in `tests`, which can be either a string or an array of
strings, using `numcores` processors. (not exported)
"""
function runtests(tests = ["all"], numcores = ceil(Int, Sys.CPU_CORES / 2))
    if isa(tests,AbstractString)
        tests = split(tests)
    end
    ENV2 = copy(ENV)
    ENV2["JULIA_CPU_CORES"] = "$numcores"
    try
        run(setenv(`$(julia_cmd()) $(joinpath(JULIA_HOME,
            Base.DATAROOTDIR, "julia", "test", "runtests.jl")) $tests`, ENV2))
    catch
        buf = PipeBuffer()
        versioninfo(buf)
        error("A test has failed. Please submit a bug report (https://github.com/JuliaLang/julia/issues)\n" *
              "including error messages above and the output of versioninfo():\n$(read(buf, String))")
    end
end

# testing

"""
    whos(io::IO=STDOUT, m::Module=Main, pattern::Regex=r"")

Print information about exported global variables in a module, optionally restricted to those matching `pattern`.

The memory consumption estimate is an approximate lower bound on the size of the internal structure of the object.
"""
function whos(io::IO=STDOUT, m::Module=Main, pattern::Regex=r"")
    maxline = displaysize(io)[2]
    line = zeros(UInt8, maxline)
    head = PipeBuffer(maxline + 1)
    for v in sort!(names(m))
        s = string(v)
        if isdefined(m, v) && ismatch(pattern, s)
            value = getfield(m, v)
            @printf head "%30s " s
            try
                if value ∈ (Base, Main, Core)
                    print(head, "              ")
                else
                    bytes = summarysize(value)
                    if bytes < 10_000
                        @printf(head, "%6d bytes  ", bytes)
                    else
                        @printf(head, "%6d KB     ", bytes ÷ (1024))
                    end
                end
                print(head, summary(value))
            catch e
                print(head, "#=ERROR: unable to show value=#")
            end
            newline = search(head, UInt8('\n')) - 1
            if newline < 0
                newline = nb_available(head)
            end
            if newline > maxline
                newline = maxline - 1 # make space for ...
            end
            line = resize!(line, newline)
            line = read!(head, line)

            write(io, line)
            if nb_available(head) > 0 # more to read? replace with ...
                print(io, '\u2026') # hdots
            end
            println(io)
            seekend(head) # skip the rest of the text
        end
    end
end
whos(m::Module, pat::Regex=r"") = whos(STDOUT, m, pat)
whos(pat::Regex) = whos(STDOUT, Main, pat)

max_stars_repo_name: visr/julia

max_issues_repo_name: visr/julia

max_forks_repo_name: visr/julia


text:
include("transvections.jl")
include("gersten_relations.jl")

function SpecialAutomorphismGroup(F::FreeGroup; ordering = KnuthBendix.LenLex, kwargs...)

    n = length(alphabet(F)) ÷ 2
    A, rels = gersten_relations(n, commutative = false)
    S = KnuthBendix.letters(A)[1:2(n^2-n)]

    maxrules = 1000*n

    rws = KnuthBendix.RewritingSystem(rels, ordering(A))
    KnuthBendix.knuthbendix!(rws; maxrules=maxrules, kwargs...)
    return AutomorphismGroup(F, S, rws, ntuple(i -> gens(F, i), n))
end

KnuthBendix.alphabet(G::AutomorphismGroup{<:FreeGroup}) = alphabet(rewriting(G))

function relations(G::AutomorphismGroup{<:FreeGroup})
    n = length(alphabet(object(G))) ÷ 2
    return last(gersten_relations(n, commutative = false))
end

ext: jl


text:
module PhaseRetrieval
# using LinearAlgebra
using Statistics
using FFTW
using SampledDomains: CartesianDomain2D, dualRange
import SampledDomains: make_centered_domain2D

import AlternatingProjections: Problem, Algorithm, solve, TwoSetsFP, IterativeAlgorithm, initial, tolerance, maxit, keephistory, snapshots # these we will need to change
using AlternatingProjections

export SHdiversity
export psf, phwrap, subpsf
export removepiston, twinphase

# export Linear

include("types.jl")
include("utils.jl")
include("ShackHartmann/SHphase.jl")
include("PSF.jl")

include("PRproblem.jl")

export AbstractPRproblem, PRproblem, appsftoPR, solve, GS, GSparam

end


text:
export amplitudeDAC, frequencyDAC, phaseDAC, modeDAC, amplitudeDACNext,
       DCSignDAC, signalTypeDAC, offsetDAC, jumpSharpnessDAC, passPDMToFastDAC,
       waveforms, DACPerformanceData

struct DACPerformanceData
  uDeltaControl::UInt8
  uDeltaSet::UInt8
  minDeltaControl::UInt8
  maxDeltaSet::UInt8
end

function passPDMToFastDAC(rp::RedPitaya, val::Bool)
  valStr = val ? "ON" : "OFF"
  send(rp, string("RP:PassPDMToFastDAC ", valStr))
end
passPDMToFastDAC(rp::RedPitaya) = occursin("ON", query(rp,"RP:PassPDMToFastDAC?"))

function amplitudeDAC(rp::RedPitaya, channel, component)
  command = string("RP:DAC:CH", Int(channel)-1, ":COMP", Int(component)-1, ":AMP?")
  return query(rp, command, Float64)
end
function amplitudeDAC(rp::RedPitaya, channel, component, value)
  if value > 1.0
    error("$value is larger than 1.0 V!")
  end
  command = string("RP:DAC:CH", Int(channel)-1, ":COMP",
                   Int(component)-1, ":AMP ", Float64(value))
  return send(rp, command)
end

function amplitudeDACNext(rp::RedPitaya, channel, component)
  command = string("RP:DAC:CH", Int(channel)-1, ":COMP", Int(component)-1, ":Next:AMP?")
  return query(rp, command, Float64)
end
function amplitudeDACNext(rp::RedPitaya, channel, component, value)
  if value > 1.0
    error("$value is larger than 1.0 V!")
  end
  command = string("RP:DAC:CH", Int(channel)-1, ":COMP",
                   Int(component)-1, ":Next:AMP ", Float64(value))
  return send(rp, command)
end

function offsetDAC(rp::RedPitaya, channel)
  command = string("RP:DAC:CH", Int(channel)-1, ":OFF?")
  return query(rp, command, Float64)
end
function offsetDAC(rp::RedPitaya, channel, value)
  if value > 1.0
    error("$value is larger than 1.0 V!")
  end
  command = string("RP:DAC:CH", Int(channel)-1, ":OFF ", Float64(value))
  return send(rp, command)
end

function frequencyDAC(rp::RedPitaya, channel, component)
  command = string("RP:DAC:CH", Int(channel)-1, ":COMP", Int(component)-1, ":FREQ?")
  return query(rp, command, Float64)
end
function frequencyDAC(rp::RedPitaya, channel, component, value)
  command = string("RP:DAC:CH", Int(channel)-1, ":COMP",
                   Int(component)-1, ":FREQ ", Float64(value))
  send(rp, command)
end

function phaseDAC(rp::RedPitaya, channel, component)
  command = string("RP:DAC:CH", Int(channel)-1, ":COMP", Int(component)-1, ":PHA?")
  return query(rp, command, Float64)
end
function phaseDAC(rp::RedPitaya, channel, component, value)
  command = string("RP:DAC:CH", Int(channel)-1, ":COMP",
                   Int(component)-1, ":PHA ", Float64(value))
  send(rp, command)
end

function jumpSharpnessDAC(rp::RedPitaya, channel)
  command = string("RP:DAC:CH", Int(channel)-1, ":JumpSharpness?")
  return query(rp, command, Float64)
end
function jumpSharpnessDAC(rp::RedPitaya, channel, value)
  command = string("RP:DAC:CH", Int(channel)-1, ":JumpSharpness ", Float64(value))
  send(rp, command)
end

#"STANDARD" or "AWG" (not yet supported)
function modeDAC(rp::RedPitaya)
  modeDAC_(rp, 1)
  modeDAC_(rp, 2)
end
function modeDAC(rp::RedPitaya, mode::String)
  modeDAC_(rp, mode, 1)
  modeDAC_(rp, mode, 2)
end

function modeDAC_(rp::RedPitaya, channel=1)
  #return query(rp, "RP:DAC:CH", Int(channel)-1,":MODe?")[2:end-1]
  return query(rp, "RP:DAC:MODe?")[2:end-1]
end
function modeDAC_(rp::RedPitaya, mode::String, channel=1)
  #send(rp, string("RP:DAC:CH", Int(channel)-1, ":MODe ", mode))
  send(rp, string("RP:DAC:MODe ", mode))
end

function signalTypeDAC(rp::RedPitaya, channel)
  command = string("RP:DAC:CH", Int(channel)-1, ":SIGnaltype?")
  @show command
  return query(rp, command)
end

waveforms() = ["SINE","SQUARE","TRIANGLE","SAWTOOTH"]

function signalTypeDAC(rp::RedPitaya, channel, sigType::String)
  if !(sigType in ["SINE","SQUARE","TRIANGLE","SAWTOOTH"] )
    error("Signal type $sigType not supported!")
  end

  command = string("RP:DAC:CH", Int(channel)-1, ":SIGnaltype ", sigType)
  return send(rp, command)
end

function DCSignDAC(rp::RedPitaya, channel)
  command = string("RP:DAC:CH", Int(channel)-1, ":SIGn?")
  return query(rp, command, Int64) == "POSITIVE" ? 1 : -1
end

function DCSignDAC(rp::RedPitaya, channel, sign::Integer)
  command = string("RP:DAC:CH", Int(channel)-1, ":SIGn ",
             sign > 0 ? "POSITIVE" : "NEGATIVE")
  return send(rp, command)
end

function readDACPerformanceData(rp::RedPitaya)
  perf = read!(rp.dataSocket, Array{UInt8}(undef, 4))
  return DACPerformanceData(perf[1], perf[2], perf[3], perf[4])
end

max_stars_repo_name: Fypsilonn/RedPitayaDAQServer

max_issues_repo_name: Fypsilonn/RedPitayaDAQServer

max_forks_repo_name: Fypsilonn/RedPitayaDAQServer


text:
# Autogenerated wrapper script for FXdiv_jll for any
export fxdiv_h

JLLWrappers.@generate_wrapper_header("FXdiv")
JLLWrappers.@declare_file_product(fxdiv_h)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_file_product(
        fxdiv_h,
        "include/fxdiv.h",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# ALL TESTS THAT CONCERN UNITCELLS

# include
using LatticePhysics
using Base.Test

# begin the Unitcell testblock
unitcell_testset = @testset "Unitcells" begin

################################################################################
#
#   TYPE DEFINITIONS OF OBJECT CLASSES IN JULIA
#
################################################################################

# begin the testset
@testset "General Unitcell functions" begin

    # get a test unitcell
    unitcell = getUnitcellHoneycomb()

    @testset "File access with JLD" begin
        # test saving
        @test typeof(saveUnitcell(unitcell))==String
        # test loading
        @test typeof(loadUnitcell(unitcell))==String
        @test typeof(loadUnitcell(unitcell.filename))==Unitcell
    end;

    @testset "Test and Information functions" begin
        # test printing information
        @test printInfo(unitcell)==nothing
        @test printInfo(unitcell, detailed=true)==nothing
        # test testing method
        @test typeof(testUnitcell(unitcell, 2, 2))==Bool
    end;

    @testset "Connectivity information functions" begin
        # testing connectivity information
        @test typeof(getConnectivityList(unitcell)) == Array{Array,1}
        # testing connection information
        @test typeof(getConnectionList(unitcell)) == Array{Array,1}
        # testing connection strength information
        @test typeof(getConnectionStrengthList(unitcell)) == Array{Any,1}
    end;

# end the testset here
end;

################################################################################
#
#   IMPLEMENTATIONS OF DIFFERENT UNITCELLS in 2D and 3D
#
################################################################################

# begin the testset
@testset "Implemented Unitcells in 2D" begin

    #   2D UNITCELLS
    #   - SQUARE / RECTANGLE
    #       - getUnitcellSquare
    @testset "2D Square unitcell" begin
        @test testUnitcell(getUnitcellSquare(), 2, 2)
        @test testUnitcell(getUnitcellSquare(1), 2, 2)
        @test testUnitcell(getUnitcellSquare(2), 2, 2)
        @test testUnitcell(getUnitcellSquare(3), 2, 2)
    end;
    #       - getUnitcellExtendedSquare
    @testset "2D Extended Square unitcell" begin
        @test testUnitcell(getUnitcellExtendedSquare(), 2, 2)
        @test testUnitcell(getUnitcellExtendedSquare(1), 2, 2)
    end;
    #       - getUnitcellCheckerboard
    @testset "2D Checkerboard unitcell" begin
        @test testUnitcell(getUnitcellCheckerboard(), 2, 2)
        @test testUnitcell(getUnitcellCheckerboard(1), 2, 2)
    end;
    #       - getUnitcellShastrySutherland
    @testset "2D Shastry-Sutherland unitcell" begin
        @test testUnitcell(getUnitcellShastrySutherland(), 2, 2)
        @test testUnitcell(getUnitcellShastrySutherland(1), 2, 2)
        @test testUnitcell(getUnitcellShastrySutherland(4), 2, 2)
    end;
    #       - getUnitcellAdvancedSquare
    @testset "2D Advanced Square unitcell" begin
        @test testUnitcell(getUnitcellAdvancedSquare(), 2, 2)
        @test testUnitcell(getUnitcellAdvancedSquare(1), 2, 2)
    end;
    #       - getUnitcellSquareOctagon
    @testset "2D Square-octagon unitcell" begin
        @test testUnitcell(getUnitcellSquareOctagon(), 2, 2)
        @test testUnitcell(getUnitcellSquareOctagon(1), 2, 2)
        @test testUnitcell(getUnitcellSquareOctagon(4), 2, 2)
        @test testUnitcell(getUnitcellSquareOctagon(2), 2, 2)
        @test testUnitcell(getUnitcellSquareOctagon(5), 2, 2)
    end;
    #       - getUnitcellBCC2D
    @testset "2D BCC (square) unitcell" begin
        @test testUnitcell(getUnitcellBCC2D(), 2, 2)
        @test testUnitcell(getUnitcellBCC2D(1), 2, 2)
    end;
    #       - getUnitcellFullyConnectedSquare
    @testset "2D Fully connected square unitcell" begin
        @test testUnitcell(getUnitcellFullyConnectedSquare(), 2, 2)
        @test testUnitcell(getUnitcellFullyConnectedSquare(1), 2, 2)
    end;

    #   - TRIANGULAR
    #       - getUnitcellTriangular
    @testset "2D Triangular unitcell" begin
        @test testUnitcell(getUnitcellTriangular(), 2, 2)
        @test testUnitcell(getUnitcellTriangular(1), 2, 2)
        @test testUnitcell(getUnitcellTriangular(3), 2, 2)
    end;
    #       - getUnitcellHoneycomb
    @testset "2D Honeycomb unitcell" begin
        @test testUnitcell(getUnitcellHoneycomb(), 2, 2)
        @test testUnitcell(getUnitcellHoneycomb(1), 2, 2)
        @test testUnitcell(getUnitcellHoneycomb(4), 2, 2)
        @test testUnitcell(getUnitcellHoneycomb(2), 2, 2)
        @test testUnitcell(getUnitcellHoneycomb(5), 2, 2)
        @test testUnitcell(getUnitcellHoneycomb(3), 2, 2)
    end;
    #       - getUnitcellKagome
    @testset "2D Kagome unitcell" begin
        @test testUnitcell(getUnitcellKagome(), 2, 2)
        @test testUnitcell(getUnitcellKagome(1), 2, 2)
    end;
    #       - getUnitcellKagomeMinus
    @testset "2D Kagome Minus unitcell" begin
        @test testUnitcell(getUnitcellKagomeMinus(), 2, 2)
        @test testUnitcell(getUnitcellKagomeMinus(1), 2, 2)
        @test testUnitcell(getUnitcellKagomeMinus(4), 2, 2)
    end;
    #       - getUnitcellHoneycombXXX
    @testset "2D Honeycomb-XXX unitcell" begin
        @test testUnitcell(getUnitcellHoneycombXXX(), 2, 2)
        @test testUnitcell(getUnitcellHoneycombXXX(1), 2, 2)
        @test testUnitcell(getUnitcellHoneycombXXX(2), 2, 2)
    end;

# end the testset here
end;

# begin the testset
@testset "Implemented Unitcells in 3D" begin

    #   3D UNITCELLS
    #   - CUBIC / FCC
    #       - getUnitcellDiamond
    @testset "3D Diamond unitcell" begin
        @test testUnitcell(getUnitcellDiamond(), 3, 3)
        @test testUnitcell(getUnitcellDiamond(1), 3, 3)
        @test testUnitcell(getUnitcellDiamond(2), 3, 3)
        @test testUnitcell(getUnitcellDiamond(3), 3, 3)
    end;
    #       - getUnitcellBCC
    @testset "3D BCC unitcell" begin
        @test testUnitcell(getUnitcellBCC(), 3, 3)
        @test testUnitcell(getUnitcellBCC(1), 3, 3)
    end;
    #       - getUnitcellPyrochlore
    @testset "3D Pyrochlore unitcell" begin
        @test testUnitcell(getUnitcellPyrochlore(), 3, 3)
        @test testUnitcell(getUnitcellPyrochlore(1), 3, 3)
    end;

    #   - (X,3)y FAMILY
    #       - getUnitcell_8_3_a
    @testset "3D (8,3)a unitcell" begin
        @test testUnitcell(getUnitcell_8_3_a(), 3, 3)
        @test testUnitcell(getUnitcell_8_3_a(1), 3, 3)
        @test testUnitcell(getUnitcell_8_3_a(4), 3, 3)
    end;
    #       - getUnitcell_8_3_b
    @testset "3D (8,3)b unitcell" begin
        @test testUnitcell(getUnitcell_8_3_b(), 3, 3)
        @test testUnitcell(getUnitcell_8_3_b(1), 3, 3)
        @test testUnitcell(getUnitcell_8_3_b(4), 3, 3)
    end;
    #       - getUnitcell_8_3_c
    @testset "3D (8,3)c unitcell" begin
        @test testUnitcell(getUnitcell_8_3_c(), 3, 3)
        @test testUnitcell(getUnitcell_8_3_c(1), 3, 3)
        @test testUnitcell(getUnitcell_8_3_c(4), 3, 3)
    end;
    #       - getUnitcell_8_3_n
    @testset "3D (8,3)n unitcell" begin
        @test testUnitcell(getUnitcell_8_3_n(), 3, 3)
        @test testUnitcell(getUnitcell_8_3_n(1), 3, 3)
        @test testUnitcell(getUnitcell_8_3_n(4), 3, 3)
    end;
    #       - getUnitcell_9_3_a
    @testset "3D (9,3)a unitcell" begin
        @test testUnitcell(getUnitcell_9_3_a(), 3, 3)
        @test testUnitcell(getUnitcell_9_3_a(1), 3, 3)
        @test testUnitcell(getUnitcell_9_3_a(4), 3, 3)
        @test testUnitcell(getUnitcell_9_3_a(2), 3, 3)
        @test testUnitcell(getUnitcell_9_3_a(5), 3, 3)
    end;
    #       - getUnitcell_10_3_a / getUnitcellHyperoctagon
    @testset "3D (10,3)a / hyperoctagon unitcell" begin
        @test testUnitcell(getUnitcell_10_3_a(), 3, 3)
        @test testUnitcell(getUnitcell_10_3_a(1), 3, 3)
        @test testUnitcell(getUnitcell_10_3_a(4), 3, 3)
        @test testUnitcell(getUnitcellHyperoctagon(), 3, 3)
        @test testUnitcell(getUnitcellHyperoctagon(1), 3, 3)
        @test testUnitcell(getUnitcellHyperoctagon(4), 3, 3)
    end;
    #       - getUnitcell_10_3_b / getUnitcellHyperhoneycomb
    @testset "3D (10,3)b / hyperhoneycomb unitcell" begin
        @test testUnitcell(getUnitcell_10_3_b(), 3, 3)
        @test testUnitcell(getUnitcell_10_3_b(1), 3, 3)
        @test testUnitcell(getUnitcell_10_3_b(4), 3, 3)
        @test testUnitcell(getUnitcell_10_3_b(2), 3, 3)
        @test testUnitcell(getUnitcell_10_3_b(5), 3, 3)
        @test testUnitcell(getUnitcellHyperhoneycomb(), 3, 3)
        @test testUnitcell(getUnitcellHyperhoneycomb(1), 3, 3)
        @test testUnitcell(getUnitcellHyperhoneycomb(4), 3, 3)
        @test testUnitcell(getUnitcellHyperhoneycomb(2), 3, 3)
        @test testUnitcell(getUnitcellHyperhoneycomb(5), 3, 3)
    end;
    #       - getUnitcell_10_3_c
    @testset "3D (10,3)c unitcell" begin
        @test testUnitcell(getUnitcell_10_3_c(), 3, 3)
        @test testUnitcell(getUnitcell_10_3_c(1), 3, 3)
        @test testUnitcell(getUnitcell_10_3_c(4), 3, 3)
    end;
    #       - getUnitcell_10_3_d
    @testset "3D (10,3)d unitcell" begin
        @test testUnitcell(getUnitcell_10_3_d(), 3, 3)
        @test testUnitcell(getUnitcell_10_3_d(1), 3, 3)
        @test testUnitcell(getUnitcell_10_3_d(4), 3, 3)
    end;

# end the testset here
end;

################################################################################
#
#   IMPLEMENTATIONS OF DIFFERENT UNITCELL CONSTRUCTION RELATED STUFF
#
################################################################################

# begin the testset
@testset "Unitcell construction related" begin

    # get test unitcells
    unitcell_2d = getUnitcellHoneycomb()
    unitcell_3d = getUnitcellDiamond()

    # test getting unitcell from sites
    @testset "Getting unitcells from collection of sites" begin
        @test testUnitcell(getUnitcellFromSites(unitcell_2d.basis, unitcell_2d.lattice_vectors), 2, 2)
        @test testUnitcell(getUnitcellFromSites(unitcell_3d.basis, unitcell_3d.lattice_vectors), 3, 3)
    end;

    # print unitcell code
    @testset "Printing unitcell generating code" begin
        @test printUnitcellGeneratingCode(unitcell_2d)==nothing
        @test printUnitcellGeneratingCode(unitcell_3d)==nothing
    end;

# end the testset here
end;

# end the unitcell test block here
end;


text:
using Test
using SparseArrays
using LinearAlgebra
using LinearAlgebraicRepresentation
using LinearAlgebraicRepresentation.Arrangement
Lar = LinearAlgebraicRepresentation

@testset "Edge fragmentation tests" begin
    V = [2 2; 4 2; 3 3.5; 1 3; 5 3; 1 2; 5 2]
    EV = SparseArrays.sparse(Array{Int8, 2}([
        [1 1 0 0 0 0 0] #1->1,2
        [0 1 1 0 0 0 0] #2->2,3
        [1 0 1 0 0 0 0] #3->1,3
        [0 0 0 1 1 0 0] #4->4,5
        [0 0 0 0 0 1 1] #5->6,7
    ]))

    @testset "intersect_edges" begin
        inters1 = Lar.Arrangement.intersect_edges(V, EV[5, :], EV[1, :])
        inters2 = Lar.Arrangement.intersect_edges(V, EV[1, :], EV[4, :])
        inters3 = Lar.Arrangement.intersect_edges(V, EV[1, :], EV[2, :])
        @test inters1 == [([2. 2.], 1/4),([4. 2.], 3/4)]
        @test inters2 == []
        @test inters3 == [([4. 2.], 1)]
    end

#    @testset "frag_edge" begin
#        rV, rEV = Lar.Arrangement.frag_edge(V, EV, 5, [1,2,3,4,5])
#        @test rV == [1.0 2.0; 5.0 2.0; 2.0 2.0; 4.0 2.0; 4.0 2.0; 2.0 2.0]
#        @test Matrix(rEV) == [1 0 0 0 0 1;
#                             0 0 0 0 1 1;
#                             0 1 0 0 1 0]
#    end
end

@testset "merge_vertices test set" begin
    n0 = 1e-12
    n1l = 1-1e-12
    n1u = 1+1e-12
    V = [ n0  n0; -n0  n0;  n0 -n0; -n0 -n0;
          n0 n1u; -n0 n1u;  n0 n1l; -n0 n1l;
         n1u n1u; n1l n1u; n1u n1l; n1l n1l;
         n1u  n0; n1l  n0; n1u -n0; n1l -n0]
    EV = Int8[1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0;
              0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0;
              0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0;
              0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0;
              0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0;
              0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0;
              0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0;
              0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0;
              0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0;
              0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0;
              0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0;
              0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1;
              1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0;
              0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0;
              0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0;
              0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1]
    EV = sparse(EV)
    V, EV = Lar.Arrangement.merge_vertices!(V, EV, [])

    @test V == [n0 n0; n0 n1u; n1u n1u; n1u n0]
    @test Matrix(EV) == [1 1 0 0;
                       0 1 1 0;
                       0 0 1 1;
                       1 0 0 1]
end

@testset "biconnected_components test set" begin
    EV = Int8[0 0 0 1 0 0 0 0 0 0 1 0; #1
              0 0 1 0 0 1 0 0 0 0 0 0; #2
              0 0 0 0 0 0 1 0 0 1 0 0; #3
              1 0 0 0 1 0 0 0 0 0 0 0; #4
              0 0 0 1 0 0 0 1 0 0 0 0; #5
              0 0 1 0 0 0 0 0 1 0 0 0; #6
              0 1 0 0 0 0 0 0 0 1 0 0; #7
              0 0 0 0 1 0 0 0 0 0 0 1; #8
              0 0 0 0 0 0 0 1 0 0 1 0; #9
              0 0 0 0 0 1 0 0 1 0 0 0; #10
              0 1 0 0 0 0 1 0 0 0 0 0; #11
              0 0 0 0 1 0 0 0 0 0 1 0; #12
              0 0 0 0 1 0 0 0 1 0 0 0] #13
    EV = sparse(EV)

    bc = Lar.Arrangement.biconnected_components(EV)
    bc = Set(map(Set, bc))

    @test bc == Set([Set([1,5,9]), Set([2,6,10]), Set([3,7,11])])
end

@testset "Face creation" begin
    @testset "External cell individuation" begin
        V = [ .5 .5;  1.5   1;  1.5  2;
             2.5  2;  2.5   1;  3.5 .5;
             3.5  3;    2 2.5;   .5  3]

        EV = Int8[-1  1  0  0  0  0  0  0  0;
                   0 -1  1  0  0  0  0  0  0;
                   0  0 -1  1  0  0  0  0  0;
                   0  0  0 -1  1  0  0  0  0;
                   0  0  0  0 -1  1  0  0  0;
                   0  0  0  0  0 -1  1  0  0;
                   0  0  0  0  0  0 -1  1  0;
                   0  0  0  0  0  0  0 -1  1;
                  -1  0  0  0  0  0  0  0  1;
                   0 -1  0  0  1  0  0  0  0]
        EV = sparse(EV)

        FE = Int8[ 0 -1 -1 -1  0  0  0  0  0  1;
                   1  1  1  1  1  1  1  1 -1  0;
                  -1  0  0  0 -1 -1 -1 -1  1 -1]
        FE = sparse(FE)

        @test Lar.Arrangement.get_external_cycle(V, EV, FE) == 3
    end

    @testset "Containment test" begin
        V = [  0   0;    4   0;    4   2;   2   4;  0 4;
              .5  .5;  2.5  .5;  2.5 2.5;  .5 2.5;
               1   1;  1.5   1;    1   2;
               2   1;    2   2;  1.5   2;
             3.5 3.5;    3 3.5;  3.5   3]
        EV1 = Int8[ 0  0  0  0  0  0  0  0  0 -1  1  0  0  0  0  0  0  0;
                    0  0  0  0  0  0  0  0  0  0 -1  1  0  0  0  0  0  0;
                    0  0  0  0  0  0  0  0  0 -1  0  1  0  0  0  0  0  0]
        EV2 = Int8[ 0  0  0  0  0  0  0  0  0  0  0  0 -1  1  0  0  0  0;
                    0  0  0  0  0  0  0  0  0  0  0  0  0 -1  1  0  0  0;
                    0  0  0  0  0  0  0  0  0  0  0  0 -1  0  1  0  0  0]
        EV3 = Int8[ 0  0  0  0  0 -1  1  0  0  0  0  0  0  0  0  0  0  0;
                    0  0  0  0  0  0 -1  1  0  0  0  0  0  0  0  0  0  0;
                    0  0  0  0  0  0  0 -1  1  0  0  0  0  0  0  0  0  0;
                    0  0  0  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  0]
        EV4 = Int8[-1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;
                    0 -1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;
                    0  0 -1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0;
                    0  0  0 -1  1  0  0  0  0  0  0  0  0  0  0  0  0  0;
                   -1  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0]
        EV5 = Int8[ 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0 -1  1  0;
                    0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0 -1  1;
                    0  0  0  0  0  0  0  0  0  0  0  0  0  0  0 -1  0  1]
        EVs = map(sparse, [EV1, EV2, EV3, EV4, EV5])

        shell1 = Int8[-1 -1  1];
        shell2 = Int8[-1 -1  1];
        shell3 = Int8[-1 -1 -1  1];
        shell4 = Int8[-1 -1 -1 -1  1];
        shell5 = Int8[-1 -1  1];
        shells = map(sparsevec, [shell1, shell2, shell3, shell4, shell5])

        shell_bboxes = []
        n = 5
        for i in 1:n
            vs_indexes = (abs.(EVs[i]')*abs.(shells[i])).nzind
            push!(shell_bboxes, Lar.bbox(V[vs_indexes, :]))
        end

        graph = Lar.Arrangement.pre_containment_test(shell_bboxes)
        @test graph == [0 0 1 1 0; 0 0 1 1 0; 0 0 0 1 0; 0 0 0 0 0; 0 0 0 1 0]

        graph = Lar.Arrangement.prune_containment_graph(n, V, EVs, shells, graph)
        @test graph == [0 0 1 1 0; 0 0 1 1 0; 0 0 0 1 0; 0 0 0 0 0; 0 0 0 0 0]
    end

    @testset "Transitive reduction" begin
        graph = [0 0 1 1 0; 0 0 1 1 0; 0 0 0 1 0; 0 0 0 0 0; 0 0 0 0 0]
        Lar.Arrangement.transitive_reduction!(graph)
        @test graph == [0 0 1 0 0; 0 0 1 0 0; 0 0 0 1 0; 0 0 0 0 0; 0 0 0 0 0]
    end

    @testset "Cell merging" begin
        graph = [0 1; 0 0]
        V = [.25 .25; .75 .25; .75 .75; .25 .75;
               0   0;   1   0;   1   1;   0   1]
        EV1 = Int8[-1  1  0  0  0  0  0  0;
                    0 -1  1  0  0  0  0  0;
                    0  0 -1  1  0  0  0  0;
                   -1  0  0  1  0  0  0  0;
                   -1  0  1  0  0  0  0  0]
        EV2 = Int8[ 0  0  0  0 -1  1  0  0;
                    0  0  0  0  0 -1  1  0;
                    0  0  0  0  0  0 -1  1;
                    0  0  0  0 -1  0  0  1]
        EVs = map(sparse, [EV1, EV2])

        shell1 = Int8[-1 -1 -1  1  0]
        shell2 = Int8[-1 -1 -1  1]
        shells = map(sparsevec, [shell1, shell2])

        boundary1 = Int8[ 1  1  0  0 -1;
                          0  0  1 -1  1]
        boundary2 = Int8[ 1  1  1 -1]
        boundaries = map(sparse, [boundary1, boundary2])

        shell_bboxes = []
        n = 2
        for i in 1:n
            vs_indexes = (abs.(EVs[i]')*abs.(shells[i])).nzind
            push!(shell_bboxes, Lar.bbox(V[vs_indexes, :]))
        end

        EV, FE = Lar.Arrangement.cell_merging(2, graph, V, EVs, boundaries, shells, shell_bboxes)

        selector = sparse(ones(Int8, 1, 3))

        @test selector*FE == [0  0  0  0  0  1  1  1 -1]
    end
end


text:
Base.@kwdef struct PackageState
    timestamp::DateTime
    load_path::Vector{String}
    id::Base.PkgId
    dir::String
    tree_hash::String
    manifest_tree_hash::Union{Nothing,String}
    project::Union{Nothing,String}
end

function Base.:(==)(s1::PackageState, s2::PackageState)
    s1.dir == s2.dir &&
    s1.id == s2.id &&
    s1.project == s2.project &&
    s1.load_path == s2.load_path &&
    s1.tree_hash == s2.tree_hash &&
    s1.manifest_tree_hash == s2.manifest_tree_hash
end

const row_names = ["Timestamp",
                   "Load path",
                   "Package ID",
                   "Source path",
                   "Tree hash",
                   "Manifest t.h.",
                   "Project"]

function tabledata(s::PackageState)
    data = [Dates.format(s.timestamp, "yyyy-mm-dd HH:MM"),
            s.load_path,
            s.id,
            s.dir,
            s.tree_hash,
            s.manifest_tree_hash,
            s.project]
    return data
end

printtable(states::Vararg{PackageState, N}; kwargs...) where N = printtable(tabledata.(states)...; kwargs...)
function printtable(datavectors::Vararg{AbstractVector, N}; kwargs...) where N
    println()
    highlighters = ()
    if N > 1
        highlighters = ( Highlighter(  (data, i, j) -> (i > 1) && length(unique(data[i,:]))> 1,
                                        crayon"fg:red"),)
    end
    pretty_table(hcat(datavectors...),
                 row_names = row_names,
                 autowrap = true,
                 linebreaks = true,
                 columns_width = min((displaysize(stdout)[2]-18-3*N)÷N,100*N),
                 alignment = :l,
                 noheader = N == 1,
                 highlighters = highlighters; kwargs...)
end

current_package_state(m::Module) = current_package_state(Base.root_module_key(m))
function current_package_state(pkg::Base.PkgId)
    date = now()
    #println("Detected loading: ", pkg)
    m = Base.root_module(pkg)
    #println("Path: ", pathof(m))
    #println("pkgdir: ", pkgdir(m))
    d = pkgdir(m)
    th = tree_hash_fmt(d)
    project, mth = project_and_tree_hash(pkg)
    return PackageState(
            timestamp = date,
            dir = d,
            id = pkg,
            project = project,
            load_path = Base.load_path(),
            tree_hash = th,
            manifest_tree_hash = mth)
end

module_states = Dict{Module, Vector{PackageState}}()

function on_load_package(pkg::Base.PkgId)
    m = Base.root_module(pkg)
    push!(module_states, m => [current_package_state(pkg)])
end

function idxstates(m::Module, ::Val{i}) where {i}
    @assert(i isa Integer, "State index must be an integer or :on_load, :newest, :current")
    return module_states[m][i]
end
idxstates(m::Module, ::Val{:on_load}) = first(module_states[m])
idxstates(m::Module, ::Val{:newest}) = last(module_states[m])
idxstates(m::Module, ::Val{:current}) = current_package_state(m)

## Pretty printing

# one-line text output
Base.show(io::IO, s::PackageState) = print(io, "PackageState($(s.id))")

# fancy multi-line output
Base.show(io::IO, ::MIME"text/plain", s::PackageState) = printtable(s; title = "PackageState of $(s.id)")


text:
ff_promote_rule(::Type{T1}, ::Type{T2}) where {T1,T2} = promote_type(T1,T2)
ff_promote_rule(::Type{T}, ::Type{<:Integer}) where {T <: Integer} = T
ff_promote_rule(::Type{T}, ::Type{<:FloatingTypes}) where {T <: FloatingTypes} = T

Base.promote_rule(::Type{V}, ::Type{T2}) where {W,T1,T2<:NativeTypes,V<:AbstractSIMDVector{W,T1}} = Vec{W,ff_promote_rule(T1,T2)}
Base.promote_rule(::Type{V1}, ::Type{V2}) where {W,T1,T2,V1<:AbstractSIMDVector{W,T1},V2<:AbstractSIMDVector{W,T2}} = Vec{W,promote_type(T1,T2)}

Base.promote_rule(::Type{VecUnroll{N,W,T1}}, ::Type{T2}) where {N,W,T1,T2<:NativeTypes} = VecUnroll{N,W,ff_promote_rule(T1,T2)}
Base.promote_rule(::Type{VecUnroll{N,W,T1}}, ::Type{Vec{W,T2}}) where {N,W,T1,T2} = VecUnroll{N,W,promote_type(T1,T2)}
Base.promote_rule(::Type{VecUnroll{N,W,T1}}, ::Type{VecUnroll{N,W,T2}}) where {N,W,T1,T2} = VecUnroll{N,W,promote_type(T1,T2)}

Base.promote_rule(::Type{Bit}, ::Type{T}) where {T <: Number} = T

# Base.promote_rule(::Type{VecTile{M,N,W,T1}}, ::Type{T2}) where {M,N,W,T1,T2<:NativeTypes} = VecTile{M,N,W,promote_rule(T1,T2)}
# Base.promote_rule(::Type{VecTile{M,N,W,T1}}, ::Type{Vec{W,T2}}) where {M,N,W,T1,T2} = VecTile{M,N,W,promote_rule(T1,T2)}
# Base.promote_rule(::Type{VecTile{M,N,W,T1}}, ::Type{VecUnroll{M,W,T2}}) where {M,N,W,T1,T2} = VecTile{M,N,W,promote_rule(T1,T2)}
# Base.promote_rule(::Type{VecTile{M,N,W,T1}}, ::Type{VecTile{M,N,W,T2}}) where {M,N,W,T1,T2} = VecTile{M,N,W,promote_rule(T1,T2)}


text:
@testset "isvalidcode" begin
    @test GADM.isvalidcode("IND")
    @test GADM.isvalidcode("BRA")
    @test GADM.isvalidcode("USA")
    @test !GADM.isvalidcode("ind")
    @test !GADM.isvalidcode("")
    @test !GADM.isvalidcode("iNd4")
end

@testset "download" begin
    # test successful download
    GADM.download("VAT")
    @test_nowarn @datadep_str "GADM_VAT"
end

@testset "dataread" begin
    path = GADM.download("VAT")
    data = GADM.dataread(path)
    @test typeof(data) === ArchGDAL.IDataset
end

@testset "getlayer" begin
    path = GADM.download("IND")
    data = GADM.dataread(path)

    # incorrect level
    @test_throws ArgumentError GADM.getlayer(data, 10)

    # correct level
    layer = GADM.getlayer(data, 0)
    @test typeof(layer) === ArchGDAL.IFeatureLayer
end

@testset "get" begin
    # invalid country codes
    @test_throws ArgumentError GADM.get("ind")
    @test_throws ArgumentError GADM.get("")
    @test_throws ArgumentError GADM.get("IND4")
    # valid country codes, children=false
    parent = GADM.get("IND")
    @test parent isa NamedTuple
    @test GeoInterface.geotype(parent.geom[1]) == :MultiPolygon
    # valid country code, children=true
    parent, children = GADM.get("IND";children=true)
    @test parent isa NamedTuple
    @test children isa NamedTuple
    @test GeoInterface.geotype(parent.geom[1]) == :MultiPolygon
    @test length(children) == 11 #number of fields in named tuple
    geometries = Tables.getcolumn(children, Symbol("geom"))
    @test geometries isa Array{ArchGDAL.IGeometry,1}
    @test length(geometries) == 36 # number of rows
    # throws error when query is invalid
    @test_throws ArgumentError GADM.get("IND", "Rio De Janerio")
    # throws argument error for supplying deeper region than available in dataset
    @test_throws ArgumentError GADM.get("VAT", "Pope")
end

@testset "basic" begin
    polygon = GADM.get("VAT")
    bounds = ArchGDAL.envelope(polygon.geom[1])
    bounds_arr = [bounds.MinX, bounds.MaxX, bounds.MinY, bounds.MaxY]

    # Vatican City bounding boxes lat min 41.9002044 lat max 41.9073912 lon min 12.4457442 lon max 12.4583653
    bounds_actual = [12.4457442, 12.4583653, 41.9002044, 41.9073912]

    for i=1:4
        diff = abs(bounds_arr[i] - bounds_actual[i])
        @test diff < 0.01
    end
end

@testset "coordinates" begin
    # invalid country code: lowercase
    @test_throws ArgumentError GADM.coordinates("ind")
    # empty country code
    @test_throws ArgumentError GADM.coordinates("")
    # invalid code other than format [A-Z]{3}
    @test_throws ArgumentError GADM.coordinates("IND4")
    # Polygon
    c = GADM.coordinates("VAT")
    @test c isa Array{Array{Array{Float64,1},1},1}
    # MultiPolygon
    c = GADM.coordinates("IND", "Gujarat")
    @test c isa Array{Array{Array{Array{Float64,1},1},1},1}
end


text:
#
# test hierarchicalsolve
#
using ApproxFun, SingularIntegralEquations, Base.Test

M = 50
L = Array(LowRankMatrix{Float64},14)
[L[i] = LowRankMatrix(0.001rand(4M,2),0.001rand(4M,2)) for i=1:2]
[L[i] = LowRankMatrix(0.01rand(2M,2),0.01rand(2M,2)) for i=[3:4;9:10]]
[L[i] = LowRankMatrix(0.1rand(M,2),0.1rand(M,2)) for i=[5:8;11:14]]
D = Array(Diagonal{Float64},8)
[D[i] = Diagonal(1./collect(1:M)) for i=1:8]

H = HierarchicalMatrix(D,L)

@test rank(H) == 400
#@test blockrank(H) == fill(2,8,8)+diagm(fill(48,8))
@test isfactored(H) == false

B = Array(Vector{Float64},8)
[B[i] = rand(M) for i=1:8]
b = rand(size(H,1))#HierarchicalVector(B)

full(H)\b
@time x = full(H)\b
H\b # includes compile time
H = HierarchicalMatrix(D,L)
@time xw = H\b # includes precomputation
@time xw = H\b # H.A's precomputed and factored :)

@test isfactored(H) == true

CH = cond(H)

@test norm(x-xw) ≤ 10CH^2*eps()
@test norm(H*x-b) ≤ 10CH*eps()
@test norm(H*xw-b) ≤ 10CH*eps()

H+H-(H-H)
@test norm(full(H-H)) < 10norm(full(H))*eps()

# Hierarchical operators case

i = 2

dom = cantor(Interval(),i)
⨍ = DefiniteLineIntegral(dom)
f = Fun(x->logabs(x-5im),dom)
sp = Space(dom)

G = GreensFun((x,y)->1/2,CauchyWeight(sp⊗sp,0);method=:Cholesky)

@time u1 = ⨍[G]\transpose(f)
println("Adaptive QR  forward error norm is: ",norm(⨍[G]*u1-f))

@test norm(⨍[G]*u1-f) < 10eps()

hdom = clustertree(dom)
hsp = Space(hdom)

G1 = GreensFun((x,y)->1/2,CauchyWeight(hsp⊗hsp,0);method=:Cholesky)
H = ⨍[G1]

@time u2 = H\f
@time u2 = H\f

@test condest(H) ≤ 20.0

import ApproxFun: ∞
@test blockrank(H) == Number[∞ 8 9 9; 8 ∞ 9 9; 9 9 ∞ 8; 9 9 8 ∞]

println("The hierarchical forward error norm is: ",norm(⨍[G]*u2-f))

@test norm(⨍[G]*u2-f) < 10eps()

# Test three domains

dom = Interval(-2.0,-1.0+0im)∪Interval(1.0,2.0+0im)∪Interval(-0.5-2im,0.5-2im)
hdom = clustertree(dom)
dom = UnionDomain(SingularIntegralEquations.collectdata(hdom))
⨍ = DefiniteLineIntegral(dom)
f = Fun(x->logabs(x-5im),dom)
sp = Space(dom)

G = GreensFun((x,y)->1/2,CauchyWeight(sp⊗sp,0);method=:Cholesky)

@time u1 = ⨍[G]\transpose(f)

println("Adaptive QR  forward error norm is: ",norm(⨍[G]*u1-f))

@test norm(⨍[G]*u1-f) < 300eps()

if VERSION ≥ v"0.5pre"
    warn("Skipping test in 0.5 due to error in type inference due to #265")
else
    hsp = Space(hdom)

    G1 = GreensFun((x,y)->1/2,CauchyWeight(hsp⊗hsp,0);method=:Cholesky)
    H = ⨍[G1]

    @time u2 = H\f
    @time u2 = H\f

    println("The hierarchical forward error norm is: ",norm(⨍[G]*u2-f))

    @test norm(⨍[G]*u2-f) < 10000eps()
end


text:
# Load the output from a previously build
include("build_output.jl")

include("builder_defs.jl")
include("builder_tools.jl")

build_jlls(ARGS, name, output, dependencies)

max_stars_repo_name: pabloferz/LAMMPSBuilder

max_issues_repo_name: pabloferz/LAMMPSBuilder

max_forks_repo_name: pabloferz/LAMMPSBuilder


text:
function walk_dom(f, x::JSCode, visited = IdDict())
    walk_dom(f, x.source, visited)
end

walk_dom(f, x, visited = IdDict()) = f(x)

walk_dom(f, x::Markdown.MD, visited = IdDict()) = walk_dom(f, x.content, visited)
walk_dom(f, x::Markdown.Header, visited = IdDict()) = walk_dom(f, x.text, visited)
walk_dom(f, x::Markdown.Paragraph, visited = IdDict()) = walk_dom(f, x.content, visited)

function walk_dom(f, x::Union{Tuple, AbstractVector, Pair}, visited = IdDict())
    get!(visited, x, nothing) !== nothing && return
    for elem in x
        walk_dom(f, elem, visited)
    end
end

function walk_dom(f, x::Node, visited = IdDict())
    get!(visited, x, nothing) !== nothing && return
    for elem in children(x)
        walk_dom(f, elem, visited)
    end
    for (name, elem) in Hyperscript.attrs(x)
        walk_dom(f, elem, visited)
    end
end

const mime_order = MIME.((
    "text/html", "text/latex", "image/svg+xml", "image/png",
    "image/jpeg", "text/markdown", "application/javascript", "text/plain"
))

function richest_mime(val)
    for mimetype in mime_order
        showable(mimetype, val) && return mimetype
    end
    error("value not writable for any mimetypes")
end

repr_richest(x) = repr(richest_mime(x), x)
repr_richest(x::String) = x
repr_richest(x::Number) = sprint(print, x)

columns(args...; class="") = DOM.div(args..., class=class * " flex flex-col")
rows(args...; class="") = DOM.div(args..., class=class * " flex flex-row")

function grid(args...; cols=4, rows=4, class="")
    class *= " grid auto-cols-max grid-cols-$(cols) grid-rows-$(rows) gap-4"
    return DOM.div(args..., class=class)
end

function styled_slider(slider, value; class="")
    return rows(slider,
                DOM.span(value, class="p-1");
                class="w-64 p-2 items-center " * class)
end


text:
using Revise

using RCall
using FastGroupBy, BenchmarkTools, DataFrames

const N = Int(1e8); const K = 100;

R"""
memory.limit(2e9)
library(data.table)
N=$N; K=$K
set.seed(1)
DT <- data.table(
  id5 = sample(K, N, TRUE),                          # large groups (int)
  id6 = sample(N/K, N, TRUE),                        # small groups (int)
  v3 =  sample(round(runif(100,max=100),4), N, TRUE) # numeric e.g. 23.5749
)
cat("GB =", round(sum(gc()[,2])/1024, 3), "\n")
pt = proc.time()
DT[, sum(v3), keyby=id5]
x = timetaken(pt); pt = proc.time()
DT[, sum(v3), keyby=id6]
y = timetaken(pt)
rm(DT); gc()
"""

@rget x;
@rget y;

srand(1);
df = DataFrame(id = rand(Int32(1):Int32(round(N/K)), N), id_small = rand(Int8(1):Int8(K),N), val = rand(round.(rand(K)*100,4), N));
x1 = @elapsed sumby!(df, :id_small, :val);

df = nothing; gc()
srand(1);
df = DataFrame(id = rand(Int32(1):Int32(round(N/K)), N), id_small = rand(Int8(1):Int8(K),N), val = rand(round.(rand(K)*100,4), N));
y1 = @elapsed sumby!(df, :id, :val);

df = nothing; gc();
srand(1);
id = rand(Int32(1):Int32(round(N/K)), N);
val = rand(round.(rand(K)*100,4), N);
y2 = @elapsed fastby!(sum, id, val);

y
y1
y2

x1/parse(Float64, x[1:end-3])
y1/parse(Float64, y[1:end-3])
y2/parse(Float64, y[1:end-3])


text:
using Pkg;
Pkg.activate("/home/users/nticea/ITensors.jl/Project.toml")
using ITensors
using JLD

## FUNCTIONS ##

struct Parameters
    # Model parameters
    N::Int
    t::Real
    U::Real
    ω::Real
    g0::Real
    g1::Real
    doping::Real

    # DMRG parameters
    sweeps::Sweeps

    # TEBD parameters
    mid::Int
    T::Int
    τ::Real
    TEBD_cutoff
    TEBD_maxdim
end

struct HubbardHolsteinModel
    # Lattice
    sites
    # Hamiltonian (in two forms)
    mpo::MPO
    gates
end

struct DMRGResults
    ground_state
    ground_state_energy
    entropy
end

struct TEBDResults
    entropy
    self_overlap
    phonon_flux
    corrs
end

function parameters(;N::Int, t::Real, U::Real=nothing, ω::Real=nothing,
                    g0::Real=nothing, g1::Real=nothing, doping::Real=0,
                    DMRG_numsweeps::Int=20, DMRG_noise=nothing, DMRG_maxdim=nothing, DMRG_cutoff=nothing,
                    T::Int=25, τ::Real=0.1, TEBD_cutoff=1E-14, TEBD_maxdim=400)

    if isnothing(U)
        U = 8*t
    end
    if isnothing(ω)
        ω = 0.5*t
    end
    if isnothing(g0)
        g0 = 0.05*t
    end
    if isnothing(g1)
        g1 = 0*g0
    end

    sweeps = Sweeps(DMRG_numsweeps)
    if isnothing(DMRG_noise)
        DMRG_noise = [1E-6,1E-6,1E-8,0]
    end
    if isnothing(DMRG_maxdim)
        DMRG_maxdim = [20,40,100,200,400]
    end

    if isnothing(DMRG_cutoff)
        DMRG_cutoff = 1E-10
    end
    setnoise!(sweeps, DMRG_noise...) # Very important to use noise for this model
    setmaxdim!(sweeps, DMRG_maxdim...)
    setcutoff!(sweeps, DMRG_cutoff...)

    Parameters(N,t,U,ω,g0,g1,doping,sweeps,ceil(Int,N/2),T,τ,TEBD_cutoff,TEBD_maxdim)
end

function HubbardHolsteinModel(p::Parameters)
    N, t, U, ω, g0, g1, τ = p.N, p.t, p.U, p.ω, p.g0, p.g1, p.τ

    # make the sites
    sites = siteinds("HubHolst", N)

    # make the hamiltonian
    ampo = OpSum()
    for j=1:N-1
        # ∑_j,σ t * (c^†_jσ c_{j+1}σ + h.c.)
        ampo += -t,"Cdagup",j,"Cup",j+1
        ampo += -t,"Cdagup",j+1,"Cup",j
        ampo += -t,"Cdagdn",j,"Cdn",j+1
        ampo += -t,"Cdagdn",j+1,"Cdn",j

        # ∑_j U * n_j↑ n_j↓
        ampo += U,"Nupdn",j,"I",j

        # ∑_j ω * nb_j
        ampo += ω,"Nb",j

        # # ∑_j g0 * nf_j (b^†_j + b_j)
        ampo += g0,"Ntot(Bd+B)",j

        # # ∑_⟨ij⟩ g1 * nf_j (b^†_i + b_i)
        ampo += g1,"Ntot",j,"Bdag+B",j+1
        ampo += g1,"Ntot",j+1,"Bdag+B",j
    end
    # Edge site
    ampo += U,"Nupdn",N
    ampo += ω,"Nb",N
    ampo += g0,"Ntot(Bd+B)",N
    H = MPO(ampo,sites)

    # make the trotter gates
    gates = ITensor[]
    for j=1:N-1
        s1 = sites[j] # site j
        s2 = sites[j+1] # site j+1

        hj_twosite = -t*(op("Cdagup*F",s1) * op("Cup",s2)  # t * (c^†_jσ c_{j+1}σ + h.c.)
                 -op("Cup*F",s1) * op("Cdagup",s2)
                 +op("Cdagdn*F",s1) * op("Cdn",s2)
                 -op("Cdn*F",s1) * op("Cdagdn",s2))
                + g1*(op("Ntot",s1) * op("Bdag+B",s2))
                + g1*(op("Bdag+B",s1) * op("Ntot",s2))

        hj_onesite = U*(op("Nupdn",s1) * op("I",s2))
                    + ω*(op("Nb",s1) * op("I",s2))
                    + g0*(op("Ntot(Bd+B)",s1) * op("I",s2))

        Gj_twosite = exp(-1.0im * τ/2 * hj_twosite)
        Gj_onesite = exp(-1.0im * τ/2 * hj_onesite)
        push!(gates,Gj_twosite)
        push!(gates,Gj_onesite)
    end
    # End site
    hn = U*op("Nupdn",sites[N])
        + ω*op("Nb",sites[N])
        + g0*op("Ntot(Bd+B)",sites[N])
    Gn = exp(-1.0im * τ/2 * hn)
    push!(gates,Gn)
    append!(gates,reverse(gates))

    HubbardHolsteinModel(sites, H, gates)
end

function initialize_wavefcn(HH::HubbardHolsteinModel, p::Parameters)
    state = [isodd(n) ? "Up,0" : "Dn,0" for n=1:p.N] # NOTE: the QN of this state is preserved through DMRG

    # Account for doping
    if p.doping > 0
        spacing = floor(1/p.doping)
        state[1:spacing:end] .= "Emp,0"
    end

    productMPS(HH.sites,state)
end

function run_DMRG(HH::HubbardHolsteinModel, p::Parameters; alg="divide_and_conquer")
    ϕ0 = initialize_wavefcn(HH,p)
    @show flux(ϕ0)
    energy, ϕ = dmrg(HH.mpo, ϕ0, p.sweeps, alg=alg)
    entropy = compute_entropy(ϕ, p.mid)
    return DMRGResults(ϕ, energy, entropy)
end

function apply_twosite_operator(ϕ::MPS, opname_1::String, opname_2::String, sites, siteidx::Int)
    ϕ = copy(ϕ)

end

"""
    function apply_onesite_operator(ϕ::MPS, opname::String, sites, siteidx::Int)
This is a wrapper function for the apply_op function, which actually applies the operator
Here, we just account for fermionic strings
"""
function apply_onesite_operator(ϕ::MPS, opname::String, sites, siteidx::Int)
    ϕ = copy(ϕ)

    ## Account for fermion sign using Jordan-Wigner strings ##
    if opname == "Cup" || opname == "Cdn"
        ϕ = apply_op(ϕ, op(opname,sites[siteidx]), siteidx)
        for i in reverse(1:(siteidx-1)) # Don't act with string on-site
            ϕ = apply_op(ϕ, op("F",sites[i]), i)
        end
        return ϕ
    elseif opname == "Cdagup" || opname == "Cdagdn"
        for i in 1:(siteidx-1) # Don't act with string on-site
            ϕ = apply_op(ϕ, op("F",sites[i]), i)
        end
        ϕ = apply_op(ϕ, op(opname,sites[siteidx]), siteidx)
        return ϕ
    end
    println("No Jordan-Wigner strings")

    # Otherwise, just apply the operator as usual
    return apply_op(ϕ, op(opname,sites[siteidx]), siteidx)
end

function apply_op(ϕ::MPS, op::ITensor, siteidx::Int)
    ϕ = copy(ϕ) # Make a copy of the original state
    orthogonalize!(ϕ, siteidx)
    new_ϕj = op * ϕ[siteidx] # Apply the local operator
    noprime!(new_ϕj)
    ϕ[siteidx] = new_ϕj
    return ϕ
end

function compute_entropy(ψ::MPS, b::Int)
    orthogonalize!(ψ, b)
    U,S,V = svd(ψ[b], (linkind(ψ, b-1), siteind(ψ,b)))
    SvN = 0.0
    for n=1:dim(S, 1)
        p = S[n,n]^2
        SvN -= p * log(p)
    end
    return SvN
end

function compute_overlap(ψ1::MPS, ψ2::MPS)
    LinearAlgebra.norm(inner(ψ1, ψ2))
end

function compute_phonon_number(ψ::MPS)
    expect(ψ,"Nb")
end

function equilibrium_correlations(dmrg_results::DMRGResults, corrtype::String, HH::HubbardHolsteinModel, p::Parameters)
    ϕ = copy(dmrg_results.ground_state)
    j = p.mid
    sites = HH.sites
    if corrtype=="spin"
        return correlation_matrix(ϕ, "Sz", "Sz")[:,j]
    elseif corrtype=="charge"
        ninj = correlation_matrix(ϕ, "Ntot", "Ntot")[:,j]
        ni = expect(ϕ, "Ntot")
        nj = ni[j]
        return ninj - nj .* ni
    elseif corrtype=="sSC"
        Σ_j = op("Cupdn", sites[j])
        ψ = apply_op(ϕ, Σ_j, j)
        function compute_corr(i::Int)
            Σ_i = op("Cdagupdn", sites[i])
            Σ_iψ = apply_op(ψ, Σ_i, i)
            return inner(ϕ,Σ_iψ)
        end
        return compute_corr.(collect(1:p.N))
    elseif corrtype=="pSC"
        Π_j = 1/sqrt(2)*(op("Cup",sites[j])*op("Cdn",sites[j+1]) + op("Cdn",sites[j])*op("Cup",sites[j+1]))
        ψ = apply_op(ϕ, Π_j, j)
        function compute_corr(i::Int)
            Π_i = 1/sqrt(2)*(op("Cdagdn",sites[i])*op("Cdagup",sites[i+1]) + op("Cdagup",sites[i])*op("Cdagdn",sites[i+1]))
            Π_iψ = apply_op(ψ, Π_i, i)
            return inner(ϕ,Π_iψ)
        end
        return compute_corr.(collect(1:p.N))
    elseif corrtype=="dSC"
        Δ_j = 1/sqrt(2)*(op("Cup",sites[j])*op("Cdn",sites[j+1]) - op("Cdn",sites[j])*op("Cup",sites[j+1]))
        ψ = apply_op(ϕ, Δ_j, j)
        function compute_corr(i::Int)
            Δ_i = 1/sqrt(2)*(op("Cdagdn",sites[i])*op("Cdagup",sites[i+1]) - op("Cdagup",sites[i])*op("Cdagdn",sites[i+1]))
            Δ_iψ = apply_op(ψ, Δ_i, i)
            return inner(ϕ,Δ_iψ)
        end
        return compute_corr.(collect(1:p.N))
    else
        @error "Unknown spin correlation type"
    end
end

function compute_correlations(dmrg_results::DMRGResults, A_t0::String, A_t::String, HH::HubbardHolsteinModel, p::Parameters)
    # Results
    corrs = []
    entropy = []
    self_overlap = []
    phonon_flux = []

    # The wavefunction being acted upon at t=0, |ψ⟩ = A_t0|ϕ⟩
    ϕ = copy(dmrg_results.ground_state)
    ψ = copy(ϕ)

    # Apply A_t0 to middle site
    ψ = apply_onesite_operator(ψ, A_t0, HH.sites, p.mid)

    # Measure correlation before time evolution
    A_tψ = apply_onesite_operator(ψ, A_t, HH.sites, p.mid)
    # Measure the correlation fcn
    @show inner(ϕ, A_tψ)

    # Parameters for time evolution
    nsteps = floor(p.T/p.τ) # Number of time steps for time evolution
    t = 0.0
    for step in 1:nsteps
        print(floor(Int,step),"-")
        ϕ = apply(HH.gates, ϕ; maxdim=p.TEBD_maxdim, cutoff=p.TEBD_cutoff)
        ψ = apply(HH.gates, ψ; maxdim=p.TEBD_maxdim, cutoff=p.TEBD_cutoff) # evolve forward

        t += p.τ

        ### SANITY CHECKS
        if step%1==0
            # Compute entropy
            push!(entropy, (compute_entropy(ϕ, p.mid),compute_entropy(ψ, p.mid)))
            println("Entropy of ϕ, ψ: ", entropy[end])

            # Compute overlap of ϕ with its original self
            push!(self_overlap, compute_overlap(ϕ,dmrg_results.ground_state))
            println("Overlap of ϕ with itself: ", self_overlap[end])

            # Phonon flux
            push!(phonon_flux, (compute_phonon_number(ϕ), compute_phonon_number(ψ)))
            println("Phonon flux ϕ: ", sum(phonon_flux[end][1])/p.N)
            println("Phonon flux ψ: ", sum(phonon_flux[end][2])/p.N)
        end

        # Calculate ⟨ϕ(t)|c_j^† c_i|ϕ(0)⟩
        function measure_corr(j::Int)
            # Apply the second measurement operator
            A_tψ = apply_onesite_operator(ψ, A_t, HH.sites, j)
            return dot(ϕ,A_tψ)
        end

        # Measure the correlation fcn
        push!(corrs,measure_corr.(collect(1:p.N)))
    end
    return TEBDResults(entropy, self_overlap, phonon_flux, hcat(corrs...))
end

function save_structs(struc, path::String)
    function Name(arg)
        string(arg)
    end
    fnames = fieldnames(typeof(struc))
    for fn in fnames
        n = Name(fn)
        d = getfield(struc, fn)
        jldopen(path, "w") do file
            write(file, n, d)
        end
    end
end

## SCRIPT ##

## PARAMETERS ##

# Model
N = 80
t = 1 ## THIS TERM IS FINE
U = 8
ω = 0*t ## THIS TERM IS FINE (by itself)
g0 = 0*t ## THIS TERM IS FINE (by itself)
g1 = 0*g0 ## THIS TERM IS FINE (by itself)

# Simulation
T = 20
τ = 0.05
DMRG_numsweeps = 40
DMRG_maxdim = 600
TEBD_maxdim = 800
TEBD_cutoff = 1E-14
DMRG_cutoff = 1E-14

# Saveout info
DO_SAVE = true
save_path = "/home/users/nticea/ITensors.jl/examples/dmrg/results/1dhubbholst_output.jld"

## TODO: Modify number of phonons on each site from script

# Specify operators of interest
A_t0 = "Cup"
A_t = "Cdagup"

# Initialize
println("Initializing...")
params = parameters(N=N, t=t, U=U, ω=ω, g0=g0, g1=g1,
                    DMRG_numsweeps=DMRG_numsweeps,
                    DMRG_maxdim=DMRG_maxdim, DMRG_cutoff=DMRG_cutoff,
                    T=T, τ=τ, TEBD_cutoff=TEBD_cutoff)
hubbholst = HubbardHolsteinModel(params)
if DO_SAVE
    save_structs(params, save_path)
end

# Run DMRG
println("Finding ground state...")
dmrg_results = run_DMRG(hubbholst, params, alg="divide_and_conquer")
if DO_SAVE
    save_structs(dmrg_results, save_path)
end

# Compute correlation functions
println("Computing correlation functions...")
tebd_results = compute_correlations(dmrg_results, A_t0, A_t, hubbholst, params)
if DO_SAVE
    save_structs(tebd_results, save_path)
end


text:
@inline eval_pow(::MinkowskiMetric, s) = abs(s)
@inline eval_pow(::Euclidean, s) = abs2(s)
@inline eval_pow(d::Minkowski, s) = abs(s)^d.p

@inline eval_diff(::MinkowskiMetric, a, b) = a - b
@inline eval_diff(::Chebyshev, ::Any, b) = b

function Distances.evaluate(d::Distances.UnionMetrics, a::AbstractVector,
                            b::AbstractVector, do_end::Bool)
    s = eval_start(d, a, b)
    @simd for i in eachindex(b)
        @inbounds ai = a[i]
        @inbounds bi = b[i]
        s = eval_reduce(d, s, eval_op(d, ai, bi))
    end
    if do_end
        return eval_end(d, s)
    else
        return s
    end
end

function Distances.evaluate(d::Distances.PreMetric, a::AbstractVector,
                            b::AbstractVector, ::Bool)
    evaluate(d, a, b)
end


text:
push!(LOAD_PATH, pwd()*"/../../src/")
push!(LOAD_PATH, pwd()*"/src/")
using Distributions
using PyPlot
using BVB

lower = -20.
upper = 20.

true_d = Cauchy(0.0, 2)
ll(x) = logpdf(true_d, x)

init_approx_d = Normal(0, 40)
x = linspace(-20,20,500)

max_iter = 30

# include("run_laplacian_gradient_boosting.jl")
include("run_laplacian_gradient_boosting_exclusive.jl")
save_trace("./cauchy/", x, init_approx_d, approx_d, true_d,
                    iter_vec, α_vec,
                    eval_vec, -elbo_vec)
show()

println("Done.")

max_stars_repo_name: richardkwo/boosting-VB

max_issues_repo_name: richardkwo/boosting-VB

max_forks_repo_name: richardkwo/boosting-VB


text:
module SpringCollabLightningDashSwitch

using ..Ahorn, Maple

sides = String[
    "left",
    "right",
    "up",
    "down"
]

textures = String[
    "default",
    "mirror"
]

@mapdef Entity "SpringCollab2020/LightningDashSwitch" LDashSwitch(x::Integer, y::Integer, side::String="up", persistent::Bool=false, sprite::String="default")

const placements = Ahorn.PlacementDict(
    "Lightning Dash Switch ($(uppercasefirst(side))) (Spring Collab 2020)" => Ahorn.EntityPlacement(
        LDashSwitch,
        "rectangle",
        Dict{String, Any}(
            "side" => side
        )
    ) for side in sides
)

Ahorn.editingOptions(entity::LDashSwitch) = Dict{String, Any}(
    "side" => sides
)

function Ahorn.selection(entity::LDashSwitch)
    x, y = Ahorn.position(entity)
    side = get(entity.data, "side", false)

    if side == "left"
        return Ahorn.Rectangle(x, y - 1, 10, 16)
    elseif side == "right"
        return Ahorn.Rectangle(x - 2, y, 10, 16)
    elseif side == "down"
        return Ahorn.Rectangle(x, y, 16, 12)
    elseif side == "up"
        return Ahorn.Rectangle(x, y - 4, 16, 12)
    end
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::LDashSwitch, room::Maple.Room)
    sprite = get(entity.data, "sprite", "default")
    side = get(entity.data, "side", "up")
    texture = sprite == "default" ? "objects/temple/dashButton00.png" : "objects/temple/dashButtonMirror00.png"

    if side == "left"
        Ahorn.drawSprite(ctx, texture, 20, 25, rot=pi)
    elseif side == "right"
        Ahorn.drawSprite(ctx, texture, 8, 8)
    elseif side == "up"
        Ahorn.drawSprite(ctx, texture, 9, 20, rot=-pi / 2)
    elseif side == "down"
        Ahorn.drawSprite(ctx, texture, 27, 7, rot=pi / 2)
    end
end

end


text:
import DataArrays: NAtype
import JSON
import Compat: Libc, unsafe_convert, parse, @compat

abstract AbstractPostgresType
type PostgresType{Name} <: AbstractPostgresType end

abstract AbstractOID
type OID{N} <: AbstractOID end

oid{T<:AbstractPostgresType}(t::Type{T}) = convert(OID, t)
pgtype(t::Type) = convert(PostgresType, t)

Base.convert{T}(::Type{Oid}, ::Type{OID{T}}) = convert(Oid, T)

function newpgtype(pgtypename, oid, jltypes)
    Base.convert(::Type{OID}, ::Type{PostgresType{pgtypename}}) = OID{oid}
    Base.convert(::Type{PostgresType}, ::Type{OID{oid}}) = PostgresType{pgtypename}

    for t in jltypes
        Base.convert(::Type{PostgresType}, ::Type{t}) = PostgresType{pgtypename}
    end
end

# To check oids in postgres, use this SQL query on a psql prompt:
#   select typname, typarray, oid from pg_type;
##
newpgtype(:bool, 16, (Bool,))
newpgtype(:bytea, 17, (Vector{UInt8},))
newpgtype(:int8, 20, (Int64,))
newpgtype(:int4, 23, (Int32,))
newpgtype(:int2, 21, (Int16,))
newpgtype(:float8, 701, (Float64,))
newpgtype(:float4, 700, (Float32,))
newpgtype(:bpchar, 1042, ())
newpgtype(:varchar, 1043, (ASCIIString,UTF8String))
newpgtype(:text, 25, ())
newpgtype(:numeric, 1700, (BigInt,BigFloat))
newpgtype(:date, 1082, ())
newpgtype(:timestamp, 1114, ())
newpgtype(:timestamptz, 1184, ())
newpgtype(:unknown, 705, (Union,NAtype))
newpgtype(:json, 114, (Dict{AbstractString,Any},))
newpgtype(:jsonb, 3802, (Dict{AbstractString,Any},))

# Support for Postgres array types, underscore indicates array
newpgtype(:_bool, 1000, (Vector{Bool},))
newpgtype(:_int8, 1016, (Vector{Int64},))
newpgtype(:_int4, 1007, (Vector{Int32},))
newpgtype(:_int2, 1005, (Vector{Int16},))
newpgtype(:_float8, 1022, (Vector{Float64},))
newpgtype(:_float4, 1021, (Vector{Float32},))
newpgtype(:_varchar, 1015, (Vector{ASCIIString}, Vector{UTF8String}))
newpgtype(:_text, 1009, (Vector{ASCIIString}, Vector{UTF8String}))

typealias PGStringTypes Union{Type{PostgresType{:bpchar}},
                              Type{PostgresType{:varchar}},
                              Type{PostgresType{:text}},
                              Type{PostgresType{:date}}}

function storestring!(ptr::Ptr{UInt8}, str::AbstractString)
    ptr = convert(Ptr{UInt8}, Libc.realloc(ptr, sizeof(str)+1))
    unsafe_copy!(ptr, unsafe_convert(Ptr{UInt8}, str), sizeof(str)+1)
    return ptr
end

# In text mode, pq returns bytea as a text string "\xAABBCCDD...", for instance
# UInt8[0x01, 0x23, 0x45] would come in as "\x012345"
function decode_bytea_hex(s::AbstractString)
    if length(s) < 2 || s[1] != '\\' || s[2] != 'x'
        error("Malformed bytea string: $s")
    end
    return hex2bytes(s[3:end])
end

jldata(::Type{PostgresType{:date}}, ptr::Ptr{UInt8}) = bytestring(ptr)

jldata(::Type{PostgresType{:timestamp}}, ptr::Ptr{UInt8}) = bytestring(ptr)

jldata(::Type{PostgresType{:timestamptz}}, ptr::Ptr{UInt8}) = bytestring(ptr)

jldata(::Type{PostgresType{:bool}}, ptr::Ptr{UInt8}) = bytestring(ptr) != "f"

jldata(::Type{PostgresType{:int8}}, ptr::Ptr{UInt8}) = parse(Int64, bytestring(ptr))

jldata(::Type{PostgresType{:int4}}, ptr::Ptr{UInt8}) = parse(Int32, bytestring(ptr))

jldata(::Type{PostgresType{:int2}}, ptr::Ptr{UInt8}) = parse(Int16, bytestring(ptr))

jldata(::Type{PostgresType{:float8}}, ptr::Ptr{UInt8}) = parse(Float64, bytestring(ptr))

jldata(::Type{PostgresType{:float4}}, ptr::Ptr{UInt8}) = parse(Float32, bytestring(ptr))

function jldata(::Type{PostgresType{:numeric}}, ptr::Ptr{UInt8})
    s = bytestring(ptr)
    return parse(search(s, '.') == 0 ? BigInt : BigFloat, s)
end

jldata(::PGStringTypes, ptr::Ptr{UInt8}) = bytestring(ptr)

jldata(::Type{PostgresType{:bytea}}, ptr::Ptr{UInt8}) = bytestring(ptr) |> decode_bytea_hex

jldata(::Type{PostgresType{:unknown}}, ptr::Ptr{UInt8}) = Union{}

jldata(::Type{PostgresType{:json}}, ptr::Ptr{UInt8}) = JSON.parse(bytestring(ptr))

jldata(::Type{PostgresType{:jsonb}}, ptr::Ptr{UInt8}) = JSON.parse(bytestring(ptr))

jldata(::Type{PostgresType{:_bool}}, ptr::Ptr{UInt8}) = map(x -> x != "f", split(bytestring(ptr)[2:end-1], ','))

jldata(::Type{PostgresType{:_int8}}, ptr::Ptr{UInt8}) = map(x -> parse(Int64, x), split(bytestring(ptr)[2:end-1], ','))

jldata(::Type{PostgresType{:_int4}}, ptr::Ptr{UInt8}) = map(x -> parse(Int32, x), split(bytestring(ptr)[2:end-1], ','))

jldata(::Type{PostgresType{:_int2}}, ptr::Ptr{UInt8}) = map(x -> parse(Int16, x), split(bytestring(ptr)[2:end-1], ','))

jldata(::Type{PostgresType{:_float8}}, ptr::Ptr{UInt8}) = map(x -> parse(Float64, x), split(bytestring(ptr)[2:end-1], ','))

jldata(::Type{PostgresType{:_float4}}, ptr::Ptr{UInt8}) = map(x -> parse(Float32, x), split(bytestring(ptr)[2:end-1], ','))

jldata(::Type{PostgresType{:_varchar}}, ptr::Ptr{UInt8}) = convert(Vector{AbstractString}, split(bytestring(ptr)[2:end-1], ','))

jldata(::Type{PostgresType{:_text}}, ptr::Ptr{UInt8}) = convert(Vector{AbstractString}, split(bytestring(ptr)[2:end-1], ','))

function pgdata(::Type{PostgresType{:bool}}, ptr::Ptr{UInt8}, data::Bool)
    ptr = data ? storestring!(ptr, "TRUE") : storestring!(ptr, "FALSE")
end

function pgdata(::Type{PostgresType{:int8}}, ptr::Ptr{UInt8}, data::Number)
    ptr = storestring!(ptr, string(convert(Int64, data)))
end

function pgdata(::Type{PostgresType{:int4}}, ptr::Ptr{UInt8}, data::Number)
    ptr = storestring!(ptr, string(convert(Int32, data)))
end

function pgdata(::Type{PostgresType{:int2}}, ptr::Ptr{UInt8}, data::Number)
    ptr = storestring!(ptr, string(convert(Int16, data)))
end

function pgdata(::Type{PostgresType{:float8}}, ptr::Ptr{UInt8}, data::Number)
    ptr = storestring!(ptr, string(convert(Float64, data)))
end

function pgdata(::Type{PostgresType{:float4}}, ptr::Ptr{UInt8}, data::Number)
    ptr = storestring!(ptr, string(convert(Float32, data)))
end

function pgdata(::Type{PostgresType{:numeric}}, ptr::Ptr{UInt8}, data::Number)
    ptr = storestring!(ptr, string(data))
end

function pgdata(::PGStringTypes, ptr::Ptr{UInt8}, data::ByteString)
    ptr = storestring!(ptr, data)
end

function pgdata(::PGStringTypes, ptr::Ptr{UInt8}, data::AbstractString)
    ptr = storestring!(ptr, bytestring(data))
end

function pgdata(::PostgresType{:date}, ptr::Ptr{UInt8}, data::AbstractString)
    ptr = storestring!(ptr, bytestring(data))
    ptr = Dates.DateFormat(ptr)
end

function pgdata(::PostgresType{:timestamp}, ptr::Ptr{UInt8}, data::AbstractString)
    ptr = storestring!(ptr, bytestring(data))
end

function pgdata(::PostgresType{:timestamptz}, ptr::Ptr{UInt8}, data::AbstractString)
    ptr = storestring!(ptr, bytestring(data))
end

function pgdata(::Type{PostgresType{:bytea}}, ptr::Ptr{UInt8}, data::Vector{UInt8})
    ptr = storestring!(ptr, bytestring("\\x", bytes2hex(data)))
end

function pgdata(::Type{PostgresType{:unknown}}, ptr::Ptr{UInt8}, data)
    ptr = storestring!(ptr, string(data))
end

function pgdata{T<:AbstractString}(::Type{PostgresType{:json}}, ptr::Ptr{UInt8}, data::Dict{T,Any})
    ptr = storestring!(ptr, bytestring(JSON.json(data)))
end

function pgdata{T<:AbstractString}(::Type{PostgresType{:jsonb}}, ptr::Ptr{UInt8}, data::Dict{T,Any})
    ptr = storestring!(ptr, bytestring(JSON.json(data)))
end

function pgdata(::Type{PostgresType{:_bool}}, ptr::Ptr{UInt8}, data::Vector{Bool})
    ptr = storestring!(ptr, string("{", join(map(x -> x ? "t" : "f", data), ','), "}"))
end

function pgdata(::Type{PostgresType{:_int8}}, ptr::Ptr{UInt8}, data::Vector{Int64})
    ptr = storestring!(ptr, string("{", join(data, ','), "}"))
end

function pgdata(::Type{PostgresType{:_int4}}, ptr::Ptr{UInt8}, data::Vector{Int32})
    ptr = storestring!(ptr, string("{", join(data, ','), "}"))
end

function pgdata(::Type{PostgresType{:_int2}}, ptr::Ptr{UInt8}, data::Vector{Int16})
    ptr = storestring!(ptr, string("{", join(data, ','), "}"))
end

function pgdata(::Type{PostgresType{:_float8}}, ptr::Ptr{UInt8}, data::Vector{Float64})
    ptr = storestring!(ptr, string("{", join(data, ','), "}"))
end

function pgdata(::Type{PostgresType{:_float4}}, ptr::Ptr{UInt8}, data::Vector{Float64})
    ptr = storestring!(ptr, string("{", join(data, ','), "}"))
end

function pgdata(::Type{PostgresType{:_varchar}}, ptr::Ptr{UInt8}, data::Vector{ASCIIString})
    ptr = storestring!(ptr, string("{", join(data, ','), "}"))
end

function pgdata(::Type{PostgresType{:_varchar}}, ptr::Ptr{UInt8}, data::Vector{UTF8String})
    ptr = storestring!(ptr, string("{", join(data, ','), "}"))
end

function pgdata(::Type{PostgresType{:_text}}, ptr::Ptr{UInt8}, data::Vector{ASCIIString})
    ptr = storestring!(ptr, string("{", join(data, ','), "}"))
end

function pgdata(::Type{PostgresType{:_text}}, ptr::Ptr{UInt8}, data::Vector{UTF8String})
    ptr = storestring!(ptr, string("{", join(data, ','), "}"))
end

# dbi
abstract Postgres <: DBI.DatabaseSystem

type PostgresDatabaseHandle <: DBI.DatabaseHandle
    ptr::Ptr{PGconn}
    status::ConnStatusType
    closed::Bool

    function PostgresDatabaseHandle(ptr::Ptr{PGconn}, status::ConnStatusType)
        new(ptr, status, false)
    end
end

type PostgresResultHandle
    ptr::Ptr{PGresult}
    types::Vector{DataType}
    nrows::Integer
    ncols::Integer
end

function PostgresResultHandle(result::Ptr{PGresult})
    status = PQresultStatus(result)
    if status == PGRES_TUPLES_OK || status == PGRES_SINGLE_TUPLE
        oids = @compat [OID{Int(PQftype(result, col))} for col in 0:(PQnfields(result)-1)]
        types = DataType[convert(PostgresType, x) for x in oids]
    else
        types = DataType[]
    end
    return PostgresResultHandle(result, types, PQntuples(result), PQnfields(result))
end

type PostgresStatementHandle <: DBI.StatementHandle
    db::PostgresDatabaseHandle
    stmt::AbstractString
    executed::Int
    paramtypes::Array{DataType}
    finished::Bool
    result::PostgresResultHandle

    function PostgresStatementHandle(db::PostgresDatabaseHandle, stmt::AbstractString, executed=0, paramtypes::Array{DataType}=DataType[])
        new(db, stmt, executed, paramtypes, false)
    end
end

function Base.copy(rh::PostgresResultHandle)
    PostgresResultHandle(PQcopyResult(result, PG_COPYRES_ATTRS | PG_COPYRES_TUPLES |
        PG_COPYRES_NOTICEHOOKS | PG_COPYRES_EVENTS), copy(rh.types), rh.ntuples, rh.nfields)
end

max_stars_repo_name: beffaCo/updatePostgresql

max_forks_repo_name: beffaCo/updatePostgresql


text:
using Sunny
using LinearAlgebra
using MPI

# MPI rank is 0-based index
rank, N_ranks = init_MPI()

# Linear schedule for external field strength
H_min = -4.0
H_max =  4.0
H = collect(range(H_min, H_max, length=N_ranks))[rank+1]

# Fixed temperature
kT = 1.0

function create_system(α::Float64)
    # Nearest-neighbor ferromagnetic Ising (z-dir) interactions
    J = -1.0
    FM = exchange(diagm([J, J, J]), Bond(1, 1, [1, 0, 0]))
    H_ext = external_field([0, 0, α])
    interactions = [FM, H_ext]

    # SC lattice -- extend lattice vector in z-direction to emulate 2D
    lvecs = Sunny.lattice_vectors(1.0, 1.0, 2.0, 90, 90 ,90)
    bvecs = [[0.0, 0.0, 0.0]]
    crystal = Crystal(lvecs, bvecs)

    # Make system and randomize spins
    extent = (120, 120, 1)
    g = 1 / Sunny.BOHR_MAGNETON
    system = SpinSystem(crystal, interactions, extent, [SiteInfo(1,1,g)])

    return system
end

α = H
system = create_system(α)
randflips!(system)

# Make replica for REMC
replica = Replica(IsingSampler(system, kT, 1), α)

# Define how the setting of α's changes
#  the sampling distribution or the system
function set_α!(replica::Replica, α::Float64)
    replica.α = α
    copy_sites = deepcopy(replica.sampler.system.sites)
    replica.sampler.system = create_system(α)
    replica.sampler.system.sites .= copy_sites
    reset_running_energy!(replica.sampler)
    reset_running_mag!(replica.sampler)
end

# Run feedback-optimized replica exchange
α_opt = run_FBO!(
    replica,
    set_α!;
    max_mcs = 200_000,
    rex_interval = 1,
    update_interval = 10_000,
    w = 0.0,
    #print_ranks_trajectory=[0],
    #trajectory_interval=10
)

# Run Replica Exchange MC (as PT) using optimized α schedule
run_REMC!(
    replica;
    therm_mcs=5000,
    measure_interval=10,
    rex_interval=1,
    max_mcs=500_000,
    bin_size=1.0,
    print_hist=true
)


text:
using SweepOperator, Test

# setup
n, p = 1000, 10
x = randn(n, p)
xtx = x'x

@testset "Sweep One By One" begin
    A = deepcopy(xtx)
    B = deepcopy(xtx)
    for j in 1:p
        sweep!(A, j)
        sweep!(A, j, true)
    end
    @test A ≈ B

    A = deepcopy(xtx)
    B = deepcopy(xtx)
    for j in 1:p
        sweep!(A, j, true)
    end
    for j in 1:p
        sweep!(A, j, true)
    end
    @test A ≈ B
end

@testset "Sweep All" begin
    A = deepcopy(xtx)
    B = deepcopy(xtx)
    sweep!(A, 1:p)
    sweep!(A, 1:p, true)
    @test A ≈ B
end

@testset "Linear Regression" begin
    y = x * collect(1.:p) + randn(n)
    xy = [x y]
    xytxy = xy'xy
    sweep!(xytxy, 1:p)
    @test xytxy[1:p, end] ≈ x\y
end