text:
function h5e_get_auto(estack_id, func, client_data)
    ccall((:H5Eget_auto2, libhdf5), Herr, (Hid, Ptr{Ptr{Nothing}}, Ptr{Ptr{Nothing}}), estack_id, func, client_data)
end

function h5e_set_auto(estack_id, func, client_data)
    ccall((:H5Eset_auto2, libhdf5), Herr, (Hid, Ptr{Nothing}, Ptr{Nothing}), estack_id, func, client_data)
end

function e_get_auto()
    errfunc = Ref(Ptr{Nothing}())
    errdata = Ref(Ptr{Nothing}())
    h5e_get_auto(HDF5.H5E_DEFAULT, errfunc, errdata)
    return errfunc[], errdata[]
end

function e_set_auto(errfunc, errdata)
    h5e_set_auto(HDF5.H5E_DEFAULT, errfunc, errdata)
end

const _errfunc = Ref(Ptr{Nothing}())
const _errdata = Ref(Ptr{Nothing}())

function disable_dag()
    if _errfunc[] == C_NULL
        _errfunc[], _errdata[] = e_get_auto()
    end
    e_set_auto(C_NULL, C_NULL)
end

function enable_dag()
    e_set_auto(_errfunc[], _errdata[])
end

ext: jl


text:
# ---
# title: 977. Squares of a Sorted Array
# id: problem977
# author: zhwang
# date: 2022-01-19
# difficulty: Easy
# categories: Array, Two Pointers
# link: <https://leetcode.com/problems/squares-of-a-sorted-array/description/>
# hidden: true
# ---
#
# Given an integer array `nums` sorted in **non-decreasing** order, return _an
# array of **the squares of each number**  sorted in non-decreasing order_.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [-4,-1,0,3,10]
#     Output: [0,1,9,16,100]
#     Explanation: After squaring, the array becomes [16,1,0,9,100].
#     After sorting, it becomes [0,1,9,16,100].
#
#
# **Example 2:**
#
#
#
#     Input: nums = [-7,-3,2,3,11]
#     Output: [4,9,9,49,121]
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 104`
#   * `-104 <= nums[i] <= 104`
#   * `nums` is sorted in **non-decreasing** order.
#
#
## @lc code=start
using LeetCode

## binary search + merge sort
function squares_of_a_sorted_array(nums::Vector{Int})::Vector{Int}
    n = length(nums)
    ## binary search
    left, right = 1, n
    while left <= right
        mid = (right + left) >> 1
        if nums[mid] < 0
            left = mid + 1
        else
            right = mid - 1
        end
    end
    ## merge sort
    i, j, res = right, left, Int[]
    while i >= 1 && j <= n
        if nums[i]^2 < nums[j]^2
            push!(res, nums[i]^2)
            i -= 1
        else
            push!(res, nums[j]^2)
            j += 1
        end
    end
    return vcat(res, (i == 0) ? nums[j:end] .^ 2 : nums[i:-1:1] .^ 2)
end

## @lc code=end


text:
@testset "proc" begin
    c = lxc(raw"\foo{bar baz}")
    @test P.proc(c) == "bar baz"
end

@testset "html" begin
    @test P.html("aaa") == "\n~~~aaa~~~\n"
end

@testset "lxargs" begin
    s = """
        :section, 1, 3, title="hello", foo="bar", a=5
        """
    a, ka = P.lxargs(s)
    @test all(a .== (:section, 1, 3))
    @test all(ka .== (
        :title => "hello",
        :foo => "bar",
        :a => 5))
    s = """
        title=5, foo="bar"
        """
    a, ka = P.lxargs(s)
    @test isempty(a)
    @test all(ka .== (:title => 5, :foo => "bar"))
    s = """
        5, 3
        """
    a, ka = P.lxargs(s)
    @test isempty(ka)
    @test all(a .== (5, 3))
    # bad ordering
    s = """
        5, a=3, 2
        """
    @test_throws ErrorException P.lxargs(s)
end

@testset "invalid" begin
    c = lxc("""\\table{a:3}""")
    @test_throws Exception lx_table(c,0)
    c = lxc("""\\table{a=5}""")
    @test_throws Exception lx_table(c,0)
    c = lxc("""\\table{\"\"b}""")
    @test_throws Exception lx_table(c,0)
end


text:
module test_bukdu_resources

using Test
using Bukdu # ApplicationController Conn Router Utils routes resources
import .Bukdu.Actions: index, show, new, edit, create, delete, update

struct ArticleController <: ApplicationController
    conn::Conn
end

index(::ArticleController)   = (:index,)
show(c::ArticleController)   = (:show, c.params.id)
new(::ArticleController)     = (:new,)
edit(c::ArticleController)   = (:edit, c.params.id)
create(::ArticleController)  = (:create,)
delete(c::ArticleController) = (:delete, c.params.id)
update(c::ArticleController) = (:update, c.params.id)

routes() do
    resources("/articles", ArticleController, :id=>Int)
end

@test Router.call(get, "/articles").got   == (:index,)
@test Router.call(get, "/articles/2").got == (:show, 2)

@test Utils.read_stdout(CLI.routes) == """
GET     /articles           ArticleController  index
GET     /articles/new       ArticleController  new
GET     /articles/:id/edit  ArticleController  edit
GET     /articles/:id       ArticleController  show
POST    /articles           ArticleController  create
DELETE  /articles/:id       ArticleController  delete
PATCH   /articles/:id       ArticleController  update
PUT     /articles/:id       ArticleController  update"""

Router.call(delete, "/articles/1")
Router.call(patch, "/articles/1")
Router.call(put, "/articles/1")

Routing.reset!()

routes() do
    resources("/articles", ArticleController, only=[index, show])
end
@test Utils.read_stdout(CLI.routes) == """
GET  /articles      ArticleController  index
GET  /articles/:id  ArticleController  show"""
Routing.reset!()

routes() do
    resources("/articles", ArticleController, :id=>Int, except=[new, edit, create, delete, update])
end
@test Utils.read_stdout(CLI.routes) == """
GET  /articles      ArticleController  index
GET  /articles/:id  ArticleController  show"""

result = Router.call(Bukdu.head, "/articles/1")
@test result.resp.status == 200
@test result.got == (:show, nothing)
@test result.route.action === show
@test result.route.param_types == Dict(:id=>Int)
@test result.route.path_params == ["id"=>"1"]

Routing.reset!()

end # module test_bukdu_resources


text:
function eval_init(nlp::Class_CUTEst, pars::Class_parameters, timer::class_advanced_timer, x)
    start_advanced_timer(timer, "INIT/evals")

    a = eval_a(nlp, x);
    if a == Float64[]
      a = zeros(1)
    end
    J = eval_jac(nlp, x)
    g = eval_grad_f(nlp, x)

    s = deepcopy(a);
    m = length(a)

    if length(SparseArrays.nonzeros(J)) > 0
      if(isbad(SparseArrays.nonzeros(J)))
          throw(Eval_NaN_error(getbad(SparseArrays.nonzeros(J)),x,"J"))
      end

      if(isbad(a))
          throw(Eval_NaN_error(getbad(a),x,"a"))
      end

      if isbad(g)
        throw(Eval_NaN_error(getbad(g),x,"g"))
      end
    end
    pause_advanced_timer(timer, "INIT/evals")

    return a, J, g
end


text:
function large03()
#
#  test polynomial suggested by Z Zeng
#
    p, z = large02()
    p = conv(p, p)
    z.mult .*= 2
    p, z
end


text:
export regOptions

type regOptions

  # overall
  parametricOnly::Bool
  matrixFree::Bool
  levels::Array{Int,1}
  interpolateToReferenceImage::Bool

  # distance
  # ngf parameter
  edgeParameterR::Float64
  edgeParameterT::Float64
  useEdgeParameterInNumerator::Bool

  # regularizer
  regularizerWeight::Float64

  # optimization
  maxIterGaussNewton::Int
  maxIterCG::Int
  tolCG::Float64
  stopping::Dict

  # additional options in dict
  additionalOptions::Dict

  # set default options
  function regOptions()
    parametricOnly = false
    matrixFree = false
    levels = [5,4,3]
    interpolateToReferenceImage = false
    edgeParameterR = 0.01
    edgeParameterT = 0.01
    useEdgeParameterInNumerator = true
    regularizerWeight = 1
    maxIterGaussNewton = 30
    maxIterCG = 2000
    tolCG = 1e-5
    stopping = Dict(
        "tolJ" => 1e-3,    # tolerance: change of the objective function
        "tolY" => 1e-2,    # tolerance: change of the variables
        "tolG" => 1e-2,    # tolerance: change of the gradient
        "tolQ" => 1e-4     # tolerance: change of quotient)
        )
    additionalOptions = Dict()
    new(parametricOnly,matrixFree,levels, interpolateToReferenceImage,
        edgeParameterR,edgeParameterT,useEdgeParameterInNumerator,
        regularizerWeight,
        maxIterGaussNewton,maxIterCG, tolCG, stopping,
        additionalOptions)
  end

end


text:
using MarketData

facts("Levels") do

    context("floor pivots") do
        @fact floorpivots(ohlc)["r3"].values[1]    --> roughly(123.310, atol=.01)              # values verified by various website calculators
        @fact floorpivots(ohlc)["r2"].values[1]    --> roughly(119.52, atol=.01)
        @fact floorpivots(ohlc)["r1"].values[1]    --> roughly(115.73, atol=.01)
        @fact floorpivots(ohlc)["pivot"].values[1] --> roughly(108.71, atol=.01)
        @fact floorpivots(ohlc)["s1"].values[1]    --> roughly(104.92, atol=.01)
        @fact floorpivots(ohlc)["s2"].values[1]    --> roughly(97.900, atol=.01)
        @fact floorpivots(ohlc)["s3"].values[1]    --> roughly(94.110, atol=.01)
        @fact floorpivots(ohlc).timestamp[end]     --> Date(2001,12,31)
    end

    context("woodiespivots") do
  #    #  @fact_approx_eq  97.37500000000001 value(wr4)[2]   # values NEED to be verified with online calculators
  #    #  @fact_approx_eq   88.62500000000001 value(ws4)[2]
        @fact woodiespivots(ohlc)["r3"].values[1]    --> roughly(124.465, atol=.01)
        @fact woodiespivots(ohlc)["r2"].values[1]    --> roughly(118.480, atol=.01)
        @fact woodiespivots(ohlc)["r1"].values[1]    --> roughly(113.655, atol=.01)
        @fact woodiespivots(ohlc)["pivot"].values[1] --> roughly(107.670, atol=.01)
        @fact woodiespivots(ohlc)["s1"].values[1]    --> roughly(102.845, atol=.01)
        @fact woodiespivots(ohlc)["s2"].values[1]    --> roughly(96.8625, atol=.01)
        @fact woodiespivots(ohlc)["s3"].values[1]    --> roughly(92.035, atol=.01)
        @fact woodiespivots(ohlc).timestamp[end]     --> Date(2001,12,31)
    end
end


text:
function filter(obj::DataFilter, arg0::NamedData)
    return jcall(obj, "filter", NamedData, (NamedData,), arg0)
end


text:
@testset "AnisotropicDiffusion" begin
    p, q = synthetic_gradient(SynthSphere(50))
    z = zeros(Float64, size(p))
    Z = AnisotropicDiffusion(z=z)(p, q)
    @test Z[38,38] ≈ -1.7507607679065658
    @test Z[1,1] ≈ -0.0005675736967002396
    @test Z[76,76] ≈ 37.33646747044026
end


text:
#=
   Simpson's paradox

   http://cplint.eu/example/inference/simpson.swinb
   """
   From "Pearl, Judea. Causality. Cambridge university press, 2009"

   Simpson's paradox ... refers to the phenomenon whereby an event C

   increases the probability of E in a given population p and, at the same time, decreases the probability
   of E in every subpopulation of p. In other words, if F and ¬F are two complementary properties describing two
   subpopulations, we might well encounter the inequalities
     P(E|C)>P(E|¬C)
     P(E|C,F)<P(E|¬C,F)
     P(E|C,¬F)<P(E|¬C,¬F)
   ...
  For example, if we associate C (connoting cause) with taking a certain drug, E (connoting effect) with
  recovery, and F with being a female, then ... the drug seems to be harmful to both males and females yet
  beneficial to the population as a whole.
  """

  cf ~/cplint/simplson.pl
     ~/webppl/simpson.wppl

=#

#=
  From the cplint model
/*
female:0.5.

recovery:0.6:- drug,\+ female.
recovery:0.7:- \+ drug,\+ female.

recovery:0.2:- drug,female.
recovery:0.3:- \+ drug,female.

drug:30/40:- \+ female.
drug:10/40:-female.

If we query for the conditional probabilities of recovery given treatment
on the whole population and on the two subpopulations, we get the results
in the tables above:
?- prob(recovery,drug,P).
P = 0.49999999999999994

?- prob(recovery,\+ drug,P).
P = 0.40000000000000013

?- prob(recovery,(drug,female),P).
P = 0.2

?- prob(recovery,(\+drug,female),P).
P = 0.3000000000000001

?- prob(recovery,(drug,\+female),P).
P = 0.6

?- prob(recovery,(\+ drug,\+female),P).
P = 0.7000000000000002

  This model give the same values:
  prob(recovery|drug): 0.5
  prob(recovery|\+ drug): 0.4000000000000001
  prob(recovery|(drug,female)): 0.20000000000000007
  prob(recovery|(\+drug,female)): 0.30000000000000004
  prob(recovery|(drug,\+female)): 0.6000000000000001
  prob(recovery|(\+ drug,\+female)): 0.7

*/

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function simpson(problem)
    female ~ flip(0.5)
    drug ~ female ? flip(10/40.0) :  flip(30/40.0)

    #=
    recovery ~ (drug == true && female==true) ? flip(0.2) :
               (drug == true && female==false) ? flip(0.6) :
               (drug == false && female==true) ? flip(0.3) :
               (drug == false && female==false) ? flip(0.7) : flip(0.0)
    =#
    # Alternative way of doing this...
    recovery ~  if     drug == true  && female==true   flip(0.2)
                elseif drug == true  && female==false  flip(0.6)
                elseif drug == false && female==true   flip(0.3)
                elseif drug == false && female==false  flip(0.7)
                else                                   flip(0.0)
                end

    # 1. prob(recovery|drug): 0.5
    # 2. prob(recovery|\+ drug): 0.4000000000000001
    # 3. prob(recovery|(drug,female)): 0.20000000000000007
    # 4. prob(recovery|(\+drug,female)): 0.30000000000000004
    # 5. prob(recovery|(drug,\+female)): 0.6000000000000001
    # 6. prob(recovery|(\+ drug,\+female)): 0.7

    if problem == 1
        true ~ Dirac(drug == true)
    elseif problem == 2
        true ~ Dirac(drug == false)
    elseif problem == 3
        true ~ Dirac(drug == true && female == true)
    elseif problem == 4
        true ~ Dirac(drug == false && female == true)
    elseif problem == 5
        true ~ Dirac(drug == true && female == false)
    else
        true ~ Dirac(drug == false && female == false)
    end

    return female, recovery, drug
end

function run_simpson(problem)

    model = simpson(problem)
    num_chains = 4

    # chains = sample(model, MH(), MCMCThreads(), 100_000, num_chains)
    chains = sample(model, MH(), MCMCThreads(), 10_000, num_chains)
    # chains = sample(model, MH(), 10_000)

    # chains = sample(model, PG(15), MCMCThreads(), 10_000, num_chains)

    # chains = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chains)
    # chains = sample(model, SMC(1000), 10_000)
    # chains = sample(model, IS(), 10_000)
    #
    display(chains)
    # display(plot(chains))

    # This the distributions of the different combinations of
    #   female, recovery, group
    genq = generated_quantities(model, chains)
    show_var_dist_pct(genq)
end

for p in 1:6
    println("\nproblem $p")
    run_simpson(p)

end

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
# Density utilities
# Author: Gabrio Rizzuti, rizzuti.gabrio@gatech.edu
# Date: August 2020

export sample_gaussmixture, negLogDensity_gaussmixture

function sample_gaussmixture(nsamples::Int64, mu::Array{Float32, 2}, sigma2::Array{Float32, 1})
    nmodes = length(sigma2)
    X = Array{Float32, 4}(undef, (1, 1, 2, nsamples))
    for i = 1:nsamples
        imode = rand(1:nmodes)
        X[1, 1, :, i] = mu[imode, :]+sqrt(sigma2[imode])*randn(Float32, 2)
    end
    return X
end

function negLogDensity_gaussmixture(X::Array{Float32, 4}, mu::Array{Float32, 2}, sigma2::Array{Float32, 1}; ε::Float32 = 1f-10)
    nmodes = length(sigma2)
    p = zeros(Float32, 1, 1, 1, size(X, 4))
    g = zeros(Float32, size(X))
    for i = 1:nmodes
        p_i = exp.(-((X[:, :, 1:1, :].-mu[i, 1]).^2f0+(X[:, :, 2:2, :].-mu[i, 2]).^2f0)/(2f0*sigma2[i]))/(nmodes*2f0*pi*sigma2[i])
        p += p_i
        g += -p_i.*(X.-reshape(mu[i, :], 1, 1, 2, 1))/sigma2[i]
    end
    return -log.(p.+ε), -g./(p.+ε)
end


text:
using jInv.Mesh
using DivSigGrad
using Multigrid

println("************************************************* Example 2D ******************************************************");

domain = [0.0, 1.0, 0.0, 1.0];
n = [100,100];
Mr = getRegularMesh(domain,n)
m  = exp(2*randn(prod(n)));
Ar = getDivSigGradMatrix(vec(m),Mr);
# Ar = Ar + 1e-4*norm(Ar,1)*speye(size(Ar,2));

levels      = 4;
numCores 	= 8;
maxIter     = 20;
relativeTol = 1e-6;
relaxType   = "SPAI";
relaxParam  = 1.0;
relaxPre 	= 2;
relaxPost   = 2;
cycleType   ='V';
coarseSolveType = "MUMPS";

MG = getMGparam(levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType);

N = size(Ar,2);

b = Ar*rand(N,20);
x = zeros(N,20);
SA_AMGsetup(Ar,MG,Float64,true,1,true);

# println("****************************** Stand-alone AMG single rhs: ******************************")
# solveMG(MG,vec(b[:,1]),vec(x[:,1]),true);

# x[:] = 0.0;
# println("****************************** Stand-alone AMG multiple rhs: ******************************")
# solveMG(MG,b,x,true);

# x[:] = 0.0;
# println("****************************** CG preconditioned with AMG single rhs: ******************************")
# solveCG_MG(Ar,MG,vec(b[:,1]),vec(x[:,1]),true)

# println("****************************** CG preconditioned with AMG multiple rhs: ******************************")
# solveCG_MG(Ar,MG,b,x,true)

println("****************************** GMRES preconditioned with AMG: single rhs: ******************************")
solveGMRES_MG(Ar,MG,vec(b[:,1]),vec(x[:,1]),true,10)

println("****************************** GMRES preconditioned with AMG: multiple rhs: ******************************")
solveGMRES_MG(Ar,MG,b,x,true,10)

println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~`")
error("Eran")
Ar = 0;
b = 0;
x = 0;
Mr = 0;
println("************************************************* Example 3D ******************************************************");

