text:
#
# Return the coordinates of the CA atoms only in a vector
#

function xCA( atoms :: Vector{Atom} )

  nCA = count( atom -> atom.name == "CA", atoms )
  x = Matrix{Float64}(undef,nCA,3)
  i = 0
  for atom in atoms
    if atom.name == "CA"
      i = i + 1
      x[i,1] = atom.x
      x[i,2] = atom.y
      x[i,3] = atom.z
    end
  end

  return x

end

max_stars_repo_name: mcubeg/JPleXL

max_issues_repo_name: mcubeg/JPleXL

max_forks_repo_name: mcubeg/JPleXL


text:
function process_likelihood_model(ρ_list::Vector,Eₘ_list::Vector)
    # Generate the A matrix used to calculate likelihoods
    # The A matrix depends on the input states and measurement operators
    dimsmatch(ρ_list,Eₘ_list)
    sum(abs,data(sum(Eₘ_list))-I)<1E-15 ||
        throw(ArgumentError("Eₘ operators do not form a valid POVM!"))
    prep = (Eₘ,ρ) -> transpose(vec(full(ρ⊗transpose(Eₘ))))
    return reduce(vcat,vec(prep.(Eₘ_list,permutedims(ρ_list))))
end

# below is the actual projected gradient descent algorithm from Knee, et al.
# Quantum process tomography via completely positive and trace-preserving
# projection. Phys. Rev. A 98, 062336 (2018).

function pgd_process_tomo(M::Matrix, A::Matrix; tol=1E-10, cptp_tol=1E-8, info=false, cbfun=nothing)
    # Choi process matrix reconstruction by maximum likelihood projected gradient descent
    size(A,1)==length(M) || throw(ArgumentError("A matrix inconsistent with number of measurements!"))
    abs(sum(M)-1)<1/4 || throw(ArgumentError("measurement counts not normalized!"))
    # infer space dimensions from A matrix
    d = isqrt(isqrt(size(A,2)))
    # initial Choi matrix guess, the identity map
    C = Matrix{ComplexF64}(I/d,d^2,d^2)
    # objective and gradient functions setup
    f = C -> loglikelihood(M,C,A)
    ∇f = C -> loglikelihood_gradient(M,C,A)
    # metaparameters & initial cost calculation
    μ = 3/2d^2; γ = 0.3
    c₁ = 1E6; c₂ = f(C)
    info && println("start cost = $c₂")
    # iterate through projected gradient descent steps, with backtracking
    h = CPTP_helpers(C)
    stp = 0
    while c₁ - c₂ > tol
        stp += 1
        c₁, ∇C = c₂, ∇f(C)
        D = project_CPTP(C .- 1/μ.*∇C, h, cptp_tol) - C
        α = 1.0; Π = γ*real(D⋅∇C)
        while (c₂ = f(C .+ α.*D)) > c₁ + α*Π
            α = α/2 # backtrack
        end
        C .= C .+ α.*D
        if cbfun !== nothing # run callback function; e.g. to calc fidelity at each step
            cbfun(stp,c₂,C)
        end
    end
    info && println("final cost = $c₂, Δ = $(c₁-c₂), number of steps: $stp")
    return C
end

function loglikelihood(M::Matrix{T}, C::Matrix, A::Matrix) where {T<:Real}
    # Binomial statistics for the measurement count probability, up to some irrelevant constant
    logP = log.(max.(real.(A*vec(C)), 1E-16))
    return -real(vec(M)⋅logP)
end

