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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 47d5b09c-c213-11eb-2472-7b23859c3caa
begin
	let
		using Pkg
		Pkg.activate(mktempdir())
		Pkg.Registry.update()
		Pkg.add("DICOM")
		Pkg.add("NIfTI")
		Pkg.add("Plots")
		Pkg.add(Pkg.PackageSpec(url="https://github.com/Dale-Black/DICOMUtils.jl"))
		Pkg.add("Images")
	end

	using PlutoUI
	using DICOM
	using NIfTI
	using Plots
	using DICOMUtils
	using LinearAlgebra
	using Images
end

# ╔═╡ e991c00b-974b-442c-93c6-55a0b12b27ed
TableOfContents()

# ╔═╡ 718b9504-f15f-41d1-b0ea-3346f17b4cca
md"""
## Load DICOMs
"""

# ╔═╡ a45f05b9-08d6-4277-9c9c-205f60b43391
volpath = "/Users/daleblack/Desktop/80.0";

# ╔═╡ 5e69faec-3f89-4f65-aab1-9429ad1e141b
volseg1 = "/Users/daleblack/Desktop/HEL_SLICER_SEG_0/S_1.20.nii";

# ╔═╡ 0399cca4-b251-451e-bf8c-96ccd207a1df
volseg2 = "/Users/daleblack/Desktop/HEL_SLICER_SEG_0/M_3.0.nii";

# ╔═╡ 20b7a9f9-4308-427f-995f-1f29f7f6df23
md"""
## Check Orientation
"""

# ╔═╡ e36c897d-07ae-4f61-8a5d-6dbdc651431a
vol = dcmdir_parse(volpath)

# ╔═╡ 8fe1a81f-be2d-4254-8029-aca2dc07a609
seg1 = niread(volseg1);

# ╔═╡ 6aea2d01-b4d5-442a-aafa-fec218cc4123
seg2 = niread(volseg2);

# ╔═╡ 644c1446-4206-4857-8aa3-ca87e0a6ff99
volraw = DICOMUtils.load_dcm_array(vol);

# ╔═╡ d729e367-61e4-4ade-b898-d66631443236
seg1raw = copy(seg1.raw);

# ╔═╡ f883ae8f-5761-471e-b80f-ec55a1fca07a
seg2raw = copy(seg2.raw);

# ╔═╡ 7d4ddb3d-0d5e-433c-a8f1-6ff77cdfd1c5
md"""
### Without Reorientation
"""

# ╔═╡ 5d6edfe9-9d47-4c8b-ac6b-ebc35e3b9915
@bind a Slider(1:size(volraw)[3], default=155, show_value=true)

# ╔═╡ b198dedc-0e5e-4c69-a755-387de3997d09
Plots.heatmap(volraw[:,:,a], c=:grays)

# ╔═╡ 4d49321a-7dbe-428b-91c4-c9561055ce0b
Plots.heatmap(seg1raw[:,:,a], c=:grays)

# ╔═╡ 55744a6c-a224-42cc-b590-bc7ec4614bc7
md"""
### Superimpose
"""

# ╔═╡ d4bc98e5-e92b-4b24-ad36-081b9c0d8507
md"""
## With Reorientation
"""

# ╔═╡ 670ec570-8a19-479d-a619-873c8fe9457a
volraw_ornt, affvol, new_affvol = DICOMUtils.orientation(volraw, (("L", "A", "S")));

# ╔═╡ 3bb89d52-da0e-42d2-9368-41c576b3ef11
seg1raw_ornt, aff1, new_aff1 = DICOMUtils.orientation(seg1raw, (("L", "A", "S")));

# ╔═╡ 16da0125-d6ad-45d8-ba30-431ff2fdc7ed
Plots.heatmap(volraw_ornt[:,:,a], c=:grays)

# ╔═╡ 8c49c4d7-8864-4a5d-9797-3e2bacb2f514
Plots.heatmap(seg1raw_ornt[:,:,a], c=:grays)

# ╔═╡ bb387fef-0034-467f-8219-232a5942c264
md"""
## Calcium Scoring
"""

# ╔═╡ 64599dbe-4f10-476b-921b-63872d7482d1
begin
	ii1 = [1 2 0 1
		   1 5 0 1]
	ii2 = [1 2 0 1
	   	   1 5 0 1]
	img = cat(ii1, ii2, dims=3)
end

# ╔═╡ 228e1866-5d50-41a6-b5f3-cce2d49cdfee
spc = [0.625, 0.625, 0.5]

# ╔═╡ 1252701e-6812-49e4-9282-7c51dd92f8cc
begin
	mm1 = [1 0 0 1
		   1 1 0 1]
	mm2 = [1 0 0 1
	   	   1 1 0 1]
	msk = cat(mm1, mm2, dims=3)
end

# ╔═╡ 288e776b-04ae-479f-8ae2-b3232037fe87
function compute_calcium_scores(image, spacing, mask, min_vol=nothing, max_vol=nothing)
	binary_mask = mask
	voxel_volume = spacing[1] * spacing[2] * spacing[3]

	agatston_score = 0
	calcium_volume = 0

	# Find individual lesions (in 3D) so that we can discard too small or too large lesions
	lesion_map = Images.label_components(binary_mask)
	n_lesions = length(unique(lesion_map))

	for lesion = 1:(n_lesions - 1)
		lesion_mask = map(x -> x == lesion, lesion_map)

		# Ignore too small or too large lesions
		lesion_volume = count(lesion_mask) * voxel_volume
		if ((min_vol != nothing) && (lesion_volume < min_vol))
			continue
		end
		if ((max_vol != nothing) && (lesion_volume > max_vol))
			continue
		end

		calcium_volume += lesion_volume

		# Calculate Agatston score for this lesion
		slices = sort(unique(lesion_mask)) .+ 1
		for z = 1:slices[end]
			fragment_mask = lesion_mask[z,:,:]
			n_pixels = count(fragment_mask)
			maximum_intensity = maximum(image[z,:,:][fragment_mask])
            if maximum_intensity < 200
                coefficient = 1
			elseif maximum_intensity < 300
				coefficient = 2
			elseif maximum_intensity < 400
				coefficient = 3
			else
				coefficient = 4
			end
			agatston_score += coefficient * n_pixels
		end
	end
	agatston_score *= spacing[1] / 3.0 * spacing[2] * spacing[3]
	return agatston_score, calcium_volume
end

# ╔═╡ ea475810-4853-4f6f-97bd-d3faf1b646d6
compute_calcium_scores(img, spc, msk)

# ╔═╡ 165084b0-4ebe-4c28-bdae-24278d7497ae
compute_calcium_scores(img, spc, msk, 0)

# ╔═╡ Cell order:
# ╠═47d5b09c-c213-11eb-2472-7b23859c3caa
# ╠═e991c00b-974b-442c-93c6-55a0b12b27ed
# ╟─718b9504-f15f-41d1-b0ea-3346f17b4cca
# ╠═a45f05b9-08d6-4277-9c9c-205f60b43391
# ╠═5e69faec-3f89-4f65-aab1-9429ad1e141b
# ╠═0399cca4-b251-451e-bf8c-96ccd207a1df
# ╟─20b7a9f9-4308-427f-995f-1f29f7f6df23
# ╠═e36c897d-07ae-4f61-8a5d-6dbdc651431a
# ╠═8fe1a81f-be2d-4254-8029-aca2dc07a609
# ╠═6aea2d01-b4d5-442a-aafa-fec218cc4123
# ╠═644c1446-4206-4857-8aa3-ca87e0a6ff99
# ╠═d729e367-61e4-4ade-b898-d66631443236
# ╠═f883ae8f-5761-471e-b80f-ec55a1fca07a
# ╟─7d4ddb3d-0d5e-433c-a8f1-6ff77cdfd1c5
# ╠═5d6edfe9-9d47-4c8b-ac6b-ebc35e3b9915
# ╠═b198dedc-0e5e-4c69-a755-387de3997d09
# ╠═4d49321a-7dbe-428b-91c4-c9561055ce0b
# ╟─55744a6c-a224-42cc-b590-bc7ec4614bc7
# ╟─d4bc98e5-e92b-4b24-ad36-081b9c0d8507
# ╠═670ec570-8a19-479d-a619-873c8fe9457a
# ╠═3bb89d52-da0e-42d2-9368-41c576b3ef11
# ╠═16da0125-d6ad-45d8-ba30-431ff2fdc7ed
# ╠═8c49c4d7-8864-4a5d-9797-3e2bacb2f514
# ╟─bb387fef-0034-467f-8219-232a5942c264
# ╠═64599dbe-4f10-476b-921b-63872d7482d1
# ╠═228e1866-5d50-41a6-b5f3-cce2d49cdfee
# ╠═1252701e-6812-49e4-9282-7c51dd92f8cc
# ╠═288e776b-04ae-479f-8ae2-b3232037fe87
# ╠═ea475810-4853-4f6f-97bd-d3faf1b646d6
# ╠═165084b0-4ebe-4c28-bdae-24278d7497ae


text:
module AVFilters
  include(joinpath(dirname(@__FILE__),"..","..","..","init.jl"))
  w(f) = joinpath(avfilter_dir, f)

  include(w("LIBAVFILTER.jl"))
end

ext: jl

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
"""
    integrate_gauss_legendre(f, a::Arb, b::Arb)

Compute an enclosure of the integral of `f` on the interval `[a, b]`
using a Gauss-Legendre quadrature of order 2.

The approximation is given by
```
(b - a) / 2 * (f(x₁) + f(x₂))
```
where
```
x₁ = (b - a) / 2 * sqrt(3) / 3 + (b + a) / 2
x₂ = -(b - a) / 2 * sqrt(3) / 3 + (b + a) / 2
```
are quadrature nodes. The remainder term is enclosed by
```
(b - a)^5 / 4320 * d4f
```
where `d4f` is an enclosure of the fourth derivative of `f` on the
interval `[a, b]`.

If it happens to be that `d4f` is not finite, for example because `f`
is not differentiable on the interval, it instead computes a naive
enclosure of the integral using the enclosure of `f` evaluated on the
interval `[a, b]`.

The function `f` should support evaluation on both `Arb` and
`ArbSeries` and should return an enclosure of the result in both
cases. For `ArbSeries` it is enough that it supports evaluation up to
degree 4 and only the zeroth and fourth order term are ever used.

- **TODO:** Optimize for performance
- **TODO:** Add option for computing integral of absolute value?
- **TODO:** Allow giving an separate function for computing derivative?
"""
function integrate_gauss_legendre(f, a::Arb, b::Arb)
    bma = b - a

    # x_series = ArbSeries((union(a, b), 1), degree = 4)
    x_series = ArbSeries(degree = 4, prec = precision(bma))
    # We need to do it in this order so that the degree of the
    # polynomial is correct when we set the constant coefficient.
    x_series[1] = 1
    Arblib.union!(Arblib.ref(x_series, 0), a, b)

    f_series = f(x_series)

    # If the fourth derivative of f is not finite compute an enclosure
    # of the integral using the enclosure of f on the interval [a, b].
    isfinite(Arblib.ref(f_series, 4)) ||
        return Arblib.mul!(bma, bma, Arblib.ref(f_series, 0))

    # remainder = (b - a)^5 / 4320 * f_series[4] * factorial(4)
    remainder = Arblib.pow!(zero(bma), bma, UInt(5))
    Arblib.mul!(remainder, remainder, Arblib.ref(f_series, 4))
    Arblib.div!(remainder, remainder, 180)

    # v = (b - a) / 2 * sqrt(Arb(3)) / 3
    v = Arblib.sqrt!(zero(bma), UInt(3))
    Arblib.mul!(v, v, bma)
    Arblib.div!(v, v, 6)

    # Set x₁ = (b + a) / 2 to begin with
    x₁ = b + a
    Arblib.mul_2exp!(x₁, x₁, -1)

    # x₂ = (b + a) / 2 - v = x₁ - v
    x₂ = Arblib.sub!(zero(bma), x₁, v)

    # x₁ = (b + a) / 2 - v = x₁ - v
    Arblib.add!(x₁, x₁, v)

    # approximation = (b - a) / 2 * (f(x₁) + f(x₂))
    approximation = f(x₁)
    Arblib.add!(approximation, approximation, f(x₂))
    Arblib.mul!(approximation, approximation, bma)
    Arblib.mul_2exp!(approximation, approximation, -1)

    # Return approximation + remainder
    return Arblib.add!(approximation, approximation, remainder)
end

"""
    integrate(f, a::Arb, b::Arb; atol, rtol, depth_start, maxevals, depth, verbose)

Compute an enclosure of the integral of `f on the interval `[a, b]`.

It uses a Guass-Legendre quadrature of order 2 through
[`integrate_gauss_legendre`](@ref) together with bisection of the
interval. On each subinterval the integral is computed using
[`integrate_gauss_legendre`](@ref) and it checks if the result
satisfies the required tolerance, if it doesn't the interval is
bisected. For more details about the integration method see
[`integrate_gauss_legendre`](@ref).

Notice that the tolerance is checked on each **subinterval** and not
on the interval as a whole. This means that even if it finishes before
reaching the maximum depth or the maximum number of evaluations the
result as a whole might not satisfy the given tolerance.

The argument `depth_start` bisect the interval using
[`bisect_interval_recursive`](@ref) before starting to compute the
integral. This can be useful if it is known beforehand that a certain
number of bisections will be necessary before the enclosures get good
enough. It defaults to `0` which corresponds to not bisecting the
interval at all before starting.

The arguments `maxevals` and `depth` can be used to limit the number
of function evaluations and the number of bisections of the interval
respectively.

If `verbose = true` then output information about the process.
"""
function integrate(
    f,
    a::Arb,
    b::Arb;
    atol = 0,
    rtol = sqrt(eps(one(a))),
    depth_start::Integer = 0,
    maxevals::Integer = 1000,
    depth::Integer = 20,
    verbose = false,
)
    isfinite(a) && isfinite(b) || return Arblib.indeterminate!(zero(a))

    intervals = bisect_interval_recursive(a, b, depth_start)

    integral = zero(a)

    verbose && @info "iteration: $(lpad(0, 2)), " *
          "starting intervals: $(lpad(length(intervals), 4))"

    iterations = 0
    evals = 0
    while !isempty(intervals)
        iterations += 1
        evals += length(intervals)

        # If we have reached the maximum number of iterations or
        # evaluations then we don't split any further.
        finished = iterations >= depth || evals >= maxevals

        # The number of intervals that had to be split further. The
        # only reason we store this in a separate variable is to
        # correctly count the number during the last iteration when
        # finished = true.
        remaining_intervals = 0

        next_intervals = empty(intervals)

        for (a, b) in intervals
            integral_part = integrate_gauss_legendre(f, a, b)

            if check_tolerance(integral_part; atol, rtol)
                integral += integral_part
            elseif finished
                remaining_intervals += 1
                integral += integral_part
            else
                remaining_intervals += 1
                push!(next_intervals, bisect_interval(a, b)...)
            end
        end

        intervals = next_intervals

        verbose && @info "iteration: $(lpad(iterations, 2)), " *
              "remaining intervals: $(lpad(remaining_intervals, 3)), " *
              "current integral: $integral"
    end

    return integral
