text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "IpoptBuilder"
version = v"3.12.10"

# Collection of sources required to build IpoptBuilder
sources = [
    "https://github.com/coin-or/Ipopt/archive/releases/3.12.10.tar.gz" =>
    "dfd29dc95ec815e1ff0a3b7dc86ecc8944b24977e40724c35dac25aa192ac3cd",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd Ipopt-releases-3.12.10/
update_configure_scripts
mkdir build
cd build/

for path in ${LD_LIBRARY_PATH//:/ }; do
    for file in $(ls $path/*.la); do
        echo "$file"
        baddir=$(sed -n "s|libdir=||p" $file)
        sed -i~ -e "s|$baddir|'$path'|g" $file
    done
done

if [ $target = "x86_64-apple-darwin14" ]; then
  # seems static linking requires apple's ar
  export AR=/opt/x86_64-apple-darwin14/bin/x86_64-apple-darwin14-ar

   # Ignore the "# Don't fix this by using the ld -exported_symbols_list flag, it doesn't exist in older darwin lds"
  # seems to work for the current version and otherwise a long list of non-Clp symbols are exported
  sed -i~ -e "s|~nmedit -s \$output_objdir/\${libname}-symbols.expsym \${lib}| -exported_symbols_list \$output_objdir/\${libname}-symbols.expsym|g" ../configure
fi

export CPPFLAGS="-DCOIN_USE_MUMPS_MPI_H"

## STATIC BUILD START
# Staticly link all dependencies and export only Ipopt symbols

# force only exporting symbols related to Ipopt
# SetIntermediateCallback is to fix https://github.com/JuliaOpt/IpoptBuilder/issues/2
sed -i~ -e 's|LT_LDFLAGS="-no-undefined"|LT_LDFLAGS="-no-undefined -export-symbols-regex \\\\"Ipopt\|SetIntermediateCallback\\\\""|g' ../configure
sed -i~ -e 's|LT_LDFLAGS="-no-undefined"|LT_LDFLAGS="-no-undefined -export-symbols-regex \\\\"Ipopt\|SetIntermediateCallback\\\\""|g' ../Ipopt/configure

../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --enable-shared --disable-static \
--enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
--with-asl-lib="-L${prefix}/lib -lasl" --with-asl-incdir="$prefix/include/asl" \
--with-blas="-L${prefix}/lib -lcoinblas -lgfortran" \
--with-lapack="-L${prefix}/lib -lcoinlapack" \
--with-metis-lib="-L${prefix}/lib -lcoinmetis" --with-metis-incdir="$prefix/include/coin/ThirdParty" \
--with-mumps-lib="-L${prefix}/lib -lcoinmumps -lcoinmetis" --with-mumps-incdir="$prefix/include/coin/ThirdParty"

## STATIC BUILD END

## DYNAMIC BUILD START
#../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --enable-shared --enable-static \
#--enable-dependency-linking lt_cv_deplibs_check_method=pass_all \
#--with-asl-lib="-L${prefix}/lib -lasl" --with-asl-incdir="$prefix/include/asl" \
#--with-blas="-L${prefix}/lib -lcoinblas -lgfortran" \
#--with-lapack="-L${prefix}/lib -lcoinlapack" \
#--with-metis-lib="-L${prefix}/lib -lcoinmetis" --with-metis-incdir="$prefix/include/coin/ThirdParty" \
#--with-mumps-lib="-L${prefix}/lib -lcoinmumps" --with-mumps-incdir="$prefix/include/coin/ThirdParty"
## DYNAMIC BUILD END

make -j${nproc}
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]
platforms = expand_gcc_versions(platforms)
# To fix gcc4 bug in Windows
#platforms = setdiff(platforms, [Windows(:x86_64, compiler_abi=CompilerABI(:gcc4)), Windows(:i686, compiler_abi=CompilerABI(:gcc4))])
push!(platforms, Windows(:i686,compiler_abi=CompilerABI(:gcc6)))
push!(platforms, Windows(:x86_64,compiler_abi=CompilerABI(:gcc6)))

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libipopt", :libipopt),
     ExecutableProduct(prefix, "ipopt", :amplexe)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/juan-pablo-vielma/ASLBuilder/releases/download/v3.1.0-1-static/build_ASLBuilder.v3.1.0.jl",
    "https://github.com/juan-pablo-vielma/COINBLASBuilder/releases/download/v1.4.6-1-static/build_COINBLASBuilder.v1.4.6.jl",
    "https://github.com/juan-pablo-vielma/COINLapackBuilder/releases/download/v1.5.6-1-static/build_COINLapackBuilder.v1.5.6.jl",
    "https://github.com/juan-pablo-vielma/COINMetisBuilder/releases/download/v1.3.5-1-static/build_COINMetisBuilder.v1.3.5.jl",
    "https://github.com/juan-pablo-vielma/COINMumpsBuilder/releases/download/v1.6.0-1-static/build_COINMumpsBuilder.v1.6.0.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: juan-pablo-vielma/IpoptBuilder

max_issues_repo_name: juan-pablo-vielma/IpoptBuilder

max_forks_repo_name: juan-pablo-vielma/IpoptBuilder


text:
using MolecularIntegrals, Profile

function profints()
    x=y=z=0.0
    xyz = [x,y,z]
    xyza = xyz + [0.1,0.05,0.025]
    xa,ya,za = xyza
    ex = 1
    A = B = xyz
    ax,ay,az = bx,by,bz = xyz
    C = D = xyza
    cx,cy,cz = dx,dy,dz = xyza
    aexpn=bexpn=cexpn=dexpn = ex
    ashell,bshell,cshell,dshell = (2,2,2,2)
    #@profview MolecularIntegrals.hrr2(ashell,bshell,cshell,dshell, aexpn,bexpn,cexpn,dexpn, A,B,C,D);
    #@profview MolecularIntegrals.hrr3(ashell,bshell,cshell,dshell, aexpn,bexpn,cexpn,dexpn, A,B,C,D);
    #@profview MolecularIntegrals.hrr1(ashell,bshell,cshell,dshell, aexpn,bexpn,cexpn,dexpn, A,B,C,D);
    @profview MolecularIntegrals.hrr(ashell,bshell,cshell,dshell, aexpn,bexpn,cexpn,dexpn, A,B,C,D);
end

function profall()
    bfs = build_basis(ethane,"cc-pvdz") #"6-31G")
    #@profview MolecularIntegrals.all_twoe_ints_chrr(bfs)
    @profview MolecularIntegrals.all_twoe_ints_rys(bfs)
end

profall()

ext: jl


text:
using Test
import BertWordPieceTokenizer as BWP

# initialize tokenizer
BWP.init("https://huggingface.co/bert-base-uncased/raw/main/vocab.txt", do_lowercase=true)

@test BWP.tokenize("I like apples") == ["[CLS]", "i", "like", "apples", "[SEP]"]
@test BWP.encode("I like apples") == ([101, 1045, 2066, 18108, 102], [0, 0, 0, 0, 0])
@test BWP.decode([101, 1045, 2123, 2102, 2066, 18108, 102]) == "i dont like apples"
@test BWP.tokenize("I like apples", second_text="I dont like apples") == ["[CLS]", "i", "like", "apples", "[SEP]", "i", "don", "##t", "like", "apples", "[SEP]"]
@test BWP.encode("I like apples", second_text="I dont like apples") == ([101, 1045, 2066, 18108, 102, 1045, 2123, 2102, 2066, 18108, 102], [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1])

BWP.enable_truncation(5)
@test BWP.tokenize("I like apples", second_text="I dont like apples") == ["[CLS]", "i", "[SEP]", "i", "[SEP]"]

BWP.enable_truncation(5, truncation="pre")
@test BWP.tokenize("I like apples", second_text="I dont like apples") == ["[CLS]", "apples", "[SEP]", "apples", "[SEP]"]


text:
function em_checkmodel(pmodel::ParametrizedSSM, params::SSMParameters)

    allzero(x::Matrix) = all(x .== 0)

    m = pmodel(params)
    C = m.C(1)

    I_x, I_u  = eye(m.nx), eye(m.nu)
    I_q0      = all(pmodel.G(1) .== 0, 2) |> vec |> float |> diagm
    I_r0      = all(pmodel.H(1) .== 0, 2) |> vec |> float |> diagm
    M         = m.A(1) .!= 0 |> float
    #I_is_q    = M^m.nx  #TODO - work out indirect stochastic row selector

    @assert kron(I_x, I_r0 * C * I_q0) * pmodel.A2.D  |> allzero
    @assert kron(I_u, I_r0 * C * I_q0) * pmodel.B2.D  |> allzero
    @assert kron(I_x, I_r0) * pmodel.C2.D             |> allzero
    @assert kron(I_u, I_r0) * pmodel.D2.D             |> allzero
    @assert kron(I_x, I_q0) * pmodel.A2.D             |> allzero
    #@assert kron(I_u, I_is_q) * pmodel.B2.D           |> allzero

end #em_modelcheck


text:
using MNIST
using Test

@testset "MNIST.jl" begin
    # Write your tests here.
end


text:
using Distributions, LinearAlgebra, Random, ProgressMeter

function test_expectation_muASAmu(D=3, seed=1, numsamples=10_000)

    rg = MersenneTwister(seed)

    μα = randn(rg, D)

    Σα = rand(rg, D, D)*0.2; Σα = Σα*Σα' .+ 0.1

    Apdf = MvNormal(μα, Σα)

    S = Diagonal(rand(rg, D) .+ 2.0) +1.0I

    μf = randn(rg, D)

    f(A) = μf'*A'*inv(S)*A*μf

    g(A) = f(Diagonal(A))[1]

    #————————————————————————————————————————————————
    # This is the result we are after
    #————————————————————————————————————————————————
    # analytical expectation
    l = 0.0
    for i in 1:length(μα)
        l += (μα[i]*μα[i] + Σα[i,i]) * μf[i]^2 / S[i,i]
    end
    #————————————————————————————————————————————————

    # numerical expectation vs analytical
    mean(@showprogress map(g, [rand(Apdf) for _ in 1:numsamples])), l

end


text:
# # [Single-rate Explicit Timestepping](@id Single-rate-Explicit-Timestepping)

# In this tutorial, we shall explore the use of explicit Runge-Kutta
# methods for the solution of nonautonomous (or non time-invariant) equations.
# For our model problem, we shall reuse the rising thermal bubble
# tutorial. See its [tutorial page](@ref Rising-Thermal-Bubble-Configuration)
# for details on the model and parameters. For the purposes of this tutorial,
# we will only run the experiment for a total of 100 simulation seconds.

using ClimateMachine
const clima_dir = dirname(dirname(pathof(ClimateMachine)));
include(joinpath(
    clima_dir,
    "tutorials",
    "Numerics",
    "TimeStepping",
    "tutorial_risingbubble_config.jl",
))

FT = Float64;

# After discretizing the spatial terms in the equation, the semi-discretization
# of the governing equations have the form:

# ``
# \begin{aligned}
#     \frac{\mathrm{d} \boldsymbol{q}}{ \mathrm{d} t} &= M^{-1}\left(M S +
#     D^{T} M (F^{adv} + F^{visc}) + \sum_{f=1}^{N_f} L^T M_f(\widehat{F}^{adv} + \widehat{F}^{visc})
#     \right) \equiv \mathcal{T}(\boldsymbol{q}).
# \end{aligned}
# ``

# Referencing the canonical form introduced in [Time integration](@ref
# Time-integration) we have that in any explicit
# formulation ``\mathcal{F}(t, \boldsymbol{q}) \equiv 0`` and, in this particular
# forumlation, ``\mathcal{T}(t, \boldsymbol{q}) \equiv \mathcal{G}(t, \boldsymbol{q})``.

# The time step restriction for an explicit method must satisfy the stable
# [Courant number](https://en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition)
# for the specific time-integrator and must be selected from the following
# constraints
#
# ``
# \Delta t_{\mathrm{explicit}} = min \left( \frac{C \Delta x_i}{u_i + a}, \frac{C \Delta x_i^2}{\nu} \right)
# ``
#
# where ``C`` is the stable Courant number, ``u_i`` denotes the velocity components,
# ``a`` the speed of sound, ``\Delta x_i`` the grid spacing (non-uniform in case of
# spectral element methods) along the direction ``(x_1,x_2,x_3)``, and ``\nu`` the
# kinematic viscosity. The first term on the right is the time step condition
# due to the non-dissipative components, while the second term to the dissipation.
# For explicit time-integrators, we have to find the minimum time step that
# satisfies this condition along all three spatial directions.
#
# ## Runge-Kutta methods
#
# A single step of an ``s``-stage Runge-Kutta (RK) method for
# solving the resulting ODE problem presented above and can be
# expressed as the following:
#
# ```math
# \begin{align}
# 	\boldsymbol{q}^{n+1} = \boldsymbol{q}^n + \Delta t \sum_{i=1}^{s} b_i \mathcal{T}(\boldsymbol{Q}^i),
# \end{align}
# ```
#
# where ``\boldsymbol{\mathcal{T}}(\boldsymbol{Q}^i)`` is the evaluation of the
# right-hand side tendency at the stage value ``\boldsymbol{Q}^i``, defined at
# each stage of the RK method:
#
# ```math
# \begin{align}
# 	\boldsymbol{Q}^i = \boldsymbol{q}^{n} +
#     \Delta t \sum_{j=1}^{s} a_{i,j}
#     \mathcal{T}(\boldsymbol{Q}^j).
# \end{align}
# ```
#
# The first stage is initialized using the field at the previous time step:
# ``\boldsymbol{Q}^{1} \leftarrow \boldsymbol{q}^n``.
#
# In the above expressions, we define
# ``\boldsymbol{A} = \lbrace a_{i,j} \rbrace \in \mathbb{R}^{s\times s}``,
# ``\boldsymbol{b} = \lbrace b_i \rbrace \in \mathbb{R}^s``, and
# ``\boldsymbol{c} = \lbrace c_i \rbrace \in \mathbb{R}^s`` as the
# characteristic coefficients of a given RK method. This means we can
# associate any RK method with its so-called *Butcher tableau*:
#
# ```math
# \begin{align}
#     \begin{array}{c|c}
#         \boldsymbol{c} &\boldsymbol{A}\\
#         \hline
#         & \boldsymbol{b}^T
#         \end{array} =
#         \begin{array}{c|c c c c}
#         c_1 & a_{1,1} & a_{1,2} & \cdots & a_{1,s}\\
#         c_2 & a_{2,1} & a_{2,2} & \cdots & a_{2,s}\\
#         \vdots & \vdots & \vdots & \ddots & \vdots\\
#         c_s & a_{s,1} & a_{s,2} & \cdots & a_{s,s}\\
#         \hline
#         & b_1 & b_2 & \cdots & b_s
#     \end{array}.
# \end{align}
# ```
#
# The vector ``\boldsymbol{c}`` is often called the *consistency vector*,
# and is typically subjected to the row-sum condition:
#
# ``
# c_i = \sum_{j=1}^{s} a_{i,j}, \quad \forall i = 1, \cdots, s.
# ``
#
# This simplifies the order conditions for higher-order RK methods.
# For more information on general RK methods, we refer the interested reader
# to Ch. 5.2 of [Atkinson2011](@cite).
#
# ### [Low-storage Runge-Kutta (LSRK) methods](@id lsrk)
# `ClimateMachine.jl` contains the following low-storage methods:
#   - Forward Euler [`LowStorageRungeKutta2N`](@ref ClimateMachine.ODESolvers.LowStorageRungeKutta2N),
#   - A 5-stage 4th-order Runge-Kutta method of Carpenter and Kennedy [`LSRK54CarpenterKennedy`](@ref ClimateMachine.ODESolvers.LSRK54CarpenterKennedy)
#   - A 14-stage 4th-order Runge-Kutta method developed by Niegemann, Diehl, and Busch [`LSRK144NiegemannDiehlBusch`](@ref ClimateMachine.ODESolvers.LSRK144NiegemannDiehlBusch).
#
# To start, let's try using the 5-stage method: `LSRK54CarpenterKennedy`.

# As is the case for all explicit methods, we are limited by the fastest
# propogating waves described by our governing equations. In our case,
# these are the acoustic waves (with approximate wave speed given by the
# speed of sound of 343 m/s).
# For the rising bubble example used here, we use 4th order polynomials in
# a discontinuous Galerkin approximation, with a domain resolution of
# 125 meters in each spatial direction. This gives an effective
# minimanl nodal distance (distance between LGL nodes) of 86 meters
# over the entire mesh. Using the equation for the explcit time-step above,
# we can determine the ``\Delta t`` by specifying the desired Courant number
# ``C`` (denoted `CFL` in the code below).
# In our case, a heuristically determined value of 0.4 is used.

timeend = FT(100)
ode_solver =
    ClimateMachine.ExplicitSolverType(solver_method = LSRK54CarpenterKennedy)
CFL = FT(0.4)
run_simulation(ode_solver, CFL, timeend);

# What if we wish to take a larger timestep size?  We could
# try to increase the target Courant number, say ``C = 1.7``, and
# re-run the simulation. But this would break! In fact, in this case the numerical
# scheme would fall outside of the stability region and the simulation would crash.
# This occurs when the time step _exceeds_ the maximal stable time-step size
# of the method. For the 5-stage method, one can typically get away with using
# time-step sizes corresponding to a Courant number of ``C \approx 0.4`` but
# typically not much larger. In contrast, we can use an LSRK method with
# a larger stability region. Let's illustrate this by using the 14-stage method
# with a ``C = 1.7`` instead.

ode_solver = ClimateMachine.ExplicitSolverType(
    solver_method = LSRK144NiegemannDiehlBusch,
)
CFL = FT(1.7)
run_simulation(ode_solver, CFL, timeend);

# And it successfully completes. Currently, the 14-stage LSRK method
# `LSRK144NiegemannDiehlBusch` contains the largest stability region of the
# low-storage methods available in `ClimateMachine.jl`.

# ### [Strong Stability Preserving Runge--Kutta (SSPRK) methods](@id ssprk)

# Just as with the LSRK methods, the SSPRK methods are self-starting,
# with ``\boldsymbol{Q}^{1} = \boldsymbol{q}^n``, and stage-values are of the form

# ```math
# \begin{align}
#     \boldsymbol{Q}^{i+1} = a_{i,1} \boldsymbol{q}^n
#     + a_{i,2} \boldsymbol{Q}^{i}
#     + \Delta t b_i\mathcal{T}(\boldsymbol{Q}^{i})
# \end{align}
# ```
#
# with the value at the next step being the ``(N+1)``-th stage value
# ``\boldsymbol{q}^{n+1} = \boldsymbol{Q}^{(N+1)}``. This allows the updates
# to be performed with only three copies of the state vector
# (storing ``\boldsymbol{q}^n``, ``\boldsymbol{Q}^{i}`` and ``\mathcal{T}(\boldsymbol{Q}^{i})``).
# We illustrate here the use of a [`SSPRK33ShuOsher`](@ref ClimateMachine.ODESolvers.SSPRK33ShuOsher) method.

ode_solver = ClimateMachine.ExplicitSolverType(solver_method = SSPRK33ShuOsher)
CFL = FT(0.2)
run_simulation(ode_solver, CFL, timeend);

# ## References
# - [Shu1988](@cite)
# - [Heun1900](@cite)


text:
using BenchmarkTools
using SimpleValueGraphs, Graphs
using SimpleWeightedGraphs, MetaGraphs
using SNAPDatasets: loadsnap, snap_graphs
using SparseArrays
using Random: MersenneTwister

function load_graph_with_rand_weights(name::Symbol)

    if name ∈ keys(snap_graphs)
        g = squash(loadsnap(name))

        rng = MersenneTwister(0)
        G = is_directed(g) ? ValDiGraph : ValGraph
        return G(g; edgeval_types=(weight=Float64, ),
                 edgeval_init=(s, d) -> (weight=rand(rng, 1:255),))
    end

    error("no graph called :$name found.")
end

function graph_with_weights(::Type{SimpleGraph{T}}, ::Type{Matrix}, name) where {T}

    gv = load_graph_with_rand_weights(name)

    g = SimpleGraph{T}(adjacency_matrix(gv))
    matrix = Matrix(weights(gv))

    return (graph=g, weights=matrix)
end

function graph_with_weights(::Type{SimpleGraph{T}}, ::Type{SparseMatrixCSC}, name) where {T}

    gv = load_graph_with_rand_weights(name)

    g = SimpleGraph{T}(adjacency_matrix(gv))
    matrix = SparseMatrixCSC(weights(gv))

    return (graph=g, weights=matrix)
end

function graph_with_weights(::Type{SimpleWeightedGraph{T}}, name) where {T}

    gv = load_graph_with_rand_weights(name)

    return SimpleWeightedGraph{T}(weights(gv))
end

function graph_with_weights(::Type{MetaGraph{T}}, name) where {T}

    T == Int || error("Currently only MetaGraph{Int} is supported")

    g, w = graph_with_weights(SimpleGraph{T}, SparseMatrixCSC, name)
    gm = MetaGraph(g)
    for e in edges(gm)
        set_prop!(gm, e, :weight, w[src(e), dst(e)])
    end

    return gm
end

function graph_with_weights(::Type{ValGraph{T}}, name) where {T}

    g, w = graph_with_weights(SimpleGraph{T}, SparseMatrixCSC, name)
    return ValGraph(g, edgeval_types=(Float64,), edgeval_init=(s, d) -> (w[s,d],))
end

suite = BenchmarkGroup()

suite["prim_mst"] = BenchmarkGroup()
suite["prim_mst"]["facebook"] = BenchmarkGroup()

suite["prim_mst"]["facebook"]["SimpleGraph{Int} + Matrix"] =
    @benchmarkable prim_mst(g, m) setup=((g, m) = graph_with_weights(SimpleGraph{Int}, Matrix, :facebook_combined))
suite["prim_mst"]["facebook"]["SimpleGraph{Int} + SparseMatrixCSC"] =
    @benchmarkable prim_mst(g, m) setup=((g, m) = graph_with_weights(SimpleGraph{Int}, SparseMatrixCSC, :facebook_combined))
suite["prim_mst"]["facebook"]["SimpleWeightedGraph{Int}"] =
    @benchmarkable prim_mst(g) setup=(g = graph_with_weights(SimpleWeightedGraph{Int}, :facebook_combined))
suite["prim_mst"]["facebook"]["MetaGraph{Int}"] =
    @benchmarkable prim_mst(g) setup=(g = graph_with_weights(MetaGraph{Int}, :facebook_combined))
suite["prim_mst"]["facebook"]["ValGraph{Int}"] =
    @benchmarkable prim_mst(g) setup=(g = graph_with_weights(ValGraph{Int}, :facebook_combined))

suite["boruvka_mst"] = BenchmarkGroup()
suite["boruvka_mst"]["facebook"] = BenchmarkGroup()

suite["boruvka_mst"]["facebook"]["SimpleGraph{Int} + Matrix"] =
    @benchmarkable boruvka_mst(g, m) setup=((g, m) = graph_with_weights(SimpleGraph{Int}, Matrix, :facebook_combined))
suite["boruvka_mst"]["facebook"]["SimpleGraph{Int} + SparseMatrixCSC"] =
    @benchmarkable boruvka_mst(g, m) setup=((g, m) = graph_with_weights(SimpleGraph{Int}, SparseMatrixCSC, :facebook_combined))
suite["boruvka_mst"]["facebook"]["SimpleWeightedGraph{Int}"] =
    @benchmarkable boruvka_mst(g) setup=(g = graph_with_weights(SimpleWeightedGraph{Int}, :facebook_combined))
suite["boruvka_mst"]["facebook"]["MetaGraph{Int}"] =
    @benchmarkable boruvka_mst(g) setup=(g = graph_with_weights(MetaGraph{Int}, :facebook_combined))
suite["boruvka_mst"]["facebook"]["ValGraph{Int}"] =
    @benchmarkable boruvka_mst(g) setup=(g = graph_with_weights(ValGraph{Int}, :facebook_combined))

suite["dijkstra_shortest_paths"] = BenchmarkGroup()
suite["dijkstra_shortest_paths"]["facebook"] = BenchmarkGroup()

suite["dijkstra_shortest_paths"]["facebook"]["SimpleGraph{Int} + Matrix"] =
    @benchmarkable dijkstra_shortest_paths(g, 1, m) setup=((g, m) = graph_with_weights(SimpleGraph{Int}, Matrix, :facebook_combined))
suite["dijkstra_shortest_paths"]["facebook"]["SimpleGraph{Int} + SparseMatrixCSC"] =
    @benchmarkable dijkstra_shortest_paths(g, 1, m) setup=((g, m) = graph_with_weights(SimpleGraph{Int}, SparseMatrixCSC, :facebook_combined))
suite["dijkstra_shortest_paths"]["facebook"]["SimpleWeightedGraph{Int}"] =
    @benchmarkable dijkstra_shortest_paths(g, 1) setup=(g = graph_with_weights(SimpleWeightedGraph{Int}, :facebook_combined))
suite["dijkstra_shortest_paths"]["facebook"]["MetaGraph{Int}"] =
    @benchmarkable dijkstra_shortest_paths(g, 1) setup=(g = graph_with_weights(MetaGraph{Int}, :facebook_combined))
suite["dijkstra_shortest_paths"]["facebook"]["ValGraph{Int}"] =
    @benchmarkable dijkstra_shortest_paths(g, 1) setup=(g = graph_with_weights(ValGraph{Int}, :facebook_combined))

suite["Experimental.dijkstra_shortest_paths"] = BenchmarkGroup()
suite["Experimental.dijkstra_shortest_paths"]["facebook"] = BenchmarkGroup()

suite["Experimental.dijkstra_shortest_paths"]["facebook"]["SimpleGraph{Int} + Matrix"] =
    @benchmarkable SimpleValueGraphs.Experimental.dijkstra_shortest_paths(g, 1, m) setup=((g, m) = graph_with_weights(SimpleGraph{Int}, Matrix, :facebook_combined))
suite["Experimental.dijkstra_shortest_paths"]["facebook"]["SimpleGraph{Int} + SparseMatrixCSC"] =
    @benchmarkable SimpleValueGraphs.Experimental.dijkstra_shortest_paths(g, 1, m) setup=((g, m) = graph_with_weights(SimpleGraph{Int}, SparseMatrixCSC, :facebook_combined))
suite["Experimental.dijkstra_shortest_paths"]["facebook"]["SimpleWeightedGraph{Int}"] =
    @benchmarkable SimpleValueGraphs.Experimental.dijkstra_shortest_paths(g, 1) setup=(g = graph_with_weights(SimpleWeightedGraph{Int}, :facebook_combined))
suite["Experimental.dijkstra_shortest_paths"]["facebook"]["MetaGraph{Int}"] =
    @benchmarkable SimpleValueGraphs.Experimental.dijkstra_shortest_paths(g, 1) setup=(g = graph_with_weights(MetaGraph{Int}, :facebook_combined))
suite["Experimental.dijkstra_shortest_paths"]["facebook"]["ValGraph{Int}"] =
    @benchmarkable SimpleValueGraphs.Experimental.dijkstra_shortest_paths(g, 1) setup=(g = graph_with_weights(ValGraph{Int}, :facebook_combined))

if !isinteractive()
    tune!(suite, verbose = true)
    results = run(suite, verbose = true)
    println(results)
end


text:
"""
    sum_gp(first_term, ratio, num_terms)

Finds sum of n terms in a geometric progression

# Input parameters

- first_term : first term of the series
- raio      : common ratio between consecutive terms -> a2/a1 or a3/a2 or a4/a3
- num_terms  : number of terms in the series till which we count sum

# Example

```julia
sum_gp(1, 2, 10)    # 1023.0
sum_gp(1, 10, 5)    # 11111.0
sum_gp(0, 2, 10)    # 0.0
sum_gp(1, 0, 10)    # 1.0
sum_gp(1, 2, 0)     # -0.0
sum_gp(-1, 2, 10)   # -1023.0
sum_gp(1, -2, 10)   # -341.0
```

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee)
"""
function sum_gp(first_term, ratio, num_terms)
    # case where ratio is 1
    if ratio == 1
        return num_terms * first_term
    end
    # ormula for finding sum of n terms of a geometric progression
    return (first_term / (1 - ratio)) * (1 - ratio^num_terms)
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
"""
    SimState
"""
mutable struct SimState
    fields::Vector{Field}
end
SimState() = SimState(Vector{Field}())

function addfield!(simstate::SimState,field::Field)
    push!(simstate.fields,field)
end


text:
mutable struct Referencer <: Rewriter
    references::Dict{String, Union{Reference, Nothing}}
    modules
    doc
    refdocs
    path
    linkfn
end

Base.show(io::IO, referencer::Referencer) = print(io, "Referencer($(Tuple(referencer.modules)), $(length(referencer.references)) references)")

function Referencer(modules; path = p"REFERENCE")
    return Referencer(
        Dict{String, Union{Reference, Nothing}}(),
        modules,
        XNode(:body),
        Dict{String, XNode}(),
        Path(path),
        ref -> "$path/$ref",
    )
end

function rewritedoc(referencer::Referencer, ::AbstractPath, doc::XNode)
    doc = populatereferences!(
        referencer.references,
        doc,
        referencer.linkfn,
        referencer.modules)
    return doc
end

function createsources!(referencer::Referencer)
    refdoc = buildreference(referencer)
    docs = Dict{AbstractPath, XNode}()

    # Conditionally update/create reference page
    if refdoc != referencer.doc
        referencer.doc = refdoc
        docs[referencer.path] =  refdoc
    end

    # Conditionally update/create binding pages
    l = ReentrantLock()
    #Threads.@threads for ref in collect(values(referencer.references))
    for ref in [v for v in values(referencer.references)]
        doc = get(referencer.refdocs, ref.fullname, nothing)
        if isnothing(doc)
            newdoc = buildreferencepage(ref)
            lock(l) do
                referencer.refdocs[ref.fullname] = newdoc
                docs[Path(referencer.linkfn(ref.fullname))] = newdoc
            end
        end
    end

    return docs
end

function reset!(referencer::Referencer)
    referencer.references = Dict{String, Union{Reference, Nothing}}()
    referencer.doc = XNode(:body)
    referencer.refdocs = Dict{String, XNode}()
    return
end

function buildreference(referencer)
    refs = referencer.references
    children = [XNode(:h1, [XLeaf("Reference")])]

    ms = sort(unique([ref.m for ref in values(referencer.references)]), by = string)

    for m in referencer.modules
        children = vcat(children,
            buildmodulereference(
                string(m),
                [ref for (k, ref) in refs if ref.m == m],
                referencer.linkfn))
    end
    return XNode(:body, children)
end

function buildmodulereference(mname, refs, linkfn)
    refs = sortrefs(refs)

    heading = XNode(:h2, [XLeaf(mname)])

    links = [
        XNode(
            :a,
            Dict(:href => linkfn(r.fullname)),
            [XNode(:code, [XLeaf(string(r.identifier))])]
        )
        for r in refs
    ]
    return [
        heading,
        XNode(:ul, [XNode(:li, [link]) for link in links])
    ]
end

function sortrefs(refs)
    return sort(refs, by = r -> (orderreftype(r.kind), string(r.identifier)))
end

orderreftype(reftype::Symbol) = ORDERREFTYPE[reftype]

const ORDERREFTYPE = Dict(
    :module => 1,
    :const => 2,
    :type => 3,
    :ptype => 4,
    :struct => 5,
    :pstruct => 6,
    :function => 7,

)

function buildreferencepage(ref::Reference)
    docs = getdocs(ref.m, ref.identifier)
    xdocs = if isnothing(docs)
        XNode(:p, XTree[
            XLeaf("No documentation found for "),
            XNode(:code, [XLeaf(ref.fullname)]),
            XLeaf(".")
            ])
    else
        convert(XTree, docs)
    end
    if !isnothing(ref.identifier) && ref.kind in (:struct, :function)
        push!(xdocs.children, XNode(:h2, [XLeaf("Methods")]))
        push!(xdocs.children, XLeaf(methods(getfield(ref.m, ref.identifier))))
    end
    return XNode(
        :body,
        [XNode(:article,
            [
                XNode(:h1, [XNode(:code, [XLeaf(ref.fullname)])]),
                xdocs,
            ],
        )]
    )
end


text:
# device management and properties

export AbstractGPUDevice

abstract type AbstractGPUDevice end

"""
    device(A::AbstractArray)

Gets the device associated to the Array `A`
"""
device(A::AbstractArray) = error("This array is not a GPU array") # COV_EXCL_LINE

"""
Hardware threads of device
"""
threads(::AbstractGPUDevice) = error("Not implemented") # COV_EXCL_LINE


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

using Test

let cmd = `$(Base.julia_cmd()) --depwarn=error --rr-detach --startup-file=no threads_exec.jl`
    for test_nthreads in (1, 2, 4, 4) # run once to try single-threaded mode, then try a couple times to trigger bad races
        new_env = copy(ENV)
        new_env["JULIA_NUM_THREADS"] = string(test_nthreads)
        run(pipeline(setenv(cmd, new_env), stdout = stdout, stderr = stderr))
    end
end

# issue #34415 - make sure external affinity settings work
if Sys.islinux()
    const SYS_rrcall_check_presence = 1008
    global running_under_rr() = 0 == ccall(:syscall, Int,
        (Int, Int, Int, Int, Int, Int, Int),
        SYS_rrcall_check_presence, 0, 0, 0, 0, 0, 0)
    if Sys.CPU_THREADS > 1 && Sys.which("taskset") !== nothing && !running_under_rr()
        run_with_affinity(spec) = readchomp(`taskset -c $spec $(Base.julia_cmd()) -e "run(\`taskset -p \$(getpid())\`)"`)
        @test endswith(run_with_affinity("1"), "2")
        @test endswith(run_with_affinity("0,1"), "3")
    end
end

# issue #34769
function idle_callback(handle)
    idle = @Base.handle_as handle UvTestIdle
    if idle.active
        idle.count += 1
        if idle.count == 1
            # We want to hit the case where we're allowing
            # the thread to go to sleep, which only happens
            # after some default amount of time (DEFAULT_THREAD_SLEEP_THRESHOLD)
            # so spend that amount of time here.
            Libc.systemsleep(0.004)
        elseif idle.count >= 10
            lock(idle.cond)
            try
                notify(idle.cond, true)
            finally
                unlock(idle.cond)
            end
            idle.active = false
        end
    end
    nothing
end

mutable struct UvTestIdle
    handle::Ptr{Cvoid}
    cond::Base.ThreadSynchronizer
    isopen::Bool
    active::Bool
    count::Int

    function UvTestIdle()
        this = new(Libc.malloc(Base._sizeof_uv_idle), Base.ThreadSynchronizer(), true, false, 0)
        Base.iolock_begin()
        Base.associate_julia_struct(this.handle, this)
        err = ccall(:uv_idle_init, Cint, (Ptr{Cvoid}, Ptr{Cvoid}),
            Base.eventloop(), this.handle)
        if err != 0
            Libc.free(this.handle)
            this.handle = C_NULL
            throw(_UVError("uv_idle_init", err))
        end
        err = ccall(:uv_idle_start, Cint, (Ptr{Cvoid}, Ptr{Cvoid}),
            this.handle, @cfunction(idle_callback, Cvoid, (Ptr{Cvoid},)))
        if err != 0
            Libc.free(this.handle)
            this.handle = C_NULL
            throw(_UVError("uv_idle_start", err))
        end
        finalizer(Base.uvfinalize, this)
        Base.iolock_end()
        return this
    end
end
Base.unsafe_convert(::Type{Ptr{Cvoid}}, idle::UvTestIdle) = idle.handle

function Base.uvfinalize(t::UvTestIdle)
    Base.iolock_begin()
    Base.lock(t.cond)
    try
        if t.handle != C_NULL
            Base.disassociate_julia_struct(t.handle) # not going to call the usual close hooks
            if t.isopen
                t.isopen = false
                ccall(:jl_close_uv, Cvoid, (Ptr{Cvoid},), t)
            end
            t.handle = C_NULL
            notify(t.cond, false)
        end
    finally
        unlock(t.cond)
    end
    Base.iolock_end()
    nothing
end

function Base.close(idle::UvTestIdle)
    Base.uvfinalize(idle)
end

function Base.wait(idle::UvTestIdle)
    Base.iolock_begin()
    Base.preserve_handle(idle)
    Base.lock(idle.cond)
    try
        idle.active = true
        Base.iolock_end()
        wait(idle.cond)
    finally
        Base.unlock(idle.cond)
        Base.iolock_begin()
        Base.unpreserve_handle(idle)
        Base.iolock_end()
    end
end

# Spawn another process as a watchdog. If this test fails, it'll unrecoverably
# hang in the event loop. Another process needs to kill it
cmd = """
    @async (Base.wait_readnb(stdin, 1); exit())
    sleep(100)
    isopen(stdin) || exit()
    println(stderr, "ERROR: Killing threads test due to watchdog expiry")
    ccall(:uv_kill, Cint, (Cint, Cint), $(getpid()), Base.SIGTERM)
"""
proc = open(pipeline(`$(Base.julia_cmd()) -e $cmd`; stderr=stderr); write=true)

let idle=UvTestIdle()
    wait(idle)
    close(idle)
end

using Base.Threads
@threads for i = 1:1
    let idle=UvTestIdle()
        wait(idle)
        close(idle)
    end
end

@test process_running(proc)

# We don't need the watchdog anymore
close(proc.in)

# https://github.com/JuliaLang/julia/pull/42973
Sys.islinux() && @testset "spawn and wait *a lot* of tasks in @profile" begin
    # Not using threads_exec.jl for better isolation, reproducibility, and a
    # tighter timeout.
    script = "profile_spawnmany_exec.jl"
    cmd = `$(Base.julia_cmd()) --depwarn=error --rr-detach --startup-file=no $script`
    @testset for n in [20000, 200000, 2000000]
        proc = run(ignorestatus(setenv(cmd, "NTASKS" => n; dir = @__DIR__)); wait = false)
        done = Threads.Atomic{Bool}(false)
        timeout = false
        timer = Timer(100) do _
            timeout = true
            for sig in [Base.SIGTERM, Base.SIGHUP, Base.SIGKILL]
                for _ in 1:1000
                    kill(proc, sig)
                    if done[]
                        if sig != Base.SIGTERM
                            @warn "Terminating `$script` required signal $sig"
                        end
                        return
                    end
                    sleep(0.001)
                end
            end
        end
        try
            wait(proc)
        finally
            done[] = true
            close(timer)
        end
        @test success(proc)
        @test !timeout
    end
end

max_stars_repo_name: Seelengrab/julia

max_issues_repo_name: Seelengrab/julia

max_forks_repo_name: Seelengrab/julia


text:
using FinancialAssets
using Test

@testset "FinancialAssets.jl" begin
    # Write your tests here.
end

max_stars_repo_name: robfs/FinancialAssets

max_issues_repo_name: robfs/FinancialAssets

max_forks_repo_name: robfs/FinancialAssets


text:
using RegistryActivity
using Test
using LibGit2
using Dates

origin = let
    general_path = RegistryActivity.general_registry()
    try
        # If the local General registry is a git repository, use it...
        LibGit2.GitRepo(general_path)
        general_path
    catch e
        if e isa LibGit2.Error.GitError && e.code == LibGit2.Error.ENOTFOUND
            # ...otherwise clone directly from upstream
            "https://github.com/JuliaRegistries/General.git"
        else
            # If there are other issues, rethrow the error
            rethrow(e)
        end
    end
end

@testset "RegistryActivity" begin
    registry = RegistryActivity.clone_registry(; origin)

    m, p, v = registry_activity(registry;
                                start_month=Date(2019, 1),
                                end_month=Date(2019, 12),
                                filter=p->(occursin(r"_jll$", last(p)["name"])),
                                )
    @test m == Date(2019, 1):Month(1):Date(2019, 12)
    @test p == [0, 0, 0, 0, 0, 0, 0, 9, 58, 143, 176, 209]
    @test v == [0, 0, 0, 0, 0, 0, 0, 12, 98, 221, 274, 343]

    m, p, v = registry_activity(registry;
                                start_month=Date(2020, 1),
                                end_month=Date(2020, 12),
                                )
    @test m == Date(2020, 1):Month(1):Date(2020, 12)
    @test p == [2774, 2885, 3013, 3132, 3280, 3389, 3541, 3680, 3827, 3981, 4091, 4215]
    @test v == [13621, 15036, 16381, 17914, 19426, 20759, 22082, 23567, 24969, 26360, 27674, 29067]
end


text:
using CLIMAParameters.Planet: R_d, cv_d, T_0, e_int_v0, e_int_i0

"""
    linearized_air_pressure(ρ, ρe_tot, ρe_pot, ρq_tot=0, ρq_liq=0, ρq_ice=0)

The air pressure, linearized around a dry rest state, from the equation of state
(ideal gas law) where:

 - `ρ` (moist-)air density
 - `ρe_tot` total energy density
 - `ρe_pot` potential energy density
and, optionally,
 - `ρq_tot` total water density
 - `ρq_liq` liquid water density
 - `ρq_ice` ice density
"""
function linearized_air_pressure(
    param_set::AbstractParameterSet,
    ρ::FT,
    ρe_tot::FT,
    ρe_pot::FT,
    ρq_tot::FT = FT(0),
    ρq_liq::FT = FT(0),
    ρq_ice::FT = FT(0),
) where {FT <: Real, PS}
    _R_d::FT = R_d(param_set)
    _cv_d::FT = cv_d(param_set)
    _T_0::FT = T_0(param_set)
    _e_int_v0::FT = e_int_v0(param_set)
    _e_int_i0::FT = e_int_i0(param_set)
    return ρ * _R_d * _T_0 +
           _R_d / _cv_d * (
        ρe_tot - ρe_pot - (ρq_tot - ρq_liq) * _e_int_v0 +
        ρq_ice * (_e_int_i0 + _e_int_v0)
    )
end

@inline linearized_pressure(atmos, state::Vars, aux::Vars) =
    linearized_pressure(
        moisture_model(atmos),
        parameter_set(atmos),
        atmos.orientation,
        state,
        aux,
    )

@inline function linearized_pressure(
    ::DryModel,
    param_set::AbstractParameterSet,
    orientation::Orientation,
    state::Vars,
    aux::Vars,
)
    ρe_pot = state.ρ * gravitational_potential(orientation, aux)
    return linearized_air_pressure(param_set, state.ρ, state.energy.ρe, ρe_pot)
end
@inline function linearized_pressure(
    ::EquilMoist,
    param_set::AbstractParameterSet,
    orientation::Orientation,
    state::Vars,
    aux::Vars,
)
    ρe_pot = state.ρ * gravitational_potential(orientation, aux)
    linearized_air_pressure(
        param_set,
        state.ρ,
        state.energy.ρe,
        ρe_pot,
        state.moisture.ρq_tot,
    )
end
@inline function linearized_pressure(
    ::NonEquilMoist,
    param_set::AbstractParameterSet,
    orientation::Orientation,
    state::Vars,
    aux::Vars,
)
    ρe_pot = state.ρ * gravitational_potential(orientation, aux)
    linearized_air_pressure(
        param_set,
        state.ρ,
        state.energy.ρe,
        ρe_pot,
        state.moisture.ρq_tot,
        state.moisture.ρq_liq,
        state.moisture.ρq_ice,
    )
end

abstract type AtmosLinearModel <: BalanceLaw end

"""
    vars_state(m::AtmosLinearModel, ::Prognostic, FT)

Conserved state variables (prognostic variables).

!!! warning

    `AtmosLinearModel` state ordering must be a contiguous subset of the initial
    state of `AtmosModel` since a shared state is used.
"""
function vars_state(lm::AtmosLinearModel, st::Prognostic, FT)
    @vars begin
        ρ::FT
        ρu::SVector{3, FT}
        energy::vars_state(energy_model(lm.atmos), st, FT)
        turbulence::vars_state(turbulence_model(lm.atmos), st, FT)
        hyperdiffusion::vars_state(hyperdiffusion_model(lm.atmos), st, FT)
        moisture::vars_state(moisture_model(lm.atmos), st, FT)
    end
end
vars_state(lm::AtmosLinearModel, st::Auxiliary, FT) =
    vars_state(lm.atmos, st, FT)
vars_state(lm::AtmosLinearModel, ::Primitive, FT) =
    vars_state(lm, Prognostic(), FT)

function update_auxiliary_state!(
    dg::DGModel,
    lm::AtmosLinearModel,
    Q::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    return false
end
function flux_second_order!(
    lm::AtmosLinearModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    aux::Vars,
    t::Real,
)
    nothing
end
integral_load_auxiliary_state!(
    lm::AtmosLinearModel,
    integ::Vars,
    state::Vars,
    aux::Vars,
) = nothing
integral_set_auxiliary_state!(lm::AtmosLinearModel, aux::Vars, integ::Vars) =
    nothing
reverse_integral_load_auxiliary_state!(
    lm::AtmosLinearModel,
    integ::Vars,
    state::Vars,
    aux::Vars,
) = nothing
reverse_integral_set_auxiliary_state!(
    lm::AtmosLinearModel,
    aux::Vars,
    integ::Vars,
) = nothing
function wavespeed(
    lm::AtmosLinearModel,
    nM,
    state::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    ref = aux.ref_state
    return soundspeed_air(parameter_set(lm.atmos), ref.T)
end

boundary_conditions(atmoslm::AtmosLinearModel) =
    (AtmosBC(atmoslm.atmos.physics), AtmosBC(atmoslm.atmos.physics))

function boundary_state!(
    nf::NumericalFluxFirstOrder,
    bc,
    atmoslm::AtmosLinearModel,
    state⁺,
    aux⁺,
    n,
    state⁻,
    aux⁻,
    t,
    state_int⁻,
    aux_int⁻,
)

    args = (; aux⁺, state⁻, aux⁻, t, n, state_int⁻, aux_int⁻)

    atmos_boundary_state!(nf, bc, atmoslm, state⁺, args)
end
function boundary_state!(
    nf::NumericalFluxSecondOrder,
    bc,
    atmoslm::AtmosLinearModel,
    _...,
)
    nothing
end
init_state_auxiliary!(
    lm::AtmosLinearModel,
    aux::MPIStateArray,
    grid,
    direction,
) = nothing
init_state_prognostic!(
    lm::AtmosLinearModel,
    state::Vars,
    aux::Vars,
    localgeo,
    t,
) = nothing

struct AtmosAcousticLinearModel{M} <: AtmosLinearModel
    atmos::M
    function AtmosAcousticLinearModel(atmos::M) where {M}
        if reference_state(atmos) === NoReferenceState()
            error("AtmosAcousticLinearModel needs a model with a reference state")
        end
        new{M}(atmos)
    end
end

function flux_first_order!(
    lm::AtmosLinearModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    tend = Flux{FirstOrder}()
    _args = (; state, aux, t, direction)

    # For some reason, we cannot call precompute, because
    # sometimes `state.ρ` is somehow `0`, which results in
    # `e_int = Inf` -> failed saturation adjustment.
    # TODO: look into this
    args = _args
    # args = merge(_args, (precomputed = precompute(lm.atmos, _args, tend),))
    flux.ρ = Σfluxes(Mass(), eq_tends(Mass(), lm, tend), lm, args)
    flux.ρu = Σfluxes(Momentum(), eq_tends(Momentum(), lm, tend), lm, args)
    flux.energy.ρe = Σfluxes(Energy(), eq_tends(Energy(), lm, tend), lm, args)
    nothing
end

struct AtmosAcousticGravityLinearModel{M} <: AtmosLinearModel
    atmos::M
    function AtmosAcousticGravityLinearModel(atmos::M) where {M}
        if reference_state(atmos) === NoReferenceState()
            error("AtmosAcousticGravityLinearModel needs a model with a reference state")
        end
        new{M}(atmos)
    end
end

function source!(
    lm::AtmosLinearModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    ::NTuple{1, Dir},
) where {Dir <: Direction}

    tend = Source()
    _args = (; state, aux, t, direction = Dir, diffusive)

    # For some reason, we cannot call precompute, because
    # sometimes `state.ρ` is somehow `0`, which results in
    # `e_int = Inf` -> failed saturation adjustment.
    # TODO: look into this
    args = _args
    # args = merge(_args, (precomputed = precompute(lm.atmos, _args, tend),))

    # Sources for the linear atmos model only appear in the momentum equation
    source.ρu = Σsources(Momentum(), eq_tends(Momentum(), lm, tend), lm, args)
    nothing
end

function numerical_flux_first_order!(
    numerical_flux::RoeNumericalFlux,
    balance_law::AtmosLinearModel,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
    direction,
) where {S, A}
    @assert moisture_model(balance_law.atmos) isa DryModel

    numerical_flux_first_order!(
        CentralNumericalFluxFirstOrder(),
        balance_law,
        fluxᵀn,
        normal_vector,
        state_prognostic⁻,
        state_auxiliary⁻,
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
        direction,
    )

    atmos = balance_law.atmos
    param_set = parameter_set(atmos)

    ρu⁻ = state_prognostic⁻.ρu

    ref_ρ⁻ = state_auxiliary⁻.ref_state.ρ
    ref_ρe⁻ = state_auxiliary⁻.ref_state.ρe
    ref_T⁻ = state_auxiliary⁻.ref_state.T
    ref_p⁻ = state_auxiliary⁻.ref_state.p
    ref_h⁻ = (ref_ρe⁻ + ref_p⁻) / ref_ρ⁻
    ref_c⁻ = soundspeed_air(param_set, ref_T⁻)

    pL⁻ = linearized_pressure(atmos, state_prognostic⁻, state_auxiliary⁻)

    ρu⁺ = state_prognostic⁺.ρu

    ref_ρ⁺ = state_auxiliary⁺.ref_state.ρ
    ref_ρe⁺ = state_auxiliary⁺.ref_state.ρe
    ref_T⁺ = state_auxiliary⁺.ref_state.T
    ref_p⁺ = state_auxiliary⁺.ref_state.p
    ref_h⁺ = (ref_ρe⁺ + ref_p⁺) / ref_ρ⁺
    ref_c⁺ = soundspeed_air(param_set, ref_T⁺)

    pL⁺ = linearized_pressure(atmos, state_prognostic⁺, state_auxiliary⁺)

    # not sure if arithmetic averages are a good idea here
    h̃ = (ref_h⁻ + ref_h⁺) / 2
    c̃ = (ref_c⁻ + ref_c⁺) / 2

    ΔpL = pL⁺ - pL⁻
    Δρuᵀn = (ρu⁺ - ρu⁻)' * normal_vector

    fluxᵀn.ρ -= ΔpL / 2c̃
    fluxᵀn.ρu -= c̃ * Δρuᵀn * normal_vector / 2
    fluxᵀn.energy.ρe -= h̃ * ΔpL / 2c̃
end

function numerical_flux_first_order!(
    ::HLLCNumericalFlux,
    balance_law::AtmosLinearModel,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
    direction,
) where {S, A}

    # There is no intermediate speed for the AtmosLinearModel.
    # As a result, HLLC simplifies to Rusanov.
    numerical_flux_first_order!(
        RusanovNumericalFlux(),
        balance_law,
        fluxᵀn,
        normal_vector,
        state_prognostic⁻,
        state_auxiliary⁻,
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
        direction,
    )
end
function numerical_flux_first_order!(
    numerical_flux::RoeNumericalFluxMoist,
    balance_law::AtmosLinearModel,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
    direction,
) where {S, A}

    numerical_flux_first_order!(
        CentralNumericalFluxFirstOrder(),
        balance_law,
        fluxᵀn,
        normal_vector,
        state_prognostic⁻,
        state_auxiliary⁻,
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
        direction,
    )

    atmos = balance_law.atmos
    param_set = parameter_set(atmos)
    FT = eltype(aux)
    ρu⁻ = state_prognostic⁻.ρu

    ref_ρ⁻ = state_auxiliary⁻.ref_state.ρ
    ref_ρe⁻ = state_auxiliary⁻.ref_state.ρe
    ref_T⁻ = state_auxiliary⁻.ref_state.T
    ref_q⁻ = state_auxiliary⁻.ref_state.ρq_tot / ref_ρ⁻
    ref_qliq⁻ = state_auxiliary⁻.ref_state.ρq_liq / ref_ρ⁻
    ref_qice⁻ = state_auxiliary⁻.ref_state.ρq_ice / ref_ρ⁻
    ref_p⁻ = state_auxiliary⁻.ref_state.p
    ref_q_pt⁻ = PhasePartition(ref_q⁻, ref_qliq⁻, ref_qice⁻)
    _R_m⁻ = gas_constant_air(param_set, ref_q_pt⁻)
    ref_h⁻ = total_specific_enthalpy(ref_ρe⁻, _R_m⁻, ref_T⁻)

    ref_c⁻ = soundspeed_air(param_set, ref_T⁻, ref_q_pt⁻)

    pL⁻ = linearized_pressure(atmos, state_prognostic⁻, state_auxiliary⁻)

    ρu⁺ = state_prognostic⁺.ρu

    ref_ρ⁺ = state_auxiliary⁺.ref_state.ρ
    ref_ρe⁺ = state_auxiliary⁺.ref_state.ρe
    ref_T⁺ = state_auxiliary⁺.ref_state.T
    ref_q⁺ = state_auxiliary⁺.ref_state.ρq_tot / ref_ρ⁺
    ref_qliq⁺ = state_auxiliary⁺.ref_state.ρq_liq / ref_ρ⁺
    ref_qice⁺ = state_auxiliary⁺.ref_state.ρq_ice / ref_ρ⁺
    ref_p⁺ = state_auxiliary⁺.ref_state.p
    ref_q_pt⁺ = PhasePartition(ref_q⁺, ref_qliq⁺, ref_qice⁺)
    _R_m⁺ = gas_constant_air(param_set, ref_q_pt⁺)
    ref_h⁺ = total_specific_enthalpy(ref_ρe⁺, _R_m⁺, ref_T⁺)
    ref_c⁺ = soundspeed_air(param_set, ref_T⁺, ref_q_pt⁺)
    pL⁺ = linearized_pressure(atmos, state_prognostic⁺, state_auxiliary⁺)

    # not sure if arithmetic averages are a good idea here
    ref_h̃ = (ref_h⁻ + ref_h⁺) / 2
    ref_c̃ = (ref_c⁻ + ref_c⁺) / 2
    ref_qt = (ref_q⁺ + ref_q⁻) / 2

    ΔpL = pL⁺ - pL⁻
    Δρuᵀn = (ρu⁺ - ρu⁻)' * normal_vector

    _T_0::FT = T_0(param_set)
    _e_int_v0::FT = e_int_v0(param_set)
    _cv_d::FT = cv_d(param_set)
    Φ = gravitational_potential(balance_law.atmos, state_auxiliary⁻)
    # guaranteed to be random
    ω = FT(π) / 3
    δ = FT(π) / 5
    random_unit_vector = SVector(sin(ω) * cos(δ), cos(ω) * cos(δ), sin(δ))
    # tangent space basis
    τ1 = random_unit_vector × normal_vector
    τ2 = τ1 × normal_vector

    ũc̃⁻ = ref_c̃ * normal_vector
    ũc̃⁺ = -ref_c̃ * normal_vector
    Λ = SDiagonal(
        abs(0 - ref_c̃),
        abs(0),
        abs(0),
        abs(0),
        abs(0 + ref_c̃),
        abs(0),
    )
    M = hcat(
        SVector(1, ũc̃⁺[1], ũc̃⁺[2], ũc̃⁺[3], ref_h̃, ref_qt),
        SVector(0, τ1[1], τ1[2], τ1[3], 0, 0),
        SVector(0, τ2[1], τ2[2], τ2[3], 0, 0),
        SVector(1, 0, 0, 0, Φ - _T_0 * _cv_d, 0),
        SVector(1, ũc̃⁻[1], ũc̃⁻[2], ũc̃⁻[3], ref_h̃, ref_qt),
        SVector(0, 0, 0, 0, _e_int_v0, 1),
    )
    Δρ = state_prognostic⁺.ρ - state_prognostic⁻.ρ
    Δρu = ρu⁺ - ρu⁻
    Δρe = state_prognostic⁺.energy.ρe - state_prognostic⁻.energy.ρe
    Δρq_tot =
        state_prognostic⁺.moisture.ρq_tot - state_prognostic⁻.moisture.ρq_tot
    Δstate = SVector(Δρ, Δρu[1], Δρu[2], Δρu[3], Δρe, Δρq_tot)

    parent(fluxᵀn) .-= M * Λ * (M \ Δstate) / 2
end

function numerical_flux_first_order!(
    ::LMARSNumericalFlux,
    balance_law::AtmosLinearModel,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
    direction,
) where {S, A}

    numerical_flux_first_order!(
        RusanovNumericalFlux(),
        balance_law,
        fluxᵀn,
        normal_vector,
        state_prognostic⁻,
        state_auxiliary⁻,
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
        direction,
    )
end


text:
using Bijectors

using Combinatorics
using DistributionsAD
using FiniteDiff
using ForwardDiff
using ReverseDiff
using Tracker
using Zygote

using Random, LinearAlgebra, Test

using Bijectors: Log, Exp, Shift, Scale, Logit, SimplexBijector, PDBijector, Permute,
    PlanarLayer, RadialLayer, Stacked, TruncatedBijector, ADBijector

using DistributionsAD: TuringUniform, TuringMvNormal, TuringMvLogNormal,
    TuringPoissonBinomial

const is_TRAVIS = haskey(ENV, "TRAVIS")
const GROUP = get(ENV, "GROUP", "All")

if GROUP == "All" || GROUP == "Interface"
    include("interface.jl")
    include("transform.jl")
    include("norm_flows.jl")
    include("bijectors/permute.jl")
end

if !is_TRAVIS && (GROUP == "All" || GROUP == "AD")
    include("ad/utils.jl")
    include("ad/distributions.jl")
end


text:
using PyPlot

N = 100
chord = 1.0
ds = .26
halfspan = 26
Vinf = 25.85 #(props are removed, so constant Vinfeff = Vinf)
Vinfeff = Vinf
Sref = 26*2

y = linspace(0,halfspan,N)
ds_ell = y[2]-y[1]
cl0 = 1.2
cl = cl0*sqrt.(1-(y/halfspan).^2)

clave_ell = mean(cl)

clopt = [1.18335, 1.18348, 1.18367, 1.1839, 1.18416, 1.18444, 1.18473, 1.18503, 1.18534, 1.18566, 1.18597, 1.18629, 1.18661, 1.18693, 1.18725, 1.18756, 1.18787, 1.18818, 1.18849, 1.18879, 1.18909, 1.18938, 1.18966, 1.18994, 1.19022, 1.19048, 1.19074, 1.191, 1.19124, 1.19148, 1.19171, 1.19193, 1.19214, 1.19234, 1.19253, 1.19271, 1.19287, 1.19303, 1.19317, 1.1933, 1.19342, 1.19352, 1.19361, 1.19368, 1.19373, 1.19376, 1.19378, 1.19377, 1.19375, 1.1937, 1.19362, 1.19352, 1.1934, 1.19324, 1.19305, 1.19283, 1.19257, 1.19228, 1.19194, 1.19156, 1.19113, 1.19065, 1.19011, 1.18951, 1.18885, 1.18811, 1.18729, 1.18639, 1.18539, 1.18429, 1.18307, 1.18172, 1.18023, 1.17857, 1.17674, 1.17471, 1.17246, 1.16994, 1.16713, 1.16399, 1.16046, 1.15648, 1.15198, 1.14687, 1.14104, 1.13435, 1.12662, 1.11765, 1.10713, 1.09472, 1.07991, 1.06206, 1.04025, 1.01321, 0.979061, 0.934965, 0.876372, 0.795405, 0.676626, 0.481338]
y2 = linspace(0,halfspan,length(clopt))
clopt_ave = mean(clopt)

CL_ell = 2*sum(cl.*Vinfeff.*chord.*ds_ell)/(Vinf.*Sref)
CL_opt = 2*sum(clopt.*Vinfeff.*chord.*ds)/(Vinf.*Sref)

PyPlot.figure()
PyPlot.plot(y,cl,"b-",label = "Elliptical cl: CL = $(round.(CL_ell,3))")
PyPlot.plot(y,ones(cl)*clave_ell,"b--",label = "Elliptical cl Average")
PyPlot.plot(y2,clopt,"r-",label = "FB results cl: CL = $(round.(CL_opt,3))")
PyPlot.plot(y2,ones(clopt)*clopt_ave,"r--",label = "FB results cl Average")
PyPlot.xlabel("Halfspan")
PyPlot.ylabel("cl")
PyPlot.legend(loc = "best")


text:
module SnakeGames

include("snake_game.jl")
include("interaction.jl")

end


text:
using Test

using MathOptInterface
const MOI = MathOptInterface

include("dummy.jl")

"""
    MutLessThan{T<:Real} <: MOI.AbstractScalarSet

A mutable `LessThan`-like set to test `copy` of indicator set
"""
mutable struct MutLessThan{T<:Real} <: MOI.AbstractScalarSet
    upper::T
    MutLessThan(v::T) where {T<:Real} = new{T}(v)
end

Base.copy(mlt::MutLessThan) = MutLessThan(Base.copy(mlt.upper))

@testset "Sets" begin
    @testset "==" begin
        # By default, `==` redirects to `===`, it works for bits type
        # but not for `BigInt`s. We define functions creating different
        # instances so that `a() !== a()`.
        a() = big(1)
        b() = big(2)
        @test a() !== a()
        @test b() !== b()
        for S in [MOI.LessThan, MOI.GreaterThan, MOI.EqualTo, MOI.PowerCone, MOI.DualPowerCone]
            @test S(a()) == S(a())
            @test S(a()) != S(b())
            @test S(b()) == S(b())
            @test S(b()) != S(a())
        end
        for S in [MOI.Interval, MOI.Semicontinuous, MOI.Semiinteger]
            @test S(a(), b()) == S(a(), b())
            @test S(a(), b()) != S(b(), a())
            @test S(a(), b()) != S(b(), b())
            @test S(a(), b()) != S(a(), a())
            @test S(a(), a()) != S(b(), b())
            @test S(a(), a()) == S(a(), a())
        end
        S = MOI.IndicatorSet
        A() = MOI.LessThan(a())
        B() = MOI.LessThan(b())
        @test S{MOI.ACTIVATE_ON_ZERO}(A()) == S{MOI.ACTIVATE_ON_ZERO}(A())
        @test S{MOI.ACTIVATE_ON_ZERO}(A()) != S{MOI.ACTIVATE_ON_ONE}(A())
        @test S{MOI.ACTIVATE_ON_ZERO}(A()) != S{MOI.ACTIVATE_ON_ZERO}(B())
        @test S{MOI.ACTIVATE_ON_ONE}(A()) != S{MOI.ACTIVATE_ON_ONE}(B())
    end
    @testset "Copy" begin
        @testset "for $S" for S in [MOI.SOS1, MOI.SOS2]
            s = S([1.0])
            s_copy = copy(s)
            s_copy.weights[1] = 2.0
            @test s.weights[1] == 1.0
        end
        @testset "IndicatorSet" begin
            s1 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(4.0))
            s2 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO}(MOI.GreaterThan(4.0))
            s1_copy = copy(s1)
            s2_copy = copy(s2)
            @test s1_copy isa MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}
            @test s1 == s1_copy
            @test s2_copy isa MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO}
            @test s2 == s2_copy
            s3 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO}(MutLessThan(4.0))
            s3_copy = copy(s3)
            @test s3.set.upper ≈ 4.0
            s3_copy.set.upper = 5.0
            @test s3.set.upper ≈ 4.0
            @test s3_copy.set.upper ≈ 5.0
        end
    end
    @testset "Broadcast" begin
        model = DummyModelWithAdd()
        x = MOI.add_variables(model, 3)
        cis = MOI.add_constraint.(model, x, MOI.EqualTo(0.0))
        @test cis isa Vector{MOI.ConstraintIndex{MOI.SingleVariable,
                                                 MOI.EqualTo{Float64}}}
        @test length(cis) == 3
    end

    @testset "dimension" begin
        @test MOI.dimension(MOI.ExponentialCone()) == 3
        @test MOI.dimension(MOI.DualExponentialCone()) == 3
        @test MOI.dimension(MOI.PowerCone(1/2)) == 3
        @test MOI.dimension(MOI.DualPowerCone(1/2)) == 3
        @test MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(5)) == 15
        @test MOI.dimension(MOI.PositiveSemidefiniteConeSquare(5)) == 25
        @test MOI.dimension(MOI.LogDetConeTriangle(5)) == 17
        @test MOI.dimension(MOI.LogDetConeSquare(5)) == 27
        @test MOI.dimension(MOI.RootDetConeTriangle(5)) == 16
        @test MOI.dimension(MOI.RootDetConeSquare(5)) == 26
        @test MOI.dimension(MOI.SOS1([1.0, 2.0])) == 2
        @test MOI.dimension(MOI.SOS2([1.0, 2.0])) == 2
        @test MOI.dimension(MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(1.0))) == 2
        @test MOI.dimension(MOI.Complements(5)) == 10
    end

    @testset "Dual Set" begin
        function dual_set_test(set1, set2)
            @test MOI.dual_set(set1) == set2
            @test MOI.dual_set_type(typeof(set1)) == typeof(set2)
            @test MOI.dual_set(set2) == set1
            @test MOI.dual_set_type(typeof(set2)) == typeof(set1)
        end
        function self_dual_set_test(set)
            @test MOI.dual_set(set) == set
            @test MOI.dual_set_type(typeof(set)) == typeof(set)
        end
        # Nonpositives
        nonpositives3 = MOI.Nonpositives(3)
        nonpositives4 = MOI.Nonpositives(4)
        self_dual_set_test(nonpositives3)
        @test MOI.dual_set(nonpositives3) != nonpositives4
        self_dual_set_test(nonpositives4)
        # Nonnegatives
        nonnegatives3 = MOI.Nonnegatives(3)
        nonnegatives4 = MOI.Nonnegatives(4)
        self_dual_set_test(nonnegatives3)
        @test MOI.dual_set(nonnegatives3) != nonnegatives4
        self_dual_set_test(nonnegatives4)
        # Zeros and Reals
        zeros3 = MOI.Zeros(3)
        zeros4 = MOI.Zeros(4)
        reals3 = MOI.Reals(3)
        reals4 = MOI.Reals(4)
        dual_set_test(zeros3, reals3)
        @test MOI.dual_set(reals3) != zeros4
        dual_set_test(zeros4, reals4)
        @test MOI.dual_set(zeros4) != reals3
        # Norm-1 and norm-∞ cones
        norminf2 = MOI.NormInfinityCone(2)
        norminf3 = MOI.NormInfinityCone(3)
        normone2 = MOI.NormOneCone(2)
        normone3 = MOI.NormOneCone(3)
        dual_set_test(norminf2, normone2)
        dual_set_test(norminf3, normone3)
        @test MOI.dual_set(norminf2) != normone3
        @test MOI.dual_set(normone2) != norminf3
        # SOC
        soc2 = MOI.SecondOrderCone(2)
        soc3 = MOI.SecondOrderCone(3)
        self_dual_set_test(soc2)
        @test MOI.dual_set(soc2) != soc3
        self_dual_set_test(soc3)
        # RSOC
        rsoc2 = MOI.RotatedSecondOrderCone(2)
        rsoc3 = MOI.RotatedSecondOrderCone(3)
        self_dual_set_test(rsoc2)
        @test MOI.dual_set(rsoc2) != rsoc3
        self_dual_set_test(rsoc3)
        # Norm-spectral and norm-nuclear cones
        normspec22 = MOI.NormSpectralCone(2, 2)
        normspec23 = MOI.NormSpectralCone(2, 3)
        normnuc22 = MOI.NormNuclearCone(2, 2)
        normnuc23 = MOI.NormNuclearCone(2, 3)
        dual_set_test(normspec23, normnuc23)
        dual_set_test(normspec22, normnuc22)
        @test MOI.dual_set(normspec22) != normnuc23
        @test MOI.dual_set(normnuc22) != normspec23
        # PSDtriangle
        psd2 = MOI.PositiveSemidefiniteConeTriangle(2)
        psd3 = MOI.PositiveSemidefiniteConeTriangle(3)
        self_dual_set_test(psd2)
        @test MOI.dual_set(psd2) != psd3
        self_dual_set_test(psd3)
        # Exponential
        exp = MOI.ExponentialCone()
        dual_exp = MOI.DualExponentialCone()
        dual_set_test(exp, dual_exp)
        @test MOI.dual_set(exp) != exp
        dual_set_test(dual_exp, exp)
        @test MOI.dual_set(dual_exp) != dual_exp
        # Power
        pow03 = MOI.PowerCone(0.3)
        pow04 = MOI.PowerCone(0.4)
        dual_pow03 = MOI.DualPowerCone(0.3)
        dual_set_test(pow03, dual_pow03)
        @test MOI.dual_set(pow03) != pow03
        dual_set_test(dual_pow03, pow03)
        @test MOI.dual_set(dual_pow03) != pow04
        @test MOI.dual_set(dual_pow03) != dual_pow03
        # PSDSquare error
        s = MOI.PositiveSemidefiniteConeSquare(4)
        err = ErrorException("""Dual of `PositiveSemidefiniteConeSquare` is not defined in MathOptInterface.
                                For more details see the comments in `src/Bridges/Constraint/square.jl`.""")
        @test_throws err MOI.dual_set(MOI.PositiveSemidefiniteConeSquare(4))
        @test_throws err MOI.dual_set_type(MOI.PositiveSemidefiniteConeSquare)
        # Not implemented
        s = MOI.LogDetConeTriangle(4)
        err = ErrorException("Dual of $s is not implemented.")
        @test_throws err MOI.dual_set(MOI.LogDetConeTriangle(4))
        err = ErrorException("Dual type of $(typeof(s)) is not implemented.")
        @test_throws err MOI.dual_set_type(MOI.LogDetConeTriangle)
    end
end


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

using Test
using Base.Threads

# threading constructs

# parallel loop with parallel atomic addition
function threaded_loop(a, r, x)
    @threads for i in r
        a[i] = 1 + atomic_add!(x, 1)
    end
end

function test_threaded_loop_and_atomic_add()
    x = Atomic()
    a = zeros(Int,10000)
    threaded_loop(a,1:10000,x)
    found = zeros(Bool,10000)
    was_inorder = true
    for i=1:length(a)
        was_inorder &= a[i]==i
        found[a[i]] = true
    end
    @test x[] == 10000
    # Next test checks that all loop iterations ran,
    # and were unique (via pigeon-hole principle).
    @test !(false in found)
    if was_inorder
        println(stderr, "Warning: threaded loop executed in order")
    end
end

test_threaded_loop_and_atomic_add()

# Helper for test_threaded_atomic_minmax that verifies sequential consistency.
function check_minmax_consistency(old::Array{T,1}, m::T, start::T, o::Base.Ordering) where T
    for v in old
        if v != start
            # Check that atomic op that installed v reported consistent old value.
            @test Base.lt(o, old[v-m+1], v)
        end
    end
end

function test_threaded_atomic_minmax(m::T,n::T) where T
    mid = m + (n-m)>>1
    x = Atomic{T}(mid)
    y = Atomic{T}(mid)
    oldx = Vector{T}(uninitialized, n-m+1)
    oldy = Vector{T}(uninitialized, n-m+1)
    @threads for i = m:n
        oldx[i-m+1] = atomic_min!(x, T(i))
        oldy[i-m+1] = atomic_max!(y, T(i))
    end
    @test x[] == m
    @test y[] == n
    check_minmax_consistency(oldy,m,mid,Base.Forward)
    check_minmax_consistency(oldx,m,mid,Base.Reverse)
end

# The ranges below verify that the correct signed/unsigned comparison is used.
test_threaded_atomic_minmax(Int16(-5000),Int16(5000))
test_threaded_atomic_minmax(UInt16(27000),UInt16(37000))

function threaded_add_locked(::Type{LockT}, x, n) where LockT
    critical = LockT()
    @threads for i = 1:n
        @test lock(critical) === nothing
        @test islocked(critical)
        x = x + 1
        @test unlock(critical) === nothing
    end
    @test !islocked(critical)
    nentered = 0
    nfailed = Atomic()
    @threads for i = 1:n
        if trylock(critical)
            @test islocked(critical)
            nentered += 1
            @test unlock(critical) === nothing
        else
            atomic_add!(nfailed, 1)
        end
    end
    @test 0 < nentered <= n
    @test nentered + nfailed[] == n
    @test !islocked(critical)
    return x
end

@test threaded_add_locked(SpinLock, 0, 10000) == 10000
@test threaded_add_locked(RecursiveSpinLock, 0, 10000) == 10000
@test threaded_add_locked(Mutex, 0, 10000) == 10000

# Check if the recursive lock can be locked and unlocked correctly.
let critical = RecursiveSpinLock()
    @test !islocked(critical)
    @test_throws AssertionError unlock(critical)
    @test lock(critical) === nothing
    @test islocked(critical)
    @test lock(critical) === nothing
    @test trylock(critical) == true
    @test islocked(critical)
    @test unlock(critical) === nothing
    @test islocked(critical)
    @test unlock(critical) === nothing
    @test islocked(critical)
    @test unlock(critical) === nothing
    @test !islocked(critical)
    @test_throws AssertionError unlock(critical)
    @test trylock(critical) == true
    @test islocked(critical)
    @test unlock(critical) === nothing
    @test !islocked(critical)
    @test_throws AssertionError unlock(critical)
    @test !islocked(critical)
end

# Make sure doing a GC while holding a lock doesn't cause dead lock
# PR 14190. (This is only meaningful for threading)
function threaded_gc_locked(::Type{LockT}) where LockT
    critical = LockT()
    @threads for i = 1:20
        @test lock(critical) === nothing
        @test islocked(critical)
        GC.gc(false)
        @test unlock(critical) === nothing
    end
    @test !islocked(critical)
end

threaded_gc_locked(SpinLock)
threaded_gc_locked(Threads.RecursiveSpinLock)
threaded_gc_locked(Mutex)

# Issue 14726
# Make sure that eval'ing in a different module doesn't mess up other threads
orig_curmodule14726 = @__MODULE__
main_var14726 = 1
module M14726
module_var14726 = 1
end

@threads for i in 1:100
    for j in 1:100
        @eval M14726 module_var14726 = $j
    end
end
@test @isdefined(orig_curmodule14726)
@test @isdefined(main_var14726)
@test @__MODULE__() == orig_curmodule14726

@threads for i in 1:100
    # Make sure current module is not null.
    # The @test might not be particularly meaningful currently since the
    # thread infrastructures swallows the error. (Same below)
    @test @__MODULE__() == orig_curmodule14726
end

module M14726_2
using Test
using Base.Threads
@threads for i in 1:100
    # Make sure current module is the same as the one on the thread that
    # pushes the work onto the threads.
    # The @test might not be particularly meaningful currently since the
    # thread infrastructures swallows the error. (See also above)
    @test @__MODULE__() == M14726_2
end
end

# Ensure only LLVM-supported types can be atomic
@test_throws TypeError Atomic{Bool}
@test_throws TypeError Atomic{BigInt}
@test_throws TypeError Atomic{ComplexF64}

# Test atomic memory ordering with load/store
mutable struct CommBuf
    var1::Atomic{Int}
    var2::Atomic{Int}
    correct_write::Bool
    correct_read::Bool
    CommBuf() = new(Atomic{Int}(0), Atomic{Int}(0), false, false)
end
function test_atomic_write(commbuf::CommBuf, n::Int)
    for i in 1:n
        # The atomic stores guarantee that var1 >= var2
        commbuf.var1[] = i
        commbuf.var2[] = i
    end
    commbuf.correct_write = true
end
function test_atomic_read(commbuf::CommBuf, n::Int)
    correct = true
    while true
        # load var2 before var1
        var2 = commbuf.var2[]
        var1 = commbuf.var1[]
        correct &= var1 >= var2
        var1 == n && break
        # Temporary solution before we have gc transition support in codegen.
        ccall(:jl_gc_safepoint, Cvoid, ())
    end
    commbuf.correct_read = correct
end
function test_atomic()
    commbuf = CommBuf()
    count = 1_000_000
    @threads for i in 1:2
        if i==1
            test_atomic_write(commbuf, count)
        else
            test_atomic_read(commbuf, count)
        end
    end
    @test commbuf.correct_write == true
    @test commbuf.correct_read == true
end
test_atomic()

# Test ordering with fences using Peterson's algorithm
# Example adapted from <https://en.wikipedia.org/wiki/Peterson%27s_algorithm>
mutable struct Peterson
    # State for Peterson's algorithm
    flag::Vector{Atomic{Int}}
    turn::Atomic{Int}
    # Collision detection
    critical::Vector{Atomic{Int}}
    correct::Vector{Bool}
    Peterson() =
        new([Atomic{Int}(0), Atomic{Int}(0)],
            Atomic{Int}(0),
            [Atomic{Int}(0), Atomic{Int}(0)],
            [false, false])
end
function test_fence(p::Peterson, id::Int, n::Int)
    @assert id == mod1(id,2)
    correct = true
    otherid = mod1(id+1,2)
    for i in 1:n
        p.flag[id][] = 1
        p.turn[] = otherid
        atomic_fence()
        while p.flag[otherid][] != 0 && p.turn[] == otherid
            # busy wait
            # Temporary solution before we have gc transition support in codegen.
            ccall(:jl_gc_safepoint, Cvoid, ())
        end
        # critical section
        p.critical[id][] = 1
        correct &= p.critical[otherid][] == 0
        p.critical[id][] = 0
        # end of critical section
        p.flag[id][] = 0
    end
    p.correct[id] = correct
end
function test_fence()
    commbuf = Peterson()
    count = 1_000_000
    @threads for i in 1:2
        test_fence(commbuf, i, count)
    end
    @test commbuf.correct[1] == true
    @test commbuf.correct[2] == true
end
test_fence()

# Test load / store with various types
let atomic_types = [Int8, Int16, Int32, Int64, Int128,
                    UInt8, UInt16, UInt32, UInt64, UInt128,
                    Float16, Float32, Float64]
    # Temporarily omit 128-bit types on 32bit x86
    # 128-bit atomics do not exist on AArch32.
    # And we don't support them yet on power, because they are lowered
    # to `__sync_lock_test_and_set_16`.
    if Sys.ARCH === :i686 || startswith(string(Sys.ARCH), "arm") ||
       Sys.ARCH === :powerpc64le || Sys.ARCH === :ppc64le
        filter!(T -> sizeof(T)<=8, atomic_types)
    end
    for T in atomic_types
        var = Atomic{T}()
        var[] = 42
        @test var[] === T(42)
        old = atomic_xchg!(var, T(13))
        @test old === T(42)
        @test var[] === T(13)
        old = atomic_cas!(var, T(13), T(14))   # this will succeed
        @test old === T(13)
        @test var[] === T(14)
        old = atomic_cas!(var, T(13), T(15))   # this will fail
        @test old === T(14)
        @test var[] === T(14)
    end
end

# Test atomic_cas! and atomic_xchg!
function test_atomic_cas!(var::Atomic{T}, range::StepRange{Int,Int}) where T
    for i in range
        while true
            old = atomic_cas!(var, T(i-1), T(i))
            old == T(i-1) && break
            # Temporary solution before we have gc transition support in codegen.
            ccall(:jl_gc_safepoint, Cvoid, ())
        end
    end
end
for T in (Int32, Int64, Float32, Float64)
    var = Atomic{T}()
    nloops = 1000
    di = nthreads()
    @threads for i in 1:di
        test_atomic_cas!(var, i:di:nloops)
    end
    @test var[] === T(nloops)
end

function test_atomic_xchg!(var::Atomic{T}, i::Int, accum::Atomic{Int}) where T
    old = atomic_xchg!(var, T(i))
    atomic_add!(accum, Int(old))
end
for T in (Int32, Int64, Float32, Float64)
    accum = Atomic{Int}()
    var = Atomic{T}()
    nloops = 1000
    @threads for i in 1:nloops
        test_atomic_xchg!(var, i, accum)
    end
    @test accum[] + Int(var[]) === sum(0:nloops)
end

function test_atomic_float(varadd::Atomic{T}, varmax::Atomic{T}, varmin::Atomic{T}, i::Int) where T
    atomic_add!(varadd, T(i))
    atomic_max!(varmax, T(i))
    atomic_min!(varmin, T(i))
end
for T in (Int32, Int64, Float32, Float64)
    varadd = Atomic{T}()
    varmax = Atomic{T}()
    varmin = Atomic{T}()
    nloops = 1000
    @threads for i in 1:nloops
        test_atomic_float(varadd, varmax, varmin, i)
    end
    @test varadd[] === T(sum(1:nloops))
    @test varmax[] === T(maximum(1:nloops))
    @test varmin[] === T(0)
end

using Dates
for period in (0.06, Dates.Millisecond(60))
    let async = Base.AsyncCondition(), t
        c = Condition()
        task = schedule(Task(function()
            notify(c)
            wait(c)
            t = Timer(period)
            wait(t)
            ccall(:uv_async_send, Cvoid, (Ptr{Cvoid},), async)
            ccall(:uv_async_send, Cvoid, (Ptr{Cvoid},), async)
            wait(c)
            sleep(period)
            ccall(:uv_async_send, Cvoid, (Ptr{Cvoid},), async)
            ccall(:uv_async_send, Cvoid, (Ptr{Cvoid},), async)
        end))
        wait(c)
        notify(c)
        delay1 = @elapsed wait(async)
        notify(c)
        delay2 = @elapsed wait(async)
        @test istaskdone(task)
        @test delay1 > 0.05
        @test delay2 > 0.05
        @test isopen(async)
        @test !isopen(t)
        close(t)
        close(async)
        @test_throws EOFError wait(async)
        @test !isopen(async)
        @test_throws EOFError wait(t)
        @test_throws EOFError wait(async)
    end
end

complex_cfunction = function(a)
    s = zero(eltype(a))
    @inbounds @simd for i in a
        s += muladd(a[i], a[i], -2)
    end
    return s
end
function test_thread_cfunction()
    @threads for i in 1:1000
        # Make sure this is not inferrable
        # and a runtime call to `jl_function_ptr` will be created
        ccall(:jl_function_ptr, Ptr{Cvoid}, (Any, Any, Any),
              complex_cfunction, Float64, Tuple{Ref{Vector{Float64}}})
    end
end
test_thread_cfunction()

# Compare the two ways of checking if threading is enabled.
# `jl_tls_states` should only be defined on non-threading build.
if ccall(:jl_threading_enabled, Cint, ()) == 0
    @test nthreads() == 1
    cglobal(:jl_tls_states) != C_NULL
else
    @test_throws ErrorException cglobal(:jl_tls_states)
end

function test_thread_range()
    a = zeros(Int, nthreads())
    @threads for i in 1:threadid()
        a[i] = 1
    end
    for i in 1:threadid()
        @test a[i] == 1
    end
    for i in (threadid() + 1):nthreads()
        @test a[i] == 0
    end
end
test_thread_range()

# Thread safety of `jl_load_and_lookup`.
function test_load_and_lookup_18020(n)
    @threads for i in 1:n
        try
            ccall(:jl_load_and_lookup,
                  Ptr{Cvoid}, (Cstring, Cstring, Ref{Ptr{Cvoid}}),
                  "$i", :f, C_NULL)
        end
    end
end
test_load_and_lookup_18020(10000)

# Nested threaded loops
# This may not be efficient/fully supported but should work without crashing.....
function test_nested_loops()
    a = zeros(Int, 100, 100)
    @threads for i in 1:100
        @threads for j in 1:100
            a[j, i] = i + j
        end
    end
    for i in 1:100
        for j in 1:100
            @test a[j, i] == i + j
        end
    end
end
test_nested_loops()

@testset "libatomic" begin
    prog = """
    using Base.Threads
    using InteractiveUtils: code_native
    function unaligned_setindex!(x::Atomic{UInt128}, v::UInt128)
        Base.llvmcall(\"\"\"
            %ptr = inttoptr i$(Sys.WORD_SIZE) %0 to i128*
            store atomic i128 %1, i128* %ptr release, align 8
            ret void
        \"\"\", Cvoid, Tuple{Ptr{UInt128}, UInt128}, unsafe_convert(Ptr{UInt128}, x), v)
    end
    code_native(stdout, unaligned_setindex!, Tuple{Atomic{UInt128}, UInt128})
    """

    mktempdir() do dir
        file = joinpath(dir, "test23901.jl")
        write(file, prog)
        run(pipeline(ignorestatus(`$(Base.julia_cmd()) --startup-file=no $file`),
                     stdout=joinpath(dir, "out.txt"),
                     stderr=joinpath(dir, "err.txt"),
                     append=false))
        out = readchomp(joinpath(dir, "out.txt"))
        err = readchomp(joinpath(dir, "err.txt"))
        @test contains(out, "libat_store") || contains(out, "atomic_store")
        @test !contains(err, "__atomic_store")
    end
end

function test_thread_too_few_iters()
    x = Atomic()
    a = zeros(Int, nthreads()+2)
    threaded_loop(a, 1:nthreads()-1, x)
    found = zeros(Bool, nthreads()+2)
    for i=1:nthreads()-1
        found[a[i]] = true
    end
    @test x[] == nthreads()-1
    # Next test checks that all loop iterations ran,
    # and were unique (via pigeon-hole principle).
    @test !(false in found[1:nthreads()-1])
    @test !(true in found[nthreads():end])
end
test_thread_too_few_iters()

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
module v2_tiny_voc

import ...flipdict, ...Settings

#20 classes on Voc dataset
# 2 dictionaries to access number<->class by O(1)
const namesdic = Dict("aeroplane"=>1,"bicycle"=>2,"bird"=>3, "boat"=>4,
            "bottle"=>5,"bus"=>6,"car"=>7,"cat"=>8,"chair"=>9,
            "cow"=>10,"diningtable"=>11,"dog"=>12,"horse"=>13,"motorbike"=>14,
            "person"=>15,"pottedplant"=>16,"sheep"=>17,"sofa"=>18,"train"=>19,"tvmonitor"=>20)

#Yolo V2-tiny pre-trained boxes width and height
const anchors = [(1.08,1.19),  (3.42,4.41),  (6.63,11.38),  (9.42,5.11),  (16.62,10.52)]

function load(;minibatch_size=1)
    return Settings(
        dataset_description = "VOC2007",
        source = "https://pjreddie.com/media/files/yolov2-tiny-voc.weights",
        weights_filepath = joinpath(@__DIR__,"v2_tiny_voc.weights"),
        image_shape = (416, 416),
        image_channels = 3,
        namesdic = namesdic,
        numsdic = flipdict(namesdic), #Flip the key value pair to help lookup by nums
        anchors = anchors,
        num_classes = length(namesdic),
        minibatch_size = minibatch_size,
        grid_x = 13,
        grid_y = 13,
        cell_bboxes = 5
    )
end

end #module


text:
@testset "filecmp" begin
    gotd()

    p1 = joinpath(td, "hello.md")
    p2 = joinpath(td, "bye.md")
    p3 = joinpath(td, "cp1.md")

    write(p1, "foo")
    write(p2, "foo")

    cp(p1, p3)

    @test F.filecmp(p1, p1)
    @test F.filecmp(p1, p2)
    @test F.filecmp(p1, p3)

    write(p2, "baz")
    @test !F.filecmp(p1, p2)

    @test !F.filecmp(p1, "foo/bar")
end


text:
# EXEMPLE DE LIGNE DE COMMANDE TESTÉE :
# ./bin/seqata.jl val -i data/ampl/01.ampl -p 3,4,5,6,7,8,9,1,10,2 -x clp
#

# ===========
# Préparation des arguments
args = Args.parse_commandline()
const TEST_LEVEL = 1
Args.set(:level, TEST_LEVEL) # sera passé temporairement à 0
Args.set(:outdir, "/tmp") # equivalent à l'option -d/tmp

# ===========
# Construction de l'instance mini
inst = instance_build_mini10()
@test inst.nb_planes == 10

# Création de la solution vide
sol = Solution(inst, algo=:earliest) # sera refait pour chaque test
# @show join(sol.planes, ",")

lg1("1. La solution initiale de alp10 doit être faisable... ")
sol = Solution(inst, algo=:earliest)
nbviols, violtxt = get_viol_description(sol)
nblines_in_violtxt = length(findall(r".+", violtxt))
@test nbviols == 0
@test nbviols == nblines_in_violtxt
@test length(violtxt) == 0
ln1("fait.")

lg1("2. La solution viole les bornes lb et ub... ")
sol = Solution(inst, algo=:earliest)
# on impose l'avion x[3] atterrit trop tot et x[5] trop tard
sol.x[3] = 0
sol.x[5] = 1000
nbviols, violtxt = get_viol_description(sol)
expected = "
p3 atterrit trop tôt  x=0 < lb=90
p5 atterrit trop tard ub=556 < x=1000
p2->p3 : écart insuffisant x[p2]=196->x[p3]=0 => sep(p2,p3)=-196 au lieu de 15
p1->p3 : écart insuffisant x[p1]=130->x[p3]=0 => sep(p1,p3)=-130 au lieu de 15
p5->p6 : écart insuffisant x[p5]=1000->x[p6]=235 => sep(p5,p6)=-765 au lieu de 8
p5->p7 : écart insuffisant x[p5]=1000->x[p7]=243 => sep(p5,p7)=-757 au lieu de 8
p5->p8 : écart insuffisant x[p5]=1000->x[p8]=251 => sep(p5,p8)=-749 au lieu de 8
p5->p9 : écart insuffisant x[p5]=1000->x[p9]=259 => sep(p5,p9)=-741 au lieu de 8
p5->p10 : écart insuffisant x[p5]=1000->x[p10]=267 => sep(p5,p10)=-733 au lieu de 8
"
if lg3()
    println("\nexpected=\n", expected)
    println("\nvioltxt=\n", violtxt)
end
nblines_in_violtxt = length(findall(r".+", violtxt))
@test nbviols == 9
@test nbviols == nblines_in_violtxt
ln1("fait.")

lg1("3. La solution viole les temps de séparation... ")

sol = Solution(inst, algo=:earliest)
# on impose l'avion x[3] à la valeur de x[5] => p4 et p5 deviennet faux
sol.x[3] = sol.x[5]
nbviols, violtxt = get_viol_description(sol)
nblines_in_violtxt = length(findall(r".+", violtxt))
@test nbviols == 2
@test nbviols == nblines_in_violtxt
if lg3()
    println("\nvioltxt=\n", violtxt)
end

ln1(" fait.\n")


text:
############################################################################################
buffer_nparameter = 100
buffer_newparam = 10
buffer_ValType = [Float32, Float64]
buffer_IType = [Int32, Int64]

############################################################################################
# ParameterBuffer
@testset "Buffer - ParameterBuffer: " begin
    for vtype in buffer_ValType
        for itype in buffer_IType
            param1 = ParameterBuffer(vtype.(1:buffer_nparameter), buffer_nparameter, itype)
            param2 = ParameterBuffer(vtype(1), buffer_nparameter, itype)
            # Check types and size of buffer
            @test eltype(param1.val) == eltype(param2.val) == vtype
            @test eltype(param1.index) == eltype(param2.index) == itype
            @test length(param1.val) == length(param1.index) ==
                length(param2.val) == length(param2.index) == buffer_nparameter
            # Resize and check again
            buffer_nparameter_new = buffer_nparameter + buffer_newparam
            _param1 = deepcopy(param1)
            _param2 = deepcopy(param2)
            resize!(_param1, buffer_nparameter_new)
            resize!(_param2, buffer_nparameter_new)
            @test length(_param1.val) == length(_param1.index) ==
                length(_param2.val) == length(_param2.index) == buffer_nparameter_new
            @test length(param1.val) + buffer_newparam == length(_param1.val)
            # Shuffle and check if correct
            θnew = collect(100.0:-1:1.0)
            __param1 = deepcopy(param1)
            shuffle!(__param1, θnew)
            @test __param1.val[begin] == θnew[begin]
        end
    end
end

############################################################################################
# ModelParameterBuffer
mutable struct Mod{T} <: AbstractModelWrapper
    val::T
end
buffer_vals = [
    (a = 1., b = 2),
    (c = zeros(10), d = zeros(2,3))
]
struct Res <: AbstractResult end
mutable struct Alg{R} <: AbstractAlgorithm
    res::R
end

function get_result(alg::Alg)
    return alg.res
end
function result!(alg::Alg, result)
    alg.res = result
end

@testset "Buffer - ModelParameterBuffer: " begin
    for val in buffer_vals
        for itype in buffer_IType
            param = ModelParameterBuffer(Mod(val), Res(), buffer_nparameter, itype)
            # Check types and size of buffer
            @test eltype(param.val) == typeof(val)
            @test length(param.val) == length(param.index) ==
                length(param.result) == length(param.weight ) == buffer_nparameter
            # Resize and check again
            buffer_nparameter_new = buffer_nparameter + buffer_newparam
            _param = deepcopy(param)
            resize!(_param, buffer_nparameter_new)
            @test length(_param.val) == length(_param.index) ==
                length(_param.result) == length(_param.weight ) == buffer_nparameter_new
            @test length(param.val) + buffer_newparam == length(_param.index)
            # Shuffle and check if correct
            __param = deepcopy(param)
            algs = [Alg(Res()) for _ in 1:buffer_nparameter]
            mods = [Mod(val) for _ in 1:buffer_nparameter]
            weights = collect(1.0:1:buffer_nparameter)
            shuffle!(__param, algs, mods, weights)
            @test __param.weight[end] == buffer_nparameter
        end
    end
end


text:
using Pkg
Pkg.activate("../Scratchy")
using Scratchy, Scratchy.Layers, Scratchy.Losses
using Scratchy.Activations, Scratchy.Models

X = [0 0; 1 1; 0 1; 1 0]
y = [0, 0, 1, 1]

layer1 = Dense("layer1", (2, 3), Sigmoid)
layer2 = Dense("layer2", (3, 1), Sigmoid)
net = Network("FFNN", [layer1, layer2], BinaryCrossEntropy, 0.2)
Scratchy.Models.display(net)

fit!(net, X, y, 1000)

logits = Scratchy.Models.forward!(net, X)
preds = Scratchy.Utils.decision(logits)
acc = Scratchy.Utils.accuracy(y, preds)

println("Logits: $(logits)")
println("Predictions: $(preds)")
println("Gold Labels: $(y)")
println("Accuracy: $(acc)")

max_stars_repo_name: jvasilakes/scratchy

max_issues_repo_name: jvasilakes/scratchy

max_forks_repo_name: jvasilakes/scratchy


text:
using SegyIO; const SIO = SegyIO
export segy_open, read_ensemble

const SEGY_DATA_TYPES = Dict(1=>SIO.IBMFloat32, 5=>Float32)
const TRACE_HDR_BYTE2SAMPLE = SIO.th_byte2sample()

#------------------------------------------------------------------------------
# Methods added by Marc.
"""
segy_open(file::String; buffer::Bool = true, warn_user::Bool = true)

Open a SEGY file and return the stream and the file header.
"""
function segy_open(file::String; buffer::Bool = true, warn_user::Bool = true)

    s = fh = nothing

    if buffer
        s = IOBuffer(read(open(file)))
    else
        s = open(file)
    end

    # Read File Header
    fh = read_fileheader(s)

    # Check fixed length trace flag
    (fh.bfh.FixedLengthTraceFlag!=1 & warn_user) &&
        @warn "Fixed length trace flag set in stream: $s"

    # Check datatype of file
    if fh.bfh.DataSampleFormat in keys(SEGY_DATA_TYPES)
        datatype = SEGY_DATA_TYPES[fh.bfh.DataSampleFormat]
    else
        @error "Data type not supported ($(fh.bfh.DataSampleFormat))"
    end

    @info("[segy_open] bfh: ns=$(fh.bfh.ns), dt=$(fh.bfh.dt), "*
            "fmt=$(fh.bfh.DataSampleFormat) ($(datatype))")

    return s, fh
end

"""
read_ensemble(s)

Read a SEGY ensemble from stream 's'. Ensemble are delimited by a change in
the values of the headers in `delimkeys`.
"""
function read_ensemble(s::IO, fh::SIO.FileHeader,
        keys::Array{String}, delimkeys::Array{String},
        maxTraces::Int, startbyte::Int, prevTrace::Tuple;
        hdrs::Union{Nothing,Array{Number,1}}=nothing,
        data::Union{Nothing,Array{Number,2}}=nothing,
        warn_user::Bool=true)

    datatype = SEGY_DATA_TYPES[fh.bfh.DataSampleFormat]

    # Make sure all delimkeys are in keys.
    # TODO: Needed for every read?
    for k in delimkeys
        in(k, keys) ? nothing : push!(keys, k)
    end

    # Move to start of trace, skipping EBCDIC header if first trace.
    startbyte = max(startbyte, 3600)
    seek(s, startbyte)

    trc_length = get_trace_length(Int(fh.bfh.ns), nbytes=sizeof(datatype))
    @debug("[read_ens] startbyte=$(startbyte), trc_length=$trc_length")

    # Allocate arrays if required.
    if hdrs == nothing
        @debug("[read_ens] allocating hdrs, ntrc=$(maxTraces).")
        hdrs = Array{BinaryTraceHeader,1}(undef, maxTraces)
    end
    if data == nothing
        @debug("[read_ens] allocating data, ns=$(fh.bfh.ns),ntrc=$(maxTraces).")
        data = Array{datatype,2}(undef, fh.bfh.ns, maxTraces)
    end

    # Read first trace in ensemble.
    SIO.read_trace!(s, fh.bfh, datatype, hdrs, data, 1, keys,
                    TRACE_HDR_BYTE2SAMPLE)

    # Initialize array of values of delimiter headers.
    prevVals = Int[]
    for k in delimkeys
        push!(prevVals, getfield(hdrs[1], Symbol(k)))
    end

    prevTrace = (prevTrace[1] == nothing) ? prevTrace = ([], []) : prevTrace

    # Read traces until delimiter headers change, or max number of traces is
    # reached.
    # TODO: Save last trace read to use as first trace in next ensemble?
    itrc = 2; ntrcEns = itrc

    while itrc <= maxTraces
        try
            SIO.read_trace!(s, fh.bfh, datatype, hdrs, data, itrc,
                            keys, TRACE_HDR_BYTE2SAMPLE)
        catch err
            if typeof(err) <: EOFError
                @debug "[read_ens] Read EOF on input stream."
                ntrcEns = itrc - 1
                break
            else
                @error "[read_ensemble] Error reading trace from input stream."
                throw(err)
            end
        end

        if is_header_changed(hdrs[itrc], delimkeys, prevVals)
            @debug("[read_ens] hdr change at trace $itrc - exiting loop")
            ntrcEns = itrc - 1
#            prevTrace = (hdrs[itrc], data[itrc])
            break
        elseif itrc == maxTraces
            @error("[read_ens] Max traces: ($maxTraces) read for ensemble "*
                   "before change in delim keys: $(delimkeys)")
            throw(BoundsError)
        end
        @debug("[read_ens] Read input trace $itrc, "*
               "cdp=$cdp, offset=$offset, byte_pos=$(position(s))")
        itrc += 1
        ntrcEns = itrc
    end

    # Set position of the next ensemble by subtracting the length of the
    # last trace read. Not needed if ntrcEns == maxTraces.
    byte_pos = position(s)
    if ntrcEns < maxTraces && !eof(s)
      byte_pos -= trc_length
      @debug("[read_ens] backing up one trace to byte_pos=$(byte_pos).")
    end
    @debug("[read_ensemble] Ensemble completed, num traces $(ntrcEns), "*
           "byte_pos=$byte_pos")

    return SIO.SeisBlock(fh, hdrs, data), byte_pos, ntrcEns, prevTrace
end

"""
Return true if the delimiter keys in hdr are different from lastVals.
"""
function is_header_changed(hdr::SIO.BinaryTraceHeader,
                           delimkeys::Array{String,1},
                           lastVals::Array{Int,1})
    changed = false
    for (i, k) in enumerate(delimkeys)
        if (lastVals[i] != getfield(hdr, Symbol(k))) return true end
    end
    return changed
end

"""
Functions to return trace length in bytes and number of traces
from start/end bytes.
"""
function get_trace_length(ns::Int; nbytes::Int=4, hdrlen::Int=240)
    return hdrlen + (ns*nbytes)
end
function get_num_traces(startbyte::Int, endbyte::Int, ns::Int;
                        hdrlen::Int=240, nbytes::Int=4)
    return (endbyte - startbyte) / get_trace_length(ns, nbytes)
end


text:
using Scallop
using Base.Test

# write your own tests here
include(joinpath(Pkg.dir("Scallop"), "src", "main.jl"))


text:
module CommonSolve

solve(args...; kwargs...) = solve!(init(args...; kwargs...))
function solve! end
function init end

end # module


text:
export Cuboid , Sphere

abstract type GeometryPrimitive end

struct Cuboid <: GeometryPrimitive
    bounds
end
function (c::Cuboid)(x)
    prod(c.bounds[1] .<= x .<= c.bounds[2])
end

struct Sphere <: GeometryPrimitive
    radius
    center
end

function (s::Sphere)(x)
    norm(x .- s.center) <= s.radius
end


text:
@testset "prior" begin

end


text:
export RejectionSampler
export InverseTransformSampler
export BoxMullerTransformSampler

@doc raw"""
    RejectionSampler(f, g, proposal_sampler, M)

A rejection sampler.

Rejection sampling involves sampling from a target distribution ``f`` using a
proposal distribution ``g`` which we are able to directly sample from. In order
for the resulting samples to be distributed according to ``f`` we require that
``\frac{f(x)}{g(x)}`` is bounded from above by some constant ``M`` (`M`) for all ``x`` in
the support of ``f``.

Under these conditions, the rejection sampling algorithm given below will
generate independent samples from ``f``:

* Sample ``X \sim g``
* Accept ``X`` with probability ``\frac{f(X)}{Mg(X)}``

# Arguments
* `f::Function`: density function for the target distribution.
* `g::Function`: density function for the proposal distribution.
* `proposal_sampler::Function`: a function that generates samples from the
    proposal distribution.
* `M::Real`: a scaling constant that bounds the ratio of the target and
    proposal density.

# Notes
* No checks are made to ensure that the scale constant is valid or that the
  proposal sampler and density match. If these conditions are not met, the
  resulting samples will not be distributed according to the target density.
* The expected acceptance rate for the sampler is the reciprocal of `M`.

# Examples
```julia
# Generating 10 standard normal samples using a Cauchy proposal
f(x) = exp(-x^2 / 2) / sqrt(2π)
g(x) = 1 / (π * (1 + x^2))
proposal_sampler = InverseTransformSampler(u -> tan(π * (u - 0.5)))
M = sqrt(2π / ℯ)
s = RejectionSampler(f, g, proposal_sampler, M)
samples = sample(s, 10)
```
"""
struct RejectionSampler <: Sampler
    f::Function
    g::Function
    proposal_sampler::Sampler
    M::Real
end

@doc raw"""
    sample(s::RejectionSampler, N)

Generate samples using a rejection sampler.

# Arguments
* `s::RejectionSampler`: a rejection sampler.
* `N::Integer`: the number of samples to generate.

# Returns
* `samples::Array{Float64}`: a 2-D array in which each column is a sample.
"""
function sample(s::RejectionSampler, N::Integer)
    # TODO: make this work for arbritary dimensions
    samples = Array{Float64}(undef, 1, N)
    i = 1
    while i <= N
        proposal = sample(s.proposal_sampler, 1)[1]

        threshold = s.f(proposal) / (s.M * s.g(proposal))
        if rand() < threshold
            samples[:, i] .= proposal
            i += 1
        end
    end

    return samples
end

@doc raw"""
    InverseTransformSampler(F_inv)

An inverse transform sampler.

Inverse transform sampling is based on the result that when
``U \sim \text{Unif}(0, 1)``, ``X = F^{-1}(U)`` will be distributed
according to the CDF ``F``. For non-decreasing CDFs we can replace ``F^{-1}``
with ``F``'s generalised inverse ``F^{-}(u) = \inf{x : F(x) > u}`` and use
the same method.

# Arguments
* `F_inv::Function`: the inverse CDF of the target distrubition.

# Notes
* Because of the symmetry of a uniform random sample, one can replace ``1-u``
  in the inverse CDF with ``u``.

# Examples
```julia
# Generating 10 standard exponential random variables
F_inv(u) = -log(u)  # using symmetry of U (see notes)
s = InverseTransformSampler(F_inv)
samples = sample(s, 10)
```
"""
struct InverseTransformSampler <: Sampler
    F_inv::Function
end

@doc raw"""
    sample(s::InverseTransformSampler, N)

Generate samples using an inverse transform sampler.

# Arguments
* `s::InverseTransformSampler`: an inverse transform sampler.
* `N::Integer`: the number of samples to generate.

# Returns
* `samples::Array{Float64}`: a 2-D array in which each column is a sample.
"""
function sample(s::InverseTransformSampler, N::Integer)
    return s.F_inv.(rand(1, N))
end

@doc raw"""
    BoxMullerTransformSampler()

A Box-Muller transform sampler.

The Box-Muller transform is a technique for generating independent, standard,
normally distributed samples. It involves sampling from a polar coordinate
system in which the angle of the sample is uniformly distrubition and the
radius is distributed as ``\text{exp}(\tfrac{1}{2})`` (equivalent to
``\chi^2_1``). It can then be shown that the cartesian coordinates of these
samples have independent standard normal distributions.

The full algorithm is given here:

* Sample ``U_1, U_2 \stackrel{\text{i.i.d.}}{\sim} \text{Unif}(0, 1)``
* Compute ``Z_1 = \sqrt{-2 \log U_1}\cos(2\pi U_2)``, ``Z_2 = \sqrt{-2 \log U_1}\sin(2\pi U_2)``

# Arguments
_None_

# Examples
```julia
# Generating 10 standard normal samples
s = BoxMullerTransformSampler()
samples = sample(s, 10)
```
"""
struct BoxMullerTransformSampler <: Sampler
end

@doc raw"""
    sample(s::BoxMullerTransformSampler, N)

Generate samples using a Box-Muller transform sampler.

# Arguments
* `s::BoxMullerTransformSampler`: a Box-Muller transform sampler.
* `N::Integer`: the number of samples to generate.

# Returns
* `samples::Array{Float64}`: a 2-D array in which each column is a sample.
"""
function sample(s::BoxMullerTransformSampler, N::Integer)
    M = N + N % 2
    U1 = rand(1, M)
    U2 = rand(1, M)
    R = sqrt.(-2 * log.(U1))
    θ = 2π * U2
    Z1 = R .* cos.(θ)
    Z2 = R .* sin.(θ)
    return hcat(Z1, Z2)[:, 1:N]
end


text:
solve_liner = Symbolics.solve_for

function solve(problem::LinearProblem)
    eqs = problem.eqs
    vars = problem.vars
    sol = solve_liner(eqs, vars)
    return Dict([(vars[i], sol[i]) for i in 1:length(vars)])
end

function one_cal(array, dict, f, ja)
    return array .- inv(Symbolics.value.(substitute.(ja, (dict,)))) * Symbolics.value.(substitute.(f, (dict,)))
end

function get_dict(array, dict, vars)
    for i in 1:length(vars)
        dict[vars[i]] = array[i]
    end
    return dict
end

function _solve(dict, f, ja, maxiter, abstol, vars)
    last_v = Inf
    new_v = collect(values(dict))
    dict = dict
    count = 0

    while count < maxiter && maximum(abs.(new_v .- last_v)) > abstol
        last_v = new_v
        new_v = one_cal(new_v, dict, f, ja)
        dict = get_dict(new_v, dict, vars)
        count += 1
    end

    if count >= maxiter
        error("Error: iterate times exceed maxiter")
    end

    return dict
end

function solve(problem::NonlinearProblem)
    eqs = problem.eqs
    dict_u = problem.vars
    vars = collect(keys(problem.vars))
    maxiters = problem.maxiters
    abstol = problem.abstol
    jac = Symbolics.jacobian(eqs, vars)
    sol = _solve(dict_u, eqs, jac, maxiters, abstol,vars)
    return sol
end


text:
module SeqUMAP

using Distances
import MultivariateStats: fit, PCA, transform
import UMAP: umap

include("encoding.jl")
include("embedding.jl")
include("projection.jl")

#encoding.jl...
export AA_DICT,
NT_DICT,
IUPACbool,
resolve_base,
resolve_seq,
string2encoding,

#embedding.jl...
collect_kmers,
kmer_count!,
get_kmer_index,
kmer_embed,

#projection.jl...
CorrectedKmer,
sequmap

end

max_stars_repo_name: alecpnkw/SeqUMAP

max_issues_repo_name: alecpnkw/SeqUMAP

max_forks_repo_name: alecpnkw/SeqUMAP