text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((2, 2),),
    ((2, 2), :ExpPSD),
    ((2, 2), :SOCExpPSD),
    ]
insts["fast"] = [
    ((3, 4),),
    ((3, 4), :ExpPSD),
    ((3, 4), :SOCExpPSD),
    ((10, 3),),
    ((5, 3), :ExpPSD),
    ((5, 3), :SOCExpPSD),
    ]
insts["various"] = [
    ((10, 20),),
    ((10, 15), :ExpPSD),
    ((10, 15), :SOCExpPSD),
    ((10, 25),),
    ((10, 20), :ExpPSD),
    ((10, 20), :SOCExpPSD),
    ((10, 30),),
    ]
return (MatrixCompletionJuMP, insts)


text:
export
    Expander,
    apply

"""
    mutable struct Expander{I,O} <: SFunc{I,O}

An Expander represents a model defined by a function that returns a
network. For a given value of an input, the conditional probability
distribution is provided by the network produced by the function
on that input.

For each such network, the expander manages a runtime to reason about it.
Expanders are lazy and do not evaluate the function until they have to.

As a result, there is state associated with Expanders. This is analysis
state rather than world state, i.e., it is the state of Scruff's
reasoning about the Expander. In keeping with Scruff design, Expanders
are immutable and all state associated with reasoning is stored in the
runtime that contains the expander. To support this, a runtime has three
fields of global state:

- `:subnets`: the expansions of all Expanders managed by the runtime
- `:subruntimes`: all the subruntimes recursively managed by this
  runtime through Expanders, keyed by the networks
- `:depth`: the depth to which Expanders in this runtime should be expanded

# Type parameters
- `I`: the input type(s) of the `Expander`
- `O`: the output type(s) of the `Expander`
"""
mutable struct Expander{I,O} <: SFunc{I,O}
    # TODO (MRH): SFunc params type parameter
    fn :: Function
    Expander(fn, I, O) = new{I,O}(memo(fn))
end

apply(expander::Expander, args...) = expander.fn(args...)

ext: jl


text:
module Newbe

__precompile__(true)

include("allmodules.jl")

### this allows to do Newbe.myColors() instal of Newbe.MyUtils.myColors()
using .MyUtils
export MyUtils

using .Behavior
export Behavior

using .Neurons
export Neurons

end

max_stars_repo_name: maalaria/Newbe

max_issues_repo_name: maalaria/Newbe

max_forks_repo_name: maalaria/Newbe


text:
module Python

using Dates,
    UnsafePointers,
    Libdl,
    Conda,
    Tables,
    TableTraits,
    IteratorInterfaceExtensions,
    Markdown,
    Requires,
    Compat,
    LinearAlgebra
using Base: @kwdef

# things not directly dependent on PyObject or libpython
include("utils.jl")

# Global configuration
# CONFIG gets populated by __init__
@kwdef mutable struct Config
    "Flags used to dlopen the Python library."
    dlopenflags::UInt32 = RTLD_LAZY | RTLD_DEEPBIND | RTLD_GLOBAL
    "Path to the Python executable."
    exepath::Union{String,Nothing} = nothing
    "Path to the Python library."
    libpath::Union{String,Nothing} = nothing
    "Handle to the open Python library."
    libptr::Ptr{Cvoid} = C_NULL
    "Used to set the Python prefix."
    pyhome::Union{String,Nothing} = nothing
    "pyhome as a Cwstring"
    pyhome_w::Vector{Cwchar_t} = []
    "Used to set the Python program name."
    pyprogname::Union{String,Nothing} = nothing
    "pyprogname as a Cwstring"
    pyprogname_w::Vector{Cwchar_t} = []
    "True if this is stackless Python."
    isstackless::Bool = false
    """True if Julia is embedded into Python (indicated by ENV["PYTHONJL_LIBPTR"] being set)."""
    isembedded::Bool = false
    "True if the Python interpreter is currently initialized."
    isinitialized::Bool = false
    "True if the Python interpreter was already initialized."
    preinitialized::Bool = false
    "The running Python version."
    version::VersionNumber = VersionNumber(0)
    "True if this is the Python in some Conda environment."
    isconda::Bool = false
    "If `isconda` is true, this is the Conda environment path."
    condaenv::String = ""
    "When true, automatically calls `pyplotshow` each time a notebook cell is evaluated."
    pyplotautoshow::Bool = true
    "When true, automatically calls `fix_qt_plugin_path` when `PyQt5`, `PyQt4`, `PySide` or `PySide2` is loaded."
    qtfix::Bool = true
    "When true, automatically sets `sys.last_traceback` when a Python exception is printed, so that `pdb.pm()` works."
    sysautolasttraceback::Bool = true
    "True if the Python input hook is currently running."
    inputhookrunning::Bool = false
    "When true, if being run inside an IPython kernel then integrate IO."
    ipythonintegration::Bool = true
    "When true, this is being run inside an IPython kernal."
    isipython::Bool = false
end
Base.show(io::IO, ::MIME"text/plain", c::Config) =
    for k in fieldnames(Config)
        println(io, k, ": ", repr(getfield(c, k)))
    end
const CONFIG = Config()

"""
    ispyreftype(::Type{T})

True if `T` is a wrapper type for a single Python reference.

Such objects must implement:
- `pyptr(o::T) =` the underlying pointer (or NULL on error)
- `unsafe_convert(::Type{CPyPtr}, o::T) = checknull(pyptr(o))`
"""
ispyreftype(::Type) = false

"""
    ispyref(x)

Equivalent to `ispyreftype(typeof(x))`.
"""
ispyref(x) = ispyreftype(typeof(x))

"""
    pyptr(o)

Retrieve the underlying Python object pointer from o.
"""
function pyptr end

convertref(::Type{T}, x) where {T} = ispyreftype(T) ? x : convert(T, x)
tryconvertref(::Type{T}, x) where {T} = ispyreftype(T) ? x : tryconvert(T, x)