end


text:
nextfloat(p::Sonum8) = reinterpret(Sonum8,reinterpret(UInt8,p)+one(UInt8))
nextfloat(p::Sonum16) = reinterpret(Sonum16,reinterpret(UInt16,p)+one(UInt16))

prevfloat(p::Sonum8) = reinterpret(Sonum8,reinterpret(UInt8,p)-one(UInt8))
prevfloat(p::Sonum16) = reinterpret(Sonum16,reinterpret(UInt16,p)-one(UInt16))

function -(x::Sonum8)
    if UInt8(x) == 0x00 || UInt8(x) == 0x80 # don't change sign for 0 and NaR
        return x
    else    # subtracting from 0x00 (two's complement def for neg)
        return Sonum8(0x00 - UInt8(x))
    end
end

function -(x::Sonum16)
    if UInt16(x) == 0x0000 || UInt16(x) == 0x8000 # don't change sign for 0 and NaR
        return x
    else    # subtracting from 0x0000 (two's complement def for neg)
        return Sonum16(0x0000 - UInt16(x))
    end
end

abs(x::Sonum8) = UInt8(x) > 0x80 ? Sonum8(0x00 - UInt8(x)) : x
abs(x::Sonum16) = UInt16(x) > 0x8000 ? Sonum8(0x0000 - UInt16(x)) : x

signbit(x::Sonum8) = UInt8(x) >= 0x80 ? true : false
signbit(x::Sonum16) = UInt16(x) >= 0x8000 ? true : false

function sign(p::T) where {T <: AbstractSonum}
    if signbit(p)       # negative and infinity case
        if isfinite(p)  # negative
            return minusone(T)
        else            # infinity
            return zero(T)
        end
    else                # positive and zero case
        if iszero(p)    # zero
            return zero(T)
        else            # positive
            return one(T)
        end
    end
end


text:
println("Caesar Cipher \n")
println("Commands: ")
println("\t - [1] Encrypt ")
println("\t - [2] Decrypt ")

println("\nWhat do you want to do?")
print("Command: ")
decision = parse(Int8, readline())

alphabet = split("abcdefghijklmnopqrstuvwxyz", "")

if decision != 1 && decision != 2
    println("Invalid Command!")

elseif decision == 1
    println("\nEnter the text you want to encrypt")
    print("Text: ")
    text = readline()
    text = split(text, "")
    encrypted_array = Int8[]
    encrypted_text = Any[]

    print("Key: ")
    key = parse(Int8,readline())

    for i in 1:length(text)

        index = findfirst(x-> x == text[i], alphabet)

        if text[i] == " "
            append!(encrypted_text, " ")
        else
            append!(encrypted_array, (index + key) % 26)
        end

    end

    for i in 1:length(encrypted_array)
        if encrypted_array[i] == 0
            index = 26
        else
            index = i
        end

        append!(encrypted_text, alphabet[index])
    end

    println("\nEncrypted Text: $(encrypted_text)")
    println(join(encrypted_text))
    println("Key $key")

end

    # for i in 1:length(text)
    #     text_array.append(char[])

max_stars_repo_name: Brymer-Meneses/starter-projects

max_issues_repo_name: brymer-meneses/starter-projects

max_forks_repo_name: brymer-meneses/starter-projects


text:
function score(x, y)
    r = sqrt(x^2 + y^2)
    if r <= 1
        return 10
    elseif r <= 5
        return 5
    elseif r <= 10
        return 1
    else
        return 0
    end
end

max_stars_repo_name: aadimator/exercism

max_issues_repo_name: aadimator/exercism

max_forks_repo_name: aadimator/exercism


text: @emit "model" ["vertices"=>[[0,0,0] [50,0,0] [0,0,50] [0,50,0]],"faces"=>[[0,1,2] [1,3,2] [0,3,1] [0,2,3]]]


text:
function run_powermodels_powerflow(json_path)
    pm = load_pm_from_json(json_path)
    model = get_model(pm["pm_model"])

    solver = get_solver(pm["pm_solver"], pm["pm_nl_solver"], pm["pm_mip_solver"],
    pm["pm_log_level"], pm["pm_time_limit"], pm["pm_nl_time_limit"], pm["pm_mip_time_limit"])

    result = _PM.run_pf(pm, model, solver,
                                setting = Dict("output" => Dict("branch_flows" => true)))
    return result
end

function run_powermodels_opf(json_path)
    pm = load_pm_from_json(json_path)
    model = get_model(pm["pm_model"])

    solver = get_solver(pm["pm_solver"], pm["pm_nl_solver"], pm["pm_mip_solver"],
    pm["pm_log_level"], pm["pm_time_limit"], pm["pm_nl_time_limit"], pm["pm_mip_time_limit"])

    if haskey(pm["branch"]["1"],"c_rating_a")
        for (key, value) in pm["gen"]
           value["pmax"] *= 0.01
           value["qmax"] *= 0.01
           value["qmin"] *= 0.01
           value["pg"] *= 0.01
           value["qg"] *= 0.01
           value["cost"] *= 100
        end

        for (key, value) in pm["branch"]
           value["c_rating_a"] *= 0.01
        end

        for (key, value) in pm["load"]
           value["pd"] *= 0.01
           value["qd"] *= 0.01
        end

        result = PowerModels._run_opf_cl(pm, model, solver,
                                        setting = Dict("output" => Dict("branch_flows" => true)))
    else
        result = PowerModels.run_opf(pm, model, solver,
                                        setting = Dict("output" => Dict("branch_flows" => true)))
    end
    return result
end

function run_powermodels_custom(json_path)
    pm = load_pm_from_json(json_path)

    _PM.correct_network_data!(pm)

    ipopt_solver = JuMP.optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)

    result = _PM.run_ac_opf(pm, ipopt_solver,
                            setting = Dict("output" => Dict("branch_flows" => true)))
    return result
end

function run_powermodels_tnep(json_path)
    pm = load_pm_from_json(json_path)
    model = get_model(pm["pm_model"])

    solver = get_solver(pm["pm_solver"], pm["pm_nl_solver"], pm["pm_mip_solver"],
    pm["pm_log_level"], pm["pm_time_limit"], pm["pm_nl_time_limit"], pm["pm_mip_time_limit"])

    result = _PM.run_tnep(pm, model, solver,
                        setting = Dict("output" => Dict("branch_flows" => true)))
    return result
end

function run_powermodels_ots(json_path)
    pm = load_pm_from_json(json_path)
    model = get_model(pm["pm_model"])

    solver = get_solver(pm["pm_solver"], pm["pm_nl_solver"], pm["pm_mip_solver"],
    pm["pm_log_level"], pm["pm_time_limit"], pm["pm_nl_time_limit"], pm["pm_mip_time_limit"])

    result = _PM.run_ots(pm, model, solver,
                        setting = Dict("output" => Dict("branch_flows" => true)))
    return result
end


text:
module SourceWalk

using Reexport

@reexport using MacroTools
using MacroTools: @forward, striplines, postwalk
using DataStructures

export sourcewalk

include("cstparser/CSTParser.jl")

module hack; using CSTParser; end

using .CSTParser

include("patch/diff.jl")
include("patch/cst.jl")

end # module


text:
module TestExamplesDGMulti2D

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "dgmulti_2d")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

@testset "DGMulti 2D" begin
  @trixi_testset "elixir_euler_weakform.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4,4),
      l2 = [0.0013463150267220328, 0.0014235793662296975, 0.0014235793662300024, 0.00472191786388071],
      linf = [0.0015248303482329195, 0.0020707330926952316, 0.002070733092696342, 0.004913455679613321],
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (SBP)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4,4),
      approximation_type = SBP(),
      l2 = [0.007465016439747669, 0.005297423547392931, 0.005297423547403158, 0.01470161132498598],
      linf = [0.021489935389522374, 0.013528869419211276, 0.013528869418737433, 0.03269072432621112],
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (Quadrilateral elements)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4,4),
      element_type = Quad(),
      l2 = [0.000290967383489527, 0.00028113809346926776, 0.0002811380934695505, 0.00102004771420306],
      linf = [0.0004970344840584673, 0.00040590009306518127, 0.0004059000930629608, 0.0014247732095258314],
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (EC) " begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4,4),
      volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
      surface_integral = SurfaceIntegralWeakForm(flux_ranocha),
      l2 = [0.008110785816182155, 0.0074686552093368745, 0.007468655209336097, 0.015986513837074563],
      linf = [0.01230954687917274, 0.013884805356942254, 0.013884805356973784, 0.040387377818142056],
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (SBP, EC)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4,4),
      volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
      surface_integral = SurfaceIntegralWeakForm(flux_ranocha),
      approximation_type = SBP(),
      l2 = [0.01285864081726596, 0.010650165503847099, 0.01065016550381281, 0.026286162111579015],
      linf = [0.037333313274372504, 0.05308320130762212, 0.05308320130841948, 0.13378665881805185],
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (Quadrilateral elements, SBP, EC)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4,4),
      element_type = Quad(),
      volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
      surface_integral = SurfaceIntegralWeakForm(flux_ranocha),
      approximation_type = SBP(),
      l2 = [0.0029319624187308896, 0.0030625695968579886, 0.003062569596855081, 0.006843948320775483],
      linf = [0.013700713240587747, 0.012810144682950497, 0.01281014468295072, 0.022412654330661308],
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (convergence)" begin
    mean_convergence = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"), 2)
    @test isapprox(mean_convergence[:l2], [4.249875508800025, 4.133727008051228, 4.133727007601049, 4.086238794189699], rtol=0.05)
  end

  @trixi_testset "elixir_euler_weakform_periodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform_periodic.jl"),
      l2 = [0.0014986508075708323, 0.001528523420746786, 0.0015285234207473158, 0.004846505183839211],
      linf = [0.0015062108658376872, 0.0019373508504645365, 0.0019373508504538783, 0.004742686826709086]
    )
  end

  @trixi_testset "elixir_euler_triangulate_pkg_mesh.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_triangulate_pkg_mesh.jl"),
      l2 = [4.664661209491976e-6, 3.7033509525940745e-6, 4.794877426562555e-6, 1.2682723101532175e-5],
      linf = [2.5099852761334418e-5, 2.2683684021362893e-5, 2.6180448559287584e-5, 5.5752932611508044e-5]
    )
  end

  @trixi_testset "elixir_euler_kelvin_helmholtz_instability.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability.jl"),
      cells_per_dimension = (8,8), tspan = (0.0, 0.2),
      l2 = [0.11140371194266245, 0.06598148257331912, 0.10448950567476646, 0.1602319547232773],
      linf = [0.2403274017090733, 0.1659961476843728, 0.12354683515032569, 0.2691592411995334],
    )
  end

  @trixi_testset "elixir_euler_rayleigh_taylor_instability.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_rayleigh_taylor_instability.jl"),
      cells_per_dimension = (8,8), tspan = (0.0, 0.2),
      l2 = [0.03548329484912729, 0.002591414376082683, 0.006916327792623457, 0.016235069002818153],
      linf = [0.4783963902824797, 0.022527207050681054, 0.040307056293369226, 0.0852365428206836],
    )
  end

  @trixi_testset "elixir_euler_brown_minion_vortex.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_brown_minion_vortex.jl"),
      num_cells_per_dimension = 4, tspan = (0.0, 0.1),
      l2 = [0.006818116685388498, 0.02123358215096057, 0.009928432022103199, 0.15364518920518527],
      linf = [0.015838346753841548, 0.09449772987646823, 0.021276761322582504, 0.3470925282340751],
    )
  end
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn isdir(outdir) && rm(outdir, recursive=true)

end # module


text:
module DataBagsTests
using Test, DataBags, Dates
using DataBags: AbstractDataBag, propertyname, withspecificvaluetype

function samevalues(A::AbstractArray, B::AbstractArray)
    @assert has_standard_indexing(A, B)
    @assert size(A) == size(B)
    for i ∈ eachindex(A, B)
        A[i] != B[i] && return false
    end
    return true
end

function maxabsdif(A::AbstractArray, B::AbstractArray)
    @assert has_standard_indexing(A, B)
    @assert size(A) == size(B)
    T = promote_type(eltype(A), eltype(B))
    res = zero(T)
    for i ∈ eachindex(A, B)
        res = max(res, abs(A[i] - B[i]))
    end
    return res
end

keywords(args...; kwds...) = kwds
arguments(args...; kwds...) = args

# UnfinishedDataBag does not extend DataBags.contents()
struct UnfinishedDataBag{K,V,D<:AbstractDict{K,V}} <: AbstractDataBag{K,V,D}
    data::D
    UnfinishedDataBag{K,V,D}(data::D) where {K,V,D<:AbstractDict{K,V}} =
        new{K,V,D}(data)
end
UnfinishedDataBag(data::D) where {K,V,D<:AbstractDict{K,V}} =
    UnfinishedDataBag{K,V,D}(data)

# Create our own data-bag with provided macro.
DataBags.@newtype SmallDB

# Check show().
show(stdout, MIME"text/plain"(),
     DataBag(a=21, b=π, c="hello", d=-1.8:0.1:14, e=rand(5), date=now()))
println()

