text:
using LinearAlgebra
using Parameters

########## Base Env ##########
# BaseEnv
struct BaseEnv <: AbstractEnv
    initial_state
end
BaseEnv(given_size::Tuple) = BaseEnv(zeros(given_size))
BaseEnv(args...) = BaseEnv(zeros(args...))
BaseEnv() = BaseEnv(0.0)
# initial_condition
NestedEnvironments.initial_condition(env::BaseEnv) = env.initial_state

########## InputAffineQuadraticCostEnv ##########
"""
An example of continuous-time nonlinear dynamical system introduced in
several studies on approximate dynamic programming.
# Reference
[1] K. G. Vamvoudakis and F. L. Lewis, “Online Actor-Critic Algorithm to Solve the Continuous-Time Infinite Horizon Optimal Control Problem,” Automatica, vol. 46, no. 5, pp. 878–888, 2010, doi: 10.1016/j.automatica.2010.02.018.
[2] V. Nevistic and J. A. Primbs, “Constrained Nonlinear Optimal Control: a Converse HJB Approach,” 1996.
"""
@with_kw struct InputAffineQuadraticCostEnv <: AbstractEnv
    Q = Matrix(I, 2, 2)
    R = 1
    P = [0.5 0; 0 1]
end

function f(env::InputAffineQuadraticCostEnv, x)
    x1 = x[1]
    x2 = x[2]
    f1 = -x1 + x2
    f2 = -0.5*x1 -0.5*x2*(1-(cos(2*x1)+2)^2)
    return  [f1, f2]
end

function g(env::InputAffineQuadraticCostEnv, x)
    x1 = x[1]
    g1 = 0
    g2 = cos(2*x1) +2
    return [g1, g2]
end

"""
    r(env, x, u)

Calculate running cost, i.e., V = ∫r dt.
"""
function r(env::InputAffineQuadraticCostEnv, x, u)
    x'*env.Q*x + u'*env.R*u
end

function V_optimal(env::InputAffineQuadraticCostEnv, x)
    return x'*env.P*x
end

function u_optimal(env::InputAffineQuadraticCostEnv, x)
    _g = g(env, x)
    return -_g'*env.P*x
end

function ẋ(env::InputAffineQuadraticCostEnv, x, t, u)
    ẋ = f(env, x) + g(env, x)*u
    return ẋ
end


text:
module GraphKernels

using Graphs
using SimpleValueGraphs
using SimpleValueGraphs: AbstractValGraph
using LinearAlgebra: dot, diag
using Statistics: mean, std
using LIBSVM
using Random: MersenneTwister, randperm
using ThreadsX

import LIBSVM: svmtrain, svmpredict

using KernelFunctions: kernelmatrix, kernelmatrix_diag
import KernelFunctions

export
    AbstractGraphKernel,

    NoGraphBaselineGraphKernel,
    ShortestPathGraphKernel,
    PyramidMatchGraphKernel,
    WeisfeilerLehmanGraphKernel,

    ConstVertexKernel,
    DiracVertexKernel,
    DotVertexKernel,

    k_fold_cross_validation,

    # reexport from KernelFunctions
    kernelmatrix,
    kernelmatrix_diag,

    # overridden methods from LIBSVM
    svmtrain,
    svmpredict

include("utils.jl")
include("replacedvertexvals.jl")
include("vertex_kernels.jl")
include("graph-kernels/abstract-graph-kernel.jl")
include("graph-kernels/baseline-graph-kernel.jl")
include("graph-kernels/pyramid-match-graph-kernel.jl")
include("graph-kernels/shortest-path-graph-kernel.jl")
include("graph-kernels/weisfeiler-lehman-graph-kernel.jl")
include("integrations/LIBSVM.jl")

# ================================================================
#     utilities
# ================================================================

"""
    k_fold_cross_validation

Simple k-fold cross validation implementation for quick testing during development.
"""
function k_fold_cross_validation(kernel::KernelFunctions.Kernel, graphs::AbstractVector{<:AbstractGraph}; k_folds=5, class_key=1, kwargs...)

    n = length(graphs)
    indices = randperm(MersenneTwister(123), n)

    acc_train = Float64[]
    acc_valid = Float64[]
    # TODO there should be better partition function that
    # ensures that not only the last partition is smaller
    for valid_indices in Iterators.partition(1:n, ceil(Int, n // k_folds))
        train_x = [graphs[indices[i]] for i in 1:n if i ∉ valid_indices]
        valid_x = [graphs[indices[i]] for i in valid_indices]

        train_y = [get_graphval(g, class_key) for g in train_x]
        valid_y = [get_graphval(g, class_key) for g in valid_x]

        model = svmtrain(train_x, train_y, kernel; kwargs...)

        train_y_pred = svmpredict(model, train_x)
        valid_y_pred = svmpredict(model, valid_x)

        push!(acc_train, mean(train_y_pred .== train_y))
        push!(acc_valid, mean(valid_y_pred .== valid_y))
    end

    return (mean_train_accuracy=mean(acc_train), std_train_accuracy=std(acc_train),
            mean_valid_accuracy=mean(acc_valid), std_valid_accuracy=std(acc_valid))
end

end

ext: jl


text:
# forward content properties
cache_key(tb::TagBlock) = cache_key(content(tb))
occupied_locs(x::TagBlock) = occupied_locs(content(x))

Base.:(==)(a::TB, b::TB) where {TB<:TagBlock} = content(a) == content(b)


text:
# This file defines a mechanism for auto-sniffing, loading and saving spectrum files

abstract type SpectrumFileType end

function savespectrum(ty::Type{<:SpectrumFileType}, fio, data)
    @error "Saving to is not implemented for $ty. Probably never will be."
end

function loadspectrum(ty::Type{<:SpectrumFileType}, filename::AbstractString; kwargs...)
    badname(sp) =
        (!haskey(sp, :Name)) ||
        startswith(sp[:Name], "Bruker") ||
        startswith(sp[:Name], "Spectrum[")
    @assert isfile(filename) "$filename is not a file in loadspectrum(...)."
    return open(filename, read = true) do ios
        spec = loadspectrum(ty, ios; kwargs...)
        spec[:Filename] = filename
        spec[:Name] = badname(spec) ? splitext(splitpath(filename)[end])[1] : spec[:Name]
        return spec
    end
end

loadspectrum(
    ty::Type{<:SpectrumFileType},
    filename::AbstractString,
    det::EDSDetector;
    kwargs...,
) = apply(loadspectrum(ty, filename; kwargs...), det)

savespectrum(ty::Type{<:SpectrumFileType}, filename::AbstractString, spec::Spectrum) =
    open(filename, write = true) do ios
        savespectrum(ty, ios, spec)
        spec[:Filename] = filename
    end

# ---------------------------------------------------------------------------------------------------#

struct BrukerSPX <: SpectrumFileType end

loadspectrum(::Type{BrukerSPX}, ios::IO) = readbrukerspx(ios)
extensions(::Type{BrukerSPX}) = (".spx",)

sniff(::Type{BrukerSPX}, ios::IO) = detectbrukerspx(ios)

# ---------------------------------------------------------------------------------------------------#

struct BrukerPDZ <: SpectrumFileType end

loadspectrum(::Type{BrukerPDZ}, ios::IO) = readbrukerpdz(ios)
extensions(::Type{BrukerPDZ}) = (".pdz",)

sniff(::Type{BrukerPDZ}, ios::IO) = detectbrukerpdz(ios)

# ---------------------------------------------------------------------------------------------------#

struct ISOEMSA <: SpectrumFileType end

loadspectrum(::Type{ISOEMSA}, ios::IO, ty::Type{<:Real} = Float64) = readEMSA(ios, ty)
savespectrum(::Type{ISOEMSA}, ios::IO, spec::Spectrum) = writeEMSA(ios, spec)

extensions(::Type{ISOEMSA}) = (".msa", ".emsa", ".ems", ".txt")

sniff(::Type{ISOEMSA}, ios::IO) = isemsa(ios)

# ---------------------------------------------------------------------------------------------------#

struct ASPEXTIFF <: SpectrumFileType end

loadspectrum(::Type{ASPEXTIFF}, ios::IO; withImgs = false, astype::Type{<:Real} = Float64) =
    readAspexTIFF(ios, withImgs = withImgs, astype = astype)
extensions(::Type{ASPEXTIFF}) = (".tif",)

sniff(::Type{ASPEXTIFF}, ios::IO) = detectAspexTIFF(ios)

# ---------------------------------------------------------------------------------------------------#

const spectrumfiletypes = (ISOEMSA, BrukerSPX, ASPEXTIFF, BrukerPDZ)

function sniffspectrum(ios::IO)
    for st in spectrumfiletypes
        try
            seekstart(ios)
            if sniff(st, ios)
                seekstart(ios)
                return st
            end
        catch
            @info "Error sniffing filetype $st"
        end
    end
    @error "This object does not represent a known spectrum file type."
end

function sniffspectrum(filename::AbstractString)
    @assert isfile(filename) "$filename is not a file in sniffspectrum(...)."
    ext = lowercase(splitext(filename)[2])
    for st in spectrumfiletypes
        if ext in extensions(st)
            try
                res = open(filename, read = true) do ios
                    sniff(st, ios) ? st : nothing
                end
                if !isnothing(res)
                    return res
                end
            catch
                @info "Error sniffing filetype $st"
            end
        end
    end
    @error "This object does not represent a known spectrum file type."
end

loadspectrum(ios::IO) = loadspectrum(sniffspectrum(ios), ios)
loadspectrum(filename::AbstractString) = loadspectrum(sniffspectrum(filename), filename)
loadspectrum(filename::AbstractString, det::EDSDetector) =
    loadspectrum(sniffspectrum(filename), filename, det)


text:
"""
    rbkiter(;input::AbstractString="", outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractString="ftt", pseudocount::Number=1f0, rowmeanlist::AbstractString="", rowvarlist::AbstractString="",colsumlist::AbstractString="", dim::Number=3, numepoch::Number=10, lower::Number=0, upper::Number=1.0f+38, expvar::Number=0.1f0, initW::Union{Nothing,AbstractString}=nothing, initV::Union{Nothing,AbstractString}=nothing, logdir::Union{Nothing,AbstractString}=nothing, perm::Bool=false, cper::Number=1f0)

Randomized Block Krylov Iteration.

Input Arguments
---------
- `input` : Julia Binary file generated by `OnlinePCA.csv2bin` function.
- `outdir` : The directory specified the directory you want to save the result.
- `scale` : {log,ftt,raw}-scaling of the value.
- `pseudocount` : The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated.
- `rowmeanlist` : The mean of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `rowvarlist` : The variance of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `colsumlist` : The sum of counts of each columns of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `dim` : The number of dimension of PCA.
- `numepoch` : The number of epoch.
- `lower` : Stopping Criteria (When the relative change of error is below this value, the calculation is terminated)
- `upper` : Stopping Criteria (When the relative change of error is above this value, the calculation is terminated)
- `initW` : The CSV file saving the initial values of eigenvectors.
- `initV` : The CSV file saving the initial values of loadings.
- `logdir` : The directory where intermediate files are saved, in every evalfreq (e.g. 5000) iteration.
- `perm` : Whether the data matrix is shuffled at random.
- `cper` : Count per X (e.g. CPM: Count per million <1e+6>)

Output Arguments
---------
- `W` : Eigen vectors of covariance matrix (No. columns of the data matrix × dim)
- `λ` : Eigen values (dim × dim)
- `V` : Loading vectors of covariance matrix (No. rows of the data matrix × dim)
- `Scores` : Principal component scores
- `ExpVar` : Explained variance by the eigenvectors
- `TotalVar` : Total variance of the data matrix
- stop : Whether the calculation is converged
"""
function rbkiter(;input::AbstractString="", outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractString="ftt", pseudocount::Number=1f0, rowmeanlist::AbstractString="", rowvarlist::AbstractString="",colsumlist::AbstractString="", dim::Number=3, numepoch::Number=10, lower::Number=0, upper::Number=1.0f+38, expvar::Number=0.1f0, initW::Union{Nothing,AbstractString}=nothing, initV::Union{Nothing,AbstractString}=nothing, logdir::Union{Nothing,AbstractString}=nothing, perm::Bool=false, cper::Number=1f0)
    # Initial Setting
    pca = RBKITER()
    N, M = nm(input)
    pseudocount, W, X, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar, lower, upper = init(input, pseudocount, dim, rowmeanlist, rowvarlist, colsumlist, initW, initV, logdir, pca, lower, upper, cper, scale)
    tmpW = zeros(Float32, M, dim)
    K = zeros(Float32, M, dim*numepoch)
    QtAt = zeros(Float32, dim*numepoch, N)
    QtAtA = zeros(Float32, dim*numepoch, M)
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    x = zeros(UInt32, M)
    normx = zeros(Float32, M)
    # If not 0 the calculation is converged
    stop = 0
    s = 1
    n = 1
    # Each epoch s
    progress = Progress(numepoch*N)
    while s <= numepoch
        open(input) do file
            stream = ZstdDecompressorStream(file)
            read!(stream, tmpN)
            read!(stream, tmpM)
            # Each step n
            while n <= N
                next!(progress)
                # Row vector of data matrix
                read!(stream, x)
                normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, cper)
                if perm
                    normx .= normx[randperm(length(normx))]
                end
                # Power iteration
                tmpW .+= normx * (normx' * W)
                n += 1
            end
            close(stream)
        end
        # Update W
        W .+= tmpW
        tmpW = zeros(Float32, M, dim)
        # NaN
        checkNaN(W, GD())
        # Retraction
        W .= Array(qr!(W).Q)
        # Check Float32
        @assert W[1,1] isa Float32
        # Save to K
        K[:, (s-1)*dim+1:s*dim] = W
        # save log file
        if logdir isa String
            stop = outputlog(s, input, dim, logdir, W, GD(), TotalVar, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, lower, upper, stop, cper)
        end
        s += 1
        if n == N + 1
            n = 1
        end
    end
    Q = Array(qr!(K).Q)
    n = 1
    println("Q' A")
    open(input) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, tmpN)
        read!(stream, tmpM)
        # Each step n
        while n <= N
            next!(progress)
            # Row vector of data matrix
            read!(stream, x)
            normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, cper)
            if perm
                normx .= normx[randperm(length(normx))]
            end
            # Power iteration
            QtAt .+= Q'*normx
            n += 1
        end
        close(stream)
    end
    n = 1
    println("Q' A' A")
    open(input) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, tmpN)
        read!(stream, tmpM)
        # Each step n
        while n <= N
            next!(progress)
            # Row vector of data matrix
            read!(stream, x)
            normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, cper)
            if perm
                normx .= normx[randperm(length(normx))]
            end
            # Power iteration
            QtAtA .+= QtAt[:,n]*normx'
            n += 1
        end
        close(stream)
    end

    # SVD
    println("svd(Q' A' A)")
    U, Σ, V = svd(QtAtA)
    println("Q*U")
    W = Q*U[:,1:dim]
    # Return, W, λ, V
    out = WλV(W, input, dim, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, TotalVar, cper)
    out = (out[1], out[2], out[3], out[4], out[5], out[6], stop)
    if outdir isa String
        output(outdir, out, expvar)
    end
    return out