function loglikelihood_gradient(M::Matrix, C::Matrix, A::Matrix)
    P = max.(real.(A*vec(C)), 1E-16)
    return unvec(-A'*(vec(M)./P))
end

function project_CPTP(C::Matrix, h::Tuple, tol::Real=1E-8)
    # generate storage objects
    X₁ = copy(vec(C)); Y₁ = zero(X₁);
    X₂ = copy(Y₁); Y₂ = copy(Y₁)
    P = copy(Y₁); Q = copy(Y₁)
    ΔP = 1.0; ΔQ = 1.0
    D, V, Mdagvec𝕀, MdagM = h
    # iterate through TP & CP projections
    while ΔP^2 + ΔQ^2 + 2*abs(P⋅X₂-P⋅X₁) + 2*abs(Q⋅Y₂-Q⋅Y₁) > tol
        project_TP!(Y₂, X₁+P, Mdagvec𝕀, MdagM)
        ΔP = norm2_diff(X₁,Y₂)
        P .= X₁ .- Y₂ .+ P
        project_CP!(X₂, Y₂+Q, D, V)
        ΔQ = norm2_diff(Y₂,X₂)
        Q .= Y₂ .- X₂ .+ Q
        X₁, X₂ = X₂, X₁
        Y₁, Y₂ = Y₂, Y₁
    end
    return unvec(X₁)
end

function project_CP!(X::Vector, vecC::Vector, D::Vector, V::Matrix)
    # Project the process onto the completely positive subspace by making the
    # Choi matrix positive semidefinite
    # We do this by taking the eigendecomposition, setting any negative
    # eigenvalues to 0, and reconstructing the Choi matrix
    C = unvec(vecC)
    hermitianize!(C)
    hermfact!(D,V,Hermitian(C))
    D .= max.(D,0)
    mul!(C,Diagonal(D),V')
    mul!(unvec(X),V,C)
end

function project_TP!(Y::Vector, vecC::Vector, Mdagvec𝕀::Vector, MdagM::Matrix)
    # Project the process onto the trace-preserving subspace
    d⁻¹ = 1/isqrt(isqrt(length(vecC)))
    mul!(Y, MdagM, vecC)
    Y .= vecC .- d⁻¹.*Y .+ d⁻¹.*Mdagvec𝕀
end

function CPTP_helpers(C::Matrix)
    D = Vector{real(eltype(C))}(undef,size(C,1))
    V = Matrix{eltype(C)}(undef,size(C))
    Mdagvec𝕀, MdagM = TP_helper_matrices(C)
    return D, V, Mdagvec𝕀, MdagM
end

function TP_helper_matrices(C::Matrix)
    d = isqrt(size(C,1))
    𝕀 = Matrix(1.0I,d,d); k = zeros(1,d)
    # this can be done more efficiently, but prob doesn't matter
    M = sum(i->(k[i]=1; k[mod1(i-1,d)]=0; (𝕀 ⊗ k) ⊗ (𝕀 ⊗ k)), 1:d)
    return M'*vec(𝕀), M'*M
end

ext: jl


text:
using Documenter
using EzXML

makedocs(
    sitename="EzXML.jl",
    modules=[EzXML],
    pages=["index.md", "manual.md", "reference.md", "devnotes.md"])

deploydocs(
    repo="github.com/bicycle1885/EzXML.jl.git",
    target="build")


text:
using Genie, MbedTLS

function configure_dev_ssl()
  cert = Genie.Assets.embedded_path(joinpath("files", "ssl", "localhost.crt")) |> MbedTLS.crt_parse_file
  key = Genie.Assets.embedded_path(joinpath("files", "ssl", "localhost.key")) |> MbedTLS.parse_keyfile

  ssl_config = MbedTLS.SSLConfig(true)
  entropy = MbedTLS.Entropy()
  rng = MbedTLS.CtrDrbg()
  MbedTLS.config_defaults!(ssl_config, endpoint=MbedTLS.MBEDTLS_SSL_IS_SERVER)

  MbedTLS.authmode!(ssl_config, MbedTLS.MBEDTLS_SSL_VERIFY_NONE)
  MbedTLS.seed!(rng, entropy)
  MbedTLS.rng!(ssl_config, rng)
  MbedTLS.own_cert!(ssl_config, cert, key)
  MbedTLS.ca_chain!(ssl_config)

  Genie.config.ssl_config = ssl_config

  nothing
end

Genie.Configuration.isdev() && Genie.config.ssl_enabled && configure_dev_ssl()


text:
using Test

@testset "test/unit/Nonlocality.jl" begin

using BellScenario

@testset "running all tests for module Nonlocality.jl" begin

  println("running Nonlocality.jl unit tests:")
  for test in readdir("./test/unit/Nonlocality/")
      # run only julia files in test directory
      if occursin(r"^.*\.jl$", test)
          println("./unit/Nonlocality/$test")
          @time include("./Nonlocality/$test")
      end
  end

end

end


text:
module JulioTests

include("utils.jl")
include("test_stack.jl")
include("test_taskgroups.jl")
include("test_channels.jl")
include("test_promise.jl")
include("test_io.jl")
include("test_locks.jl")
include("test_select.jl")
include("test_cancellation.jl")
include("test_introduction.jl")
include("test_structured_concurrency.jl")
include("test_search3.jl")
include("test_happy_eyeballs.jl")
include("test_custom_select.jl")
include("test_contextmanagers.jl")
include("test_dying_philosopher.jl")
include("test_doctest.jl")

end  # module JulioTests


text:
################################################################################
#
#          NfOrd/Elem.jl : Elements of orders of number fields
#
# This file is part of hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016 Tommy Hofmann
#
################################################################################

export ==, +, -, *, ^, add!, conjugates_arb, conjugates_arb_log, discriminant,
       divexact, elem_in_nf, coordinates, isone, iszero, minkowski_map, mod,
       mul!, norm, one, parent, powermod, rand, rand!, representation_matrix,
       show, trace, t2, zero

################################################################################
#
#  Deepcopy
#
################################################################################

function Base.deepcopy_internal(x::NfAbsOrdElem{S, T}, dict::IdDict) where {S, T}
  z = parent(x)()
  z.elem_in_nf = Base.deepcopy_internal(x.elem_in_nf, dict)
  if x.has_coord
    z.has_coord = true
    z.coordinates = Base.deepcopy_internal(x.coordinates, dict)::Vector{fmpz}
  end
  return z
end

################################################################################
#
#  Conversion from matrix
#
################################################################################

function elem_from_mat_row(O::NfAbsOrd, M::fmpz_mat, i::Int, d::fmpz = fmpz(1))
  return O(fmpz[M[i, j] for j=1:degree(O)])
end

################################################################################
#
#  Parent object overloading
#
################################################################################

@doc Markdown.doc"""
      (O::NumFieldOrd)(a::NumFieldElem, check::Bool = true) -> NumFieldOrdElem

Given an element $a$ of the ambient number field of $\mathcal O$, this
function coerces the element into $\mathcal O$. It will be checked that $a$
is contained in $\mathcal O$ if and only if `check` is `true`.
"""
(O::NfAbsOrd{S, T})(a::T, check::Bool = true) where {S, T} = begin
  if nf(O) !== parent(a)
    throw(error("Underlying number fields not equal"))
  end
  if check
    (x, y) = _check_elem_in_order(a,O)
    !x && error("Number field element not in the order")
    return NfAbsOrdElem(O, deepcopy(a), y)
  else
    return NfAbsOrdElem(O, deepcopy(a))
  end
end

@doc Markdown.doc"""
      (O::NumFieldOrd)(a::NumFieldOrdElem, check::Bool = true) -> NumFieldOrdElem

Given an element $a$ of some order in the ambient number field of
$\mathcal O$, this function coerces the element into $\mathcal O$. It
will be checked that $a$ is contained in $\mathcal O$ if and only if
`check` is `true`.
"""
(O::NfAbsOrd{S, T})(a::NfAbsOrdElem{S, T}, check::Bool = true) where {S, T} = begin
  b = nf(parent(a))(a)
  if check
    (x, y) = _check_elem_in_order(b,O)
    !x && error("Number field element not in the order")
    return NfAbsOrdElem(O, deepcopy(b), y)
  else
    return NfAbsOrdElem(O, deepcopy(b))
  end
end

(O::NfAbsOrd{S, T})(a::T, arr::Vector{fmpz}, check::Bool = false) where {S, T} = begin
  if check
    (x, y) = _check_elem_in_order(a,O)
    (!x || arr != y ) && error("Number field element not in the order")
    return NfAbsOrdElem(O, deepcopy(a), y)
  else
    return NfAbsOrdElem(O, deepcopy(a), deepcopy(arr))
  end
end

(O::NfAbsOrd{S, T})(a::T, arr::fmpz_mat, check::Bool = false) where {S, T} = begin
  if check
    (x, y) = _check_elem_in_order(a,O)
    (!x || arr != y ) && error("Number field element not in the order")
    return NfAbsOrdElem(O, deepcopy(a), y)
  else
    return NfAbsOrdElem(O, deepcopy(a), deepcopy(arr))
  end
end

@doc Markdown.doc"""
      (O::NumFieldOrd)(a::IntegerUnion) -> NumFieldOrdElem

Given an element $a$ of type `fmpz` or `Integer`, this
function coerces the element into $\mathcal O$.
"""
(O::NfAbsOrd)(a::IntegerUnion) = begin
  return NfAbsOrdElem(O, nf(O)(a))
end

@doc Markdown.doc"""
      (O::NfAbsOrd)(arr::Vector{fmpz})

Returns the element of $\mathcal O$ with coefficient vector `arr`.
"""
(O::NfAbsOrd)(arr::Vector{fmpz}) = begin
  return NfAbsOrdElem(O, deepcopy(arr))
end

(O::NfAbsOrd)(arr::fmpz_mat) = begin
  return NfAbsOrdElem(O, arr)
end

@doc Markdown.doc"""
      (O::NfAbsOrd)(arr::Vector{Integer})

Returns the element of $\mathcal O$ with coefficient vector `arr`.
"""
(O::NfAbsOrd)(arr::Vector{S}) where {S <: Integer} = begin
  return NfAbsOrdElem(O, deepcopy(arr))
end

(O::NfAbsOrd)() = NfAbsOrdElem(O)

################################################################################
#
#  Parent check
#
################################################################################

function check_parent(x::NfAbsOrdElem{S, T}, y::NfAbsOrdElem{S, T}) where {S, T}
  return parent(x) === parent(y)
end

################################################################################
#
#  "Assure" functions for fields
#
################################################################################

function assure_has_coord(a::NfAbsOrdElem)
  if a.has_coord
    return nothing
  else
    (x, y) = _check_elem_in_order(a.elem_in_nf, parent(a))
    !x && error("Not a valid order element")
    a.coordinates = y
    a.has_coord = true
    return nothing
  end
end

################################################################################
#
#  Coordinates
#
################################################################################

@doc Markdown.doc"""
    coordinates(a::NfAbsOrdElem) -> Vector{fmpz}

Returns the coefficient vector of $a$ with respect to the basis of the order.
"""
function coordinates(a::NfAbsOrdElem; copy::Bool = true)
  assure_has_coord(a)
  @hassert :NfOrd 2 a == dot(a.coordinates, basis(parent(a), copy = false))
  if copy
    return deepcopy(a.coordinates)
  else
    return a.coordinates
  end
end

################################################################################
#
#  Characteristic and minimal polynomial
#
################################################################################

@doc Markdown.doc"""
    charpoly(a::NfAbsOrdElem) -> fmpz_poly
    charpoly(a::NfAbsOrdElem, FlintZZ) -> fmpz_poly

The characteristic polynomial of $a$.
"""
function charpoly(a::NfAbsOrdElem, Zx::FmpzPolyRing = FmpzPolyRing(FlintZZ, :x, false))
  return Zx(charpoly(elem_in_nf(a)))
end

@doc Markdown.doc"""
    minpoly(a::NfAbsOrdElem) -> fmpz_poly

The minimal polynomial of $a$.
"""
function minpoly(a::NfAbsOrdElem, Zx::FmpzPolyRing = FmpzPolyRing(FlintZZ, :x, false))
  return Zx(minpoly(elem_in_nf(a)))
end

################################################################################
#
#  String I/O
#
################################################################################

function AbstractAlgebra.expressify(a::NfAbsOrdElem; context = nothing)
  return AbstractAlgebra.expressify(a.elem_in_nf, context = context)
end

function show(io::IO, x::NfAbsOrdElem)
  print(io, AbstractAlgebra.obj_to_string(x, context = io))
end

################################################################################
#
#  Adhoc division
#
################################################################################

function //(x::nf_elem, y::NfOrdElem)
  check_parent(x, y.elem_in_nf)
  return x//y.elem_in_nf
end

function //(y::NfOrdElem, x::nf_elem)
  check_parent(x, y.elem_in_nf)
  return y.elem_in_nf//x
end

################################################################################
#
#  Modular reduction
#
################################################################################

function mod(a::NfAbsOrdElem, m::Union{fmpz, Int})
  d = degree(parent(a))
  ar = coordinates(a)
  for i in 1:d
    ar[i] = mod(a.coordinates[i], m)
  end
  return NfAbsOrdElem(parent(a), ar) # avoid making a copy of ar
end

################################################################################
#
#  Modular exponentiation
#
################################################################################

function powermod(a::NfAbsOrdElem, i::fmpz, p::fmpz)

  #if contains_equation_order(parent(a))#This doesn't work!
  if is_defining_polynomial_nice(nf(parent(a)))
    return powermod_fast(a, i, p)
  else
    return powermod_gen(a, i, p)
  end
end

function powermod_gen(a::NfAbsOrdElem, i::fmpz, p::fmpz)
  if i == 0
    return one(parent(a))
  end
  if i == 1
    b = mod(a,p)
    return b
  end
  if mod(i,2) == 0
    j = div(i, 2)
    b = powermod(a, j, p)
    b = b^2
    b = mod(b,p)
    return b
  end
  b = mod(a * powermod(a, i - 1, p), p)
  return b
end

function powermod_fast(a::NfAbsOrdElem, i::fmpz, p::fmpz)
  if i == 0
    return one(parent(a))
  end
  if i == 1
    bb = mod(a, p)
    return bb
  end

  b = a.elem_in_nf
  d = denominator(b)
  f, e = ppio(d, p) # f = p^k, and e is a unit mod p
  b *= e
  e = invmod(e, p)
  e = powermod(e, i, p)

  y = one(parent(b))
  while i > 1
    if iszero(b)
      return zero(parent(a))
    end
    if iseven(i)
      b = mod(b*b, p)
    else
      y = mod(b*y, p)
      b = mod(b*b, p)
    end
    i = div(i, 2)
  end
  b = mod(b*y, p)

  return mod(parent(a)(b)*e, p)
end

function powermod_fast(a::NfAbsOrdElem{NfAbsNS, NfAbsNSElem}, i::fmpz, p::fmpz)
  if i == 0
    return one(parent(a))
  end
  if i == 1
    bb = mod(a, p)
    return bb
  end

  b = a.elem_in_nf
  d = denominator(b)
  f, e = ppio(d, p) # f = p^k, and e is a unit mod p
  b *= e
  e = invmod(e, p)
  e = powermod(e, i, p)

  y = one(parent(b))
  while i > 1
    if iszero(b)
      return zero(parent(a))
    end
    if iseven(i)
      mul!(b, b, b)
      mod!(b, p)
      #b = mod(b*b, p)
    else
      mul!(y, b, y)
      mod!(y, p)
      #y = mod(b*y, p)
      mul!(b, b, b)
      mod!(b, p)
      #b = mod(b*b, p)
    end
    i = div(i, 2)
  end
  mul!(b, b, y)
  mod!(b, p)
  #b = mod(b*y, p)

  return mod(parent(a)(b*e), p)
end

function powermod(a::NfOrdElem, i::fmpz, I::NfOrdIdl)
  if i == 0
    return one(parent(a))
  end
  b = mod(a, I)
  if i == 1
    return b
  end
  if iszero(b)
    return b
  end
  y = one(parent(b))
  while i > 1
    if iseven(i)
      b = mod(b*b, I)
    else
      y = mod(b*y, I)
      b = mod(b*b, I)
    end
    i = div(i, 2)
  end
  b = mod(b*y, I)

  return mod(b, I)
end

powermod(a::NfAbsOrdElem, i::Integer, m::Integer) = powermod(a, fmpz(i), fmpz(m))

powermod(a::NfAbsOrdElem, i::fmpz, m::Integer)  = powermod(a, i, fmpz(m))

powermod(a::NfAbsOrdElem, i::Integer, m::fmpz)  = powermod(a, fmpz(i), m)

################################################################################
#
#  Representation matrices
#
################################################################################

@doc Markdown.doc"""
    representation_matrix(a::NfAbsOrdElem) -> fmpz_mat

Returns the representation matrix of the element $a$.
"""
function representation_matrix(a::NfAbsOrdElem)
  O = parent(a)
  assure_has_basis_matrix(O)
  assure_has_basis_mat_inv(O)
  A = representation_matrix(a, nf(O))
  mul!(A, O.basis_matrix, A)
  mul!(A, A, O.basis_mat_inv)
  !(A.den == 1) && error("Element not in order")
  return A.num
end

@doc Markdown.doc"""
    representation_matrix(a::NfAbsOrdElem, K::AnticNumberField) -> FakeFmpqMat

Returns the representation matrix of the element $a$ considered as an element
of the ambient number field $K$. It is assumed that $K$ is the ambient number
field of the order of $a$.
"""
function representation_matrix(a::NfAbsOrdElem{S, T}, K::S) where {S, T}
  nf(parent(a)) != K && error("Element not in this field")
  A, d = Nemo.representation_matrix_q(a.elem_in_nf)
  A.base_ring = FlintZZ
  z = FakeFmpqMat(A, d)
  return z
end

@doc Markdown.doc"""
    representation_matrix_mod(a::NfAbsOrdElem, d::fmpz) -> fmpz_mat

Returns the representation matrix of the element $a$ with entries reduced mod $d$.
"""
function representation_matrix_mod(a::NfAbsOrdElem, d::fmpz)
  O = parent(a)
  A, den = representation_matrix_q(elem_in_nf(a))
  BM = basis_matrix(O, copy = false)
  BMinv = basis_mat_inv(O, copy = false)
  d2 = BM.den * BMinv.den * den
  d2c, d2nc = ppio(d2, d)
  d1 = d * d2c
  if fits(Int, d1)
    R = ResidueRing(FlintZZ, Int(d1), cached = false)
    AR = map_entries(R, A)
    BMR = map_entries(R, BM.num)
    BMinvR = map_entries(R, BMinv.num)
    mul!(AR, BMR, AR)
    mul!(AR, AR, BMinvR)
    inver = invmod(d2nc, d1)
    mul!(AR, AR, R(inver))
    dv = R(d2c)
    for i = 1:nrows(AR)
      for j = 1:ncols(AR)
        if !iszero(AR[i, j])
          AR[i, j] = divexact(AR[i, j], dv)
        end
      end
    end
    res = lift!(A, AR)
    mod!(res, d)
    return res
  else
    RR = ResidueRing(FlintZZ, d1, cached = false)
    ARR = map_entries(RR, A)
    BMRR = map_entries(RR, BM.num)
    mul!(ARR, BMRR, ARR)
    dv = RR(d2c)
    BMinvRR = map_entries(RR, BMinv.num)
    mul!(ARR, ARR, BMinvRR)
    inver = invmod(d2nc, d1)
    ARR *= inver
    for i = 1:nrows(ARR)
      for j = 1:ncols(ARR)
        if !iszero(ARR[i, j])
          ARR[i, j] = divexact(ARR[i, j], dv)
        end
      end
    end
    res1 = lift(ARR)
    mod!(res1, d)
    return res1
  end
end

################################################################################
#
#  Random element generation
#
################################################################################

# TODO: Make this faster, don't allocate the ar array ...
function rand!(z::NfAbsOrdElem{S, T}, B::Vector{NfAbsOrdElem{S, T}}, R) where {S, T}
  O = parent(z)
  y = O()
  for i in 1:degree(O)
    y = mul!(y, rand(R), B[i])
    z = add!(z, z, y)
  end
  return z
end

function rand!(z::NfAbsOrdElem, O::NfAbsOrd, R::UnitRange{T}) where T <: Integer
  y = O()
  B = basis(O, copy = false)
  for i in 1:degree(O)
    y = mul!(y, rand(R), B[i])
    z = add!(z, z, y)
  end
  return z
end

@doc Markdown.doc"""
    rand(O::NfOrd, R::UnitRange{Integer}) -> NfAbsOrdElem

Computes a coefficient vector with entries uniformly distributed in `R` and returns
the corresponding element of the order.
"""
function rand(O::NfOrd, R::UnitRange{T}) where T <: Integer
  z = O()
  rand!(z, O, R)
  return z
end

function rand!(z::NfAbsOrdElem, O::NfOrd, n::IntegerUnion)
  return rand!(z, O, -n:n)
end

@doc Markdown.doc"""
    rand(O::NfOrd, n::IntegerUnion) -> NfAbsOrdElem

Computes a coefficient vector with entries uniformly distributed in
$\{-n,\dotsc,-1,0,1,\dotsc,n\}$ and returns the corresponding element of the
order $\mathcal O$.
"""
function rand(O::NfOrd, n::Integer)
  return rand(O, -n:n)
end

function rand(O::NfOrd, n::fmpz)
  z = O()
  rand!(z, O, BigInt(n))
  return z
end

function rand!(z::NfAbsOrdElem, O::NfOrd, n::fmpz)
  return rand!(z, O, BigInt(n))
end

################################################################################
#
#  Conversion
#
################################################################################

(K::AnticNumberField)(x::NfAbsOrdElem{AnticNumberField, nf_elem}) = elem_in_nf(x)

(K::NfAbsNS)(x::NfAbsOrdElem{NfAbsNS, NfAbsNSElem}) = elem_in_nf(x)

################################################################################
#
#  Factorization
#
################################################################################

@doc Markdown.doc"""
    factor(a::NfOrdElem) -> Fac{NfOrdElem}

Computes a factorization of $a$ into irreducible elements. The return value
is a factorization `fac`, which satisfies `a = unit(fac) * prod(p^e for (p, e)
in fac)`.

The function requires that $a$ is non-zero and that all prime ideals containing
$a$ are principal, which is for example satisfied if class group of the order
of $a$ is trivial.
"""
function factor(a::NfOrdElem)
  iszero(a) && throw(error("Element must be non-zero"))
  OK = parent(a)
  I = a * OK
  D = Dict{NfOrdElem, Int}()
  u = a
  for (p, e) in factor(I)
    b, c = is_principal(p)
    !b && throw(error("Prime ideal dividing the element not principal"))
    D[c] = e
    u = divexact(u, c^e)
  end
  return Nemo.Fac(u, D)
end


text:
module SeisDvv

using DSP, GLM, DataFrames, Statistics, FFTW, Interpolations, SeisNoise

include("tools.jl")
include("WCC.jl")
include("Stretching.jl")
include("DTW.jl")
include("MWCS.jl")
include("Wavelets.jl")

end


text:
using Test
using Pluto
import Pluto: samplenotebook, Notebook, Client, Cell, load_notebook, load_notebook_nobackup, save_notebook, run_reactive!, WorkspaceManager

function notebook_inputs_equal(nbA, nbB)
    x = normpath(nbA.path) == normpath(nbB.path)

    to_compare(cell) = (cell.uuid, cell.code)
    y = to_compare.(nbA.cells) == to_compare.(nbB.cells)

    x && y
end

"The converse of Julia's `Base.sprint`."
function sread(f::Function, input::String, args...)
    io = IOBuffer(input)
    output = f(io, args...)
    close(io)
    return output
end

@testset "Notebook File I/O" begin
    notebookA = samplenotebook()
    notebookB = Notebook(joinpath(tempdir(), "é🧡💛.jl"), [
        Cell("z = y"),
        Cell("v = u"),
        Cell("y = x"),
        Cell("x = w"),
        Cell("using Base"),
        Cell("t = 1"),
        Cell("w = v"),
        Cell("u = t"),
    ])
    for nb in [notebookA, notebookB]
        client = Client(Symbol("client", rand(UInt16)), nothing)
        client.connected_notebook = nb
        Pluto.connectedclients[client.id] = client
    end

    @testset "I/O" begin
        @test let
            nb = notebookB
            save_notebook(nb)
            result = load_notebook_nobackup(nb.path)
            notebook_inputs_equal(nb, result)
        end
    end
    @testset "Line endings" for nb in [notebookA, notebookB]
        file_contents = sprint(save_notebook, nb)

        @test let
            result = sread(load_notebook_nobackup, file_contents, nb.path)
            notebook_inputs_equal(nb, result)
        end

        @test let
            file_contents_windowsed = replace(file_contents, "\n" => "\r\n")
            result_windowsed = sread(load_notebook_nobackup, file_contents_windowsed, nb.path)
            notebook_inputs_equal(nb, result_windowsed)
        end
    end
    @testset "Bad code" begin
        nb = Notebook(joinpath(tempdir(), "é🧡💛.jl"), [
            Cell("z = y"),
            Cell("y = z")
        ])
        @test let
            file_contents = sprint(save_notebook, nb)
            result = sread(load_notebook_nobackup, sprint(save_notebook, nb), nb.path)
            notebook_inputs_equal(nb, result)
        end
    end
    @testset "Backup" begin
        notebookA.path = tempname() * ".jl"
        notebookB.path = tempname() * ".jl"
        backup_path = notebookB.path * ".backup"
        isfile(backup_path) && rm(backup_path)
        run_reactive!(notebookB, notebookB.cells)
        save_notebook(notebookB)
        notebookC = load_notebook(notebookB.path)
        @test !isfile(backup_path)
        WorkspaceManager.unmake_workspace(notebookB)
    end

    # TODO: tests for things mentioned in https://github.com/fonsp/Pluto.jl/issues/9
    # TODO: test bad dirs, filenames, permissions
end


text:
module YaoHIR

export GenericRoutine, Routine,
    IntrinsicRoutine,
    Operation,
    AdjointOperation,
    Chain, Gate, Ctrl,
    BlockIR

using MLStyle
using Expronicon
using Core: CodeInfo
using Core.Compiler: IRCode

include("types.jl")
include("intrinsic.jl")
include("printing.jl")

end


text:
using Measures
using Statistics
import Images

using FITSIO

"""
    ds9show(images..., [lock=true], [pad=nothing])

Open one or more images in a SAO DS9 window.
By default, the scales and zooms of the images are locked.
Pass lock=false to disable.

If you pass `pad=true`, the images are padded with NaN
so that they have equal axes. By default (pad=nothing),
padding will be applied when all images are 2D and
locked.
"""
function ds9show(
    # Accept any number of images
    imgs...;
    # Flags
    lock=true,
    pad=nothing,
    setscale=nothing,
    setcmap=nothing,
    loadcmap=nothing,
    τ=nothing,
    regions=String[]
)
    # See this link for DS9 Command reference
    # http://ds9.si.edu/doc/ref/command.html#fits

    # If lock is true (deafult) and the images have different sizes (and are 2D), and
    # padding has not been disabled, turn padding on.
    if length(imgs) > 1 && isnothing(pad) && lock && all(==(2), length.(size.(imgs))) && !all(i -> size(i) == size(first(imgs)), imgs)
        @warn "Padding images so that locked axes work correctly. Disable with either `pad=false` or `lock=false`"
        pad = true
    else
        pad = false
    end

    # If pad is set, pad out the images with NaN to be the same size.
    if pad
        # Sort of an ungly line. paddedviews returns a tuple of views.
        # We need to convert to a vector of arrays.
        imgs = [collect.(Images.paddedviews(NaN, imgs...))...]
    end

    # Detect images with complex values, and just show the power
    imgs = map(imgs) do img
        if eltype(img) <: Complex
            return abs.(img)
        end
        if eltype(img) <: Bool
            return Int8.(img)
        end
        return img
    end

    # For each image, write a temporary file.
    # We clean these up after DS9 exits.
    fnames = String[]
    for img in imgs
        tempfile = tempname()
        writefits(tempfile, img)
        push!(fnames, tempfile)
    end

    # Decide where to look to open DS9 depending on the platform
    if Sys.iswindows()
        cmd = `C:\\SAOImageDS9\\ds9.exe $fnames`
    elseif Sys.isapple() && isdir("/Applications/SAOImageDS9.app")
        cmd = `open -W /Applications/SAOImageDS9.app --args $fnames`
    elseif Sys.isapple() && isdir("/Applications/SAOImage DS9.app")
        cmd = `open -W /Applications/SAOImage\ DS9.app --args $fnames`
    else
        @warn "Untested system for ds9show. Assuming ds9 is in PATH." maxlog=1
        cmd = `ds9 $fnames`
    end

    if !isnothing(setscale)
        cmd =  `$cmd -scale mode $setscale`
    end

    if !isnothing(loadcmap)
        cmd =  `$cmd -cmap load $loadcmap`
    elseif !isnothing(setcmap)
        cmd =  `$cmd -cmap $setcmap`
    end

    if !isnothing(τ)
        σ = std(filter(isfinite, view(last(imgs),:,:,1,1,1)))
        l1 = -0.5τ*σ
        l2 = +1τ*σ
        cmd = `$cmd -scale limits $l1 $l2`
    end

    # If lock=true, then add almost all possible lock flags to the command
    if lock
        cmd = `$cmd -lock frame image -lock crosshair image -lock crop image -lock slice image -lock bin yes -lock axes yes -lock scale yes -lock scalelimits yes -lock colorbar yes -lock block yes -lock smooth yes `
    end

    for region in regions
        cmd = `$cmd -regions command "$region"`
    end

    # Open DS9 asyncronously.
    # When it finishes or errors, delete the temporary FITS files.
    task = @async begin
        try
            # Open DS9
            run(cmd)
        finally
            # Then delete files whenever it's done
            println()
            @info "cleaning up tmp. files from DS9" fnames
            rm.(fnames)
        end
        return nothing
    end

    # Return the async task, in case the user wants to run something when they are done.
    return task
end

# If called with a vector of images, expand and call the function above
ds9show(imgs::AbstractArray{<:AbstractArray}; kwargs...) = ds9show(imgs...;kwargs...)
export ds9show

using ImageShow, Images, ColorSchemes

function imshow2(
    img;
    clims=extrema(skipmissing(filter(isfinite, img))),
    τ=nothing,
    cmap=:Grey
)
    # c = convert(Array{Float32}, img)
    c =  img
    # Fast path: pure grayscale image
    if cmap == :Grey
        c = copy(c)
        min, max = clims
        c .-= min
        c ./= max - min
        if any(!isfinite, c)
            c = collect(view(c, :,reverse(axes(c,2)))')
            mask = .! isfinite.(c)
            c[mask] .= 0
            coloured = RGBA.(Gray.(c))
            coloured[mask] .= RGBA(0,0,0,0)
            coloured
        else
            Gray.(c[:,end:-1:begin]')
        end
    else
        cscheme = ColorSchemes.colorschemes[cmap]
        if any(!isfinite, c)
            c = collect(view(c, :,reverse(axes(c,2)))')
            mask = .! isfinite.(c)
            c[mask] .= 0
            coloured = RGBA.(get(cscheme::ColorScheme, c, clims))
            coloured[mask] .= RGBA(0,0,0,0)
            coloured
        else
            get(cscheme::ColorScheme, c[:,end:-1:begin]', clims)
        end
    end

end
export imshow2

# This is a Plots recipe for how to display our image type.
# If the user runs `using Plots; plot(img)` this will be called.
using RecipesBase
@recipe function f(img::DirectImage)

    τ = get(plotattributes, :τ, nothing)
    lims = get(plotattributes, :lims, nothing)

    skyconvention = get(plotattributes, :skyconvention, false)

    clims = nothing
    if !(isnothing(τ) || ismissing(τ) || τ==:none)
        σ = std(filter(isfinite, view(img,:,:,1,1,1)))
        clims = (-0.5τ*σ, +1τ*σ)
    end

    delete!(plotattributes, :add_marker)

    if !isnothing(lims)
        if length(lims) == 1
            lims = (-lims, lims)
        end
        xlims := lims
        ylims := lims
    end

    unit = "px"
    platescale = 1.0
    if hasproperty(img, :PLATESCALE)
        platescale = img.PLATESCALE
        unit = "mas"
    elseif hasproperty(img, :PIXSCALE)
        platescale = img.PIXSCALE*3.6e6
        unit = "mas"
    elseif hasproperty(img, :PLATESC)
        platescale = img.PLATESC*1e3
        unit = "mas"
    elseif hasproperty(img, :CDELT1)
        platescale = img.CDELT1*1e3
        unit = "mas"
    end

    colorbar_title = ""
    if hasproperty(img, :UNIT)
        colorbar_title = img.UNIT
    end

    seriestype := :heatmap
    xguide --> "X - $unit"
    yguide --> "Y - $unit"
    label --> ""
    # grid --> false
    # legend --> false
    # background_color_outside --> "#000000"
    # background_color_inside --> "#000000"
    # foreground_color --> "#eee"
    aspect_ratio --> 1
    clims --> clims
    colorbar_title --> colorbar_title

    if skyconvention
        xflip --> true
        return (
            UnitRange(axes(img,1)).*platescale,
            UnitRange(axes(img,2)).*platescale,
            reverse(collect(transpose(img[:,:,1,1,1])), dims=2)
        )
    else
        return (
                UnitRange(axes(img,1)).*platescale,
                UnitRange(axes(img,2)).*platescale,
                collect(transpose(img[:,:,1,1,1]))
        )
    end
end

"""
imshow(img; τ=5, lims=1200)

^ displays an image with +5 std, -2.5std, xlims and ylims set to 1200 mas.

Display an image with convenient formatting.
Looks into the FITS headers for the platescale, if available.

!! NOTE: You must install the Plots package to use imshow.
Run `using Plots` before `using DirectImages` to enable.
"""
function imshow end

function imshow(img; τ=5, lims=nothing, kwargs... )
    error("The Plots package is not active. Run `using Plots` before `using DirectImages` to enable.")
end
function imshow!(img; τ=5, lims=nothing, kwargs... )
    error("The Plots package is not active. Run `using Plots` before `using DirectImages` to enable.")
end

# Optionally depend on Plots. If the user imports it, this code will be run to set up
# our `imshow` function.
using Requires
function __init__()
    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin

        # Imshow just ensures that the argument is a DirectImage and then
        # dispatches to the Plots recipe above with a default value for τ.
        function imshow(img::DirectImage; τ=5, lims=nothing, kwargs...)
            Plots.plot(img; τ, lims, kwargs...)
        end
        function imshow(img; τ=5, lims=nothing, kwargs...)
            Plots.plot(DirectImage(img); τ, lims, kwargs...)
        end
        function imshow!(img::DirectImage; τ=5, lims=nothing, kwargs...)
            Plots.plot!(img; τ, lims, kwargs...)
        end
        function imshow!(img; τ=5, lims=nothing, kwargs...)
            Plots.plot!(DirectImage(img); τ, lims, kwargs...)
        end

    end

    nothing
end
export imshow
export imshow!

max_stars_repo_name: sefffal/SpaceImages


text:
module DistributionsAD

using PDMats,
      LinearAlgebra,
      Distributions,
      Random,
      SpecialFunctions,
      StatsFuns,
      Compat,
      Requires,
      ZygoteRules,
      ChainRules,  # needed for `ChainRules.chol_blocked_rev`
      ChainRulesCore,
      FillArrays,
      Adapt

using SpecialFunctions: logabsgamma, digamma
using LinearAlgebra: copytri!, AbstractTriangular
using Distributions: AbstractMvLogNormal,
                     ContinuousMultivariateDistribution
using Base.Iterators: drop

import StatsBase
import StatsFuns: logsumexp,
                  binomlogpdf,
                  nbinomlogpdf,
                  poislogpdf,
                  nbetalogpdf
import Distributions: MvNormal,
                      MvLogNormal,
                      logpdf,
                      quantile,
                      PoissonBinomial,
                      Binomial,
                      BetaBinomial,
                      Erlang

export TuringScalMvNormal,
       TuringDiagMvNormal,
       TuringDenseMvNormal,
       TuringMvLogNormal,
       TuringPoissonBinomial,
       TuringWishart,
       TuringInverseWishart,
       arraydist,
       filldist

include("common.jl")
include("arraydist.jl")
include("filldist.jl")
include("univariate.jl")
include("multivariate.jl")
include("matrixvariate.jl")
include("flatten.jl")

include("chainrules.jl")
include("zygote.jl")

@init begin
    @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" begin
        using .ForwardDiff
        using .ForwardDiff: @define_binary_dual_op # Needed for `eval`ing diffrules here
        include("forwarddiff.jl")
    end

    @require ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" begin
        # ensures that we can load ForwardDiff without depending on it
        # (it is a dependency of ReverseDiff and therefore always available)
        @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" begin
            include("reversediff.jl")
        end
    end

    @require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" begin
        using DiffRules
        using SpecialFunctions
        using LinearAlgebra: AbstractTriangular
        using .Tracker: Tracker, TrackedReal, TrackedVector, TrackedMatrix,
                        TrackedArray, TrackedVecOrMat, track, @grad, data
        include("tracker.jl")
    end

    @require LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02" begin
        using .LazyArrays: BroadcastArray, BroadcastVector, LazyArray

        const LazyVectorOfUnivariate{
            S<:ValueSupport,
            T<:UnivariateDistribution{S},
            Tdists<:BroadcastVector{T},
        } = VectorOfUnivariate{S,T,Tdists}

        function Distributions._logpdf(
            dist::LazyVectorOfUnivariate,
            x::AbstractVector{<:Real},
        )
            return sum(copy(logpdf.(dist.v, x)))
        end

        function Distributions.logpdf(
            dist::LazyVectorOfUnivariate,
            x::AbstractMatrix{<:Real},
        )
            size(x, 1) == length(dist) ||
                throw(DimensionMismatch("Inconsistent array dimensions."))
            return vec(sum(copy(logpdf.(dists, x)), dims = 1))
        end

        const LazyMatrixOfUnivariate{
            S<:ValueSupport,
            T<:UnivariateDistribution{S},
            Tdists<:BroadcastArray{T,2},
        } = MatrixOfUnivariate{S,T,Tdists}

        function Distributions._logpdf(
            dist::LazyMatrixOfUnivariate,
            x::AbstractMatrix{<:Real},
        )
            return sum(copy(logpdf.(dist.dists, x)))
        end

        lazyarray(f, x...) = LazyArray(Base.broadcasted(f, x...))
        export lazyarray
    end
end

end


text:
# Include files

include("quadratic.jl");

f = (x) -> 2x^2 + 4x - 4;

println(quadratic(f))

max_stars_repo_name: neilpanchal/Julia-Algorithms

max_issues_repo_name: neilpanchal/Julia-Algorithms

max_forks_repo_name: neilpanchal/Julia-Algorithms


text:
function ExceptionalDataContext()
    return ExceptionalDataContext(())
end

function get_celestial_bodies(obj::ExceptionalDataContext)
    return jcall(obj, "getCelestialBodies", LazyLoadedCelestialBodies, ())
end

function get_data_providers_manager(obj::LazyLoadedDataContext)
    return jcall(obj, "getDataProvidersManager", DataProvidersManager, ())
end

function get_frames(obj::ExceptionalDataContext)
    return jcall(obj, "getFrames", LazyLoadedFrames, ())
end

function get_geo_magnetic_fields(obj::ExceptionalDataContext)
    return jcall(obj, "getGeoMagneticFields", LazyLoadedGeoMagneticFields, ())
end

function get_gravity_fields(obj::ExceptionalDataContext)
    return jcall(obj, "getGravityFields", LazyLoadedGravityFields, ())
end

function get_time_scales(obj::ExceptionalDataContext)
    return jcall(obj, "getTimeScales", LazyLoadedTimeScales, ())
end


text:
to_int(s) = parse(Int, s)

using DataStructures

function pass_day(fish)
    zero_count = length(filter(==(0), fish))
    fish = map(grow, fish)
    for _ = 1:zero_count
        append!(fish, 8)
    end
    fish
end

function solve(filename, num_days)
    fish = parse_file(filename)
    for day = 1:num_days
        fish = pass_day(fish)
    end
    length(fish)
end

solve("6eg.txt", 80) == 5934
solve("6.txt", 80)

grow(n) = n == 0 ? 6 : n - 1
parse_file(filename) = parse.(Int, split(read(filename, String), ","))

function solve2(filename, num_days)
    fishes = parse_file(filename)
    counts = counter(fishes)
    for _ in 1:num_days
        new_counts = counter(Int)
        for (fish, count) in counts
            new_counts[grow(fish)] = count
        end
        new_counts[8] += new_counts[0]
        counts = new_counts
    end
    sumcounts
end
solve2("6eg.txt", 256)
solve2("6eg.txt", 256) == 26984457539
solve2("6.txt", 80) == 5934


text:
using SeqLoggers
using Logging
using Test

# const SERVER_URL = "http://localhost:5341/"
const SERVER_URL = "http://subdn215:5341/"

logger = SeqLogger(SERVER_URL; App="Trialrun")
Logging.global_logger(logger)
@info "Test"
@info "Test invalid string\n, \r, \\ \""
@test length(logger.event_batch) == 2
SeqLoggers.flush_global_logger()
@test isempty(logger.event_batch)

unbatched_logger = SeqLogger(SERVER_URL; batch_size=1, App="Trialrun", Env="Test")
Logging.global_logger(unbatched_logger)
@info "Batch=1 test; doesn't need to flush"
@test isempty(unbatched_logger.event_batch)

# write TeeLogger
tee_logger = SeqLoggers.LoggingExtras.TeeLogger(
        SeqLogger(SERVER_URL; App="Trialrun"),
        Logging.ConsoleLogger(stderr)
)

Logging.with_logger(tee_logger) do
    @info "Test teelogger, plus does not need to flush, because run with `with_logger`"
end
@test isempty(tee_logger.loggers[1].event_batch)

SeqLoggers.Logging.global_logger(tee_logger)
@test SeqLoggers.Logging.global_logger() == tee_logger
@info "In Atom main module, the current logger is an atomlogger"
SeqLoggers.flush_current_logger() # flushes the current_logger()
@info "And to flush global logger, we need the flush_global_logger function"
SeqLoggers.flush_global_logger() # flushes the global_logger()

logger1 = Logging.ConsoleLogger(stderr)
logger2 = SeqLoggers.LoggingExtras.TeeLogger(
    SeqLogger(SERVER_URL; App="Trialrun", recursion1="yeeees"),
    Logging.ConsoleLogger(stderr)
)
logger3 = SeqLogger(SERVER_URL, ; App="Trialrun", recursion2="yes")
Logging.global_logger(logger1)
@info "Test logger recursion 0"
Logging.with_logger(logger2) do
    @info "Test logger recursion 1 {recursion1}"
    @test Logging.current_logger() == logger2
    @test Logging.global_logger() == logger1
    Logging.with_logger(logger3) do
        @info "Test logger recursion 2 {recursion2}"
        @test Logging.current_logger() == logger3
        @test Logging.global_logger() == logger1
    end
end

@testset "run_with_logger SeqLoggers Flush event batch" begin
logger = SeqLogger(SERVER_URL; App="Trialrun")

@test_throws DomainError run_with_logger(logger, -1) do x
    @info "Run stuff"
    sqrt(x)
end

@test isempty(logger.event_batch)
end

@testset "Throw exception in all available Loggers" begin

console_logger = global_logger()
@test_throws ArgumentError run_with_logger(console_logger) do
    throw(ArgumentError("Tester, indeed!"))
end

seq_logger = SeqLogger(SERVER_URL; App="TrialRun")
@test_throws ArgumentError run_with_logger(seq_logger) do
    @info "sfd"
    @warn "sfd"
    @error "Logged error with backtrace" back_trace="This is a backtrace"
    @error "Logged error without backtrace"
    throw(ArgumentError("Real exception thrown"))
end

# Run this test after load logger from config was successful
file_path = joinpath(@__DIR__, "data", "config.json")
tee_logger = SeqLoggers.load_logger_from_config(file_path)
@test_throws ArgumentError run_with_logger(tee_logger) do
    @info "sfd"
    @warn "sfd"
    @error "Logged error with backtrace" back_trace="This is a backtrace"
    @error "Logged error without backtrace"
    throw(ArgumentError("Real exception thrown"))
end

end


text:
module MonteCarlo

using Random
using Statistics

# European Option Calculation from
# https://groups.google.com/forum/?hl=en&fromgroups=#!topic/julia-dev/ImhGsqX_IHc

function perf_euro_option_devec(npaths)
    steps = 250
    r = 0.05
    sigma = .4
    T = 1
    dt = T/steps
    K = 100
    S = fill(100.0, npaths)
    t1 = (r - 0.5*sigma^2)*dt
    t2 = sigma*sqrt(dt)

    for i=1:steps
        for j=1:npaths
            S[j] *= exp(t1 + t2 * randn())
        end
    end

    return mean(exp(-r*T) .* max.(K .- S, 0))
end

function perf_euro_option_vec(npaths)
    steps = 250
    r = 0.05
    sigma = .4
    T = 1
    dt = T/steps
    K = 100
    S = fill(100.0, npaths)
    t1 = (r - 0.5*sigma^2)*dt
    t2 = sigma*sqrt(dt)
    R = Array{Float64}(undef, npaths)

    for i=1:steps
        S .*= exp.(t2 .* randn!(R) .+ t1)
    end

    return mean(exp(-r*T) .* max.(K .- S, 0))
end

end # module


text:
push!(LOAD_PATH, pwd())

using OffsetArrays
using DataStructures
using Circular

function Base.insert!(a::OffsetArray{T, 1, Array{T, 1}}, i::Integer, v::T) where T
    j = 1 + i - firstindex(a)
    insert!(a.parent, j, v)
end

function Base.splice!(a::OffsetArray{T, 1, Array{T, 1}}, i::Integer) where T
    j = 1 + i - firstindex(a)
    splice!(a.parent, j)
end

function solve011(filename)
    println(filename)
    basefreq = 0
    visited = Set([basefreq])
    open(filename) do file
        for line in eachline(file)
            change = parse(Int, line)
            basefreq += change
        end
    end

    return basefreq
end

function solve012(filename)
    nstep = 0
    basefreq = 0
    visited = Set([basefreq])
    changes = readlines(filename)
    while true
        for change_str in changes
            change = parse(Int, change_str)
            basefreq += change
            if in(basefreq, visited)
                return basefreq
            end
            union!(visited, Set([basefreq]))
        end
    end
end

function solve021(filename)
    count2, count3 = open(filename) do input
        count2 = 0
        count3 = 0
        counts = zeros(Int, 26)
        offset = Int('a')
        for id in eachline(input)
            # println(id)
            for char in id
                loc = 1 + Int(char) - offset
                counts[loc] += 1
            end
            # println(counts)
            n2 = any(counts .== 2)
            n3 = any(counts .== 3)
            # println("$n2 $n3")
            count2 += Int(n2)
            count3 += Int(n3)
            # println("$count2 $count3")
            # Reset and check next id.
            counts .= 0
        end
        return count2, count3
    end
    return count2 * count3
end

function solve022(filename)
    ids = readlines(filename)
    dist = id_dist.(reshape(ids, :, 1), reshape(ids, 1, :))
    # println(dist)

    id_indx = findfirst(dist) do d
        d == 1
    end
    correct_id1 = ids[id_indx[1]]
    correct_id2 = ids[id_indx[2]]
    println(correct_id1)
    println(correct_id2)
    matching = Array{Char}(undef, length(correct_id1) - 1)
    i = 1
    for cc in zip(correct_id1, correct_id2)
        c1, c2 = cc
        if c1 == c2
            matching[i] = c1
            i += 1
        end
    end
    return string(matching...)
end

function id_dist(id1, id2)
    diff = 0
    # println(id1)
    # println(id2)
    for cc in zip(split(id1, ""), split(id2, ""))
        c1, c2 = cc
        diff += Int(c1 != c2)
    end
    # println("dist $diff")
    return diff
end

function solve031(filename)
    claim_count = count_claims(readlines(filename))
    # display(claim_count[1:9, 1:9]); println()
    return sum(claim_count .>= 2)
end

function count_claims(claims)
    claim_count = zeros(Int, 1000, 1000)
    for claim in claims
        id, left, top, width, height = parse_claim(claim)
        claim_count[1+top:top+height, 1+left:left+width] .+= 1
    end
    return claim_count
end

function parse_claim(claim_str)
    claim_regex = r"#([0-9]*) @ ([0-9]*),([0-9]*): ([0-9]*)x([0-9]*)"
    id, left, top, width, height = match(claim_regex, claim_str).captures
    left, top, width, height = parse.(Int, [left, top, width, height])
    return id, left, top, width, height
end

function solve032(filename)
    claims = readlines(filename)
    claim_count = count_claims(claims)
    # display(claim_count[1:9, 1:9]); println()
    for claim in claims
        id, left, top, width, height = parse_claim(claim)
        claimed_section = claim_count[1+top:top+height, 1+left:left+width]
        # display(claimed_section); println()
        if all(claimed_section .== 1)
            return id
        end
    end
end

function solve041(filename)
    records = readlines(filename)
    minutes_asleep = create_sleep_table(records)
    most_sleep = first(sort(
        collect(minutes_asleep), by=entry->sum(entry[2]), rev=true))
    most_sleep_id = first(most_sleep)
    id_number = parse(Int, most_sleep_id[2:lastindex(most_sleep_id)])
    asleep = last(most_sleep)
    most_sleep_minute = argmax(asleep) - 1
    return id_number * most_sleep_minute
end

function create_sleep_table(records)
    sort!(records)
    n = length(records)
    minutes_asleep = Dict{String, Array{Int}}()
    record_regex = r"\[.*:([0-9]*)\].*(#[0-9]*|falls|wakes).*"
    asleep = Array{Int}(undef, 60)
    sleep_start::Int = 0
    sleep_end::Int = 0
    for record in records
        println(record)
        m = match(record_regex, record)
        if m[2][1] == '#'
            id = m[2]
            if haskey(minutes_asleep, id)
                asleep = minutes_asleep[id]
            else
                asleep = zeros(Int, 60)
                minutes_asleep[id] = asleep
            end
        elseif m[2] == "falls"
            sleep_start = parse(Int, m[1])
        else
            sleep_end = parse(Int, m[1])
            # Array is 1 to 60, minutes 0 to 59, so offset by 1.
            asleep[1+sleep_start:sleep_end] .+= 1
        end
    end
    return minutes_asleep
end

function print_sleep_table(sleep_table)
    print("    ")
    for i in 1:60
        print(div(i-1, 10))
    end
    print("\n    ")
    for i in 1:60
        print(mod(i-1, 10))
    end
    println()

    for entry in sleep_table
        id, asleep = entry
        print("$id ")
        for a in asleep
            print(a)
        end
        println()
    end
end

function solve042(filename)
    records = readlines(filename)
    minutes_asleep = create_sleep_table(records)
    print_sleep_table(minutes_asleep)
    most_sleep = first(sort(
        collect(minutes_asleep), by=entry->maximum(entry[2]), rev=true))
    most_sleep_id = first(most_sleep)
    id_number = parse(Int, most_sleep_id[2:lastindex(most_sleep_id)])
    asleep = last(most_sleep)
    most_sleep_minute = argmax(asleep) - 1
    return id_number, most_sleep_minute, id_number * most_sleep_minute
end

function solve051(filename)
    polymer = readlines(filename)[1]
    remaining = react_polymer(polymer)
    return length(remaining)
end

function react_polymer(polymer)
    # println(polymer)
    remaining = Stack{Char}()
    current = 'a'
    for current in polymer
        # println(current)
        if length(remaining) > 0
            lastc = pop!(remaining)
            if lowercase(lastc) == lowercase(current) && lastc != current
                # lastc and current cancel, so do nothing and they are gone.
            else
                push!(remaining, lastc)
                push!(remaining, current)
            end
        else
            push!(remaining, current)
        end
        # println(remaining)
    end
    return remaining
end

function solve052(filename)
    polymer = readlines(filename)[1]
    polymer_lengths = Array{Int}(undef, 26)
    for entry in enumerate("abcdefghijklmnopqrstuvwxyz")
        i, c = entry
        stripped = strip_polymer(polymer, c)
        # println(reverse(collect(stripped)))
        remaining = react_polymer(reverse(collect(stripped)))
        n = length(remaining)
        println(entry..., ": ", n)
        polymer_lengths[i] = n
    end
    return minimum(polymer_lengths)
end

function strip_polymer(polymer, c)
    remaining = Stack{Char}()
    for current in polymer
        if lowercase(c) != lowercase(current)
            push!(remaining, current)
        end
    end
    # println(c)
    # println(remaining)
    return remaining
end

function solve061(filename)
    point_strs = readlines(filename)
    points = transpose(parse.(Int, hcat(split.(point_strs, ",")...)))
    # display(points)
    dist = Array{Int}(undef, size(points, 1))
    grid = zeros(Int, maximum(points[:,1])+2, maximum(points[:,2])+2)
    grid_x_axis, grid_y_axis = axes(grid)
    for ix in grid_x_axis, iy in grid_y_axis
        x, y = ix-1, iy-1
        for i in axes(points,1)
            dist[i] = manhattan_dist([x,y], points[i,:])
        end
        # Leave ties unclaimed.
        claimed = argmin(dist)
        if sum(dist .== dist[claimed]) > 1
            continue
        end
        grid[ix, iy] = claimed
    end
    # display(grid)

    # Look at boundary to find unbounded point regions.
    bounded = fill(true, size(points, 1))
    for ix in grid_x_axis
        iymin, iymax = firstindex(grid, 2), lastindex(grid, 2)
        if grid[ix, iymin] > 0
            bounded[grid[ix, iymin]] = false
        end
        if grid[ix, iymax] > 0
            bounded[grid[ix, iymax]] = false
        end
    end
    for iy in grid_y_axis
        ixmin, ixmax = firstindex(grid, 1), lastindex(grid, 1)
        if grid[ixmin, iy] > 0
            bounded[grid[ixmin, iy]] = false
        end
        if grid[ixmax, iy] > 0
            bounded[grid[ixmax, iy]] = false
        end
    end
    # display(bounded)

    areas = Array{Int}(undef, size(points, 1))
    for i in axes(points, 1)
        areas[i] = sum(grid .== i)
    end
    largest_bounded_area = maximum(areas[bounded])
    return largest_bounded_area
end

function manhattan_dist(x1, x2)
    abs(x1[1]-x2[1]) + abs(x1[2]-x2[2])
end

function solve062(filename)
    MAX_DIST = 10000
    point_strs = readlines(filename)
    points = transpose(parse.(Int, hcat(split.(point_strs, ",")...)))
    # display(points)
    dist = Array{Int}(undef, size(points, 1))
    grid = zeros(Int, maximum(points[:,1])+2, maximum(points[:,2])+2)
    grid_x_axis, grid_y_axis = axes(grid)
    for iy in grid_y_axis, ix in grid_x_axis
        x, y = ix-1, iy-1
        for i in axes(points,1)
            dist[i] = manhattan_dist([x,y], points[i,:])
        end
        if sum(dist) < MAX_DIST
            grid[ix, iy] = 1
        end
    end
    # display(grid)
    return sum(grid)
end

function solve071(filename)
    nodes, pred_dag, succ_dag = open(build_dag, filename)
    next_nodes = binary_minheap(Char)
    remaining_nodes = Set(nodes)
    completed_nodes = fill('0', length(remaining_nodes))
    schedule_nodes!(next_nodes, remaining_nodes, pred_dag)
    # display(pred_dag)
    # display(next_nodes)
    # display(remaining_nodes)
    # display(completed_nodes)
    for i in axes(completed_nodes, 1)
        n = pop!(next_nodes)
        completed_nodes[i] = n
        clear_pred!(pred_dag, n)
        schedule_nodes!(next_nodes, remaining_nodes, pred_dag)
        # display(pred_dag)
        # display(next_nodes)
        # display(remaining_nodes)
        # display(completed_nodes)
    end
    return string(completed_nodes...)
end

function schedule_nodes!(next_nodes, remaining_nodes, pred_dag)
    # Find workable nodes and queue them up.
    for n in copy(remaining_nodes)
        if length(pred_dag[n]) == 0
            push!(next_nodes, n)
            pop!(remaining_nodes, n)
        end
    end
end

function clear_pred!(pred_dag, n)
    for arcset in pred_dag
        orig, pred = arcset
        if n in pred
            pop!(pred, n)
        end
    end
end

function build_dag(input)
    nodes = Set{Char}()
    pred_dag = Dict{Char, Set{Char}}()
    succ_dag = Dict{Char, Set{Char}}()
    for line = eachline(input)
        step_re = r"Step ([A-Z]).*step ([A-Z]).*"
        pred, succ = match(step_re, line).captures
        pred, succ = pred[1], succ[1]  # String to Char.
        if haskey(pred_dag, succ)
            pred_dag[succ] = push!(pred_dag[succ], pred)
        else
            pred_dag[succ] = Set{Char}([pred])
        end
        if haskey(succ_dag, pred)
            succ_dag[pred] = push!(succ_dag[pred], succ)
        else
            succ_dag[pred] = Set{Char}([succ])
        end
        push!(nodes, pred)
        push!(nodes, succ)
    end
    for n in nodes
        if !haskey(pred_dag, n); pred_dag[n] = Set{Char}(); end
        if !haskey(succ_dag, n); succ_dag[n] = Set{Char}(); end
    end
    return nodes, pred_dag, succ_dag
end

function solve072(filename)
    N_WORKERS = 5
    TIME_OFFSET = 60
    CHAR_OFFSET = Int('A') - 1
    time_required = [Int(c) - CHAR_OFFSET + TIME_OFFSET for c in 'A':'Z']
    nodes, pred_dag, succ_dag = open(build_dag, filename)
    next_nodes = binary_minheap(Char)
    available_workers = Set(1:N_WORKERS)
    next_workers = PriorityQueue{Tuple{Int, Char, Int}, Int}()
    remaining_nodes = Set(nodes)
    completed_nodes = fill('0', length(remaining_nodes))
    schedule_nodes!(next_nodes, remaining_nodes, pred_dag)
    t = 0
    schedule_workers!(next_workers, next_nodes, available_workers, t, time_required)
    # display(pred_dag)
    # display(next_workers)
    # display(remaining_nodes)
    # display(completed_nodes)
    for i in axes(completed_nodes, 1)
        w, n, tf = dequeue!(next_workers)
        completed_nodes[i] = n
        push!(available_workers, w)
        t = tf
        clear_pred!(pred_dag, n)
        schedule_nodes!(next_nodes, remaining_nodes, pred_dag)
        schedule_workers!(next_workers, next_nodes, available_workers, t, time_required)
        # display(pred_dag)
        # display(next_workers)
        # display(remaining_nodes)
        # display(completed_nodes)
    end
    return string(completed_nodes...), t
end

function schedule_workers!(next_workers, next_nodes, available_workers, t, time_required)
    CHAR_OFFSET = Int('A') - 1
    while length(next_nodes) > 0 && length(available_workers) > 0
        n = pop!(next_nodes)
        w = pop!(available_workers)
        i = Int(n) - CHAR_OFFSET
        dt = time_required[i]
        enqueue!(next_workers, (w, n, t+dt)=>t+dt)
    end
end

function solve081(filename)
    CHILD = 1
    METADATA = 2
    HEADER1 = 3
    HEADER2 = 4
    inputs = readline(filename)
    tot::Int = 0
    next_mode = Stack{Int}()
    push!(next_mode, HEADER1)
    local mode::Int = 0
    n_child::Int = 0
    n_metadata::Int = 0
    for c in split(inputs)
        mode = pop!(next_mode)
        val = parse(Int, c)
        # println("Pre ", mode, " ", val, " ", next_mode)
        if mode == HEADER1
            # Reading number of children
            n_child = val
            push!(next_mode, HEADER2)
        elseif mode == HEADER2
            # Reading number of metadata
            n_metadata = val
            if n_metadata > 0
                push!(next_mode, repeat([METADATA], n_metadata)...)
            end
            if n_child > 0
                push!(next_mode, repeat([HEADER1], n_child)...)
            end
        elseif mode == METADATA
            tot += val
        end
        # println("Post ", tot, " ", next_mode)
    end
    return tot
end

mutable struct Node{T}
    val::T
    children::Array{T,1}
    metadata::Array{T,1}
end

function solve082(filename)
    CHILD = 1
    METADATA = 2
    HEADER1 = 3
    HEADER2 = 4
    inputs = readline(filename)
    tot::Int = 0
    tree = Node{Int}[]
    push!(tree, Node{Int}(-1, [], []))
    max_node_idx::Int = 1
    next_mode = Stack{Tuple{Int, Int, Int}}()
    push!(next_mode, (HEADER1, 0, max_node_idx))
    local mode::Int = 0
    n_child::Int = 0
    n_metadata::Int = 0
    for c in split(inputs)
        mode, parent, node_idx = pop!(next_mode)
        val = parse(Int, c)
        # println("Pre ", mode, " ", parent, " ", next_mode)
        # display(tree); println()
        if mode == HEADER1
            # Reading number of children
            n_child = val
            push!(next_mode, (HEADER2, parent, node_idx))
        elseif mode == HEADER2
            # Reading number of metadata
            n_metadata = val
            if n_metadata > 0
                push!(next_mode, repeat([(METADATA, node_idx, 0)], n_metadata)...)
            end
            if n_child > 0
                push!(next_mode, [(HEADER1, node_idx, max_node_idx + i) for i in n_child:-1:1]...)
                for i in 1:n_child
                    push!(tree, Node{Int}(-1, [], []))
                end
                max_node_idx += n_child
            end
            if parent > 0
                parent_node = tree[parent]
                push!(parent_node.children, node_idx)
            end
        elseif mode == METADATA
            parent_node = tree[parent]
            push!(parent_node.metadata, val)
        end
        # println("Post ", parent, " ", next_mode)
        # display(tree); println()
    end

    setval!(tree, tree[1])
    display(tree); println()
    return tree[1].val
end

function setval!(tree::Array{Node{T},1}, n::Node{T}) where T<:Integer
    if n.val >= 0
        return
    elseif length(n.children) == 0
        n.val = sum(n.metadata)
    else
        vals = similar(n.children)
        for tup in enumerate(n.children)
            idx, i = tup
            c = tree[i]
            if c.val < 0
                setval!(tree, c)
            end
            vals[idx] = c.val
        end
        n.val = 0
        for m in n.metadata
            if m <= length(n.children)
                n.val += vals[m]
            end
        end
    end
end

function solve091(args...)
    if length(args) == 1
        re = r"([0-9]*) players; last marble is worth ([0-9]*) points"
        input_str = string(open(readchomp, args[1]))
        println(input_str)
        input_m = match(re, input_str)
        show(input_m); println()
        if input_m !== nothing
            input = input_m.captures
        else
            error("Invalid input file.")
        end
    elseif length(args) == 2
        input = args
    else
        error("Must provide at least two extra arguments: N_PLAYERS and MAX_MARBLE")
    end

    n_players, max_marble = map(x -> parse(Int, x), input)
    # Increasing index is clockwise.
    marbles = OffsetArray{Int}(undef, 0:0)
    marbles[0] = 0
    current_marble_loc = 0
    current_player = 1
    scores = zeros(Int, n_players)
    for m = 1:max_marble
        if mod(m, 23) ≠ 0
            # Insert marble between marbles 1 and 2 clockwise from current.
            i = mod(current_marble_loc + 2, length(marbles))
            insert!(marbles, i, m)
            current_marble_loc = i
        else
            # Score marble and remove.
            scores[current_player] += m
            i_remove = mod(current_marble_loc - 7, length(marbles))
            scores[current_player] += marbles[i_remove]
            splice!(marbles, i_remove)
            current_marble_loc = i_remove
        end
        # println(current_player, ", ", m, ": ", circshift(marbles, 1))
        current_player = mod(current_player + 1, n_players+1)
        (current_player == 0) && (current_player = 1)
    end

    return maximum(scores)
end

function solve092(args...)
    println("solve092 ", args)
    if length(args) == 1
        re = r"([0-9]*) players; last marble is worth ([0-9]*) points"
        input_str = string(open(readchomp, args[1]))
        println(input_str)
        input_m = match(re, input_str)
        show(input_m); println()
        if input_m !== nothing
            input = input_m.captures
        else
            error("Invalid input file.")
        end
    elseif length(args) == 2
        input = args
    else
        error("Must provide at least two extra arguments: N_PLAYERS and MAX_MARBLE")
    end

    n_players, max_marble = map(x -> parse(Int, x), input)
    # Increasing index is clockwise.
    marbles = circular_list(0)
    current_player = 1
    scores = zeros(Int, n_players)
    println("Progress: (. = 10000 marbles)")
    for m = 1:max_marble
        mod(m, 10_000) == 0 && print(".")
        if mod(m, 23) ≠ 0
            # Insert marble between marbles 1 and 2 clockwise from current.
            marbles = insertafter!(next(marbles), m)
        else
            # Score marble and remove.
            scores[current_player] += m
            for i = 1:7; marbles = prev(marbles); end
            scores[current_player] += value(marbles)
            removeleftshift!(marbles)
        end
        # println(current_player, ", ", m, ": ", circshift(marbles, 1))
        current_player = mod(current_player + 1, n_players+1)
        (current_player == 0) && (current_player = 1)
    end
    println()
    return maximum(scores)
end


text:
using LinearAlgebra

# Fourier Types {Parameter type, full or half flag, Differentiable Flag}
mutable struct ConcatBasis{T,TB1,TB2} <: AbstractFuncApprox where {T<:Real,TB1,TB2}
    ϕ1::TB1
    ϕ2::TB2

    function ConcatBasis(::Type{T}, ϕ1, ϕ2) where {T}
        new{T,typeof(ϕ1),typeof(ϕ2)}(clone(ϕ1), clone(ϕ2))
    end
end

function call(ϕ::ConcatBasis{T}, x) where {T<:Real}
    N = get_num_outputs(ϕ)
    N1 = get_num_outputs(ϕ.ϕ1)
    feats = zeros(T, N)
    call!(view(feats, 1:N1), ϕ.ϕ1, x)
    # call!(view(feats, N1:N), ϕ.ϕ2, x)

    return feats
end

function call!(out, ϕ::ConcatBasis{T}, x) where {T<:Real}
    N = get_num_outputs(ϕ)
    N1 = get_num_outputs(ϕ.ϕ1)
    call!(view(out, 1:N1), ϕ.ϕ1, x)
    # call!(view(out, N1:N), ϕ.ϕ2, x)
end

function get_num_outputs(ϕ::ConcatBasis{T})::Int where {T}
    return get_num_outputs(ϕ.ϕ1) + get_num_outputs(ϕ.ϕ2)
end

function get_output_range(ϕ::ConcatBasis{T}) where {T}
    N = get_num_outputs(ϕ)
    N1 = get_num_outputs(ϕ.ϕ1)
    ranges = zeros(T, (N, 2))
    ranges[1:N1, :] .= get_output_range(ϕ.ϕ1)
    ranges[N1:N, :] .= get_output_range(ϕ.ϕ2)

    return ranges
end

function clone(ϕ::ConcatBasis{T})::ConcatBasis{T} where {T<:Real}
    ϕ2 = ConcatBasis(T, ϕ.ϕ1, ϕ.ϕ2)
    return ϕ2
end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
using SolidStateDetectors
using IntervalSets
using SpecialFunctions
T=Float32

struct CustomChargeDensity{T, N} <: SolidStateDetectors.AbstractChargeDensity{T}
	G1::SolidStateDetectors.AbstractGeometry{T,N}
	G2::SolidStateDetectors.AbstractGeometry{T,N}
	L1::SolidStateDetectors.AbstractChargeDensity{T}
	L2::SolidStateDetectors.AbstractChargeDensity{T}
	G3::SolidStateDetectors.AbstractGeometry{T,N}
	L3::SolidStateDetectors.AbstractChargeDensity{T}
	G4::SolidStateDetectors.AbstractGeometry{T,N}
	L4::SolidStateDetectors.AbstractChargeDensity{T}
	G5::SolidStateDetectors.AbstractGeometry{T,N}
	L5::SolidStateDetectors.AbstractChargeDensity{T}
	G6::SolidStateDetectors.AbstractGeometry{T,N}
	L6::SolidStateDetectors.AbstractChargeDensity{T}
	G7::SolidStateDetectors.AbstractGeometry{T,N}
	L7::SolidStateDetectors.AbstractChargeDensity{T}
	G8::SolidStateDetectors.AbstractGeometry{T,N}
	L8::SolidStateDetectors.AbstractChargeDensity{T}
	G9::SolidStateDetectors.AbstractGeometry{T,N}
	L9::SolidStateDetectors.AbstractChargeDensity{T}
	G10::SolidStateDetectors.AbstractGeometry{T,N}
	L10::SolidStateDetectors.AbstractChargeDensity{T}
	G11::SolidStateDetectors.AbstractGeometry{T,N}
	L11::SolidStateDetectors.AbstractChargeDensity{T}
	G12::SolidStateDetectors.AbstractGeometry{T,N}
	L12::SolidStateDetectors.AbstractChargeDensity{T}
	G13::SolidStateDetectors.AbstractGeometry{T,N}
	L13::SolidStateDetectors.AbstractChargeDensity{T}
	G14::SolidStateDetectors.AbstractGeometry{T,N}
	L14::SolidStateDetectors.AbstractChargeDensity{T}
	G15::SolidStateDetectors.AbstractGeometry{T,N}
	L15::SolidStateDetectors.AbstractChargeDensity{T}
	G16::SolidStateDetectors.AbstractGeometry{T,N}
	L16::SolidStateDetectors.AbstractChargeDensity{T}
	G17::SolidStateDetectors.AbstractGeometry{T,N}
	L17::SolidStateDetectors.AbstractChargeDensity{T}
end

function SolidStateDetectors.get_charge_density(cdm::CustomChargeDensity{T}, pt::SolidStateDetectors.AbstractCoordinatePoint{T})::T where {T}
	(pt in cdm.G1) * SolidStateDetectors.get_charge_density(cdm.L1, pt) + 800000000000.0 * (1 + erf((sqrt(pt[1]^2+pt[2]^2)-(0.7))/(2*0.1)))/2 +
	(pt in cdm.G2) * SolidStateDetectors.get_charge_density(cdm.L2, pt) +
	(pt in cdm.G3) * SolidStateDetectors.get_charge_density(cdm.L3, pt) +
	(pt in cdm.G4) * SolidStateDetectors.get_charge_density(cdm.L4, pt) +
	(pt in cdm.G5) * SolidStateDetectors.get_charge_density(cdm.L5, pt) +
	(pt in cdm.G6) * SolidStateDetectors.get_charge_density(cdm.L6, pt) +
	(pt in cdm.G7) * SolidStateDetectors.get_charge_density(cdm.L7, pt) +
	(pt in cdm.G8) * SolidStateDetectors.get_charge_density(cdm.L8, pt) +
	(pt in cdm.G9) * SolidStateDetectors.get_charge_density(cdm.L9, pt) +
	(pt in cdm.G10) * SolidStateDetectors.get_charge_density(cdm.L10, pt) * erfc(-(pt[3]-(1.0074e-05))/(2*3.7e-08))/2 +
	(pt in cdm.G11) * SolidStateDetectors.get_charge_density(cdm.L11, pt) * (1 + erf(-(pt[3]-(0.0019797316718427))/(2*1.3416407864998738e-07)))/2 +
	(pt in cdm.G12) * SolidStateDetectors.get_charge_density(cdm.L12, pt) * (1 + erf(-(pt[3]-(0.0019797316718427))/(2*1.3416407864998738e-07)))/2 +
	(pt in cdm.G13) * SolidStateDetectors.get_charge_density(cdm.L13, pt) * (1 + erf(-(pt[3]-(0.0019797316718427))/(2*1.3416407864998738e-07)))/2 +
	(pt in cdm.G14) * SolidStateDetectors.get_charge_density(cdm.L14, pt) * (1 + erf(-(pt[3]-(0.0019797316718427))/(2*1.3416407864998738e-07)))/2 +
	(pt in cdm.G15) * SolidStateDetectors.get_charge_density(cdm.L15, pt) * (1 + erf(-(pt[3]-(0.0019797316718427))/(2*1.3416407864998738e-07)))/2 +
	(pt in cdm.G16) * SolidStateDetectors.get_charge_density(cdm.L16, pt) * (1 + erf(-(pt[3]-(0.0019797316718427))/(2*1.3416407864998738e-07)))/2 +
	(pt in cdm.G17) * SolidStateDetectors.get_charge_density(cdm.L17, pt) * (1 + erf(-(pt[3]-(0.0019797316718427))/(2*1.3416407864998738e-07)))/2
end

G1 = SolidStateDetectors.Tube{T}(0..0.01627006165897972,0..2π,0..0.002,CartesianPoint{T}(-1.2390881971262907e-19, 9.912705577010326e-19, 0))
G2 = SolidStateDetectors.Tube{T}(0..0.01627006165897972,0..2π,0..1e-05,CartesianPoint{T}(-1.2390881971262907e-19, 9.912705577010326e-19, 0))
CD1 = SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,50000000000.0),(0.0f0,0.0f0,0.0f0))
CD2 = SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,-1.9e+16),(0.0f0,0.0f0,0.0f0))
G3 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 2e-05,CartesianPoint{T}(0, 0, 0.00199),0)
CD3= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))
G4 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 2e-05,CartesianPoint{T}(0.009020061658979717, 0.0, 0.00199),0)
CD4= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))
G5 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 2e-05,CartesianPoint{T}(0.0045100308294898596, 0.0078116025403784425, 0.00199),0)
CD5= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))
G6 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 2e-05,CartesianPoint{T}(-0.004510030829489857, 0.007811602540378445, 0.00199),0)
CD6= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))
G7 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 2e-05,CartesianPoint{T}(-0.009020061658979717, 1.1046389638781276e-18, 0.00199),0)
CD7= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))
G8 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 2e-05,CartesianPoint{T}(-0.004510030829489863, -0.007811602540378441, 0.00199),0)
CD8= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))
G9 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 2e-05,CartesianPoint{T}(0.004510030829489859, -0.0078116025403784425, 0.00199),0)
CD9= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))
G10 = SolidStateDetectors.Tube{T}(0..0.01627006165897972,0..2π,0..1.48e-07,CartesianPoint{T}(-1.2390881971262907e-19, 9.912705577010326e-19, 1.0001000000000001e-05))
CD10= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,-1.9e+16),(0.0f0,0.0f0,0.0f0))
G11 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 5.366563145999495e-07,CartesianPoint{T}(0, 0, 0.0019797316718427),0)
CD11= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))
G12 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 5.366563145999495e-07,CartesianPoint{T}(0.009020061658979717, 0.0, 0.0019797316718427),0)
CD12= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))
G13 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 5.366563145999495e-07,CartesianPoint{T}(0.0045100308294898596, 0.0078116025403784425, 0.0019797316718427),0)
CD13= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))
G14 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 5.366563145999495e-07,CartesianPoint{T}(-0.004510030829489857, 0.007811602540378445, 0.0019797316718427),0)
CD14= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))
G15 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 5.366563145999495e-07,CartesianPoint{T}(-0.009020061658979717, 1.1046389638781276e-18, 0.0019797316718427),0)
CD15= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))
G16 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 5.366563145999495e-07,CartesianPoint{T}(-0.004510030829489863, -0.007811602540378441, 0.0019797316718427),0)
CD16= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))
G17 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 5.366563145999495e-07,CartesianPoint{T}(0.004510030829489859, -0.0078116025403784425, 0.0019797316718427),0)
CD17= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))