# C API
include("cpython/CPython.jl")

const C = CPython
const CPyPtr = C.PyPtr

include("gil.jl")
include("eval.jl")
include("PyRef.jl")
include("PyCode.jl")
include("PyInternedString.jl")
include("builtins.jl")
include("PyException.jl")
include("PyObject.jl")
include("PyDict.jl")
include("PyList.jl")
include("PySet.jl")
include("PyIterable.jl")
include("PyIO.jl")
include("PyBuffer.jl")
include("PyArray.jl")
include("PyObjectArray.jl")
include("PyPandasDataFrame.jl")

include("julia.jl")
include("gui.jl")
include("matplotlib.jl")
include("ipython.jl")

include("init.jl")

end # module


text:
export gpuradixsort, gpuradixsort!

using CUDA
# CUDA.allowscalar(false)

const RADIX_SIZE = 11
const RADIX_MASK = UInt16(2^RADIX_SIZE-1)

# The winner of bencharmks/bencharmks-countmap
function radixhist!(buffer, v, number_of_radix)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = gridDim().x * blockDim().x

    local_buffer = CUDA.zeros(2^RADIX_SIZE, number_of_radix)

    for k in 1:number_of_radix
        for j = i:stride:length(v)
            bits_to_shift = (k-1)*RADIX_SIZE
            vj_shifted = v[j] >> bits_to_shift
            b = Int(vj_shifted & RADIX_MASK) + 1

            local_buffer[b, k] += 1
        end
    end

    for k in 1:number_of_radix
        for i in 1:2^RADIX_SIZE
            CUDA.@atomic buffer[i, k] = buffer[i, k] + local_buffer[i, k]
        end
    end

    return
end

# count
function radixhist(v::CuArray{T}; threads = 256, blocks = 1024) where {T}
    number_of_radix = ceil(Int, 8sizeof(T) / RADIX_SIZE)

    buffer = CUDA.zeros(Int, 2^RADIX_SIZE, number_of_radix)
    CUDA.@sync @cuda threads = threads blocks = blocks radixhist!(buffer, v, number_of_radix)

    for i in 1:number_of_radix
        buffer[:, i] = cumsum(buffer[:, i])
    end
    buffer
end

gpuradixsort(v) = gpuradixsort!(copy(v))

function gpuradixsort!(vs::AbstractVector{T}) where T
    bin = collect(radixhist(cu(vs)))
    ts = similar(vs)

    # use the histogram to sort the data
    hi = length(vs)
    lo = 1
    len = hi-lo+1
    iters = ceil(Int, 8sizeof(T) / RADIX_SIZE)
    for j = 1:iters
        # Unroll first data iteration, check for degenerate case
        v = vs[hi]
        idx = Int((v >> ((j-1)*RADIX_SIZE)) & RADIX_MASK) + 1

        cbin = bin[:,j]
        ci = cbin[idx]
        ts[ci] = vs[hi]
        cbin[idx] -= 1

        # Finish the loop...
        @inbounds for i in hi-1:-1:lo
            v = vs[i]
            idx = Int((v >> ((j-1)*RADIX_SIZE)) & RADIX_MASK) + 1
            ci = cbin[idx]
            ts[ci] = vs[i]
            cbin[idx] -= 1
        end
        vs,ts = ts,vs
    end

    vs
end


text:
using LibGit2

const PKGDIR = Pkg.dir("LibGit2")
const TESTDIR = joinpath(PKGDIR, "test")
const FIXTURE_DIR = joinpath(TESTDIR, "fixtures")
const LIBGIT2_FIXTURE_DIR = joinpath(PKGDIR, "vendor", "libgit2", "tests", "resources")

context(f::Function) = f()
context(f::Function, s::AbstractString) = (println(s); context(f))

cleanup_dir(p) = begin
    if isdir(p)
        rm(p, recursive=true)
    end
end

function cleanup!(repo::GitRepo, dir::AbstractString)
    close(repo)
    LibGit2.free!(repo)
    Base.gc()
    rm(dir, recursive=true)
end

function remote_transport_test(f::Function)
    tmp_dir = mktempdir()
    test_repo = init_repo(tmp_dir, bare=false)
    test_repo_dir = joinpath(TESTDIR, joinpath("fixtures", "testrepo.git", "."))
    try
        test_remote = remote_add!(test_repo, "origin", test_repo_dir)
        f(test_repo, test_remote)
    finally
        cleanup!(test_repo, tmp_dir)
    end
end
remote_transport_test(f::Function, s::AbstractString) =
    (println(s); remote_transport_test(f))

copy_recur(path::AbstractString, dest::AbstractString) = begin
    apath, adest = abspath(path), abspath(dest)
    if isfile(apath)
        cp(apath, adest)
        return
    end
    contents = readdir(path)
    for c in contents
        pc, dc = joinpath(apath, c), joinpath(adest, c)
        if isfile(pc)
            cp(pc, dc)
        elseif isdir(pc)
            isdir(dc) || mkdir(dc)
            copy_recur(pc, dc)
        end
    end
end

function repo_clone_test(f::Function)
    tmp_dir = mktempdir()
    source = joinpath(TESTDIR, "fixtures", "testrepo.git")
    try
        f(source, tmp_dir)
    finally
        cleanup_dir(tmp_dir)
    end
end
repo_clone_test(f::Function, s::AbstractString) = (println(s); repo_clone_test(f))

function with_test_index(f::Function)
    try
        test_index_path = joinpath(TESTDIR, "fixtures", "testrepo.git", "index")
        test_index = GitIndex(test_index_path)
        f(test_index, test_index_path)
    finally
    end
end
with_test_index(f::Function, s::AbstractString) =
    (println(s); with_test_index(f))

