text:
module CudaUtilsModule

import CUDAdrv
import CUDAdrv: CuModule, CuModuleFile, unsafe_unload!,
                CuFunction, cudacall, CuDevice, CuDim

const ptxdict = Dict()
const mdlist = CuModule[]

function mdinit(devlist, ptxfile)
  global ptxdict
  global mdlist

  # isempty(mdlist) || error("mdlist is not empty")
  for dev in devlist
    # device(dev)
    CuDevice(dev)
    md = CuModuleFile(ptxfile)

    # Utils functions
    ptxdict["setup_kernel"] = CuFunction(md, "setup_kernel")
    ptxdict["perturb"]      = CuFunction(md, "perturb")
    ptxdict["veccost2"]     = CuFunction(md, "veccost2")

    # Add vector to matrix
    ptxdict["vec_add"]         = CuFunction(md, "vec_add")

    # Forward pass of the viterbi algorithm
    ptxdict["viterbi_forward"] = CuFunction(md, "viterbi_forward")

    # ICM functions
    ptxdict["condition_icm3"] = CuFunction(md, "condition_icm3")

    push!(mdlist, md)
  end
end

# mdclose() = (for md in mdlist; unsafe_unload!(md); end; empty!(mdlist); empty!(ptxdict))
mdclose() = (empty!(mdlist); empty!(ptxdict))

function finit()
  mdclose()
end

function init(devlist, ptxfile)
  mdinit(devlist, ptxfile)
end

# CUDAdrv.Mem.Buffer
function condition_icm3(
  nblocks::Integer, nthreads::CuDim,
  d_ub::CUDAdrv.Mem.Buffer,     # in/out. h-by-n. where we accumulate
  d_bbs::CUDAdrv.Mem.Buffer,     # in. binaries
  d_codek::CUDAdrv.Mem.Buffer,  # in. n-long. indices into d_bb
  conditioning::Cint,
  m::Cint,
  n::Cint)                        # in. size( d_ub, 2 )

  fun = ptxdict["condition_icm3"]
  cudacall(fun,
    (Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cuchar}, Cint, Cint, Cint),
    d_ub, d_bbs, d_codek, conditioning, m, n; blocks=nblocks, threads=nthreads)

  return nothing
end

function vec_add(
  nblocks::CuDim, nthreads::CuDim,
  d_matrix::CUDAdrv.Mem.Buffer,
  d_vec::CUDAdrv.Mem.Buffer,
  n::Cint,
  h::Cint)

  fun = ptxdict["vec_add"]
  cudacall(fun,
    Tuple{Ptr{Cfloat}, Ptr{Cfloat}, Cint, Cint},
    d_matrix, d_vec, n, h; blocks=nblocks, threads=nthreads)
  end

function viterbi_forward(
  nblocks::CuDim, nthreads::CuDim,
  d_matrix::CUDAdrv.Mem.Buffer,
  d_vec::CUDAdrv.Mem.Buffer,
  d_minv::CUDAdrv.Mem.Buffer,
  d_mini::CUDAdrv.Mem.Buffer,
  n::Cint, j::Cint)

  fun = ptxdict["viterbi_forward"]
  cudacall(fun,
    (Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cint}, Cint, Cint),
    d_matrix, d_vec, d_minv, d_mini, n, j; blocks=nblocks, threads=nthreads)
end

function veccost2(
  nblocks::CuDim, nthreads::CuDim,
  d_rx::CUDAdrv.Mem.Buffer,
  d_codebooks::CUDAdrv.Mem.Buffer,
  d_codes::CUDAdrv.Mem.Buffer,
  d_veccost::CUDAdrv.Mem.Buffer, # out.
  d::Cint,
  m::Cint,
  n::Cint)

  fun = ptxdict["veccost2"]
  cudacall(fun,
    (Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cuchar}, Ptr{Cfloat}, Cint, Cint, Cint),
    d_rx, d_codebooks, d_codes, d_veccost, d, m, n; blocks=nblocks, threads=nthreads,
    shmem=Int(d*sizeof(Cfloat)))
end

function perturb(
  nblocks::CuDim, nthreads::CuDim,
  state::CUDAdrv.Mem.Buffer,
  codes::CUDAdrv.Mem.Buffer,
  n::Cint,
  m::Cint,
  k::Cint)

  fun = ptxdict["perturb"]
  cudacall(fun,
    (Ptr{Nothing}, Ptr{Cuchar}, Cint, Cint, Cint),
    state, codes, n, m, k; blocks=nblocks, threads=nthreads)
end

function setup_kernel(
  nblocks::CuDim, nthreads::CuDim,
  n::Cint,
  state::CUDAdrv.Mem.Buffer)

  fun = ptxdict["setup_kernel"]
  cudacall(fun,
    (Cint, Ptr{Nothing}),
    n, state; blocks=nblocks, threads=nthreads)
end

end #cudaUtilsModule


