text:
export ITensor,
       norm,
       combiner,
       combinedindex,
       delta,
       δ,
       exp,
       replaceindex!,
       inds,
       isNull,
       scale!,
       normalize!,
       randn!,
       multSiteOps,
       order,
       permute,
       randomITensor,
       diagITensor,
       tensor,
       array,
       matrix,
       vector,
       scalar,
       store,
       dense

"""
An ITensor is a tensor whose interface is
independent of its memory layout. Therefore
it is not necessary to know the ordering
of an ITensor's indices, only which indices
an ITensor has. Operations like contraction
and addition of ITensors automatically
handle any memory permutations.
"""
mutable struct ITensor{N}
  store::TensorStorage
  inds::IndexSet{N}
  #TODO: check that the storage is consistent with the
  #total dimension of the indices (possibly only in debug mode);
  ITensor{N}(st,is::IndexSet{N}) where {N} = new{N}(st,is)
end
ITensor(st,is::IndexSet{N}) where {N} = ITensor{N}(st,is)

inds(T::ITensor) = T.inds
store(T::ITensor) = T.store

# TODO: do we need these? I think yes, for add!(::ITensor,::ITensor)
setinds!(T::ITensor,is...) = (T.inds = IndexSet(is...))
setstore!(T::ITensor,st::TensorStorage) = (T.store = st)

#
# Dense ITensor constructors
#

ITensor(T::Tensor{<:Number,N}) where {N} = ITensor{N}(store(T),inds(T))
ITensor{N}(T::Tensor{<:Number,N}) where {N} = ITensor{N}(store(T),inds(T))

# Convert the ITensor to a Tensor that shares the same
# data and indices as the ITensor
# TODO: should we define a `convert(::Type{<:Tensor},::ITensor)`
# method?
tensor(A::ITensor) = Tensor(store(A),inds(A))

"""
    ITensor(iset::IndexSet)

Construct an ITensor having indices
given by the IndexSet `iset`
"""
ITensor(is::IndexSet) = ITensor(Float64,is)
ITensor() = ITensor{0}(Dense{Nothing}(),IndexSet())
ITensor(inds::Vararg{Index,N}) where {N} = ITensor(IndexSet{N}(inds...))

function ITensor(::Type{T},
                 inds::IndexSet{N}) where {T<:Number,N}
  return ITensor{N}(Dense{float(T)}(zeros(float(T),dim(inds))),inds)
end
ITensor(::Type{T},inds::Index...) where {T<:Number} = ITensor(T,IndexSet(inds...))

function ITensor(::UndefInitializer,
                 inds::IndexSet{N}) where {N}
  return ITensor{N}(Dense{Float64}(Vector{Float64}(undef,dim(inds))),inds)
end
ITensor(x::UndefInitializer,inds::Index...) = ITensor(x,IndexSet(inds...))

function ITensor(x::S,inds::IndexSet{N}) where {S<:Number,N}
  return ITensor{N}(Dense{float(S)}(fill(float(x),dim(inds))),inds)
end

"""
    ITensor(x)

Construct a scalar ITensor with value `x`.

    ITensor(x,i,j,...)

Construct an ITensor with indices `i`,`j`,...
and all elements set to `x`.
"""
ITensor(x::S,inds::Index...) where {S<:Number} = ITensor(x,IndexSet(inds...))

function ITensor(A::Array{S},inds::IndexSet{N}) where {S<:Number,N}
  length(A) ≠ dim(inds) && throw(DimensionMismatch("In ITensor(Array,IndexSet), length of Array ($(length(A))) must match total dimension of IndexSet ($(dim(inds)))"))
  return ITensor{N}(Dense{float(S)}(float(vec(A))),inds)
end
ITensor(A::Array{S},inds::Index...) where {S<:Number} = ITensor(A,IndexSet(inds...))

#
# Diag ITensor constructors
#

"""
diagITensor(::Type{T}, is::IndexSet)

Make a sparse ITensor of element type T with non-zero elements
only along the diagonal. Defaults to having `zero(T)` along the diagonal.
The storage will have Diag type.
"""
function diagITensor(::Type{T},
                     is::IndexSet{N}) where {T<:Number,N}
  return ITensor{N}(Diag{Vector{T}}(zeros(T,minDim(is))),is)
end

"""
diagITensor(::Type{T}, is::Index...)

Make a sparse ITensor of element type T with non-zero elements
only along the diagonal. Defaults to having `zero(T)` along the diagonal.
The storage will have Diag type.
"""
diagITensor(::Type{T},inds::Index...) where {T<:Number} = diagITensor(T,IndexSet(inds...))

"""
diagITensor(v::Vector{T}, is::IndexSet)

Make a sparse ITensor with non-zero elements only along the diagonal.
The diagonal elements will be set to the values stored in `v` and
the ITensor will have element type `float(T)`.
The storage will have Diag type.
"""
function diagITensor(v::Vector{T},
                     is::IndexSet) where {T<:Number}
  length(v) ≠ minDim(is) && error("Length of vector for diagonal must equal minimum of the dimension of the input indices")
  return ITensor(Diag{Vector{float(T)}}(v),is)
end

"""
diagITensor(v::Vector{T}, is::Index...)

Make a sparse ITensor with non-zero elements only along the diagonal.
The diagonal elements will be set to the values stored in `v` and
the ITensor will have element type `float(T)`.
The storage will have Diag type.
"""
function diagITensor(v::Vector{T},
                     is::Index...) where {T<:Number}
  return diagITensor(v,IndexSet(is...))
end

"""
diagITensor(is::IndexSet)

Make a sparse ITensor of element type Float64 with non-zero elements
only along the diagonal. Defaults to storing zeros along the diagonal.
The storage will have Diag type.
"""
diagITensor(is::IndexSet) = diagITensor(Float64,is)

"""
diagITensor(is::Index...)

Make a sparse ITensor of element type Float64 with non-zero elements
only along the diagonal. Defaults to storing zeros along the diagonal.
The storage will have Diag type.
"""
diagITensor(inds::Index...) = diagITensor(IndexSet(inds...))

"""
diagITensor(x::T, is::IndexSet) where {T<:Number}

Make a sparse ITensor with non-zero elements only along the diagonal.
The diagonal elements will be set to the value `x` and
the ITensor will have element type `float(T)`.
The storage will have Diag type.
"""
function diagITensor(x::T,
                     is::IndexSet) where {T<:Number}
  return ITensor(Diag{Vector{float(T)}}(fill(float(x),minDim(is))),is)
end

"""
diagITensor(x::T, is::Index...) where {T<:Number}

Make a sparse ITensor with non-zero elements only along the diagonal.
The diagonal elements will be set to the value `x` and
the ITensor will have element type `float(T)`.
The storage will have Diag type.
"""
function diagITensor(x::T,
                     is::Index...) where {T<:Number}
  return diagITensor(x,IndexSet(is...))
end

"""
    delta(::Type{T},inds::IndexSet)

Make a diagonal ITensor with all diagonal elements 1.
"""
function delta(::Type{T},is::IndexSet) where {T}
  return ITensor(Diag{float(T)}(one(T)),is)
end

"""
    delta(::Type{T},inds::Index...)

Make a diagonal ITensor with all diagonal elements 1.
"""
function delta(::Type{T},is::Index...) where {T}
  return delta(T,IndexSet(is...))
end

delta(is::IndexSet) = delta(Float64,is)

delta(is::Index...) = delta(IndexSet(is...))
const δ = delta

"""
dense(T::ITensor)

Make a copy of the ITensor where the storage is the dense version.
For example, an ITensor with Diag storage will become Dense storage.
"""
function dense(T::ITensor)
  ITensor(dense(Tensor(store(T),inds(T))))
end

"""
complex(T::ITensor)

Convert to the complex version of the storage.
"""
Base.complex(T::ITensor) = ITensor(complex(Tensor(store(T),inds(T))))

# This constructor allows many IndexSet
# set operations to work with ITensors
IndexSet(T::ITensor) = inds(T)

eltype(T::ITensor) = eltype(store(T))

"""
    order(A::ITensor) = ndims(A)

The number of indices, `length(inds(A))`.
"""
order(T::ITensor) = order(inds(T))
ndims(T::ITensor) = order(inds(T))

"""
    dim(A::ITensor) = length(A)

The total number of entries, `prod(size(A))`.
"""
dim(T::ITensor) = dim(inds(T))

"""
    dims(A::ITensor) = size(A)

Tuple containing `size(A,d) == dim(inds(A)[d]) for d in 1:ndims(A)`.
"""
dims(T::ITensor) = dims(inds(T))
Base.size(A::ITensor) = dims(inds(A))
Base.size(A::ITensor{N}, d::Int) where {N} = d in 1:N ? dim(inds(A)[d]) :
  d>0 ? 1 : error("arraysize: dimension out of range")

isNull(T::ITensor) = (eltype(T) === Nothing)

Base.copy(T::ITensor{N}) where {N} = ITensor{N}(copy(tensor(T)))

# TODO: make versions where the element type can be specified
# Should this be called `array`? (Version that makes a view, if dense)
array(T::ITensor) = array(tensor(T))

matrix(T::ITensor{N}) where {N} = (N==2 ? array(tensor(T)) : throw(DimensionMismatch("ITensor must be order 2 to convert to a Matrix")))

vector(T::ITensor{N}) where {N} = (N==1 ? array(tensor(T)) : throw(DimensionMismatch("ITensor must be order 1 to convert to a Vector")))

scalar(T::ITensor) = T[]

function array(T::ITensor{N},is::Vararg{Index,N}) where {N}
  perm = getperm(inds(T),is)
  return array(permutedims(tensor(T),perm))
end

function matrix(T::ITensor{N},i1::Index,i2::Index) where {N}
  N≠2 && throw(DimensionMismatch("ITensor must be order 2 to convert to a Matrix"))
  return array(T,i1,i2)
end

Base.getindex(T::ITensor{N},vals::Vararg{Int,N}) where {N} = tensor(T)[vals...]::Number

function Base.getindex(T::ITensor{N},
                       ivs::Vararg{IndexVal,N}) where {N}
  p = getperm(inds(T),ivs)
  vals = permute(val.(ivs),p)
  return T[vals...]
end

# TODO: we should figure out if this is how we want to do
# slicing
#function getindex(T::ITensor,ivs::AbstractVector{IndexVal}...)
#  p = getperm(inds(T),map(x->x[1], ivs))
#  vals = map(x->val.(x), ivs[[p...]])
#  return Tensor(store(T),inds(T))[vals...]
#end

Base.getindex(T::ITensor) = tensor(T)[]

Base.setindex!(T::ITensor{N},x::Number,vals::Vararg{Int,N}) where {N} = (tensor(T)[vals...] = x)

function Base.setindex!(T::ITensor,x::Number,ivs::IndexVal...)
  p = getperm(inds(T),ivs)
  vals = permute(val.(ivs),p)
  return T[vals...] = x
end

# TODO: what was this doing?
#function setindex!(T::ITensor,
#                   x::Union{<:Number, AbstractArray{<:Number}},
#                   ivs::Union{IndexVal, AbstractVector{IndexVal}}...)
#  remap_ivs = map(x->x isa IndexVal ? x : x[1], ivs)
#  p = getperm(inds(T),remap_ivs)
#  vals = map(x->x isa IndexVal ? val(x) : val.(x), ivs[p])
#  storage_setindex!(store(T),inds(T),x,vals...)
#  return T
#end

function Base.fill!(T::ITensor,
                    x::Number)
  # TODO: automatically switch storage type if needed?
  fill!(tensor(T),x)
  return T
end

# TODO: implement in terms of delta tensors (better for QNs)
function replaceindex!(A::ITensor,i::Index,j::Index)
  pos = indexpositions(A,i)
  inds(A)[pos[1]] = j
  return A
end

function replaceinds!(A::ITensor,inds1,inds2)
  is1 = IndexSet(inds1)
  is2 = IndexSet(inds2)
  pos = indexpositions(A,is1)
  for (j,p) ∈ enumerate(pos)
    inds(A)[p] = is2[j]
  end
  return A
end

# TODO: can we turn these into a macro?
# for f ∈ (prime,setprime,noprime,...)
#   f(A::ITensor,vargs...) = ITensor(store(A),f(inds(A),vargs...))
#   f!(A::ITensor,vargs...) = ( f!(inds(A),vargs...); return A )
# end

# TODO: implement more in-place versions

