text:
struct Available
    header::String
    files::Vector{NamedTuple}
end
function Available()
    header = "YYYYMMDD HHMMSS   name of the file(up to 80 characters)"
    new(header, Vector{NamedTuple}())
end
function Available(path::String)
    lines = readlines(path)
    ioc = findfirst(x -> occursin("YYYYMMDD HHMMSS", x), lines)
    headerlines = isnothing(ioc) ? [] : lines[1:ioc[1]]
    header = join(headerlines, "\n")
    filelines = isnothing(ioc) ? lines : lines[ioc+1:end]
    files = NamedTuple[]
    for l in filelines
        sl = split(l)
        date = DateTime(sl[1]*sl[2], "yyyymmddHHMMSS")
        push!(files, (date=date, filename=sl[3]))
    end
    Available(header, files)
end
Available(fpdir::FlexpartDir) = Available(getdir(fpdir, :available))
Base.show(io::IO, available::Available) = Base.display(available.files)

function write(available::Available, path::String)
    (tmppath, tmpio) = mktemp()

    for line in format(available) Base.write(tmpio, line*"\n") end

    close(tmpio)
    dest = path

    mv(tmppath, dest, force=true)
end

function write(fpdir::FlexpartDir, available::Available)
    write(available, getdir(fpdir, :available))
end

function dateYY(d)
    y = Dates.year(d)
    if 80 <= y <= 99
        d+Dates.Year(1900)
    elseif 0 <= y <= 79
        d+Dates.Year(2000)
    else
        error("don't know what to do with year $d")
    end
end

function format(available::Available)
    strs = [available.header]
    for elem in available.files
        str = Dates.format(elem.date, "yyyymmdd")*" "*Dates.format(elem.date, "HHMMSS")*"      "*elem.filename*"      "*"ON DISK"
        push!(strs, str)
    end
    strs
end

function update(available::Available, inputdir::String)
    files = readdir(inputdir)
    newfiles = NamedTuple[]

    for file in files
        m = match(r"^([A-Z]*)(\d{8,10})$", file)
        if !isnothing(m)
            x = m.captures[2]
            m_sep = parse.(Int, [x[1:2], x[3:4], x[5:6], x[7:8]])
            formated_date = DateTime(m_sep...)
            push!(newfiles, (date=dateYY.(formated_date), filename=file))
        end
    end
    Available(available.header, newfiles)
end

function update(available::Available, dates::Vector{DateTime}, prefix::String)
    fnames = prefix .* Dates.format.(dates, "yymmddHH")
    Available(available.header, [(date=date, filename=fname) for (date, fname) in zip(dates, fnames)])
end

function updated_available(fpdir::FlexpartDir)
    update(Available(fpdir), getdir(fpdir, :input))
end

function formatinput(fp_input_path::String)
    fp_input_files = readdir(fp_input_path)
    m = [match(r"\d{8}$", file).match for file in fp_input_files]
    m_sep = [parse.(Int, [x[1:2], x[3:4], x[5:6], x[7:8]]) for x in m]
    formated_date = [DateTime(y, m, d, h) for (y, m, d, h) in m_sep]
    dates = dateYY.(formated_date)
    [Dates.format(d, "yyyymmdd")*" "*Dates.format(d, "HHMMSS")*"      "*fn*"      "*"ON DISK" for (d, fn) in zip(dates, fp_input_files)]

    # write_formated_av(formated_av)
end

function formatinput(dates:: Array{DateTime, 1}, prefix::String)
    fname = prefix .* Dates.format.(dates, "yymmddHH")
    [Dates.format(d, "yyyymmdd")*" "*Dates.format(d, "HHMMSS")*"      "*fn*"      "*"ON DISK" for (d, fn) in zip(dates, fname)]
end

function update_available(fpdir::FlexpartDir, dates:: Array{DateTime, 1}, prefix::String)
    formated = formatinput(dates, prefix)
    avpath = getdir(fpdir, :available)
    update_available(avpath, formated)
end

function update_available(fpdir::FlexpartDir)
    formated = formatinput(getdir(fpdir, :input))
    avpath = getdir(fpdir, :available)
    update_available(avpath, formated)
end

function update_available(avpath::String, formated::Vector{<:String})
    av_file = readlines(avpath)
    ioc = findfirst(x -> occursin("YYYYMMDD HHMMSS", x), av_file)
    new_av = isnothing(ioc) ? [] : av_file[1:ioc[1]]
    for l in formated
        push!(new_av, l)
    end
    open(avpath, "w") do f
        for l in new_av
            Base.write(f, l*"\n")
        end
    end
end


text:
## Converison

#ensure that COnversion is called
coefficients(cfs::AbstractVector,A::Fourier{DD,R1},B::Laurent{DD,R2}) where {DD,R1,R2} =
    mul_coefficients(Conversion(A,B),cfs)
coefficients(cfs::AbstractVector,A::Laurent{DD,R1},B::Fourier{DD,R2}) where {DD,R1,R2} =
    mul_coefficients(Conversion(A,B),cfs)

hasconversion(::Fourier{DD,R1},::Laurent{DD,R2}) where {DD,R1,R2} = true
hasconversion(::Laurent{DD,R1},::Fourier{DD,R2}) where {DD,R1,R2} = true

Conversion(a::Laurent{DD,R1},b::Fourier{DD,R2}) where {DD,R1,R2} = ConcreteConversion(a,b)
Conversion(a::Fourier{DD,R1},b::Laurent{DD,R2}) where {DD,R1,R2} = ConcreteConversion(a,b)

function getindex(C::ConcreteConversion{Laurent{DD,R1},Fourier{DD,R2},T},k::Integer,j::Integer) where {DD,R1,R2,T}
    if k==j==1
        one(T)
    elseif iseven(k) && k==j
        -one(T)*im
    elseif iseven(k) && k+1==j
        one(T)*im
    elseif isodd(k) && (k==j || k-1==j )
        one(T)
    else
        zero(T)
    end
end

function getindex(C::ConcreteConversion{Fourier{DD,R1},Laurent{DD,R2},T},k::Integer,j::Integer) where {DD,R1,R2,T}
    if k==j==1
        one(T)
    elseif iseven(k) && k==j
        one(T)/2*im
    elseif iseven(k) && k+1==j
        one(T)/2
    elseif isodd(k) && k==j
        one(T)/2
    elseif isodd(k) && j==k-1
        -one(T)*im/2
    else
        zero(T)
    end
end

bandwidths(::ConcreteConversion{Laurent{DD,R1},Fourier{DD,R2}}) where {DD,R1,R2} = 1,1
bandwidths(::ConcreteConversion{Fourier{DD,R1},Laurent{DD,R2}}) where {DD,R1,R2} = 1,1

for RULE in (:conversion_rule,:maxspace_rule,:union_rule)
    @eval begin
        # override both to avoid SumSpace overrides
        function $RULE(A::Laurent{DD,R1},B::Fourier{DD,R2}) where {DD,R1,R2}
            @assert domainscompatible(A,B)
            B
        end
        function $RULE(A::Fourier{DD,R1},B::Laurent{DD,R2}) where {DD,R1,R2}
            @assert domainscompatible(A,B)
            A
        end
    end
end

conversion_type(A::Fourier{DD,R1},B::Fourier{DD,R2}) where {DD<:Circle,R1,R2} = domain(A).orientation ? A : B

hasconversion(A::Fourier{DD,R1},B::Fourier{DD,R2}) where {DD,R1,R2} = domain(A) == reverseorientation(domain(B))
function Conversion(A::Fourier{DD,R1},B::Fourier{DD,R2}) where {DD,R1,R2}
    if A==B
        ConversionWrapper(Operator(I,A))
    else
        @assert domain(A) == reverseorientation(domain(B))
        ConcreteConversion(A,B)
    end
end
bandwidths(::ConcreteConversion{Fourier{DD,R1},Fourier{DD,R2}}) where {DD,R1,R2}=0,0

getindex(C::ConcreteConversion{Fourier{DD,R1},Fourier{DD,R2},T},k::Integer,j::Integer) where {DD,R1,R2,T} =
    k==j ? (iseven(k) ? (-one(T)) : one(T)) : zero(T)

### Cos/Sine

function Derivative(S::Union{CosSpace,SinSpace},order)
    @assert isa(domain(S),PeriodicSegment)
    ConcreteDerivative(S,order)
end

bandwidths(D::ConcreteDerivative{CS}) where {CS<:CosSpace} = iseven(D.order) ? (0,0) : (0,1)
bandwidths(D::ConcreteDerivative{S}) where {S<:SinSpace} = iseven(D.order) ? (0,0) : (1,0)
rangespace(D::ConcreteDerivative{S}) where {S<:CosSpace} = iseven(D.order) ? D.space : SinSpace(domain(D))
rangespace(D::ConcreteDerivative{S}) where {S<:SinSpace} = iseven(D.order) ? D.space : CosSpace(domain(D))

function getindex(D::ConcreteDerivative{CS,OT,T},k::Integer,j::Integer) where {CS<:CosSpace,OT,T}
    d=domain(D)
    m=D.order
    C=convert(T,2/complexlength(d)*π)

    if k==j && mod(m,4)==0
        (C*(k-1))^m
    elseif k==j && mod(m,4)==2
        -(C*(k-1))^m
    elseif j==k+1 && mod(m,4)==1
        -(C*k)^m
    elseif j==k+1 && mod(m,4)==3
        (C*k)^m
    else
        zero(T)
    end
end

function getindex(D::ConcreteDerivative{CS,OT,T},k::Integer,j::Integer) where {CS<:SinSpace,OT,T}
    d=domain(D)
    m=D.order
    C=convert(T,2/complexlength(d)*π)

    if k==j && mod(m,4)==0
        (C*k)^m
    elseif k==j && mod(m,4)==2
        -(C*k)^m
    elseif j==k-1 && mod(m,4)==1
        (C*j)^m
    elseif j==k-1 && mod(m,4)==3
        -(C*j)^m
    else
        zero(T)
    end
end

# Use Laurent derivative
Derivative(S::Fourier{DD,RR},k::Integer) where {DD<:Circle,RR} =
    DerivativeWrapper(Derivative(Laurent(S),k)*Conversion(S,Laurent(S)),k)

Integral(::CosSpace,m::Integer) =
    error("Integral not defined for CosSpace.  Use Integral(CosSpace()|(2:∞)) if first coefficient vanishes.")

Integral(sp::SinSpace{<:PeriodicSegment}, m::Integer) = ConcreteIntegral(sp,m)

bandwidths(D::ConcreteIntegral{<:SinSpace}) = iseven(D.order) ? (0,0) : (1,0)
rangespace(D::ConcreteIntegral{<:CosSpace}) = iseven(D.order) ? D.space : SinSpace(domain(D))
rangespace(D::ConcreteIntegral{<:SinSpace}) = iseven(D.order) ? D.space : CosSpace(domain(D))

function getindex(D::ConcreteIntegral{CS,OT,T},k::Integer,j::Integer) where {CS<:SinSpace,OT,T}
    d=domain(D)
    @assert isa(d,PeriodicSegment)
    m=D.order
    C=convert(T,2/complexlength(d)*π)

    if k==j && mod(m,4)==0
        (C*k)^(-m)
    elseif k==j && mod(m,4)==2
        -(C*k)^(-m)
    elseif j==k-1 && mod(m,4)==1
        -(C*j)^(-m)
    elseif j==k-1 && mod(m,4)==3
        (C*j)^(-m)
    else
        zero(T)
    end
end

function Integral(S::SubSpace{<:CosSpace,<:AbstractInfUnitRange{Int},<:PeriodicSegment},k::Integer)
    @assert first(S.indexes)==2
    ConcreteIntegral(S,k)
end

bandwidths(D::ConcreteIntegral{<:SubSpace{<:CosSpace,<:AbstractInfUnitRange{Int},<:PeriodicSegment}}) =
    (0,0)
rangespace(D::ConcreteIntegral{<:SubSpace{<:CosSpace,<:AbstractInfUnitRange{Int},<:PeriodicSegment}}) =
    iseven(D.order) ? D.space : SinSpace(domain(D))

function getindex(D::ConcreteIntegral{<:SubSpace{<:CosSpace,<:AbstractInfUnitRange{Int},<:PeriodicSegment}},
                  k::Integer,j::Integer)
    d=domain(D)
    m=D.order
    T=eltype(D)
    C=convert(T,2/complexlength(d)*π)

    if k==j
        if mod(m,4)==0
            (C*k)^(-m)
        elseif mod(m,4)==2
            -(C*k)^(-m)
        elseif mod(m,4)==1
        (C*k)^(-m)
        else   # mod(m,4)==3
            -(C*k)^(-m)
        end
    else
        zero(T)
    end
end

# CosSpace Multiplicaiton is the same as Chebyshev

Multiplication(f::Fun{CS},sp::CS) where {CS<:CosSpace} = ConcreteMultiplication(f,sp)
Multiplication(f::Fun{SS},sp::SS) where {SS<:SinSpace} = ConcreteMultiplication(f,sp)
Multiplication(f::Fun{CS},sp::SinSpace) where {CS<:CosSpace} = ConcreteMultiplication(f,sp)
function Multiplication(f::Fun{SS},sp::CosSpace) where SS<:SinSpace
    @assert domain(f) == domain(sp)
    a=f.coefficients/2
    if length(a) == 0
        A=ZeroOperator(sp,SinSpace(domain(sp)))
        MultiplicationWrapper(f,A)
    elseif length(a) == 1
        A=ToeplitzOperator([0.],[a[1];0.;-a]) + HankelOperator(a)
        MultiplicationWrapper(f,SpaceOperator(A,sp,SinSpace(domain(sp))))
    else
        A=ToeplitzOperator(a[2:end],[a[1];0.;-a]) + HankelOperator(a)
        MultiplicationWrapper(f,SpaceOperator(A,sp,SinSpace(domain(sp))))
    end
end

bandwidths(M::ConcreteMultiplication{CS,CS}) where {CS<:CosSpace} =
    (ncoefficients(M.f)-1,ncoefficients(M.f)-1)
