text:
__precompile__(true)

module QuNet

# TODO: Update what we should be using vs importing. Be discriminating!
using LightGraphs, SimpleWeightedGraphs, GraphPlot, MetaGraphs, GraphRecipes
using LinearAlgebra, StatsBase, Statistics
using Documenter, Colors, Plots, LaTeXStrings

# using SatelliteToolbox
# using QuadGK

import SparseArrays: dropzeros!
import Base: *, print, string
import GraphPlot: gplot
import QuadGK: quadgk
import SatelliteToolbox: expatmosphere

abstract type QObject end
abstract type QNode <: QObject end
abstract type QChannel <: QObject end

TIME_STEP = 0.01

# WARNING The order of these is fairly important.
# Don't change them willy-nilly unless you like segfaults and screaming
include("Network.jl")
include("TemporalGraphs.jl")
include("CostVector.jl")
include("Node.jl")
include("Channel.jl")
include("Percolation.jl")
include("Routing.jl")
include("Plot.jl")
include("Utilities.jl")
include("Benchmarking.jl")
include("GraphInterface.jl")
include("TypeTree.jl")

export
# Abstract Classes
QObject, QNode, QChannel,

# Benchmarking.jl
percolation_bench, dict_average, dict_err, make_user_pairs, net_performance,

# Channel.jl
BasicChannel, AirChannel,

# CostVector.jl
zero_costvector, unit_costvector, convert_costs, get_pathcv,

# Network.jl
QNetwork, GridNetwork, add, update,
refresh_graph!,

# Node.jl
Coords, BasicNode, Velocity, PlanSatNode,

# Percolation.jl

# Plot.jl
gplot,

# GraphInterface.jl
hard_rem_edge!,

# Routing.jl
shortest_path,

# TemporalGraphs.jl

# Utilities.jl
dB_to_P, P_to_dB, Z_to_dB, dB_to_Z, purify

# TypeTree.jl
qunet_type_tree

end

max_stars_repo_name: drpeterrohde/QuNet

max_issues_repo_name: SamuelMarks/QuNet

max_forks_repo_name: SamuelMarks/QuNet


text:
abstract type Preconditioner end

struct IdentityPreconditioner <: Preconditioner
end

struct JacobiPreconditioner{T <: Number} <: Preconditioner
    diagonal_list::Array{Array{T, 1}, 1}
end

function apply!(precond::IdentityPreconditioner, baselines) where {T <: Number}
    # Do nothing
end

function apply!(precond::JacobiPreconditioner{T}, baselines) where {T <: Number}
    @assert length(baselines) == length(precond.diagonal_list)

    for idx in eachindex(baselines)
        baselines[idx] .*= precond.diagonal_list[idx]
    end
end

function jacobi_preconditioner(baseline_length_list, T)
    # TODO: once we specify a mask or flags, we need to take them into account here!
    diagonals = [Array{T}(undef, length(x)) for x in baseline_length_list]

    for cur_diag_idx in eachindex(diagonals)
        diagonals[cur_diag_idx] .= 1 ./ baseline_length_list[cur_diag_idx]
    end

    JacobiPreconditioner{T}(diagonals)
end

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ a4eba47e-3cde-11eb-2a1c-7f5102c9a9e1
begin
	using HDF5
	using Plots
	using PlutoUI
	using LaTeXStrings
	using Printf
	import Pda
end

# ╔═╡ b0c18c94-3cde-11eb-0c12-731d453f428b
md"""## Exponential RTP in free space"""

# ╔═╡ f2663b88-3cdf-11eb-2869-5f3d5ff56948
md"""The MSD of exponential RTPs with complete reorientation and initial condition $$P(\mathbf{x}, \theta,t) = \delta(\mathbf{x})/2\pi$$ is

$$\left<x^2\right> =\left<y^2\right>=  U^2_0\tau^2_R\left( \frac{t}{\tau_R} + \exp(-t/\tau_R) -1\right)$$

Ref: [Run-and-tumble bacteria slowly approaching the diffusive regime](https://doi.org/10.1103/PhysRevE.101.062607), Villa-Torrealba *et. al*, PRE, 2020
"""

# ╔═╡ 3635fca6-3ce0-11eb-0d4b-1310b1716b86
function msd_exponential(t::Real, U0::Real, tauR::Real)
	return U0^2*tauR^2*(t/tauR+exp(-t/tauR) -1)
end

# ╔═╡ 8b42ee14-3ce2-11eb-0702-5da8b9887bd9
function get_var(file::String)
	t = Pda.get_h5data(file, "x/t")[2:end]
	v = Pda.get_h5data(file, "x/v")[2:end]
	tauR = Pda.get_h5attr(file, "tauR")
	U0 = Pda.get_h5attr(file, "U0")
	ell = U0*tauR
	# non-dimensionalized by tauR and ell
	t ./= tauR
	v ./= ell^2
	return t,v
end

# ╔═╡ d0f3bf32-3cde-11eb-13c0-6d1d3bd0b070
begin
	file = "U1.000tauR1.200.h5"
	t,v = get_var(file)
	plot(t,v,
		framestyle=:box,
		label="Langevin",
		lw=2)
	plot!(t, msd_exponential.(t, 1, 1)/1^2,
		lw=2,
		linestyle=:dash,
		label="theory",
		legend=:topleft,
		xscale=:log10,
		yscale=:log10)
	xlabel!(L"t/\tau_R")
	ylabel!(L"\mathrm{var}(x)/\ell^2")
	title!("Exponential RTP in free space")
end

# ╔═╡ 48f91100-3ce2-11eb-1331-930d15c04e16
md"""## Pareto RTP in free space"""

# ╔═╡ 53054a04-3ce2-11eb-2869-973ecb35edeb
begin
	file1= "U1.000alpha1.200free.h5"
	t1,v1 = get_var(file1)
	plot(t1,v1,
		framestyle=:box,
		label="Langevin",
		lw=2)
	alpha = Pda.get_h5attr(file1, "alpha")
	# #theory is <x^2> \sim t^(3-\alpha)
	i = 3000
	plot!(t1[end-i:end], (t1[end-i:end]/t1[end]).^(3-alpha)*v1[end],
		lw=2,
		linestyle=:dash,
		label="theory",
		legend=:topleft,
		xscale=:log10,
		yscale=:log10)
	xlabel!(L"t/\tau_R")
	ylabel!(L"\mathrm{var}(x)/\ell^2")
	title!("Pareto RTP in free space")
end

# ╔═╡ Cell order:
# ╠═a4eba47e-3cde-11eb-2a1c-7f5102c9a9e1
# ╟─b0c18c94-3cde-11eb-0c12-731d453f428b
# ╟─f2663b88-3cdf-11eb-2869-5f3d5ff56948
# ╠═3635fca6-3ce0-11eb-0d4b-1310b1716b86
# ╠═8b42ee14-3ce2-11eb-0702-5da8b9887bd9
# ╠═d0f3bf32-3cde-11eb-13c0-6d1d3bd0b070
# ╟─48f91100-3ce2-11eb-1331-930d15c04e16
# ╠═53054a04-3ce2-11eb-2869-973ecb35edeb


text:
macro wandbconfig(exs...)
    kwargs = Any[]
    for ex in exs
        if ex isa Expr && ex.head === :(=)
            k, v = ex.args
            if !(k isa Symbol)
                k = Symbol(k)
            end
            push!(kwargs, Expr(:kw, k, esc(v)))
        elseif ex isa Expr && ex.head === :... # Keyword splatting
            push!(kwargs, esc(ex))
        else
            k, v = ex, ex
            if !(k isa Symbol)
                k = Symbol(k)
            end
            push!(kwargs, Expr(:kw, k, esc(v)))
        end
    end
    return quote
        if !("JULIA_NO_WANDB" in keys(ENV)) || ENV["JULIA_NO_WANDB"] == "" || ENV["JULIA_NO_WANDB"] == "false" || ENV["JULIA_NO_WANDB"] == false
            kw = Dict(pairs((;$(kwargs...))))
            wandb.config.update(kw, allow_val_change=true)
        end
    end
end


text:
default_partitioner(::Type{<:SparseMatrixVBC{U, W, Tv, Ti}}) where {U, W, Tv, Ti} =
    AlternatePacker(
        EquiChunker(),
        EquiChunker(),
        DynamicTotalChunker(ConstrainedCost(model_SparseMatrixVBC_memory(Tv, Ti), WidthCost(), W)),
        DynamicTotalChunker(ConstrainedCost(permutedims(model_SparseMatrixVBC_memory(Tv, Ti)), WidthCost(), U)),
        DynamicTotalChunker(ConstrainedCost(model_SparseMatrixVBC_memory(Tv, Ti), WidthCost(), W)),
    )

function SparseMatrixVBC{U, W}(A::SparseMatrixCSC{Tv, Ti}, method=default_partitioner(SparseMatrixVBC{U, W, Tv, Ti})) where {U, W, Tv, Ti}
    Π, Φ = pack_plaid(A, method)
    return SparseMatrixVBC{U, W}(A, convert(SplitPartition, Π), convert(SplitPartition, Φ))
end

function SparseMatrixVBC{U, W}(A::SparseMatrixCSC{Tv, Ti}, Π::SplitPartition{Ti}, Φ::SplitPartition{Ti}) where {U, W, Tv, Ti}
    @inbounds begin
        (m, n) = size(A)
        K = length(Π)
        L = length(Φ)
        Π_spl = Π.spl
        Φ_spl = Φ.spl
        Π_asg = convert(MapPartition, Π).asg
        hst = zeros(Ti, K)
        A_pos = A.colptr
        A_idx = A.rowval
        A_val = A.nzval
        pos = undefs(Ti, L + 1)
        ofs = undefs(Ti, L + 1)
        pos[1] = 1
        ofs[1] = 1
        for l = 1:L
            ofs[l + 1] = ofs[l]
            pos[l + 1] = pos[l]
            j = Φ_spl[l]
            j′ = Φ_spl[l + 1]
            w = j′ - j
            for q = A_pos[j]:A_pos[j′] - 1
                i = A_idx[q]
                k = Π_asg[i]
                if hst[k] < l
                    u = Π_spl[k + 1] - Π_spl[k]
                    pos[l + 1] += 1
                    ofs[l + 1] += u * w
                end
                hst[k] = l
            end
        end

        idx = Vector{Ti}(undef, pos[end] - 1)
        Δw = fld(CpuId.simdbytes(), sizeof(Tv)) #TODO need to define an upper limit on this one
        val = Vector{Tv}(undef, ofs[end] - 1 + U * (Δw * cld(W, Δw)))
        for Q = ofs[end] : ofs[end]  - 1 + U * (Δw * cld(W, Δw)) #extra stuff to keep vector access in bounds
            val[Q] = zero(Tv)
        end

        A_q = ones(Int, W)

        for k = 1:K
            @assert Π_spl[k + 1] - Π_spl[k] <= U
        end

        for l = 1:L
            j = Φ_spl[l]
            w = Φ_spl[l + 1] - j
            @assert w <= W
            if w == 1
                Q = pos[l]
                q = ofs[l]
                A_q_1 = A_pos[j]
                while A_q_1 < A_pos[j + 1]
                    k = Π_asg[A_idx[A_q_1]]
                    for i = Π_spl[k] : Π_spl[k + 1] - 1
                        if A_q_1 < A_pos[j + 1] && A_idx[A_q_1] == i
                            val[q] = A_val[A_q_1]
                            A_q_1 += 1
                        else
                            val[q] = zero(Tv)
                        end
                        q += 1
                    end
                    idx[Q] = k
                    Q += 1
                end
            else
                k = K + 1
                for Δj = 1:w
                    A_q[Δj] = A_pos[j + Δj - 1]
                    if A_q[Δj] < A_pos[j + Δj]
                        k = min(k, Π_asg[A_idx[A_q[Δj]]])
                    end
                end
                Q = pos[l]
                q = ofs[l]
                while k != K + 1
                    for i = Π_spl[k] : Π_spl[k + 1] - 2
                        for Δj = 1:w
                            if A_q[Δj] < A_pos[j + Δj] && A_idx[A_q[Δj]] == i
                                val[q] = A_val[A_q[Δj]]
                                A_q[Δj] += 1
                            else
                                val[q] = zero(Tv)
                            end
                            q += 1
                        end
                    end
                    k′ = K + 1
                    let i = Π_spl[k + 1] - 1
                        for Δj = 1:w
                            if A_q[Δj] < A_pos[j + Δj]
                                if A_idx[A_q[Δj]] == i
                                    val[q] = A_val[A_q[Δj]]
                                    A_q[Δj] += 1
                                else
                                    val[q] = zero(Tv)
                                end
                                if A_q[Δj] < A_pos[j + Δj]
                                    k′ = min(k′, Π_asg[A_idx[A_q[Δj]]])
                                end
                            else
                                val[q] = zero(Tv)
                            end
                            q += 1
                        end
                    end
                    idx[Q] = k
                    Q += 1
                    k = k′
                end
            end
        end
        return SparseMatrixVBC{U, W, Tv, Ti}(m, n, Π, Φ, pos, idx, ofs, val)
    end
end


text:
#!/usr/bin/env julia

using DataFrames

import Celeste.AccuracyBenchmark
import Celeste.ArgumentParse

const OUTPUT_DIRECTORY = joinpath(splitdir(Base.source_path())[1], "output")
const CELESTE_PRIOR_CATALOG_CSV = joinpath(OUTPUT_DIRECTORY, "celeste_prior_catalog.csv")
const COADD_CATALOG_FITS = joinpath(Pkg.dir("Celeste"), "test", "data", "coadd_for_4263_5_119.fit")

GALAXY_ONLY_COLUMNS = [
    :de_vaucouleurs_mixture_weight,
    :minor_major_axis_ratio,
    :half_light_radius_px,
    :angle_deg,
]

parser = Celeste.ArgumentParse.ArgumentParser()
ArgumentParse.add_argument(
    parser,
    "action",
    help=string(
        "coadd: convert a Stripe82 coadd catalog to CSV form; ",
        "prior: draw sources at random from the Celeste prior",
    ),
)
ArgumentParse.add_argument(
    parser,
    "coadd_fits",
    help="Stripe82 catalog FITS file, for 'coadd'",
    default=COADD_CATALOG_FITS,
    required=false,
)
parsed_args = ArgumentParse.parse_args(parser, ARGS)