prime!(A::ITensor,vargs...)= ( prime!(inds(A),vargs...); return A )
prime(A::ITensor,vargs...)= ITensor(store(A),prime(inds(A),vargs...))
Base.adjoint(A::ITensor) = prime(A)

setprime(A::ITensor,vargs...) = ITensor(store(A),setprime(inds(A),vargs...))

noprime(A::ITensor,vargs...) = ITensor(store(A),noprime(inds(A),vargs...))

mapprime(A::ITensor,vargs...) = ITensor(store(A),mapprime(inds(A),vargs...))

swapprime(A::ITensor,vargs...) = ITensor(store(A),swapprime(inds(A),vargs...))

addtags(A::ITensor,vargs...) = ITensor(store(A),addtags(inds(A),vargs...))

removetags(A::ITensor,vargs...) = ITensor(store(A),removetags(inds(A),vargs...))

replacetags!(A::ITensor,vargs...) = ( replacetags!(inds(A),vargs...); return A )
replacetags(A::ITensor,vargs...) = ITensor(store(A),replacetags(inds(A),vargs...))

settags!(A::ITensor,vargs...) = ( settags!(inds(A),vargs...); return A )
settags(A::ITensor,vargs...) = ITensor(store(A),settags(inds(A),vargs...))

swaptags(A::ITensor,vargs...) = ITensor(store(A),swaptags(inds(A),vargs...))

# TODO: implement in a better way (more generically for other storage)
Base.:(==)(A::ITensor,B::ITensor) = (norm(A-B) == zero(promote_type(eltype(A),eltype(B))))

function Base.isapprox(A::ITensor,
                       B::ITensor;
                       atol::Real=0.0,
                       rtol::Real=Base.rtoldefault(eltype(A),eltype(B),atol))
    return norm(A-B) <= atol + rtol*max(norm(A),norm(B))
end

# TODO: bring this back or just use T[]?
# I think T[] may not be generic, since it may only work if order(T)==0
# so it may be nice to have a seperate scalar(T) for when dim(T)==1
#function scalar(T::ITensor)
#  !(order(T)==0 || dim(T)==1) && throw(ArgumentError("ITensor with inds $(inds(T)) is not a scalar"))
#  return scalar(tensor(store(T),inds(T)))
#end

function Random.randn!(T::ITensor)
  return randn!(tensor(T))
end

const Indices = Union{IndexSet,Tuple{Vararg{Index}}}

"""
    randomITensor([S,] inds)

Construct an ITensor with type S (default Float64) and indices inds, whose elements are normally distributed random numbers.
"""
function randomITensor(::Type{S},
                       inds::Indices) where {S<:Number}
  T = ITensor(S,IndexSet(inds))
  randn!(T)
  return T
end
function randomITensor(::Type{S},
                       inds::Index...) where {S<:Number}
  return randomITensor(S,IndexSet(inds...))
end
randomITensor(inds::Indices) = randomITensor(Float64,
                                             IndexSet(inds))
randomITensor(inds::Index...) = randomITensor(Float64,
                                              IndexSet(inds...))

function combiner(inds::IndexSet; kwargs...)
  tags = get(kwargs, :tags, "CMB,Link")
  new_ind = Index(prod(dims(inds)), tags)
  new_is = IndexSet(new_ind, inds)
  return ITensor(Combiner(),new_is),new_ind
end
combiner(inds::Index...; kwargs...) = combiner(IndexSet(inds...); kwargs...)

combinedindex(T::ITensor) = store(T) isa Combiner ? store(T).ci : nothing

LinearAlgebra.norm(T::ITensor) = norm(tensor(T))

function dag(T::ITensor)
  TT = conj(tensor(T))
  return ITensor(store(TT),dag(inds(T)))
end

function permute(T::ITensor,new_inds)
  perm = getperm(new_inds,inds(T))
  Tp = permutedims(tensor(T),perm)
  return ITensor(Tp)
end
permute(T::ITensor,inds::Index...) = permute(T,IndexSet(inds...))

function Base.:*(T::ITensor,x::Number)
  return ITensor(x*tensor(T))
end
Base.:*(x::Number,T::ITensor) = T*x
#TODO: make a proper element-wise division
Base.:/(A::ITensor,x::Number) = A*(1.0/x)

Base.:-(A::ITensor) = ITensor(-tensor(A))
function Base.:+(A::ITensor,B::ITensor)
  C = copy(A)
  C = add!(C,B)
  return C
end
function Base.:-(A::ITensor,B::ITensor)
  C = copy(A)
  C = add!(C,-1,B)
  return C
end

function Base.:*(A::ITensor,B::ITensor)
  (Alabels,Blabels) = compute_contraction_labels(inds(A),inds(B))
  CT = contract(tensor(A),Alabels,tensor(B),Blabels)
  C = ITensor(CT)
  if warnTensorOrder > 0 && order(C) >= warnTensorOrder
    println("Warning: contraction resulted in ITensor with $(order(C)) indices")
  end
  return C
end

LinearAlgebra.dot(A::ITensor,B::ITensor) = (dag(A)*B)[]

"""
    exp(A::ITensor, Lis::IndexSet; hermitian = false)

Compute the exponential of the tensor `A` by treating it as a matrix ``A_{lr}`` with
the left index `l` running over all indices in `Lis` and `r` running over all
indices not in `Lis`. Must have `dim(Lis) == dim(inds(A))/dim(Lis)` for the exponentiation to
be defined.
When `hermitian=true` the exponential of `Hermitian(reshape(A, dim(Lis), :))` is
computed internally.
"""
function LinearAlgebra.exp(A::ITensor,
                           Linds,
                           Rinds = prime(IndexSet(Linds));
                           ishermitian = false)
  Lis,Ris = IndexSet(Linds),IndexSet(Rinds)
  Lpos,Rpos = getperms(inds(A),Lis,Ris)
  expAT = exp(tensor(A),Lpos,Rpos;ishermitian=ishermitian)
  return ITensor(expAT)
end

expHermitian(A::ITensor,Linds,Rinds) = exp(A,Linds,Rinds;ishermitian=true)

#######################################################################
#
# In-place operations
#

"""
    normalize!(T::ITensor)

Normalize an ITensor in-place, such that norm(T)==1.
"""
normalize!(T::ITensor) = scale!(T,1/norm(T))

"""
    copyto!(B::ITensor, A::ITensor)

Copy the contents of ITensor A into ITensor B.
```
B .= A
```
"""
function Base.copyto!(R::ITensor{N},T::ITensor{N}) where {N}
  perm = getperm(inds(R),inds(T))
  TR = permutedims!(tensor(R),tensor(T),perm)
  return ITensor(TR)
end

"""
    add!(B::ITensor, A::ITensor)
    add!(B::ITensor, α::Number, A::ITensor)

Add ITensors B and A (or α*A) and store the result in B.
```
B .+= A
B .+= α .* A
```
"""
add!(R::ITensor,T::ITensor) = add!(R,1,T)

function add!(R::ITensor{N},α::Number,T::ITensor{N}) where {N}
  return apply!(R,T,(r,t)->r+α*t)
end

"""
    add!(A::ITensor, α::Number, β::Number, B::ITensor)

Add ITensors α*A and β*B and store the result in A.
```
A .= α .* A .+ β .* B
```
"""
function add!(R::ITensor{N},αr::Number,αt::Number,T::ITensor{N}) where {N}
  return apply!(R,T,(r,t)->αr*r+αt*t)
end

function apply!(R::ITensor{N},T::ITensor{N},f::Function) where {N}
  perm = getperm(inds(R),inds(T))
  TR,TT = tensor(R),tensor(T)

  # TODO: Include type promotion from α
  TR = convert(promote_type(typeof(TR),typeof(TT)),TR)
  TR = permutedims!!(TR,TT,perm,f)

  setstore!(R,store(TR))
  setinds!(R,inds(TR))
  return R
end

"""
    axpy!(a::Number, v::ITensor, w::ITensor)
```
w .+= a .* v
```
"""
axpy!(a::Number,v::ITensor,w::ITensor) = add!(w,a,v)

# This is not implemented correctly
#"""
#w .= a .* v + b .* w
#"""
#axpby!(a::Number,v::ITensor,b::Number,w::ITensor) = add!(w,b,w,a,v)

"""
    scale!(A::ITensor,x::Number) = rmul!(A,x)

Scale the ITensor A by x in-place. May also be written `rmul!`.
```
A .*= x
```
"""
function scale!(T::ITensor,x::Number)
  TT = tensor(T)
  scale!(TT,x)
  return T
end
rmul!(T::ITensor,fac::Number) = scale!(T,fac)

"""
    mul!(A::ITensor,x::Number,B::ITensor)

Scalar multiplication of ITensor B with x, and store the result in A.
Like `A .= x .* B`, and equivalent to `add!(A, 0, x, B)`.
"""
mul!(R::ITensor,α::Number,T::ITensor) = add!(R,0,α,T)
mul!(R::ITensor,T::ITensor,α::Number) = mul!(R,α,T)

function summary(io::IO,
                   T::ITensor)
  print(io,"ITensor ord=$(order(T))")
  for i = 1:order(T)
    print(io," ",inds(T)[i])
  end
  print(io," \n",typeof(store(T)))
end

# TODO: make a specialized printing from Diag
# that emphasizes the missing elements
function Base.show(io::IO,T::ITensor)
  summary(io,T)
  print(io,"\n")
  if !isNull(T)
    Base.print_array(io,array(T))
  end
end

function Base.show(io::IO,
                   mime::MIME"text/plain",
                   T::ITensor)
  summary(io,T)
end

function Base.similar(T::ITensor,
                      element_type=eltype(T))
  return ITensor(similar(tensor(T),element_type))
end

function multSiteOps(A::ITensor,
                     B::ITensor)
  R = prime(A,"Site")
  R *= B
  return mapprime(R,2,1)
end

function readCpp!(io::IO,T::ITensor;kwargs...)
  T.inds = read(io,IndexSet;kwargs...)
  read(io,12) # ignore scale factor
  storage_type = read(io,Int32) # see StorageType enum above
  if storage_type==0 # Null
    T.store = Dense{Nothing}()
  elseif storage_type==1  # DenseReal
    T.store = read(io,Dense{Float64};kwargs...)
  elseif storage_type==2  # DenseCplx
    T.store = read(io,Dense{ComplexF64};kwargs...)
  elseif storage_type==3  # Combiner
    T.store = CombinerStorage(T.inds[1])
  #elseif storage_type==4  # DiagReal
  #elseif storage_type==5  # DiagCplx
  #elseif storage_type==6  # QDenseReal
  #elseif storage_type==7  # QDenseCplx
  #elseif storage_type==8  # QCombiner
  #elseif storage_type==9  # QDiagReal
  #elseif storage_type==10 # QDiagCplx
  #elseif storage_type==11 # ScalarReal
  #elseif storage_type==12 # ScalarCplx
  else
    throw(ErrorException("C++ ITensor storage type $storage_type not yet supported"))
  end
end

function Base.read(io::IO,::Type{ITensor};kwargs...)
  format = get(kwargs,:format,"hdf5")
  T = ITensor()
  if format=="cpp"
    readCpp!(io,T;kwargs...)
  else
    throw(ArgumentError("read ITensor: format=$format not supported"))
  end
  return T
end


text:
module Decoder

using Base.Iterators: take

using ...Client

include("field.jl")
include("process.jl")
include("ticktype.jl")

# Make a shortcut
const pop = popfirst!

function decode(msg, w, ver)

  it = Iterators.Stateful(Field.(msg))

  # The first field is the message ID
  id::Int = pop(it)

  # The second field (version) is ignored for id < 75 and != 3, 5, 11, 17, 21
  if id  < 75 && id ∉ (3, 5, 11, 17, 21)             ||
     id == 21 && ver < Client.PRICE_BASED_VOLATILITY
    pop(it)
  end

  f = get(process, id, nothing)

  if isnothing(f)
    @error "decoder: unknown message" id

  else
    try
      f(it, w, ver)
    catch e
      if e isa EOFError
        @error "decoder: reached end of message" id
      else
        rethrow()
      end
    end

    isempty(it) || @error "decoder: message not fully parsed" id ignored=collect(String, it)
  end
end

end

ext: jl


text:
export rule

