text:
using DSGE, ModelConstructors
using Distributions, Test, LinearAlgebra

m = AnSchorfheide()
# not moved below here
# prior
priordensity = exp(DSGE.prior(m))
@testset "Ensure prior density is a density" begin
    @test 0 <= priordensity <= 1
end

# Pseudo-measurement equation matrices in Systems
system = compute_system(m)
system[:ZZ_pseudo]
system[:DD_pseudo]

nothing

# Below have all been moved to ModelConstructors.jl
# Test Parameter type

# UnscaledParameter, fixed=false
#=α =  parameter(:α, 0.1596, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Normal(0.30, 0.05), fixed=false)
@testset "Test non-fixed UnscaledParameter" begin
    @test isa(α, UnscaledParameter)
    @test α.key == :α
    @test isa(α.prior.value, Normal)
    @test α.prior.value.μ == 0.3
    @test α.description == "No description available."
    @test α.tex_label == ""
    @test isa(α.transform, ModelConstructors.SquareRoot)
end

# UnscaledParameter, fixed = true
α_fixed =  parameter(:α_fixed, 0.1596, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.Untransformed(), Normal(0.30, 0.05), fixed=true)
@testset "Test fixed UnscaledParameter" begin
    @test α_fixed.transform_parameterization == (0.1596,0.1596)
    @test isa(α_fixed.transform, ModelConstructors.Untransformed)
end

# UnscaledParameter, fixed = true, transform should be overwritten given fixed
α_fixed =  parameter(:α_fixed, 0.1596, (1e-5, 0.999), (1e-5, 0.999), ModelConstructors.SquareRoot(), Normal(0.30, 0.05), fixed=true)
@testset "Test fixed UnscaledParameter, ensuring transform is overwritten" begin
    @test isa(α_fixed.transform, ModelConstructors.Untransformed)
end

# Fixed UnscaledParameter, minimal constructor
δ = parameter(:δ, 0.025)
@testset "Test fixed UnscaledParameter minimal constructor" begin
    @test δ.fixed
    @test δ.transform_parameterization == (0.025, 0.025)
    @test δ.valuebounds == (0.025, 0.025)
end

# Scaled parameter
β = parameter(:β, 0.1402, (1e-5, 10.), (1e-5, 10.), ModelConstructors.Exponential(), GammaAlt(0.25, 0.1), fixed=false,  scaling = x -> (1 + x/100)\1, description="β: Discount rate.", tex_label="\\beta ")
@testset "Test ScaledParameter constructor" begin
    @test isa(β, ScaledParameter)
    @test isa(β.prior.value, Gamma)
    @test isa(β.transform, ModelConstructors.Exponential)
end

# Invalid transform
@testset "Ensure error thrown on invalid transform" begin
    @test_throws UndefVarError α_bad = parameter(:α, 0.1596, (1e-5, 0.999), (1e-5, 0.999),
                                                  InvalidTransform(), Normal(0.30, 0.05), fixed=false)
end

# Arithmetic with parameters
@testset "Check arithmetic with parameters" begin
    @test promote_type(AbstractParameter{Float64}, Float16) == Float64
    @test promote_type(AbstractParameter{Float64}, Int8) == Float64
    ## @test promote_rule(AbstractParameter{Float64}, Float16) == Float64
    ## @test promote_rule(AbstractParameter{Float64}, Int8) == Float64
    @test δ + δ == 0.05
    @test δ^2 == 0.025^2
    @test -δ == -0.025
    @test log(δ) == log(0.025)
end

# transform_to_real_line and transform_to_model_space
cx = 2 * (α - 1/2)
@testset "Check parameter transformations for optimization" begin
    @test abs(transform_to_real_line(α) - cx / sqrt(1 - cx^2)) < .001
    @test transform_to_real_line(δ) == 0.025
end

m = AnSchorfheide()
let lastparam = parameter(:p, 0.0)
    for θ in m.parameters
        isa(θ, Parameter) && (lastparam = θ)
    end
    @testset "Check AnSchorfheide last parameter" begin
        @test isa(lastparam, Parameter)
        @test lastparam.value == 0.20*2.237937
    end
end
# transform_to_real_line and transform_to_model_space, acting on the entire parameter vector. they should be inverses!
pvec = m.parameters
vals = transform_to_real_line(pvec)
transform_to_model_space!(m, vals)
@testset "Check parameter transformations for optimization part 2" begin
    @test pvec == m.parameters
end

# all fixed parameters should be unchanged by both transform_to_real_line and transform_to_model_space
@testset "Check fixed parameters are unchanged by optimization transformations" begin
    for θ in m.parameters
        if θ.fixed
            @test θ.value == transform_to_real_line(θ)
            @test θ.value == transform_to_model_space(θ, θ.value)
        end
    end
end

# settings
# settings - boolean, string, and number. adding to model. overwriting. filestrings. testing/not testing.
n_mh_blocks = Setting(:n_mh_blocks, 22) # short constructor
reoptimize = Setting(:reoptimize, false)
vint = Setting(:data_vintage, "REF", true, "vint", "Date of data") # full constructor

@testset "Check settings corresponding to parameters" begin
    @test promote_rule(Setting{Float64}, Float16) == Float64
    @test promote_rule(Setting{Bool}, Bool) == Bool
    @test promote_rule(Setting{String}, String) == String
    @test convert(Int64, n_mh_blocks) == 22
    @test convert(String, vint) == "REF"

    @test get_setting(m, :n_mh_blocks) == m.settings[:n_mh_blocks].value
    m.testing = true
    @test get_setting(m, :n_mh_blocks) == m.test_settings[:n_mh_blocks].value
    @test ModelConstructors.filestring(m) == "_test"

    m.testing = false
    m <= Setting(:n_mh_blocks, 5, true, "mhbk", "Number of blocks for Metropolis-Hastings")
    @test m.settings[:n_mh_blocks].value == 5
    @test occursin(r"^\s*_mhbk=5_vint=(\d{6})", ModelConstructors.filestring(m))
    ModelConstructors.filestring(m, "key=val")
    ModelConstructors.filestring(m, ["key=val", "foo=bar"])
    m.testing = true

    # Overwriting settings
    a = gensym() # unlikely to clash
    b = gensym()
    m <= Setting(a, 0, true, "abcd", "a")
    m <= Setting(a, 1)
    @test m.test_settings[a].value == 1
    @test m.test_settings[a].print == true
    @test m.test_settings[a].code == "abcd"
    @test m.test_settings[a].description == "a"
    m <= Setting(b, 2, false, "", "b")
    m <= Setting(b, 3, true, "abcd", "b1")
    @test m.test_settings[b].value == 3
    @test m.test_settings[b].print == true
    @test m.test_settings[b].code == "abcd"
    @test m.test_settings[b].description == "b1"
end

# model paths. all this should work without errors
m.testing = true
addl_strings = ["foo=bar", "hat=head", "city=newyork"]
@testset "Check proper model paths" begin
    for fn in [:rawpath, :workpath, :tablespath, :figurespath]
        @eval $(fn)(m, "test")
        @eval $(fn)(m, "test", "temp")
        @eval $(fn)(m, "test", "temp", addl_strings)
    end
end=#


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1VolumeProjection <: SwaggerModel
    configMap::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ConfigMapProjection } # spec name: configMap
    downwardAPI::Any # spec type: Union{ Nothing, IoK8sApiCoreV1DownwardAPIProjection } # spec name: downwardAPI
    secret::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SecretProjection } # spec name: secret
    serviceAccountToken::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ServiceAccountTokenProjection } # spec name: serviceAccountToken

    function IoK8sApiCoreV1VolumeProjection(;configMap=nothing, downwardAPI=nothing, secret=nothing, serviceAccountToken=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1VolumeProjection, Symbol("configMap"), configMap)
        setfield!(o, Symbol("configMap"), configMap)
        validate_property(IoK8sApiCoreV1VolumeProjection, Symbol("downwardAPI"), downwardAPI)
        setfield!(o, Symbol("downwardAPI"), downwardAPI)
        validate_property(IoK8sApiCoreV1VolumeProjection, Symbol("secret"), secret)
        setfield!(o, Symbol("secret"), secret)
        validate_property(IoK8sApiCoreV1VolumeProjection, Symbol("serviceAccountToken"), serviceAccountToken)
        setfield!(o, Symbol("serviceAccountToken"), serviceAccountToken)
        o
    end
end # type IoK8sApiCoreV1VolumeProjection