ccdm = CustomChargeDensity(G1,G2,CD1,CD2,G3,CD3,G4,CD4,G5,CD5,G6,CD6,G7,CD7,G8,CD8,G9,CD9,G10,CD10,G11,CD11,G12,CD12,G13,CD13,G14,CD14,G15,CD15,G16,CD16,G17,CD17)

max_stars_repo_name: leenderthayen/NDP

max_issues_repo_name: leenderthayen/NDP

max_forks_repo_name: leenderthayen/NDP


text:
@enum FCType fcgiven fclevel fcslope fcnatural
const global fcrate = fcslope
export FCType, fcgiven, fclevel, fcslope, fcrate, fcnatural

"""
    FCType

An enumeration type representing the possible types of final conditions.

  * `fcgiven` - the final conditions are provided by the user.
  * `fclevel` - the final conditions are prescribed from the levels of the
    steady state solution.
  * `fcrate` or `fcslope` - the final conditions are such that the slope of the dynamic
    solution matches the slope of the steady state solution.
  * `fcnatural` - the final condition sets the second difference equal to zero.

"""
FCType

"""
    AbstractSolverData

An abstract type for the internal data structure of a dynamic solver.

### Implementation (for developers)

To implement a derived class, one must specialize the following methods. Note
that in the signatures of most of these methods, the `sd` argument of the new
type is the last argument..

  * Constructor - normally it would take a model and a plan together with a
    final condition type.
  * `global_RJ(x, exog, sd)` - compute the residual and the Jacobian (return a
    2-tuple) of the system at the given point `x` and exogenous data `exog`.
  * `global_R!(res, x, exog, sd)` (optional) - compute the residual in place.
  * `assign_exog_data!(x, exog, sd)` - this function applies the exogenous
    constraints as per the plan and the given exogenous data.
  * `assign_final_condition!(x, exog, sd, ::Val{T})` where T is one of the
    FCType constants. This method applies the given type of final condition. If
    T is fcgiven, the values of the final conditions are taken from exog.
  * `assign_update_step(x, lambda, dx, sd)` - this method must implement x = x +
    lambda * dx
  * `update_plan!(sd, model, plan)` - update the solver data to reflect the given
    simulation plan.

"""
abstract type AbstractSolverData end