rangespace(M::ConcreteMultiplication{CS,CS}) where {CS<:CosSpace} = domainspace(M)
getindex(M::ConcreteMultiplication{CS,CS},k::Integer,j::Integer) where {CS<:CosSpace} =
    chebmult_getindex(M.f.coefficients,k,j)

function getindex(M::ConcreteMultiplication{SS,SS},k::Integer,j::Integer) where SS<:SinSpace
    a=M.f.coefficients
    ret=toeplitz_getindex([zero(eltype(a));-a],a,k,j)/2
    if k ≥ 2
        ret+=hankel_getindex(a,k,j)/2
    end
    ret
end

bandwidths(M::ConcreteMultiplication{SS,SS}) where {SS<:SinSpace}=ncoefficients(M.f)+1,ncoefficients(M.f)-1
rangespace(M::ConcreteMultiplication{SS,SS}) where {SS<:SinSpace}=CosSpace(domain(M))

function getindex(M::ConcreteMultiplication{Cs,SS},k::Integer,j::Integer) where {SS<:SinSpace,Cs<:CosSpace}
    a=M.f.coefficients
    ret=toeplitz_getindex(a,k,j)/2
    if length(a)>=3
        ret-=hankel_getindex(view(a,3:length(a)),k,j)/2
    end
    ret
end

bandwidths(M::ConcreteMultiplication{Cs,SS}) where {SS<:SinSpace,Cs<:CosSpace} =
    (ncoefficients(M.f)-1,ncoefficients(M.f)-1)
rangespace(M::ConcreteMultiplication{Cs,SS}) where {SS<:SinSpace,Cs<:CosSpace} =
    SinSpace(domain(M))

function Multiplication(a::Fun{Fourier{D,R},T},sp::Fourier{D,R}) where {T,D,R}
    d=domain(a)
    c,s=components(a)
    O=Operator{float(T)}[Multiplication(c,CosSpace(d)) Multiplication(s,SinSpace(d));
                        Multiplication(s,CosSpace(d)) Multiplication(c,SinSpace(d))]
    MultiplicationWrapper(a,InterlaceOperator(O,space(a),sp))
end

coefficienttimes(f::Fun{Fourier{DD,RR}},g::Fun{Fourier{DD,RR}}) where {DD,RR} = Multiplication(f,space(g))*g

## Definite integral

for SP in (:CosSpace,:SinSpace,:Fourier)
    @eval begin
        DefiniteIntegral(S::$SP{D,R}) where {D,R} =
            ConcreteDefiniteIntegral{typeof(S),prectype(S)}(S)
        DefiniteLineIntegral(S::$SP{D,R}) where {D,R} =
            ConcreteDefiniteLineIntegral{typeof(S),real(prectype(S))}(S)
    end
end

getindex(Σ::ConcreteDefiniteIntegral{CosSpace{D,R},T},k::Integer) where {T,D<:PeriodicSegment,R} =
    k == 1 ? convert(T,complexlength(domain(Σ))) : zero(T)

getindex(Σ::ConcreteDefiniteIntegral{SinSpace{D,R},T},k::Integer) where {T,D<:PeriodicSegment,R} =
    zero(T)

getindex(Σ::ConcreteDefiniteIntegral{Fourier{D,R},T},k::Integer) where {T,D<:PeriodicSegment,R} =
    k == 1 ? convert(T,complexlength(domain(Σ))) : zero(T)

getindex(Σ::ConcreteDefiniteIntegral{CosSpace{D,R},T},k::Integer) where {T,D<:Circle,R} =
    k==2 ? convert(T,complexlength(domain(Σ))/2) : zero(T)

getindex(Σ::ConcreteDefiniteIntegral{SinSpace{D,R},T},k::Integer) where {T,D<:Circle,R} =
    k == 1 ? convert(T,0.5im*complexlength(domain(Σ))) : zero(T)

getindex(Σ::ConcreteDefiniteIntegral{Fourier{D,R},T},k::Integer) where {T,D<:Circle,R} =
    k == 2 ? convert(T,0.5im*complexlength(domain(Σ))) : (k==3 ? convert(T,complexlength(domain(Σ))/2) : zero(T))

getindex(Σ::ConcreteDefiniteLineIntegral{CosSpace{D,R},T},k::Integer) where {T,D<:PeriodicSegment,R} =
    k==1 ? convert(T,arclength(domain(Σ))) : zero(T)

getindex(Σ::ConcreteDefiniteLineIntegral{SinSpace{D,R},T},k::Integer) where {T,D<:PeriodicSegment,R} =
    zero(T)

getindex(Σ::ConcreteDefiniteLineIntegral{Fourier{D,R},T},k::Integer) where {T,D<:PeriodicSegment,R} =
    k==1 ? convert(T,arclength(domain(Σ))) : zero(T)

getindex(Σ::ConcreteDefiniteLineIntegral{CosSpace{D,R},T},k::Integer) where {T,D<:Circle,R} =
    k==1 ? convert(T,arclength(domain(Σ))) : zero(T)

getindex(Σ::ConcreteDefiniteLineIntegral{SinSpace{D,R},T},k::Integer) where {T,D<:Circle,R} =
    zero(T)

getindex(Σ::ConcreteDefiniteLineIntegral{Fourier{D,R},T},k::Integer) where {T,D<:Circle,R} =
    k==1 ? convert(T,arclength(domain(Σ))) : zero(T)

bandwidths(Σ::ConcreteDefiniteIntegral{CosSpace{D,R}}) where {D<:PeriodicSegment,R} = 0,0
bandwidths(Σ::ConcreteDefiniteIntegral{SinSpace{D,R}}) where {D<:PeriodicSegment,R} = 0,0
bandwidths(Σ::ConcreteDefiniteIntegral{Fourier{D,R}}) where {D<:PeriodicSegment,R} = 0,0
bandwidths(Σ::ConcreteDefiniteIntegral{CosSpace{D,R}}) where {D<:Circle,R} = 0,1
bandwidths(Σ::ConcreteDefiniteIntegral{SinSpace{D,R}}) where {D<:Circle,R} = 0,0
bandwidths(Σ::ConcreteDefiniteIntegral{Fourier{D,R}}) where {D<:Circle,R} = 0,2
bandwidths(Σ::ConcreteDefiniteLineIntegral{CosSpace{D,R}}) where {D,R} = 0,0
bandwidths(Σ::ConcreteDefiniteLineIntegral{SinSpace{D,R}}) where {D,R} = 0,0
bandwidths(Σ::ConcreteDefiniteLineIntegral{Fourier{D,R}}) where {D,R} = 0,0

transformtimes(f::Fun{CS},g::Fun{Fourier{D,R}}) where {CS<:CosSpace,D,R} =
    transformtimes(Fun(Fourier(domain(f)),interlace(f.coefficients,zeros(cfstype(f),ncoefficients(f)-1))),g)
transformtimes(f::Fun{SS},g::Fun{Fourier{D,R}}) where {SS<:SinSpace,D,R} =
    transformtimes(Fun(Fourier(domain(f)),interlace(zeros(cfstype(f),ncoefficients(f)+1),f.coefficients)),g)
transformtimes(f::Fun{CS},g::Fun{SS}) where {CS<:CosSpace,SS<:SinSpace} =
    transformtimes(Fun(Fourier(domain(f)),interlace(f.coefficients,zeros(cfstype(f),ncoefficients(f)-1))),
                    Fun(Fourier(domain(g)),interlace(zeros(cfstype(g),ncoefficients(g)+1),g.coefficients)))
transformtimes(f::Fun{Fourier{D,R}},g::Fun{CS}) where {CS<:CosSpace,D,R} = transformtimes(g,f)
transformtimes(f::Fun{Fourier{D,R}},g::Fun{SS}) where {SS<:SinSpace,D,R} = transformtimes(g,f)
transformtimes(f::Fun{SS},g::Fun{CS}) where {SS<:SinSpace,CS<:CosSpace} = transformtimes(g,f)

ReverseOrientation(S::Fourier{D}) where {D} = ReverseOrientationWrapper(NegateEven(S,reverseorientation(S)))
Reverse(S::Fourier{D}) where {D} = ReverseWrapper(NegateEven(S,S))

## Multivariate

for TYP in (:Fourier,:Laurent,:CosSpace,:SinSpace,:Taylor)
    @eval begin
        function Dirichlet(S::TensorSpace{Tuple{$TYP{PeriodicSegment{T},R},PS}},k) where {PS,T,R}
            op = [Operator(I,S.spaces[1])⊗Evaluation(S.spaces[2],leftendpoint,k);
                            ReverseOrientation(S.spaces[1])⊗Evaluation(S.spaces[2],rightendpoint,k) ]
            DirichletWrapper(SpaceOperator(op,S,PiecewiseSpace(rangespace(op).spaces)),k)
        end
        function Dirichlet(S::TensorSpace{Tuple{PS,$TYP{PeriodicSegment{T},R}}},k) where {PS,T,R}
            op = [Evaluation(S.spaces[1],leftendpoint,k)⊗Operator(I,S.spaces[2]);
                            Evaluation(S.spaces[1],rightendpoint,k)⊗ReverseOrientation(S[2]) ]
            DirichletWrapper(SpaceOperator(op,S,PiecewiseSpace(rangespace(op).spaces)),k)
        end
    end
end

ext: jl


text:
using QuadGK
using Roots

abstract type Distribution end

abstract type DiscreteDistribution <: Distribution end

"""
BinomialDistribution creates a struct for the binomial distribution.
See also: [`binomial_mean`](@ref), [`binomial_variance`](@ref), [`binomial_probability_mass_funcion`](@ref).
"""
struct BinomialDistribution <: DiscreteDistribution
    n::Integer
    p::AbstractFloat
    q::AbstractFloat
    function BinomialDistribution(n::Integer, p::AbstractFloat)
        q = 1 - p
        new(n, p, q)
    end
end

"""
Calculates the binomial mean for a binomial distribution.
"""
function binomial_mean(distribution::BinomialDistribution)
    mean = distribution.n * distribution.p
    return mean
end

"""
Calculates the variance for a binomial distribution.
"""
function binomial_variance(distribution::BinomialDistribution)
    variance = distribution.n * distribution.p * distribution.q
    return variance
end

"""
Calculates the probability that ``P(Y = k)`` for a binomial distribution ``Y``.
"""
function binomial_probability_mass_funcion(distribution::BinomialDistribution, k::Integer)
    pᵧ = binomial(distribution.n, k) * distribution.p^k * distribution.q^(distribution.n - k)
    return pᵧ
end

abstract type ContiniousDistribution <: Distribution end

struct NormalDistribution <: ContiniousDistribution
    μ::AbstractFloat
    σ²::AbstractFloat
end

function normal_probability_density_function(distribution::NormalDistribution, x::AbstractFloat)
    μ = distribution.μ
    σ = √(distribution.σ²)
    f = 1 / (σ * √(2 * π)) * ℯ^(-1 / 2 * ((x - μ) / σ)^2)
    return f
end

function normal_cumulative_distribution_function(distribution::NormalDistribution, x::AbstractFloat)
    X = distribution
    F, error = quadgk(t -> normal_probability_density_function(X, t), -Inf, x)
    return F
end

abstract type GammaDistributions <: ContiniousDistribution end

struct GammaDistribution <: GammaDistributions
    α::AbstractFloat
    β::AbstractFloat
end

function gamma_probability_density_function(distribution::GammaDistribution, x::AbstractFloat)
    α = distribution.α
    β = distribution.β
    f = β^α / gamma_function(α) * x^(α - 1) * ℯ^(-β * x)
end

function gamma_cumulative_distribution_function(distribution::GammaDistribution, x::AbstractFloat)
    X = distribution
    F, error = quadgk(t -> gamma_probability_density_function(X, t), -Inf, x)
    return F
end

function gamma_function(α::AbstractFloat)
    α > 0 || error("α must be in ℝ > 0")
    Γ, error = quadgk(x -> x^(α - 1) * exp(-x), 0, Inf)
    return Γ
end

function gamma_function(α::Integer)
    α > 0 || error("α must be in ℝ > 0")
    Γ = factorial(α - 1)
    return Γ
end

struct ExponentialDistribution <: GammaDistributions
    α::Integer
    β::AbstractFloat
    function ExponentialDistribution(β::AbstractFloat)
        α = 1
        return new(α, β)
    end
end

function exponential_probability_density_function(distribution::ExponentialDistribution, x::AbstractFloat)
    f = gamma_probability_density_function(distribution, x)
end

function exponential_cumulative_distribution_function(distribution::ExponentialDistribution, x::AbstractFloat)
    F = gamma_cumulative_distribution_function(distribution, x)
end

struct ChiSquaredDistribution <: GammaDistributions
    df::Integer
    α::AbstractFloat
    β::AbstractFloat
    function ChiSquaredDistribution(df::Integer)
        α = df / 2
        β = 1 / 2
        return new(df, α, β)
    end
end

function chi_squared_probability_density_function(distribution::ChiSquaredDistribution, x::AbstractFloat)
    df = distribution.df
    C = 1 / (gamma_function(df / 2) * 2^(df / 2))
    f = C * x^(df / 2 - 1) * ℯ^(-x / 2)
end

function chi_squared_cumulative_distribution_function(distribution::ChiSquaredDistribution, x::AbstractFloat)
    X = distribution
    F, error = quadgk(t -> chi_squared_probability_density_function(X, t), 0, x)
    return F
end

abstract type tDistributions <: ContiniousDistribution end

struct tDistribution <: tDistributions
    df::Integer
end

function t_probability_density_function(distribution::tDistributions, x::AbstractFloat)
    df = distribution.df
    C = 1 / √(df * π) * gamma_function((df + 1) / 2) / gamma_function(df / 2)
    f = C * (1 + x^2 / df)^((-df - 1) / 2)
end

function t_cumulative_distribution_function(distribution::tDistributions, x::AbstractFloat)
    X = distribution
    F, error = quadgk(t -> t_probability_density_function(X, t), -Inf, x)
    return F
end

struct CauchyDistribution <: tDistributions
    df::Integer
    function CauchyDistribution()
        df = 1
        new(df)
    end
end