function main()
    if parsed_args["action"] == "coadd"
        catalog = AccuracyBenchmark.load_coadd_catalog(parsed_args["coadd_fits"])
        catalog_label = splitext(basename(parsed_args["coadd_fits"]))[1]
        output_filename = joinpath(OUTPUT_DIRECTORY, "$(catalog_label).csv")
    elseif parsed_args["action"] == "prior"
        catalog = AccuracyBenchmark.generate_catalog_from_celeste_prior(500, 12345)
        output_filename = CELESTE_PRIOR_CATALOG_CSV
    else
        @printf("Invalid action: '%s'\n\n", action)
        print_usage()
        return
    end

    # for stars, ensure galaxy-only fields are NA
    for column_symbol in GALAXY_ONLY_COLUMNS
        catalog[catalog[:is_star], column_symbol] = NA
    end

    if !isdir(OUTPUT_DIRECTORY)
        mkdir(OUTPUT_DIRECTORY)
    end
    AccuracyBenchmark.write_catalog(output_filename, catalog)
    AccuracyBenchmark.append_hash_to_file(output_filename)
end

main()


text:
#!/bin/bash
export JULIA_NUM_THREADS=4
julia -e 'idx_day_to_use=1; include("examples/calc_ccfs.jl")';
julia -e 'idx_day_to_use=2; include("examples/calc_ccfs.jl")';
julia -e 'idx_day_to_use=3; include("examples/calc_ccfs.jl")';
julia -e 'idx_day_to_use=4; include("examples/calc_ccfs.jl")';
julia -e 'idx_day_to_use=5; include("examples/calc_ccfs.jl")';
julia -e 'idx_day_to_use=6; include("examples/calc_ccfs.jl")';
julia -e 'idx_day_to_use=7; include("examples/calc_ccfs.jl")';
julia -e 'idx_day_to_use=8; include("examples/calc_ccfs.jl")';
julia -e 'idx_day_to_use=9; include("examples/calc_ccfs.jl")';
julia -e 'idx_day_to_use=10; include("examples/calc_ccfs.jl")';
julia -e 'idx_day_to_use=11; include("examples/calc_ccfs.jl")';
julia -e 'idx_day_to_use=12; include("examples/calc_ccfs.jl")';
#julia -e 'idx_day_to_use=13; include("examples/calc_ccfs.jl")';
#julia -e 'idx_day_to_use=14; include("examples/calc_ccfs.jl")';


text:
using Bridge
using Bridge.Models: ℝ
import Bridge: b, σ, B, β, a, constdiff

struct SinDiffusion{T} <: ContinuousTimeProcess{ℝ{1, T}}
    a::T
    b::T
    c::T
    σ::T
end

b(t, x, P::SinDiffusion) = ℝ{1}(P.a + P.b*sin.(P.c * x))
σ(t, x, P::SinDiffusion) = ℝ{1}(P.σ)
constdiff(::SinDiffusion) = true
clone(::SinDiffusion, θ) = SinDiffusion(θ...)
params(P::SinDiffusion) = [P.a, P.b, P.c, P.σ]

struct SinDiffusionAux{T,S} <: ContinuousTimeProcess{ℝ{1, T}}
    a::T
    b::T
    c::T
    σ::T
    t::Float64
    u::S
    T::Float64
    v::S
end

#B(t, P::SinDiffusionAux) = @SMatrix [0.0]
#β(t, P::SinDiffusionAux) = ℝ{1}(0.0)

#B(t, P::SinDiffusionAux) = @SMatrix [(P.b*P.c*cos.(P.c * P.v)) * t/P.T + (P.b*P.c*cos.(P.c * P.u)) * (1-t/P.T)]
#β(t, P::SinDiffusionAux) = ℝ{1}((P.a + P.b*sin.(P.c * P.v) - P.b*P.c*cos.(P.c * P.v) * P.v) * t/P.T + (P.a + P.b*sin.(P.c * P.u) - P.b*P.c*cos.(P.c * P.u) * P.u) * (1-t/P.T))

#B(t, P::SinDiffusionAux) = @SMatrix [P.b*P.c*cos.(P.c * P.v)]
#β(t, P::SinDiffusionAux) = ℝ{1}(P.a + P.b*sin.(P.c * P.v) - P.b*P.c*cos.(P.c * P.v) * P.v)

B(t, P::SinDiffusionAux) = @SMatrix [0 + (t/P.T)/5]#+ (t/P.T)/5
β(t, P::SinDiffusionAux) = ℝ{1}((P.v-P.u)/P.T)

σ(t, P::SinDiffusionAux) = ℝ{1}(P.σ)
dependsOnParams(::SinDiffusionAux) = (4,)
constdiff(::SinDiffusionAux) = true
b(t, x, P::SinDiffusionAux) = B(t, P) * x + β(t, P)
a(t, P::SinDiffusionAux) = σ(t, P) * σ(t, P)'
clone(P::SinDiffusionAux, θ) = SinDiffusionAux(θ..., P.t, P.u, P.T, P.v)


text:
module TestParameters

using Test
using LightGBM

NSAMPLES_TRAIN = 2000
NSAMPLES_TEST = 500
NFEATURES = 20

X_train = randn(NSAMPLES_TRAIN, NFEATURES)
y_train_binary = rand(0:1, NSAMPLES_TRAIN)
y_train_regression = rand(NSAMPLES_TRAIN)
X_test = randn(NSAMPLES_TEST, NFEATURES)

norm(x) = sqrt(sum(x .^ 2))

@testset "parameters -- boosting" begin

    # The scheme is to set up an old style gbdt and a new DART and show that the learned
    # models are "diffferent", i.e. the parameter passed through and made a difference
    # Check for both classification and regression
    # Why do we use num_iterations 12? well, GOSS doesn't "work" until after 1/learning_rate
    # and default learning rate is 0.1, so we need to make sure models are going for at least
    # enough time to be different.

    classifier_gdbt = LightGBM.LGBMClassification(boosting = "gbdt", num_iterations=12)
    classifier_dart = LightGBM.LGBMClassification(boosting = "dart", num_iterations=12)
    classifier_goss = LightGBM.LGBMClassification(boosting = "goss", num_iterations=12)
    regressor_gbdt = LightGBM.LGBMRegression(boosting = "gbdt", num_iterations=12)
    regressor_dart = LightGBM.LGBMRegression(boosting = "dart", num_iterations=12)
    regressor_goss = LightGBM.LGBMRegression(boosting = "goss", num_iterations=12)

    LightGBM.fit!(classifier_dart, X_train, y_train_binary, verbosity = -1)
    LightGBM.fit!(classifier_gdbt, X_train, y_train_binary, verbosity = -1)
    LightGBM.fit!(classifier_goss, X_train, y_train_binary, verbosity = -1)
    LightGBM.fit!(regressor_gbdt, X_train, y_train_regression, verbosity = -1)
    LightGBM.fit!(regressor_dart, X_train, y_train_regression, verbosity = -1)
    LightGBM.fit!(regressor_goss, X_train, y_train_regression, verbosity = -1)

    p_gbdt = LightGBM.predict(classifier_gdbt, X_test, verbosity = -1)
    p_dart = LightGBM.predict(classifier_dart, X_test, verbosity = -1)
    p_goss = LightGBM.predict(classifier_goss, X_test, verbosity = -1)
    r_gbdt = LightGBM.predict(regressor_gbdt, X_test, verbosity = -1)
    r_dart = LightGBM.predict(regressor_dart, X_test, verbosity = -1)
    r_goss = LightGBM.predict(regressor_goss, X_test, verbosity = -1)

    # checking the models are sufficiently different by making sure
    # that the distance between the vectors is relatively large
    @test norm(p_dart .- p_gbdt) >= sqrt(NSAMPLES_TEST) * 1e-5
    @test norm(p_goss .- p_gbdt) >= sqrt(NSAMPLES_TEST) * 1e-5
    @test norm(p_dart .- p_goss) >= sqrt(NSAMPLES_TEST) * 1e-5
    # Check also for regression
    @test norm(r_dart .- r_gbdt) >= sqrt(NSAMPLES_TEST) * 1e-5
    @test norm(r_goss .- r_gbdt) >= sqrt(NSAMPLES_TEST) * 1e-5
    @test norm(r_dart .- r_goss) >= sqrt(NSAMPLES_TEST) * 1e-5

end

@testset "parameters -- dart classifcation" begin

    # the scheme is to set up a DART classifier and show that changing each one of
    # the parameters individually results in a difference to the model

    dart_default = LightGBM.LGBMClassification(boosting = "dart")
    dart_drop_rate = LightGBM.LGBMClassification(boosting = "dart", drop_rate = 0.9)
    dart_max_drop = LightGBM.LGBMClassification(boosting = "dart", max_drop = 1)
    dart_skip_drop = LightGBM.LGBMClassification(boosting = "dart", skip_drop = 1.0)
    dart_xdm = LightGBM.LGBMClassification(boosting = "dart", xgboost_dart_mode = true)
    dart_uniform_drop = LightGBM.LGBMClassification(boosting = "dart", uniform_drop = true)
    dart_drop_seed = LightGBM.LGBMClassification(boosting = "dart", drop_seed = 20)

    LightGBM.fit!(dart_default, X_train, y_train_binary, verbosity = -1)
    LightGBM.fit!(dart_drop_rate, X_train, y_train_binary, verbosity = -1)
    LightGBM.fit!(dart_max_drop, X_train, y_train_binary, verbosity = -1)
    LightGBM.fit!(dart_skip_drop, X_train, y_train_binary, verbosity = -1)
    LightGBM.fit!(dart_xdm, X_train, y_train_binary, verbosity = -1)
    LightGBM.fit!(dart_uniform_drop, X_train, y_train_binary, verbosity = -1)
    LightGBM.fit!(dart_drop_seed, X_train, y_train_binary, verbosity = -1)

    p_dart_default = LightGBM.predict(dart_default, X_test, verbosity = -1)
    p_dart_drop_rate = LightGBM.predict(dart_drop_rate, X_test, verbosity = -1)
    p_dart_max_drop = LightGBM.predict(dart_max_drop, X_test, verbosity = -1)
    p_dart_skip_drop = LightGBM.predict(dart_skip_drop, X_test, verbosity = -1)
    p_dart_xdm = LightGBM.predict(dart_xdm, X_test, verbosity = -1)
    p_dart_uniform_drop = LightGBM.predict(dart_uniform_drop, X_test, verbosity = -1)
    p_dart_drop_seed = LightGBM.predict(dart_drop_seed, X_test, verbosity = -1)

    @test norm(p_dart_default .- p_dart_drop_rate) >= sqrt(NSAMPLES_TEST) * 1e-5
    @test norm(p_dart_default .- p_dart_max_drop) >= sqrt(NSAMPLES_TEST) * 1e-10
    @test norm(p_dart_default .- p_dart_skip_drop) >= sqrt(NSAMPLES_TEST) * 1e-5
    @test norm(p_dart_default .- p_dart_xdm) >= sqrt(NSAMPLES_TEST) * 1e-10
    @test_broken norm(p_dart_default .- p_dart_uniform_drop) >= sqrt(NSAMPLES_TEST) * 1e-5
    @test norm(p_dart_default .- p_dart_drop_seed) >= sqrt(NSAMPLES_TEST) * 1e-5

end

@testset "parameters -- dart regression" begin

    # the scheme is to set up a DART regressor and show that changing each one of
    # the parameters individually results in a difference to the model

    dart_default = LightGBM.LGBMRegression(boosting = "dart")
    dart_drop_rate = LightGBM.LGBMRegression(boosting = "dart", drop_rate = 0.9)
    dart_max_drop = LightGBM.LGBMRegression(boosting = "dart", max_drop = 1)
    dart_skip_drop = LightGBM.LGBMRegression(boosting = "dart", skip_drop = 1.0)
    dart_xdm = LightGBM.LGBMRegression(boosting = "dart", xgboost_dart_mode = true)
    dart_uniform_drop = LightGBM.LGBMRegression(boosting = "dart", uniform_drop = true)
    dart_drop_seed = LightGBM.LGBMRegression(boosting = "dart", drop_seed = 20)

    LightGBM.fit!(dart_default, X_train, y_train_regression, verbosity = -1)
    LightGBM.fit!(dart_drop_rate, X_train, y_train_regression, verbosity = -1)
    LightGBM.fit!(dart_max_drop, X_train, y_train_regression, verbosity = -1)
    LightGBM.fit!(dart_skip_drop, X_train, y_train_regression, verbosity = -1)
    LightGBM.fit!(dart_xdm, X_train, y_train_regression, verbosity = -1)
    LightGBM.fit!(dart_uniform_drop, X_train, y_train_regression, verbosity = -1)
    LightGBM.fit!(dart_drop_seed, X_train, y_train_regression, verbosity = -1)

    r_dart_default = LightGBM.predict(dart_default, X_test, verbosity = -1)
    r_dart_drop_rate = LightGBM.predict(dart_drop_rate, X_test, verbosity = -1)
    r_dart_max_drop = LightGBM.predict(dart_max_drop, X_test, verbosity = -1)
    r_dart_skip_drop = LightGBM.predict(dart_skip_drop, X_test, verbosity = -1)
    r_dart_xdm = LightGBM.predict(dart_xdm, X_test, verbosity = -1)
    r_dart_uniform_drop = LightGBM.predict(dart_uniform_drop, X_test, verbosity = -1)
    r_dart_drop_seed = LightGBM.predict(dart_drop_seed, X_test, verbosity = -1)

    @test norm(r_dart_default .- r_dart_drop_rate) >= sqrt(NSAMPLES_TEST) * 1e-5
    @test norm(r_dart_default .- r_dart_max_drop) >= sqrt(NSAMPLES_TEST) * 1e-10
    @test norm(r_dart_default .- r_dart_skip_drop) >= sqrt(NSAMPLES_TEST) * 1e-5
    @test norm(r_dart_default .- r_dart_xdm) >= sqrt(NSAMPLES_TEST) * 1e-10
    @test_broken norm(r_dart_default .- r_dart_uniform_drop) >= sqrt(NSAMPLES_TEST) * 1e-5
    @test norm(r_dart_default .- r_dart_drop_seed) >= sqrt(NSAMPLES_TEST) * 1e-5