text:
function compress_system(	particle_type::String,
						R::Array{Float64, 2},
						Lx::Float64,
						Ly::Float64,
						Lz::Float64,
						phi_target::Float64,
						X::Array{Float64, 1},
						Y::Array{Float64, 1},
						Z::Array{Float64, 1},
						Q0::Array{Float64, 1},
						Q1::Array{Float64, 1},
						Q2::Array{Float64, 1},
						Q3::Array{Float64, 1},
						A11::Array{Float64, 1},
						A12::Array{Float64, 1},
						A13::Array{Float64, 1},
						A21::Array{Float64, 1},
						A22::Array{Float64, 1},
						A23::Array{Float64, 1},
						A31::Array{Float64, 1},
						A32::Array{Float64, 1},
						A33::Array{Float64, 1},
						position_constraint_axis::String,
						position_constraint_lower::Float64,
						position_constraint_upper::Float64,
						orientation_axis::Array{Float64, 1},
						orientation_constraint_axis::Array{Float64, 1},
						orientation_constraint_lower::Float64,
						orientation_constraint_upper::Float64,
						sigma_translation::Float64,
						sigma_translation_max::Float64,
						sigma_rotation::Float64,
						sigma_rotation_max::Float64,
						sigma_ratio::Float64,
						delta_phi::Float64,
						number_of_relaxation_sweeps_max::Int64)

	number_of_particles::Int64 = size(R, 1)

	phi_initial::Float64 = 0.0
	if particle_type == "sphere"
		phi_initial = sum(4.0 * pi / 3.0 * R[:, 1].^3) / (Lx * Ly * Lz)
	elseif particle_type == "ellipsoid"
		phi_initial = sum(4.0 * pi / 3.0 * R[:, 1] .* R[:, 2] .* R[:, 3]) / (Lx * Ly * Lz)
	elseif particle_type == "cuboid"
		phi_initial = sum(8.0 * R[:, 1] .* R[:, 2] .* R[:, 3]) / (Lx * Ly * Lz)
	end

	is_converged::Bool = false
	is_relaxed::Bool = false
	phi::Float64 = phi_initial
	phi_prim::Float64 = 0.0
	Lx_prim::Float64 = 0.0
	Ly_prim::Float64 = 0.0
	Lz_prim::Float64 = 0.0
	current_compression::Int64 = 0

	X_prim::Array{Float64, 1} = zeros(number_of_particles)
	Y_prim::Array{Float64, 1} = zeros(number_of_particles)
	Z_prim::Array{Float64, 1} = zeros(number_of_particles)
	Q0_prim::Array{Float64, 1} = zeros(number_of_particles)
	Q1_prim::Array{Float64, 1} = zeros(number_of_particles)
	Q2_prim::Array{Float64, 1} = zeros(number_of_particles)
	Q3_prim::Array{Float64, 1} = zeros(number_of_particles)
	A11_prim::Array{Float64, 1} = zeros(number_of_particles)
	A12_prim::Array{Float64, 1} = zeros(number_of_particles)
	A13_prim::Array{Float64, 1} = zeros(number_of_particles)
	A21_prim::Array{Float64, 1} = zeros(number_of_particles)
	A22_prim::Array{Float64, 1} = zeros(number_of_particles)
	A23_prim::Array{Float64, 1} = zeros(number_of_particles)
	A31_prim::Array{Float64, 1} = zeros(number_of_particles)
	A32_prim::Array{Float64, 1} = zeros(number_of_particles)
	A33_prim::Array{Float64, 1} = zeros(number_of_particles)

	while !is_converged && phi < phi_target
		current_compression += 1
		phi_prim = phi_initial + convert(Float64, current_compression) * delta_phi

		Lx_prim = (phi / phi_prim)^(1/3) * Lx
		Ly_prim = (phi / phi_prim)^(1/3) * Ly
		Lz_prim = (phi / phi_prim)^(1/3) * Lz

		for current_particle = 1:number_of_particles
			X_prim[current_particle] = Lx_prim / Lx * X[current_particle]
			Y_prim[current_particle] = Ly_prim / Ly * Y[current_particle]
			Z_prim[current_particle] = Lz_prim / Lz * Z[current_particle]

			Q0_prim[current_particle] = Q0[current_particle]
			Q1_prim[current_particle] = Q1[current_particle]
			Q2_prim[current_particle] = Q2[current_particle]
			Q3_prim[current_particle] = Q3[current_particle]
			A11_prim[current_particle] = A11[current_particle]
			A12_prim[current_particle] = A12[current_particle]
			A13_prim[current_particle] = A13[current_particle]
			A21_prim[current_particle] = A21[current_particle]
			A22_prim[current_particle] = A22[current_particle]
			A23_prim[current_particle] = A23[current_particle]
			A31_prim[current_particle] = A31[current_particle]
			A32_prim[current_particle] = A32[current_particle]
			A33_prim[current_particle] = A33[current_particle]
		end

		(	X_prim,
			Y_prim,
			Z_prim,
			Q0_prim,
			Q1_prim,
			Q2_prim,
			Q3_prim,
			A11_prim,
			A12_prim,
			A13_prim,
			A21_prim,
			A22_prim,
			A23_prim,
			A31_prim,
			A32_prim,
			A33_prim,
			sigma_translation,
			sigma_rotation,
			is_relaxed) = evolve_system(particle_type,
										R,
										Lx_prim,
										Ly_prim,
										Lz_prim,
										X_prim,
										Y_prim,
										Z_prim,
										Q0_prim,
										Q1_prim,
										Q2_prim,
										Q3_prim,
										A11_prim,
										A12_prim,
										A13_prim,
										A21_prim,
										A22_prim,
										A23_prim,
										A31_prim,
										A32_prim,
										A33_prim,
										position_constraint_axis,
										position_constraint_lower,
										position_constraint_upper,
										orientation_axis,
										orientation_constraint_axis,
										orientation_constraint_lower,
										orientation_constraint_upper,
										sigma_translation,
										sigma_translation_max,
										sigma_rotation,
										sigma_rotation_max,
										sigma_ratio,
										number_of_relaxation_sweeps_max,
										0)

		if is_relaxed # If system has been successfully relaxed.
			phi = phi_prim
			Lx = Lx_prim
			Ly = Ly_prim
			Lz = Lz_prim

			for current_particle = 1:number_of_particles
				X[current_particle] = X_prim[current_particle]
				Y[current_particle] = Y_prim[current_particle]
				Z[current_particle] = Z_prim[current_particle]

				Q0[current_particle] = Q0_prim[current_particle]
				Q1[current_particle] = Q1_prim[current_particle]
				Q2[current_particle] = Q2_prim[current_particle]
				Q3[current_particle] = Q3_prim[current_particle]
				A11[current_particle] = A11_prim[current_particle]
				A12[current_particle] = A12_prim[current_particle]
				A13[current_particle] = A13_prim[current_particle]
				A21[current_particle] = A21_prim[current_particle]
				A22[current_particle] = A22_prim[current_particle]
				A23[current_particle] = A23_prim[current_particle]
				A31[current_particle] = A31_prim[current_particle]
				A32[current_particle] = A32_prim[current_particle]
				A33[current_particle] = A33_prim[current_particle]
			end

			println(join(("   Current volume fraction: ", string(round(phi, convert(Int64, ceil(-log10(delta_phi))))))))
		else
			is_converged = true
		end
	end

	# Due to rounding errors, it happens that phi > phi_target. Compensate by
	# slightly expanding the system.
	if phi > phi_target
		linear_expansion_factor::Float64 = (phi / phi_target)^(1/3)
		phi = phi_target
		Lx *= linear_expansion_factor
		Ly *= linear_expansion_factor
		Lz *= linear_expansion_factor
		X *= linear_expansion_factor
		Y *= linear_expansion_factor
		Z *= linear_expansion_factor
	end

	return (Lx, Ly, Lz, phi, X, Y, Z, Q0, Q1, Q2, Q3, A11, A12, A13, A21, A22, A23, A31, A32, A33, sigma_translation, sigma_rotation)
end

ext: jl

max_stars_repo_name: roding/whitefish

max_issues_repo_name: roding/whitefish

max_forks_repo_name: roding/whitefish


text:
# using CPLEX
# using GLPK
# using Cbc
using Clp
using JuMP

include("seqata_test_model_util.jl")

# ===========
# Préparartion des arguments
args = Args.parse_commandline()
const TEST_LEVEL = 1
Args.set(:level, TEST_LEVEL) # sera passé temporairement à 0
Args.set(:outdir, "/tmp") # eqivalent à  -d/tmp
Args.set(:itermax, 300)

# ===========
inst = instance_build_mini10()
@test inst.nb_planes == 10

# Création d'une struct encapsulant la spécifiation d'un test LP
specs = [
   LpSolverSpec(:lp,  :clp, LpTimingSolver,  Clp.Optimizer),
#    LpSolverSpec(:lp2, :clp, Lp2TimingSolver, Clp.Optimizer),
#    LpSolverSpec(:lp3, :clp, Lp3TimingSolver, Clp.Optimizer),
#    LpSolverSpec(:lp4, :clp, Lp4TimingSolver, Clp.Optimizer),
]

for spec in specs
    test_name = "Test annealing LP $(spec.algo):$(spec.external_mip_solver)"
    lg1("$test_name ... ")
    @testset "$test_name" begin
        Args.set(:level, 0)
        test_one_lp_annealing(spec)
        Args.set(:level, TEST_LEVEL)
    end
    ln1(" => fait.")
end


text:
using Distributed
addprocs(1)
@everywhere using BLG_SQUID
# using MKL
# BLAS.set_num_threads(2)
@everywhere using MKL, SharedArrays
@everywhere using Quantica, StaticArrays, Parameters, LinearAlgebra, StatsBase
@everywhere using Baselet, Arpack, StaticArrays, CSV, DataFrames, Dates
@everywhere using Optim, ProgressMeter, Distributed
@everywhere using PhysicalConstants.CODATA2018: ustrip, @u_str, ħ, k_B, m_e, e, μ_B

using BLG_SQUID, Quantica, StaticArrays, Parameters, LinearAlgebra, StatsBase, SharedArrays
using Baselet, Arpack, StaticArrays, CSV, DataFrames, Dates
using PhysicalConstants.CODATA2018: ustrip, @u_str, ħ, k_B, m_e, e, μ_B
using Optim, ProgressMeter

using ElectronDisplay, VegaLite, CairoMakie

max_stars_repo_path: src/using Distributed.jl

max_issues_repo_path: src/using Distributed.jl

max_forks_repo_path: src/using Distributed.jl


text:
using NamedDims
using NamedDims:
    names,
    unify_names,
    unify_names_longest,
    unify_names_shortest,
    dim_noerror,
    tuple_issubset,
    tuple_cat,
    order_named_inds,
    permute_dimnames,
    remaining_dimnames_from_indexing,
    remaining_dimnames_after_dropping
using Test