@rule GCV(:z, Marginalisation) (q_y_x::Any, q_κ::Any, q_ω::Any, meta::GCVMetadata) = begin

    m, v = mean(q_y_x), cov(q_y_x)
    psi = (m[1] - m[2]) ^ 2 + v[1,1]+ v[2,2] -v[1,2] -v[2,1]
    A = exp(-mean(q_ω)+var(q_ω) / 2)

    a = mean(q_κ)
    b = psi * A
    c = -mean(q_κ)
    d = var(q_κ)

    return ExponentialLinearQuadratic(get_approximation(meta), a, b, c, d)
end

@rule GCV(:z, Marginalisation) (q_y::Any,q_x::Any, q_κ::Any, q_ω::Any, meta::GCVMetadata) = begin

    my, vy = mean(q_y), cov(q_y)
    mx, vx = mean(q_x), cov(q_x)
    psi = (my - mx) ^ 2 + vy+ vx
    A = exp(-mean(q_ω)+var(q_ω) / 2)

    a = mean(q_κ)
    b = psi * A
    c = -mean(q_κ)
    d = var(q_κ)

    return ExponentialLinearQuadratic(get_approximation(meta), a, b, c, d)
end


text:
using ConnectFourSolver
using Test

@testset "ConnectFourSolver.jl" begin
    solver = Solver()
    book = get_book()
    load_book(solver, book)

    open(joinpath(@__DIR__, "Test_L1_R1")) do f
        for line in eachline(f)
            p = Position()
            moves, score = split(line)
            for x in moves
                p(parse(Int, x) - 1)  # starts from 0
            end
            @test solver(p) == parse(Int, score)
        end
    end
end


text:
function mandelbrot(x, y)
    z = c = x + y*im
    for i in 1:100.0; abs(z) > 2 && return i; z = z^2 + c; end; 0
end

x = Node(range(-2, 1, length = 400))
y = Node(range(-1, 1, length = 300))

fig, ax, img = heatmap(x, y, mandelbrot, colormap = Reverse(:deep),
    figure = (resolution = (400, 300),))
hidedecorations!(ax)

record(fig, "mandelbrot.mp4", 1:200) do frame
    x.val = x[] .+ ((-0.562 .- x[]) .* 0.05)
    y[] = y[] .+ ((0.645 .- y[]) .* 0.05)
    autolimits!(ax)
end


text:
#######################################
#
# Plot spectrum of signal
#
#######################################
"""
    plot_spectrum(signal::Vector, fs::Real; kwargs...)

Plot the spectrum of a signal.
"""
function plot_spectrum(
    signal::Vector,
    fs::Real;
    titletext::S = "",
    Fmin::Number = 0,
    Fmax::Number = 90,
    targetFreq::F = 0.0,
    dBPlot::Bool = true,
    noise_level::Number = 0,
    signal_level::Number = 0,
) where {S<:AbstractString,F<:AbstractFloat}

    # Determine fft frequencies
    signal_length = length(signal)
    frequencies = range(0, stop = 1, length = Int(signal_length / 2 + 1)) * fs / 2

    # Calculate fft and convert to power
    fftSweep = 2 / signal_length * fft(signal)
    spectrum = abs.(fftSweep[1:div(signal_length, 2)+1])  # Amplitude
    spectrum = spectrum .^ 2

    valid_idx = (frequencies .<= Fmax) .& (frequencies .>= Fmin)
    spectrum = spectrum[valid_idx]
    frequencies = frequencies[valid_idx]

    # Want a log plot?
    if dBPlot
        spectrum = 10 * log10.(spectrum)
        ylabel = "Response Power (dB)"
    else
        ylabel = "Response Power (uV^2)"
    end

    # Plot signal
    Plots.plot(frequencies, spectrum, lab = "Spectrum", color = :black)
    xlims!(Fmin, Fmax)
    xlabel!("Frequency (Hz)")
    ylabel!(ylabel)
    p = title!(titletext)

    # Plot the noise level if requested
    if noise_level != 0
        if dBPlot
            noise_level = 10 * log10.(noise_level)
        end
        p = Plots.plot!(
            [Fmin, targetFreq + 2],
            [noise_level, noise_level],
            lab = "Noise",
            color = :red,
        )
    end

    # Plot the signal level if requested
    if signal_level != 0
        if dBPlot
            signal_level = 10 * log10.(signal_level)
        end
        Plots.plot!(
            [Fmin, targetFreq],
            [signal_level, signal_level],
            lab = "Signal",
            color = :green,
        )

        targetFreqIdx = something(
            findfirst(
                isequal(minimum(abs.(frequencies .- targetFreq))),
                abs.(frequencies .- targetFreq),
            ),
            0,
        )
        targetFreq = frequencies[targetFreqIdx]
        targetResults = spectrum[targetFreqIdx]
        #TODO remove label for circle rather than make it empty
        p = Plots.plot!(
            [targetFreq],
            [targetResults],
            marker = (:circle, 5, 0.1, :green),
            markerstrokecolor = :green,
            lab = "",
        )
    end
    return p
end

#######################################
#
# Filter response
#
#######################################

# Plot filter response
function plot_filter_response(
    zpk_filter::FilterCoefficients,
    fs::Integer;
    lower::Number = 1,
    upper::Number = 30,
    sample_points::Int = 1024,
)

    frequencies = range(lower, stop = upper, length = 1024)
    h = freqresp(zpk_filter, frequencies * ((2pi) / fs))
    magnitude_dB = 20 * log10.(convert(Array{Float64}, abs.(h)))
    phase_response = (360 / (2 * pi)) * unwrap(convert(Array{Float64}, angle.(h)))

    p1 = Plots.plot(frequencies, magnitude_dB, lab = "")
    p2 = Plots.plot(frequencies, phase_response, lab = "")

    p = Plots.plot(
        p1,
        p2,
        ylabel = ["Magnitude (dB)" "Phase (degrees)"],
        xlabel = "Frequency (Hz)",
        layout = @layout([a; b]),
    )
end


text:
const FluxMPI_initialized = Ref(false)

"""
    Initialized()

Has FluxMPI been initialized?
"""
Initialized() = FluxMPI_initialized[]

"""
    Init(; gpu_devices::Union{Nothing,Vector{Int}} = nothing, verbose::Bool = false)

Setup `FluxMPI`. If GPUs are available and CUDA is functional, each rank is allocated a GPU in a
round-robin fashion.

If calling this function, no need to call `MPI.Init` first.
"""
function Init(; gpu_devices::Union{Nothing,Vector{Int}}=nothing, verbose::Bool=false)
    if Initialized()
        verbose && clean_println("FluxMPI already initialized; Skipping...")
        return
    end

    !MPI.Initialized() && MPI.Init()

    rank = local_rank()

    if CUDA.functional()
        gpu_device = if gpu_devices === nothing
            device_count = length(CUDA.devices())
            (rank + 1) % device_count
        else
            gpu_devices[rank + 1]
        end
        verbose && clean_println("Using GPU $gpu_device")
        CUDA.device!(gpu_device)
    else
        verbose && clean_println("Using CPU")
    end
    FluxMPI_initialized[] = true
    return
end

"""
    local_rank()

Get the rank of the process.
"""
@inline local_rank() = MPI.Comm_rank(MPI.COMM_WORLD)

"""
    total_workers()

Get the total number of workers.
"""
@inline total_workers() = MPI.Comm_size(MPI.COMM_WORLD)

"""
    clean_println(args...; kwargs...)

Add `rank` and `size` information to the printed statement
"""
function clean_println(args...; kwargs...)
    if !Initialized()
        println(args...; kwargs...)
        return
    end
    rank = local_rank()
    size = total_workers()
    if size == 1
        println(args...; kwargs...)
        return
    end
    for r in 0:(size - 1)
        r == rank && println("$(now()) [$(rank) / $(size)] ", args...; kwargs...)
        MPI.Barrier(MPI.COMM_WORLD)
    end
    return
end

@nograd clean_println

"""
    clean_print(args...; kwargs...)

Add `rank` and `size` information to the printed statement
"""
function clean_print(args...; kwargs...)
    if !Initialized()
        print(args...; kwargs...)
        return
    end
    rank = local_rank()
    size = total_workers()
    if size == 1
        print(args...; kwargs...)
        return
    end
    for r in 0:(size - 1)
        r == rank && print("$(now()) [$(rank) / $(size)] ", args...; kwargs...)
        MPI.Barrier(MPI.COMM_WORLD)
    end
    return
end

@nograd clean_print


text:
include("../utils.jl")
include("group_utils.jl")

function _to_affine(p, is_vector = false)
    d = length(p)
    N = Int((sqrt(4d + 1) - 1) / 2)
    P = Matrix{eltype(p)}(undef, N + 1, N + 1)
    setindex!(P, p.parts[1], 1:N, N + 1)
    setindex!(P, p.parts[2], 1:N, 1:N)
    fill!(view(P, N + 1, 1:N), 0)
    if is_vector
        P[N+1, N+1] = 0
    else
        P[N+1, N+1] = 1
    end
    return P
end

@testset "Special Euclidean group" begin
    G = SpecialEuclidean(3)
    @test repr(G) == "SpecialEuclidean(3)"
    M = base_manifold(G)
    @test M === TranslationGroup(3) × SpecialOrthogonal(3)
    Rn = Rotations(3)
    x = Matrix(I, 3, 3)

    t = Vector{Float64}.([1:3, 2:4, 4:6])
    ω = [[1.0, 2.0, 3.0], [3.0, 2.0, 1.0], [1.0, 3.0, 2.0]]
    tuple_pts = [(ti, exp(Rn, x, hat(Rn, x, ωi))) for (ti, ωi) in zip(t, ω)]
    tuple_v = ([-1.0, 2.0, 1.0], hat(Rn, x, [1.0, 0.5, -0.5]))

    @testset "product point" begin
        reshapers = (Manifolds.ArrayReshaper(), Manifolds.StaticReshaper())
        for reshaper in reshapers
            shape_se = Manifolds.ShapeSpecification(reshaper, M.manifolds...)
            pts = [Manifolds.prod_point(shape_se, tp...) for tp in tuple_pts]
            v_pts = [Manifolds.prod_point(shape_se, tuple_v...)]

            g1, g2 = pts[1:2]
            t1, R1 = g1.parts
            t2, R2 = g2.parts
            g1g2 = Manifolds.prod_point(shape_se, R1 * t2 + t1, R1 * R2)
            @test compose(G, g1, g2) ≈ g1g2
            @test _to_affine(g1g2) ≈ _to_affine(g1) * _to_affine(g2)

            w = translate_diff(G, pts[1], Identity(G), v_pts[1])
            w2 = allocate(w)
            w2.parts[1] .= w.parts[1]
            w2.parts[2] .= pts[1].parts[2] * w.parts[2]
            @test _to_affine(w2, true) ≈ _to_affine(pts[1]) * _to_affine(v_pts[1], true)

            test_group(G, pts, v_pts; test_diff = true, diff_convs = [(), (LeftAction(),)])
        end
    end

    @testset "affine matrix" begin
        shape_se = Manifolds.ShapeSpecification(Manifolds.ArrayReshaper(), M.manifolds...)
        pts = [_to_affine(Manifolds.prod_point(shape_se, tp...)) for tp in tuple_pts]
        v_pts = [_to_affine(Manifolds.prod_point(shape_se, tuple_v...), true)]
        test_group(G, pts, v_pts; test_diff = true, diff_convs = [(), (LeftAction(),)])
    end

    @testset "hat/vee" begin
        shape_se = Manifolds.ShapeSpecification(Manifolds.ArrayReshaper(), M.manifolds...)
        x = Manifolds.prod_point(shape_se, tuple_pts[1]...)
        V = Manifolds.prod_point(shape_se, tuple_v...)
        vexp = [V.parts[1]; vee(Rn, x.parts[2], V.parts[2])]
        v = vee(G, x, V)
        @test v ≈ vexp
        @test hat(G, x, v) ≈ V

        v = vee(G, _to_affine(x), _to_affine(V, true))
        @test v ≈ vexp
        @test hat(G, _to_affine(x), v) ≈ _to_affine(V, true)
    end
end