"""
    global_RJ(x::Array, exog::Array, sd::AbstractSolverData)

Compute the residual and the Jacobian of the dynamic system of equations at the
given point `x` using the given exogenous data `exog`. Returns a tuple of two
elements. The first is a vector of residuals and the second is the Jacobian
matrix.

"""
global_RJ(point::AbstractArray{Float64}, exog_data::AbstractArray{Float64}, g::AbstractSolverData) = error("Not implemented for $(typeof(g)).")
# export global_RJ

"""
    global_R!(res::Vector, x::Array, exog::Array, sd::AbstractSolverData)

Compute the residual of the dynamic system of equations at the given point `x` using the given exogenous data `exog`.
The residual is updated in place.
"""
global_R!(res::AbstractArray{Float64}, point::AbstractArray{Float64}, exog_data::AbstractArray{Float64}, g::AbstractSolverData) = error("Not implemented for $(typeof(g)).")
# export global_R!

"""
    assign_exog_data(x::Array, exog::Array, sd::AbstractSolverData)

Apply the exogenous constraints according to the simulation plan (passed to the
constructor of sd) with exogenous values taken from `exog`.

"""
assign_exog_data!(x::AbstractArray{Float64,2}, exog::AbstractArray{Float64,2}, g::AbstractSolverData) = error("Not implemented for $(typeof(g)).")
# export assign_exog_data!