end


text:
for P in (:Year, :Month, :Week, :Day)
    @eval begin
        @inline $P(td::TimeDate) = $P(on_date(td))
        @inline $P(tdz::TimeDateZone) = $P(on_date(ZonedDateTime(tdz)))
    end
end
for P in (:Hour, :Minute, :Second, :Millisecond)
    @eval begin
        @inline $P(td::TimeDate) = $P(at_time(td))
        @inline $P(tdz::TimeDateZone) = $P(at_time(ZonedDateTime(tdz)))
    end
end
for P in (:Microsecond, :Nanosecond)
    @eval begin
        @inline $P(td::TimeDate) = $P(at_time(td))
        @inline $P(tdz::TimeDateZone) = $P(TimeDate(tdz))
    end
end

for P in (:Years, :Months)
    @eval begin
        @inline $P(td::TimeDate) = $P(on_date(td))
        @inline $P(tdz::TimeDateZone) = $P(on_date(ZonedDateTime(tdz)))
    end
end

Weeks(td::TimeDate) = Weeks(DateTime(td))
Weeks(tdz::TimeDateZone) = Weeks(DateTime(td))
Days(td::TimeDate) = Days(DateTime(td))
Days(tdz::TimeDateZone) = Days(DateTime(td))

for P in (:Hours, :Minutes, :Seconds, :Milliseconds)
    @eval begin
        @inline $P(td::TimeDate) = $P(at_time(td))
        @inline $P(tdz::TimeDateZone) = $P(at_time(ZonedDateTime(tdz)))
    end
end
for P in (:Microseconds, :Nanoseconds)
    @eval begin
        @inline $P(td::TimeDate) = $P(at_time(td))
        @inline $P(tdz::TimeDateZone) = $P(TimeDate(tdz))
    end
end

for (P, Q) in ((:Year,:year), (:Month,:month), (:Week,:week), (:Day,:day),
               (:Hour,:hour), (:Minute,:minute), (:Second,:second),
               (:Millisecond,:millisecond), (:Microsecond,:microsecond),
               (:Nanosecond,:nanosecond))
    @eval begin
        @inline $Q(td::TimeDate) = $P(td).value
        @inline $Q(tdz::TimeDateZone) = $P(tdz).value
    end
end

for (P, Q) in ((:Years,:years), (:Months,:months), (:Weeks,:weeks), (:Days,:days),
               (:Hours,:hours), (:Minutes,:minutes), (:Seconds,:seconds),
               (:Milliseconds,:milliseconds), (:Microseconds,:microseconds),
               (:Nanoseconds,:nanoseconds))
    @eval begin
        @inline $Q(td::TimeDate) = $P(td).value
        @inline $Q(tdz::TimeDateZone) = $P(tdz).value
    end
end


text:
"""
    readXYZ(io::IO) -> Vector{AtomPosition{3}}

Reads an XYZ file into a `Vector{AtomPosition{3}}`.
"""
function readXYZ(io::IO)
    # Loop through each line in the file
    itr = eachline(io)
    # Get the number of atoms from the first line
    n_at = parse(Int, iterate(itr)[1])
    v = Vector{AtomPosition{3}}(undef, n_at)
    # Skip the second line
    iterate(itr)
    # Loop through the rest
    for (n, ln) in enumerate(itr)
        # If the line is empty, skip it
        isempty(ln) && continue
        # Split up the current line
        entries = split(ln)
        # Get the name (first entry)
        name = entries[1]
        # Get the position (next three entries)
        pos = parse.(Float64, entries[2:4])
        # Add to vector
        v[n] = AtomPosition{3}(name, pos)
    end
    return v
end

readXYZ(filename::AbstractString) = open(readXYZ, filename)

"""
    readXSF3D(
        io::IO;
        spgrp::Integer = 0,
        origin::AbstractVector{<:Real} = [0, 0, 0]
        ctr::Symbol = :P
    ) -> CrystalWithDatasets{3}

Reads in an XCrysDen XSF file from an input stream and returns a `CrystalWithDatasets{3}` with all
datasets that have been included within the file.

Space group and origin information are not supplied in XSF files, but they can be supplied using
the `spgrp` and `origin` keyword arguments. Centering information can be provided using the `ctr`
argument, but is overridden by a space group assignment.
"""
function readXSF3D(
    io::IO;
    spgrp::Integer = 0,
    origin::AbstractVector{<:Real} = [0, 0, 0],
    ctr::Symbol = :P
)
    # Function for getting lattice basis vectors
    function getlattice!(itr)
        vecs = [parse.(Float64, split(iterate(itr)[1])) for n in 1:3]
        @debug string("Found vectors:\n", vecs)
        return BasisVectors{3}(hcat(vecs...))
    end
    # Function for getting 3D lists of atoms
    function getatoms!(itr, basis, natom)
        # Vector of AtomPositions
        apos = Vector{AtomPosition{3}}(undef, natom)
        for n in 1:natom
            # Get next line
            entries = split(iterate(itr)[1])
            # Get atomic number (first entry)
            num = parse(Int, entries[1])
            # Get atomic position (next 3)
            pos = SVector{3,Float64}(parse.(Float64, entries[s]) for s in 2:4)
            # assign atom position struct
            apos[n] = AtomPosition{3}(num, basis\pos)
        end
        return AtomList{3}(basis, apos)
    end
    # Function to get grid data
    # By default, trim the edges of the grid (repeated points)
    function getgrid!(itr, sz; trim=true)
        # Store everything in this vector (to be reshaped)
        v = Vector{Float64}(undef, prod(sz))
        # Entry counter
        entry = 0
        # Line to be processed
        ln = ""
        while entry < prod(sz) || !contains(ln, "END_DATAGRID")
            # Get the next line
            ln = iterate(itr)[1]
            # Get the data entries
            dt = filter(!isnothing, [tryparse(Float64, s) for s in split(ln)])
            # Append the entries to the vector
            v[entry .+ (1:length(dt))] = dt
            # Move the entry counter over
            entry = entry + length(dt)
        end
        # Reshape and trim the edges of the grid, if needed
        return reshape(v, sz[1], sz[2], sz[3])[1:end - trim, 1:end - trim, 1:end - trim]
    end
    # Manually iterate through the file
    # This allows lines to be skipped
    iter = eachline(io)
    # Counter for debugging
    count = 0
    # Preallocated variables
    # Basis vectors
    prim = zeros(BasisVectors{3})
    conv = zeros(BasisVectors{3})
    # Atomic positions
    local atom_list::AtomList{3}
    data = Dict{String,RealSpaceDataGrid{3,Float64}}()
    # There is a break statement to get of out of this loop
    while true
        # Get the line contents
        ln = let i = iterate(iter)
            # Iteration returns nothing once it's out of lines, so break from loop
            isnothing(i) ? break : i[1]
        end
        count += 1
        # Check for empty lines or comment lines; advance if present
        # Comments are only allowed between blocks in a valid XSF file
        # No need to check within the loops below
        (isempty(ln) || startswith(ln, "#") || all(isspace, ln)) && continue
        # The block is supposed to be given as "BEGIN_BLOCK_DATAGRID_3D"
        # But bin2xsf seems to write "BEGIN_BLOCK_DATAGRID3D" (missing underscore)
        if contains(ln, "BEGIN_BLOCK_DATAGRID") && contains(ln, "3D")
            # Skip this line
            ln = iterate(iter)[1]
            @debug "Line contents: \"" * ln * "\""
            # Loop until the end of the block
            while !contains(ln, "END_BLOCK")
                # Check for a datagrid entry
                if contains(ln, "DATAGRID_3D")
                    # Get the name of the datagrid
                    @debug "Line contents: \"" * ln * "\""
                    name = split(ln, "DATAGRID_3D_")[2]
                    # Get the dimensions of the datagrid
                    ln = iterate(iter)[1]
                    dims = parse.(Int, split(ln))
                    grid = Array{Float64,3}(undef, dims...)
                    # Get the shift off the origin
                    ln = iterate(iter)[1]
                    orig = parse.(Float64, split(ln))
                    # Get the lattice
                    latt = getlattice!(iter)
                    # Get the datagrid
                    # TODO: Determine if trimming is needed based on periodicity
                    # Periodic structures should get trimmed
                    grid = getgrid!(iter, dims, trim=true)
                    # Export the data in a RealSpaceDataGrid
                    data[name] = RealSpaceDataGrid(latt, orig, grid)
                end
                ln = iterate(iter)[1]
            end
        # Get the primitive cell vectors
        elseif contains(ln, "PRIMVEC")
            prim = getlattice!(iter)
        # Get the conventional cell vectors (if present)
        elseif contains(ln, "CONVVEC")
            conv = getlattice!(iter)
        # Get primitive cell atomic coordinates
        elseif contains(ln, "PRIMCOORD")
            # Get the number of atoms
            ln = iterate(iter)[1]
            natom = parse(Int, split(ln)[1])
            @debug string("natom = ", natom)
            # The next lines have all of the atoms
            atom_list = getatoms!(iter, prim, natom)
        elseif contains(ln, "CONVCOORD") && isempty(atom_list)
            # TODO: complete this, but it's low priority
            # The PRIMCOORD block gets everything we need right now
        elseif contains(ln, "ATOMS") && isempty(atom_list)
            # TODO: complete this, but it's low priority
            # The PRIMCOORD block gets everything we need right now
        elseif contains(ln, "DIM-GROUP")
            # Check the line below and see if the structure is periodic
            ln = iterate(iter)[1]
            if first(strip(ln)) == '3' && spgrp == 0
                # If the structure is a crystal, set spgrp to 1
                # spgrp 0 implies a non-periodic structure
                spgrp = 1
            end
        elseif contains(ln, "CRYSTAL") && spgrp == 0
            # If the structure is a crystal, set spgrp to 1
            # spgrp 0 implies a non-periodic structure
            spgrp = 1
        end
    end
    # If the conventional cell hasn't been defined, generate it
    if iszero(conv)
        conv = BasisVectors{3}(REDUCTION_MATRIX_3D[ctr] \ matrix(prim))
    end
    # Generate the real space lattice
    latt = RealLattice{3}(prim, conv)
    xtal = Crystal{3}(latt, spgrp, origin, atom_list, atom_list)
    return CrystalWithDatasets{3,String,RealSpaceDataGrid{3}}(xtal, data)
end