function quantile_finder(distribution::ChiSquaredDistribution, α::AbstractFloat)
    f(x) = chi_squared_cumulative_distribution_function(distribution, x) - 1 + α
    quantile = find_zero(f, (0, 200))
    return quantile
end

function quantile_finder(distribution::NormalDistribution, α::AbstractFloat)
    f(x) = normal_cumulative_distribution_function(distribution, x) - 1 + α
    quantile = find_zero(f, (-5, 5))
    return quantile
end

function quantile_finder(distribution::tDistributions, α::AbstractFloat)
    f(x) = t_cumulative_distribution_function(distribution, x) - 1 + α
    quantile = find_zero(f, (1, 1000))
    return quantile
end

max_issues_repo_name: persson-io/Stochastics

max_forks_repo_name: persson-io/Stochastics


text:
windows = false
cd("/remote/ceph/group/gedet/data/sim/2020/2020-02-06_8380701d_st_ffischer")
include("gerda_data_processing/init.jl");
include("gerda_data_processing/fct.jl");

# All
ch_table = Table(ch = [0, 1, 2, 3, 4, 5, 6, 7, 7, 7, 8, 9, 10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36],
                run = [86,86,86,86,86,78,78,67,66,65,86,86,86,86,86,87,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,83,86,86,86,86]);
#= # Just outliers
ch_table = Table(ch = [5, 6, 7, 7, 7, 13,32],
                run = [78,78,67,66,65,87,83]);
# Just run0086
ch_table = Table(ch = [0, 1, 2, 3, 4, 8, 9, 10,11,12,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,33,34,35,36],
                run = [86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86,86]); =#
#= ch_table = Table(ch = [0],
                run = [86]); =#
current_dir = pwd()
cal = JSON.parsefile("data_pulses/FEP_pulses/calibration/cal.json");

pre_cut = 50        # keV
multiplicity = 1;
std_intersect = 2;

for s in unique(ch_table.run)

    set = findall(x->x==s, ch_table.run)[1]
    int_run = ch_table[set].run;
    str_run = lpad(int_run, 4, "0");

    ddir = "/remote/ceph/group/gerda/data/phase2/blind/v04.00/gen"
    cd(ddir*"/tier4/all/phy/run"*str_run*"/")
    files = glob("*.root")
    run_keys  = []
    for meta_key in files
        push!(run_keys, split(meta_key, "-phy-all-tier")[1])
    end
    cd(current_dir)
    @info("Run keys for "* str_run *" are gathered! Number of files: "*string(length(run_keys)))
    #= pro = Progress(length(run_keys), dt=0.5,
                        barglyphs=BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
                        barlen=10) =#
    for id in run_keys
        filtered_dict = Dict();
        for ch in unique(ch_table.ch)
            filtered_dict[ch] = Table(  energy       = [],
                                        multiplicity = [],
                                        timestamp    = [],
                                        #det          = [],
                                        channel      = [],
                                        waveform     = []);
        end
        channel_index   = findall(x->x==ch_table[set].run, ch_table.run);
        unique_channels = []
        new_id = true
        for eval_channel in unique(ch_table.ch[channel_index])
            ch_str = lpad(eval_channel, 2, "0");
            filename = "data_pulses/baselines_raw/run"*str_run*"/"*ch_str*"-"*channel_to_name[eval_channel]*"/"
            filename *= id*"-phy-bl-raw.h5"
            if !isfile(filename)
                push!(unique_channels, eval_channel)
            else
                @info("File already exist: "*filename)
            end
        end
        unique_channels = unique(unique_channels)

        filename1 = ddir*"/tier1/ged/phy/run"*str_run*"/"*id*"-phy-ged-tier1.root";
        filename4 = ddir*"/tier4/all/phy/run"*str_run*"/"*id*"-phy-all-tier4.root"

        #>------- Tier4 ---------------------------------------<#
        #>-----------------------------------------------------<#
        tree = TFile(filename4)["tier4"];

        Base.run(`clear`)
        current_file_number = findall(x->x == id, run_keys)[1]
        @info("File "*string(current_file_number)*" of "*string(length(run_keys)))
        @info(">---------------------------------------------------------------------<")

        for eval_channel in unique_channels
            index = [];
            #push!(index, findall(y->sum(y) >= E - ΔE && sum(y) <= E + ΔE, tree.energy[1:end]));     # Select energy range
            push!(index, findall(x->x[eval_channel + 1] == 0, tree.energy[1:end]));                  # Select channel
            #push!(index, findall(x->x <= 1, tree.multiplicity[1:end]));                             # Select multiplicity
            #push!(index, findall(x->x==0, tree.isBL[1:end]));                                       # Select "is baseline"
            #push!(index, findall(x->x==0, tree.isTP[1:end]));                                       # Select "no test pulse"
            # Add additional filters here
            # push!(index, findall(..., ...))

            filtered_index = index[1]
            #=
            for j in 2:1:length(index)-1
                for r in 1:1:length(filtered_index)
                    if !(filtered_index[r] in index[j])
                        filtered_index[r] = 0
                    end
                end
            end
            =#
            #filtered_index = filtered_index[findall(x->x != 0, filtered_index)];

            if length(filtered_index) > 0
                energy       = sum.(tree.energy[filtered_index])
                multiplicity = tree.multiplicity[filtered_index]
                timestamp    = tree.timestamp[filtered_index]
                eventNumber  = tree.eventNumber[filtered_index]
                #ch           = findall.(x-> x > 5, tree.energy[filtered_index]);
                eventChannelNumber = tree.eventChannelNumber[filtered_index]
                channel  = zeros(length(filtered_index)) .+ eval_channel
                #detector = []
                #for c in ch
                #    push!(channel, c[1]-1)
                #    push!(detector, channel_to_name[c[1]-1])
                #end
                events    = TypedTables.Table(raw2mgtevent.(TFile(filename1)["MGTree"].event[filtered_index]))
                waveforms = events.fAuxWaveforms;
                waveform  = []
                for w in waveforms
                    wf_time  = w[eval_channel+1].wf.time
                    #bl       = sum(w[eval_channel+1].wf.value) / length(w[eval_channel+1].wf.value)
                    wf_value = w[eval_channel+1].wf.value #.- bl
                    #wf_value .*= cal[string(eval_channel)]["value"]
                    push!(waveform, RDWaveform(wf_time, wf_value))
                end
                events   = nothing;
                append!(filtered_dict[eval_channel], Table( energy       = energy,
                                                            multiplicity = multiplicity,
                                                            timestamp    = timestamp,
                                                            #det          = detector,
                                                            channel      = channel,
                                                            waveform     = waveform));
            end
            @info("Ch"*string(eval_channel)*"-"*channel_to_name[eval_channel]*" | Events: "*string(length(filtered_dict[eval_channel].energy)))
            #= next!(pro) =#
        end

        @info(">---------------------------------------------------------------------<")
        @info("Start storing the raw baselines")

        for eval_channel in unique_channels
            if length(filtered_dict[eval_channel].energy) > 0

                #>------- Store channel -------------------------------<#
                #>-----------------------------------------------------<#
                ch_str = lpad(eval_channel, 2, "0");
                filename = "data_pulses/baselines_raw/run"*str_run*"/"*ch_str*"-"*channel_to_name[eval_channel]*"/"
                if !isdir(filename)
                    mkpath(filename)
                end

                #= str_start = lpad(1, 7, "0")
                str_end   = lpad(length(bl_index), 7, "0")
                filename *= str_start*"-"*str_end*".h5" =#
                bl_index = 1:1:length(filtered_dict[eval_channel].energy)
                filename *= id*"-phy-bl-raw.h5"
                HDF5.h5open(filename, "w") do h5f
                    LegendHDF5IO.writedata( h5f, "data", Table( energy       = float.(filtered_dict[eval_channel].energy[bl_index]),
                                                                multiplicity = float.(filtered_dict[eval_channel].multiplicity[bl_index]),
                                                                timestamp    = float.(filtered_dict[eval_channel].timestamp[bl_index]),
                                                                #det          = string.(channel_data.det),
                                                                channel      = float.(filtered_dict[eval_channel].channel[bl_index]),
                                                                waveform     = StructArray{RDWaveform}(Array{RDWaveform,1}(filtered_dict[eval_channel].waveform[bl_index])),
                                                                #wf_cut       = StructArray{RDWaveform}(Array{RDWaveform,1}(channel_data.wf_cut))
                                                                ))
                end

                @info("Ch"*string(eval_channel)*"-"*channel_to_name[eval_channel]*" | "*string(length(bl_index))*" baselines stored!")
            else
                @info("Ch"*string(eval_channel)*"-"*channel_to_name[eval_channel]*" | No baselines found!")
            end
        end
    end

end

include("Filter_by_BL_converter.jl")


text:
using TropicalTensors.Reversible, TropicalTensors
using Yao, TropicalYao
using Test
using Viznet
using NiLang, NiLang.AD

@testset "spinglass" begin
    sg = Spinglass(SquareLattice(8, 10), ones(Int32, 142), zeros(Int32, 80))
    Lx, Ly = size(sg.lattice)
    reg = ArrayReg(ones(Tropical{Int32}, 1<<Ly))
    A = stack4reg(reg, Ly)
    B = stack4reg(reg, Lx-1)
    @test check_inv(isolve, (Int32(0), sg, reg, A, B))
    reg = ArrayReg(ones(Tropical{Int32}, 1<<Ly))
    A = stack4reg(reg, Ly)
    B = stack4reg(reg, Lx-1)
    @test isolve(Int32(0), sg, reg, A, B)[1] === Int32(142)
end

@testset "optconfig" begin
    function opt_config_largemem(sg::Spinglass{LT,T}) where {LT,T}
        reg = ArrayReg(ones(Tropical{T}, 1<<sg.lattice.Ny))
        A = stack4reg(reg, length(sg.lattice))
        eng, sg, reg, A = isolve_largemem(T(0.0), sg, reg, A)
        sgg = Spinglass(sg.lattice, GVar.(sg.Js, zero(sg.Js)), GVar.(sg.hs, zero(sg.hs)))
        gres = (~isolve_largemem)(GVar(eng, T(1)), sgg, GVar(reg), GVar(A))
        return TropicalTensors.SpinglassOptConfig(sg, eng, grad.(gres[2].Js), grad.(gres[2].hs))
    end
    sg = rand_spinglass(Int32, SquareLattice(9, 7); jt=Randpm(), ht=Zero(), seed=2)
    optc = opt_config(sg)
    optc2 = opt_config_largemem(sg)
    @test optc.grad_Js == optc2.grad_Js
    @test optc.grad_hs == optc2.grad_hs
    @test optc.eng == optc2.eng
end


text:
module TrussVisualization

using Requires

function __init__()
    @require Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" @eval begin
        include("makie.jl")
        export visualize
    end
end

end


text: abstract type StereoReader end


text:
"""
Module CommonUtils

General purpose utilities usable by all element types

"""

module ExplicitJacobians

using LinearAlgebra # for I matrix in geometricFactors
using ForwardDiff
using StaticArrays
using SparseArrays
using UnPack
using SetupDG

export init_jacobian_matrices
export hadamard_jacobian,accum_hadamard_jacobian!
export reduce_jacobian!
export hadamard_sum, hadamard_sum!
export banded_matrix_function
export columnize

# for constructing DG matrices
export build_rhs_matrix, assemble_global_SBP_matrices_2D

# use w/reshape to convert from matrices to arrays of arrays
columnize(A) = SVector{size(A,2)}([A[:,i] for i in 1:size(A,2)])

# sparse matrix assembly version
# can only deal with one coordinate component at a time in higher dimensions
function hadamard_jacobian(Q::SparseMatrixCSC,dF,U::AbstractArray,scale = -1)

    Nfields = length(U)
    NpK = size(Q,2)
    blockIds = repeat([NpK],Nfields)
    A = spzeros(NpK*Nfields,NpK*Nfields)

    accum_hadamard_jacobian!(A,Q,dF,U,scale)

    return A
end

# compute and accumulate contributions from a Jacobian function dF
function accum_hadamard_jacobian!(A::SparseMatrixCSC,Q::SparseMatrixCSC,
    dF,U::AbstractArray,scale = -1)

    Nfields = length(U)

    num_pts = size(Q,1)
    ids(m) = (1:num_pts) .+ (m-1)*num_pts
    Block(m,n) = CartesianIndices((ids(m),ids(n)))

    # loop over non-zero indices in Q
    Qnz = zip(findnz(Q)...)
    for (i,j,Qij) in Qnz
        Ui = getindex.(U,i)
        Uj = getindex.(U,j)
        dFdU = dF(Ui,Uj)
        for n = 1:length(U), m=1:length(U)
            A[Block(m,n)[i,j]] += dFdU[m,n]*Qij
        end
    end

    # add diagonal entry assuming Q = +/- Q^T
    for m = 1:Nfields, n = 1:Nfields
        #Asum = scale * vec(sum(A[Block(m,n)],dims=1)) # CartesianIndices broken
        Asum = scale * vec(sum(A[ids(m),ids(n)],dims=1)) # TODO: fix slowness (maybe just alloc issue?)
        # A[Block(m,n)] .+= diagm(Asum) # can't index all at once - bug related to CartesianIndices?
        for i = 1:num_pts
            A[Block(m,n)[i,i]] += Asum[i]
        end
    end
end

# computes block-banded matrix whose bands are entries of matrix-valued
# function evals (e.g., a Jacobian function).
function banded_matrix_function(mat_fun,U::AbstractArray)
    Nfields = length(U)
    num_pts = length(U[1])

    A = spzeros(Nfields*num_pts,Nfields*num_pts)
    ids(m) = (1:num_pts) .+ (m-1)*num_pts
    Block(m,n) = CartesianIndices((ids(m),ids(n)))

    for i = 1:num_pts
        mat_i = mat_fun(getindex.(U,i))
        for n = 1:Nfields, m = 1:Nfields
            A[Block(m,n)[i,i]] = mat_i[m,n] # TODO: replace with fast sparse constructor
        end
    end
    return A
end