domain = [0.0, 1.0, 0.0, 1.0, 0.0, 1.0];
n = [100,100,100];
Mr = getRegularMesh(domain,n)
m  = exp(randn(prod(n)));
Ar = getDivSigGradMatrix(vec(m),Mr);
Ar = Ar + 1e-6*norm(Ar,1)*speye(size(Ar,2));

MG = getMGparam(levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType);

N = size(Ar,2);
b = Ar*rand(N,3);
x = zeros(N,3);

SA_AMGsetup(Ar,MG,Float64,true,size(b,2),true);
println("Stand-alone AMG:")
solveMG(MG,b,x,true);

println("CG preconditioned with AMG:")
x[:] = 0.0;
solveCG_MG(Ar,MG,b,x,true)


text:
@testset "Indices" begin
    inds = Indices{Int}((1,3,5))

    @test length(inds) == 3
    @test 3 ∈ inds
    @test 4 ∉ inds
    @test all(isodd, inds)
    @test keys(inds) === inds
end

@testset "VectorIndices" begin
    @test HashIndices() isa HashIndices{Any}

    inds = Indices{Int64}()

    @test isinsertable(inds)
    @test length(inds) == 0
    @test keys(inds) === inds
    @test isempty(inds)
    @test isequal(copy(inds), inds)
    @test_throws IndexError inds[10]
    @test length(unset!(inds, 10)) == 0
    io = IOBuffer(); print(io, inds); @test String(take!(io)) == "{}"
    io = IOBuffer(); show(io, MIME"text/plain"(), inds); @test String(take!(io)) == "0-element Indices{Int64,Array{Int64,1}}"
    @test_throws IndexError delete!(inds, 10)

    insert!(inds, 10)

    @test length(inds) == 1
    @test keys(inds) === inds
    @test !isempty(inds)
    @test isequal(copy(inds), inds)
    @test inds[10] == 10
    @test_throws IndexError insert!(inds, 10)
    @test length(set!(inds, 10)) == 1
    @test_throws IndexError insert!(inds, 10)
    io = IOBuffer(); print(io, inds); @test String(take!(io)) == "{10}"
    io = IOBuffer(); show(io, MIME"text/plain"(), inds); @test String(take!(io)) == "1-element Indices{Int64,Array{Int64,1}}\n 10"
    @test !isequal(inds, empty(inds))
    @test isequal(inds, copy(inds))
    @test isempty(empty(inds))

    delete!(inds, 10)

    @test isequal(inds, Indices{Int64}())

    for i = 2:2:100
        insert!(inds, i)
    end
    @test issetequal(inds, HashIndices(2:2:100))
    @test all(in(i, inds) == iseven(i) for i in 2:100)
    @test isempty(empty!(inds))
end


text:
type GEChoiceModel <: ChoiceModel
	genome::Vector{Int}
	activegene::Int
	function GEChoiceModel(genomelength::Int)
		genomelength > 0 || error("genome must have a length of 1 or more")
		new(fill(0.0, genomelength), 1)
	end
end

function GEChoiceModel(params::Vector)
	cm = GEChoiceModel(length(params))
	setparams!(cm, params)
	cm
end

function godelnumber(cm::GEChoiceModel, cc::ChoiceContext)

	# wrap to beginning of genome if necessary
	if cm.activegene > length(cm.genome)
		cm.activegene = 1
	end

	allele = cm.genome[cm.activegene]
	@assert allele >= 0

	# println("Choice Context: $(cc.lowerbound), $(cc.upperbound)")
	# println("Allele: $(allele)")

	gn = if (cc.datatype <: Bool) || (cc.datatype <: Integer ? (cc.lowerbound > typemin(cc.datatype)) : isfinite(cc.lowerbound))
			if (cc.datatype <: Bool) || (cc.datatype <: Integer ? (cc.upperbound < typemax(cc.datatype)) : isfinite(cc.upperbound))
				# println("both finite")
				cc.lowerbound + mod(allele, cc.upperbound - cc.lowerbound + ((cc.datatype <: Integer) ? 1 : 0.0))
			else
				# println("lower finite; upper infinite")
				cc.lowerbound + allele
			end
		else
			if (cc.datatype <: Integer ? (cc.upperbound < typemax(cc.datatype)) : isfinite(cc.upperbound))
				# println("lower infinite; upper finite")
				cc.upperbound - allele
			else
				# println("lower infinite; upper infinite")
				allele
			end
		end

	cptrace = Dict{Symbol, Any}(:all => allele, :gen => cm.activegene)

	cm.activegene += 1

	# println("GN: $(gn)")

	convert(cc.datatype, gn), cptrace

end

function resetstate!(cm::GEChoiceModel)
	cm.activegene = 1
end

getparams(cm::GEChoiceModel) = convert(Vector{Float64}, cm.genome)

function setparams!(cm::GEChoiceModel, params)
	length(params) == length(cm.genome) || error("Expected $(length(cm.genome)) parameters, but got $(length(params))")
	all(p -> 0<=p<=Float64(typemax(UInt32)), params) || error("All parameters must be between 0 and $(typemax(UInt32))")
	cm.genome = map(p->Int(round(p)), params)
end

paramranges(cm::GEChoiceModel) = fill((Float64(0),Float64(typemax(UInt32))), length(cm.genome))

show(io::IO, cm::GEChoiceModel) = print(io, "Grammmatical Evolution choice model (genome: ", cm.genome, ")")


text:
using Test
using NeXLCore