text:
#=

    This is a port of the SPPL model trueskill-poisson-binomial.pynb

    This model:
    * Prior()
    Summary Statistics
  parameters      mean       std   naive_se      mcse          ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64      Float64   Float64       Float64

      perfB1   79.3452    0.9921     0.0099    0.0102    9771.5110    0.9999   105070.0105
      skillA   99.9423    9.4506     0.0945    0.0981    9873.7861    0.9999   106169.7428
      perfA1   89.9122    9.0726     0.0907    0.0948   10013.5210    0.9999   107672.2688
      result    0.8594    0.3476     0.0035    0.0032   10370.4888    1.0000   111510.6318

    * With observation result == 0
    Summary Statistics
    parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
        Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

        perfB1   79.4940    0.8266     0.0083    0.0249   1315.1090    1.0011      171.1045
        skillA   86.2807    4.2628     0.0426    0.1226   1134.3992    0.9999      147.5929
        perfA1   75.9473    3.2281     0.0323    0.0936   1280.7372    1.0003      166.6325
        result    0.0003    0.0173     0.0002    0.0003   3133.8451    1.0002      407.7342

=#

using Turing, StatsPlots, Distributions, StatsBase
using CSV
include("jl_utils.jl")

@model function trueskill_poisson_binomial()
    perfB1 ~ truncated(Binomial(100, 0.9), 0,80)
    skillA ~ truncated(Poisson(100), 77,125)
    perfA1 ~ Binomial(skillA, 0.9)

    if perfA1 > perfB1
       result ~ Dirac(1)
    else
        result ~ Dirac(0)
    end

    # Observation
    true ~ Dirac(result == 0)
end

model = trueskill_poisson_binomial()

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5),  10_000)
# chns = sample(model, PG(5),  MCMCThreads(), 1_000, 4)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, HMC(0.1,6), 1_000)
# chns = sample(model, NUTS(), 1_000)

display(chns)
# display(plot(chns))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Test
using DFTK
include("testcases.jl")

if !Sys.iswindows() && mpi_nprocs() == 1
@testset "Test run_wannier90" begin
    using wannier90_jll

    Si = ElementPsp(silicon.atnum, psp=load_psp("hgh/lda/Si-q4"))
    model  = model_LDA(silicon.lattice, [Si => silicon.positions])
    basis  = PlaneWaveBasis(model; Ecut=5, kgrid=[4, 4, 4])
    scfres = self_consistent_field(basis, tol=1e-12, n_bands=12)

    fileprefix = "wannier90_outputs/Si"
    run_wannier90(scfres; fileprefix,
                  n_wannier=8, bands_plot=true,
                  num_print_cycles=50, num_iter=500,
                  dis_win_max=17.185257,
                  dis_froz_max=6.8318033,
                  dis_num_iter=120,
                  dis_mix_ratio=1.0,
                  wannier_plot=true)

    @test  isfile("wannier90_outputs/Si.amn")
    @test  isfile("wannier90_outputs/Si.chk")
    @test  isfile("wannier90_outputs/Si.eig")
    @test  isfile("wannier90_outputs/Si.mmn")
    @test  isfile("wannier90_outputs/Si.nnkp")
    @test  isfile("wannier90_outputs/Si.win")
    @test  isfile("wannier90_outputs/Si.wout")
    @test !isfile("wannier90_outputs/Si.werr")

    # remove produced files
    rm("wannier90_outputs", recursive=true)
end
end


text:
using MarketRiskPrototype
using Test
using Curves: Curve, Tenor, @t_str
using CurrencyAmounts: @currencies, CurrencyAmount, Currency
using Dates

@testset "MarketRiskPrototype.jl" begin

    @currencies EUR, USD, GBP

    dates = [Date("2020-06-01"), Date("2020-06-02"), Date("2020-06-03"), Date("2020-06-04")]

    fx_curve_1 = Curve([t"2D", t"1W", t"1M", t"3M", t"6M", t"1Y"], [5.5, 5.7, 5.8, 6.1, 6.5, 7.0].*USD./EUR)
    fx_curve_2 = Curve([t"2D", t"1M", t"2M", t"3M", t"6M", t"1Y"], [5.2, 5.4, 5.5, 6.3, 6.9, 7.4].*USD./EUR)

    usd_eur_fx_time_series = FXForwardTimeSeries(EUR, USD, dates, [fx_curve_1, fx_curve_2, fx_curve_1, fx_curve_2])
    @test usd_eur_fx_time_series isa MarketRiskPrototype.AbstractTimeSeries

    disc_curve_1 = Curve([t"2D", t"1M", t"2M", t"3M", t"6M", t"1Y"], [0.99, 0.98, 0.96, 0.95, 0.93, 0.9])
    disc_curve_2 = Curve([t"2D", t"1M", t"2M", t"3M", t"6M", t"1Y"], [0.98, 0.97, 0.955, 0.945, 0.925, 0.91])

    usd_disc_time_series = DiscountTimeSeries(USD, dates, [disc_curve_1, disc_curve_2, disc_curve_1, disc_curve_2])
    @test usd_disc_time_series isa MarketRiskPrototype.AbstractTimeSeries

    fx_curve_1b = Curve([t"2D", t"1W", t"1M", t"3M", t"6M", t"1Y"], 0.5.*[5.5, 5.7, 5.8, 6.1, 6.5, 7.0].*EUR./GBP)
    fx_curve_2b = Curve([t"2D", t"1M", t"2M", t"3M", t"6M", t"1Y"], 0.5.*[5.2, 5.4, 5.5, 6.3, 6.9, 7.4].*EUR./GBP)

    eur_gbp_fx_time_series = FXForwardTimeSeries(GBP, EUR, dates, [fx_curve_1b, fx_curve_2b, fx_curve_1b, fx_curve_2b])

    disc_curve_1b = Curve([t"2D", t"1M", t"2M", t"3M", t"6M", t"1Y"], 0.95.*[0.99, 0.98, 0.96, 0.95, 0.93, 0.9])
    disc_curve_2b = Curve([t"2D", t"1M", t"2M", t"3M", t"6M", t"1Y"], 0.95.*[0.98, 0.97, 0.955, 0.945, 0.925, 0.91])

    eur_disc_time_series = DiscountTimeSeries(EUR, dates, [disc_curve_1b, disc_curve_2b, disc_curve_1b, disc_curve_2b])

    usd_eur_specific_market_data = MarketRiskPrototype.SpecificFXForwardMarketData(EUR, USD, usd_eur_fx_time_series, usd_disc_time_series)
    @test usd_eur_specific_market_data isa MarketRiskPrototype.AbstractProductSpecificMarketData

    market_data_container = FXForwardMarketDataContainer([usd_eur_fx_time_series, eur_gbp_fx_time_series], [usd_disc_time_series, eur_disc_time_series])
    @test market_data_container isa MarketRiskPrototype.AbstractMarketDataContainer
    @test length(market_data_container.fx_forward_curves) == 2 && length(market_data_container.discount_curves) == 2

    fx_forward_1 = FXForward(100000.0EUR, USD, 2.0USD/EUR, t"3M")
    @test fx_forward_1 isa MarketRiskPrototype.AbstractProduct

    fx_forward_2 = FXForward(100000.0GBP, EUR, 1.5EUR/GBP, t"2M")
    fx_forward_3 = FXForward(100000.0USD, EUR, 0.7EUR/USD, t"3M")

    market_data_recoverd = get_specific_market_data(fx_forward_1, market_data_container)
    @test usd_eur_specific_market_data == market_data_recoverd

    pv = price(fx_forward_1, market_data_container, dates[1])
    @test pv < 0USD

    pv2 = price(fx_forward_2, market_data_container, dates[2])
    @test pv2 < 0EUR

    pv3 = price(fx_forward_3, market_data_container, dates[1])
    @test pv != 0EUR

    all_trades = [fx_forward_1, fx_forward_2, fx_forward_3]

    # to calculate pv for multiple trades, just broadcast:
    pv_all_trades = price.(all_trades, market_data_container, dates[1])
    @test length(pv_all_trades) == 3

    pvs = price_time_series(fx_forward_1, market_data_container)
    @test length(pvs) == 4

    # to calculate time series over multiple trades, just broadcast:
    pvs_all_trades = price_time_series.(all_trades, market_data_container)
    @test length(pvs_all_trades) == 3

    var = value_at_risk(pvs .- pv, 0.9)
    @test var isa CurrencyAmount

    shifted_returns = create_historical_scenarios(market_data_container; time_horizon=1, lookback=2)
    @test shifted_returns isa FXForwardMarketDataContainer

    pvs_shifted = price_time_series(fx_forward_1, shifted_returns)
    @test length(pvs_shifted) == 3
    @test pvs_shifted[1] isa CurrencyAmount

    pls = calculate_scenario_pls(fx_forward_1, shifted_returns)
    @test pls isa Vector{CurrencyAmount{Float64, Currency{:USD}}}
    @test length(pls) == 2

    shifted_returns2 = create_historical_scenarios(market_data_container; time_horizon=1, lookback=1, valuation_date=Date("2020-06-03"))
    pvs_shifted2 = price_time_series(fx_forward_1, shifted_returns2)
    @test length(pvs_shifted2) == 2

    var2 = value_at_risk.(all_trades, market_data_container, Date("2020-06-03"); time_horizon=1, lookback=1)
    @test var2 isa Vector{<: CurrencyAmount}
    @test length(var2) == 3

    spot_rates = get_spot_rates(market_data_container)
    var3 = value_at_risk(all_trades, market_data_container, spot_rates, Date("2020-06-03"); time_horizon=1, lookback=1)

end

max_issues_repo_name: lungben/MarketRiskPrototype

max_forks_repo_name: lungben/MarketRiskPrototype


text:
# HACK: work around Pkg.jl#2500
test_project = Base.active_project()
preferences_file = joinpath(dirname(@__DIR__), "LocalPreferences.toml")
test_preferences_file = joinpath(dirname(test_project), "LocalPreferences.toml")
if isfile(preferences_file) && !isfile(test_preferences_file)
    cp(preferences_file, test_preferences_file)
end

using LULESH
using Test
using MPI

if parse(Bool, get(ENV, "ENZYME_CI", "false"))
    using Enzyme_jll
    @info "Testing against" Enzyme_jll.libEnzyme
    args = `--enzyme`
else
    args = ``
end

using JET
MPI.Init()

@testset "JET" begin
    domain = Domain(LuleshProblem(1, true, 45, 1, 1, 1, Vector, Float64, MPI.COMM_WORLD))
    test_call(lagrangeLeapFrog, (typeof(domain),))
    test_opt(lagrangeLeapFrog, (typeof(domain),))
end

MPI.install_mpiexecjl(; destdir = ".", force=true)

function run_example(name, args; nranks=nothing)
    example = joinpath(@__DIR__, "..", "examples", name)

    @testset "$(basename(example)) $args $nranks" begin
        cmd = `$(Base.julia_cmd()) --startup-file=no $(example) $(args)`
        if nranks !== nothing
            cmd = `./mpiexecjl --project=$(test_project) -n $(nranks) $(cmd)`
        end
        @debug "Testing $example $nranks" cmd
        @test success(pipeline(cmd, stderr=stderr))
    end
end

@testset "examples" begin
    @testset "benchmark.jl" begin
        run_example("benchmark.jl", `-s 8 --mpi $(args)`; nranks=1)
    end
end


text:
# Uses new formulation from limbdark paper.
include("cel_bulirsch.jl")

function s2(r::T,b::T) where {T <: Real}
# For now, just compute linear component:
Lambda1 = zero(typeof(b))
if b >= 1.0+r ||  r == 0.0
  # No occultation:
  Lambda1 = zero(b)  # Case 1
elseif b <= r-1.0
  # Full occultation:
  Lambda1 = zero(b)  # Case 11