function create_test_repo(test_path)
    if isdir(abspath(test_path))
        rm(test_path, recursive=true)
    end
    repo = init_repo(test_path)
    open(joinpath(test_path, "README"), "w") do fh
        write(fh, "foo\n")
    end
    return repo
end

function seed_test_repo(repo)
    sig = Signature("test", "test@test.com")
    idx = GitIndex(repo)
    push!(idx, "README")
    tree_id = write_tree!(idx)

    msg = "This is a commit\n"
    tree = lookup(GitTree, repo, tree_id)
    commit_id = commit(repo, "HEAD", msg, sig, sig, tree)
    return commit_id, tree_id
end

function clone_sandbox(repo_path::AbstractString)
        tmp_dir = mktempdir()
        repo_path = abspath(repo_path)
        run(`git clone --quiet -- $repo_path $tmp_dir`)
        return GitRepo(tmp_dir), tmp_dir
end

function rename_git_files!(tmp_dir::AbstractString)
    # rename git specific files and folders in cloned test dir
    for d in [".gitted", "dot_git"]
        d = joinpath(tmp_dir, d)
        isdir(d) && mv(d, joinpath(tmp_dir, ".git"))
    end
    for f in ["gitattributes", "gitignore", "gitmodules"]
        p = joinpath(tmp_dir, f)
        isfile(p) && mv(p, joinpath(tmp_dir, "."*f))
    end
end

function setup_sandbox(repo_name::AbstractString)
    tmp_dir = mktempdir()
    fixture_dir = joinpath(LIBGIT2_FIXTURE_DIR, repo_name)
    copy_recur(fixture_dir, tmp_dir)
    rename_git_files!(tmp_dir)
    return GitRepo(tmp_dir), tmp_dir
end

function sandboxed_test(f::Function, reponame::AbstractString)
    repo, tmp_dir = setup_sandbox(reponame)
    try
        f(repo, tmp_dir)
    finally
        cleanup!(repo, tmp_dir)
    end
end
sandboxed_test(f::Function, reponame::AbstractString, s::AbstractString) = (println(s);
                                                            sandboxed_test(f, reponame))

function sandboxed_clone_test(f::Function, reponame::AbstractString)
    remote, tmp_dir1 = setup_sandbox(reponame)
    GitConfig(remote)["core.bare"] = true
    repo, tmp_dir2 = clone_sandbox(tmp_dir1)
    try
        f(repo, remote, tmp_dir1)
    finally
        cleanup!(repo, tmp_dir1)
        cleanup!(remote, tmp_dir2)
    end
end
sandboxed_clone_test(f::Function, reponame::AbstractString, s::AbstractString) =
    (println(s); sandboxed_clone_test(f, reponame))

function sandboxed_checkout_test(f::Function)
    test_repo,  test_repo_dir  = setup_sandbox("testrepo")
    test_clone, test_clone_dir = clone_sandbox(test_repo_dir)
    test_bare,  test_bare_dir  = setup_sandbox("testrepo.git")
    GitConfig(test_bare)["core.bare"] = true
    try
        f(test_repo, test_clone, test_bare)
    finally
        cleanup!(test_repo, test_repo_dir)
        cleanup!(test_clone, test_clone_dir)
        cleanup!(test_bare, test_bare_dir)
    end
end
sandboxed_checkout_test(f::Function, s::AbstractString) =
    (println(s); sandboxed_checkout_test(f))

function with_repo_access(f::Function)
    path = joinpath(TESTDIR, "fixtures", "testrepo.git")
    repo = GitRepo(path)
    try
        f(repo, path)
    finally
        close(repo)
        LibGit2.free!(repo)
        Base.gc()
    end
end
with_repo_access(f::Function, s::AbstractString) =
    (println(s); with_repo_access(f))

function with_tmp_repo_access(f::Function)
    tmp_dir  = mktempdir()
    repo_dir = joinpath(TESTDIR, "fixtures", "testrepo.git", ".")
    run(`git clone --quiet -- $repo_dir $tmp_dir`)
    repo = GitRepo(tmp_dir)
    try
        f(repo, tmp_dir)
    finally
        cleanup!(repo, tmp_dir)
    end
end
with_tmp_repo_access(f::Function, s::AbstractString) =
    (println(s); with_tmp_repo_access(f))

function with_new_repo(f::Function; bare=false)
    tmp_dir = mktempdir()
    repo = LibGit2.init_repo(tmp_dir; bare=bare)
    try
        f(repo, tmp_dir)
    finally
        cleanup!(repo, tmp_dir)
    end
end
with_new_repo(f::Function, s::AbstractString; bare::Bool=false) =
    (println(s); with_new_repo(f; bare=bare))

touch_test(parent::AbstractString, file::AbstractString, content=nothing) = begin
    fp = joinpath(parent, file)
    dr = dirname(fp)
    isdir(dr) ||  mkdir(dr)
    open(fp, "w") do fh
        if content != nothing
            write(fh, content)
        end
    end
    return fp
end

with_test_repo(f::Function, reponame::AbstractString, clone::Bool=false) = begin
    tmp_dir = mktempdir()
    fixture_dir = joinpath(FIXTURE_DIR, reponame)
    if clone
        run(`git clone --quiet -- $fixture_dir $tmp_dir`)
    else
        copy_recur(fixture_dir, tmp_dir)
    end
    rename_git_files!(tmp_dir)
    repo = GitRepo(tmp_dir)
    try
        f(repo, tmp_dir)
    finally
        cleanup!(repo, tmp_dir)
    end
end

with_standard_test_repo(f::Function, clone::Bool) = with_test_repo(f, "testrepo_wd", clone)
with_standard_test_repo(f::Function, s::AbstractString) =
    (println(s); with_standard_test_repo(f, false))

with_merged_test_repo(f::Function, clone::Bool) = with_test_repo(f, "mergedrepo_wd", clone)
with_merged_test_repo(f::Function, s::AbstractString) =
    (println(s); with_merged_test_repo(f, false))