@testset "X-ray" begin

    e1, e2, e3 = element(20), element(84), element(92)

    @testset "Elements" begin

        @test z(e1) == 20
        @test z(e2) == 84
        @test z(e3) == 92

        @test isapprox(a(e1), 40.08, atol = 0.01)
        @test isapprox(a(e2), 209.00, atol = 0.01)
        @test isapprox(a(e3), 238.03, atol = 0.01)

        @test symbol(e1) == "Ca"
        @test symbol(e2) == "Po"
        @test symbol(e3) == "U"

        @test name(e1) == "Calcium"
        @test name(e2) == "Polonium"
        @test name(e3) == "Uranium"

        @test e1 ≠ e2
        @test e1 == element(20)
        @test e2 == element(84)
        @test e3 == element(92)

        @test e1 < e2
        @test e1 < e3
        @test e2 < e3

        @test e2 >= e1
        @test e2 >= e2
        @test e3 >= e2
        @test e3 >= e1

        @test e2 > e1
        @test e3 > e2
        @test e3 > e1

        @test n"Os" == elements[76]
        @test n"Osmium" == elements[76]
        @test n"76" == elements[76]
    end

    @testset "Shells" begin
        @test length(ksubshells) == 1
        @test length(lsubshells) == 3
        @test length(msubshells) == 5
        @test length(nsubshells) == 7
        @test length(osubshells) == 9

        @test ksubshells == map(n->parse(SubShell,n), ("K",))
        @test lsubshells == map(n->parse(SubShell,n), ("L1", "L2", "L3"))
        @test msubshells == map(n->parse(SubShell,n), ("M1", "M2", "M3", "M4", "M5"))
        @test nsubshells == map(n->parse(SubShell,n),("N1", "N2", "N3", "N4", "N5", "N6", "N7"))
        @test osubshells ==
              map(n->parse(SubShell,n), ("O1", "O2", "O3", "O4", "O5", "O6", "O7", "O8", "O9"))

        @test shell(parse(SubShell,"K")) == Shell(1)
        @test shell(parse(SubShell,"L2")) == Shell(2)
        @test shell(parse(SubShell,"M4")) == Shell(3)
        @test shell(parse(SubShell,"N7")) == Shell(4)
        @test shell(parse(SubShell,"O1")) == Shell(5)

        @test capacity(n"K1") == 2
        @test capacity(n"L2") == 2
        @test capacity(n"M5") == 6
        @test capacity(n"N6") == 6
        @test capacity(n"O9") == 10

        as1, as2, as3 = AtomicSubShell(e1.number, subshell("L3")),
        AtomicSubShell(e2.number, subshell("M5")),
        AtomicSubShell(e3.number, subshell("K"))

        @test repr(as1) == "Ca L3"
        @test repr(as2) == "Po M5"
        @test repr(as3) == "U K"

        @test as1 == n"Ca L3"
        @test as2 == n"Po M5"
        @test as3 == n"U K"

        @test shell(as1) == Shell(2)
        @test shell(as2) == Shell(3)
        @test shell(as3) == Shell(1)

        @test n"Ca K" > n"Ca L1"
        @test n"Ca K" < n"Ti K"
        @test as1 < as2
        @test !(as2 < as1)
        @test !(as1 < as1)

        @test n"Ca K" == n"Ca K"
        @test !(n"Ca K" == n"Ca L1")

        @test isless(Shell(3), Shell(1))  # Binding energy order (not N)
        @test !isless(Shell(3), Shell(3))
        @test !isless(n"Ca K", n"Ca L1") # Binding energy order (not N)
        @test isless(n"Ca L3", n"Ca L1")
        @test isless(n"Ca K", n"Fe K")
        @test !isless(n"Fe K", n"Fe K")

        # Check that these funtions are defined for all available AtomicSubShell
        @test all(energy.(atomicsubshells(el)) ≠ 0.0 for el in eachelement())
        @test all(fluorescenceyield.(atomicsubshells(el), CullenEADL) ≠ 0.0 for el in eachelement())
        @test all(map(sh->ionizationcrosssection(sh,2.0*energy(sh)),atomicsubshells(el)) ≠ 0.0 for el in eachelement())
        @test all(map(sh->relativeionizationcrosssection(sh,2.0*energy(sh)),atomicsubshells(el)) ≠ 0.0 for el in eachelement())
        @test all(all(map(ass->parse(AtomicSubShell,repr(ass)) == ass, atomicsubshells(el))) for el in eachelement())
    end

    @testset "Transitions" begin
        @test n"K-L3" in alltransitions
        @test n"L3-M5" in alltransitions
        #@test_throws ArgumentError n"M5-L3"
        #@test_throws ArgumentError n"K-N7"

        @test shell(n"K-L3") == Shell(1)
        @test shell(n"L3-M5") == Shell(2)
        @test shell(n"M5-N7") == Shell(3)

        @test all(tr -> shell(tr) == Shell(1), ktransitions)
        @test all(tr -> shell(tr) == Shell(2), ltransitions)
        @test all(tr -> shell(tr) == Shell(3), mtransitions)
        @test all(tr -> shell(tr) == Shell(4), ntransitions)

        @test length(alltransitions) == sum(length.((
            ktransitions,
            ltransitions,
            mtransitions,
            ntransitions,
            otransitions,
        )))

        @test n"K-L3" == n"K-L3"
        @test n"K-L3" ≠ n"K-L2"

        @test transition(n"Fe L3-M5") == n"L3-M5"
        @test element(n"Fe L3-M5") == n"Fe"

        @test isapprox(energy(n"Fe", n"L3-M5"), 705.0, atol = 1.0)
        @test isapprox(energy(n"U", n"L3-M5"), 13614.6, atol = 1.0)
        @test energy(n"Fe", n"L3-M5") == energy(n"Fe L3-M5")
        @test isapprox(energy(n"C K-L2"), 277.4, atol = 1.0)

        # from https://www.physics.nist.gov/cgi-bin/XrayTrans/search.pl?element=Fe&trans=KL3&lower=&upper=&units=A
        @test isapprox(λ(n"Fe K-L3"), 1.936e-8, atol = 1e-10) #
        @test isapprox(λ(n"Cu L3-M5"), 13.336e-8, atol = 1e-10)

        x = n"Ba L3-M5"
        @test isapprox(NeXLCore.speedOfLight, λ(x) * ν(x), atol = 1.0e7) # C = λν
        @test isapprox(energy(x), NeXLCore.plancksConstant * ν(x), atol = 1.0) # E = hν
        @test isapprox(energy(x), NeXLCore.plancksConstant / (2π) * ω(x), atol = 1.0) # E = ħω

        x = n"U M5-N7"
        @test isapprox(NeXLCore.speedOfLight, λ(x) * ν(x), atol = 1.0e7) # C = λν
        @test isapprox(energy(x), NeXLCore.plancksConstant * ν(x), atol = 1.0) # E = hν
        @test isapprox(energy(x), NeXLCore.plancksConstant / (2π) * ω(x), atol = 1.0) # E = ħω

        @test isapprox(energy(n"Fe L3"), 708.0, atol = 1.0)
        @test energy(n"Fe L3") == energy(n"Fe", n"L3")

        @test has(n"C", n"K-L2")
        @test !has(n"Li", n"K-L3")
        @test !has(n"C", n"L3-M5")
        @test has(n"U", n"L3-M5")
        @test has(n"U", n"M5-N7")

        @test all(tr -> shell(tr) == Shell(2), characteristic(n"Fe", ltransitions))
        @test all(tr -> shell(tr) == Shell(1), characteristic(n"Fe", ktransitions))

        @test length(characteristic(n"Fe", ltransitions, 0.0)) == 14
        @test length(characteristic(n"Fe", ltransitions, 0.1)) == 6
        @test length(characteristic(n"Fe", ltransitions, 0.01)) == 9

        @test isless(n"Fe K-L3", n"Fe K-L2")
        @test !isless(n"Fe K-L3", n"Fe K-L3")
        @test !isless(n"Fe K-L2", n"Fe K-L3")

        # Check that these funtions are defined for all available transitions
        @test all(energy.(characteristic(el,alltransitions)) ≠ 0.0 for el in eachelement())
        @test all(weight.(NormalizeToUnity, characteristic(el,alltransitions)) ≠ 0.0 for el in eachelement())
        @test all(weight.(NormalizeByShell, characteristic(el,alltransitions)) ≠ 0.0 for el in eachelement())
        @test all(weight.(NormalizeBySubShell, characteristic(el,alltransitions)) ≠ 0.0 for el in eachelement())
        @test all(all(map(cxr->parse(CharXRay,repr(cxr)) == cxr, characteristic(el,alltransitions))) for el in eachelement())
    end
    @testset "Other" begin
        @test firstsubshell(Shell(3)) == n"M1"
        @test lastsubshell(Shell(3)) == n"M5"
        @test firstsubshell(Shell(1)) == n"K1"
        @test lastsubshell(Shell(1)) == n"K1"
        @test firstsubshell(Shell(2)) == n"L1"
        @test lastsubshell(Shell(2)) == n"L3"
        @test length(eachelement()) == 92
        @test isequal(eachelement()[1], n"H")
        @test isequal(eachelement()[end], n"U")
        @test eachelement()[26]==n"Fe"
        @test eachelement()[41]==n"Nb"

        @test all(ze -> ze[1]==z(ze[2]), zip(1:92, eachelement()))
        @test NeXLCore.fluorescenceyield(20, 1, 4, 9, CullenEADL) == 0.0
        @test isapprox(
            NeXLCore.fluorescenceyield(25, 1, 4, 9, CullenEADL),
            0.0007828,
            atol = 0.0000001,
        )
        @test configuration(n"Fe") == "1s² 2s² 2p⁶ 3s² 3p⁶ 4s² 3d⁶"
        @test density(n"Au") == 19.3
        @test density(pure(n"Au")) == density(n"Au")
        @test get(mat"Fe2O3", :Drip, "Drip") == "Drip"
        m = parse(
            Material,
            "Fe2O3",
            name = "Hematite",
            properties = Dict{Symbol,Any}(:Density => 5.25, :Drizzle => 99),
        )
        m[:Drop] = "Drop"
        @test get(m, :Density, "Drip") == 5.25
        @test get(m, :Drop, "Drip") == "Drop"
        @test m[:Drizzle] == 99
        @test typeof(m[:Drizzle]) == Int64
        lm = labeled(m)
        @test all(e -> lm[MassFractionLabel("Hematite", e)] == m[e], keys(m))
        m = asnormalized(mat"0.8*Fe+0.15*Ni+0.04*Cr", 1.0)
        @test isapprox(m[n"Fe"], 0.8081, atol = 0.0001)
        @test isapprox(m[n"Ni"], 0.1515, atol = 0.0001)
        @test isapprox(m[n"Cr"], 0.0404, atol = 0.0001)
        @test analyticaltotal(m) == 1.0
        lm = labeled(m)
        @test all(
            el -> lm[MassFractionLabel("N[0.8⋅Fe+0.15⋅Ni+0.04⋅Cr,1.0]", el)] == m[el],
            keys(m),
        )
        @test haskey(m, n"Fe")
        @test !haskey(m, n"Li")
        m = parsedtsa2comp("K411,(O:42.36),(Mg:8.85),(Si:25.38),(Ca:11.06),(Fe:11.21),4.0")
        @test name(m) == "K411"
        @test m[n"Fe"] == 0.1121
        @test density(m) == 4.0
        @test todtsa2comp(m) == "K411,(O:42.36),(Mg:8.85),(Si:25.38),(Ca:11.06),(Fe:11.21),4.0"
        cl = NeXLCore.compositionlibrary()
        @test cl["SRM93a"][n"Na"] == 0.0295
        @test isapprox(cl["NIST K2789"][n"U"], 0.2120, atol = 0.0001)
    end
    @testset "Film" begin
        f = Film(pure(n"C"), 1.0e-6)
        @test repr(f) == "10.0 nm of Pure C"
        @test isapprox(
            0.9673838,
            transmission(f, n"O K-L3", deg2rad(40.0)),
            atol = 0.000001,
        )

        @test material(f)[n"C"] == 1.0
        @test thickness(f) == 1.0e-6 # cm
    end
    @testset "nex" begin
        @test enx"Fe K" == energy(n"Fe K")
        @test enx"Fe K-L3" == energy(n"Fe K-L3")
    end
end


text:
# Interval
# ========
#
# Base interval types and utilities.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# Note, just to be clear: this shadows IntervalTrees.Interval
"A genomic interval specifies interval with some associated metadata"
struct Interval{T} <: IntervalTrees.AbstractInterval{Int64}
    seqname::String
    first::Int64
    last::Int64
    strand::Strand
    metadata::T
end

function Interval(seqname::AbstractString, first::Integer, last::Integer, strand::Union{Strand,Char}=STRAND_BOTH, metadata=nothing)
    return Interval{typeof(metadata)}(seqname, first, last, strand, metadata)
end

function Interval(seqname::AbstractString, range::UnitRange{T}, strand::Union{Strand,Char}=STRAND_BOTH, metadata=nothing) where T<:Integer
    return Interval{typeof(metadata)}(seqname, first(range), last(range), strand, metadata)
end

"""
    Interval{T}(data)

The returned data is converted to Interval{T} if there is an implemented [`Base.convert`](https://docs.julialang.org/en/v1/base/base/#Base.convert) function for the type of data.
This method provides a useful hook for converting custom types to Interval{T}.
"""
function Interval{T}(data) :: Interval{T} where T
    return data #Note: the returned data is converted to Interval{T}.
end

function BioGenerics.seqname(i::Interval)
    return i.seqname
end

function BioGenerics.metadata(i::Interval)
    return i.metadata
end

function strand(i::Interval)
    return i.strand
end

"""
    leftposition(i::Interval)

Return the leftmost position of `i`.
"""
function BioGenerics.leftposition(i::Interval)
    return i.first
end

"""
    rightposition(i::Interval)

Return the rightmost position of `i`.
"""
function BioGenerics.rightposition(i::Interval)
    return i.last
end

IntervalTrees.first(i::Interval) = leftposition(i)
IntervalTrees.last(i::Interval) = rightposition(i)

function Base.isless(a::Interval{T}, b::Interval{T}, seqname_isless::Function=isless) where T
    a_seqname = seqname(a)
    b_seqname = seqname(b)

    if a_seqname != b_seqname
        return seqname_isless(a_seqname, b_seqname)
    end

    a_leftposition = leftposition(a)
    b_leftposition = leftposition(b)

    if a_leftposition != b_leftposition
        return a_leftposition < b_leftposition
    end

    a_rightposition = rightposition(a)
    b_rightposition = rightposition(b)

    if a_rightposition != b_rightposition
        return a_rightposition < b_rightposition
    end

    a_strand = strand(a)
    b_strand = strand(b)

    if a_strand != b_strand
        return a_strand < b_strand
    end

    return false
end

"""
Check if two intervals are well ordered.

Intervals are considered well ordered if seqname(a) <= seqname(b) and leftposition(a) <= leftposition(b).
"""
function isordered(a::Interval{T}, b::Interval{T}, seqname_isless::Function=isless) where T

    a_seqname = seqname(a)
    b_seqname = seqname(b)

    if a_seqname != b_seqname
        return seqname_isless(a_seqname, b_seqname)
    end

    a_leftposition = leftposition(a)
    b_leftposition = leftposition(b)

    if a_leftposition != b_leftposition
        return a_leftposition < b_leftposition
    end

    return true