"""
    assign_final_condition!(x::Array, exog::Array, sd::AbstractSolverData, ::Val{T})

Apply final conditions of the given type. The type parameter `T` must be one of
the constants of type `FCType`. If `T` is fcgiven, the final condition values
are taken from `exog`. If `T` is one of `fclevel` or `fcslope`, then the final
conditions are taken from the steady state solution in the model instance passed
to the constructor of `sd`.

"""
assign_final_condition!(x::AbstractArray{Float64,2}, exog::AbstractArray{Float64,2}, g::AbstractSolverData, ::Val{T}) where T = error("Not implemented for $(typeof(g)) with final condition $(T).")
# export assign_final_condition!

"""
    assign_update_step!(x::Array, lambda::Real, dx::Array, sd::AbstractSolverData)

Assign the iterative update. Equivalent to `x = x + lambda * dx`.
"""
function assign_update_step! end
@inline assign_update_step!(x::AbstractArray{Float64}, lambda::Float64, Δx::AbstractArray{Float64}, ::AbstractSolverData) = BLAS.axpy!(lambda, Δx, x)  # axpy!(a,x,y) <=> y .= a .* x .+ y
# export assign_update_step!

"""
    update_plan!(sd, model, plan)

Update the solver data structure for the given simulation plan. The `model` must
be the same instance that was used to create `sd`. Typically, the plan should
span the same range as the one used when `sd` was created, although it might
have different exogenous/endogenous assignments.

"""
update_plan!(sd::AbstractSolverData, m::Model, p::Plan) = error("Not implemented for $(typeof(sd)).")
# export update_plan!


text:
using ClimateMachine
const clima_dir = dirname(dirname(pathof(ClimateMachine)));

if parse(Bool, get(ENV, "CLIMATEMACHINE_PLOT_EDMF_COMPARISON", "false"))
    plot_dir = joinpath(clima_dir, "output", "bomex_edmf", "pycles_comparison")
else
    plot_dir = nothing
end

include(joinpath(@__DIR__, "compute_mse.jl"))

#! format: off
best_mse = Dict()
best_mse[:Bomex] = Dict()
best_mse[:Bomex]["ρ"] = 3.4943021267397123e-02
best_mse[:Bomex]["ρu[1]"] = 3.0714039084256679e+03
best_mse[:Bomex]["ρu[2]"] = 1.3375796498101822e-03
best_mse[:Bomex]["moisture.ρq_tot"] = 4.8463531712319707e-02
best_mse[:Bomex]["turbconv.environment.ρatke"] = 6.6626422991098286e+02
best_mse[:Bomex]["turbconv.environment.ρaθ_liq_cv"] = 8.5667099987715503e+01
best_mse[:Bomex]["turbconv.environment.ρaq_tot_cv"] = 1.6439116552012851e+02
best_mse[:Bomex]["turbconv.updraft[1].ρa"] = 7.9577348413012515e+01
best_mse[:Bomex]["turbconv.updraft[1].ρaw"] = 8.4352188057391225e-02
best_mse[:Bomex]["turbconv.updraft[1].ρaθ_liq"] = 9.0101464252959325e+00
best_mse[:Bomex]["turbconv.updraft[1].ρaq_tot"] = 1.0768120864370509e+01
#! format: on

sufficient_mse(computed_mse, best_mse) = computed_mse <= best_mse + eps()

function test_mse(computed_mse, best_mse, key)
    mse_not_regressed = sufficient_mse(computed_mse[key], best_mse[key])
    @test mse_not_regressed
    mse_not_regressed || @show key
end

computed_mse = Dict(
    k => compute_mse(
        solver_config.dg.grid,
        solver_config.dg.balance_law,
        time_data,
        dons_arr,
        data_files[k],
        k,
        best_mse[k],
        plot_dir,
    ) for k in keys(data_files)
)

@testset "BOMEX EDMF Solution Quality Assurance (QA) tests" begin
    #! format: off
    test_mse(computed_mse[:Bomex], best_mse[:Bomex], "ρ")
    test_mse(computed_mse[:Bomex], best_mse[:Bomex], "ρu[1]")
    test_mse(computed_mse[:Bomex], best_mse[:Bomex], "moisture.ρq_tot")
    test_mse(computed_mse[:Bomex], best_mse[:Bomex], "turbconv.updraft[1].ρa")
    test_mse(computed_mse[:Bomex], best_mse[:Bomex], "turbconv.updraft[1].ρaw")
    test_mse(computed_mse[:Bomex], best_mse[:Bomex], "turbconv.updraft[1].ρaθ_liq")
    test_mse(computed_mse[:Bomex], best_mse[:Bomex], "turbconv.updraft[1].ρaq_tot")
    test_mse(computed_mse[:Bomex], best_mse[:Bomex], "turbconv.environment.ρatke")
    test_mse(computed_mse[:Bomex], best_mse[:Bomex], "turbconv.environment.ρaθ_liq_cv")
    test_mse(computed_mse[:Bomex], best_mse[:Bomex], "turbconv.environment.ρaq_tot_cv")
    #! format: on
end

max_stars_repo_name: jm-c/CLIMA

max_issues_repo_name: jm-c/CLIMA

max_forks_repo_name: jm-c/CLIMA


text:
module adapters_postgresql_structured_test

using Test # @test
using Octo.Adapters.PostgreSQL # Repo.connect Schema from SELECT FROM WHERE INNER JOIN ON GROUP BY

struct Weather
end

struct City
end

Schema.model(Weather, table_name="weather")
Schema.model(City, table_name="cities")

w = from(Weather)
c = from(City)

@test to_sql([SELECT * FROM (w, c) WHERE w.city == c.name]) == "SELECT * FROM weather, cities WHERE city = name"
@test to_sql([SELECT * FROM w INNER JOIN c ON w.city == c.name]) == "SELECT * FROM weather INNER JOIN cities ON city = name"

##

struct Distributor
end

struct Film
end

Schema.model(Distributor, table_name="distributors")
Schema.model(Film, table_name="films")

d = from(Distributor, :d)
f = from(Film, :f)

@test to_sql([SELECT (f.title,) FROM (d, f)]) ==
             "SELECT f.title FROM distributors d, films f"

@test to_sql([SELECT (f.title, f.did, d.name, f.date_prod, f.kind) FROM (d, f) WHERE f.did == d.did]) ==
             "SELECT f.title, f.did, d.name, f.date_prod, f.kind FROM distributors d, films f WHERE f.did = d.did"

d = from(Distributor)
f = from(Film)
total = as(SUM(f.len), :total)
@test to_sql([SELECT (f.kind, total) FROM f GROUP BY :kind]) ==
             "SELECT kind, SUM(len) AS total FROM films GROUP BY kind"
@test to_sql([SELECT d.name FROM d WHERE d.name LIKE "W%"]) ==
             "SELECT name FROM distributors WHERE name LIKE 'W%'"
@test to_sql([SELECT * INTO :films_recent FROM f WHERE f.date_prod >= "2002-01-01"]) ==
             "SELECT * INTO films_recent FROM films WHERE date_prod >= '2002-01-01'"

# issue 10
struct Here
end
Schema.model(Here, table_name="heres", primary_key="id")
heres = from(Here, :heres)
@test to_sql([INSERT INTO heres VALUES ("test1", "SR")]) == "INSERT INTO heres VALUES ('test1', 'SR')"

end # module adapters_postgresql_structured_test


text:
using LightXML
using Dates
using Statistics
using Plots

const datadir = joinpath(dirname(pathof(Glucose)), "..", "data", "OhioT1DM-training")
const ohio5559 = joinpath(datadir, "559-ws-training.xml")

function average_over_bins(data::AbstractArray, bin_size::Int64)
  size = min(length(data), bin_size)
  bins = Iterators.partition(data, size) |> collect
  map(row -> mean(row), bins)
end

function normalize(data::AbstractArray)
  (data .- mean(data)) ./ std(data)
end

function normalize_binned_data(data::AbstractArray, bin_size::Int64)
  normalize(average_over_bins(data, bin_size))
end

function prepare_data(var::String, bin_size::Int64; file_location = ohio5559)
  if var == "basis_steps"
    ode_data = prepare_steps_data(file_location = file_location)
  elseif var == "bolus"
    ode_data = prepare_bolus_data(file_location = file_location)
  end

  println(length(ode_data[1, :]))
  println(length(ode_data[2, :]))

  glucose_data = normalize_binned_data(ode_data[1, :], bin_size)
  other_data = normalize_binned_data(ode_data[2, :], bin_size)

  println(length(glucose_data))
  println(length(other_data))

  u0 = [glucose_data[1]; other_data[1]]

  (u0, Array(transpose(hcat([glucose_data, other_data]...))))
end

function prepare_steps_data(; file_location = ohio5559)
  full_data_glucose = parse_data_from_xml("glucose_level"; round_time = true, file_location = file_location)
  full_data_steps = parse_data_from_xml("basis_steps"; round_time = true, file_location = file_location)

  times_glucose, ode_data_glucose = full_data_glucose[1, :], full_data_glucose[2, :] # 1000 is arbitrary
  times_steps, ode_data_steps = full_data_steps[1, :], full_data_steps[2, :] # 1000 is arbitrary

  sorted_segments = sorted_glucose_steps_segments()
  times = sorted_segments[1]
  start_time = minimum(times)
  times = map(time -> time - start_time, times)

  ode_data_glucose = map(i -> ode_data_glucose[i], findall(time -> (time - start_time) in times, unique(times_glucose)))
  ode_data_steps = map(i -> ode_data_steps[i], findall(time -> (time - start_time) in times, unique(times_steps)))

  Array(transpose(hcat([ode_data_glucose, ode_data_steps]...)))
end

function prepare_all_data(bin_size::Int64; file_location = ohio5559)
  full_data_glucose = parse_data_from_xml("glucose_level"; round_time = true, file_location = file_location)
  full_data_bolus = parse_data_from_xml("bolus"; round_time = true, file_location = file_location)
  full_data_steps = parse_data_from_xml("basis_steps"; round_time = true, file_location = file_location)

  times_glucose, ode_data_glucose = full_data_glucose[1, :], full_data_glucose[2, :]
  times_bolus, ode_data_bolus = full_data_bolus[1, :], full_data_bolus[2, :]
  times_steps, ode_data_steps = full_data_steps[1, :], full_data_steps[2, :] # 1000 is arbitrary

  sorted_segments = sorted_glucose_steps_segments()

  sorted_segment_intersections = map(segment -> length(intersect(segment, times_bolus)), sorted_segments)
  idx = findmax(sorted_segment_intersections)[2]
  times = sorted_segments[idx]

  start_time = minimum(times)
  times = map(time -> time - start_time, times)

  ode_data_glucose = normalize_binned_data(map(i -> ode_data_glucose[i], findall(time -> (time - start_time) in times, times_glucose)), bin_size)
  ode_data_steps = normalize_binned_data(map(i -> ode_data_steps[i], findall(time -> (time - start_time) in times, times_steps)), bin_size)

  ode_data_bolus_new = []
  for time in times
    if (time + start_time) in times_bolus
      idx = findall(x -> x == (time + start_time), unique(times_bolus))[1]
      push!(ode_data_bolus_new, ode_data_bolus[idx])
    else
      push!(ode_data_bolus_new, 0.0)
    end
  end
  ode_data_bolus = normalize_binned_data(ode_data_bolus_new, bin_size)

  ode_data = Array(transpose(hcat([ode_data_glucose, ode_data_steps, ode_data_bolus]...)))
  u0 = ode_data[:, 1]

  u0, ode_data
end

function prepare_all_data_meals_hypo(bin_size::Int64; file_location = ohio5559, hypo_id::Int = 2, normal = true)
  full_data_glucose = parse_data_from_xml("glucose_level"; round_time = true, file_location = file_location)
  full_data_bolus = parse_data_from_xml("bolus"; round_time = true, file_location = file_location)
  full_data_steps = parse_data_from_xml("basis_steps"; round_time = true, file_location = file_location)
  full_data_meals = parse_data_from_xml("meal"; round_time = true, file_location = file_location)
  full_data_hypoevents = parse_data_from_xml("hypo_event"; round_time = true, file_location = file_location)

  # set start time to zero
  offset_start_time = minimum(vcat(full_data_glucose[1, :],
    full_data_bolus[1, :],
    full_data_steps[1, :],
    full_data_meals[1, :],
    full_data_hypoevents[1, :]))

  full_data_glucose[1, :] = full_data_glucose[1, :] .- offset_start_time
  full_data_bolus[1, :] = full_data_bolus[1, :] .- offset_start_time
  full_data_steps[1, :] = full_data_steps[1, :] .- offset_start_time
  full_data_meals[1, :] = full_data_meals[1, :] .- offset_start_time
  full_data_hypoevents[1, :] = full_data_hypoevents[1, :] .- offset_start_time

  times_glucose, ode_data_glucose = full_data_glucose[1, :], full_data_glucose[2, :]
  times_bolus, ode_data_bolus = full_data_bolus[1, :], full_data_bolus[2, :]
  times_steps, ode_data_steps = full_data_steps[1, :], full_data_steps[2, :] # 1000 is arbitrary
  times_meals, ode_data_meals = full_data_meals[1, :], full_data_meals[2, :]
  times_hypoevents, ode_data_hypoevents = full_data_hypoevents[1, :], full_data_hypoevents[2, :]

  sorted_segments = sorted_glucose_steps_segments()
  intersections = map(x -> (length(intersect(x, times_bolus)), length(intersect(x, times_meals)), length(intersect(x, times_hypoevents))), sorted_segments)
  at_least_one_hypoevent_indices = findall(x -> x[3] != 0, intersections)

  times = sorted_segments[at_least_one_hypoevent_indices[hypo_id]]

  start_time = minimum(times)
  times = map(time -> time - start_time, times)
  println("START_TIME: ", start_time)
  if normal
    ode_data_glucose = normalize_binned_data(map(i -> ode_data_glucose[i], findall(time -> (time - start_time) in times, times_glucose)), bin_size)
    ode_data_steps = normalize_binned_data(map(i -> ode_data_steps[i], findall(time -> (time - start_time) in times, times_steps)), bin_size)
  else
    ode_data_glucose = map(i -> ode_data_glucose[i], findall(time -> (time - start_time) in times, times_glucose))
    ode_data_steps = map(i -> ode_data_steps[i], findall(time -> (time - start_time) in times, times_steps))
  end

  ode_data_bolus_new = []
  ode_data_meals_new = []
  ode_data_hypoevents_new = []
  for time in times
    if (time + start_time) in times_bolus
      idx = findall(x -> x == (time + start_time), times_bolus)[1]
      idx2 = findall(x -> x == (time + start_time), unique(times_bolus))[1]

      if idx != idx2
        println("HERE")
      end

      push!(ode_data_bolus_new, ode_data_bolus[idx])
    else
      push!(ode_data_bolus_new, 0.0)
    end

    if (time + start_time) in times_meals
      idx = findall(x -> x == (time + start_time), times_meals)[1]
      push!(ode_data_meals_new, ode_data_meals[idx])
    else
      push!(ode_data_meals_new, 0.0)
    end

    if (time + start_time) in times_hypoevents
      idx = findall(x -> x == (time + start_time), times_hypoevents)[1]
      push!(ode_data_hypoevents_new, ode_data_hypoevents[idx])
    else
      push!(ode_data_hypoevents_new, 0.0)
    end
  end

  if normal
    ode_data_bolus = normalize_binned_data(ode_data_bolus_new, bin_size)
    ode_data_meals = normalize_binned_data(ode_data_meals_new, bin_size)
    ode_data_hypoevents = normalize_binned_data(ode_data_hypoevents_new, bin_size)
  else
    ode_data_bolus = ode_data_bolus_new
    ode_data_meals = ode_data_meals_new
    ode_data_hypoevents = ode_data_hypoevents_new
  end

  ode_data = Array(transpose(hcat([ode_data_glucose, ode_data_steps, ode_data_bolus, ode_data_meals, ode_data_hypoevents]...)))
  u0 = ode_data[:, 1]

  u0, ode_data
end

function prepare_bolus_data(; file_location = ohio5559)
  full_data_glucose = parse_data_from_xml("glucose_level"; round_time = true, file_location = file_location)
  full_data_bolus = parse_data_from_xml("bolus"; round_time = true, file_location = file_location)

  times_glucose, ode_data_glucose = full_data_glucose[1, :], full_data_glucose[2, :]
  times_bolus, ode_data_bolus = full_data_bolus[1, :], full_data_bolus[2, :]

  times = times_glucose[5183:5999]
  start_time = minimum(times)
  times = map(time -> time - start_time, times)

  ode_data_glucose = ode_data_glucose[5183:5999]
  ode_data_bolus_new = []

  for time in times
    if (time + start_time) in times_bolus
      idx = findall(x -> x == (time + start_time), unique(times_bolus))[1]
      push!(ode_data_bolus_new, ode_data_bolus[idx])
    else
      push!(ode_data_bolus_new, 0.0)
    end
  end

  ode_data_bolus = ode_data_bolus_new

  println(length(ode_data_glucose))
  println(length(ode_data_bolus))

  Array(transpose(hcat([ode_data_glucose, ode_data_bolus]...)))

end

function sorted_glucose_steps_segments(; file_location = ohio5559)
  full_data_glucose = parse_data_from_xml("glucose_level"; round_time = true, file_location = file_location)
  full_data_bolus = parse_data_from_xml("bolus"; round_time = true, file_location = file_location)
  full_data_steps = parse_data_from_xml("basis_steps"; round_time = true, file_location = file_location)
  full_data_meals = parse_data_from_xml("meal"; round_time = true, file_location = file_location)
  full_data_hypoevents = parse_data_from_xml("hypo_event"; round_time = true, file_location = file_location)

  # set start time to zero
  offset_start_time = minimum(vcat(full_data_glucose[1, :],
    full_data_bolus[1, :],
    full_data_steps[1, :],
    full_data_meals[1, :],
    full_data_hypoevents[1, :]))

  full_data_glucose[1, :] = full_data_glucose[1, :] .- offset_start_time
  full_data_bolus[1, :] = full_data_bolus[1, :] .- offset_start_time
  full_data_steps[1, :] = full_data_steps[1, :] .- offset_start_time
  full_data_meals[1, :] = full_data_meals[1, :] .- offset_start_time
  full_data_hypoevents[1, :] = full_data_hypoevents[1, :] .- offset_start_time

  times_glucose, ode_data_glucose = full_data_glucose[1, :], full_data_glucose[2, :]
  times_bolus, ode_data_bolus = full_data_bolus[1, :], full_data_bolus[2, :]
  times_steps, ode_data_steps = full_data_steps[1, :], full_data_steps[2, :] # 1000 is arbitrary

  times_gs = intersect(times_glucose, times_steps)
  shifted_times_gs = [0.0, times_gs[1:length(times_gs)-1]...][1:length(times_gs)]
  consec_time_differences = (times_gs.-shifted_times_gs)[2:length(times_gs)]
  jump_indices = findall(x -> !(x in [0.0, 5.0]), consec_time_differences)

  segments = []
  for i in 1:length(jump_indices)
    if i == 1
      push!(segments, times_gs[1:jump_indices[i]])
    else
      push!(segments, times_gs[(jump_indices[i-1]+1):jump_indices[i]])
    end
  end

  reverse(sort(segments, by = length))
end

function parse_data_from_xml(data_name::String; round_time::Bool = false, file_location = ohio5559)
  # xdoc is an instance of XMLDocument, which maintains a tree structure
  xdoc = parse_file(file_location)

  # get the root element
  xroot = root(xdoc)  # an instance of XMLElement

  data_element = find_element(xroot, data_name)

  parsed_data = [[], []]

  for event in child_nodes(data_element)
    if is_elementnode(event)
      # println(event)
      if (data_name in ["glucose_level", "basis_steps"])
        push!(parsed_data[1], ceil(DateTime(attribute(XMLElement(event), "ts"), "dd-mm-yyyy HH:MM:SS"), Dates.Minute(round_time ? 5 : 1)).instant.periods.value / (1000.0 * 60))
        push!(parsed_data[2], parse(Float64, attribute(XMLElement(event), "value")))
      elseif (data_name == "meal")
        push!(parsed_data[1], ceil(DateTime(attribute(XMLElement(event), "ts"), "dd-mm-yyyy HH:MM:SS"), Dates.Minute(round_time ? 5 : 1)).instant.periods.value / (1000.0 * 60))
        push!(parsed_data[2], parse(Float64, attribute(XMLElement(event), "carbs")))
      elseif (data_name == "bolus")
        push!(parsed_data[1], ceil(DateTime(attribute(XMLElement(event), "ts_begin"), "dd-mm-yyyy HH:MM:SS"), Dates.Minute(round_time ? 5 : 1)).instant.periods.value / (1000.0 * 60))
        push!(parsed_data[2], parse(Float64, attribute(XMLElement(event), "dose")))
      elseif (data_name == "hypo_event")
        push!(parsed_data[1], ceil(DateTime(attribute(XMLElement(event), "ts"), "dd-mm-yyyy HH:MM:SS"), Dates.Minute(round_time ? 5 : 1)).instant.periods.value / (1000.0 * 60))
        push!(parsed_data[2], 1.0)
      end
    end
  end

  free(xdoc)

  println(string("length(parsed_data): ", length(parsed_data[1])))

  # shift times to begin at t=0
  start_time = minimum(parsed_data[1])
  # parsed_data = (map(time -> time - start_time, parsed_data[1]), parsed_data[2])

  # reshape into matrix
  Array(transpose(hcat(parsed_data...)))
end

function plot_data(ode_data::AbstractArray, var1::String, var2::String)
  tspan = (0.0, Float64(length(ode_data[1, :]) - 1))
  tsteps = range(tspan[1], tspan[2], length = length(ode_data[1, :]))

  pl = plot(tsteps, ode_data[1, :], color = :red, label = "Data: $(var1)", xlabel = "t", title = "$(var1) and $(var2)")
  plot!(tsteps, ode_data[2, :], color = :blue, label = "Data: $(var2)")
  pl
end

function plot_data_all(ode_data::AbstractArray, var1::String, var2::String)
  tspan = (0.0, Float64(length(ode_data[1, :]) - 1))
  tsteps = range(tspan[1], tspan[2], length = length(ode_data[1, :]))

  pl = plot(tsteps, ode_data[1, :], color = :red, label = "Data: Glucose", xlabel = "t", title = "Glucose, $(var1), and $(var2)")
  plot!(tsteps, ode_data[2, :], color = :blue, label = "Data: $(var1)")
  plot!(tsteps, ode_data[3, :], color = :green, label = "Data: $(var2)")
  pl
end

function plot_data_all_meals_hypo(ode_data::AbstractArray)
  tspan = (0.0, Float64(length(ode_data[1, :]) - 1))
  tsteps = range(tspan[1], tspan[2], length = length(ode_data[1, :]))

  pl = plot(tsteps, ode_data[1, :], color = :red, label = "Data: Glucose", xlabel = "t", title = "Glucose, Steps, Bolus, Meals, Hypo Events")
  plot!(tsteps, ode_data[2, :], color = :blue, label = "Data: Steps")
  plot!(tsteps, ode_data[3, :], color = :green, label = "Data: Bolus")
  plot!(tsteps, ode_data[4, :], color = :purple, label = "Data: Meals")
  plot!(tsteps, ode_data[5, :], color = :orange, label = "Data: Hypo Events")
  pl
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect