text:
using SeisIO, SeisNoise, JLD2, Statistics, FFTW, DSP, Plots

include("utils.jl")
include("partition.jl")

"""
    selective_stacking(data::CorrData, reference::CorrData; threshold::Float64=0.0, slice::Union{Bool, Float64, Array{Float64,1}}=false, metric::String="cc", win_len::Float64=10.0, win_step::Float64=5.0)

Stack the windows in a CorrData object that exceed a correlation-coefficient threshold with respect to a reference.

# Arguments
- `data::CorrData,`    : Input CorrData matrix used to define the reference and each window
- `reference::CorrData`    : Reference cross-correlation function to use in computing the correlation coefficient
- `threshold::Float64,`    : Value of correlation-coefficient below which we zero windows
- `slice::Union{Bool, Float64, Array{Float64,1}}`    : whether to slice the cross-correlations before computing convergence. If so, a single Float64 specifies the
                                                       lag at which to slice to keep just the ballistic wave and an array of Float64 specifies the start and end lag
                                                       to keep only the coda.

# Output
- `stackedData::CorrData,`    : Slectively stacked data
- `cList::Array{Float64,1}`    : Correlation coefficient of each window with respect to the reference
"""
function selective_stacking(data::CorrData, reference::CorrData; threshold::Float64=0.0, slice::Union{Bool, Float64, Array{Float64,1}}=false, metric::String="cc", coh_win_len::Float64=10.0, coh_win_step::Float64=5.0, filter::Union{Bool, Array{Float64,1}}=false, phase_smoothing::Float64=0.0)
    # slice data if given a time window
    # TODO: find a better way to pass arguments to this function that only apply to coh, or only to cc. e.g., win_len has no meaning if metric=="cc"
    if typeof(slice) != Bool
        ref = copy(reference)
        d = copy(data)
        if typeof(slice)==Float64
            # time vector
            tvec = -reference.maxlag:1/reference.fs:reference.maxlag
            # find all times within [-slice,slice] time
            t_inds = findall(x->(x>=-slice && x<=slice), tvec)

            # slice reference (maintain 2d array)
            ref.corr = reference.corr[t_inds, :]
            # slice data (maintain 2d array)
            d.corr = data.corr[t_inds, :]
        elseif typeof(slice)==Array{Float64,1}
            # convert startlag/endlag[s] to startlag/windowlength[samples]
            win_len = Int(diff(slice)[1] * reference.fs)
            startlag = Int(slice[1] * reference.fs)

            # partition reference to keep only the coda
            ref.corr = partition(reference.corr, startlag, win_len)
            # partition data to keep only the coda
            d.corr = partition(data.corr, startlag, win_len)
        else
            println("Please choose an allowable slicing operation. Exiting.")
            exit()
        end
    else
        # default to full reference cross-correlation
        ref = reference
        d = data
    end

    if metric=="cc"
        # compute correlation coefficient for each window in data with respect to reference
        cList = get_cc(d, ref)
    elseif metric=="coh"
        # compute coherence for each window as a function of frequency
        cohList = get_coh(d, ref, coh_win_len, coh_win_step)

        # use only coherencies for requested frequencies in the mean
        # default to using all frequencies
        if filter==false
            filter=[0.0, ref.fs/2]
        end
        # set acceptable frequencies to use in mean coh calculation
        freqs = rfftfreq(Int(coh_win_len*ref.fs), ref.fs)
        freq_inds = findall(x->(x>=filter[1] && x<=filter[2]), freqs)

        # compute correlation coefficient for each window in data with respect to reference
        ccList=get_cc(d, ref)

        # compute mean coherence
        mcohList = mean(cohList[freq_inds, :], dims=1)[1,:]

        # scale coherence by sign of the correlation coefficient to make extraction of constructively coherent traces possible
        cList = mcohList .* sign.(ccList)
    end

    # find windows that exceed the threshold
    if typeof(slice)==Array{Float64,1}
        # keep only xcorrs coherent in positive && negative coda
        # for now, we take only the windowed cross-correaltions that improve both the positive and negative coda
        good_fit_neg = findall(x->(x>=threshold), cList[1:Int(length(cList)/2)])
        good_fit_pos = findall(x->(x>=threshold), cList[Int(length(cList)/2+1):end])
        good_fit = intersect(good_fit_neg, good_fit_pos)
    else
        # find cross-correlations that fall below the cc threshold
        good_fit = findall(x->(x>=threshold), cList)
    end

    # copy data to extract well-fitting windows
    tempData = copy(data)
    tempData.corr = tempData.corr[:, good_fit]

    # linearly stack all data that exceeds the correlation-coefficient threshold
    stackedData = stack(tempData, allstack=true, phase_smoothing=phase_smoothing)

    return stackedData, cList
end

"""
    get_cc(data::CorrData, reference::CorrData)

Compute the correlation coefficient for each window in a CorrData object with respect to a reference

# Arguments
- `data::CorrData`    : Input CorrData matrix used to define the reference and each window
- `reference::CorrData`    : Reference cross-correlation function to use in computing the correlation coefficient

# Output
- `cList::Array{Float64,1}`    :  Correlation coefficient for each window with respect to the reference
"""
function get_cc(data::CorrData, reference::CorrData)
    # array of correlation coefficients for each constituent cross-correlation window
    cList = Array{Float32,1}(undef,size(data.corr)[2])

    # iterate over each windowed cross-correlation and compute correlation-coefficient
    for j=1:size(data.corr)[2]
        cList[j] = cor(reference.corr[:,1], data.corr[:,j])
    end

    return cList
end

"""
    get_coh(data::CorrData, reference::CorrData, win_len::Float64, win_step::Float64)

Compute the coherence between each window of a cross-correlation CorrData object with respect to a reference.

# Arguments
- `data::CorrData`    : Cross-correlation data to compute coherency
- `reference::CorrData`    : Reference cross-correlation function
- `win_len::Float64`    : Window length [s] when computing the psd using Welch's method. This parameter decides the resolvable frequencies in the coherence calculation
- `win_step::Float64`    : Step size [s] for windowing of each cross-correlation when computing the psd using Welch's method

# Outputs
- `coh::Array{Float32,2}`    : Matrix of size (nFreqs, nWindows) containing coherency values for each window and frequency
"""
function get_coh(data::CorrData, reference::CorrData, win_len::Float64, win_step::Float64)
    # convert window length and window step from time to samples
    win_len = Int(win_len * data.fs)
    win_step = Int(win_step * data.fs)

    # normalize signals
    data.corr = normalize(data.corr)
    reference.corr = normalize(reference.corr)

    # compute the coherence for each window in CorrData dim=2
    Ryy = welch_psd(reference.corr[:, 1], reference.corr[:, 1], data.fs, win_len, win_step)
    # initialize arrays to hold psd for each window
    Rxy = Array{eltype(data.corr), 2}(undef, size(Ryy,1), size(data.corr,2))
    Rxx = Array{eltype(data.corr), 2}(undef, size(Ryy,1), size(data.corr,2))
    # iterate over windows of CorrData
    for i=1:size(data.corr,2)
        # estimate the power spectral density using Welch's method
        Rxy[:, i] = welch_psd(data.corr[:, i], reference.corr[:, 1], data.fs, win_len, win_step, cross=true)
        Rxx[:, i] = welch_psd(data.corr[:, i], data.corr[:, i], data.fs, win_len, win_step)
    end

    # compute the sample coherence function
    coh = (Rxy.^2) ./ (Rxx .* Ryy)

    return coh
end

"""
    welch_psd(A::Array{R,1}, B::Array{R,1}, fs::Float64, win_len::Int64, win_step::Int64; cross=false) where R<:Real

Compute the cross power spectral density of two given arrays using Welch's method

Reference: Welch, P. D. (1967), "The use of Fast Fourier Transform for the estimation of power spectra:
A method based on time averaging over short, modified periodograms" (PDF),
IEEE Transactions on Audio and Electroacoustics, AU-15 (2): 70–73, doi:10.1109/TAU.1967.1161901

# Arguments
- `A::Array{R,1} where R<:Real`    : Signal 1
- `B::Array{R,1} where R<:Real`    : Signal 1
- `fs::Float64`    : sampling frequency
- `win_len::Int64`    : window length [samples] for FFT computation
- `win_step::Int64`    : step size [samples] for window advancement
- `cross::Bool`    : whether or not this is the psd of a cross-correlation

# Outputs
- `real.(cpsd)::Array{R,2}`    : Estimated cross power spectral density of two given time series
"""
function welch_psd(A::Array{R,1}, B::Array{R,1}, fs::Float64, win_len::Int64, win_step::Int64; cross=false) where R<:Real
    # partition the data into m windows
    A = window(A, win_len, win_step)
    B = window(B, win_len, win_step)

    # compute the FFT on the windowed data without bandpassing
    A_fft = process_fft(A, 0.0001, fs/2 - 0.0001, fs, corners=1)
    B_fft = process_fft(B, 0.0001, fs/2 - 0.0001, fs, corners=1)

    # compute the squared-magnitude of each signal block
    AB_fft2 = A_fft .* conj.(B_fft)

    # average over signal blocks to obtain an estimate of the cross power spectral density
    cpsd = mean(AB_fft2, dims=2)

    if cross
        power!(cpsd)
    end

    return real.(cpsd)
end

"""
    power!(A::AbstractArray)

Compute the power spectrum of a given FFT

# Arguments
- `A::AbstractArray`    : Input FFT

# Outputs
-  `power::AbstractArray`    : Power spectrum
"""
function power!(A::AbstractArray)
    return (abs.(A)).^2
end
power(A::AbstractArray) = (U = deepcopy(A);power!(U);return U)

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation


text:
function Documented(def::StructDefinition{false})
    spec = def.spec
    ext = extension(spec)
    doc = string(
        """
        Minimalistic wrapper for $(spec.name).
        """,
        !isnothing(ext) ? """

        Extension: $(ext.name)
        """ : "",
        """

        API documentation: https://www.khronos.org/registry/vulkan/specs/1.2-extensions/man/html/$(spec.name).html
        """,
    )
    Documented(def, doc)
end

function Documented(def::StructDefinition{true})
    spec = def.spec
    ext = extension(spec)
    doc = string(
        """
        High-level wrapper for $(spec.name).
        """,
        !isnothing(ext) ? """

        Extension: $(replace(string(ext.name), '_' => "\\\\_"))
        """ : "",
        """

        [API documentation](https://www.khronos.org/registry/vulkan/specs/1.2-extensions/man/html/$(spec.name).html)
        """,
    )
    Documented(def, doc)
end

function Documented(def::Constructor{HandleDefinition})
    doc = string(' '^4, reconstruct_call(def.p), '\n'^2)
    Documented(def, doc)