end

"""
Return true if interval `a` entirely precedes `b`.
"""
function precedes(a::Interval{T}, b::Interval{T}, seqname_isless::Function=isless) where T
    return (rightposition(a) < leftposition(b) && seqname(a) == seqname(b)) || seqname_isless(seqname(a), seqname(b))::Bool
end

function Base.:(==)(a::Interval{T}, b::Interval{T}) where T
    return seqname(a)       == seqname(b) &&
           leftposition(a)  == leftposition(b) &&
           rightposition(a) == rightposition(b) &&
           strand(a)        == strand(b) &&
           metadata(a)      == metadata(b)
end

"Return true if interval `a` overlaps interval `b`, with no consideration to strand"
function BioGenerics.isoverlapping(a::Interval{S}, b::Interval{T}) where {S, T}
    return leftposition(a) <= rightposition(b) &&
           leftposition(b) <= rightposition(a) &&
           seqname(a)      == seqname(b)
end

function Base.show(io::IO, i::Interval)
    if get(io, :compact, false)
        print(io, seqname(i), ":", leftposition(i), "-", rightposition(i), "  ", strand(i), "  ", metadata(i) === nothing ? "nothing" : metadata(i))
    else
        println(io, summary(i), ':')
        println(io, "  sequence name: ", seqname(i))
        println(io, "  leftmost position: ", leftposition(i))
        println(io, "  rightmost position: ", rightposition(i))
        println(io, "  strand: ", strand(i))
          print(io, "  metadata: ", metadata(i) === nothing ? "nothing" : metadata(i))
    end
end

function metadatatype(::Type{T}) where T
    return _metadatatype(eltype(T))
end

function metadatatype(x::Any)
    return metadatatype(typeof(x))
end

function _metadatatype(::Type{Interval{T}}) where T
    return T
end

function _metadatatype(::Type{T}) where T
    return T
end


text:
#
# Optimal data formatter for a fixed width screen.
# (see <https://research.google.com/pubs/pub44667.html>)
#

module PrettyPrinting

export
    pprint,
    pprintln

import Base:
    IndexStyle,
    getindex,
    size,
    show,
    *, /, |

const DEFAULT_LINE_WIDTH = 79
const DEFAULT_BREAK_COST = 1
const DEFAULT_SPILL_COST = 2

include("tile.jl")
include("expr.jl")
include("fit.jl")

"""
    pprint([io::IO], data)

Display the data so that it fits the width of the output screen.
"""
pprint(@nospecialize data) =
    pprint(stdout, data)

pprint(io::IO, data) =
    pprint(io, tile(data))

function pprint(io::IO, lt::Layout)
    render(io, best_fit(io, lt))
    nothing
end

"""
    pprintln([io::IO], data)

Display the data using `pprint` and print a line break.
"""
pprintln(@nospecialize data) =
    pprintln(stdout, data)

function pprintln(io::IO, data)
    pprint(io, data)
    println(io)
end

"""
    quoteof(obj)

Convert an object to its representation in Julia AST.

Implement this method to customize [`pprint`](@ref) on a user-defined type.
"""
function quoteof
end

"""
    tile(obj)

Convert an object to the corresponding layout expression.

Implement this method to customize [`pprint`](@ref) on a user-defined type.
"""
function tile
end

end


text:
module Nauty

using LightGraphs
using Permutations

import LightGraphs

export nauty, NautyIsomorphism

const LIB_FILE = "$(@__DIR__)" * "/../deps/nauty"

const WORDSIZE = ccall((:wordsize, LIB_FILE), Int, ())

include("dense.jl")
include("sparse.jl")

# data structure copied from https://github.com/bovine3dom/Nauty.jl

const Nbool = Cint
mutable struct optionblk
    getcanon::Cint	# make canong and canonlab?
    digraph::Nbool	# multiple edges or loops?
    writeautoms::Nbool	# write automorphisms?
    writemarkers::Nbool	# write stats on pts fixed, etc.?
    defaultptn::Nbool	# set lab,ptn,active for single cell?
    cartesian::Nbool	# use cartesian rep for writing automs?
    linelength::Cint	# max chars/line (excl. '\n') for output
    outfile::Ptr{Cvoid}	# FILE *outfile;
		# file for output, if any
    userrefproc::Ptr{Cvoid}	# void (*userrefproc)
		# replacement for usual refine procedure
		# (graph*,int*,int*,int,int*,int*,set*,int*,int,int);
    userautomproc::Ptr{Cvoid}	# void (*userautomproc)
		# procedure called for each automorphism
		# (int,int*,int*,int,int,int);
    userlevelproc::Ptr{Cvoid}	# void (*userlevelproc)
		# procedure called for each level
		# (int*,int*,int,int*,statsblk*,int,int,int,int,int,int);
    usernodeproc::Ptr{Cvoid}	# void (*usernodeproc)
		# procedure called for each node
		# (graph*,int*,int*,int,int,int,int,int,int);
    usercanonproc::Ptr{Cvoid}	# Cint  (*usercanonproc)
		# procedure called for better labellings
		# (graph*,int*,graph*,int,int,int,int);
    invarproc::Ptr{Cvoid}	# void (*invarproc)
		# procedure to compute vertex-invariant
		# (graph*,int*,int*,int,int,int,int*,int,Nbool,int,int);
    tc_level::Cint	# max level for smart target cell choosing
    mininvarlevel::Cint	# min level for invariant computation
    maxinvarlevel::Cint	# max level for invariant computation
    invararg::Cint	# value passed to (*invarproc)()
    dispatch::Ptr{Cvoid}	# dispatchvec *dispatch;
		# vector of object-specific routines
    schreier::Nbool	# use random schreier method
    extra_options::Ptr{Cvoid}	# void *extra_options;
		# arbitrary extra options
end

const CONSOLWIDTH = 78
const DISPATCH_GRAPH = Ref{Ptr{Cvoid}}(0)
const ADJACENCIES = Ref{Ptr{Cvoid}}(0)
function __init__()
    DISPATCH_GRAPH[] = cglobal((:dispatch_graph,LIB_FILE),Nothing)
    ADJACENCIES[] = cglobal((:adjacencies,LIB_FILE),Nothing)
    nothing
end
DEFAULTOPTIONS_GRAPH() = optionblk(0,false,false,false,true,false,CONSOLWIDTH,
                    C_NULL,C_NULL,C_NULL,C_NULL,C_NULL,C_NULL,C_NULL,
                    100,0,1,0,DISPATCH_GRAPH[],false,C_NULL)
DEFAULTOPTIONS_DIGRAPH() = optionblk(0,true,false,false,true,false,CONSOLWIDTH,
                    C_NULL,C_NULL,C_NULL,C_NULL,C_NULL,C_NULL,ADJACENCIES[],
                    100,0,999,0,DISPATCH_GRAPH[],false,C_NULL)
DEFAULTOPTIONS(::Type{DenseNautyGraph}) = DEFAULTOPTIONS_GRAPH()
DEFAULTOPTIONS(::Type{DenseNautyDiGraph}) = DEFAULTOPTIONS_DIGRAPH()
DEFAULTOPTIONS(::DenseNautyGraph) = DEFAULTOPTIONS_GRAPH()
DEFAULTOPTIONS(::DenseNautyDiGraph) = DEFAULTOPTIONS_DIGRAPH()

function pprintobject(name, obj, sep=", ")
  print("$name(")
  print(join(map(fn -> "$fn=$(getfield(obj, fn))", fieldnames(typeof(obj))), sep))
  print(")")
end

function Base.show(io::IO, ::MIME"text/plain", options::Nauty.optionblk)
    pprintobject("optionblk", options)
end

mutable struct statsblk
    grpsize1::Cdouble	# /* size of group is */
    grpsize2::Cint	# /* grpsize1 * 10^grpsize2 */
    numorbits::Cint	# /* number of orbits in group */
    numgenerators::Cint	# /* number of generators found */
    errstatus::Cint	# /* if non-zero : an error code */
    numnodes::Culong	# /* total number of nodes */
    numbadleaves::Culong	# /* number of leaves of no use */
    maxlevel::Cint	# /* maximum depth of search */
    tctotal::Culong	# /* total size of all target cells */
    canupdates::Culong	# /* number of updates of best label */
    invapplics::Culong	# /* number of applications of invarproc */
    invsuccesses::Culong	# /* number of successful uses of invarproc() */
    invarsuclevel::Cint	# /* least level where invarproc worked */
end

const statsblk() = statsblk(zeros(13)...)

function Base.show(io::IO, ::MIME"text/plain", stats::Nauty.statsblk)
    pprintobject("statsblk", stats)
end

"""densenauty is the low-level interface to nauty.

Its arguments match closely the call to densenauty: a graph, an option block, an optionally a partition, in the form of a pair of lists (lab,ptn).

The return value is (stats, orbits, ...) a stats block, a description of the minimal vertices (numbered from 0) in each orbit, and (if options.getcanon) a list (also numbered from 0) matching the graph's vertices with those in the canonical graph, and the canonical graph itself.
"""
function densenauty(g::DenseNautyGraph,
               options = DEFAULTOPTIONS_GRAPH()::optionblk,
                    partition = nothing::Union{Nothing,Tuple{Vector{Cint},Vector{Cint}}})
    __densenauty(g, options, partition)
end
function densenauty(g::DenseNautyDiGraph,
                    options = DEFAULTOPTIONS_DIGRAPH()::optionblk,
                    partition = nothing::Union{Nothing,Tuple{Vector{Cint},Vector{Cint}}})
    __densenauty(g, options, partition)
end

function __densenauty(g::DenseNautyXGraph, options::optionblk, partition)
    n = nv(g)
    stats = statsblk()

    if options.defaultptn == 1
        @assert partition == nothing
        lab = zeros(Cint, n)
        ptn = zeros(Cint, n)
    else
        @assert length(partition[1]) == length(partition[2]) == n
        lab = partition[1]
        ptn = partition[2]
    end

    if options.getcanon == 1
        outgraph = typeof(g)(n)
    else
        outgraph = nothing
    end

    orbits = zeros(Cint, n)

    ccall((:densenauty, LIB_FILE), Cvoid,
          (Ptr{UInt64}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ref{optionblk}, Ref{statsblk}, Cint, Cint, Ptr{UInt64}),
          g.data.chunks, lab, ptn, orbits, options, stats, num_setwords(n), n, outgraph == nothing ? C_NULL : outgraph.data.chunks)

    if options.getcanon == 1
        result = (lab, outgraph)
    else
        result = ()
    end

    (stats, orbits, result...)
end

# ⚠ not thread-safe!
# it doesn't seem possible to put this function inside the nauty function,
# because the @cfunction call can't see it.
const __gens = Permutation[]
function userautomproc_jl_global(count::Cint, permptr::Ptr{Cint}, orbitsptr::Ptr{Cint}, numorbits::Cint, stabvertex::Cint, n::Cint)
    perm = unsafe_wrap(Array, permptr, n)
    orbits = unsafe_wrap(Array, orbitsptr, n)

    push!(__gens, Permutation(perm.+1))
    nothing
end
"""nauty is a higher-level interface to nauty, which relies on densenauty.

The arguments are a graph g and optional named arguments getcanon::Bool, automgroup::Bool and partition, which is either "nothing" or a list of lists of vertices (numbered from 1).

The return value is a dictionary with entries
:orbits (orbits numbered from 1), an Int[]
:grpsize (possibly with rounding errors), a BigInt
:numorbits (the number of orbits under the automorphism group), an Int
:numgenerators (the number of group generators), an Int
if automgroup,
:generators (the generators of the automorphism group), a Permutation[]
if getcanon,
:lab (the correspondence between old vertices and new ones), a Permutation
:canong (the canonically labelled graph), a DenseNautyGraph
"""
function nauty(g::DenseNautyXGraph;
               getcanon = false::Bool,
               automgroup = false::Bool,
               partition = nothing::Union{Nothing,Tuple{Vector{Cint},Vector{Cint}},Vector{Vector{Int}}})
    options = DEFAULTOPTIONS(g)

    if getcanon
        options.getcanon = true
    end

    if partition == nothing
        labptn = nothing
    elseif isa(partition,Tuple)
        options.defaultptn = false
        labptn = partition
    else
        options.defaultptn = false
        n = nv(g)
        lab = Vector{Cint}(undef,n)
        ptn = ones(Cint,n)
        i = 0
        for p=partition
            @assert !isempty(p)
            lab[i+1:i+length(p)] = p.-1
            i += length(p)
            ptn[i] = 0
        end
        labptn = (lab,ptn)
    end

    local generators
    function userautomproc_jl(count::Cint, permptr::Ptr{Cint}, orbitsptr::Ptr{Cint}, numorbits::Cint, stabvertex::Cint, n::Cint)
        perm = unsafe_wrap(Array, permptr, n)
        orbits = unsafe_wrap(Array, orbitsptr, n)

        push!(generators, Permutation(perm.+1))
        nothing
    end

    if automgroup
        generators = Permutation[]
        userautomproc_c = @cfunction($userautomproc_jl, Nothing, (Cint, Ptr{Cint}, Ptr{Cint}, Cint, Cint, Cint))
        options.userautomproc = unsafe_convert(Ptr{Cvoid},userautomproc_c)
        GC.@preserve userautomproc_c rv = densenauty(g, options, labptn)

    else
        rv = densenauty(g, options, labptn)
    end

    stats = rv[1]
    stats.errstatus == 0 || error("densenauty: error $(stats.errstatus)")

    result = Dict()
    push!(result, :orbits => rv[2].+1)
    push!(result, :grpsize => div(BigInt(ldexp(significand(stats.grpsize1),precision(Float64)))*BigInt(10)^stats.grpsize2*BigInt(2)^exponent(stats.grpsize1)+BigInt(2)^(precision(Float64)-1),BigInt(2)^precision(Float64)))
    push!(result, :numorbits => Int(stats.numorbits))
    push!(result, :numgenerators => Int(stats.numgenerators))
    if automgroup
        push!(result, :generators => generators)
    end
    if getcanon
        push!(result, :lab => rv[3].+1)
        push!(result, :canong => rv[4])
    end
    result