const _property_map_IoK8sApiCoreV1VolumeProjection = Dict{Symbol,Symbol}(Symbol("configMap")=>Symbol("configMap"), Symbol("downwardAPI")=>Symbol("downwardAPI"), Symbol("secret")=>Symbol("secret"), Symbol("serviceAccountToken")=>Symbol("serviceAccountToken"))
const _property_types_IoK8sApiCoreV1VolumeProjection = Dict{Symbol,String}(Symbol("configMap")=>"IoK8sApiCoreV1ConfigMapProjection", Symbol("downwardAPI")=>"IoK8sApiCoreV1DownwardAPIProjection", Symbol("secret")=>"IoK8sApiCoreV1SecretProjection", Symbol("serviceAccountToken")=>"IoK8sApiCoreV1ServiceAccountTokenProjection")
Base.propertynames(::Type{ IoK8sApiCoreV1VolumeProjection }) = collect(keys(_property_map_IoK8sApiCoreV1VolumeProjection))
Swagger.property_type(::Type{ IoK8sApiCoreV1VolumeProjection }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1VolumeProjection[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1VolumeProjection }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1VolumeProjection[property_name]

function check_required(o::IoK8sApiCoreV1VolumeProjection)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1VolumeProjection }, name::Symbol, val)
end

ext: jl


text:
const PokemonIconSheetFile = datapath("sprites", "swsh-pokemon-icons.png")
const ItemIconSheetFile = datapath("sprites", "swsh-item-icons.png")
const GenderIconSheetFile = datapath("sprites", "gender-icons.png")
const TypeIconSheetFile = datapath("sprites", "type-icons.png")

macro _sheet_func(func_name, sheet_name, has_mask = false, preprocess = nothing)
    gray_name = Symbol(sheet_name, :Gray)
    file_name = Symbol(sheet_name, :File)
    TMask = has_mask ? Matrix{Float32} : Nothing
    TRGB = :(Nullable{SpriteSheet{RGB{Float32}, $TMask}})
    TGray = :(Nullable{SpriteSheet{Gray{Float32}, $TMask}})
    quote
        const $sheet_name = Ref{$TRGB}(nothing)
        const $gray_name = Ref{$TGray}(nothing)
        function $func_name(; force = false, gray = false)
            if gray
                @_load_data $gray_name force begin
                    sheet = SpriteSheet($file_name, gray = true)
                    $preprocess
                    sheet
                end
            else
                @_load_data $sheet_name force begin
                    sheet = SpriteSheet($file_name, gray = false)
                    $preprocess
                    sheet
                end
            end
        end
    end |> esc
end

@_sheet_func pokemon_icons PokemonIconSheet
@_sheet_func item_icons ItemIconSheet false begin
    sheet.image .= blur(sheet.image, 1f0)
end
@_sheet_func gender_icons GenderIconSheet true
@_sheet_func type_icons TypeIconSheet

function initialize_icons()
    pokemon_icons()
    item_icons()
    gender_icons()
    type_icons()
    nothing
end


text:
"""
    node_degrees(parent, [root = 1])

Return a vector telling the degree of each node.
"""
function node_degrees(parent::Vector{Int}, root = 1)
    @assert parent[root] == root
    idx = similar(parent)
    idx .= 0
    for p in parent
        idx[p] += 1
    end
    idx[root] -= 1      # root does not have itself as children
    return idx
end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
#!/usr/bin/env julia

#=
# data should be sorted and works by divide and conquer
# improved variant of binary search for equally distributed data
# it tries probing the position
#
# Interpolation search finds particular item by computing probe position.
# Start from middle of list. If not a match, divide list using probe position and find new middle.
# Do binary search in 2 blocks.
# Repeat
#
# run-time complexity is O(log log n)
=#

function interpolation_search(search_for, data_array, low=1, high=-1)
  if high == -1
    high = length(data_array)
  end
  if high < low
    return -1
  end

  ## calculation of probe position
  mid = Int(floor(
          (
           low +
           ((high - low) / (data_array[high] - data_array[low])) *
           (search_for - data_array[low])
          )
        ))
  if mid < low
    return -1
  end

  if data_array[mid] == search_for
    return mid
  end

  if data_array[mid] < search_for
    return interpolation_search(search_for, data_array, mid+1, high)
  else
    return interpolation_search(search_for, data_array, low, mid-1)
  end

  -1
end

     ## 1   2   3   4   5   6   7   8   9   10
data = [10, 11, 12, 22, 24, 27, 34, 47, 56, 67]

@time println("at index: ", interpolation_search(47, data))
@time println("at index: ", interpolation_search(12, data))
@time println("at index: ", interpolation_search(7, data))
@time println("at index: ", interpolation_search(67, data))
@time println("at index: ", interpolation_search(34, data))


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/09/2014

module AbstractPilotResponseImpl

export AbstractPilotResponse

abstract AbstractPilotResponse

end

module AbstractPilotResponseInterfaces

export updatePilotResponse

function updatePilotResponse()
end

end


text:
using Convex, SCS, ToeplitzMatrices, LinearAlgebra

function structHMirls(y,N,k,omega,sigma)
    # Line spectral estimator based on matrix completion via an iteratively
    # reweighted least-squares approach.
    # See: Kümmerle, Verdun '19: "Completion of Structured Low-Rank Matrices via
    # Iteratively Reweighted Least Squares".
    #
    # Inputs:
    # Y: Data matrix
    # N: ambient dimension \geq M
    # k: number of sinusoids
    # omega: indices of the observed entries of the data in the ambient space
    # sigma: noise variance of the additive white Gaussian noise with mean zero
    #
    # M = size(Y,1): number of measurements per snapshot
    # L = size(Y,2): number of snapshots
    #
    # Output:
    # f_hat: estimate of the frequencies of the K sinusoids in [rad]

    maxiters = 100
    #lambda = 1e-3
    lambda = sqrt(N*log(length(omega))*sigma)
    #lambda = 1e-1
    alpha = 0
    d1::Int = ceil(N/2)
    d2 = N + 1 - d1

    phi(y) = y[omega]

    function phiStar(z,N,omega)
        x = zeros(Complex{Float64},N)
        x[omega] = z
    end

    function truncSVD(A,k)
        svdA = svd(A)
        return svdA.U[:,1:k]*Diagonal(svdA.S[1:k])*svdA.Vt[1:k,:]
    end

    function kronSum(A,B)
        d1 = size(A)[1]
        d2 = size(B)[1]
        return kron(Matrix{Float64}(I, d2, d2),A) + kron(B,Matrix{Float64}(I, d1, d1))
    end

# Optimization variables
z = ComplexVariable(N)
H = ComplexVariable(d1,d2)
Hvec = ComplexVariable(d1*d2)

# initialize
z1 = phiStar(y,N,omega)
epsilon = svdvals(Hankel(z1[1:d1],z1[d1:end]))[1]
W = Matrix{Float64}(I, d1*d2, d1*d2)
error = Inf
ii = 0
z_new = Inf*ones(N)