end

backquoted(arg) = string('`', arg, '`')

function document_return_codes(spec::SpecFunc)
    res = ""
    if !isempty(spec.success_codes) || !isempty(spec.error_codes)
        res *= "Return codes:"
    end
    if !isempty(spec.success_codes) && must_return_success_code(spec)
        res *= string("\n- Success:\n  - ", join(backquoted.(remove_vk_prefix.(spec.success_codes)), "\n  - "))
    end
    if !isempty(spec.error_codes)
        res *= string("\n- Error:\n  - ", join(backquoted.(remove_vk_prefix.(spec.error_codes)), "\n  - "))
    end

    if !isempty(res)
        res = '\n' * res * '\n'
    end

    res
end

document_return_codes(spec::SpecHandle) = document_return_codes(add_constructor(spec))
document_return_codes(spec::SpecStruct) = ""

function append_to_argdoc!(argdocs, spec, str)
    i = findfirst(==(wrap_identifier(spec)), first.(argdocs))
    argdocs[i] = concat_right(argdocs[i], str)
end

function Documented(def::APIFunction{SpecFunc})
    p = def.p
    spec = def.spec
    params = children(spec)
    externsync_params = filter(x -> x.is_externsync, params)
    argdocs = document_arguments(p)
    if !isempty(externsync_params)
        for param ∈ externsync_params
            append_to_argdoc!(argdocs, param, " (externsync)")
        end
    end
    extra = ""
    if spec isa SpecFunc
        if any(is_data_with_retrievable_size, params)
            extra = """

            !!! warning
                The pointer returned by this function holds memory owned by Julia. It is therefore **your** responsibility to free it after use (e.g. with `Libc.free`).
            """
        elseif any(is_opaque_data, params)
            opaque_params = params[findall(is_opaque_data, params)]
            map(opaque_params) do param
                append_to_argdoc!(
                    argdocs,
                    param,
                    " (must be a valid pointer with `$(wrap_identifier(len(param)))` bytes)",
                )
            end
        end
    end

    extra *= """

    [API documentation](https://www.khronos.org/registry/vulkan/specs/1.2-extensions/man/html/$(spec.name).html)
    """

    ext = extension(spec)

    doc = string(
        """
            $(reconstruct_call(p))
        """,
        !isnothing(ext) ? """

        Extension: $(replace(string(ext.name), '_' => "\\\\_"))
        """ : "",
        document_return_codes(spec),
        isempty(argdocs) ? "" : """

        Arguments:
        - $(join(last.(argdocs), "\n- "))
        """,
        extra,
    )
    Documented(def, doc)
end

function document_arguments(p)
    args, kwargs = p[:args], p[:kwargs]
    if isempty(args) && isempty(kwargs)
        []
    else
        backquoted_args = map(x -> name(x) => backquoted(x), args)
        documented_kwargs = map(kwargs) do kwarg
            identifier, value = kwarg.args
            identifier => string(backquoted(identifier), ": defaults to ", backquoted(value))
        end
        [backquoted_args; documented_kwargs]
    end
end

function docstring(ex, docstring)
    Dict(
        :category => :doc,
        :ex => ex,
        :docstring => docstring,
    )
end

concat_right(pair, val) = pair.first => (pair.second * val)

ext: jl


text:
"This is a single line comment."
"""
This is a multi line comment.
"""

max_stars_repo_name: smenjas/programming-languages-compared

max_issues_repo_name: smenjas/programming-languages-compared

max_forks_repo_name: smenjas/programming-languages-compared


text:
l="\n"
△(n,k)=((i,j)->(" "^(k-j+1)*"*"^i*" "^(k-j))^(n==k ? 2 : 1)).(1:2:2n,1:n)
~n=print(join(△(n,2n),l),l,join(△(n,n),l))

max_stars_repo_name: pgrepds/Codegolf

max_issues_repo_name: pgrepds/Codegolf

max_forks_repo_name: pgrepds/Codegolf


text:
module AReg
using PolynomialRoots
using StatsBase
using Statistics

export ARp, ARp_est, fit , simulate,simulateAR1, unc_mean
##
"Autoregressive Structure (only normal errors are considered for the moment)"
mutable struct ARp
    vecpar::Vector{Real}
    sigma::Real
    start_vals::Vector{Real}
end

mutable struct ARp_est
    ar::ARp
    obs::Vector{Real}
    estimates::Vector{Real}
    sigma::Real
end

function unc_mean(arp::ARp)
    "compute the unconditional mean of an ARp type"
    p = length(arp.vecpar)
    p==0 ? (return 0) : (return arp.vecpar[1]/(1-sum(arp.vecpar[2:end])))
end
function is_ar_stationary(arp::ARp)
    coeff = arp.vecpar
    #coefficients of the characteristic polynomial
    char_pol_coeff = [coeff[1]; - coeff[2:end]]
    #compute the roots of the char poly
    tmp =abs2.(roots(Array{Float64}(char_pol_coeff)))
    #if all the roots lie outside the unit circle then the process is stationary
    all(tmp.>1) ? (return true,unc_mean(arp)) : (return false)
end

##
"fit an AR(p) with predetermined p, the starting values are not estimated
To Do : need to add a method for the observations only.. it should return the
ARp_est with optimal p
"
function fitARp(obs::Vector{<:Real}, order;demean=true)
    p = order # number of autoregressive terms (constant excluded)

    ## Ols estimation of the AR(p) coefficients
    Y = obs[p+1:end]
    T = length(Y)
    X = ones(T,p+1)
    #println(p)
    for i=1:p
        X[:,p+2-i] = obs[p+1-i:end-i]
    end

    est_par = (X'*X)\(X'Y)
    #standard deviation of residuals
    sigma = std(X*est_par - Y)
    return est_par, sigma
end

function fitARp(obs::Vector{<:Real},arp::ARp;demean=true)
    p = length(arp.vecpar)-1
    est_par, sigma = fitARp(obs, p)

    return ARp_est(arp,obs,est_par, sigma)
end

"Sample an AR process of lenght T and  parameters specified in vecpar starting
from phi_0 to phi_p where p is the order of the process"
function simulate(ar::ARp,T::Int=1000)
    sigma = ar.sigma
    vecpar = ar.vecpar
    start_vals = ar.start_vals
    #Y_T = zeros(T)
    p = length(vecpar)-1
    epsilon_T = sigma.* randn(T)
    Y_T = epsilon_T

    if isempty(start_vals) # if no start values are given pick them at ranUtilities
        if p>0
            start_vals = sigma.*randn(p)
        else
            start_vals = []
        end
    end
    lagged_vals = start_vals
    for t = 1:T
        Y_T[t] +=  vecpar'*[1;lagged_vals ]
        if p>0
            next_lagged_vals = circshift(lagged_vals,-1)
            next_lagged_vals[end] = Y_T[t]
        else
            next_lagged_vals = Y_T[t]
        end
        lagged_vals = next_lagged_vals
    end
    return Y_T
end

function simulateAR1(vecPar::Vector{<:Real},sigma::Float64;T::Int=1000)
    #start on the unconditional mean
    μ = vecPar[1]./(1-vecPar[2])
    #Y_T = zeros(T)
    epsilon_T = sigma.* randn(T)
    Y_T = epsilon_T
    lagged_val = μ
    for t = 1:T
        Y_T[t] +=  vecPar[1] + vecPar[2] * lagged_val
        lagged_val = Y_T[t]
    end
    return Y_T
end

function ar1_variance(x)
    try
        return AReg.fitARp(Float64.(x), 1)[2]
    catch
        return NaN
    end
end
export ar1_variance

end


text:
using Test
using MiniLogging
using MiniLogging: @debug, @info, @warn, @error, @critical
using MiniLogging.Hierarchy
using MiniLogging.Hierarchy: Ancestors

@test collect(Ancestors("")) == []
@test collect(Ancestors("a")) == [""]
@test collect(Ancestors("a.bb.ccc")) == ["a.bb", "a", ""]
@test collect(Ancestors("❤.❤❤.❤❤❤")) == ["❤.❤❤", "❤" , ""]
@test_throws ErrorException collect(Ancestors(".a.bb.ccc"))
@test_throws ErrorException collect(Ancestors("a.b.cc."))
@test_throws ErrorException collect(Ancestors("a.b..cc"))

t = Tree()
push!(t, "a.b.c")
@test parent_node(t, "a.b.c") == ""
push!(t, "a.b.c")
@test parent_node(t, "a.b.c") == ""
push!(t, "a.b")
@test parent_node(t, "a.b") == ""
@test parent_node(t, "a.b.c") == "a.b"
push!(t, "a.b.c.d")
@test parent_node(t, "a.b.c.d") == "a.b.c"
push!(t, "a.b.c.d2")
@test parent_node(t, "a.b.c.d2") == "a.b.c"
push!(t, "a.b.c.d.❤")
@test parent_node(t, "a.b.c.d.❤") == "a.b.c.d"
push!(t, "a.b.❤.d")
@test parent_node(t, "a.b.❤.d") == "a.b"
push!(t, "a")
@test parent_node(t, "a") == ""
push!(t, "")
@test parent_node(t, "") == ""

open("stdout.out", "w") do f1
    open("stderr.out", "w") do f2
        redirect_stdout(f1) do
            redirect_stderr(f2) do
                include("test_log.jl")
            end
        end
    end
end

out1 = readlines("stdout.log")
err1 = readlines("stderr.log")
out2 = readlines("stdout.out")
err2 = readlines("stderr.out")

@test out1 == out2
@test err1 == err2

max_stars_repo_name: colinfang/MiniLogging

max_issues_repo_name: colinfang/MiniLogging

max_forks_repo_name: colinfang/MiniLogging


text:
export BlockSparse,
       BlockSparseTensor,
       Block,
       block,
       BlockOffset,
       BlockOffsets,
       blockoffsets,
       blockview,
       nnzblocks,
       nnz,
       findblock,
       isblocknz

#
# BlockSparse storage
#

struct BlockSparse{ElT,VecT,N} <: TensorStorage{ElT}
  data::VecT
  blockoffsets::BlockOffsets{N}  # Block number-offset pairs
  function BlockSparse(data::VecT,
                       blockoffsets::BlockOffsets{N};
                       sorted=true) where {VecT<:AbstractVector{ElT},
                                           N} where {ElT}
    sorted && check_blocks_sorted(blockoffsets)
    new{ElT,VecT,N}(data, blockoffsets)
  end
end

function BlockSparse(::Type{ElT},
                     blockoffsets::BlockOffsets,
                     dim::Integer;
                     vargs...) where {ElT<:Number}
  return BlockSparse(zeros(ElT, dim), blockoffsets; vargs...)