end

"""NautyIsomorphism() may be supplied as "alg" argument to LightGraphs.Experimental.has_isomorphism
"""
struct NautyIsomorphism <: LightGraphs.Experimental.IsomorphismAlgorithm
end

function LightGraphs.Experimental.has_isomorph(g1::AbstractGraph, g2::AbstractGraph, ::NautyIsomorphism;
                         vertex_relation::Union{Nothing, Function}=nothing,
                         edge_relation::Union{Nothing, Function}=nothing)::Bool
    !LightGraphs.Experimental.could_have_isomorph(g1, g2) && return false

    vertex_relation == nothing || error("didn't code vertex relations yet")
    edge_relation == nothing || error("didn't code edge relations yet")

    return nauty(DenseNautyGraph(g1),getcanon=true)[:canong] ==
        nauty(DenseNautyGraph(g2),getcanon=true)[:canong];
end

#!!! sparse version

end


text:
using OrdinaryDiffEq, ODE, ODEInterfaceDiffEq, LSODA, Sundials, DiffEqDevTools, Plots; gr()

## Define the ThreeBody Problem
const threebody_μ = parse(Float64,"0.012277471")
const threebody_μ′ = 1 - threebody_μ

f = (du,u,p,t) -> begin
  @inbounds begin
  # 1 = y₁
  # 2 = y₂
  # 3 = y₁'
  # 4 = y₂'
  D₁ = ((u[1]+threebody_μ)^2 + u[2]^2)^(3/2)
  D₂ = ((u[1]-threebody_μ′)^2 + u[2]^2)^(3/2)
  du[1] = u[3]
  du[2] = u[4]
  du[3] = u[1] + 2u[4] - threebody_μ′*(u[1]+threebody_μ)/D₁ - threebody_μ*(u[1]-threebody_μ′)/D₂
  du[4] = u[2] - 2u[3] - threebody_μ′*u[2]/D₁ - threebody_μ*u[2]/D₂
  end
end

t₀ = 0.0; T = parse(Float64,"17.0652165601579625588917206249")
tspan = (t₀,2T)

prob = ODEProblem(f,[0.994, 0.0, 0.0, parse(Float64,"-2.00158510637908252240537862224")],tspan)

test_sol = TestSolution(T,[prob.u0])
abstols = 1.0 ./ 10.0 .^ (3:13); reltols = 1.0 ./ 10.0 .^ (0:10);

sol = solve(prob,Vern9(),abstol=1e-14,reltol=1e-14)
@show sol[1] - sol[end]
@show sol[end] - prob.u0;

apr = appxtrue(sol,test_sol)
@show sol[end]
@show apr.u[end]
@show apr.errors

setups = [Dict(:alg=>DP5())
          #Dict(:alg=>ode45()) #fails
          Dict(:alg=>BS5())
          Dict(:alg=>Tsit5())
          Dict(:alg=>dopri5())];
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=100)
plot(wp)

setups = [Dict(:alg=>DP5(),:dense=>false)
          #Dict(:alg=>ode45()) # Fails
          Dict(:alg=>BS5(),:dense=>false)
          Dict(:alg=>Tsit5(),:dense=>false)
          Dict(:alg=>dopri5())];
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,numruns=100)
plot(wp)

setups = [Dict(:alg=>DP5())
          #Dict(:alg=>ode45()) #fails
          Dict(:alg=>BS5())
          Dict(:alg=>Tsit5())
          Dict(:alg=>dopri5())];
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,numruns=100)
plot(wp)

setups = [Dict(:alg=>DP5())
          Dict(:alg=>Vern6())
          Dict(:alg=>Vern7())
          Dict(:alg=>TanYam7())
          Dict(:alg=>Vern8())
          Dict(:alg=>DP8())
          Dict(:alg=>dop853())
          Dict(:alg=>Vern9())];
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=100)
plot(wp)

setups = [Dict(:alg=>DP5())
          Dict(:alg=>Vern6())
          Dict(:alg=>Vern7())
          Dict(:alg=>TanYam7())
          Dict(:alg=>Vern8())
          Dict(:alg=>DP8())
          Dict(:alg=>dop853())
          Dict(:alg=>Vern9())];
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,dense=false,numruns=100,verbose=false)
plot(wp)

setups = [Dict(:alg=>DP5())
          Dict(:alg=>Vern6())
          Dict(:alg=>Vern7())
          Dict(:alg=>TanYam7())
          Dict(:alg=>Vern8())
          Dict(:alg=>DP8())
          Dict(:alg=>dop853())
          Dict(:alg=>Vern9())];
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,numruns=100)
plot(wp)

#setups = [Dict(:alg=>ode78())
#          Dict(:alg=>VCABM())
#          Dict(:alg=>CVODE_Adams())];
#wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,dense=false,numruns=100)

setups = [Dict(:alg=>DP5())
          Dict(:alg=>lsoda())
          Dict(:alg=>Vern8())
          Dict(:alg=>ddeabm())
          Dict(:alg=>odex())
          Dict(:alg=>ARKODE(Sundials.Explicit(),order=6))
    ];
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=100)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (3:13); reltols = 1.0 ./ 10.0 .^ (0:10);
setups = [Dict(:alg=>Tsit5())
          Dict(:alg=>Vern9())
          Dict(:alg=>AitkenNeville(min_order=1, max_order=9, init_order=4, threading=true))
          Dict(:alg=>ExtrapolationMidpointDeuflhard(min_order=1, max_order=9, init_order=4, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, threading=true))]
solnames = ["Tsit5","Vern9","AitkenNeville","Midpoint Deuflhard","Midpoint Hairer Wanner"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,names=solnames,
                      save_everystep=false,verbose=false,numruns=100)
plot(wp)

setups = [Dict(:alg=>ExtrapolationMidpointDeuflhard(min_order=1, max_order=9, init_order=9, threading=false))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, threading=false))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, sequence = :romberg, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, sequence = :bulirsch, threading=true))]
solnames = ["Deuflhard","No threads","standard","Romberg","Bulirsch"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,names=solnames,
                      save_everystep=false,verbose=false,numruns=100)
plot(wp)

setups = [Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=10, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=5, max_order=11, init_order=10, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=15, init_order=10, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=5, max_order=7, init_order=6, threading=true))]
solnames = ["1","2","3","4","5"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,names=solnames,
                      save_everystep=false,verbose=false,numruns=100)
plot(wp)

using SciMLBenchmarks
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

max_stars_repo_name: Gregliest/SciMLBenchmarksOutput

max_issues_repo_name: Gregliest/SciMLBenchmarksOutput

max_forks_repo_name: Gregliest/SciMLBenchmarksOutput


text:
#structs & enums
@enum SLOPE begin
  decrease = -1
  undefined = 0
  increase = 1
end

mutable struct DEPTH_LOG
  index::Int
  length::Int
  slope::Int
  end

mutable struct DEPTH_LOG_WINDOW
  index::Int
  value::Int
end

println("#######################")
println("### AoC2021 - Day 1 ###")
println("#######################")

println(">> @arg1 : Input file")
println()

if ( size(ARGS)[1] != 1 )
  println("FATAL ERROR : Invalid number of args. Input file is required as @arg1.")
  quit()
end

window_size = 5

# open input file
fd = open(ARGS[1], "r")
vals = readlines(ARGS[1])
if (size(vals)[1] <= 1)
  println("INVALID INPUT : Not enough data for meaningful answer.")
  quit()
end

# convert to int array
measures = [parse(Int,x) for x in vals]

# count depth changes
res = DEPTH_LOG[]
depth_log = DEPTH_LOG( 0, 0, 0)
n = size(measures)[1]

for i =1:n

  if (i==1)
    depth_log.length += 1
    continue;
  end

  if measures[i] > measures[i-1] # increase
    if (depth_log.slope == 0)
      depth_log.slope = 1
    end
    if (depth_log.slope != 1)
      push!(res, depth_log)
      global depth_log = DEPTH_LOG( i, 0, 1)
    end
    depth_log.length+=1
  elseif measures[i] < measures[i-1] # decrease
    if (depth_log.slope == 0)
      depth_log.slope = -1
    end
    if (depth_log.slope != -1)
      push!(res, depth_log)
      global depth_log = DEPTH_LOG( i, 0, -1)
    end
    depth_log.length+=1
  end

  if (i==n)
    push!(res, depth_log)
  end;

end #!for

# prints logs
# println("~~~ DEPTH LOGS ~~~")
# for result in res
#   println(result)
# end
# println()

# print AoC result
cpt = 0
for log in res
  if log.slope == 1
    global cpt += log.length
  end
end
println("Number of larger than previous measures : $cpt")

# - part 2
res2 = DEPTH_LOG_WINDOW[]
for i=window_size:n
  depth_log_wdw = DEPTH_LOG_WINDOW(i, 0)
  for j=0:window_size-1
    depth_log_wdw.value += measures[i-j]
  end
  push!(res2, depth_log_wdw)
end

# Find how many sliding windows values are superior than their precedecessors
cpt2=0
for k=1:size(res2)[1]
  if (k==1)
    continue
  end
  if res2[k].value > res2[k-1].value
    global cpt2 += 1
  end
end
println("Number of larger than previous sliding windows of size $window_size : $cpt2")

# close file
close(fd)


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.6.3
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://twitter.com/genkuroki/status/1459220790444920833

# %%
using Distributions, StatsPlots

# %%
mixnormal = MixtureModel([Normal(0, 1), Normal(1, 2)], [1/3, 2/3])

# %%
X = rand(mixnormal, 10^5)
f(x) = (1/3)*pdf(Normal(0, 1), x) + (2/3)*pdf(Normal(1, 2), x)

histogram(X; norm=true, alpha=0.3, label="sample of mixnormal")
plot!(f; label="(1/3)N(x|0,1)+(2/3)N(x|1,2))", lw=2)

# %%

max_stars_repo_path: 0024/mixed normal dist.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0024/mixed normal dist.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0024/mixed normal dist.jl

max_forks_repo_name: genkuroki/public


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Julia 1.3.1
#     language: julia
#     name: julia-1.3
# ---

# # Example 4
#
# Two three-dimensional variables are read from the netcdf files generated in `Example3`, combined into a single data structure, and then re-written together into a new netcdf file.

# +
using NCTiles,NCDatasets

inputs = "input/"
outputs = "output/"
nt="0003"

file_in1=outputs*"ex3/THETA/THETA.$nt.nc"
file_in2=outputs*"ex3/SALT/SALT.$nt.nc"
~isfile(file_in1) ? error("Running Example3 first is needed to run Example4") : nothing

file_out=outputs*"ex4/TS.$nt.nc"
if ~ispath(outputs*"ex4/"); mkpath(outputs*"ex4/"); end

README = readlines(joinpath(inputs,"README"))

# Get the first 3D variable
ncvars,ncdims,fileatts = readncfile(file_in1)

# Add a second 3D variable
tmp,_,_ = readncfile(file_in2)
T=ncvars["THETA"]
S=tmp["SALT"]
ncvars["SALT"]=NCvar(S.name,S.units,T.dims,S.values,S.atts,T.backend)

# Rewrite to a new file
write(ncvars,file_out,README=README,globalattribs=fileatts)
# -


text:
#=
Script to contract a given tensor network with a given contraction plan
=#

using PicoQuant
using ArgParse
using JSON

"""
    parse_commandline()

Parse command line options and return argument dictionary
"""
function parse_commandline(ARGS)
    s = ArgParseSettings()
        @add_arg_table! s begin
            "--dsl", "-f"
                help = "Name of the dsl file to load"
                required = true
            "--dsl_data", "-d"
                help = "H5 file containing the tensor data"
                required = true
            "--output", "-o"
                help = "H5 file to output data to"
                arg_type = String
                required = true
        end
        return parse_args(ARGS, s)
end

function main(ARGS)
    parsed_args = parse_commandline(ARGS)

    dsl_filename = parsed_args["dsl"]
    data_filename = parsed_args["dsl_data"]
    output_filename = parsed_args["output"]

    DSLBackend(dsl_filename, data_filename, output_filename)

    execute_dsl_file(dsl_filename, data_filename, output_filename)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main(ARGS)
end


text:
function Maneuver(arg0::AttitudeProvider, arg1::ManeuverTriggers, arg2::PropulsionModel)
    return Maneuver((AttitudeProvider, ManeuverTriggers, PropulsionModel), arg0, arg1, arg2)
end