end

@testset "parameters -- goss classifcation" begin

    # the scheme is to set up a goss classifier and show that changing each one of
    # the parameters individually results in a difference to the model

    goss_default = LightGBM.LGBMClassification(boosting = "goss", num_iterations=12)
    goss_top_rate = LightGBM.LGBMClassification(boosting = "goss", num_iterations=12, top_rate = 0.01)
    goss_other_rate = LightGBM.LGBMClassification(boosting = "goss", num_iterations=12, other_rate = 0.01)

    LightGBM.fit!(goss_default, X_train, y_train_binary, verbosity = -1)
    LightGBM.fit!(goss_top_rate, X_train, y_train_binary, verbosity = -1)
    LightGBM.fit!(goss_other_rate, X_train, y_train_binary, verbosity = -1)

    p_goss_default = LightGBM.predict(goss_default, X_test, verbosity = -1)
    p_goss_top_rate = LightGBM.predict(goss_top_rate, X_test, verbosity = -1)
    p_goss_other_rate = LightGBM.predict(goss_other_rate, X_test, verbosity = -1)

    @test norm(p_goss_default .- p_goss_top_rate) >= sqrt(NSAMPLES_TEST) * 1e-5
    @test norm(p_goss_default .- p_goss_other_rate) >= sqrt(NSAMPLES_TEST) * 1e-5

end

@testset "parameters -- goss regression" begin

    # the scheme is to set up a goss regressor and show that changing each one of
    # the parameters individually results in a difference to the model

    goss_default = LightGBM.LGBMRegression(boosting = "goss", num_iterations=12)
    goss_top_rate = LightGBM.LGBMRegression(boosting = "goss", num_iterations=12, top_rate = 0.01)
    goss_other_rate = LightGBM.LGBMRegression(boosting = "goss", num_iterations=12, other_rate = 0.01)

    LightGBM.fit!(goss_default, X_train, y_train_regression, verbosity = -1)
    LightGBM.fit!(goss_top_rate, X_train, y_train_regression, verbosity = -1)
    LightGBM.fit!(goss_other_rate, X_train, y_train_regression, verbosity = -1)

    r_goss_default = LightGBM.predict(goss_default, X_test, verbosity = -1)
    r_goss_top_rate = LightGBM.predict(goss_top_rate, X_test, verbosity = -1)
    r_goss_other_rate = LightGBM.predict(goss_other_rate, X_test, verbosity = -1)

    @test norm(r_goss_default .- r_goss_top_rate) >= sqrt(NSAMPLES_TEST) * 1e-5
    @test norm(r_goss_default .- r_goss_other_rate) >= sqrt(NSAMPLES_TEST) * 1e-5

end

end # end module


text:
struct Event
	id_node::Int
	position_node::Int
	position_type::Int
	time::DateTime
end

struct EventWriter
	file::IOStream
	queue::Queue{Event}

	function EventWriter(filename::String)
		file = open(filename,"w+")
		queue = Queue{Event}()
		_write_header(file)
		new(file,queue)
	end
end

function _write_header(file)
	write(file,"Time,Id,position\n")
end

function add_event(event_writer::EventWriter,event::Event)
	#FIXME: Just a way to handle null event. Must be substituted by a better structure
	if event.id_node == -1
		return
	end
	enqueue!(event_writer.queue,event)
end

function write_events(event_writer::EventWriter)
	for event in event_writer.queue
		write(event_writer.file,"$(event.time),$(event.id_node),$(event.position_node)\n")
		dequeue!(event_writer.queue)
	end
end

function close_writer(event_writer::EventWriter)
	write_events(event_writer)
	close(event_writer.file)
end

max_stars_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_issues_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_forks_repo_name: Luigi-Crisci/Hypergraphs-Epidemics


text:
type SGsd <: StepData
	getStep::Function
	stepString::AbstractString
	function SGsd(stepString::AbstractString)
		new(computeSGStep,stepString)
	end
	function SGsd()
		new(computeSGStep,"SG")
	end
end

function computeSGStep(x, k, gnum, sd::SGsd, problem::Problem; outputLevel =0)

	func = consume(problem.getNextSampleFunction)

	(f, g) = func(x)

	gnum += 1

	(g, gnum)

end


text:
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 5.
#

module blocksys

include("types.jl")
import .types.FL, .types.SquareSparseMatrix

using SparseArrays

# constants
macheps = eps(FL)

"""
### Asserts provided arguments match in size.

### Arguments:
- `A` — the `A` matrix
- `n` — the size of the matrix
- `b` — the `b` vector

### Returns: `Nothing`

"""
function assertMatrixSize(A, n::Int, b::Vector=[0 for _ in 1:n])::Nothing

    # if size(A, 1) == size(A, 2) == n == length(b)
    if n == length(b)
        return
    end

    throw(ErrorException("provided arguments do not match in sizes"))

end

"""
### Asserts `n` is divisible by `l`.

### Arguments:
- `A` — the `A` matrix
- `n` — the size of the matrix
- `b` — the `b` vector

### Returns: `Nothing`

"""
function assertBlockSize(n::Int, l::Int)::Nothing

    if n%l != 0
        throw(ErrorException("n is not divisible by l"))
    end

end

"""
### Calculates the permutation vector of given matrix.
This function implements the *partial pivoting* algorithm that assures that the diagonal of the matrix contains numbers that are not near zero or equal to zero.

### Arguments:
- `A` — the `A` matrix
- `n` — size of the matrix
- `l` — size of each chunk (submatrix)
- `limited` — this flag limits the bounds of the searching loop; it may be used when the function gives a permutation that leads to having a value near zero on the diagonal of the matrix; default value is set to `false`

### Returns
— the pivot vector that simulates all the matrix rows swaps that would have occurred

"""
function partialPivot(A, n::Int, l::Int, limited::Bool = false)::Vector

    assertMatrixSize(A, n)
    assertBlockSize(n, l)

    # save the new order of the rows
    pivot = [x for x in 1 : n]

    for k in 1:Int(n/l)

        for p in (k-1)*l+1 : k*l
            # get current value
            curr = abs(A[pivot[p], p])

            # check if it’s the biggest value
            max_value = curr
            max_index = pivot[p]
            search_end = limited ? k*l : k*l + l
            for i in (p+1) : search_end

                # special case for k=n/l
                # because this is the last submatrix there are no submatrices underneath;
                # thus,the internal loop used for searching max value needs to step earlier
                if i > n
                    break
                end
                # get another value from that column
                potential = A[pivot[i], p]
                if abs(potential) > max_value
                    # its absolute value is bigger — mark as ready for swap
                    max_value = abs(potential)
                    max_index = pivot[i]
                end
            end
            if max_value < macheps
                throw(ErrorException("a value near or equal zero was selected as the main value (diagonal value) — try to run `partialPivot` with `limited` set to `true` or do not use pivoting for this particular case at all"))
            end
            # swap the rows if needed
            if max_value > curr
                pivot[p], pivot[max_index] = pivot[max_index], pivot[p]
            end

        end

    end

    return pivot

end

"""
### Gives LU decomposition of given matrix.

### Arguments:
- `A` — the matrix to decompose
- `b` — the `b` vector (the vector that is on the right side of the `Ax = b` equation)
- `n` — the size of the matrix
- `l` — size of each chunk (submatrix)
- `pivot` — optional pivoting vector that simulates matrix row swap (can be generated using the `partialPivot` function); default value is set to a linear vector [1,…,n]

### Returns a tuple:
- `[1]` — matrices `L` and `U` contained in one matrix for optimization
- `[2]` — altered `b' = b^{(k)}` vector

"""
function decomposeIntoLU(A, b::Vector, n::Int, l::Int, pivot::Vector = [x for x in 1:n])#::Tuple{types.SquareSparseMatrix, Vector}

    assertMatrixSize(A, n)
    assertBlockSize(n, l)

    # this a more optimized way of storing both L and U matrices
    # we copy `A` at the beginning because that’s the start of the sequence:
    # A^{(1)} → A^{(2)} → … → A^{(n)} = U
    # in the end both matrices are triangular, so there is no problem with storing them this way
    # and the `L` matrix does not store any information on the diagonal of the matrix so there is no collision between `U` and `L`
    LU = copy(A)

    bprime = copy(b)

    # these two adjacent `for` loops are just one loop split into `l`-sized chunks
    # — its time complexity is still O(n)
    for k in 1:Int(n/l)
        for p in (k-1)*l+1 : k*l

            # update all rows beneath
            for i in (p+1) : k*l + l
                # special case k = n/l
                if i > n
                    break
                end

                lᵢₚ = LU[pivot[i], p] / LU[pivot[p], p]

                # process the row with new coefficients
                for j in p : k*l + l
                    if j > n
                        break
                    end
                    LU[pivot[i], j] = LU[pivot[i], j] - lᵢₚ * LU[pivot[p], j]
                end

                # adjust the `b` matrix accordingly
                bprime[pivot[i]] = bprime[pivot[i]] - lᵢₚ * bprime[pivot[p]]

                # save the `L` matrix part
                LU[pivot[i], p] = lᵢₚ
            end
        end
    end

    return LU, bprime

end

"""
### Solves the standard equation `Ux = y` where `L` is an upper-triangular matrix.

### Arguments:
- `U` — the upper-triangular matrix
- `y` — the right side of the equation
- `n` — the size of the matrix
- `l` — size of each chunk (submatrix)
- `pivot` — optional pivoting vector that simulates matrix row swap (can be generated using the `partialPivot` function); default value is set to a linear vector [1,…,n]

### Returns
— the `x` output vector (solution to the `Ux = y` equation)

"""
function solveUpperTriangularMatrix(U, y::Vector, n::Int, l::Int, pivot::Vector = [x for x in 1:n])::Array

    # init the output vector
    x = [0.0 for _ in 1:n]

    # calc the values of the output vector
    for k in Int(n/l) : -1 : 1
        for p in k*l : -1 : (k-1)*l+1

            tmp = y[pivot[p]]
            for i in (p+1) : (k*l + l)
                if i > n
                    break
                end
                tmp -= x[pivot[i]] * U[pivot[p], i]
            end
            x[pivot[p]] = tmp/U[pivot[p], p]

        end
    end

    return x

end

"""
### Solves the `LUx = b` equation.

### Arguments:
- `LU` — LU decomposition of a given matrix
- `b` — the `b` vector (the right side of the `LUx = b` equation)
- `n` — the size of the matrix
- `l` — size of each chunk (submatrix)
- `pivot` — optional pivoting vector that simulates matrix row swap (can be generated using the `partialPivot` function); default value is set to a linear vector [1,…,n]

### Returns
— the `x` output vector (solution to the `LUx = b` equation)

"""
function solveFromLU(LU, b::Vector, n::Int, l::Int, pivot::Vector = [x for x in 1:n])::Array

    # Ly = b

    y = zeros(n)

    # calc the values of the intermediary vector
    for k in 1:Int(n/l)
        for p in (k-1)*l+1 : k*l

            tmp = b[pivot[p]]
            for i in (p-1) : -1 : (k-1)*l - l
                if i <= 0
                    break
                end
                tmp -= y[pivot[i]] * LU[pivot[p], i]
            end

            y[pivot[p]] = tmp#/1

        end
    end

    # Ux = y
    return solveUpperTriangularMatrix(LU, y, n, l, pivot)

end

"""
### Solves the `Ax = b` equation.

### Arguments:
- `A` — the input matrix
- `b` — the `b` vector (the right side of the `LUx = b` equation)
- `n` — the size of the matrix
- `l` — size of each chunk (submatrix)
- `pivot` — optional pivoting vector that simulates matrix row swap (can be generated using the `partialPivot` function); default value is set to a linear vector [1,…,n]

### Returns
— the `x` output vector (solution to the `Ax = b` equation)

"""
function gaussianElimination(A, b::Vector, n::Int, l::Int, pivot::Vector = [x for x in 1:n])::Array

    # get the `U` matrix and `b'`
    # here we don’t actually need the `L` matrix, but it is generated either way — but that doesn’t affect performance
    LU, bprime = decomposeIntoLU(A, b, n, l, pivot)

    return solveUpperTriangularMatrix(LU, bprime, n, l, pivot)

end

"""
### Generates the `b` vector with respect to given `A` matrix.
*(Assuming `x` is a vector of ones.)*

### Arguments:
- `A` — the input matrix
- `n` — the size of the matrix
- `l` — size of each chunk (submatrix)

### Returns
— the `b` vector

"""
function generateB(A, n::Int, l::Int)::Vector

    b = zeros(n)

    for k in 1:Int(n/l)
        for p in (k-1)*l+1 : k*l

            for j in k : k*l + l
                if j > n
                    break
                end
                b[p] += A[p, j]
            end

        end
    end

    return b

end

end

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
"""
dyn_back is code for the traceback. The mxmy, mxiy, mxdy type code simply does the work of taking a HMM state for both
	x/y (MxMy) is "M"atch in X, "M"atch in Y, with "I"nsert, "D"elete other options, calculating the cost of transitioning
	from this state at this node to the next one and thereby either choosing the best scoring transition if we are beneath
	some threshold or otherwise making the choice stochastically through a softmax decision. These functions are what we
	call in the recursive traceback after our tensors have been constructed.
	The numbers 1-7 are used to index the states, which are MxMy, MxIy, MxDy, etc.
"""
module dyn_back

using StatsBase

#This is same as other softmax code but we'll have a local copy.
function back_softmax(vec, num_samples::Int=1, temperature::Float64=1.0)
	vec = vec ./ temperature
	max_value = maximum(vec)

	log_partition = log.(sum(exp.(vec .- max_value))) .+ max_value

	prob_vec = exp.(vec .- log_partition)
	wv = Weights(prob_vec)
	result = sample(1:length(prob_vec), wv, num_samples)
	return result
end