@testset "dim" begin
    @testset "small case" begin
        @test dim((:x, :y), :x) == 1
        @test dim((:x, :y), :y) == 2
        @test_throws ArgumentError dim((:x, :y), :z)  # not found
    end
    @testset "large case" begin
        @test dim((:x, :y, :a, :b, :c, :d), :x) == 1
        @test dim((:x, :y, :a, :b, :c, :d), :a) == 3
        @test dim((:x, :y, :a, :b, :c, :d), :d) == 6
        @test_throws ArgumentError dim((:x, :y, :a, :b, :c, :d), :z) # not found
    end
end
@testset "allocations: dim" begin
    @test 0 == @allocated (()->dim((:a, :b), :b))()
    @test 0 == @allocated (()->dim_noerror((:a, :b, :c), :c))()
    if VERSION >= v"1.1"
        @test 0 == @allocated (()->dim((:a,:b), (:a,:b)))()
    else
        @test_broken 0 == @allocated (()->dim((:a,:b), (:a,:b)))()
    end
end

@testset "unify_names_*" begin
    @test_throws DimensionMismatch unify_names((:a,), (:a, :b,))

    @test unify_names_longest((:a,), (:a, :b,)) == (:a, :b)
    @test unify_names_longest((:a,), (:a, :_)) == (:a, :_)
    @test unify_names_longest((:a, :b), (:a, :_, :c)) == (:a, :b, :c)

    @test unify_names_shortest((:a,), (:a, :b,)) == (:a,)
    @test unify_names_shortest((:a,), (:a, :_)) == (:a,)
    @test unify_names_shortest((:a, :b), (:a, :_, :c)) == (:a, :b)

    @test_throws DimensionMismatch unify_names_longest((:a, :b, :c), (:b, :a))

    for unify in (unify_names, unify_names_longest, unify_names_shortest)
        @test unify((:a,), (:a,)) == (:a,)
        @test unify((:a, :b), (:a, :b)) == (:a, :b)
        @test unify((:a, :_), (:a, :b)) == (:a, :b)
        @test unify((:a, :_), (:a, :_)) == (:a, :_)

        @test unify((:a, :b, :c), (:_, :_, :_)) == (:a, :b, :c)
        @test unify((:a, :_, :c), (:_, :b, :_)) == (:a, :b, :c)
        @test unify((:_, :_, :_), (:_, :_, :_)) == (:_, :_, :_)

        @test_throws DimensionMismatch unify((:a,), (:b,))
        @test_throws DimensionMismatch unify((:a,:b), (:b, :a))
        @test_throws DimensionMismatch unify((:a, :b, :c), (:_, :_, :d))
    end

    # vararg version
    @test unify_names((:a, :_), (:a, :b,), (:_, :b)) == (:a, :b)
    @test unify_names((:a, :b,)) == (:a, :b)
    @test_throws DimensionMismatch unify_names((:a, :_), (:a, :b,), (:_, :c))
end
@testset "allocations: unify_names_*" begin
    for unify in (unify_names, unify_names_longest, unify_names_shortest)
        if VERSION >= v"1.2"
            @test 0 == @allocated (()->unify((:a, :b), (:a, :_)))()
        else
            @test_broken 0 == @allocated (()->unify((:a, :b), (:a, :_)))()
        end
    end
    if VERSION >= v"1.1"
        @test 0 == @allocated (()->unify_names_longest((:a, :b), (:a, :_, :c)))()
        @test 0 == @allocated (()->unify_names_shortest((:a, :b), (:a, :_, :c)))()
        @test 0 == @allocated (()->unify_names((:a, :b), (:a, :_), (:_, :b)))()
    else
        @test_broken 0 == @allocated (()->unify_names_longest((:a, :b), (:a, :_, :c)))()
        @test_broken 0 == @allocated (()->unify_names_shortest((:a, :b), (:a, :_, :c)))()
        @test_broken 0 == @allocated (()->unify_names((:a, :b), (:a, :_), (:_, :b)))()
    end

end

@testset "order_named_inds" begin
    @test order_named_inds((:x,)) == (:,)
    @test order_named_inds((:x,); x=2) == (2,)

    @test order_named_inds((:x, :y,)) == (:, :)
    @test order_named_inds((:x, :y); x=2) == (2, :)
    @test order_named_inds((:x, :y); y=2, ) == (:, 2)
    @test order_named_inds((:x, :y); y=20, x=30) == (30, 20)
    @test order_named_inds((:x, :y); x=30, y=20) == (30, 20)
end
@testset "allocations: order_named_inds" begin
    if v"1.1-" <= VERSION <= v"1.2-" # test passes on 1.1, including 1.1.1
        @test 0 == @allocated (()->order_named_inds((:a, :b, :c); b=1, c=2))()
        @test 0 == @allocated (()->order_named_inds((:a, :b, :c), (b=1, c=2)))()
    else
        @test_broken 0 == @allocated (()->order_named_inds((:a, :b, :c); b=1, c=2))()
        @test_broken 0 == @allocated (()->order_named_inds((:a, :b, :c), (b=1, c=2)))()
    end
end

@testset "remaining_dimnames_from_indexing" begin
    @test remaining_dimnames_from_indexing((:a, :b, :c), (10,20,30)) == tuple()
    @test remaining_dimnames_from_indexing((:a, :b, :c), (10, :,30)) == (:b,)
    @test remaining_dimnames_from_indexing((:a, :b, :c), (1:1, [true], [20])) == (:a, :b, :c)

    @test remaining_dimnames_from_indexing((:a, :b, :c), (:, [CartesianIndex()], :, :)) == (:a, :_, :b, :c)
    @test remaining_dimnames_from_indexing((:a, :b, :c), (1, [CartesianIndex()], 2, :)) == (:_, :c)
    @test remaining_dimnames_from_indexing((:a, :b, :c), (CartesianIndex(1,1), :)) == (:c,)
end
@testset "allocations: remaining_dimnames_from_indexing" begin
    @test 0 == @allocated (()->remaining_dimnames_from_indexing((:a, :b, :c), (:,390,:)))()
end

@testset "remaining_dimnames_after_dropping" begin
    @test remaining_dimnames_after_dropping((:a, :b, :c), 1) == (:b, :c)
    @test remaining_dimnames_after_dropping((:a, :b, :c), 3) == (:a, :b)
    @test remaining_dimnames_after_dropping((:a, :b, :c), (1,3)) == (:b,)
    @test remaining_dimnames_after_dropping((:a, :b, :c), (3,1)) == (:b,)
end
@testset "allocations: remaining_dimnames_after_dropping" begin
    @test 0 == @allocated remaining_dimnames_after_dropping((:a,:b,:c,:d,:e), 4)
    @test 0 == @allocated (()->remaining_dimnames_after_dropping((:a,:b,:c,:d,:e), (1,3)))()
end

@testset "permute_dimnames" begin
    @test permute_dimnames((:a, :b, :c), (1, 2, 3)) == (:a ,:b, :c)
    @test permute_dimnames((:a, :b, :c), (3, 2, 1)) == (:c ,:b, :a)

    # permute_dimnames allows non-bijective "permutations"
    @test permute_dimnames((:a, :b, :c), (3, 3, 1)) == (:c ,:c, :a)
    @test permute_dimnames((:a, :b, :c), (3, 1))== (:c, :a)

    @test_throws BoundsError permute_dimnames((:a, :b, :c), (30, 30, 30))
    @test_throws BoundsError permute_dimnames((:a, :b), (1, 0))
end
@testset "allocations: permute_dimnames" begin
    if VERSION >= v"1.1"
        @test 0 == @allocated permute_dimnames((:a,:b,:c), (1,3,2))
    end
end

@testset "tuple_issubset" begin
    @test tuple_issubset((:a, :c), (:a, :b, :c)) == true
    @test tuple_issubset((:a, :b, :c), (:a, :c)) == false
end
@testset "allocations: tuple_issubset" begin
    @test 0 == @allocated tuple_issubset((:a, :c), (:a, :b, :c))
    @test 0 == @allocated tuple_issubset((:a, :b, :c), (:a, :c))