text:
function next(side, kol)
    return HorizonSide((Int(side) + kol) % 4)
end

function goAround(side)
    kol = 0
    while isborder(robot, side)
        if !isborder(robot, next(side, 1))
            move!(robot, next(side, 1))
            kol += 1
        else
            break
        end
    end
    if isborder(robot, side)
        for i = 1:kol
            if !isborder(robot, next(side, 3))
                move!(robot, next(side, 3))
            end
        end
        return false
    else
        move!(robot, side)
        for i = 1:kol
            if !isborder(robot, next(side, 3))
                move!(robot, next(side, 3))
            end
        end
        return true
    end
end

function moveAndKol(side, kol)
    flag = false
    ok = true
    count = 0
    while true
        if (isborder(robot, next(side, kol)) && !flag)
            flag = true
            count += 1
        elseif !isborder(robot, next(side, kol))
            flag = false
        end
        if (isborder(robot, side))
            ok = goAround(side)
        end
        if ok
            if !isborder(robot, side)
                move!(robot, side)
            end
        else
            break
        end
    end
    return count
end

function kol(side, move_side, ok)
    count = 0
    n = 1
    if ok
        n = 1
    else
        n = 3
    end
    while !isborder(robot, side)
        count += moveAndKol(move_side, n)
        move_side = next(move_side, 2)
        n = (n + 2) % 4
        move!(robot, side)
    end
    return count
end

function f()
    count = 0
    count += kol(Nord, Ost, true)
    while !isborder(robot, Sud)
        move!(robot, Sud)
    end
    count += kol(Ost, Nord, false)
    println(count)
end

f()

max_stars_repo_name: TheXlebushek/mirea-progs

max_issues_repo_name: TheXlebushek/mirea-progs

max_forks_repo_name: TheXlebushek/mirea-progs


text:
using SeisReconstruction
using SeisProcessing
using LinearAlgebra
using Test

d = SeisLinearEvents()

#randomly decimate 50% of traces

mask = rand(1,size(d,2));
perc=50
mask[(LinearIndices(mask .< perc/100))[findall(mask .< perc/100)]] .= 0;
mask[(LinearIndices(mask .>= perc/100))[findall(mask .>= perc/100)] ] .= 1;

ddec = Array{Float64}(undef,size(d));

for it = 1 : size(d,1)
    ddec[it,:] = d[it:it,:].*mask;
end

dmwni = SeisMWNI(ddec,dt=0.004,fmax=100,Niter_external=5,Niter_internal=10)

# test that quality factor is greater than 10 Decibels
quality_factor = 10*log10(norm(d[:],2)/norm(dmwni[:]-d[:],2))
println("Quality factor = ",quality_factor)
@test quality_factor > 5


text:
module p1

# Return boundary conditions values and indices, ϕb and β, for problem 1,
# with N subdivisions
function boundaryconditions1(N::Integer)
    ϕb = zeros(N+1,N+1)                 # boundary conditions
    β  = zeros(Int,N+1,N+1)             # boundary condition locations

    # Conditions on the perimeter
    ϕb[:,1] = ϕb[:,N+1] = ϕb[1,:] = ϕb[N+1,:] = 0.0
    β[:,1]  = β[:,N+1]  = β[1,:]  = β[N+1,:]  = 1

    # Conditions on the interior
    # xmin < x < xmax  →  N*xmin + 1 < i < N*xmax + 1
    xmin = 0.5
    xmax = 0.75
    ymin = 0.625
    ymax = 0.875

    xs = int([ceil(N*xmin+1):floor(N*xmax+1)])
    ys = int([ceil(N*ymin+1):floor(N*ymax+1)])

    ϕb[xs,ys] = 10.0
    β[xs,ys]  = 1

    return ϕb, β
end

# Return the charge distribution for problem 1
function chargedensity(N::Integer, x::Real, y::Real, q::Real)
    ρ = zeros(N+1,N+1)

    x <= 1.0 || error("x must be less than 1")
    y <= 1.0 || error("y must be less than 1")

    # Spread the point charge across four points; treat it like uniform density
    ρ[floor(x*N+1),floor(y*N+1)] = N*N*q/4
    ρ[floor(x*N+2),floor(y*N+1)] = N*N*q/4
    ρ[floor(x*N+1),floor(y*N+2)] = N*N*q/4
    ρ[floor(x*N+2),floor(y*N+2)] = N*N*q/4

    return ρ
end

chargedensity1(N::Integer) = chargedensity(N, 0.5, 0.25, 20)

# Get the grid location nearest to a given x,y value
function nearestpoint(x::Real, y::Real, N::Integer)

    i = int(x*N+1)
    j = int(y*N+1)

    return i, j

end

end


text:
import Base: ==

# Basic julia sugar

function desugar(ex)
  MacroTools.prewalk(ex) do ex
    @capture(ex, (xs__,)) ? :(tuple($(xs...))) :
    @capture(ex, xs_[i__]) ? :(getindex($xs, $(i...))) :
    ex
  end
end

# Constants

immutable Constant{T}
  value::T
end

tocall(c::Constant) = c.value

isconstant(v::Vertex) = isa(value(v), Constant)

mapconst(f, g) = map(x -> isa(x, Constant) ? Constant(f(x.value)) : f(x), g)

a::Constant == b::Constant = a.value == b.value

Base.hash(c::Constant, h::UInt = UInt(0)) = hash((Constant, c.value), h)

for (c, v) in [(:constant, :vertex), (:dconstant, :dvertex)]
  @eval $c(x) = $v(Constant(x))
  @eval $c(v::Vertex) = $v(v)
end

type Do end

tocall(::Do, a...) = :($(a...);)

# Line Numbers

immutable Line
  file::String
  line::Int
end

const noline = Line("", -1)

function Line(ex::Expr)
  @assert ex.head == :line
  Line(string(ex.args[2]), ex.args[1])
end

function normlines(ex)
  line = noline
  ex′ = :(;)
  for ex in ex.args
    isline(ex) && (line = Line(ex); continue)
    line == noline && (push!(ex′.args, ex); continue)
    @assert @capture(ex, var_ = val_)
    push!(ex′.args, :($var = $line($val)))
  end
  return ex′
end

function applylines(ex)
  ex′ = :(;)
  for ex in ex.args
    @capture(ex, (var_ = val_) | val_)
    val = MacroTools.postwalk(val) do ex
      @capture(ex, l_Frame(x_)) && return x # Ignore frames for now
      @capture(ex, l_Line(x_)) || return ex
      push!(ex′.args, Expr(:line, l.line, symbol(l.file)))
      @gensym edge
      push!(ex′.args, :($edge = $x))
      return edge
    end
    isexpr(val, Symbol) ? (ex′.args[end].args[1] = val) :
      push!(ex′.args, var == nothing ? :($var = $val) : val)
  end
  return ex′
end

immutable Frame{T}
  f::T
end

immutable SkipFrame end

# Static tuples

# TODO: just use `getindex` and `tuple` to represent these?
immutable Split
  n::Int
end

# TODO: printing
function normsplits(ex)
  MacroTools.prewalk(ex) do ex
    @capture(ex, (xs__,) = y_) || return ex
    @gensym edge
    quote
      $edge = $y
      $((:($(xs[i]) = $(Split(i))($edge)) for i = 1:length(xs))...)
    end
  end |> MacroTools.flatten |> block
end

tocall(::typeof(tuple), args...) = :($(args...),)

tocall(s::Split, x) = :($x[$(s.n)])

group(xs...) = vertex(tuple, xs...)

function detuple(v::IVertex)
  postwalk(v) do v
    if isa(value(v), Split) && value(v[1]) == tuple
      v[1][value(v).n]
    else
      v
    end
  end
end

# Bindings

immutable Bind
  name::Symbol
end

# TODO: printing
function insertbinds(ex)
  ls = map(ex.args) do l
    @capture(l, x_ = y_) || return l
    :($x = $(Bind(x))($y))
  end
  :($(ls...);)
end

# Inputs

immutable Input end

splitnode(v, n) = vertex(Split(n), v)

inputnode(n) = splitnode(constant(Input()), n)

isinput(v::IVertex) = isa(value(v), Split) && value(v[1]) == Constant(Input())

function bumpinputs(v::IVertex)
  prewalk(v) do v
    isinput(v) ?
      inputnode(value(v).n + 1) :
      v
  end
end

function spliceinput(v::IVertex, input::IVertex)
  postwalk(v) do v
    value(v) == Constant(Input()) ? input : v
  end
end

spliceinputs(v::IVertex, inputs::Vertex...) =
  spliceinput(v, group(inputs...))

function graphinputs(v::IVertex)
  n = 0
  prewalk(v) do v
    isinput(v) && (n = max(n, value(v).n))
    v
  end
  return n
end

# Closures

immutable Flosure end
immutable LooseEnd end

# TODO: scope
function normclosures(ex)
  bs = bindings(ex)
  MacroTools.prewalk(shortdef(ex)) do ex
    @capture(ex, (args__,) -> body_) || return ex
    @assert all(arg -> isa(arg, Symbol), args)
    closed = filter(x -> inexpr(body, x), bs)
    vars = vcat(closed, args)
    body = MacroTools.prewalk(body) do ex
      ex in vars ?
        Expr(:call, Split(findfirst(x->x==ex, vars)), LooseEnd()) :
        ex
    end
    :($(Flosure())($body, $(closed...)))
  end |> MacroTools.flatten |> block
end

flopen(v::IVertex) = mapconst(x->x==LooseEnd()?Input():x,v)


text:
using InplaceOps
using Test, LinearAlgebra

n = 10

X = randn(n,n)

v = randn(n)
v_orig = v

Z = randn(n,n)
Z_orig = Z