with_merge_test_repo(f::Function, clone::Bool) =
    with_test_repo(f, "merge_testrepo_wd", false)
with_merge_test_repo(f::Function, s::AbstractString) =
    (println(s); with_merge_test_repo(f, false))


text:
# All processes related to survival:
# - growth
# - density-independent mortality
# - competition

"""
    survive!(patch)

Density independent survival of individuals in a patch. The actual mortality
probability is calculated in one of three ways: (1) with a metabolic formula,
modified by the passed `mortality` variable and an individual's temperature
adaptation; (2) via the individual's habitat (i.e. precipitation) adaptation;
(3) as a global constant probability defined by `mortality`.
Which mode is selected depends on `mortalitytype`.
"""
function survive!(patch::Patch)
    mtype = setting("mortalitytype")
    temp = patch.temp
    idx = 1
    while idx <= size(patch.community,1)
        (patch.community[idx].marked) && continue
        dies = false
        if mtype == "metabolic"
            mass = patch.community[idx].size
            deathrate = setting("mortality") * mass^(-1/4) * exp(-act/(boltz*temp))
            dieprob = (1 - exp(-deathrate))
            (rand() * patch.community[idx].tempadaptation < dieprob) && (dies = true)
        elseif mtype == "habitat"
            precadaptation = pdf(Normal(patch.community[idx].traits["precopt"], patch.community[idx].traits["prectol"]), patch.prec)
            maxprecadapt = pdf(Normal(patch.prec, setting("perfecttol")), patch.prec)
            fitness = precadaptation / maxprecadapt
            fitness > 1 && (fitness = 1)
            fitness < 0 && (fitness = 0)
            (rand() < setting("mortality") + (1 - fitness)) && (dies = true)
            #(rand() < setting("mortality")+(1-patch.community[idx].precadaptation)) && (dies = true)
        elseif mtype == "global"
            (rand() < setting("mortality")) && (dies = true)
        else
            simlog("Invalid mortality type $mtype", 'e')
        end
        if dies
            @simlog("$(idstring(patch.community[idx])) has died.", 'd')
            splice!(patch.community, idx)
            continue
        end
        idx += 1
    end
end

"""
    survive!(world)

World-wide mortality. Sounds apocalyptic, but is just a fact of life.
"""
function survive!(world::Array{Patch,1})
    for patch in world
        (patch.isisland || !setting("static")) && survive!(patch) # pmap(!,patch) ???
    end
end

"""
    grow!(patch, growthrate, capgrowth)

Growth of individuals in the given patch. The actual growthrate is calculated
with a metabolic formula, modified by the passed `growthrate` variable.
"""
function grow!(patch::Patch, growthrate::Float64, capgrowth::Bool)
    temp = patch.temp
    idx = 1
    while idx <= size(patch.community,1)
        if !patch.community[idx].marked
            repsize = patch.community[idx].traits["repsize"]
            mass = patch.community[idx].size
            if mass <= repsize # stop growth if reached repsize
                growth = growthrate * mass^(3/4) * exp(-act/(boltz*temp))
                newmass = mass + growth
                if capgrowth && newmass > patch.community[idx].traits["repsize"]
                    newmass = patch.community[idx].traits["repsize"]
                end
                if newmass > 0 && mass > 0
                    patch.community[idx].size = newmass
                else
                    splice!(patch.community, idx)
                    continue
                end
            end
        end
        idx += 1
    end
end

"""
    grow!(world)

Carry out growth for all patches.
"""
function grow!(world::Array{Patch,1})
    for patch in world
        # pmap(!,patch) ???
        if patch.isisland || !setting("static")
            grow!(patch, setting("growthrate"), setting("capgrowth"))
        end
    end
end

"""
    compete!(patch)

While there are too many organisms in a patch, pick two at random and kill the
one that is less adapted to the local precipitation levels.
"""
function compete!(patch::Patch)
    totalmass = sum(map(x -> x.size, patch.community))
    while totalmass >= patch.capacity # occupied capacity larger than available
        firstind, secondind = rand(eachindex(patch.community), 2)
        firstind == secondind && length(eachindex(patch.community)) > 1 && continue
        if patch.community[firstind].precadaptation < patch.community[secondind].precadaptation
            totalmass -= patch.community[firstind].size
            splice!(patch.community, firstind) # profiling: expensive!
        elseif patch.community[firstind].precadaptation > patch.community[secondind].precadaptation
            totalmass -= patch.community[secondind].size
            splice!(patch.community, secondind) # profiling: expensive!
        else
            victim = rand([firstind, secondind])
            totalmass -= patch.community[victim].size
            splice!(patch.community, victim) # profiling: expensive!
        end
    end
end

"""
    compete!(world)

Carry out competition on all patches.
"""
function compete!(world::Array{Patch,1})
    for patch in world
        (patch.isisland || !setting("static")) && compete!(patch) # pmap(!,patch) ???
    end
end

max_stars_repo_name: CCTB-Ecomods/gemm

max_issues_repo_name: CCTB-Ecomods/gemm

max_forks_repo_name: CCTB-Ecomods/gemm


text:
# This file is a part of ForwardDiffPullbacks.jl, licensed under the MIT License (MIT).

__precompile__(true)

"""
    ForwardDiffPullbacks

ChainRulesCore compatible pullbacks using ForwardDiff.
"""
module ForwardDiffPullbacks

using LinearAlgebra

import ChainRulesCore
import ForwardDiff
import StaticArrays

using ChainRulesCore: AbstractTangent, Tangent, NoTangent, ZeroTangent, AbstractThunk, unthunk

# using Requires

include("dual_numbers.jl")
include("chainrules_types.jl")
include("fwd_back.jl")
include("broadcasting.jl")
include("with_forwarddiff.jl")
include("rrules.jl")