end

function BlockSparse(::Type{ElT},
                     ::UndefInitializer,
                     blockoffsets::BlockOffsets,
                     dim::Integer; vargs...) where {ElT<:Number}
  return BlockSparse(Vector{ElT}(undef,dim),blockoffsets; vargs...)
end

BlockSparse(blockoffsets::BlockOffsets,
            dim::Integer; vargs...) = BlockSparse(Float64,blockoffsets,dim; vargs...)

BlockSparse(::UndefInitializer,
            blockoffsets::BlockOffsets,
            dim::Integer; vargs...) = BlockSparse(Float64,undef,blockoffsets,dim; vargs...)

#
# Random
#

function Base.randn(::Type{ <: BlockSparse{ElT}},
                    blockoffsets::BlockOffsets,
                    dim::Integer) where {ElT <: Number}
  return BlockSparse(randn(ElT, dim), blockoffsets)
end

#function BlockSparse{ElR}(data::VecT,offsets) where {ElR,VecT<:AbstractVector{ElT}} where {ElT}
#  ElT == ElR ? BlockSparse(data,offsets) : BlockSparse(ElR.(data),offsets)
#end
#BlockSparse{ElT}() where {ElT} = BlockSparse(ElT[],BlockOffsets())

function Base.similar(D::BlockSparse)
  return BlockSparse(similar(data(D)),blockoffsets(D))
end

# TODO: test this function
Base.similar(D::BlockSparse,
             ::Type{ElT}) where {ElT} = BlockSparse(similar(data(D),ElT),
                                                    copy(blockoffsets(D)))
Base.copy(D::BlockSparse) = BlockSparse(copy(data(D)),
                                        copy(blockoffsets(D)))

# TODO: check the offsets are the same?
function Base.copyto!(D1::BlockSparse,D2::BlockSparse)
  blockoffsets(D1) ≠ blockoffsets(D1) && error("Cannot copy between BlockSparse storages with different offsets")
  copyto!(data(D1),data(D2))
  return D1
end

# convert to complex
# TODO: this could be a generic TensorStorage function
Base.complex(D::BlockSparse{T}) where {T} = BlockSparse{complex(T)}(complex(data(D)),
                                                                    blockoffsets(D))

function Base.conj(D::BlockSparse{<: Real}; always_copy = false)
  if always_copy
    return copy(D)
  end
  return D
end

function Base.conj(D::BlockSparse; always_copy = false)
  if always_copy
    return conj!(copy(D))
  end
  return BlockSparse(conj(data(D)), blockoffsets(D))
end

function scale!(D::BlockSparse,α::Number)
  scale!(data(D),α)
  return D
end

Base.ndims(::BlockSparse{T,V,N}) where {T,V,N} = N

Base.eltype(::BlockSparse{T}) where {T} = eltype(T)
# This is necessary since for some reason inference doesn't work
# with the more general definition (eltype(Nothing) === Any)
Base.eltype(::BlockSparse{Nothing}) = Nothing
Base.eltype(::Type{BlockSparse{T}}) where {T} = eltype(T)

dense(::Type{<:BlockSparse{ElT,VecT}}) where {ElT,VecT} = Dense{ElT,VecT}

function Base.promote_rule(::Type{<:BlockSparse{ElT1,VecT1,N}},
                           ::Type{<:BlockSparse{ElT2,VecT2,N}}) where {ElT1,ElT2,VecT1,VecT2,N}
  return BlockSparse{promote_type(ElT1,ElT2),promote_type(VecT1,VecT2),N}
end

function Base.promote_rule(::Type{<:BlockSparse{ElT1,VecT1,N1}},
                           ::Type{<:BlockSparse{ElT2,VecT2,N2}}) where {ElT1,ElT2,VecT1,VecT2,N1,N2}
  return BlockSparse{promote_type(ElT1,ElT2),promote_type(VecT1,VecT2),NR} where {NR}
end

function Base.promote_rule(::Type{<:BlockSparse{ElT1,Vector{ElT1},N1}},
                           ::Type{ElT2}) where {ElT1,ElT2<:Number,N1}
  ElR = promote_type(ElT1,ElT2)
  VecR = Vector{ElR}
  return BlockSparse{ElR,VecR,N1}
end

function Base.convert(::Type{<:BlockSparse{ElR,VecR,N}},
                      D::BlockSparse{ElD,VecD,N}) where {ElR,VecR,N,ElD,VecD}
  return BlockSparse(convert(VecR,data(D)),
                     blockoffsets(D))
end

function Base.:*(D::BlockSparse,x::Number)
  return BlockSparse(x*data(D),blockoffsets(D))
end
Base.:*(x::Number,D::BlockSparse) = D*x

"""
blockdim(T::BlockSparse,pos::Int)

Get the block dimension of the block at position pos.
"""
blockdim(D::BlockSparse,
         pos::Int) = blockdim(blockoffsets(D),
                              nnz(D),
                              pos)

"""
blockdim(T::BlockSparse,block::Block)

Get the block dimension of the block.
"""
function blockdim(D::BlockSparse,
                  block::Block)
  pos = findblock(D, block)
  return blockdim(D, pos)
end

findblock(D::BlockSparse{<:Number,
                         <:AbstractVector,
                         N},
          block::Block{N};
          vargs...) where {N} = findblock(blockoffsets(D),
                                          block;
                                          vargs...)

"""
isblocknz(T::BlockSparse,
          block::Block)

Check if the specified block is non-zero.
"""
function isblocknz(T::BlockSparse{ElT,VecT,N},
                   block::Block{N}) where {ElT,VecT,N}
  isnothing(findblock(T,block)) && return false
  return true
end

# Given a specified block, return a Dense storage that is a view to the data
# in that block. Return nothing if the block is structurally zero
function blockview(T::BlockSparse,
                   block)
  #error("Block must be structurally non-zero to get a view")
  !isblocknz(T,block) && return nothing
  blockoffsetT = offset(T,block)
  blockdimT = blockdim(T,block)
  dataTslice = @view data(T)[blockoffsetT+1:blockoffsetT+blockdimT]
  return Dense(dataTslice)
end

function Base.:+(D1::BlockSparse,D2::BlockSparse)
  # This could be of order nnzblocks, avoid?
  if blockoffsets(D1) == blockoffsets(D2)
    return BlockSparse(data(D1)+data(D2),blockoffsets(D1))
  end
  blockoffsetsR,nnzR = union(blockoffsets(D1),nnz(D1),
                             blockoffsets(D2),nnz(D2))
  R = BlockSparse(undef,blockoffsetsR,nnzR)
  for (blockR,offsetR) in blockoffsets(R)
    blockview1 = blockview(D1,blockR)
    blockview2 = blockview(D2,blockR)
    blockviewR = blockview(R,blockR)
    if isnothing(blockview1)
      copyto!(blockviewR,blockview2)
    elseif isnothing(blockview2)
      copyto!(blockviewR,blockview1)
    else
      # TODO: is this fast?
      blockviewR .= blockview1 .+ blockview2
    end
  end
  return R
end

# Helper function for HDF5 write/read of BlockSparse
function offsets_to_array(boff::BlockOffsets{N}) where {N}
  nblocks = length(boff)
  asize = (N+1)*nblocks
  n = 1
  a = Vector{Int}(undef,asize)
  for bo in boff
    for j=1:N
      a[n] = bo[1][j]
      n += 1
    end
    a[n] = bo[2]
    n += 1
  end
  return a
end

# Helper function for HDF5 write/read of BlockSparse
function array_to_offsets(a,N::Int)
  asize = length(a)
  nblocks = div(asize,N+1)
  boff = BlockOffsets{N}(undef,nblocks)
  j = 0
  for b=1:nblocks
    boff[b] = Pair( ntuple(i->(a[j+i]),N) , a[j+N+1] )
    j += (N+1)
  end
  return boff
end

function HDF5.write(parent::Union{HDF5File, HDF5Group},
                    name::String,
                    B::BlockSparse)
  g = g_create(parent,name)
  attrs(g)["type"] = "BlockSparse{$(eltype(B))}"
  attrs(g)["version"] = 1
  if eltype(B) != Nothing
    write(g,"ndims",ndims(B))
    write(g,"data",data(B))
    off_array = offsets_to_array(blockoffsets(B))
    write(g,"offsets",off_array)
  end
end

function HDF5.read(parent::Union{HDF5File,HDF5Group},
                   name::AbstractString,
                   ::Type{Store}) where {Store <: BlockSparse}
  g = g_open(parent,name)
  ElT = eltype(Store)
  typestr = "BlockSparse{$ElT}"
  if read(attrs(g)["type"]) != typestr
    error("HDF5 group or file does not contain $typestr data")
  end
  N = read(g,"ndims")
  data = read(g,"data")
  off_array = read(g,"offsets")
  boff = array_to_offsets(off_array,N)
  return BlockSparse(data,boff)
end


text:
using DustExtinction
using DustExtinction: bounds, checkbounds
using Test, Measurements, Unitful, UnitfulAstro, Random

Random.seed!(9994445781)

include("deprecate.jl")
include("color_laws.jl")
include("dust_maps.jl")
include("fittable_laws.jl")
include("mixture_laws.jl")

@testset "interfaces" begin
    for LAW in [CCM89, OD94, CAL00, GCC09, VCG04, FM90, G16, F99, F04, F19, M14]
        @test bounds(LAW) == bounds(LAW())
        @test checkbounds(LAW, 1000) == checkbounds(LAW(), 1000)
        low, high = bounds(LAW)
        @test all(checkbounds.(LAW, rand(low:high, 1000)))
        @test all(isfinite.(LAW().(bounds(LAW))))
    end
    @test bounds(DustExtinction.ExtinctionLaw) == (0, Inf)
end