@testset "DataBags" begin
    # Check missing specialization of DataBags.contents().
    Q = UnfinishedDataBag(Dict{Symbol,Any}())
    @test_throws ErrorException Q.units

    # Check impossible indexing of weird key type.
    R = DataBag(Dict{Int,Int}())
    @test_throws ErrorException propertyname(typeof(R), Symbol(42))

    # DataBag with string keys.
    D1 = Dict("units" => "km", "Δx" => 0.20, "Δy" => 0.15)
    A1 = DataBag(D1)
    T1 = typeof(A1)
    B1 = wrap(DataBag, D1)
    @test length(A1) == length(D1)
    @test keytype(T1) == keytype(D1)
    @test valtype(T1) == valtype(D1)
    @test keytype(A1) == keytype(D1)
    @test valtype(A1) == valtype(D1)
    for (k,v) in zip(keys(A1), values(A1))
        @test D1[k] == v
    end
    for (k,v) in pairs(A1)
        @test D1[k] == v
    end
    # A1 and D1 should be different, B1 abd D1 should be the same:
    @test A1.units == A1["units"] == D1["units"]
    A1.units = 50
    @test contents(D1) === D1
    @test contents(B1) === D1
    @test contents(A1) !== D1
    @test propertyname(T1, :gizmo) == "gizmo"
    @test A1.units == A1["units"] != D1["units"]
    @test B1.units == B1["units"] == D1["units"]
    @test haskey(A1, "Δx") == true
    delete!(A1, "Δx")
    @test haskey(A1, "Δx") == false
    merge!(A1, D1)
    @test haskey(A1, "Δx") == true
    @test A1.Δx == A1["Δx"] == D1["Δx"]
    @test getkey(A1, :foo, π) == π
    @test getkey(A1, "units", π) == "units"
    @test get(A1, :foo, π) == π
    @test get(A1, "units", π) == "km"
    @test (!haskey(A1, "foo") && get!(A1, "foo", π) == π &&
           haskey(A1, "foo") && A1.foo == π)
    @test pop!(A1, "foo", π) == π && !haskey(A1, "foo")
    @test pop!(A1, "foo", :bar) == :bar
    @test pop!(A1, "units") == "km"
    @test_throws KeyError pop!(A1, "units")
    C1 = convert(DataBag, A1) # should yield A1
    C1.bar = "foo"
    @test haskey(A1, "bar") && A1.bar == "foo"
    delete!(delete!(A1, "units"), "bar")
    C1 = convert(DataBag{String,Float32}, A1) # should not yield A1
    C1.bar = π
    @test !haskey(A1, "bar") && C1.bar ≈ π

    # DataBag with symbolic keys.
    D2 = Dict(:units => "µm", :Δx => 0.20, :Δy => 0.15)
    D3 = Dict(:name => "Mr. Doe", :Δx => 0.30, :Δy => 0.20)
    A2 = DataBag(D2)
    T2 = typeof(A2)
    @test length(A2) == length(D2)
    @test keytype(T2) == keytype(D2)
    @test keytype(A2) == keytype(D2)
    @test keytype(A2) == keytype(D2)
    @test valtype(A2) == valtype(D2)
    @test propertyname(T2, :gizmo) == :gizmo
    @test propertynames(A2) == keys(A2)
    @test A2.units == A2[:units] == D2[:units]
    A2.units = 50
    @test A2.units == A2[:units] != D2[:units]
    @test haskey(A2, :Δx) == true
    delete!(A2, :Δx)
    @test haskey(A2, :Δx) == false
    merge!(A2, D2)
    @test haskey(A2, :Δx) == true
    @test A2.units == A2[:units] == D2[:units]
    A3 = empty(A2)
    @test length(A3) == 0 && length(A2) > 0
    @test keytype(A3) == keytype(A2)
    @test valtype(A3) == valtype(A2)
    A4 = empty(A2, String)
    @test length(A4) == 0
    @test keytype(A4) == String
    @test valtype(A4) == valtype(A2)
    A5 = empty(A2, String, Int32)
    @test length(A5) == 0
    @test keytype(A5) == String
    @test valtype(A5) == Int32
    empty!(A2)
    @test length(A2) == 0
    merge!(A2, D2)
    merge!(A3, D3)
    A6 = empty(A3)
    merge!(A6, D2, D3)
    @test length(A6) == 4
    @test A6.units == A2.units
    @test A6.name == A3.name
    @test A6.Δx == A3.Δx && A6.Δy == A3.Δy
    merge!(empty!(A6), A3, A2)
    @test length(A6) == 4
    @test A6.units == A2.units
    @test A6.name == A3.name
    @test A6.Δx == A2.Δx && A6.Δy == A2.Δy

    # Test different constructors and combination of values.
    B0 = DataBag()
    @test keytype(B0) == Symbol && valtype(B0) == Any
    B1 = DataBag( a  =  1,  b  =  2,  c  =  3,  d  =  "µm",  e  =  :dif)
    @test keytype(B1) == Symbol && valtype(B1) == Any
    B2 = DataBag(:a  => 1, :b  => 2, :c  => 3, :d  => "µm", :e  => :dif)
    @test keytype(B2) == Symbol && valtype(B2) == Any
    B3 = DataBag("a" => 1, "b" => 2, "c" => 3, "d" => "µm", "e" => :dif)
    @test keytype(B3) == String && valtype(B3) == Any
    N1 = Dict{Symbol,Number}() # to store the numerical values of A1
    N2 = Dict{Symbol,Number}() # to store computed values
    for (k,v) in B1
        if isa(v, Number)
            N1[k] = v
            N2[k] = 1 - v
        end
    end
    B4 = DataBag(N1)
    @test keytype(B4) == keytype(N1) && valtype(B4) == Any
    B5 = DataBag{Symbol}(N2)
    @test keytype(B5) == Symbol && valtype(B5) == Any
    B6 = DataBag{Symbol,Float64}(N1)
    @test keytype(B6) == Symbol && valtype(B6) == Float64
    @test all(x -> x == 1, values(merge(+, B5, B4)))
    @test all(x -> x == 1, values(merge(+, N2, B4)))
    @test all(x -> x == 1, values(merge!(+, B6, B5)))

    # Explictely check wrap() and check that the 2 objects share the same data.
    W1 = wrap(DataBag, D1)
    @test keytype(W1) == keytype(D1) && valtype(W1) == valtype(D1)
    @test !haskey(W1, "newkey") && !haskey(D1, "newkey")
    W1.newkey = 12345
    @test haskey(W1, "newkey") && haskey(D1, "newkey")
    @test W1.newkey == pop!(D1, "newkey")
    @test !haskey(W1, "newkey") && !haskey(D1, "newkey")
    #
    W2 = wrap(DataBag, D2)
    @test keytype(W2) == keytype(D2) && valtype(W2) == valtype(D2)
    @test !haskey(W2, :newkey) && !haskey(D2, :newkey)
    W2.newkey = 12345
    @test haskey(W2, :newkey) && haskey(D2, :newkey)
    @test W2.newkey == pop!(D2, :newkey)
    @test !haskey(W2, :newkey) && !haskey(D2, :newkey)
    #
    W3 = wrap(DataBag{keytype(D2)}, D2)
    @test keytype(W3) == keytype(D2) && valtype(W3) == valtype(D2)
    @test !haskey(W3, :newkey) && !haskey(D2, :newkey)
    W3.newkey = 12345
    @test haskey(W3, :newkey) && haskey(D2, :newkey)
    @test W3.newkey == pop!(D2, :newkey)
    @test !haskey(W3, :newkey) && !haskey(D2, :newkey)
    #
    W4 = wrap(DataBag{keytype(N2),valtype(N2)}, N2)
    @test keytype(W4) == keytype(N2) && valtype(W4) == valtype(N2)
    @test !haskey(W4, :newkey) && !haskey(N2, :newkey)
    W4.newkey = 12345
    @test haskey(W4, :newkey) && haskey(N2, :newkey)
    @test W4.newkey == pop!(N2, :newkey)
    @test !haskey(W4, :newkey) && !haskey(N2, :newkey)
    #
    W5 = wrap(DataBag{keytype(N1),valtype(N1),typeof(N1)}, N1)
    @test keytype(W5) == keytype(N1) && valtype(W5) == valtype(N1)
    @test !haskey(W5, :newkey) && !haskey(N1, :newkey)
    W5.newkey = 12345
    @test haskey(W5, :newkey) && haskey(N1, :newkey)
    @test W5.newkey == pop!(N1, :newkey)
    @test !haskey(W5, :newkey) && !haskey(N1, :newkey)

    # Check contents(Dict{...}, ...)
    @test isa(contents(Dict{Any,Any},     D1), Dict{String,Any})
    @test isa(contents(Dict{String,Any},  D1), Dict{String,Any})
    @test isa(contents(Dict{Any,Any},     N1), Dict{Symbol,Any})
    @test isa(contents(Dict{Any,Number},  N1), Dict{Symbol,Number})
    @test isa(contents(Dict{Any,Integer}, N1), Dict{Symbol,Integer})
    @test isa(contents(Dict{Any,Int16},   N1), Dict{Symbol,Int16})
    kwds = keywords(i8=Int8(1), i16=Int16(2), i32=Int32(3), i64=Int64(4))
    @test isa(contents(Dict{Any,Any};      kwds...), Dict{Symbol,Any})
    @test isa(contents(Dict{Symbol,Any};   kwds...), Dict{Symbol,Any})
    @test isa(contents(Dict{Any,Integer};  kwds...), Dict{Symbol,Integer})
    @test isa(contents(Dict{Symbol,Int16}; kwds...), Dict{Symbol,Int16})
    @test_throws MethodError contents(Dict{String,Int16}; kwds...)
    args1 = (:i8=>Int8(1), :i16=>Int16(2), :i32=>Int32(3), :i64=>Int64(4))
    @test isa(contents(Dict{Any,Any},      args1...), Dict{Symbol,Any})
    @test isa(contents(Dict{Symbol,Any},   args1...), Dict{Symbol,Any})
    @test isa(contents(Dict{Any,Integer},  args1...), Dict{Symbol,Integer})
    @test isa(contents(Dict{Symbol,Int16}, args1...), Dict{Symbol,Int16})
    args2 = map(kv -> String(kv[1]) => kv[2], args1)
    @test isa(contents(Dict{Any,Any},      args2...), Dict{String,Any})
    @test isa(contents(Dict{String,Any},   args2...), Dict{String,Any})
    @test isa(contents(Dict{Any,Integer},  args2...), Dict{String,Integer})
    @test isa(contents(Dict{String,Int16}, args2...), Dict{String,Int16})
    args3 = ("i8"=>Int8(1), :i16=>Int16(2), :"32"=>Int32(3), :i64=>Int64(4))
    @test isa(contents(Dict{Any,Any},      args3...), Dict{Any,Any})
    @test isa(contents(Dict{Any,Integer},  args3...), Dict{Any,Integer})
    @test isa(contents(Dict{Any,Int16},    args3...), Dict{Any,Int16})

    # Test data-bag type created by DataBags.@newtype.
    X = SmallDB(date=now(), value=π, arr=rand(3))
    @test isa(X, SmallDB)
    @test isa(X, DataBags.AbstractDataBag{Symbol,Any})
    @test isa(X, AbstractDict)
    @test length(X) == 3
    X.more = "some thing else"
    @test length(X) == 4
    merge!(X, D2)
    @test X.Δx == D2[:Δx]

    # Miscellaneaous.
    @test withspecificvaluetype(valtype(D2), D2) == D2
    @test withspecificvaluetype(Real, N2) == Dict{keytype(N2),Real}(N2)
end

end # module


text:
function Base.show(io::IO, m::SingleCellExpr )
  println(string( "And object of class " , typeof( m) , " with " , size( m.data)[2], " Cells and ", size(m.data)[1] , " Genes" )  )
  if ( typeof(m.raw) != Int64 )
    println(string( "The data has been normalized" ))
  end
  if ( typeof(m.zscored) != Int64 )
    println(string( "The data has been z-scored" ))
  end
end


text:
struct telescope
    f::Float64 # Focal Length, meters
    D::Float64 # Aperture, meters
    σ::Float64 # Sensor Noise, pixels
    p::Float64 # pixel pitech, m
    λ::Float64 # wavelength, m
    𝓂::AbstractArray{Int,1} # Image Size (m1 x m2), pixels
    σ_s::Float64 # Std Dev of signal derived from telescope parameters
    telescope(f,D,σ,p,λ,𝓂) = new(f,D,σ,p,λ,𝓂,0.45*λ*f/D/p)
end

struct target
    SNR::Float64 # Signal to Noise Ratio
    vel::Float64 # Target Velocity Magnitude, pixels/image
end

# --------
# A Detection is the result of extracting the point sources from an SSA image
struct Detection
    ID::Int
    Time::Float64
    Pd::Float64
    Centroid::SVector{2,Float64}
    CentroidUncertainty::SVector{2,Float64}
    Signal::SVector{1,Float64}
    SignalUncertainty::SVector{1,Float64}
    TrackingWindow::SVector{4,Float64}
end
# --------

#=struct DataPipe{T1, T2, T3}
    BackgroundMethod::T1
    ClusterMethod::T2
    Threshold::T3
end =#

# --------
# Object to hold information about a track (result of Kalman Filter)
struct Track
    Measurements::Array{Detection,1}
    State::Array{SMatrix{4,1,Float64},1}
    Covariance::Array{SMatrix{4,4,Float64},1}
    Residuals::Array{SMatrix{2,1,Float64},1}
end
# --------

# --------
# TrackHistory holds an array of Tracks to keep the history
struct TrackHistory
    Tracks::Array{Track,1}
end
# --------

# --------
# An AssociationMatrix is an TxS matrix of integers where each row, T
# is a unique track and each Column, S, is a time history containing
# the index of the Detection that is associated to that track at time, ts.
struct AssociationMatrix
    TOMHT::Array{Int,2}
end
# --------

# --------
# Clusters holds a single Track and all Detections within some measurement
# gating distance away (ex: mahalanobis distance)
struct Cluster
    Track::Track
    Detections::Array{Detection,1}
end
# --------

# --------
# HypothesisCost is the same size as AssociationMatrix and has a running value
# of the cost function for each Hypothesis. Col(1) = Cost(1), Col(2) = Cost(1) +
# Cost(2), Col(3) = Cost(2) + Cost(3) = Cost(1) + Cost(2) + Cost(3),... etc.
struct HypothesisCost
    CostMatrix::Array{Float64,2}
end
# --------

# --------
# Tuning object to hold all of the parameters I dont want to recompile to change
struct TuningParams
    ProcessNoise::Float64
    M::Int64
    ScanningWindow::Int64
    Pfa::Float64
    Inflator::Float64
    TrackingWindow::Int64
end
# --------
# --------
# Tuning object to hold all of the parameters I dont want to recompile to change
struct Frame
    Raw::Array{Float64,2}
    Background::Array{Float64,2}
    Detections::Vector{Detection}
    Tracks::Vector{Track}
    Associations::AssociationMatrix
    Cost::HypothesisCost
end
# --------

# --------
# Catalog - filled with all tracks that leave the scene and persist at the end of the processing
struct Catalog
    EndFrames::Vector{Int64} # Frame number where track ends
    Tracks::Vector{Track}    # Track
end
# --------

export telescope, target, Detection, Track, TrackHistory, AssociationMatrix, Cluster, HypothesisCost, TuningParams, Frame, Catalog


text:
using Documenter, Dimensionless

makedocs(sitename="Dimensionless.jl",
pages = ["Home" => "index.md",
        "example.md",
        "library.md"
        ],
format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true")
)

deploydocs(
    repo = "github.com/martinkosch/Dimensionless.jl.git",
)


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 9c5d8972-7485-11eb-25f7-155fbbd40d98
using Revise

# ╔═╡ 2c24aaca-7483-11eb-05a9-6f847b7d09b7
using Distributions

# ╔═╡ 42bb4640-7483-11eb-0661-677c524c62b0
using Plots

# ╔═╡ eeb0be1c-7483-11eb-341f-b1996e9a53a6
using OmegaMH

# ╔═╡ 9d2f2714-749b-11eb-1403-e3a045df62ca
using ReplicaExchange

# ╔═╡ ebd11c9a-7484-11eb-2941-c78e3e22d260
using Random

# ╔═╡ 0464d3ec-7489-11eb-1081-d761307263db
using Pkg; Pkg.add("PyPlot")

# ╔═╡ 94d6f4e6-74ac-11eb-0871-a31b2df5526d
using PlutoUI

# ╔═╡ 128974e2-7483-11eb-03a0-15b2abc75d66
md"# MCMC Test"

# ╔═╡ 318e4a98-7483-11eb-1916-c9630ecf07e2
md"First, let's setup our target density"

# ╔═╡ 27ced392-7483-11eb-3d36-7b8277daadff
dist = MvNormal([0.0, 0.0], [1.0, 1.0])

# ╔═╡ 7e4e66ce-7483-11eb-238e-d15279086f3b
target_density(x, y) = pdf(dist, [x, y])