function mxmy(trinucs, last_letter::AbstractString, short_node::Int64, far_node::Int64, cumulative_score::Float64, cur_score::Float64, passed_values::Array{Float64, 1}, bbt, thresh, temperature)
	short_mm, short_md, short_im, short_dm, short_dd, short_mi, short_ii, far_mm, far_md, far_im, far_dm, far_dd, far_mi, far_ii = passed_values;
	mxmy_back = bbt[last_letter][short_node - 1, far_node - 1, 1:7, 1:64]

	mxmy_array = zeros(7 * 64)

	#Similar to other functions, we just look at the appropriate next "cell" in tensor and assess scores with the transition before making this "next stop"
	mxmy_array[1:64] = -vec(mxmy_back[1, 1:end]) .+ (short_mm + far_mm)
	mxmy_array[65:128] = -vec(mxmy_back[2, 1:end]) .+ (short_im + far_mm)
	mxmy_array[129:192] = -vec(mxmy_back[3, 1:end]) .+ (short_dm + far_mm)
	mxmy_array[193:256] = -vec(mxmy_back[4, 1:end]) .+ (short_mm + far_im)
	mxmy_array[257:320] = -vec(mxmy_back[5, 1:end]) .+ (short_dm + far_im)
	mxmy_array[321:384] = -vec(mxmy_back[6, 1:end]) .+ (short_mm + far_dm)
	mxmy_array[385:448] = -vec(mxmy_back[7, 1:end]) .+ (short_im + far_dm)

	if rand() < thresh
		next_stop = findmax(-mxmy_array)[2]
	else
		next_stop = back_softmax(-mxmy_array, 1, temperature)[1]
	end
	#some logic about handling this next stop...
	next_state = div(next_stop-1, 64) + 1
	next_letter = div(mod(next_stop-1, 64), 16) + 1
	next_codon = trinucs[mod(next_stop-1, 64) + 1]
	cumulative_score += mxmy_array[next_stop]
	cur_score = mxmy_array[next_stop]
	return next_state, next_letter, next_codon, cur_score, cumulative_score
end

function mxiy(trinucs, last_letter::AbstractString, short_node::Int64, far_node::Int64, cumulative_score::Float64, cur_score::Float64, passed_values::Array{Float64, 1}, bbt, thresh, temperature)
	short_mm, short_md, short_im, short_dm, short_dd, short_mi, short_ii, far_mm, far_md, far_im, far_dm, far_dd, far_mi, far_ii = passed_values;
	mxiy_back = bbt[last_letter][short_node - 1, far_node, 1:5, 1:64]

	mxiy_array = zeros(5 * 64)

	mxiy_array[1:64] = -vec(mxiy_back[1, 1:end]) .+ (short_mm + far_mm)
	mxiy_array[65:128] = -vec(mxiy_back[2, 1:end]) .+ (short_im + far_mm)
	mxiy_array[129:192] = -vec(mxiy_back[3, 1:end]) .+ (short_dm + far_mm)
	mxiy_array[193:256] = -vec(mxiy_back[4, 1:end]) .+ (short_mm + far_im)
	mxiy_array[257:320] = -vec(mxiy_back[5, 1:end]) .+ (short_dm + far_im)

	if rand() < thresh
		next_stop = findmax(-mxiy_array)[2]
	else
		next_stop = back_softmax(-mxiy_array, 1, temperature)[1]
	end

	next_state = div(next_stop-1, 64) + 1
	next_letter = div(mod(next_stop-1, 64), 16) + 1
	next_codon = trinucs[mod(next_stop-1, 64) + 1]
	cumulative_score += mxiy_array[next_stop]
	cur_score = mxiy_array[next_stop]
	return next_state, next_letter, next_codon, cur_score, cumulative_score
end

function mxdy(trinucs, last_letter::AbstractString, short_node::Int64, far_node::Int64, cumulative_score::Float64, cur_score::Float64, passed_values::Array{Float64, 1}, bbt, thresh, temperature)
	short_mm, short_md, short_im, short_dm, short_dd, short_mi, short_ii, far_mm, far_md, far_im, far_dm, far_dd, far_mi, far_ii = passed_values;
	mxdy_back = bbt[last_letter][short_node - 1, far_node - 1, [1, 2, 3, 6, 7], 1:64]

	mxdy_array = zeros(5 * 64)

	mxdy_array[1:64] = -vec(mxdy_back[1, 1:end]) .+ (short_mm + far_mm)
	mxdy_array[65:128] = -vec(mxdy_back[2, 1:end]) .+ (short_im + far_mm)
	mxdy_array[129:192] = -vec(mxdy_back[3, 1:end]) .+ (short_dm + far_mm)
	mxdy_array[193:256] = -vec(mxdy_back[4, 1:end]) .+ (short_mm + far_dm)
	mxdy_array[257:320] = -vec(mxdy_back[5, 1:end]) .+ (short_im + far_dm)

	if rand() < thresh
		next_stop = findmax(-mxdy_array)[2]
	else
		next_stop = back_softmax(-mxdy_array, 1, temperature)[1]
	end

	next_state = [1, 2, 3, 6, 7][div(next_stop-1, 64) + 1][1]
	next_letter = div(mod(next_stop-1, 64), 16) + 1
	next_codon = trinucs[mod(next_stop-1, 64) + 1]
	cumulative_score += mxdy_array[next_stop]
	cur_score = mxdy_array[next_stop]
	return next_state, next_letter, next_codon, cur_score, cumulative_score
end

function ixmy(trinucs, last_letter::AbstractString, short_node::Int64, far_node::Int64, cumulative_score::Float64, cur_score::Float64, passed_values::Array{Float64, 1}, bbt, thresh, temperature)
	short_mm, short_md, short_im, short_dm, short_dd, short_mi, short_ii, far_mm, far_md, far_im, far_dm, far_dd, far_mi, far_ii = passed_values;
	ixmy_back = bbt[last_letter][short_node, far_node - 1, [1, 2, 4, 6, 7], 1:64]

	ixmy_array = zeros(5 * 64)

	ixmy_array[1:64] = -vec(ixmy_back[1, 1:end]) .+ (short_mm + far_mm)
	ixmy_array[65:128] = -vec(ixmy_back[2, 1:end]) .+ (short_im + far_mm)
	ixmy_array[129:192] = -vec(ixmy_back[3, 1:end]) .+ (short_dm + far_mm)
	ixmy_array[193:256] = -vec(ixmy_back[4, 1:end]) .+ (short_mm + far_dm)
	ixmy_array[257:320] = -vec(ixmy_back[5, 1:end]) .+ (short_im + far_dm)

	if rand() < thresh
		next_stop = findmax(-ixmy_array)[2]
	else
		next_stop = back_softmax(-ixmy_array, 1, temperature)[1]
	end

	next_state = [1, 2, 4, 6, 7][div(next_stop-1, 64) + 1][1]
	next_letter = div(mod(next_stop-1, 64), 16) + 1
	next_codon = trinucs[mod(next_stop-1, 64) + 1]
	cumulative_score += ixmy_array[next_stop]
	cur_score = ixmy_array[next_stop]
	return next_state, next_letter, next_codon, cur_score, cumulative_score
end

function dxmy(trinucs, last_letter::AbstractString, short_node::Int64, far_node::Int64, cumulative_score::Float64, cur_score::Float64, passed_values::Array{Float64, 1}, bbt, thresh, temperature)
	short_mm, short_md, short_im, short_dm, short_dd, short_mi, short_ii, far_mm, far_md, far_im, far_dm, far_dd, far_mi, far_ii = passed_values;
	dxmy_back = bbt[last_letter][short_node - 1, far_node - 1, [1, 3, 4, 5, 6], 1:64]
	dxmy_array = zeros(5 * 64)

	dxmy_array[1:64] = -vec(dxmy_back[1, 1:end]) .+ (short_mm + far_mm)
	dxmy_array[65:128] = -vec(dxmy_back[2, 1:end]) .+ (short_dm + far_mm)
	dxmy_array[129:192] = -vec(dxmy_back[3, 1:end]) .+ (short_mm + far_im)
	dxmy_array[193:256] = -vec(dxmy_back[4, 1:end]) .+ (short_dm + far_im)
	dxmy_array[257:320] = -vec(dxmy_back[5, 1:end]) .+ (short_mm + far_dm)

	if rand() < thresh
		next_stop = findmax(-dxmy_array)[2]
	else
		next_stop = back_softmax(-dxmy_array, 1, temperature)[1]
	end

	next_state = [1, 3, 4, 5, 6][div(next_stop - 1, 64) + 1][1]
	next_letter = div(mod(next_stop-1, 64), 16) + 1
	next_codon = trinucs[mod(next_stop-1, 64) + 1]
	cumulative_score += dxmy_array[next_stop]
	cur_score = dxmy_array[next_stop]
	return next_state, next_letter, next_codon, cur_score, cumulative_score
end

function ixdy(trinucs, last_letter::AbstractString, short_node::Int64, far_node::Int64, cumulative_score::Float64, cur_score::Float64, passed_values::Array{Float64, 1}, bbt, thresh, temperature)
	short_mm, short_md, short_im, short_dm, short_dd, short_mi, short_ii, far_mm, far_md, far_im, far_dm, far_dd, far_mi, far_ii = passed_values;
	ixdy_back = bbt[last_letter][short_node, far_node - 1, [1, 2, 6, 7], 1:64]

	ixdy_array = zeros(4 * 64)

	ixdy_array[1:64] = -vec(ixdy_back[1, 1:end]) .+ (short_mm + far_mm)
	ixdy_array[65:128] = -vec(ixdy_back[2, 1:end]) .+ (short_im + far_mm)
	ixdy_array[129:192] = -vec(ixdy_back[3, 1:end]) .+ (short_mm + far_dm)
	ixdy_array[193:256] = -vec(ixdy_back[4, 1:end]) .+ (short_im + far_dm)

	if rand() < thresh
		next_stop = findmax(-ixdy_array)[2]
	else
		next_stop = back_softmax(-ixdy_array, 1, temperature)[1]
	end

	next_state = [1, 2, 6, 7][div(next_stop - 1, 64) + 1][1]
	next_letter = div(mod(next_stop-1, 64), 16) + 1
	next_codon = trinucs[mod(next_stop-1, 64) + 1]
	cumulative_score += ixdy_array[next_stop]
	cur_score = ixdy_array[next_stop]
	return next_state, next_letter, next_codon, cur_score, cumulative_score
end

function dxiy(trinucs, last_letter::AbstractString, short_node::Int64, far_node::Int64, cumulative_score::Float64, cur_score::Float64, passed_values::Array{Float64, 1}, bbt, thresh, temperature)
	short_mm, short_md, short_im, short_dm, short_dd, short_mi, short_ii, far_mm, far_md, far_im, far_dm, far_dd, far_mi, far_ii = passed_values;
	dxiy_back = bbt[last_letter][short_node - 1, far_node, [1, 3, 4, 5], 1:64]

	dxiy_array = zeros(4 * 64)

	dxiy_array[1:64] = -vec(dxiy_back[1, 1:end]) .+ (short_mm + far_mm)
	dxiy_array[65:128] = -vec(dxiy_back[2, 1:end]) .+ (short_dm + far_mm)
	dxiy_array[129:192] = -vec(dxiy_back[3, 1:end]) .+ (short_mm + far_im)
	dxiy_array[193:256] = -vec(dxiy_back[4, 1:end]) .+ (short_dm + far_im)

	if rand() < thresh
		next_stop = findmax(-dxiy_array)[2]
	else
		next_stop = back_softmax(-dxiy_array, 1, temperature)[1]
	end

	next_state = [1, 3, 4, 5][div(next_stop - 1, 64) + 1][1]
	next_letter = div(mod(next_stop-1, 64), 16) + 1
	next_codon = trinucs[mod(next_stop-1, 64) + 1]

	cumulative_score += dxiy_array[next_stop]
	cur_score = dxiy_array[next_stop]
	return next_state, next_letter, next_codon, cur_score, cumulative_score
end

end


text:
function factor(n)
    answer = []
    hw = []
    for i in eratosphen(n)
        if n%i == 0
            push!(answer,i)
            push!(hw,1)
            n = n/i
            while n%i == 0
                hw[end] += 1
                n = n/i
            end
        end
    end
    return answer,hw
end

function eratosphen(n)
    ser = fill(true,n)
    ser[1] = false
    k = 2
    while k !== nothing && k<n
        ser[k^2:k:end] .= false
        k=findnext(ser, k+1)
    end
    return findall(ser)
end


text:
module Sequences

import MutableArithmetics
const MA = MutableArithmetics

export shortest_common_supersequence

"""
Given arrays (or strings) a and b, returns a table of size
(length(a) + 1) x (length(b) + 1)
Each entry i,j in the table is the length of the longest
common subsequence among a[1:(i-1)] and b[1:(j-1)]
Adapted from https://github.com/WestleyArgentum/Subsequences.jl
Copyright 2015 WestleyArentum and other contributors, released
under the MIT "Expat" License.
"""
function longest_common_subsequence_table(a, b)
    lengths = zeros(Int, length(a) + 1, length(b) + 1)
    for i in 1:length(a)
        ai = a[i]
        for j in 1:length(b)
            if ai == b[j]
                lengths[i+1, j+1] = lengths[i, j] + 1
            else
                lengths[i+1, j+1] = max(lengths[i+1, j], lengths[i, j+1])
            end
        end
    end
    lengths
end

"""
Given arrays (or strings) a and b, returns a table of size
(length(a) + 1) x (length(b) + 1)
Each entry i,j in the table is the length of the shortest
common supersequence among a[1:(i-1)] and b[1:(j-1)]
"""
function shortest_common_supersequence_table(a, b)
    table = longest_common_subsequence_table(a, b)
    for i in 1:size(table, 1)
        for j in 1:size(table, 2)
            table[i, j] = (i - 1) + (j - 1) - table[i, j]
        end
    end
    table
end