@testset "redden/deredden" begin
    wave = [3372., 4404., 5428., 6509., 8090.,
            3683., 4393., 5519., 6602., 8046.,
            12660., 16732., 22152.,
            5244., 6707., 7985., 9055.,
            6993.,
            3502., 4676., 4127.,
            4861., 5479.,
            3546., 4925., 6335., 7799., 9294.,
            3047., 4711., 5498.,
            6042., 7068., 8066.,
            4814., 6571., 8183.]

    # Av=0.3, Rv=3.1, law=ccm89
    ref_values = [0.6376288197244566, 0.6937993158943373, 0.7556003129054399, 0.7959216081971665, 0.8508142213648981, 0.6517368195297566, 0.6930246485833793, 0.7596467718783736, 0.7990835963390542, 0.8493553321733043, 0.9264835715958143, 0.9524307890106202, 0.9694546832913796, 0.7467375901587614, 0.8026674093639407, 0.847310247751628, 0.8771971841412943, 0.8125473475954554, 0.6442645437827309, 0.7127773090414787, 0.6749407318003154, 0.7249166081334651, 0.7578933014813821, 0.6461535564599767, 0.7288900774623704, 0.7900092679188825, 0.8409405763429691, 0.881972061390043, 0.6091986958411011, 0.7151404275193919, 0.7587309863781261, 0.7798780926274524, 0.8151697825084553, 0.8500202342098636, 0.7219200730427752, 0.7980286586611023, 0.8538476253948861]

    flux = ones(length(wave))
    output = @inferred broadcast((l, w, f)->redden(l, w, f; Av = 0.3), CCM89, wave, flux)
    @test output ≈ ref_values
    @test @inferred(broadcast((l, w, f)->deredden(l, w, f; Av = 0.3), CCM89, wave, output)) ≈ flux
    flux .= redden.(CCM89, wave, flux, Av = 0.3)
    @test flux ≈ ref_values

    # interfaces
    @test redden(CCM89, wave[1], flux[1]) == redden(CCM89(), wave[1], flux[1])
    @test redden(CCM89, wave[1], flux[1], Rv = 2.8) == redden(CCM89(Rv = 2.8), wave[1], flux[1])
    @test redden.(CCM89, wave, flux) == redden.(CCM89(), wave, flux)
    @test redden.(CCM89, wave, flux, Rv = 2.8) == redden.(CCM89(Rv = 2.8), wave, flux)

    # Measurements
    flux = ones(length(wave)) .± 0.1
    output = @inferred broadcast((l, w, f)->redden(l, w, f; Av = 0.3), CCM89, wave, flux)
    @test Measurements.value.(output) ≈ ref_values
    @test @inferred(broadcast((l, w, f)->deredden(l, w, f; Av = 0.3), CCM89, wave, output)) ≈ flux
    flux .= redden.(CCM89, wave, flux, Av = 0.3)
    @test flux ≈ ref_values

    # Unitful
    flux = ones(length(wave))u"Jy"
    wave = wave * u"angstrom"
    output = @inferred broadcast((l, w, f)->redden(l, w, f; Av = 0.3), CCM89, wave, flux)
    @test redden.(CCM89, wave, ustrip.(u"Jy", flux), Av = 0.3) ≈ ustrip.(u"Jy", output)
    @test ustrip.(output) ≈ ref_values
    @test @inferred(broadcast((l, w, f)->deredden(l, w, f; Av = 0.3), CCM89, wave, output)) ≈ flux
    @test deredden.(CCM89, wave, ustrip.(u"Jy", output), Av = 0.3) ≈ ustrip.(u"Jy", flux)
    flux .= redden.(CCM89, wave, flux, Av = 0.3)
    @test flux ≈ output
end


text:
module VulkanGen

using StructArrays: StructVector
using Accessors: @set
using Graphs
using MLStyle
using DocStringExtensions
using Reexport
using Dictionaries
using .Meta: isexpr

@template (FUNCTIONS, METHODS, MACROS) = """
                                         $(DOCSTRING)
                                         $(TYPEDSIGNATURES)
                                         """

@template TYPES = """
                  $(DOCSTRING)
                  $(TYPEDEF)
                  $(TYPEDSIGNATURES)
                  $(TYPEDFIELDS)
                  $(SIGNATURES)
                  """

include("types.jl")

include("spec/VulkanSpec.jl")
@reexport using .VulkanSpec

include("wrapper.jl")

end # module VulkanGen

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
using NCDatasets
using Dates
using TOML
using CFTime
using Random
using UnPack

tomlpath = joinpath(@__DIR__, "sbm_config.toml")
parsed_toml = TOML.parsefile(tomlpath)
config = Wflow.Config(tomlpath)

@testset "configuration file" begin
    @test parsed_toml isa Dict{String,Any}
    @test config isa Wflow.Config
    @test Dict(config) == parsed_toml
    @test pathof(config) == tomlpath
    @test dirname(config) == dirname(tomlpath)

    # test if the values are parsed as expected
    @test config.casename == "testcase"
    @test config.starttime === DateTime(2000)
    @test config.endtime === DateTime(2000, 2)
    @test config.output.path == "data/output_moselle.nc"
    @test config.output isa Wflow.Config
    @test collect(keys(config.output)) == ["lateral", "vertical", "path"]
end

@testset "checkdims" begin
    @test_throws AssertionError Wflow.checkdims(("z", "lat", "time"))
    @test_throws AssertionError Wflow.checkdims(("z", "lat"))
    @test Wflow.checkdims(("lon", "lat", "time")) == ("lon", "lat", "time")
    @test Wflow.checkdims(("time", "lon", "lat")) == ("time", "lon", "lat")
    @test Wflow.checkdims(("lat", "lon", "time")) == ("lat", "lon", "time")
    @test Wflow.checkdims(("time", "lat", "lon")) == ("time", "lat", "lon")
end

@testset "timecycles" begin
    @test Wflow.timecycles([Date(2020, 4, 21), Date(2020, 10, 21)]) == [(4, 21), (10, 21)]
    @test_throws ErrorException Wflow.timecycles([Date(2020, 4, 21), Date(2021, 10, 21)])
    @test_throws ErrorException Wflow.timecycles(collect(1:400))
    @test Wflow.timecycles(collect(1:12)) == collect(zip(1:12, fill(1, 12)))
    @test Wflow.timecycles(collect(1:366)) ==
          monthday.(Date(2000, 1, 1):Day(1):Date(2000, 12, 31))

    @test Wflow.monthday_passed((1,1), (1,1))  # same day
    @test Wflow.monthday_passed((1,2), (1,1))  # day later
    @test Wflow.monthday_passed((2,1), (1,1))  # month later
    @test Wflow.monthday_passed((2,2), (1,1))  # month and day later
    @test !Wflow.monthday_passed((2,1), (2,2))  # day before
    @test !Wflow.monthday_passed((1,2), (2,2))  # month before
    @test !Wflow.monthday_passed((1,1), (2,2))  # day and month before
end

# test reading and setting of warm states (reinit=false)
# modify existing config and initialize model with warm states
@test config.model.reinit
config["model"]["reinit"] = false
@test !config.model.reinit
model = Wflow.initialize_sbm_model(config)

@unpack vertical, clock, reader, writer = model

@testset "output and state names" begin
    ncdims = ("lon", "lat", "layer", "time")
    @test dimnames(writer.dataset["ustorelayerdepth"]) == ncdims
    ncvars = [k for k in keys(writer.dataset) if !in(k, ncdims)]
    @test "snow" in ncvars
    @test "q_river" in ncvars
    @test "q_land" in ncvars
    @test length(writer.state_parameters) == 14
end

@testset "warm states" begin
    @test Wflow.param(model, "lateral.river.reservoir.volume")[2] ≈ 2.7393418e7
    @test Wflow.param(model, "vertical.satwaterdepth")[41120] ≈ 201.51429748535156
    @test Wflow.param(model, "vertical.snow")[41120] ≈ 4.21874475479126
    @test Wflow.param(model, "vertical.tsoil")[41120] ≈ -1.9285825490951538
    @test Wflow.param(model, "vertical.ustorelayerdepth")[1][1] ≈ 16.73013687133789
    @test Wflow.param(model, "vertical.snowwater")[41120] ≈ 0.42188167572021484
    @test Wflow.param(model, "vertical.canopystorage")[1] ≈ 0.0
    @test Wflow.param(model, "lateral.subsurface.ssf")[39308] ≈ 1.1614302208e11
    @test Wflow.param(model, "lateral.river.q")[5501] ≈ 111.46229553222656
    @test Wflow.param(model, "lateral.river.h")[5501] ≈ 8.555977821350098
    @test Wflow.param(model, "lateral.land.q")[39626] ≈ 1.0575231313705444
    @test Wflow.param(model, "lateral.land.h")[39626] ≈ 0.03456519544124603
end

@testset "reducer" begin
    V = [6, 5, 4, 1]
    @test Wflow.reducerfunction("maximum")(V) == 6
    @test Wflow.reducerfunction("mean")(V) == 4
    @test Wflow.reducerfunction("median")(V) == 4.5
    @test Wflow.reducerfunction("first")(V) == 6
    @test Wflow.reducerfunction("last")(V) == 1
    @test_throws ErrorException Wflow.reducerfunction("other")
end

@testset "network" begin
    @unpack network = model
    @unpack indices, reverse_indices = model.network.land
    # test if the reverse index reverses the index
    linear_index = 100
    cartesian_index = indices[linear_index]
    @test cartesian_index === CartesianIndex(115, 6)
    @test reverse_indices[cartesian_index] === linear_index
end

Wflow.close_files(model, delete_output = false)

@testset "NetCDF creation" begin
    path = Base.Filesystem.tempname()
    _ = Wflow.create_tracked_netcdf(path)
    # safe to open the same path twice
    ds = Wflow.create_tracked_netcdf(path)
    close(ds)  # path is removed on process exit
end


text:
function merge_sort!(A, start, stop)
    if start ≥ stop
        return
    end

    middle = floor(Int, (start + stop) / 2)
    merge_sort!(A, start, middle)
    merge_sort!(A, middle + 1, stop)
    merge!(A, start, middle, stop)
end

function merge!(A, start, middle, stop)
    n₁ = middle - start + 1
    n₂ = stop - middle

    L = zeros(Float64, (n₁+1,))
    R = zeros(Float64, (n₂+1,))

    for i in 1:n₁
        L[i] = A[start + i - 1]
    end

    for j in 1:n₂
        R[j] = A[middle + j]
    end

    L[n₁ + 1] = Inf
    R[n₂ + 1] = Inf

    i = j = 1

    for k in start:stop
        if L[i] ≤ R[j]
            A[k] = L[i]
            i += 1
        else
            A[k] = R[j]
            j += 1
        end
    end
end

max_stars_repo_name: blallen/CodePractice

max_issues_repo_name: blallen/CodePractice

max_forks_repo_name: blallen/CodePractice


text:
#!/usr/bin/env julia

Base.eval(:(have_color=true))

function die(msg)
    print_with_color(:red, STDERR, "ERROR: " * msg * "\n")
    exit(1)
end

if length(ARGS) != 1
    println("Usage: ", PROGRAM_FILE, "<folder>")
    println("    where <folder> is the output of a run in Janus.")
    println("    The program will output in stdout a csv with the")
    println("    temporal correlation of the configurations.")
    quit()
end

folder = ARGS[1]

if folder[end] == '/'
    thefolder = folder |> dirname |> basename
else
    thefolder = folder |> basename
end

# Older format
# parsedfoldername = match(r"T1_(0.\d)_T2_(0.\d)_tw1_(\d+)_replica_(\d+)", thefolder)