"""
    readXSF3D(
        filename::AbstractString;
        spgrp::Integer = 0,
        origin::AbstractVector{<:Real} = [0, 0, 0]
        ctr::Symbol = :P
    ) -> CrystalWithDatasets{3}

Reads an XSF file at path `filename`.
"""
function readXSF3D(filename::AbstractString; kwargs...)
    return open(filename) do io
        readXSF3D(io; kwargs...)
    end
end

"""
    writeXSF(io, xtal::Crystal{D})

Writes the crystal component of an XCrysDen XSF file.
"""
function writeXSF(io::IO, xtal::Crystal{D}) where D
    println(io, "# Written by Xtal.jl")
    # Dimension information
    println(io, "DIM_GROUP")
    println(io, D, "  1")
    # Primitive cell vectors
    println(io, "PRIMVEC")
    for (n, x) in enumerate(matrix(prim(xtal)))
        @printf(io, "%20.14f  ", x)
        n % D == 0 && println(io)
    end
    # Coordinates of the atoms in the primitive cell
    println(io, "PRIMCOORD")
    # Throw warning if the space group number is greater than 1
    # This probably means the generating set of atoms is incomplete
    xtal.sgno > 1 && @warn "Symmetrization of the unit cell has not been implemented yet."
    # Print the number of atoms in the structure
    println(io, lpad(string(natom(xtal.gen)), 6), "    1")
    # Print all of the atoms
    for a in cartesian(xtal.gen)
        @printf(io, "%6i", atomicno(a))
        for x in coord(a)
            @printf(io, "%20.14f  ", x)
        end
        println(io)
    end
end

"""
    writeXSF(io::IO, key, data::RealSpaceDataGrid{D,T}; periodic=true)

Writes the crystal component of an XCrysDen XSF file. By default, automatic wrapping of the
datagrid occurs (values are repeated at the end of each dimension).
"""
function writeXSF(io::IO, key, data::RealSpaceDataGrid{D,T}; periodic=true) where {D,T}
    println(io, "    DATAGRID_", D, "D_", key)
    # Print the grid size (+1)
    println(io, " "^8, join([rpad(string(n + 1), 8) for n in gridsize(data)]))
    # Print the grid offset
    print(io, " "^4)
    for x in data.orig
        @printf(io, "%20.14f", x)
    end
    print(io, "\n" * " "^4)
    # Print the basis vectors for the grid
    for (n, x) in enumerate(matrix(basis(data)))
        @printf(io, "%20.14f", x)
        n % D == 0 && print(io, "\n" * " "^4)
    end
    # This should generate a view that wraps around for a "general grid"
    if periodic
        g = view(grid(data), (1 .+ (0:n) .% n for n in gridsize(data))...)
    else
        g = view(M, (1:n for n in size(M))...)
    end
    # Print the actual data in the grid
    for (n,x) in enumerate(g)
        @printf(io, "%20.14f", x)
        n % 4 == 0 && print(io, "\n" *  " "^4)
    end
    # End the datagrid
    length(g) % 4 == 0 || println(io)
    println(io, "    END_DATAGRID_", D, "D")
end

function writeXSF(
    io::IO,
    data::Dict{K, RealSpaceDataGrid{D,T}};
    periodic=true
) where {K,D,T}
    println(io, "BEGIN_BLOCK_DATAGRID_", D, "D")
    println(io, "Written by Xtal.jl")
    for (k,d) in data
        writeXSF(io::IO, k, d, periodic=periodic)
    end
    println(io, "END_BLOCK_DATAGRID_", D, "D")
end

"""
    writeXSF(io::IO, xtaldata::CrystalWithDatasets{D,K,V})

Writes a `CrystalWithDatasets` to an XCrysDen XSF file.
"""
function writeXSF(
    io::IO,
    xtaldata::CrystalWithDatasets;
    periodic=true
)
    writeXSF(io, Crystal(xtaldata))
    writeXSF(io, data(xtaldata), periodic=periodic)
end

function writeXSF(filename::AbstractString, xtal::Crystal)
    open(filename, write=true) do io
        writeXSF(io, xtal)
    end
end

function writeXSF(filename::AbstractString, xtaldata::CrystalWithDatasets; periodic=true)
    open(filename, write=true) do io
        writeXSF(io, xtaldata, periodic=periodic)
    end
end

"""
    readCPcoeff(io::IO, Lmax::Val{L}) -> SphericalComponents{L}

Reads in the spherical harmonic projection coefficients from a CPpackage2 calculation.

By default, CPpackage2 gives the coefficients for spherical harmonics up to a maximum l value of 6.
"""
function readCPcoeff(io::IO, Lmax::Val{L}=Val(6)) where L
    # All the data should be in a Vector{Float64} with this one line
    data = parse.(Float64, [v[2] for v in split.(readlines(io))])
    @debug "$(length(data)) lines in file"
    # Number of spherical harmonic coefficients
    ncoeff = (L + 1)^2
    natom = div(length(data), ncoeff)
    @debug "ncoeff = $ncoeff, natom = $natom"
    return [SphericalComponents{L}(data[(n - 1)*ncoeff .+ (1:ncoeff)]) for n in 1:natom]
end

readCPcoeff(filename::AbstractString) = open(readCPcoeff, filename)

"""
    readCPgeo(io::IO) -> Vector{AtomPosition{3}}

Reads the atomic positions used for a CPpackage2 calculation.
"""
function readCPgeo(io::IO)
    lns = readlines(io)
    atomnames = [v[1] for v in split.(lns)]
    positions = [parse.(Float64, v[2:4]) for v in split.(lns)]
    return AtomPosition{3}.(atomnames, positions)
end

readCPgeo(filename::AbstractString) = open(readCPgeo, filename)

"""
    readCPcell(io::IO) -> BasisVectors{3}

Reads the basis vectors of the unit cell used for a CPpackage2 calculation.
"""
function readCPcell(io::IO)
    matrix = hcat([parse.(Float64, v) for v in split.(readlines(io))]...)
    return BasisVectors{3}(matrix)
end

readCPcell(filename::AbstractString) = open(readCPcell, filename)

# Include code for processing specific types of files
include("software/abinit.jl")
include("software/vasp.jl")


text:
using Elly

function main()
    dfs = HDFSClient("h6-hdfs-master", 8020)
    readdir(dfs)

end


text:
function ChangesOfVariables.with_logabsdet_jacobian(::typeof(log1pexp), x::Real)
    y = log1pexp(x)
    return y, x - y
end

function ChangesOfVariables.with_logabsdet_jacobian(::typeof(logexpm1), x::Real)
    y = logexpm1(x)
    return y, x - y
end

function ChangesOfVariables.with_logabsdet_jacobian(::typeof(log1mexp), x::Real)
    y = log1mexp(x)
    return y, x - y
end

function ChangesOfVariables.with_logabsdet_jacobian(::typeof(log2mexp), x::Real)
    y = log2mexp(x)
    return y, x - y
end

function ChangesOfVariables.with_logabsdet_jacobian(::typeof(logit), x::Real)
    y = logit(x)
    y, -log(x * (1 - x))
end

function ChangesOfVariables.with_logabsdet_jacobian(::typeof(logistic), x::Real)
    y = logistic(x)
    y, log(y * (1 - y))
end

function ChangesOfVariables.with_logabsdet_jacobian(::typeof(logcosh), x::Real)
    abs_x = abs(x)
    a = - 2 * abs_x
    z = log1pexp(a)
    y = abs_x + z - IrrationalConstants.logtwo
    return y, log1mexp(a) - z
end


text:
using SampleProjectTest
using SampleProject
using Test

@testset "Trigonometry" begin
    @bench @testset "Sin" begin
        s = sin_sum(10_000_000)
        @test isapprox(s, 0; atol = 1e-6)
    end

    @bench @testset "Cos" begin
        s = cos_sum(10_000_000)
        @test isapprox(s, 1; atol = 1e-6)
    end
end

max_stars_repo_name: maxbennedich/julia-regression-analysis

max_issues_repo_name: maxbennedich/julia-regression-analysis

max_forks_repo_name: maxbennedich/julia-regression-analysis


text:
using DeformableRegistration.ImageProcessing
using DeformableRegistration.Types
using Images
using Base.Test
using MicroLogging

function checkImageProperties(img::regImage)
    # check if image has two dimensions
    @test ndims(img.data) == 2

    # check if image values are between 0 and 1
    @test float64(minfinite(img.data))>=0.0
    @test float64(maxfinite(img.data))<=1.0
    # check if colorspace is float
    @test eltype(img.data) <: Float64
end

@testset "ImageProcessing" begin
@testset "load or create image" begin
    # test loadImage
    testimage = dirname(Base.source_path()) * "/testdata/luebeck.jpg"
    img = loadImage(testimage)
    checkImageProperties(img)
    # check if spatial domain is the one that was specified
    img = loadImage(testimage)
    @test img.voxelsize == [1.0,1.0]
    @test img.shift == [0.0,0.0]

    # test createImage
    imgdata = rand(43,123)
    img = createImage(imgdata)
    checkImageProperties(img)
    @test size(imgdata) == (43,123)
end

@testset "restriction" begin
    # test restrictResolutionToLevel
    imgdata = rand(64,128)
    img = createImage(imgdata)
    checkImageProperties(img)
    restrictedImage = restrictResolutionToLevel(img,1)
    @test img.voxelsize[1] == ceil(restrictedImage.voxelsize[1] / 2)
    @test img.voxelsize[2] == ceil(restrictedImage.voxelsize[2] / 2)

    # test restrictResolutionToLevel
    imgdata = rand(65,128)
    img = createImage(imgdata)
    # check image properties for each level of the created image
    restrictedImage = restrictResolutionToLevel(img,2)
    @test img.voxelsize[1] == ceil(ceil(restrictedImage.voxelsize[1] / 2) / 2)
    @test img.voxelsize[2] == ceil(ceil(restrictedImage.voxelsize[2] / 2) / 2)

end

@testset "write image" begin
    using FileIO, QuartzImageIO

    imgdata = rand(64,128) # x \in [0, 1]
    img = createImage(imgdata)
    # load image, write image, load the written image and compare it
    testimagewrite = dirname(Base.source_path()) * "/testdata/testimage_write.jpg"
    save(testimagewrite, img.data)
    imgwritten = loadImage(testimagewrite)
    A = convert(Array{Float64,2},img.data)
    B = convert(Array{Float64,2},imgwritten.data)
    @test norm(A .- B)/length(A) < 0.0001
end
end


text:
"""
    publish(channel, message)

Post a message to a channel.
"""
publish(channel, message; client=get_global_client()) = execute(["PUBLISH", channel, message], client)

"""
    subscribe(fn::Function,
              channel,
              channels...;
              stop_fn::Function=(msg) -> false,
              err_cb::Function=(err) -> log_error(err))

Listen for messages published to the given channels in a do block. Optionally provide a stop
function `stop_fn(msg)` which gets run as a callback everytime a subscription message is received,
the subscription loop breaks if the `stop_fn` returns `true`. Optionally provide `err_cb(err)`
function which gets run on encountering an exception in the main subscription loop.

# Examples
```julia-repl
julia> channels = ["first", "second"];

julia> publisher = Client();

julia> subscriber = Client();

julia> stop_fn(msg) = msg[end] == "close subscription";  # stop the subscription loop if the message matches

julia> messages = [];

julia> @async subscribe(channels...; stop_fn=stop_fn, client=subscriber) do msg
           push!(messages, msg)
       end;  # Without @async this function will block, alternatively use Thread.@spawn

julia> wait_until_subscribed(subscriber);

julia> subscriber.is_subscribed
true

julia> subscriber.subscriptions
Set{String} with 2 elements:
  "second"
  "first"

julia> publish("first", "hello"; client=publisher);

julia> publish("second", "world"; client=publisher);

julia> println(messages)
Any[["message", "first", "hello"], ["message", "second", "world"]]  # message has the format [<message type>, <channel>, <actual message>]

julia> unsubscribe("first"; client=subscriber);

julia> wait_until_channel_unsubscribed(subscriber, "first");

julia> subscriber.subscriptions
Set{String} with 1 element:
  "second"

julia> unsubscribe(; client=subscriber);  # unsubscribe from all channels

julia> wait_until_unsubscribed(subscriber);

julia> subscriber.is_subscribed
false

julia> subscriber.subscriptions
Set{String}()
```
"""
function subscribe(fn::Function, channel, channels...; stop_fn::Function=(msg) -> false, err_cb::Function=(err) -> log_error(err), client::Client=get_global_client())
    if client.is_subscribed
        throw(RedisError("SUBERROR", "Cannot open multiple subscriptions in the same Client instance"))
    end

    @lock client.lock client.subscriptions = Set([channel, channels...])
    execute(["SUBSCRIBE", client.subscriptions...], client)
    @lock client.lock set_subscribed!(client)
    yield()
    try
        while true
            type, chnl = msg = recv(client.socket)

            if type == "message" && chnl in client.subscriptions
                fn(msg)
                stop_fn(msg) && break

            elseif type == "unsubscribe"
                if isnothing(chnl)
                    @lock client.lock client.subscriptions = Set{String}()
                elseif chnl in client.subscriptions
                    @lock client.lock delete!(client.subscriptions, chnl)
                end

                isempty(client.subscriptions) && break
            end
        end
    catch err
        err_cb(err)
    finally
        if !isempty(client.subscriptions)
            unsubscribe(client.subscriptions...; client=client)
            @lock client.lock client.subscriptions = Set{String}()
            @lock client.lock flush!(client)
        end
        @lock client.lock set_unsubscribed!(client)
    end