function __init__()
    # Possible extensions:
    # @require Nabla = "49c96f43-aa6d-5a04-a506-44c7070ebe78" include("nabla_support.jl")
    # @require ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" include("reversediff_support.jl")
    # @require Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" include("tracker_support.jl")
end

end # module


text:
using Images
using ImageView
using EMIRT

println("compare image and segmentation")
println("Usage: julia compare_img_seg.jl /path/of/img.h5 /path/of/seg.h5")
@assert length(ARGS) >= 2

fimg = ARGS[1]
fseg = ARGS[2]

img = readimg(fimg)[:,:,1:64]
seg = readseg(fseg)[:,:,1:64]

include(joinpath(Pkg.dir(), "EMIRT/plugins/show.jl"))

imgc, imgslice = show(img, seg)

#If we are not in a REPL
if (!isinteractive())

    # Create a condition object
    c = Condition()

    # Get the main window (A Tk toplevel object)
    win = toplevel(imgc)

    # Notify the condition object when the window closes
    bind(win, "<Destroy>", e->notify(c))

    # Wait for the notification before proceeding ...
    wait(c)
end


text:
"""
    OrderedDither(mat::AbstractMatrix)

Generalized ordered dithering algorithm using a threshold map.
Takes a normalized threshold matrix `mat`.

When applying the algorithm to an image, the threshold matrix is repeatedly tiled
to match the size of the image. It is then applied as a per-pixel threshold map.
Optionally, this final threshold map can be inverted by selecting `invert_map=true`.
"""
struct OrderedDither{T<:AbstractMatrix{<:Rational},R<:Real} <: AbstractDither
    mat::T
    color_error_multiplier::R
end
function OrderedDither(mat; invert_map=false, color_error_multiplier=0.5)
    if invert_map
        return OrderedDither(1 .- mat, color_error_multiplier)
    end
    return OrderedDither(mat, color_error_multiplier)
end

function binarydither!(alg::OrderedDither, out::GenericGrayImage, img::GenericGrayImage)
    # eagerly promote to the same eltype to make for-loop faster
    FT = floattype(eltype(img))
    mat = FT.(alg.mat)

    T = eltype(out)
    black, white = T(0), T(1)

    # Precompute lookup tables for modulo indexing of threshold matrix
    matsize = size(mat)
    rlookup = [mod1(i, matsize[1]) for i in 1:size(img)[1]]
    clookup = [mod1(i, matsize[2]) for i in 1:size(img)[2]]

    @inbounds @simd for i in CartesianIndices(img)
        r, c = Tuple(i)
        out[i] = ifelse(img[i] > mat[rlookup[r], clookup[c]], white, black)
    end
    return out
end

# Using Pattern Dithering by Thomas Knoll / Adobe Inc. Patent expired in 2019.
# https://patents.google.com/patent/US6606166B1/en
# Implemented according to Joel Yliluoma's pseudocode
# https://bisqwit.iki.fi/story/howto/dither/jy/
function colordither(
    alg::OrderedDither,
    img::GenericImage,
    cs::AbstractVector{<:Pixel},
    metric::DifferenceMetric,
)
    cs_lab = Lab.(cs)
    cs_xyz = XYZ.(cs)

    # Precompute lookup tables for modulo indexing of threshold matrix
    mat = numerator.(alg.mat)
    nmax = maximum(mat)
    matsize = size(mat)
    rlookup = [mod1(i, matsize[1]) for i in 1:size(img)[1]]
    clookup = [mod1(i, matsize[2]) for i in 1:size(img)[2]]

    # Allocate matrices
    candidates = Array{UInt8}(undef, nmax)
    index = Matrix{UInt8}(undef, size(img)...)

    @inbounds for I in CartesianIndices(img)
        r, c = Tuple(I)
        err = zero(XYZ)
        px = XYZ(img[I])

        for j in 1:nmax
            col = px + alg.color_error_multiplier * err
            idx = _closest_color_idx(col, cs_lab, metric)

            # We are in loop (idx ↔ err) if we already computed `idx` as the closest color
            # after the first iteration. Breaking out of this loops lets us avoid calling
            # the expensive closest color computation.
            if (j > 2) && (idx in candidates[2:(j - 1)])
                # Find the range of the loop in `candidates`:
                looprange = (findfirst(i -> i == idx, candidates[2:(j - 1)]) + 1):(j - 1)
                # Fill the rest of `candidates` with this loop:
                candidates[j:end] .= Iterators.take(
                    Iterators.cycle(candidates[looprange]), nmax - j + 1
                )
                break
            else
                candidates[j] = idx
                err = px - cs_xyz[idx]
            end
        end
        # Sort candidates by luminance (dark to bright)
        index[I] = partialsort!(candidates, mat[rlookup[r], clookup[c]]; by=i -> cs_lab[i].l)
    end
    return index
end