parsedfoldername = match(r"T1_(0.\d)_T2_(0.\d)_tw1_(\d+)_\w+_replica_(\d+)", thefolder)
if parsedfoldername == nothing
    die("Couldn't parse $folder name.")
else
    T1, T2, tw1, replica = parsedfoldername.captures
    T1                   = parse(Float64 , T1)
    T2                   = parse(Float64 , T2)
    tw1                  = parse(Int64   , tw1)
    replica              = parse(Int64   , replica)
end

if !isfile(folder * "/S00001/run.log")
    die("No run.log in folder.")
end

Jhash = @sprintf("%d",read(folder * "/S00001/conf_Js") |> hash)

cd(folder * "/S00001/confs")

confs = readdir(".")

isconf(f) = ismatch(r"conf_[0-9A-F]{10}", f)

if !all(map(isconf, confs))
    warn("Some files don't match conf_\\x{10}, filtering them.")
    filter!(isconf, confs)
end

function uncompact(load)
    N = length(load)
    data = BitArray(N*8)
    for ib in 1:N
        is = 8(ib-1)
        for j in 0:7
            data[1+is+j] = (load[ib]>>j) & 0x1
        end
    end
    return data
end

function conf2data(file)
    data = Dict()
    open(file,"r") do f
        data["itime"]           = read(f,Cint)
        data["Dt"]              = read(f,Clonglong)
        data["t"]               = read(f,Clonglong) # The MC step
        data["timewindow"]      = read(f,Clonglong)
        data["flagJ"]           = read(f,Cint)
        data["flagS"]           = read(f,Cint)
        data["flagh"]           = read(f,Cint)

        data["L"]               = read(f,Cint)
        data["β"]               = read(f,Cdouble)
        data["H"]               = read(f,Cdouble)
        data["tHon"]            = read(f,Culonglong)

        data["ir"]              = read(f,Cint)
        data["nfiles_dir"]      = read(f,Cint)

        data["energy_sigma"]    = read(f,Cdouble)
        data["energy_tau"]      = read(f,Cdouble)
        data["energymag_sigma"] = read(f,Cdouble)
        data["energymag_tau"]   = read(f,Cdouble)
        data["mag_sigma"]       = read(f,Cdouble)
        data["mag_tau"]         = read(f,Cdouble)

        num                     = data["L"]^3 ÷ 8
        data["spins_σ"]         = uncompact(read(f,Cchar,num))
        data["spins_τ"]         = uncompact(read(f,Cchar,num))
        @assert eof(f)
    end
    return data
end

function getC(bconfA, bconfB)
    spinsA_σ = bconfA["spins_σ"]
    spinsA_τ = bconfA["spins_τ"]
    spinsB_σ = bconfB["spins_σ"]
    spinsB_τ = bconfB["spins_τ"]
    N = length(spinsA_σ)
    differing_spins_σ = map(≠,spinsA_σ, spinsB_σ) |> sum
    differing_spins_τ = map(≠,spinsA_τ, spinsB_τ) |> sum
    # C = 1*N↑↑ + (-1)*N↑↓
    Cσ = (N-differing_spins_σ) - differing_spins_σ
    Cσ /= N
    Cτ = (N-differing_spins_τ) - differing_spins_τ
    Cτ /= N
    return (Cσ+Cτ)/2
end

function rectify_times!(array)
    L = length(array)
    for i in 2:L
        if array[i] == array[i-1]
            for j in i:L
                array[j] += 1
            end
        end
    end
end

# This definition was used when the list of times didn't include all
# the tw1 combinations. Not, it has been superseeded by one that
# ignores its argument.
# function generatetimes(tw1)
#     tws_A = logspace(log10(1),log10(tw1),100)
#     tws_A = map(x->floor(Int64,x), tws_A)
#     rectify_times!(tws_A)

#     tws_B = logspace(log10(tw1+1),log10(2e8),100)
#     tws_B = map(x->floor(Int64,x), tws_B)
#     rectify_times!(tws_B)

#     tws_C = logspace(log10(2e8+1),log10(3e8+10),100)
#     tws_C = map(x->floor(Int64,x), tws_C)
#     rectify_times!(tws_C)

#     tws = unique([tws_A;tws_B;tws_C]) |> sort
#     ts = logspace(4,8,10)
#     ts = map(x->floor(Int64,x), ts) |> sort

#     return tws, ts
# end

function generatetimes(tw1)
    const tw1list = map(Int64,[1e4,1e5,1e6,1e7,1e8])
    tws = logspace(log10(1),log10(3e8+10),100)

    for tw1 in tw1list
        newlogs = logspace(log10(tw1), log10(3e8), 50)
        tws = vcat(tws, newlogs)
    end

    endlogs = logspace(log10(2e8), log10(3e8), 50)
    tws = vcat(tws, endlogs)

    tws = map(x->floor(Int64,x), tws) |> sort
    rectify_times!(tws)

    tws = vcat(tw1list, tws) |> unique |> sort
    ts = logspace(4,8,10)
    ts = map(x->floor(Int64,x), ts) |> sort

    return tws, ts
end

function search_E_flukes(mc_list, E_list, tw1)
    crop1 = findfirst(x->x>tw1, mc_list)
    crop2 = findfirst(x->x>=2e8, mc_list)
    mc = mc_list[crop1+50:crop2-1]
    E = E_list[crop1+50:crop2-1]
    threshold = mean(E)
    if any(E .> mean(E) + 5std(E)) || any(E .< mean(E) - 5std(E))
        warn("Possible FPGA corruption in $folder.")
    end
end

function measureandprint(confs, tw1, Jhash)
    confs = sort(confs)

    twlist, tlist = generatetimes(tw1)

    ℓ   = length(confs)
    ntw = length(twlist)
    nt  = length(tlist)

    bconfs = map(conf2data, confs)

    L = bconfs[1]["L"]
    Tlist = [1/c["β"] for c in bconfs]
    readedT1, readedT2 = maximum(Tlist), minimum(Tlist)
    if readedT1 ≉ T1
        die("Readed T1 ($readedT1) is not ≈ folder T1 ($T1)")
    end
    if readedT2 ≉ T2
        die("Readed T2 ($T2) is not ≈ folder T2 ($T2)")
    end

    timeslist = [c["t"] for c in bconfs]
    realmag_σ_list = [c["mag_sigma"] for c in bconfs]
    realmag_τ_list = [c["mag_tau"] for c in bconfs]
    E_σ_list = [c["energy_sigma"] for c in bconfs]
    E_τ_list = [c["energy_tau"] for c in bconfs]

    E_list = 0.5*(E_σ_list + E_τ_list)
    M_list = 0.5*(realmag_σ_list + realmag_τ_list)

    # Consistency check
    measmag_σ_list = [(2sum(c["spins_σ"])-L^3)/L^3 for c in bconfs]
    measmag_τ_list = [(2sum(c["spins_τ"])-L^3)/L^3 for c in bconfs]
    if measmag_σ_list ≉ realmag_σ_list
        die("Readed magnetization is not ≈ calculated magnetization (σ spins)")
    end
    if measmag_τ_list ≉ realmag_τ_list
        die("Readed magnetization is not ≈ calculated magnetization (τ spins)")
    end

    # Search for FPGA flukes
    search_E_flukes(timeslist, E_list, tw1)

    for i in 1:ℓ
        tw = bconfs[i]["t"]
        tw ∉ twlist && continue
        for j in i+1:ℓ
            tw_plus_t = bconfs[j]["t"]
            t = tw_plus_t - tw
            t ∉ tlist && continue
            # Wohoo! measure!
            T = 1/bconfs[i]["β"]
            C = getC(bconfs[i],bconfs[j])
            println("$L,$tw1,$tw,$T,\"$Jhash\",\"Energy\",0,$(E_list[i]),$T1,$T2")
            println("$L,$tw1,$tw,$T,\"$Jhash\",\"Magnetization\",0,$(M_list[i]),$T1,$T2")
            println("$L,$tw1,$tw,$T,\"$Jhash\",\"Correlation\",$t,$C,$T1,$T2")
        end
    end

end

println("\"L\",\"tw1\",\"mc\",\"T\",\"Jhash\",\"observable\",\"parameter\",\"value\",\"T1\",\"T2\"")
measureandprint(confs, tw1, Jhash)

max_stars_repo_name: RedPointyJackson/tfg

max_issues_repo_name: RedPointyJackson/tfg

max_forks_repo_name: RedPointyJackson/tfg


text:
using CodecZlib, Requires

abstract type ArchiveFormat{S} end

decompressor_stream(::Type{ArchiveFormat{:GZIP}}, io::IOStream) = GzipDecompressorStream(io)

# Handles the case when we are loading a supported archive format but the codec is not loaded
function decompressor_stream(::Type{ArchiveFormat{S}}, io::IOStream) where S
    throw(CodecMissingError(S))
end

# called by __init__() to enable support for optional codecs
function define_optional_decompressor_streams()
    @require CodecBzip2="523fee87-0ab8-5b00-afb7-3ecf72e48cfd" begin
        decompressor_stream(::Type{ArchiveFormat{:BZIP2}}, io::IOStream) = CodecBzip2.Bzip2DecompressorStream(io)
    end
    @require CodecXz="ba30903b-d9e8-5048-a5ec-d1f5b0d4b47b" begin
        decompressor_stream(::Type{ArchiveFormat{:XZ}}, io::IOStream) = CodecXz.XzDecompressorStream(io)
    end
end

# decompress the stream if it's compressed (by a supported codec)
function decompress(io)
    format = FileIO.detect_compressor(io, formats=["GZIP", "BZIP2", "XZ"])
    if format !== nothing
        return decompressor_stream(ArchiveFormat{Symbol(format)}, io)
    else
        return io # not compressed
    end
end


text:
"random number between a and b"
function randin(a::Real, b::Real)
    -Inf < a ≤ b < Inf || DomainError()
    a + (b - a) * rand()
end

"""
Weighted random integer.
Given weights = (w1, w2, ...), returns a random integer i with probability proportional to w_i.
The weights should be 0 ≤ w_i < ∞, and do not need to sum to 1.
"""
function wrand(weights)
    wr = rand() * sum(w for w in weights)
    for (i,w) in enumerate(weights)
        @assert 0 ≤ w < Inf
        wr -= w
        if wr ≤ 0
            return i
        end
    end
end

"""
Returns a random vector x = (x1, x2, ..., xn), such that
xi ≥ 0 and sum(xi) = 1
"""
function randsimplex(n::Int)
    # https://cs.stackexchange.com/questions/3227/uniform-sampling-from-a-simplex
    @assert n > 0
    diff(sort([0; rand(n - 1); 1]))