# ╔═╡ 9a7dc4f2-7483-11eb-36be-251be000793e
target_density(0.2, 0.4)

# ╔═╡ 68100d38-7483-11eb-07d5-45547ad45143
Plots.surface(-4:0.1:4, -4:0.1:4, target_density, c = :blues)

# ╔═╡ f744509a-7483-11eb-3168-9bff07e7b245
md"Now let's sample from the distribution using Metropolis Hastings";

# ╔═╡ 09fe4874-7484-11eb-28be-9f5df0586c87
state_init = rand(2)

# ╔═╡ 20537202-7484-11eb-1d70-a1bcbaa3c265
target_logdensity(state) = logpdf(dist, state);

# ╔═╡ 3b0ac92e-7484-11eb-2e90-cf94eb2e9731
num_samples = 100001;

# ╔═╡ 4343d3b0-7484-11eb-388b-576f50c83335
function propose_and_logratio(rng, state)
	prop_dist = MvNormal(state, [0.1, 0.1])
	sample = rand(rng, prop_dist)
	(sample, 0)
end

# ╔═╡ e84726c8-7484-11eb-1512-c5cee451507a
rng = Random.MersenneTwister(0)

# ╔═╡ f8e56030-7484-11eb-2bb6-afa2f952b085
samples = OmegaMH.mh(rng, target_logdensity, num_samples, state_init, propose_and_logratio)

# ╔═╡ 747691fe-7487-11eb-292e-2975ef75e7f0
Xs = [i[1] for i in samples]

# ╔═╡ 814cf4ae-7487-11eb-1b18-c11f124dd18a
Ys = [i[2] for i in samples]

# ╔═╡ ff22c6e2-7487-11eb-3bbe-db691feb921b
Plots.histogram2d(Xs, Ys, nbins = 100)

# ╔═╡ 22183a3e-7487-11eb-192c-4370f0b71496
md"Now let's create an example that breaks MH"

# ╔═╡ 30131820-7487-11eb-11de-9561b77139e8
dist_hard = MixtureModel([MvNormal([-2.0, -2.0], [0.5, 0.5]),
		                  MvNormal([2.0, 2.0], [1.5, 0.1])],
						[0.5, 0.5])

# ╔═╡ b9c69864-7487-11eb-3942-e146b6ae5ea4
Plots.surface(-4:0.1:4, -4:0.1:4, (x, y) -> pdf(dist_hard, [x, y]), c = :blues)

# ╔═╡ 1c2dedee-7489-11eb-3a2c-519ab256505f
state_init_hard = [-1.2, -1.3]

# ╔═╡ b0030bf6-74a0-11eb-14c8-4f09bcb0327c
target_density_hard(state) = pdf(dist_hard, state);

# ╔═╡ d51078da-7488-11eb-334f-c1a0ce33384c
target_logdensity_hard(state) = logpdf(dist_hard, state);

# ╔═╡ ff1e6132-7488-11eb-3a82-d1e9fa247b42
samples_hard = OmegaMH.mh(rng,
						  target_logdensity_hard,
						  num_samples,
						  state_init,
						  propose_and_logratio)

# ╔═╡ 519dc07e-7489-11eb-295b-3769fa285c29
Xs_hard = [i[1] for i in samples_hard]

# ╔═╡ 5a8e043c-7489-11eb-23a5-9ba3521900bd
Ys_hard = [i[2] for i in samples_hard]

# ╔═╡ 4d3c151c-7489-11eb-2352-65f743d197d0
Plots.histogram2d(Xs_hard, Ys_hard, nbins = 100)

# ╔═╡ dd7f08f4-748a-11eb-15ae-1f79ca5c1a0c
md"Unfortunately, this doesn;t work in this new space; the MCMC chain fails to jump from mode to mode"

# ╔═╡ 92926c20-749c-11eb-37f7-e7ec05a5ee41
# relax(density, temp) = state -> 1 - exp(- density(state)/temp)
relax(density, temp) = state -> density(state)^(1/temp)

# ╔═╡ 98c5076a-74a1-11eb-0f34-99ffa441dee6
logrelax(logdensity, temp) = state -> (1/temp) * logdensity(state)

# ╔═╡ d1910068-74a0-11eb-21e6-e1d2410c1646
target_density_hard_rlx = relax(target_density_hard, 2.0)

# ╔═╡ a2b3e42c-74a0-11eb-1678-7148d7776b6e
Plots.surface(-4:0.1:4, -4:0.1:4, (x, y) -> target_density_hard_rlx([x,y]), c = :blues)

# ╔═╡ cd57df24-74a2-11eb-166f-31b691e72e12
temps = [1.0, 2.0, 20.0, 200.0]

# ╔═╡ e5df8ac4-74a2-11eb-3a36-31d7ac95fba8
function simulate_n(temp, state, samples_per_swap)
	relaxed_logdensity = logrelax(target_logdensity_hard, temp)
	OmegaMH.mh(rng,
			   relaxed_logdensity,
			   samples_per_swap,
			   state,
			   propose_and_logratio)
end

# ╔═╡ b75ae81e-74a3-11eb-2ac0-efd1a5c21f21
simulate_1(temp, state, samples_per_swap) =
	simulate_n(temp, state, samples_per_swap)[end]

# ╔═╡ e208a218-74a3-11eb-28ac-0dec5c478e15
function evaluate(temp, state)
	relaxed_logdensity = logrelax(target_logdensity_hard, temp)
	relaxed_logdensity(state)
end

# ╔═╡ ee2bb440-748a-11eb-046b-e5a4741ffe99
samples_hard_wow = re!(rng, temps, 1, 10000, [deepcopy(state_init) for i = 1:length(temps)], simulate_n, simulate_1, evaluate)

# ╔═╡ 59576824-74b3-11eb-1a98-d773404b49e0
unpack(s) = [i[1] for i in s], [i[2] for i in s]

# ╔═╡ b870ea20-74a4-11eb-2ab7-6972d5d2a105
Plots.histogram2d(unpack(samples_hard_wow)..., nbins = 100)

# ╔═╡ 25be5166-74a8-11eb-1255-bf738e624e97
Plots.contour(-4:0.1:4, -4:0.1:4, (x, y) -> target_density_hard_rlx([x,y]), c = :blues)

# ╔═╡ de1d4ab8-74a7-11eb-16a1-c728c0e288b8
function mcmc_dynamics(pdfs, traces, lbub, xrng = -4:0.1:4, yrng = -4:0.1:4)
	# plt = Plots.contoour(pdfs[1])
	plots = []
	for i = 1:length(pdfs)
		plt = Plots.contour(xrng, yrng,
							  (x, y) -> pdfs[i]([x, y]),
							   c = :blues,
							   aspect_ratio = 1.0)
		Xs = [j[1] for j in traces[i]]
	    Ys = [j[2] for j in traces[i]]
		plot!(plt, Xs[lbub], Ys[lbub])
			# , alpha = range(0.0; stop = 1.0, length = length(lbub)))
		push!(plots, plt)
	end
	Plots.plot(plots..., layour = (4, 1))
end

# ╔═╡ 745f702a-74a8-11eb-1ada-cf83599e0b17
pdfs = [relax(target_density_hard, temp) for temp in temps]

# ╔═╡ 0603ea0c-74ad-11eb-29ef-d58f5a17aba4
swap_every = 10

# ╔═╡ 0a36105a-74ad-11eb-27c8-ed6919e3c5c8
samples_per_swap = 10000

# ╔═╡ 17ace2ae-74ad-11eb-34be-3be5f9eb6f28
niters = swap_every * samples_per_swap

# ╔═╡ b82be6b2-74a8-11eb-1df2-9d26ac88ccbb
traces = re_all!(rng, temps, swap_every, samples_per_swap, [deepcopy(rand(dist_hard)) for i = 1:length(temps)], simulate_n, simulate_1, evaluate)

# ╔═╡ 4eb016a8-74b3-11eb-0fb3-e78e6117666b
Plots.histogram2d(unpack(traces[1])..., nbins = 100)

# ╔═╡ 9e771a08-74ac-11eb-1260-ad7e2af51b1b
@bind lb Slider(1:niters)

# ╔═╡ b8ac427c-74ac-11eb-0eca-a7234b0ab59a
@bind ub Slider(lb:niters)

# ╔═╡ 98226b20-74a8-11eb-28d4-c5d50971fcd5
mcmc_dynamics(pdfs, traces, lb:ub)

# ╔═╡ bd699f26-74ac-11eb-019e-bb6ef2c119d2

# ╔═╡ Cell order:
# ╠═9c5d8972-7485-11eb-25f7-155fbbd40d98
# ╟─128974e2-7483-11eb-03a0-15b2abc75d66
# ╠═2c24aaca-7483-11eb-05a9-6f847b7d09b7
# ╠═42bb4640-7483-11eb-0661-677c524c62b0
# ╠═eeb0be1c-7483-11eb-341f-b1996e9a53a6
# ╠═9d2f2714-749b-11eb-1403-e3a045df62ca
# ╠═ebd11c9a-7484-11eb-2941-c78e3e22d260
# ╠═0464d3ec-7489-11eb-1081-d761307263db
# ╠═94d6f4e6-74ac-11eb-0871-a31b2df5526d
# ╟─318e4a98-7483-11eb-1916-c9630ecf07e2
# ╠═27ced392-7483-11eb-3d36-7b8277daadff
# ╠═7e4e66ce-7483-11eb-238e-d15279086f3b
# ╠═9a7dc4f2-7483-11eb-36be-251be000793e
# ╠═68100d38-7483-11eb-07d5-45547ad45143
# ╠═f744509a-7483-11eb-3168-9bff07e7b245
# ╠═09fe4874-7484-11eb-28be-9f5df0586c87
# ╠═20537202-7484-11eb-1d70-a1bcbaa3c265
# ╠═3b0ac92e-7484-11eb-2e90-cf94eb2e9731
# ╠═4343d3b0-7484-11eb-388b-576f50c83335
# ╠═e84726c8-7484-11eb-1512-c5cee451507a
# ╠═f8e56030-7484-11eb-2bb6-afa2f952b085
# ╠═747691fe-7487-11eb-292e-2975ef75e7f0
# ╠═814cf4ae-7487-11eb-1b18-c11f124dd18a
# ╠═ff22c6e2-7487-11eb-3bbe-db691feb921b
# ╠═22183a3e-7487-11eb-192c-4370f0b71496
# ╠═30131820-7487-11eb-11de-9561b77139e8
# ╠═b9c69864-7487-11eb-3942-e146b6ae5ea4
# ╠═1c2dedee-7489-11eb-3a2c-519ab256505f
# ╠═b0030bf6-74a0-11eb-14c8-4f09bcb0327c
# ╠═d51078da-7488-11eb-334f-c1a0ce33384c
# ╠═ff1e6132-7488-11eb-3a82-d1e9fa247b42
# ╠═519dc07e-7489-11eb-295b-3769fa285c29
# ╠═5a8e043c-7489-11eb-23a5-9ba3521900bd
# ╠═4d3c151c-7489-11eb-2352-65f743d197d0
# ╠═dd7f08f4-748a-11eb-15ae-1f79ca5c1a0c
# ╠═92926c20-749c-11eb-37f7-e7ec05a5ee41
# ╠═98c5076a-74a1-11eb-0f34-99ffa441dee6
# ╠═d1910068-74a0-11eb-21e6-e1d2410c1646
# ╠═a2b3e42c-74a0-11eb-1678-7148d7776b6e
# ╠═cd57df24-74a2-11eb-166f-31b691e72e12
# ╠═e5df8ac4-74a2-11eb-3a36-31d7ac95fba8
# ╠═b75ae81e-74a3-11eb-2ac0-efd1a5c21f21
# ╠═e208a218-74a3-11eb-28ac-0dec5c478e15
# ╠═ee2bb440-748a-11eb-046b-e5a4741ffe99
# ╠═59576824-74b3-11eb-1a98-d773404b49e0
# ╠═b870ea20-74a4-11eb-2ab7-6972d5d2a105
# ╠═25be5166-74a8-11eb-1255-bf738e624e97
# ╠═de1d4ab8-74a7-11eb-16a1-c728c0e288b8
# ╠═745f702a-74a8-11eb-1ada-cf83599e0b17
# ╠═0603ea0c-74ad-11eb-29ef-d58f5a17aba4
# ╠═0a36105a-74ad-11eb-27c8-ed6919e3c5c8
# ╠═17ace2ae-74ad-11eb-34be-3be5f9eb6f28
# ╠═b82be6b2-74a8-11eb-1df2-9d26ac88ccbb
# ╠═4eb016a8-74b3-11eb-0fb3-e78e6117666b
# ╠═9e771a08-74ac-11eb-1260-ad7e2af51b1b
# ╠═b8ac427c-74ac-11eb-0eca-a7234b0ab59a
# ╠═98226b20-74a8-11eb-28d4-c5d50971fcd5
# ╠═bd699f26-74ac-11eb-019e-bb6ef2c119d2

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
#tcm_writeTreeToLatex

function tcm_writeTreeToLatex(outpath, outdir, tree_left, tree_right, tree_product, tree_isLeaf, var_names, doPDF)

	# temp_fn = x -> x ? "\\textbf{stop}" : "\\textbf{go}";
	# tree_action = map(temp_fn, tree_action)

	numNodes = length(tree_left);
	isclosed = zeros(Bool, numNodes)
	cn = 1;

	outhandle = open(outpath, "w")
	print(outhandle, "\\documentclass{standalone}\n")
	print(outhandle, "\\usepackage{forest}\n\n")

	print(outhandle, "\\begin{document}\n")
	print(outhandle, "\\begin{forest}\n")
	print(outhandle, "for tree={draw=none, l sep=40pt, s sep = 15pt, edge={->,thick}}\n")
	writeNodeToLatex(outhandle, 1, tree_left, tree_right, tree_product, tree_isLeaf, var_names)
	print(outhandle, "\\end{forest}\n")
	print(outhandle, "\\end{document}\n")
	close(outhandle)

	if (doPDF)
		run(`pdflatex -output-directory $outdir $outpath`)
	end

end

function writeNodeToLatex(outhandle, cn, tree_left, tree_right, tree_product, tree_isLeaf, var_names)

	if (tree_isLeaf[cn])
		if (!isempty(var_names))
			single_var_name = var_names[tree_product[cn]];
		else
			single_var_name = string(tree_product[cn]);
		end

		temp = string("[", single_var_name, ", draw")
		if (cn == 2)
			temp = string(temp, ", edge label={node[midway,left, yshift = 0.1cm]{\\emph{yes}}},")
		end

		if (cn == 3)
			temp = string(temp, ", edge label={node[midway,right, yshift = 0.1cm]{\\emph{no}}}")
		end
		temp = string(temp, "] \n");
		print(outhandle, temp);
	else

		if (!isempty(var_names))
			single_var_name = var_names[tree_product[cn]];
		else
			single_var_name = string(tree_product[cn]);
		end

		temp = string("[", single_var_name, ", ")
		if (cn == 2)
			temp = string(temp, "edge label={node[midway,left, yshift = 0.1cm]{\\emph{yes}}},")
		end

		if (cn == 3)
			temp = string(temp, "edge label={node[midway,right, yshift = 0.1cm]{\\emph{no}}}")
		end

		temp = string(temp, "\n");

		print(outhandle, temp)

		left = tree_left[cn];
		writeNodeToLatex(outhandle, left, tree_left, tree_right, tree_product, tree_isLeaf,  var_names)

		right = tree_right[cn];
		writeNodeToLatex(outhandle, right, tree_left, tree_right, tree_product, tree_isLeaf,  var_names)

		print(outhandle, "]\n")
	end