"""
Given arrays a and b, returns the shortest
common supersequence of a and b.
See: https://www.ics.uci.edu/~eppstein/161/960229.html
"""
function shortest_common_supersequence(a::Union{NTuple, AbstractVector}, b::Union{NTuple, AbstractVector})
    table = shortest_common_supersequence_table(a, b)
    result = promote_type(eltype(a), eltype(b))[]
    i = size(table, 1) - 1
    j = size(table, 2) - 1
    while i > 0 || j > 0
        if i == 0
            push!(result, b[j])
            j -= 1
        elseif j == 0
            push!(result, a[i])
            i -= 1
        elseif a[i] == b[j]
            push!(result, a[i])
            i -= 1
            j -= 1
        elseif table[i, j + 1] > table[i + 1, j]
            push!(result, b[j])
            j -= 1
        else
            push!(result, a[i])
            i -= 1
        end
    end
    reverse!(result)
    result
end

function shortest_common_supersequence(a::AbstractString, b::AbstractString)
    scs = shortest_common_supersequence(collect(a), collect(b))
    join(scs)
end

function mergesorted!(result, v1::AbstractArray, v2::AbstractArray, isless, combine, filter=x -> !iszero(x))
    i = 1
    i1 = 1
    i2 = 1
    while i1 <= length(v1) && i2 <= length(v2)
        x1 = v1[i1]
        x2 = v2[i2]
        if isless(x1, x2)
            if filter(x1)
                result[i] = x1
                i += 1
            end
            i1 += 1
        elseif isless(x2, x1)
            if filter(x2)
                result[i] = x2
                i += 1
            end
            i2 += 1
        else
            c = combine(x1, x2)
            if filter(c)
                result[i] = combine(x1, x2)
                i += 1
            end
            i1 += 1
            i2 += 1
        end
    end
    for j in i1:length(v1)
        if filter(v1[j])
            result[i] = v1[j]
            i += 1
        end
    end
    for j in i2:length(v2)
        if filter(v2[j])
            result[i] = v2[j]
            i += 1
        end
    end
    resize!(result, i - 1)
    result
end

function mergesorted(v1::AbstractArray, v2::AbstractArray, isless=Base.isless,
                     combine=Base.:+, filter=x -> !iszero(x))
    T = MA.promote_operation(combine, eltype(v1), eltype(v2))
    result = Vector{T}(undef, length(v1) + length(v2))
    mergesorted!(result, v1, v2, isless, combine, filter)
end

end


text:
module OSQP

export OSQPMathProgBaseInterface
using SparseArrays
using LinearAlgebra

using OSQP_jll

include("constants.jl")
include("types.jl")
include("interface.jl")
include("MOI_wrapper.jl")
const Optimizer = MathOptInterfaceOSQP.Optimizer

end # module


text:
using StrTables, Unicode_Entities

@static VERSION < v"0.7.0-DEV" ? (using Base.Test) : (using Test)

# Test the functions lookupname, matches, longestmatches, completions
# Check that characters from all 3 tables (BMP, non-BMP, 2 character) are tested

const def = Unicode_Entities.default

const datapath = "../data"
const dpath = "ftp://ftp.unicode.org/Public/UNIDATA/"
const fname = "UnicodeData.txt"

const symnam = Vector{String}()
const symval = Vector{Char}()

"""Load up all names and characters from original data file"""
function load_unicode_data()
    lname = joinpath(datapath, fname)
    isfile(lname) || download(string(dpath, fname), lname)
    aliasnam = Vector{String}()
    aliasval = Vector{Char}()
    open(lname, "r") do f
        while (l = chomp(readline(f))) != ""
            flds = split(l, ";")
            str = flds[2]
            alias = flds[11]
            ch = Char(parse_hex(UInt32, flds[1]))
            if str[1] == '<'
                str != "<control>" && continue
                str = ""
            end
            if str != ""
                push!(symnam, str)
                push!(symval, ch)
            end
            if alias != ""
                push!(aliasnam, alias)
                push!(aliasval, ch)
            end
        end
    end
    # Check for duplicates
    names = Set{String}(symnam)
    for (str,ch) in zip(aliasnam, aliasval)
        if !(str in names)
            push!(symnam, str)
            push!(symval, ch)
        end
    end
end

load_unicode_data()

@testset "Unicode_Entities" begin

    @testset "matches data file" begin
        for (i, ch) in enumerate(symval)
            list = matchchar(def, ch)
            if !isempty(list)
                @test symnam[i] in list
            end
        end
        for (i, nam) in enumerate(symnam)
            str = lookupname(def, nam)
            if str != ""
                @test symval[i] == str[1]
            end
        end
    end

@testset "lookupname" begin
    @test lookupname(def, "foobar")   == ""
    @test lookupname(def, SubString("My name is Spock", 12)) == ""
    @test lookupname(def, "end of text") == "\x03" # \3
    @test lookupname(def, "TIBETAN LETTER -A") == "\u0f60"
    @test lookupname(def, "LESS-THAN OR SLANTED EQUAL TO") == "\u2a7d"
    @test lookupname(def, "REVERSED HAND WITH MIDDLE FINGER EXTENDED") == "\U1f595"
end

@testset "matches" begin
    @test isempty(matches(def, ""))
    @test isempty(matches(def, "\uf900"))
    @test isempty(matches(def, SubString("This is \uf900", 9)))
    for (chrs, exp) in (("\U1f596", ["RAISED HAND WITH PART BETWEEN MIDDLE AND RING FINGERS"]),
                        ("\u0f4a", ["TIBETAN LETTER REVERSED TA"]),
                        (".", ["FULL STOP", "PERIOD"]))
        res = matches(def, chrs)
        @test length(res) >= length(exp)
        @test intersect(res, exp) == exp
    end
end

@testset "longestmatches" begin
    @test isempty(longestmatches(def, "\uf900 abcd"))
    @test isempty(longestmatches(def, SubString("This is \uf900 abcd", 9)))
    for (chrs, exp) in (("\U1f596 abcd", ["RAISED HAND WITH PART BETWEEN MIDDLE AND RING FINGERS"]),
                        (".abcd", ["FULL STOP", "PERIOD"]),
                        ("\u0f4a#123", ["TIBETAN LETTER REVERSED TA", "TIBETAN LETTER TTA"]))
        res = longestmatches(def, chrs)
        @test length(res) >= length(exp)
        @test intersect(res, exp) == exp
    end
end

@testset "completions" begin
    @test isempty(completions(def, "ScottPaulJones"))
    @test isempty(completions(def, SubString("My name is Scott", 12)))
    for (chrs, exp) in (("ZERO", ["ZERO WIDTH JOINER", "ZERO WIDTH NO-BREAK SPACE",
                                  "ZERO WIDTH NON-JOINER", "ZERO WIDTH SPACE"]),
                        ("BACK OF", ["BACK OF ENVELOPE"]))
        res = completions(def, chrs)
        @test length(res) >= length(exp)
        @test intersect(res, exp) == exp
    end
end
end


text:
# Autogenerated wrapper script for CyrusSASL_jll for i686-w64-mingw32
export libanonymous, libcrammd5, libdigestmd5, libntlm, libotp, libplan, libsasl2, libscram