"""
    Bayer(; level)

Ordered dithering using the Bayer matrix as a threshold matrix.
The Bayer matrix is of dimension ``2^{n+1} \\times 2^{n+1}``, where ``n`` is the `level`,
which defaults to `1`.

[1]  Bayer, B.E., "An Optimum Method for Two-Level Rendition of Continuous
     Tone Pictures," IEEE International Conference on Communications,
     Conference Records, 1973, pp. 26-11 to 26-15.
"""
function Bayer(; level=1, kwargs...)
    bayer = bayer_matrix(level) .+ 1
    return OrderedDither(bayer//(2^(2 * level + 2) + 1); kwargs...)
end

"""
    bayer_matrix(n::Int)::AbstractMatrix{Int}

Contruct (un-normalized) Bayer matrices of level `n` through recursive definition.
"""
function bayer_matrix(n::Int)::AbstractMatrix{Int}
    n < 0 && throw(DomainError(n, "Bayer matrix only defined for n ≥ 0."))

    if n == 0
        b = [0 2; 3 1]
    else
        bₙ₋₁ = 4 * bayer_matrix(n - 1)
        b = [(bₙ₋₁) (bₙ₋₁.+2); (bₙ₋₁.+3) (bₙ₋₁.+1)]
    end
    return b
end

"""
    ClusteredDots()

Clustered dots ordered dithering.
Uses ``6 \\times 6`` threshold matrix `CLUSTERED_DOTS_MAT`.
"""
ClusteredDots(; kwargs...) = OrderedDither(CLUSTERED_DOTS_MAT; kwargs...)
const CLUSTERED_DOTS_MAT =
    [
        35 30 18 22 31 36
        29 15 10 17 21 32
        14 9 5 6 16 20
        13 4 1 2 11 19
        28 8 3 7 24 25
        34 27 12 23 26 33
    ]//37

"""
    CentralWhitePoint()

Central white point ordered dithering.
Uses ``6 \\times 6`` threshold matrix `CENTRAL_WHITE_POINT_MAT`.
"""
CentralWhitePoint(; kwargs...) = OrderedDither(CENTRAL_WHITE_POINT_MAT; kwargs...)
const CENTRAL_WHITE_POINT_MAT =
    [
        35 26 22 18 30 34
        31 14 10 6 13 25
        19 7 2 1 9 21
        23 11 3 4 5 17
        27 15 8 12 16 29
        36 32 20 24 28 33
    ]//37

"""
    BalancedCenteredPoint()

Balanced centered point ordered dithering.
Uses ``6 \\times 6`` threshold matrix `BALANCED_CENTERED_POINT_MAT`.
"""
BalancedCenteredPoint(; kwargs...) = OrderedDither(BALANCED_CENTERED_POINT_MAT; kwargs...)
const BALANCED_CENTERED_POINT_MAT =
    [
        31 23 17 22 34 36
        25 12 8 10 27 29
        14 6 1 3 15 20
        16 4 2 5 13 19
        28 9 7 11 26 30
        33 21 18 24 32 35
    ]//37

"""
    Rhombus()

Diagonal ordered matrix with balanced centered points.
Uses ``8 \\times 8`` threshold matrix `RHOMBUS_MAT`.
"""
Rhombus(; kwargs...) = OrderedDither(RHOMBUS_MAT; kwargs...)
const RHOMBUS_MAT =
    [
        14 10 6 13 19 23 27 20
        7 2 1 9 26 31 32 24
        11 3 4 5 22 30 29 28
        15 8 12 16 18 25 21 17
        19 23 27 20 14 10 6 13
        26 31 32 24 7 2 1 9
        22 30 29 28 11 3 4 5
        18 25 21 17 15 8 12 16
    ]//33


text:
#!/usr/bin/env julia

export pe5

pe5()=lcm(1:20)

if !haskey(ENV,"PROJECT_EULER_WITHOUT_RESULT")
    pe5()|>print
end

max_stars_repo_path: 5. Smallest multiple/5.jl

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_path: 5. Smallest multiple/5.jl

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_path: 5. Smallest multiple/5.jl

max_forks_repo_name: miRoox/ProjectEuler


text:
function qsfit(source::QSO{TRecipe}) where TRecipe <: DefaultRecipe
    elapsed = time()
    @assert length(source.specs) == 1
    pspec = PreparedSpectrum(source, id=1)

    model = Model(pspec.domain)
    model[:Continuum] = SumReducer([])
    model[:main] = SumReducer([])
    push!(model[:main].list, :Continuum)
    select_reducer!(model, :main)
    delete!(model.revals, :default_sum)

    # TODO if source.options[:instr_broadening]
    # TODO     GFit.set_instr_response!(model[1], (l, f) -> instrumental_broadening(l, f, source.spectra[id].resolution))
    # TODO end

    println(logio(source), "\nFit continuum components...")
    QSFit.add_qso_continuum!(source, pspec, model)
    QSFit.add_host_galaxy!(source, pspec, model)
    QSFit.add_balmer_cont!(source, pspec, model)
    fitres = fit!(source, model, pspec)
    QSFit.renorm_cont!(source, pspec, model)
    freeze(model, :qso_cont)
    haskey(model, :galaxy)  &&  freeze(model, :galaxy)
    haskey(model, :balmer)  &&  freeze(model, :balmer)
    evaluate!(model)

    println(logio(source), "\nFit iron templates...")
    model[:Iron] = SumReducer([])
    push!(model[:main].list, :Iron)
    QSFit.add_iron_uv!( source, pspec, model)
    QSFit.add_iron_opt!(source, pspec, model)

    if length(model[:Iron].list) > 0
        fitres = fit!(source, model, pspec)
        haskey(model, :ironuv   )  &&  freeze(model, :ironuv)
        haskey(model, :ironoptbr)  &&  freeze(model, :ironoptbr)
        haskey(model, :ironoptna)  &&  freeze(model, :ironoptna)
    end
    evaluate!(model)

    println(logio(source), "\nFit known emission lines...")
    QSFit.add_emission_lines!(source, pspec, model)
    QSFit.guess_emission_lines!(source, pspec, model)
    QSFit.add_patch_functs!(source, pspec, model)

    fitres = fit!(source, model, pspec)
    for lname in keys(pspec.lcs)
        freeze(model, lname)
    end

    println(logio(source), "\nFit unknown emission lines...")
    QSFit.add_unknown_lines!(source, pspec, model)

    println(logio(source), "\nLast run with all parameters free...")
    thaw(model, :qso_cont)
    haskey(model, :galaxy   )  &&  thaw(model, :galaxy)
    haskey(model, :balmer   )  &&  thaw(model, :balmer)
    haskey(model, :ironuv   )  &&  thaw(model, :ironuv)
    haskey(model, :ironoptbr)  &&  thaw(model, :ironoptbr)
    haskey(model, :ironoptna)  &&  thaw(model, :ironoptna)
    for lname in keys(pspec.lcs)
        thaw(model, lname)
    end
    for j in 1:source.options[:n_unk]
        cname = Symbol(:unk, j)
        if model[cname].norm.val > 0
            thaw(model, cname)
        else
            freeze(model, cname)
        end
    end
    fitres = fit!(source, model, pspec)

    if QSFit.neglect_weak_features!(source, pspec, model, fitres)
        println(logio(source), "\nRe-run fit...")
        fitres = fit!(source, model, pspec)
    end

    println(logio(source))
    show(logio(source), fitres)

    out = QSFitResults(source, pspec, model, fitres)
    elapsed = time() - elapsed
    println(logio(source), "\nElapsed time: $elapsed s")
    close_logio(source)
    return out
end


text:
using TuringModels, StatsFuns

Turing.setadbackend(:reversediff);

d = DataFrame(CSV.read(joinpath(@__DIR__, "..", "..", "data", "chimpanzees.csv"),
  delim=';'));
size(d) # Should be 504x8

# pulled_left, condition, prosoc_left
@model m10_3(y, x₁, x₂) = begin
    α ~ Normal(0, 10)
    βp ~ Normal(0, 10)
    βpC ~ Normal(0, 10)

    logits = α .+ (βp .+ βpC * x₁) .* x₂
    y .~ BinomialLogit.(1, logits)
end;

chns = sample(m10_3(d[:,:pulled_left], d[:,:condition], d[:,:prosoc_left]),
  Turing.NUTS(0.65), 2000);

# Rethinking result

m_10_03t_result = "
      Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
 a    0.05   0.13      -0.15       0.25  3284    1
 bp   0.62   0.22       0.28       0.98  3032    1
 bpC -0.11   0.26      -0.53       0.29  3184    1
";

# Describe the draws

describe(chns)

# End of m10.03t.jl


text:
export parse

#
# Parse given `ms` into a Dict().
#

function parse (ms::Int)
  qr = divrem(ms, 1000)
  ms = qr[2]

  qr = divrem(qr[1], 60)
  s = qr[2]

  qr = divrem(qr[1], 60)
  m = qr[2]

  qr = divrem(qr[1], 24)
  h = qr[2]
  d = qr[1]

  return [
    "days" => d,
    "hours" => h,
    "minutes" => m,
    "seconds" => s,
    "milliseconds" => ms
  ]
end

max_stars_repo_name: queckezz/ms-parse

max_issues_repo_name: queckezz/ms-parse

max_forks_repo_name: queckezz/ms-parse


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using Base: @propagate_inbounds

@doc doc"""
    OnlineUvMean{T<:AbstractFloat}

*BAT-internal, not part of stable public API.*

Univariate mean implemented via Kahan-Babuška-Neumaier summation.
"""
mutable struct OnlineUvMean{T<:AbstractFloat}
    sum_v::DoubleFloat{T}
    sum_w::DoubleFloat{T}

    OnlineUvMean{T}() where {T<:AbstractFloat} = new{T}(zero(DoubleFloat{T}), zero(DoubleFloat{T}))

    OnlineUvMean{T}(sum_v::Real, sum_w::Real) where {T<:AbstractFloat} = new{T}(sum_v, sum_w)
end

OnlineUvMean() = OnlineUvMean{Float64}()

@inline Base.getindex(omn::OnlineUvMean{T}) where {T<:AbstractFloat} = T(omn.sum_v / omn.sum_w)

function Base.empty!(target::OnlineUvMean{T}) where T
    target.sum_v = zero(DoubleFloat{T})
    target.sum_w = zero(DoubleFloat{T})

    target
end

function Base.merge!(target::OnlineUvMean{T}, others::OnlineUvMean...) where {T}
    sum_v = target.sum_v
    sum_w = target.sum_w

    @inbounds @simd for x in others
        sum_w += x.sum_w
        sum_v += x.sum_v
    end

    target.sum_w = sum_w
    target.sum_v = sum_v

    target
end

@inline function Base.push!(omn::OnlineUvMean{T}, data::Real, weight::Real = 1) where {T}
    # Workaround for lack of promotion between, e.g., Float32 and DoubleFloat{Float64}
    weight_conv = T(weight)

    omn.sum_v += DoubleFloat{T}(weight_conv * data)
    omn.sum_w += DoubleFloat{T}(weight_conv)
    omn
end

function reweight_relative!(target::OnlineUvMean{T}, reweighting_factor::Real) where T
    target.sum_w *= reweighting_factor
    target
end

@doc doc"""
    OnlineUvVar{T<:AbstractFloat,W}

*BAT-internal, not part of stable public API.*

Implementation based on variance calculation Algorithms of Welford and West.

`W` must either be `Weights` (no bias correction) or one of `AnalyticWeights`,
`FrequencyWeights` or `ProbabilityWeights` to specify the desired bias
correction method.
"""
mutable struct OnlineUvVar{T<:AbstractFloat,W}
    n::Int64
    sum_w::DoubleFloat{T}
    sum_w2::DoubleFloat{T}
    mean_x::T
    s::T

    OnlineUvVar{T,W}() where {T<:AbstractFloat,W} =
        new{T,W}(
            zero(Int64), zero(DoubleFloat{T}), zero(DoubleFloat{T}),
            zero(T), zero(T)
        )

    OnlineUvVar{T,W}(n::Int64, sum_w::DoubleFloat{T}, sum_w2::DoubleFloat{T}, mean_x::T, s::T) where {T<:AbstractFloat,W} =
        new{T,W}(
            n, sum_w, sum_w2, mean_x, s
        )
end

OnlineUvVar() = OnlineUvVar{Float64, ProbabilityWeights}()

@propagate_inbounds Base.getindex(ocv::OnlineUvVar{T, Weights}) where {T} =
    ifelse(ocv.sum_w > zero(T), T(ocv.s / ocv.sum_w), T(NaN))

@propagate_inbounds function Base.getindex(ocv::OnlineUvVar{T, AnalyticWeights}) where {T}
    d = ocv.sum_w - ocv.sum_w2 / ocv.sum_w
    ifelse(ocv.sum_w > zero(T) && d > zero(T), T(ocv.s / d), T(NaN))
end

@propagate_inbounds Base.getindex(ocv::OnlineUvVar{T, FrequencyWeights}) where {T} =
    ifelse(ocv.sum_w > one(T), T(ocv.s / (ocv.sum_w - one(T))), T(NaN))

@propagate_inbounds Base.getindex(ocv::OnlineUvVar{T, ProbabilityWeights}) where {T} =
    ifelse(ocv.n > one(T) && ocv.sum_w > zero(T), T(ocv.s * ocv.n / ((ocv.n - one(T)) * ocv.sum_w)), T(NaN))

function Base.empty!(target::OnlineUvVar{T,W}) where {T,W}
    target.sum_w = zero(DoubleFloat{T})
    target.sum_w2 = zero(DoubleFloat{T})
    target.mean_x = zero(T)
    target.s = zero(T)

    target
end

function Base.merge!(target::OnlineUvVar{T,W}, others::OnlineUvVar...) where {T,W}
    n = target.n
    sum_w = target.sum_w
    sum_w2 = target.sum_w2
    mean_x = target.mean_x
    s = target.s

    @inbounds @simd for x in others
        n += x.n

        dx = mean_x - x.mean_x

        new_sum_w = (sum_w + x.sum_w)
        mean_x = (sum_w * DoubleFloat(mean_x) + x.sum_w * DoubleFloat(x.mean_x)) / new_sum_w

        s += x.s + sum_w * x.sum_w / new_sum_w * dx * dx

        sum_w = new_sum_w
        sum_w2 += x.sum_w2

    end

    target.n = n
    target.sum_w = sum_w
    target.sum_w2 = sum_w2
    target.mean_x = mean_x
    target.s = s

    target
end

@inline function Base.push!(ocv::OnlineUvVar{T}, data::Real, weight::Real = 1) where T
    # Ignore zero weights (can't be handled)
    if weight ≈ 0
        return ocv
    end

    # Workaround for lack of promotion between, e.g., Float32 and DoubleFloat{Float64}
    weight_conv = T(weight)

    n = ocv.n
    sum_w = ocv.sum_w
    sum_w2 = ocv.sum_w2
    mean_x = ocv.mean_x
    s = ocv.s

    n += one(n)
    sum_w += DoubleFloat{T}(weight_conv)
    sum_w2 += DoubleFloat{T}(weight_conv^2)
    dx = data - mean_x
    new_mean_x = mean_x + dx * weight_conv / sum_w
    new_dx = data - new_mean_x

    s = muladd(dx, weight_conv * new_dx, s)
    mean_x = new_mean_x

    ocv.n = n
    ocv.sum_w = sum_w
    ocv.sum_w2 = sum_w2
    ocv.mean_x = mean_x
    ocv.s = s

    ocv
end

function reweight_relative!(target::OnlineUvVar{T,W}, reweighting_factor::Real) where {T,W}
    target.sum_w *= reweighting_factor
    target.sum_w2 *= reweighting_factor
    target.s *= reweighting_factor

    target
end

mutable struct BasicUvStatistics{T<:Real,W}
    mean::OnlineUvMean{T}
    var::OnlineUvVar{T,W}
    maximum::T
    minimum::T

    BasicUvStatistics{T,W}() where {T<:Real,W} =
        new(OnlineUvMean{T}(), OnlineUvVar{T,W}(), typemin(T), typemax(T))
    BasicUvStatistics{T,W}(mean::OnlineUvMean{T}, var::OnlineUvVar{T,W}, maximum::T, minimum::T) where {T<:Real,W} =
        new(mean, var, maximum, minimum)
end

function Base.empty!(stats::BasicUvStatistics{T,W}) where {T,W}
    empty!(stats.mean)
    empty!(stats.var)
    stats.maximum = typemin(T)
    stats.minimum = typemax(T)

    stats
end

@inline function Base.push!(stats::BasicUvStatistics{T,W}, data::Real, weight::Real = 1) where {T,W}
    nmaximum = stats.maximum
    nminimum = stats.minimum

    push!(stats.mean, data, weight)
    push!(stats.var, data, weight)
    nmaximum = max(stats.maximum, maximum(data))
    nminimum = min(stats.minimum, minimum(data))

    stats.maximum = nmaximum
    stats.minimum = nminimum

    stats
end

function Base.merge!(target::BasicUvStatistics, others::BasicUvStatistics...)
    t_mean = target.mean
    t_var = target.var
    t_maximum = target.maximum
    t_minimum = target.minimum

    for x in others
        t_mean = merge!(t_mean, x.mean)
        t_var = merge!(t_var, x.var)
        t_maximum = max(t_maximum, x.maximum)
        t_minimum = min(t_minimum, x.minimum)
    end

    target.mean = t_mean
    target.var = t_var
    target.maximum = t_maximum
    target.minimum = t_minimum

    target
end

function reweight_relative!(stats::BasicUvStatistics{T,W}, reweighting_factor::Real) where {T,W}
    reweight_relative!(stats.mean, reweighting_factor)
    reweight_relative!(stats.var, reweighting_factor)

    stats
end

const OnlineUvStatistic = Union{BAT.OnlineUvMean, BAT.OnlineUvVar, BAT.BasicUvStatistics}

@inline function Base.append!(stats::OnlineUvStatistic, data::AbstractVector{<:Real})
    push!.(Ref(stats), data)
    stats
end

@inline function Base.append!(stats::OnlineUvStatistic, data::AbstractVector{<:Real}, weights::AbstractVector{<:Real})
    push!.(Ref(stats), data, weights)
    stats
end

Base.merge(target::S, others::S...) where{S <: OnlineUvStatistic} = merge!(deepcopy(target), others...)