else
  if b == 0
    Lambda1 = -2/3*sqrt(1.0-r^2)^3 # Case 10
  elseif b==r
    if r == 0.5
      Lambda1 = 1/3-4/(9pi) # Case 6
    elseif r < 0.5
      m = 4r^2
      Lambda1 = 1/3+2/(9pi)*cel_bulirsch(m,one(r),m-3,(1-m)*(2m-3)) # Case 5
    else
      m = 4r^2; minv = inv(m)
      Lambda1 = 1/3+1/(9pi*r)*cel_bulirsch(minv,one(r),m-3,1-m) # Case 7
    end
  else
    onembpr2 = 1-(b+r)^2; onembmr2=1-(b-r)^2; fourbr = 4b*r
    k2 = onembpr2/fourbr+1
    if (b+r) > 1.0 # k^2 < 1, Case 2, Case 8
      k2c = -onembpr2/fourbr; kc = sqrt(k2c)
      Eofk = cel_bulirsch(k2,kc,one(b),one(b),k2c) # Complete elliptic integral of second kind
      Lambda1 = onembmr2*(k2c*(b-r)*(b+r)*cel_bulirsch(k2,kc,(b-r)^2*k2c,zero(b),3one(b))
          -(3-6*r^2-2*b*r)*cel_bulirsch(k2,kc,one(b),one(b),zero(b))-fourbr*Eofk)/(9*pi*sqrt(b*r))
    elseif (b+r) < 1.0  # k^2 > 1, Case 3, Case 9
      k2inv = inv(k2); k2c =onembpr2/onembmr2; kc = sqrt(k2c)
      Eofk = cel_bulirsch(k2inv,kc,one(b),one(b),k2c) # Complete elliptic integral of second kind
      bmrdbpr = (b-r)/(b+r);
      mu = 3bmrdbpr/onembmr2
      p = bmrdbpr^2*onembpr2/onembmr2
      Lambda1 = 2*sqrt(onembmr2)*(onembpr2*cel_bulirsch(k2inv,kc,p,1.0+mu,p+mu)
             -(4-7r^2-b^2)*Eofk)/(9*pi)
    else
      # b+r = 1 or k^2=1, Case 4 (extending r up to 1)
      Lambda1 = 2/(3pi)*acos(1.-2.*r)-4/(9pi)*(3+2r-8r^2)*sqrt(r*b)-2/3*convert(typeof(b),r>.5)
    end
  end
end
flux = 1.0-1.5*Lambda1-convert(typeof(b),r>b)
return flux
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
#=
    unconditional.jl : General representation of sfuncs that don't depend on anything.
=#

export Unconditional

"""
    abstract type Unconditional{T} <: SFunc{Tuple{}, T}

`Unconditional` is a general representation of an *sfunc* that does not depend
on anything.  It has no input.

# Type parameters
- `T`: the input type(s) of the `Unconditional`
- `P`: the parameter type(s) of the `Unconditional`
"""
abstract type Unconditional{T} <: SFunc{Tuple{}, T} end

@impl begin
    struct UnconditionalSample end
    function sample(sf::Unconditional{T}, i::Tuple{})::T
        sample(sf, x) # should never be called
    end
end

@impl begin
    struct UnconditionalLogcpdf end
    function logcpdf(sf::Unconditional{T}, i::Tuple{}, o::T)
        logcpdf(sf, i, o) # should never be called
    end
end

@impl begin
    struct UnconditionalMakeFactors
        numpartitions::Int64 = 100
    end

    function make_factors(sf::Unconditional{T},
            range::Vector{T},
            parranges::NTuple{N,Vector},
            id,
            parids::Tuple)::Tuple{Vector{<:Scruff.Utils.Factor}, Vector{<:Scruff.Utils.Factor}} where {T,N}

        (lowers, uppers) = bounded_probs(sf, range, ())
        keys = (id,)
        dims = (length(range),)
        return ([Factor(dims, keys, lowers)], [Factor(dims, keys, uppers)])
    end
end

# compute_pi(u::Unconditional, range, parranges, incoming_pis) = compute_pi(u, range)

# unconditional sfuncs do not have parents, so send_lambda trivially sends an empty message
@impl begin
    struct UnconditionalSendLambda end
    function send_lambda(sf::Unconditional{T},
                    lambda::Score{T},
                    range::Vector{T},
                    parranges::NTuple{N,Vector},
                    incoming_pis::Tuple,
                    parent_idx::Integer)::Score where {T,N}
        SoftScore(Float64[], Float64[])
    end
end


text:
using PatModules
@auxmodule utils () begin

    function do_nothing()
    end
end


text:
# Code below is hard-copied from
# https://github.com/ahwillia/NonNegLeastSquares.jl
# See that repository for the license of this file.
# (Temporary stop-gap until the above package is registered)

using Distributed

function nonneg_lsq(A, B; alg::Symbol = :pivot, gram::Bool = false, kwargs...)

    if alg == :fnnls
        return fnnls(A, B; gram = gram, kwargs...)
    else
        error("Specified algorithm :", alg, " not recognized.")
    end
end

# If second input is a vector, convert it to a matrix
nonneg_lsq(A, b::AbstractVector; kwargs...) = nonneg_lsq(A, b[:, :]; kwargs...)

"""
x = fnnls(AtA, Atb; ...)
Returns x that solves A*x = b in the least-squares sense, subject to x >=0. The
inputs are the cross-products AtA = A'*A and Atb = A'*b. Uses the modified
active set method of Bro and De Jong (1997).
Optional arguments:
    tol: tolerance for nonnegativity constraints
    max_iter: maximum number of iterations (counts inner loop iterations)
References:
    Bro R, De Jong S. A fast non-negativitity-constrained least squares
    algorithm. Journal of Chemometrics. 11, 393–401 (1997)
"""
function fnnls(
    AtA,
    Atb::AbstractVector{T};
    tol::Float64 = 1e-8,
    max_iter = 30 * size(AtA, 2),
) where {T}

    n = size(AtA, 1)
    x = zeros(T, n)
    s = zeros(T, n)

    # P is a bool array storing positive elements of x
    # i.e., x[P] > 0 and x[~P] == 0
    P = x .> tol
    w = Atb - AtA * x

    # We have reached an optimum when either:
    #   (a) all elements of x are positive (no nonneg constraints activated)
    #   (b) ∂f/∂x = A' * (b - A*x) > 0 for all nonpositive elements of x
    iter = 0
    while sum(P) < n && any(w[(!).(P)] .> tol) && iter < max_iter

        # find i that maximizes w, restricting i to indices not in P
        # Note: the while loop condition guarantees at least one w[~P]>0
        i = argmax(w .* (!).(P))

        # Move i to P
        P[i] = true

        # Solve least-squares problem, with zeros for columns/elements not in P
        s[P] = AtA[P, P] \ Atb[P]
        s[(!).(P)] .= zero(eltype(s)) # zero out elements not in P

        # Inner loop: deal with negative elements of s
        while any(s[P] .<= tol)
            iter += 1

            # find indices in P where s is negative
            ind = @__dot__ (s <= tol) & P

            # calculate step size, α, to prevent any xᵢ from going negative
            α = minimum(x[ind] ./ (x[ind] - s[ind]))

            # update solution (pushes some xᵢ to zero)
            x += α * (s - x)

            # Remove all i in P where x[i] == 0
            for i = 1:n
                if P[i] && abs(x[i]) < tol
                    P[i] = false # remove i from P
                end
            end

            # Solve least-squares problem again, zeroing nonpositive columns
            s[P] = AtA[P, P] \ Atb[P]
            s[(!).(P)] .= zero(eltype(s)) # zero out elements not in P
        end

        # update solution
        x = deepcopy(s)
        w .= Atb - AtA * x
    end
    return x
end

function fnnls(
    A,
    B::AbstractMatrix;
    gram::Bool = false,
    use_parallel::Bool = false,
    kwargs...,
)

    n = size(A, 2)
    k = size(B, 2)

    if gram
        # A,B are actually Gram matrices
        AtA = A
        AtB = B
    else
        # cache matrix computations
        AtA = A' * A
        AtB = A' * B
    end

    if use_parallel && nprocs() > 1
        X = @distributed (hcat) for i = 1:k
            fnnls(AtA, AtB[:, i]; kwargs...)
        end
    else
        X = Array{eltype(B)}(undef, n, k)
        for i = 1:k
            X[:, i] = fnnls(AtA, AtB[:, i]; kwargs...)
        end
    end

    return X
end


text:
using IterativeSolvers
using Base.Test
using LinearMaps

srand(1234321)

@testset "LSQR" begin

@testset "Small dense matrix" for T = (Float32, Float64)
    A = rand(T, 10, 5)
    b = rand(T, 10)
    x, history = lsqr(A, b, log = true)
    @test isa(history, ConvergenceHistory)
    @test norm(x - A\b) ≤ √eps(T)
    @test history.isconverged
    @test last(history[:resnorm]) ≈ norm(b - A * x) atol=√eps(T)
end

function sol_matrix(m, n)
    mn = min(m, n)
    spdiagm((1.0 : mn - 1, 1.0 : mn), (-1, 0), m, n)
end

@testset "SOL test" for (m, n) = ((10, 10), (20, 10), (20, 10))
    # Test adapted from the BSD-licensed Matlab implementation at
    #    http://www.stanford.edu/group/SOL/software/lsqr.html
    #              Michael Saunders, Systems Optimization Laboratory,
    #              Dept of MS&E, Stanford University.
    #-----------------------------------------------------------------------
    A = LinearMap(sol_matrix(m, n))
    x = float.(n : -1 : 1)
    b = float(A*x)
    x_lsqr = lsqr(A, b, atol=1e-6, btol=1e-6, conlim=1e10, maxiter=10n)
    @test norm(b - A * x_lsqr) ≤ 1e-4
end
end


text:
using DashBootstrapComponents

list_group = dbc_listgroup([
    dbc_listgroupitem("Item 1"),
    dbc_listgroupitem("Item 2"),
    dbc_listgroupitem("Item 3"),
]);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
"""
    inner_update_cd!(i, j, B, resid, X, Z, norms, lambda, reg)

Updates a single coefficient estimate in place (to be called by
`update_cd_cyclic!`, `update_cd_random!`, `update_cd_active_cyclic!`, or
`update_cd_active_random!`) when `X` and `Z` are both standardized.

# Arguments

- i = row index of the coefficient to update
- j = column index of the coefficient to update
- B = 2d array of floats consisting of coefficient estimates
- resid = 2d array of floats consisting of the residuals
- X = 2d array of floats consisting of the row covariates, with all
  categorical variables coded in appropriate contrasts
- Z = 2d array of floats consisting of the column covariates, with all
  categorical variables coded in appropriate contrasts
- norms = `nothing`
- lambda = lambda penalty, a floating scalar
- reg = 2d array of bits, indicating whether or not to regularize each of the
  coefficients

# Value

None; updates a coefficient in place

"""
function inner_update_cd!(i::Int64, j::Int64, B::AbstractArray{Float64,2},
                          resid::AbstractArray{Float64,2},
                          X::AbstractArray{Float64,2},
                          Z::AbstractArray{Float64,2},
                          norms::Nothing, lambda::Float64, reg::BitArray{2})

    gradient = calc_grad(X[:,i], Z[:,j], resid) # Calculate gradient
    b = B[i,j] # Store the current coefficient.

    b2update = b - gradient # L2 update
    b2sign = sign(b2update) # L2 update sign

    # Update the current coefficient
    b = reg[i,j] ? prox(b, gradient, b2sign, lambda, norms) : b2update

    if b != B[i,j]
        ger!(B[i,j]-b, X[:,i], Z[:,j], resid) # Update residuals
        B[i,j] = b # Update coefficient estimates in the current copy
    end
end

"""
    inner_update_cd!(i, j, B, resid, X, Z, norms, lambda, reg)

Updates a single coefficient estimate in place (to be called by
`update_cd_cyclic!`, `update_cd_random!`, `update_cd_active_cyclic!`, or
`update_cd_active_random!`) when `X` and `Z` are not standardized.

# Arguments

- i = row index of the coefficient to update
- j = column index of the coefficient to update
- B = 2d array of floats consisting of coefficient estimates
- resid = 2d array of floats consisting of the residuals
- X = 2d array of floats consisting of the row covariates, with all
  categorical variables coded in appropriate contrasts
- Z = 2d array of floats consisting of the column covariates, with all
  categorical variables coded in appropriate contrasts
- norms = 2d array of floats consisting of the norms corresponding to each
  coefficient
- lambda = lambda penalty, a floating scalar
- reg = 2d array of bits, indicating whether or not to regularize each of the
  coefficients

# Value

None; updates a coefficient in place

"""
function inner_update_cd!(i::Int64, j::Int64, B::AbstractArray{Float64,2},
                          resid::AbstractArray{Float64,2},
                          X::AbstractArray{Float64,2},
                          Z::AbstractArray{Float64,2},
                          norms::AbstractArray{Float64,2}, lambda::Float64,
                          reg::BitArray{2})

    gradient = calc_grad(X[:,i], Z[:,j], resid)  # Calculate gradient
    b = B[i,j] # Store the current coefficient.

    b2update = b - gradient/norms[i,j] # L2 update
    b2sign = sign(b2update) # L2 update sign

    # Update the current coefficient
    b = reg[i,j] ? prox(b, gradient, b2sign, lambda, norms[i,j]) : b2update

    if b != B[i,j]
        ger!(B[i,j]-b, X[:,i], Z[:,j], resid) # Update residuals
        B[i,j] = b # Update coefficient estimates in the current copy
    end