end

@testset "tuple_cat" begin
    @test tuple_cat((1, 2), (3, 4, 5), (6,)) == (1, 2, 3, 4, 5, 6)
    @test tuple_cat((1, 2)) == (1, 2)
end
@testset "allocations: tuple_cat" begin
    @test 0 == @allocated tuple_cat((1, 2), (3, 4, 5), (6,))
end


text:
module TestEAGO

import EAGO
using MathOptInterface
using Test

const MOI = MathOptInterface
const OPTIMIZER = MOI.instantiate(MOI.OptimizerWithAttributes(EAGO.Optimizer, MOI.Silent() => true))
const BRIDGED = MOI.instantiate(MOI.OptimizerWithAttributes(EAGO.Optimizer, MOI.Silent() => true), with_bridge_type = Float64)
const CONFIG = MOI.Test.Config(atol = 1e-3, rtol = 1e-3, optimal_status = MOI.OPTIMAL,
                               exclude = Any[MOI.DualObjectiveValue, MOI.ConstraintBasisStatus, MOI.VariableName, MOI.ConstraintName, MOI.delete,
                                             MOI.ConstraintDual, MOI.ListOfModelAttributesSet, MOI.add_constrained_variables])

"""
    runtests()

This function runs all functions in the this Module starting with `test_`.
"""
function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
end

"""
    test_runtests()

This function runs all the tests in MathOptInterface.Test.

Pass arguments to `exclude` to skip tests for functionality that is not
implemented or that your solver doesn't support.
"""
function test_runtests()
    MOI.Test.runtests(BRIDGED, CONFIG,
                      exclude = [# IPOPT Inherited test exclusions
                                "test_model_ScalarFunctionConstantNotZero",
                                "test_solve_TerminationStatus_DUAL_INFEASIBLE",
                                "test_linear_VectorAffineFunction_empty_row",
                                "test_solve_DualStatus_INFEASIBILITY_CERTIFICATE_",
                                "test_model_LowerBoundAlreadySet",
                                "test_model_UpperBoundAlreadySet",
                                "test_model_copy_to_UnsupportedAttribute",
                                "test_model_copy_to_UnsupportedConstraint",
                                "test_objective_set_via_modify",
                                "test_conic_linear_VectorOfVariables_2",

                                # Cbc default test exclusions
                                "test_linear_Indicator_",
                                "test_linear_SOS1_integration",
                                "test_linear_SOS2_integration",
                                "test_solve_SOS2_add_and_delete",
                                "test_conic_NormInfinityCone_INFEASIBLE",
                                "test_conic_NormOneCone_INFEASIBLE",
                                "test_solve_TerminationStatus_DUAL_INFEASIBLE",

                                # EAGO test exclusions
                                "test_attribute_NumberOfThreads",
                                "test_modification_",
                                "test_linear_integration_delete_variables",

                                # EAGO Exclusions to Resolve (by adding conic support later and fixing twice solve issues)
                                "conic_NormOneCone_VectorAffineFunction",
                                "conic_NormOneCone_VectorOfVariables",
                                "conic_NormInfinityCone_VectorOfVariables",
                                "conic_NormInfinityCone_VectorAffineFunction",
                                "test_conic_NormInfinityCone_3",
                                "test_conic_NormOneCone",
                                "test_conic_linear_VectorOfVariables",
                                "conic_linear_VectorAffineFunction",
                                "conic_linear_VectorAffineFunction_2",

                                "linear_integer_solve_twice",
                                "linear_integration",

                                "test_quadratic_SecondOrderCone_basic",
                                "test_quadratic_constraint_GreaterThan",
                                "test_quadratic_constraint_LessThan",
                                "test_quadratic_constraint_minimize",
                                "test_quadratic_duplicate_terms",
                                "test_quadratic_integration",
                                "test_quadratic_nonconvex_constraint_integration",
                                "test_quadratic_homogeneous",
                                "test_quadratic_nonhomogeneous",
                                "test_quadratic_constraint_integration",
                                "test_quadratic_constraint_basic",

                                "test_constraint_qcp_duplicate_diagonal",
                                "test_constraint_qcp_duplicate_off_diagonal",
                                "test_objective_get_ObjectiveFunction_ScalarAffineFunction",
                                "test_objective_qp_ObjectiveFunction_edge_cases",
                                "test_objective_qp_ObjectiveFunction_zero_ofdiag",

                                "test_model_ModelFilter_ListOfConstraintIndices",
                                "test_model_ModelFilter_ListOfConstraintTypesPresent"

                                ],
                      exclude_tests_after = v"0.10.5")
end

end

@testset "MOI" begin
    TestEAGO.runtests()
end

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
using Flux.Tracker

# This replicates the housing data example from the Knet.jl readme. Although we
# could have reused more of Flux (see the mnist example), the library's
# abstractions are very lightweight and don't force you into any particular
# strategy.

cd(@__DIR__)

isfile("housing.data") ||
  download("https://raw.githubusercontent.com/MikeInnes/notebooks/master/housing.data",
           "housing.data")

rawdata = readdlm("housing.data")'

# The last feature is our target -- the price of the house.

x = rawdata[1:13,:]
y = rawdata[14:14,:]

# Normalise the data
x = (x .- mean(x,2)) ./ std(x,2)

# The model

W = param(randn(1,13)/10)
b = param([0.])

# using CuArrays
# W, b, x, y = cu.((W, b, x, y))

predict(x) = W*x .+ b
meansquarederror(ŷ, y) = sum((ŷ .- y).^2)/size(y, 2)
loss(x, y) = meansquarederror(predict(x), y)

function update!(ps, η = .1)
  for w in ps
    w.data .-= w.grad .* η
    w.grad .= 0
  end
end

for i = 1:10
  back!(loss(x, y))
  update!((W, b))
  @show loss(x, y)
end

predict(x[:,1]) / y[1]

max_stars_repo_name: domluna/model-zoo

max_issues_repo_name: domluna/model-zoo