# init sparse matrices based on elem-to-elem connectivity
# dims = tuple of matrix dimensions for each matrix (assumes square matrices)
function init_jacobian_matrices(md::MeshData, dims, Nfields=1)
    @unpack FToF = md
    Nfaces,K = size(FToF)
    EToE = @. (FToF.-1) ÷ Nfaces + 1
    A = spzeros.(K.*dims,K.*dims)

    ids(e,block_size) = (@. (1:block_size) + (e-1)*block_size)

    for e = 1:K, f = 1:Nfaces
        enbr = EToE[f,e]

        for (i,block_size) = enumerate(dims)
            id,id_nbr = ids.((e,enbr), block_size)

            # init to non-zeros to retain sparsity pattern
            fill!(A[i][id, id], 1e-16)
            fill!(A[i][id, id_nbr], 1e-16)
        end
    end

    return repeat.(A,Nfields,Nfields)
end

# accumulate VhTr*A*Vh on-the-fly
# TODO: speed up via SparseMatrixSimpleBSR
function reduce_jacobian!(A_reduced::SparseMatrixCSC, A::SparseMatrixCSC,
                          md::MeshData, Vh, Nfields=1)

    @unpack FToF = md
    Nfaces,K = size(FToF)
    EToE = @. (FToF.-1) ÷ Nfaces + 1
    Nh,Np = size(Vh)

    ids(e,npts) = (@. (1:npts) + (e-1)*npts)

    for i = 1:Nfields, j = 1:Nfields
        for e = 1:K, f = 1:Nfaces
            enbr = EToE[f,e]
            id,id_nbr   = ids.((e,enbr), Nh)
            idr,idr_nbr = ids.((e,enbr), Np)

            offset(i,num_pts) = (i-1)*num_pts*K

            # convert to matrix for fast NLA
            A_reduced[idr .+ offset(i,Np),idr .+ offset(j,Np)]     .= transpose(Vh)*Matrix(A[id.+offset(i,Nh),id.+offset(j,Nh)])*Vh
            A_reduced[idr .+ offset(i,Np),idr_nbr .+ offset(j,Np)] .= transpose(Vh)*Matrix(A[id.+offset(i,Nh),id_nbr.+offset(j,Nh)])*Vh
        end
    end
end

# =============== for residual evaluation ================

# use ATr for faster col access of sparse CSC matrices
function hadamard_sum(ATr::SparseMatrixCSC,F,u::AbstractArray)
    m, n = size(ATr)
    # rhs = [zeros(n) for i in eachindex(u)]
    rhs = MVector{length(u)}([zeros(n) for i in eachindex(u)]) # probably faster w/StaticArrays?
    hadamard_sum!(rhs,ATr,F,u)
    return rhs
end

# computes ∑ A_ij * F(u_i,u_j) = (A∘F)*1 for flux differencing
function hadamard_sum!(rhs::AbstractArray,ATr::SparseMatrixCSC,F,u::AbstractArray)
    cols = rowvals(ATr)
    vals = nonzeros(ATr)
    m, n = size(ATr)
    for i = 1:n
        ui = getindex.(u,i)
        val_i = zeros(length(u))
        for j in nzrange(ATr, i) # column-major: extracts ith col of ATr = ith row of A
            col = cols[j]
            Aij = vals[j]
            uj = getindex.(u,col)
            val_i += Aij * F(ui,uj)
        end
        setindex!.(rhs,val_i,i)
    end
end

# ============== Constructing DG matrices =================

# flexible but slow function to construct global matrices based on rhs evals
# Np,K = number dofs and elements
# vargs = other args for applyRHS
function build_rhs_matrix(applyRHS,Np,K,vargs...)
    u = zeros(Np,K)
    A = spzeros(Np*K,Np*K)
    for i in eachindex(u)
        u[i] = 1
        r_i = applyRHS(u,vargs...)
        A[:,i] = droptol!(sparse(r_i[:]),1e-12)
        u[i] = 0
    end
    return A
end

# inputs = ref elem and mesh data
# Qrhskew,Qshskew = skew symmetric hybridized SBP operators
# note that md::MeshData needs FToF to also be periodic
function assemble_global_SBP_matrices_2D(rd::RefElemData, md::MeshData,
    Qrhskew, Qshskew, dtol=1e-12)

    @unpack rxJ,sxJ,ryJ,syJ,nxJ,nyJ,mapP,FToF = md
    @unpack Nfaces,wf,wq = rd
    Nh = size(Qrhskew,1)
    NfqNfaces,K = size(nxJ)
    Nfq = convert(Int,NfqNfaces/Nfaces)

    EToE = @. (FToF.-1) ÷ Nfaces + 1
    mapPerm = ((mapP.-1) .% NfqNfaces) .+ 1 # mod out face offset
    fids = length(wq)+1:Nh # last indices correspond to face nodes

    ids(e) = (1:Nh) .+ (e-1)*Nh # block offsets
    Block(e1,e2) = CartesianIndices((ids(e1),ids(e2))) # emulating BlockArrays, but faster
    face_ids = (x,f)->reshape(x,Nfq,Nfaces)[:,f]

    Ax,Ay,Bx,By = ntuple(x->spzeros(Nh*K,Nh*K),4)
    for e = 1:K # loop over elements

        # self-contributions - assume affine for now
        Ax_local = rxJ[1,e]*Qrhskew + sxJ[1,e]*Qshskew
        Ay_local = ryJ[1,e]*Qrhskew + syJ[1,e]*Qshskew
        Ax[Block(e,e)] .= sparse(Ax_local)
        Ay[Block(e,e)] .= sparse(Ay_local)

        # neighboring contributions (ignore self-neighbors)
        for (f,enbr) in enumerate(EToE[:,e])
            if enbr!=e
                fperm = face_ids(reshape(mapPerm[:,e],Nfq,Nfaces),f)
                Axnbr = spdiagm(0 => face_ids(.5*wf.*nxJ[:,e],f))
                Aynbr = spdiagm(0 => face_ids(.5*wf.*nyJ[:,e],f))
                Bx[Block(e,enbr)[face_ids(fids,f),fids[fperm]]] .= Axnbr # TODO: switch to block access pattern
                By[Block(e,enbr)[face_ids(fids,f),fids[fperm]]] .= Aynbr
            end
        end
    end
    return droptol!.((Ax,Ay,Bx,By),dtol)
end

end

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
"""
    chinese_remainder_theorem(array::Array{Int64}, remainder::Array{Int64}))

The Chinese remainder theorem is based on simultaneous congruence. In more detail, if `x`
knows the `remainder` of the Euclidean division of an `array[i]` by several integers, then
`x` can determine the remainder of the division of `array[i]` by the product of these
integers uniquely, following the condition that the divisors are pairwise coprime. For more
information see: [https://en.wikipedia.org/wiki/Chinese_remainder_theorem](https://en.wikipedia.org/wiki/Chinese_remainder_theorem)

# Arguments
- `array::Array{Int64}`: Array with elements.
- `remainder::Array{Int64}`: Number of elements

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> num = [3, 4, 5]
julia> rem = [2, 3, 1]
julia> ClassicAlgorithmsCollections.chinese_remainder_theorem(num, rem)
11
```
"""
function chinese_remainder_theorem(array::Array{Int64}, remainder::Array{Int64})
    x = 1
    n = length(array)
    while true
        i = 1
        while i < n
            if x % array[i] != remainder[i]
                break
            end
            i += 1
        end
        if i == n
            return x
        end
        x += 1
    end
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
using Test
using DataFlowTasks
using DataFlowTasks: memory_overlap
using LinearAlgebra

@testset "Dependencies" begin
    @testset "Arrays and subarrays" begin
        m,n = 20,20
        A   = rand(m,n)
        C   = A
        B   = copy(A)
        @test memory_overlap(A,C) == true
        @test memory_overlap(A,B) == false
        A11 = view(A,1:10,1:10)
        A12 = view(A,11:20,1:10)
        @test memory_overlap(A,A11) == true
        @test memory_overlap(A11,A11) == true
        @test memory_overlap(A11,A12) == false
    end
end


text:
warn("The Datetime.jl package is deprecated, but will remain maintained for the 0.3 release. Please consider using the Dates.jl package instead, which mirrors the Dates module in Base Julia in 0.4")

const FILEPATH = joinpath(dirname(@__FILE__),"tzdata")

#create _leaps and _leaps1
function make_leaps(fname)
  temp = Int[]
  open(joinpath(FILEPATH, fname*".dat")) do fin
    #readdlm was giving an InexactError, so I read the file by hand
    for x in readlines(fin)
      push!(temp, int(x))
    end
  end
  open(joinpath(FILEPATH, fname), "w") do fout
    serialize(fout, temp)
  end
end

for fname in ["_leaps", "_leaps1"]
  make_leaps(fname)
end

#create Zone***DATA
function make_zones(fname)
  #There is probably a better way to read this in
  temp = {}
  open(joinpath(FILEPATH, fname*".dat")) do fin
    for x in readlines(fin)
      ss = split(x, " ")
      push!(temp, ss[1], int(ss[2]), int(ss[3]))
    end
  end

  #reshape temp to matrix for serializing
  mat = transpose(reshape(temp, 3, div(length(temp), 3)))
  open(joinpath(FILEPATH, fname), "w") do fout
    serialize(fout, mat)
  end
end

tzs = map(chomp, collect(open(readlines, joinpath(FILEPATH, "timezonenames"))))

map(make_zones, tzs)


text:
import SymbolicUtils: symtype
struct Parameter{T} end

isparameter(x) = false
isparameter(::Sym{<:Parameter}) = true
isparameter(::Sym{<:FnType{<:Any, <:Parameter}}) = true

SymbolicUtils.symtype(s::Symbolic{Parameter{T}}) where T = T
SymbolicUtils.similarterm(t::Term{T}, f, args) where {T<:Parameter} = Term{T}(f, args)

Base.convert(::Type{Num}, x::Symbolic{Parameter{T}}) where {T<:Number} = Num(x)

"""
$(SIGNATURES)

Define one or more known variables.
"""
macro parameters(xs...)
    esc(_parse_vars(:parameters, Parameter{Real}, xs))
end


text:
# Iterative Loop
for number in 10:20
    println(number)
end

# while loop
number = 1
while number <= 5
    number += 1
    print(number)
end

# Recursive Loop
function fact(n)
    if n == 1
        return(1)
    else
        return n * fact(n-1)
        println(n)
    end
end

fact(5)

array = [5, 10, 15, 20]

newarray = []
for element in array
    if element == 10
        append!(newarray, element)
    elseif element == 20
        append!(newarray, element)
    end
end
println(newarray)

println(array)

# What all can we loop?

# zip
array2 = [30, 40, 50, 60]

for (i, w) in zip(array, array2)
    println(i, w)
end

Hello = "hello world"

function filterstring(input)
    validchars = ['H', 'W', 'l', 'r', 'o', 'd']
    for character in input
        if character in validchars
            print(character)
        end
    end
end

filterstring(Hello)

max_stars_repo_path: 3/Comprehensive Julia Tutorial 3.jl

max_stars_repo_name: Mechachleopteryx/JuliaLessons

max_issues_repo_path: 3/Comprehensive Julia Tutorial 3.jl

max_issues_repo_name: Mechachleopteryx/JuliaLessons

max_forks_repo_path: 3/Comprehensive Julia Tutorial 3.jl

max_forks_repo_name: Mechachleopteryx/JuliaLessons


text:
module TestCode
using PlotlyLocal, Distributions, Colors

data = [Dict(:y => randn(100), :type => "bar"), Dict(:x => 1:100, :y => rand(100))]
layout = Dict(:title => "My Title")
p = plot(data, layout)
data!(p, Dict(:y => cumsum(randn(100)), :x => 1:100))
plot!(p)

data = [Dict(:z => randn(100,100), :type => :surface)]
p = plot(data, Dict())
layout!(p, title = "New Title!")
plot!(p)

end #module