end

"""
    update_cd_cyclic!(B, resid, X, Z, norms, lambda, reg)

Cyclically updates coefficients in place for each coordinate descent iteration.

# Arguments

- B = 2d array of floats consisting of coefficient estimates
- resid = 2d array of floats consisting of the residuals
- X = 2d array of floats consisting of the row covariates, with all
  categorical variables coded in appropriate contrasts
- Z = 2d array of floats consisting of the column covariates, with all
  categorical variables coded in appropriate contrasts
- norms = 2d array of floats consisting of the norms corresponding to each
  coefficient or `nothing`
- lambda = lambda penalty, a floating scalar
- reg = 2d array of bits, indicating whether or not to regularize each of the
  coefficients

# Value

None; updates coefficients in place

"""
function update_cd_cyclic!(B::AbstractArray{Float64,2},
                           resid::AbstractArray{Float64,2},
                           X::AbstractArray{Float64,2},
                           Z::AbstractArray{Float64,2},
                           norms, lambda::Float64, reg::BitArray{2})

    # Cycle through all the coefficients
    for j = 1:size(B,2), i = 1:size(B,1)
        inner_update_cd!(i, j, B, resid, X, Z, norms, lambda, reg)
    end
end

"""
    update_cd_random!(B, resid, X, Z, norms, lambda, reg)

Randomly updates coefficients in place for each coordinate descent iteration.

# Arguments

- B = 2d array of floats consisting of coefficient estimates
- resid = 2d array of floats consisting of the residuals
- X = 2d array of floats consisting of the row covariates, with all
  categorical variables coded in appropriate contrasts
- Z = 2d array of floats consisting of the column covariates, with all
  categorical variables coded in appropriate contrasts
- norms = 2d array of floats consisting of the norms corresponding to each
  coefficient or `nothing`
- lambda = lambda penalty, a floating scalar
- reg = 2d array of bits, indicating whether or not to regularize each of the
  coefficients

# Value

None; updates coefficients in place

"""
function update_cd_random!(B::AbstractArray{Float64,2},
                           resid::AbstractArray{Float64,2},
                           X::AbstractArray{Float64,2},
                           Z::AbstractArray{Float64,2},
                           norms, lambda::Float64, reg::BitArray{2})

    # Randomly shuffle the coefficient indices
    rand_idx = Random.shuffle(collect(zip(repeat(1:size(B,1), inner=size(B,2)),
                                      repeat(1:size(B,2), outer=size(B,1)))))

    # Update the coeffcients randomly
    for (i, j) in rand_idx
        inner_update_cd!(i, j, B, resid, X, Z, norms, lambda, reg)
    end
end

"""
    update_cd_active_cyclic!(B, resid, X, Z, norms, lambda, reg, nonreg_idx,
                             active_idx)

Cyclically updates active set of coefficients in place for each coordinate
descent iteration.

# Arguments

- B = 2d array of floats consisting of coefficient estimates
- resid = 2d array of floats consisting of the residuals
- X = 2d array of floats consisting of the row covariates, with all
  categorical variables coded in appropriate contrasts
- Z = 2d array of floats consisting of the column covariates, with all
  categorical variables coded in appropriate contrasts
- norms = 2d array of floats consisting of the norms corresponding to each
  coefficient or `nothing`
- lambda = lambda penalty, a floating scalar
- reg = 2d array of bits, indicating whether or not to regularize each of the
  coefficients
- nonreg_idx = tuple with the 2d indices of the non-regularized coefficients
  as two 1d arrays of integers
- active_idx = tuple with the 2d indices of the active coefficients as two 1d
  arrays of integers

# Value

None; updates coefficients in place

# Some notes

Given that you pass in the indices for the non-regularized and active
(regularized) coefficients separately, this function can be further optimized
so that you don't check for regularization when updating each coefficient
with `inner_update!`.

"""
function update_cd_active_cyclic!(B::AbstractArray{Float64,2},
                                  resid::AbstractArray{Float64,2},
                                  X::AbstractArray{Float64,2},
                                  Z::AbstractArray{Float64,2},
                                  norms, lambda::Float64, reg::BitArray{2},
                                  nonreg_idx::Tuple{AbstractArray{Int64,1},
                                                    AbstractArray{Int64,1}},
                                  active_idx::Tuple{AbstractArray{Int64,1},
                                                    AbstractArray{Int64,1}})

  # Cycle through the non-regularized coefficients
  for (i,j) in zip(nonreg_idx[1], nonreg_idx[2])
    inner_update_cd!(i, j, B, resid, X, Z, norms, lambda, reg)
  end

  # Cycle through the active set
  for (i,j) in zip(active_idx[1], active_idx[2])
    inner_update_cd!(i, j, B, resid, X, Z, norms, lambda, reg)
  end
end

"""
    update_cd_active_random!(B, resid, X, Z, norms, lambda, reg, nonreg_idx,
                             active_idx)

Randomly updates active set of coefficients in place for each coordinate
descent iteration.

# Arguments

- B = 2d array of floats consisting of coefficient estimates
- resid = 2d array of floats consisting of the residuals
- X = 2d array of floats consisting of the row covariates, with all
  categorical variables coded in appropriate contrasts
- Z = 2d array of floats consisting of the column covariates, with all
  categorical variables coded in appropriate contrasts
- norms = 2d array of floats consisting of the norms corresponding to each
  coefficient or `nothing`
- lambda = lambda penalty, a floating scalar
- reg = 2d array of bits, indicating whether or not to regularize each of the
  coefficients
- nonreg_idx = tuple with the 2d indices of the non-regularized coefficients
  as two 1d arrays of integers
- active_idx = tuple with the 2d indices of the active coefficients as two 1d
  arrays of integers

# Value

None; updates coefficients in place

"""
function update_cd_active_random!(B::AbstractArray{Float64,2},
                                  resid::AbstractArray{Float64,2},
                                  X::AbstractArray{Float64,2},
                                  Z::AbstractArray{Float64,2},
                                  norms, lambda::Float64, reg::BitArray{2},
                                  nonreg_idx::Tuple{AbstractArray{Int64,1},
                                                    AbstractArray{Int64,1}},
                                  active_idx::Tuple{AbstractArray{Int64,1},
                                                    AbstractArray{Int64,1}})

    # Randomly shuffle the coefficient indices for non-regularized
    # coefficients and those in the active set
    rand_idx = Random.shuffle(collect(zip(vcat(nonreg_idx[1], active_idx[1]),
                                          vcat(nonreg_idx[2], active_idx[2]))))

    # Update the non-regularized coefficients and active set randomly
    for (i, j) in rand_idx
        inner_update_cd!(i, j, B, resid, X, Z, norms, lambda, reg)
    end
end

"""
    cd!(X, Y, Z, lambda, B, regXidx, regZidx, reg, norms;
        isVerbose, stepsize, isRandom, thresh, maxiter)

Performs coordinate descent using either random or cyclic updates. Does NOT
take advantage of the active set; see `cd_active!`.

# Arguments

- X = 2d array of floats consisting of the row covariates, with all
  categorical variables coded in appropriate contrasts
- Y = 2d array of floats consisting of the multivariate response
- Z = 2d array of floats consisting of the column covariates, with all
  categorical variables coded in appropriate contrasts
- lambda = lambda penalty, a floating scalar
- B = 2d array of floats consisting of starting coefficient estimates
- regXidx = 1d array of indices corresponding to regularized X covariates
- regZidx = 1d array of indices corresponding to regularized Z covariates
- reg = 2d array of bits, indicating whether or not to regularize each of the
  coefficients
- norms = 2d array of floats consisting of the norms corresponding to each
  coefficient or `nothing`

# Keyword arguments

- isVerbose = boolean flag indicating whether or not to print messages.
  Defaults to `true`.
- stepsize = float; step size for updates
  (irrelevant for coordinate descent). Defaults to `0.01`.
- isRandom = boolean flag indicating whether to use random or cyclic updates.
  Defaults to `true`.
- thresh = threshold at which the coefficients are considered to have
  converged, a floating scalar. Defaults to `10^(-7)`.
- maxiter = maximum number of update iterations. Defaults to `10^10`.

# Value

None; updates coefficients in place

# Some notes

Convergence is determined as when the log ratio of the current and previous
criteria is less than the threshold `thresh`.

"""
function cd!(X::AbstractArray{Float64,2}, Y::AbstractArray{Float64,2},
             Z::AbstractArray{Float64,2}, lambda::Float64,
             B::AbstractArray{Float64,2},
             regXidx::AbstractArray{Int64,1},
             regZidx::AbstractArray{Int64,1}, reg::BitArray{2}, norms;
             isVerbose::Bool=true, stepsize::Float64=0.01,
             isRandom::Bool=true, thresh::Float64=10.0^(-7),
             maxiter::Int=10^10)

    # Set the updates to random or cyclic.
    if isRandom == true
        update! = update_cd_random!
        update_active! = update_cd_active_random!
    else
        update! = update_cd_cyclic!
        update_active! = update_cd_active_cyclic!
    end

    # Calculate residuals
    resid = calc_resid(X, Y, Z, B)
    # Denominators of criterion
    crit_denom = [size(X,1)*size(Z,1), size(X,2)*size(Z,2)]
    # Placeholder to store the old criterion
    oldcrit = 1.0
    # Calculate the current criterion
    crit = criterion(B[regXidx, regZidx], resid, lambda, crit_denom)

    iter = 0
    # Iterate until coefficients converge or maximum iterations have been
    # reached.
    while (abs(log(crit/oldcrit)) > thresh) && (iter < maxiter)
        # Store the current criterion
        oldcrit = crit
        # Update the coefficient estimates
        update!(B, resid, X, Z, norms, lambda, reg)
        # Calculate the criterion after updating
        crit = criterion(B[regXidx, regZidx], resid, lambda, crit_denom)

        iter += 1 # Increment the number of iterations
        # Print warning message if coefficient estimates do not converge.
        if iter==(maxiter)
            println("Estimates did not converge.")
        end
    end

    println_verbose(string("Criterion: ", crit), isVerbose)
    println_verbose(string("Number of iterations: ", iter), isVerbose)
end