max_forks_repo_name: domluna/model-zoo


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>10.0, :alg=>TMJets(abstol=1e-15, orderT=10, orderQ=2, maxsteps=50_000)),
    "X0"   => [(1, Hyperrectangle(low=[0.0, 0.0, 0.0], high=[0.0, 0.0, 0.0]))],
    "info" => Dict(
        "name" => "cardiac_cell",
        "dim" => 3,
        "linear" => false,
        "hybrid" => true
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
#!/usr/bin/env julia

using FromFile

@from "Helper.jl" import Helper

const OXYGEN = 1
const CO2 = 2

function get_column(i::Int, lines::Vector{String})::Vector{Char}
    return [line[i] for line in lines]
end

function get_value(data::Vector{String}; which=OXYGEN)::Int
    lines = copy(data)

    col_idx = 1
    while length(lines) != 1
        column::Vector{Char} = get_column(col_idx, lines)
        zeros = count(c -> c == '0', column)
        ones = length(column) - zeros
        to_keep::Char = begin
            if which == OXYGEN
                result = '1'    # if ones >= zeros
                if zeros > ones
                    result = '0'
                end
            else    # CO2
                result = '0'    # if zeros <= ones
                if ones < zeros
                    result = '1'
                end
            end
            result
        end
        lines = [line for line in lines if line[col_idx] == to_keep]
        col_idx += 1
    end
    return parse(Int, lines[1], base=2)
end

function main()
    # fname = "example.txt"
    fname = "input.txt"

    data = Helper.read_lines(fname)

    oxygen = get_value(data, which=OXYGEN)
    println(oxygen)
    co2 = get_value(data, which=CO2)
    println(co2)
    println("---")
    println(oxygen * co2)
end

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


text:
module Exceptions

using SearchLight

export NotPersistedException, MissingDatabaseConfigurationException, DatabaseAdapterException
export UnretrievedModelException, InvalidModelException

abstract type SearchLightException <: Exception end

Base.showerror(io::IO, ex::T) where {T<:SearchLightException} = print(io, ex.msg)

struct NotPersistedException <: Exception
  model
  msg
end
NotPersistedException(model) = NotPersistedException(model, "NotPersistedException: Model is not persisted \n$model")

struct MissingDatabaseConfigurationException <: Exception
  msg::String
end
MissingDatabaseConfigurationException() = MissingDatabaseConfigurationException("The database configuration can not be found")

struct NotConnectedException <: Exception
  msg::String
end
NotConnectedException() = NotConnectedException("SearchLight is not connected to the database")

struct DatabaseAdapterException <: Exception
  msg::String
end
DatabaseAdapterException() = DatabaseAdapterException("The SearchLight database adapter has thrown an unexpected exception")

struct UnretrievedModelException <: Exception
  model
  id
  msg
end
UnretrievedModelException(model, id) = UnretrievedModelException(model, id, "UnretrievedModelException: the $(typeof(model)) data could not be retrieved for id $id. \nModel: \n$model")

struct InvalidModelException <: Exception
  model
  errors
  msg
end
InvalidModelException(model, errors) = InvalidModelException("The $(typeof(model)) model has validation errors:\n$errors")

end


text:
using TensorTrains
using Test
using Random

function ttcheckerr(Y1, Y2)
    return norm(add(1.0, Y1, -1.0, Y2)) < length(Y1) * eps(Float64)*max(norm(Y1),norm(Y2))
end

@testset "TensorTrains" begin
    @testset "dot prod" begin
        T1 = ttones(3); T2 = ttones(3)
        @test dot(T1,T2) ≈ 8.0
        T1 = ttrandn(5); T2 = ttrandn(5)
        @test dot(T1,T2) ≈ vec(decompress(T1))'*vec(decompress(T2))
        @test dot(T1,T2) ≈ vec(unfold(T1))'*vec(unfold(T2))
    # @test dot(ttones)
    end;
    @testset "norm" begin
        T1 = ttones(3);
        @test norm(T1) ≈ √8
        T1 = ttrandn(5);
        @test norm(T1)  ≈ √dot(T1,T1)
        @test norm!(T1) ≈ √dot(T1,T1)
    end;
    @testset "add" begin
        T1 = ttones(3);
        add!(T1, ttones(3));
        @test norm(T1) ≈ 2*√8;
        add!(T1, -1.0, ttones(3));
        @test norm(T1) ≈ √8;
        T2 = add(2.0, T1, -1.0, T1);
        @test norm(T2) ≈ norm(T1);
    end
    @testset "orthogonality" begin
        T1 = ttones(4);
        add!(T1, ttones(4));
        add!(T1, ttones(4));
        rightOrth!(T1);
        @test norm(T1[end]) ≈ norm(T1);
        leftOrth!(T1);
        @test norm(T1[1]) ≈ norm(T1);
    end
    @testset "svd" begin
        T1 = ttones(4);
        add!(T1, ttones(4));
        add!(T1, ttones(4));
        svdtrunc!(T1, svd!(T1), 1e-16);
        @test ranks(T1) == [1,1,1];
        Random.seed!(1234);
        T1 = ttrandn(8);
        add!(T1, ttrandn(8));
        add!(T1, ttrandn(8));
        T2 = deepcopy(T1);
        svdtrunc!(T1, svd!(T1), 1e-1);
        @test norm(add(1.0, T1, -1.0, T2)) < 2e-1;
        T1 = deepcopy(T2);
        svdtrunc!(T1, svd!(T1), 2e-1);
        @test norm(add(1.0, T1, -1.0, T2)) < 4e-1;
    end
    @testset "matvec" begin
        Random.seed!(1234);
        d = 8;
        Id = tteye(d);
        T1 = ttrandn(d);
        @test norm(T1) ≈ norm(matvec(Id,T1))
        Y1 = matvec(add(1.0,Id,1.0,Id),T1);
        Y2 = matvec(scale(Id,2.0),T1);
        @test ttcheckerr(Y1,Y2);
        A1 = deepcopy(Id);
        add!(A1, 2.5, Id);
        Y1 = matvec(A1, T1);
        A2 = deepcopy(Id);
        scale!(A2, 3.5);
        Y2 = matvec(A2, T1);
        @test ttcheckerr(Y1,Y2);
    end
    @testset "matmat" begin
        Random.seed!(1234);
        d = 8;
        Id = tteye(d);
        A1 = deepcopy(Id);
        scale!(A1,Float64(π));
        T1 = ttrandn(d);
        A2 = matmat(Id,A1);
        Y1 = matvec(A1,T1);
        Y2 = matvec(A2,T1);
        @test ttcheckerr(Y1,Y2);
        A2 = matmat(A2,Id);
        Y2 = matvec(A2,T1);
        @test ttcheckerr(Y1,Y2);
    end
    @testset "hadamard" begin
        Random.seed!(1234);
        T1 = ttones(6);
        T2 = ttrandn(6);
        @test norm(hadamard(T1,T2)) ≈ norm(T2)
    end
    @testset "kron" begin
        Random.seed!(1234);
        A1 = tteye(8);
        A2 = tteye(8);
        T1 = ttrandn(8);
        T2 = ttrandn(8);
        scale!(A1, 0.5);
        scale!(A2, Float64(π));
        Y1 = matvec(ttkron(A1,A2),ttkron(T1,T2));
        Y2 = ttkron(matvec(A1,T1),matvec(A2,T2));
        @test ttcheckerr(Y1,Y2);
    end
end;


text:
# Estimates the Klein consumption equation by GMM
using Econometrics, DelimitedFiles, Statistics
function KleinMoments(θ, y, x, z)
	e = y - x*θ
	m = e.*z
end

function main()
data = readdlm("klein.data")
# construct missing lags, and drop first row that has missing data
profits = data[:,3]
output = data[:,7]
data = [data lag(profits,1) lag(output,1)]
data = data[2:end,:]
n = size(data,1)
# define instruments
exogs = [1, 6, 8, 9, 10, 11, 12]
exogs = data[:,exogs]
exogs = [ones(n,1) exogs]
# CONSUMPTION
println("CONSUMPTION EQUATION")
# define variables in consumption equation
y = data[:,2]
profits = data[:,3]
lagprofits = data[:,11]
wp = data[:,4]
wg = data[:,8]
wages = wp + wg
# regressors in consumption equation
x = [profits lagprofits wages]
x = [ones(n,1) x]
# GMM estimation using CUE
θstart = x\y # ols start values
names = ["Constant", "Profits", "Profits-1", "Wages"]
moments = θ -> KleinMoments(θ, y, x, exogs)
# estimation results using CUE
gmmtitle = "Klein model 1 GMM example CUE"
gmmresults(moments, θstart, "",  gmmtitle, names)
nothing
end
main()

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
module DOM

using Hyperscript

const global_unique_dom_id_counter = Ref(0)

"""
    get_unique_dom_id()

We could use a unique ID like uuid4, but since every dom element gets
such an id, I prefer to keep the id as short as possible, so we just use a counter.
"""
function get_unique_dom_id()
    global_unique_dom_id_counter[] += 1
    return string(global_unique_dom_id_counter[])
end

"""
Dome node with unique ID, to make it easier to interpolate it.
"""
function um(tag, args...; kw...)
    m(tag, args..., dataJscallId = get_unique_dom_id(); kw...)
end

function m_unesc(tag, args...; kw...)
    m(Hyperscript.NOESCAPE_HTMLSVG_CONTEXT, tag, args..., dataJscallId = get_unique_dom_id(); kw...)
end

for node in [:a, :abbr, :address, :area, :article, :aside, :audio, :b,
             :base, :bdi, :bdo, :blockquote, :body, :br, :button, :canvas,
             :code, :col, :colgroup, :data, :datalist, :caption, :cite,
             :dd, :del, :details, :dfn, :dialog, :div, :dl, :dt, :em, :embed,
             :figcaption, :figure, :footer, :form, :h1, :h2, :fieldset,
             :h3, :h4, :h5, :h6, :head, :header, :hgroup, :hr, :html, :i, :iframe,
             :input, :ins, :kbd, :label, :legend, :li, :link, :img,
             :main, :map, :mark, :math, :menu, :menuitem, :meta, :meter, :nav,
             :object, :ol, :optgroup, :option, :output, :p, :param, :noscript,
             :picture, :pre, :progress, :q, :rb, :rp, :rt, :rtc, :ruby, :s, :samp,
             :section, :select, :slot, :small, :source, :span,
             :strong, :sub, :summary, :sup, :svg, :table, :tbody, :td, :template,
             :textarea, :tfoot, :th, :thead, :time, :title, :tr,
             :track, :u, :ul, :var, :video, :wbr, :font]

    node_name = string(node)
    unesc = Symbol(node_name * "_unesc")
    @eval $(node)(args...; kw...) = um($(node_name), args...; kw...)
    @eval $(unesc)(args...; kw...) = m_unesc($(node_name), args...; kw...)
end

style(args...; kw...) = m_unesc("style", args...; kw...)
script(args...; kw...) = m_unesc("script", args...; kw...)

end

using .DOM

function selector(node)
    query_string = "[data-jscall-id=$(repr(uuid(node)))]"
    return js"(document.querySelector($(query_string)))"
end

# default turn attributes into strings
attribute_render(session::Session, parent, attribute::String, x) = string(x)
attribute_render(session::Session, parent, attribute::String, x::Nothing) = x
attribute_render(session::Session, parent, attribute::String, x::Bool) = x

function attribute_render(session::Session, parent, attribute::String, obs::Observable)
    onjs(session, obs, js"value=> JSServe.update_node_attribute($(parent), $attribute, value)")
    return attribute_render(session, parent, attribute, obs[])
end

struct DontEscape
    x::Any
end

function Hyperscript.printescaped(io::IO, x::DontEscape, escapes)
    print(io, x.x)
end

function attribute_render(session::Session, parent, attribute::String, jss::JSCode)
    register_resource!(session, jss)
    return string(jss)
end

function attribute_render(session::Session, parent, attribute::String, jss::Asset)
    if parent isa Hyperscript.Node{Hyperscript.CSS}
        # css seems to require an url object
        return "url($(url(jss, session.url_serializer)))"
    else
        return "$(url(jss, session.url_serializer))"
    end
end

render_node(session::Session, x) = x

const BOOLEAN_ATTRIUTES = Set([
    "allowfullscreen",
    "allowpaymentrequest",
    "async",
    "autofocus",
    "autoplay",
    "checked",
    "controls",
    "default",
    "defer",
    "disabled",
    "formnovalidate",
    "hidden",
    "ismap",
    "itemscope",
    "loop",
    "multiple",
    "muted",
    "nomodule",
    "novalidate",
    "open",
    "readonly",
    "required",
    "reversed",
    "selected",
    "typemustmatch"
])

is_boolean_attribute(attribute::String) = attribute in BOOLEAN_ATTRIUTES

function render_node(session::Session, node::Node)
    # give each node a unique id inside the dom
    new_attributes = Dict{String, Any}()
    newchildren = map(children(node)) do elem
        return jsrender(session, elem)
    end
    for (k, v) in Hyperscript.attrs(node)
        rendered = attribute_render(session, node, k, v)
        # We code nothing to mean omitting the attribute!
        if is_boolean_attribute(k)
            if rendered isa Bool
                if rendered
                    # only add attribute if true!
                    new_attributes[k] = true
                end
            else
                error("Boolean attribute $(k) expects a boolean! Found: $(typeof(rendered))")
            end
        else
            new_attributes[k] = rendered
        end
    end
    return Node(
        Hyperscript.context(node),
        Hyperscript.tag(node),
        newchildren,
        new_attributes
    )
end

# jsrender(session, x) will be called anywhere...
# if there is nothing sessions specific in the dom, fallback to jsrender without session
function jsrender(session::Session, node::Node)
    render_node(session, node)
end

function uuid(node::Node)
    get(Hyperscript.attrs(node), "data-jscall-id") do
        error("Node $(node) doesn't have a unique id. Make sure to use DOM.$(Hyperscript.tag(node))")
    end
end

"""
    jsrender([::Session], x::Any)
Internal render method to create a valid dom. Registers used observables with a session
And makes sure the dom only contains valid elements. Overload jsrender(::YourType)
To enable putting YourType into a dom element/div.
You can also overload it to take a session as first argument, to register
messages with the current web session (e.g. via onjs).
"""
jsrender(::Session, x::Any) = jsrender(x)
jsrender(::Session, x::Symbol) = DOM.p(string(x))
jsrender(::Session, x::Hyperscript.Styled) = x
jsrender(x) = x


text:
"""
    PyramidMatchGraphKernel <: AbstractGraphKernel

A graph kernel that compares vertex embeddings of graph using pyramid match.

For each graph, the vertices are mapped into a hypercube of dimension `embedding_dim`.
The hypercube is the split into buckets of certain granularity, controlled by `histogram_levels`,
and the number of vertices in a bucket for both graphs are compared.
If vertex labels are provided, only vertices with the same label in a bucket are compared.

# Keywords
- `embedding_dim=6`: The number of dimensions of the space where the vertices are embedded.
    This is an upper bound - in some cases only an embedding of lower dimension is created.
- `histogram_levels=4`: The number of levels for which histograms are compared. For
    each level `l ∈ 0:histogram_levels`, the embeddings are split into 2^l buckets and
    the number of vertices in that bucket are compared separately.
- `vertex_labels`: By which vertex labels the vertices are grouped together. The embedding
    is created on all vertices, but then the histograms are separately compared for each
    group of vertex labels. Finally the results of each group are added together.
    Either `Auto`, `:` (all), or a tuple of vertex value keys.

# References
<https://ojs.aaai.org/index.php/AAAI/article/view/10839>
"""
struct PyramidMatchGraphKernel{VL <: LabelsType} <: AbstractGraphKernel

    embedding_dim::Int
    histogram_levels::Int
    vertex_labels::VL

    function PyramidMatchGraphKernel(;embedding_dim::Integer=6, histogram_levels::Integer=4, vertex_labels=Auto())

        embedding_dim >= 1 || throw(DomainError(embedding_dim, "embedding_dim must be >= 1"))
        histogram_levels >= 0 || throw(DomainError(embedding_dim, "histogram_levels must be >= 0"))

        return new{typeof(vertex_labels)}(Int(embedding_dim), Int(histogram_levels), vertex_labels)
    end
end

function preprocessed_form(kernel::PyramidMatchGraphKernel, g::AbstractGraph)

    embedding = _embedding(g, kernel.embedding_dim)
    d = size(embedding, 2) # d might sometimes be smaller than embedding_dim

    vertex_labels = _labels(kernel.vertex_labels, vertexvals_type(g))

    # TODO, we should deduce the correct type instead of using Any
    vertex_class_sizes = counter(Any)
    for v in vertices(g)
        class_key = tuple((get_vertexval(g, v, i) for i in vertex_labels)...)
        inc!(vertex_class_sizes, class_key)
    end

    # The rows of embedding correspond to vertices. We group vertices together
    # by their labels and associate a sub matrix of the corresponding rows
    # with each label class
    insert_row_index = Dict{Any, Int}()
    points_by_vertex_class = Dict{Any, Matrix{Float64}}()
    for (k, n_rows) in vertex_class_sizes
        insert_row_index[k] = 0
        points_by_vertex_class[k] = Matrix{Float64}(undef, n_rows, d)
    end

    for (v, row) in enumerate(eachrow(embedding))
        k = tuple((get_vertexval(g, v, i) for i in vertex_labels)...)
        row_index = (insert_row_index[k] += 1)
        points_by_vertex_class[k][row_index, :] = row
    end

    for points in values(points_by_vertex_class)
        # the j-th column of this matrix contains the j-th coordinates
        # of the embedding vectors sorted in ascending order
        sort!(points; dims=1)
    end

    return points_by_vertex_class
end

function apply_preprocessed(kernel::PyramidMatchGraphKernel, points_by_class_1, points_by_class_2)

    L = kernel.histogram_levels
    result = 0.0

    for (class, points1) in points_by_class_1

        points2 = get(points_by_class_2, class, nothing)
        points2 == nothing && continue

        # TODO the embedding dim might actually be smaller for some graph
        # we should maybe consider some kind of scaling in such a case
        d = min(size(points1, 2), size(points2, 2))
        n1 = size(points1, 1)
        n2 = size(points2, 1)

        hist_intersect = zeros(Int64, L + 1)

        # TODO ensure no int overflow, same result on 32 bit platform

        # TODO we don't need to store hist_intersect as a vector
        # TODO add @inbounds
        hist_intersect[1] = d * min(n1, n2)
        for l in 1:L
            cell_boundaries = range(0.0, 1.0, length=2^l + 1)
            for j in 1:d
                i1 = 1
                i2 = 1
                while i1 <= n1
                    # TODO is is possible that searchsortedlast is not implemented
                    # efficiently on a StepRangeLen
                    cell_num = searchsortedlast(cell_boundaries, points1[i1, j])
                    # TODO maybe we should verify here, that cell_num is a valid index
                    cell_lower = cell_boundaries[cell_num]
                    cell_upper = cell_boundaries[cell_num + 1]
                    # TODO maybe we need some correction for rounding errors here

                    i1 += 1
                    num1 = 1
                    # the first loop is just for safety we probably don't need it
                    while i1 <= n1 && points1[i1, j] < cell_lower
                        i1 += 1
                    end
                    # count number of vertices from g1 that fall into the bucket
                    # specified [cell_lower, cell_upper) along the j-th dimension
                    # of the hypercube
                    while i1 <= n1 && points1[i1, j] < cell_upper
                        num1 += 1
                        i1 += 1
                    end

                    # count number of vertices from g2 that fall into that bucket
                    # We could also use binary search here
                    while i2 <= n2 && points2[i2, j] < cell_lower
                        i2 += 1
                    end
                    num2 = 0
                    while i2 <= n2 && points2[i2, j] < cell_upper
                        num2 += 1
                        i2 += 1
                    end
                    hist_intersect[l + 1] += min(num1, num2)
                end
            end
        end

        result += (hist_intersect[L + 1] + sum(l -> (hist_intersect[l+1] - hist_intersect[l+2]) / 2^(L - l), 0:L-1))

    end
    return result
end

function _embedding(g::AbstractGraph, embedding_dim::Int)

    evs = eigvecs(adjacency_matrix(g); sortby=x -> -abs(x))
    embedding = abs.(@view evs[:, 1:min(size(evs, 2), embedding_dim)])
    # theoretically clamping is not necessary, this is just a precaution
    # against rounding errors. Clamping to prevfloat(-1.0) should make some calculations easier
    clamp!(embedding, 0.0, prevfloat(1.0))
    return embedding
end


text:
using RegularizedRegression: LassoRegression, RidgeRegression, learn, coefficients

@testset "Lasso regression (via GLMNet.jl)" begin
    for _ in 1:NumTestRepetitions
        test_learner(LassoRegression, 5.0) # Seems this needs to be higher for SLOPE than for SCAD and MCP
    end
end

@testset "Ridge regression (via GLMNet.jl)" begin
    for _ in 1:NumTestRepetitions
        test_learner(RidgeRegression, 5.0) # Seems this needs to be higher for SLOPE than for SCAD and MCP
    end
end


text:
__precompile__(true)
module WeakRefStrings

export WeakRefString, WeakRefStringArray

using Nulls

"""
A custom "weakref" string type that only points to external string data.
Allows for the creation of a "string" instance without copying data,
which allows for more efficient string parsing/movement in certain data processing tasks.

**Please note that no original reference is kept to the parent string/memory, so `WeakRefString` becomes unsafe
once the parent object goes out of scope (i.e. loses a reference to it)**

Internally, a `WeakRefString{T}` holds:

  * `ptr::Ptr{T}`: a pointer to the string data (code unit size is parameterized on `T`)
  * `len::Int`: the number of code units in the string data
  * `ind::Int`: a field that can be used to store an integer, like an index into an array; this can be helpful
                in certain cases when the underlying source may need to move around (which would invalidate
                the WeakRefString's `ptr` field), a new WeakRefString can created using the same offset into
                the parent data as the old one.
"""
struct WeakRefString{T} <: AbstractString
    ptr::Ptr{T}
    len::Int # of code units
    ind::Int # used to keep track of a string data index
end

WeakRefString(ptr::Ptr{T}, len) where {T} = WeakRefString(ptr, Int(len), 0)
WeakRefString(t::Tuple{Ptr{T}, Int, Int}) where {T} = WeakRefString(t[1], t[2], t[3])

const NULLSTRING = WeakRefString(Ptr{UInt8}(0), 0)
const NULLSTRING16 = WeakRefString(Ptr{UInt16}(0), 0)
const NULLSTRING32 = WeakRefString(Ptr{UInt32}(0), 0)
Base.zero(::Type{WeakRefString{T}}) where {T} = WeakRefString(Ptr{T}(0), 0, 0)
Base.endof(x::WeakRefString) = x.len
Base.length(x::WeakRefString) = x.len
Base.next(x::WeakRefString, i::Int) = (Char(unsafe_load(x.ptr, i)), i + 1)

import Base: ==
function ==(x::WeakRefString{T}, y::WeakRefString{T}) where {T}
    x.len == y.len && (x.ptr == y.ptr || ccall(:memcmp, Cint, (Ptr{T}, Ptr{T}, Csize_t),
                                           x.ptr, y.ptr, x.len) == 0)
end

function Base.hash(s::WeakRefString{T}, h::UInt) where {T}
    h += Base.memhash_seed
    ccall(Base.memhash, UInt, (Ptr{T}, Csize_t, UInt32), s.ptr, s.len, h % UInt32) + h
end

Base.show(io::IO, ::Type{WeakRefString{T}}) where {T} = print(io, "WeakRefString{$T}")
function Base.show(io::IO, x::WeakRefString{T}) where {T}
    print(io, '"')
    for c in x
        print(io, c)
    end
    print(io, '"')
    return
end

chompnull(x::WeakRefString{T}) where {T} = unsafe_load(x.ptr, x.len) == T(0) ? x.len - 1 : x.len

Base.string(x::WeakRefString{UInt16}) = x == NULLSTRING16 ? "" : String(transcode(UInt8, unsafe_wrap(Array, x.ptr, chompnull(x))))
Base.string(x::WeakRefString{UInt32}) = x == NULLSTRING32 ? "" : String(transcode(UInt8, unsafe_wrap(Array, x.ptr, chompnull(x))))

Base.convert(::Type{WeakRefString{UInt8}}, x::String) = WeakRefString(pointer(x), length(x))
Base.convert(::Type{String}, x::WeakRefString) = convert(String, string(x))
Base.string(x::WeakRefString) = x == NULLSTRING ? "" : unsafe_string(x.ptr, x.len)
Base.String(x::WeakRefString) = string(x)

struct WeakRefStringArray{T <: ?WeakRefString, N} <: AbstractArray{T, N}
    data::Vector{Vector{UInt8}}
    elements::Array{T, N}
end

WeakRefStringArray(data::Vector{UInt8}, ::Type{T}, rows::Integer) where {T <: ?WeakRefString} = WeakRefStringArray([data], Vector{T}(zeros(Nulls.T(T), rows)))
WeakRefStringArray(data::Vector{UInt8}, A::Array{T}) where {T <: ?WeakRefString} = WeakRefStringArray([data], A)

Base.size(A::WeakRefStringArray) = size(A.elements)
Base.getindex(A::WeakRefStringArray, i::Int) = A.elements[i]
Base.getindex(A::WeakRefStringArray{T, N}, I::Vararg{Int, N}) where {T, N} = A.elements[I...]
Base.setindex!(A::WeakRefStringArray{T, N}, v, i::Int) where {T, N} = setindex!(A.elements, v, i)
Base.setindex!(A::WeakRefStringArray{T, N}, v, I::Vararg{Int, N}) where {T, N} = setindex!(A.elements, v, I...)
Base.resize!(A::WeakRefStringArray, i) = resize!(A.elements, i)

function Base.append!(a::WeakRefStringArray{T, 1}, b::WeakRefStringArray{T, 1}) where {T}
    append!(a.data, b.data)
    append!(a.elements, b.elements)
    return a
end

end # module


text:
using KROME
using Printf

function av_slab(io::IO=stdout)
  nmols = krome_nmols()[] # read Fortran module variable
  n = zeros(nmols) # default abundances (number density)

  idx_H = krome_idx_H()[] + 1
  idx_H2 = krome_idx_H2()[] + 1
  idx_Hj = krome_idx_Hj()[] + 1
  idx_Cj = krome_idx_Cj()[] + 1
  idx_CO = krome_idx_CO()[] + 1
  idx_O = krome_idx_O()[] + 1
  idx_E = krome_idx_E()[] + 1

  # == Variable for NL97 network
  dust_to_gas_ratio = 1.
  deff = 1.

  tmax = 3.1536e16
  dt = 3.1536e12

  rho = 3.841e-22

  abundc = 1.e-4
  abundo = 3.e-4
  abundhe = 0.1
  abundS = 0.#0.6e-6
  abundSi = 3.37e-6

  AV_conversion_factor =  6.289e-22 #5.348e-22

  krome_init() # init krome (mandatory)

  abar = 1 + abundhe*4. + abundc*12. + abundo*16. + abundSi*28. + abundS*32.
  # == CHANGE HERE, F1: n=1e3 with G0=1e1 or F4: n=10^5.5 with G0=1e5 ==
  nH = 10.0^3 #10.^5.5#(=rho/abar/krome_p_mass()[])
  Ghab = 1.e1*1.69/2. # 1.e5*1.69/2.
  crate = 5.e-17

  opt_depth_max = 100.
  opt_depth = 1e-6

  NH2 = 0.
  H2sh = 1.

  println(io, "opt_depth_1,n(idx_H2)_2,n(idx_H)_3,n(idx_Hj)_4,n(idx_CO)_5,n(idx_Cj)_6,n(idx_E)_7,n(idx_O)_8,T_9,H2sh_10")

  while opt_depth <= opt_depth_max
    n .= 1.e-40
    n[idx_H] = nH
    n[idx_Cj] = abundc*nH
    n[idx_O] = abundo*nH
    n[idx_E] = n[idx_Cj]
    T = 50.
    Tdust = 10.

    # == Calculations for reactions 3 and 7
    fshield_dust    = exp(-2.5e0 * opt_depth)
    G_dust  = Ghab * fshield_dust

    ch34 = 5.087e2 * T^1.586e-2
    ch35 = -0.4723e0 - 1.102e-5 * log(T)

    if n[idx_E] == 0e0
      # If the fractional ionization is zero, then there won't be any recombination,
      # so the value we use for phi doesn't matter too much -- 1d20 is simply an
      # arbitrary large number

      phi = 1e20
    else
      phi = G_dust * sqrt(T) / (n[idx_E])
    end

    if phi < 1e-6
      h_gr = 1.225e-13 * dust_to_gas_ratio
    else
      hgrvar1  = 8.074e-6 * phi^1.378e0
      hgrvar2  = (1e0 + ch34 * phi^ch35)
      h_gr     = 1.225e-13 * dust_to_gas_ratio / (1e0 + hgrvar1 * hgrvar2)
    end

    uv_ion = 0.
    # == End: Calculations for reactions 3 and 7

    COsh = 0.5

    krome_set_user_ghab(Ghab)
    krome_set_user_crate(crate)
    krome_set_user_av(opt_depth)
    krome_set_user_h2self(H2sh)
    krome_set_user_tdust(Tdust)
    krome_set_user_coself(COsh)
    krome_set_user_uv_ion(uv_ion)
    krome_set_user_h_gr(h_gr)
    krome_set_user_dust_to_gas_ratio(dust_to_gas_ratio)
    krome_set_user_deff(deff)

    i=0
    ttot=0.

    while ttot < tmax

      dtC=dt
      ttot=ttot+dtC

      if  ttot > tmax
        dtC = dtC - (ttot-tmax)
        ttot = tmax
      end

      # == Calculations for reactions 5, this needs to be done here in this loop as it depends on the abundance which change over time

      NHtot = opt_depth / (dust_to_gas_ratio * AV_conversion_factor)

      p1 = log10(NHtot / 1e18)
      if p1 < 0e0
        p1 = 0e0
      elseif p1 > 5e0
        p1 = 5e0
      end

      if n[idx_E]/nH < 1e-4
        p2 = -4.
      elseif n[idx_E]/nH > 0.1
        p2 = -1.
      else
        p2 = log10(n[idx_E]/nH)
      end

      f4 = 1.06  + 4.08e-2 * p2 + 6.51e-3 * p2^2e0
      f5 = 1.90  + 0.678   * p2 + 0.113   * p2^2e0
      f6 = 0.990 - 2.74e-3 * p2 + 1.13e-3 * p2^2e0

      ion = f4 * (-15.6 - 1.10 * p1 + 9.13e-2 * p1^2e0) + f5 * 0.87 * exp(-((p1 - 0.41) / 0.84)^2e0)

      heat = f6 * (-26.5 - 0.920 * p1 + 5.89e-2 * p1^2e0) + f6 * 0.96 * exp(-((p1 - 0.38) / 0.87)^2e0)

      ion = 1e1^ion
      heat = 1e1^heat

      krome_set_user_xr_ion(ion)

      # == End: Calculations for reactions 5

      # Convert to modifiable variables and then back
      T_ = fill(T)
      dtC_ = fill(dtC)
      # Original call: krome(n,T,dtC)
      krome(n,T_,dtC_)
      T = T_[]
      dtC = dtC_[]

      i=i+1

    end

    # write(*,*) opt_depth,n(idx_H2),n(idx_H),n(idx_Hj),n(idx_CO),n(idx_Cj),n(idx_E),n(idx_O),T,H2sh
    println(io,
            join(map(fsn, [opt_depth,
                            n[idx_H2],n[idx_H],n[idx_Hj],n[idx_CO],n[idx_Cj],n[idx_E],n[idx_O],
                            T,H2sh]), " "))

    NH_old = opt_depth/AV_conversion_factor

    opt_depth = opt_depth*1.1

    NH2 = NH2 + (opt_depth/AV_conversion_factor-NH_old)*n[idx_H2]/nH
    b5 = sqrt(1.38065e-16*T/krome_p_mass()[])/1e5
    x = NH2/5e14
    H2sh = 0.965/(1+x/b5)^2. + 0.035/sqrt(1+x)*exp(-8.5e-4*sqrt(1+x))
  end
end

# Return formatted value in Fortran scientific notation
function fsn(value, width=25, precision=16)
  @assert precision <= 64 "precision must be <= 64"

  # Extract exponent and base
  if value == 0.0
    exponent = -1
    base = 0.0
  else
    exponent = floor(Int, log10(abs(value)))
    base = value / float(10)^exponent
  end

  # Make base have a leading zero
  base /= 10
  exponent += 1

  # Get formatted values
  base_formatted = (@sprintf "%.64f" base)[1:(precision + 2 + (base < 0))]
  exponent_formatted = @sprintf "E%+04d" exponent

  padding = width - length(base_formatted) - length(exponent_formatted)

  return " "^padding * base_formatted * exponent_formatted
end


text:
using AdventOfCodeSolutions
using Test

@testset "AdventOfCodeSolutions.jl" begin
    foreach(include, readdir("../2020/src/", join=true))
end