text:
#=

  The BUGS Book
  Example 2.7.1 Surgery (continued): prediction

  Cf ~/jags/bugs_book_2_7_1.R
  """
  model {
    theta ~ dbeta(3, 27) # prior distribution
    Y ~ dbin(theta, 20) # sampling distribution
    P6 <- step(Y - 5.5) # =1 if y >= 6, 0 otherwise
  }

Summary Statistics
  parameters      mean       std   naive_se      mcse         ess      rhat
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64

       theta    0.0987    0.0524     0.0003    0.0011   1799.2552    1.0049
           y    1.9620    1.6057     0.0080    0.0367   1271.8025    1.0081

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

       theta    0.0193    0.0588    0.0900    0.1299    0.2219
           y    0.0000    1.0000    2.0000    3.0000    6.0000

Distributions of variable (num:40)
0.00000 =>   38931  (0.973275)
1.00000 =>    1069  (0.026725)
  4.176741 seconds (17.21 M allocations: 923.018 MiB, 28.68% gc time)

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function bugs_book_2_7_1()

    theta ~ Beta(3,27) # prior distribution
    y ~ Binomial(20, theta) # sampling distribution
    p6 = y >= 6

    return p6
end

model = bugs_book_2_7_1()

num_chains = 4

# chains = sample(model, Prior(), MCMCThreads(), 10_000, num_chains)
chains = sample(model, MH(), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, MH(
#                        # :alpha => Normal(2,sqrt(2)),
#                        # :beta => Normal(2,sqrt(2)),
#                        # :sigma => Gamma(2,2)
#                        ), MCMCThreads(), 40_000, num_chains)

# chains = sample(model, PG(15), MCMCThreads(), 1_000, num_chains)
# chains = sample(model, SMC(10_000), MCMCThreads(), 20_000, num_chains)
# chains = sample(model, IS(), MCMCThreads(), 10_000, num_chains)

# Both HMC and NUTS give the following error:
# ERROR: LoadError: TaskFailedException:
# TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing,Float64,3}
# chains = sample(model, HMC(0.1,5), MCMCThreads(), 10_000, num_chains) # Error
# chains = sample(model, NUTS(0.65), MCMCThreads(), 10_000, num_chains) # Error

display(chains)
# display(plot(chains))

gen = generated_quantities(model, chains)
show_var_dist_pct(gen, 40)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
module TuringCallbacks

using Reexport

using LinearAlgebra
using Logging
using DocStringExtensions
import Turing
@reexport using OnlineStats # used to compute different statistics on-the-fly

using TensorBoardLogger
const TBL = TensorBoardLogger

import DataStructures: DefaultDict

export TensorBoardCallback, DefaultDict, WindowStat, Thin, Skip

include("stats.jl")
include("tensorboardlogger.jl")
include("callbacks/tensorboard.jl")

end


text:
using Luxor

struct Cell <: Drawable
    center::Union{Nothing,Drawable}
    top_left::Union{Nothing,Label}
    top_right::Union{Nothing,Label}
    bottom_left::Union{Nothing,Label}
    bottom_right::Union{Nothing,Label}
    bg_color::Union{Nothing,Color}
    line_thickness::Integer
    line_color::Color
end

function Cell(
    center;
    tl=nothing,
    tr=nothing,
    bl=nothing,
    br=nothing,
    bg_color=nothing,
    line_thickness=1,
    line_color="#000",
    )
    Cell(
        center,
        Label(tl),
        Label(tr),
        Label(bl),
        Label(br),
        Color(bg_color),
        line_thickness,
        Color(line_color),
    )
end

function Cell(c::Cell; kwargs...)
    function copyfield(name)
        if haskey(kwargs, name)
            return kwargs[name]
        end
        getfield(c, name)
    end
    Cell(
        copyfield(:center),
        tl=copyfield(:top_left),
        tr=copyfield(:top_right),
        bl=copyfield(:bottom_left),
        br=copyfield(:bottom_right),
        bg_color=copyfield(:bg_color),
        line_thickness=copyfield(:line_thickness),
        line_color=copyfield(:line_color)
    )
end

hasaftereffect(c::Cell) = c.line_thickness > 0

function draw(c::Cell, pos, cell_size)
    if !isnothing(c.bg_color)
        @layer begin
            sethue(c.bg_color.name)
            box(pos, cell_size, cell_size, :fill)
        end
    end
    if !isnothing(c.center)
        draw(c.center, pos, cell_size; halign=:center, valign=:middle)
    end
end

function drawaftereffect(c::Cell, pos, cell_size)
    line_thickness = c.line_thickness
    line_color = c.line_color
    if line_thickness > 0
        @layer begin
            setline(line_thickness)
            sethue(line_color.name)
            box(pos, cell_size, cell_size, :stroke)
        end
    end
end


text:
# Example script for basic frame connection class usage.

using SMLMFrameConnection
using SMLMSim
using SMLMData
using ImageView

## Simulate some data using SMLMSim.
# Define simulation parameters.
# NOTE: Care must be taken to maintain unit conventions - SMLMSim works in
#       physical units, but here we're maintaining units of pixels and frames.
γ = 1e3 # emission rate, photons/frame
q = [0.0 0.5
    1.0e-3 0.0] # kinetic transmissions (1/frame): [on->on, on->off; off->on, off->off]
fluor = SMLMSim.GenericFluor(; γ = γ, q = q)
pattern = SMLMSim.Point2D()
ρ = 5.0 # emitters / pixel^2
xsize = 8 # pixels
ysize = 8
pixelsize = 1.0 # 1.0 to keep units of pixels
framerate = 1.0 # 1.0 to keep units of frames
nframes = 10000
ndatasets = 1 # for now, must be 1 due to unresolved bug in defineidealFC()

# Generate the simulation
smld_true, smld_model, smld_noisy = SMLMSim.sim(;
    ρ = ρ,
    σ_PSF = 1.3,
    minphotons = 50,
    ndatasets = ndatasets,
    nframes = nframes,
    framerate = framerate,
    pattern = pattern,
    molecule = fluor,
    camera = SMLMSim.IdealCamera(; xpixels = xsize, ypixels = ysize, pixelsize = pixelsize)
)

# Populate some missing fields not added during the simulation.
smld_noisy.bg = zeros(Float64, length(smld_noisy.framenum))
smld_noisy.σ_bg = fill(Inf64, length(smld_noisy.framenum))
smld_noisy.σ_photons = fill(Inf64, length(smld_noisy.framenum))

# Perform frame connection.
smld_connected, smld_preclustered, smld_combined, params = SMLMFrameConnection.frameconnect(smld_noisy;
    nnearestclusters = 2, nsigmadev = 5.0,
    maxframegap = 5, nmaxnn = 2)

## Make some circle images of the results (circle radii indicate localization
## precision).
# Plot the combined results overlain with the original data:
#   Original data shown in magenta, combined results in green.
mag = 100.0 # image magnification
circleim_original = SMLMData.makecircleim(smld_noisy, mag)
circleim_combined = SMLMData.makecircleim(smld_combined, mag)
ImageView.imshow(ImageView.RGB.(circleim_original, circleim_combined, circleim_original))

# Plot the "ideal" result (all localizations of same emitter w/in 5 frames are
# combined) overlain with the LAP-FC result:
#   Ideal result shown in magenta, LAP-FC results in green
#   -> white circles indicate "ideal" performance of LAP-FC
smld_idealconnected, smld_idealcombined =
    SMLMFrameConnection.defineidealFC(smld_noisy; maxframegap = 5)
circleim_ideal = SMLMData.makecircleim(smld_idealcombined, mag)
ImageView.imshow(ImageView.RGB.(circleim_ideal, circleim_combined, circleim_ideal))


text:
using AbstractGPs, KernelFunctions
using TemporalGPs
using SpecialFunctions
using StaticArrays
using LinearAlgebra
using Zygote

#----------------------------------------------------------------------------
# periodic kernel
struct ApproxPeriodicKernel{J,T} <: KernelFunctions.SimpleKernel
    l::T
end

# assuming length( k.r ) == 1
ApproxPeriodicKernel(k::PeriodicKernel{T}; degree_of_approx = 3) where {T} = ApproxPeriodicKernel{degree_of_approx,T}( 4*only(k.r) )
ApproxPeriodicKernel(; r=1/4, degree_of_approx = 3 ) = ApproxPeriodicKernel(PeriodicKernel(; r=[r]); degree_of_approx=degree_of_approx)
# k = ApproxPeriodicKernel()

function TemporalGPs.to_sde(k::ApproxPeriodicKernel{J,T}, s::SArrayStorage{T}) where {J,T<:Real}
    D  = 2*(J+1)
    w0 = 2π
    F  = kron( Diagonal(0:J), [0 -w0; w0 0])
    H  = kron( ones(1, J+1),[1 0])
    Fs = SMatrix{D,D,T}(F)
    q  = T(0)
    Hs = SVector{D,T}(H)#SMatrix{1,D,T}(H)
    return Fs, q, Hs
end

function TemporalGPs.stationary_distribution(k::ApproxPeriodicKernel{J,T}, s::SArrayStorage{T}) where {J,T<:Real}
    D = 2*(J+1)
    q = [ j>0 ? 2*besseli(j, 1/k.l^2)*exp(-k.l^2) : besseli(0, 1/k.l^2)*exp(-k.l^2) for j in 0:J]
    P = kron(Diagonal(q), Diagonal([1, 1]))
    m = SVector{D, T}( zeros(T,D) )
    Ps = SMatrix{D,D,T}(P)
    return TemporalGPs.Gaussian(m, Ps)
end

Zygote.@adjoint function to_sde(k::ApproxPeriodicKernel, storage_type)
    return to_sde(k, storage_type), Δ->(nothing, nothing)
end

Zygote.@adjoint function TemporalGPs.stationary_distribution(k::ApproxPeriodicKernel, storage_type)
    return TemporalGPs.stationary_distribution(k, storage_type), Δ->(nothing, nothing)
end

#------------------------------------------------------------------------
# Kernel sum (deep)
KernelFunctions.KernelSum(kernel::Kernel) = kernel

function TemporalGPs.lgssm_components(k::KernelSum, ts::AbstractVector, storage_type::TemporalGPs.StorageType)
    As_l, as_l, Qs_l, emission_proj_l, x0_l = TemporalGPs.lgssm_components(k.kernels[1], ts, storage_type)
    As_r, as_r, Qs_r, emission_proj_r, x0_r = TemporalGPs.lgssm_components( KernelSum( k.kernels[2:end]... ), ts, storage_type)

    As = map(TemporalGPs.blk_diag, As_l, As_r)
    as = map(vcat, as_l, as_r)
    Qs = map(TemporalGPs.blk_diag, Qs_l, Qs_r)
    emission_projections = TemporalGPs._sum_emission_projections(emission_proj_l, emission_proj_r)
    x0 = TemporalGPs.Gaussian(vcat(x0_l.m, x0_r.m), TemporalGPs.blk_diag(x0_l.P, x0_r.P))

    return As, as, Qs, emission_projections, x0
end

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
"""
**Student-T likelihood**