"""
    cd_active!(X, Y, Z, lambda, B, regXidx, regZidx, reg, norms;
               isVerbose, stepsize, isRandom, thresh, maxiter)

Performs coordinate descent, taking advantage of the active set, using either
random or cyclic updates.

# Arguments

- X = 2d array of floats consisting of the row covariates, with all
  categorical variables coded in appropriate contrasts
- Y = 2d array of floats consisting of the multivariate response
- Z = 2d array of floats consisting of the column covariates, with all
  categorical variables coded in appropriate contrasts
- lambda = lambda penalty, a floating scalar
- B = 2d array of floats consisting of starting coefficient estimates
- regXidx = 1d array of indices corresponding to regularized X covariates
- regZidx = 1d array of indices corresponding to regularized Z covariates
- reg = 2d array of bits, indicating whether or not to regularize each of the
  coefficients
- norms = 2d array of floats consisting of the norms corresponding to each
  coefficient or `nothing`

# Keyword arguments

- isVerbose = boolean flag indicating whether or not to print messages.
  Defaults to `true`.
- stepsize = float; step size for updates
  (irrelevant for coordinate descent). Defaults to `0.01`.
- isRandom = boolean flag indicating whether to use random or cyclic updates.
  Defaults to `true`.
- thresh = threshold at which the coefficients are considered to have
  converged, a floating scalar. Defaults to `10^(-7)`.
- maxiter = maximum number of update iterations. Defaults to `10^10`.

# Value

None; updates coefficients in place

# Some notes

Convergence is determined as when the log ratio of the current and previous
criteria is less than the threshold `thresh`.

"""
function cd_active!(X::AbstractArray{Float64,2}, Y::AbstractArray{Float64,2},
                    Z::AbstractArray{Float64,2}, lambda::Float64,
                    B::AbstractArray{Float64,2},
                    regXidx::AbstractArray{Int64,1},
                    regZidx::AbstractArray{Int64,1}, reg::BitArray{2}, norms;
                    isVerbose::Bool=true, stepsize::Float64=0.01,
                    isRandom::Bool=true, thresh::Float64=10.0^(-7),
                    maxiter::Int=10^10)

    # Set the updates to random or cyclic.
    if isRandom == true
        update! = update_cd_random!
        update_active! = update_cd_active_random!
    else
        update! = update_cd_cyclic!
        update_active! = update_cd_active_cyclic!
    end

    # Calculate residuals
    resid = calc_resid(X, Y, Z, B)
    # Denominators of criterion
    crit_denom = [size(X,1)*size(Z,1), size(X,2)*size(Z,2)]
    # Placeholder to store the old criterion
    oldcrit = 1.0

    # Initial iteration of updates
    iter = 1
    # Update the coefficient estimates
    update!(B, resid, X, Z, norms, lambda, reg)
    # Calculate the criterion after the first update
    crit = criterion(B[regXidx, regZidx], resid, lambda, crit_denom)

    # Indices of non-regularized coefficients
    nonreg_cart_idx = findall(reg.==false)
    nonreg_idx = Tuple([map(i->i[1], nonreg_cart_idx),
                           map(i->i[2], nonreg_cart_idx)])
    # Identify active set
    active_cart_idx = findall((B.!=0.0) .& (reg.==true))
    active_idx = Tuple([map(i->i[1], active_cart_idx),
                           map(i->i[2], active_cart_idx)])

    # Iterate until coefficients converge or maximum iterations have been
    # reached.
    while (abs(log(crit/oldcrit)) > thresh) && (iter < maxiter)
        println_verbose(string("Iterating through active set at iteration ",
                               iter, "."), isVerbose)
        # Iterate over active set
        while (abs(log(crit/oldcrit)) > thresh) && (iter < maxiter)
            # Store the current criterion
            oldcrit = crit
            # Update the non-regularized and active coefficients
            update_active!(B, resid, X, Z, norms, lambda, reg, nonreg_idx,
                           active_idx)
            # Calculate the criterion after updating
            crit = criterion(B[regXidx, regZidx], resid, lambda, crit_denom)

            iter += 1 # Increment the number of iterations
            # Warning message if coefficient estimates do not converge.
            if iter==(maxiter)
                println("Estimates did not converge.")
            end
        end

        # Store the current criterion
        oldcrit = crit
        # Update the coefficient estimates
        update!(B, resid, X, Z, norms, lambda, reg)
        # Calculate the criterion after updating
        crit = criterion(B[regXidx, regZidx], resid, lambda, crit_denom)

        # Re-identify active set
        active_cart_idx = findall((B.!=0.0) .& (reg.==true))
        active_idx = Tuple([map(i->i[1], active_cart_idx),
                           map(i->i[2], active_cart_idx)])
        iter += 1 # Increment the number of iterations
        # Warning message if coefficient estimates do not converge.
        if iter==(maxiter)
            println("Estimates did not converge.")
        end
    end

    println_verbose(string("Criterion: ", crit), isVerbose)
    println_verbose(string("Number of iterations: ", iter), isVerbose)
end


text:
println("including tssm_schroedinger.jl for type ", T)
for (METHOD, SUF, COMPLEX_METHOD, DIM) in (
                 (:Schroedinger1D, :_schroedinger_1d, true, 1 ),
                 (:Schroedinger2D, :_schroedinger_2d, true, 2 ),
                 (:Schroedinger3D, :_schroedinger_3d, true, 3 ),
                 (:SchroedingerReal1D, :_schroedinger_real_1d, false, 1 ),
                 (:SchroedingerReal2D, :_schroedinger_real_2d, false, 2 ),
                 (:SchroedingerReal3D, :_schroedinger_real_3d, false, 3 ),
                 (:SchroedingerRotating2D, :_schroedinger_rotating_2d, true, 2 ),
                 (:SchroedingerRotating3D, :_schroedinger_rotating_3d, true, 3 ),
                 (:SchroedingerRotatingReal2D, :_schroedinger_rotating_real_2d, false, 2 ),
                 (:SchroedingerRotatingReal3D, :_schroedinger_rotating_real_3d, false, 3 ),
                 (:SchroedingerHermite1D, :_schroedinger_hermite_1d, true, 1 ),
                 (:SchroedingerHermite2D, :_schroedinger_hermite_2d, true, 2 ),
                 (:SchroedingerHermite3D, :_schroedinger_hermite_3d, true, 3 ),
                 (:SchroedingerHermiteReal1D, :_schroedinger_hermite_real_1d, false, 1 ),
                 (:SchroedingerHermiteReal2D, :_schroedinger_hermite_real_2d, false, 2 ),
                 (:SchroedingerHermiteReal3D, :_schroedinger_hermite_real_3d, false, 3 ),
                 (:SchroedingerGeneralizedLaguerre2D, :_schroedinger_gen_laguerre_2d, true, 2 ),
                 (:SchroedingerGeneralizedLaguerreHermite3D, :_schroedinger_gen_laguerre_hermite_3d, true, 3 ),
                 (:SchroedingerGeneralizedLaguerreReal2D, :_schroedinger_gen_laguerre_real_2d, false, 2 ),
                 (:SchroedingerGeneralizedLaguerreHermiteReal3D, :_schroedinger_gen_laguerre_hermite_real_3d, false, 3 ),
                )
    println("    ", METHOD)
    if T == :Float128
        SUF = Symbol(SUF, "_wf128")
    end
    WF = Symbol(:Wf,METHOD)
    @eval begin

        function get_hbar(m::($METHOD){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_hbar",SUF))), ($T),
                 (Ptr{Nothing}, ), m.m )
        end

        function get_cubic_coupling(m::($METHOD){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_cubic_coupling",SUF))), ($T),
                 (Ptr{Nothing}, ), m.m )
        end

        function set_cubic_coupling!(m::($METHOD){$T}, cc::Real)
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_cubic_coupling",SUF))), Nothing,
                 (Ptr{Nothing}, ($T)), m.m, cc )
        end

        function load_potential!(m::($METHOD){$T}, filename::AbstractString)
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"load_potential",SUF))), Nothing,
                 (Ptr{Nothing}, Cstring, Int32,), m.m, filename, length(filename))
        end

        function save_potential(m::($METHOD){$T}, filename::AbstractString)
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"save_potential",SUF))), Nothing,
                 (Ptr{Nothing}, Cstring, Int32,), m.m, filename, length(filename))
        end

        function kinetic_energy(psi::($WF){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"kinetic_energy_wf",SUF))), ($T),
                 (Ptr{Nothing}, ), psi.p )
        end

        function potential_energy(psi::($WF){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"potential_energy_wf",SUF))), ($T),
                 (Ptr{Nothing}, ), psi.p )
        end

        function interaction_energy(psi::($WF){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"interaction_energy_wf",SUF))), ($T),
                 (Ptr{Nothing}, ), psi.p )
        end

        function get_energy_expectation_deviation(psi::($WF){$T})
           ans = zeros(($T), 2)
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_energy_expectation_deviation_wf",SUF))), Nothing,
                (Ptr{Nothing}, Ptr{$T}), psi.p, ans )
           tuple(ans...)
        end

        function get_realspace_observables(psi::($WF){$T})
           ans = zeros(($T), 2+2*($DIM))
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_realspace_observables_wf",SUF))), Nothing,
                (Ptr{Nothing}, Ptr{$T}), psi.p, ans )
           tuple(ans...)
        end

    end #eval

    if COMPLEX_METHOD
        @eval begin
            function propagate_B!(psi::($WF){$T}, dt::Number)
                ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"propagate_B_wf",SUF))), Nothing,
                        (Ptr{Nothing}, Complex{$T},), psi.p, dt)
            end

            function propagate_B_derivative!(this::($WF){$T}, other::($WF){$T}, dt::Number)
               if this.m ≠ other.m
                   error("this and other must belong to the same method")
               end
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"propagate_B_derivative_wf",SUF))), Nothing,
                                (Ptr{Nothing}, Ptr{Nothing}, Complex{$T}),
                                 this.p, other.p, dt)
            end

            function imaginary_time_propagate_A!(psi::($WF){$T}, dt::Number)
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"imaginary_time_propagate_A_wf",SUF))), Nothing,
                                (Ptr{Nothing}, Complex{$T},), psi.p, dt)
            end

            function imaginary_time_propagate_B!(psi::($WF){$T}, dt::Number, method_for_B::Integer=0)
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"imaginary_time_propagate_B_wf",SUF))), Nothing,
                                (Ptr{Nothing}, Complex{$T}, Int32), psi.p, dt, method_for_B)
            end

            function add_apply_B!(this::($WF){$T}, other::($WF){$T},
                                 coefficient::Number=1.0)
               if this.m ≠ other.m
                   error("this and other must belong to the same method")
               end
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"add_apply_B_wf",SUF))), Nothing,
                                (Ptr{Nothing}, Ptr{Nothing}, Complex{$T}),
                                 this.p, other.p, coefficient)
            end

            function selfconsistent_nonlinear_step!(psi::($WF){$T}, dt::Number,
                           dt1::Number, eps::Number=100.0_prec*eps(($T)), max_iters::Integer=30)
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"selfconsistent_nonlinear_step_wf",SUF))), Nothing,
                          (Ptr{Nothing}, Complex{$T}, Complex{$T}, ($T), Int32), psi.p, dt, dt1, eps, max_itesr)
            end

            function kinetic_matrix_element(psi1::($WF){$T}, psi2::($WF){$T})
                if psi1.m ≠ psi2.m
                    error("psi1 and psi2 must belong to the same method")
                end
                ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"kinetic_matrix_element_wf",SUF))), (Complex{$T}),
                    (Ptr{Nothing}, Ptr{Nothing}), psi1.p, psi2.p )
            end

            function potential_matrix_element(psi1::($WF){$T}, psi2::($WF){$T})
                if psi1.m ≠ psi2.m
                    error("psi1 and psi2 must belong to the same method")
                end
                ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"potential_matrix_element_wf",SUF))), (Complex{$T}),
                    (Ptr{Nothing}, Ptr{Nothing}), psi1.p, psi2.p )
            end

         end #eval
    else
        @eval begin
            function imaginary_time_propagate_A!(psi::($WF){$T}, dt::Real)
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"imaginary_time_propagate_A_wf",SUF))), Nothing,
                                (Ptr{Nothing}, ($T),), psi.p, dt)
            end

            function imaginary_time_propagate_B!(psi::($WF){$T}, dt::Real, method_for_B::Integer=0)
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"imaginary_time_propagate_B_wf",SUF))), Nothing,
                                (Ptr{Nothing}, ($T), Int32), psi.p, dt, method_for_B)
            end

            function add_apply_B!(this::($WF){$T}, other::($WF){$T},
                                 coefficient::Number=1.0)
               if this.m ≠ other.m
                   error("this and other must belong to the same method")
               end
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"add_apply_B_wf",SUF))), Nothing,
                                (Ptr{Nothing}, Ptr{Nothing}, ($T)),
                                 this.p, other.p, coefficient)
            end

            function selfconsistent_nonlinear_step!(psi::($WF){$T}, dt::Number,
                           dt1::Number, eps::Number=100.0_prec*eps(($T)), max_iters::Integer=30)
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"selfconsistent_nonlinear_step_wf",SUF))), Nothing,
                          (Ptr{Nothing}, ($T), ($T), ($T), Int32), psi.p, dt, dt1, eps, max_itesr)
            end

            function kinetic_matrix_element(psi1::($WF){$T}, psi2::($WF){$T})
                if psi1.m ≠ psi2.m
                    error("psi1 and psi2 must belong to the same method")
                end
                ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"kinetic_matrix_element_wf",SUF))), ($T),
                    (Ptr{Nothing}, Ptr{Nothing}), psi1.p, psi2.p )
            end

            function potential_matrix_element(psi1::($WF){$T}, psi2::($WF){$T})
                if psi1.m ≠ psi2.m
                    error("psi1 and psi2 must belong to the same method")
                end
                ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"potential_matrix_element_wf",SUF))), ($T),
                    (Ptr{Nothing}, Ptr{Nothing}), psi1.p, psi2.p )
            end

        end #eval
    end #if

    if DIM==1
        @eval begin
            function set_potential!(m::($METHOD){$T}, f::Function)
               f_c = cfunction_check_return_type(f, ($T), (($T),))
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_potential",SUF))), Nothing,
                     (Ptr{Nothing}, Ptr{Nothing}), m.m, f_c )
            end

            function get_potential(m::($METHOD){$T}, unsafe_access::Bool=false)
               dims = zeros(Int32, 1)
               Vp = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_potential",SUF))), Ptr{$T},
                     (Ptr{Nothing}, Ptr{Int32}), m.m, dims )
               if Vp==C_NULL
                   return zeros(($T), dims[1])
               end
               V = unsafe_wrap(Array, Vp, dims[1], own=false)
               if unsafe_access
                  return V
               else
                  return copy(V)
               end
            end

            function observable(psi::($WF){$T}, f::Function)
               f_c = cfunction_check_return_type(f, ($T), (($T), ))
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"observable_wf",SUF))), ($T),
                     (Ptr{Nothing}, Ptr{Nothing} ), psi.p, f_c )
            end

        end #eval
        if COMPLEX_METHOD
        @eval begin
            function set_potential_t!(m::($METHOD){$T}, f::Function)
               f_c = cfunction_check_return_type(f, ($T), (($T),($T)))
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_potential_t",SUF))), Nothing,
                     (Ptr{Nothing}, Ptr{Nothing}, Bool), m.m, f_c, false )
            end

            function set_potential_t_derivative!(m::($METHOD){$T}, f::Function)
               f_c = cfunction_check_return_type(f, ($T), (($T),($T)))
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_potential_t",SUF))), Nothing,
                     (Ptr{Nothing}, Ptr{Nothing}, Bool), m.m, f_c, true )
            end

            function get_potential_t(m::($METHOD){$T}, t::Real)
               dims = zeros(Int32, 1)
               Vp = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_potential_t",SUF))), Ptr{$T},
                     (Ptr{Nothing}, ($T), Bool, Ptr{Int32}), m.m, t, false, dims )
               if Vp==C_NULL
                   return zeros(($T), dims[1])
               end
               V = unsafe_wrap(Array, Vp, dims[1], own=false)
               return copy(V)
            end

            function get_potential_t_derivative(m::($METHOD){$T}, t::Real)
               dims = zeros(Int32, 1)
               Vp = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_potential_t",SUF))), Ptr{$T},
                     (Ptr{Nothing}, ($T), Bool, Ptr{Int32}), m.m, t, true, dims )
               if Vp==C_NULL
                   return zeros(($T), dims[1])
               end
               V = unsafe_wrap(Array, Vp, dims[1], own=false)
               return copy(V)
            end

        end #eval
        end #if
    elseif DIM==2
        @eval begin
            function set_potential!(m::($METHOD){$T}, f::Function)
               f_c = cfunction_check_return_type(f, ($T), (($T), ($T)))
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_potential",SUF))), Nothing,
                     (Ptr{Nothing}, Ptr{Nothing}), m.m, f_c )
            end

            function get_potential(m::($METHOD){$T}, unsafe_access::Bool=false)
               dims = zeros(Int32, 2)
               Vp = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_potential",SUF))), Ptr{$T},
                     (Ptr{Nothing}, Ptr{Int32}), m.m, dims )
               if Vp==C_NULL
                   return zeros(($T), dims[1], dims[2])
               end
               V = unsafe_wrap(Array, Vp, dims[1], dims[2], own=false)
               if unsafe_access
                  return V
               else
                  return copy(V)
               end
            end

            function observable(psi::($WF){$T}, f::Function)
               f_c = cfunction_check_return_type(f, ($T), (($T), ($T)))
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"observable_wf",SUF))), ($T),
                     (Ptr{Nothing}, Ptr{Nothing} ), psi.p, f_c)
            end

        end #eval
        if COMPLEX_METHOD
        @eval begin
            function set_potential_t!(m::($METHOD){$T}, f::Function)
               f_c = cfunction_check_return_type(f, ($T), (($T),($T),($T)))
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_potential_t",SUF))), Nothing,
                     (Ptr{Nothing}, Ptr{Nothing}, Bool), m.m, f_c, false )
            end

            function set_potential_t_derivative!(m::($METHOD){$T}, f::Function)
               f_c = cfunction_check_return_type(f, ($T), (($T),($T),($T)))
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_potential_t",SUF))), Nothing,
                     (Ptr{Nothing}, Ptr{Nothing}, Bool), m.m, f_c, true )
            end

            function get_potential_t(m::($METHOD){$T}, t::Real)
               dims = zeros(Int32, 2)
               Vp = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_potential_t",SUF))), Ptr{$T},
                     (Ptr{Nothing}, ($T), Bool, Ptr{Int32}), m.m, t, false, dims )
               if Vp==C_NULL
                   return zeros(($T), dims[1], dims[2])
               end
               V = unsafe_wrap(Array, Vp, dims[2], own=false)
               return copy(V)
            end

            function get_potential_t_derivative(m::($METHOD){$T}, t::Real)
               dims = zeros(Int32, 2)
               Vp = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_potential_t",SUF))), Ptr{$T},
                     (Ptr{Nothing}, ($T), Bool, Ptr{Int32}), m.m, t, true, dims )
               if Vp==C_NULL
                   return zeros(($T), dims[1], dims[2])
               end
               V = unsafe_wrap(Array, Vp, dims[2], own=false)
               return copy(V)
            end

        end #eval
        end #if
    elseif DIM==3
        @eval begin
            function set_potential!(m::($METHOD){$T}, f::Function)
               f_c = cfunction_check_return_type(f, ($T), (($T), ($T), ($T)))
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_potential",SUF))), Nothing,
                     (Ptr{Nothing}, Ptr{Nothing}), m.m, f_c )
            end

            function get_potential(m::($METHOD){$T}, unsafe_access::Bool=false)
               dims = zeros(Int32, 3)
               Vp = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_potential",SUF))), Ptr{$T},
                     (Ptr{Nothing}, Ptr{Int32}), m.m, dims )
               if Vp==C_NULL
                   return zeros(($T), dims[1], dims[2], dims[3])
               end
               V = unsafe_wrap(Array, Vp, dims[1], dims[2], dims[3], own=false)
               if unsafe_access
                  return V
               else
                  return copy(V)
               end
            end

            function observable(psi::($WF){$T}, f::Function)
               f_c = cfunction_check_return_type(f, ($T), (($T), ($T), ($T)))
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"observable_wf",SUF))), ($T),
                     (Ptr{Nothing}, Ptr{Nothing} ), psi.p, f_c )
            end

        end #eval
        if COMPLEX_METHOD
        @eval begin
            function set_potential_t!(m::($METHOD){$T}, f::Function)
               f_c = cfunction_check_return_type(f, ($T), (($T),($T),($T)($T)))
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_potential_t",SUF))), Nothing,
                     (Ptr{Nothing}, Ptr{Nothing}, Bool), m.m, f_c, false )
            end

            function set_potential_t_derivative!(m::($METHOD){$T}, f::Function)
               f_c = cfunction_check_return_type(f, ($T), (($T),($T),($T)($T)))
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_potential_t",SUF))), Nothing,
                     (Ptr{Nothing}, Ptr{Nothing}, Bool), m.m, f_c, true )
            end

            function get_potential_t(m::($METHOD){$T}, t::Real)
               dims = zeros(Int32, 3)
               Vp = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_potential_t",SUF))), Ptr{$T},
                     (Ptr{Nothing}, ($T), Bool, Ptr{Int32}), m.m, t, false, dims )
               if Vp==C_NULL
                   return zeros(($T), dims[1], dims[2], dims[3])
               end
               V = unsafe_wrap(Array, Vp, dims[3], own=false)
               return copy(V)
            end

            function get_potential_t_derivative(m::($METHOD){$T}, t::Real)
               dims = zeros(Int32, 3)
               Vp = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_potential_t",SUF))), Ptr{$T},
                     (Ptr{Nothing}, ($T), Bool, Ptr{Int32}), m.m, t, true, dims )
               if Vp==C_NULL
                   return zeros(($T), dims[1], dims[2], dims[3])
               end
               V = unsafe_wrap(Array, Vp, dims[3], own=false)
               return copy(V)
            end

        end #eval
        end #if
    end #if