end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
function solve_stidsen(model ; showplot = false, docovercuts = true , global_branch = false, use_nsga = true, global_nsga = true, lift_covers = true, time_limit=300, fill_triangles=true, preprocess=true)

	tic()
	vm = copy(model)
	vd = getvOptData(vm)
	@assert length(vd.objs) == 2
	@assert all(isonlybinary.(vm, vd.objs))
	@assert vd.objSenses[1] == vd.objSenses[2]

	sense = vd.objSenses[1] == :Max ? Max : Min
	JuMP.setobjectivesense(vm, vd.objSenses[1])
	cstrData = ConstraintData(vm)

	z1, z2 = vd.objs
	status = global_branch ? solve(vm, method=:lexico, verbose=false) : solve(vm, method=:dicho, round_results=true)
	status == :Optimal || return status

	YN_convex = map(x->round.(Int, x), getY_N(vm))
	XE_convex = [[round(getvalue(JuMP.Variable(vm, i), j)) for i = 1:vm.numCols] for j = 1:length(YN_convex)]

	if length(YN_convex) == 1 || YN_convex[1] == YN_convex[end]
		return @NT(YN = YN_convex, XE = XE_convex, nodes=0)
	end

	if !issorted(YN_convex, by = first)
		s = sortperm(YN_convex, by = first)
        YN_convex = YN_convex[s]
        XE_convex = XE_convex[s]
	end

	if use_nsga
		modelnsga = copy(vm)
		modelnsga.ext[:vOpt].objs = [z1, z2]

		if global_nsga
			ns = nsga_binary(100, 500, modelnsga, seed=XE_convex, pmut=0.3, showprogress=false)
		else
			ns = union((nsga_binary(50, 500, modelnsga, seed=[XE_convex[i], XE_convex[i+1]], pmut=0.3, showprogress=false) for i = 1:length(XE_convex)-1)...)
			#if we do one nsga per triangle, we can have dominated solutions here.
			NSGAII.fast_non_dominated_sort!(ns, sense==Max ? NSGAII.Max() : NSGAII.Min())
			filter!(x->x.rank == 1, ns)
		end

		ns = unique(x->x.y, ns)
		sort!(ns, by=x->x.y[1])
		LN_NSGA = NonDomPoints(sense, map(x->x.pheno, ns), map(x->Tuple(x.y), ns))

		if fill_triangles
			diffs = [YN_convex[i] .- YN_convex[i+1] for i = 1:length(YN_convex)-1]
			areas = abs.(map(prod, diffs))
			while maximum(areas) > 2*mean(areas)
				i = indmax(areas)
				list_to_add = filter!(x-> YN_convex[i][1] < x.y[1] < YN_convex[i+1][1], ns)
				isempty(list_to_add) && break
				sol_to_add = list_to_add[length(list_to_add)÷2 + 1]
				insert!(YN_convex, i+1, sol_to_add.y)
				insert!(XE_convex, i+1, sol_to_add.pheno)
				diffs = [YN_convex[i] .- YN_convex[i+1] for i = 1:length(YN_convex)-1]
				areas = abs.(map(prod, diffs))
			end
		end
	end

	nbNodesTotal = 0
	LNGlobal = use_nsga ? LN_NSGA :  NonDomPoints(sense, XE_convex, Tuple.(YN_convex))

	t = toq()

	for i = 1:length(YN_convex)-1

		t > time_limit && break

		tic()
		#println("exploring triangle ◬ ($(YN_convex[i]) - $(YN_convex[i+1])) area : $(areas[i])")

		m = copy(vm)
		LN = NonDomPoints(sense, [XE_convex[i], XE_convex[i+1]], [Tuple(YN_convex[i]), Tuple(YN_convex[i+1])])

		if use_nsga
			for j = 1:length(LN_NSGA.yn)
				if LN.yn[1][1] < LN_NSGA.yn[j][1] < LN.yn[end][1]
					push!(LN, LN_NSGA.xe[j], LN_NSGA.yn[j])
				end
			end
		end

		Ƶ1, Ƶ2 = copy(z1, m), copy(z2, m)
		Ƶ = LN.λ[1]*Ƶ1 + LN.λ[2]*Ƶ2
		JuMP.setobjective(m, vd.objSenses[1], Ƶ.aff)

		if sense == Max
			bound1 = LN.yn[1][1] + 1
			bound2 = LN.yn[end][2] + 1
			cstr1 = @constraint(m, Ƶ1.aff >= bound1)
			cstr2 = @constraint(m, Ƶ2.aff >= bound2)
		else
			bound1 = LN.yn[end][1] - 1
			bound2 = LN.yn[1][2] - 1
			cstr1 = @constraint(m, Ƶ1.aff <= bound1)
			cstr2 = @constraint(m, Ƶ2.aff <= bound2)
		end

		#Stack of nodes to evaluate
		S = [Node(m, bound1, bound2, cstr1, cstr2)]

		if preprocess
			preprocess!(S[1], LN, Ƶ1, Ƶ2)
		end

		if docovercuts
			apply_cuts!(S[1], cstrData, lift_covers)
		end

		t += toq()
		#Solve while there are nodes to process
		cpt = 0
		while !isempty(S) && t < time_limit
			tic()
			sort!(S, by = x->x.zparent, rev=(sense==Min))
			process_node_stidsen(pop!(S), S, sense, LN, Ƶ1, Ƶ2, LNGlobal, cstrData, showplot)
			nbNodesTotal += 1
			cpt += 1
			t+=toq()
		end

		for i = 1:length(LN.xe)
			push!(LNGlobal, LN.xe[i], LN.yn[i])
		end

	end

	return @NT(YN = LNGlobal.yn	, XE = LNGlobal.xe, nodes = nbNodesTotal, timeout = t > time_limit)
end

function setRHS!(n::Union{Node, NodeParragh})
	JuMP.setRHS(n.cstr1, n.bound1)
	JuMP.setRHS(n.cstr2, n.bound2)
end

function fix_variables!(n::Union{Node, NodeParragh})
	for i in n.f1
		setlowerbound(JuMP.Variable(n.m, i), 1.)
	end
	for i in n.f0
		setupperbound(JuMP.Variable(n.m, i), 0.)
	end
end

function unfix_variables!(n::Union{Node, NodeParragh})
	for i in n.f1
		setlowerbound(JuMP.Variable(n.m, i), 0.)
	end
	for i in n.f0
		setupperbound(JuMP.Variable(n.m, i), 1.)
	end
end

function apply_cuts!(n::Node, cstrData, lift_covers, nb_try = 2)
	for i = 1:nb_try
		res = @suppress solve(n.m, ignore_solve_hook=true, relaxation=true)
		res != :Optimal && return
		n.z = round(getobjectivevalue(n.m), 8)
		n.x = round.(n.m.colVal, 8)
		success = find_cover_cuts(n, cstrData, lift_covers)
		!success && return
	end
	#println(n.m)
	return
end

function preprocess!(n, LN, obj1, obj2)
	cpt = 0
	for i = 1:n.m.numCols
		setlowerbound(JuMP.Variable(n.m, i), 1.)
		res = solve(n.m, ignore_solve_hook=true, relaxation=true, suppress_warnings=true)
		x = round.(n.m.colVal, 8)
		z = round(getobjectivevalue(n.m), 8)
		z1, z2 = evaluate(x, obj1), evaluate(x, obj2)
		setlowerbound(JuMP.Variable(n.m, i), 0.)
		setupperbound(JuMP.Variable(n.m, i), 0.)
		if isfathomable(z, LN) || res != :Optimal
			#println("fixed variable $i to 0 $res, $z    |||   $x")
			push!(n.f0, i)
			cpt += 1
		else
			res = solve(n.m, ignore_solve_hook=true, relaxation=true, suppress_warnings=true)
			x = round.(n.m.colVal, 8)
			z = round(getobjectivevalue(n.m), 8)
			z1, z2 = evaluate(x, obj1), evaluate(x, obj2)
			setupperbound(JuMP.Variable(n.m, i), 1.)
			if isfathomable(z, LN) || res != :Optimal
				#println("fixed variable $i to 1")
				setlowerbound(JuMP.Variable(n.m, i), 1.)
				push!(n.f1, i)
				cpt += 1
			else
			end
		end
	end
end

function process_node_stidsen(n::Node, S, sense, LN, obj1, obj2, LNGlobal, cstrData, showplot)

	n.zparent != Inf && isfathomable(n.zparent, LN) && return

	fix_variables!(n)
	setRHS!(n)
	res = solve(n.m, ignore_solve_hook=true, relaxation=true, suppress_warnings=true)
	unfix_variables!(n)

	res != :Optimal && return
	n.z = round(getobjectivevalue(n.m), 8)
	n.x = round.(n.m.colVal, 8)

	z1, z2 = evaluate(n.x, obj1), evaluate(n.x, obj2)

	if isfathomable(n.z, LN)
		showplot && plot_int_found(LN, LNGlobal, z1, z2, marker="k.", sleeptime=0.01)
		return
	end

	if isbinary(n.x)
		if isweaklydominated(z1, z2, LN)
			# println("int, dominated but not fathomable : paretobranch")
			showplot && plot_int_found(LN, LNGlobal, z1, z2, marker = "kx")
			append!(S, paretobranch(sense, n, z1, z2, LN, obj1, obj2, showplot))
		else
			# println("new int solution found")
			push!(LN, n.x, (z1, z2))
			append!(S, paretobranch(sense, n, z1, z2, LN, obj1, obj2, showplot))
			showplot && plot_int_found(LN, LNGlobal, z1, z2, sleeptime=1)
		end
	else
		# println("not binary")
		if isweaklydominated(z1, z2, LN)
			# println("not binary, dominated : paretobranch")
			showplot && plot_int_found(LN, LNGlobal, z1, z2, marker="r.")
			for node in paretobranch(sense, n, z1, z2, LN, obj1, obj2, showplot)
				push!(S, node)
			end
		else
			# println("not binary, not dominated : basic branch")
			showplot && plot_int_found(LN, LNGlobal, z1, z2, sleeptime=0.01, marker="k.")
			n1, n2 = basicbranch(n)
			push!(S, n1)
			push!(S, n2)
		end
	end
end

function paretobranch(sense::Type{Min}, n, z1, z2, LN, obj1, obj2, showplot)
	λ = LN.λ
	WS = sum((z1,z2).*λ)
	i = findfirst(x -> weakly_dominates(sense, x, (z1,z2)), LN.yn)
	#Find the first dominated nadir going left from (z1,z2)
	left = i-1
	while left > 0 && sum(nadirs(LN)[left].*λ) <= WS
		left -= 1
	end

	#Find the first dominated nadir going right from (z1,z2)
	right = i
	while right < length(LN) && sum(nadirs(LN)[right].*λ) <= WS
		right += 1
	end

	boundz1 = left > 0 ? nadirs(LN)[left][1] - 1. : NaN
	boundz2 = right < length(LN) ? nadirs(LN)[right][2] - 1. : NaN

	res = Node[]

	if !isnan(boundz1)
		n1 = shallowcopy(n)
		n1.bound1 = boundz1
		# @constraint(n1.m, copy(obj1, n1.m).aff <= boundz1)
		push!(res, n1)
	end
	if !isnan(boundz2)
		n2 = shallowcopy(n)
		n2.bound2 = boundz2
		# @constraint(n2.m, copy(obj2, n2.m).aff <= boundz2)
		push!(res, n2)
	end
	showplot && plot_pareto_branch(LN, z1, z2, boundz1, boundz2, sleeptime=0.01)
	#@assert length(res) >= 1
	res
end

function paretobranch(sense::Type{Max}, n, z1, z2, LN, obj1, obj2, showplot)
	λ = LN.λ
	WS = sum((z1,z2).*λ)
	i = findfirst(x -> weakly_dominates(sense, x, (z1,z2)), LN.yn)
	#Find the first dominated nadir going left from (z1,z2)
	left = i-1
	while left > 0 && sum(nadirs(LN)[left].*λ) >= WS
		left -= 1
	end

	#Find the first dominated nadir going right from (z1,z2)
	right = i
	while right < length(LN) && sum(nadirs(LN)[right].*λ) >= WS
		right += 1
	end

	boundz1 = right < length(LN) ? nadirs(LN)[right][1] + 0.5 : NaN
	boundz2 = left > 0 ? nadirs(LN)[left][2] + 0.5 : NaN

	res = Node[]

	if !isnan(boundz1)
		n1 = shallowcopy(n)
		n1.bound1 = boundz1
		# @constraint(n1.m, copy(obj1, n1.m).aff >= boundz1)
		push!(res, n1)
	end
	if !isnan(boundz2)
		n2 = shallowcopy(n)
		n2.bound2 = boundz2
		# @constraint(n2.m, copy(obj2, n2.m).aff >= boundz2)
		push!(res, n2)
	end
	showplot && plot_pareto_branch(LN, z1, z2, boundz1, boundz2, sleeptime=0.5)
	res
end

function basicbranch(n)
	i = indmax(map(x->x>0.5 ? 1-x : x, n.x))
	n1, n2 = shallowcopy(n), shallowcopy(n)
	push!(n1.f0, i)
	push!(n2.f1, i)
	n1, n2
end


text:
using OpenQuantumBase, Test

@test q_translate("ZZ+0.5ZI-XZ") == σz⊗σz + 0.5σz⊗σi - σx⊗σz
@test single_clause(["x"], [2], 0.5, 2) == 0.5σi⊗σx
@test single_clause(["x"], [2], 0.5, 2, sp = true) == 0.5spσi⊗spσx
@test single_clause(["z","z"], [2,3], -2, 4) == -2σi⊗σz⊗σz⊗σi
@test single_clause([σ₊], [3], 0.1, 3) == 0.1σi⊗σi⊗σ₊
@test single_clause([σ₊], [3], 0.1, 3, sp=true) == 0.1spσi⊗spσi⊗σ₊
@test standard_driver(2) == σx⊗σi + σi⊗σx
@test standard_driver(2, sp = true) == spσx⊗spσi + spσi⊗spσx
@test collective_operator("z", 3) ≈ σz⊗σi⊗σi + σi⊗σz⊗σi + σi⊗σi⊗σz
@test collective_operator("z", 3, sp = true) ≈ spσz⊗spσi⊗spσi + spσi⊗spσz⊗spσi + spσi⊗spσi⊗spσz
@test local_field_term([-1.0, 0.5], [1,3], 3) ≈ -1.0σz⊗σi⊗σi + 0.5σi⊗σi⊗σz
@test local_field_term([-1.0, 0.5], [1,3], 3, sp = true) ≈ -1.0spσz⊗spσi⊗spσi + 0.5spσi⊗spσi⊗spσz
@test two_local_term([-1.0, 0.5], [[1,3],[1,2]], 3) ≈ -1.0σz⊗σi⊗σz + 0.5σz⊗σz⊗σi