# Define OP
constraints = [H[i,j] == z[i+j-1] for i = 1:d1, j = 1:d2]
constraints += [Hvec == vec(H)]
objective = abs(vec(H)'*W*vec(H)) + 1/(2*lambda)*sumsquares(z[omega]-y)
problem = minimize(objective, constraints)

while error > 1e-3 && ii < maxiters
ii += 1
z_old = z_new
if ii == 1
    problem1 = minimize(epsilon^2*sumsquares(z) + 1/(2*lambda)*sumsquares(z[omega]-y), constraints)
    solve!(problem1, SCS.Optimizer)
else
    objective = abs(2*Hvec'*W*Hvec) + 1/(2*lambda)*sumsquares(z[omega]-y)
    problem = minimize(objective, constraints)
    solve!(problem, SCS.Optimizer)
end
z_new = z.value
epsilon = min(epsilon, norm(z_new - z_old)+alpha^(2*k))
error = norm(z_new-z_old)/norm(z_new)
Hz = Hankel(z_new[1:d1],z_new[d1:end])
#Hz = Array{Complex{Float64},2}(undef,d1,d2)
#for i = 1:d1, j = 1:d2
#Hz[i,j] = z_new[i+j-1]
#end
W = inv(kronSum(truncSVD(Hz*Hz',k),truncSVD(Hz'*Hz,k)) + (epsilon^2)*Matrix{Float64}(I, d1*d2, d1*d2))
#println(error)
end

# Retrieve Frequencies
L = 1
clmnsY0 = setdiff(1:d2*L,(1:L).*d2)
clmnsY1 = clmnsY0.+1

Y0 = H.value[:,clmnsY0]
Y1 = H.value[:,clmnsY1]

F = svd(Y0)
p = sortperm(F.S,rev=true)[1:k]

f = eigvals(Diagonal(ones(K)./F.S[p])*F.U[:,p]'*Y1*F.V[:,p])
f_hat = (angle.(f)/(2*pi) .+ 1).%1
#println(f_hat)
#println(svdvals(H.value))
return f_hat

end


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

"The `Future` module implements future behavior of already existing functions,
which will replace the current version in a future release of Julia."
module Future

using Random

## copy!

# This has now been moved to Base (#29178), and should be deprecated in the
# next "deprecation phase".

"""
    Future.copy!(dst, src) -> dst

Copy `src` into `dst`.
This function has now been moved into `Base`, consider using `copy!(dst, src)` instead.
"""
copy!(dst::AbstractSet, src::AbstractSet) = Base.copy!(dst, src)
copy!(dst::AbstractDict, src::AbstractDict) = Base.copy!(dst, src)
copy!(dst::AbstractArray, src::AbstractArray) = Base.copy!(dst, src)

## randjump

"""
    randjump(r::MersenneTwister, steps::Integer) -> MersenneTwister

Create an initialized `MersenneTwister` object, whose state is moved forward
(without generating numbers) from `r` by `steps` steps.
One such step corresponds to the generation of two `Float64` numbers.
For each different value of `steps`, a large polynomial has to be generated internally.
One is already pre-computed for `steps=big(10)^20`.
"""
randjump(r::MersenneTwister, steps::Integer) =
    Random._randjump(r, Random.DSFMT.calc_jump(steps))

end # module Future

max_stars_repo_name: skumagai/julia

max_issues_repo_name: skumagai/julia

max_forks_repo_name: skumagai/julia


text:
struct EnvironmentSampler{E<:AbstractEnvironment,B<:TrajectoryBuffer}
    environments::Vector{E}
    buffers::Vector{B}
    function EnvironmentSampler(env_tconstructor; dtype::Maybe{Type} = nothing)
        nt = Threads.nthreads()
        envs = [env_tconstructor(nt)...]
        bufs = [TrajectoryBuffer(first(envs), dtype = dtype) for _ = 1:nt]
        new{eltype(envs),eltype(bufs)}(envs, bufs)
    end
end

function rollout!(
    policy!,
    B::TrajectoryBuffer,
    sampler::EnvironmentSampler,
    n::Integer;
    reset! = randreset!,
    Hmax::Integer = n,
    nthreads::Integer = Threads.nthreads(),
    truncate::Bool = true,
)
    n > 0 || argerror("n must be > 0")
    0 < Hmax <= n || argerror("Hmax must be in range (0, nsamples]")
    if !(0 < nthreads <= Threads.nthreads())
        argerror("nthreads must be in range (0, Threads.nthreads()]")
    end

    foreach(empty!, sampler.buffers)

    if nthreads == 1 # short circuit to avoid threading overhead
        _rollout(sampler, policy!, reset!, n, Hmax)
    else
        _threaded_rollout(sampler, policy!, reset!, n, Hmax, nthreads)
    end

    # TODO test
    ns = truncate ? n : sum(nsamples, sampler.buffers)
    collate!(B, sampler.buffers, ns)
    return B
end

function rollout(
    policy!,
    sampler::EnvironmentSampler,
    n::Integer;
    dtype::Maybe{DataType} = nothing,
    kwargs...,
)
    rollout!(
        policy!,
        TrajectoryBuffer(first(sampler.environments), dtype = dtype),
        sampler,
        n;
        kwargs...,
    )
end

function _rollout(
    sampler::EnvironmentSampler,
    policy!,
    reset!::R,
    n::Integer,
    Hmax::Integer,
) where {R}
    tid = Threads.threadid()
    env = sampler.environments[tid]
    B = sampler.buffers[tid]
    togo = n
    while togo > 0
        reset!(env)
        len = _rollout!(policy!, B, env, max(1, min(Hmax, togo)))
        togo -= len
    end
    return nothing
end

function _threaded_rollout(
    sampler::EnvironmentSampler,
    policy!,
    reset!,
    n::Integer,
    Hmax::Integer,
    nthreads::Integer,
)
    alive = Atomic{Bool}(true)
    iters = [Atomic{Int}(0) for _ = 1:nthreads]
    togo = Atomic{Int}(n)

    # TODO remove this once 1.3 is dropped
    @static if VERSION < v"1.4"
        Threads.@sync for i = 1:nthreads
            Threads.@spawn _thread_worker(
                sampler,
                policy!,
                reset!,
                togo,
                convert(Int, Hmax),
                iters,
                iters[i],
            )
        end
    else
        Threads.@sync for i = 1:nthreads
            Threads.@spawn _thread_worker(
                $sampler,
                $policy!,
                $reset!,
                $togo,
                $(convert(Int, Hmax)),
                $iters,
                $(iters[i]),
            )
        end
    end
    return nothing
end

function _thread_worker(
    sampler::EnvironmentSampler,
    policy!,
    @specialize(reset!),
    togo::Atomic{Int},
    Hmax::Int,
    iters::Vector{Atomic{Int}},
    iter::Atomic{Int},
)
    tid = Threads.threadid()
    env = sampler.environments[tid]
    B = sampler.buffers[tid]
    stopcb = () -> togo[] <= 0
    t = time()
    while togo[] > 0
        if iter[] <= minimum(getindex, iters)
            atomic_add!(iter, 1)
            reset!(env)
            len = _rollout!(policy!, B, env, max(1, min(togo[], Hmax)), stopcb)
            atomic_sub!(togo, len)
            t = time()
        else
            time() - t > 5 && internalerror("Timeout on thread $tid.")
            # See: https://github.com/JuliaLang/julia/issues/33097
            ccall(:jl_gc_safepoint, Cvoid, ())
        end
    end
    return nothing
end


text:
using DiffEqMonteCarlo, StochasticDiffEq, DiffEqBase,
      DiffEqProblemLibrary, OrdinaryDiffEq
using Test

using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
import DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_linear, prob_sde_2Dlinear

prob = prob_sde_linear
prob2 = MonteCarloProblem(prob)
sim = solve(prob2,SRIW1(),dt=1//2^(3),num_monte=10,adaptive=false)

m = timestep_mean(sim,3)
m2,v = timestep_meanvar(sim,3)
med = timestep_median(sim,3)
quant = timestep_quantile(sim,0.5,3)
@test quant ≈ med
@test m ≈ m2
m3,m4,c = timestep_meancov(sim,3,3)
@test m ≈ m3
@test v ≈ c
m3,m4,c = timestep_meancor(sim,3,3)
@test c ≈ one(c)
vecarr = timeseries_steps_mean(sim)
vecarr = timeseries_steps_median(sim)
vecarr = timeseries_steps_quantile(sim,0.5)
m_series,v_series = timeseries_steps_meanvar(sim)
summ = MonteCarloSummary(sim)

m4,v4 = m_series[3],v_series[3]
covar_mat = timeseries_steps_meancov(sim)[3,3]
@test m ≈ m4
@test v ≈ v4
@test m ≈ covar_mat[1]
@test m ≈ covar_mat[2]
@test v ≈ covar_mat[3]

@test (get_timestep(sim,1)...,) == (get_timepoint(sim,0.0)...,)
@test (get_timestep(sim,2)...,) == (get_timepoint(sim,1//2^(3))...,)
@test (get_timestep(sim,3)...,) == (get_timepoint(sim,1//2^(2))...,)

sim = solve(prob2,SRIW1(),dt=1//2^(3),num_monte=10)
m = timepoint_mean(sim,0.5)
med = timepoint_median(sim,0.6)
quant = timepoint_quantile(sim,0.5,0.6)
@test quant ≈ med
m2,v = timepoint_meanvar(sim,0.5)
@test m ≈ m2
m3,m4,c = timepoint_meancov(sim,0.5,0.5)
@test m ≈ m3
@test v ≈ c
m3,m4,c = timepoint_meancor(sim,0.5,0.5)
@test c ≈ one(c)
m_series = timeseries_point_mean(sim,0:1//2^(3):1)
m_series = timeseries_point_median(sim,0:1//2^(3):1)
m_series = timeseries_point_quantile(sim,0.5,0:1//2^(3):1)
m_series,v_series = timeseries_point_meanvar(sim,0:1//2^(3):1)
summ = MonteCarloSummary(sim,0:1//2^(3):1)
m5,v5 = m_series[5],v_series[5]
@test m ≈ m5
@test v ≈ v5
m6,m7,v6 = timeseries_point_meancov(sim,0:1//2^(3):1,0:1//2^(3):1)[5,5]
@test m ≈ m6
@test m ≈ m7
@test v ≈ v6

prob = prob_sde_2Dlinear
prob2 = MonteCarloProblem(prob)
sim = solve(prob2,SRIW1(),dt=1//2^(3),num_monte=10,adaptive=false)

m = timestep_mean(sim,3)
m2,v = timestep_meanvar(sim,3)
med = timestep_median(sim,3)
quant = timestep_quantile(sim,0.5,3)
@test quant ≈ med
@test m ≈ m2
m3,m4,c = timestep_meancov(sim,3,3)
@test m ≈ m3
@test v ≈ c
m3,m4,c = timestep_meancor(sim,3,3)
@test c ≈ ones(size(c)...)
vecarr = timeseries_steps_mean(sim)
vecarr = timeseries_steps_median(sim)
vecarr = timeseries_steps_quantile(sim,0.5)
m_series,v_series = timeseries_steps_meanvar(sim)
summ = MonteCarloSummary(sim)

m4,v4 = m_series[3],v_series[3]
covar_mat = timeseries_steps_meancov(sim)[3,3]
@test m ≈ m4
@test v ≈ v4
@test m ≈ covar_mat[1]
@test m ≈ covar_mat[2]
@test v ≈ covar_mat[3]

@test (get_timestep(sim,1)...,) == (get_timepoint(sim,0.0)...,)
@test (get_timestep(sim,2)...,) == (get_timepoint(sim,1//2^(3))...,)
@test (get_timestep(sim,3)...,) == (get_timepoint(sim,1//2^(2))...,)

sim = solve(prob2,SRIW1(),dt=1//2^(3),num_monte=10)
m = timepoint_mean(sim,0.5)
med = timepoint_median(sim,0.6)
quant = timepoint_quantile(sim,0.5,0.6)
@test quant ≈ med
m2,v = timepoint_meanvar(sim,0.5)
@test m ≈ m2
m3,m4,c = timepoint_meancov(sim,0.5,0.5)
@test m ≈ m3
@test v ≈ c
m3,m4,c = timepoint_meancor(sim,0.5,0.5)
@test c ≈ ones(size(c)...)
m_series = timeseries_point_mean(sim,0:1//2^(3):1)
m_series = timeseries_point_median(sim,0:1//2^(3):1)
m_series = timeseries_point_quantile(sim,0.5,0:1//2^(3):1)
m_series,v_series = timeseries_point_meanvar(sim,0:1//2^(3):1)
summ = MonteCarloSummary(sim,0:1//2^(3):1)
m5,v5 = m_series[5],v_series[5]
@test m ≈ m5
@test v ≈ v5
m6,m7,v6 = timeseries_point_meancov(sim,0:1//2^(3):1,0:1//2^(3):1)[5,5]
@test m ≈ m6
@test m ≈ m7
@test v ≈ v6


text:
include("population.jl")

function geneticAlgorithm(data, fitValue, repeatSize, selectSize, crossoverPoint, mutationPercentage)
    calcuatePopulationFitness!(data, fitValue)
    bestFit = data[1].fitness
    reapeatCounter = 1
    popGen = 0

    while (bestFit > 0 && reapeatCounter <= repeatSize)
        data = selectPopulation(data, selectSize)
        data = crossoverPopulation(data,crossoverPoint)
        mutatePopulation!(data,mutationPercentage)
        calculatePopulationFitness!(data,fitValue)

        if bestFit == data[1].fitness
            reapeatCounter +=1

        else
            bestFit = data[1].fitness
            reapeatCounter = 1
        end
        popGen +=1
    end
    return popGen, reapeatCounter, data
end

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
import CairoMakie as MK
import Clustering
import Colors
import Random
import Plots
import SparseArrays

using LinearAlgebra
using ProgressMeter
using Statistics
using StatsPlots

function assignment_summary_df(assignments::Pair{Symbol, Vector{Int}}...; min_molecules_per_cell::Int)
    assignments = Dict(assignments)
    return DataFrame(Dict(
        "name" => collect(keys(assignments)),
        "noise_frac" => [round(mean(x .== 0), sigdigits=3) for x in values(assignments)],
        "n_cells" => [sum(B.count_array(x, drop_zero=true) .>= min_molecules_per_cell) for x in values(assignments)]
    ))[:, [:name, :n_cells, :noise_frac]]
end

function plot_subset(df_spatial::DataFrame, dapi_arr::Union{Matrix, Nothing}, (xs, xe), (ys, ye); polygons::Union{Bool, Vector{Matrix{Float64}}}=true,
        markersize=2.0, alpha=0.2, min_molecules_per_cell::Int=1, grid_step::Float64=5.0, bandwidth::Float64=grid_step, min_border_length=3,
        cell_col::Symbol=:cell, dapi_alpha=0.9, polygon_line_width::T1 where T1 <: Real=2, dens_threshold::Float64=1e-5, noise::Bool=true,
        size_mult=1/3, plot_raw_dapi::Bool=true, color_col::Symbol=:ncv_color, annotation_col::Union{Symbol, Nothing}=nothing, build_panel::Bool=true,
        grid_alpha::Float64=0.5, ticks=false, swap_plots::Bool=false, dapi_color::Symbol=:twilight, clims=nothing, polygon_line_color="black",
        plot_bg_dapi::Bool=true, title::String="", kwargs...)
    df_subs = @where(df_spatial, :x .>= xs, :x .<= xe, :y .>= ys, :y .<= ye);

    if (typeof(polygons) == Bool)
        if polygons
            polygons = B.boundary_polygons(df_subs, df_subs[!, cell_col], grid_step=grid_step, min_molecules_per_cell=min_molecules_per_cell,
                bandwidth=bandwidth, min_border_length=min_border_length, dens_threshold=dens_threshold)
        else
            polygons = Matrix{Float64}[]
        end
    end

    # xs, xe, ys, ye = round.(Int, [minimum(df_subs.x), maximum(df_subs.x), minimum(df_subs.y), maximum(df_subs.y)]);
    xs, xe, ys, ye = round.(Int, [xs, xe, ys, ye]);

    xticks_vals = range(0, xe-xs, length=5)
    yticks_vals = range(0, ye-ys, length=5)

    yticks = xticks = Float64[]
    if ticks
        xticks = (xticks_vals, ["$f" for f in round.(Int, range(xs, xe, length=5))])
        yticks = (yticks_vals, ["$f" for f in round.(Int, range(ys, ye, length=5))])
    end

    plot_size = ((xe-xs), ye-ys) .* size_mult
    if plot_raw_dapi
        plot_size = (2 * plot_size[1], plot_size[2])
    end

    is_noise = noise ? (df_subs[!, cell_col] .== 0) : nothing
    annotation = (annotation_col === nothing) ? nothing : df_subs[!, annotation_col]

    polygon_kwargs = Dict(:strokewidth => polygon_line_width, :strokecolor => polygon_line_color)
    if dapi_arr === nothing
        return B.plot_molecules(df_subs, polygons; color=color_col, markersize=markersize, alpha=alpha, offset=(-xs, -ys), size=plot_size,
            is_noise=is_noise, polygon_kwargs=polygon_kwargs, noise_kwargs=Dict(:markersize => markersize), annotation=annotation, kwargs...)
    end

    dapi_subs = copy(dapi_arr[ys:ye, xs:xe]')
    plot_args = Dict(:size => plot_size, :xticks => xticks, :yticks => yticks)

    fig = MK.Figure(resolution=plot_size)
    get_axis() = ticks ?
        MK.Axis(fig, xticks=xticks, yticks=yticks, title=title) :
        MK.Axis(fig, xticklabelsvisible=false, xticksvisible=false, yticklabelsvisible=false, yticksvisible=false, title=title)

    pi1, pi2 = swap_plots ? (2, 1) : (1, 2)
    fig[1, pi1] = get_axis()

    if plot_bg_dapi
        MK.image!(Colors.GrayA.(1 .- dapi_subs ./ maximum(dapi_subs), dapi_alpha), interpolate=false)
    end

    B.plot_molecules!(df_subs, polygons; color=color_col, markersize=markersize, alpha=alpha, offset=(-xs, -ys),
        polygon_kwargs=polygon_kwargs, is_noise=is_noise,
        noise_kwargs=Dict(:markersize => markersize), annotation=annotation, kwargs...)

    if !plot_raw_dapi
        return fig
    end

    # Possible colorschemes: :delta, :twilight, :PuBuGn_9, :PuBu_9, :dense, :tofino, :berlin, :delta
    fig[1, pi2] = get_axis()
    hm = MK.heatmap!(dapi_subs, colormap=dapi_color)
    MK.poly!([MK.Point2.(eachrow(p .- [xs ys])) for p in polygons]; color="transparent", polygon_kwargs...)

    if clims !== nothing
        hm.colorrange = clims
    end

    if !build_panel
        error("build_panel is not supported")
        # return plt1, plt2
    end

    return fig
end

function cell_coord_frame(df_spatial::DataFrame, cell_id::Int; cell_col::Symbol=:cell, offset::Float64=0.1)
    sample_df = df_spatial[df_spatial[!, cell_col] .== cell_id,:];
    return [round.(Int, ((1. + offset) * s - offset * e, (1. + offset) * e - offset * s)) for (s,e) in [B.val_range(sample_df.x), B.val_range(sample_df.y)]];
end

function plot_comparison_for_cell(df_spatial::DataFrame, cell_id::Int, args...; cell1_col::Symbol=:cell, cell2_col::Symbol=:cell_paper, offset::Float64=0.1, kwargs...)
    sample_df = df_spatial[df_spatial[!, cell1_col] .== cell_id,:];
    paper_ids = setdiff(unique(sample_df[!, cell2_col]), [0]);

    if !isempty(paper_ids)
        sample_df = df_spatial[(df_spatial[!, cell1_col] .== cell_id) .| in.(df_spatial[!, cell2_col], Ref(paper_ids)),:];
    end;

    xc, yc = median(sample_df.x[sample_df[!, cell1_col] .== cell_id]), median(sample_df.y[sample_df[!, cell1_col] .== cell_id])
    xls, yls = [round.(Int, ((1. + offset) * s - offset * e, (1. + offset) * e - offset * s)) for (s,e) in [B.val_range(sample_df.x), B.val_range(sample_df.y)]];
    return plot_comparison_for_cell(df_spatial, xls, yls, args...; xc=xc, yc=yc, cell_col=cell1_col, kwargs...)
end

function plot_comparison_for_cell(df_spatial::DataFrame, xls::Tuple{T, T}, yls::Tuple{T, T}, seg_arr::Union{Matrix{<:Integer}, Nothing},
        dapi_arr::Union{Matrix{<:Real}, Nothing}; paper_polys::Array{Matrix{Float64}, 1}=Matrix{Float64}[], polygon_line_width::Float64=2.0, polygon_alpha::Float64=0.5,
        size_mult::Float64=1.0, grid_alpha::Float64=0.0, markersize::Float64=2.0, title="", center_mult::Float64=3.0, noise::Bool=false, paper_poly_color="darkred", paper_line_mult::Float64=1.0,
        plot_raw_dapi::Bool=true, kwargs...) where T <: Real

    xls, yls = max.(xls, 1), max.(yls, 1)

    if dapi_arr !== nothing
        xls = min.(xls, size(dapi_arr, 2))
        yls = min.(yls, size(dapi_arr, 1))
    end

    if seg_arr !== nothing
        paper_polys = extract_polygons_from_label_grid(copy(seg_arr[yls[1]:yls[2], xls[1]:xls[2]]'))
        # paper_polys = [Plots.Shape(pg[:,1], pg[:,2]) for pg in paper_polys]
    else
        paper_polys = [pg .- [xls[1] yls[1]] for pg in paper_polys]
    end

    if dapi_arr === nothing
        plot_raw_dapi = false
    end

    fig = plot_subset(df_spatial, dapi_arr, xls, yls; size_mult=size_mult, grid_alpha=grid_alpha, markersize=markersize, noise=noise,
        polygon_line_width=polygon_line_width, polygon_alpha=polygon_alpha, plot_raw_dapi=plot_raw_dapi, kwargs...);

    if length(paper_polys) > 0
        for i in 1:length(fig.scene.children)
            MK.poly!(fig[1, i], [MK.Point2.(eachrow(p)) for p in paper_polys]; color="transparent", strokecolor=paper_poly_color, strokewidth=polygon_line_width*paper_line_mult)
        end;
    end

    return fig
end

function joint_ordering(matrices::Matrix{T}...) where T <: Real
    joint_mtx = hcat(matrices...)
    cell_dists = Symmetric(1 .- cor(joint_mtx));
    cell_ord = Clustering.hclust(cell_dists, linkage=:ward).order;

    cell_ords = [Int[] for m in matrices]
    cum_sum = 0
    for i in 1:length(matrices)
        cell_ords[i] = cell_ord[(cell_ord .<= (cum_sum + size(matrices[i], 2))) .& (cell_ord .> cum_sum)] .- cum_sum
        cum_sum += size(matrices[i], 2)
    end

    gene_dists = Symmetric(1 .- cor(joint_mtx'));
    gene_ord = Clustering.hclust(gene_dists, linkage=:ward).order;

    return cell_ords, gene_ord
end

## Comparison of segmentations

prepare_qc_df(df_spatial::DataFrame, cell_col::Symbol=:cell; kwargs...) =
    prepare_qc_df(df_spatial, Int.(df_spatial[!, cell_col]); kwargs...)

function prepare_qc_df(df_spatial::DataFrame, assignment::Vector{<:Int}; min_area::T where T<:Real, min_molecules_per_cell::T2 where T2 <: Real,
        max_molecules_per_cell::Int=0, max_elongation::T3 where T3 <: Real = 30.0, dapi_arr::Union{Matrix{<:Real}, Nothing}=nothing)
    qc_per_cell = B.get_cell_qc_df(df_spatial, assignment, dapi_arr=dapi_arr);
    qc_per_cell[!, :cell_id] = 1:size(qc_per_cell, 1)
    qc_per_cell[!, :sqr_area] = sqrt.(qc_per_cell.area)

    if max_molecules_per_cell == 0
        max_molecules_per_cell = maximum(qc_per_cell.n_transcripts)
    end

    return @where(qc_per_cell, :n_transcripts .>= min_molecules_per_cell, :n_transcripts .<= max_molecules_per_cell,
        :sqr_area .>= min_area, :elongation .<= max_elongation)
end

hist_bins(vals::Vector{<:Real}...; n_bins::Int=100, min_val::T where T<: Real=0.0, m_quantile::T2 where T2<:Real=0.99, max_val_mult::Float64=1.0) =
    range(min_val, maximum([quantile(v, m_quantile) / m_quantile for v in vals if length(v) > 0])*max_val_mult, length=n_bins)

function plot_qc_comparison(qc_per_cell_dfs::Union{Array{DataFrame, 1}, Tuple{DataFrame, DataFrame}}; labels::Vector{String}=["Baysor", "DAPI"],
        max_quants::Vector{<:Real}=[0.9999, 0.99, 0.99, 0.999, 0.999], n_bins::Vector{<:Real}=[75, 100, 100, 100, 100], kwargs...)
    m_names = [:n_transcripts, :density, :elongation, :sqr_area, :mean_dapi_brightness]
    plot_titles = ["Num. of transcripts", "Density", "Elongation", "sqrt(Area)", "Mean DAPI brightness"]
    plots = Plots.Plot[]
    for (cs, xlab, mq, nb) in zip(m_names, plot_titles, max_quants, n_bins)
        if any(!(cs in propertynames(qdf)) for qdf in qc_per_cell_dfs)
            continue
        end
        t_bins = hist_bins([qdf[!, cs] for qdf in qc_per_cell_dfs]..., m_quantile=mq, n_bins=nb)
        plt = Plots.plot(widen=false, xlabel=xlab, ylabel="Num. of cells", xlims=B.val_range(t_bins))
        for (qdf, lab) in zip(qc_per_cell_dfs, labels)
            if length(unique(qdf[!, cs])) < 2
                continue
            end
            Plots.histogram!(qdf[!, cs], label=lab, bins=t_bins, alpha=0.6, kwargs...)
        end
        push!(plots, plt)
    end

    return Plots.plot(plots..., size=(400 * ceil(Int, length(plots) / 2), 600))
end

function match_assignments(qc_per_cell_dfs::Union{Array{DataFrame, 1}, Tuple{DataFrame, DataFrame}}, df_spatial::DataFrame, cell_cols::Vector{Symbol}; kwargs...)
    assignments = [denserank(ifelse.(in.(df_spatial[!,s], Ref(Set(qdf.cell_id))), df_spatial[!,s], 0)) .- 1 for (s, qdf) in zip(cell_cols, qc_per_cell_dfs)];
    return match_assignments(assignments...; kwargs...)
end

function match_assignments(assignment1::Vector{<:Integer}, assignment2::Vector{<:Integer}; bin_match_frac::Float64=0.05)
    # Basic statistics
    contingency_table = sparse_counts(assignment1 .+ 1, assignment2 .+ 1);
    ci,ri,v = SparseArrays.findnz(contingency_table)
    m_norm1 = sparse(ci, ri, v ./ sum(contingency_table, dims=1)[ri])
    m_norm2 = sparse(ci, ri, v ./ sum(contingency_table, dims=2)[ci])

    noise_fracs = Vector.((m_norm2[2:end,1], m_norm1[1,2:end]));
    max_overlaps = Vector.((maximum(m_norm2[2:end,:], dims=2)[:], maximum(m_norm1[:,2:end], dims=1)[:]));

    bin_match = dropzeros!(sparse(ci, ri, m_norm1.nzval .>= bin_match_frac))[2:end, 2:end];
    n_overlaps = (sum(bin_match, dims=2)[:], sum(bin_match, dims=1)[:])

    ctn_min = sparse(ci, ri, min.(m_norm1.nzval, m_norm2.nzval))
    match_fracs = range(0.0, 1.0, length=30)
    match_nums = [sum(any(ctn_min .>= mt, dims=1)) for mt in match_fracs];

    # Derivative statistics
    match_noise = [nf .> 0.5 for nf in noise_fracs]
    multiple_overlap = [((mo .< 0.7) .& (nf .< 0.25)) for (mo, nf) in zip(max_overlaps, noise_fracs)];
    @assert all([!any(mn .& mo) for (mn, mo) in zip(match_noise, multiple_overlap)])

    return (contingency=contingency_table, noise_fracs=noise_fracs, max_overlaps=max_overlaps, n_overlaps=n_overlaps, match_fracs=match_fracs,
        match_nums=match_nums, match_noise=match_noise, multiple_overlap=multiple_overlap)
end

function plot_matching_comparison(match_results::NamedTuple; labels::Vector{String}=["Baysor", "DAPI"], n_bins::Vector{<:Real}=[50, 30, 30], kwargs...)
    m_names = [:noise_fracs, :max_overlaps]
    plot_titles = ["Fraction of molecules, not assigned in other segmentation", "Fraction of molecules, matching to a single cell"]
    plots = Plots.Plot[]
    for (n, xlab, nb) in zip(m_names, plot_titles, n_bins)
        vals = match_results[n]
        t_bins = hist_bins(vals..., m_quantile=1.0, n_bins=nb, max_val_mult=1.02)
        plt = Plots.histogram(vals[1], label=labels[1], bins=t_bins, xlims=B.val_range(t_bins), widen=false,
            xlabel=xlab, ylabel="Num. of cells", kwargs...)
        Plots.histogram!(vals[2], label=labels[2], bins=t_bins, alpha=0.6)
        push!(plots, plt)
    end

    max_bin = maximum(maximum.(match_results.n_overlaps))
    plt = Plots.histogram(match_results.n_overlaps[1], label=labels[1], bins=0:0.5:max_bin, xticks=0:max_bin, widen=false,
        xlabel="Number of overlapping cells", ylabel="Num. of cells", kwargs...)
    Plots.histogram!(match_results.n_overlaps[2], label=labels[2], bins=0:0.5:max_bin, alpha=0.6)
    push!(plots, plt)

    plt = Plots.plot(match_results.match_fracs, match_results.match_nums, xlabel="Minimal fraction of matching molecules", ylabel="Num. of matching cells",
        xlims=(0.0, 1.0), ylims=(0, maximum(match_results.match_nums) * 1.05), widen=false, legend=:none, lw=2.0)
    push!(plots, plt)

    return Plots.plot(plots..., layout=(2, 2), size=(900, 600))
end

function build_statistics_df(qc_per_cell_dfs::Union{Array{DataFrame, 1}, Tuple{DataFrame, DataFrame}}, match_res::NamedTuple, df_spatial::DataFrame; labels::Vector{Symbol}=[:Baysor, :DAPI], sigdigits::Int=3)
    stat_df = DataFrame(Dict(:Metric => ["Num. cells", "Total num. molecules", "Fraction of noise", "Fraction of matching cells", "Fraction of matching to noise", "Fraction of multiple overlap"]))
    for i in 1:length(qc_per_cell_dfs)
        stat_df[!, labels[i]] = Any[
            size(qc_per_cell_dfs[i], 1), sum(qc_per_cell_dfs[i].n_transcripts), round(1 .- sum(qc_per_cell_dfs[i].n_transcripts) / size(df_spatial, 1), sigdigits=sigdigits),
            round(mean(.!match_res.match_noise[i] .& (.!match_res.multiple_overlap[i])), sigdigits=sigdigits), round(mean(match_res.match_noise[i]), sigdigits=sigdigits),
            round(mean(match_res.multiple_overlap[i]), sigdigits=sigdigits)
        ]
    end

    return stat_df
end

estimate_embedding(df_spatial::DataFrame, qc_per_cell_dfs::Union{Array{DataFrame, 1}, Tuple{DataFrame, DataFrame}}, cell_cols::Vector{Symbol}; kwargs...) =
    estimate_embedding([convert_segmentation_to_counts(df_spatial.gene, df_spatial[!, cq])[:, qdf.cell_id] for (cq, qdf) in zip(cell_cols, qc_per_cell_dfs)]...; kwargs...)

function estimate_embedding(count_matrices::Matrix{<:Real}...; joint::Bool=true, kwargs...)
    cm_merged = hcat(count_matrices...);
    cm_merged = cm_merged ./ sum(cm_merged, dims=1)
    ids_per_mat = B.split_ids(vcat([repeat([i], inner=size(cm, 2)) for (i,cm) in enumerate(count_matrices)]...))
    umap_merged = fit(B.UmapFit, cm_merged; kwargs...);
    return [umap_merged.embedding[:, ids] for ids in ids_per_mat], umap_merged, cm_merged
end

function plot_qc_embeddings(qc_per_cell_dfs::Union{Array{DataFrame, 1}, Tuple{DataFrame, DataFrame}}, match_res::NamedTuple, embeddings::Array{Matrix{Float64}, 1}; log_colors::Bool=true, size=(1000, 1000),
        labels::Vector{String}=["Baysor", "DAPI"], legend::Symbol=:best, ms::Float64=3.0, alpha1::Float64=0.25, plot_num_transcripts::Bool=true, layout=nothing,
        plot_partial::Bool=false, palette=Plots.distinguishable_colors(4, cchoices=[80, 100]), build_figure::Bool=true, kwargs...)
    plts = Plots.Plot[]
    c_vals = [qdf.n_transcripts for qdf in qc_per_cell_dfs]
    if log_colors
        c_vals = [log.(cv) for cv in c_vals]
    end

    factor_per_cell = String[]

    if plot_partial
        factor_per_cell = [ifelse.(match_res.multiple_overlap[i], "Multiple", ifelse.(match_res.match_noise[i], "No",
            ifelse.(match_res.max_overlaps[i] .< 0.7, "Partial", "Full"))) for i in 1:length(qc_per_cell_dfs)];
    else
        factor_per_cell = [ifelse.(match_res.multiple_overlap[i], "Multiple", ifelse.(match_res.match_noise[i], "No", "Yes")) for i in 1:length(qc_per_cell_dfs)];
    end

    for i in 1:length(qc_per_cell_dfs)
        plt = Plots.plot(format=:png, size=(600, 600), legend=legend, bg_legend=Colors.RGBA(1.0, 1.0, 1.0, 0.5), title="$(labels[i]), matching", legend_title="Match")
        for (ci,ids) in enumerate(B.split_ids(denserank(factor_per_cell[i])))
            Plots.scatter!(embeddings[i][1,ids], embeddings[i][2,ids], label=factor_per_cell[i][ids[1]], ms=ms, alpha=alpha1, markerstrokewidth=0; color=palette[ci], kwargs...)
        end
        push!(plts, plt)

        plot_num_transcripts || continue

        colors = map_to_colors(c_vals[i], lims=B.val_range(vcat(c_vals...)))[:colors]
        plt = Plots.scatter(embeddings[i][1,:], embeddings[i][2,:], color=colors, ms=ms, markerstrokewidth=0, format=:png, size=(600, 600), legend=false,
            title="$(labels[i]), num. transcripts")
        push!(plts, plt)
    end;

    if !build_figure
        return plts
    end

    return Plots.plot(plts..., size=size, format=:png, layout=something(layout, length(plts)))
end

plot_expression_vec_comparison(df_spatial::DataFrame, qc_per_cell_dfs::Dict{Symbol, DataFrame}, gene_names; labels=["Baysor", "DAPI"], kwargs...) =
    B.plot_expression_vectors([B.count_array(vcat(B.split(df_spatial.gene, df_spatial[!, cs] .+ 1)[2:end][qdf.cell_id]...)) for (cs, qdf) in qc_per_cell_dfs]...; gene_names=gene_names, labels=labels, kwargs...)

plot_expression_vec_comparison(df_spatial::DataFrame, qc_per_cell_dfs::Union{Array{DataFrame, 1}, Tuple{DataFrame, DataFrame}}, cell_cols::Vector{Symbol}, gene_names; labels=["Baysor", "DAPI"], kwargs...) =
    B.plot_expression_vectors([B.count_array(vcat(B.split(df_spatial.gene, df_spatial[!, cs] .+ 1)[2:end][qdf.cell_id]...)) for (cs, qdf) in zip(cell_cols, qc_per_cell_dfs)]...; gene_names=gene_names, labels=labels, kwargs...)

function estimate_non_matching_part_correlation(df_spatial::DataFrame, qc_per_cell_dfs::Union{Array{DataFrame, 1}, Tuple{DataFrame, DataFrame}},
        match_res::NamedTuple; cell_cols::Vector{Symbol}=[:cell, :cell_paper], rev::Bool=false, max_overlap_borders::Tuple{Float64, Float64}=(0.25, 0.75))
    data_id = rev ? 2 : 1
    if rev
        cell_cols = reverse(cell_cols)
    end
    part_overlap_mask = (match_res.max_overlaps[data_id] .> max_overlap_borders[1]) .& (match_res.max_overlaps[data_id] .< max_overlap_borders[2]);
    cur_overlaps = match_res.max_overlaps[data_id][part_overlap_mask]
    partially_matching_ids = qc_per_cell_dfs[data_id].cell_id[part_overlap_mask];
    match_cors = Float64[]

    n_genes = maximum(df_spatial.gene)
    for t_cell_id in partially_matching_ids
        t_df = df_spatial[df_spatial[!, cell_cols[1]] .== t_cell_id,:];
        t_mcp = mode(t_df[!, cell_cols[2]][t_df[!, cell_cols[2]] .!= 0]);

        match_expr = B.prob_array(t_df.gene[t_df[!, cell_cols[2]] .== t_mcp], max_value=n_genes);
        non_match_expr = B.prob_array(t_df.gene[t_df[!, cell_cols[2]] .!= t_mcp], max_value=n_genes);

        push!(match_cors, cor(non_match_expr, match_expr))
    end

    return match_cors, partially_matching_ids, cur_overlaps, part_overlap_mask
end

function correlation_effect_size_df(datasets::NamedTuple, part_cor_key::Symbol, cell_col_labels::Vector{String}, cell_col_names::Vector{Symbol})
    dfs = DataFrame[]
    for i in 1:2
        for d in datasets
            (part_cor_key in keys(d)) || continue

            corr_res = d[part_cor_key][i]
            qc_df = d[:qc_per_cell_dfs][cell_col_names[i]]
            corr = corr_res[1];
            ord = sortperm(corr);
            corr = corr[ord]

            mismatch_frac = 1 .- corr_res[3][ord];
            n_trans_per_cell = qc_df.n_transcripts[corr_res[4]][ord];

            push!(dfs, DataFrame(
                :Correlation => corr,
                :MolFrac => cumsum(mismatch_frac .* n_trans_per_cell) ./ sum(qc_df.n_transcripts),
                :CellFrac => (1:length(corr)) ./ size(qc_df, 1),
                :Dataset => d[:name],
                :Segmentation => cell_col_labels[i]
            ))
        end
    end
    return vcat(dfs...)
end

function filter_assignment(assignment::Vector{Int}, min_mols_per_cell::Int)
    passed_ids = Set(findall(B.count_array(assignment, drop_zero=true) .>= min_mols_per_cell))
    return denserank([((a in passed_ids) ? a : 0) for a in assignment]) .- 1
end

function append_matching_statistics!(d::Dict{Symbol, Any}, cell_col_names::Vector{Symbol};
        target_cell_col_names::Union{Vector{Symbol}, Nothing}=nothing)
    cur_names = sort(intersect(cell_col_names, propertynames(d[:df])))
    if length(setdiff(cell_col_names, cur_names)) > 0
        @warn "The following col names weren't found in the data: $(setdiff(cell_col_names, cur_names))"
    end

    if target_cell_col_names === nothing
        target_cell_col_names = cur_names
    end

    d[:qc_per_cell_dfs] = Dict(k => prepare_qc_df(d[:df], k; min_area=d[:min_area],
        min_molecules_per_cell=d[:min_mols_per_cell], max_molecules_per_cell=get(d, :max_mols_per_cell, 0), max_elongation=15) for k in cur_names);

    for cs in cur_names # filter cells
        d[:df][!, cs][.!in.(d[:df][!, cs], Ref(Set(d[:qc_per_cell_dfs][cs].cell_id)))] .= 0
        any(d[:df][!, cs] .!= 0) || error("All cells were filtered out for '$cs'")
    end

    d[:assignment_filt] = Dict(cs => denserank(d[:df][!, cs]) .- 1 for cs in cur_names)

    # for (i1,c1) in enumerate(cur_names)
        # for c2 in cur_names[(i1 + 1):length(cur_names)]
    @showprogress for c1 in cur_names
        for c2 in target_cell_col_names
            (c1 != c2) || continue
            parts = [get(split(String(c1), '_'), 2, "cell"), get(split(String(c2), '_'), 2, "cell")]
            parts = [p for p in parts if length(p) > 0]
            sm = Symbol(join(vcat(["match_res"], parts), "_"))
            d[sm] = match_assignments(d[:assignment_filt][c1], d[:assignment_filt][c2]);

            sel_cells = [c1, c2]
            spc = Symbol(join(vcat(["part_cors"], parts), "_"))
            qc_dfs = [d[:qc_per_cell_dfs][k] for k in sel_cells]
            d[spc] = [estimate_non_matching_part_correlation(d[:df], qc_dfs, d[sm], rev=r, cell_cols=sel_cells) for r in [false, true]]

            ssd = Symbol(join(vcat(["stat_df"], parts), "_"))
            d[ssd] = build_statistics_df([d[:qc_per_cell_dfs][k] for k in sel_cells], d[sm], d[:df]; labels=sel_cells)

            sfv = Symbol(join(vcat(["frac_vals"], parts), "_"))
            overlap_mask = d[sm].n_overlaps[1] .== 1
            overlap_nums, uniq_overlap_ids = vec.(collect.(findmax(d[sm].contingency[2:end,:][overlap_mask, 2:end], dims=2)));
            d[sfv] = (
                overlap = overlap_nums ./ d[:qc_per_cell_dfs][c2].n_transcripts[getindex.(uniq_overlap_ids, 2)],
                assigned = d[:qc_per_cell_dfs][c1].n_transcripts[overlap_mask] ./ d[:qc_per_cell_dfs][c2].n_transcripts[getindex.(uniq_overlap_ids, 2)],
                overlap_ids = findall(overlap_mask)
            )
        end
    end
end

## Quick visualizations

function plot_correlation(d::Dict{Symbol, Any}, part_cor_key::Symbol)
    pdfs1 = DataFrame(:Correlation => d[part_cor_key][1][1], :Type => d[:name]);
    pdfs2 = DataFrame(:Correlation => d[part_cor_key][2][1], :Type => d[:name]);
    t_heights = (0.9 .* (size(pdfs1, 1), size(pdfs2, 1)) ./ max(size(pdfs1, 1), size(pdfs2, 1)))

    display(("Baysor", sum(pdfs1.Correlation .> 0.5), sum(pdfs1.Correlation .< 0.5)))
    display(("Paper", sum(pdfs2.Correlation .> 0.5), sum(pdfs2.Correlation .< 0.5)))

    plt = @df pdfs1 violin(:Type, :Correlation, side=:right, label="Baysor", ylabel="Correlation between parts", legendtitle="Source",
        bg_legend=Colors.RGBA(1.0, 1.0, 1.0, 0.9), bar_width=t_heights[1], size=(300, 240), legend=:bottomright, xgrid=false) # , color=color_per_col_label["Baysor"]
    @df pdfs2 violin!(:Type, :Correlation, side=:left, label="Paper", bar_width=t_heights[2]) # , color=color_per_col_label["Paper"]

    return plt
end

function plot_correlation_lines(d::Dict{Symbol, Any}, part_cor_key::Symbol; lw=3, max_cor=0.5)
    vals = sort(d[part_cor_key][1][1])
    vals = vals[vals .<= max_cor]
    plt = Plots.plot(vals, 1:length(vals), label="Baysor", size=(300, 200), # , color=color_per_col_label["Baysor"]
        xlabel="Correlation", ylabel="Num. cells", legend=:topleft, lw=lw)
    vals = sort(d[part_cor_key][2][1])
    vals = vals[vals .<= max_cor]
    Plots.plot!(vals, 1:length(vals), label="Paper", lw=lw) # , color=color_per_col_label["Paper"]
    return plt
end

function plot_matching_hists(frac_vals::NamedTuple, match_res::NamedTuple)
    bins = 0.0:0.02:1.0
    plt1 = Plots.histogram(frac_vals.overlap, bins=bins, label=false,
        xlabel="Overlap size / poly-a cell size", ylabel="Cell density", normalize=true, lw=0, xlim=(0, 1))
    Plots.stephist!(frac_vals.overlap, bins=bins, label="", normalize=true, color="black")
    Plots.vline!([0.9975], label="", color="black", lw=2.5);

    bins = 0.0:0.04:2.0
    plt2 = Plots.histogram(min.(frac_vals.assigned, 1.99), bins=bins, normalize=true, lw=0,
        legend=false, xlabel="Target cell size / poly-A cell size", ylabel="Cell density")
    Plots.stephist!(min.(frac_vals.assigned, 1.99), bins=bins, label="", normalize=true, color="black")
    Plots.vline!([1], label="", color="black", lw=2.5);

    vals = B.prob_array(match_res.n_overlaps[1] .+ 1);
    xvals = 0:(length(vals) - 1)
    plt3 = Plots.bar(xvals, vals, xlabel="Num. overlaps", ylabel="Fraction", label=false, xticks=xvals, xgrid=false)

    plt = Plots.plot(plt1, plt2, plt3, size=(800, 200), layout=(1, 3))
    return plt
end

## Matching benchmarks

function estimate_jaccard_per_molecule(df::DataFrame, cell_col1::Symbol, cell_col2::Symbol; min_mols_per_cell::Int)
    mol_ids_per_cell = [Set.(B.split_ids(df[!, k] .+ 1)[2:end]) for k in [cell_col1, cell_col2]]
    cms = [B.convert_segmentation_to_counts(df.gene, df[!, k]) for k in [cell_col1, cell_col2]]

    jac_dists = Dict{Tuple{Int, Int}, Float64}();
    exp_cors = Dict{Tuple{Int, Int}, Float64}();
    dists = Float64[]
    for i in 1:size(df, 1)
        c1, c2 = df[i, cell_col1], df[i, cell_col2]
        if (c1 == 0) || (c2 == 0)
            push!(dists, 0.0)
            continue
        end

        mpc1, mpc2 = mol_ids_per_cell[1][c1], mol_ids_per_cell[2][c2]
        if (length(mpc1) < min_mols_per_cell) || (length(mpc2) < min_mols_per_cell)
            push!(dists, 0.0)
            continue
        end

        if !((c1, c2) in keys(jac_dists))
            jac_dists[(c1, c2)] = length(intersect(mpc1, mpc2)) / length(union(mpc1, mpc2))

            exp_cors[(c1, c2)] = cor(cms[1][:,c1], cms[2][:,c2])
        end
        push!(dists, jac_dists[(c1, c2)])
    end

    return (dists=dists, jac_dists=jac_dists, exp_cors=exp_cors, cms=cms)
end

function precision_recall_df(contingency::SparseArrays.SparseMatrixCSC{Int64, Int64})
    overlap_sizes, overlap_ids = findmax(contingency[2:end, 2:end], dims=1);
    overlap_sizes = collect(overlap_sizes)[:];
    overlap_ids = getindex.(overlap_ids, 1)[:];

    n_mols_per_cell_src = sum(contingency[2:end,:], dims=2)[:];
    n_mols_per_cell_targ = sum(contingency[:, 2:end], dims=1)[:];

    n_overlaps = sum(contingency[2:end, 2:end] ./ n_mols_per_cell_targ' .> 0.25, dims=1)[:]
    recall = overlap_sizes ./ n_mols_per_cell_targ;
    precision = overlap_sizes ./ n_mols_per_cell_src[overlap_ids];
    return DataFrame(Dict(:precision => precision, :recall => recall, :n_overlaps => n_overlaps, :overlap_sizes => overlap_sizes))
end

max_stars_repo_name: hiraksarkar/BaysorAnalysis

max_issues_repo_name: hiraksarkar/BaysorAnalysis

max_forks_repo_name: hiraksarkar/BaysorAnalysis


text:
using BenchmarkTools
using Printf
using Transducers

benchmark = false

exampleFilename = "src/advent/day13/example.txt"
example1Filename = "src/advent/day13/example1.txt"
example2Filename = "src/advent/day13/example2.txt"
example3Filename = "src/advent/day13/example3.txt"
example4Filename = "src/advent/day13/example4.txt"
example5Filename = "src/advent/day13/example5.txt"
inputFilename = "src/advent/day13/input.txt"

struct Notes
    t::Int64
    buses::Array{Union{Missing, Int64}}
end

function readData(filename)
    f = open(filename)
    t = parse(Int64, readline(f))
    buses = split(readline(f), ",") |>
            Map(s -> s == "x" ? Missing() : parse(Int64, s)) |>
            collect
    close(f)
    Notes(t, buses)
end

exampleNotes = readData(exampleFilename)
example1Notes = readData(example1Filename)
example2Notes = readData(example2Filename)
example3Notes = readData(example3Filename)
example4Notes = readData(example4Filename)
example5Notes = readData(example5Filename)
inputNotes = readData(inputFilename)

#### part1

function part1(notes)
    result = notes.buses |>
            Filter(bus -> !ismissing(bus)) |>
            Map(bus -> (bus = bus, t = notes.t + (bus - (notes.t % bus)) % bus)) |>
            foldxl((a, b) -> a.t < b.t ? a : b)
    result.bus * (result.t - notes.t)
end

@assert part1(exampleNotes) == 295

@show part1(inputNotes)
benchmark && @btime part1(inputNotes);

#### part2

function part2(notes)
    buses = NamedTuple{(:number, :index), Tuple{Int64, Int64}}[]
    for i = 1:length(notes.buses)
        if !ismissing(notes.buses[i])
            push!(buses, (number = notes.buses[i], index = i))
        end
    end
    sort!(buses; rev = true, by = bus -> bus.number)
    referenceBus = buses[1]
    otherBuses = @view buses[2:length(buses)]

    t = referenceBus.number
    while true
        matched = true
        for bus in otherBuses
            t2 = t + bus.index - referenceBus.index
            if t2 % bus.number != 0
                matched = false
                break
            end
        end
        matched && break
        t += referenceBus.number
    end
    t - (referenceBus.index - 1)
end

@assert part2(exampleNotes) == 1068781
@assert part2(example1Notes) == 3417
@assert part2(example2Notes) == 754018
@assert part2(example3Notes) == 779210
@assert part2(example4Notes) == 1261476
@assert part2(example5Notes) == 1202161486

#@show part2(inputNotes)
#benchmark && @btime part2(inputNotes);


text:
using MethodAnalysis: methodinstances
using PrecompileMacro
using Test

@testset "PrecompileMacro" begin
    @precompile g(path::String, id::Int) = string(path, id)

    # Verify that one methodinstance is ready even though we didn't call the function yet.
    @test length(methodinstances(g)) == 1

    # To be sure that the method works.
    @test g("lorem", 1) == "lorem1"

    @precompile h(a::Int, b::Int; c=3) = string(a, b, c)
    @test length(methodinstances(h)) == 1

    @test_throws LoadError eval(:(@precompile f(path::AbstractString) = 1))

    "Some doc"
    @precompile k() = 3
    @test length(methodinstances(k)) == 1
    @test contains(string(@doc k), "Some doc")
end


text:
@discord_object struct MessageApplication
    id::Snowflake
    cover_image::String
    description::String
    icon::String
    name::String
end


text:
using Documenter, AllanDeviations

makedocs(
	format = :html,
    sitename = "AllanDeviations.jl",
    authors = "Julien Kluge",
    pages = [
        "Home" => "index.md",
        "Installation Guide" => "installation.md",
        "Quick Start Guide" => "quickstart.md",
        "Library" => Any[
			"API guide" => "lib/apiguide.md",
			"Allan Deviation" => "lib/allandev.md",
			"Modified Allan Deviation" => "lib/mallandev.md",
			"Hadamard Deviation" => "lib/hadamarddev.md",
			"Time Deviation" => "lib/timedev.md",
			"Total Deviation" => "lib/totaldev.md",
			"Max. time interval error" => "lib/mtie.md"
		]
    ]
)

deploydocs(
    repo   = "github.com/JulienKluge/AllanDeviations.jl.git",
    target = "build",
    osname = "linux",
    julia  = "nightly",
    deps = nothing,
    make = nothing
)