end

"""
Returns a random vector x = x0 * v[0] + x1 * v[1] + ... + xn * v[n]
where xi ≥ 0 and sum(xi) = 1
"""
function randsimplex(V::Vector{Vector{Float64}})
    @assert allsame(length.(V))
    X = randsimplex(length(V))
    sum(x*v for (x,v) in zip(X,V))
end


text:
export vect_to_ind, ind_to_vect, power_flow, cost, markovdecision

using MAT
using JuMP
using CPLEX
# include("savecosts.jl")

"""
    vect_to_ind(vec, m, n)

This function converts the vector vec into a single integer betwen 1 and n^m,
where `m` is the length of the vector, `n` is the number of possibilities
in each entry of the vector and each entry of vec is between 1 and `n`.
"""
function vect_to_ind(vec, m, n)
	@assert (length(vec) == m)
	vec=vec-ones(m)
	return sum(vec[i]*n^(m-i) for i = 1:m)+1
end

"""
    ind_to_vect(ind, m, n)

This function converts the integer `ind` into the corresponding vector of length `m`
such that `ind_to_vect(vect_to_ind(vec, m, n), m, n) = vec`.
Here `m` is the length of the vector, `n` is the number of possibilities
in each entry of the vector and `ind` is an integer between 1 and n^m.
"""
function ind_to_vect(ind, m, n)
	ind = ind-1
	vect = zeros(m)
	for i = 1:m
		to_div = ind - sum(vect[j]*n^(m-j) for j = 1:m)
		vect[i] = div(to_div, n^(m-i))
	end
	return vect+ones(m)
end

"""
    power_flow(p, r, x, I_up, Vmin, Vmax)

Computes the values of branch and nodal electrical quantities for given injection `p`
and returns the cost associated to the unfeasibility of the voltages and currents.
"""
function power_flow(power, r, x, I_up, Vmin, Vmax)
	Cunfeas = 8300
	solver = CplexSolver(
    	CPX_PARAM_SCRIND=0,   # Verbose solver output
    )

	m = Model(solver=solver)

	@variable(m, v[1:5]>=0)
	JuMP.fix(v[1], 0.0529)
	@variable(m, l[2:5]>=0)
	@variable(m, P[1:5])
	@variable(m, Q[1:5])
	@variable(m, q[1:5])
	@variable(m, p[1:5])
	for i = 1:5
		setlowerbound(q[i], -0.01)
		setupperbound(q[i], 0.01)
	end
	for i = 2:5
		JuMP.fix(p[i], power[i-1])
	end

	@constraint(m, P[1] == 0)
	@constraint(m, Q[1] == 0)
	@constraint(m, [i=1:5], P[i] == p[i] + sum(P[j]-r[j-1]*l[j] for j = i+1:5))
	@constraint(m, [i=1:5], Q[i] == q[i] + sum(Q[j]-x[j-1]*l[j] for j = i+1:5))
	@constraint(m, [i=2:5], v[i-1] == v[i] -2*(r[i-1]*P[i]+x[i-1]*Q[i])+(x[i-1]^2+r[i-1]^2)*l[i])
	@constraint(m, [i=2:5], v[i]*l[i] >= P[i]^2 + Q[i]^2)

	@objective(m, Min, 0)

	status = solve(m)
	if status == :Optimal
		v_val = getvalue(v)
		l_val = getvalue(l)
		cost_unfeas = 0
		for j = 2:4
			if l_val[j] > I_up[j-1]
				cost_unfeas = cost_unfeas + Cunfeas*(l_val[j]-I_up[j-1])
			end
		end
		for i = 1:4
			if v_val[i] > Vmax[i]
				cost_unfeas = cost_unfeas + Cunfeas*(v_val[i]-Vmax[i])
			end
			if v_val[i] < Vmin[i]
				cost_unfeas = cost_unfeas + Cunfeas*(Vmin[i]-v_val[i])
			end
		end
	else
		cost_unfeas = 1e9
	end

	return cost_unfeas
end

"""
    cost(a, k, t, data)

Computes the cost for a given action, stage and time step.
`a` is an integer from 1:54 representing the action (+10%, +0% or -10% for each household
 and the activation level for the flexible load).
`k` is an integer from 1:9 representing the state (demand or solar, +10%, +0% or -10%).
`t` is an integer from 1:169 representing the time step (one week, hourly resolution).
"""
function cost(a, k, t, data)
	Ccurt = 300
	Cflex = 1
	forecast_demand = data["pDemand"][t, :]
	forecast_solar = data["pSolarMax"][t, :]

	# Apply uncertainty given by k to demand and solar production
	state = ind_to_vect(k, 3, 3)
	uncertainty = state[2:3]
	uncertainty[find(iszero, uncertainty-1)] = 1.1
	uncertainty[find(iszero, uncertainty-2)] = 1
	uncertainty[find(iszero, uncertainty-3)] = 0.9
	forecast_demand = uncertainty[1]*forecast_demand
	forecast_solar = uncertainty[2]*forecast_solar

    # Apply curtailment action given by k
	actions = ind_to_vect(a, 4, 3)
	curtailment = actions[2:4]
	curtailment[find(iszero, curtailment-1)] = 0
	curtailment[find(iszero, curtailment-2)] = -0.5
	curtailment[find(iszero, curtailment-3)] = -1
	power = zeros(4)
	power[1:3] = forecast_solar[1:3] + forecast_solar[1:3].*curtailment - forecast_demand[1:3]
	power[4] = actions[1] == 1 ? 0 : forecast_demand[4]
	cost_flex = actions[1] == 1 ? 0 : Cflex
	cost_curt = sum(-Ccurt*forecast_solar[1:3].*curtailment)

	# Get data
	r = data["r"]
	x = data["x"]
	I_up = data["I_up"]
	Vmin = data["Vmin"]
	Vmax = data["Vmax"]

	cost_unfeas = power_flow(power, r, x, I_up, Vmin, Vmax)
	tot_cost = cost_curt + cost_unfeas + cost_flex
end

"""
    markovdecision(data)

Computes the optimal action and the value function using the value iteration algorithm
data is the dictionnary given by the call `matread(\"data_4nodes.mat\")`.
"""
function markovdecision(data)
	pDemand = data["pDemand"]
	nb_steps = size(pDemand)[1]
	nb_nodes = size(pDemand)[2]
	nb_stages = 18
	cost_akt = zeros(2*3^3, nb_stages, nb_steps)
	big_P = zeros(18,18)
	P = zeros(9,9)
	for i = 1:9
		for j = 1:9
			val = 1
			vec1 = ind_to_vect(i, 3, 3)[2:3]
			vec2 = ind_to_vect(j, 3, 3)[2:3]
			for k in abs.(vec1-vec2)
				val = val*(k == 0 ? 0.5 : 0.25)
			end
			P[i,j] = val
		end
	end
	V = zeros(nb_stages, nb_steps)
	opt_actions = zeros(nb_stages, nb_steps)
	for t = nb_steps-1:-1:1
		for k = 1:nb_stages
			state = ind_to_vect(k,3,3)[1]
			nb_actions = state[1] == 1 ? 2*3^3 : 3^3
			all_possibilities = zeros(nb_actions)
			for a = 1:nb_actions
				actions = ind_to_vect(a,4,3)
				if actions[1] == 1
					big_P = [P zeros(9,9);zeros(9,9) P]
				elseif actions[1] == 2
					big_P = [zeros(9,9) P;zeros(9,9) zeros(9,9)]
				end
				cost_a = cost(a, k, t, data)
				cost_akt[a,k,t] = cost_a
				all_possibilities[a] = cost_a + big_P[k,:]'*V[:,t+1]
			end
			V[k,t] = minimum(all_possibilities)
			opt_actions[k,t] = indmin(all_possibilities)
		end
		println("Iteration t = ", t)
	end
	savecost(cost_akt)
	return (V, opt_actions)
end

function run()
    data = matread("./data/data_4nodes.mat")
    (V, opt_actions) = markovdecision(data)
    @show V
    @show opt_actions
end

max_stars_repo_name: qlete/ANManagement

max_issues_repo_name: qlete/ANManagement

max_forks_repo_name: qlete/ANManagement


text:
#=
This file is auto-generated. Do not edit.
=#

"""Data Structure for the technical parameters of thermal generation technologies."""
mutable struct TechThermal <: TechnicalParams
    rating::Float64  # Thermal limited MVA Power Output of the unit. <= Capacity
    primemover::PrimeMovers  # PrimeMover Technology according to EIA 923
    fuel::ThermalFuels  # PrimeMover Fuel according to EIA 923
    activepowerlimits::Min_Max
    reactivepowerlimits::Union{Nothing, Min_Max}
    ramplimits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    timelimits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    internal::InfrastructureSystemsInternal
end

function TechThermal(rating, primemover, fuel, activepowerlimits, reactivepowerlimits, ramplimits, timelimits, )
    TechThermal(rating, primemover, fuel, activepowerlimits, reactivepowerlimits, ramplimits, timelimits, InfrastructureSystemsInternal())
end

function TechThermal(; rating, primemover, fuel, activepowerlimits, reactivepowerlimits, ramplimits, timelimits, )
    TechThermal(rating, primemover, fuel, activepowerlimits, reactivepowerlimits, ramplimits, timelimits, )
end

# Constructor for demo purposes; non-functional.

function TechThermal(::Nothing)
    TechThermal(;
        rating=0.0,
        primemover=OT::PrimeMovers,
        fuel=OTHER::ThermalFuels,
        activepowerlimits=(min=0.0, max=0.0),
        reactivepowerlimits=nothing,
        ramplimits=nothing,
        timelimits=nothing,
    )
end

"""Get TechThermal rating."""
get_rating(value::TechThermal) = value.rating
"""Get TechThermal primemover."""
get_primemover(value::TechThermal) = value.primemover
"""Get TechThermal fuel."""
get_fuel(value::TechThermal) = value.fuel
"""Get TechThermal activepowerlimits."""
get_activepowerlimits(value::TechThermal) = value.activepowerlimits
"""Get TechThermal reactivepowerlimits."""
get_reactivepowerlimits(value::TechThermal) = value.reactivepowerlimits
"""Get TechThermal ramplimits."""
get_ramplimits(value::TechThermal) = value.ramplimits
"""Get TechThermal timelimits."""
get_timelimits(value::TechThermal) = value.timelimits
"""Get TechThermal internal."""
get_internal(value::TechThermal) = value.internal


text:
# simulates data from prior, trains and tests the net, and returns
# the trained net and the information for transforming the inputs

using Statistics, Flux
using Base.Iterators