v0 = [1.0+0.0im, 0]
v1 = [0, 1.0+0.0im]
@test q_translate_state("0011") == v0⊗v0⊗v1⊗v1
@test q_translate_state("0.1(111)+(011)") == 0.1*v1⊗v1⊗v1 + v0⊗v1⊗v1
@test q_translate_state("(111)+(011)"; normal=true) == (v1⊗v1⊗v1 + v0⊗v1⊗v1)/sqrt(2)


text:
# This file is a part of SimilaritySearch.jl

function create_vectors(n, dim, normalize=false)
    D = [rand(Float32, dim) for i in 1:n]
    if normalize
        for u in D
            normalize!(u)
        end
    end

    D
end

function create_sequence(dim, sort)
    s = rand(1:10, dim)
    if sort
        sort!(s)
        s = unique(s)
    end

    s
end

include("testresults.jl")
include("testseq.jl")
include("testpivots.jl")
include("testsearchgraph.jl")


text:
"""
    Column

    Data structure for holding columns and related information.
"""
struct Column
    cost::Float64  # Native cost

    spind::Int  # Index of convexity constraint

    # Column coefficients, in sparse form.
    # These do not include the convexity constraint
    rind::Vector{Int}
    rval::Vector{Float64}

    is_vertex::Bool # Whether column is vertex (true) or ray (false)
end


text:
###############################################################################
#
# Make g_min and g_max control
#
###############################################################################
"""
    gsw_control!(
                photo_set::AbstractPhotoModelParaSet{FT},
                canopyi::CanopyLayer{FT},
                envir::AirLayer{FT},
                ind::Int
    ) where {FT<:AbstractFloat}
    gsw_control!(
                photo_set::AbstractPhotoModelParaSet{FT},
                canopyi::CanopyLayer{FT},
                envir::AirLayer{FT}
    ) where {FT<:AbstractFloat}

make sure g_sw is in its physiological range limited by diffusion, given
- `photo_set` [`C3ParaSet`] or [`C4ParaSet`] type parameter set
- `canopyi` [`CanopyLayer`](@ref) type struct
- `envir` [`AirLayer`] type struct
- `ind` Nth leaf

Note that this function is meant to use jointly with gas_exchange! when
    computing optimal stomtal conductance.
"""
function gsw_control!(
            photo_set::AbstractPhotoModelParaSet{FT},
            canopyi::CanopyLayer{FT},
            envir::AirLayer{FT},
            ind::Int
) where {FT<:AbstractFloat}
    # if g_sw is low than g_min
    if canopyi.g_sw[ind] < canopyi.g_min
        gas_exchange!(photo_set, canopyi, envir, GswDrive(), ind,
                      canopyi.g_min);
    elseif canopyi.g_sw[ind] > canopyi.g_max
        gas_exchange!(photo_set, canopyi, envir, GswDrive(), ind,
                      canopyi.g_max);
    end

    return nothing
end

function gsw_control!(
            photo_set::AbstractPhotoModelParaSet{FT},
            canopyi::CanopyLayer{FT},
            envir::AirLayer{FT}
) where {FT<:AbstractFloat}
    # if g_sw is low than g_min
    for i in eachindex(canopyi.g_sw)
        canopyi.ps.APAR = canopyi.APAR[i];
        leaf_ETR!(photo_set, canopyi.ps);
        gsw_control!(photo_set, canopyi, envir, i);
    end

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
solve(i::MSetInstance, ::DynamicProgramming; kwargs...) = msets_dp(i; kwargs...)

function msets_dp(i::MSetInstance)
  # V[µ, δ]: µ items, up to δ.
  V = Matrix{Float64}(undef, m(i), dimension(i))
  S = Dict{Tuple{Int, Int}, Vector{Int}}()

  # Initialise: µ == 1, just take the best element among [δ+1, d]; δ == d, take nothing.
  V[1, dimension(i)] = 0.0
  S[1, dimension(i)] = Int[]

  for δ in (dimension(i) - 1):-1:1
     v, x = findmax(values(i, (δ + 1):dimension(i)))
     V[1, δ] = v
     S[1, δ] = [x]
  end

  for µ in 2:m(i)
    V[µ, dimension(i)] = 0.0
    S[µ, dimension(i)] = Int[]
  end

  # Dynamic part.
  for µ in 2:m(i)
    for δ in (dimension(i) - 1):-1:0
      obj_idx = δ + 1
      take_δ = value(i, obj_idx) + V[µ - 1, δ + 1]
      dont_take_δ = V[µ, δ + 1]

      if take_δ > dont_take_δ
        V[µ, δ + 1] = take_δ
        S[µ, δ + 1] = vcat(δ, S[µ - 1, δ + 1])
      else
        V[µ, δ + 1] = dont_take_δ
        S[µ, δ + 1] = S[µ, δ + 1]
      end
    end
  end

  return MSetSolution(i, S[m(i), 1] .+ 1)
end


text:
module StateSpaceStudy

    using LinearAlgebra, Random, Distributions, PyPlot

    include("filteringrecursion.jl")

    include("smoothingrecursion.jl")

    include("simulatedata.jl")

    include("Example_4_3.jl")

    include("Example_4_2.jl")

    export simulatedata, filteringrecursion, smoothingrecursion

end


text:
import Base.Dates: Day, Hour

warsaw = resolve("Europe/Warsaw", tzdata["europe"]...)

normal = DateTime(2015, 1, 1, 0)   # 24 hour day in warsaw
spring = DateTime(2015, 3, 29, 0)  # 23 hour day in warsaw
fall = DateTime(2015, 10, 25, 0)   # 25 hour day in warsaw

# Unary plus
@test +ZonedDateTime(normal, warsaw) == ZonedDateTime(normal, warsaw)

# Period arithmetic
@test ZonedDateTime(normal, warsaw) + Day(1) == ZonedDateTime(normal + Day(1), warsaw)
@test ZonedDateTime(spring, warsaw) + Day(1) == ZonedDateTime(spring + Day(1), warsaw)
@test ZonedDateTime(fall, warsaw) + Day(1) == ZonedDateTime(fall + Day(1), warsaw)

@test ZonedDateTime(normal, warsaw) + Hour(24) == ZonedDateTime(normal + Hour(24), warsaw)
@test ZonedDateTime(spring, warsaw) + Hour(24) == ZonedDateTime(spring + Hour(25), warsaw)
@test ZonedDateTime(fall, warsaw) + Hour(24) == ZonedDateTime(fall + Hour(23), warsaw)

# Do the same calculations but backwards over the transitions.
@test ZonedDateTime(normal + Day(1), warsaw) - Day(1) == ZonedDateTime(normal, warsaw)
@test ZonedDateTime(spring + Day(1), warsaw) - Day(1) == ZonedDateTime(spring, warsaw)
@test ZonedDateTime(fall + Day(1), warsaw) - Day(1) == ZonedDateTime(fall, warsaw)

@test ZonedDateTime(normal + Day(1), warsaw) - Hour(24) == ZonedDateTime(normal, warsaw)
@test ZonedDateTime(spring + Day(1), warsaw) - Hour(23) == ZonedDateTime(spring, warsaw)
@test ZonedDateTime(fall + Day(1), warsaw) - Hour(25) == ZonedDateTime(fall, warsaw)

# Ensure that arithmetic around transitions works.
@test ZonedDateTime(spring, warsaw) + Hour(1) == ZonedDateTime(spring + Hour(1), warsaw)
@test ZonedDateTime(spring, warsaw) + Hour(2) == ZonedDateTime(spring + Hour(3), warsaw)
@test ZonedDateTime(fall, warsaw) + Hour(2) == ZonedDateTime(fall + Hour(2), warsaw, 1)
@test ZonedDateTime(fall, warsaw) + Hour(3) == ZonedDateTime(fall + Hour(2), warsaw, 2)

# Non-Associativity
hour_day = (ZonedDateTime(spring, warsaw) + Hour(24)) + Day(1)
day_hour = (ZonedDateTime(spring, warsaw) + Day(1)) + Hour(24)

@test hour_day - day_hour == Hour(1)

# Arithmetic with a StepRange should always work even when the start/stop lands on
# ambiguous or non-existent DateTimes.
ambiguous = DateTime(2015, 10, 25, 2)  # Ambiguous hour in Warsaw
nonexistent = DateTime(2014, 3, 30, 2)  # Non-existent hour in Warsaw

range = ZonedDateTime(ambiguous - Day(1), warsaw):Hour(1):ZonedDateTime(ambiguous - Day(1) + Hour(1), warsaw)
@test range .+ Day(1) == ZonedDateTime(ambiguous, warsaw, 1):Hour(1):ZonedDateTime(ambiguous + Hour(1), warsaw)

range = ZonedDateTime(ambiguous - Day(1) - Hour(1), warsaw):Hour(1):ZonedDateTime(ambiguous - Day(1), warsaw)
@test range .+ Day(1) == ZonedDateTime(ambiguous - Hour(1), warsaw, 1):Hour(1):ZonedDateTime(ambiguous, warsaw, 2)

range = ZonedDateTime(nonexistent - Day(1), warsaw):Hour(1):ZonedDateTime(nonexistent - Day(1) + Hour(1), warsaw)
@test range .+ Day(1) == ZonedDateTime(nonexistent + Hour(1), warsaw):Hour(1):ZonedDateTime(nonexistent + Hour(1), warsaw)

range = ZonedDateTime(nonexistent - Day(1) - Hour(1), warsaw):Hour(1):ZonedDateTime(nonexistent - Day(1), warsaw)
@test range .+ Day(1) == ZonedDateTime(nonexistent - Hour(1), warsaw):Hour(1):ZonedDateTime(nonexistent - Hour(1), warsaw)


text:
# change dir
cd("G:\\My Drive\\18. Github Repo\\DataScience");

# versioninfo
versioninfo()

# load some analysis and data handling libs
using BenchmarkTools, DataFrames, DelimitedFiles, CSV, XLSX

# get the data
# download it as .csv file
P = download("https://raw.githubusercontent.com/nassarhuda/easy_data/master/programming_languages.csv", "programminglanguages.csv");

# read in the csv file
# P is Matrix{Any} very hard to work with
P, H = readdlm("programminglanguages.csv", ','; header=true);

# using CSV into DataFrame
C = CSV.read("programminglanguages.csv", DataFrame);

# look at and access data frame
C[1:10, :]
C.year
names(C)

# short summary
describe(C)

# test the performances
@btime P, H = readdlm("programminglanguages.csv", '.', header=true);
@btime C = CSV.read("programminglanguages.csv", DataFrame);

# wwrite csv
CSV.write("test.csv", C);

# Switch to XLSX
T = XLSX.readdata(".\\data\\zillow_data_download_april2020.xlsx", "Sale_counts_city", "A1:F9");
typeof(T)

# Matrix Any is hard to work with
df = DataFrame(T[2:end, :]);
rename!(df, names(df) .=> Symbol.(T[1,:]))
# equiv to
df = DataFrame(T[2:end, :], Symbol.(T[1, :]))

# convert everything the reasonable data type
nm = names(df);
for name in nm
    type = typeof(df[2, Symbol(name)]);
    try
        df[:, Symbol(name)] = convert(Array{type, 1}, df[:, Symbol(name)]);
    catch er;
    end
end
df

# read the whole table
@time G = XLSX.readtable("data/zillow_data_download_april2020.xlsx","Sale_counts_city");
dg = DataFrame(G...)

# convert the dg into proper type
nm = names(dg);
for name in nm
    type = typeof(dg[2, Symbol(name)]);
    try
        dg[!, Symbol(name)] = convert(Array{Union{Missing, type}, 1}, dg[!, Symbol(name)]);
    catch er;
        println("Column $name missing value detected ... ");
        throw(er);
    end
end
dg

# using data frame:C to find the language
function find_lang(data::DataFrame, lang::String)
    idx = findfirst(data[!, :language] .== lang);
    if !isnothing(idx)
        return data[idx, :year];
    else
        @warn("-- $lang not found... ");
    end
end

find_lang(C, "Julia")

# using the dataframe to find number of languages
function number_lang(data::DataFrame, year::Int64)
    ans = length(findall(data[!, :year] .== year));
    return ans;
end

number_lang(C, 2011)

max_stars_repo_path: 01. Data.jl

max_issues_repo_path: 01. Data.jl

max_forks_repo_path: 01. Data.jl


text:
module Examples

include("example_base.jl")

include("arnold_cat_map.jl")
include("bakers_map.jl")
include("henon_map.jl")
include("linz_sprott_99.jl")
include("lorenz_63.jl")
include("lozi_map.jl")
include("rnn.jl")
include("standard_map.jl")
include("van_der_pol.jl")

end

using .Examples.ExampleBase: LEDemo, LEExample,
    ContinuousExample, DiscreteExample

# Make LyapunovExponents.lorenz_63 etc. available:
using .Examples.ArnoldCatMap
using .Examples.BakersMap
using .Examples.HenonMap
using .Examples.LinzSprott99
using .Examples.Lorenz63
using .Examples.LoziMap
using .Examples.RNN
using .Examples.StandardMap
using .Examples.VanDerPol

const DEMOS = [
    lorenz_63,
    linz_sprott_99,
    van_der_pol,
    beer_95,
    henon_map,
    standard_map,
    bakers_map,
    arnold_cat_map,
    lozi_map,
]