end

"""
    unsubscribe([channels...]) -> nothing

Unsubscribes the client from the given channels, or from all of them if none is given.
"""
unsubscribe(channels...; client=get_global_client()) = execute_without_recv(["UNSUBSCRIBE", channels...], client)

"""
    psubscribe(fn::Function,
               pattern,
               patterns...;
               stop_fn::Function=(msg) -> false,
               err_cb::Function=(err) -> log_error(err))

Listen for messages published to the given channels matching ghe given patterns in a do block.
Optionally provide a stop function `stop_fn(msg)` which gets run as a callback everytime a
subscription message is received, the subscription loop breaks if the `stop_fn` returns `true`.
Optionally provide `err_cb(err)` function which gets run on encountering an exception in the main
subscription loop.

# Examples
```julia-repl
julia> patterns = ["first*", "second*"];

julia> publisher = Client();

julia> subscriber = Client();

julia> stop_fn(msg) = msg[end] == "close subscription";  # stop the subscription loop if the message matches

julia> messages = [];

julia> @async psubscribe(patterns...; stop_fn=stop_fn, client=subscriber) do msg
           push!(messages, msg)
       end;  # Without @async this function will block, alternatively use Thread.@spawn

julia> wait_until_subscribed(subscriber);

julia> subscriber.is_subscribed
true

julia> subscriber.psubscriptions
Set{String} with 2 elements:
  "first*"
  "second*"

julia> publish("first_pattern", "hello"; client=publisher);

julia> publish("second_pattern", "world"; client=publisher);

julia> println(messages)
Any[["pmessage", "first*", "first_pattern", "hello"], ["pmessage", "second*", "second_pattern", "world"]]  # message has the format [<message type>, <pattern>, <channel>, <actual message>]

julia> punsubscribe("first*"; client=subscriber);

julia> wait_until_pattern_unsubscribed(subscriber, "first*");

julia> subscriber.psubscriptions
Set{String} with 1 element:
  "second*"

julia> punsubscribe(; client=subscriber);  # unsubscribe from all patterns

julia> wait_until_unsubscribed(subscriber);

julia> subscriber.is_subscribed
false

julia> subscriber.psubscriptions
Set{String}()
```
"""
function psubscribe(fn::Function, pattern, patterns...; stop_fn::Function=(msg) -> false, err_cb::Function=(err) -> log_error(err), client::Client=get_global_client())
    if client.is_subscribed
        throw(RedisError("SUBERROR", "Cannot open multiple subscriptions in the same Client instance"))
    end

    @lock client.lock client.psubscriptions = Set([pattern, patterns...])
    execute(["PSUBSCRIBE", client.psubscriptions...], client)
    @lock client.lock set_subscribed!(client)

    try
        while true
            type, pttrn = msg = recv(client.socket)

            if type == "pmessage" && pttrn in client.psubscriptions
                fn(msg)
                stop_fn(msg) && break

            elseif type == "punsubscribe"
                if isnothing(pttrn)
                    @lock client.lock client.psubscriptions = Set{String}()
                elseif pttrn in client.psubscriptions
                    @lock client.lock delete!(client.psubscriptions, pttrn)
                end

                isempty(client.psubscriptions) && break
            end
        end
    catch err
        err_cb(err)
    finally
        if !isempty(client.psubscriptions)
            punsubscribe(client.psubscriptions...; client=client)
            @lock client.lock client.psubscriptions = Set{String}()
            @lock client.lock flush!(client)
        end
        @lock client.lock set_unsubscribed!(client)
    end
end

"""
    unsubscribe([channels...]) -> nothing

Unsubscribes the client from the given patterns, or from all of them if none is given.
"""
punsubscribe(patterns...; client=get_global_client()) = execute_without_recv(["PUNSUBSCRIBE", patterns...], client)


text:
const Target = AddressMap{<:Union{Empty, Select}}

function (t::Target)(am::A) where A <: AddressMap
    d = Dict()
    for (k, v) in shallow_iterator(am)
        if haskey(t, k)
            d[k] = get_value(v)
        else
            t(d, get_sub(v))
        end
    end
    d
end

@inline function Base.in(addr, tg::Target)
    get_sub(tg, addr) === SelectAll()
end

@inline Base.getindex(tg::AddressMap{S}, addr) where S <: Select = get_sub(tg, addr)
@inline get_sub(tg::Target, addr) = get_sub(tg, addr)
@inline Base.merge(::SelectAll, ::Target) = SelectAll()
@inline Base.merge(::Target, ::SelectAll) = SelectAll()
@inline Base.merge(::SelectAll, ::SelectAll) = SelectAll()
@inline Base.merge(::SelectAll, ::Empty) = SelectAll()
@inline Base.merge(::Empty, ::SelectAll) = SelectAll()

const DynamicTarget = DynamicMap{Select}

@inline haskey(dm::DynamicTarget, addr::A) where A <: Address = haskey(dm.tree, addr) && !isempty(get_sub(dm, addr))

Base.merge!(dt::DynamicTarget, s::SelectAll) = dt

@inline Base.push!(tg::DynamicTarget, addr) = set_sub!(tg, addr, SelectAll())
@inline function Base.push!(tg::DynamicTarget, addr::Tuple{}) end
@inline function Base.push!(tg::DynamicTarget, addr::Tuple{T}) where T
    Base.push!(tg, addr[1])
end
@inline function Base.push!(tg::DynamicTarget, addr::Tuple) where T
    hd, tl = addr[1], addr[2 : end]
    sub = get_sub(tg, hd)
    if sub isa DynamicTarget
        push!(sub, tl)
    else
        new = target(tl)
        merge!(new, sub)
        set_sub!(tg, hd, new)
    end
end

function target(v::Vector{T}) where T <: Tuple
    tg = DynamicTarget()
    for k in v
        push!(tg, k)
    end
    tg
end

function target(k::A) where A <: Address
    tg = DynamicTarget()
    push!(tg, k)
    tg
end

function target(k::Tuple)
    tg = DynamicTarget()
    push!(tg, k)
    tg
end

function flatten(t::Target)
    addrs, chd, _ = collect(t)
    addrs, Float64[], chd
end


text:
# # N queens

using Convex, GLPK, LinearAlgebra, SparseArrays, Test
aux(str) = joinpath(@__DIR__, "aux_files", str) # path to auxiliary files
include(aux("antidiag.jl"))

n = 8
# We encode the locations of the queens with a matrix of binary random variables.
x = Variable((n, n), :Bin)

# Now we impose the constraints: at most one queen on any anti-diagonal, at most one queen on any diagonal, and we must have exactly one queen per row and per column.
## At most one queen on any anti-diagonal
constr = Constraint[sum(antidiag(x, k)) <= 1 for k = -n+2:n-2]
## At most one queen on any diagonal
constr += Constraint[sum(diag(x, k)) <= 1 for k = -n+2:n-2]
## Exactly one queen per row and one queen per column
constr += Constraint[sum(x, dims=1) == 1, sum(x, dims=2) == 1]
p = satisfy(constr)
solve!(p, GLPK.Optimizer)

# Let us test the results:
for k = -n+2:n-2
	@test evaluate(sum(antidiag(x, k))) <= 1
	@test evaluate(sum(diag(x, k))) <= 1
end
@test all(evaluate(sum(x, dims=1)) .≈ 1)
@test all(evaluate(sum(x, dims=2)) .≈ 1)


text:
if length(ARGS) < 1
    println("Please provide input file -- bed format")
    exit(0)
else
    bedf = ARGS[1]
end

f = open(bedf, "r")
n = countlines(f)
seekstart(f)

for i = 1:n
    line = readline(f)
    if !(line[1:5]=="track")

        s = split(line, "\t")
        p_start = parse(Int64, s[2])
        p_end   = parse(Int64, s[3])

        println(string(p_start , "\t" , p_end , "\t" ,(p_end - p_start ),"\t", s[4]))

    end
end


text:
#src # Shock capturing with flux differencing and stage limiter

# This tutorial contains a short summary of the idea of shock capturing for DGSEM with flux differencing
# and its implementation in [Trixi.jl](https://github.com/trixi-framework/Trixi.jl).
# In the second part, an implementation of a positivity preserving limiter is added to the simulation.

# # Shock capturing with flux differencing

# The following rough explanation is on a very basic level. More information about an entropy stable
# shock-capturing strategy for DGSEM discretizations of advection dominated problems, such as the
# compressible Euler equations or the compressible Navier-Stokes equations, can be found in
# [Hennemann et al. (2021)](https://doi.org/10.1016/j.jcp.2020.109935). In
# [Rueda-Ramírez et al. (2021)](https://doi.org/10.1016/j.jcp.2021.110580) you find the extension to
# the systems with non-conservative terms, such as the compressible magnetohydrodynamics (MHD) equations.

# The strategy for a shock-capturing method presented by Hennemann et al. is based on a hybrid blending
# of a high-order DG method with a low-order variant. The low-order subcell finite volume (FV) method is created
# directly with the Legendre-Gauss-Lobatto (LGL) nodes already used for the high-order DGSEM.
# Then, the final method is a convex combination with regulating indicator $\alpha$ of these two methods.

# Since the surface integral is equal for both the DG and the subcell FV method, only the volume integral divides
# between the two methods.

# This strategy for the volume integral is implemented in Trixi under the name of
# [`VolumeIntegralShockCapturingHG`](@ref) with the three parameters of the indicator and the volume fluxes for
# the DG and the subcell FV method.

# Note, that the DG method is based on the flux differencing formulation. Hence, you have to use a
# two-point flux, such as [`flux_ranocha`](@ref), [`flux_shima_etal`](@ref), [`flux_chandrashekar`](@ref) or [`flux_kennedy_gruber`](@ref),
# for the DG volume flux. We would recommend to use the entropy conserving flux `flux_ranocha` by
# [Ranocha (2018)](https://cuvillier.de/en/shop/publications/7743) for the compressible Euler equations.
# ````julia
# volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
#                                                  volume_flux_dg=volume_flux_dg,
#                                                  volume_flux_fv=volume_flux_fv)
# ````

# We now focus on a choice of the shock capturing indicator `indicator_sc`.
# A possible indicator is $\alpha_{HG}$ presented by Hennemann et al. (p.10), which depends on the
# current approximation with modal coefficients $\{m_j\}_{j=0}^N$ of a given `variable`.

# The indicator is calculated for every DG element by itself. First, we calculate a smooth $\alpha$ by
# ```math
# \alpha = \frac{1}{1+\exp(-\frac{-s}{\mathbb{T}}(\mathbb{E}-\mathbb{T}))}
# ```
# with the total energy $\mathbb{E}=\max\big(\frac{m_N^2}{\sum_{j=0}^N m_j^2}, \frac{m_{N-1}^2}{\sum_{j=0}^{N-1} m_j^2}\big)$,
# threshold $\mathbb{T}= 0.5 * 10^{-1.8*(N+1)^{1/4}}$ and parameter $s=ln\big(\frac{1-0.0001}{0.0001}\big)\approx 9.21024$.

# For computational efficiency, $\alpha_{min}$ is introduced und used for
# ```math
# \tilde{\alpha} = \begin{cases}
# 0, & \text{if } \alpha<\alpha_{min}\\
# \alpha, & \text{if } \alpha_{min}\leq \alpha \leq 1- \alpha_{min}\\
# 1, & \text{if } 1-\alpha_{min}<\alpha.
# \end{cases}
# ```

# Moreover, the parameter $\alpha_{max}$ sets a maximal value for $\alpha$ by
# ```math
# \alpha = \min\{\tilde{\alpha}, \alpha_{max}\}.
# ```
# This allows to control the maximal dissipation.

# To remove numerical artifact the final indicator is smoothed with all the neighboring elements'
# indicators. This is activated with `alpha_smooth=true`.
# ```math
# \alpha_{HG} = \max_E \{ \alpha, 0.5 * \alpha_E\},
# ```
# where $E$ are all elements sharing a face with the current element.

# Furthermore, you can specify the variable used for the calculation. For instance you can choose
# `density`, `pressure` or both with `density_pressure` for the compressible Euler equations.
# For every equation there is also the option to use the first conservation variable with `first`.