function MakeNeuralMoments(model::SNMmodel;TrainTestSize=1, Epochs=1000)
    data = 0.0
    datadesign = 0.0
    nParams = size(model.lb,1)
    # training and testing
    if (TrainTestSize == 1) TrainTestSize = Int64(2*nParams*1e4); end # use a default size if none provided
    params = zeros(TrainTestSize,nParams)
    statistics = zeros(TrainTestSize,size(model.auxstat(model.lb,1)[1],1))
    Threads.@threads for s = 1:TrainTestSize
        ok = false
        θ = model.priordraw()
        W = (model.auxstat(θ,1))[1]
        # repeat draw if necessary
        while any(isnan.(W))
            θ = model.priordraw()
            W = model.auxstat(θ,1)[1]
        end
        params[s,:] = θ
        statistics[s,:] = W
    end
    # transform stats to robustify against outliers
    q50 = zeros(size(statistics,2))
    q01 = similar(q50)
    q99 = similar(q50)
    iqr = similar(q50)
    for i = 1:size(statistics,2)
        q = quantile(statistics[:,i],[0.01, 0.25, 0.5, 0.75, 0.99])
        q01[i] = q[1]
        q50[i] = q[3]
        q99[i] = q[5]
        iqr[i] = q[4] - q[2]
    end
    nninfo = (q01, q50, q99, iqr)
    statistics = TransformStats(statistics, nninfo)
    # train net
    TrainingProportion = 0.5 # size of training/testing
    params = Float32.(params)
    s = Float32.(std(params, dims=1)')
    statistics = Float32.(statistics)
    trainsize = Int(TrainingProportion*TrainTestSize)
    yin = params[1:trainsize, :]'
    yout = params[trainsize+1:end, :]'
    xin = statistics[1:trainsize, :]'
    xout = statistics[trainsize+1:end, :]'
    # define the neural net
    nStats = size(xin,1)
    NNmodel = Chain(
        Dense(nStats, 10*nParams, tanh),
        Dense(10*nParams, 3*nParams, tanh),
        Dense(3*nParams, nParams)
    )
    loss(x,y) = Flux.huber_loss(NNmodel(x)./s, y./s; δ=0.1) # Define the loss function
    # monitor training
    function monitor(e)
        println("epoch $(lpad(e, 4)): (training) loss = $(round(loss(xin,yin); digits=4)) (testing) loss = $(round(loss(xout,yout); digits=4))| ")
    end
    # do the training
    bestsofar = 1.0e10
    pred = 0.0 # define it here to have it outside the for loop
    batches = [(xin[:,ind],yin[:,ind])  for ind in partition(1:size(yin,2), 50)]
    bestmodel = 0.0
    for i = 1:Epochs
        if i < 20
            opt = Momentum() # the optimizer
        else
            opt = ADAMW() # the optimizer
        end
        Flux.train!(loss, Flux.params(NNmodel), batches, opt)
        current = loss(xout,yout)
        if current < bestsofar
            bestsofar = current
            bestmodel = NNmodel
            xx = xout
            yy = yout
            println("________________________________________________________________________________________________")
            monitor(i)
            pred = NNmodel(xx)
            error = yy .- pred
            results = [pred;error]
            rmse = sqrt.(mean(error.^Float32(2.0),dims=2))
            println(" ")
            println("RMSE for model parameters ")
            prettyprint(reshape(round.(rmse,digits=3),1,nParams))
            println(" ")
            println("dstats prediction of parameters:")
            dstats(pred'; short=true);
            println(" ")
            println("dstats prediction error:")
            dstats(error'; short=true);
        end
    end
    bestmodel, nninfo
end


text:
struct MSTState{C <: Codon}
    parents::Vector{C}
    ranks::Vector{Int}
end

function MSTState{C}() where C <: Codon
    MSTState{C}(collect((C(i) for i in UInt64(0):UInt64(63))), zeros(Int, 64))
end

@inline function reset!(x::MSTState{C}) where C <: Codon
    fill!(x.ranks, 0)
    copy!(x.parents, (C(i) for i in UInt64(0):UInt64(63)))
end

cdn2i(c::Codon) = UInt64(c) + 1

@inline function parentof(vertex::C, state::MSTState{C}) where C <: Codon
    @inbounds return state.parents[cdn2i(vertex)]
end

@inline function setparentof(vertex::C, parent::C, state::MSTState{C}) where C <: Codon
    @inbounds state.parents[cdn2i(vertex)] = parent
end

@inline function rankof(vertex::C, state::MSTState{C}) where C <: Codon
    @inbounds return state.ranks[cdn2i(vertex)]
end

@inline function incrank(vertex::C, state::MSTState{C}) where C <: Codon
    @inbounds state.ranks[cdn2i(vertex)] += 1
end

function find(vertex::C, state::MSTState{C}) where C <: Codon
    p = parentof(vertex, state)
    if p != vertex
        setparentof(vertex, find(p, state), state)
    end
    return parentof(vertex, state)
end

function union(vertexa::C, vertexb::C, state::MSTState{C}) where C <: Codon
    roota = find(vertexa, state)
    rootb = find(vertexb, state)
    if roota != rootb
        ranka = rankof(roota, state)
        rankb = rankof(rootb, state)
        if ranka > rankb
            setparentof(rootb, roota, state)
        else
            setparentof(roota, rootb, state)
            if ranka == rankb
                incrank(rootb, state)
            end
        end
    end
end

function mst(codongraph::CodonGraph, state::MSTState)
    DS = DN = 0
    reset!(state)
    for (v1, v2, DS_i, DN_i) in codongraph
        if find(v1, state) != find(v2, state)
            union(v1, v2, state)
            DS += DS_i
            DN += DN_i
        end
    end
    return DS, DN
end


text:
module PassiveTracerFlows

using
  CUDA,
  Reexport

@reexport using FourierFlows

include("traceradvectiondiffusion.jl")

@reexport using PassiveTracerFlows.TracerAdvectionDiffusion

end # module

max_stars_repo_name: FourierFlows/PassiveTracerFlows

max_issues_repo_name: FourierFlows/PassiveTracerFlows

max_forks_repo_name: FourierFlows/PassiveTracerFlows


text:
using Rewrite
using Test

include("patterns.jl")
include("rules.jl")


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, Pkg

name = "ruby"
version = v"2.7.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource(
        "https://cache.ruby-lang.org/pub/ruby/2.7/ruby-2.7.1.tar.gz",
        "d418483bdd0000576c1370571121a6eb24582116db0b7bb2005e90e250eae418",
    ),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/ruby-*/
# otherwise readline and therefore ruby is held back
apk del python
# living on the edge...
sed -i -e 's/v[[:digit:]]\..*\//edge\//g' /etc/apk/repositories
apk update
# these return an error code, but work anyways ¯\_(ツ)_/¯
apk upgrade --update-cache --available --latest || true
apk add ruby-full || true
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --with-baseruby=/usr/bin/ruby --enable-shared
make -j${nproc}
make install
${bindir}/ruby -e 'puts $:' | sed "s:${prefix}/::g" > ${prefix}/RUBYLIB
for bin in $(grep -rl "${bindir}/ruby" ${bindir})
do
    sed -i "s:${bindir}/ruby:/usr/bin/env ruby:" ${bin}
done
"""

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

# Unfortunately Windows and Mac don't seem to cross-compile out of the box.
# It also seems like some of the libraries don't support case-insensitive file systems,
# so this might be a problem when trying to get those to work.
platforms = filter(p -> Sys.islinux(p) || Sys.isfreebsd(p), supported_platforms())
# TODO: fix armv7l musl. Probably an upstream issue though
filter!(!=(Platform("armv7l", "linux"; libc="musl")), platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libruby", :libruby),
    ExecutableProduct("ruby", :ruby),
    # These are actually ruby scripts
    ExecutableProduct("bundle", :bundle),
    ExecutableProduct("bundler", :bundler),
    ExecutableProduct("erb", :erb),
    ExecutableProduct("gem", :gem),
    ExecutableProduct("irb", :irb),
    ExecutableProduct("racc", :racc),
    ExecutableProduct("racc2y", :racc2y),
    ExecutableProduct("rake", :rake),
    ExecutableProduct("rdoc", :rdoc),
    ExecutableProduct("ri", :ri),
    # these are all subdirs of ${prefix} that need to be in the env var RUBYLIB for ruby to
    # work properly
    FileProduct("RUBYLIB", :RUBYLIB),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
    Dependency("Libiconv_jll"),
    Dependency("OpenSSL_jll"),
    Dependency("Readline_jll"),
    Dependency("Zlib_jll"),
    Dependency("Gdbm_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
	struct ExtractScalar{T}
		c::T
		s::T
	end

Extracts a numerical value, centred by subtracting `c` and scaled by multiplying by `s`.
Strings are converted to numbers.

The extractor returns `ArrayNode{Matrix{Union{Missing, Int64}},Nothing}` or it subtypes.
If passed `missing`, it extracts missing values which Mill understands and can work with.

It can be created also using `extractscalar(Float32, 5, 2)`

# Example
```jldoctest
julia> ExtractScalar(Float32, 2, 3)(1)
1×1 Mill.ArrayNode{Array{Float32,2},Nothing}:
 -3.0

julia> ExtractScalar(Float32, 2, 3)(missing)
1×1 Mill.ArrayNode{Array{Missing,2},Nothing}:
 missing
```
"""
struct ExtractScalar{T} <: AbstractExtractor
	c::T
	s::T
end

ExtractScalar(T = Float32) = ExtractScalar(zero(T), one(T))
ExtractScalar(T, c, s) = ExtractScalar(T(c), T(s))

extractscalar(::Type{T}, m = zero(T), s = one(T)) where {T<:Number} = ExtractScalar(T, m, s)
function extractscalar(::Type{T}, e::Entry) where {T<:Number}
	if unify_types(e) <: AbstractString
		values = parse.(T, keys(e.counts))
	else
		values = keys(e.counts)
	end
	min_val = minimum(values)
	max_val = maximum(values)
	c = min_val
	s = max_val == min_val ? 1 : 1 / (max_val - min_val)
	ExtractScalar(Float32(c), Float32(s))
end

(s::ExtractScalar{T})(v::W) where {T,W<:Union{Missing, Nothing}} = ArrayNode(fill(missing,(1,1)))
(s::ExtractScalar{T})(v::W) where {T,W<:ExtractEmpty} = ArrayNode(fill(zero(T),(1,0)))
(s::ExtractScalar{T})(v::Number) where {T} = ArrayNode(s.s .* (fill(T(v),1,1) .- s.c))
(s::ExtractScalar{T})(v::AbstractString) where{T} = s((tryparse(T,v)))
(s::ExtractScalar)(v) = s(missing)

Base.length(e::ExtractScalar) = 1

# data type has different hashes for each patch version of julia
# see https://discourse.julialang.org/t/datatype-hash-differs-per-patch-version/48827
Base.hash(e::ExtractScalar{T}, h::UInt) where {T} = hash((e.c, e.s), h)
Base.:(==)(e1::ExtractScalar{T}, e2::ExtractScalar{T}) where {T} = e1.c === e2.c && e1.s === e2.s


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, Pkg

name = "Antlr4CppRuntime"
version = v"4.9.3"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/antlr/antlr4/archive/refs/tags/$(version).tar.gz",
                  "efe4057d75ab48145d4683100fec7f77d7f87fa258707330cadd1f8e6f7eecae"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/antlr4*/runtime/Cpp
atomic_patch -p1 ../../../patches/forward-cmake-toolchain.patch
mkdir build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release
make -j${nproc}
make install

# Remove large static library
rm "${prefix}/lib/libantlr4-runtime.a"
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct("libantlr4-runtime", :libantlr4_cpp_runtime)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Libuuid_jll", uuid="38a345b3-de98-5d2b-a5d3-14cd9215e700"); platforms=filter(p -> !(Sys.iswindows(p) || Sys.isapple(p)), platforms))
]

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

max_stars_repo_name: azreika/Yggdrasil

max_issues_repo_name: azreika/Yggdrasil

max_forks_repo_name: azreika/Yggdrasil


text:
# select algorithm used to estimate Sep for triangular generalized problems
abstract type TrGSepAlgo end
# look-ahead is standard in LAPACK
struct TrGSepLookAhead <: TrGSepAlgo end
# alternative is Higham's iterative 1-norm estimator
struct TrGSepIter <: TrGSepAlgo end

function _trgsep(A::StridedMatrix{T},B::StridedMatrix{T},
                 D::StridedMatrix{T},E::StridedMatrix{T},
                 algo::TrGSepAlgo = TrGSepLookAhead()) where {T}
    m = checksquare(A)
    n = checksquare(B)

    # C and F are used to hold intermediate results
    C = zeros(T,m,n)
    F = zeros(T,m,n)

    m1 = checksquare(D)
    n1 = checksquare(E)
    (m1 == m && n1 == n) || throw(DimensionMismatch(
        "dimensions of D and E must match A and B"))

    rtyone = one(real(T))
    rtyzero = zero(real(T))
    scale = rtyone
    scaloc = rtyone
    dif = rtyzero
    cone = one(T)
    rdsum = rtyone
    rdscale = rtyzero
    for j=1:n
        for i=m:-1:1
            # build 2x2 problem
            Z = [A[i,i] -B[j,j]; D[i,i] -E[j,j]]
            fZ = lu(Z) # should be _getc2
            rhs = [C[i,j]; F[i,j]]
            if algo isa TrGSepLookAhead
                # from zlatdf
                # set up a RHS b s.t. solution x = Z \ b has large norm
                x = rhs[fZ.p] # ipiv if _getc2
                pmone = -one(T)
                bp = x[1] + cone
                bm = x[1] - cone
                # L-part
                splus = rtyone
                L21 = fZ.L[2,1]
                splus += real(L21' * L21)
                sminu = real(L21' * x[2])
                splus *= real(x[1])
                if splus > sminu
                    x[1] = bp
                elseif sminu > splus
                    x[1] = bm
                else
                    x[1] += pmone
                    pmone = cone
                end
                x[2] -= x[1] * L21
                # U-part
                w = [x[1]; x[2] + cone]
                x[2] -= cone
                splus = zero(real(T))
                sminu = zero(real(T))
                for i=2:-1:1
                    t = cone / fZ.U[i,i]
                    w[i] *= t
                    x[i] *= t
                    if i==1
                        k=2
                        U12 = fZ.U[i,k]
                        w[i] -= w[k] * (U12*t)
                        x[i] -= x[k] * (U12*t)
                    end
                    splus += abs(w[i])
                    sminu += abs(x[i])
                end
                if splus > sminu
                    x .= w
                end
                x = x[fZ.p] # jpiv if _getc2
                rdscale, rdsum = _ssq(x, rdscale, rdsum)
                if rdscale != 0
                    dif = sqrt(2*rtyone*m*n) / (rdscale * sqrt(rdsum))
                end
            else
                error("missing logic for $algo")
            end

            C[i,j] = x[1]
            F[i,j] = x[2]
            # substitute R[i,j], L[i,j] into remaining eq
            if i > 1
                α = -x[1]
                C[1:i-1,j] .+= α * A[1:i-1,i]
                F[1:i-1,j] .+= α * D[1:i-1,i]
            end
            if j < n
                α = x[2]
                C[i,j+1:n] .+= α * B[j,j+1:n]
                F[i,j+1:n] .+= α * E[j,j+1:n]
            end
        end
    end
    return dif
end

# based on ztgsen
"""
    eigvalscond(S::GeneralizedSchur, nsub) => pl, pr

compute approx. reciprocal norms of projectors on left/right subspaces
associated w/ leading `nsub`×`nsub` block of `S`. Use `ordschur` to select
eigenvalues of interest.

An approximate bound on avg. absolute error of associated eigenvalues
is `ϵ * norm(vcat(A,B)) / pl`. See LAPACK documentation for further details.
"""
function eigvalscond(S::GeneralizedSchur{Ty}, nsub) where Ty
    # nsub is m in ztgsen
    n = size(S.T,1)
    rtyone = one(real(Ty))
    rtyzero = zero(real(Ty))
    if nsub == 0 || nsub == n
        return rtyone, rtyone
    end
    # solve S₁₁ R - L S₂₂ = σ S₁₂
    #       T₁₁ R - L T₂₂ = σ T₁₂
    n1 = nsub
    n2 = n - nsub
    R = S.S[1:n1,n1+1:n]
    L = S.T[1:n1,n1+1:n]
    R, L, scale = trsylvester!(S.S[1:n1,1:n1],S.S[n1+1:n,n1+1:n],R,
                               S.T[1:n1,1:n1],S.T[n1+1:n,n1+1:n],L)

    #=
    ## DEVNOTE: this has been taking MINUTES to compile on v1.1.1+
    ## so we live with the allocations for now

    R, L, scale = trsylvester!(view(S.S,1:n1,1:n1),view(S.S,n1+1:n,n1+1:n),R,
                               view(S.T,1:n1,1:n1),view(S.T,n1+1:n,n1+1:n),L)
    =#

    rdscale = rtyzero
    sumsq = rtyone
    rdscale, sumsq = _ssq(vec(R), rdscale, sumsq)
    pl = rdscale * sqrt(sumsq)
    if pl == 0
        pl = rtyone
    else
        pl = scale / (sqrt(scale*scale / pl + pl) * sqrt(pl))
    end
    rdscale = rtyzero
    sumsq = rtyone
    rdscale, sumsq = _ssq(vec(L), rdscale, sumsq)
    pr = rdscale * sqrt(sumsq)
    if pr == 0
        pr = rtyone
    else
        pr = scale / (sqrt(scale*scale / pr + pr) * sqrt(pr))
    end
    pl, pr
end

function subspacesep(S::GeneralizedSchur{Ty}, nsub::Integer;
                     pnorm::Union{Int,Symbol}=:F) where Ty
    n = size(S.S,1)
    rtyone = one(real(Ty))
    n2 = n - nsub
    sep = zeros(real(Ty),2)
    # these could use views for all, but for now that takes minutes to compile
    function fu(X)
        # solve T₁₁ R - R T₂₂ = σ X
        R = reshape(X[1:nsub*n2],nsub,n2)
        L = reshape(X[nsub*n2+1:2*nsub*n2],nsub,n2)
        R, L, s = trsylvester!(S.S[1:nsub,1:nsub],
                               S.S[nsub+1:n,nsub+1:n],
                               R,
                               S.T[1:nsub,1:nsub],
                               S.T[nsub+1:n,nsub+1:n],
                               L)
        scale = s
        return vcat(vec(R),vec(L))
    end
    function fuH(X)
        # solve T₁₁ᴴ R - R T₂₂ᴴ = σ X
        R = reshape(X[1:nsub*n2],nsub,n2)
        L = reshape(X[nsub*n2+1:2*nsub*n2],nsub,n2)
        R, L, s = adjtrsylvester!(S.S[1:nsub,1:nsub],
                                  S.S[nsub+1:n,nsub+1:n],
                                  R,
                                  S.T[1:nsub,1:nsub],
                                  S.T[nsub+1:n,nsub+1:n],
                                  L)
        scale = s
        return vcat(vec(R),vec(L))
    end
    function fl(X)
        # solve T₁₁ R - R T₂₂ = σ X
        R = reshape(X[1:nsub*n2],n2,nsub)
        L = reshape(X[nsub*n2+1:2*nsub*n2],n2,nsub)
        R, L, s = trsylvester!(S.S[nsub+1:n,nsub+1:n],
                               S.S[1:nsub,1:nsub],
                               R,
                               S.T[nsub+1:n,nsub+1:n],
                               S.T[1:nsub,1:nsub],
                               L)
        scale = s
        return vcat(vec(R),vec(L))
    end
    function flH(X)
        # solve T₁₁ᴴ R - R T₂₂ᴴ = σ X
        R = reshape(X[1:nsub*n2],n2,nsub)
        L = reshape(X[nsub*n2+1:2*nsub*n2],n2,nsub)
        R, L, s = adjtrsylvester!(S.S[nsub+1:n,nsub+1:n],
                                  S.S[1:nsub,1:nsub],
                                  R,
                                  S.T[nsub+1:n,nsub+1:n],
                                  S.T[1:nsub,1:nsub],
                                  L)
        scale = s
        return vcat(vec(R),vec(L))
    end
    if pnorm == :F || pnorm == :Frobenius
        sep[1] = _trgsep(view(S.S,1:nsub,1:nsub),
                          view(S.S,nsub+1:n,nsub+1:n),
                          view(S.T,1:nsub,1:nsub),
                          view(S.T,nsub+1:n,nsub+1:n))
        sep[2] = _trgsep(view(S.S,nsub+1:n,nsub+1:n),
                         view(S.S,1:nsub,1:nsub),
                         view(S.T,nsub+1:n,nsub+1:n),
                         view(S.T,1:nsub,1:nsub))
    elseif pnorm == 1
        scale = rtyone
        estu = norm1est!(fu,fuH,zeros(Ty,2*nsub*n2))
        sep[1] = scale / estu
        scale = rtyone
        estl = norm1est!(fl,flH,zeros(Ty,2*nsub*n2))
        sep[2] = scale / estl
    else
        throw(ArgumentError("pnorm must be :Frobenius or 1"))
    end
    return sep
end