text:
@testset "Algebra" begin
    Mod.@polyvar x y
    @test 2 .- ((1 .+ (-x)) .* 4) ./ 2 == x.^2 .* (1 ./ x) .* 2
    @test dot(0, x^2 - 2*x^2) == dot((x^2 - x)', x^2 - x^2)
    @test dot(x + 1, 2) == 2x + 2
    @test 2 .* x .+ 2 == (x + 3) .+ (x .- 1)
    @test ((x + y) .- y) ./ y == x / y
    @test (-2)*x == -(2*x)
    @test x * x == x^2
    @test 4x == 2.0(2x)
    @test x * (2x) == 2x^2
    @inferred 2 * (2x)
    @inferred 2.0 * (2x)
    #@test eltype(2.0 * (2x)) == Float64
    t = @inferred 1.5 * (2x)
    @test t isa AbstractTerm{Float64}
    @test t == 3x

    p = @inferred y^2 * (2.0x^2 + 3y)
    @test p == 2x^2*y^2 + 3y^3
    @test p isa AbstractPolynomial{Float64}

    p = @inferred (2x^2 + 3.0y) * y^2
    @test p == 2x^2*y^2 + 3y^3
    @test p isa AbstractPolynomial{Float64}

    p = @inferred (2.0y^2) * (2x^2 + 3y)
    @test p == 4x^2*y^2 + 6y^3
    @test p isa AbstractPolynomial{Float64}

    p = @inferred (2x^2 + 3y) * (2.0y^2)
    @test p == 4x^2*y^2 + 6y^3
    @test p isa AbstractPolynomial{Float64}

    p = @inferred (y^2 + 2.0) * (2x + y)
    @test p == 2x*y^2 + y^3 + 4x + 2y
    @test p isa AbstractPolynomial{Float64}

    p = @inferred CustomPoly(x + 2y) * CustomTerms(x^2 + 1.0)
    @test p == x^3 + 2x^2*y + x + 2y
    @test p isa AbstractPolynomial{Float64}

    p = @inferred CustomPoly(x + 2y) + 1
    @test p == x + 2y + 1
    @test p isa AbstractPolynomial{Int}

    p = @inferred CustomPoly(x + 2y) + 1.0
    @test p == x + 2y + 1
    @test p isa AbstractPolynomial{Float64}

    @testset "Inference" begin
        @inferred x^2
        @inferred x^2-2x
        @inferred x+2x
        @inferred -2x+2
        @inferred x^2-2x+2
        @inferred x^2-2x
        @inferred 2x+2
        @inferred (1+x)^0
        @inferred (1+x)^1
        @inferred (1+x)^2
        @inferred (1+x)^3
        @inferred (x^2)^0
        @inferred (x^2)^1
        @inferred (x^2)^2
        @inferred (x^2)^3
        @inferred (2x) * (3x)
        @inferred (2x)^0
        @inferred (2x)^1
        @inferred (2x)^2
        @inferred (2x)^3
    end

    @test iszero((x+x-2*x) * (x * (x^2 + y^2)))
    @test iszero((0*x) * (x*y * (x^2 + y^2)))

    @testset "Scalar - Array" begin
        Mod.@polyvar x y
        @test x + [x^2+y y; y x*y] == [x+x^2+y x+y; x+y x+x*y]
        @test [x^2+y y; y x*y] + x == [x^2+x+y y+x; y+x x*y+x]

        @test x * [1+y y] == [x*(1+y) x*y]
        @test [1+y y] * x == [(1+y)*x y*x]
        @test [1+y y] / x == [(1+y)/x y/x]
    end
end


text:
using Catlab.CategoricalAlgebra.CSets
using Catlab.Theories
using AutoHashEquals

"""
Data for an individual component (each vector corresponds to its elements)
1.) how many of each color (for each in-arrow) targets this point
2.) what color this point targets (for each out arrow)

This could be extended to add extra automorphism-invariant properties.
"""
@auto_hash_equals struct CDataPoint
  indata::Vector{Vector{Int}}
  outdata::Vector{Int}
end

"""Data required to color a CSet (each element of each component)"""
const CData = Dict{Symbol, Vector{CDataPoint}}

"""
Computes colors for a CSet, distinguishing nodes by their immediate
connectivity. It is not sufficient to compute the automorphism group, but it is
a good starting point.

This does not generalize to ACSets. We cannot naively throw the attributes as
raw data into the color data. It will make indistinguishable elements (e.g. two
elements that map to different data but otherwise can be permuted) as
distinguishable.
"""
function compute_color_data(g::StructACSet{S}, color::CDict)::CData where {S}
  tabs, arrs, srcs, tgts = ob(S), hom(S), dom(S), codom(S)
  res = CData()
  for tab in tabs # compute colordata for each tab
    subres = []  # vector of data for each in-arrow
    for (arr, src, _) in filter(x->x[3]==tab, collect(zip(arrs,srcs,tgts))) # each incoming arrow
      color_src = color[src]
      subsubres = zeros(Int, nparts(g, tab), max0(color_src))
      for (colorsrc, arrtgt) in zip(color_src, g[arr])
        subsubres[arrtgt, colorsrc] += 1
      end
      push!(subres, subsubres)
    end

    # Also compute per-element data for table `tgt` (now, regard as a src)
    out_subres = Vector{Int}[color[tgt][g[oga]]
                             for (oga, src, tgt) in zip(arrs, srcs, tgts)
                             if src==tab]

    # Combine the two pieces of data for each elmeent in tgt, store in res
    res[tab] = [CDataPoint([ssr[i,:] for ssr in subres],
                           [osr[i] for osr in out_subres])
                for i in parts(g,tab)]
  end
  return res
end

nocolor(g::StructACSet{S}) where {S} = CDict([k => ones(Int, nparts(g, k)) for k in ob(S)])

"""
Iterative color refinement based on the number (and identity) of incoming and
outgoing arrows.
Inputs:
 - g: CSet we are color saturating
 - init_color: initial coloring, if any (default: uniform)
 - history: record each iteration (if true)
Returns:
- trajectory of colorings
- permutation-invariant hash value for final coloring (indicator fn)
"""
function color_saturate(g::StructACSet{S};
                      init_color::Union{Nothing,CDict}=nothing,
                      history::Bool=false
                     )::Pair{Vector{CDict},UInt64} where {S}
  res = CDict[]
  # Default: uniform coloring
  new_color = init_color === nothing ? nocolor(g) : init_color

  prev_n, curr_n, iter = 0, 1, 0
  hashes = Dict{Symbol, Vector{UInt}}()
  while prev_n != curr_n
    if history
      push!(res, new_color)
    end
    iter += 1
    prev_color = new_color
    # All that matters about newdata's type is that it is hashable
    newdata = compute_color_data(g, prev_color)
    # Distinguish by both color AND newly computed color data
    new_datahash = Dict{Symbol, Vector{UInt}}(
      [k=>map(hash, zip(prev_color[k],v)) for (k, v) in collect(newdata)])
    # Identify set of new colors for each component
    hashes = Dict{Symbol, Vector{UInt}}(
      [k=>sort(collect(Set(v))) for (k, v) in new_datahash])
    # Assign new colors by hash value of color+newdata
    new_color = CDict([
      k=>[findfirst(==(new_datahash[k][i]), hashes[k])
          for i in 1:nparts(g, k)]
      for (k, v) in new_datahash])
    prev_n = sum(map(max0, values(prev_color)))
    curr_n = sum(map(max0, values(new_color)))
  end
  push!(res, new_color)
  indicator = hash(Dict([k=>Set(v) for (k,v) in hashes]))
  return res => indicator
end


text:
module P4est_wrapper

################################################################
# Check correct installation
################################################################
# Load in `deps.jl`, complaining if it does not exist
const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("P4est_wrapper was not build properly. Please run Pkg.build(\"p4est.jl\").")
end
include(depsjl_path)
# Check if p4est library was found
if(!P4EST_FOUND)
    error("P4est_wrapper was not build properly. Please run Pkg.build(\"p4est.jl\").")
end

# Use library path from build process
const p4est_lib = P4EST_LIB

################################################################
# Use dependencies
################################################################
using Libdl
using MPI
using CEnum

################################################################
# Include p4est wrappers
################################################################
include("ctypes.jl")
export Ctm, Ctime_t, Cclock_t

# Include P4EST bindings

# The following file fix and/or complete
# the automated generation of types and data structures
include(joinpath(@__DIR__, "common_fixes.jl"))

# The following bindings are automatically generated.
# Read more info in gen/README.md
bindings_dir = joinpath(@__DIR__, "bindings")
include(joinpath(bindings_dir, "sc_common.jl"))
include(joinpath(bindings_dir, "sc_api.jl"))
include(joinpath(bindings_dir, "p4est_common.jl"))
include(joinpath(bindings_dir, "p4est_api.jl"))
include(joinpath(bindings_dir, "p8est_common.jl"))
include(joinpath(bindings_dir, "p8est_api.jl"))
include(joinpath(bindings_dir, "p6est_common.jl"))
include(joinpath(bindings_dir, "p6est_api.jl"))

# The following file fix and/or complete
# the automated generation of functions
include(joinpath(@__DIR__, "api_fixes.jl"))

################################################################
# Export everything starting with ...
################################################################
foreach(names(@__MODULE__, all=true)) do s
    if startswith(string(s), "sc") || startswith(string(s), "p4est") || startswith(string(s), "p6est") || startswith(string(s), "p8est")
        @eval export $s
    end
end

end # module


text:
using DashBootstrapComponents

items = [
    dbc_dropdownmenuitem("Action"),
    dbc_dropdownmenuitem("Another action"),
    dbc_dropdownmenuitem("Something else here"),
    dbc_dropdownmenuitem(divider = true),
    dbc_dropdownmenuitem("Something else here after the divider"),
];

dropdown = dbc_row([
    dbc_col(
        dbc_dropdownmenu(
            label = "Left-aligned menu (default)",
            children = items,
            align_end = false,
        ),
    ),
    dbc_col(
        dbc_dropdownmenu(label = "Right-aligned menu", children = items, align_end = true),
    ),
]);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
"""
    relerror(beam1::Particles, beam2::Particles)

Compute the relative error of particle field from `beam1` comparing to `beam2`.
"""
function relerror(beam1::Particles, beam2::Particles)
    beam1.npar == beam2.npar || error("npar mismatches")
    npar = beam1.npar
    efield_error = relerror(beam1.self_efields, beam2.self_efields, npar)
    bfield_error = relerror(beam1.self_bfields, beam2.self_bfields, npar)
    return max(efield_error, bfield_error)
end

"""
    relerror(A::AbstractVector{SVector{3,T}}, B::AbstractVector{SVector{3,T}}, n) where {T}

Compute the relative error of an array of Vector3 `A` comparing to another array of Vector3
`B`.
"""
function relerror(A::AbstractVector{SVector{3,T}}, B::AbstractVector{SVector{3,T}}, n) where {T}
    denominator = 0.0
    numerator = 0.0
    for i in 1:n
        a = A[i]
        b = B[i]
        numerator += norm(a - b)^2
        denominator += norm(b)^2
    end
    return sqrt(numerator / denominator)
end


text:
export split_dims

# Split the dimension of a ND-tensor into 3 parts:
#   (dim_pre, dim_mid, dim_post)
function split_dims{T}(tensor::T, dim::Int)
  dims = size(tensor)
  dim_pre  ::Int = prod(dims[1:dim-1])
  dim_mid  ::Int = dims[dim]
  dim_post ::Int = prod(dims[dim+1:end])

  (dim_pre, dim_mid, dim_post)
end


text:
################################################################################
#
#  Triangularization
#
################################################################################

function _pivot(A, start_row, col)
  if !iszero(A[start_row, col])
    return 1;
  end

  for j in start_row + 1:nrows(A)
    if !iszero(A[j, col])
      swap_rows!(A, j, start_row)
      return -1
    end
  end

  return 0
end

function _strong_echelon_form(A::Generic.Mat{NfOrdQuoRingElem}, strategy)
  B = deepcopy(A)

  if nrows(B) < ncols(B)
    B = vcat(B, zero_matrix(base_ring(B), ncols(B) - nrows(B), ncols(B)))
  end

  if strategy == :split
    q, w = z_split(ideal(base_ring(A)))
    R = order(ideal(base_ring(A)))
    ideals = q
    if length(w) != 0
      push!(ideals, prod(w))
    end
    C = _strong_echelon_form_split(B, ideals)
    return C
  elseif strategy == :no_split
    C = _strong_echelon_form_nonsplit(B)
    return C
  else
    error("Invalid strategy")
  end
end

function strong_echelon_form(A::Generic.Mat{NfOrdQuoRingElem}, shape::Symbol = :upperright, strategy::Symbol = :split)
  if shape == :lowerleft
    h = _strong_echelon_form(invert_cols(A), strategy)
    invert_cols!(h)
    invert_rows!(h)
    return h
  elseif shape == :upperright
    h = _strong_echelon_form(A, strategy)
    return h
  else
    error("Not yet implemented")
  end
end

function triangularize!(A::Generic.Mat{NfOrdQuoRingElem})
  n = nrows(A)
  m = ncols(A)
  d = one(base_ring(A))

  t_isdiv = 0.0
  t_xxgcd = 0.0
  t_arith = 0.0

  row = 1
  col = 1
  while row <= nrows(A) && col <= ncols(A)
    #println("doing row $row")
    t = _pivot(A, row, col)
    if t == 0
      col = col + 1
      continue
    end
    d = d*t
    for i in (row + 1):nrows(A)
      if iszero(A[i, col])
        continue
      end

      t_isdiv += @elapsed b, q = isdivisible(A[i, col], A[row, col])

      if b
        for k in col:m
          t_arith += @elapsed A[i, k] = A[i, k] - q*A[row, k]
        end
        @hassert :NfOrdQuoRing 1 A[i, col] == zero(base_ring(A))
      else
        t_xxgcd += @elapsed g,s,t,u,v = xxgcd(A[row, col], A[i, col])
        @hassert :NfOrdQuoRing 1 isone(s*v - t*u)

        for k in col:m
          t_arith += @elapsed t1 = s*A[row, k] + t*A[i, k]
          t_arith += @elapsed t2 = u*A[row, k] + v*A[i, k]
          A[row, k] = t1
          A[i, k] = t2
        end
      end
    end
    row = row + 1;
    col = col + 1;
  end
  #println("  === Time triangularization")
  #println("    isdivisbible: $t_isdiv")
  #println("    xxgcd       : $t_xxgcd")
  #println("    arith       : $t_arith")
  #println("    total time  : $(toc())")
  return d
end

function triangularize(A::Generic.Mat{NfOrdQuoRingElem})
  #println("copying ...")
  B = deepcopy(A)
  #println("done")
  triangularize!(B)
  return B
end

################################################################################
#
#  Strong echelon form
#
################################################################################

# Naive version of inplace strong echelon form
# It is assumed that A has more rows then columns.
function strong_echelon_form_naive!(A::Generic.Mat{NfOrdQuoRingElem})
  #A = deepcopy(B)
  n = nrows(A)
  m = ncols(A)

  @assert n >= m

  #print("triangularizing ... ")
  triangularize!(A)
  #println("done")

  T = zero_matrix(base_ring(A), 1, ncols(A))

  # We do not normalize!
  for j in 1:m
    if !iszero(A[j,j]) != 0
      # This is the reduction
      for i in 1:j-1
        if iszero(A[i, j])
          continue
        else
          q, r = divrem(A[i, j], A[j, j])
          for l in i:m
            A[i, l] = A[i, l] - q*A[j, l]
          end
        end
      end

      a = annihilator(A[j, j])

      for k in 1:m
        T[1, k] = a*A[j, k]
      end
    else
      for k in 1:m
        T[1, k] = A[j, k]
      end
    end

    for i in j+1:m

      if iszero(T[1, i])
        continue
      end

      if iszero(A[i, i])
        for k in i:m
          T[1, k], A[i, k] = A[i, k], T[1, k]
        end
      else
        b, q = isdivisible(T[1, i], A[i, i])
        if b
          for k in i:m
            T[1, k] = T[1, k] - q*A[i, k]
          end
          @hassert :NfOrdQuoRing 1 T[1, i] == zero(base_ring(A))
        else
          g,s,t,u,v = xxgcd(A[i, i], T[1, i])

          for k in i:m
            t1 = s*A[i, k] + t*T[1, k]
            t2 = u*A[i, k] + v*T[1, k]
            A[i, k] = t1
            T[1, k] = t2
          end
        end
      end
    end
  end
  return A
end

################################################################################
#
#  Howell form
#
################################################################################

function howell_form!(A::Generic.Mat{NfOrdQuoRingElem})
  @assert nrows(A) >= ncols(A)

  k = nrows(A)

  strong_echelon_form_naive!(A)

  for i in 1:nrows(A)
    if iszero_row(A, i)
      k = k - 1

      for j in (i + 1):nrows(A)
        if !iszero_row(A, j)
          swap_rows!(A, i, j)
          j = nrows(A)
          k = k + 1
        end
      end
    end
  end
  return k
end

function howell_form(A::Generic.Mat{NfOrdQuoRingElem})
  B = deepcopy(A)

  if nrows(B) < ncols(B)
    B = vcat(B, zero_matrix(base_ring(B), ncols(B) - nrows(B), ncols(B)))
  end

  howell_form!(B)

  return B
end

################################################################################
#
#  Determinant
#
################################################################################

function det(M::Generic.Mat{NfOrdQuoRingElem})
  nrows(M) != ncols(M) && error("Matrix must be square matrix")
  N = deepcopy(M)
  d = triangularize!(N)
  z = one(base_ring(M))
  for i in 1:nrows(N)
    z = z * N[i, i]
  end
  return z*d
  q, r = divrem(z, d)
  @hassert :NfOrdQuoRing 1 iszero(r)
  return divexact(z, d)
end

################################################################################
#
#  Z Split
#
################################################################################

function z_split(I::NfOrdIdl)
  b = basis_matrix(I)
  O = order(I)
  n = degree(O)
  c = coprime_base([b[i, i] for i in 1:n])
  nI = norm(I)
  if isone(nI)
    return NfOrdIdl[I], NfOrdIdl[]
  end
  val = Vector{Int}(undef, length(c))
  for i in 1:length(c)
    val[i] = valuation(nI, c[i])
  end
  if n == 1
    nz = one(FlintZZ)
  else
    nz = prod(b[i, i] for i in 2:n)
  end

  A = NfOrdIdl[]
  B = NfOrdIdl[]

  for i in 1:length(c)
    a = ideal(O, c[i]^val[i]) + I
    if iszero(mod(nz, c[i]))
      push!(B, a)
    else
      push!(A, a)
    end
  end
  return A, B
end

function can_map_into_integer_quotient(Q::NfOrdQuoRing)
  B = basis_matrix(ideal(Q))
  for i in 2:ncols(B)
    if !isone(B[i, i])
      return false
    end
  end
  return true
end

function map_into_integer_quotient(Q::NfOrdQuoRing)
  B = basis_matrix(ideal(Q))
  m = B[1, 1]
  R = ResidueRing(FlintZZ, m, cached = false)
  f = (x -> R(x.elem.coordinates[1]))
  g = (y -> Q(y.data)::NfOrdQuoRingElem)
  return R, f, g
end

function can_make_small(Q::Generic.ResRing{fmpz})
  if nbits(modulus(Q)) < Sys.WORD_SIZE - 1
    return true
  else
    return false
  end
end

function make_small(Q::Generic.ResRing{fmpz})
  R = ResidueRing(FlintZZ, Int(modulus(Q)), cached = false)
  f = (x -> R(x.data)::nmod)
  g = (x -> Q(x.data)::Generic.Res{fmpz})
  return R, f, g
end

function _strong_echelon_form_split(M::MatElem{NfOrdQuoRingElem}, ideals)
  Q = base_ring(M)
  R = base_ring(Q)
  modulus = ideal(Q)

  n = nrows(M)
  m = ncols(M)

  M_cur = zero_matrix(Q, n, m)

  if length(ideals) == 1
    return _strong_echelon_form_nonsplit(M)
  end

  I = ideals[1]

  RmodI, mRmodI = quo(R, I)
  MmodI = zero_matrix(RmodI, n, m)
  for i in 1:n
    for j in 1:m
      MmodI[i, j] = RmodI(lift(R, M[i, j]))
    end
  end
  echelon_modI = _strong_echelon_form_nonsplit(MmodI)

  for i in 1:n
    for j in 1:m
      M_cur[i, j] = Q(lift(R, echelon_modI[i, j]))
    end
  end

  _assure_weakly_normal_presentation(I)
  gI = gcd(Q(I.gen_one), Q(I.gen_two))

  @assert ideal(R, lift(R, gI)) + modulus == I
#    assert R*R!gi + Modulus(Rd) eq i;

  r = M_cur
  l = gI

  for i in 2:length(ideals)
    I = ideals[i]

    m_cur = zero_matrix(Q, n, m)

    RmodI, mRmodI = quo(R, I)
    MmodI = zero_matrix(RmodI, n, m)

    for i in 1:n
      for j in 1:m
        MmodI[i, j] = RmodI(lift(R, M[i, j]))
      end
    end

    echelon_modI = _strong_echelon_form_nonsplit(MmodI)

    for i in 1:n
      for j in 1:m
        m_cur[i, j] = Q(lift(R, echelon_modI[i, j]))
      end
    end

    _assure_weakly_normal_presentation(I)
    gI = gcd(Q(I.gen_one), Q(I.gen_two))

    @assert ideal(R, lift(R, gI)) + modulus == I

    g, a, b, e, f = xxgcd(l, gI)
    gg = g
    ginv = inv(g)
    g = g * ginv
    a = a * ginv
    b = b * ginv
    e = e * gg
    f = f * gg
    @assert g == a * l + b * gI
    @assert 0 == e * l + f * gI
    @assert 1 == a * f - b * e
    #@show inv(g)
    a = a * l
    b = b * gI
    #@assert isone(g)
    r = r * b + m_cur * a
    l = l * gI
  end
  return r
end

#    if l cmpeq 1 then
#      r := m;
#      l := gi;
#      L := i;
#    else
#      f, a, b := Idempotents(L, i);
#      if f then
#        a := Rd!a;
#        b := Rd!b;
#      else
#        g, a,b := Xgcd(l, gi);
#        a *:= l;
#        b *:= gi;
#        assert g eq 1;
#      end if;
#      assert 1 eq a+b;
#      r := r*b+m*a;
#      l *:= gi;
#      L *:= i;

function _strong_echelon_form_nonsplit(M)
  Q = base_ring(M)
  I = ideal(Q)

  n = nrows(M)
  m = ncols(M)

  if can_map_into_integer_quotient(Q)
    RmodIZ, f, g = map_into_integer_quotient(Q)
    M_cur = zero_matrix(Q, n, m)
    if can_make_small(RmodIZ)
      RmodIZsmall, ff, gg = make_small(RmodIZ)
      M_temp = zero_matrix(RmodIZsmall, n, m)
      for i in 1:n
        for j in 1:m
          M_temp[i, j] = ff(f(M[i, j]))
        end
      end
      strong_echelon_form!(M_temp)
      for i in 1:n
        for j in 1:m
          M_cur[i, j] = g(gg(M_temp[i, j]))
        end
      end
    else
      forflint = zero_matrix(FlintZZ, n, m)
      for i in 1:n
        for j in 1:m
          forflint[i, j] = f(M[i, j]).data
        end
      end
      ccall((:fmpz_mat_strong_echelon_form_mod, :libflint), Nothing, (Ref{fmpz_mat}, Ref{fmpz}), forflint, RmodIZ.modulus)
      for i in 1:n
        for j in 1:m
          M_cur[i, j] = Q(forflint[i, j])
        end
      end
      #error("dasds")
    end
    return M_cur
  else
    N = deepcopy(M)
    strong_echelon_form_naive!(N)
    return N
  end
end

function test_pseudohnf()
  Qx, x = FlintQQ["x"]
  for i in 2:15
    K, a = NumberField(x^i - 10, "a")
    O = maximal_order(K)
    lp = NfOrdFracIdl[]
    for p in [2, 3, 5, 7, 11, 13]
      pp = prime_decomposition(O, p)
      for P in pp
        push!(lp, frac_ideal(O, P[1]))
      end
    end

    pm = PseudoMatrix(matrix(K, 5, 5, [ rand(-2^10:2^10) for i in 1:25]), rand(lp, 5))

    @time d = numerator(det(pm))

    if iszero(norm(d))
      continue
    end

    q, w = z_split(d)
    R = order(d)
    ideals = q

    if length(w) != 0
      push!(ideals, prod(w))
    end

    @show length(ideals)

    gc()
    @time pseudo_hnf_mod(pm, d, strategy = :split)
    gc()
    @time pseudo_hnf_mod(pm, d, strategy = :no_split)
    gc()
    @time pseudo_hnf_kb(pm)
    gc()

    hnf_new = pseudo_hnf_mod(pm, d)
    hnf_old = pseudo_hnf_kb(pm)

    @assert det(hnf_new) == det(hnf_old)

    @assert Hecke._spans_subset_of_pseudohnf(hnf_new, hnf_old, :upperright)
  end
  println("PASS")
end


text:
# module TestRclNodeBase

using Rcl
import Rcl.internal: rcutils_error_is_set, rcutils_reset_error
using Test

@testset "node base" begin
    # initialize
    @test !rcl.ok()
    rcl.init()
    @test rcl.ok()

    @test !rcutils_error_is_set()
    nb = rcl.NodeBase("test_node")
    @test rcl.isvalid(nb)

    @test name(nb) == "test_node"
    @test namespace(nb) == "/"

    # finalizing it shouldn't cause errors
    finalize(nb)
    @test !rcutils_error_is_set()

    # But checking if the node is still valid apparently does set the error
    # TODO: file a bug about this with the rcl folks?
    @test !rcl.isvalid(nb)
    @test rcutils_error_is_set()
    # The error was expected so reset it
    rcutils_reset_error()
    @test !rcutils_error_is_set()

    @testset "name and namspace" begin
        foo_node = rcl.NodeBase("foo", "foo_ns")
        @test name(foo_node) == "foo"
        @test namespace(foo_node) == "/foo_ns"

        bar_node = rcl.NodeBase("bar", "/biz/baz")
        @test name(bar_node) == "bar"
        @test namespace(bar_node) == "/biz/baz"

        finalize(foo_node)
        finalize(bar_node)
    end

    # shutdown
    @test rcl.ok()
    rcl.shutdown()
    @test !rcl.ok()
end

# end # module TestRclBase
nothing


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    ContrationView, edgecontraction

struct ContractionView{T<:UndirectedGraph} <: UndirectedGraph
    graph::T
    collapsed::Set{Int}
    vnode::Int
    vneighbormap::Dict{Int, Dict{Int,Int}}
end

"""
    edgecontraction(graph::UndirectedGraph, nodes) -> ContractionView

Generate edge-contraction view.
"""
function edgecontraction(graph::UndirectedGraph, nodes)
    # Determine outgoing edges to be contracted
    incs = Set{Int}()
    for n in nodes
        union!(incs, incidences(graph, n))
    end
    inedges = edgeset(nodesubgraph(graph, nodes))
    outedges = setdiff(incs, inedges)

    # Virtual node
    vnode = maximum(nodeset(graph)) + 1

    # Neighbormap
    vnbrmap = Dict(vnode => Dict{Int,Int}())
    for e in outedges
        (u, v) = getedge(graph, e)
        if u in nodes
            if !(v in keys(vnbrmap))
                vnbrmap[v] = copy(neighbors(graph, v))
            end
            vnbrmap[v][e] = vnode
            vnbrmap[vnode][e] = v
        elseif v in nodes
            if !(u in keys(vnbrmap))
                vnbrmap[u] = copy(neighbors(graph, u))
            end
            vnbrmap[u][e] = vnode
            vnbrmap[vnode][e] = u
        else
            @assert false
        end
    end
    return ContractionView(graph, Set(nodes), vnode, vnbrmap)
end

function neighbors(view::ContractionView, idx::Int)
    idx in view.collapsed && throw(KeyError(idx))
    if idx in keys(view.vneighbormap)
        return view.vneighbormap[idx]
    else
        return neighbors(view.graph, idx)
    end
end

function getedge(view::ContractionView, idx::Int)
    (u, v) = getedge(view.graph, idx)
    if u in view.collapsed
        return (view.vnode, v)
    elseif v in view.collapsed
        return (u, view.vnode)
    else
        return (u, v)
    end
end

function hasedge(view::ContractionView, u::Int, v::Int)
    u in view.collapsed && throw(KeyError(idx))
    v in view.collapsed && throw(KeyError(idx))
    return findedgekey(view.graph, u, v) !== nothing
end

nodeset(view::ContractionView) = union(setdiff(nodeset(view.graph), view.collapsed), view.vnode)
edgeset(view::ContractionView) = edgeset(view.graph)

nodecount(view::ContractionView) = nodecount(view.graph) - length(view.collapsed) + 1
edgecount(view::ContractionView) = edgecount(view.graph)


text:
struct RandomBatchQs <: BatchPQs
    k::Int

    function RandomBatchQs(; k::Int)::RandomBatchQs
        (k < 1) && throw(ArgumentError("Invalid batch size k=$(k)."))
        return new(k)
    end
end

"""
Select batch by selecting batch_size indices from candidate_indices at random
"""
function select_batch(qs::RandomBatchQs, x::Array{T, 2}, labels::Dict{Symbol, Vector{Int}}, candidate_indices::Vector{Int})::Vector{Int} where T <: Real
    num_observations = length(candidate_indices)
    if num_observations <= qs.k
        return candidate_indices
    end

    return sample(candidate_indices, qs.k, replace=false)
end


text:
@everywhere function removeAllSingleSignatures(exposures_init, allSignatures, genome,)

    accr_first, frob_rel_div_first, norm_one_dif_first, exposures = evalSingleSample(exposures_init, allSignatures, genome)

    exposuresOutput = exposures;
    accr = accr_first;
    frob_rel_div = frob_rel_div_first;
    norm_one_dif =  norm_one_dif_first;

    numSig = count(!iszero, exposures);

    for j = 1 : ( numSig - 1)

        accr_temp, frob_rel_div_temp, norm_one_dif_temp, exposuresSample_temp, fID = removeOneSig(exposuresOutput, allSignatures, genome);

        if ((accr_first - accr_temp ) < 0.01)

            exposuresOutput = exposuresSample_temp;
            accr = accr_temp;
            frob_rel_div = frob_rel_div_temp;
            norm_one_dif =  norm_one_dif_temp;
        else
            break
        end
    end

    return exposuresOutput, accr, frob_rel_div, norm_one_dif
end

max_stars_repo_name: ypradat/sigprofilerjulia

max_issues_repo_name: ypradat/sigprofilerjulia

max_forks_repo_name: ypradat/sigprofilerjulia


text:
struct Snowflake
    n::UInt64
end

Snowflake(s::AbstractString) = Snowflake(parse(UInt64, s))
Base.show(io::IO, s::Snowflake) = print(io, string(s.n; base=10))
Base.:(==)(s::Snowflake, n::Integer) = s.n == n
Base.:(==)(n::Integer, s::Snowflake) = n == s.n
StructTypes.StructType(::Type{Snowflake}) = StructTypes.StringType()
HTTP.escapeuri(s::Snowflake) = escapeuri(s.n)

snowflake2datetime(s::Snowflake) = unix2datetime(((s.n >> 22) + 1420070400000) / 1000)
worker_id(s::Snowflake) = (s.n & 0x3e0000) >> 17
process_id(s::Snowflake) = (s.n & 0x1f000) >> 12
increment(s::Snowflake) = s.n & 0xfff