C = cholesky(X'X)
v_check = C \ v
@! v = C \ v
@test v == v_check
@test v === v_orig

A = qr(X)
v_check = A.Q * v
@! v = A.Q * v
@test v == v_check
@test v === v_orig

@test_throws MethodError (u = ones(10); @! u = A*u)

Z_check = X * X
@! Z = X * X
@test Z == Z_check
@test Z === Z_orig

@test_throws MethodError (U = ones(10,10); @! U = U*X)
@test_throws MethodError (U = ones(10,10); @! U = X*U)
@test_throws MethodError (U = ones(10,10); @! U *= X)

Z_check = X' * X
@! Z = X' * X
@test Z == Z_check
@test Z === Z_orig

Z_check = X * X'
@! Z = X * X'
@test Z == Z_check
@test Z === Z_orig


text:
module AliasedSeabed

export blackwell_asbmask, m1, m2

using Statistics
using ImageFiltering

function m1(ntheta; Ttheta=702, dtheta=28)
    meanfilter(ntheta, dtheta, dtheta).^2 .> Ttheta
end

function m2(nphi; Tphi=282, dphi=52)
    meanfilter(nphi, dphi, dphi).^2  .> Tphi
end

function myfindn(x)
    I = findall(!iszero, x)
    (getindex.(I, 1), getindex.(I, 2))
end

function meanfilter(A, height, width)
    kernel  = Array{Float64}(undef, height, width)
    kernel .= 1 / (width * height)
    kernel = centered(kernel)
    imfilter(A, kernel, Fill(0))
end

"""
    bwselect(BW, c , r)
Select objects in a binary image.
Similar to the MATLAB function of the same name.
"""
function bwselect(BW, c, r)
    # constants
    north = CartesianIndex(-1,  0)
    south = CartesianIndex( 1,  0)
    east  = CartesianIndex( 0,  1)
    west  = CartesianIndex( 0, -1)

    queue = CartesianIndex.(r,c)

    m,n = size(BW)
    out = falses(m,n)

    while !isempty(queue)
        node = pop!(queue)

        if BW[node]
            wnode = node
            enode = node + east

             # Move west until node is false
            while checkbounds(Bool, BW, wnode) && BW[wnode]
                out[wnode] = true
                if checkbounds(Bool, BW, wnode + north) && !out[wnode + north]
                    push!(queue, wnode + north)
                end
                if checkbounds(Bool, BW, wnode + south) && !out[wnode + south]
                    push!(queue, wnode + south)
                end
                wnode += west
            end

            # Move east until node is false
            while checkbounds(Bool, BW, enode) && BW[enode]
                out[enode] = true
                if checkbounds(Bool, BW, enode + north) && !out[enode + north]
                    push!(queue, enode + north)
                end
                if checkbounds(Bool, BW, enode + south) && !out[enode + south]
                    push!(queue, enode + south)
                end
                enode += east
            end
        end

    end
    return out
end

"""

Aliased seabed detection according to [Blackwell et al. (2019)](https://arxiv.org/abs/1904.10736)

"""
function blackwell_asbmask(Sv, ntheta, nphi; Ttheta=702, Tphi=282, dtheta=28, dphi=52, minSv=nothing)

    _m1 = m1(ntheta, Ttheta=Ttheta, dtheta=dtheta)
    _m2 = m2(nphi, Tphi=Tphi, dphi=dphi)

    m = (_m1 .| _m2)

    a = Sv[m]

    b = filter(x -> x !=-999, a) # remove 'missing' cells

    if length(b) < 1
        # No ASB found
        return map(x -> false, Sv)
    end

    md = median(b)

    @info("medianSv = $md")

    if minSv == nothing
        T = md
    else
        T = max(md, minSv)
    end

    @info("T = $T")

    aboveT = Sv.> T # Points above lower threshold.
    rows, cols = myfindn(m)

    bw = bwselect(aboveT, cols, rows)

    bw .| m

end

end # module


text:
@testset "Upgrading" begin
    @testset "JLSO.upgrade_jlso" begin
        @testset "no change for current version" begin
            d = Dict("metadata" => Dict("version" => v"3.0"))
            new_d = @suppress_out upgrade_jlso(d)
            @test new_d == Dict("metadata" => Dict("version" => v"3.0"))
        end

        @testset "upgrade 1.0" begin
            d = Dict("metadata" => Dict("version" => v"1.0", "format"=>:serialize))
            new_d = @suppress_out upgrade_jlso(d)
            @test new_d["metadata"]["version"] == "3"
            @test new_d["metadata"]["format"] == "julia_serialize"
            @test new_d["metadata"]["compression"] == "none"
            @test isempty(new_d["metadata"]["project"])
            @test isempty(new_d["metadata"]["manifest"])
            @test isempty(new_d["objects"])

            @testset "Don't rename bson format" begin
                d = Dict("metadata" => Dict("version" => v"1.0", "format"=>:bson))
                new_d = @suppress_out upgrade_jlso(d)
                @test new_d["metadata"]["version"] == "3"
                @test new_d["metadata"]["format"] == "bson"
                @test new_d["metadata"]["compression"] == "none"
                @test isempty(new_d["metadata"]["project"])
                @test isempty(new_d["metadata"]["manifest"])
                @test isempty(new_d["objects"])
            end

            @testset "Invalid VersionNumber" begin
                d = Dict(
                    "metadata" => Dict(
                        "version" => v"1.0", "format" => :bson,
                        # We don't intend to go back an make a 0.3.14 release at any point.
                        "pkgs" => merge(pkgs, Dict("JLSO" => v"0.3.14")),
                    )
                )

                new_d = @suppress_out upgrade_jlso(d)
                @test !isempty(new_d["metadata"]["project"])
                @test !isempty(new_d["metadata"]["manifest"])

                # Test that the project.toml has the JLSO dep, but not a specific version
                # (fallback when package versions don't exist in the registry)
                _project = Pkg.TOML.parse(new_d["metadata"]["project"])
                @test haskey(_project, "deps")
                @test haskey(_project["deps"], "JLSO")
                @test haskey(_project, "compat")
                # Check that the compat section doesn't have our package version
                @test !haskey(_project["compat"], "JLSO")
            end
        end
    end

    # The below is how we saves the specimens for compat testing
    # JLSO.save("specimens/v1_serialize.jlso", datas; format=:serialize)
    # JLSO.save("specimens/v1_bson.jlso", datas; format=:bson)

    @testset "Can still load old files" begin
        dir = joinpath(@__DIR__, "specimens")
        @testset "$fn" for fn in readdir(dir)
            if Int === Int32 && occursin("_serialize", fn)
                # Julia 64-bit serializations won't load on 32-bit systems
                @test_warn(JLSO.LOGGER, r"MethodError*", JLSO.load(joinpath(dir, fn)))
            else
                jlso_data = @suppress_out JLSO.load(joinpath(dir, fn))
                @test jlso_data == datas
            end
        end
    end

    # Test upgrading files to avoid repeated upgrading
    @testset "JLSO.upgrade" begin
        src_dir = joinpath(@__DIR__, "specimens")
        @testset "upgrade" begin
            mktempdir() do dest_dir
                @testset "$fn" for fn in readdir(src_dir)
                    src = joinpath(src_dir, fn)
                    dest = joinpath(dest_dir, fn)
                    @suppress_out JLSO.upgrade(src, dest)

                    if Int === Int32 && occursin("_serialize", fn)
                        # Julia 64-bit serializations won't load on 32-bit systems
                        @test_warn(JLSO.LOGGER, r"MethodError*", JLSO.load(dest))
                    else
                        jlso_data = JLSO.load(dest)
                        @test jlso_data == datas
                    end
                end
            end
        end
        @testset "batch upgrade" begin
            mktempdir() do dest_dir
                fns = readdir(src_dir)
                src = joinpath.(src_dir, fns)
                dest = joinpath.(dest_dir, fns)

                @suppress_out JLSO.upgrade(src, dest)
                jlso_data = JLSO.load(first(dest))
                @test jlso_data == datas
            end
        end
        @testset "upgrade w/ manifest" begin
            src = joinpath(src_dir, "v3_bson_none.jlso")
            orig_jlso = read(src, JLSOFile)

            # Extract the project & manifest to modify
            new_project = deepcopy(orig_jlso.project)
            new_manifest = deepcopy(orig_jlso.manifest)

            # Check the current state is what we expect
            @test new_project["compat"]["Memento"] == "0.10, 0.11, 0.12"
            @test new_manifest["Memento"][1]["version"] == "0.12.1"

            # Update those values
            new_project["compat"]["Memento"] = "0.11, 0.12"
            new_manifest["Memento"][1]["version"] = "0.12.0"

            mktempdir() do dest_dir
                dest = joinpath(dest_dir, "v3_bson_none.jlso")
                @suppress_out JLSO.upgrade(src, dest, new_project, new_manifest)
                new_jlso = read(dest, JLSOFile)

                @test new_project["compat"]["Memento"] == "0.11, 0.12"
                @test new_manifest["Memento"][1]["version"] == "0.12.0"
            end
        end
    end
end

max_stars_repo_name: invenia/JLSO

max_issues_repo_name: invenia/JLSO

max_forks_repo_name: invenia/JLSO


text:
using DelimitedFiles

function step_fish(input_fishnum)
    if input_fishnum == 0
        return 6
    else
        return input_fishnum - 1
    end
end

function step_day(input_list)
    num_fish_to_add = count(x->(x==0), input_list)
    output_list = map(x->step_fish(x), input_list)
    append!(output_list, ones(Int, num_fish_to_add)*8)
    return output_list
end

function test_step_day()
    test_list = [3,2,1]
    @assert step_day(test_list) == [2,1,0]
    test_list = [0,2,1]
    @assert step_day(test_list) == [6,1,0,8]
    test
end

function step_multiple_days(fish_list, num_days)
    for i in 1:num_days
        println("Processing day ", i, " num_fish is ", length(fish_list))
        fish_list = step_day(fish_list)
    end
    return fish_list
end

# read input file, note single quotes around the delimiter
input = readdlm("input.txt", ',',Int)
# convert to a column vector
input = vec(input)

# step 80 days
output = step_multiple_days(input, 80)
print("number of fish is ", length(output))


text:
#
# A simple and flexible atom implementation
#
export Atom, atomic_system, periodic_system, isolated_system

# Valid types for atom identifiers
const AtomId = Union{Symbol,AbstractString,Integer}

struct Atom{D, L<:Unitful.Length, V<:Unitful.Velocity, M<:Unitful.Mass}
    position::SVector{D, L}
    velocity::SVector{D, V}
    atomic_symbol::Symbol
    atomic_number::Int
    atomic_mass::M
    data::Dict{Symbol, Any}  # Store arbitrary data about the atom.
end
velocity(atom::Atom)      = atom.velocity
position(atom::Atom)      = atom.position
atomic_mass(atom::Atom)   = atom.atomic_mass
atomic_symbol(atom::Atom) = atom.atomic_symbol
atomic_number(atom::Atom) = atom.atomic_number
element(atom::Atom)       = elements[atomic_symbol(atom)]
n_dimensions(atom::Atom{D}) where {D} = D

Base.hasproperty(at::Atom, x::Symbol) = hasfield(Atom, x) || haskey(at.data, x)
Base.getproperty(at::Atom, x::Symbol) = hasfield(Atom, x) ? getfield(at, x) : getindex(at.data, x)
function Base.propertynames(at::Atom, private::Bool=false)
    if private
        (fieldnames(Atom)..., keys(at.data)...)
    else
        (filter(!isequal(:data), fieldnames(Atom))..., keys(at.data)...)
    end
end

function Atom(identifier::AtomId,
              position::AbstractVector{L},
              velocity::AbstractVector{V}=zeros(3)u"bohr/s";
              atomic_symbol=Symbol(elements[identifier].symbol),
              atomic_number=elements[identifier].number,
              atomic_mass::M=elements[identifier].atomic_mass,
              kwargs...) where {L <: Unitful.Length, V <: Unitful.Velocity, M <: Unitful.Mass}
    Atom{length(position), L, V, M}(position, velocity, atomic_symbol,
                                    atomic_number, atomic_mass, Dict(kwargs...))
end

# Update constructor: Amend any atom by extra data.
function Atom(;atom, kwargs...)
    extra = atom isa Atom ? atom.data : (; )
    Atom(atomic_symbol(atom), position(atom), velocity(atom);
         atomic_symbol=atomic_symbol(atom),
         atomic_number=atomic_number(atom),
         atomic_mass=atomic_mass(atom),
         extra..., kwargs...)
end
Atom(atom::Atom; kwargs...) = Atom(; atom, kwargs...)

function Base.convert(::Type{Atom}, id_pos::Pair{<:AtomId,<:AbstractVector{<:Unitful.Length}})
    Atom(id_pos...)
end

function Base.show(io::IO, at::Atom{D, L}) where {D, L}
    pos  = ustrip.(at.position)
    print(io, "Atom($(at.atomic_symbol), [", join(pos, ", "), "]u\"$(unit(L))\")")
end

#
# Special high-level functions to construct atomic systems
#
atomic_system(atoms::AbstractVector{<:Atom}, box, bcs) = FlexibleSystem(atoms, box, bcs)
atomic_system(atoms::AbstractVector, box, bcs) = FlexibleSystem(convert.(Atom, atoms), box, bcs)

function isolated_system(atoms::AbstractVector{<:Atom})
    # Use dummy box and boundary conditions
    D = n_dimensions(first(atoms))
    atomic_system(atoms, infinite_box(D), fill(DirichletZero(), D))
end
isolated_system(atoms::AbstractVector) = isolated_system(convert.(Atom, atoms))

function periodic_system(atoms::AbstractVector,
                         box::AbstractVector{<:AbstractVector},
                         boundary_conditions=fill(Periodic(), length(box));
                         fractional=false, kwargs...)
    lattice = hcat(box...)
    !fractional && return atomic_system(atoms, box, boundary_conditions)

    parse_fractional(atom::Atom) = atom
    function parse_fractional(atom::Pair)::Atom
        id, pos_fractional = atom
        Atom(id, lattice * pos_fractional)
    end
    atomic_system(parse_fractional.(atoms), box, boundary_conditions)
end


text:
using JuMP, Cbc, NPZ

function Solve_func(year, tol)

	println("Solving weights for $year ...\n\n")

	array = npzread(string(year, "_input.npz"))

	A1 = array["A1"]
	A2 = array["A2"]
	b = array["b"]

	model = Model(Cbc.Optimizer)
	set_optimizer_attribute(model, "logLevel", 1)
	N = size(A1)[2]

	@variable(model, r[1:N] >= 0)
	@variable(model, s[1:N] >= 0)

	@objective(model, Min, sum(r[i] + s[i] for i in 1:N))

	# bound on top by tolerance
	@constraint(model, [i in 1:N], r[i] + s[i] <= tol)

	# Ax = b
	@constraint(model, [i in 1:length(b)], sum(A1[i,j] * r[j] + A2[i,j] * s[j]
		                          for j in 1:N) == b[i])

	optimize!(model)
	termination_status(model)

	r_vec = value.(r)
	s_vec = value.(s)

	npzwrite(string(year, "_output.npz"), Dict("r" => r_vec, "s" => s_vec))

	println("\n")

end

year_list = [x for x in 2012:2030]
tol_list = [0.40, 0.38, 0.35, 0.33, 0.30, 0.37, 0.38,
			0.38, 0.39, 0.39, 0.38, 0.40, 0.39, 0.41,
			0.41, 0.42, 0.42, 0.42, 0.42]

# Run solver function for all years and tolerances (in order)
for i in zip(year_list, tol_list)
	Solve_func(i[1], i[2])
end

max_stars_repo_name: chusloj/taxdata

max_issues_repo_name: chusloj/taxdata

max_forks_repo_name: chusloj/taxdata


text:
#for undef and occursin
using Compat
using Test

#examples_path = joinpath(dirname(dirname(@__FILE__)), "examples")
#for example in ["simple_mopec2.jl"]
#    include(joinpath(examples_path, example))
#end

# This doesn't appears to be in Compat?
@static if VERSION >= v"0.7.0"
    function isdef(symbol)
        return isdefined(@__MODULE__, symbol)
    end
else
    function isdef(symbol)
        return isdefined(symbol)
    end
end

# On the new macros are really robust
if VERSION >= v"1.0.0"
    include("variable.jl")
end

include("simple_mopec2.jl")
include("gnep_river_basin.jl")
include("ovf_loss_nl.jl")
include("ovf_loss_equil.jl")

# write your own tests here


text:
module TestUtils

using Test
using LightGBM

X_train = randn(1000, 20)
y_train_binary = rand(0:1, 1000)
y_train_regression = rand(1000)

@testset "stringifyparams -- convert to zero-based" begin
    indices = [1, 3, 5, 7, 9]
    classifier = LightGBM.LGBMClassification(categorical_feature = indices)
    ds_parameters = LightGBM.stringifyparams(classifier, LightGBM.DATASETPARAMS)

    expected = "categorical_feature=0,2,4,6,8"
    @test occursin(expected, ds_parameters)
end

@testset "loadmodel predicts same as original model -- regression" begin
    # Arrange
    estimator = LightGBM.LGBMRegression(objective = "regression")
    LightGBM.fit!(estimator, X_train, y_train_regression)
    expected_prediction = predict(estimator, X_train)
    # Save the fitted model.
    model_filename = joinpath(@__DIR__, "fixture.model")
    savemodel(estimator, model_filename)

    # Act
    estimator_from_file = LGBMRegression()
    loadmodel!(estimator_from_file, model_filename)
    actual_prediction = predict(estimator_from_file, X_train)

    # Assert
    @test expected_prediction == actual_prediction

    # Teardown
    rm(model_filename)

end

@testset "loadmodel predicts same as original model -- binary" begin
    # Arrange
    estimator = LightGBM.LGBMClassification(objective = "binary", num_class = 1)
    LightGBM.fit!(estimator, X_train, y_train_binary)
    expected_prediction = predict(estimator, X_train)
    # Save the fitted model.
    model_filename = joinpath(@__DIR__, "fixture.model")
    savemodel(estimator, model_filename)

    # Act
    estimator_from_file = LGBMClassification()
    loadmodel!(estimator_from_file, model_filename)
    actual_prediction = predict(estimator_from_file, X_train)

    # Assert
    @test expected_prediction == actual_prediction

    # Teardown
    rm(model_filename)
end

@testset "loadmodel predicts same as original model with custom pararms -- regression" begin
    # Arrange
    estimator = estimator = LGBMRegression(
        objective = "regression",
        num_iterations = 100,
        learning_rate = .2,
        early_stopping_round = 3,
        feature_fraction = .5,
        bagging_fraction = .6,
        bagging_freq = 2,
        num_leaves = 100,
        metric = ["auc", "binary_logloss"]
    )
    LightGBM.fit!(estimator, X_train, y_train_regression)
    expected_prediction = predict(estimator, X_train)
    # Save the fitted model.
    model_filename = joinpath(@__DIR__, "fixture.model")
    savemodel(estimator, model_filename)

    # Act
    estimator_from_file = LGBMRegression()
    loadmodel!(estimator_from_file, model_filename)
    actual_prediction = predict(estimator_from_file, X_train)

    # Assert
    @test expected_prediction == actual_prediction

    # Teardown
    rm(model_filename)

end

@testset "loadmodel predicts same as original model with custom params -- binary" begin
    # Arrange
    estimator = estimator = LGBMClassification(
        objective = "binary",
        num_iterations = 100,
        learning_rate = .2,
        early_stopping_round = 3,
        feature_fraction = .5,
        bagging_fraction = .6,
        bagging_freq = 2,
        num_leaves = 100,
        num_class = 1,
        metric = ["auc", "binary_logloss"]
    )
    LightGBM.fit!(estimator, X_train, y_train_binary)
    expected_prediction = predict(estimator, X_train)
    # Save the fitted model.
    model_filename = joinpath(@__DIR__, "fixture.model")
    savemodel(estimator, model_filename)

    # Act
    estimator_from_file = LGBMClassification()
    loadmodel!(estimator_from_file, model_filename)
    actual_prediction = predict(estimator_from_file, X_train)

    # Assert
    @test expected_prediction == actual_prediction

    # Teardown
    rm(model_filename)
end

end # Module


text:
using ArgParse2
using DataFrames
using DataStructures
using BioSequences
using CodecZlib
using BioAlignments
using DataFramesMeta
using Statistics
using st16SseqJuliaTools
using ProgressMeter
using XAM

struct Containment
    template:: String #id of the first sequences
    reference:: String #id of other sequence
    template_from::Int64
    template_to::Int64
    identity:: Float64 #Aligned part identuty
    coverage:: Float64 #Coverage fraction in trms of template length
    tlength:: Int64 #template length
    rlength:: Int64 #reference length
    endfragid_shorter::Float64
    endfragid_longer::Float64
    aln_fulllen::PairwiseAlignment{LongSequence{DNAAlphabet{4}},LongSequence{DNAAlphabet{4}}}
    aln_3end_shorter::PairwiseAlignment{LongSequence{DNAAlphabet{4}},LongSequence{DNAAlphabet{4}}}
    aln_3end_longer::PairwiseAlignment{LongSequence{DNAAlphabet{4}},LongSequence{DNAAlphabet{4}}}
end

function main()
    parser = ArgumentParser(prog = "analyse_alnself.jl",
                        description = "Detect chimeras",
                        epilog = "Removing potential chimeras and repair snps")
    add_argument!(parser, "--ref", "-r", help = "reference fasta",type = String)
    add_argument!(parser, "--salmon-rez", "-s", help = "Salmon rezults",type = String)
    add_argument!(parser, "--out", "-o", help = "Fasta file to write out cleaned up sequences",default="chimerasfree.fasta",type = String)
    add_argument!(parser, "--fused", "-f", help = "Fasta file to write out cleaned up sequences where contigs of each cluster are joined in one sequence",default="chimerasfree.fasta",type = String)

    args = parse_args(parser)
    salmon_data = read_salmon(args.salmon_rez)
    ref = fasta_to_dict(args.ref)
    length_data=Dict{String,Int64}()
    for k in keys(ref)
        length_data[k] = length(ref[k])
    end

    #Limits for real match ... maybe should be options
    endlength1 = 15 #15 #length from end to rtest separately
    endlength2 = 50 #50 #length from end to rtest separately
    endfrag_id_limit = 0.8 #if sequences match general cut-off but the end fragment is dissimilar they are not merged
    contained_idf_limit_c = 0.96 # 0.96 Limid for identity for containment 0.97 #z1 - 7_416;size=1", "4_2136;size=1. Salmonela and coli
    contained_covf_limit_c = 0.99 # Limid for coverage for containment

    costmodel = BioAlignments.CostModel(match=0, mismatch=1, insertion=1, deletion=1);
    #known_chimeras = ["1_1868;size=92", "1_6093;size=54", "1_2188;size=22","1_1760;size=8"]
    pairs_4contain = Array{Containment,1}()
    aln_costs = AffineGapScoreModel(
        match=5,
        mismatch=-4,
        gap_open=-10,
        gap_extend=-5
    )
    aln_costs_wogaps = AffineGapScoreModel(
        match=5,
        mismatch=-4,
        gap_open=-10,
        gap_extend=-5
    )

    ref_keys = collect(keys(ref))
    @showprogress 1 "Pairwise alignments..."  for k1 in ref_keys
        #TO DO - below there are three consequtive alignments with some similar steps - should be a function there
        pairs_part_4contain = Array{Containment,1}(undef,length(ref_keys))
        s1 = split(k1,"@")[end]
        Threads.@threads for i in 1:length(ref_keys)
            k2 = ref_keys[i]
            s2 = split(k2,"@")[end]
            #Limit containments analysis to the same species
            if (k1 != k2) && (s1 == s2)
                rn = k1
                tn = tempname = k2
                rs = ref[rn]
                ts = ref[tn]
                tlen = length(ts)
                rlen = length(rs)

                # Calculate overall alignment

                aln_fulllen = alignment(pairalign(LocalAlignment(), rs, ts, aln_costs))
                matches_fulllen = BioAlignments.count_matches(aln_fulllen)
                mismatches_fulllen = BioAlignments.count_mismatches(aln_fulllen)
                ins_fulllen = count_insertions(aln_fulllen)
                del_fulllen = count_deletions(aln_fulllen)
                s_fulllen = pos_fulllen = first(aln_fulllen.a.aln.anchors).seqpos
                e_fulllen = rpos_fulllen = last(aln_fulllen.a.aln.anchors).seqpos
                t_pos_fulllen=first(aln_fulllen.a.aln.anchors).refpos+1
                t_rpos_fulllen=last(aln_fulllen.a.aln.anchors).refpos
                taln_fulllen = abs(t_pos_fulllen-t_rpos_fulllen)+1
                idfrac_fulllen = matches_fulllen/taln_fulllen
                #calculate identity at 3'
                tendfrag1 = ts[t_rpos_fulllen-endlength1+1:t_rpos_fulllen]
                rendfrag1 = rs[e_fulllen-endlength1+1:e_fulllen]
                tendfrag2 = ts[t_rpos_fulllen-endlength2+1:t_rpos_fulllen]
                rendfrag2 = rs[e_fulllen-endlength2+1:e_fulllen]
                aln_3end1 = alignment(pairalign(GlobalAlignment(), rendfrag1, tendfrag1, aln_costs_wogaps))
                aln_3end2 = alignment(pairalign(GlobalAlignment(), rendfrag2, tendfrag2, aln_costs_wogaps))
                matches_3end1 = BioAlignments.count_matches(aln_3end1)
                matches_3end2 = BioAlignments.count_matches(aln_3end2)

                t_pos_3end1=first(aln_3end1.a.aln.anchors).refpos+1
                t_rpos_3end1=last(aln_3end1.a.aln.anchors).refpos
                taln_3end1 = abs(t_pos_3end1-t_rpos_3end1)+1
                idfrac_3end1 = matches_3end1/taln_3end1

                t_pos_3end2=first(aln_3end2.a.aln.anchors).refpos+1
                t_rpos_3end2=last(aln_3end2.a.aln.anchors).refpos
                taln_3end2 = abs(t_pos_3end2-t_rpos_3end2)+1
                idfrac_3end2 = matches_3end2/taln_3end2

                if taln_fulllen > 0
                    pairs_part_4contain[i] = Containment(tempname,
                    rn,
                    t_pos_fulllen,
                    t_rpos_fulllen,
                    idfrac_fulllen,
                    taln_fulllen/length(ts),
                    length(ts),
                    length(rs),
                    idfrac_3end1,
                    idfrac_3end2,
                    aln_fulllen,
                    aln_3end1,
                    aln_3end2)
                end
            end
        end
        for ii in 1:length(ref_keys)
            if isassigned(pairs_part_4contain,ii)
                push!(pairs_4contain,pairs_part_4contain[ii])
            end
        end
    end

    #Collect data for  contained detection

    contained = Array{String,1}()
    #reformat collected data into dictionary
    tn_rn_map = Dict{String,Array{Containment,1}}()
    for pair in pairs_4contain
        #tn,rn,idf,covf,rl =  pair[1],pair[2],pair[5],pair[6],pair[7]
        if !(pair.template in keys(tn_rn_map))
            tn_rn_map[pair.template] = Array{Containment,1}()
        end
        push!(tn_rn_map[pair.template],pair)
    end

    #detect contained

    for tn in keys(tn_rn_map)
        tl = length(ref[tn])
        cont = false
        for rd in tn_rn_map[tn]
            t_s = salmon_data[split(tn,";")[1]]
            r_s = salmon_data[split(rd.reference,";")[1]]
            t_clnr = split(tn,"_")[1]
            r_clnr = split(rd.reference,"_")[1]
            tl_shorter_fraction = abs(rd.rlength-tl)/tl
            if tl <= rd.rlength #&& (t_clnr == r_clnr)
                if (round(rd.identity,digits=2) >= contained_idf_limit_c &&
                    round(rd.coverage,digits=2) >= contained_covf_limit_c &&
                    ( ( (round(rd.endfragid_longer,digits=2) >= endfrag_id_limit) && (round(rd.endfragid_shorter,digits=2) >= endfrag_id_limit)  ) || (tl_shorter_fraction > 0.2) ) )
                    if tl == rd.rlength
                        #if lengths matches add abundance criteria
                        if r_s > t_s
                            cont = true
                        end
                    else
                        cont = true
                    end
                end
            end
        end
        if cont
            push!(contained,tn)
        end
    end
    filtered_out = OrderedDict{String,LongSequence{BioSequences.DNAAlphabet{4}}}()
    ct_contained = 0
    ct_good = 0
    for k in keys(ref)
        iscont = (k in contained)
        if  !iscont
            ct_good += 1
            filtered_out[k]= ref[k]
        else
            if iscont
                ct_contained += 1
            end
        end
    end
    dict_to_fasta(filtered_out,args.out)
    #Join from the same cluster
    fused_seqs = OrderedDict{String,Array{String,1}}()
    fused_ids = OrderedDict{String,Array{String,1}}()
    for k in keys(filtered_out)
        cid = split(k,"_")[1]
        if !(cid in keys(fused_seqs))
            fused_seqs[cid] = Array{String,1}()
            fused_ids[cid] = Array{String,1}()
        end
        push!(fused_seqs[cid],String(filtered_out[k]))
        push!(fused_ids[cid],k)
    end
    fused_seqnames = OrderedDict{String,LongSequence{DNAAlphabet{4}}}()
    spacer = String(fill('N',300))
    for k in keys(fused_ids)
        seqname=join(fused_ids[k],";")
        seqname = replace(seqname,";size=1@"=>"|")
        seqname = replace(seqname,k*"_"=>"")
        prefix = k*"@"*string(length(fused_ids[k]))*":"
        fused_seqnames[prefix*seqname] = LongSequence{DNAAlphabet{4}}(join(fused_seqs[k],spacer))
    end
    dict_to_fasta(fused_seqnames,args.fused)

    println(stderr,"$ct_contained contained sequences were detected, $ct_good sequences were  written to ",args.out)
    println(stderr,"Fused sequences are written to ",args.fused, ", name format: <clusterid>@<number of contigs>: ....<contig id|species based on blast>...separated by ;")
end

main()