Student-t likelihood for regression: ``\\frac{\\Gamma((\\nu+1)/2)}{\\sqrt{\\nu\\pi}\\sigma\\Gamma(\\nu/2)}\\left(1+(y-f)^2/(\\sigma^2\\nu)\\right)^{(-(\\nu+1)/2)}``
see [wiki page](https://en.wikipedia.org/wiki/Student%27s_t-distribution)

```julia
StudentTLikelihood(ν::T,σ::Real=one(T)) #ν is the number of degrees of freedom
#σ is the variance for local scale of the data.
```

---

For the analytical solution, it is augmented via:
```math
p(y|f,\\omega) = \\mathcal{N}(y|f,\\sigma^2\\omega)
```
Where ``\\omega \\sim \\mathcal{IG}(\\frac{\\nu}{2},\\frac{\\nu}{2})`` where ``\\mathcal{IG}`` is the inverse gamma distribution
See paper [Robust Gaussian Process Regression with a Student-t Likelihood](http://www.jmlr.org/papers/volume12/jylanki11a/jylanki11a.pdf)
"""
struct StudentTLikelihood{T<:Real} <: RegressionLikelihood{T}
    ν::T
    α::T
    σ::T
    β::LatentArray{Vector{T}}
    θ::LatentArray{Vector{T}}
    function StudentTLikelihood{T}(ν::T,σ::T=one(T)) where {T<:Real}
        new{T}(ν,(ν+one(T))/2.0,σ)
    end
    function StudentTLikelihood{T}(ν::T,σ::T,β::AbstractVector{<:AbstractVector{T}},θ::AbstractVector{<:AbstractVector{T}}) where {T<:Real}
        new{T}(ν,(ν+one(T))/2.0,σ,β,θ)
    end
end

function StudentTLikelihood(ν::T,σ::T=one(T)) where {T<:Real}
    StudentTLikelihood{T}(ν,σ)
end

function init_likelihood(likelihood::StudentTLikelihood{T},inference::Inference{T},nLatent::Int,nSamplesUsed::Int,nFeatures::Int) where T
    if inference isa AnalyticVI || inference isa GibbsSampling
        StudentTLikelihood{T}(
        likelihood.ν,likelihood.σ,
        [abs2.(T.(rand(T,nSamplesUsed))) for _ in 1:nLatent],
        [zeros(T,nSamplesUsed) for _ in 1:nLatent])
    else
        StudentTLikelihood{T}(likelihood.ν,likelihood.σ)
    end
end

function pdf(l::StudentTLikelihood,y::Real,f::Real)
    tdistpdf(l.ν,(y-f)/l.σ)
end

function Base.show(io::IO,model::StudentTLikelihood{T}) where T
    print(io,"Student-t likelihood")
end

function compute_proba(l::StudentTLikelihood{T},μ::AbstractVector{T},σ²::AbstractVector{T}) where {T<:Real}
    return μ,max.(σ²,0.0).+0.5*l.ν*l.σ^2/(0.5*l.ν-1)
end

## Local Updates ##

function local_updates!(model::VGP{T,<:StudentTLikelihood,<:AnalyticVI}) where {T}
    model.likelihood.β .= broadcast((Σ,μ,y)->0.5*(Σ+abs2.(μ-y).+model.likelihood.σ^2*model.likelihood.ν),diag.(model.Σ),model.μ,model.y)
    model.likelihood.θ .= broadcast(β->model.likelihood.α./β,model.likelihood.β)
end

function local_updates!(model::SVGP{T,<:StudentTLikelihood,<:AnalyticVI}) where {T}
    model.likelihood.β .= broadcast((K̃,κ,Σ,μ,y)->0.5*(K̃ + opt_diag(κ*Σ,κ) + abs2.(κ*μ-y).+model.likelihood.σ^2*model.likelihood.ν),model.K̃,model.κ,model.Σ,model.μ,model.inference.y)
    model.likelihood.θ .= broadcast(β->model.likelihood.α./β,model.likelihood.β)
end

function sample_local!(model::VGP{T,<:StudentTLikelihood,<:GibbsSampling}) where {T}
    model.likelihood.β .= broadcast((μ::AbstractVector{<:Real},y)->rand.(InverseGamma.(model.likelihood.α,0.5*(abs2.(μ-y).+model.likelihood.σ^2*model.likelihood.ν))),model.μ,model.y)
    model.likelihood.θ .= broadcast(β->1.0./β,model.likelihood.β)
    return nothing
end

## Global Gradients ##

@inline ∇E_μ(model::AbstractGP{T,<:StudentTLikelihood,<:GibbsorVI}) where {T} = hadamard.(model.likelihood.θ,model.inference.y)
@inline ∇E_μ(model::AbstractGP{T,<:StudentTLikelihood,<:GibbsorVI},i::Int) where {T} = hadamard(model.likelihood.θ[i],model.inference.y[i])

@inline ∇E_Σ(model::AbstractGP{T,<:StudentTLikelihood,<:GibbsorVI}) where {T} = 0.5.*model.likelihood.θ
@inline ∇E_Σ(model::AbstractGP{T,<:StudentTLikelihood,<:GibbsorVI},i::Int) where {T} = 0.5.*model.likelihood.θ[i]

## ELBO Section ##

function ELBO(model::AbstractGP{T,<:StudentTLikelihood,<:AnalyticVI}) where {T}
    return expecLogLikelihood(model) - InverseGammaKL(model) - GaussianKL(model)
end

function expecLogLikelihood(model::VGP{T,<:StudentTLikelihood,<:AnalyticVI}) where {T}
    # model.likelihood.β .= broadcast((Σ,μ,y)->0.5*(Σ+abs2.(μ-y).+model.likelihood.σ^2*model.likelihood.ν),diag.(model.Σ),model.μ,model.y)
    # model.likelihood.θ .= broadcast(β->model.likelihood.α./β,model.likelihood.β)
    tot = -0.5*model.nLatent*model.nSample*(log(twoπ*model.likelihood.σ^2))
    tot += -sum(broadcast((α,β)->sum(log.(β).-digamma(α)), model.likelihood.α,model.likelihood.β))
    tot += -0.5.*sum(broadcast((θ,Σ,μ,y)->dot(θ,Σ)+dot(θ,abs2.(μ))-2.0*dot(θ,μ.*y)+dot(θ,abs2.(y)),model.likelihood.θ,diag.(model.Σ),model.μ,model.y))
    return tot
end

function expecLogLikelihood(model::SVGP{T,<:StudentTLikelihood,<:AnalyticVI}) where {T}
    tot = -0.5*model.nLatent*model.inference.nSamplesUsed*(log(twoπ*model.likelihood.σ^2))
    tot += -sum(broadcast(β->sum(digamma(model.likelihood.α).-log.(β)),model.likelihood.β))
    tot += -0.5.*sum(broadcast((θ,K̃,κ,Σ,κμ,y)->dot(θ,(K̃+opt_diag(κ*Σ,κ)+abs2.(κμ)-2.0*(κμ).*y+abs2.(y))),model.likelihood.θ,model.K̃,model.κ,model.Σ,model.κ.*model.μ,model.inference.y))
    return model.inference.ρ*tot
end

function InverseGammaKL(model::AbstractGP{T,<:StudentTLikelihood}) where {T}
    α_p = model.likelihood.ν/2; β_p= α_p*model.likelihood.σ^2
    model.inference.ρ*sum(broadcast(InverseGammaKL,model.likelihood.α,model.likelihood.β,α_p,β_p))
end

## PDF and Log PDF Gradients ## (verified gradients)

function grad_log_pdf(l::StudentTLikelihood{T},y::Real,f::Real) where {T<:Real}
    (one(T)+l.ν) * (y-f) / ((f-y)^2 + l.σ^2*l.ν)
end

function hessian_log_pdf(l::StudentTLikelihood{T},y::Real,f::Real) where {T<:Real}
    v = l.ν * l.σ^2; Δ² = (f-y)^2
    (one(T)+l.ν) * (-v + Δ²) / (v+Δ²)^2
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "WaveFD"
version = v"0.3.0"

# Collection of sources required to build AzStorage
sources = [
    GitSource(
        "https://github.com/ChevronETC/WaveFD.jl.git",
        "d0aa5dda65904317b4c7b30e6fd90a90b6c368c7"
    )
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/WaveFD.jl/src

CXXFLAGS="-funroll-loops"
if [[ "${target}" == i686-linux-gnu || "${target}" == x86_64-linux-gnu || "${target}" == powerpc6rle-linux-gnu ]]; then
    CXXFLAGS+=" -D__FUNCTION_CLONES__"
fi

echo "target=$target, CXXFLAGS=$CXXFLAGS"

cmake . -DCMAKE_CXX_FLAGS="${CXXFLAGS}" -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN%.*}_gcc.cmake -DCMAKE_BUILD_TYPE=Release
make -j${nproc} VERBOSE=1
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libprop2DAcoIsoDenQ_DEO2_FDTD", :libprop2DAcoIsoDenQ_DEO2_FDTD),
    LibraryProduct("libprop2DAcoVTIDenQ_DEO2_FDTD", :libprop2DAcoVTIDenQ_DEO2_FDTD),
    LibraryProduct("libprop2DAcoTTIDenQ_DEO2_FDTD", :libprop2DAcoTTIDenQ_DEO2_FDTD),
    LibraryProduct("libprop3DAcoIsoDenQ_DEO2_FDTD", :libprop3DAcoIsoDenQ_DEO2_FDTD),
    LibraryProduct("libprop3DAcoVTIDenQ_DEO2_FDTD", :libprop3DAcoVTIDenQ_DEO2_FDTD),
    LibraryProduct("libprop3DAcoTTIDenQ_DEO2_FDTD", :libprop3DAcoTTIDenQ_DEO2_FDTD),
    LibraryProduct("libillumination", :libillumination),
    LibraryProduct("libspacetime", :libspacetime)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("FFTW_jll")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"9", allow_unsafe_flags=true)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@testset "Utilities" begin
    @testset "Impute.dim" begin
        X = rand(10, 5)
        KA = KeyedArray(X; A=1:10, B=collect("abcde"))

        @test Impute.dim(X, 1) == Impute.dim(X, :rows) == Impute.dim(KA, :A)
        @test first(eachslice(X, dims=1)) == first(eachslice(KA, dims=1)) == first(eachslice(KA, dims=:A))

        @test Impute.dim(X, 2) == Impute.dim(X, :cols) == Impute.dim(KA, :B)
        @test first(eachslice(X, dims=2)) == first(eachslice(KA, dims=2)) == first(eachslice(KA, dims=:B))
    end
end


text:
"""
    kernelmatrix!(K::AbstractMatrix, κ::Kernel, x::AbstractVector)
    kernelmatrix!(K::AbstractMatrix, κ::Kernel, x::AbstractVector, y::AbstractVector)

In-place version of [`kernelmatrix`](@ref) where pre-allocated matrix `K` will be
overwritten with the kernel matrix.

    kernelmatrix!(K::AbstractMatrix, κ::Kernel, X::AbstractMatrix; obsdim::Integer=2)
    kernelmatrix!(
        K::AbstractMatrix,
        κ::Kernel,
        X::AbstractMatrix,
        Y::AbstractMatrix;
        obsdim::Integer=2,
    )

Equivalent to `kernelmatrix!(K, κ, ColVecs(X))` and
`kernelmatrix(K, κ, ColVecs(X), ColVecs(Y))` respectively.
Set `obsdim=1` to get `RowVecs`.

See also: [`ColVecs`](@ref), [`RowVecs`](@ref)
"""
kernelmatrix!

"""
    kernelmatrix(κ::Kernel, x::AbstractVector)

Compute the kernel `κ` for each pair of inputs in `x`.
Returns a matrix of size `(length(x), length(x))` satisfying
`kernelmatrix(κ, x)[p, q] == κ(x[p], x[q])`.

    kernelmatrix(κ::Kernel, x::AbstractVector, y::AbstractVector)

Compute the kernel `κ` for each pair of inputs in `x` and `y`.
Returns a matrix of size `(length(x), length(y))` satisfying
`kernelmatrix(κ, x, y)[p, q] == κ(x[p], y[q])`.

    kernelmatrix(κ::Kernel, X::AbstractMatrix; obsdim::Int=2)
    kernelmatrix(κ::Kernel, X::AbstractMatrix, Y::AbstractMatrix; obsdim::Int=2)

Equivalent to `kernelmatrix(κ, ColVecs(X))` and `kernelmatrix(κ, ColVecs(X), ColVecs(Y))`
respectively.
Set `obsdim=1` to get `RowVecs`.

See also: [`ColVecs`](@ref), [`RowVecs`](@ref)
"""
kernelmatrix

"""
    kernelmatrix_diag!(K::AbstractVector, κ::Kernel, x::AbstractVector)
    kernelmatrix_diag!(K::AbstractVector, κ::Kernel, x::AbstractVector, y::AbstractVector)

In place version of [`kernelmatrix_diag`](@ref).

    kernelmatrix_diag!(K::AbstractVector, κ::Kernel, X::AbstractMatrix; obsdim::Int=2)
    kernelmatrix_diag!(
        K::AbstractVector,
        κ::Kernel,
        X::AbstractMatrix,
        Y::AbstractMatrix;
        obsdim::Int=2,
    )

Equivalent to `kernelmatrix_diag!(K, κ, ColVecs(X))` and
`kernelmatrix_diag!(K, κ, ColVecs(X), ColVecs(Y))` respectively.
Set `obsdim=1` to get `RowVecs`.

See also: [`ColVecs`](@ref), [`RowVecs`](@ref)
"""
kernelmatrix_diag!

"""
    kernelmatrix_diag(κ::Kernel, x::AbstractVector)

Compute the diagonal of `kernelmatrix(κ, x)` efficiently.

    kernelmatrix_diag(κ::Kernel, x::AbstractVector, y::AbstractVector)

Compute the diagonal of `kernelmatrix(κ, x, y)` efficiently.
Requires that `x` and `y` are the same length.

    kernelmatrix_diag(κ::Kernel, X::AbstractMatrix; obsdim::Int=2)
    kernelmatrix_diag(κ::Kernel, X::AbstractMatrix, Y::AbstractMatrix; obsdim::Int=2)

Equivalent to `kernelmatrix_diag(κ, ColVecs(X))` and
`kernelmatrix_diag(κ, ColVecs(X), ColVecs(Y))` respectively.

See also: [`ColVecs`](@ref), [`RowVecs`](@ref)
"""
kernelmatrix_diag

#
# Kernel implementations. Generic fallbacks that depend only on kernel evaluation.
#

kernelmatrix!(K::AbstractMatrix, κ::Kernel, x::AbstractVector) = kernelmatrix!(K, κ, x, x)

function kernelmatrix!(K::AbstractMatrix, κ::Kernel, x::AbstractVector, y::AbstractVector)
    validate_inplace_dims(K, x, y)
    K .= κ.(x, permutedims(y))
    return K
end

kernelmatrix(κ::Kernel, x::AbstractVector) = kernelmatrix(κ, x, x)

function kernelmatrix(κ::Kernel, x::AbstractVector, y::AbstractVector)
    validate_inputs(x, y)
    return κ.(x, permutedims(y))
end

function kernelmatrix_diag!(K::AbstractVector, κ::Kernel, x::AbstractVector)
    validate_inplace_dims(K, x)
    return map!(x -> κ(x, x), K, x)
end

function kernelmatrix_diag!(
    K::AbstractVector, κ::Kernel, x::AbstractVector, y::AbstractVector
)
    return map!(κ, K, x, y)
end

kernelmatrix_diag(κ::Kernel, x::AbstractVector) = map(x -> κ(x, x), x)

kernelmatrix_diag(κ::Kernel, x::AbstractVector, y::AbstractVector) = map(κ, x, y)

#
# SimpleKernel optimisations.
#

function kernelmatrix!(K::AbstractMatrix, κ::SimpleKernel, x::AbstractVector)
    validate_inplace_dims(K, x)
    pairwise!(K, metric(κ), x)
    return map!(x -> kappa(κ, x), K, K)
end

function kernelmatrix!(
    K::AbstractMatrix, κ::SimpleKernel, x::AbstractVector, y::AbstractVector
)
    validate_inplace_dims(K, x, y)
    pairwise!(K, metric(κ), x, y)
    return map!(x -> kappa(κ, x), K, K)
end

function kernelmatrix(κ::SimpleKernel, x::AbstractVector)
    return map(x -> kappa(κ, x), pairwise(metric(κ), x))
end

function kernelmatrix(κ::SimpleKernel, x::AbstractVector, y::AbstractVector)
    validate_inputs(x, y)
    return map(x -> kappa(κ, x), pairwise(metric(κ), x, y))
end

function kernelmatrix_diag(κ::SimpleKernel, x::AbstractVector)
    return map(x -> kappa(κ, x), colwise(metric(κ), x))
end

function kernelmatrix_diag(κ::SimpleKernel, x::AbstractVector, y::AbstractVector)
    return map(x -> kappa(κ, x), colwise(metric(κ), x, y))
end

#
# Wrapper methods for AbstractMatrix inputs to maintain obsdim interface.
#

const defaultobs = 2

function kernelmatrix!(
    K::AbstractMatrix, κ::Kernel, X::AbstractMatrix; obsdim::Int=defaultobs
)
    return kernelmatrix!(K, κ, vec_of_vecs(X; obsdim=obsdim))
end

function kernelmatrix!(
    K::AbstractMatrix,
    κ::Kernel,
    X::AbstractMatrix,
    Y::AbstractMatrix;
    obsdim::Int=defaultobs,
)
    return kernelmatrix!(K, κ, vec_of_vecs(X; obsdim=obsdim), vec_of_vecs(Y; obsdim=obsdim))
end

function kernelmatrix(κ::Kernel, X::AbstractMatrix; obsdim::Int=defaultobs)
    return kernelmatrix(κ, vec_of_vecs(X; obsdim=obsdim))
end

function kernelmatrix(κ::Kernel, X::AbstractMatrix, Y::AbstractMatrix; obsdim=defaultobs)
    return kernelmatrix(κ, vec_of_vecs(X; obsdim=obsdim), vec_of_vecs(Y; obsdim=obsdim))
end

function kernelmatrix_diag!(
    K::AbstractVector, κ::Kernel, X::AbstractMatrix; obsdim::Int=defaultobs
)
    return kernelmatrix_diag!(K, κ, vec_of_vecs(X; obsdim=obsdim))
end

function kernelmatrix_diag!(
    K::AbstractVector,
    κ::Kernel,
    X::AbstractMatrix,
    Y::AbstractMatrix;
    obsdim::Int=defaultobs,
)
    return kernelmatrix_diag!(
        K, κ, vec_of_vecs(X; obsdim=obsdim), vec_of_vecs(Y; obsdim=obsdim)
    )
end

function kernelmatrix_diag(κ::Kernel, X::AbstractMatrix; obsdim::Int=defaultobs)
    return kernelmatrix_diag(κ, vec_of_vecs(X; obsdim=obsdim))
end

function kernelmatrix_diag(
    κ::Kernel, X::AbstractMatrix, Y::AbstractMatrix; obsdim::Int=defaultobs
)
    return kernelmatrix_diag(
        κ, vec_of_vecs(X; obsdim=obsdim), vec_of_vecs(Y; obsdim=obsdim)
    )
end


text:
using ADCME, LinearAlgebra, PyCall
using DelimitedFiles
using PyPlot

# read data
elem = readdlm("meshdata/elem.txt", Int64)
node = readdlm("meshdata/nodes.txt")
dof = readdlm("meshdata/dof.txt", Int64)[:]
elem_ = constant(elem)
ne = size(elem,1)
nv = size(node, 1)

# precompute
localcoef = zeros(ne, 3, 3)
areas = zeros(ne)
for e = 1:ne
    el = elem[e,:]
    x1, y1 = node[el[1],:]
    x2, y2 = node[el[2],:]
    x3, y3 = node[el[3],:]
    A = [x1 y1 1.0; x2 y2 1.0; x3 y3 1.0]
    localcoef[e,:,:] = inv(A)
    areas[e] = 0.5*abs(det(A))
end

# compute right hand side using midpoint rule
rhs = zeros(nv)
for i = 1:ne
    el = elem[i,:]
    rhs[el] .+= 4*areas[i]/3
end

areas = constant(areas)
localcoef = constant(localcoef)
D = constant(diagm(0=>ones(2)))
# D = Variable(2.0) .* [1.0 0.0;0.0 1.0]

function body(i, tai, taj, tav)
    el = elem_[i-1]
    a = areas[i-1]
    L = localcoef[i-1]
    LocalStiff = Array{PyObject}(undef, 3, 3)
    for i = 1:3
        for j = 1:3
            LocalStiff[i,j] = a*[L[1,i] L[2,i]]*D*[L[1,j];L[2,j]]|>squeeze
        end
    end
    ii = reshape([el el el], (-1,))
    jj = reshape([el;el;el], (-1,))
    tai = write(tai, i, ii)
    taj = write(taj, i, jj)
    # op = tf.print(el)
    # i = bind(i, op)
    tav = write(tav, i, vcat(LocalStiff[:]...))
    return i+1, tai, taj, tav
end

i = constant(2, dtype=Int32)
tai = TensorArray(ne+1, dtype=Int64)
taj = TensorArray(ne+1, dtype=Int64)
tav = TensorArray(ne+1)
tai = write(tai, 1, constant(ones(Int64,9)))
taj = write(taj, 1, constant(ones(Int64,9)))
tav = write(tav, 1, constant(zeros(9)))
_, ii, jj, vv = while_loop((i, tas...)->i<=ne+1, body, [i, tai, taj, tav])
ii = reshape(stack(ii),(-1,)); jj = reshape(stack(jj),(-1,)); vv = reshape(stack(vv),(-1,))

A = SparseTensor(ii, jj, vv, nv, nv)

ndof = [x for x in setdiff(Set(1:nv), Set(dof))]
A = scatter_update(A, dof, ndof, spzero(length(dof), length(ndof)))
A = scatter_update(A, ndof, dof, spzero(length(ndof), length(dof)))
A = scatter_update(A, dof, dof, spdiag(length(dof)))
rhs[dof] .= 0.0
sol = A\rhs

# loss = sum((sol - (@. 1-node[:,1]^2-node[:,2]^2))^2)
sess = Session(); init(sess)
S = run(sess, sol)
close("all")
scatter3D(node[:,1], node[:,2], S, marker="^", label = "FEM")
scatter3D(node[:,1], node[:,2], (@. 1-node[:,1]^2-node[:,2]^2), marker = "+", label = "Exact")
legend()


text:
# This script contains commands to calculate basic summary
# statistics on the dataset

import MLJ: schema, std, mean, median, coerce, coerce!, scitype
using DataFrames
using CSV
using UrlDownload
using PyPlot

# Load files - this assumes that the repository has been cloned locally,
# that jurisdiction-level files have been concatenated (using the Python script
# in `extract.py`) and saved locally.

path_nat = "/Users/GD/Desktop/wcpdb_nat.csv"
path_subnat = "/Users/GD/Desktop/wcpdb_subnat.csv"

data_nat = CSV.read(path_nat, DataFrame)
data_subnat = CSV.read(path_subnat, DataFrame)

data_nat_2018 = data_nat[data_nat[:Year] .== 2018, :]
data_subnat_2018 = data_subnat[data_subnat[:Year] .== 2018, :]

dropmissing!(data_nat_2018, [:Tax_dummy, :ETS_dummy])
dropmissing!(data_subnat_2018, [:Tax_dummy, :ETS_dummy])

gdf_nat = groupby(data_nat_2018, [:Jurisdiction, :Year])
gdf_subnat = groupby(data_subnat_2018, [:Jurisdiction, :Year])

# Total number of schemes
temp = combine(gdf_nat, :Tax_dummy => sum => :Sum_tax, :ETS_dummy => sum => :Sum_ets)
temp = combine(gdf_subnat, :Tax_dummy => sum => :Sum_tax, :ETS_dummy => sum => :Sum_ets)

temp.Sum_tax_bin = [x > 0 ? 1 : 0 for x in temp.Sum_tax]
temp.Sum_tax_ets = [x > 0 ? 1 : 0 for x in temp.Sum_ets]

gtemp = groupby(temp, :Year)
gtemp = combine(gtemp, :Sum_tax_bin => sum => :tax_tot,
                :Sum_tax_ets => sum => :ets_tot)

# Sector level stats
sectors = r"ABFLOW003|ABFLOW012|ABFLOW028|ABFLOW029|ABFLOW034|ABFLOW035"

data_2018 = vcat(data_nat_2018, data_subnat_2018)
dropmissing!(data_2018, :IEA_CODE)
data_2018_sec = data_2018[occursin.(sectors, data_2018.IEA_CODE), :]

gdata_2018_sec = groupby(data_2018_sec, [:Jurisdiction, :Year, :IEA_CODE])

gsec_temp = combine(gdata_2018_sec, :Tax_dummy => sum => :Sum_tax, :ETS_dummy => sum => :Sum_ets)

gsec_temp.Sum_tax_bin = [x > 0 ? 1 : 0 for x in gsec_temp.Sum_tax]
gsec_temp.Sum_tax_ets = [x > 0 ? 1 : 0 for x in gsec_temp.Sum_ets]

ggsec_temp = groupby(gsec_temp, [:Year, :IEA_CODE])
ggsec_temp = combine(ggsec_temp, :Sum_tax_bin => sum => :tax_tot,
                :Sum_tax_ets => sum => :ets_tot)

# Sector-level carbon prices
path_prices = "/Users/GD/Sync/PE/Data/ECP/ECP_calculation/Price/Sector_level/EW_sector_prices.csv"
prices = CSV.read(path_prices, DataFrame)

prices_2017 = prices[prices[:, :Year] .== 2017, :]
prices_2017 = prices_2017[occursin.(sectors, prices_2017.Flow), :]

gprices_2017 = groupby(prices_2017, [:Year, :Flow])

sumstat = combine(gprices_2017, :Total_price_sector_2019USD => describe)

max_stars_repo_name: g-dolphin/WorldCarbonPricingDatabase

max_issues_repo_name: g-dolphin/WorldCarbonPricingDatabase

max_forks_repo_name: g-dolphin/WorldCarbonPricingDatabase


text:
# joint neighbor probability for symmetric SNE
# function joint_neighbor(x::VecOfVec{<:Number}, kernel = cauchy, distance = euclidean)
function joint_neighbor(X::AbstractMatrix, kernel = cauchy, distance = euclidean)
    n = size(X, 2)
    joint_neighbor!(zeros(n, n), X, kernel, distance)
end
function joint_neighbor!(K::AbstractMatrix, X::AbstractMatrix,
                         kernel = cauchy, distance = euclidean)
    unnormalized_joint_neighbor!(K, X, kernel, distance)
    K ./= sum(K) # entire matrix is normalized probability distribution now
end

# sum over all elements of Q via matrix-vector products to work with fast transforms
function normalization_constant(Q::AbstractMatOrFac)
    n, m = size(Q)
    x, y = ones(n), ones(m)
    return dot(x, Q*y)
end

function unnormalized_joint_neighbor!(K::AbstractMatrix, X::AbstractMatrix,
                         kernel = cauchy, distance = euclidean)
    x = [c for c in eachcol(X)]
    xt = permutedims(x)
    @. K = kernel(distance(x, xt)) # IDEA: lazy: LazyMatrixSum(gramian(kernel, x), kernel(0)*I(length(x)))
    for i in diagind(K) # "Because we are only interested in modeling pairwise similarities, we set the value of pi|i to zero."
        K[i] = 0
    end
end

# lazy representation
function unnormalized_joint_neighbor(X::AbstractMatrix, kernel = cauchy, distance = euclidean, dofast::Bool = false)
    d, n = size(X)
    k = Cauchy()
    k² = k^2
    Q = gramian(k, X)
    Q² = gramian(k², X) # squared kernel
    D = -1I(n) # to subtract diagonal
    Q = Matrix(Q) + D # IDEA: Lazy
    Q² = Matrix(Q²) + D
    return Q, Q²
end


text:
using UnivariateFunctions
using StochasticIntegrals
using Dates
using LinearAlgebra
using Distributions: Normal, MersenneTwister, quantile
using Statistics: var, mean, cov
using Sobol
using Random
const tol = 10*eps()

USD_mean_reversion  = 0.15
USD_vol = 0.01
GBP_mean_reversion = 0.1
GBP_vol = 0.001
GBP_FX_Vol = 0.2

today = Date(2016,12,1)
date_2020 = Date(2020,12,1)
later_date = Date(2035,12,1)
later_later_date = Date(2080,12,1)

USD_hw_a_curve = PE_Function(USD_vol, USD_mean_reversion, today,0)
USD_hw_aB_curve = USD_hw_a_curve * (1/USD_mean_reversion) * (1 - PE_Function(1.0, -USD_mean_reversion, today,0))
GBP_hw_a_curve = PE_Function(GBP_vol, GBP_mean_reversion, today,0)
GBP_hw_aB_curve = GBP_hw_a_curve * (1/GBP_mean_reversion) * (1 - PE_Function(1.0, -GBP_mean_reversion, today,0))

brownian_corr_matrix = Symmetric([1.0 0.75 0.5 0.0;
                                  0.0 1.0 0.5 0.25;
                                  0.0 0.0 1.0 0.25;
                                  0.0 0.0 0.0 1.0])
brownian_ids = [:USD_IR, :GBP_IR, :GBP_FX, :BARC]
USD_IR_a_ito  = ItoIntegral(:USD_IR, USD_hw_a_curve)
USD_IR_aB_ito = ItoIntegral(:USD_IR, USD_hw_aB_curve)
GBP_IR_a_ito  = ItoIntegral(:GBP_IR, GBP_hw_a_curve)
GBP_IR_aB_ito = ItoIntegral(:GBP_IR, GBP_hw_aB_curve)
GBP_FX_ito    = ItoIntegral(:GBP_FX, GBP_FX_Vol)
ito_integrals = Dict([:USD_IR_a, :USD_IR_aB, :GBP_IR_a, :GBP_IR_aB, :GBP_FX] .=> [USD_IR_a_ito, USD_IR_aB_ito, GBP_IR_a_ito, GBP_IR_aB_ito, GBP_FX_ito])

ito_set_ = ItoSet(brownian_corr_matrix, brownian_ids, ito_integrals)
# The next ito integral should have constant vol
abs(get_volatility(ito_set_,  :GBP_FX, Date(2020,1,1)) - GBP_FX_Vol) < tol
abs(get_volatility(ito_set_,  :GBP_FX, Date(2022,1,1)) - GBP_FX_Vol) < tol
# The next ito integral has changing vol
abs(get_volatility(ito_set_,  :USD_IR_a, Date(2020,1,1)) - get_volatility(ito_set_,  :USD_IR_a, Date(2022,1,1)) ) > 0.005
abs(get_volatility(ito_set_,  :USD_IR_a, Date(2020,1,1)) - evaluate(USD_hw_a_curve, Date(2020,1,1)) ) < tol

covar = CovarianceAtDate(ito_set_, years_from_global_base(today), years_from_global_base(date_2020))
abs(covar.covariance_[5,5] - GBP_FX_Vol^2 * (years_from_global_base(date_2020) - years_from_global_base(today))) < tol
abs(covar.covariance_[1,1] - get_variance(USD_IR_a_ito, years_from_global_base(today), years_from_global_base(date_2020))) < tol

cov_date = CovarianceAtDate(ito_set_, today, later_date)
abs(cov_date.covariance_[5,5] - GBP_FX_Vol^2 * (years_from_global_base(later_date) - years_from_global_base(today))) < tol

abs(get_volatility(ito_set_,  :USD_IR_a, Date(2020,1,1)) - get_volatility(cov_date,  :USD_IR_a, Date(2020,1,1)) ) < tol
abs(cov_date.covariance_[1,3] - get_covariance(cov_date, :USD_IR_a, :GBP_IR_a)) < tol
abs(get_variance(cov_date, :USD_IR_a) - get_covariance(cov_date, :USD_IR_a, :USD_IR_a)) < tol

# Test correlation.
abs(get_correlation(cov_date, :USD_IR_a, :GBP_IR_a) - (get_covariance(cov_date, :USD_IR_a, :GBP_IR_a)/sqrt(get_variance(cov_date, :USD_IR_a) * get_variance(cov_date, :GBP_IR_a)))) < tol
abs(get_correlation(cov_date, :GBP_IR_aB, :GBP_IR_a) - (get_covariance(cov_date, :GBP_IR_aB, :GBP_IR_a)/sqrt(get_variance(cov_date, :GBP_IR_aB) * get_variance(cov_date, :GBP_IR_a)))) < tol

Random.seed!(1234)
## Test random draws
function abs_value_of_dict_differences(dictA::Dict, dictB::Dict)
   differences = merge(-, dictA, dictB)
   return sum(abs.(values(differences)))
end
abs_value_of_dict_differences(get_draws(cov_date), get_draws(cov_date)) > tol
abs_value_of_dict_differences(get_draws(cov_date;  uniform_draw = rand(5)), get_draws(cov_date; uniform_draw =  rand(5))) > tol

function test_random_points_pdf(covar::CovarianceAtDate)
    draws = get_draws(covar)
    pdf_val = pdf(covar, draws)
    return (pdf_val >= 0.0)
end
all([test_random_points_pdf(cov_date) for i in 1:1000])

# Distribution Testing
function SplitDicts(dictarray)
    return get.(dictarray, :USD_IR_a, 0), get.(dictarray, :USD_IR_aB, 0), get.(dictarray, :GBP_IR_a, 0), get.(dictarray, :GBP_IR_aB, 0), get.(dictarray, :GBP_FX, 0)
end
defaults = get_draws(cov_date,100000)
normal_twister = Mersenne(MersenneTwister(123), length(cov_date.covariance_labels_))
normals = get_draws(cov_date,100000; number_generator = normal_twister)
normal_samples = SplitDicts(normals)

s = SobolGen(SobolSeq(length(cov_date.ito_set_.ito_integrals_)))
sobols = get_draws(cov_date, 100000; number_generator = s)
sobol_samples = SplitDicts(sobols)
zero_draws = get_zero_draws(cov_date,2)

abs(var(normal_samples[1]) - get_variance(cov_date, :USD_IR_a))  < 0.001
abs(var(normal_samples[2]) - get_variance(cov_date, :USD_IR_aB)) < 0.11
abs(var(normal_samples[3]) - get_variance(cov_date, :GBP_IR_a))  < 1e-06
abs(var(normal_samples[4]) - get_variance(cov_date, :GBP_IR_aB)) < 1e-04
abs(var(normal_samples[5]) - get_variance(cov_date, :GBP_FX))    < 0.02

abs(var(sobol_samples[1])  - get_variance(cov_date, :USD_IR_a))  < 1e-04
abs(var(sobol_samples[2])  - get_variance(cov_date, :USD_IR_aB)) < 0.001
abs(var(sobol_samples[3])  - get_variance(cov_date, :GBP_IR_a))  < 1e-05
abs(var(sobol_samples[4])  - get_variance(cov_date, :GBP_IR_aB)) < 1e-05
abs(var(sobol_samples[5])  - get_variance(cov_date, :GBP_FX))    < 0.02

abs(cov(sobol_samples[5], sobol_samples[1]) - get_covariance(cov_date, :GBP_FX, :USD_IR_a)) < 0.001

#  Test likelihood
function test_random_points_loglikelihood(covar::CovarianceAtDate)
    draws = get_draws(covar)
    log_likelihood_val = log_likelihood(covar, draws)
    return log_likelihood_val > -Inf
end
all([test_random_points_loglikelihood(cov_date) for i in 1:1000])

# Testing data conversions - From draws
draws = get_draws(cov_date, 7)
arr, labs = to_array(draws; labels = [:USD_IR_a, :USD_IR_aB, :GBP_IR_a, :GBP_IR_aB])
size(arr) == (7,4)
labs ==  [:USD_IR_a, :USD_IR_aB, :GBP_IR_a, :GBP_IR_aB]
arr2, labs2 = to_array(draws)
size(arr2) == (7,5)
Set(labs2) == Set([:USD_IR_a, :USD_IR_aB, :GBP_IR_a, :GBP_IR_aB, :GBP_FX])
dd = to_dataframe(draws; labels = [:USD_IR_a, :USD_IR_aB, :GBP_IR_a, :GBP_IR_aB])
size(dd) == (7,4)
dd2 = to_dataframe(draws)
size(dd2) == (7,5)
Set(names(dd2)) == Set([:USD_IR_a, :USD_IR_aB, :GBP_IR_a, :GBP_IR_aB, :GBP_FX])

# Testing data conversions - From array
draws = to_draws(arr; labels = labs)
length(draws) == 7
Set(collect(keys(draws[1]))) == Set([:USD_IR_a, :USD_IR_aB, :GBP_IR_a, :GBP_IR_aB])
draws = to_draws(arr)
length(draws) == 7
Set(collect(keys(draws[1]))) == Set([:x1, :x2, :x3, :x4])
dd = to_dataframe(arr; labels = labs)
size(dd)[1] == 7
Set(collect(names(dd))) == Set([:USD_IR_a, :USD_IR_aB, :GBP_IR_a, :GBP_IR_aB])
dd2 = to_dataframe(arr)
size(dd2)[1] == 7
Set(names(dd2)) == Set([:x1, :x2, :x3, :x4])

# Testing data conversions - From dataframe
draws = to_draws(dd; labels = labs)
length(draws) == 7
Set(collect(keys(draws[1]))) == Set([:USD_IR_a, :USD_IR_aB, :GBP_IR_a, :GBP_IR_aB])
draws = to_draws(dd)
length(draws) == 7
Set(collect(keys(draws[1]))) == Set([:USD_IR_a, :USD_IR_aB, :GBP_IR_a, :GBP_IR_aB])
X, labs = to_array(dd; labels = labs)
size(X) == (7, 4)
Set(labs) == Set([:USD_IR_a, :USD_IR_aB, :GBP_IR_a, :GBP_IR_aB])
X, labs = to_array(dd)
size(X) == (7, 4)
Set(labs) == Set([:USD_IR_a, :USD_IR_aB, :GBP_IR_a, :GBP_IR_aB])

# Testing hypercube generation
confidence_level = 0.95
num = 1000000
confidence_hc = get_confidence_hypercube(covar, confidence_level, num)
# Now each edge should be same number of standard deviations away from the mean:
devs = confidence_hc[:USD_IR_a][2]/sqrt(covar.covariance_[1,1])
abs(devs - confidence_hc[:GBP_IR_aB][2]/sqrt(covar.covariance_[2,2])) < tol
abs(devs - confidence_hc[:GBP_IR_a][2]/sqrt(covar.covariance_[3,3]))  < tol
abs(devs - confidence_hc[:USD_IR_aB][2]/sqrt(covar.covariance_[4,4])) < tol
abs(devs - confidence_hc[:GBP_FX][2]/sqrt(covar.covariance_[5,5]))    < tol
# And the confidence hypercube should contain confidence_level of the distribution.
function _sobols(chol, num::Int, sob_seq::SobolGen)
    dims = size(chol)[1]
    array = Array{Float64,2}(undef, num, dims)
    for i in 1:num
        sobs = next!(sob_seq)
        normal_draw = quantile.(Ref(Normal()), sobs)
        scaled_draw = chol * normal_draw
        array[i,:] = scaled_draw
    end
    return array
end
function estimate_mass_in_hypercube()
    dims = length(covar.covariance_labels_)
    data = _sobols(covar.chol_, num, SobolGen(SobolSeq(dims)))
    number_of_draws = size(data)[1]
    in_confidence_area = 0
    cutoffs = devs .* sqrt.(diag(covar.covariance_))
    for i in 1:number_of_draws
        in_confidence_area += all(abs.(data[i,:]) .< cutoffs)
    end
    mass_in_area = in_confidence_area/number_of_draws
    return mass_in_area
end
abs(estimate_mass_in_hypercube() - confidence_level) < 0.01

# Testing hypercube generation
confidence_level = 0.5
num = 1000
confidence_hc = get_confidence_hypercube(covar, confidence_level, 500)
# Now each edge should be same number of standard deviations away from the mean:
devs = confidence_hc[:USD_IR_a][2]/sqrt(covar.covariance_[1,1])
abs(devs - confidence_hc[:GBP_IR_aB][2]/sqrt(covar.covariance_[2,2])) < tol
abs(devs - confidence_hc[:GBP_IR_a][2]/sqrt(covar.covariance_[3,3]))  < tol
abs(devs - confidence_hc[:USD_IR_aB][2]/sqrt(covar.covariance_[4,4])) < tol
abs(devs - confidence_hc[:GBP_FX][2]/sqrt(covar.covariance_[5,5]))    < tol
# And the confidence hypercube should contain confidence_level of the distribution.
abs(estimate_mass_in_hypercube() - confidence_level) < 0.03

# Testing Antithetic Variates
normals_antithetic = get_draws(cov_date,100000; number_generator = normal_twister, antithetic_variates = true)
normal_samples_antithetic = SplitDicts(normals_antithetic)
## Antithetic should have lower variance due to some data copying.
var(normal_samples[1]) > var(normal_samples_antithetic[1])
## Antithetic each pair should sum to 0. So whole vector should sum to zero.
sum(normal_samples_antithetic[1]) < tol
sum(normal_samples_antithetic[2]) < 100*tol
sum(normal_samples_antithetic[3]) < tol
sum(normal_samples_antithetic[4]) < tol
sum(normal_samples_antithetic[5]) < tol

# Now we will do a test for finding the expectation of a lognormal.
theoretical_expectation = exp(get_variance(cov_date, :USD_IR_a)/2)
log_normal_samples = exp.(normal_samples[1])
log_normal_expectation_estimate = mean(log_normal_samples)
log_normal_antithetic_samples = exp.(normal_samples_antithetic[1])
log_normal_antithetic_expectation_estimate = mean(log_normal_antithetic_samples)
abs(log_normal_antithetic_expectation_estimate - theoretical_expectation) < abs(log_normal_expectation_estimate - theoretical_expectation)


text:
using GridapFSI

output = main(
  Um=0.2,
  Re=20,
  rho_s=1.0e-3,
  strategy="biharmonic",
  alpha_m=1.0e-6,
  alpha_m_weight="volume",
  model="../models/elasticFlagFine.json",
  dt=1.0,
  tf=25.0,
  theta=0.6
)


text:
const DOC_OBSERVATIONS =
    "on predictions `ŷ`, "*
    "given ground truth observations `y`. "
const DOC_WEIGHTS =
    "Optionally specify per-sample weights, `w`. "
const DOC_CLASS_WEIGHTS =
    "An optional `AbstractDict`, denoted `class_w` above, "*
    "keyed on `levels(y)`, specifies class weights. "

macro create_aliases(M_ex)
    esc(quote
        M = $M_ex
        for alias in Symbol.(instances(M))
        # isdefined(parentmodule(M), alias) || eval(:(const $alias = $M()))
        eval(:(const $alias = $M()))
        end
        end)
end

function detailed_doc_string(M; typename="", body="", footer="", scitype="")

    _instances = _decorate(instances(M))
    human_name = MLJBase.human_name(M)
    if isempty(scitype)
        scitype = target_scitype(M) |> string
    end

    if isempty(typename)
        ret = "    $M\n\n"
    else
        ret = "    MLJBase.$typename\n\n"
    end

    ret *= "A measure type for $(human_name)"
    isempty(_instances) ||
        (ret  *= ", which includes the instance(s): "*
         "$_instances")
    ret *= ".\n\n"
    ret *= "    $(name(M))()(ŷ, y)\n"
    supports_weights(M) &&
        (ret *= "    $(name(M))()(ŷ, y, w)\n")
    supports_class_weights(M) &&
        (ret *= "    $(name(M))()(ŷ, y, class_w)\n")
    ret *= "\n"
    if isempty(fieldnames(M))
            ret *= "Evaluate the $(human_name) "
    else
        ret *= "Evaluate the default instance of $(name(M)) "
    end
    ret *= "$DOC_OBSERVATIONS"
    supports_weights(M) &&
        (ret *= DOC_WEIGHTS)
    supports_class_weights(M) &&
        (ret *= DOC_CLASS_WEIGHTS)
    ret *= "\n\n"
    isempty(body) || (ret *= "$body\n\n")
    ret *= "Requires `scitype(y)` to be a subtype of $scitype; "
    ret *= "`ŷ` must be an array of $(prediction_type(M)) predictions. "
    isempty(footer) ||(ret *= "\n\n$footer")
    ret *= "\n\n"
    ret *= "For more information, run `info($(name(M)))`. "
    return ret
end

_err_create_docs() = error(
    "@create_docs syntax error. Usage: \n"*
    "@create_docs(MeasureType, typename=..., body=..., scitype=..., footer=...")
macro create_docs(M_ex, exs...)
    M_ex isa Symbol || _err_create_docs()
    t = ""
    b = ""
    s = ""
    f = ""
    for ex in exs
        ex.head == :(=) || _err_create_docs()
        ex.args[1] == :typename && (t = ex.args[2])
        ex.args[1] == :body &&     (b = ex.args[2])
        ex.args[1] == :scitype &&  (s = ex.args[2])
        ex.args[1] == :footer &&   (f = ex.args[2])
    end
    esc(quote
        "$(detailed_doc_string($M_ex, typename=$t, body=$b, scitype=$s, footer=$f))"
        function $M_ex end
        end)
end

# TODO: I wonder why this is not a macro?

"""
    metadata_measure(T; kw...)

Helper function to write the metadata for a single measure.
"""
function metadata_measure(T; name::String="",
                          human_name="",
                          instances::Vector{String}=String[],
                          target_scitype=Unknown,
                          prediction_type::Symbol=:unknown,
                          orientation::Symbol=:unknown,
                          aggregation=Mean(),
                          is_feature_dependent::Bool=false,
                          supports_weights::Bool=true,
                          supports_class_weights::Bool=false,
                          docstring::String="",
                          distribution_type=Unknown)
    pred_str        = "$prediction_type"
    orientation_str = "$orientation"
#    dist = ifelse(ismissing(distribution_type), missing, "$distribution_type")
    ex = quote

        # traits common with models:
        if !isempty($name)
            StatisticalTraits.name(::Type{<:$T}) = $name
        end
        if !isempty($docstring)
            StatisticalTraits.docstring(::Type{<:$T}) = $docstring
        end
        StatisticalTraits.target_scitype(::Type{<:$T}) = $target_scitype
        StatisticalTraits.prediction_type(::Type{<:$T}) = Symbol($pred_str)
        StatisticalTraits.supports_weights(::Type{<:$T}) = $supports_weights

        # traits specific to measures:
        if !isempty($instances)
            StatisticalTraits.instances(::Type{<:$T}) = $instances
        end
        if !isempty($human_name)
            StatisticalTraits.human_name(::Type{<:$T}) = $human_name
        end
        StatisticalTraits.orientation(::Type{<:$T}) = Symbol($orientation_str)
        StatisticalTraits.aggregation(::Type{<:$T}) = $aggregation
        StatisticalTraits.is_feature_dependent(::Type{<:$T}) =
            $is_feature_dependent
        StatisticalTraits.supports_class_weights(::Type{<:$T}) =
            $supports_class_weights
        StatisticalTraits.distribution_type(::Type{<:$T}) = $distribution_type

    end
    parentmodule(T).eval(ex)
end