using OpenSSL_jll
using SQLite_jll
using Gdbm_jll
JLLWrappers.@generate_wrapper_header("CyrusSASL")
JLLWrappers.@declare_library_product(libanonymous, "libanonymous-3.dll")
JLLWrappers.@declare_library_product(libcrammd5, "libcrammd5-3.dll")
JLLWrappers.@declare_library_product(libdigestmd5, "libdigestmd5-3.dll")
JLLWrappers.@declare_library_product(libntlm, "libntlm-3.dll")
JLLWrappers.@declare_library_product(libotp, "libotp-3.dll")
JLLWrappers.@declare_library_product(libplan, "libplain-3.dll")
JLLWrappers.@declare_library_product(libsasl2, "libsasl2-3.dll")
JLLWrappers.@declare_library_product(libscram, "libscram-3.dll")
function __init__()
    JLLWrappers.@generate_init_header(OpenSSL_jll, SQLite_jll, Gdbm_jll)
    JLLWrappers.@init_library_product(
        libanonymous,
        "lib/sasl2/libanonymous-3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcrammd5,
        "lib/sasl2/libcrammd5-3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libdigestmd5,
        "lib/sasl2/libdigestmd5-3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libntlm,
        "lib/sasl2/libntlm-3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libotp,
        "lib/sasl2/libotp-3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libplan,
        "lib/sasl2/libplain-3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsasl2,
        "bin/libsasl2-3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libscram,
        "lib/sasl2/libscram-3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
fname = :ST_ARCsparse
c = Combi(hessian_sparse,PDataST,solve_modelST_ARC,preprocessST,decreaseGen,TparamsST())
include("Template.jl")


text:
board = [[" ", " ", " "],
         [" ", " ", " "],
         [" ", " ", " "]]

player1 = "X"
player2 = "O"

player1_win = -1
tie = 0
player2_win = 1

function check_win(board)

    # Rows
    for i in 1:3
        if (board[i][1] == board[i][2] == board[i][3]) && (board[i][1] in [player1, player2])
            if (board[i][1] == player1)
                return player1_win
            end
            return player2_win
        end
    end

    # Column
    for i in 1:3
        if (board[1][i] == board[2][i] == board[3][i]) && (board[1][i] in [player1, player2])
            if (board[1][i] == player1)
                return player1_win
            end
            return player2_win
        end
    end

    # Principal diagonal
    if (board[1][1] == board[2][2] == board[3][3]) && (board[1][1] in [player1, player2])
        if (board[1][1] == player1)
            return player1_win
        end
        return player2_win
    end

    # Other diag
    if (board[1][3] == board[2][2] == board[3][1]) && (board[1][3] in [player1, player2])
        if (board[1][3] == player1)
            return player1_win
        end
        return player2_win
    end

    # Tie
    if length(check_available(board)) == 0
        return tie
    end

    # No one wins, game still goes on
    return nothing
end

function minimax(board, depth, is_maximising)
    result = check_win(board)

    if result !== nothing
        return result
    elseif is_maximising

        best_score = -Inf

        for i in 1:3, j in 1:3
            if board[i][j] == " "

                board[i][j] = player2

                score = minimax(board, depth + 1, false)

                board[i][j] = " "

                best_score = max(score, best_score)
            end
        end

        return best_score
    else
        best_score = Inf

        for i in 1:3, j in 1:3
            if board[i][j] == " "

                board[i][j] = player1

                score = minimax(board, depth + 1, true)

                board[i][j] = " "

                best_score = min(score, best_score)
            end
        end

        return best_score
    end

end

function best_move(board)
    best_score = -Inf

    best_pos = (0, 0)

    for i in 1:3, j in 1:3
        if board[i][j] == " "

            board[i][j] = player2

            score = minimax(board, 0, false)

            board[i][j] = " "

            if score > best_score
                best_score = score
                best_pos = (i, j)
            end
        end
    end
    return best_pos
end

function print_board(board)
    println("    1   2   3")
    println(" 1  ", join(board[1], " | "))
    println("   ———┼———┼———")
    println(" 2  ", join(board[2], " | "))
    println("   ———┼———┼———")
    println(" 3  ", join(board[3], " | "))
end

function make_move!(board, player, (x, y))
    board[x][y] = player
end

function check_available(board)
    available_cells = []
    for i in 1:3, j in 1:3
        if board[i][j] == " "
            push!(available_cells, (i, j))
        end
    end
    return available_cells
end

function ai_move!(board)
    make_move!(board, player2, best_move(board))
end

print_board(board)

while true
    print_board(board)

    print("Enter: ")
    txt = split(readline())
    make_move!(board, player1, (parse(Int, txt[1]), parse(Int, txt[2])))

    if check_win(board) == player1_win
        println(1)
        break
    elseif check_win(board) == player2_win
        println(2)
        break
    elseif check_win(board) == tie
        println("Tie!")
        break
    end

    ai_move!(board)
    print_board(board)

    if check_win(board) == player1_win
        println(1)
        break
    elseif check_win(board) == player2_win
        println(2)
        break
    elseif check_win(board) == tie
        println("Tie!")
        break
    end
end

max_stars_repo_name: soumitradev/youtube-tutorial

max_issues_repo_name: soumitradev/youtube-tutorial

max_forks_repo_name: soumitradev/youtube-tutorial


text:
function time_unit(t)
    units = ["hr", "min", "s", "ms", "us", "ns", "ps"]
    sizes = [60.0 * 60.0, 60.0, 1.0, 1e-3, 1e-6, 1e-9, 1e-12]

    for (unit, size) in zip(units, sizes)
        if t / 10 > size
            return unit, size
        end
    end

    return units[end], sizes[end]
end

function fmt_time(ts)
    avg = sum(ts) / length(ts)
    unit, size = time_unit(avg)

    return @sprintf "%d %s" (avg / size) unit
end

function bench_graph_union(nodes :: Int, edges :: Int)
    us = [convert(Int, floor(nodes * rand()) + 1) for e in 1:edges]
    vs = [convert(Int, floor(nodes * rand()) + 1) for e in 1:edges]

    uf = UnionFinder(nodes)
    t = @elapsed union!(uf, us, vs)

    return uf, t
end

function bench_grid_union(nodes :: Int, edges :: Int)
    width = convert(Int, ceil(sqrt(nodes)))

    us = [convert(Int, floor(nodes * rand()) + 1) for e in 1:edges]

    dirs = [1, nodes - 1, width, nodes - width]
    vs = [(dirs[convert(Int, floor(4 * rand()) + 1)] + us[e]) % nodes + 1
          for e in 1:edges]

    uf = UnionFinder(nodes)
    t = @elapsed union!(uf, us, vs)

    return uf, t
end

function bench_avg(nodes :: Int, frac :: Float64, sweeps :: Int, f :: Function)
    uts = Vector{Float64}(undef, sweeps)
    cts = Vector{Float64}(undef, sweeps)

    edges = convert(Int, ceil(nodes * frac))

    for i in 1:sweeps
        uf, uts[i] = f(nodes, edges)
        cts[i] = @elapsed CompressedFinder(uf)

        uts[i] = edges == 0 ? uts[i] : uts[i] / edges
        cts[i] = cts[i] / nodes
    end

    @printf "%12s/edge %12s/node\n" fmt_time(uts) fmt_time(cts)
end

function benchmark_main()
    @printf "%30s%17s %17s\n" " " "UnionFinder" "CompressedFinder"

    @printf "%30s" "Sparse 1,000 node graph"
    bench_avg(1000, 0.1, 1000, bench_graph_union)
    @printf "%30s" "Dense 1,000 node graph"
    bench_avg(1000, 0.8, 1000, bench_graph_union)
    @printf "%30s" "Sparse 1,000,000 node graph"
    bench_avg(1000 * 1000, 0.1, 10, bench_graph_union)
    @printf "%30s" "Dense 1,000,000 node graph"
    bench_avg(1000 * 1000, 0.8, 10, bench_graph_union)
    println()
    @printf "%30s" "Sparse 1,000 node grid"
    bench_avg(1000, 0.1, 1000, bench_grid_union)
    @printf "%30s" "Dense 1,000 node grid"
    bench_avg(1000, 0.8, 1000, bench_grid_union)
    @printf "%30s" "Sparse 1,000,000 node grid"
    bench_avg(1000 * 1000, 0.1, 10, bench_grid_union)
    @printf "%30s" "Dense 1,000,000 node grid"
    bench_avg(1000 * 1000, 0.8, 10, bench_grid_union)
end

benchmark_main()


text:
# Julia wrapper for header: /usr/include/libavformat/avformat.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_get_packet,
    av_append_packet,
    #avformat_version,
    avformat_configuration,
    avformat_license,
    av_register_all,
    av_register_input_format,
    av_register_output_format,
    avformat_network_init,
    avformat_network_deinit,
    av_iformat_next,
    av_oformat_next,
    avformat_alloc_context,
    avformat_free_context,
    avformat_get_class,
    avformat_new_stream,
    av_new_program,
    av_find_input_format,
    av_probe_input_format,
    av_probe_input_format2,
    av_probe_input_buffer,
    avformat_open_input,
    avformat_find_stream_info,
    av_find_best_stream,
    av_read_packet,
    av_read_frame,
    av_seek_frame,
    avformat_seek_file,
    av_read_play,
    av_read_pause,
    av_close_input_file,
    avformat_close_input,
    avformat_write_header,
    av_write_frame,
    av_interleaved_write_frame,
    av_interleave_packet_per_dts,
    av_write_trailer,
    av_guess_format,
    av_guess_codec,
    av_hex_dump,
    av_hex_dump_log,
    av_pkt_dump2,
    av_pkt_dump_log2,
    av_codec_get_id,
    av_codec_get_tag,
    av_find_default_stream_index,
    av_index_search_timestamp,
    av_add_index_entry,
    av_url_split,
    av_dump_format,
    av_get_frame_filename,
    av_filename_number_test,
    av_sdp_create,
    av_match_ext,
    avformat_query_codec,
    avformat_get_riff_video_tags,
    avformat_get_riff_audio_tags

function av_get_packet(s,pkt,size::Integer)
    ccall((:av_get_packet,libavformat),Cint,(Ptr{AVIOContext},Ptr{AVPacket},Cint),s,pkt,size)
end

function av_append_packet(s,pkt,size::Integer)
    ccall((:av_append_packet,libavformat),Cint,(Ptr{AVIOContext},Ptr{AVPacket},Cint),s,pkt,size)
end

function avformat_version()
    ccall((:avformat_version,libavformat),UInt32,())
end

function avformat_configuration()
    ccall((:avformat_configuration,libavformat),Ptr{UInt8},())
end

function avformat_license()
    ccall((:avformat_license,libavformat),Ptr{UInt8},())
end

function av_register_all()
    ccall((:av_register_all,libavformat),Void,())
end

function av_register_input_format(format)
    ccall((:av_register_input_format,libavformat),Void,(Ptr{AVInputFormat},),format)
end

function av_register_output_format(format)
    ccall((:av_register_output_format,libavformat),Void,(Ptr{AVOutputFormat},),format)
end

function avformat_network_init()
    ccall((:avformat_network_init,libavformat),Cint,())
end

function avformat_network_deinit()
    ccall((:avformat_network_deinit,libavformat),Cint,())
end

function av_iformat_next(f)
    ccall((:av_iformat_next,libavformat),Ptr{AVInputFormat},(Ptr{AVInputFormat},),f)
end

function av_oformat_next(f)
    ccall((:av_oformat_next,libavformat),Ptr{AVOutputFormat},(Ptr{AVOutputFormat},),f)
end

function avformat_alloc_context()
    ccall((:avformat_alloc_context,libavformat),Ptr{AVFormatContext},())
end

function avformat_free_context(s)
    ccall((:avformat_free_context,libavformat),Void,(Ptr{AVFormatContext},),s)
end

function avformat_get_class()
    ccall((:avformat_get_class,libavformat),Ptr{AVClass},())
end

function avformat_new_stream(s,c)
    ccall((:avformat_new_stream,libavformat),Ptr{AVStream},(Ptr{AVFormatContext},Ptr{AVCodec}),s,c)
end

function av_new_program(s,id::Integer)
    ccall((:av_new_program,libavformat),Ptr{AVProgram},(Ptr{AVFormatContext},Cint),s,id)
end

function av_find_input_format(short_name)
    ccall((:av_find_input_format,libavformat),Ptr{AVInputFormat},(Ptr{UInt8},),short_name)
end

function av_probe_input_format(pd,is_opened::Integer)
    ccall((:av_probe_input_format,libavformat),Ptr{AVInputFormat},(Ptr{AVProbeData},Cint),pd,is_opened)
end

function av_probe_input_format2(pd,is_opened::Integer,score_max)
    ccall((:av_probe_input_format2,libavformat),Ptr{AVInputFormat},(Ptr{AVProbeData},Cint,Ptr{Cint}),pd,is_opened,score_max)
end

function av_probe_input_buffer(pb,fmt,filename,logctx,offset::Integer,max_probe_size::Integer)
    ccall((:av_probe_input_buffer,libavformat),Cint,(Ptr{AVIOContext},Ptr{Ptr{AVInputFormat}},Ptr{UInt8},Ptr{Void},UInt32,UInt32),pb,fmt,filename,logctx,offset,max_probe_size)
end

function avformat_open_input(ps,filename,fmt,options)
    ccall((:avformat_open_input,libavformat),Cint,(Ptr{Ptr{AVFormatContext}},Ptr{UInt8},Ptr{AVInputFormat},Ptr{Ptr{AVDictionary}}),ps,filename,fmt,options)
end

function avformat_find_stream_info(ic,options)
    ccall((:avformat_find_stream_info,libavformat),Cint,(Ptr{AVFormatContext},Ptr{Ptr{AVDictionary}}),ic,options)
end

function av_find_best_stream(ic,_type::AVMediaType,wanted_stream_nb::Integer,related_stream::Integer,decoder_ret,flags::Integer)
    ccall((:av_find_best_stream,libavformat),Cint,(Ptr{AVFormatContext},AVMediaType,Cint,Cint,Ptr{Ptr{AVCodec}},Cint),ic,_type,wanted_stream_nb,related_stream,decoder_ret,flags)
end

function av_read_packet(s,pkt)
    ccall((:av_read_packet,libavformat),Cint,(Ptr{AVFormatContext},Ptr{AVPacket}),s,pkt)
end

function av_read_frame(s,pkt)
    ccall((:av_read_frame,libavformat),Cint,(Ptr{AVFormatContext},Ptr{AVPacket}),s,pkt)
end

function av_seek_frame(s,stream_index::Integer,timestamp::Int64,flags::Integer)
    ccall((:av_seek_frame,libavformat),Cint,(Ptr{AVFormatContext},Cint,Int64,Cint),s,stream_index,timestamp,flags)
end

function avformat_seek_file(s,stream_index::Integer,min_ts::Int64,ts::Int64,max_ts::Int64,flags::Integer)
    ccall((:avformat_seek_file,libavformat),Cint,(Ptr{AVFormatContext},Cint,Int64,Int64,Int64,Cint),s,stream_index,min_ts,ts,max_ts,flags)
end

function av_read_play(s)
    ccall((:av_read_play,libavformat),Cint,(Ptr{AVFormatContext},),s)
end

function av_read_pause(s)
    ccall((:av_read_pause,libavformat),Cint,(Ptr{AVFormatContext},),s)
end

function av_close_input_file(s)
    ccall((:av_close_input_file,libavformat),Void,(Ptr{AVFormatContext},),s)
end

function avformat_close_input(s)
    ccall((:avformat_close_input,libavformat),Void,(Ptr{Ptr{AVFormatContext}},),s)
end

function avformat_write_header(s,options)
    ccall((:avformat_write_header,libavformat),Cint,(Ptr{AVFormatContext},Ptr{Ptr{AVDictionary}}),s,options)
end

function av_write_frame(s,pkt)
    ccall((:av_write_frame,libavformat),Cint,(Ptr{AVFormatContext},Ptr{AVPacket}),s,pkt)
end

function av_interleaved_write_frame(s,pkt)
    ccall((:av_interleaved_write_frame,libavformat),Cint,(Ptr{AVFormatContext},Ptr{AVPacket}),s,pkt)
end

function av_interleave_packet_per_dts(s,out,pkt,flush::Integer)
    ccall((:av_interleave_packet_per_dts,libavformat),Cint,(Ptr{AVFormatContext},Ptr{AVPacket},Ptr{AVPacket},Cint),s,out,pkt,flush)
end

function av_write_trailer(s)
    ccall((:av_write_trailer,libavformat),Cint,(Ptr{AVFormatContext},),s)
end

function av_guess_format(short_name,filename,mime_type)
    ccall((:av_guess_format,libavformat),Ptr{AVOutputFormat},(Ptr{UInt8},Ptr{UInt8},Ptr{UInt8}),short_name,filename,mime_type)
end

function av_guess_codec(fmt,short_name,filename,mime_type,_type::AVMediaType)
    ccall((:av_guess_codec,libavformat),Cint,(Ptr{AVOutputFormat},Ptr{UInt8},Ptr{UInt8},Ptr{UInt8},AVMediaType),fmt,short_name,filename,mime_type,_type)
end

function av_hex_dump(f,buf,size::Integer)
    ccall((:av_hex_dump,libavformat),Void,(Ptr{FILE},Ptr{UInt8},Cint),f,buf,size)
end

function av_hex_dump_log(avcl,level::Integer,buf,size::Integer)
    ccall((:av_hex_dump_log,libavformat),Void,(Ptr{Void},Cint,Ptr{UInt8},Cint),avcl,level,buf,size)
end

function av_pkt_dump2(f,pkt,dump_payload::Integer,st)
    ccall((:av_pkt_dump2,libavformat),Void,(Ptr{FILE},Ptr{AVPacket},Cint,Ptr{AVStream}),f,pkt,dump_payload,st)
end

function av_pkt_dump_log2(avcl,level::Integer,pkt,dump_payload::Integer,st)
    ccall((:av_pkt_dump_log2,libavformat),Void,(Ptr{Void},Cint,Ptr{AVPacket},Cint,Ptr{AVStream}),avcl,level,pkt,dump_payload,st)
end

function av_codec_get_id(tags,tag::Integer)
    ccall((:av_codec_get_id,libavformat),Cint,(Ptr{Ptr{AVCodecTag}},UInt32),tags,tag)
end

function av_codec_get_tag(tags,id::AVCodecID)
    ccall((:av_codec_get_tag,libavformat),UInt32,(Ptr{Ptr{AVCodecTag}},AVCodecID),tags,id)
end

function av_find_default_stream_index(s)
    ccall((:av_find_default_stream_index,libavformat),Cint,(Ptr{AVFormatContext},),s)
end

function av_index_search_timestamp(st,timestamp::Int64,flags::Integer)
    ccall((:av_index_search_timestamp,libavformat),Cint,(Ptr{AVStream},Int64,Cint),st,timestamp,flags)
end

function av_add_index_entry(st,pos::Int64,timestamp::Int64,size::Integer,distance::Integer,flags::Integer)
    ccall((:av_add_index_entry,libavformat),Cint,(Ptr{AVStream},Int64,Int64,Cint,Cint,Cint),st,pos,timestamp,size,distance,flags)
end

function av_url_split(proto,proto_size::Integer,authorization,authorization_size::Integer,hostname,hostname_size::Integer,port_ptr,path,path_size::Integer,url)
    ccall((:av_url_split,libavformat),Void,(Ptr{UInt8},Cint,Ptr{UInt8},Cint,Ptr{UInt8},Cint,Ptr{Cint},Ptr{UInt8},Cint,Ptr{UInt8}),proto,proto_size,authorization,authorization_size,hostname,hostname_size,port_ptr,path,path_size,url)
end

function av_dump_format(ic,index::Integer,url,is_output::Integer)
    ccall((:av_dump_format,libavformat),Void,(Ptr{AVFormatContext},Cint,Ptr{UInt8},Cint),ic,index,url,is_output)
end

function av_get_frame_filename(buf,buf_size::Integer,path,number::Integer)
    ccall((:av_get_frame_filename,libavformat),Cint,(Ptr{UInt8},Cint,Ptr{UInt8},Cint),buf,buf_size,path,number)
end

function av_filename_number_test(filename)
    ccall((:av_filename_number_test,libavformat),Cint,(Ptr{UInt8},),filename)
end

function av_sdp_create(ac,n_files::Integer,buf,size::Integer)
    ccall((:av_sdp_create,libavformat),Cint,(Ptr{Ptr{AVFormatContext}},Cint,Ptr{UInt8},Cint),ac,n_files,buf,size)
end

function av_match_ext(filename,extensions)
    ccall((:av_match_ext,libavformat),Cint,(Ptr{UInt8},Ptr{UInt8}),filename,extensions)
end

function avformat_query_codec(ofmt,codec_id::AVCodecID,std_compliance::Integer)
    ccall((:avformat_query_codec,libavformat),Cint,(Ptr{AVOutputFormat},AVCodecID,Cint),ofmt,codec_id,std_compliance)
end

function avformat_get_riff_video_tags()
    ccall((:avformat_get_riff_video_tags,libavformat),Ptr{AVCodecTag},())
end

function avformat_get_riff_audio_tags()
    ccall((:avformat_get_riff_audio_tags,libavformat),Ptr{AVCodecTag},())
end


text:
struct Difference{PA, PB}
    a::PA
    b::PB
end

Base.:*(n::Number, d::Difference) = Difference(*(n, d.a), *(n, d.b))
Base.:+(d1::Difference, d2::Difference) = Difference(d1.a + d2.a, d1.b + d2.b)
any_inside(d::Difference) = Difference(any_inside(d.a), any_inside(d.b))

struct CollisionCache{GeomA, GeomB, M, D <: Difference}
    bodyA::GeomA
    bodyB::GeomB
    simplex_points::MVector{M, D}
end

function CollisionCache(geomA, geomB)
    N = dimension(typeof(geomA))
    @assert dimension(typeof(geomA)) == dimension(typeof(geomB))
    CollisionCache(N, geomA, geomB)
end

function CollisionCache(::Val{N}, geomA, geomB) where N
    interior_point = any_inside(Difference(geomA, geomB))
    simplex = MVector(ntuple(i -> interior_point, Val(N + 1)))
    CollisionCache(geomA, geomB, simplex)
end

function reset!(cache::CollisionCache)
    interior_point = any_inside(Difference(cache.bodyA, cache.bodyB))
    simplex = cache.simplex_points
    @inbounds for i in eachindex(simplex)
        simplex[i] = interior_point
    end
    cache
end

dimension(::Type{CollisionCache{G1, G2, M, D}}) where {G1,G2,M,D} = dimension(G1)

function support_vector_max(geometry::gt.GeometryPrimitive, direction, initial_guess::Tagged)
    best_pt, score = gt.support_vector_max(geometry, gt.Vec(direction))
    Tagged(SVector(best_pt))
end

function support_vector_max(pt::gt.Vec{N, T}, direction, initial_guess::Tagged) where {N,T}
    Tagged(SVector(pt))
end

function support_vector_max(pt::gt.Point{N, T}, direction, initial_guess::Tagged) where {N,T}
    Tagged(SVector(pt))
end

function support_vector_max(simplex::Union{gt.AbstractSimplex, gt.AbstractFlexibleGeometry}, direction, initial_guess::Tagged)
    best_pt, score = gt.support_vector_max(simplex, gt.Vec(direction))
    Tagged(SVector(best_pt))
end

function support_vector_max(mesh::gt.HomogenousMesh{gt.Point{N, T}}, direction, initial_guess::Tagged) where {N,T}
    best_arg, best_value = gt.argmax(x-> dot(SVector(x), direction), gt.vertices(mesh))
    best_vec = SVector(best_arg)
    Tagged(best_vec)
end

any_inside(pt::SVector) = Tagged(pt)
support_vector_max(pt::SVector, direction, initial_guess::Tagged) = Tagged(pt)

function transform_simplex(cache::CollisionCache, poseA, poseB)
    transform_simplex(dimension(typeof(cache)), cache, poseA, poseB)
end

@generated function transform_simplex(::Val{N}, cache::CollisionCache, poseA, poseB) where {N}
    transform_simplex_impl(N, cache, poseA, poseB)
end

function transform_simplex_impl(N, cache, poseA, poseB)
    Expr(:call, :(SVector),
        [:((poseA(value(cache.simplex_points[$i].a)) -
            poseB(value(cache.simplex_points[$i].b)))) for i in 1:(N + 1)]...)
end

# Note: it looks like this can be replaced with transpose(weights) * points in Julia 1.3 (before that, it's a lot slower)
@generated function linear_combination(weights::StaticVector{N}, points::StaticVector{N}) where {N}
    expr = :(weights[1] * points[1])
    for i = 2 : N
        expr = :($expr + weights[$i] * points[$i])
    end
    return quote
        Base.@_inline_meta
        $expr
    end
end

struct GJKResult{M, N, T}
    simplex::SVector{M, SVector{N, T}}
    weights::SVector{M, T}
    in_collision::Bool
    closest_point_in_body::Difference{SVector{N, T}, SVector{N, T}}
    iterations::Int
end

function Base.getproperty(result::GJKResult, sym::Symbol)
    if sym === :signed_distance
        @warn "The `signed_distance` field was removed. Please use the `separation_distance` and `simplex_penetration_distance` functions."
        if result.in_collision
            return -simplex_penetration_distance(result)
        else
            return separation_distance(result)
        end
    else
        return Core.getfield(result, sym)
    end
end

closest_point_in_world(result::GJKResult) = linear_combination(result.simplex, result.weights)
closest_point_in_body(result::GJKResult) = result.closest_point_in_body # just for symmetry with the world function
separation_distance(result::GJKResult) = (@assert !result.in_collision; norm(closest_point_in_world(result)))
simplex_penetration_distance(result::GJKResult) = penetration_distance(result.simplex)

function gjk!(cache::CollisionCache,
              poseA::Transformation,
              poseB::Transformation,
              max_iter=100,
              atol=1e-6)
    rotAinv = transform_deriv(inv(poseA), 0)
    rotBinv = transform_deriv(inv(poseB), 0)
    simplex = transform_simplex(cache, poseA, poseB)
    iter = 1

    while true
        weights = projection_weights(simplex)
        min_weight, index_to_replace = findmin(weights)
        best_point = linear_combination(weights, simplex)
        if min_weight > 0
            # Nominally in collision; but check for numerical issues
            separation_squared = best_point ⋅ best_point
            @assert separation_squared < sqrt(1_000_000 * eps(typeof(separation_squared)))
            closest_point_in_body = linear_combination(weights, cache.simplex_points)
            return GJKResult(simplex, weights, true, closest_point_in_body, iter)
        end

        direction = -best_point
        direction_in_A = rotAinv * direction
        direction_in_B = rotBinv * direction

        starting_vertex_index = 1
        starting_vertex = cache.simplex_points[starting_vertex_index]
        starting_vertex_score =
            dot(value(starting_vertex.a), direction_in_A) +
            dot(value(starting_vertex.b), direction_in_B)
        for j in 2:length(cache.simplex_points)
            candidate = cache.simplex_points[j]
            candidate_score =
                dot(value(candidate.a), direction_in_A) +
                dot(value(candidate.b), direction_in_B)
            if candidate_score > starting_vertex_score
                starting_vertex_score = candidate_score
                starting_vertex = candidate
            end
        end

        improved_vertex = Difference(
            support_vector_max(cache.bodyA, direction_in_A, starting_vertex.a),
            support_vector_max(cache.bodyB, -direction_in_B, starting_vertex.b))
        improved_point = poseA(value(improved_vertex.a)) - poseB(value(improved_vertex.b))
        score = dot(improved_point, direction)
        if score <= dot(best_point, direction) + atol || iter >= max_iter
            closest_point_in_body = linear_combination(weights, cache.simplex_points)
            return GJKResult(simplex, weights, false, closest_point_in_body, iter)
        else
            cache.simplex_points[index_to_replace] = improved_vertex
            simplex = setindex(simplex, improved_point, index_to_replace)
            # simplex[index_to_replace] = improved_point
        end
        iter += 1
    end
end

function penetration_distance(simplex)
    best_dist_squared = nothing
    for i in eachindex(simplex)
        face = simplex_face(simplex, i)
        weights = projection_weights(face)
        closest_point = linear_combination(weights, face)
        dist_squared = closest_point ⋅ closest_point
        if best_dist_squared === nothing || dist_squared < best_dist_squared
            best_dist_squared = dist_squared
        end
    end
    return sqrt(best_dist_squared)
end

function gjk(geomA, geomB,
             poseA::Transformation=IdentityTransformation(),
             poseB::Transformation=IdentityTransformation())
    cache = CollisionCache(geomA, geomB)
    gjk!(cache, poseA, poseB)
end


text:
F, z = CyclotomicField(2)
@test Hecke.force_coerce_cyclo(F, z) == z


text:
struct PlotDisplay <: AbstractDisplay end

abstract type AbstractBackend end
function backend_display end

"""
Currently available displays by backend
"""
const available_backends = AbstractBackend[]
const current_backend = Ref{Union{Missing,AbstractBackend}}(missing)
const use_display = Ref{Bool}(true)

function inline!(inline = true)
    use_display[] = !inline
end

function register_backend!(backend::AbstractBackend)
    push!(available_backends, backend)
    if(length(available_backends) == 1)
        current_backend[] = backend
    end
    nothing
end

function __init__()
    pushdisplay(PlotDisplay())
end

function Base.display(d::PlotDisplay, scene::Scene)
    use_display[] || throw(MethodError(display, (d, scene)))
    try
        backend_display(current_backend[], scene)
    catch ex
        if ex isa MethodError && ex.f in (backend_display, backend_show)
            throw(MethodError(display, (d, scene)))
        else
            rethrow()
        end
    end
end

Base.showable(mime::MIME, scene::Scene) = backend_showable(current_backend[], mime, scene)

# have to be explicit with mimetypes to avoid ambiguity

function backend_show end
for M in (MIME"text/plain", MIME)
    @eval function Base.show(io::IO, m::$M, scene::Scene)
        res = get(io, :juno_plotsize, size(scene))
        resize!(scene, res...)
        update!(scene)
        AbstractPlotting.backend_show(current_backend[], io, m, scene)
    end
end

function backend_showable(backend, m::MIME, scene::Scene)
    hasmethod(backend_show, Tuple{typeof(backend), IO, typeof(m), typeof(scene)})
end

# fallback show when no backend is selected
function backend_show(backend, io::IO, ::MIME"text/plain", scene::Scene)
    println(io, "Scene ($(size(scene, 1))px, $(size(scene, 2))px):")
    println(io, "events:")
    for field in fieldnames(Events)
        println(io, "    ", field, ": ", to_value(getfield(scene.events, field)))
    end
    println(io, "plots:")
    for plot in scene.plots
        println(io, "   *", typeof(plot))
    end
    print(io, "subscenes:")
    for subscene in scene.children
        print(io, "\n   *scene($(size(subscene, 1))px, $(size(subscene, 2))px)")
    end
    return
end

function backend_show(backend, io::IO, ::MIME"text/plain", plot::Combined)
    println(io, typeof(plot))
    println(io, "plots:")
    for p in plot.plots
        println(io, "   *", typeof(p))
    end
    print(io, "attributes:")
    for (k, v) in theme(plot)
        print(io, "\n  $k : $(typeof(v))")
    end
end

function Base.show(io::IO, ::MIME"text/plain", plot::Atomic)
    println(io, typeof(plot))
    print(io, "attributes:")
    for (k, v) in theme(plot)
        print(io, "\n  $k : $(typeof(to_value(v)))")
    end
end

# Stepper for generating progressive plot examples
mutable struct Stepper
    scene::Scene
    folder::String
    step::Int
end

function Stepper(scene, path)
    ispath(path) || mkpath(path)
    Stepper(scene, path, 1)
end

function FileIO.save(filename::String, scene::Scene)
    open(filename, "w") do s
        show(IOContext(s, :full_fidelity => true), MIME"image/png"(), scene)
    end
end
"""
    step!(s::Stepper)
steps through a `Makie.Stepper` and outputs a file with filename `filename-step.jpg`.
This is useful for generating progressive plot examples.
"""
function step!(s::Stepper)
    FileIO.save(joinpath(s.folder, basename(s.folder) * "-$(s.step).jpg"), s.scene)
    s.step += 1
    return s
end

"""
Record all window events that happen while executing function `f`
for `scene` and serializes them to `path`.
"""
function record_events(f, scene::Scene, path::String)
    display(scene)
    result = Vector{Pair{Float64, Pair{Symbol, Any}}}()
    for field in fieldnames(Events)
        on(getfield(scene.events, field)) do value
            value = isa(value, Set) ? copy(value) : value
            push!(result, time() => (field => value))
        end
    end
    f()
    open(path, "w") do io
        serialize(io, result)
    end
end

"""
Replays the serialized events recorded with `record_events` in `path` in `scene`.
"""
replay_events(scene::Scene, path::String) = replay_events(()-> nothing, scene, path)
function replay_events(f, scene::Scene, path::String)
    display(scene)
    events = open(io-> deserialize(io), path)
    sort!(events, by = first)
    for i in 1:length(events)
        t1, (field, value) = events[i]
        field == :mousedrag && continue
        if field == :mousebuttons
            Base.invokelatest() do
                getfield(scene.events, field)[] = value
            end
        else
            Base.invokelatest() do
                getfield(scene.events, field)[] = value
            end
        end
        f()
        yield()
        if i < length(events)
            t2, (field, value) = events[i + 1]
            # min sleep time 0.001
            (t2 - t1 > 0.001) && sleep(t2 - t1)
        end
    end
end

struct RecordEvents
    scene::Scene
    path::String
end

function Base.display(d::PlotDisplay, re::RecordEvents)
    display(d, re.scene)
end

export Stepper, step!, replay_events, record_events, RecordEvents


text:
# Autogenerated wrapper script for BerkeleyDB_jll for aarch64-linux-gnu
export db_archive, db_checkpoint, db_convert, db_deadlock, db_dump, db_hotbackup, db_load, db_log_verify, db_printlog, db_recover, db_replicate, db_stat, db_tuner, db_upgrade, db_verify, libdb

using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("BerkeleyDB")
JLLWrappers.@declare_executable_product(db_archive)
JLLWrappers.@declare_executable_product(db_checkpoint)
JLLWrappers.@declare_executable_product(db_convert)
JLLWrappers.@declare_executable_product(db_deadlock)
JLLWrappers.@declare_executable_product(db_dump)
JLLWrappers.@declare_executable_product(db_hotbackup)
JLLWrappers.@declare_executable_product(db_load)
JLLWrappers.@declare_executable_product(db_log_verify)
JLLWrappers.@declare_executable_product(db_printlog)
JLLWrappers.@declare_executable_product(db_recover)
JLLWrappers.@declare_executable_product(db_replicate)
JLLWrappers.@declare_executable_product(db_stat)
JLLWrappers.@declare_executable_product(db_tuner)
JLLWrappers.@declare_executable_product(db_upgrade)
JLLWrappers.@declare_executable_product(db_verify)
JLLWrappers.@declare_library_product(libdb, "libdb-18.1.so")
function __init__()
    JLLWrappers.@generate_init_header(OpenSSL_jll)
    JLLWrappers.@init_executable_product(
        db_archive,
        "bin/db_archive",
    )

    JLLWrappers.@init_executable_product(
        db_checkpoint,
        "bin/db_checkpoint",
    )

    JLLWrappers.@init_executable_product(
        db_convert,
        "bin/db_convert",
    )

    JLLWrappers.@init_executable_product(
        db_deadlock,
        "bin/db_deadlock",
    )

    JLLWrappers.@init_executable_product(
        db_dump,
        "bin/db_dump",
    )

    JLLWrappers.@init_executable_product(
        db_hotbackup,
        "bin/db_hotbackup",
    )

    JLLWrappers.@init_executable_product(
        db_load,
        "bin/db_load",
    )

    JLLWrappers.@init_executable_product(
        db_log_verify,
        "bin/db_log_verify",
    )

    JLLWrappers.@init_executable_product(
        db_printlog,
        "bin/db_printlog",
    )

    JLLWrappers.@init_executable_product(
        db_recover,
        "bin/db_recover",
    )

    JLLWrappers.@init_executable_product(
        db_replicate,
        "bin/db_replicate",
    )

    JLLWrappers.@init_executable_product(
        db_stat,
        "bin/db_stat",
    )

    JLLWrappers.@init_executable_product(
        db_tuner,
        "bin/db_tuner",
    )

    JLLWrappers.@init_executable_product(
        db_upgrade,
        "bin/db_upgrade",
    )

    JLLWrappers.@init_executable_product(
        db_verify,
        "bin/db_verify",
    )

    JLLWrappers.@init_library_product(
        libdb,
        "lib/libdb.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
const jlAState{T} = JLArray{T,1}
const jlAStateBatch{T} = JLArray{T,2}
const jlAStateBatchVec{T} = JLArray{T,3}

function build_rng_generator_T(arrT::GPUArray, seed)
    c = similar(arrT)
    fill!(c, seed)
    return GPUArrays.global_rng(c)
end

@inline state_uview(σ::gpuAStateBatch, i)    = view(σ, :, i)
@inline state_uview(σ::gpuAStateBatchVec, i) = view(σ, :, :, i)
@inline state_uview(σ::gpuAStateBatchVec, batch, el) = view(σ, :, batch, el)


text:
#import COMMON\common everywhere you want to import/"using" modules, as well as using globals
import Base.showerror

export
    printfields, capitalize, TooManyIterationsException

macro inherit(name, base, fields)
    base_type = Core.eval(@__MODULE__, base)
    base_fieldnames = fieldnames(base_type)
    base_types = [t for t in base_type.types]
    base_fields = [:($f::$T) for (f, T) in zip(base_fieldnames, base_types)]
    res = :(mutable struct $name end)
    push!(res.args[end].args, base_fields...)
    push!(res.args[end].args, fields.args...)
    return res
end

"""
    @printfields
Prints names and contents of `object`'s fields.
"""
macro printfields(object)
    return :(for name in fieldnames(typeof($object))
        println(name," ",typeof(getfield($object,name))," ", getfield($object,name))
    end)
end

"Capitalizes a word"
capitalize(str::AbstractString) = begin
    return string(uppercase(str[1]), lowercase(str[2:end]))
end

"""
    TooManyIterationsException
Exception for a routine with a limited number of Iterations.
"""
struct TooManyIterationsException
    max_iterations::Int64
end
Base.showerror(io::IO, e::TooManyIterationsException) = print(io, "rached maximum number ($(e.max_iterations)) of iterations")

#=
common:
- Julia version:
- Author: Dan
- Date: 2021-06-02
=#


text:
function transformcoords(a,b,c,A,B,C,X)

R = [a[:] b[:] c[:]]*inv([A[:] B[:] C[:]]);

d = size(X,1)
n = size(X,2)
x = R*reshape(X, (d, n));

return x, R

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
export PolarizedModel, coherencymatrix

import ComradeBase: AbstractPolarizedModel, m̆, evpa, CoherencyMatrix, StokesVector

"""
    $(TYPEDEF)
Wrapped model for a polarized model. This uses the stokes representation of the image.
"""
struct PolarizedModel{TI,TQ,TU,TV} <: AbstractPolarizedModel
    I::TI
    Q::TQ
    U::TU
    V::TV
end

@inline visanalytic(::Type{PolarizedModel{I,Q,U,V}}) where {I,Q,U,V} = visanalytic(I)*visanalytic(Q)*visanalytic(U)*visanalytic(V)
@inline imanalytic(::Type{PolarizedModel{I,Q,U,V}}) where {I,Q,U,V} = imanalytic(I)*imanalytic(Q)*imanalytic(U)*imanalytic(V)

@inline function intensity_point(pmodel::PolarizedModel, u, v)
    I = intensity_point(pmodel.I, u, v)
    Q = intensity_point(pmodel.Q, u, v)
    U = intensity_point(pmodel.U, u, v)
    V = intensity_point(pmodel.V, u, v)
    return StokesVector(I,Q,U,V)
end

"""
    $(SIGNATURES)
Computes the visibility in the stokes basis of the polarized model
"""
@inline function visibility(pimg::PolarizedModel, u, v)
    si = visibility(pimg.I, u, v)
    sq = visibility(pimg.Q, u, v)
    su = visibility(pimg.U, u, v)
    sv = visibility(pimg.V, u, v)
    return StokesVector(si, sq, su, sv)
end

"""
    $(SIGNATURES)
Finds the polarized intensity map of the polarized model pmodel.
"""
function intensitymap!(pimg::IntensityMap{<:StokesVector}, pmodel::PolarizedModel)
    intensitymap!(stokes(pimg, :I), pmodel.I)
    intensitymap!(stokes(pimg, :Q), pmodel.Q)
    intensitymap!(stokes(pimg, :U), pmodel.U)
    intensitymap!(stokes(pimg, :V), pmodel.V)
end

function intensitymap(pmodel::PolarizedModel, fovx::Real, fovy::Real, nx::Int, ny::Int; pulse=DeltaPulse())
    imgI = intensitymap(pmodel.I, fovx, fovy, nx, ny)
    imgQ = intensitymap(pmodel.Q, fovx, fovy, nx, ny)
    imgU = intensitymap(pmodel.U, fovx, fovy, nx, ny)
    imgV = intensitymap(pmodel.V, fovx, fovy, nx, ny)
    pimg = StructArray{StokesVector{eltype(imgI)}}((imgI.im, imgQ.im, imgU.im, imgV.im))
    return IntensityMap(pimg, fovx, fovy, pulse)
end

"""
    $(SIGNATURES)
Computes the coherency matrix of the polarized model.
"""
@inline function coherencymatrix(pimg::PolarizedModel, u, v)
    si = visibility(pimg, u, v)
    return convert(CoherencyMatrix, si)
end

"""
    $(SIGNATURES)
electric vector position angle or EVPA of the polarized model
"""
@inline function evpa(pimg, u, v)
    sq = visibility(pimg.Q, u, v)
    su = visibility(pimg.U, u, v)
    return 1/2*angle(su/sq)
end

"""
    $(SIGNATURES)
Computes the fractional linear polarization in the visibility domain
```math
    \\breve{m} = \\frac{Q + iU}{I}
```
"""
@inline function m̆(pimg, u, v)
    Q = visibility(pimg.Q, u, v)
    U = visibility(pimg.U, u, v)
    I = visibility(pimg.I, u, v)
    return (Q+1im*U)/I
end


text:
module PlotsTests

include("add_packages.jl")
include("imgcomp.jl")

# don't actually show the plots
Random.seed!(1234)
default(show=false, reuse=true)
img_tol = isinteractive() ? 1e-2 : 10e-2

@testset "GR" begin
    ENV["PLOTS_TEST"] = "true"
    ENV["GKSwstype"] = "100"
    @test gr() == Plots.GRBackend()
    @test backend() == Plots.GRBackend()

    @static if Sys.islinux()
        image_comparison_facts(:gr, tol=img_tol, skip = [25, 30])
    end
end

# @static if isinteractive()
#     @testset "PyPlot" begin
#         @test pyplot() == Plots.PyPlotBackend()
#         @test backend() == Plots.PyPlotBackend()
#
#         image_comparison_facts(:pyplot, tol=img_tol, skip = [2, 25, 30, 31])
#     end
# end

@testset "UnicodePlots" begin
    @test unicodeplots() == Plots.UnicodePlotsBackend()
    @test backend() == Plots.UnicodePlotsBackend()

    # lets just make sure it runs without error
    @test isa(plot(rand(10)), Plots.Plot) == true
end

# The plotlyjs testimages return a connection error on travis:
# connect: connection refused (ECONNREFUSED)

# @static if isinteractive()
#     @testset "PlotlyJS" begin
#         @test plotlyjs() == Plots.PlotlyJSBackend()
#         @test backend() == Plots.PlotlyJSBackend()
#
#         image_comparison_facts(:plotlyjs,
#             skip=[
#                 2,  # animation (skipped for speed)
#                 25,
#                 27, # (polar plots) takes very long / not working
#                 30,
#                 31, # animation (skipped for speed)
#             ],
#             tol=img_tol)
#     end
# end

# InspectDR returns that error on travis:
# ERROR: LoadError: InitError: Cannot open display:
#  in Gtk.GLib.GError(::Gtk.##229#230) at /home/travis/.julia/v0.5/Gtk/src/GLib/gerror.jl:17

# @testset "InspectDR" begin
#     @test inspectdr() == Plots.InspectDRBackend()
#     @test backend() == Plots.InspectDRBackend()
#
#     image_comparison_facts(:inspectdr,
#         skip=[
#             2,  # animation
#             6,  # heatmap not defined
#             10, # heatmap not defined
#             22, # contour not defined
#             23, # pie not defined
#             27, # polar plot not working
#             28, # heatmap not defined
#             31, # animation
#         ],
#         tol=img_tol)
# end

# @testset "Plotly" begin
#     @test plotly() == Plots.PlotlyBackend()
#     @test backend() == Plots.PlotlyBackend()
#
#     # # until png generation is reliable on OSX, just test on linux
#     # @static Sys.islinux() && image_comparison_facts(:plotly, only=[1,3,4,7,8,9,10,11,12,14,15,20,22,23,27], tol=img_tol)
# end

# @testset "Immerse" begin
#     @test immerse() == Plots.ImmerseBackend()
#     @test backend() == Plots.ImmerseBackend()
#
#     # as long as we can plot anything without error, it should be the same as Gadfly
#     image_comparison_facts(:immerse, only=[1], tol=img_tol)
# end

# @testset "PlotlyJS" begin
#     @test plotlyjs() == Plots.PlotlyJSBackend()
#     @test backend() == Plots.PlotlyJSBackend()
#
#     # as long as we can plot anything without error, it should be the same as Plotly
#     image_comparison_facts(:plotlyjs, only=[1], tol=img_tol)
# end

# @testset "Gadfly" begin
#     @test gadfly() == Plots.GadflyBackend()
#     @test backend() == Plots.GadflyBackend()
#
#     @test typeof(plot(1:10)) == Plots.Plot{Plots.GadflyBackend}
#     @test plot(Int[1,2,3], rand(3)) == not(nothing)
#     @test plot(sort(rand(10)), rand(Int, 10, 3)) == not(nothing)
#     @test plot!(rand(10,3), rand(10,3)) == not(nothing)
#
#     image_comparison_facts(:gadfly, skip=[4,6,23,24,27], tol=img_tol)
# end

@testset "Axes" begin
    p = plot()
    axis = p.subplots[1][:xaxis]
    @test typeof(axis) == Plots.Axis
    @test Plots.discrete_value!(axis, "HI") == (0.5, 1)
    @test Plots.discrete_value!(axis, :yo) == (1.5, 2)
    @test Plots.ignorenan_extrema(axis) == (0.5,1.5)
    @test axis[:discrete_map] == Dict{Any,Any}(:yo  => 2, "HI" => 1)

    Plots.discrete_value!(axis, ["x$i" for i=1:5])
    Plots.discrete_value!(axis, ["x$i" for i=0:2])
    @test Plots.ignorenan_extrema(axis) == (0.5, 7.5)
end

@testset "NoFail" begin
    histogram([1, 0, 0, 0, 0, 0])
end

# tests for preprocessing recipes

# @testset "recipes" begin

    # user recipe

    # type T end
    # @recipe function f(::T)
    #     line := (3,0.3,:red)
    #     marker := (20,0.5,:blue,:o)
    #     bg := :yellow
    #     rand(10)
    # end
    # plot(T())

    # plot recipe

    # @recipe function f(::Type{Val{:hiplt}},plt::Plot)
    #     line := (3,0.3,:red)
    #     marker := (20,0.5,:blue,:o)
    #     t := :path
    #     bg:=:green
    #     ()
    # end
    # plot(rand(10),t=:hiplt)

    # series recipe

    # @recipe function f(::Type{Val{:hi}},x,y,z)
    #     line := (3,0.3,:red)
    #     marker := (20,0.5,:blue,:o)
    #     t := :path
    #     ()
    # end
    # plot(rand(10),t=:hiplt)

# end

end # module


text:
"The Omega Lanaguage (Not Library)"
module Lang

include("oexpr.jl")       # Abstract Syntax tree for Omega
# include("parse.jl")     # Parse Omega from tex
# include("repl.jl")      # Parse Omega from tex
include("compile.jl")     # Compile OmegaLang Expression into shit
include("unparse.jl")     # Convert expression tree into pretty printed syntax
include("expr.jl")        # Convert from Julia Expression to Omegea Expression

end


text:
module MambaModels

using Reexport

@reexport using StatisticalRethinking, Distributions, RDatasets, DataFrames
@reexport using StatsBase, StatsPlots, StatsFuns, Plots
@reexport using CSV, DelimitedFiles, Serialization
@reexport using MCMCChains, Mamba

using DataStructures

const src_path_m = @__DIR__

"""

# rel_path_m

Relative path using the MambaModels src/ directory.

### Example to get access to the data subdirectory
```julia
rel_path_m("..", "data")
```
"""
rel_path_m(parts...) = normpath(joinpath(src_path_m, parts...))

include("scriptentry_m.jl")
include("generate_m.jl")

export
  rel_path_m,
  script_dict_m,
  generate_m

end # module