# This indicator is implemented in Trixi and called [`IndicatorHennemannGassner`](@ref) with the parameters
# `equations`, `basis`, `alpha_max`, `alpha_min`, `alpha_smooth` and `variable`.
# ````julia
# indicator_sc = IndicatorHennemannGassner(equations, basis,
#                                          alpha_max=0.5,
#                                          alpha_min=0.001,
#                                          alpha_smooth=true,
#                                          variable=variable)
# ````

# # Positivity preserving limiter

# Some numerical solutions are physically meaningless, for instance negative values of pressure
# or density for the compressible Euler equations. This often results in crashed simulations since
# the calculation of numerical fluxes or stable time steps uses mathematical operations like roots or
# logarithms. One option to avoid these cases are a-posteriori positivity preserving limiters.
# Trixi provides the fully-discrete positivity-preserving limiter of
# [Zhang, Shu (2011)](https://doi.org/10.1098/rspa.2011.0153).

# It works the following way. For every passed (scalar) variable and for every DG element we calculate
# the minimal value $value_{min}$. If this value falls below the given threshold $\varepsilon$,
# the approximation is slightly adapted such that the minimal value of the relevant variable lies
# now above the threshold.
# ```math
# \underline{u}^{new} = \theta * \underline{u} + (1-\theta) * u_{mean}
# ```
# where $\underline{u}$ are the collected pointwise evaluation coefficients in element $e$ and
# $u_{mean}$ the integral mean of the quantity in $e$. The new coefficients are a convex combination
# of these two values with factor
# ```math
# \theta = \frac{value_{mean} - \varepsilon}{value_{mean} - value_{min}},
# ```
# where $value_{mean}$ is the relevant variable evaluated for the mean value $u_{mean}$.

# The adapted approximation keeps the exact same mean value, but the relevant variable is now greater
# or equal the threshold $\varepsilon$ at every node in every element.

# We specify the variables the way we did before for the shock capturing variables. For the
# compressible Euler equations `density`, `pressure` or the combined variable `density_pressure`
# are a reasonable choice.

# You can implement the limiter in Trixi using [`PositivityPreservingLimiterZhangShu`](@ref) with parameters
# `threshold` and `variables`.
# ````julia
# stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=thresholds,
#                                                      variables=variables)
# ````
# Then, the limiter is added to the time integration method in the `solve` function. For instance, like
# ````julia
# CarpenterKennedy2N54(stage_limiter!, williamson_condition=false)
# ````
# or
# ````julia
# SSPRK43(stage_limiter!).
# ````

# # Simulation with shock capturing and positivity preserving

# Now, we can run a simulation using the described methods of shock capturing and positivity
# preserving limiters. We want to give an example for the 2D compressible Euler equations.
using OrdinaryDiffEq, Trixi

equations = CompressibleEulerEquations2D(1.4)

# As our initial condition we use the Sedov blast wave setup.
function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D)
  ## Set up polar coordinates
  inicenter = SVector(0.0, 0.0)
  x_norm = x[1] - inicenter[1]
  y_norm = x[2] - inicenter[2]
  r = sqrt(x_norm^2 + y_norm^2)

  r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
  ## r0 = 0.5 # = more reasonable setup
  E = 1.0
  p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2)
  p0_outer = 1.0e-5 # = true Sedov setup
  ## p0_outer = 1.0e-3 # = more reasonable setup

  ## Calculate primitive variables
  rho = 1.0
  v1  = 0.0
  v2  = 0.0
  p   = r > r0 ? p0_outer : p0_inner

  return prim2cons(SVector(rho, v1, v2, p), equations)
end
initial_condition = initial_condition_sedov_blast_wave
#-
basis = LobattoLegendreBasis(3)

# We set the numerical fluxes and divide between the surface flux and the two volume fluxes for the DG
# and FV method. Here, we are using [`flux_lax_friedrichs`](@ref) and [`flux_ranocha`](@ref).
surface_flux = flux_lax_friedrichs
volume_flux  = flux_ranocha

# Now, we specify the shock capturing indicator $\alpha$.

# We implement the described indicator of Hennemann, Gassner as explained above with parameters
# `equations`, `basis`, `alpha_max`, `alpha_min`, `alpha_smooth` and `variable`.
# Since density and pressure are the critical variables in this example, we use
# `density_pressure = density * pressure = rho * p` as indicator variable.
indicator_sc = IndicatorHennemannGassner(equations, basis,
                                         alpha_max=0.5,
                                         alpha_min=0.001,
                                         alpha_smooth=true,
                                         variable=density_pressure)

# Now, we can use the defined fluxes and the indicator to implement the volume integral using shock
# capturing.
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)

# We finalize the discretization by implementing Trixi's `solver`, `mesh`, `semi` and `ode`,
# while `solver` now has the extra parameter `volume_integral`.
solver = DGSEM(basis, surface_flux, volume_integral)

coordinates_min = (-2.0, -2.0)
coordinates_max = ( 2.0,  2.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=6,
                n_cells_max=10_000)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan);

# We add some callbacks to get an solution analysis and use a CFL-based time step size calculation.
analysis_callback = AnalysisCallback(semi, interval=100)

stepsize_callback = StepsizeCallback(cfl=0.8)

callbacks = CallbackSet(analysis_callback, stepsize_callback);

# We now run the simulation using the positivity preserving limiter of Zhang and Shu for the variables
# density and pressure.
stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(5.0e-6, 5.0e-6),
                                                     variables=(Trixi.density, pressure))

sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);

using Plots
plot(sol)


text:
# TO DO:
# As soon as someone sends me a file with waveform traces that have the
# attributes event_id, magnitude_id, focal_mechanism_id, I can trivially
# output SeisEvents. No one has done that yet.

@doc """
    (H,R) = asdf_rqml(fpat::String)

Read QuakeXML (qml) from ASDF file(s) matching file string pattern `fpat`. Returns:
* `H`, Array{SeisHdr,1}
* `R`, Array{SeisSrc,1}

""" asdf_rqml
function asdf_rqml(io::HDF5File)
  EvCat = Array{SeisHdr,1}()
  EvSrc = Array{SeisSrc,1}()

  if has(io, "QuakeML")
    qml = parse_string(String(UInt8.(read(io["QuakeML"]))))
    event_xml!(EvCat, EvSrc, qml)
    free(qml)
  end
  return EvCat, EvSrc
end

function asdf_rqml(fpat::String)
  files = safe_isfile(fpat) ? [fpat] : ls(fpat)
  EvCat = Array{SeisHdr,1}()
  EvSrc = Array{SeisSrc,1}()
  for file in files
    io = h5open(file, "r")
    if has(io, "QuakeML")
      qml = parse_string(String(UInt8.(read(io["QuakeML"]))))
      event_xml!(EvCat, EvSrc, qml)
      free(qml)
    end
    close(io)
  end
  return(EvCat, EvSrc)
end

function asdf_wqml!(hdf::HDF5File, HDR::Array{SeisHdr,1}, SRC::Array{SeisSrc,1},
  ovr::Bool, v::Integer)
  hq = Int8[]
  io = IOBuffer()
  H0 = SeisHdr[]
  R0 = SeisSrc[]
  if has(hdf, "QuakeML")
    if ovr
      (v > 0) && println("ovewriting QuakeML...")
      new_qml!(io)
    else
      hq = read(hdf["QuakeML"])

      # check whether we can append directly
      nq = length(hq)
      si = nq-29
      test_str = String(UInt8.(hq[si:nq]))
      if test_str == "</eventParameters>\n</quakeml>\n"
        # behavior for SeisIO-compatible QuakeML
        deleteat!(hq, si:nq)
        append!(hq, ones(Int8, 4).*Int8(32))
      else
        # behavior for other QuakeML
        qml = parse_string(String(UInt8.(hq)))
        event_xml!(H0, R0, qml)
        free(qml)
        new_qml!(io)
        hq = Int8[]
        GC.gc()
      end
    end

    # delete hdf["QuakeML"]. I fear one can't resize HDF5 arrays in-place.
    # One certainly can't in the Julia HDF5 interface...
    o_delete(hdf, "QuakeML")
  else
    new_qml!(io)
  end
  if isempty(H0) && isempty(R0)
    write_qml!(io, HDR, SRC, v)
  else
    append!(H0, HDR)
    append!(R0, SRC)
    write_qml!(io, H0, R0, v)
  end
  buf = vcat(hq, take!(io))
  hdf["QuakeML"] = buf
  close(io)
  return nothing
end

@doc """
    asdf_wqml(fname, SHDR::Array{SeisHdr,1}, SSRC::Array{SeisSrc,1}[, KWs])
    asdf_wqml(fname, H::SeisHdr, R::SeisSrc[, KWs])

Write QuakeXML (qml) to "QuakeML/" dataset in ASDF file `fname` from `SHDR` and
`SSRC`.

    asdf_wqml(fname, evt::SeisEvent[, KWs])
    asdf_wqml(fname, evt::Array{SeisEvent,1}[, KWs])

As above, for the `:hdr` and `:source` fields of `evt`.

|KW     | Type      | Default   | Meaning                                     |
|:---   |:---       |:---       |:---                                         |
| ovr   | Bool      | false     | overwrite QML in existing ASDF file? [^1]   |
| v     | Integer   | 0         | verbosity                                   |

[^1] By default, data are appended to the existing contents of "QuakeML/".

!!! warning

    To write data from `R ∈ SSRC`, it must be true that `R.eid == H.id` for some `H ∈ SHDR`.

See also: `write_qml`
""" asdf_wqml
function asdf_wqml(hdf_out::String, HDR::Array{SeisHdr,1}, SRC::Array{SeisSrc,1};
  ovr::Bool=false,
  v::Integer=0)
  if isfile(hdf_out)
    hdf = h5open(hdf_out, "r+")
    fmt = read(attrs(hdf)["file_format"])
    (fmt == "ASDF") || (close(hdf); error("invalid ASDF file!"))
  else
    hdf = h5open(hdf_out, "cw")
    attrs(hdf)["file_format"] = "ASDF"
    attrs(hdf)["file_format_version"] = "1.0.2"
  end
  asdf_wqml!(hdf, HDR, SRC, ovr, v)
  close(hdf)

  # logging
  opts = string(", ovr=\"", ovr, "\", v=", v)
  for h in HDR
    fwrite_note_quake!(h, "asdf_wqml", hdf_out, opts)
  end
  for r in SRC
    fwrite_note_quake!(r, "asdf_wqml", hdf_out, opts)
  end
  return nothing
end
asdf_wqml(hdf_out::String, H::SeisHdr, R::SeisSrc; ovr::Bool=false, v::Integer=0) = asdf_wqml(hdf_out, [H], [R], ovr=ovr, v=v)
asdf_wqml(hdf_out::String, W::SeisEvent; ovr::Bool=false, v::Integer=0) = asdf_wqml(hdf_out, [W.hdr], [W.source], ovr=ovr, v=v)

function asdf_wqml(hdf_out::String, events::Array{SeisEvent,1};
  ovr::Bool=false,
  v::Integer=0)

  N = length(events)
  H = Array{SeisHdr, 1}(undef, N)
  R = Array{SeisSrc, 1}(undef, N)
  for i in 1:N
    H[i] = events[i].hdr
    R[i] = events[i].source
  end
  asdf_wqml(hdf_out, H, R, ovr=ovr, v=v)
  return nothing
end


text:
"""
    ET0calc_read(fileETo)
Read output ETo output file of EToCalc.exe
See http://www.fao.org/land-water/databases-and-software/eto-calculator/en/

**Input**
* fileETo: full file name (ETo extension)

**Output**
* datafame containing time and reference evaporanspitation

**Example**
```
et_calc = ET0calc_read(joinpath(dirname(@__DIR__),"test","input","ET0calc.ETo"));
```
"""
function ET0calc_read(filein::String)::DataFrame
    # read header
    parsehead(fi) = fi |> readline |> x-> split(x,':') |> x -> Meta.parse(x[1])
    resol,day,month,year = Dates.Day(1),1,1,1901
    open(filein,"r") do fid
        foreach(x->readline(fid),1:1)
        record_type = parsehead(fid)
        if record_type == 2
            resol = Dates.Day(10);
        elseif record_type == 3
            resol = Dates.Month(1);
        end
        day,month,year = parsehead(fid),parsehead(fid),parsehead(fid);
    end
    # read data
    temp = readdlm(filein,skipstart=8);
    return DataFrame(datetime=collect(DateTime(year,month,day):resol:DateTime(year,month,day)+(size(temp,1)-1)*resol),
                    ET0 = temp[:,1])
end

"""
    ET0calc_dsc(fileout;name,lat,lon,H,from_date,to_date,wind_height,psych,
                miss_wind,miss_rad,allow_dev,temp_range,humi_range,wind_range,
                Angstrom)
Write ET0calc input settings (dsc) file
See http://www.fao.org/land-water/databases-and-software/eto-calculator/en/

**Input**
* `fileout`: output file name
* `name`: ["Site", "Country"] vector
* `lat,lon,H`: Latitude, longitude, height (deg,deg,m)
* `from_date,to_date`: starting to end date (DateTime)
* `wind_height`: heigth of wind measurement (m)
* `psych`: coefficient psychrometer (depending on type of ventilation)
* `miss_wind`: default value for missing wind data (m/s)
* `alow_dev`: Allowable deviation (%) from theoretical maximum radiation
* `temp_range`: allowed ragne [minimum, maximum] (deg)
* `humi_range`: allowed humidity range [minimum, maximum] (%)
* `wind_range`: allowed wind speed range [minimum, maximum] (m/s)
* `Angstrom`: Adjustment for station elevation (Angstrom formula)

**Example**
```
out_file = joinpath(dirname(@__DIR__),"test","output","ET0calc.DSC");
latitude = -34.873305; # degrees
longitude = -58.139995; # degrees
altitude = 27.95; # meters
utc_time = Dates.Hour(3); # will be SUBTRACTED
wind_height = 2.; # meters
ET0calc_dsc(out_file,name = ["LaPlata","Argentina"],
            lat = latitude,lon=longitude,H=altitude,
            from_date=DateTime(2010,1,1),to_date=DateTime(2011,2,3),
            wind_height=wind_height,allow_dev=10,
            miss_wind=3.,miss_rad=0.2,psych=0.009,
            temp_rang=[-10.,40.],humi_range=[1.0,99.],wind_range=[1.5,10.7]);
```
"""
function ET0calc_dsc(fileout; name::Vector{String}=["Site","Country"],
                    lat::Float64=0.0,lon::Float64=0.0,H::Float64=0.0,
                    from_date=DateTime(2010,1,1),to_date=DateTime(2010,1,2),
                    wind_height::Float64=2.0,
                    psych::Float64=0.000800,miss_wind::Float64=2.0,
                    miss_rad::Float64=0.19,allow_dev::Int=5,
                    temp_rang=[-20.,50.],humi_range=[0.0,100.],wind_range=[0.,10],
                    Angstrom=[0.25,0.5,1])::Nothing
    open(fileout,"w") do fid
        @printf(fid,"%s\r\n%s\r\n",name[1],name[2]);
        @printf(fid,"%10.2f                   : Latitude (degrees)\r\n",lat);
        @printf(fid,"%10.2f                   : Longitude (degrees)\r\n",lon);
        @printf(fid,"%7i                      : Altitude (meters above sea level)\r\n",H);
        @printf(fid,"      1                      : daily data\r\n");
        @printf(fid,"%4i%4i%6i               : From Date (Day - Month - Year)\r\n",
                Dates.day(from_date),Dates.month(from_date),Dates.year(from_date));
        @printf(fid,"%4i%4i%6i               : To Date (Day - Month - Year)\r\n",
                Dates.day(to_date),Dates.month(to_date),Dates.year(to_date));
        @printf(fid,"%14.6f               : coefficient psychrometer (depending on type of ventilation)\r\n",psych);
        @printf(fid,"%9.1f                    : height wind speed measurement (meter)\r\n",wind_height);
        @printf(fid,"%9.1f                    : Missing air humidity data: difference between Tdew and Tmin (°C)\r\n",0.0)
        @printf(fid,"%10.2f                   : Missing radiation data: Coefficient of Hargreaves Equation (MJ/m2.day)\r\n",miss_rad);
        @printf(fid,"%9.1f                    : Missing wind speed: estimate for U2 (m/sec)\r\n",miss_wind);
        @printf(fid,"%10.2f                   : Coefficient: a Angstrom formula\r\n",Angstrom[1]);
        @printf(fid,"%10.2f                   : Coefficient: b Angstrom formula\r\n",Angstrom[2]);
        @printf(fid,"%7i                      : Option: Adjustment for station elevation (Angstrom formula)\r\n",Angstrom[3]);
        @printf(fid,"%7i                      : Climatic data range: lower limit for Minimal temperature (°C)\r\n",temp_rang[1]);
        @printf(fid,"%7i                      : Climatic data range: upper limit for Maximal temperature (°C)\r\n",temp_rang[2]);
        @printf(fid,"%7i                      : Climatic data range: lower limit for Minimal Relative Humidity (%%)\r\n",humi_range[1]);
        @printf(fid,"%7i                      : Climatic data range: upper limit for Maximal wind speed (m/sec)\r\n",wind_range[2]);
        @printf(fid,"%7i                      : Climatic data range: Allowable deviation (%%) from theoretical maximum radiation\r\n",allow_dev);
    end
end

"""
    ET0calc_dta(datawrite,filewrite;radiation_units)
Write ET0calc input data (dta) file
See http://www.fao.org/land-water/databases-and-software/eto-calculator/en/
See `test/input/ET0calc_codes.txt` for output parameter codes

**Input**
* `datawrite`: dataframe containing datetime,Tmin,Tmax(|Tmean),RHmin,RHmax(|RHmean),ux,Rn|Rs (net/solar or shortwave radiation) columns
* `filewrite`: output file name
* `radiation_units`: radiation units ("W.m^-2" or "MJ.m^-2/day")

**Example**
```
out_file = joinpath(dirname(@__DIR__),"test","output","ET0calc_out.DTA");
datause = DataFrame(
        datetime=collect(DateTime(2010,1,1):Dates.Day(1):DateTime(2010,1,10)),
        Tmin = zeros(10).-5.0, Tmax = zeros(10).+40.0,
        RHmin = zeros(10).+1.0, RHmax = zeros(10).+99.0,
        ux = zeros(10).+1.0, Rn = ones(10).+30.0)
# Optional: get minim/max values from hourly samples using ResampleAndFit Pkg
# min_data = ResampleAndFit.aggregate2(datause,resol=Dates.Day(1),fce=minimum);
# max_data = ResampleAndFit.aggregate2(datause,resol=Dates.Day(1),fce=maximum);
# mean_data = ResampleAndFit.aggregate2(datause,resol=Dates.Day(1),fce=mean);
ET0calc_dta(datause,out_file,radiation_units="W.m^-2");
```
"""
function ET0calc_dta(datawrite,filewrite;radiation_units="W.m^-2")::Nothing
    open(filewrite,"w") do fid
        ET0calc_dta_head(fid,names(datawrite),radiation_units);
        for i in 1:size(datawrite,1)
            @printf(fid,"%10i",i)
            for j in names(datawrite)
                if j != :datetime
                    @printf(fid,"%10.2f",
                        (ismissing(datawrite[j][i]) || isnan(datawrite[j][i])) ?
                        -999. : datawrite[j][i])
                end
            end
            @printf(fid,"    -999.0\r\n");
        end
    end
end

"""
Auxiliary function to write header for ET0calc DTA file
"""
function ET0calc_dta_head(fid::IOStream,col_names,radiation_units)::Nothing
    @printf(fid,"        Nr");
    for j in col_names
        if j == :Tmin
            @printf(fid,"       103")
        elseif j == :Tmax
            @printf(fid,"       101")
        elseif j == :Tmean
            @printf(fid,"       102")
        elseif j == :RHmax
            @printf(fid,"       201")
        elseif j == :RHmin
            @printf(fid,"       203")
        elseif j == :RHmean
            @printf(fid,"       202")
        elseif j == :ux
            @printf(fid,"       301")
        elseif j == :Rn
            if radiation_units == "W.m^-2"
                @printf(fid,"       432")
            else
                @printf(fid,"       431")
            end
        elseif j == :Rs
            if radiation_units == "W.m^-2"
                @printf(fid,"       422")
            else
                @printf(fid,"       421")
            end
        end
    end
    @printf(fid,"       501\r\n")
end


text:
module Poker

include("cards.jl")
include("odds.jl")
include("game.jl")

export +, -, >, <, fresh_deck, deal, cards, isless, isflush, to_string,
    calculate_odds, get_preflop_odds, bet, get_return, get_returns,
    get_returns_table

end


text:
using .Unitful

struct QuantityConstructor{D,U} end
QuantityConstructor{D,U}(x::T) where {T,D,U} = Quantity{T,D,U}(x)

constructor_of(::Type{Quantity{T,D,U}}) where {T,D,U} = QuantityConstructor{D,U}


text:
using Pkg
Pkg.activate("Search")
using Search
using Distributions
using LinearAlgebra
using Plots
initDist = MvNormal(Array{Float64}([1000, 1000, 100, 100]))

searchDist = MvNormal((Array{Float64}([1000,1000])))

searchSchedule = [0,1];
function update!(targets,Δt)
    targets[:,1:2] = targets[:,1:2] .+ targets[:,3:4] .* Δt
end
Search.initializeScenario(initDist,searchDist,searchSchedule,Search.linearTargetModel!)
N = 10;
targets = Array{Float64}(undef,(N,4))
for n = 1:N
    targets[n,:] = rand(initDist)
end
lim = 10000
xrange = (-lim,lim);
yrange = (-lim,lim);
scatter(targets[:,1],targets[:,2];xlim=xrange,ylim=yrange)

function initTargets(initDist::MvNormal;N=100)
    targets = Array{Float64}(undef,(N,4))
    for n = 1:N
        targets[n,:] = rand(initDist)
    end
    targets = [targets zeros(N,1)]
    return targets
end

function isFound(x1,x2,threshold)
    norm(x2.-x1)<=threshold ? true : false
end

Δt = 1/1;
targets = initTargets(initDist;N);
scatter(targets[:,1],targets[:,2];xlim=xrange,ylim=yrange)
for t = 1:Δt:50
    update!(targets,Δt)
    p = scatter(targets[:,1],targets[:,2];xlim=xrange,ylim=yrange)
    for i in 1:N
        if isFound(targets[i,1:2],[0 0],10000)
            targets[i,:] .= NaN;
            targets[i,end] = 1;
        end
    end
    if all(i -> i==1,targets[:,5])
        println("All targets found")
        break
    end
    p = scatter(targets[:,1],targets[:,2];xlim=xrange,ylim=yrange)
    display(p)
    sleep(.1)
end

dims(x) = length(size(x))

function gridDist(d::MvNormal;δ = .1)
    Δ = 2*sqrt(det(cov(d)))
    i = -Δ:δ:Δ;
    j = -Δ:δ:Δ;
    z = [pdf(d,[i,j]) for i in i, j in j]
    return i,j,z
end

function addGrid(i,j,A::AbstractArray,k,l,B::AbstractArray,Δ)
    x = minimum([i...,k...]):Δ :maximum([i...,k...]);
    y = minimum([j...,l...]):Δ :maximum([j...,l...]);
    _,_,pad_A = padGrid(x,y,i,j,A,Δ)
    _,_,pad_B = padGrid(x,y,k,l,B,Δ)
    G = pad_A .+ pad_B
    return x,y,G
end

function multGrid(i,j,A::AbstractArray,k,l,B::AbstractArray,Δ)
    x,y,overlap_A = overlapGrid(k,l,i,j,A,Δ)
    _,_,overlap_B = overlapGrid(i,j,k,l,B,Δ)
    G = overlap_A .* overlap_B
    return x,y,G
end

function padGrid(x,y,i,j,A,Δ)
    a = abs(minimum(i) - minimum(x))/Δ
    b = length(y)
    pad_upper = zeros(convert.(Int,[a,b])...)
    println(typeof(a))
    println(typeof(b))
    a = abs(maximum(i) - maximum(x))/Δ
    pad_lower = zeros(convert.(Int,[a,b])...)
    println(typeof(a))
    println(typeof(b))
    a = length(i)
    b = abs(minimum(j) - minimum(y))/Δ
    pad_left = zeros(convert.(Int,[a,b])...)
    println(typeof(a))
    println(typeof(b))
    b = abs(maximum(j) - maximum(y))/Δ
    pad_right = zeros(convert.(Int,[a,b])...)
    println(typeof(a))
    println(typeof(b))
    pad_A = [pad_upper;pad_left A pad_right;pad_lower]
    return x,y,pad_A
end

function overlapGrid(k,l,i,j,A,Δ)
    x = maximum( [minimum(k),minimum(i)] ):Δ:minimum( [maximum(k),maximum(i)] )
    y = maximum( [minimum(l),minimum(j)] ):Δ:minimum( [maximum(l),maximum(j)] )
    xi = indexin([minimum(x),maximum(x)],i)
    yi = indexin([minimum(y),maximum(y)],j)
    overlap_A = A[xi[1]:xi[2],yi[1]:yi[2]]
    return x,y,overlap_A
end


text:
using BinDeps

@BinDeps.setup

## Start defining dependency
const libhypre = library_dependency("libhypre", aliases=["libhypre", "libHYPRE"])
const hyprever = "2.11.2"
const hyprefilebase = "hypre-$(hyprever)"

provides(Sources,
         URI("https://computation.llnl.gov/projects/hypre-scalable-linear-solvers-multigrid-methods/download/$(hyprefilebase).tar.gz"),
         libhypre)

## Define paths used later.
const prefix = BinDeps.usrdir(libhypre)
const srcrootdir = BinDeps.srcdir(libhypre)
const srcdir = joinpath(srcrootdir, hyprefilebase, "src")
const cmakebuilddir = joinpath(srcdir, "cmbuild")

## Add support for basic "make flags"
if "clean" in ARGS || "clean-all" in ARGS
    dirsToRm = [BinDeps.srcdir(libhypre)]
    "clean-all" in ARGS && push!(dirsToRm, prefix, BinDeps.downloadsdir(libhypre))
    for dir in dirsToRm
        rm(dir, recursive = true, force = true)
    end
    exit(0)
end

mpicc = get(ENV, "JULIA_MPI_C_COMPILER", chomp(readstring(`which mpicc`)))
mpicxx = get(ENV, "JULIA_MPI_CXX_COMPILER", chomp(readstring(`which mpicxx`)))

## Build steps
provides(SimpleBuild,
(@build_steps begin
    GetSources(libhypre)
    CreateDirectory(prefix)
    @build_steps begin
        ChangeDirectory(cmakebuilddir)
        `cmake -DHYPRE_SHARED:BOOL=ON -DHYPRE_INSTALL_PREFIX:PATH=$prefix -DCMAKE_C_COMPILER:FILEPATH=$(mpicc) -DCMAKE_CXX_COMPILER:FILEPATH=$(mpicxx) ..`
        `cmake -L ..`
        `make`
        `make install`
    end
end),
libhypre,
os = :Unix)

## Install
@BinDeps.install Dict(:libhypre => :libhypre)


text:
# ## InflationFixedExclusionCPI.jl - Función de inflación de exclusión fija de gastos básicos utilizando fórmula del IPC

"""
    InflationFixedExclusionCPI{N} <: InflationFunction

Función de inflación para computar la inflación de exclusión fija de gastos
básicos para `N` bases del IPC.

## Utilización

    function (inflfn::InflationFixedExclusionCPI)(base::VarCPIBase{T}, i::Int) where T

Define cómo opera InflationFixedExclusionCPI sobre un objeto de tipo
`VarCPIBase`, con `i` listas de exclusión de gastos básicos (posiciones en el
IPC) para las bases 2000 y 2010.

    function (inflfn::InflationFixedExclusionCPI)(cs::CountryStructure, ::CPIVarInterm)

Define cómo opera InflationFixedExclusionCPI sobre un CountryStructure

## Ejemplos

Instanciamos la función y le damos dos listas de exclusión, una por cada base en
el CountryStructure. Se exclyen los gastos básicos con números de columnas
`[25,30,54,88]` del primer `VarCPIBase` y los gastos básicos `[65,95,85]` del
segundo `VarCPIBase`.
```julia-repl
julia> InfExc = InflationFixedExclusionCPI([25,30,54,88], [65,95,85])
(::InflationFixedExclusionCPI{2}) (generic function with 6 methods)
```

Cuando `gt00` es un `VarCPIBase` para la base 2000 y su vector de exclusión es
el primero de la lista:
```julia-repl
julia> InfExc(gt00, 1)
120-element Vector{Float32}:
 1.4066696
 0.7496357
 0.5115509
 0.39703846
 ⋮
 0.23813248
 0.27743578
 0.46629906
```

Cuando `gt10` es un `VarCPIBase` para la base 2010 y su vector de exclusión es
el segundo de la lista:
```julia-repl
julia> InfExc(gt10, 2)
```

Cuando `gtdata` es un `CountryStructure`:
```julia-repl
julia> InfExc(gtdata)
231-element Vector{Float32}:
 8.727622
 8.73872
 8.972383
 9.031701
 ⋮
 4.8354387
 5.24472
 6.0075045
```
"""
Base.@kwdef struct InflationFixedExclusionCPI{N} <: InflationFunction
    # Tupla con vectores de gastos básicos a excluir en cada base (tantos vectores como bases)
    v_exc::NTuple{N,Vector{Int}}
end

# Ampliación para que reciba solo los vectores, no necesariamente como una tupla.
InflationFixedExclusionCPI(v_exc...) = InflationFixedExclusionCPI(v_exc)

# Extender el método de nombre y de tag
measure_name(inflfn::InflationFixedExclusionCPI) = "Exclusión fija de gastos básicos IPC " * string(map(length, inflfn.v_exc))

# Método para obtener parámetros
params(inflfn::InflationFixedExclusionCPI) = inflfn.v_exc

# Cómputo del resumen intermensual utilizando la lista de exclusión i
function (inflfn::InflationFixedExclusionCPI)(base::VarCPIBase{T}, i::Int) where T
    # Elección del vector de exclusión a utlizar dependiendo de que base se está tomando
    # 1 -> base 2000, 2-> base 2010 (se puede generalizar para más bases)
    exc = inflfn.v_exc[i]
    # Capitalizar los índices de precios a partir del objeto base::VarCPIBase
    base_ipc= capitalize(base.v, base.baseindex)
    # creación de una copia de la lista original de pesos desde base.w
    w_exc = copy(base.w)
    # Asignación de peso cero a los gastos básicos de la lista de exclusión (exc = inflfn.v_exc[i])
    # (j itera sobre los elementos de la lista de exclusión)
        for j in exc w_exc[j] = 0 end
    # Renormalización de pesos
    w_exc = w_exc / sum(w_exc)
    # Obtener Ipc con exclusión
    cpi_exc = base_ipc*w_exc
    # Obtener variación intermensual
    varm_cpi_exc =  varinterm(cpi_exc)
end

# Redefinición del método para obtener variaciones intermensuales del CountryStructure
function (inflfn::InflationFixedExclusionCPI)(cs::CountryStructure, ::CPIVarInterm)
    # Acá se llama a inflfn(base, i), en donde base es de tipo VarCPIBase e i es la posición del vector de exclusión.
    # Esta es la función que debe definirse para cualquier medida de inflación.
    l = length(cs.base)
    vm = mapfoldl(i -> inflfn(cs.base[i],i), vcat, 1:l)
    vm
end


text:
using Test
using NQCDynamics
using NQCDynamics.Calculators
using LinearAlgebra: tr, Diagonal, eigvecs, eigvals
using RingPolymerArrays: RingPolymerArrays

@testset "General constructors" begin
    model = NQCModels.DoubleWell()
    @test Calculators.Calculator(model, 1, Float64) isa Calculators.DiabaticCalculator
    @test Calculators.Calculator(model, 1, 2, Float64) isa Calculators.RingPolymerDiabaticCalculator
    model = NQCModels.Free()
    @test Calculators.Calculator(model, 1, Float64) isa Calculators.AdiabaticCalculator
    @test Calculators.Calculator(model, 1, 2, Float64) isa Calculators.RingPolymerAdiabaticCalculator
    model = NQCModels.MiaoSubotnik()
    @test Calculators.Calculator(model, 1, Float64) isa Calculators.LargeDiabaticCalculator
end

@testset "AdiabaticCalculator" begin
    model = NQCModels.Free(3)
    calc = Calculators.AdiabaticCalculator{Float64}(model, 5)
    r = rand(3, 5)

    @test @allocated(Calculators.evaluate_potential!(calc, r)) == 0
    @test @allocated(Calculators.evaluate_derivative!(calc, r)) == 0

    @test @allocated(Calculators.get_potential(calc, r)) == 0
    @test @allocated(Calculators.get_derivative(calc, r)) == 0
end

@testset "RingPolymerAdiabaticCalculator" begin
    model = NQCModels.Free(3)
    calc = Calculators.RingPolymerAdiabaticCalculator{Float64}(model, 5, 10)
    r = rand(3, 5, 10)

    @test @allocated(Calculators.evaluate_potential!(calc, r)) == 0
    @test @allocated(Calculators.evaluate_derivative!(calc, r)) == 0

    @test @allocated(Calculators.get_potential(calc, r)) == 0
    @test @allocated(Calculators.get_derivative(calc, r)) == 0
end

@testset "DiabaticCalculator" begin
    model = NQCModels.DoubleWell()
    calc = Calculators.DiabaticCalculator{Float64}(model, 1)
    r = rand(1,1)
    quantities = (:potential, :derivative, :eigen, :adiabatic_derivative, :nonadiabatic_coupling)

    @testset "Potential evaluation" begin
        r = [0.1;;]
        true_potential = potential(model, r)

        @test calc.potential != true_potential
        @test getfield(calc, :potential).position != r

        Calculators.get_potential(calc, r)

        @test calc.potential == true_potential

        # Check that the position for potential is the only one that has been updated.
        for name in quantities
            if name !== :potential
                @test getfield(calc, name).position != r
            else
                @test getfield(calc, name).position == r
            end
        end
    end

    @testset "Dependent evaluation" begin
        r = [0.0;;]

        # Check the position fields are different
        for name in quantities
            @test getfield(calc, name).position != r
        end

        Calculators.get_nonadiabatic_coupling(calc, r)

        # Check the position fields are correct
        for name in quantities
            @test getfield(calc, name).position == r
        end

        # Check all the results
        @test calc.potential ≈ potential(model, r)
        @test calc.derivative ≈ derivative(model, r)
        @test calc.eigen.values ≈ eigvals(potential(model, r))
        @test abs.(calc.eigen.vectors) ≈ abs.(eigvecs(potential(model, r)))
        @test isapprox(
            calc.adiabatic_derivative[1],
            calc.eigen.vectors' * derivative(model, r)[1] * calc.eigen.vectors
        )
    end

    @testset "Zero allocations" begin
        @test @allocated(Calculators.evaluate_potential!(calc, r)) == 0
        @test @allocated(Calculators.evaluate_derivative!(calc, r)) == 0
        @test @allocated(Calculators.evaluate_eigen!(calc, r)) == 0
        @test @allocated(Calculators.evaluate_adiabatic_derivative!(calc, r)) == 0
        @test @allocated(Calculators.evaluate_nonadiabatic_coupling!(calc, r)) == 0
    end
end

@testset "RingPolymerDiabaticCalculator" begin
    model = NQCModels.DoubleWell()

    reset_calc() = Calculators.RingPolymerDiabaticCalculator{Float64}(model, 1, 10)

    r = rand(1,1,10)
    r_centroid = rand(1,1)
    standard_quantities = (:potential, :derivative, :eigen, :adiabatic_derivative, :nonadiabatic_coupling)
    centroid_quantities = (:centroid, :centroid_potential, :centroid_derivative,
        :centroid_eigen, :centroid_adiabatic_derivative, :centroid_nonadiabatic_coupling)

    @testset "Potential evaluation" begin
        calc = reset_calc()
        true_potential = [potential(model, r[:,:,i]) for i=1:10]

        @test calc.potential != true_potential
        @test getfield(calc, :potential).position != r

        Calculators.get_potential(calc, r)

        @test calc.potential == true_potential

        # Check that the position for potential is the only one that has been updated.
        for name in standard_quantities
            if name !== :potential
                @test getfield(calc, name).position != r
            else
                @test getfield(calc, name).position == r
            end
        end
    end

    @testset "Dependent evaluation" begin
        calc = reset_calc()

        # Reset the position fields
        for name in standard_quantities
            getfield(calc, name).position .= NaN
        end

        Calculators.get_nonadiabatic_coupling(calc, r)

        # Check the position fields are correct
        for name in standard_quantities
            @test getfield(calc, name).position == r
        end

        # Check all the results
        @test calc.potential ≈ [potential(model, r[:,:,i]) for i=1:10]
        @test calc.derivative ≈ [derivative(model, r[k,j,i]) for k=1:1, j=1:1, i=1:10]
        @test [calc.eigen[i].values for i=1:10] ≈ eigvals.(calc.potential)
        @test [abs.(calc.eigen[i].vectors) for i=1:10] ≈ [abs.(eigvecs(calc.potential[i])) for i=1:10]
        @test isapprox(
            calc.adiabatic_derivative[1],
            calc.eigen[1].vectors' * derivative(model, r[:,:,1])[1] * calc.eigen[1].vectors
        )
    end

    @testset "Centroid dependent evaluation" begin
        calc = reset_calc()
        r_centroid = RingPolymerArrays.get_centroid(r)

        # Check the position fields are different
        for name in centroid_quantities
            @test getfield(calc, name).position != r
        end

        Calculators.get_centroid_nonadiabatic_coupling(calc, r)

        # Check the position fields are correct
        for name in centroid_quantities
            @test getfield(calc, name).position == r
        end

        # Check all the results
        @test calc.centroid_potential ≈ potential(model, r_centroid)
        @test calc.centroid_derivative ≈ derivative(model, r_centroid)
        @test calc.centroid_eigen.values ≈ eigvals(potential(model, r_centroid))
        @test abs.(calc.centroid_eigen.vectors) ≈ abs.(eigvecs(potential(model, r_centroid)))
        @test isapprox(
            calc.centroid_adiabatic_derivative[1],
            calc.centroid_eigen.vectors' * derivative(model, r_centroid)[1] * calc.centroid_eigen.vectors
        )
    end

    @testset "Extra quantities" begin
        calc = reset_calc()

        Calculators.get_traceless_potential(calc, r)
        Calculators.get_traceless_derivative(calc, r)
        Calculators.get_traceless_adiabatic_derivative(calc, r)

        @test calc.V̄[1] ≈ tr(calc.potential[1]) / nstates(model)
        @test tr(calc.traceless_potential[1]) ≈ 0 atol=eps(Float64)
        @test calc.D̄[1] ≈ tr(calc.derivative[1]) / nstates(model)
        @test tr(calc.traceless_derivative[1]) ≈ 0 atol=eps(Float64)
        @test calc.traceless_adiabatic_derivative[1] ≈ calc.adiabatic_derivative[1] .- Diagonal(fill(calc.D̄[1], nstates(model)))
    end

    @testset "Zero allocations" begin

        calc = reset_calc()

        # Call all functions to ensure they have been compiled
        Calculators.evaluate_potential!(calc, r)
        Calculators.evaluate_derivative!(calc, r)
        Calculators.evaluate_eigen!(calc, r)
        Calculators.evaluate_adiabatic_derivative!(calc, r)
        Calculators.evaluate_nonadiabatic_coupling!(calc, r)

        Calculators.evaluate_traceless_potential!(calc, r)
        Calculators.evaluate_V̄!(calc, r)
        Calculators.evaluate_traceless_derivative!(calc, r)
        Calculators.evaluate_D̄!(calc, r)
        Calculators.evaluate_traceless_adiabatic_derivative!(calc, r)

        Calculators.evaluate_centroid!(calc, r)
        Calculators.evaluate_centroid_potential!(calc, r)
        Calculators.evaluate_centroid_derivative!(calc, r)
        Calculators.evaluate_centroid_eigen!(calc, r)
        Calculators.evaluate_centroid_adiabatic_derivative!(calc, r)
        Calculators.evaluate_centroid_nonadiabatic_coupling!(calc, r)

        calc = reset_calc()

        # Check for no allocations
        @test @allocated(Calculators.evaluate_potential!(calc, r)) == 0
        @test @allocated(Calculators.evaluate_derivative!(calc, r)) == 0
        @test @allocated(Calculators.evaluate_eigen!(calc, r)) == 0
        @test @allocated(Calculators.evaluate_adiabatic_derivative!(calc, r)) == 0
        @test @allocated(Calculators.evaluate_nonadiabatic_coupling!(calc, r)) == 0

        @test @allocated(Calculators.evaluate_traceless_potential!(calc, r)) == 0
        @test @allocated(Calculators.evaluate_V̄!(calc, r)) == 0
        @test @allocated(Calculators.evaluate_traceless_derivative!(calc, r)) == 0
        @test @allocated(Calculators.evaluate_D̄!(calc, r)) == 0
        @test @allocated(Calculators.evaluate_traceless_adiabatic_derivative!(calc, r)) == 0

        @test @allocated(Calculators.evaluate_centroid!(calc, r)) == 0
        @test @allocated(Calculators.evaluate_centroid_potential!(calc, r)) == 48 # not sure why this isn't zero
        @test @allocated(Calculators.evaluate_centroid_derivative!(calc, r)) == 0
        @test @allocated(Calculators.evaluate_centroid_eigen!(calc, r)) == 0
        @test @allocated(Calculators.evaluate_centroid_adiabatic_derivative!(calc, r)) == 0
        @test @allocated(Calculators.evaluate_centroid_nonadiabatic_coupling!(calc, r)) == 0
    end
end

@testset "LargeDiabaticCalculator" begin
    model = NQCModels.MiaoSubotnik()
    calc = Calculators.LargeDiabaticCalculator{Float64}(model, 1)
    r = rand(1,1)

    Calculators.get_nonadiabatic_coupling(calc, r)
    @test all(x->x==1, values(calc.stats))

    Calculators.evaluate_potential!(calc, r)
    Calculators.evaluate_derivative!(calc, r)
    Calculators.evaluate_eigen!(calc, r)
    Calculators.evaluate_adiabatic_derivative!(calc, r)
    Calculators.evaluate_nonadiabatic_coupling!(calc, r)

    @test @allocated(Calculators.evaluate_potential!(calc, r)) == 0
    @test @allocated(Calculators.evaluate_derivative!(calc, r)) == 0
    @test @allocated(Calculators.evaluate_eigen!(calc, r)) == 56896 # nonzero due to eigenroutines
    @test @allocated(Calculators.evaluate_adiabatic_derivative!(calc, r)) == 0
    @test @allocated(Calculators.evaluate_nonadiabatic_coupling!(calc, r)) == 0
end

@testset "FrictionCalculator" begin
    model = CompositeFrictionModel(Free(), RandomFriction(1))
    calc = Calculators.FrictionCalculator{Float64}(model, 1)
    r = rand(1,1)

    Calculators.get_potential(calc, r)
    Calculators.get_derivative(calc, r)
    Calculators.get_friction(calc, r)
    @test all(x->x==1, values(calc.stats))

    Calculators.evaluate_potential!(calc, r)
    Calculators.evaluate_derivative!(calc, r)
    Calculators.evaluate_friction!(calc, r)
    @test @allocated(Calculators.evaluate_potential!(calc, r)) == 0
    @test @allocated(Calculators.evaluate_derivative!(calc, r)) == 0
    @test @allocated(Calculators.evaluate_friction!(calc, r)) == 192
end

@testset "RingPolymerFrictionCalculator" begin
    model = CompositeFrictionModel(Free(), RandomFriction(1))
    calc = Calculators.RingPolymerFrictionCalculator{Float64}(model, 1, 10)
    r = rand(1,1,10)

    Calculators.get_potential(calc, r)
    Calculators.get_derivative(calc, r)
    Calculators.get_friction(calc, r)
    @test all(x->x==1, values(calc.stats))

    Calculators.evaluate_potential!(calc, r)
    Calculators.evaluate_derivative!(calc, r)
    Calculators.evaluate_friction!(calc, r)
    @test @allocated(Calculators.evaluate_potential!(calc, r)) == 0
    @test @allocated(Calculators.evaluate_derivative!(calc, r)) == 0
    @test @allocated(Calculators.evaluate_friction!(calc, r)) == 1920
end

@testset "DiabaticFrictionCalculator" begin
    model = NQCModels.MiaoSubotnik()
    calc = Calculators.DiabaticFrictionCalculator{Float64}(model, 1)
    r = rand(1,1)

    Calculators.get_nonadiabatic_coupling(calc, r)

    Calculators.evaluate_potential!(calc, r)
    Calculators.evaluate_derivative!(calc, r)
    Calculators.evaluate_eigen!(calc, r)
    Calculators.evaluate_adiabatic_derivative!(calc, r)
    Calculators.evaluate_nonadiabatic_coupling!(calc, r)

    @test @allocated(Calculators.evaluate_potential!(calc, r)) == 0
    @test @allocated(Calculators.evaluate_derivative!(calc, r)) == 0
    @test @allocated(Calculators.evaluate_eigen!(calc, r)) == 56896 # nonzero due to eigenroutines
    @test @allocated(Calculators.evaluate_adiabatic_derivative!(calc, r)) == 0
    @test @allocated(Calculators.evaluate_nonadiabatic_coupling!(calc, r)) == 0
end


text:
using Dates

# create a Date / dateTime variable
d1 = Date(2021, 1, 25)
dt1 = DateTime(2021, 4, 12) # if the 'time' elements are not specified this will detail to 00:00:00
dt2 = DateTime(2021, 9, 16, 14, 22, 1)
println("Date and DateTime : $d1 / $dt1")
println("DateTime (2) : $dt2")

# parse a date using a DateFormat
df1 = DateFormat("y-m-d")
d2 = Date("2021-02-03", df1) # parse a date 'string'
println("Date (2) is : $d2 ( $(typeof(d2) ) )")

# calculate the difference between two Dates
d3 = Date(2021, 2, 10)
d4 = Date(2021, 4, 9)
dt3 = DateTime(2021, 5, 3, 10, 19, 11)
dt4 = DateTime(2021, 5, 3, 12, 9, 2)
println("Is $d3 < $d4 : $(d3 < d4)")
println("The difference between $d4 and $d3 is : $(d4 - d3)")
println("The difference between $dt4 and $dt3 is : $(dt4 - dt3)") # note answer is in milliseconds (Millisecond type)

# date query methods
println("The day of the week of $d3 is : $(Dates.dayofweek(d3))")
println("The day name of the week of $d3 is : $(Dates.dayname(d3))")

# there are also a number of 'adjuster' methods
println("The first day of the week where the date is $d4 is : $(Dates.firstdayofweek(d4))")
println("The last day of the month (Date $d4) is : $(Dates.lastdayofmonth(d4))")

max_stars_repo_name: rpillar/JuliaExamples

max_issues_repo_name: rpillar/JuliaExamples

max_forks_repo_name: rpillar/JuliaExamples


text:
# libpurex:
# Exponential families paramterised by their means

include("binary_search.jl");

struct Gaussian
    σ2;
end

# convenience
Gaussian() = Gaussian(1);

struct Bernoulli
end

struct Exponential
end

struct Poisson
end

# KL divergence

rel_entr(x, y) = x==0 ? 0. : x*log(x/y);

d(expfam::Gaussian,    μ, λ) = (μ-λ)^2/(2*expfam.σ2);
d(expfam::Bernoulli,   μ, λ) = rel_entr(μ, λ) + rel_entr(1-μ, 1-λ);
d(expfam::Exponential, μ, λ) = λ == 0 ? Inf : μ/λ - log(μ/λ) - 1;
d(expfam::Poisson,     μ, λ) = rel_entr(μ, λ) - μ + λ;

# solution for λ in \min_λ d(μ, λ) - λ x
# i.e. λ satisfying (λ - μ)ϕ''(λ) == x

invh(expfam::Gaussian,    μ, x) = μ + x*expfam.σ2;
invh(expfam::Bernoulli,   μ, x) = 2μ/(1-x+sqrt((x-1)^2 + 4*x*μ));
invh(expfam::Exponential, μ, x) = x > 0 ? Inf : 2μ/(1 + sqrt(1 - 4*x*μ));
invh(expfam::Poisson,     μ, x) = μ/(1-x);

sample(rng, expfam::Gaussian,    μ) = μ + sqrt(expfam.σ2)*randn(rng);
sample(rng, expfam::Bernoulli,   μ) = rand(rng) ≤ μ;
sample(rng, expfam::Exponential, μ) = randexp(rng)*μ;

# upward and downward confidence intervals
dup(expfam::Gaussian, μ, v) = μ + sqrt(2*expfam.σ2*v);
ddn(expfam::Gaussian, μ, v) = μ - sqrt(2*expfam.σ2*v);

function dup(expfam::Bernoulli, μ, v)
    μ == 1 ? 1. : binary_search(λ -> d(expfam, μ, λ) - v, μ, 1);
end

function ddn(expfam::Bernoulli, μ, v)
    μ == 0 ? 0. : binary_search(λ -> v - d(expfam, μ, λ), 0, μ);
end

function dup(expfam::Exponential, μ, v)
    binary_search(λ -> d(expfam, μ, λ) - v, μ, μ*exp(v+1));
end

function ddn(expfam::Exponential, μ, v)
    binary_search(λ -> v - d(expfam, μ, λ), 0, μ);
end


text:
# just for testing/debugging purposes
function junky(y,hrel)
    ℓ₀ = 3
    if true
        @assert -1 ≤ y ≤ 1
        @assert hrel ≥ 0.0
    end
    ξL = y * (ℓ₀-1) + ℓ₀ #  1 ≤ ξL = position relative to mass on left  (dim'less)
    ξR = y * (ℓ₀-1) - ℓ₀ # -1 ≥ ξR = position relative to mass on right (dim'less)
    if hrel ≈ 0.0
        return ξL^-4 - ξR^-4
    else
        aa = 0.711557; bb = -7.58607
        function g(hox)
            min(1/2.5,abs(hox))
        end
        function f(hox)
            (1+ aa*g(hox) + bb*g(hox)^2)^(-1)
        end
        function ww(ξ)
            sign(ξ)/ξ^4 * f(hrel/ξ)
        end
        wL = ww(ξL)
        wR = ww(ξR)
        return wL + wR
    end
end

max_stars_repo_name: diakosmos/razor

max_issues_repo_name: diakosmos/razor

max_forks_repo_name: diakosmos/razor