text:
# This example requires Cloudy to be installed.
#using Pkg; Pkg.add(PackageSpec(name="Cloudy", version="0.1.0"))
using Cloudy
const PDistributions = Cloudy.ParticleDistributions

# Import modules
using Distributions  # probability distributions and associated functions
using StatsBase
using LinearAlgebra
using StatsPlots
using Plots
using JLD2
using Random

# Import Calibrate-Emulate-Sample modules
using EnsembleKalmanProcesses.EnsembleKalmanProcessModule
using EnsembleKalmanProcesses.Observations
using EnsembleKalmanProcesses.ParameterDistributionStorage
using EnsembleKalmanProcesses.DataStorage

# Import the module that runs Cloudy
include("GModel.jl")
using .GModel

################################################################################
#                                                                              #
#                      Cloudy Calibrate-Emulate-Sample Example                 #
#                                                                              #
#                                                                              #
#     This example uses Cloudy, a microphysics model that simulates the        #
#     coalescence of cloud droplets into bigger drops, to demonstrate how      #
#     the full Calibrate-Emulate-Sample pipeline can be used for Bayesian      #
#     learning and uncertainty quantification of parameters, given some        #
#     observations.                                                            #
#                                                                              #
#     Specifically, this examples shows how to learn parameters of the         #
#     initial cloud droplet mass distribution, given observations of some      #
#     moments of that mass distribution at a later time, after some of the     #
#     droplets have collided and become bigger drops.                          #
#                                                                              #
#     In this example, Cloudy is used in a "perfect model" (aka "known         #
#     truth") setting, which means that the "observations" are generated by    #
#     Cloudy itself, by running it with the true parameter values. In more     #
#     realistic applications, the observations will come from some external    #
#     measurement system.                                                      #
#                                                                              #
#     The purpose is to show how to do parameter learning using                #
#     Calibrate-Emulate-Sample in a simple (and highly artificial) setting.    #
#                                                                              #
#     For more information on Cloudy, see                                      #
#              https://github.com/CliMA/Cloudy.jl.git                          #
#                                                                              #
################################################################################

rng_seed = 41
Random.seed!(rng_seed)

homedir = pwd()
println(homedir)
figure_save_directory = homedir*"/output/"
data_save_directory = homedir*"/output/"
if ~isdir(figure_save_directory)
    mkdir(figure_save_directory)
end
if ~isdir(data_save_directory)
    mkdir(data_save_directory)
end

###
###  Define the (true) parameters and their priors
###

# Define the parameters that we want to learn
# We assume that the true particle mass distribution is a Gamma distribution
# with parameters N0_true, θ_true, k_true
param_names = ["N0", "θ", "k"]
n_param = length(param_names)

N0_true = 300.0  # number of particles (scaling factor for Gamma distribution)
θ_true = 1.5597  # scale parameter of Gamma distribution
k_true = 0.0817  # shape parameter of Gamma distribution
params_true = [N0_true, θ_true, k_true]
# Note that dist_true is a Cloudy distribution, not a Distributions.jl
# distribution
dist_true = PDistributions.Gamma(N0_true, θ_true, k_true)

###
###  Define priors for the parameters we want to learn
###

# Define constraints
lbound_N0 = 0.4 * N0_true
lbound_θ = 1.0e-1
lbound_k = 1.0e-4
c1 = bounded_below(lbound_N0)
c2 = bounded_below(lbound_θ)
c3 = bounded_below(lbound_k)
constraints = [[c1], [c2], [c3]]

# We choose to use normal distributions to represent the prior distributions of
# the parameters in the transformed (unconstrained) space. i.e log coordinates
d1 = Parameterized(Normal(4.5, 1.0)) #truth is 5.19
d2 = Parameterized(Normal(0.0, 2.0)) #truth is 0.378
d3 = Parameterized(Normal(-1.0, 1.0))#truth is -2.51
distributions = [d1, d2, d3]

param_names = ["N0", "θ", "k"]

priors = ParameterDistribution(distributions, constraints, param_names)

###
###  Define the data from which we want to learn the parameters
###

data_names = ["M0", "M1", "M2"]
moments = [0.0, 1.0, 2.0]
n_moments = length(moments)

###
###  Model settings
###

# Collision-coalescence kernel to be used in Cloudy
coalescence_coeff = 1/3.14/4/100
kernel_func = x -> coalescence_coeff
kernel = Cloudy.KernelTensors.CoalescenceTensor(kernel_func, 0, 100.0)

# Time period over which to run Cloudy
tspan = (0., 1.0)

###
###  Generate (artificial) truth samples
###  Note: The observables y are related to the parameters θ by:
###        y = G(x1, x2) + η
###

g_settings_true = GModel.GSettings(kernel, dist_true, moments, tspan)
gt = GModel.run_G(params_true, g_settings_true, PDistributions.update_params,
                  PDistributions.moment, Cloudy.Sources.get_int_coalescence)
n_samples = 100
yt = zeros(length(gt),n_samples)
# In a perfect model setting, the "observational noise" represent the internal
# model variability. Since Cloudy is a purely deterministic model, there is no
# straightforward way of coming up with a covariance structure for this internal
# model variability. We decide to use a diagonal covariance, with entries
# (variances) largely proportional to their corresponding data values, gt.
Γy = convert(Array, Diagonal([100.0, 5.0, 30.0]))
μ = zeros(length(gt))

# Add noise
for i in 1:n_samples
    yt[:, i] = gt .+ rand(MvNormal(μ, Γy))
end

truth = Observations.Obs(yt, Γy, data_names)
truth_sample = truth.mean
###
###  Calibrate: Ensemble Kalman Inversion
###

N_ens = 50 # number of ensemble members
N_iter = 8 # number of EKI iterations
# initial parameters: N_params x N_ens
initial_params = EnsembleKalmanProcessModule.construct_initial_ensemble(priors, N_ens; rng_seed=6)
ekiobj = EnsembleKalmanProcessModule.EnsembleKalmanProcess(initial_params, truth_sample, truth.obs_noise_cov,
                   Inversion(), Δt=0.1)

# Initialize a ParticleDistribution with dummy parameters. The parameters
# will then be set in run_G_ensemble
dummy = 1.0
dist_type = PDistributions.Gamma(dummy, dummy, dummy)
g_settings = GModel.GSettings(kernel, dist_type, moments, tspan)

# EKI iterations
for i in 1:N_iter

    params_i = mapslices(x -> transform_unconstrained_to_constrained(priors, x),
                         get_u_final(ekiobj); dims=1)
    g_ens = GModel.run_G_ensemble(params_i, g_settings,
                                  PDistributions.update_params,
                                  PDistributions.moment,
                                  Cloudy.Sources.get_int_coalescence)
    EnsembleKalmanProcessModule.update_ensemble!(ekiobj, g_ens)
end

# EKI results: Has the ensemble collapsed toward the truth?
transformed_params_true = transform_constrained_to_unconstrained(priors,
                                                                 params_true)
println("True parameters (transformed): ")
println(transformed_params_true)

println("\nEKI results:")
println(mean(get_u_final(ekiobj), dims=2))

u_stored= get_u(ekiobj,return_array=false)
g_stored= get_g(ekiobj,return_array=false)
@save data_save_directory*"parameter_storage.jld2" u_stored
@save data_save_directory*"data_storage.jld2" g_stored

#plots
gr(size=(1800,600))

u_init = get_u_prior(ekiobj)
for i in 1:N_iter
    u_i = get_u(ekiobj,i)
    p1 = plot(u_i[1,:], u_i[2,:], seriestype=:scatter, xlims = extrema(u_init[1,:]), ylims = extrema(u_init[2,:]))
    plot!(p1,[transformed_params_true[1]], xaxis="u1", yaxis="u2", seriestype="vline",
        linestyle=:dash, linecolor=:red, label = false,
        title = "EKI iteration = " * string(i)
        )
    plot!(p1,[transformed_params_true[2]], seriestype="hline", linestyle=:dash, linecolor=:red, label = "optimum")

    p2 = plot(u_i[2,:], u_i[3,:], seriestype=:scatter, xlims = extrema(u_init[2,:]), ylims = extrema(u_init[3,:]))
    plot!(p2,[transformed_params_true[2]], xaxis="u1", yaxis="u2", seriestype="vline",
        linestyle=:dash, linecolor=:red, label = false,
        title = "EKI iteration = " * string(i)
        )
    plot!(p2,[transformed_params_true[3]], seriestype="hline", linestyle=:dash, linecolor=:red, label = "optimum")

    p3 = plot(u_i[3,:], u_i[1,:], seriestype=:scatter, xlims = extrema(u_init[3,:]), ylims = extrema(u_init[1,:]))
    plot!(p3,[transformed_params_true[3]], xaxis="u1", yaxis="u2", seriestype="vline",
        linestyle=:dash, linecolor=:red, label = false,
        title = "EKI iteration = " * string(i)
        )
    plot!(p3,[transformed_params_true[1]], seriestype="hline", linestyle=:dash, linecolor=:red, label = "optimum")

    p=plot(p1, p2, p3, layout=(1,3))
    display(p)
    sleep(0.5)
end


text:
using Test, NiLang, NiLang.AD

@testset "complex GVar" begin
    a = 1.0+ 2im
    @test GVar(a) == Complex(GVar(1.0), GVar(2.0))
    @test GVar(a, a) == Complex(GVar(1.0, 1.0), GVar(2.0, 2.0))

    gx = GVar(1.0 + 1.0im)
    gx2 = chfield(gx, grad, 1.0+0.0im)
    @test gx2 == Complex(GVar(1.0, 1.0), GVar(1.0, 0.0))
end

@i function fr(f, loss, args...; il)
    f(args...)
    loss += (args |> tget(il)).re
end
@i function fi(f, loss, args...; il)
    f(args...)
    loss += (args |> tget(il)).im
end
function ccheck_grad(f, args; verbose=true, iloss=1)
    check_grad(fr, (f, 0.0, args...); verbose=verbose, iloss=2, il=1) &&
     check_grad(fi, (f, 0.0, args...); verbose=verbose, iloss=2, il=1)
end

@testset "check grad" begin
    x = 1.0 - 4.0im
    y = 2.0 - 2.3im
    z = 3.0 + 1.0im
    r = 4.0
    for opm in [⊕, ⊖]
        @test check_inv(opm(complex), (1+2.0im, 2.0, 3.0); verbose=true)
        @test ccheck_grad(opm(complex), (1+2.0im, 2.0, 3.0); verbose=true, iloss=1)
        for (subop, args) in [
            (opm(identity), (x,y)), (opm(+), (x, y, z)),
            (opm(-), (x, y, z)), (opm(*), (x, y, z)),
            (opm(/), (x, y, z)), (opm(^), (x, y, r)),
            (opm(exp), (x, y)), (opm(log), (x, y)),
            (opm(inv), (x, y))
            ]
            @test ccheck_grad(subop, args; verbose=true, iloss=1)
            r1 = subop(args...)
            r2 = [(opm == (⊕) ? Base.:+ : Base.:-)(args[1], subop.f(args[2:end]...)), args[2:end]...]
            @test all(r1 .≈ r2)
        end

        for (subop, args) in [
            (opm(angle), (r, y)), (opm(abs), (r, y)), (opm(abs), (r, 0.0im)),
            (opm(abs2), (r, y))
            ]
            @show subop, args
            r1 = [subop(args...)...]
            r2 = [(opm == (⊕) ? Base.:+ : Base.:-)(args[1], subop.f(args[2:end]...)), args[2:end]...]
            @test r1 ≈ r2
            @test check_grad(subop, args; verbose=true, iloss=1)
        end
    end
    for op in [-, NEG]
        @test check_inv(op, (x,); verbose=true)
        @test ccheck_grad(op, (x,); verbose=true, iloss=1)
    end
end

ext: jl


text:
# this is likely to be slower but easy to see that it's correct

function move_combine_compact!(x)
    pts_scored = 0
    updated = false

    if all(==(0), x)
        return x, pts_scored, updated
    end

    # combine phase
    @inbounds for i in 1:3
        if 0 < x[i] == x[i+1]
            updated = true
            x[i] += 1
            x[i+1] = 0
            pts_scored += 1 << x[i]
        end
    end

    # compact phase i.e. remove 0 in between
    @inbounds if x[1] == 0
        updated = true
        x[1], x[2], x[3], x[4] = x[2], x[3], x[4], 0
    end

    @inbounds if x[2] == 0
        updated = true
        x[2], x[3], x[4] = x[3], x[4], 0
    end

    @inbounds if x[3] == 0
        updated = true
        x[3], x[4] =  x[4], 0
    end
    x, pts_scored, updated
end


text:
using MatrixCore
using Test

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


text:
using VideoPlayer
using Test

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


text:
function wavenumbers_refine(ω::T, medium::Medium{T}, species::Species{T};
        tol::T = 1e-7,
        k_effs::Vector{Complex{T}} = Complex{T}[],
        kws...) where T<:Number

    # the dispersion equation is given by: `dispersion(k1,k2) = 0` where k_eff = k1 + im*k2.
    dispersion = dispersion_equation(ω, medium, species; tol = tol, kws...)

    # add any specified keffs
    k_vecs = [[real(kp),imag(kp)] for kp in k_effs]

    k_vecs = map(k_vecs) do k_vec
       res = optimize(dispersion, k_vec, Optim.Options(iterations = 10 * Int(abs(round(log(tol)))), g_tol = tol^3.0, x_tol=tol))
       if res.minimum > T(100)*tol
           [zero(T),-one(T)]
       else
           res.minimizer
       end
    end
    k_vecs = reduce_kvecs(k_vecs, tol)

    # Delete unphysical waves, including waves travelling backwards with almost no attenuation. This only is important in the limit of very low frequency or very weak scatterers.
    deleteat!(k_vecs, findall([k_vec[2] < -sqrt(tol) for k_vec in k_vecs]))

    k_effs = [ k_vec[1] + k_vec[2]*im for k_vec in k_vecs]
    k_effs = sort(k_effs, by=imag)

    return k_effs
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
# common directories and files for all tests
const fake_data = "../data/foo.conllu"
const real_tdata = "/state/partition1/tmp/okirnap/UD_English/en-ud-train.conllu"
const chmodel = "../data/english_chmodel.jld" # converted chmodelfile

# to work with biyofiz version of the data
const real_tdata3 = "/mnt/ai/data/nlp/conll17/ud-treebanks-v2.0/UD_English-LinES/en_lines-ud-train.conllu"
const real_ddata3 = "/mnt/ai/data/nlp/conll17/ud-treebanks-v2.0/UD_English-LinES/en_lines-ud-dev.conllu"

##### Dead directories #####
# const real_ddata = "/state/partition1/tmp/okirnap/UD_English/en-ud-dev.conllu"
# const real_tdata2 = "/state/partition1/tmp/okirnap/UD_English-LinES/en_lines-ud-train.conllu"
# const real_ddata2 = "/state/partition1/tmp/okirnap/UD_English-LinES/en_lines-ud-dev.conllu"
# const sc_tdata = "/ai/data/nlp/conll17/ud-treebanks-v2.0/UD_English/en-ud-train.conllu"
# const sc_ddata = "/ai/data/nlp/conll17/ud-treebanks-v2.0/UD_English/en-ud-dev.conllu"
# const sc_chmodel = "/Users/okirnap/charNL/trained_models/converged/english_chmodel.jld"


text:
# generic fallback behaviour
DiscreteFunction(dom, cod) = error("not implemented")
# specific constructions
DiscreteFunction(dom :: Segment, cod :: Segment) = ResFun(dom, cod)
DiscreteFunction(dom :: AbstractProduct, cod :: Segment) = ExtFun(dom, cod)
DiscreteFunction(dom :: AbstractProduct, cod :: AbstractProduct) = FTuple(dom, cod)

# generic fallback
#is_composable(dom, cod) = false
#is_composable(rng :: Segment, dom :: Segment) = true
#is_composable(rng :: Segment, dom :: AbstractProduct) = false
#is_composable(rng :: DProd, dom :: DProd) =

function ∘(f1 :: AbstractDiscreteFunction, f2 :: AbstractDiscreteFunction)
    dom = domain(f2)
    cod = codomain(f1)
    g = DiscreteFunction(dom, cod)
    for x in dom
        g[x] = f1(f2(x))
    end
    return g
end

max_stars_repo_name: kirtsar/DiscreteFunctions

max_issues_repo_name: kirtsar/DiscreteFunctions


text:
# Stricker Klaus 2020 - License is MIT: http://julialang.org/license
# ReadWriteDlm2 - rwd2_write.jl - https://github.com/strickek/ReadWriteDlm2.jl

"""

    writedlm2(f, A; opts...)
    writedlm2(f, A, delim; opts...)

Write `A` (a vector, matrix, or an iterable collection of iterable rows, a
`Tables` source) as text to `f` (either a filename or an IO stream). The columns
are separated by `';'`, another `delim` (Char or String) can be defined.

By default, a pre-processing of values takes place. Before writing as strings,
decimal marks are changed from `'.'` to `','`. With the keyword argument
`decimal=` another decimal mark can be defined.
To switch off this pre-processing set: `decimal='.'`.

In `writedlm2()` the output format for `Date` and `DateTime` data can be
defined with format strings. Defaults are the ISO formats. Day (`E`, `e`)
and month (`U`, `u`) names are written in the `locale` language. For writing
`Complex` numbers the imaginary component suffix can be selected with the
`imsuffix=` keyword argument.

# Additional Keyword Arguments

* `decimal=','`: Character for writing decimal marks
* `dtfs=\"yyyy-mm-ddTHH:MM:SS.s\"`: DateTime write format
* `dfs=\"yyyy-mm-dd\"`: Date write format
* `locale=\"english\"`: Language for DateTime writing
* `imsuffix=\"im\"`: Complex Imag suffix `\"im\"`, `\"i\"` or `\"j\"`
* `missingstring=\"na\"`: How missing values are written

# Code Example
```jldoctest
julia> using ReadWriteDlm2, Dates

julia> A = Any[1 1.2; "text" Date(2017)];

julia> writedlm2("test.csv", A)

julia> read("test.csv", String)
"1;1,2\\ntext;2017-01-01\\n"
```
"""
writedlm2(io::IO, a; opts...) =
    writedlm2auto(io, a, ';'; opts...)

writedlm2(io::IO, a, dlm; opts...) =
    writedlm2auto(io, a, dlm; opts...)

writedlm2(f::AbstractString, a; opts...) =
    writedlm2auto(f, a, ';'; opts...)

writedlm2(f::AbstractString, a, dlm; opts...) =
    writedlm2auto(f, a, dlm; opts...)

function writedlm2auto(f, a, dlm;
    decimal::AbstractChar=',',
    dtfs::AbstractString="yyyy-mm-ddTHH:MM:SS.s",
    dfs::AbstractString="yyyy-mm-dd",
    locale::AbstractString="english",
    imsuffix::AbstractString="im",
    missingstring::AbstractString="na",
    opts...)

    ((!isempty(dtfs) && !occursin(Regex("[^YymdHMSs]"), dtfs)) ||
    (!isempty(dfs) && !occursin(Regex("[^YymdHMSs]"), dfs))) && info(
    """
    Format string for DateTime(`$dtfs`) or Date(`$dfs`)
    contains numeric code elements only. At least one non-numeric
    code element or character is needed for parsing dates.
    """)

    (string(dlm) == string(decimal)) && error(
        "Error: decimal = delim = ´$(dlm)´ - change decimal or delim!")

    ((imsuffix != "im") && (imsuffix != "i") && (imsuffix != "j")) && error(
        "Only `\"im\"`, `\"i\"` or `\"j\"` are valid for `imsuffix`.")

    if isa(a, Union{Nothing, Missing, Number, TimeType})
         a = [a]  # create 1 element Array
    elseif Tables.istable(a) == true # Tables interface
        a = matrix2(a)  # Matrix{Any}, first row columnnames, row 2:end -> data
    end

    if isa(a, AbstractArray)
        fdt = !isempty(dtfs)  # Bool: format DateTime
        dtdf = DateFormat(dtfs, locale)
        fd = !isempty(dfs)    # Bool: format Date
        ddf = DateFormat(dfs, locale)
        ft = (decimal != '.')   # Bool: format Time (change decimal)

        # create b for manipulation/write - keep a unchanged
        b = similar(a, Any)
        for i in eachindex(a)
            b[i] =
            isa(a[i], AbstractFloat) ? floatformat(a[i], decimal) :
            isa(a[i], Missing) ? missingstring :
            isa(a[i], Nothing) ? "nothing" :
            isa(a[i], DateTime) && fdt ? Dates.format(a[i], dtdf) :
            isa(a[i], Date) && fd ? Dates.format(a[i], ddf) :
            isa(a[i], Time) && ft ? timeformat(a[i], decimal) :
            isa(a[i], Complex) ? complexformat(a[i], decimal, imsuffix) :
            string(a[i])
        end
        else  # a is not a Number, TimeType or Array -> no preprocessing
        b = a
    end

    writedlm(f, b, dlm; opts...)

end # End function writedlm2auto()


text:
############################################################################################
"""
$(TYPEDEF)

Utility struct to help assign boundaries to parameter - keeps parameter unconstrained. Useful for assigning buffer values for functions of parameter.

# Fields
$(TYPEDFIELDS)
"""
struct Unconstrained <: AbstractConstraint
    function Unconstrained()
        return new()
    end
end

############################################################################################
function construct_transform(constrained::Unconstrained, val)
    return Bijectors.Identity{0}(), Bijectors.inverse(Bijectors.Identity{0}())
end

############################################################################################
function _check(
    _rng::Random.AbstractRNG,
    constraint::Unconstrained,
    val::Union{R,Array{R},AbstractArray}
) where {R<:Real}
    return true
end

############################################################################################
# Unconstrained constraint -> just flatten
function construct_flatten(
    output::Type{T},
    flattentype::F,
    unflattentype::U,
    constraint::Unconstrained,
    val
) where {T<:AbstractFloat,F<:FlattenTypes,U<:UnflattenTypes}
    return construct_flatten(T, flattentype, unflattentype, val)
end

############################################################################################
#=
!NOTE: For this constraint, we use a Bijector as Transformer, so we do not need to add a new functors
1. MyTransformer <: AbstractTransformer, MyInverseTransformer <: AbstractTransformer
2. define a function construct_transform(MyConstraint, val) -> MyTransformer, MyInverseTransformer
3. overload unconstrain and log_abs_det_jac on MyTransformer.
4. overload constrain on MyInverseTransformer.
=#

############################################################################################
#Export
export
    Unconstrained,
    construct_flatten,
    construct_transform,
    constrain,
    unconstrain,
    log_abs_det_jac


text:
#i64o-shifts.jl

#bitshifting operations on superints
#iterative leftshift and rightshift operations on Array SuperInts

lsh(a::UInt64, b::Int64) = a << b
lsh(a::UInt64, b::UInt16) = a << b

function lsh!{FSS}(a::ArrayNum{FSS}, b::Int64)
  (b < 0) && return rsh!(a, -b)
  lsh!(a, to16(b))
end
#destructive version which clobbers the existing array
function lsh!{FSS}(a::ArrayNum{FSS}, b::UInt16)
  #kick it back to right shift if it's negative
  l = __cell_length(FSS)
  #calculate how many cells apart our two ints shall be.
  celldiff::UInt16 = b >> 6
  #calculate how much we have to shift
  shift::UInt16 = b & 0x003F
  shift == 0 && @goto cellmove          #skip it, if it's more efficient.
  c_shift::UInt16 = 0x0040 - shift

  #go ahead and shift all blocks.
  for idx = 1:l-1
    @inbounds a.a[idx] = (a.a[idx] << shift) | a.a[idx + 1] >> c_shift
  end
  #finish the last block.
  @inbounds a.a[l] = a.a[l] << shift
  #cut in here if we skipped the cell movement.

  @label cellmove
  #move the cells
  (celldiff == 0) && return
  splitdex::Int64 = l - celldiff
  for idx = 1:l
    @inbounds a.a[idx] = (idx <= splitdex) ? a.a[idx + celldiff] : 0
  end
  return a
end

@fracproc lsh shft

rsh(a::UInt64, b::Int64) = a >> b
rsh(a::UInt64, b::UInt16) = a >> b
rsh!(a::UInt64, b::UInt16) = rsh!(a, Int64(b))

function rsh!{FSS}(a::ArrayNum{FSS}, b::Int64)
  (b < 0) && return lsh!(a, -b)
  rsh!(a, to16(b))
end
function rsh!{FSS}(a::ArrayNum{FSS}, b::UInt16)
  #kick it back to right shift if it's negative

  l = __cell_length(FSS)
  #calculate how many cells apart our two ints shall be.
  celldiff::UInt16 = b >> 6
  #calculate how much we have to shift
  shift::UInt16 = b & 0x003F
  shift == 0 && @goto cellmove          #skip it, if it's more efficient.
  c_shift::UInt16 = 0x0040 - shift

  #go ahead and shift all blocks.
  for idx = l:-1:2
    @inbounds a.a[idx] = (a.a[idx] >> shift) | a.a[idx - 1] << c_shift
  end

  #finish the last block.
  @inbounds a.a[1] = a.a[1] >> shift
  #cut in here if we skipped the cell movement.
  @label cellmove
  #move the cells
  (celldiff == 0) && return
  splitdex::Int64 = celldiff

  for idx = l:-1:1
    @inbounds a.a[idx] = (idx > splitdex) ? a.a[idx - celldiff] : 0
  end
  a
end

@fracproc rsh shft

################################################################################
## rightshift with underflow check.

@generated function rsh_underflow_check{FSS}(num::UInt64, shift::UInt16, ::Type{Val{FSS}})
  mfsize = Unums.max_fsize(FSS)
  quote
    mask = ((shift > $mfsize) * f64) | mask_bot($mfsize - shift)
    ubit_needed = ((num & mask) != z64) * UNUM_UBIT_MASK
    num = rsh(num, shift)
    (num, ubit_needed)
  end
end

@generated function underflow_check{FSS}(num::ArrayNum{FSS}, shift::UInt16)
  l = __cell_length(FSS)
  mfsize = Unums.max_fsize(FSS)
  quote
    shift = (shift > $mfsize) ? ($mfsize + o16) : shift
    middle_cell = $l - div(shift, 0x0040)
    middle_mask = 0x0040 - shift % 0x0040
    for idx = __cell_length(FSS):-1:(middle_cell + 1)
      @inbounds (num.a[idx] == 0) || return UNUM_UBIT_MASK
    end
    @inbounds res = ((num.a[middle_cell] & mask_bot(middle_mask)) != z64) * UNUM_UBIT_MASK
    return res
  end
end

function rsh_underflow_check!{FSS}(num::ArrayNum{FSS}, shift::UInt16)
  ubit_needed = underflow_check(num, shift)
  rsh!(num, shift)
  ubit_needed
end

################################################################################
## frac versions of this function

function frac_rsh_underflow_check!{ESS, FSS}(x::UnumSmall{ESS,FSS}, shift::UInt16)
  (x.fraction, ubit_needed) = rsh_underflow_check(x.fraction, shift, Val{FSS})
  x.flags |= ubit_needed
  return x
end

function frac_rsh_underflow_check!{ESS, FSS}(x::UnumLarge{ESS,FSS}, shift::UInt16)
  ubit_needed = rsh_underflow_check!(x.fraction, shift)
  x.flags |= ubit_needed
  return x
end

################################################################################
## shift with carry
doc"""
  `shift_with_carry(carry::UInt64, fraction::UInt64, shift::UInt16)`
  `shift_with_carry(carry::UInt64, fraction::ArrayNum, shift::UInt16)`

  shifts a fraction by shift spaces to the left, and copies that number of
  bits into the fraction at the shift position.
"""
function shift_with_carry(carry::UInt64, fraction::UInt64, shift::UInt16)
  fraction = rsh(fraction, shift)
  carry_xfer = carry << (0x0040 - shift)
  carry >>= shift
  old_fraction = fraction
  fraction += carry_xfer
  (old_fraction > fraction) && (carry += o64)
  return (carry, fraction)
end
function shift_with_carry{FSS}(carry::UInt64, fraction::ArrayNum{FSS}, shift::UInt16)
  rsh!(fraction, shift)
  shiftbits = shift % 0x0040
  shiftcell = (shift ÷ 0x0040) + 1
  carry_xfer = carry << (0x0040 - shiftbits)
  carry >>= shift
  @inbounds begin
    oldvalue = fraction.a[shiftcell]
    fraction.a[shiftcell] += carry_xfer
    (oldvalue > fraction.a[shiftcell]) && (carry += o64)
  end
  if (shiftcell > 1)
    @inbounds fraction.a[shiftcell - 1] = carry
    carry = 0
  end
  return (carry, fraction)
end

@universal function frac_shift_with_carry!(carry::UInt64, x::Unum, shift::UInt16)
  (carry, x.fraction) = shift_with_carry(carry, x.fraction, shift)
  return carry
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
int main() {

  int[] a = new int[10];
  int j=0;
  while (j<a.length) {
     a[j] = j;
     j++;
  }

  for (int x : a)
     printInt(x);

  int x = 45;
  printInt(x);
  return 0;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
# using Brillouin: irrfbz_path
using Compat: eachslice
using Counters: counter
using LinearAlgebra: cross
using Spglib: get_ir_reciprocal_mesh, get_spacegroup_type

export ReciprocalPoint,
    ReciprocalLattice,
    MonkhorstPackGrid,
    # ReciprocalPath,
    reciprocal_mesh,
    coordinates,
    weights

"""
    ReciprocalLattice(mat::SMatrix)

Construct a `ReciprocalLattice`.

!!! warning
    You should not use this function directly, always use `inv` of a `Lattice`.
"""
struct ReciprocalLattice{T} <: AbstractLattice{T}
    data::SMatrix{3,3,T,9}
end
@functor ReciprocalLattice

"""
    inv(lattice::Lattice)
    inv(lattice::ReciprocalLattice)

Get the reciprocal of a `Lattice` or a `ReciprocalLattice`.
"""
function Base.inv(lattice::Lattice)
    Ω = det(lattice.data)  # Cannot use `cellvolume`, it takes the absolute value!
    𝐚, 𝐛, 𝐜 = basis_vectors(lattice)
    return ReciprocalLattice(
        inv(Ω) * transpose(hcat(cross(𝐛, 𝐜), cross(𝐜, 𝐚), cross(𝐚, 𝐛))),
    )
end
function Base.inv(lattice::ReciprocalLattice)
    Ω⁻¹ = det(lattice.data)  # Cannot use `cellvolume`, it takes the absolute value!
    𝐚⁻¹, 𝐛⁻¹, 𝐜⁻¹ = basis_vectors(lattice)
    return Lattice(inv(Ω⁻¹) * hcat(cross(𝐛⁻¹, 𝐜⁻¹), cross(𝐜⁻¹, 𝐚⁻¹), cross(𝐚⁻¹, 𝐛⁻¹)))
end

"""
    basis_vectors(lattice::ReciprocalLattice)

Get the three basis vectors from a `ReciprocalLattice`.
"""
basis_vectors(lattice::ReciprocalLattice) = lattice[1, :], lattice[2, :], lattice[3, :]

"""
    ReciprocalPoint(x, y, z, w)

Represent a special point of the 3D Brillouin zone. Each of them has a weight `w`.
"""
struct ReciprocalPoint{T}
    coord::SVector{3,T}
    weight::Float64
end
ReciprocalPoint(coord::AbstractVector{T}, weight) where {T} =
    ReciprocalPoint{T}(SVector{3}(coord), weight)
ReciprocalPoint(x, y, z, w) = ReciprocalPoint(SVector(x, y, z), w)
@functor ReciprocalPoint (coord,)

"""
    MonkhorstPackGrid(mesh, is_shift)

Represent the Monkhorst--Pack grid.

# Arguments
- `mesh`: A length-three vector specifying the k-point grid (``nk_1 × nk_2 × nk_3``) as in Monkhorst--Pack grids.
- `is_shift`: A length-three vector specifying whether the grid is displaced by half a grid step in the corresponding directions.
"""
struct MonkhorstPackGrid
    mesh::SVector{3,UInt}
    is_shift::SVector{3,Bool}
    function MonkhorstPackGrid(mesh, is_shift)
        @assert all(mesh .>= 1)
        if eltype(is_shift) != Bool
            is_shift = Bool.(is_shift)
        end
        return new(mesh, is_shift)
    end
end

# See example in https://spglib.github.io/spglib/python-spglib.html#get-ir-reciprocal-mesh
"""
    reciprocal_mesh(cell::Cell, mesh, is_shift; kwargs...)

List the `ReciprocalPoint`s from the mesh of the reciprocal space of a `Cell`.

# Arguments
- `cell::Cell`: the cell.
- `mesh`: a vector of three integers which specify the mesh numbers along reciprocal primitive axis.
- `is_shift=falses(3)`: a vector of three elements specifying whether the mesh is shifted along the axis in half of adjacent mesh points irrespective of the mesh numbers. The allowed values are `0`, `1`, `true`, and `false`.
- `is_time_reversal=true`: whether to impose the time reversal symmetry on the mesh.
- `symprec=1e-5`: distance tolerance in Cartesian coordinates to find crystal symmetry.
- `cartesian=false`: whether to return the reciprocal points in Cartesian coordinates.
- `ir_only=true`: whether to return the reciprocal points only in the irreducible Brillouin zone.
"""
function reciprocal_mesh(
    cell::Cell,
    mesh,
    is_shift = falses(3);
    is_time_reversal = true,
    symprec = 1e-5,
    cartesian = false,
    ir_only = true,
)
    _, mapping, grid = get_ir_reciprocal_mesh(
        cell,
        mesh,
        is_shift;
        is_time_reversal = is_time_reversal,
        symprec = symprec,
    )
    shift = is_shift ./ 2  # true / 2 = 0.5, false / 2 = 0
    mapping = convert(Vector{Int}, mapping)
    weights = counter(mapping)
    total_number = length(mapping)  # Number of all k-points, not only the irreducible ones
    crystal_coord = if ir_only
        map(unique(mapping)) do index
            x, y, z = (grid[:, index] .+ shift) ./ mesh
            weight = weights[index] / total_number
            ReciprocalPoint(x, y, z, weight)
        end
    else
        map(eachslice(grid; dims = 2)) do point
            x, y, z = (point .+ shift) ./ mesh
            weight = 1 / total_number
            ReciprocalPoint(x, y, z, weight)
        end
    end
    if cartesian
        t = CartesianFromFractional(inv(Lattice(cell)))
        return map(t, crystal_coord)
    else
        return crystal_coord
    end
end
reciprocal_mesh(cell::Cell, mp::MonkhorstPackGrid; kwargs...) =
    reciprocal_mesh(cell, mp.mesh, mp.is_shift; kwargs...)

"""
    coordinates(arr::AbstractArray{<:ReciprocalPoint})

Get the coordinates of an array of `ReciprocalPoint`s.
"""
coordinates(arr::AbstractArray{<:ReciprocalPoint}) = map(x -> x.coord, arr)

"""
    weights(arr::AbstractArray{<:ReciprocalPoint})

Get the weights of an array of `ReciprocalPoint`s.
"""
weights(arr::AbstractArray{<:ReciprocalPoint}) = map(x -> x.weight, arr)

function Base.show(io::IO, x::ReciprocalPoint)
    if get(io, :compact, false) || get(io, :typeinfo, nothing) == typeof(x)
        Base.show_default(IOContext(io, :limit => true), x)  # From https://github.com/mauro3/Parameters.jl/blob/ecbf8df/src/Parameters.jl#L556
    else
        println(io, string(typeof(x)))
        print(io, " coord = ", x.coord, ", weight = ", x.weight)
    end
end

# struct ReciprocalPath{N}
#     special_points::Dict{Symbol,SVector{N,Float64}}
#     suggested_paths::Vector{Vector{Symbol}}
#     lattice::Lattice
# end
# """
#     ReciprocalPath(lattice::Lattice, spgnum::Integer)
#     ReciprocalPath(cell::Cell)

# Construct a `ReciprocalPath` from a `Lattice` or a `Cell`.
# """
# function ReciprocalPath(lattice::Lattice, spgnum::Integer)
#     kpath = irrfbz_path(spgnum, collect(basis_vectors(lattice)))
#     return ReciprocalPath(kpath.points, kpath.paths, lattice)
# end
# function ReciprocalPath(cell::Cell)
#     spg = get_spacegroup_type(cell)
#     return ReciprocalPath(Lattice(cell), spg.number)
# end

"""
    coordinates(path::ReciprocalPath, cartesian = false)

If `cartesian` is `true`, return the coordinates in the Cartesian coordinate system.
"""
# coordinates(path::ReciprocalPath, cartesian = false) =
#     cartesian ?
#     Dict(
#         key => CartesianFromFractional(inv(path.lattice))(value) for
#         (key, value) in path.special_points
#     ) : path.special_points


text:
using Random
using LinearAlgebra
"""
dir_sample(dir_type::String, d::Int, m::Int, seed::Int)

"dir_sample"
samples "m" directions in dimension "d" of type "dir_type"

Inputs:
(1) dir_type: a string, the direction type
(2) d: an integer, the dimension in which the directions are sampled
(3) m: an integer, the number of directions to generate
(4) seed: an integer, used for random sampling

Output:
an m x d matrix, where each row is a direction in R^d
"""
function dir_sample(dir_type::String, d::Int, m::Int, seed::Int,
    theta = pi/4)
    ## set the seed for random number generation
    Random.seed!(seed)
    if dir_type=="pos"
        directions = zeros(0,d)
        num_row = 0
        while num_row<m
            direction = randn(1,d)
            if sum(direction.>0)==d
                directions = [directions; direction]
                num_row+=1
            end
        end
        return directions
    elseif dir_type=="unif"
        directions = randn(m,d)
        return directions
    elseif dir_type=="angle"
        directions = zeros(0,d)
        num_row = 0
        vec_1 = ones(d)/sqrt(d) ## the vector [1,1,...,1]\in\R^d
        num_row = 0
        while num_row<m
            direction = randn(1,d)
            direction = direction/norm(direction)
            if acos(dot(vec_1, direction'))<theta
                directions = [directions; direction]
                num_row+=1
            end
        end
        return directions
    end
end

# ## sample usage code
# dir_type = "pos" ## "unif"
# d = 2
# m = 100
# seed = 1
# A = dir_sample(dir_type::String, d::Int, m::Int, seed::Int)

# using Plots
# pyplot()

# plot()
# for i = 1:size(A,1)
#     plot!([0, A[i,1]], [0,A[i,2]], legend = false)
# end
# plot!()

"""
"UnifRandBall"
generates a d x N matrix, whose columns represent points,
randomly uniformly sampled in the unit ball of R^d
"""
function UnifRandBall(d::Int,N::Int)
    Collect = zeros(d,0) ## Set an (d x 0) empty matrix for storing
    i = 0 ## number of collected points
    while i<N
        pt = 2*rand(d).-1 ## (uniformly) generate a random point inside [-1,1]^d
        if norm(pt)<1 ## if the point is inside the unit ball, collect it
            i = i+1
            Collect = [Collect pt]
        end
    end
    return Collect
end

"""
pt_sample(pt_type::String, d::Int, n::Int, seed::Int)

"pt_sample"
samples "n" points in R^d of type "pt_type"

Inputs:
(1) pt_type: a string, the point type
(2) d: an integer, the dimension in which the directions are sampled
(3) n: an integer, the number of directions to generate
(4) seed: an integer, used for random sampling
(5) r(optional): the reciprocal of the curvature of the de-sitter space, between 0 and 1, exclusive

Output:
a d x n matrix, where each row is a direction in R^d
"""
function pt_sample(pt_type::String, d::Int, n::Int, seed::Int,
        r=0.1::Float64)
    ## set the seed for random number generation
    Random.seed!(seed)
    if pt_type=="uniform"
        pts = UnifRandBall(d,n)
        return pts
    elseif pt_type=="sphere"
        pts = UnifRandBall(d,n)
        for a = 1:n
            pts[:,a] = pts[:,a]/norm(pts[:,a])
        end
        return pts
    elseif pt_type=="de-sitter"
        c = 1/r ## the constant/curvature c in x_1^2+...+x_(d-1)^2-x_d^2 = 1/c^2
        pts = zeros(d,0)
        a = 0
        while a<n
            new_pt = UnifRandBall(d-1,1)
            x_d_a2 = sum(new_pt.^2)-1/c^2
            if x_d_a2>0
                new_pt = [new_pt;(-1)^a*sqrt(x_d_a2)] ## (-1)^a is used to cover both the upper and lower piece
                pts = [pts new_pt]
                a+=1
            end
        end
        return pts
    end
end

# ## sample usage code
# pt_type = "sphere" ## "uniform"
# d = 2
# n = 100
# seed = 1
# X = pt_sample(pt_type::String, d::Int, n::Int, seed::Int)

# using Plots
# pyplot()

# scatter(X[1,:], X[2,:], legend = false)


text:
module VectorizationBase

import ArrayInterface, LinearAlgebra, Libdl, Hwloc, IfElse
using ArrayInterface:
    StaticInt, Zero, One, StaticBool, True, False,
    contiguous_axis, contiguous_axis_indicator, contiguous_batch_size, stride_rank,
    device, CPUPointer, CPUIndex, eq, ne, lt, le, gt, ge,
    known_length, known_first, known_last, strides, offsets,
    static_first, static_last, static_length
import IfElse: ifelse

asbool(::Type{True}) = true
asbool(::Type{False}) = false
# TODO: see if `@inline` is good enough.
# @inline asvalbool(r) = Val(map(Bool, r))
# @inline asvalint(r) = Val(map(Int, r))
@generated function asvalint(r::T) where {T<:Tuple{Vararg{StaticInt}}}
    t = Expr(:tuple)
    for s ∈ T.parameters
        push!(t.args, s.parameters[1])
    end
    Expr(:call, Expr(:curly, :Val, t))
end
@generated function asvalbool(r::T) where {T<:Tuple{Vararg{StaticBool}}}
    t = Expr(:tuple)
    for b ∈ T.parameters
        push!(t.args, b === True)
    end
    Expr(:call, Expr(:curly, :Val, t))
end
@inline val_stride_rank(A) = asvalint(stride_rank(A))
@inline val_dense_dims(A) = asvalbool(ArrayInterface.dense_dims(A))

# doesn't export `Zero` and `One` by default, as these names could conflict with an AD library
export Vec, Mask, EVLMask, MM, stridedpointer, vload, vstore!, StaticInt, True, False, Bit,
    vbroadcast, mask, vfmadd, vfmsub, vfnmadd, vfnmsub,
    VecUnroll, Unroll, pick_vector_width

using Base: llvmcall, VecElement, HWReal, tail
const LLVMCALL = GlobalRef(Base, :llvmcall)

const FloatingTypes = Union{Float32, Float64} # Float16

const SignedHW = Union{Int8,Int16,Int32,Int64}
const UnsignedHW = Union{UInt8,UInt16,UInt32,UInt64}
const IntegerTypesHW = Union{SignedHW,UnsignedHW}
const IntegerTypes = Union{StaticInt,IntegerTypesHW}

struct Bit; data::Bool; end # Dummy for Ptr
const Boolean = Union{Bit,Bool}
const NativeTypesExceptBit = Union{Bool,HWReal}
const NativeTypes = Union{NativeTypesExceptBit, Bit}

const _Vec{W,T<:Number} = NTuple{W,Core.VecElement{T}}

abstract type AbstractSIMD{W,T <: Union{<:StaticInt,NativeTypes}} <: Real end
abstract type AbstractSIMDVector{W,T} <: AbstractSIMD{W,T} end
struct VecUnroll{N,W,T,V<:Union{NativeTypes,AbstractSIMD{W,T}}} <: AbstractSIMD{W,T}
    data::Tuple{V,Vararg{V,N}}
    @inline (VecUnroll(data::Tuple{V,Vararg{V,N}})::VecUnroll{N,W,T,V}) where {N,W,T,V<:AbstractSIMD{W,T}} = new{N,W,T,V}(data)
    @inline (VecUnroll(data::Tuple{T,Vararg{T,N}})::VecUnroll{N,T,T}) where {N,T<:NativeTypes} = new{N,1,T,T}(data)
    # # following two definitions are for checking that you aren't accidentally creating `VecUnroll{0}`s.
    # @inline (VecUnroll(data::Tuple{V,Vararg{V,N}})::VecUnroll{N,W,T,V}) where {N,W,T,V<:AbstractSIMD{W,T}} = (@assert(N > 0); new{N,W,T,V}(data))
    # @inline (VecUnroll(data::Tuple{T,Vararg{T,N}})::VecUnroll{N,T,T}) where {N,T<:NativeTypes} = (@assert(N > 0); new{N,1,T,T}(data))

    # @inline VecUnroll{N,W,T,V}(data::Tuple{V,Vararg{V,N}}) where {N,W,T,V<:AbstractSIMDVector{W,T}} = new{N,W,T,V}(data)
    # @inline (VecUnroll(data::Tuple{V,Vararg{V,N}})::VecUnroll{N,W,T,Vec{W,T}}) where {N,W,T,V<:AbstractSIMDVector{W,T}} = new{N,W,T,V}(data)
    # @inline (VecUnroll(data::Tuple{V,Vararg{V,N}})::VecUnroll{N,W,T,V}) where {N,W,T,V<:AbstractSIMDVector{W,T}} = new{N,W,T,V}(data)
end
# @inline VecUnroll(data::Tuple) = VecUnroll(promote(data))
# const AbstractSIMD{W,T} = Union{AbstractSIMDVector{W,T},VecUnroll{<:Any,W,T}}
const VecOrScalar = Union{AbstractSIMDVector,NativeTypes}
const NativeTypesV = Union{AbstractSIMD,NativeTypes,StaticInt}
# const NativeTypesV = Union{AbstractSIMD,NativeTypes,StaticInt}
const IntegerTypesV = Union{AbstractSIMD{<:Any,<:IntegerTypes},IntegerTypesHW}
struct Vec{W,T} <: AbstractSIMDVector{W,T}
    data::NTuple{W,Core.VecElement{T}}
    @inline Vec{W,T}(x::NTuple{W,Core.VecElement{T}}) where {W,T<:NativeTypes} = new{W,T}(x)
    @generated function Vec(x::Tuple{Core.VecElement{T},Vararg{Core.VecElement{T},_W}}) where {_W,T<:NativeTypes}
        W = _W + 1
        # @assert W === pick_vector_width(W, T)# || W === 8
        vtyp = Expr(:curly, :Vec, W, T)
        Expr(:block, Expr(:meta,:inline), Expr(:(::), Expr(:call, vtyp, :x), vtyp))
    end
    # @inline function Vec(x::NTuple{W,<:Core.VecElement}) where {W}
    #     T = eltype(x)
    #     @assert W === pick_vector_width(W, T)
    #     # @assert ispow2(W) && (W ≤ max(pick_vector_width(W, T), 8))
    #     new{W,T}(x)
    # end
end

@inline Base.copy(v::AbstractSIMDVector) = v
@inline asvec(x::_Vec) = Vec(x)
@inline asvec(x) = x
@inline data(vu::VecUnroll) = getfield(vu, :data)

@inline unrolleddata(x) = x
@inline unrolleddata(x::VecUnroll) = getfield(x, :data)

@inline _demoteint(::Type{T}) where {T} = T
@inline _demoteint(::Type{Int64}) = Int32
@inline _demoteint(::Type{UInt64}) = UInt32

abstract type AbstractMask{W,U<:UnsignedHW} <: AbstractSIMDVector{W,Bit} end
struct Mask{W,U} <: AbstractMask{W,U}
    u::U
    @inline function Mask{W,U}(u::Unsigned) where {W,U} # ignores U...
        U2 = mask_type(StaticInt{W}())
        new{W,U2}(u % U2)
    end
end
struct EVLMask{W,U} <: AbstractMask{W,U}
    u::U
    evl::UInt32
    @inline function EVLMask{W,U}(u::Unsigned, evl) where {W,U} # ignores U...
        U2 = mask_type(Val{W}())
        new{W,U2}(u % U2, evl % UInt32)
    end
end
@inline Mask{W}(u::U) where {W,U<:Unsigned} = Mask{W,U}(u)
@inline EVLMask{W}(u::U, i) where {W,U<:Unsigned} = EVLMask{W,U}(u, i)
@inline Mask(m::EVLMask{W,U}) where {W,U} = Mask{W,U}(getfield(m,:u))
# Const prop is good enough; added an @inferred test to make sure.
# Removed because confusion can cause more harm than good.

@inline Base.broadcastable(v::AbstractSIMDVector) = Ref(v)

Vec{W,T}(x::Vararg{NativeTypes,W}) where {W,T<:NativeTypes} = Vec(ntuple(w -> Core.VecElement{T}(x[w]), Val{W}()))
Vec{1,T}(x::Union{Float32,Float64}) where {T<:NativeTypes} = T(x)
Vec{1,T}(x::Union{Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Bool}) where {T<:NativeTypes} = T(x)

@inline Base.length(::AbstractSIMDVector{W}) where W = W
@inline Base.size(::AbstractSIMDVector{W}) where W = (W,)
@inline Base.eltype(::AbstractSIMD{W,T}) where {W,T} = T
@inline Base.conj(v::AbstractSIMDVector) = v # so that things like dot products work.
@inline Base.adjoint(v::AbstractSIMDVector) = v # so that things like dot products work.
@inline Base.transpose(v::AbstractSIMDVector) = v # so that things like dot products work.

# Not using getindex/setindex as names to emphasize that these are generally treated as single objects, not collections.
@generated function extractelement(v::Vec{W,T}, i::I) where {W,I <: IntegerTypesHW,T}
    typ = LLVM_TYPES[T]
    instrs = """
        %res = extractelement <$W x $typ> %0, i$(8sizeof(I)) %1
        ret $typ %res
    """
    call = :($LLVMCALL($instrs, $T, Tuple{_Vec{$W,$T},$I}, data(v), i))
    Expr(:block, Expr(:meta, :inline), call)
end
@generated function insertelement(v::Vec{W,T}, x::T, i::I) where {W,I <: IntegerTypesHW,T}
    typ = LLVM_TYPES[T]
    instrs = """
        %res = insertelement <$W x $typ> %0, $typ %1, i$(8sizeof(I)) %2
        ret <$W x $typ> %res
    """
    call = :(Vec($LLVMCALL($instrs, _Vec{$W,$T}, Tuple{_Vec{$W,$T},$T,$I}, data(v), x, i)))
    Expr(:block, Expr(:meta, :inline), call)
end
@inline (v::AbstractSIMDVector)(i::IntegerTypesHW) = extractelement(v, i - one(i))
@inline (v::AbstractSIMDVector)(i::Integer) = extractelement(v, Int(i) - 1)
Base.@propagate_inbounds (vu::VecUnroll)(i::Integer, j::Integer) = getfield(vu, :data)[j](i)

@inline Base.Tuple(v::Vec{W}) where {W} = ntuple(v, Val{W}())

# Use with care in function signatures; try to avoid the `T` to stay clean on Test.detect_unbound_args

@inline data(v) = v
@inline data(v::Vec) = getfield(v, :data)

function Base.show(io::IO, v::AbstractSIMDVector{W,T}) where {W,T}
    name = typeof(v)
    print(io, "$(name)<")
    for w ∈ 1:W
        print(io, repr(extractelement(v, w-1)))
        w < W && print(io, ", ")
    end
    print(io, ">")
end
Base.bitstring(m::AbstractMask{W}) where {W} = bitstring(data(m))[end-W+1:end]
function Base.show(io::IO, m::AbstractMask{W}) where {W}
    bits = data(m)
    if m isa EVLMask
        print(io, "EVLMask{$W,Bit}<")
    else
        print(io, "Mask{$W,Bit}<")
    end
    for w ∈ 0:W-1
        print(io, bits & 1)
        bits >>= 1
        w < W-1 && print(io, ", ")
    end
    print(io, ">")
end
function Base.show(io::IO, vu::VecUnroll{N,W,T,V}) where {N,W,T,V}
    println(io, "$(N+1) x $V")
    d = data(vu)
    for n in 1:N+1
        show(io, d[n]);
        n > N || println(io)
    end
end

"""
The name `MM` type refers to _MM registers such as `XMM`, `YMM`, and `ZMM`.
`MMX` from the original MMX SIMD instruction set is a [meaningless initialism](https://en.wikipedia.org/wiki/MMX_(instruction_set)#Naming).

The `MM{W,X}` type is used to represent SIMD indexes of width `W` with stride `X`.
"""
struct MM{W,X,I<:Union{HWReal,StaticInt}} <: AbstractSIMDVector{W,I}
    i::I
    @inline MM{W,X}(i::T) where {W,X,T<:Union{HWReal,StaticInt}} = new{W,X::Int,T}(i)
end
@inline MM(i::MM{W,X}) where {W,X} = MM{W,X}(getfield(i, :i))
@inline MM{W}(i::Union{HWReal,StaticInt}) where {W} = MM{W,1}(i)
@inline MM{W}(i::Union{HWReal,StaticInt}, ::StaticInt{X}) where {W,X} = MM{W,X}(i)
@inline data(i::MM) = getfield(i, :i)

@inline extractelement(i::MM{W,X,I}, j) where {W,X,I} = getfield(i, :i) + (X % I) * (j % I)

function Base.getproperty(::AbstractSIMD, ::Symbol)
    throw(ErrorException("""
`Base.getproperty` not defined on AbstractSIMD.
If you wish to work with the data as a tuple, it is recommended to use `Tuple(v)`. Once you have an ordinary tuple, you can access
individual elements normally. Alternatively, you can index using parenthesis, e.g. `v(1)` indexes the first element.
Parenthesis are used instead of `getindex`/square brackets because `AbstractSIMD` objects represent a single number, and
for `x::Number`, `x[1] === x`.

If you wish to perform a reduction on the collection, the naming convention is prepending the base function with a `v`. These functions
are not overloaded, because for `x::Number`, `sum(x) === x`. Functions include `vsum`, `vprod`, `vmaximum`, `vminimum`, `vany`, and `vall`.

If you wish to define a new operation applied to the entire vector, do not define it in terms of operations on the individual eleemnts.
This will often lead to bad code generation -- bad in terms of both performance, and often silently producing incorrect results!
Instead, implement them in terms of existing functions defined on `::AbstractSIMD`. Please feel free to file an issue if you would like
clarification, and especially if you think the function may be useful for others and should be included in `VectorizationBase.jl`.
"""))
end

"""
  pause()

For use in spin-and-wait loops, like spinlocks.
"""
@inline pause() = ccall(:jl_cpu_pause, Cvoid, ())

# notinthreadedregion() = iszero(ccall(:jl_in_threaded_region, Cint, ()))
# function assert_init_has_finished()
#     _init_has_finished[] || throw(ErrorException("bad stuff happened"))
#     return nothing
# end

include("static.jl")
include("cartesianvindex.jl")
include("topology.jl")
include("cpu_info.jl")
# include("cache_inclusivity.jl")
include("early_definitions.jl")
include("promotion.jl")
include("llvm_types.jl")
include("lazymul.jl")
include("strided_pointers/stridedpointers.jl")
# include("strided_pointers/bitpointers.jl")
include("strided_pointers/cartesian_indexing.jl")
include("strided_pointers/grouped_strided_pointers.jl")
include("strided_pointers/cse_stridemultiples.jl")
include("llvm_intrin/binary_ops.jl")
include("vector_width.jl")
include("ranges.jl")
include("llvm_intrin/conversion.jl")
include("llvm_intrin/masks.jl")
include("llvm_intrin/intrin_funcs.jl")
include("llvm_intrin/memory_addr.jl")
include("llvm_intrin/unary_ops.jl")
include("llvm_intrin/vbroadcast.jl")
(VERSION ≥ v"1.6.0-DEV.674") && include("llvm_intrin/verf.jl")
include("llvm_intrin/vector_ops.jl")
include("llvm_intrin/nonbroadcastingops.jl")
include("llvm_intrin/integer_fma.jl")
include("vecunroll/memory.jl")
include("vecunroll/fmap.jl")
include("base_defs.jl")
include("alignment.jl")
include("special/misc.jl")
include("special/double.jl")
include("special/exp.jl")
# include("special/log.jl")

demoteint(::Type{T}, W) where {T} = False()
demoteint(::Type{UInt64}, W::StaticInt) = gt(W, pick_vector_width(UInt64))
demoteint(::Type{Int64}, W::StaticInt) = gt(W, pick_vector_width(Int64))

@generated function simd_vec(::DemoteInt, y::_T, x::Vararg{_T,_W}) where {DemoteInt<:StaticBool,_T,_W}
    W = 1 + _W
    T = DemoteInt === True ? _demoteint(_T) : _T
    trunc = T !== _T
    Wfull = nextpow2(W)
    ty = LLVM_TYPES[T]
    init = W == Wfull ? "undef" : "zeroinitializer"
    instrs = ["%v0 = insertelement <$Wfull x $ty> $init, $ty %0, i32 0"]
    Tup = Expr(:curly, :Tuple, T)
    for w ∈ 1:_W
        push!(instrs, "%v$w = insertelement <$Wfull x $ty> %v$(w-1), $ty %$w, i32 $w")
        push!(Tup.args, T)
    end
    push!(instrs, "ret <$Wfull x $ty> %v$_W")
    llvmc = :($LLVMCALL($(join(instrs,"\n")), _Vec{$Wfull,$T}, $Tup))
    trunc ? push!(llvmc.args, :(y % $T)) : push!(llvmc.args, :y)
    for w ∈ 1:_W
        ref = Expr(:ref, :x, w)
        trunc && (ref = Expr(:call, :%, ref, T))
        push!(llvmc.args, ref)
    end
    quote
        $(Expr(:meta,:inline))
        Vec($llvmc)
    end
end
function vec_quote(demote, W, Wpow2, offset::Int = 0)
    call = Expr(:call, :simd_vec, Expr(:call, demote ? :True : :False)); Wpow2 += offset
    iszero(offset) && push!(call.args, :y)
    foreach(w -> push!(call.args, Expr(:ref, :x, w)), max(1,offset):min(W,Wpow2)-1)
    # foreach(w -> push!(call.args, Expr(:call, :VecElement, Expr(:call, :zero, :T))), W+1:Wpow2)
    call
end
@generated function _vec(::StaticInt{_Wpow2}, ::DemoteInt, y::T, x::Vararg{T,_W}) where {DemoteInt<:StaticBool, _Wpow2, _W, T <: NativeTypes}
    W = _W + 1
    demote = DemoteInt === True
    Wpow2 = demote ? 2_Wpow2 : _Wpow2
    if W ≤ Wpow2
        vec_quote(demote, W, Wpow2)
    else
        tup = Expr(:tuple)
        offset = 0
        while offset < W
            push!(tup.args, vec_quote(demote, W, Wpow2, offset)); offset += Wpow2
        end
        Expr(:call, :VecUnroll, tup)
    end
end
@inline function Vec(y::T, x::Vararg{T,_W}) where {_W, T <: NativeTypes}
    W = StaticInt{_W}() + One()
    _vec(pick_vector_width(W, T), demoteint(T, W), y, x...)
end

@inline reduce_to_onevec(f::F, vu::VecUnroll) where {F} = ArrayInterface.reduce_tup(f, data(vu))

include("precompile.jl")
_precompile_()

function __init__()
    ccall(:jl_generating_output, Cint, ()) == 1 && return
    reset_features!()
    if unwrap(cpu_name()) !== Symbol(Sys.CPU_NAME::String)
        @info "Defining CPU name."
        define_cpu_name()
    end
    safe_topology_load!()
    redefine_attr_count()
    foreach(redefine_cache, 1:4)
    return nothing
end

end # module


text:
# This file was generated, do not modify it. # hide
using MLJ, PyPlot
import Statistics

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
function advance(x::DiseaseState, ::Type{S}) where {S <: DiseaseStateSequence}
  states = convert(DiseaseStates, S)
  current_index = findfirst(Ref(x) .== states)
  return states[current_index + 1]
end

function _count_by_state(events::Events{T},
                         state::DiseaseState,
                         time) where T <: DiseaseStateSequence
  n_ids = 0
  if state ∉ convert(DiseaseStates, T)
    error("Invalid state specified")
  end
  if state == convert(DiseaseStates, T)[1] # S
    nextstate = advance(state, T) # Either E or I
    @simd for i = 1:individuals(events)
      n_ids += (events[nextstate][i] > time) || (events[nextstate][i] === NaN)
    end
  elseif state in convert(DiseaseStates, T)[2:end-1]
    nextstate = advance(state, T) # Either I or R
    # E/I at or before time and I/R after time, or I/R never
    @simd for i = 1:individuals(events)
      n_ids += (events[state][i] <= time) && ((events[nextstate][i] > time) || (events[nextstate][i] === NaN))
    end
  elseif state == convert(DiseaseStates, T)[end] # I or R
    @simd for i = 1:individuals(events)
      n_ids += events[state][i] <= time # I/R at or before time
    end
  end
  @debug "$n_ids individual(s) in state $state at t = $time"
  return n_ids
end

function _state_color(x::DiseaseState)
  return [:purple, :lightblue4, :lightgreen, :yellow][findfirst(Ref(x) .== Pathogen._DiseaseStateVector)]
end

function _epidemic_curve(events::Events{T},
                         state::DiseaseState,
                         min,
                         max) where T <: DiseaseStateSequence
  if min >= max
    @error "Minimum time must be less than maximum time"
  end
  local times
  if state == convert(DiseaseStates, T)[1]
    nextstate = convert(DiseaseStates, T)[2]
    times = events[nextstate]
  elseif state in convert(DiseaseStates, T)[2:end-1]
    nextstate = advance(state, T)
    times = events[[state; nextstate]][:]
  elseif state == convert(DiseaseStates, T)[end]
    times = events[state]
  else
    @error "Invalid state specified"
  end
  times = times[Ref(min) .< times .< Ref(max)]
  sort!(times)
  insert!(times, 1, min)
  push!(times, max)
  counts = [_count_by_state(events, state, t) for t in times]
  return times, counts
end

function _epidemic_curve(events::Events{T},
                         state::DiseaseState,
                         times) where T <: DiseaseStateSequence
  counts = Array{Int64, 1}(undef, length(times))
  @simd for t in eachindex(times)
    counts[t] = _count_by_state(events, state, t)
  end
  return counts
end

@recipe function f(events::Events{T},
                   state::DiseaseState,
                   min, max) where {
                   T <: DiseaseStateSequence}
  xguide --> "Time"
  yguide --> "N"
  label --> convert(Char, state)
  linecolor --> _state_color(state)
  seriestype --> :steppre
  _epidemic_curve(events, state, min, max)
end

@recipe function f(events::Events{T},
                   state::DiseaseState) where T <: DiseaseStateSequence
  events, state, 0.0, maximum(events)
end

@recipe function f(
  events::Events{S},
  min,
  max) where {
  S <: DiseaseStateSequence}

  for s in convert(Tuple, S)
    @series begin
      events, s, min, max
    end
  end
end

@recipe function f(events::Events{T}) where T <: DiseaseStateSequence
  events, 0.0, maximum(events)
end

@recipe function f(events::Events{T},
                   state::DiseaseState,
                   times) where {
                   T <: DiseaseStateSequence}
  xguide --> "Time"
  yguide --> "N"
  label --> convert(Char, state)
  linecolor --> _state_color(state)
  seriestype --> :steppre
  _epidemic_curve(events, state, times)
end

@recipe function f(
  events::Events{S},
  times) where {
  S <: DiseaseStateSequence}
  for s in convert(Tuple, S)
    @series begin
      events, s, times
    end
  end
end


text:
# power analysis functions for proportions

"""
    powerPTest(h::Real = 0, n::Real = 0, alpha::Float64 = 0.05, alternative = "two.sided")

Compute power of a sample of `n` observations to test the effect size `h`
for one sample proportion (P) compared to a constant (c) at type I error = `alpha` (default: 0.05).
The effect size `h` = ϕ₁ - ϕ₂ where ϕ₁ = 2*asin(sqrt(P)) and ϕ₂ = 2*asin(sqrt(c)). The option
`alternative` is `two.sided` for H₀: ϕ₁ = ϕ₂ vs H₁: ϕ₁ ≠ ϕ₂, `less` for H₁: ϕ₁ < ϕ₂, and `greater`
for H₁: ϕ₁ > ϕ₂.
"""
function powerPTest(;
    h::Real = 0,
    n::Real = 0,
    alpha::Float64 = 0.05,
    alternative::String = "two.sided"
    )

    check_args(h=h,n=n,alpha=alpha)

    if alternative == "less"
        tside = 1
    elseif alternative in ("two","two.sided","two-sided","two sided")
        tside = 2
        h = abs(h)
    elseif alternative == "greater"
        tside = 3
    else
        error("`alternative` must be `less`, `two.sided`, or `greater`.")
    end

    if tside == 1
        return cdf(Normal(),quantile(Normal(),alpha) - h*sqrt(n))
    elseif tside == 2
        return ccdf(Normal(),cquantile(Normal(),alpha/2) - h*sqrt(n)) + cdf(Normal(),quantile(Normal(),alpha/2) - h*sqrt(n))
    elseif tside == 3
        return ccdf(Normal(),cquantile(Normal(),alpha) - h*sqrt(n))
    end
    error("internal error")
end

"""
    samplesizePTest(h::Real = 0, alpha::Float64 = 0.05, power::Float64 = 0.8, alternative = "two.sided")

Compute the sample size to test the effect size `h` for one sample proportion (P) compared to a constant (c)
with power > `power` (default: 0.8) and type I error = `alpha` (default: 0.05).
The effect size `h` = ϕ₁ - ϕ₂ where ϕ₁ = 2*asin(sqrt(P)) and ϕ₂ = 2*asin(sqrt(c)). The option
`alternative` is `two.sided` for H₀: ϕ₁ = ϕ₂ vs H₁: ϕ₁ ≠ ϕ₂, `less` for H₁: ϕ₁ < ϕ₂, and `greater`
for H₁: ϕ₁ > ϕ₂.
"""
function samplesizePTest(;
    h::Real = 0,
    alpha::Float64 = 0.05,
    power::Float64 = 0.8,
    alternative::String = "two"
    )

    check_args(h=h,alpha=alpha,power=power)

    return ceil(Int64,fzero(x->powerPTest(h = h, n = x, alpha = alpha, alternative = alternative) - power, 2 + 1e-10, 1e+09))
end

"""
    effectsizePTest(n::Real = 0, alpha::Float64 = 0.05, power::Float64 = 0.8, alternative = "two.sided")

Compute the effect size that a sample size `n` can detect for one sample proportion (P) compared to a constant (c)
with power > `power` (default: 0.8) and type I error = `alpha` (default: 0.05).
The effect size `h` = ϕ₁ - ϕ₂ where ϕ₁ = 2*asin(sqrt(P)) and ϕ₂ = 2*asin(sqrt(c)). The option
`alternative` is `two.sided` for H₀: ϕ₁ = ϕ₂ vs H₁: ϕ₁ ≠ ϕ₂, `less` for H₁: ϕ₁ < ϕ₂, and `greater`
for H₁: ϕ₁ > ϕ₂.
"""
function effectsizePTest(;
    n::Real = 0,
    alpha::Real = 0.05,
    power::Real = 0.0,
    alternative::String = "two"
    )

    check_args(n=n,alpha=alpha,power=power)

    return fzero(x->powerPTest(h = x, n = n, alpha = alpha, alternative = alternative) - power, 1e-10, 1 - 1e-10)
end

"""
    alphaPTest(h::Real = 0, n::Real = 0, power::Float64 = 0.8, alternative = "two.sided")

Compute the probability of type I error that a sample size `n` has in detecting the effect size `h`
for one sample proportion (P) compared to a constant (c) with power > `power` (default: 0.8)
and type I error = `alpha` (default: 0.05). The effect size `h` = ϕ₁ - ϕ₂ where ϕ₁ = 2*asin(sqrt(P)) and ϕ₂ = 2*asin(sqrt(c)).
The option `alternative` is `two.sided` for H₀: ϕ₁ = ϕ₂ vs H₁: ϕ₁ ≠ ϕ₂, `less` for H₁: ϕ₁ < ϕ₂, and `greater`
for H₁: ϕ₁ > ϕ₂.
"""
function alphaPTest(;
    h::Real = 0,
    n::Real = 0,
    power::Float64 = 0.8,
    alternative::String = "two"
    )

    check_args(h=h,n=n,power=power)

    return fzero(x->powerPTest(h = h, n = n, alpha = x, alternative = alternative) - power, 1e-10, 1 - 1e-10)
end

"""
    pwr.PTest(h::Real = 0, n::Real = 0, alpha::Float64 = 0.0, power::Float64 = 0.0, alternative = "two.sided")

Compute one of the test parameters such as sample size (`n`),
effect size (`h`), type I error (`alpha`), and power (`power`).
The parameter to be estimated must be set to zero (default).
The effect size `h` = ϕ₁ - ϕ₂ where ϕ₁ = 2*asin(sqrt(P)) and ϕ₂ = 2*asin(sqrt(c)).
The option `alternative` is `two.sided` for H₀: ϕ₁ = ϕ₂ vs H₁: ϕ₁ ≠ ϕ₂, `less` for H₁: ϕ₁ < ϕ₂, and `greater`
for H₁: ϕ₁ > ϕ₂.
"""
function PTest(;
    h::Real = 0,
    n::Real = 0,
    alpha::Float64 = 0.0,
    power::Float64 = 0.0,
    alternative::String = "two.sided"
    )
    if sum([x == 0 for x in (h,n,alpha,power)]) != 1
        error("exactly one of `h`, `n`, `power`, and `alpha` must be zero")
    end

    if power == 0.0
        power = powerPTest(h = h, n = n, alpha = alpha, alternative = alternative)
    elseif alpha == 0.0
        alpha = alphaPTest(h = h, n = n, power = power, alternative = alternative)
    elseif h == 0
        h = effectsizePTest(n = n, alpha = alpha, power = power, alternative = alternative)
    elseif n == 0
        n = samplesizePTest(h = h, alpha = alpha, power = power, alternative = alternative)
    end

    alt = Dict("two" => "two-sided", "two.sided" => "two-sided","two sided" => "two-sided", "less" => "less", "greater" => "greater")

    return htest(
        string("Proportion power calculation for binomial distribution (arcsine transformation)"),
        OrderedDict(
            "h" => h,
            "n" => n,
            "alpha" => alpha,
            "power" => power,
            "alternative" => alt[alternative]
            )
        )
end

max_stars_repo_name: mwsohn/pwr

max_issues_repo_name: mwsohn/pwr

max_forks_repo_name: mwsohn/pwr


text:
function PlanetConstants{FT}(::Type{EarthConstants}; kwargs...) where {FT}

  defaults = (
      k_Karman =          FT(0.4)                             , # Von Karman constant (1)
      molmass_dryair =    FT(28.97e-3)                        , # Molecular weight dry air (kg/mol)
      kappa_d =           FT(2//7)                            , # Adiabatic exponent dry air
      ρ_liq =             FT(1e3)                             , # Density of water (kg/m^3)
      ρ_ice =             FT(916.7)                           , # Density of ice (kg/m^3)
      molmass_water =     FT(18.01528e-3)                     , # Molecular weight (kg/mol)
      cp_v =              FT(1859)                            , # Isobaric specific heat vapor (J/kg/K)
      cp_l =              FT(4181)                            , # Isobaric specific heat liquid (J/kg/K)
      cp_i =              FT(2100)                            , # Isobaric specific heat ice (J/kg/K)
      T_freeze =          FT(273.15)                          , # Freezing point temperature (K)
      T_min =             FT(150.0)                           , # Minimum temperature guess in saturation adjustment (K)
      T_max =             FT(1000.0)                          , # Maximum temperature guess in saturation adjustment (K)
      T_icenuc =          FT(233.00)                          , # Homogeneous nucleation temperature (K)
      T_triple =          FT(273.16)                          , # Triple point temperature (K)
      LH_v0 =             FT(2.5008e6)                        , # Latent heat vaporization at T_0 (J/kg)
      LH_s0 =             FT(2.8344e6)                        , # Latent heat sublimation at T_0 (J/kg)
      press_triple =      FT(611.657)                         , # Triple point vapor pressure (Pa)
      ρ_ocean =           FT(1.035e3)                         , # Reference density sea water (kg/m^3)
      cp_ocean =          FT(3989.25)                         , # Specific heat sea water (J/kg/K)
      planet_radius =     FT(6.371e6)                         , # Mean planetary radius (m)
      day =               FT(86400)                           , # Length of day (s)
      Omega =             FT(7.2921159e-5)                    , # Ang. velocity planetary rotation (1/s)
      grav =              FT(9.81)                            , # Gravitational acceleration (m/s^2)
      orbit_semimaj =     FT(1*astro_unit(FT))                , # Length of semimajor orbital axis (m)
      TSI =               FT(1362)                            , # Total solar irradiance (W/m^2)
      MSLP =              FT(1.01325e5)                       , # Mean sea level pressure (Pa)
  )
  return PlanetConstants{FT}(;pairs(merge(defaults, kwargs.data))...)

end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
#
# PRIM-related types and functions
#

struct PrimBox <: Wrapper
    pyo::PyObject

    function PrimBox(pyo::PyObject)
        new(pyo)
    end
end

function Base.show(io::IO, box::PrimBox)
    print(io, box.pyo.__str__())
end

function show_tradeoff(box::PrimBox)
    fig = pycall(box.pyo.show_tradeoff, PyObject)
    return fig
end

function show_details(box::PrimBox, fig=nothing)
    pycall(box.pyo[:show_details], PyObject, fig=fig)
    nothing
end

stats(box::PrimBox) = box.pyo[:stats]

limits(box::PrimBox) = box.pyo[:limits]

Base.length(box::PrimBox) = box.pyo[:__len__]()

# From original python doc:
# x : a matrix-like object (pandas.DataFrame, numpy.recarray, etc.)
# the independent variables [can be anything convertible to a DataFrame]
#
# y : a list-like object, the column name (str), or callable
# the dependent variable either provided as a list-like object
# classifying the data into cases of interest (e.g., False/True),
# a list-like object storing the raw variable value (in which case
# a threshold must be given), a string identifying the dependent
# variable in x, or a function called on each row of x to compute the
# dependent variable
struct Prim <: Wrapper
    pyo::PyObject

    function Prim(x::DataSet, y::Vector;
                  threshold=nothing, threshold_type=">",
                  obj_func=prim.lenient1,
                  peel_alpha=0.05, paste_alpha=0.05, mass_min=0.05,
                  include=nothing, exclude=nothing, coi=nothing)

        pandasDF = pandas_dataframe(x; include=include, exclude=exclude)

        # Convert y into Vector{Bool} by matching category of interest
        # Note that classification and coi can be strings or symbols,
        # as long as they're consistent (i.e., 'in' and '==' work.)
        if coi != nothing
            if coi isa AbstractArray
                y = [value in coi for value in y]
            else
                y = (y .== coi)
            end
        end

        prim = rhodium.Prim(pandasDF, y;
                            threshold=threshold, threshold_type=threshold_type,
                            obj_func=obj_func, peel_alpha=peel_alpha, paste_alpha=paste_alpha,
                            mass_min=mass_min) # include=include, exclude=exclude, coi=coi)
        return new(prim)
    end
end

function find_box(p::Prim)
    pyo = pycall(p.pyo.find_box, PyObject)
    return PrimBox(pyo)
end

function find_all(p::Prim)
    boxes = pycall(p.pyo[:find_all], PyVector)
    return boxes
end

# Python code says this needs to be tested...
"""
    perform_pca(p::Prim, subsets::Dict=nothing, exclude=nothing)

Pre-process the data by performing a pca based rotation on it.

Arguments:
    subsets: optional dict with group name as key and a list of
             uncertainty names as values. If this is used, a
             constrained PCA-PRIM is executed
             (N.B. the list of uncertainties should not contain
             categorical uncertainties.
    exclude: optional list of str, the uncertainties that should
             be excluded from the rotation
"""
function perform_pca(p::Prim, subsets::Dict=nothing, exclude=nothing)
    p.pyo[:perform_pca](;subsets=subsets, exclude=exclude)
end


text:
"""
    absorb_fromleft(A::Array{T, 2}, B::Array{T, 3}) where T<:Number

Absorb the matrix A in B where A is at the left of B.

 -- A -- B -- = -- C --
         |         |

"""
function absorb_fromleft(A::Array{T, 2}, B::Array{T, 3}) where T<:Number
    a1, a2 = size(A)
    b1, b2, b3 = size(B)
    B = reshape(B, b1, b2*b3)
    return reshape(A*B, a1, b2, b3)
end

"""
    absorb_fromright(A::Array{T, 3}, B::Array{T, 2}) where T<:Number

Absorb the matrix B in A where B is at the right of A.

 -- A -- B -- = -- C --
    |              |

"""
function absorb_fromright(A::Array{T, 3}, B::Array{T, 2}) where T<:Number
    a1, a2, a3 = size(A)
    b1, b2 = size(B)
    A = reshape(A, a1*a2, a3)
    return reshape(A*B, a1, a2, b2)
end

"""
    prop_right2(L::Array{T, 2}, A::Array{T, 3}, B::Array{T, 3}) where T<:Number

Propagate the tensor L through the tensors A, conj(B).

 ---- A ---     ---
 |    |         |
 L    |     =  L2
 |    |         |
 ---- B ---     ---

"""
function prop_right2(L::Array{T, 2}, A::Array{T, 3},
                     B::Array{T, 3}) where T<:Number
    Ar = reshape(A, size(A, 1), size(A, 2)*size(A, 3))
    T1 = BLAS.gemm('T', 'N', L, Ar)
    T1 = reshape(T1, size(T1, 1)*size(B, 2), size(A, 3))
    Br = reshape(B, size(B, 1)*size(B, 2), size(B, 3))
    L2 = BLAS.gemm('T', 'N', T1, conj(Br))
    return L2
end

"""
    prop_right3(L::Array{T, 3}, A::Array{T, 3}, M::Array{T, 4},
                B::Array{T, 3}) where T<:Number

Propagate the tensor L through the tensors A, M, conj(B).

 ---- A ---         ---
 |    |             |
 L--- M ---  =  new_L
 |    |             |
 ---- B ---         ---

"""
function prop_right3(L::Array{T, 3}, A::Array{T, 3}, M::Array{T, 4},
                     B::Array{T, 3}) where T<:Number
    @tensoropt L1[l1, l2, s2, c] := L[l1, l2, l3]*conj(B[l3, s2, c])
    @tensoropt L2[l1, s1, b, c] := L1[l1, l2, s2, c]*M[l2, s1, s2, b]
    @tensoropt new_L[a, b, c] := L2[l1, s1, b, c]*A[l1, s1, a]
    return new_L
end

"""
    prop_right4(L::Array{T, 4}, A::Array{T, 3}, M1::Array{T, 4},
                M2::Array{T, 4}, B::Array{T, 3}) where T<:Number

Propagate the tensor L through the tensors A, M1, M2, conj(B).

 ---- A ---         ---
 |    |             |
 |--- M1---         |--
 L    |      =  new_L
 |----M2---         |--
 |    |             |
 ---- B ---         ---

"""
function prop_right4(L::Array{T, 4}, A::Array{T, 3}, M1::Array{T, 4},
                     M2::Array{T, 4}, B::Array{T, 3}) where T<:Number
    @tensoropt L1[l1, l2, l3, s3, d] := L[l1, l2, l3, l4]*conj(B[l4, s3, d])
    @tensoropt L2[l1, l2, s2, c, d] := L1[l1, l2, l3, s3, d]*M2[l3, s2, s3, c]
    @tensoropt L3[l1, s1, b, c, d] := L2[l1, l2, s2, c, d]*M1[l2, s1, s2, b]
    @tensoropt new_L[a, b, c, d] := L3[l1, s1, b, c, d]*A[l1, s1, a]
    return new_L
end

"""
    prop_left2(A::Array{T, 3}, B::Array{T, 3}, R::Array{T, 2}) where T<:Number

Propagate the tensor R through the tensors A, conj(B).

 --- A ----      ---
     |    |        |
     |    R   =   R2
     |    |        |
 --- B ----      ---

"""
function prop_left2(A::Array{T, 3}, B::Array{T, 3},
                    R::Array{T, 2}) where T<:Number
    Ar = reshape(A, size(A, 1)*size(A, 2), size(A, 3))
    T1 = BLAS.gemm('N', 'N', Ar, R)
    T1 = reshape(T1, size(A, 1), size(A, 2)*size(T1, 2))
    Br = reshape(B, size(B, 1), size(B, 2)*size(B, 3))
    L2 = BLAS.gemm('N', 'C', T1, Br)
    return L2
end

"""
    prop_left3(A::Array{T, 3}, M::Array{T, 4}, B::Array{T, 3},
               R::Array{T, 3}) where T<:Number

Propagate the tensor R through the tensors A, M, conj(B).

 --- A ----      ---
     |    |        |
 --- M -- R   =    new_R
     |    |        |
 --- B ----      ---

"""
function prop_left3(A::Array{T, 3}, M::Array{T, 4}, B::Array{T, 3},
                    R::Array{T, 3}) where T<:Number
    @tensoropt new_R[a, b, c] := (
        A[a, s1, r1]*M[b, s1, s2, r2]*conj(B[c, s2, r3])*R[r1, r2, r3]
    )
    return new_R
end

"""
    prop_left4(A::Array{T, 3}, M1::Array{T, 4}, M2::Array{T, 4},
               B::Array{T, 3}, R::Array{T, 4}) where T<:Number

Propagate the tensor R through the tensors A, M1, M2, conj(B).

 --- A -----     ----
     |     |        |
 --- M1 ---|     ---|
     |     R   =    new_R
 --- M2 ---|     ---|
     |     |        |
 --- B -----     ----

"""
function prop_left4(A::Array{T, 3}, M1::Array{T, 4}, M2::Array{T, 4},
                    B::Array{T, 3}, R::Array{T, 4}) where T<:Number
    @tensoropt R1[a, s1, r2, r3, r4] := A[a, s1, r1]*R[r1, r2, r3, r4]
    @tensoropt R2[a, b, s2, r3, r4] := R1[a, s1, r2, r3, r4]*M1[b, s1, s2, r2]
    @tensoropt R3[a, b, c, s3, r4] := R2[a, b, s2, r3, r4]*M2[c, s2, s3, r3]
    @tensoropt new_R[a, b, c, d] := R3[a, b, c, s3, r4]*conj(B[d, s3, r4])
    return new_R
end

"""
    prop_right_subexp(L::Array{T, 3}, W::Array{T, 4}, M::Array{T, 3}) where T<:Number

Propagate the tensor L through the tensors W, conj(M) for the DMRG3S algorithm.

 ----               ---
 |                  |
 |    |             |----
 L -- W ---     =   P
 |    |             |----
 |    |             |
 ---- M ---         ---

"""
function prop_right_subexp(L::Array{T, 3}, W::Array{T, 4},
                           M::Array{T, 3}) where T<:Number
    @tensoropt P[l1, s1, r2, r3] := (
        L[l1, l2, l3]*conj(M[l3, s2, r3])*W[l2, s1, s2, r2]
        )
    return P
end

"""
    prop_left_subexp(W::Array{T, 4}, M::Array{T, 3},
                     R::Array{T, 3}) where T<:Number

Propagate the tensor R through the tensors W, conj(M) for the DMRG3S algorithm.

        ----     ----
           |        |
           |     ---|
     |     |   =    P
 --- W --- R     ---|
     |     |        |
     |     |        |
 --- B -----     ----

"""
function prop_left_subexp(W::Array{T, 4}, M::Array{T, 3},
                          R::Array{T, 3}) where T<:Number
    @tensoropt R1[r1, r2, s2, l3] := R[r1, r2, r3]*conj(M[l3, s2, r3])
    @tensoropt P[s1, r1, l2, l3] := R1[r1, r2, s2, l3]*W[l2, s1, s2, r2]
    return P
end

"""
    absorb_Le(Le::Array{T, 3}, M::Array{T, 2}) where T<:Number

Absorb M and conj(M) into the left environment Le.

 ----- M ------     ----
 |                  |
Le ------------ =  Le --
 |                  |
 --- conj(M) --     ----

"""
function absorb_Le(Le::Array{T, 3}, M::Matrix{T}) where T<:Number
    if size(M, 1) == size(M, 2)
        @tensoropt Le[l1, l2, l3] = Le[a, l2, b]*M[a, l1]*conj(M[b, l3])
    else
        @tensoropt Le[l1, l2, l3] := Le[a, l2, b]*M[a, l1]*conj(M[b, l3])
    end
    return Le
end

"""
    absorb_Re(Re::Array{T, 3}, M::Array{T, 2}) where T<:Number

Absorb M and conj(M) into the right environment Re.

 ---- M ------     ----
             |        |
 ---------- Re  =  -- Re
             |        |
 - conj(M) ---     ----

"""
function absorb_Re(Re::Array{T, 3}, M::Matrix{T}) where T<:Number
    if size(M, 1) == size(M, 2)
        @tensoropt Re[r1, r2, r3] = M[r1, a]*conj(M[r3, b])*Re[a, r2, b]
    else
        @tensoropt Re[r1, r2, r3] := M[r1, a]*conj(M[r3, b])*Re[a, r2, b]
    end
    return Re
end

"""
    build_local_hamiltonian(Le::Array{T, 3}, W::Array{T, 4},
                            Re::Array{T, 3}) where T<:Number

Build the local Hamiltonian with the left and right environments.

 ---      ----
 |     |     |     |
Le --- W --- Re =  Hi
 |     |     |     |
 ---      ----

"""
function build_local_hamiltonian(Le::Array{T, 3}, W::Array{T, 4},
                                 Re::Array{T, 3})  where T<:Number
    @tensoropt Hi[l1, s1, r1, l3, s2, r3] := (
        Le[l1, l2, l3]*W[l2, s1, s2, r2]*Re[r1, r2, r3]
        )
    Hi = reshape(Hi, (size(Le, 1)*size(W, 2)*size(Re, 1),
                      size(Le, 3)*size(W, 3)*size(Re, 3)))
    return Hi
end

"""
    build_local_hamiltonian(Le::Array{T, 3}, W::Array{T, 4}, Re::Array{T, 3},
                            cache::Cache{T}) where T<:Number

Build the local Hamiltonian with the left and right environments. Pass a
preallocated cache that can be used to store the final tensor.
"""
function build_local_hamiltonian(
    Le::Array{T, 3}, W::Array{T, 4}, Re::Array{T, 3}, cache::Cache{T}
    )  where T<:Number
    # Check if any cache element of appropriate size can be used.
    needed_space = (size(Le, 1)*size(W, 2)*size(Re, 1)
                    *size(Le, 3)*size(W, 3)*size(Re, 3))
    loc_in_cache = is_in_cache(cache, needed_space)
    if loc_in_cache > 0
        Hi = reshape(cache.elts[loc_in_cache],
                     (size(Le, 1), size(W, 2), size(Re, 1),
                      size(Le, 3), size(W, 3), size(Re, 3)))
        @tensoropt Hi[l1, s1, r1, l3, s2, r3] = (
            Le[l1, l2, l3]*W[l2, s1, s2, r2]*Re[r1, r2, r3]
            )
    else
        @tensoropt Hi[l1, s1, r1, l3, s2, r3] := (
            Le[l1, l2, l3]*W[l2, s1, s2, r2]*Re[r1, r2, r3]
            )
    end
    Hi = reshape(Hi, (size(Le, 1)*size(W, 2)*size(Re, 1),
                      size(Le, 3)*size(W, 3)*size(Re, 3)))
    update_cache!(cache, Hi)
    return Hi
end

"""
    build_local_hamiltonian_2(Le::Array{T, 3}, W1::Array{T, 4},
                              W2::Array{T, 4}, Re::Array{T, 3},
                              cache::Cache{T}) where T<:Number

Build the local Hamiltonian with the left and right environments for DMRG2.

 ---            ---
 |    |     |     |    |
Le -- W1 -- W2 -- Re = Hi
 |    |     |     |    |
 ---            ---

"""
function build_local_hamiltonian_2(Le::Array{T, 3}, W1::Array{T, 4},
                                   W2::Array{T, 4}, Re::Array{T, 3},
                                   cache::Cache{T}) where T<:Number
    @tensoropt Hi[l1, s1, s3, r1, l3, s2, s4, r3] := (Le[l1, l2, l3]
                                                      *W1[l2, s1, s2, a]
                                                      *W2[a, s3, s4, r2]
                                                      *Re[r1, r2, r3])
    Hi = reshape(Hi, (size(Le, 1)*size(W1, 2)*size(W2, 2)*size(Re, 1),
                      size(Le, 3)*size(W1, 3)*size(W2, 3)*size(Re, 3)))
    # Check if any cache element of appropriate size can be used.
    needed_space = (size(Le, 1)*size(W1, 2)*size(W2, 2)*size(Re, 1)
                    *size(Le, 3)*size(W1, 3)*size(W2, 3)*size(Re, 3))
    loc_in_cache = is_in_cache(cache, needed_space)
    if loc_in_cache > 0
        Hi = reshape(cache.elts[loc_in_cache],
                     ((size(Le, 1), size(W1, 2), size(W2, 2), size(Re, 1),
                      size(Le, 3),size(W1, 3),size(W2, 3),size(Re, 3))))
        @tensoropt Hi[l1, s1, s3, r1, l3, s2, s4, r3] = (
            Le[l1, l2, l3]*W1[l2, s1, s2, a]*W2[a, s3, s4, r2]*Re[r1, r2, r3]
            )
    else
        @tensoropt Hi[l1, s1, s3, r1, l3, s2, s4, r3] := (
            Le[l1, l2, l3]*W1[l2, s1, s2, a]*W2[a, s3, s4, r2]*Re[r1, r2, r3]
            )
    end
    Hi = reshape(Hi, ((size(Le, 1)*size(W1, 2)*size(W2, 2)*size(Re, 1),
                       size(Le, 3)*size(W1, 3)*size(W2, 3)*size(Re, 3))))
    update_cache!(cache, Hi)
    return Hi
end


text:
module Reachable

module ModuleA
export RchA
struct RchA end
end # ModuleA

module ModuleB
using Reachable.ModuleA
export rchb
rchb(::RchA) = "hello"
f(a) = 1
end # ModuleB

end # Reachable


text:
using Pkg

Pkg.add("SymPy")


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

mutable struct IoK8sApiStorageV1CSINodeSpec <: SwaggerModel
    drivers::Any # spec type: Union{ Nothing, Vector{IoK8sApiStorageV1CSINodeDriver} } # spec name: drivers

    function IoK8sApiStorageV1CSINodeSpec(;drivers=nothing)
        o = new()
        validate_property(IoK8sApiStorageV1CSINodeSpec, Symbol("drivers"), drivers)
        setfield!(o, Symbol("drivers"), drivers)
        o
    end
end # type IoK8sApiStorageV1CSINodeSpec

const _property_map_IoK8sApiStorageV1CSINodeSpec = Dict{Symbol,Symbol}(Symbol("drivers")=>Symbol("drivers"))
const _property_types_IoK8sApiStorageV1CSINodeSpec = Dict{Symbol,String}(Symbol("drivers")=>"Vector{IoK8sApiStorageV1CSINodeDriver}")
Base.propertynames(::Type{ IoK8sApiStorageV1CSINodeSpec }) = collect(keys(_property_map_IoK8sApiStorageV1CSINodeSpec))
Swagger.property_type(::Type{ IoK8sApiStorageV1CSINodeSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiStorageV1CSINodeSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiStorageV1CSINodeSpec }, property_name::Symbol) =  _property_map_IoK8sApiStorageV1CSINodeSpec[property_name]

function check_required(o::IoK8sApiStorageV1CSINodeSpec)
    (getproperty(o, Symbol("drivers")) === nothing) && (return false)
    true
end

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


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

# %%
_foo(expr) = expr
_foo(expr::Symbol) = :foo
function _foo(expr::Expr)
    if expr.head === :macrocall
        Expr(expr.head, expr.args[1], _foo.(expr.args[2:end])...)
    else
        Expr(expr.head, _foo.(expr.args)...)
    end
end

macro foo(expr)
    _foo(expr) |> esc
end

# %%
_bar(expr) = expr
_bar(expr::Symbol) = :bar
function _bar(expr::Expr)
    if expr.head === :macrocall
        Expr(expr.head, expr.args[1], _bar.(expr.args[2:end])...)
    else
        Expr(expr.head, _bar.(expr.args)...)
    end
end

macro bar(expr)
    _bar(expr) |> esc
end

# %%
@macroexpand @foo hoge(x)

# %%
@macroexpand @bar hoge(x)

# %%
@macroexpand @bar @foo hoge(x)

# %%
@macroexpand hoge(@foo @bar x)

# %%
foo(x) = "foo"

# %%
bar(x) = "bar"

# %%
@foo hoge(x)

# %%
@bar hoge(x)

# %%
@bar @foo hoge(x)

# %%
@foo @bar hoge(x)

# %%
@macroexpand1 @foo @bar hoge(x)

# %%
@macroexpand @foo @bar hoge(x)

# %%

max_stars_repo_path: 0020/order of the actions of macros.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0020/order of the actions of macros.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0020/order of the actions of macros.jl

max_forks_repo_name: genkuroki/public


text:
abstract type Model end

predict(model::Model, pref::Preference, user::Int, item::Int) = error("predict is not defined for $(typeof(model)).")

objective(model::Model, dataset::Dataset) = error("objective is not defined for $(typeof(model)).")

update!(model::Model, dataset::Dataset) = error("update! is not defined for $(typeof(model)).")


text:
using QuantumPropagators
using Test
using SafeTestsets

# Note: comment outer @testset to stop after first @safetestset failure
@time @testset verbose=true "QuantumPropagators" begin

    print("\n* Spec.Rad. (test_specrad.jl):")
    @time @safetestset "Spec.Rad." begin include("test_specrad.jl") end

    print("\n* Cheby (test_cheby.jl):")
    @time @safetestset "Cheby" begin include("test_cheby.jl") end

    print("\n* Newton (test_newton.jl):")
    @time @safetestset "Newton" begin include("test_newton.jl") end

    print("\n* Exp (test_expprop.jl):")
    @time @safetestset "Exp" begin include("test_expprop.jl") end

    print("\n* Propagate (test_propagate.jl):")
    @time @safetestset "Propagate" begin include("test_propagate.jl") end

    print("\n")

end
;


text:
using TimerOutputs

mutable struct Report
    reason::Union{ReportReasons.ReportReason, Nothing}
    egraph::EGraph
    iterations::Int
    to::TimerOutput
end

Report() = Report(nothing, EGraph(), 0, TimerOutput())

Report(g::EGraph) = Report(nothing, g, 0, TimerOutput())

# string representation of timedata
function Base.show(io::IO, x::Report)
    g = x.egraph
    println(io, "Equality Saturation Report")
    println(io, "=================")
    println(io, "\tStop Reason: $(x.reason)")
    println(io, "\tIterations: $(x.iterations)")
    # println(io, "\tRules applied: $(g.age)")
    println(io, "\tEGraph Size: $(g.numclasses) eclasses, $(length(g.memo)) nodes")
    print_timer(io, x.to)
end


text:
using Knet

"""
Smooth L1 Loss.
- As stated in the Fast R-CNN paper, it is calculated as:

    for |x| < beta --> 0.5 * x^2
    otherwise      --> |x| - 0.5

    where x = gt - pred
"""

function smooth_l1(x; beta=1)
    low_idx = findall(Array(value(x)) .< beta)
    high_idx = findall(Array(value(x)) .>= beta)
    loss_sum = sum(x[high_idx] .- 0.5)
    loss_sum += sum(0.5 .* (x[low_idx].^2))
    return loss_sum
end

max_stars_repo_name: barisbatuhan/retinaface

max_issues_repo_name: barisbatuhan/retinaface

max_forks_repo_name: barisbatuhan/retinaface


text:
using Conda, VersionParsing, Test

exe = Sys.iswindows() ? ".exe" : ""

Conda.update()

env = :test_conda_jl
rm(Conda.prefix(env); force=true, recursive=true)

@test Conda.exists("curl", env)
Conda.add("curl", env)

@testset "Install Python package" begin
    Conda.add("python", env)
    pythonpath = joinpath(Conda.python_dir(env), "python" * exe)
    @test isfile(pythonpath)

    cmd = Conda._set_conda_env(`$pythonpath -c "import zmq"`, env)
    @test_throws Exception run(cmd)
    Conda.add("pyzmq", env)
    run(cmd)
end

curlvers = Conda.version("curl",env)
@test curlvers >= v"5.0"
@test Conda.exists("curl==$curlvers", env)

curl_path = joinpath(Conda.bin_dir(env), "curl" * exe)
@test isfile(curl_path)

@test "curl" in Conda.search("cu*", env)

Conda.rm("curl", env)
@test !isfile(curl_path)

pythonpath = joinpath(Conda.PYTHONDIR, "python" * exe)
@test isfile(pythonpath)
pyversion = read(`$pythonpath -c "import sys; print(sys.version)"`, String)
@test pyversion[1:1] == Conda.MINICONDA_VERSION

condarc_path = joinpath(homedir(), ".condarc")
last_modified = mtime(condarc_path)

Conda.add_channel("foo", env)
@test Conda.channels(env) == ["foo", "defaults"]

# Testing that calling the function twice do not fail
Conda.add_channel("foo", env)

# Validate that only the Conda.jl RC file was modified
@test occursin("foo", read(Conda.conda_rc(env), String))
@test !isfile(condarc_path) || !occursin("foo", read(condarc_path, String))

Conda.rm_channel("foo", env)
@test mtime(condarc_path) == last_modified

@test Conda.channels(env) == ["defaults"]

# Add a package from a specific channel
Conda.add("requests", env; channel="conda-forge")

@testset "Batch install and uninstall" begin
    Conda.add(["affine", "ansi2html"], env)
    installed = Conda._installed_packages(env)
    @test "affine" ∈ installed
    @test "ansi2html" ∈ installed

    Conda.rm(["affine", "ansi2html"], env)
    installed = Conda._installed_packages(env)
    @test "affine" ∉ installed
    @test "ansi2html" ∉ installed
end

# Run conda clean
Conda.clean(; debug=true)

@testset "Exporting and creating environments" begin
    new_env = :test_conda_jl_2
    Conda.add("curl", env)
    Conda.export_list("conda-pkg.txt", env)

    # Create a new environment
    rm(Conda.prefix(new_env); force=true, recursive=true)
    Conda.import_list(
        IOBuffer(read("conda-pkg.txt")), new_env; channels=["defaults", "conda-forge"]
    )

    # Ensure that our new environment has our channels and package installed.
    Conda.channels(new_env) == ["defaults", "conda-forge"]
    installed = Conda._installed_packages(new_env)
    @test "curl" ∈ installed
    rm("conda-pkg.txt")
end


text:
export rayleigh_criterion

"""
    rayleigh_criterion(focal_length, diameter; λ=500e-9)

Calculates the resolution of a microscope according to Rayleigh.
`focal_length` is the focal length of the lens.
`λ` is the wavelength and `diameter` is the diameter of the aperture of the lens.
See [Wikipedia](https://en.wikipedia.org/wiki/Angular_resolution#The_Rayleigh_criterion).

 # Examples
```julia-repl
julia> rayleigh_criterion(100, 200)
3.05e-7

julia> rayleigh_criterion(100, 200, λ=100e-9)
6.099999999999999e-8
```
"""
function rayleigh_criterion(focal_length, diameter; λ=500e-9)
    return rayleigh_criterion(calc_NA(focal_length, diameter), λ=λ)
end

"""
    rayleigh_criterion(NA; λ=500e-9)

Calculates the resolution of a microscope according to Rayleigh.
`λ` is the wavelength and `NA` the numerical aperture of the system.
See [Wikipedia](https://en.wikipedia.org/wiki/Angular_resolution#The_Rayleigh_criterion).

 # Examples
```julia-repl
julia> rayleigh_criterion(1.0)
3.05e-7

julia> rayleigh_criterion(1.0, λ=100e-9)
6.099999999999999e-8
```
"""
function rayleigh_criterion(NA; λ=500e-9)
    return 0.61 * λ / NA
end

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics