text:
using Test
using SimplePosets
using LinearAlgebra
using SimpleGraphs

@test 1==1

P = Chain(5) + Chain(3)
@test height(P) == 5

P = StandardExample(4)
A = minimals(P)
@test length(A) == 4
B = maximals(P)
@test length(B) == 4
@test length(union(A,B))==8

P = PartitionLattice(5)
@test card(P) == 52
@test height(P) == 5
@test inv(P') == P

P = Divisors(2*3*5)
M = mobius_matrix(P)
Z = zeta_matrix(P)
n = card(P)
@test M*Z == Matrix{Int}(I,n,n)

P = BooleanLattice(3)
@test element_type(P) == String

P = RandomPoset(10,2)
@test card(P) == 10

P = Divisors(2*3*5)
@test elements(relabel(P)) == collect(1:8)
@test length(random_linear_extension(P)) == 8

d = random_average_height(P,20)
@test d[1]==0.0


text:
using DBGraphs
using Base.Test
using JuliaDB
using LightGraphs
# write your own tests here

path = Pkg.dir("DBGraphs", "test", "data")
# g = DBGraph(loadfiles(path, indexcols=["src","dst"]))
g = DBGraph(loadfiles(path))

@show g
@test g.src == :src
@test g.dst == :dst
# @test 2 ∈ neighbors(g, 1)

@test outdegree(g, 1) == 2
@test indegree(g, 1) == 0
@test indegree(g, 3) == 0
@test outdegree(g, 3) == 2

@test ne(bfs_tree(g, 1)) == 1
@test ne(bfs_tree(g, 3)) == 10
for e in edges(bfs_tree(g, 3))
  println(e)
end

ext: jl


text:
using CSV
using DataFrames
using SeasonalTrendLoess
using Test

airp = CSV.File("airpassengers.csv") |> DataFrame

@testset "Test errors" begin
    @test_throws AssertionError stl(rand(100),10; ns=5)
    @test_throws AssertionError stl(rand(100),10; ns=8)
end

@testset "Air passengers" begin
    # TODO
    # no proper tests yet
    s = stl(float(airp.passengers), 12; ns = 7, no = 3)
    s = stl(float(airp.passengers), 12; ns = 7, robust = true, spm = true)
end


text:
using PowerModelsACDC, PowerModels, Ipopt, Juniper, JuMP, Cbc, Gurobi, Cbc
casename = "case9"
bt = 100 ## constraint tightening setting 95, 90, 85, 80
    include("./test/data/tnep/PSCC/conv_spec.jl")
    file = "./test/data/tnep/PSCC/$casename.m"
    data1 = PowerModels.parse_file(file)
    PowerModelsACDC.process_additional_data!(data1)
    # PowerModelsACDCInv.process_additional_data!(data1)
    data = deepcopy(data1)

    #increase load and generation by 3 times (except 6 and 24 bus)
    for i = 1: length(keys(data["load"]))
        data["load"]["$i"]["pd"] = 3* data1["load"]["$i"]["pd"]
        data["load"]["$i"]["qd"] = 3* data1["load"]["$i"]["qd"]
    end
    for j = 1: length(keys(data["gen"]))
        data["gen"]["$j"]["pmax"] = 3*data1["gen"]["$j"]["pmax"]
        data["gen"]["$j"]["qmax"] = 3*data1["gen"]["$j"]["qmax"]
        data["gen"]["$j"]["qmin"] = 3*data1["gen"]["$j"]["qmin"]
    end
    #constraint tightening
    for (b,branch) in data["branch"]
        branch["rate_a"] = bt/100* branch["rate_a"]
    end
    for (b,branchdc) in data["branchdc_ne"]
        branchdc["rateA"] = bt/100* branchdc["rateA"]
    end

    ipopt = JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol" => 1e-4, "print_level" => 1)
    cbc_solver = JuMP.optimizer_with_attributes(Cbc.Optimizer)
    gurobi = JuMP.optimizer_with_attributes(Gurobi.Optimizer)
    juniper = JuMP.optimizer_with_attributes(Juniper.Optimizer, "mip_solver" => cbc_solver, "nl_solver" => ipopt)

    s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => false, "process_data_internally" => false)

    converter_parameters_rxb(data)
    # converter_cost(data)      # only to use with 6 and 24 bus systems

    data_pf1 = deepcopy(data)

##
resultDC = run_tnepopf(data, DCPPowerModel, gurobi, setting = s)
resultDC["built_cv"], resultDC["built_br"] = display_results(data_pf1, resultDC)

resultSOCBF = run_tnepopf_bf(data, SOCBFPowerModel, gurobi, setting = s)
resultSOCBF["built_cv"], resultSOCBF["built_br"] = display_results(data_pf1, resultSOCBF)

resultSOCWR = run_tnepopf(data, SOCWRPowerModel, gurobi, setting = s)
resultSOCWR["built_cv"], resultSOCWR["built_br"] = display_results(data_pf1, resultSOCWR)

resultQC     =  run_tnepopf(data, QCRMPowerModel, gurobi; setting = s)
resultQC["built_cv"], resultQC["built_br"] = display_results(data_pf1, resultQC)

resultAC_juni     =  run_tnepopf(data, ACPPowerModel, juniper; setting = s)
resultAC_juni["built_cv"], resultAC_juni["built_br"] = display_results(data_pf1, resultAC_juni)

resultLPAC     =  run_tnepopf(data, LPACCPowerModel, gurobi; setting = s)
resultLPAC["built_cv"], resultLPAC["built_br"] = display_results(data_pf1, resultLPAC)

function converter_cost(data)
    for(c,conv) in data["convdc_ne"]
        conv["cost"] = conv["Pacmax"]*0.083 *100+ 28
        display(conv["cost"])
    end
end

function converter_parameters_rxb(data)
for (c,conv) in data["convdc_ne"]

    bus = conv["busac_i"]
    display(bus)
    base_kV = data["bus"]["$bus"]["base_kv"]
    base_S = sqrt((100*conv["Pacmax"])^2+(100*conv["Qacmax"])^2) #base MVA = 100
    base_Z = base_kV^2/base_S # L-L votlage/3 phase power
    base_Y= 1/base_Z
    display("baseS:$base_S")
    conv["xtf"] = 0.10*100/base_S #new X =old X *(100MVA/old Sbase)
    conv["rtf"] = conv["xtf"]/100
    conv["bf"] = 0.08*base_S/100
    conv["xc"] = 0.07*100/base_S #new X =old X *(100MVA/old Zbase)
    conv["rc"] = conv["xc"]/100 #new X =old X *(100MVA/old Zbase)
    rtf = conv["rtf"]
    xtf = conv["xtf"]
    bf = conv["bf"]
    Pmax = conv["Pacmax"]
    Pmin =  conv["Pacmin"]
    Qmax = conv["Qacmax"]
    Qmin =  conv["Qacmin"]

    conv["Imax"] = sqrt(Pmax^2+Qmax^2)
    xc = conv["xc"]
    rc = conv["rc"]
    Imax = conv["Imax"]

    println("rtf:$rtf","     ", "xtf:$xtf","     ", "bf:$bf", "     ","rc:$rc", "     ","xc:$xc", "     ","Imax:$Imax","     ","Pmin: $Pmin", "     ","Pmax: $Pmax" )
    if xtf > 0.1 || xc > 0.1
    end
end
end


text:
module model
using Distributions
include("data_repeatrate.jl")
include("tree_number_estimate.jl")
include("kmean.jl")
include("sigma_estimate.jl")

function fit(df)
    x1 = 30
    x2 = 38
    y1 = 65
    y2 = 75

    one_tree,df_close = data_repeatrate.Onetree(df,x1,x2,y1,y2)
    #println(one_tree)

    fruit_number = data_repeatrate.Areafruit(one_tree)
    #println(fruit_number)

    df_c = df[df.Close.>0,:]

    tree_number = round(Int, tree_number_estimate.number_estimate(df_c, fruit_number))
    println(tree_number)
    df_k = [df_c.X df_c.Y]

    #df_k = DataFrames(X = df_c.X, Y = df_c.Y)
    #println(size(df_k))
    centers = kmean.K_mean(df_k', tree_number, df_c.Close .^(1/2))
    # println("tree centers")
    # println(centers)
    clustInd = kmean.k_mean_predict(df_k', centers)
    # println(clustInd)
    df_centers = centers[:, clustInd]
    loc = [df.X df.Y]
    esigma_x, esigma_y = sigma_estimate.sigma_est(loc, df_centers, df.Close')
    esigma = [esigma_x 0; 0 esigma_y]
    return fruit_number, centers, esigma
end

function predict(loc::AbstractMatrix{<:Real},
    fruit_number::Float64,
    centers::AbstractMatrix{<:Real},
    esigma::AbstractMatrix{<:Real})
    d, n = size(loc)
    d2, k = size(centers)
    (d == 2) || throw(ArgumentError(
    "dims of data $d should be 2"))
    (d == d2) || throw(ArgumentError(
    "dims of data $d and center $d2 should be equal"))
    pred_fruit_arr = Array{Float64, 2}(undef, n, k)

    for j in 1:k
        center_t = centers[:, j]
        distriNormal = MvNormal(center_t, esigma)
        for i in 1:n
            loc_t = loc[:, i]
            pred_fruit_arr[i, j] = fruit_number * pdf(distriNormal, loc_t)
        end
    end
    pred_fruit = sum(pred_fruit_arr, dims=2)[:,1]
    return pred_fruit

end # function

function mse(pred_fruit, gt_fruit)
    len1 = length(pred_fruit)
    len2 = length(gt_fruit)
    (len1 == len2) || throw(ArgumentError(
    "length of pred_fruit $len1 and gt_fruit $len2 should be equal"))
    return mean((pred_fruit - gt_fruit).^ 2)
end # function

function accuarcy_soft(pred_fruit, gt_fruit)
    len1 = length(pred_fruit)
    len2 = length(gt_fruit)
    (len1 == len2) || throw(ArgumentError(
    "length of pred_fruit $len1 and gt_fruit $len2 should be equal"))
    correct_soft = sum(abs.(pred_fruit - gt_fruit) .< 3)
    return correct_soft/len1
end

end  # module model


text:
module TestLift

using StructuredQueries
const SQ = StructuredQueries
using Base.Test

##############################################################################
##
## Test standard lifting semantics
##
##############################################################################

types = [
    Float16,
    Float32,
    Float64,
    Int128,
    Int16,
    Int32,
    Int64,
    Int8,
    UInt16,
    UInt32,
    UInt64,
    UInt8,
]

f(x::Number) = 5 * x
f(x::Number, y::Number) = x + y
f(x::Number, y::Number, z::Number) = x + y * z

for T in types
    a = one(T)
    x = Nullable{T}(a)
    y = Nullable{T}()

    U1 = Core.Inference.return_type(f, Tuple{T})
    @test isequal(SQ.lift(f, x), Nullable(f(a)))
    @test isequal(SQ.lift(f, y), Nullable{U1}())

    U2 = Core.Inference.return_type(f, Tuple{T, T})
    @test isequal(SQ.lift(f, x, x), Nullable(f(a, a)))
    @test isequal(SQ.lift(f, x, y), Nullable{U2}())

    U3 = Core.Inference.return_type(f, Tuple{T, T, T})
    @test isequal(SQ.lift(f, x, x, x), Nullable(f(a, a, a)))
    @test isequal(SQ.lift(f, x, y, x), Nullable{U3}())
end

##############################################################################
##
## Test non-standard lifting semantics
##
##############################################################################

# three-valued logic

# & truth table
v1 = SQ.lift(&, Nullable(true), Nullable(true))
v2 = SQ.lift(&, Nullable(true), Nullable(false))
v3 = SQ.lift(&, Nullable(true), Nullable{Bool}())
v4 = SQ.lift(&, Nullable(false), Nullable(true))
v5 = SQ.lift(&, Nullable(false), Nullable(false))
v6 = SQ.lift(&, Nullable(false), Nullable{Bool}())
v7 = SQ.lift(&, Nullable{Bool}(), Nullable(true))
v8 = SQ.lift(&, Nullable{Bool}(), Nullable(false))
v9 = SQ.lift(&, Nullable{Bool}(), Nullable{Bool}())

@test isequal(v1, Nullable(true))
@test isequal(v2, Nullable(false))
@test isequal(v3, Nullable{Bool}())
@test isequal(v4, Nullable(false))
@test isequal(v5, Nullable(false))
@test isequal(v6, Nullable(false))
@test isequal(v7, Nullable{Bool}())
@test isequal(v8, Nullable(false))
@test isequal(v9, Nullable{Bool}())

# | truth table
u1 = SQ.lift(|, Nullable(true), Nullable(true))
u2 = SQ.lift(|, Nullable(true), Nullable(false))
u3 = SQ.lift(|, Nullable(true), Nullable{Bool}())
u4 = SQ.lift(|, Nullable(false), Nullable(true))
u5 = SQ.lift(|, Nullable(false), Nullable(false))
u6 = SQ.lift(|, Nullable(false), Nullable{Bool}())
u7 = SQ.lift(|, Nullable{Bool}(), Nullable(true))
u8 = SQ.lift(|, Nullable{Bool}(), Nullable(false))
u9 = SQ.lift(|, Nullable{Bool}(), Nullable{Bool}())

@test isequal(u1, Nullable(true))
@test isequal(u2, Nullable(true))
@test isequal(u3, Nullable(true))
@test isequal(u4, Nullable(true))
@test isequal(u5, Nullable(false))
@test isequal(u6, Nullable{Bool}())
@test isequal(u7, Nullable(true))
@test isequal(u8, Nullable{Bool}())
@test isequal(u9, Nullable{Bool}())

# others

x1 = Nullable(1)
x2 = Nullable(2)
y = Nullable{Int}()
z1 = 1
z2 = 2

@test SQ.lift(isnull, x1) == false
@test SQ.lift(isnull, y) == true

@test SQ.lift(isless, x1, y) == true
@test SQ.lift(isless, y, x1) == false
@test SQ.lift(isless, x1, x2) == true
@test SQ.lift(isless, x2, x1) == false
@test SQ.lift(isless, y, y) == false
@test SQ.lift(isless, x1, z2) == true
@test SQ.lift(isless, x2, z1) == false
@test SQ.lift(isless, z1, x2) == true
@test SQ.lift(isless, z2, x1) == false

@test SQ.lift(get, x1) == 1
@test_throws NullException SQ.lift(get, y)
@test SQ.lift(get, y, 1) == 1

end


text:
"""
    struct SourceManager <: Any
Hold a pointer to a `clang::SourceManager` object.
"""
struct SourceManager
    ptr::CXSourceManager
end

Base.unsafe_convert(::Type{CXSourceManager}, x::SourceManager) = x.ptr
Base.cconvert(::Type{CXSourceManager}, x::SourceManager) = x


text:
function ties{T<:Real}(a::Array{T,1})
    unique(a[2:end][a[2:end] .== a[1:end-1]])
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    VarCoefDerivativeOperator

A dissipation operator on a nonperiodic finite difference grid.
"""
struct VarCoefDerivativeOperator{T,CoefficientCache,Parallel,SourceOfCoefficients,Grid} <: AbstractVariableCoefficientDerivativeOperator{T}
    coefficients::VarCoefDerivativeCoefficients{T,CoefficientCache,Parallel,SourceOfCoefficients}
    grid::Grid
    Δx::T
    factor::T
    b::Vector{T}

    function VarCoefDerivativeOperator(coefficients::VarCoefDerivativeCoefficients{T,CoefficientCache,Parallel,SourceOfCoefficients},
                                       grid::Grid, b::Vector{T}) where {T,CoefficientCache,Parallel,SourceOfCoefficients,Grid}
        @argcheck checkbounds(Bool, grid, coefficients.coefficient_cache) DimensionMismatch
        @argcheck length(grid) == length(b)

        Δx = step(grid)
        factor = inv(Δx^derivative_order(coefficients))
        new{T,CoefficientCache,Parallel,SourceOfCoefficients,Grid}(
            coefficients, grid, Δx, factor, b)
    end
end

function Base.show(io::IO, D::VarCoefDerivativeOperator{T}) where {T}
    print(io, "SBP variable coefficient ")
    if  derivative_order(D) == 1
        print(io, "1st")
    elseif  derivative_order(D) == 2
        print(io, "2nd")
    elseif  derivative_order(D) == 3
        print(io, "3rd")
    else
        print(io, derivative_order(D), "th")
    end
    print(io, "derivative operator with order of accuracy ")
    print(io, accuracy_order(D), " {T=", T, ", Parallel=", typeof(D.coefficients.parallel), "} \n")
    print(io, "on a grid in [", first(grid(D)), ", ", last(grid(D)),
                "] using ", length(grid(D)), " nodes \n")
    print(io, "and coefficients given in \n")
    print(io, source_of_coeffcients(D))
end

"""
    mul!(dest::AbstractVector, D::VarCoefDerivativeOperator, u::AbstractVector, α, β)

Compute `α*D*u + β*dest` and store the result in `dest`.
"""
Base.@propagate_inbounds function mul!(dest::AbstractVector, D::VarCoefDerivativeOperator,
                                        u::AbstractVector, α, β)
    @boundscheck begin
        @argcheck size(D, 2) == length(u) DimensionMismatch
        @argcheck size(D, 1) == length(dest) DimensionMismatch
    end
    @inbounds mul!(dest, D.coefficients, u, D.b, D.factor*α, β)
end

"""
    mul!(dest::AbstractVector, D::VarCoefDerivativeOperator, u::AbstractVector, α)

Compute `α*D*u` and store the result in `dest`.
"""
Base.@propagate_inbounds function mul!(dest::AbstractVector, D::VarCoefDerivativeOperator,
                                        u::AbstractVector, α)
    @boundscheck begin
        @argcheck size(D, 2) == length(u) DimensionMismatch
        @argcheck size(D, 1) == length(dest) DimensionMismatch
    end
    @inbounds mul!(dest, D.coefficients, u, D.b, D.factor*α)
end

"""
    var_coef_derivative_operator(source_of_coefficients, derivative_order, accuracy_order, xmin, xmax, N, left_weights, right_weights, bfunc, parallel=Val{:serial}())

Create a `VarCoefDerivativeOperator` approximating a `derivative_order`-th
derivative with variable coefficients `bfunc` on a grid between `xmin` and
`xmax` with `N` grid points up to order of accuracy `accuracy_order` with
coefficients given by `source_of_coefficients`.
The evaluation of the derivative can be parallised using threads by chosing
`parallel=Val{:threads}())`.
"""
function var_coef_derivative_operator(source_of_coefficients, derivative_order, accuracy_order,
                                        xmin, xmax, N, bfunc, parallel=Val{:serial}())
    grid = construct_grid(source_of_coefficients, accuracy_order, xmin, xmax, N)
    coefficients = var_coef_derivative_coefficients(source_of_coefficients, derivative_order, accuracy_order, grid, parallel)
    VarCoefDerivativeOperator(coefficients, grid, bfunc.(grid))
end

"""
    mass_matrix(D::Union{DerivativeOperator,VarCoefDerivativeOperator})

Create the diagonal mass matrix for the SBP derivative operator `D`.
"""
function mass_matrix(D::Union{DerivativeOperator,VarCoefDerivativeOperator})
    m = fill(one(eltype(D)), length(grid(D)))
    @unpack left_weights, right_weights = D.coefficients

    m[1:length(left_weights)] = left_weights
    m[end:-1:end-length(right_weights)+1] = right_weights
    Diagonal(D.Δx * m)
end


text:
abstract type MappingFunctions end

struct Sigmoid   <: MappingFunctions end
struct HardLimit <: MappingFunctions end
struct ReLU      <: MappingFunctions end
struct Cosine    <: MappingFunctions end
struct Qubit     <: MappingFunctions end
struct HyperbolicTangent <: MappingFunctions end

function mapping(W::AbstractMatrix{T}, b::AbstractMatrix{T}, X::AbstractMatrix{T}, mf::MappingFunctions) where {T<:BaseType}
    return map(x->mapfunc(x, mf), W*X .+ b)
end

function mapping(W::AbstractVector{T}, b::T, X::AbstractVector{T}, mf::MappingFunctions) where {T<:BaseType}
    return map(x-mapfunc(x, mf), dot(conj(W),X) + b)
end

function mapfunc(x::T, mf::Sigmoid) where {T<:BaseType}
    return one(T)/(one(T)+exp(-x))
end

function mapfunc(x::T, mf::HardLimit) where {T<:BaseType}
    return x >= zero(T) ? one(T) : zero(T)
end

function mapfunc(x::T, mf::ReLU) where {T<:BaseType}
    return x >= zero(T) ? x : zero(T)
end

function mapfunc(x::T, mf::Cosine) where {T<:BaseType}
    return cos(x)
end

function mapfunc(x::T, mf::Qubit) where {T<:BaseType}
    return sin(0.5*pi*x)^2
end

function mapfunc(x::T, mf::HyperbolicTangent) where {T<:BaseType}
    v = exp(-x)
    return (one(T)-v)/(one(T)+v)
end

# split-type functions
for mf in (:Sigmoid, :HardLimit, :ReLU, :Cosine, :Qubit, :HyperbolicTangent)
    @eval begin
        struct $(Symbol(mf,:Split)) <: MappingFunctions end
        export $(Symbol(mf,:Split))

        function mapfunc(x::T, mf::$(Symbol(mf,:Split))) where {T<:AbstractFloat}
            return mapfunc(x, $(mf)())
        end

        function mapfunc(x::Complex{T}, mf::$(Symbol(mf,:Split))) where {T<:AbstractFloat}
            return complex(mapfunc(real(x), $(mf)()),
                           mapfunc(imag(x), $(mf)()))
        end

        function mapfunc(x::Quaternion{T},  mf::$(Symbol(mf,:Split))) where {T<:AbstractFloat}
            return quaternion(mapfunc(real(x),  $(mf)()),
                              mapfunc(imagi(x), $(mf)()),
                              mapfunc(imagj(x), $(mf)()),
                              mapfunc(imagk(x), $(mf)()))
        end
    end
end


text:
"""
   lenardcoefficient(e0::Float64, ashell::AtomicSubShell)

Computes the Lenard coefficient according to the algorithm of Heinrich.
"Heinrich K. F. J. (1967) EPASA 2, paper no. 7"
"""
lenardcoefficient(e0::Float64, ashell::AtomicSubShell) = 4.5e5 / ((0.001 * e0)^1.65 - (0.001 * energy(ashell))^1.65)

"""
    ionizationdepthratio(primary::AtomicSubShell, secondary::AtomicSubShell, e0::Float64)

Ionization depth ratio from "Reed S.J.B. (1990) Microbeam Analysis, p.109"
"""
function ionizationdepthratio(primary::AtomicSubShell, secondary::AtomicSubShell, e0::Float64)
    uA, uB = e0 / energy(secondary), e0 / energy(primary)
    return ((((uB * log(uB)) - uB) + 1.0) / (((uA * log(uA)) - uA) + 1.0))
end

"""
    familyfactor(shellA::AtomicSubShell, shellB::AtomicSubShell)::Float64

Accounts for the differences in ionization cross section between K , L and M shells
"""
function familyfactor(shellA::AtomicSubShell, shellB::AtomicSubShell)::Float64
    fA, fB = shell(shellA), shell(shellB)
    if fA == fB
        res = 1.0 # fA==fB
    else
        if (n(fA) == 1) && (n(fB) == 2)
            res = 1.0 / 0.24
        elseif (n(fA) == 2) && (n(fB) == 1)
            res = 0.24
        elseif (n(fA) == 3) && ((n(fB) == 1) || (n(fB) == 2))
            res = 0.02
        else
            res = 0.0
        end
    end
    return res
end

"""
    ionizationfraction(shell::AtomicSubShell)

The fraction of the ionizations to attribute to the specified shell.  Computed
from the jump ratio.
"""
function ionizationfraction(shell::AtomicSubShell)
    r = jumpratio(shell)
    return r >= 1.0 ? (r - 1.0) / r : 0.0
end

"""
The `ReedInternal` structure optimizes the calculation of the Reed correction algorithm.
"""
struct ReedInternal
    primary::CharXRay
    kk::Float64
    lenard::Float64
    muB::Float64

    function ReedInternal(comp::Material, primary::CharXRay, secondary::AtomicSubShell, e0::Float64)
        bElm, aElm = element(primary), element(secondary)
        @assert aElm in keys(comp) "Element $(aElm.symbol) not in $(comp). (aElm)"
        @assert bElm in keys(comp) "Element $(bElm.symbol) not in $(comp). (bElm)"
        cB = comp[bElm]
        muB_A, muB = mac(aElm, primary), mac(comp, primary)
        ionizeF = ionizationfraction(secondary)
        fluorB = meanfluorescenceyield(element(primary), shell(inner(primary)))
        v = lenardcoefficient(e0, secondary) / muB # keV
        ss = ionizationdepthratio(inner(primary), secondary, e0)
        f = familyfactor(secondary, inner(primary))
        k = f * 0.5 * cB * (muB_A / muB) * ionizeF * fluorB * (a(aElm) / a(bElm)) * ss
        return new(primary, k, v, muB)
    end
end

Base.show(io::IO, ri::ReedInternal) = print(io, repr(ri.primary))

"""
The `ReedFluorescence` structure implements `FluorescenceCorrection` for the Reed fluorescence model.
"""
struct ReedFluorescence <: FluorescenceCorrection
    comp::Material
    secondary::AtomicSubShell # Emitter
    e0::Float64 # Incident beam energy
    # Cached partial calculations
    exciters::Vector{ReedInternal}
end

Base.show(io::IO, reed::ReedFluorescence) =
    print(io, "Reed[$(reed.secondary) due to $(reed.exciters) for $(name(reed.comp)) at $(reed.e0/1000.0) keV]")

"""
   F(reed::ReedFluorescence, secondary::CharXRay, toa::Float64)

Compute the enhancement of the secondary characteristic X-ray due to the
primaries specified in reed.
"""
function F(reed::ReedFluorescence, secondary::CharXRay, toa::Float64)
    function finternal(ri::ReedInternal, secondary::CharXRay, toa::Float64, comp::Material)
        @assert ri.kk >= 0.0 "ri.kk = $(ri.kk) for $(ri) and $(secondary) in $(comp)"
        u = mac(comp, secondary) / (sin(toa) * ri.muB)
        # TODO: Evaluate whether weight(ri.primary) is necessary/correct???
        return normweight(ri.primary) * ri.kk * ((log(1.0 + u) / u) + (log(1.0 + ri.lenard) / ri.lenard))
    end
    return isempty(reed.exciters) ? 1.0 :
           1.0 + mapreduce(ex -> finternal(ex, secondary, toa, reed.comp), +, reed.exciters)
end
"""
    fluorescencecorrection(::Type{ReedFluorescence}, comp::Material, primary::Vector{CharXRay}, secondary::AtomicSubShell, e0::Float64)

Construct an instance of a ReedFluorescence correction structure to compute the
secondary fluorescence due to a primary characteristic X-ray in the specified
material and beam energy.
"""
function fluorescencecorrection(
    ::Type{ReedFluorescence},
    comp::Material,
    primarys::Vector{CharXRay},
    secondary::AtomicSubShell,
    e0::Float64,
)
    ris = Vector{ReedInternal}()
    if element(secondary) in keys(comp)
        for primary in primarys
            if (energy(primary) >= energy(secondary)) && (element(primary) in keys(comp))
                push!(ris, ReedInternal(comp, primary, secondary, e0))
            end
        end
    end
    return ReedFluorescence(comp, secondary, e0, ris)
end

NeXLCore.minproperties(::Type{ReedFluorescence}) = (:BeamEnergy, :TakeOffAngle)


text:
module TestBasics
using ConfigEnv
using Test

const dir = dirname(@__FILE__)

# There is no "USER" variable on windows.
initial_value = haskey(ENV, "USER") ? ENV["USER"] : "WINDOWS"
ENV["USER"] = initial_value

@testset "Locality" begin
    cfg = dotenv()
    @test haskey(cfg.dict, "QWERTY")
    @test cfg["QWERTY"] == "ZXC"
end

@testset "basic" begin
    # basic input
    str = "BASIC=basic"
    file = joinpath(dir, ".env.override")
    file2 = joinpath(dir, ".testenv")

    # can turn off override of ENV vars
    cfg = dotenv(file, overwrite = false)
    @test ENV["USER"] != cfg["USER"]
    @test ENV["USER"] == initial_value

    # iobuffer, string, file
    @test ConfigEnv.parse(str) == Dict("BASIC"=>"basic")
    @test ConfigEnv.parse(read(file)) == Dict("CUSTOMVAL123"=>"yes","USER"=>"replaced value")
    @test dotenv(file).dict == Dict("CUSTOMVAL123"=>"yes","USER"=>"replaced value")

    @test isempty(dotenv("inexistentfile.env"))

    # length of returned values
    @test length(dotenv(file2).dict) == 10

    # appropiately loaded into ENV if CUSTOM_VAL is non existent
    @test ENV["CUSTOMVAL123"] == "yes"

    # Test that EnvProxyDict is reading from ENV
    ENV["SOME_RANDOM_KEY"] = "abc"
    cfg = dotenv(file)
    @test !haskey(cfg.dict, "SOME_RANDOM_KEY")
    @test cfg["SOME_RANDOM_KEY"] == "abc"
    @test get(cfg, "OTHER_RANDOM_KEY", "zxc") == "zxc"
end

@testset "Overwrite dotenv" begin
    # basic input
    file = joinpath(dir, ".env.override")

    # By default force override
    cfg = dotenv(file)
    @test ENV["USER"] == cfg["USER"]
    @test ENV["USER"] == "replaced value"

    # Restore previous environment
    ENV["USER"] = initial_value
end

@testset "No overwrite dotenvx" begin
    # basic input
    file = joinpath(dir, ".env1")

    # By default force overwrite
    ENV["FOO"] = "FUU"
    cfg = dotenvx(file)
    @test ENV["FOO"] == "FUU"
    @test cfg["FOO"] == "BAR"
end

@testset "Other means to load data" begin
    txt = "FOO = BAR"
    cfg = dotenv(txt)
    @test cfg.dict["FOO"] == "BAR"

    iob = IOBuffer()
    print(iob, "FOO = BAR")
    cfg = dotenv(iob)
    @test cfg.dict["FOO"] == "BAR"
end

end # module


text:
module TestCron

using ZulipReminderBot
using ZulipReminderBot: cron_worker, heartbeat!, TimedMessage, Message, curts
using Test

function mockmessage(content; sender_id = 1, stream = "Stream1", topic = "Topic1", type = "stream")
    return Message(stream, topic, type, sender_id, content)
end

@testset "Two incoming messages after current time" begin
    tmsg1 = TimedMessage(1, 1, 100, mockmessage("Ping1"))
    tmsg2 = TimedMessage(2, 2, 200, mockmessage("Ping2"))

    cin = Channel{Tuple{TimedMessage, Int}}(Inf)
    cout = Channel{TimedMessage}(Inf)
    put!(cin, (tmsg1, 1))
    put!(cin, (tmsg2, 1))

    sched = TimedMessage[]
    heartbeat!(sched, 10, cin, cout)
    @test !isready(cin)
    @test !isready(cout)

    @test length(sched) == 2
    @test sched[1].id == 2
    @test sched[2].id == 1
end

@testset "Two message, one should be executed" begin
    tmsg1 = TimedMessage(1, 1, 100, mockmessage("Ping1"))
    tmsg2 = TimedMessage(2, 2, 200, mockmessage("Ping2"))

    cin = Channel{Tuple{TimedMessage, Int}}(Inf)
    cout = Channel{TimedMessage}(Inf)
    put!(cin, (tmsg1, 1))
    put!(cin, (tmsg2, 1))

    sched = TimedMessage[]
    heartbeat!(sched, 150, cin, cout)

    @test !isready(cin)

    @test length(sched) == 1
    @test sched[1].id == 2

    @test isready(cout)
    tmsg = take!(cout)
    @test !isready(cout)
    @test tmsg.id == 1
end

@testset "Message and remove, same sender" begin
    tmsg1 = TimedMessage(1, 1, 100, mockmessage("Ping1", sender_id = 1))
    tmsg2 = TimedMessage(1, 2, 0, mockmessage("Ping2", sender_id = 1))

    cin = Channel{Tuple{TimedMessage, Int}}(Inf)
    cout = Channel{TimedMessage}(Inf)
    put!(cin, (tmsg1, 1))

    sched = TimedMessage[]
    heartbeat!(sched, 50, cin, cout)

    put!(cin, (tmsg2, -1))
    heartbeat!(sched, 51, cin, cout)

    @test isempty(sched)
    @test !isready(cin)
    @test !isready(cout)
end

@testset "Message and remove, different senders" begin
    tmsg1 = TimedMessage(1, 1, 100, mockmessage("Ping1", sender_id = 1))
    tmsg2 = TimedMessage(1, 2, 0, mockmessage("Ping2", sender_id = 2))

    cin = Channel{Tuple{TimedMessage, Int}}(Inf)
    cout = Channel{TimedMessage}(Inf)
    put!(cin, (tmsg1, 1))

    sched = TimedMessage[]
    heartbeat!(sched, 50, cin, cout)

    put!(cin, (tmsg2, -1))
    heartbeat!(sched, 51, cin, cout)

    @test !isempty(sched)
    @test length(sched) == 1
    @test sched[1].id == 1
    @test !isready(cin)
    @test !isready(cout)
end

@testset "cronworker" begin
    ts = curts()
    tmsg1 = TimedMessage(1, 1, ts - 10_000, mockmessage("Ping1"))
    tmsg2 = TimedMessage(2, 2, ts + 10_000, mockmessage("Ping2"))

    cin = Channel{Tuple{TimedMessage, Int}}(Inf)
    cout = Channel{TimedMessage}(Inf)
    put!(cin, (tmsg1, 1))
    put!(cin, (tmsg2, 1))

    sched = TimedMessage[]
    # I can't quite get the logic of tasks, but ok...
    task = @async cron_worker(cin, cout, sched, 0.05)
    sleep(0.5)

    @test !isready(cin)
    @test isready(cout)
    tmsg = take!(cout)
    @test tmsg.id == 1
    @test length(sched) == 1
    @test sched[1].id == 2
end

end


text:
eltype(::Type{R}) where {T, R<:AbstractSemiringElement{T}} = T

==(l::R, r::S) where {R<:AbstractSemiringElement, S<:AbstractSemiringElement} = issamesemiring(R, S) ? val(l) == val(r) : throw(DomainError("incompatable semirings"))
≈(l::R, r::S) where {R<:AbstractSemiringElement, S<:AbstractSemiringElement} = issamesemiring(R, S) ? val(l) ≈ val(r) : throw(DomainError("incompatable semirings"))

zero(::R) where R<:AbstractSemiringElement = zero(R)
one(::R) where R<:AbstractSemiringElement = one(R)
iszero(x::R) where {R<:AbstractSemiringElement} = x == zero(R)
isone(x::R) where {R<:AbstractSemiringElement} = x == one(R)

convert(::Type{R}, x) where R<:AbstractSemiringElement = R(x)
convert(::Type{R}, x::Number) where R<:AbstractSemiringElement = R(x) # resolving ambiguity
convert(::Type{R}, x::S) where {R<:AbstractSemiringElement, S<:AbstractSemiringElement} = issamesemiring(R, S) ? R(val(x)) : throws(DomainError("incompatable semirings"))
convert(::Type{R}, x::R) where R<:AbstractSemiringElement = x

promote_rule(::Type{R}, ::Type{R}) where R<:AbstractSemiringElement = R
function promote_rule(::Type{R}, ::Type{S}) where {T, U, R<:AbstractSemiringElement{T}, S<:AbstractSemiringElement{U}}
    Z = typejoin(R, S)
    isabstracttype(Z) ? throw(DomainError("incompatable semiring types")) : Z{promote_type(T, U)}
end

# \oplus and \otimes alias for + and *

⊕(l::AbstractSemiringElement, r::AbstractSemiringElement) = +(l, r)
⊙(l::AbstractSemiringElement, r::AbstractSemiringElement) = *(l, r)

(::Type{R})(v::T) where {T, R<:AbstractSemiringElement} = R(convert(T, v))

# some operations on the semiring type

semiring(::S) where S<: AbstractSemiringElement = S
issamesemiring(::Type{R}, ::Type{S}) where {R<:AbstractSemiringElement, S<: AbstractSemiringElement} = !isabstracttype(typejoin(R, S))
issamesemiring(x::R, y::S) where {S<:AbstractSemiringElement, R<:AbstractSemiringElement} = issamesemiring(R, S)

# implementations stemming from the optional interface

-(l::AbstractSemiringElement, r::AbstractSemiringElement) = l+(-r)
inv(x::R) where R<:AbstractSemiringElement = one(R)/x
mulinv(x::AbstractSemiringElement) = inv(x)
addinv(x::R) where R<:AbstractSemiringElement = -x


text:
#####
# some common components for Prioritized Experience Replay based methods
#####

const PERLearners = Union{PrioritizedDQNLearner,RainbowLearner,IQNLearner}

function RLBase.update!(learner::Union{DQNLearner, PERLearners}, t::AbstractTrajectory)
    length(t[:terminal]) < learner.min_replay_history && return

    learner.update_step += 1

    if learner.update_step % learner.target_update_freq == 0
        copyto!(learner.target_approximator, learner.approximator)
    end

    learner.update_step % learner.update_freq == 0 || return

    inds, batch = sample(learner.rng, t, learner.sampler)

    if t isa PrioritizedTrajectory
        priorities = update!(learner, batch)
        t[:priority][inds] .= priorities
    else
        update!(learner,batch)
    end
end

function RLBase.update!(trajectory::PrioritizedTrajectory, p::QBasedPolicy{<:PERLearners}, env::AbstractEnv, ::PostActStage)
    push!(trajectory[:reward], get_reward(env))
    push!(trajectory[:terminal], get_terminal(env))
    push!(trajectory[:priority], p.learner.default_priority)
end


text:
export tiar

using IterativeSolvers
using LinearAlgebra
using Random

"""
    tiar(nep,[maxit=30,][σ=0,][γ=1,][linsolvecreator=DefaultLinSolverCreator(),][tolerance=eps()*10000,][neigs=6,][errmeasure,][v=rand(size(nep,1),1),][logger=0,][check_error_every=1,][orthmethod=DGKS,][proj_solve=false,][inner_solver_method=DefaultInnerSolver(),][inner_logger=0])

Run the tensor infinite Arnoldi method on the nonlinear eigenvalue problem stored in `nep`. This is equivalent to `iar`, but handles orthogonalization with
a tensor representation.

The target `σ` is the center around which eiganvalues are computed.
The value `γ` corresponds
to scaling and specifying a shift and scaling is effectively the same as the
transformation `λ=γs+σ` where `s`
is now the eigenvalue parameter. If you want eigenvalues in a disk centered,
select `σ` as the center
of the disk and `γ` as the radius.
The vector `v` is the starting vector for constructing the Krylov space. The orthogonalization
method, used in contructing the orthogonal basis of the Krylov space, is specified by
`orthmethod`, see the package `IterativeSolvers.jl`.
The iteration is continued until `neigs` Ritz pairs have converged.
This function throws a `NoConvergenceException` if the wanted eigenpairs are not computed after `maxit` iterations.
However, if `neigs` is set to `Inf` the iteration is continued until `maxit` iterations without an error being thrown.
The parameter `proj_solve` determines if the Ritz paris are extracted using the Hessenberg matrix (false),
or as the solution to a projected problem (true). If true, the method is descided by `inner_solver_method`, and the
logging of the inner solvers are descided by `inner_logger`, which works in the same way as `logger`.

See [`augnewton`](@ref) for other parameters.

# Example
```julia-repl
julia> using NonlinearEigenproblems, LinearAlgebra
julia> nep=nep_gallery("dep0",100);
julia> v0=ones(size(nep,1));
julia> λ,v=tiar(nep;v=v0,tol=1e-5,neigs=3);
julia> norm(compute_Mlincomb!(nep,λ[1],v[:,1])) # Is it an eigenvalue?
julia> λ    # print the computed eigenvalues
3-element Array{Complex{Float64},1}:
 0.050462487743188206 - 3.4059600538119376e-18im
 -0.07708769561361105 + 8.611006691570004e-19im
   0.1503916927814904 + 9.388210527944734e-18im
```

# References
* Algorithm 2 in Jarlebring, Mele, Runborg, The Waveguide Eigenvalue Problem and the Tensor Infinite Arnoldi Method, SIAM J. Scient. computing, 39 (3), A1062-A1088, 2017

"""
tiar(nep::NEP;params...)=tiar(ComplexF64,nep;params...)
function tiar(
    ::Type{T},
    nep::NEP;
    orthmethod::Type{T_orth}=DGKS,
    maxit=30,
    linsolvercreator=DefaultLinSolverCreator(),
    tol=eps(real(T))*10000,
    neigs=6,
    errmeasure::ErrmeasureType = DefaultErrmeasure(nep),
    σ=zero(T),
    γ=one(T),
    v=randn(real(T),size(nep,1)),
    logger=0,
    check_error_every=1,
    proj_solve=false,
    inner_solver_method=DefaultInnerSolver(),
    inner_logger=0)where{T,T_orth<:IterativeSolvers.OrthogonalizationMethod}

    @parse_logger_param!(logger)
    @parse_logger_param!(inner_logger)

    # Ensure types σ and v are of type T
    σ=T(σ)
    v=Array{T,1}(v)

    # initialization
    n = size(nep,1); m = maxit;

    if n<m
        msg="Loss of orthogonality in the matrix Z. The problem size is too small, use iar instead.";
        throw(LostOrthogonalityException(msg))
    end

    # initialize variables
    a  = zeros(T,m+1,m+1,m+1);
    Z  = zeros(T,n,m+1);
    t  = zeros(T,m+1);
    tt = zeros(T,m+1);
    g  = zeros(T,m+1,m+1);
    f  = zeros(T,m+1,m+1);
    ff = zeros(T,m+1,m+1);
    H  = zeros(T,m+1,m);
    h  = zeros(T,m+1);
    hh = zeros(T,m+1);
    y  = zeros(T,n,m+1);
    α=Array{T,1}(γ.^(0:m)); α[1]=zero(T);
    local M0inv::LinSolver=create_linsolver(linsolvercreator,nep,σ)
    err = NaN*ones(m+1,m+1);
    λ=zeros(T,m+1); Q=zeros(T,n,m+1);
    Z[:,1]=v; Z[:,1]=Z[:,1]/norm(Z[:,1]);
    a[1,1,1]=one(T);

    # temp var for plot
    conv_eig_hist=zeros(Int,m+1)

    local pnep::NEP;
    if (proj_solve)
        pnep=create_proj_NEP(nep,maxit,T);
    end

    k=1; conv_eig=0;
    while (k <= m)&(conv_eig<neigs)

        # computation of y[:,2], ..., y[:,k+1]
        y[:,2:k+1]=Z[:,1:k]*transpose(a[1:k,k,1:k])
        broadcast!(/,view(y,:,2:k+1),view(y,:,2:k+1),(1:k)')

        # computation of y[:,1]
        y[:,1] = compute_Mlincomb!(nep,σ,y[:,1:k+1],α[1:k+1]);
        y[:,1] = -lin_solve(M0inv,y[:,1]);

        # Gram–Schmidt orthogonalization in Z
        Z[:,k+1]=y[:,1];
        t[k+1] = orthogonalize_and_normalize!(view(Z,:,1:k), view(Z,:,k+1), view(t,1:k), orthmethod)

        # compute the matrix G
        for l=1:k+1
            for i=2:k+1
                g[i,l]=a[i-1,k,l]/(i-1);
            end
         g[1,l]=t[l];
        end

        # compute h (orthogonalization with tensors factorization)
        h = zero(h)
        Ag = zero(h[1:k])
        for l=1:k
            mul!(Ag,a[1:k,1:k,l]',g[1:k,l])
            h[1:k] .+= Ag;
        end

        # compute the matrix F
        f=g;
        Ah = zero(f[1:k+1,1])
        for l=1:k
            mul!(Ah,a[1:k+1,1:k,l],h[1:k])
            f[1:k+1,l] .-= Ah;
        end

        # re-orthogonalization
        # compute hh (re-orthogonalization with tensors factorization)
        hh = zero(hh)
        Af = zero(hh[1:k])
        for l=1:k
            mul!(Af,a[1:k,1:k,l]',f[1:k,l])
            hh[1:k] .+= Af;
        end

        # compute the matrix FF
        ff=f;
        Ah=zero(ff[1:k+1,1])
        for l=1:k
            mul!(Ah,a[1:k+1,1:k,l],hh[1:k])
            ff[1:k+1,l] .-= Ah;
        end

        # update the orthogonalization coefficients
        h=h+hh; f=ff;
        β=norm(view(f,1:k+1,1:k+1)); # equivalent to Frobenius norm

        # extend the matrix H
        H[1:k,k]=h[1:k]; H[k+1,k]=β;

        # extend the tensor
        for i=1:k+1
            for l=1:k+1
                a[i,k+1,l]=f[i,l]/β;
            end
        end

        # compute Ritz pairs (every p iterations)
        if (rem(k,check_error_every)==0)||(k==m)
            D,W = eigen(H[1:k,1:k])
            VV=Z[:,1:k]*transpose(a[1,1:k,1:k])	# extract proper subarray
            Q = VV*W
            λ = σ .+ γ ./ D

            if (proj_solve)  # Projected solve to extract eigenvalues (otw hessenberg matrix)
                set_projectmatrices!(pnep,Z[:,1:k],Z[:,1:k]);
                # Make a call to the inner solve method
                λproj,Qproj=inner_solve(inner_solver_method,T,pnep,
                                        λv=copy(λ),
                                        neigs=size(λ,1)+3,
                                        σ=σ,
                                        tol=tol/10,
                                        innerlogger=inner_logger);

                II=sortperm(abs.(λproj .- σ));
                λproj=λproj[II]; Qproj=Qproj[:,II];
                Q=Z[:,1:k]*Qproj;
                λ=λproj;
             end

            conv_eig=0;
            err[k,1:size(λ,1)]=
              map(s-> estimate_error(errmeasure,λ[s],Q[:,s]), 1:size(λ,1))
            # Log them and compute the converged
            push_iteration_info!(logger,2, k,err=err[k,1:size(λ,1)], λ=λ,
                                 continues=true);
            for s=1:size(λ,1)
                if err[k,s]<tol;
                    conv_eig=conv_eig+1;
                    push_info!(logger,"+", continues=true);
                elseif err[k,s]<tol*10
                    push_info!(logger,"=", continues=true);
                else
                    push_info!(logger,"-", continues=true);
                end
            end
            push_info!(logger,"");

            idx=sortperm(err[k,1:k]); # sort the error
            err[k,1:k]=err[k,idx];
            # extract the converged Ritzpairs
            if (k==m)||(conv_eig>=neigs)
                nrof_eigs = Int(min(length(λ),neigs))
                λ=λ[idx[1:nrof_eigs]]
                Q=Q[:,idx[1:nrof_eigs]]
            end
            conv_eig_hist[k]=conv_eig
        end
        k=k+1;
    end
    k=k-1
    # NoConvergenceException
    if conv_eig<neigs && neigs != Inf
        err=err[end,1:neigs];
        idx=sortperm(err); # sort the error
        λ=λ[idx];  Q=Q[:,idx]; err=err[idx];
        msg="Number of iterations exceeded. maxit=$(maxit)."
        if conv_eig<3
            msg=string(msg, " Check that σ is not an eigenvalue.")
        end
        throw(NoConvergenceException(λ,Q,err,msg))
    end

    # extract the converged Ritzpairs
    λ=λ[1:min(length(λ),conv_eig)];
    Q=Q[:,1:min(size(Q,2),conv_eig)];
    return λ,Q,Z[:,1:k],conv_eig_hist
end


text:
using Test
using Dorothy

@DorothyAll()

@testset "Geometry" begin

	@testset "Distances" begin
		P1 = [0.0, 0.0, 0.0]
		P2 = [1.0, 1.0, 1.0]
		P3 = [2.0, 3.0, 4.0]
		@test sqnorm(P1) == 0.0
		@test sqnorm(P2) == 3.0
		@test sqnorm(P3) == 29.0
		@test dist(P1, P1) == 0.0
		@test dist(P1, P2) == dist(P2, P1) == √3.0
		@test dist(P1, P3) == dist(P3, P1) == √29.0
		@test sqdist(P1, P1) == 0.0
		@test sqdist(P1, P2) == sqdist(P2, P1) == 3.0
		@test sqdist(P1, P3) == sqdist(P3, P1) == 29.0
	end

	@testset "Angles" begin
		P1 = [1.0, 0.0, 0.0]
		P2 = [0.0, 0.0, 0.0]
		P3 = [0.0, 1.0, 0.0]
		P4 = [1.0, 1.0, 1.0]
		@test angle(P1, P2, P1) == 0.0
		@test angle(P1, P2, P3) ≈ angle(P3, P2, P1) ≈ τ/4
		@test angle(P2, P1, P3) ≈ angle(P2, P3, P1) ≈ τ/8
		@test angle(P1, P2, P4) ≈ angle(P4, P2, P1) ≈ 0.9553166181245092
	end

	@testset "Dihedrals" begin
		P1 = [1.0, 0.0, 0.0]
		P2 = [0.0, 0.0, 0.0]
		P3 = [0.0, 0.0, 1.0]
		P4 = [1.0, 0.0, 1.0]
		P5 = [-1.0, 0.0, 1.0]
		P6 = [0.0, 1.0, 1.0]
		P7 = [0.0, -1.0, 1.0]
		@test dihedral(P1, P2, P3, P4) == dihedral(P4, P3, P2, P1) == 0.0
		@test dihedral(P1, P2, P3, P5) == dihedral(P5, P3, P2, P1) == τ/2
		@test dihedral(P1, P2, P3, P6) == dihedral(P6, P3, P2, P1) == τ/4
		@test dihedral(P1, P2, P3, P7) == dihedral(P7, P3, P2, P1) == -τ/4
	end

	@testset "COM" begin
		P1 = [1.0, 0.0, 0.0]
		P2 = [0.0, 1.0, 0.0]
		P3 = [0.0, 0.0, 1.0]
		R = Vector3D[P1, P2, P3]
		@test cog(R) == [1/3, 1/3, 1/3]
	end

	@testset "Transformations" begin
		P1 = [1.0, 0.0, 0.0]
		@test translation(x = 1.0) * P1 == [2.0, 0.0, 0.0]
		rotatez4 = rotation(z = τ/4)
		@test rotatez4 * P1 ≈ [0.0, 1.0, 0.0]
		rotatez4 = rotation(τ/4, [0.0, 0.0, 1.0])
		@test rotatez4 * P1 ≈ [0.0, 1.0, 0.0]
		rotatez4 = rotation(τ/4, [0.0, 0.0, 1.0], O = [0.0, 0.0, 0.0])
		@test rotatez4 * P1 ≈ [0.0, 1.0, 0.0]
		rotatez4 = rotation(τ/4, [0.0, 0.0, 10.0] - [0.0, 0.0, 1.0])
		@test rotatez4 * P1 ≈ [0.0, 1.0, 0.0]
		scale2 = scaling(2.0)
		@test scale2 * P1 == [2.0, 0.0, 0.0]
		@test translation([1.0, 0.0, 1.0]) * scaling(2.0) * rotation(z = τ/4) *
				P1 ≈ [1.0, 2.0, 1.0]
	end

	@testset "Superposition" begin
		m1 = readf("$(Dorothy.datapath)/1BTL.pdb")
		m2 = MolecularModel(m1)
		@test rmsd(m1.R, m2.R) ≈ 0.0
		m1.R .= translation([25.0, 30.0, -1.0]) * rotation(z = τ/4) *
				translation([-5.0, 3.0, 1.0]) .* m2.R
		@test rmsd(m1.R, m2.R) ≈ 36.167398917591804
		T = superposition(m1.R, m2.R)
		m1.R .= T .* m1.R
		@test rmsd(m1.R, m2.R) + 1.0 ≈ 1.0
	end

	@testset "Fit line" begin
		m = readf("$(Dorothy.datapath)/1BTL.pdb")
		infermasses!(m)
		h1 = view(m, 1:120) # First α helix
		V = fitline(h1.R)
		@test V ≈ [0.776711550200754, 0.3946883305177638, 0.4908566894092842]
		V = fitline(h1.R, h1.masses)
		@test V ≈ [0.783371080605647, 0.3869484031172630, 0.4864161627616218]
	end

	@testset "Fit plane" begin
		m = readf("$(Dorothy.datapath)/1BTL.pdb")
		infermasses!(m)
		h18 = view(m, [1:120; 1866:2032]) # First and last α helices
		V = fitplane(h18.R)
		@test V ≈ [0.02187509375361663, -0.815568953125186, 0.578246282280793]
		V = fitplane(h18.R, h18.masses)
		@test V ≈ [0.02428341812381785, -0.828435699705839, 0.559557510053364]
	end

end # @testset


text:
module sp

using StaticArrays: SVector
import StaticArrays: push, setindex
import Base: *, +, ^, promote_rule, convert, show, isless, size, getindex

immutable QuasiStatic{N, T} <: AbstractVector{T}
    data::SVector{N, T}
    length::Int

    function QuasiStatic{N, T}() where {N, T}
        data = SVector{N, T}(ntuple(_ -> T(), Val{N}))
        new{N, T}(data, 0)
    end
    QuasiStatic{N, T}(data::SVector{N, T}, length=N) where {N, T} = new{N, T}(data, length)
end

function QuasiStatic{N}(data::SVector{M, T}) where {N, M, T}
    q = QuasiStatic{N, T}()
    for d in data
        q = push(q, d)
    end
    q
end

QuasiStatic(data::SVector{M, T}, length=M) where {M, T} = QuasiStatic{M, T}(data, length)

@inline getindex(v::QuasiStatic, i) = (@boundscheck (i >= 1 && i <= v.length) || throw(BoundsError()); v.data[i])
size(v::QuasiStatic) = (v.length,)
function push(v::QuasiStatic{N}, x) where {N}
    if v.length >= N
        error("Capacity reached")
    end
    QuasiStatic(setindex(v.data, x, v.length + 1), v.length + 1)
end

abstract type PolynomialLike end
abstract type TermLike <: PolynomialLike end
abstract type MonomialLike <: TermLike end
abstract type VariableLike <: MonomialLike end

immutable Variable{Name} <: VariableLike
end

show(io::IO, v::Variable{Name}) where Name = print(io, Name)

immutable Monomial{N, V} <: MonomialLike
    exponents::NTuple{N, Int64}
end

exponents(m::Monomial) = m.exponents
isless(m1::Monomial, m2::Monomial) = variables(m1) < variables(m2) && exponents(m1) < exponents(m2)
format_exponent(e) = e == 1 ? "" : "^$e"

function show(io::IO, m::Monomial)
    for (i, v) in enumerate(variables(m))
        if m.exponents[i] != 0
            print(io, v, format_exponent(m.exponents[i]))
        end
    end
end

variables(::Type{Monomial{N, V}}) where {N, V} = V
variables(m::Monomial) = variables(typeof(m))

convert(::Type{<:Monomial}, v::Variable{Name}) where {Name} = Monomial{1, (Name,)}((1,))

*(v1::Variable{N1}, v2::Variable{N2}) where {N1, N2} = Monomial{2, (N1, N2)}((1, 1))
*(v1::Variable{N}, v2::Variable{N}) where N = Monomial{1, (N,)}((2,))

immutable Term{T, MonomialType <: Monomial} <: TermLike
    coefficient::T
    monomial::MonomialType
end

(::Type{TermType})() where {T, N, V, TermType <: Term{T, Monomial{N, V}}} = TermType(0, Monomial{N, V}(ntuple(_ -> 0, Val{N})))

format_coefficient(c) = c == 1 ? "" : string(c)
show(io::IO, t::Term) = print(io, format_coefficient(t.coefficient), t.monomial)

convert(T::Type{<:Term}, m::Monomial) = T(1, m)
convert(T::Type{<:Term}, v::Variable) = T(Monomial(v))

@generated function convert(::Type{Term{T, Mono1}}, t::Term{T, Mono2}) where {T, Mono1, Mono2}
    args = Any[0 for v in variables(Mono1)]
    for (j, var) in enumerate(variables(Mono2))
        I = find(v -> v == var, variables(Mono1))
        if isempty(I)
            throw(InexactError())
        elseif length(I) > 1
            error("Duplicate variables in destination $Mono1")
        end
        args[I[1]] = :(t.monomial.exponents[$j])
    end
    quote
        Term{$T, $Mono1}(t.coefficient,
            Monomial{$(length(args)), $(variables(Mono1))}($(Expr(:tuple, args...))))
    end
end

immutable Polynomial{T <: Term, V <: AbstractVector{T}} <: PolynomialLike
    terms::V
end

Polynomial(terms::V) where {T <: Term, V <: AbstractVector{T}} = Polynomial{T, V}(terms)
convert(T::Type{<:Polynomial}, t::Term) = T(QuasiStatic(SVector(t)))
convert(T::Type{<:Polynomial}, m::Monomial) = T(Term(m))
convert(T::Type{<:Polynomial}, v::Variable) = T(Term(v))

# convert(T::Type{Polynomial{T, V1}}, p::Polynomial{T, V2}) where {T, V1, V2} = T(convert(V1, p.terms))

convert(T::Type{Polynomial{T1, V1}}, p::Polynomial{T2, V2}) where {T1, T2, V1, V2} = T(convert(V1, p.terms))

function show(io::IO, p::Polynomial)
    if !isempty(p.terms)
        print(io, p.terms[1])
        for i in 2:length(p.terms)
            print(io, " + ", p.terms[i])
        end
    end
end

function promote_rule(::Type{<:MonomialLike}, ::Type{<:MonomialLike})
    Monomial
end

function promote_rule(::Type{<:TermLike}, ::Type{<:TermLike})
    Term
end

@generated function promote_rule(::Type{Term{T, Mono1}}, ::Type{Term{T, Mono2}}) where {T, Mono1, Mono2}
    vars = Tuple(sort(collect(union(Set(variables(Mono1)),
                                    Set(variables(Mono2))))))
    quote
        Term{T, Monomial{$(length(vars)), $(vars)}}
    end
end

function promote_rule(::Type{<:PolynomialLike}, ::Type{<:PolynomialLike})
    Polynomial
end

@generated function promote_rule(::Type{Polynomial{T1, V1}}, ::Type{Polynomial{T2, V2}}) where {T1, T2, V1, V2}
    termtype = promote_rule(T1, T2)
    quote
        Polynomial{$termtype, Vector{$termtype}}
    end
end

function (+)(v1::Variable{Name}, v2::Variable{Name}) where {Name}
    Term(v1) + Term(v2)
end

function (+)(t1::Term{T, Mono}, t2::Term{T, Mono}) where {T, Mono}
    if t1.monomial < t2.monomial
        Polynomial(QuasiStatic(SVector(t1, t2)))
    elseif t1.monomial > t2.monomial
        Polynomial(QuasiStatic(SVector(t2, t1)))
    else
        terms = QuasiStatic{2, Term{T, Mono}}()
        terms = push(terms, Term{T, Mono}(t1.coefficient + t2.coefficient, t1.monomial))
        Polynomial(terms)
        # Polynomial([Term{T, Mono}(t1.coefficient + t2.coefficient, t1.monomial)])
    end
end

@generated function (+)(p1::Polynomial{T1, QuasiStatic{N1, T1}}, p2::Polynomial{T2, QuasiStatic{N2, T2}}) where {T1, N1, T2, N2}
    T = promote_type(T1, T2)
    Q = QuasiStatic{N1 + N2, T}
    quote
        terms = $Q()
        i1 = 1
        i2 = 1
        while i1 <= length(p1.terms) && i2 <= length(p2.terms)
            t1 = convert($T, p1.terms[i1])
            t2 = convert($T, p2.terms[i2])
            if t1.monomial < t2.monomial
                terms = push(terms, t1)
                i1 += 1
            elseif t1.monomial > t2.monomial
                terms = push(terms, t2)
                i2 += 1
            else
                terms = push(terms, $T(t1.coefficient + t2.coefficient,
                             t1.monomial))
                i1 += 1
                i2 += 1
            end
        end
        for i in i1:length(p1.terms)
            terms = push(terms, p1.terms[i])
        end
        for i in i2:length(p2.terms)
            terms = push(terms, p2.terms[i])
        end
        Polynomial(terms)
    end
end

@generated function (*)(m1::Monomial{N1, V1}, m2::Monomial{N2, V2}) where {N1, V1, N2, V2}
    vars = Tuple(sort(collect(union(Set(V1), Set(V2)))))
    args = []
    for (i, v) in enumerate(vars)
        i1 = findfirst(V1, v)
        i2 = findfirst(V2, v)
        if i1 != 0 && i2 != 0
            push!(args, :(m1.exponents[$i1] + m2.exponents[$i2]))
        elseif i1 != 0
            push!(args, :(m1.exponents[$i1]))
        else
            @assert i2 != 0
            push!(args, :(m2.exponents[$i2]))
        end
    end
    Expr(:call, :(Monomial{$(length(args)), $vars}), Expr(:tuple, args...))
end

(*)(t1::Term, t2::Term) = Term(t1.coefficient * t2.coefficient, t1.monomial * t2.monomial)
(+)(m1::Monomial, m2::Monomial) = Term(m1) + Term(m2)

(+)(t1::PolynomialLike, t2::PolynomialLike) = +(promote(t1, t2)...)
(*)(t1::TermLike, t2::TermLike) = *(promote(t1, t2)...)
*(x::Number, m::Monomial) = Term(x, m)
*(x::Number, v::Variable) = x * Monomial(v)
^(v::Variable{Name}, x::Integer) where Name = Monomial{1, (Name,)}((x,))

end


text:
using CompScienceMeshes

Γ = meshsphere(1.0, 0.3)
γ = bisecting_refinement(Γ)

include(Pkg.dir("CompScienceMeshes","examples","plotlyjs_patches.jl"))
p = patch(γ)
PlotlyJS.plot(p)


text:
module JungleHelperBouncyShroom

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/BouncyShroomUp" BouncyShroomUp(x::Integer, y::Integer, yeetx::Integer = 200, yeety::Integer = -290, spriteDirectory::String="JungleHelper/BouncyShroom")
@mapdef Entity "JungleHelper/BouncyShroomLeft" BouncyShroomLeft(x::Integer, y::Integer, yeetx::Integer = 200, yeety::Integer = -290, spriteDirectory::String="JungleHelper/BouncyShroom")
@mapdef Entity "JungleHelper/BouncyShroomRight" BouncyShroomRight(x::Integer, y::Integer, yeetx::Integer = 200, yeety::Integer = -290, spriteDirectory::String="JungleHelper/BouncyShroom")

const placements = Ahorn.PlacementDict()

BouncyShroom = Dict{String, Type}(
    "up" => BouncyShroomUp,
    "left" => BouncyShroomLeft,
    "right" => BouncyShroomRight
)

BouncyShroomUnion = Union{BouncyShroomUp, BouncyShroomLeft, BouncyShroomRight}

for (dir, entity) in BouncyShroom
	key = "Bouncy Shroom ($(uppercasefirst(dir))) (Jungle Helper)"
	placements[key] = Ahorn.EntityPlacement(
		entity
	)
end

directions = Dict{String, String}(
    "JungleHelper/BouncyShroomUp" => "up",
    "JungleHelper/BouncyShroomLeft" => "left",
    "JungleHelper/BouncyShroomRight" => "right",
)

function Ahorn.selection(entity::BouncyShroomUnion)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 16, y - 16, 24, 24)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::BouncyShroomUp, room::Maple.Room) = Ahorn.drawSprite(ctx, get(entity, "spriteDirectory", "JungleHelper/BouncyShroom") * "/mushroom00", -4, -4)
Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::BouncyShroomLeft, room::Maple.Room) = Ahorn.drawSprite(ctx, get(entity, "spriteDirectory", "JungleHelper/BouncyShroom") * "/mushroom_ld_00", -4, -4)
Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::BouncyShroomRight, room::Maple.Room) = Ahorn.drawSprite(ctx, get(entity, "spriteDirectory", "JungleHelper/BouncyShroom") * "/mushroom_rd_00", -4, -4)

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
@testset "Gpt" begin
  using Transformers.GenerativePreTrain
  for f ∈ readdir("./gpt/")
    include("./gpt/$f")
  end
end


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

mutable struct IoK8sApiMetricsV1beta1ContainerMetrics <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    usage::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: usage

    function IoK8sApiMetricsV1beta1ContainerMetrics(;name=nothing, usage=nothing)
        o = new()
        validate_property(IoK8sApiMetricsV1beta1ContainerMetrics, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiMetricsV1beta1ContainerMetrics, Symbol("usage"), usage)
        setfield!(o, Symbol("usage"), usage)
        o
    end
end # type IoK8sApiMetricsV1beta1ContainerMetrics

const _property_map_IoK8sApiMetricsV1beta1ContainerMetrics = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("usage")=>Symbol("usage"))
const _property_types_IoK8sApiMetricsV1beta1ContainerMetrics = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("usage")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}")
Base.propertynames(::Type{ IoK8sApiMetricsV1beta1ContainerMetrics }) = collect(keys(_property_map_IoK8sApiMetricsV1beta1ContainerMetrics))
Swagger.property_type(::Type{ IoK8sApiMetricsV1beta1ContainerMetrics }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiMetricsV1beta1ContainerMetrics[name]))}
Swagger.field_name(::Type{ IoK8sApiMetricsV1beta1ContainerMetrics }, property_name::Symbol) =  _property_map_IoK8sApiMetricsV1beta1ContainerMetrics[property_name]

function check_required(o::IoK8sApiMetricsV1beta1ContainerMetrics)
    true
end

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


text:
using OrdinaryDiffEq, DiffEqProblemLibrary.ODEProblemLibrary
using Test

ODEProblemLibrary.importodeproblems()

@testset "init" begin
    prob = ODEProblemLibrary.prob_ode_linear
    prob2D = ODEProblemLibrary.prob_ode_2Dlinear

    inferred = [BS3(), Tsit5(), RK4(), Vern6()]
    for alg in inferred
        @inferred init(prob, alg)
        @inferred init(prob2D, alg)
    end

    notinferred = [SDIRK2(), TRBDF2(), KenCarp4(), Rosenbrock23(), Rodas4()]
    for alg in notinferred
        @test_broken @inferred init(prob, alg)
        @test_broken @inferred init(prob2D, alg)
    end
end


text:
#!/usr/bin/env julia

p(n) = mapreduce(x -> (4 / (1 + (x * x))), +, (0.5 / n):(1 / n):1) / n

function picalc(numsteps)

  # Get number of processes.
  # np = nworkers()
  np = 1

  println("Calculating PI using:")
  println("  ", numsteps, " slices")
  println("  ", np, " worker(s)")
  start = time()

  mypi = p(numsteps)

  stop = time()

  elapsed = (stop - start)

  println("Obtained value of PI: ", mypi)
  println("Time taken: ", round(elapsed; digits=3), " seconds")

end

const numsteps = if length(ARGS) > 0
    parse(Int, ARGS[1])
else
    1_000_000_000
end

# Warm up kernel
print("  Warming up...")
warms = time()
p(10)
warmt = time() - warms
println("done. [", round(warmt; digits=3), "s]\n")

# Run the full example
picalc(numsteps)


text:
# This file is a part of EmpiricalDistributions.jl, licensed under the MIT License (MIT).

_bin_left(edge::AbstractVector{<:Real}) = view(edge, firstindex(edge):lastindex(edge)-1)
_bin_right(edge::AbstractVector{<:Real}) = view(edge, firstindex(edge)+1:lastindex(edge))
_bin_centers(edge::AbstractVector{<:Real}) = (_bin_left(edge) .+ _bin_right(edge)) ./ 2

_bin_widths(edge::AbstractVector{<:Real}) = _bin_right(edge) .- _bin_left(edge)

# Ranges need special treatment, since subtraction of integer ranges results in "step cannot be zero" error
function _bin_widths(edge::AbstractRange{<:Real})
    ref = step(edge)
    T = typeof(ref)
    stp = zero(T)
    n = length(eachindex(edge)) - 1
    StepRangeLen(ref, stp, n)
end

function _ratio(a::Real, b::Real)
    R = float(promote_type(typeof(a), typeof(b)))
    a == b ?  one(R) : convert(R, a / b)
end

function _linear_interpol(x_lo::Real, x_hi::Real, y_lo::Real, y_hi::Real, x::Real)
    T = promote_type(typeof(y_lo), typeof(y_hi), typeof(x))
    w_hi = T(_ratio(x - x_lo, x_hi - x_lo))
    w_lo = one(w_hi) - w_hi
    y = y_lo * w_lo + y_hi * w_hi
    @assert y_lo <= y <= y_hi
    y
end

function _linear_interpol(X::AbstractVector{<:Real}, Y::AbstractVector{<:Real}, x::Real)
    i_lo, i_hi = _find_idxs_lohi(X, x)
    x_lo, x_hi = X[i_lo], X[i_hi]
    y_lo, y_hi = Y[i_lo], Y[i_hi]
    _linear_interpol(x_lo, x_hi, y_lo, y_hi, x)
end

function _find_bin(edge::AbstractVector{<:Real}, closed_left::Bool, x::Real)
    if closed_left == false
        searchsortedfirst(edge, x) - 1
    else
        searchsortedlast(edge, x)
    end
end

function _find_bin(edges::NTuple{N,AbstractVector{<:Real}}, closed_left::Bool, x::NTuple{N,Real}) where N
    map(_find_bin, edges, map(_ -> closed_left, x), x)
end

function _find_idxs_lohi(X::AbstractVector{<:Real}, x::Real)
    x_min, x_max = first(X), last(X)
    x_min <= x <= x_max || throw(DomainError(x, "require value between x_min and x_max"))
    idxs = searchsorted(X, x)
    i1, i2 = first(idxs), last(idxs)
    i1, i2 = first(idxs), last(idxs)
    from, to = min(i1, i2), max(i1, i2)
    @assert checkbounds(Bool, X, from:to)
    (from, to)
end

function _mean(h::StatsBase.Histogram{<:Real, N}) where {N}
    @assert h.isdensity # Implementation requires normalized histogram

    Y = h.weights
    X = _bin_centers.(h.edges)
    W = _bin_widths.(h.edges)

    mean_est = mapreduce(
        idx -> map(getindex, X, idx.I) .* Y[idx] .* prod(map(getindex, W, idx.I)),
        (a, b) -> a .+ b,
        CartesianIndices(Y)
    ) ./ sum(idx -> Y[idx] .* prod(map(getindex, W, idx.I)), CartesianIndices(Y))
end

_findmaxidx_tuple_or_int(A::AbstractVector{<:Real}) = findmax(A)[2]
_findmaxidx_tuple_or_int(A::AbstractArray{<:Real}) = findmax(A)[2].I

function _mode(h::StatsBase.Histogram)
    @assert h.isdensity # Implementation requires normalized histogram

    maxidx = _findmaxidx_tuple_or_int(h.weights)
    mode_corner1 = map(getindex, h.edges, maxidx)
    mode_corner2 = map(getindex, h.edges, maxidx .+ 1)
    (mode_corner1 .+ mode_corner2) ./ 2
end

function _var(h::StatsBase.Histogram{<:Real,N}, mean_est::NTuple{N,Real}) where {N}
    @assert h.isdensity # Implementation requires normalized histogram

    Y = h.weights
    X = _bin_centers.(h.edges)
    W = _bin_widths.(h.edges)

    mean_est = mapreduce(
        idx -> (map(getindex, X, idx.I) .- mean_est).^2 .* Y[idx] .* prod(map(getindex, W, idx.I)),
        (a, b) -> a .+ b,
        CartesianIndices(Y)
    ) ./ sum(idx -> Y[idx] .* prod(map(getindex, W, idx.I)), CartesianIndices(Y))
end

function _cov(h::StatsBase.Histogram{<:Real,N}, mean_est::NTuple{N,Real}) where {N}
    @assert h.isdensity # Implementation requires normalized histogram

    Y = h.weights
    X = _bin_centers.(h.edges)
    W = _bin_widths.(h.edges)

    dim_idxs = 1:N

    [
        mapreduce(
            idx -> (X[i][idx.I[i]] - mean_est[i]) * (X[j][idx.I[j]] - mean_est[j]) .* Y[idx] .* prod(map(getindex, W, idx.I)),
            (a, b) -> a .+ b,
            CartesianIndices(Y)
        ) ./ sum(idx -> Y[idx] .* prod(map(getindex, W, idx.I)), CartesianIndices(Y))
        for i in dim_idxs, j in dim_idxs
    ]
end

_rand_uniform(rng::AbstractRNG, T::Type{<:Real}, lo::Real, hi::Real) = fma(hi - lo, rand(rng, T), lo)

# For tests of adapt()
struct _AdaptToArray end
Adapt.adapt_storage(::Val{_AdaptToArray}, xs::AbstractArray) = convert(Array, xs)


text:
X = [1 2; 3 4]
@test mat_embed(X, (2,3), zero_pos=[2, 3], tpin="rec", tpout="rec") == [1 0 2; 0 3 4]
@test mat_embed(X, (3,2), zero_pos=[3, 4], tpin="rec", tpout="rec") == [1 0; 3 2; 0 4]
@test mat_embed(X, (4,4), zero_pos=[3, 4], tpin="rec", tpout="her") == [0 1.0 3.0 0.0; 1.0 0 0.0 2.0; 3.0 0.0 0 4.0; 0.0 2.0 4.0 0]

X = [2 3; 3 2]
@test mat_embed(X, (2,3), zero_pos=[1,2,3], tpin="her", tpout="rec") == [0 0 3; 0 2 2]
@test mat_embed(X, (2,2), zero_pos=[4], tpin="her", tpout="rec") == [2 2; 3 0]
@test mat_embed(X, (3,3), zero_pos=[], tpin="her", tpout="her") == [0 2 3; 2 0 2; 3 2 0]


text:
# Autogenerated wrapper script for TMatrix_jll for i686-linux-gnu-libgfortran3
export tmatrix_fixed_orientation, tmatrix_random_orientation

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("TMatrix")
JLLWrappers.@declare_library_product(tmatrix_fixed_orientation, "libtmatrixfo.so")
JLLWrappers.@declare_library_product(tmatrix_random_orientation, "libtmatrixro.so")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        tmatrix_fixed_orientation,
        "lib/libtmatrixfo.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

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

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Catalyst
using Polyurethanes

#
# Reaction without evaporation
#
reaction_noevap = @reaction_network begin
    k1, U --> NCO + DIPA_l
    k2, NCO + OH --> POL
end k1 k2

#
# Reaction with evaporation
#
vap(OH,A,b) = A * exp(OH / b)
reaction_evap = @reaction_network begin
    (k1, km1), U <--> NCO + DIPA_l
    k2, NCO + OH --> POL
    vap(OH, A, b), DIPA_l --> DIPA_v
end k1 km1 k2 A b

# Directory of the data
dir = "./data"

#
# BD systems, without evaporation
#
bd_noevap = [
    System(
        title="BD IPDI 110C 1",
        label="NCO/DIPA = 1.00",
        experimental_data=readdlm("$dir/bd_ipdi_110C_1.0.dat"),
        c0=(U = 3.4675, NCO = 0., DIPA_l = 0., OH = 3.4675, POL = 0.),
        p0=(k1 = 6.8e-6, k2 = 2.2e-4),
        lower=(k1 = 0., k2 = 0.),
        upper=(k1 = +Inf, k2 = +Inf)
    )
    System(
        title="BD IPDI 110C 0.5",
        label="NCO/DIPA = 0.50",
        experimental_data=readdlm("$dir/bd_ipdi_110C_0.5.dat"),
        c0=(U = 2.246154, NCO = 2.246154, DIPA_l = 0.,	OH = 4.49231, POL = 0.),
        p0=(k1 = 6.8e-6, k2 = 7.7e-5),
        lower=(k1 = 0., k2 = 0.),
        upper=(k1 = +Inf, k2 = +Inf)
    )
    System(
        title="BD IPDI 110C 0.25",
        label="NCO/DIPA = 0.25",
        experimental_data=readdlm("$dir/bd_ipdi_110C_0.25.dat"),
        c0=(U = 1.603468, NCO = 4.810405, DIPA_l = 0.,	OH = 6.41387, POL = 0.),
        p0=(k1 = 6.8e-6, k2 = 7.7e-5),
        lower=(k1 = 0., k2 = 0.),
        upper=(k1 = +Inf, k2 = +Inf)
    )
]

#
# BD systems with evaporation
#
bd_evap = [
    System(
        title="BD IPDI 110C 1",
        label="NCO/DIPA = 1.00",
        experimental_data=readdlm("$dir/bd_ipdi_110C_1.0.dat"),
        c0=(U = 3.4675, NCO = 0., DIPA_l = 0., OH = 3.4675, POL = 0., DIPA_v = 0.),
        p0=(k1 = 6.8e-6, km1 = 4.5e-4, k2 = 2.2e-4, A = 1.13e-6, b = 0.5),
        lower=(k1 = 0., km1 = 0., k2 = 0., A = 0., b = -Inf),
        upper=(k1 = +Inf, km1 = +Inf, k2 = +Inf, A = +Inf, b = +Inf)
    )
    System(
        title="BD IPDI 110C 0.5",
        label="NCO/DIPA = 0.50",
        experimental_data=readdlm("$dir/bd_ipdi_110C_0.5.dat"),
        c0=(U = 2.246154, NCO = 2.246154, DIPA_l = 0.,	OH = 4.49231, POL = 0., DIPA_v = 0.),
        p0=(k1 = 6.8e-6, km1 = 4.5e-4, k2 = 7.7e-5, A = 1.13e-6, b = 0.5),
        lower=(k1 = 0., km1 = 0., k2 = 0., A = 0., b = -Inf),
        upper=(k1 = +Inf, km1 = +Inf, k2 = +Inf, A = +Inf, b = +Inf)
    )
    System(
        title="BD IPDI 110C 0.25",
        label="NCO/DIPA = 0.25",
        experimental_data=readdlm("$dir/bd_ipdi_110C_0.25.dat"),
        c0=(U = 1.603468, NCO = 4.810405, DIPA_l = 0.,	OH = 6.41387, POL = 0., DIPA_v = 0.),
        p0=(k1 = 6.8e-6, km1 = 4.5e-4, k2 = 7.7e-5, A = 1.13e-6, b = 0.5),
        lower=(k1 = 0., km1 = 0., k2 = 0., A = 0., b = -Inf),
        upper=(k1 = +Inf, km1 = +Inf, k2 = +Inf, A = +Inf, b = +Inf)
    )
]

#
# PEG at 100oC, with evaporation
#
peg_100_evap = [
    System(
        title="PEG IPDI 100C 1",
        label="NCO/DIPA = 1.00",
        experimental_data=readdlm("$dir/peg_ipdi_100C_1.0.dat"),
        c0=(U = 2.121362, NCO = 0., DIPA_l = 0., OH = 2.12136, POL = 0., DIPA_v = 0.),
        p0=(k1 = 6.8e-6, km1 = 4.5e-4, k2 = 2.2e-4, A = 1.13e-6, b = 0.5),
        lower=(k1 = 0., km1 = 0., k2 = 0., A = 0., b = -Inf),
        upper=(k1 = +Inf, km1 = +Inf, k2 = +Inf, A = +Inf, b = +Inf)
    )
    System(
        title="PEG IPDI 100C 0.25",
        label="NCO/DIPA = 0.25",
        experimental_data=readdlm("$dir/peg_ipdi_100C_0.25.dat"),
        c0=(U = 0.627473, NCO = 1.882418, DIPA_l = 0., OH = 2.50989, POL = 0., DIPA_v = 0.),
        p0=(k1 = 6.8e-6, km1 = 4.5e-4, k2 = 2.2e-4, A = 1.13e-6, b = 0.5),
        lower=(k1 = 0., km1 = 0., k2 = 0., A = 0., b = -Inf),
        upper=(k1 = +Inf, km1 = +Inf, k2 = +Inf, A = +Inf, b = +Inf)
    )
    System(
        title="PEG IPDI 100C 1 100h",
        label="NCO/DIPA = 1. (100h)",
        experimental_data=readdlm("$dir/peg_ipdi_100C_1.0_100h.dat"),
        c0=(U = 2.121362, NCO = 0., DIPA_l = 0., OH = 2.12136, POL = 0., DIPA_v = 0.),
        p0=(k1 = 6.8e-6, km1 = 4.5e-4, k2 = 2.2e-4, A = 1.13e-6, b = 0.5),
        lower=(k1 = 0., km1 = 0., k2 = 0., A = 0., b = -Inf),
        upper=(k1 = +Inf, km1 = +Inf, k2 = +Inf, A = +Inf, b = +Inf)
    )
]

#
# PEG at 110oC, with evaporation
#
peg_110_evap = [
    System(
        title="PEG IPDI 110C 1",
        label="NCO/DIPA = 1.00",
        experimental_data=readdlm("$dir/peg_ipdi_110C_1.0.dat"),
        c0=(U = 2.121362, NCO = 0., DIPA_l = 0., OH = 2.12136, POL = 0., DIPA_v = 0.),
        p0=(k1 = 6.8e-6, km1 = 4.5e-4, k2 = 2.2e-4, A = 1.13e-6, b = 0.5),
        lower=(k1 = 0., km1 = 0., k2 = 0., A = 0., b = -Inf),
        upper=(k1 = +Inf, km1 = +Inf, k2 = +Inf, A = +Inf, b = +Inf)
    )
    System(
        title="PEG IPDI 110C 0.25",
        label="NCO/DIPA = 0.25",
        experimental_data=readdlm("$dir/peg_ipdi_110C_0.25.dat"),
        c0=(U = 0.627473, NCO = 1.882418, DIPA_l = 0., OH = 2.50989, POL = 0., DIPA_v = 0.),
        p0=(k1 = 6.8e-6, km1 = 4.5e-4, k2 = 2.2e-4, A = 1.13e-6, b = 0.5),
        lower=(k1 = 0., km1 = 0., k2 = 0., A = 0., b = -Inf),
        upper=(k1 = +Inf, km1 = +Inf, k2 = +Inf, A = +Inf, b = +Inf)
    )
]

#
# PEG 100oC, without evaporation
#
peg_100_noevap = [
    System(
        title="PEG IPDI 100C 1",
        label="NCO/DIPA = 1.00",
        experimental_data=readdlm("$dir/peg_ipdi_100C_1.0.dat"),
        c0=(U = 2.121362, NCO = 0., DIPA_l = 0., OH = 2.12136, POL = 0.),
        p0=(k1 = 6.8e-6, k2 = 2.2e-4),
        lower=(k1 = 0., k2 = 0.),
        upper=(k1 = +Inf, k2 = +Inf)
    )
    System(
        title="PEG IPDI 100C 0.25",
        label="NCO/DIPA = 0.25",
        experimental_data=readdlm("$dir/peg_ipdi_100C_0.25.dat"),
        c0=(U = 0.627473, NCO = 1.882418, DIPA_l = 0., OH = 2.50989, POL = 0.),
        p0=(k1 = 6.8e-6, k2 = 2.2e-4),
        lower=(k1 = 0., k2 = 0.),
        upper=(k1 = +Inf, k2 = +Inf)
    )
    System(
        title="PEG IPDI 100C 1 100h",
        label="NCO/DIPA = 1. (100h)",
        experimental_data=readdlm("$dir/peg_ipdi_100C_1.0_100h.dat"),
        c0=(U = 2.121362, NCO = 0., DIPA_l = 0., OH = 2.12136, POL = 0.),
        p0=(k1 = 6.8e-6, k2 = 2.2e-4),
        lower=(k1 = 0., k2 = 0.),
        upper=(k1 = +Inf, k2 = +Inf)
    )
]

#
# PEG 110oC without evaporation
#
peg_110_noevap = [
    System(
        title="PEG IPDI 110C 1",
        label="NCO/DIPA = 1.00",
        experimental_data=readdlm("$dir/peg_ipdi_110C_1.0.dat"),
        c0=(U = 2.121362, NCO = 0., DIPA_l = 0., OH = 2.12136, POL = 0.),
        p0=(k1 = 6.8e-6, k2 = 2.2e-4),
        lower=(k1 = 0., k2 = 0.),
        upper=(k1 = +Inf, k2 = +Inf)
    )
    System(
        title="PEG IPDI 110C 0.25",
        label="NCO/DIPA = 0.25",
        experimental_data=readdlm("$dir/peg_ipdi_110C_0.25.dat"),
        c0=(U = 0.627473, NCO = 1.882418, DIPA_l = 0., OH = 2.50989, POL = 0.),
        p0=(k1 = 6.8e-6, k2 = 2.2e-4),
        lower=(k1 = 0., k2 = 0.),
        upper=(k1 = +Inf, k2 = +Inf)
    )
]

function run_system(basename,systems,reaction)
    # Fit each set independently
    for (isys,system) in pairs(systems)
        println("Running: ", system.title)
        result = fit(reaction,system)
        plt = plot(result)
        savefig(plt,"./plots/$(basename)_$isys.png")
        writecsv(result,"./simulations/$(basename)_$isys.csv")
    end
    # Obtain best overall fit
    println("Running fit for all sets...")
    result = fitall(reaction,systems)
    plt = plot(result)
    savefig(plt,"./plots/$basename.png")
    writecsv(result,"./simulations/$(basename)_all.csv")
end

for (basename, system, reaction) in [
                            ("bd_evap", bd_evap, reaction_evap),
                            ("bd_noevap", bd_noevap, reaction_noevap),
                            ("peg_100_evap", peg_100_evap, reaction_evap),
                            ("peg_100_noevap", peg_100_noevap, reaction_noevap),
                            ("peg_110_evap", peg_110_evap, reaction_evap),
                            ("peg_110_noevap", peg_110_noevap, reaction_noevap),
                          ]
    run_system(basename, system, reaction)
end


text:
using LinearAlgebra, ForwardDiff, NLPModels

export Bisseccao, Newton_rc_bissec

function Bisseccao(g, a, b, max_bissec; λ = 0)
    ϵ = 1e-4
    status= :resolvido
    iter = 0
    while abs(g(λ)) > ϵ
        if g(a)*g(b)==0 && g(b)==0
            λ=b
        else
            λ=a
        end
        if g(a) * g(b) < 0
            while abs(b-a) > ϵ
                λ = (b + a) / 2
                if g(λ) * g(a) < 0
                    b = λ
                elseif g(λ) * g(b) < 0
                    a = λ
                end
            end
        end
        if g(a) * g(b) > 0
            if b < 1000
                b = b*5+1e-3
            elseif a < 8000
                a = a + 40
            else
                status= :bisseccao_falhou
                break
            end
        end

        if iter > max_bissec
            println("Maximum bissection.")
            break
        end
        iter +=1
    end
    return λ, status
end

function Newton_rc_bissec(
    nlp::AbstractNLPModel;
    η1::Float64 = 1e-2,
    η2::Float64 = 0.75,
    Δ::Float64 = 5.0,
    max_bissec = 1000,
    a::Float64 = 0.0,
    b::Float64 = 1.0,
    atol::Real = 1e-6,
    rtol::Real = 1e-6,
    max_eval::Int = 1000,
    max_iter::Int = 0,
    max_time::Float64 = 10.0
    )

    if !unconstrained(nlp)
        error("Problem is not unconstrained")
    end

    x = copy(nlp.meta.x0)
    n = nlp.meta.nvar

    f(x) = obj(nlp, x)
    ∇f(x) = grad(nlp, x)
    H(x) = hess(nlp, x)

    fx = f(x)
    ∇fx = ∇f(x)
    Hx = Matrix(H(x))
    g(λ) = norm(inv(Hx + λ * I) * ∇f(x)) - Δ
    d = - Hx \ ∇fx

    ϵ = atol + rtol * norm(∇fx)
    t₀ = time()

    iter = 0
    Δt = time() - t₀
    solved = norm(∇fx) < ϵ # First order stationary
    tired = neval_obj(nlp) ≥ max_eval > 0|| iter ≥ max_iter > 0 || Δt ≥ max_time > 0 # Excess time, iteration, evaluations

    status = :unknown

    @info log_header(
    [:iter, :fx, :ngx, :nf, :Δt],
    [Int, Float64, Float64, Int, Float64],
    hdr_override=Dict(:fx => "f(x)", :ngx => "‖∇f(x)‖", :nf => "#f")
    )

    @info log_row(
    Any[iter, fx, norm(∇fx), neval_obj(nlp), Δt]
    )

    while !(solved || tired)
        if norm(d) < Δ
            d = - Hx \ ∇f(x)
        else
            λ, stat = Bisseccao(g, a, b, max_bissec)
            if stat == :bisseccao_falhou
                @warn("Bissection fail")
                break
            end
            d = - (Hx + λ * I) \ ∇f(x)
        end

        Ared = f(x) - f(x + d)
        Pred = f(x) - (f(x) + dot(d, ∇f(x)) + dot(d, Hx * d) /  2)
        ρ = Ared / Pred
        if ρ < η1
            Δ = Δ / 2
        elseif ρ < η2
            x = x + d
        else
            x = x + d
            Δ = 2Δ
        end

        fx = f(x)
        ∇fx = ∇f(x)
        Hx = Matrix(H(x))

        iter += 1
        Δt = time() - t₀
        solved = norm(∇fx) < ϵ # First order stationary
        tired = neval_obj(nlp) ≥ max_eval > 0|| iter ≥ max_iter > 0 || Δt ≥ max_time > 0 # Excess time, iteration, evaluations

        @info log_row(
        Any[iter, fx, norm(∇fx), neval_obj(nlp), Δt]
        )
    end

    if solved
        status = :first_order
    elseif tired
        if neval_obj(nlp) ≥ max_eval > 0
            status = :max_eval
        elseif iter ≥ max_iter > 0
            status = :max_iter
        elseif Δt ≥ max_time > 0
            status = :max_time
        end
    end

    return GenericExecutionStats(
    status,
    nlp,
    solution=x,
    objective=f(x),
    dual_feas=norm(∇fx),
    elapsed_time=Δt,
    iter=iter
    )
end


text:
# low-level wrappers

function cublasCreate_v2()
  handle = Ref{cublasHandle_t}()
  cublasCreate_v2(handle)
  handle[]
end

function cublasXtCreate()
  handle = Ref{cublasXtHandle_t}()
  cublasXtCreate(handle)
  handle[]
end

function cublasXtGetBlockDim(handle)
  bd = Ref{Int}()
  cublasXtGetBlockDim(handle, bd)
  bd[]
end

function cublasXtGetPinningMemMode(handle)
  mm = Ref{cublasXtPinningMemMode_t}()
  cublasXtGetPinningMemMode(handle, mm)
  mm[]
end

function cublasGetVersion_v2(handle)
  version = Ref{Cint}()
  cublasGetVersion_v2(handle, version)
  version[]
end

function cublasGetProperty(property::CUDAapi.libraryPropertyType)
  value_ref = Ref{Cint}()
  cublasGetProperty(property, value_ref)
  value_ref[]
end


text:
using JLD
using ProgressMeter
using Distributed

# See Wiersig paper for details/definitions of these functions. There is an overall sign change
# in the Green's function to match definitions in Heider's paper.
"""
    ∂νG(b::Boundary, r′::Point2D, t::Float64, n::Complex{Float64}, k::Complex{Float64})

Normal derivative of Green's function at point r′ (external or internal to boundary),
point parameterized by t on the boundary, and specified index of refraction and complex
wave vector k.
"""
function ∂νG(b::Boundary, r′::Point2D, t::Float64, n::Complex{Float64}, k::Complex{Float64})
  r = boundaryPoint(b, t)
  diff = r-r′
  absDiff = norm(diff)
  nv = normVector(b, t)
  # Normalize this normal vector because Wiersig paper defines it as unit vector.
  ν = nv/norm(nv)
  cosα = dot(Point2D(ν), diff)/absDiff
  -im*n*k/4*cosα*BoundaryIntegralMethod.h1(n*k*absDiff)
end

"""
    G(b::Boundary, r′::Point2D, t::Float64, n::Complex{Float64}, k::Complex{Float64})

Green's function at point r′ (external or internal to boundary),
point parameterized by t on the boundary, and specified index of refraction and complex
wave vector k.

"""
function G(b::Boundary, r′::Point2D, t::Float64, n::Complex{Float64}, k::Complex{Float64})
  r = boundaryPoint(b, t)
  absDiff = norm(r-r′)
  1/4*im*BoundaryIntegralMethod.h0(n*k*absDiff)
end

function ds(b::Boundary, j::Int, N::Int)
  arcLen(b, j*π/N, 2N)
end

"""
    ψField(b::Boundary, N::Int, r′::Point2D, n::Complex{Float64}, k::Complex{Float64}, ψ_Γ::Array{Complex{Float64},1}, ∂νψ_Γ::Array{Complex{Float64},1})

Calculate the field at a point given
the field and its normal derivative on the boundary.
"""
function ψField(b::Boundary, N::Int, ptSpacing::Float64, r′::Point2D, n::Complex{Float64}, k::Complex{Float64}, ψ_Γ::Array{Complex{Float64},1}, ∂νψ_Γ::Array{Complex{Float64},1}, dsList::Array{Float64})
  # Note that due to sign differences in the Green's function in Heider's paper, there is a plus in the discrete integral below
  # instead of the customary minus
  fieldVal = try
    sum([dsList[j+1]*(ψ_Γ[j+1]*∂νG(b, r′, j*π/N, n, k) + G(b, r′, j*π/N, n, k)*∂νψ_Γ[j+1]) for j in 0:2N-1])
  catch err
    @error exception=err
    foundPointIdx = -1
    for j in 0:2N-1
      pt = boundaryPoint(b, j*π/N)
      if isapprox(norm(pt-r′), 0.0, atol=ptSpacing/4)
        foundPointIdx = j
        break
      end
    end
    @info "Singular field at r′ = $r′, foundPointIdx = $foundPointIdx, t = $(foundPointIdx*π/N), boundary point: $(boundaryPoint(b, foundPointIdx*π/N)), using boundary value $(ψ_Γ[foundPointIdx+1])"
    ψ_Γ[foundPointIdx+1]
  end

  # fieldVal
  converted = Base.convert(Complex{Float32}, fieldVal)
end

function exportTotalField(b::Boundary, computedBoundary::Array{Complex{Float64},1}, k::Complex{Float64};
    outFile::String = nothing, xRange=(-1.5,1.5), yRange=(-1.5,1.5), dx=.01, dy=.01, boundaryPoints = 200)
  @info "Calculating field in range ($xRange, $yRange) dx=$dx, dy=$dy..."
  ys = yRange[1]:dy:yRange[2]
  xs = xRange[1]:dx:xRange[2]
  #field = Array(Complex{Float64},(length(ys),length(xs)))
  N = Base.convert(Int,length(computedBoundary)/4)
  ψ_Γ = computedBoundary[1:2N]
  ∂νψ_Γ = computedBoundary[2N+1:end]
  dr = sqrt(dx^2+dy^2)
  discreteB = DiscretizedBoundary(b, boundaryPoints)

  dsList = [ds(b, j, N) for j in 0:2N-1]
  fieldRows = []
  @showprogress 1 for i in 1:length(ys)
    row = @distributed (hcat) for j in 1:length(xs)
      r′ = Point2D(xs[j], ys[i])
      inside = isInside(discreteB, r′)
      n = inside ? nInt(b) : nExt(b)
      ψField(b, N, dr, r′, n, k, ψ_Γ, ∂νψ_Γ, dsList)
    end
    push!(fieldRows, row)
  end
  field = vcat(fieldRows...)

  if outFile != nothing
    JLD.jldopen(outFile, "w") do file
      write(file, "meta/k/full", k)
      write(file, "meta/k/Re", real(k))
      write(file, "meta/k/Im", imag(k))
      write(file, "meta/xRange", xRange)
      write(file, "meta/yRange", yRange)
      write(file, "meta/boundaryObject", b)
      write(file, "boundary/full", computedBoundary)
      write(file, "boundary/Re", real(computedBoundary))
      write(file, "boundary/Im", imag(computedBoundary))
      write(file, "field/full", field)
      write(file, "field/Re", real(field))
      write(file, "field/Im", imag(field))
    end
  end
  field
end


text:
@testset "LinearSolvers" begin
    @testset "testIterativeWrapper.jl" begin
        include("testIterativeWrapper.jl")
    end
    @testset "testBlockIterativeWrapper.jl" begin
        include("testBlockIterativeWrapper.jl")
    end
    @testset "testJuliaWrapper.jl" begin
        include("testJuliaWrapper.jl")
    end
end


text:
struct Card
    suit::Int64
    rank::Int64
    function Card(suit::Int64, rank::Int64)
        @assert(1 <= suit <= 4, "suite is not between 1 and 4")
        @assert(1 <= rank <= 13, "rank is not between 1 and 13")
        new(suit, rank)
    end
end

# Global Variables
const suit_names = ["♣", "♦", "♥", "♠"]
const rank_names = ["A", "2", "3", "4", "5", "6", "7", "8", "9", "10", "J", "Q", "K"]

function Base.show(io::IO, card::Card)
    print(io, rank_names[card.rank], suit_names[card.suit])
end

import Base.< # Must use import
# this won't work:
# function base.<(...)
# end
function <(c1::Card, c2::Card)
    # suit is compared first
    # than rank
    (c1.suit, c1.rank) < (c2.suit, c2.rank)
end

# Abstract Types and Subtyping
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_abstract_types_and_subtyping
abstract type CardSet # <: inheritance
end

struct Deck <: CardSet
    cards :: Array{Card, 1} # 1-D array of Card
end

function Deck()
    deck = Deck(Card[])
    for suit = 1:4
        for rank = 1:13
            push!(deck.cards, Card(suit, rank))
        end
    end

    deck
end

function Base.show(io::IO, cs::CardSet)
    for card in cs.cards
        print(io, card, " ")
    end
    println()
end

function Base.pop!(cs::CardSet)
    pop!(cs.cards)
end

function Base.push!(cs::CardSet, card::Card)
    push!(cs.cards, card)
    nothing
end

import Random: shuffle!
function shuffle!(deck::Deck)
    shuffle!(deck.cards)
    deck
end

# Exercise 18-2
function Base.isless(x::Card, y::Card)
    x < y
end

function Base.sort!(deck::Deck)
    sort!(deck.cards)
    deck
end

# Abstract Types and Subtyping
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_abstract_types_and_subtyping
struct Hand <: CardSet
    cards :: Array{Card, 1}
    label :: String
end

function Hand(label::String="")
    Hand(Card[], label)
end

# Abstract Types and Functions
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_abstract_types_and_functions
"""
    move!(source::CardSet, target::CardSet, n::Int)

Move `n` number of cards from `source` to `target`.
"""
function move!(source::CardSet, target::CardSet, n::Int)
    @assert 1<= n <= length(source.cards)

    for i = 1:n
        card = pop!(source)
        push!(target, card)
    end

    nothing
end


text:
using Documenter, ContinuousWavelets
ENV["PLOTS_TEST"] = "true"
ENV["GKSwstype"] = "100"
makedocs(sitename = "ContinuousWavelets.jl",
    format = Documenter.HTML(),
    authors = "David Weber",
    clean = true,
    pages = [
        "basic usage" => "index.md",
        "Install" => "installation.md",
        "CWT" => [
            "Available Wavelet Families" => "coreType.md",
            "CWT Construction" => "CWTConstruction.md",
            "Wavelet Frequency Spacing" => "spacing.md",
            "Boundary Conditions" => "bound.md",
            "Inversion" => "inverse.md"
        ],
    ])

deploydocs(
    repo = "github.com/dsweber2/ContinuousWavelets.jl.git",
)


text:
using TestFunctionRunner
TestFunctionRunner.@run(paths = ["../lib/TryExperimental"])


text:
using Distributed
using CSV
using DataFrames
using BenchmarkTools
using Plots
using FFTW
using Dates

addprocs(16)
include("1DPackage.jl")
#this package is not installed make sure you are in right working dic
#using .OneDim
@everywhere using .OneDim
@everywhere using SharedArrays

function main(
    lambda,
    mu,
    g,
    T,
    Nt,
    times_saved,
    dt = 0.1,
    L = 501,
    N = 511,
    re = 1.25,
    x0 = 0.0,
)
    M = 1.0
    #for spectral galerkin chose uneven grid size (symmetric k vector)

    q_array = [0.1 * i for i = 1:300]

    D = 1.0

    xstep = L / N
    xarray = [-L / 2 + i * xstep for i = 1:N]
    kstep = 2.0 * pi / L
    k = vcat(0:(N/2), -(N + 1)/2+1:-1) * kstep
    #for phi in equilibrium at t = 0 (note r has to finite if field is non interacting)
    # only valid for r >0
    phi = sampler((N,)) .* sqrt.(T ./ (mu .+ k .^ 2))
    phi[1] = 0
    phi = fft(phi) / N

    parameter_phi = [lambda, mu, g, k, xstep, 2 * D * T, D]
    parameter_impurity = [M, lambda, xstep]
    parameter_V = [re, L]
    x = [x0]
    x_dot = [0.0]

    o = solver(
        Nt,
        dt,
        L,
        N,
        phi,
        x,
        [0.0],
        parameter_phi,
        parameter_impurity,
        parameter_V,
        phi_galerk_modelA,
        NewtonImp,
        V_fct,
        times_saved,
    )

    Ẏ = [i[1] for i in o[3]]
    Y = [i[1] for i in o[2]]
    Ẏ_iẎ_j = Ẏ * Ẏ'
    eq = (1im * (q_array * Y'))
    return Ẏ_iẎ_j, eq, o
end

function average(
    lambda,
    mu,
    g,
    T,
    Nt,
    times_saved,
    dt,
    N_monte_carlo = 300,
    L = 511,
    N = 511,
)
    ans = main(lambda, mu, g, T, Nt, times_saved, dt, L, N)
    Ẏ_iẎ_j = 1.0 / N * ans[1]
    eq = 1.0 / N * ans[2]
    for _ = 1:N_monte_carlo-1
        ans = main(lambda, mu, g, T, Nt, times_saved, dt, L, N)
        Ẏ_iẎ_j += 1.0 / N_monte_carlo * ans[1]
        eq += 1.0 / N_monte_carlo * ans[2]
    end
    folder = string(Dates.today()) * "lambda-$lambda"
    mkdir(folder)
    cd(folder)
    df_y = DataFrame(Ẏ_iẎ_j, :auto)
    CSV.write("y_sy_t.csv", df_y)
    df_para = DataFrame(
        lambda = lambda,
        mu = mu,
        g = g,
        T = T,
        Nt = Nt,
        dt = dt,
        N_monte_carlo = N_monte_carlo,
        N = N,
        L = L,
    )
    CSV.write("parameters.csv", df_para)
    df_eq = DataFrame(eq, :auto)
    CSV.write("eq.csv", df_eq)
    cd(pwd, "..")
end

max_stars_repo_name: phyjonas/ImpuritySGPE

max_issues_repo_name: phyjonas/ImpuritySGPE

max_forks_repo_name: phyjonas/ImpuritySGPE


text:
"""
    bin(data<:Vector{<AbstractFloat}, ticks:Vector{AbstractFloat})

Makes low-to-high histogram of `data` split into bins demarked by `ticks` spanning the range of `data`. Returns the weight of each bin (counts of data per bin), a list of data per bin, and the indices into input `data` set for each bin.
"""
function bin(data::Vector{<:AbstractFloat}, ticks::Vector{<:AbstractFloat})
    sortI = sortperm(data)
    sort!(ticks)
    n = length(ticks)

    weights = zeros(n-1)
    binvals = fill([], n-1)
    bindices = fill([], n-1)
    for i = 1:n-1
        inbin = findall((data[sortI] .≥ ticks[i]) .& (data[sortI] .< ticks[i + 1]))
        weights[i] = length(inbin)
        bindices[i] = sortI[inbin]  # index into original data set, not ordered set
        binvals[i] = data[inbin]
    end

    return (weights, binvals, bindices)
end

bin(data::Vector{<:AbstractFloat}, ticks::LinRange) = bin(data, Vector(ticks))
bin(data::Vector{<:AbstractFloat}, n::Integer) = bin(data, LinRange(min(data...), max(data...), n))