end # for


text:
# automatically generated by the FlatBuffers compiler, do not modify

MyGame.Example.eval(quote

FlatBuffers.@with_kw mutable struct TestSimpleTableWithEnum
    color::Color = 2
end
FlatBuffers.@ALIGN(TestSimpleTableWithEnum, 1)
FlatBuffers.slot_offsets(::Type{T}) where {T<:TestSimpleTableWithEnum} = [
    0x00000004
]

TestSimpleTableWithEnum(buf::AbstractVector{UInt8}) = FlatBuffers.read(TestSimpleTableWithEnum, buf)
TestSimpleTableWithEnum(io::IO) = FlatBuffers.deserialize(io, TestSimpleTableWithEnum)

end)

max_stars_repo_name: rjkat/flatbuffers-julia

max_issues_repo_name: rjkat/flatbuffers-julia

max_forks_repo_name: rjkat/flatbuffers-julia


text:
# Autogenerated wrapper script for PCRE2_jll for aarch64-linux-musl
export libpcre2_16, libpcre2_32, libpcre2_8

JLLWrappers.@generate_wrapper_header("PCRE2")
JLLWrappers.@declare_library_product(libpcre2_16, "libpcre2-16.so.0")
JLLWrappers.@declare_library_product(libpcre2_32, "libpcre2-32.so.0")
JLLWrappers.@declare_library_product(libpcre2_8, "libpcre2-8.so.0")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libpcre2_16,
        "lib/libpcre2-16.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcre2_32,
        "lib/libpcre2-32.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpcre2_8,
        "lib/libpcre2-8.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Statistics:mean, std
using SpecialFunctions:erf

"""
    scale_minmax(scores)

Transform an array of scores into a range between [0,1] using min-max scaling.

Parameters
----------
      scores::Tuple{Scores, Scores}
A tuple consisting of two vectors representing training and test scores.

Returns
----------
    normalized_scores::Tuple{Scores, Scores}
The normalized train and test scores.

Examples
----------
    scores_train, scores_test = ([1, 2, 3], [4, 3, 2, 1, 0])
    scale_minmax(scores_train, scores_test) # ([0.0, 0.5, 1.0], [1.0, 1.0, 0.5, 0.0, 0.0])
"""
function scale_minmax(scores::Tuple{Scores, Scores})::Tuple{Scores, Scores}
    scores_train, scores_test = scores
    minTrain, maxTrain = extrema(scores_train)
    @assert minTrain < maxTrain "Cannot normalize scores if they are all equal"
    f = scores -> clamp.((scores .- minTrain) ./ (maxTrain - minTrain), 0, 1)
    f(scores_train), f(scores_test)
end

"""
    scale_unify(scores)

Transform an array of scores into a range between [0,1] using unifying scores as described in [1].

Parameters
----------
    scores::Tuple{Scores, Scores}
A tuple consisting of two vectors representing training and test scores.

Returns
----------
    unified_scores::Tuple{Scores, Scores}
The unified train and test scores.

Examples
----------
    scores_train, scores_test = ([1, 2, 3], [4, 3, 2, 1, 0])
    unify(scores_train, scores_test) # ([0.0, 0.0, 0.68..], [0.95.., 0.68.., 0.0, 0.0, 0.0])

References
----------
Kriegel, Hans-Peter; Kroger, Peer; Schubert, Erich; Zimek, Arthur (2011): Interpreting and Unifying Outlier Scores.
"""
function scale_unify(scores::Tuple{Scores, Scores})::Tuple{Scores, Scores}
    scores_train, scores_test = scores
    μ, σ = mean(scores_train), std(scores_train)
    @assert σ > 0 "Cannot normalize scores if they are all equal"
    f = scores -> clamp.(erf.((scores .- μ) ./ (σ * √2)), 0, 1)
    f(scores_train), f(scores_test)
end


text:
## PART 1

println("~~ QUESTION 4 ~~");
println("\nPart 1");

# Read in values from previous simulation (but only after thermalization):
println("\n\tReading in values from MD simulation...");

temps_1069 = readdlm("q4/temps-rho-0.75-temp-1.069.dat")[401:end]
temps_1304 = readdlm("q4/temps-rho-0.75-temp-1.304.dat")[401:end]

pe_1069 = readdlm("q4/pe-rho-0.75-temp-1.069.dat")[401:end]
pe_1304 = readdlm("q4/pe-rho-0.75-temp-1.304.dat")[401:end]

vir_1069 = readdlm("q4/vir-rho-0.75-temp-1.069.dat")[401:end]
vir_1304 = readdlm("q4/vir-rho-0.75-temp-1.304.dat")[401:end]

# Slap it all into a single array so that we can reuuse
#   old code
md_1069 = [temps_1069 pe_1069 vir_1069];
md_1304 = [temps_1304 pe_1304 vir_1304];