text:
module TDMSReader

import BitOperations: bget
import DataStructures: OrderedDict

include("types.jl")

const _example_tdms=(@__DIR__) * "\\..\\test\\example_files\\reference_file.tdms"
const _example_incremental=[(@__DIR__) * "\\..\\test\\example_files\\incremental_test_$i.tdms" for i=1:6]
const _example_DAQmx=(@__DIR__) * "\\..\\test\\example_files\\DAQmx example.tdms"

export readtdms

readtdms() = readtdms(_example_tdms)
function readtdms(fn::AbstractString)
    s = open(fn)
    f=File()
    objdict=ObjDict()
    while !eof(s)
        (toc,nextsegmentoffset,rawdataoffset)=readleadin(s)
        if toc.kTocNewObjList
            empty!(objdict.current)
        end
        if toc.kTocMetaData
            readmetadata!(f, objdict, s)
        end
        if toc.kTocRawData
            readrawdata!(objdict, nextsegmentoffset-rawdataoffset, s)
        end
    end
    close(s)
    f
end

function readseginfo(fn::AbstractString)
    #STILL COULD USE SOME WORK HERE TO MAKE IT CLEAN AND NEAT
    s = open(fn)
    f=File()
    objdict=ObjDict()
    tdmsSeg = OrderedDict()
    lead_size = 28 #lead in is 28 bytes
    segCnt = 0
    while !eof(s)
        startPos = position(s)
        (toc,nextsegmentoffset,rawdataoffset)=readleadin(s)
        if toc.kTocNewObjList
            empty!(objdict.current)
        end

        nobj = read(s, UInt32)
        seek(s,startPos+lead_size)
        if toc.kTocMetaData
            readmetadata!(f, objdict, s)
        end
        if toc.kTocRawData
            #readrawdata!(objdict, nextsegmentoffset-rawdataoffset, s)
        end
        nextsegmentPos = Int64(startPos+nextsegmentoffset+lead_size)
        dataPos = Int64(startPos+lead_size+rawdataoffset)
        rawdatasize = nextsegmentoffset-rawdataoffset
        tdmsTmp = SegInfo(startPos,toc,nextsegmentPos,dataPos,nobj,rawdatasize)
        segCnt +=1
        tdmsSeg[string("seg",segCnt)] = tdmsTmp
        seek(s,nextsegmentPos)
    end
    close(s)
    return f,objdict,tdmsSeg
end

function readleadin(s::IO)
    @assert ntoh(read(s, UInt32)) == 0x54_44_53_6D
    toc=ToC(ltoh(read(s, UInt32)))
    toc.kTocBigEndian && throw(ErrorException("Big Endian files not supported"))
    read(s, UInt32) == 4713 || throw(ErrorException("File not recongnized as TDMS formatted file"))

    return (toc=toc,nextsegmentoffset=read(s, UInt64),rawdataoffset=read(s, UInt64),)
end

function readmetadata!(f::File, objdict::ObjDict, s::IO)
    n = read(s, UInt32)
    for i=1:n
        readobj!(f, objdict, s)
    end
end

function readobj!(f::File, objdict::ObjDict, s::IO)
    b=UInt8[]
    readbytes!(s, b, read(s, UInt32))
    objpath = String(b); empty!(b)
    # @info "Read @ position $(hexstring(position(s)))"
    rawdata=read(s, UInt32)
    hasrawdata=false
    hasnewchunk=false
    if rawdata==0xFF_FF_FF_FF #No Raw Data
    elseif rawdata==zero(UInt32) #Keep Chunk layout
        haskey(objdict.full, objpath) || throw(ErrorException("Previous Segment Missing"))
        hasrawdata=true
    elseif rawdata==0x00_00_12_69 || rawdata==0x00_00_13_69
        readDAQmx(s::IO, rawdata)
        throw(ErrorException("Not Implemented"))
    else
        hasrawdata=true
        hasnewchunk=true
        T=tdsTypes[read(s, UInt32)]
        ( T <: TDMSUnimplementedType ) && throw(ErrorException("TDMS Data Type of $T is not supported"))
        read(s,UInt32)==1 || throw(ErrorException("TDMS Array Dimension is not 1"))
        n=read(s,UInt64)
        if T == String
            throw(ErrorException("Need Functionality to Read String as raw data"))
        end
    end

    if objpath=="/"
        hasrawdata && throw(ErrorException("TDMS root should not have raw data"))
        props=f.props
        readprop!(props,s)
    else
        m=match(r"\/'(.+?)'(?:\/'(.+)')?", objpath)
        isnothing(m.captures) && throw(ErrorException("Object Path $objpath is malformed"))
        group,channel=m.captures[1:2]
        g=if haskey(f.groups, group)
            f[group]
        else
            get!(f.groups,group,Group())
        end
        if isnothing(channel) # Is a Group
            hasrawdata && throw(ErrorException("TDMS Group should not have raw data"))
            readprop!(g.props,s)
        else # Is  a Channel
            if haskey(g.channels,channel)
                if hasnewchunk && eltype(g[channel].data)==Nothing #need to replace existing channel if created without having data
                    setindex!(g.channels,TDMSReader.Channel{T}(g[channel].props),channel) #convert existing channel
                end
                chan = g[channel]
            else
                chan=if hasrawdata
                    get!(g.channels, channel, TDMSReader.Channel{T}())
                else
                    get!(g.channels, channel, TDMSReader.Channel{Nothing}())
                end
            end
            readprop!(chan.props, s)
        end
    end

    if hasrawdata
        if hasnewchunk
            objdict.full[objpath]=Chunk{T}(chan.data,n)
        end
        objdict.current[objpath]=objdict.full[objpath]
    end
end

function readprop!(props, s::IO)
    b=UInt8[]
    n=read(s,UInt32)
    for i=1:n
        readbytes!(s, b, read(s, UInt32))
        propname = String(b); empty!(b)
        T = tdsTypes[read(s,UInt32)]
        if T==String
            readbytes!(s, b, read(s, UInt32))
            props[propname] = String(b)
            empty!(b)
        else
            propval=read(s,T)
            props[propname]=propval
        end
    end
end

function readrawdata!(objects::NTuple{N,Chunk}, nbytes::Integer, s::IO) where N
    n = 0
    while n < nbytes && !eof(s)
        for x in objects
            T = eltype(x.data)
            for i=1:x.nsamples
                push!(x.data, read(s,T))
                n += sizeof(T)
            end
        end
    end
    return nothing
end

function readrawdata!(objdict::ObjDict, nbytes::Integer, s::IO)
    n = 0
    while n < nbytes && !eof(s)
        for (key,val) in objdict.current
            n += readchunk!(val.data, val.nsamples, s)
        end
    end
    return nothing
end

function readchunk!(v::Vector{T}, n::Integer, s::IO) where {T}
    for i = 1:n
        push!(v, read(s,T))
    end
    n*sizeof(T)
end

function seekalign(s::IO)
    x = position(s)
    mask = typeof(x)(0b11)
    seek(s, ifelse(x & mask > 0,x  & ~mask + 4, x))
end

function readDAQmx(s::IO, id)
    @info "Read DAQmc Raw Data @ $(hexstring(id))"
    T = tdsTypes[read(s, UInt32)]
    @info "Data type $T"
    read(s, UInt32)==1 || throw(ErrorException("TDMS Array Dimension is not 1"))
    chunksize = read(s, UInt64)
    @info "Chunk size = $chunksize"

    @info "Read vector of format change scalers"
    vectorsize = read(s, UInt32)
    V = tdsTypes[read(s, UInt32)]
    rawbufferindex = read(s, UInt32)
    rawbyteoffset = read(s, UInt32)
    sampleformatbitmap = read(s, UInt32)
    scaleid = read(s,UInt32)
    @info "Vector Size = $vectorsize"
    @info "DAQmx data type = $V"
    @info "Raw Buffer Index = $rawbufferindex"
    @info "Raw byte offset = $rawbyteoffset"
    @info "Sample Format Bitmatp = $sampleformatbitmap"
    @info "Scale ID = $scaleid"

    for i in 1:n

    end

end

function hexstring(x::Integer)
    "0x$(lpad(string(x,base=16),8,'0'))"
end

end # module

max_stars_repo_name: DanBerge/TDMSReader

max_issues_repo_name: DanBerge/TDMSReader


text:
module A

value = 1

end

ext: jl


text:
using Dates
using BusinessDays

export BusinessDayAdjustment, NoAdjustment, Following, Preceding
export ModifiedFollowing, ModifiedPreceding
export EOMNoAdjustment, EOMPreceding,EOMFollowing

abstract type BusinessDayAdjustment end
abstract type NoAdjustment <: BusinessDayAdjustment end
abstract type Following <: BusinessDayAdjustment end
abstract type ModifiedFollowing <: BusinessDayAdjustment end
abstract type Preceding <: BusinessDayAdjustment end
abstract type ModifiedPreceding <: BusinessDayAdjustment end
abstract type EOMNoAdjustment <: BusinessDayAdjustment end
abstract type EOMPreceding <: BusinessDayAdjustment end
abstract type EOMFollowing <: BusinessDayAdjustment end

BusinessDays.tobday(hc::HolidayCalendar, dt::Date, ::Type{NoAdjustment}) = dt
BusinessDays.tobday(hc::HolidayCalendar, dt::Date, ::Type{Following}) = tobday(hc, dt; forward=true)
BusinessDays.tobday(hc::HolidayCalendar, dt::Date, ::Type{Preceding}) = tobday(hc, dt; forward=false)
BusinessDays.tobday(hc::HolidayCalendar, df::Date, ::Type{EOMNoAdjustment}) = Dates.lastdayofmonth(dt)
BusinessDays.tobday(hc::HolidayCalendar, df::Date, ::Type{EOMPreceding}) = tobday(hc, Dates.lastdayofmonth(dt), Preceding)
BusinessDays.tobday(hc::HolidayCalendar, df::Date, ::Type{EOMFollowing}) = tobday(hc, Dates.lastdayofmonth(dt), Following)

BusinessDays.tobday(hc::HolidayCalendar, dt::Date, ::Type{ModifiedFollowing}) = begin
    day_adj = tobday(hc, dt, Following)
    if month(day_adj) ≠ month(dt)
        day_adj = tobday(hc, dt, Preceding)
    end
    day_adj
end
BusinessDays.tobday(hc::HolidayCalendar, dt::Date, ::Type{ModifiedPreceding}) = begin
    day_adj = tobday(hc, dt, Preceding)
    if month(day_adj) ≠ month(dt)
        day_adj = tobday(hc, dt, Following)
    end
    day_adj

end

max_stars_repo_name: binnisb/JuliaLearn

max_issues_repo_name: binnisb/JuliaLearn

max_forks_repo_name: binnisb/JuliaLearn


text:
# Linear discriminant
module LinDiscClassifier
	export LinDiscModel, LinDisc_train, LinDisc_exec
	using StatsFuns: softmax!

	# LinDiscModel struct
	struct LinDiscModel
		mean::Vector{Vector{Float64}}		# mean vector
		icm::Matrix{Float64}			# inverse covariance matrix covariance matrix
		priors::Dict{Int,Float64}		# class priors
		classes::Vector{Int}			# classes
		r1::Float64				# diagonal regularization weight
		r2::Float64				# trace regularization weight
	end

	# Printers
	Base.show(io::IO, m::LinDiscModel) = print(io, "Linear Discriminant model, $(length(m.mean)) classes, r1=$(m.r1), r2=$(m.r2)")

	# Train methods (classification)
	LinDisc_train(x::Matrix{Float64}, y::Vector{Int}, r1::Float64=0, r2::Float64=0) = begin

		m = size(x,1) 							# number of variables
		classes = sort(unique(y))					# classes
		priors = Dict(yi=>sum(yi.==y)/length(y) for yi in classes)	# priors

		mv = [vec(mean(x[:,y.==yi],2)) for yi in classes]		# mean vectors

		cv = zeros(m,m)
		@simd for yi in classes # iterate through the classes
			cv += priors[yi]*cov(x[:,y.==yi],2)
		end

		cv = (1.0-r1-r2)*cv + r1*Matrix(Diagonal(diag(cv))) + r2*trace(cv)*eye(m)

		return LinDiscModel(mv, inv(cv), priors, classes, r1, r2)
	end

	# Execution methods (classification)
	LinDisc_exec(m::LinDiscModel, x::Matrix{Float64})=begin
		C = length(m.classes)
		out = zeros(Float64, C, size(x,2))
		@simd for c in 1:C
		 	# 			Linear term		Constant term log(P(ωᵢ)) - ...
			@inbounds out[c,:] = (m.icm*m.mean[c])'*x .+ log(m.priors[m.classes[c]]).-1/2*m.mean[c]'*m.icm*m.mean[c]
		end

		# Apply softmax
		for j in 1:size(out,2)
			@inbounds softmax!(view(out,:,j))
		end

		return out
	end
end

##########################
# FunctionCell Interface #
##########################
"""
	lindisc(r1=0.0, r2=0.0)

Constructs an untrained cell that when piped data inside, returns a linear liscriminant classifier trained
function cell based on the input data and labels.

# Arguments
  * `r1`::Float64 represents the diagonal covariance matrix regularization parameter
  * `r2`::Float64 represents the covariance matrix trace regularization parameter

For more information:
	[1] L. Kuncheva "Combining Pattern Classifiers 2'nd Ed." 2014, ISBN 978-1-118-31523-1
"""
lindisc(r1::Float64=0.0, r2::Float64=0.0) = FunctionCell(lindisc, (r1,r2), ModelProperties(), "Linear discriminant classifier: r1=$r1, r2=$r2")

############################
# DataCell/Array Interface #
############################
"""
	lindisc(x, r1=0.0, r2=0.0)

Trains a linear discriminant classifcation model that using the data `x` and regularization parameters `r1`, `r2`.
"""
# Training
lindisc(x::T where T<:CellDataL, r1::Float64=0.0, r2::Float64=0.0) = lindisc((getx!(x), gety(x)), r1, r2)
lindisc(x::Tuple{T,S} where T<:AbstractVector where S<:AbstractVector, r1::Float64=0.0, r2::Float64=0.0) = lindisc((mat(x[1], LearnBase.ObsDim.Constant{2}()), x[2]), r1, r2)
lindisc(x::Tuple{T,S} where T<:AbstractMatrix where S<:AbstractVector, r1::Float64=0.0, r2::Float64=0.0) = begin

	@assert nobs(x[1]) == nobs(x[2]) "[lindisc] Expected $(nobs(x[1])) labels/values, got $(nobs(x[2]))."

	# Transform labels first
	enc = labelencn(x[2])
	yenc = label2ind.(x[2],enc)

	# Train model
	lddata = LinDiscClassifier.LinDisc_train(getobs(x[1]), yenc, r1, r2)

	# Build model properties
	modelprops = ModelProperties(nvars(x[1]), length(enc.label), enc)

	FunctionCell(lindisc, Model(lddata, modelprops), "Linear discriminant classifier: r1=$r1, r2=$r2")

end

# Execution
lindisc(x::T where T<:CellData, model::Model{<:LinDiscClassifier.LinDiscModel}) =
	datacell(lindisc(getx!(x), model), gety(x))
lindisc(x::T where T<:AbstractVector, model::Model{<:LinDiscClassifier.LinDiscModel}) =
	lindisc(mat(x, LearnBase.ObsDim.Constant{2}()), model)
lindisc(x::T where T<:AbstractMatrix, model::Model{<:LinDiscClassifier.LinDiscModel}) =
	LinDiscClassifier.LinDisc_exec(model.data, getobs(x))


text:
using GLM, QuantileRegression, Winston
using Requests

# Load data
url = "http://vincentarelbundock.github.io/Rdatasets/csv/quantreg/engel.csv"
df = readtable(Requests.get_streaming(url))
ResultQR = qreg(foodexp~income, df, .5)

# Fit quantile regression for a bunch of different quantiles
QNum = 35; # Number of equally spaced quantiles
DataPlot = reduce(vcat,[[coeftable(qreg(foodexp ~ income, df, i/(QNum+1); method = :ip)).cols[2][2] coeftable(qreg(foodexp ~ income, df, i/(QNum+1); method = :ip)).cols[3][2]]  for i in 1:QNum])

# Fit OLS model to compare
ResultLM = lm(foodexp~income, df) # Fit the model using OLS from the GLM-package
ols = vcat([[coeftable(ResultLM).cols[1][2] coeftable(ResultLM).cols[2][2]] for i =1:QNum]...)

PlotDF =  DataFrame(X = linspace(1,QNum,QNum)/(QNum+1),
                         Y = DataPlot[:,1],
                         ols = ols[:,1],
                         Ymin = DataPlot[:,1] - 1.96 * DataPlot[:,2],
                         Ymax = DataPlot[:,1] + 1.96 * DataPlot[:,2])

x = PlotDF[:X]
p = FramedPlot()
βτ   = Curve(x, PlotDF[:Y])
βols = Curve(x, PlotDF[:ols] , "color", "red")
ci_low  = Curve(x, PlotDF[:Ymin], "type" , "dash")
ci_high = Curve(x, PlotDF[:Ymax], "type" , "dash")
setattr(βτ,"label","β(τ)")
setattr(βols, "label", "β(ols)")
lgnd = Legend(.8,.2,[βτ, βols])
add(p, βτ, βols, ci_low, ci_high, lgnd)
setattr(p, "title" , "Quantile regression: Food Expenditure ~ Income")
setattr(p, "xlabel", "Quantiles")
setattr(p, "ylabel", "Coefficient on Income")
savefig(p, Pkg.dir("QuantileRegression")"/examples/qreg_example_plot.png")


text:
"""
```
function impulse_responses(m::AbstractDSGEVECMModel{S}, data::AbstractArray{S},
                           coint_mat::AbstractMatrix{S}, method::Symbol,
                           n_obs_shock::Int; horizon::Int = 0,
                           flip_shocks::Bool = false, verbose::Symbol = :none) where {S <: Real}

function impulse_responses(m::AbstractDSGEVECMModel{S}, coint_mat::AbstractMatrix{S}, method::Symbol,
                           n_obs_shock::Int; horizon::Int = 0, use_intercept::Bool = false,
                           flip_shocks::Bool = false, verbose::Symbol = :none) where {S <: Real}
```
computes the VECM impulse responses identified by the DSGE
```
sₜ = TTT × sₜ₋₁ + RRR × impact[:, i],
yₜ = ZZ × sₜ + DD + MM × impact[:, i],
```
where `impact[:, i]` is a linear combination of
(orthogonal) structural shocks `ϵₜ ∼ 𝒩 (0, I)`, and
`MM × impact[:, i]` are the correlated measurement errors.

We draw a β and Σᵤ from the posterior implied by the DSGE
and data, and we then compute normal VECM impulse responses given those
coefficients, innovations variance-covariance matrix, and
the matrix specifying cointegrating relationships in observables. The weight placed
on the DSGE is encoded by the field `λ` of the DSGEVECM object `m`.

Given β, Σᵤ, we compute impulse responses using one of the
available identifiction strategies to the VECM system
```
Δŷₜ₊₁ = eₜβₑ + X̂ₜ₊₁βᵥ + uₜ₊₁,
```
where `βₑ` are the coefficients for the error correction terms;
`eₜ₊₁` are the error correction terms specifying the cointegrating relationships;
`βᵥ` are the coefficients for the VAR terms;
`X̂ₜ₊₁` are the lags of observables in period `t + 1`, i.e. `yₜ, yₜ₋₁, ..., yₜ₋ₚ`,
and `uₜ₊₁ ∼ 𝒩 (0, Σ)`.

If the second function is used (where `data` is not an input), then we assume
the user wants to compute the VECM approximation of the DSGE,
regardless of the `λ` value in `m`. Note that this function will not
update the value of `λ` in `m` (even though we are computing the DSGE-VECM(∞) approximation).

### Inputs
* `coint_mat::AbstractMatrix{S}`: matrix specifying the cointegrating relationships
    in observables. Given a matrix `data` with dimensions `n_observables × T`,
    multiplying `coint_mat * data` should yield a `n_coint × T` matrix, where
    `n_coint` is the number of cointegrating relationships and `T` is
    the number of periods of data.
* `method::Symbol`: The available methods are `:cholesky`, `:maxBC`, and `:choleskyLR`.
    See the docstrings `impulse_responses` for VECMs specifically.
* `n_obs_shock::Int`: The index of the observable corresponding to the orthogonalized shock
    causing the impulse response.

### Keyword Arguments
* `horizon::Int`: the desired horizon of the impulse responses.
* `use_intercept::Bool`: use an intercept term for the VECM approximation
* `flip_shocks::Bool`: default is a "negative" impulse response on impact.
    Set to `true` for the positive impulse response.
"""
function impulse_responses(m::AbstractDSGEVECMModel{S}, data::AbstractArray{S},
                           coint_mat::AbstractMatrix{S}, method::Symbol,
                           n_obs_shock::Int; horizon::Int = impulse_response_horizons(m),
                           flip_shocks::Bool = false, verbose::Symbol = :none) where {S <: Real}
    β, Σ = compute_system(m, data; verbose = verbose)
    Σ += Σ'
    Σ ./= 2

    return impulse_responses(β, Σ, coint_mat, n_obs_shock, horizon;
                             method = method, flip_shocks = flip_shocks,
                             use_intercept = true)
end

function impulse_responses(m::AbstractDSGEVECMModel{S}, coint_mat::AbstractMatrix{S}, method::Symbol,
                           n_obs_shock::Int; horizon::Int = 0, use_intercept::Bool = false,
                           flip_shocks::Bool = false, verbose::Symbol = :none) where {S <: Real}
    β, Σ = compute_system(m; verbose = verbose, use_intercept = use_intercept)
    Σ += Σ'
    Σ ./= 2

    return impulse_responses(β, Σ, coint_mat, n_obs_shock,
                             horizon > 0 ? horizon : impulse_response_horizons(m);
                             method = method, use_intercept = use_intercept,
                             flip_shocks = flip_shocks)
end

"""
```
function impulse_responses(m::AbstractDSGEVECMModel{S}, data::AbstractArray{S},
    X̂::Matrix{S} = Matrix{S}(undef, 0, 0);
    horizon::Int = 0, MM::Matrix{S} = Matrix{S}(undef, 0, 0),
    flip_shocks::Bool = false, draw_shocks::Bool = false,
    verbose::Symbol = :none) where {S <: Real}
```
computes the VECM impulse responses identified by the DSGE
```
sₜ = TTT × sₜ₋₁ + RRR × impact[:, i],
yₜ = ZZ × sₜ + DD + MM × impact[:, i],
```
where `impact[:, i]` is a linear combination of
(orthogonal) structural shocks `ϵₜ ∼ 𝒩 (0, I)`, and
`MM × impact[:, i]` are the correlated measurement errors.

The VECM impulse responses are computed according to
```
Δŷₜ₊₁ = eₜβₑ + X̂ₜ₊₁βᵥ + uₜ₊₁,
```
where `βₑ` are the coefficients for the error correction terms;
`eₜ` are the error correction terms specifying the cointegrating relationships;
`βᵥ` are the coefficients for the VAR terms;
`X̂ₜ₊₁` are the lags of observables in period `t + 1`, i.e. `yₜ, yₜ₋₁, ..., yₜ₋ₚ`,
and `uₜ₊₁ ∼ 𝒩 (0, Σ)`. Note these impulses responses are *not*
computed in deviations from the baseline forecast `Δŷₜ₊₁ = eₜ₊₁βₑ + X̂ₜ₊₁βᵥ`.
To compute these impulse responses, use the keyword `deviations`.

The shock `uₜ₊₁` is identified by assuming
```
Σᵤ = 𝔼[u × u'] = chol(Σᵤ) × Ω × ϵₜ,
```
where the rotation matrix `Ω` is the `Q` matrix from a QR decomposition
of the impact response matrix corresponding to the state space system, i.e.
```
Ω, _ = qr(∂yₜ / ∂ϵₜ').
```
The impact response matrix is constructed using only the stationary component of the
state space system and ignores the cointegration components of `ZZ` and `DD`.

For reference, see Del Negro and Schorfheide (2004),
Del Negro and Schorfheide (2006), and Del Negro and Schorfheide (2009).

### Inputs
* `coint_mat::Matrix{S}`: matrix specifying the cointegrating relationships
    in the actual `data` matrix. Evaluating `coint_mat * data` should yield
    a time series of the cointegrating relationships.
* `X̂::Matrix{S}`: covariates for the first "forecast" period
    of the impulse response, i.e. if we have a VECM with `p` lags, then
```
X̂ = [eₜ, 1, ŷₜ, ŷₜ₋₁, ..., ŷₜ₋ₚ₊₁]
```
so that, when β is the vector of VECM coefficients, then
```
𝔼[ŷₜ₊₁] = kron(I, X̂') * β.
```
Internally, we do equivalent matrix operations to avoid allocating
the Kronecker product.

****
NOTE: this function generally involves taking random draws from
probability distributions, so seeds need to be set
to achieve reproducibility.
****

### Keywords
* `horizon::Int`: horizon of impulse responses
* `flip_shocks::Bool`: impulse response shocks are negative by default. Set to `true` for
    a positive signed shock.
* `draw_shocks::Bool`: true if you want to draw shocks along the entire horizon
* `deviations::Bool`: set true to compute the impulse response in deviations
    rather than as a forecast. Mechnically, we ignore `X̂` (treated as zeros)
    and the intercept term.
* `verbose::Symbol`: quantity of output desired
"""
function impulse_responses(m::AbstractDSGEVECMModel{S}, data::AbstractArray{S},
                           coint_mat::Matrix{S}, X̂::Vector{S} = Vector{S}(undef, 0);
                           horizon::Int = 0, MM::Matrix{S} = Matrix{S}(undef, 0, 0),
                           flip_shocks::Bool = false, draw_shocks::Bool = false,
                           deviations::Bool = false, verbose::Symbol = :none) where {S <: Real}

    # Prepare X̂
    n_obs = size(data, 1)
    n_coint = size(coint_mat, 1)
    k = n_coint + n_obs * get_lags(m) + 1
    if isempty(X̂)
        XX = lag_data(data, get_lags(m); use_intercept = true)
        addcoint = coint_mat * data[:, end]
        X̂ = vcat(addcoint, 1, data[:, end], XX[end, 1+1:k - n_obs])
    end
    h = (horizon > 0) ? horizon : impulse_response_horizons(m) # grab horizons

    # Compute underlying state space system
    system = compute_system(m; get_system = true, use_intercept = true)

    # Compute VECM coefficients
    β, Σ = compute_system(m, data; verbose = verbose)
    Σ += Σ'
    Σ ./= 2.

    # Use rotation identification
    if isempty(MM)
        if hasmethod(measurement_error, (typeof(m),))
            _, MM = measurement_error(m)
        else
            MM = zeros(S, n_observables(m), n_shocks(m))
        end
    end
    return impulse_responses(system[:TTT], system[:RRR], system[:ZZ], system[:DD], MM,
                             system[:QQ], k, n_obs, n_coint, β, Σ,
                             coint_mat, X̂, h; flip_shocks = flip_shocks,
                             draw_shocks = draw_shocks, deviations = deviations)
end

"""
```
function impulse_responses(TTT::Matrix{S}, RRR::Matrix{S}, ZZ::Matrix{S},
                           DD::Vector{S}, MM::Matrix{S}, QQ::Matrix{S},
                           k::Int, n_obs::Int, n_coint::Int, β::Matrix{S}, Σ::Matrix{S},
                           coint_mat::Matrix{S}, horizon::Int, X̂::Matrix{S} = zeros(S, k);
                           flip_shocks::Bool = false, draw_shocks::Bool = false,
                           deviations::Bool = false,
                           test_shocks::Matrix{S} =
                           Matrix{S}(undef, 0, 0)) where {S<:Real}
```
computes the VECM impulse responses identified by the state space system
```
sₜ = TTT × sₜ₋₁ + RRR × ϵₜ
yₜ = ZZ × sₜ + DD + MM × ϵₜ
```
where `ϵₜ ∼ 𝒩 (0, QQ)` and `MM × ϵₜ` are the correlated measurement errors.

Consider the VECM
```
Δŷₜ₊₁ = eₜβₑ + X̂ₜ₊₁βᵥ + uₜ₊₁,
```
where `βₑ` are the coefficients for the error correction terms;
`eₜ` are the error correction terms specifying the cointegrating relationships;
`βᵥ` are the coefficients for the VAR terms (including the intecept)o;
`X̂ₜ₊₁` are the lags of observables in period `t + 1`, i.e. `yₜ, yₜ₋₁, ..., yₜ₋ₚ₊₁`;
and `uₜ₊₁ ∼ 𝒩 (0, Σ)`. Note these impulses responses are *not*
computed in deviations from the baseline forecast `Δŷₜ₊₁ = eₜ₊₁βₑ + X̂ₜ₊₁βᵥ`.
To compute these impulse responses, set the keyword `deviations = true`.

The shock `uₜ₊₁` is identified via
```
Σᵤ = 𝔼[u × u'] = chol(Σᵤ) × Ω × ϵₜ,
```
where the rotation matrix `Ω` is the `Q` matrix from a QR decomposition
of the impact response matrix corresponding to the state space system, i.e.
```
Ω, _ = qr(∂yₜ / ∂ϵₜ').
```
The impact response matrix is constructed using only the stationary component of the
state space system and ignores the cointegration components of `ZZ` and `DD`.

The data are assumed to have dimensions `n_obs × T`, and
the cointegration relationships in the data are given by `coint_mat * data`, where
`coint_mat` has dimensions `n_coint × n_obs`. The variable `k` is the
number of total regressors in the VECM, including cointegration terms.

For reference, see Del Negro and Schorfheide (2004), Del Negro and Schorfheide (2006),
and Del Negro, Schorfheide, Smets, and Wouters (2007).
"""
function impulse_responses(TTT::Matrix{S}, RRR::Matrix{S}, ZZ::Matrix{S},
                           DD::Vector{S}, MM::Matrix{S}, QQ::Matrix{S},
                           k::Int, n_obs::Int, n_coint::Int, β::Matrix{S}, Σ::Matrix{S},
                           coint_mat::Matrix{S}, horizon::Int, X̂::Vector{S} = zeros(S, k);
                           flip_shocks::Bool = false, draw_shocks::Bool = false,
                           deviations::Bool = false,
                           test_shocks::Matrix{S} =
                           Matrix{S}(undef, 0, 0)) where {S <: Real}

    # Grab stationary components
    ZZ = ZZ[1:n_obs, :]
    DD = DD[1:n_obs]
    MM = MM[1:n_obs, :]

    # Compute impulse responses of predicted values for each β, Σ, and rotation
    a0_m = convert(Matrix{S},
                   dropdims(impulse_responses(TTT, RRR, ZZ, DD, MM,
                                              sqrt.(QQ), 1, # 1 b/c just want impact and sqrt -> 1 sd shock
                                              accumulate = false); dims = 2)')
    rotation, r_a = qr(a0_m)
    rotation = sign.(diag(r_a)) .* rotation' # normalizes each row (change signs) so that lower diagonal (r_a') has positive diagonal elements
    Σ_chol = cholesky(Σ).L * rotation # mapping from structural shocks to innovations in VECM

    # Remove intercept term if in deviations
    if deviations
        β = β[vcat(1:n_coint, n_coint + 2:size(β, 1)), :]
        X̂ = zeros(S, size(β, 1))
        k -= 1
    end

    if draw_shocks || !isempty(test_shocks)
        ŷ = Matrix{S}(undef, n_obs, horizon)

        if isempty(test_shocks)
            shocks = randn(size(RRR, 2), horizon) # shocks getting drawn are structural shocks
            if flip_shocks
                @warn "Shocks are being drawn, so flipping shocks does not make sense in this context. Ignoring `flip_shocks` . . ."
            end
        else
            shocks = test_shocks
        end
        for t = 1:horizon
            out      = vec(X̂' * β) + Σ_chol * shocks[:, t] # X̂ normally would be [X̂ 0 0; 0 X̂ 0; 0 0 X̂] if n_obs = 3,
            ŷ[:, t]  = out                                 # but this way of coding it results in less memory storage
            addcoint = X̂[1:n_coint] + coint_mat * out      # Predicted cointegration terms

            X̂ = deviations ? vcat(addcoint,  out, X̂[n_coint + 1:k - n_obs]) :
                vcat(addcoint, 1.,  out, X̂[n_coint + 1 + 1:k - n_obs]) # XXl = X̂[n_coint + 1 + 1:k - n_obs]
        end
    else
        nshocks = size(RRR, 2)
        ŷ       = Array{S, 3}(undef, n_obs, horizon, nshocks)
        old_X̂   = X̂
        shocks  = zeros(S, nshocks)

        for i = 1:nshocks
            X̂ = copy(old_X̂)
            shocks[i] = flip_shocks ? sqrt(QQ[i, i]) :
                -sqrt(QQ[i, i]) # a negative 1 s.d. shock by default
            out        = vec(X̂' * β) + Σ_chol * shocks # do impact separately
            shocks[i]  = 0. # set back to zero
            ŷ[:, 1, i] = out
            addcoint   = X̂[1:n_coint] + coint_mat * out
            X̂          = deviations ? vcat(addcoint,  out, X̂[n_coint + 1:k - n_obs]) :
                vcat(addcoint, 1.,  out, X̂[n_coint + 1 + 1:k - n_obs]) # XXl = X̂[n_coint + 1 + 1:k - n_obs]
            for t = 2:horizon
                out        = vec(X̂' * β)
                ŷ[:, t, i] = out
                addcoint   = X̂[1:n_coint] + coint_mat * out
                X̂          = deviations ? vcat(addcoint,  out, X̂[n_coint + 1:k - n_obs]) :
                    vcat(addcoint, 1.,  out, X̂[n_coint + 1 + 1:k - n_obs]) # XXl = X̂[n_coint + 1 + 1:k - n_obs]
            end
        end
    end

    return ŷ
end


text:
"""
PkgFiles encodes information about the current location of a package.
Fields:
- `id`: the `PkgId` of the package
- `basedir`: the current base directory of the package
- `files`: a list of files (relative path to `basedir`) that define the package.

Note that `basedir` may be subsequently updated by Pkg operations such as `add` and `dev`.
"""
mutable struct PkgFiles
    id::PkgId
    basedir::String
    files::Vector{String}
end

PkgFiles(id::PkgId, path::AbstractString) = PkgFiles(id, path, String[])
PkgFiles(id::PkgId, ::Nothing) = PkgFiles(id, "")
PkgFiles(id::PkgId) = PkgFiles(id, normpath(basepath(id)))
PkgFiles(id::PkgId, files::AbstractVector{<:AbstractString}) =
    PkgFiles(id, normpath(basepath(id)), files)

# Abstraction interface
Base.PkgId(info::PkgFiles) = info.id
srcfiles(info::PkgFiles) = info.files
basedir(info::PkgFiles) = info.basedir

function Base.show(io::IO, info::PkgFiles)
    compact = get(io, :compact, false)
    if compact
        print(io, "PkgFiles(", info.id.name, ", ", info.basedir, ", ")
        show(io, info.files)
        print(io, ')')
    else
        println(io, "PkgFiles(", info.id, "):")
        println(io, "  basedir: \"", info.basedir, '"')
        print(io, "  files: ")
        show(io, info.files)
    end
end


text:
"""
Simulations to estimate coverage of confidence bands with Blasques methods and Buccheri's enhancement
"""

#region import and models

begin
using DrWatson

using ScoreDrivenERGM
import ScoreDrivenERGM:StaticNets, DynNets
import ScoreDrivenERGM.DynNets:SdErgm,SdErgmDirBin0Rec0, sample_ergm_sequence, stats_from_mat, array_2_vec_all_par, unrestrict_all_par, conf_bands_par_uncertainty, avg_grad_and_hess_obj_SD_filter_time_seq, conf_bands_par_uncertainty, number_ergm_par, estimate_filter_and_conf_bands, conf_bands_coverage, estimate, mle_distrib_filtered_par, plot_filtered_and_conf_bands
using ScoreDrivenERGM.Utilities

using PyPlot
pygui(true)

using ForwardDiff
using StatsBase
using LinearAlgebra
using Distributions
using Statistics

using JLD2

model_mle = DynNets.SdErgmDirBin0Rec0_mle()
model_pmle = DynNets.SdErgmDirBin0Rec0_pmle()
indTvPar = trues(2)

end
#endregion

# #region coverage simulations
begin
model = model_pmle
N=100
T=100
dgpType =  "AR"

nVals = [ 50, 100, 300]
tVals = [100, 300]
models = [model_pmle]
nSampleCoverage=50

dgpOptions = (minAlpha = 0.2, maxAlpha = 0.3, nCycles=1.5, phaseAlpha = 0.1π, phaseshift = 0.1, plotFlag=false, B =0.98, sigma = 0.0005, A = 0.3)
quantilesVals = [[0.975, 0.025]]

parDgpT = DynNets.sample_time_var_par_from_dgp(model_mle, dgpType, N, T;  dgpOptions...)

# @elapsed allCoverBuccheri, allCoverBlasques, allvEstSdResPar, allfVecT_filt, allConfBandsBuccheri,allConfBandsBlasques, allErrFlags = conf_bands_coverage(model, dgpType, dgpOptions, T, N, 2, quantilesVals; plotFlag=false)

# nNVals = length(nVals)
# nTVals = length(tVals)
# nModels = length(models)

# allCoverBuccheriVarN = Array{typeof(allCoverBuccheri),3}(undef, nNVals, nTVals, nModels)
# allCoverBlasquesVarN = Array{typeof(allCoverBlasques),3}(undef, nNVals, nTVals, nModels)
# allvEstSdResParVarN = Array{typeof(allvEstSdResPar),3}(undef, nNVals, nTVals, nModels)
# allfVecT_filtVarN = Array{typeof(allfVecT_filt),3}(undef, nNVals, nTVals, nModels)
# allConfBandsBuccheriVarN = Array{typeof(allConfBandsBuccheri),3}(undef, nNVals, nTVals, nModels)
# allConfBandsBlasquesVarN = Array{typeof(allConfBandsBlasques),3}(undef, nNVals, nTVals, nModels)
# fractErrVarN = Array{typeof(allErrFlags),3}(undef, nNVals, nTVals, nModels)

# parDgpTvarN = Array{Array{Float64,2},2}(undef, nNVals, nTVals)

# for (indT, T) in Iterators.enumerate(tVals)
#     for (indN, N) in Iterators.enumerate(nVals)

#         parDgpTvarN[indN, indT] = parDgpT
#         for (indM, model) in Iterators.enumerate(models)

#             @elapsed allCoverBuccheri, allCoverBlasques, allvEstSdResPar, allfVecT_filt, allConfBandsBuccheri, allConfBandsBlasques, allErrFlags = conf_bands_coverage(model, dgpType, dgpOptions, T, N,  nSampleCoverage, quantilesVals; plotFlag = false)

#             allCoverBuccheriVarN[indN, indT, indM] = allCoverBuccheri
#             allCoverBlasquesVarN[indN, indT, indM] = allCoverBlasques
#             allvEstSdResParVarN[indN, indT, indM] = allvEstSdResPar
#             allfVecT_filtVarN[indN, indT, indM] = allfVecT_filt
#             allConfBandsBuccheriVarN[indN, indT, indM] = allConfBandsBuccheri
#             allConfBandsBlasquesVarN[indN, indT, indM] = allConfBandsBlasques
#             fractErrVarN[indN, indT, indM] = allErrFlags

#         end
#     end
# end
# nErgmPar=2

# @save("./data/confBands_$(dgpType)_B_$(dgpOptions.B)_A_$(dgpOptions.A)_sig_$(dgpOptions.sigma)_$(nVals)_$(tVals)_nSample_$(nSampleCoverage)_pmle.jld", allCoverBuccheriVarN, allCoverBlasquesVarN, allfVecT_filtVarN, allConfBandsBuccheriVarN, allConfBandsBlasquesVarN, fractErrVarN, parDgpTvarN)
end

begin
using JLD2
nSampleCoverage=50
dgpType = "SD"
models = [model_pmle]

nNVals = length(nVals)
nTVals = length(tVals)
nModels = length(models)

@load("$(datadir())\\old_pre_drWatson_git_ignored\\confBands_$(dgpType)_B_$(dgpOptions.B)_A_$(dgpOptions.A)_sig_$(dgpOptions.sigma)_$(nVals)_$(tVals)_nSample_$(nSampleCoverage)_pmle.jld", allCoverBuccheriVarN, allCoverBlasquesVarN, allfVecT_filtVarN, allConfBandsBuccheriVarN, allConfBandsBlasquesVarN, fractErrVarN, parDgpTvarN)

avgCover =zeros(2,nNVals, nTVals,nModels, 2, nSampleCoverage)
for (indT, T) in Iterators.enumerate(tVals)
    for (indN, N) in Iterators.enumerate(nVals)
        for (indM, model) in Iterators.enumerate(models)
            for indPar in 1:2
                for n=1:nSampleCoverage
                    avgCover[indPar, indN, indT, indM, 1, n] = mean(allCoverBuccheriVarN[indN, indT, indM][indPar,:,1,n])
                    avgCover[indPar, indN, indT, indM, 2, n] = mean(allCoverBlasquesVarN[indN, indT, indM][indPar,:,1,n])
                end
            end
        end
    end
end

indM =1
indB =2

nominalLevel = 0.95
parNames = ["θ", "η"]
BandNames = ["Parameters + Filtering Uncertainty", "Parameters Uncertainty"]

fig, ax1 = plt.subplots(2, length(tVals),figsize=(12, 6), sharey =true)
fig.canvas.set_window_title("Confidence Bands' Coverages $(BandNames[indB])")
fig.subplots_adjust(left=0.075, right=0.95, top=0.9, bottom=0.25)
fig.suptitle("Confidence Bands' Coverages $(BandNames[indB]) DGP = $(dgpType)")

for (indT, T) in Iterators.enumerate(tVals)
    for indPar in 1:2

        data = avgCover[indPar,:,indT,indM,indB,:]'
        bp = ax1[indPar, indT].boxplot(data, notch=0, sym="+", vert=1, whis=1.5)

        ax1[indPar, indT].yaxis.grid(true, linestyle="-", which="major", color="lightgrey", alpha=0.5)

        # Hide these grid behind plot objects
        xlims = ax1[indPar, indT].get_xlim()
        ax1[indPar, indT].hlines(nominalLevel, xlims[1], xlims[2], linestyle=":" , colors = "r")
        ax1[indPar, indT].set_axisbelow(true)
        ax1[indPar, indT].set_title("T = $T")
        ax1[indPar, indT].set_xlabel("Network Size")
        ax1[indPar, indT].set_ylabel("Coverages for $(parNames[indPar])")
        ax1[indPar, indT].set_xticklabels(nVals, rotation=45, fontsize=8)
    end
end
tight_layout()

end

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
# store perf info of each variable
mutable struct SplitInfo_gpu{T<:AbstractFloat, S<:Int}
    gain::T
    ∑δL::Vector{T}
    ∑δ²L::Vector{T}
    ∑𝑤L::T
    ∑δR::Vector{T}
    ∑δ²R::Vector{T}
    ∑𝑤R::T
    gainL::T
    gainR::T
    𝑖::S
    feat::S
    cond::T
end

struct TreeNode_gpu{T<:AbstractFloat, S<:Int, B<:Bool}
    left::S
    right::S
    feat::S
    cond::T
    gain::T
    pred::Vector{T}
    split::B
end

TreeNode_gpu(left::S, right::S, feat::S, cond::T, gain::T, K) where {T<:AbstractFloat, S<:Int} = TreeNode_gpu(left, right, feat, cond, gain, zeros(T,K), true)
TreeNode_gpu(pred::Vector{T}) where {T} = TreeNode_gpu(0, 0, 0, zero(T), zero(T), pred, false)

# single tree is made of a root node that containes nested nodes and leafs
struct TrainNode_gpu{T<:AbstractFloat, S<:Int}
    parent::S
    depth::S
    ∑δ::Vector{T}
    ∑δ²::Vector{T}
    ∑𝑤::T
    gain::T
    𝑖::Vector{S}
    𝑗::Vector{S}
end

# single tree is made of a root node that containes nested nodes and leafs
struct Tree_gpu{T<:AbstractFloat, S<:Int}
    nodes::Vector{TreeNode_gpu{T,S,Bool}}
end

# gradient-boosted tree is formed by a vector of trees
struct GBTree_gpu{T<:AbstractFloat, S<:Int}
    trees::Vector{Tree_gpu{T,S}}
    params::EvoTypes
    metric::Metric
    K::Int
    levels
end


text:
# Utility functions for field tracing.

export select_seeds

"""
	 select_seeds(x, y; nsegment=(5, 5))

Generate uniform seeding points in the grid range `x` and `y` in `nsegment`. If `nsegment`
specified, use the keyword input, otherwise it will be overloaded by the 3D version seed
generation function!
"""
function select_seeds(x, y; nsegment=(5, 5))
   xmin, xmax = extrema(x)
   ymin, ymax = extrema(y)

   dx = (xmax - xmin) / nsegment[1]
   dy = (ymax - ymin) / nsegment[2]

   xrange = xmin .+ dx*((1:nsegment[1]) .- 0.5)
   yrange = ymin .+ dy*((1:nsegment[2]) .- 0.5)

   seeds = zeros(eltype(x[1]), 2, prod(nsegment))
   for i in 0:prod(nsegment)-1
      seeds[1,i+1] = xrange[i % nsegment[1] + 1]
      seeds[2,i+1] = yrange[i ÷ nsegment[2] + 1]
   end
   seeds
end

"""
	 select_seeds(x, y, z; nsegment=(5, 5, 5))

Generate uniform seeding points in the grid range `x`, `y` and `z` in `nsegment`.
"""
function select_seeds(x, y, z; nsegment=(5, 5, 5))
   xmin, xmax = extrema(x)
   ymin, ymax = extrema(y)
   zmin, zmax = extrema(z)

   dx = (xmax - xmin) / nsegment[1]
   dy = (ymax - ymin) / nsegment[2]
   dz = (zmax - zmin) / nsegment[3]

   xrange = xmin .+ dx*((1:nsegment[1]) .- 0.5)
   yrange = ymin .+ dy*((1:nsegment[2]) .- 0.5)
   zrange = zmin .+ dz*((1:nsegment[3]) .- 0.5)

   seeds = zeros(eltype(x[1]), 3, prod(nsegment))
   for k in 1:nsegment[3], j in 1:nsegment[2], i in 1:nsegment[1]
      seeds[1,i+(j-1)*nsegment[2]+(k-1)*nsegment[2]*nsegment[3]] = xrange[i]
      seeds[2,i+(j-1)*nsegment[2]+(k-1)*nsegment[2]*nsegment[3]] = yrange[j]
      seeds[3,i+(j-1)*nsegment[2]+(k-1)*nsegment[2]*nsegment[3]] = zrange[k]
   end
   seeds
end


text:
# Invertible multiscale HINT network from Kruse et. al (2020)
# Author: Gabrio Rizzuti, grizzuti3@gatech.edu
# Date: October 2020

export NetworkMultiScaleHINT, NetworkMultiScaleHINT3D

"""
    H = NetworkMultiScaleHINT(n_in, n_hidden, L, K; split_scales=false, k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, ndims=2)

    H = NetworkMultiScaleHINT3D(n_in, n_hidden, L, K; split_scales=false, k1=3, k2=3, p1=1, p2=1, s1=1, s2=1)

 Create a multiscale HINT network for data-driven generative modeling based
 on the change of variables formula.

 *Input*:

 - 'n_in': number of input channels

 - `n_hidden`: number of hidden units in residual blocks

 - `L`: number of scales (outer loop)

 - `K`: number of flow steps per scale (inner loop)

 - `split_scales`: if true, split output in half along channel dimension after each scale. Feed one half through the next layers,
    while saving the remaining channels for the output.

 - `k1`, `k2`: kernel size for first and third residual layer (`k1`) and second layer (`k2`)

 - `p1`, `p2`: respective padding sizes for residual block layers

 - `s1`, `s2`: respective strides for residual block layers

 - `ndims` : number of dimensions

 *Output*:

 - `H`: multiscale HINT network

 *Usage:*

 - Forward mode: `Z, logdet = H.forward(X)`

 - Inverse mode: `X = H.inverse(Z)`

 - Backward mode: `ΔX, X = H.backward(ΔZ, Z)`

 *Trainable parameters:*

 - None in `H` itself

 - Trainable parameters in activation normalizations `H.AN[i]`,
 and in coupling layers `H.CL[i]`, where `i` ranges from `1` to `depth`.

 See also: [`ActNorm`](@ref), [`CouplingLayerHINT!`](@ref), [`get_params`](@ref), [`clear_grad!`](@ref)
"""
struct NetworkMultiScaleHINT <: InvertibleNetwork
    AN::AbstractArray{ActNorm, 2}
    CL::AbstractArray{CouplingLayerHINT, 2}
    X_dims::Union{Array{Tuple, 1}, Nothing}
    L::Int64
    K::Int64
    split_scales::Bool
end

@Flux.functor NetworkMultiScaleHINT

# Constructor
function NetworkMultiScaleHINT(n_in::Int64, n_hidden::Int64, L::Int64, K::Int64;
                               split_scales=false, k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, ndims=2)

    AN = Array{ActNorm}(undef, L, K)
    CL = Array{CouplingLayerHINT}(undef, L, K)
    if split_scales
        X_dims = Array{Tuple}(undef, L-1)
        channel_factor = 2
    else
        X_dims = nothing
        channel_factor = 4
    end

    # Create layers
    for i=1:L
        for j=1:K
            AN[i, j] = ActNorm(n_in*4; logdet=true)
            CL[i, j] = CouplingLayerHINT(n_in*4, n_hidden; permute="full", k1=k1, k2=k2, p1=p1, p2=p2,
                                         s1=s1, s2=s2, logdet=true, ndims=ndims)
        end
        n_in *= channel_factor
    end

    return NetworkMultiScaleHINT(AN, CL, X_dims, L, K, split_scales)
end

NetworkMultiScaleHINT3D(args...; kw...) = NetworkMultiScaleHINT(args...; kw..., ndims=3)

# Concatenate states Zi and final output
function cat_states(X_save::AbstractArray{Array, 1}, X::AbstractArray{Float32, 4})
    X_full = []
    for j=1:size(X_save, 1)
        X_full = cat(X_full, vec(X_save[j]); dims=1)
    end
    X_full = cat(X_full, vec(X); dims=1)
    return Float32.(X_full)  # convert to Array{Float32, 1}
end

# Split 1D vector in latent space back to states Zi
function split_states(X_dims::AbstractArray{Tuple, 1}, X_full::AbstractArray{Float32, 1})
    L = length(X_dims) + 1
    X_save = Array{Array}(undef, L-1)
    count = 1
    for j=1:L-1
        X_save[j] = reshape(X_full[count: count + prod(XY_dims[j])-1], X_dims[j])
        count += prod(X_dims[j])
    end
    X = reshape(X_full[count: count + prod(X_dims[end])-1], Int.(X_dims[end].*(.5, .5, 4, 1)))
    return X_save, X
end

# Forward pass and compute logdet
function forward(X, H::NetworkMultiScaleHINT)
    H.split_scales && (X_save = Array{Array}(undef, H.L-1))
    logdet = 0f0
    for i=1:H.L
        X = wavelet_squeeze(X)
        for j=1:H.K
            X_, logdet1 = H.AN[i, j].forward(X)
            X, logdet2 = H.CL[i, j].forward(X_)
            logdet += (logdet1 + logdet2)
        end
        if H.split_scales && i < H.L    # don't split after last iteration
            X, Z = tensor_split(X)
            X_save[i] = Z
            H.X_dims[i] = size(Z)
        end
    end
    H.split_scales && (X = cat_states(X_save, X))
    return X, logdet
end

# Inverse pass and compute gradients
function inverse(Z, H::NetworkMultiScaleHINT)
    H.split_scales && ((X_save, Z) = split_states(H.X_dims, Z))
    for i=H.L:-1:1
        if H.split_scales && i < H.L
            Z = tensor_cat(Z, X_save[i])
        end
        for j=H.K:-1:1
            Z_ = H.CL[i, j].inverse(Z)
            Z = H.AN[i, j].inverse(Z_; logdet=false)
        end
        Z = wavelet_unsqueeze(Z)
    end
    return Z
end

# Backward pass and compute gradients
function backward(ΔZ, Z, H::NetworkMultiScaleHINT; set_grad::Bool=true)

    # Split data and gradients
    if H.split_scales
        ΔX_save, ΔZ = split_states(H.X_dims, ΔZ)
        X_save, Z = split_states(H.X_dims, Z)
    end

    if ~set_grad
        Δθ = Array{Parameter, 1}(undef, 0)
        ∇logdet = Array{Parameter, 1}(undef, 0)
    end

    for i=H.L:-1:1
        if H.split_scales && i < H.L
            ΔZ = tensor_cat(ΔZ, ΔX_save[i])
            Z = tensor_cat(Z, X_save[i])
        end
        for j=H.K:-1:1
            if set_grad
                ΔZ_, Z_ = H.CL[i, j].backward(ΔZ, Z)
                ΔZ, Z = H.AN[i, j].backward(ΔZ_, Z_)
            else
                ΔZ_, Δθcl, Z_, ∇logdet_cl = H.CL[i, j].backward(ΔZ, Z; set_grad=set_grad)
                ΔZ, Δθx, Z, ∇logdet_x = H.AN[i, j].backward(ΔZ_, Z_; set_grad=set_grad)
                Δθ = cat(Δθx, Δθcl, Δθ; dims=1)
                ∇logdet = cat(∇logdet_x, ∇logdet_cl, ∇logdet; dims=1)
            end
        end
        ΔZ = wavelet_unsqueeze(ΔZ)
        Z = wavelet_unsqueeze(Z)
    end
    set_grad ? (return ΔZ, Z) : (return ΔZ, Δθ, Z, ∇logdet)
end

## Jacobian-related utils

function jacobian(ΔX, Δθ::Array{Parameter, 1}, X, H::NetworkMultiScaleHINT)
    if H.split_scales
        X_save = Array{Array}(undef, H.L-1, 2)
        ΔX_save = Array{Array}(undef, H.L-1, 2)
    end
    logdet = 0f0
    GNΔθ = Array{Parameter, 1}(undef, 0)
    idxblk = 0
    for i=1:H.L
        X = wavelet_squeeze(X)
        ΔX = wavelet_squeeze(ΔX)
        for j=1:H.K
            npars_ij = 2+length(get_params(H.CL[i, j]))
            Δθij = Δθ[idxblk+1:idxblk+npars_ij]
            ΔX_, X_, logdet1, GNΔθ1 = H.AN[i, j].jacobian(ΔX, Δθij[1:2], X)
            ΔX, X, logdet2, GNΔθ2 = H.CL[i, j].jacobian(ΔX_, Δθij[3:end], X_)
            logdet += (logdet1 + logdet2)
            GNΔθ = cat(GNΔθ, GNΔθ1, GNΔθ2; dims=1)
            idxblk += npars_ij
        end
        if H.split_scales && i < H.L    # don't split after last iteration
            X, Z = tensor_split(X)
            ΔX, ΔZ = tensor_split(ΔX)
            X_save[i] = Z
            ΔX_save[i] = ΔZ
            H.X_dims[i] = size(Z)
        end
    end
    if H.split_scales
        X = cat_states(X_save, X)
        ΔX = cat_states(ΔX_save, ΔX)
    end
    return ΔX, X, logdet, GNΔθ
end

adjointJacobian(ΔZ, Z, H::NetworkMultiScaleHINT) = backward(ΔZ, Z, H; set_grad=false)

## Other utils

# Clear gradients
function clear_grad!(H::NetworkMultiScaleHINT)
    depth = length(H.CL)
    for j=1:depth
        clear_grad!(H.AN[j])
        clear_grad!(H.CL[j])
    end
end

# Get parameters
function get_params(H::NetworkMultiScaleHINT)
    L, K = size(H.CL)
    p = Array{Parameter, 1}(undef, 0)
    for i = 1:L
        for j = 1:K
            p = cat(p, get_params(H.AN[i, j]); dims=1)
            p = cat(p, get_params(H.CL[i, j]); dims=1)
        end
    end
    return p
end


text:
using GroupSlices
using Test

# Tests for groupslices, groupinds, firstinds, lastinds
A = [1 2 3 ; 4 5 6 ; 7 8 9]
B = [11 12 13 ; 14 15 16 ; 17 18 19]
C = [21 22 23 ; 24 25 26 ; 27 28 29]
D = cat(A, B, C, C, B, C, A; dims=3)

ic = [1;2;3;3;2;3;1]
ib = Vector{Int}[]
push!(ib,[1;7])
push!(ib,[2;5])
push!(ib,[3;4;6])
ia = [1;2;3]
ia1 = [1;2;3]
ia2 = [7;5;6]

ic_test = groupslices(D,3)
ib_test = groupinds(ic_test)
ia_test = firstinds(ic_test)
ia1_test = firstinds(ib_test)
ia2_test = lastinds(ib_test)

@test isequal(ic, ic_test)
@test isequal(ib, ib_test)
@test isequal(ia, ia_test)
@test isequal(ia1, ia1_test)
@test isequal(ia2, ia2_test)

E =  [1.0  1.0  2.0  2.0  1.0  1.0  1.0  2.0  2.0  2.0
      1.0  1.0  2.0  2.0  2.0  1.0  2.0  2.0  2.0  1.0
      1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  2.0]

ib = Vector{Int}[[1,2,6],[3,4,8,9],[5,7],[10]]
ia = [1; 3; 5; 10]
ic = [1; 1; 3; 3; 5; 1; 5; 3; 3; 10]
ia1 = [1; 3; 5; 10]
ia2 = [6; 9; 7; 10]

ic_test = groupslices(E,2)
ib_test = groupinds(ic_test)
ia_test = firstinds(ic_test)
ia1_test = firstinds(ib_test)
ia2_test = lastinds(ib_test)

@test isequal(ic, ic_test)
@test isequal(ib, ib_test)
@test isequal(ia, ia_test)
@test isequal(ia1, ia1_test)
@test isequal(ia2, ia2_test)

# Tests when adding keyword dims
V = [1,1,5,6,1,5,9]
@test all(x == V[i] for (x,i) in zip(V, groupslices(V)))

@test groupslices(E; dims=1) == groupslices(E, 1) == groupslices(E)
@test groupslices(E; dims=2) == groupslices(E, 2)


text:
opt = Adagrad()

A = [1.0 2.0; 3.0 4.0]
B = [5. 6. 7.; 8. 9. 0.]
g_t = [A, B]

@testset "Optimizer List" begin
    @test_nowarn update(opt, g_t)
end


text:
"A measure associated with a parameterization "
struct ParameterizationMeasure{P,T} <: DomainIntegrals.Weight{T}
    param   ::  P
end

ParameterizationMeasure(param) = ParameterizationMeasure{typeof(param),eltype(domain(param))}(param)

DomainIntegrals.unsafe_weightfun(p::ParameterizationMeasure, t) = gradientnorm(p.param, t)
DomainIntegrals.support(p::ParameterizationMeasure) = domain(p.param)

export BoundaryIntegralOperator
"A boundary integral operator with a parameterized boundary."
struct BoundaryIntegralOperator{K,P,T} <: IntegralOperator
    kernel      ::  K
    boundary    ::  Domain{T}
    param       ::  P
end

BasisFunctions.name(op::BoundaryIntegralOperator) = "Boundary integral operator"

kernel(op::BoundaryIntegralOperator) = op.kernel
boundary(op::BoundaryIntegralOperator) = op.boundary
parameterization(op::BoundaryIntegralOperator) = op.param
parameterdomain(op::BoundaryIntegralOperator) = domain(parameterization(op))

eval_kernel(op::BoundaryIntegralOperator, x, y) =
    kernel(op)(x, y, parameterization(op))

measure(op::BoundaryIntegralOperator) = ParameterizationMeasure(parameterization(op))

singularity(op::BoundaryIntegralOperator) = singularity(kernel(op))

export op_apply
"Apply the integral operator to the given density function and evaluate at the point t."
function op_apply(op::BoundaryIntegralOperator, density, t)
    μ = measure(op)
    integrand = y -> eval_kernel(op, t, y) * density(y) * unsafe_weightfun(μ, y)
    DomainIntegrals.integral(QuadAdaptive(), integrand, parameterdomain(op), LogSingPoint(t))
end


text:
module test_jive_onlyonce_evaluated

using Test
using Jive
empty!(Jive.onlyonce_evaluated)
empty!(Jive.onlyonce_called)

include(normpath(@__DIR__, "heavy.jl"))
@test val == 0

include(normpath(@__DIR__, "heavy.jl"))
@test val == 42

include(normpath(@__DIR__, "heavy.jl"))
@test val == 42

end # module test_jive_onlyonce_evaluated

module test_jive_onlyonce_called_block

using Test
using Jive
empty!(Jive.onlyonce_evaluated)
empty!(Jive.onlyonce_called)

function f(x)
    @onlyonce begin
        x += 2
    end
    x
end # function f(x)

@test f(1) == 3
@test f(1) == 1
@test f(1) == 1

end # module test_jive_onlyonce_called_block

module test_jive_onlyonce_called_call

using Test
using Jive
empty!(Jive.onlyonce_evaluated)
empty!(Jive.onlyonce_called)

function f(x)
    @onlyonce(x += 2)
    x
end # function f(x)

@test f(1) == 3
@test f(1) == 1
@test f(1) == 1

end # module test_jive_onlyonce_called_call


text:
import Base: getindex, keys, read

struct Board
    values::Matrix{Int}
    marked::Set{CartesianIndex{2}}

    Board(v) = new(v, Set())
end

getindex(board::Board, idx) = getindex(board.values, idx)
keys(board::Board) = keys(board.values)

function read_input(f::AbstractString)
    return open(f, "r") do io
        numbers = parse.(Int, split(readline(io), ','))
        boards = []
        while !eof(io)
            # skip empty line
            @assert readline(io) == ""
            push!(boards, read(io, Board))
        end
        return ([boards...], numbers)
    end
end

function read(io::IO, ::Type{Board}; dims = (5, 5))
    board = vcat([let line = split(readline(io), r"\s+"; keepempty = false)
                      parse.(Int, line)'
                  end for _ in 1:dims[1]]...)
    @assert size(board) == dims
    return Board(board)
end

function mark!(board::Board, number::Int)::Board
    union!(board.marked, findall(==(number), board.values))
    return board
end

const ROWS_AND_COLUMNS = vcat([[[CartesianIndex(i, j) for j in 1:5],
                                [CartesianIndex(j, i) for j in 1:5]]
                               for i in 1:5]...)

function winning_score(board::Board)::Union{Nothing, Int}
    return any(fs -> fs ⊆ board.marked, ROWS_AND_COLUMNS) ?
        sum(board[setdiff(keys(board), board.marked)]) :
        nothing
end

function find_first_winner!(boards::AbstractVector{<: Board},
                            numbers::AbstractVector{<: Integer}
                           )::Union{Nothing, Int}
    for n in numbers, b in boards
        mark!(b, n)
        score = winning_score(b)
        if score !== nothing
            return score * n
        end
    end
    return nothing
end

function find_last_winner!(boards::AbstractVector{<: Board},
                           numbers::AbstractVector{<: Integer}
                          )::Union{Nothing, Int}
    winners = findall(b -> winning_score(b) !== nothing, boards)
    playing = setdiff(keys(boards), winners)
    for n in numbers, p in copy(playing)
        b = boards[p]
        mark!(b, n)
        score = winning_score(b)
        if score !== nothing
            deleteat!(playing, findfirst(==(p), playing))
            if isempty(playing)
                return score * n
            end
        end
    end
    return nothing
end

let (bs, ns) = read_input("day04-test.in")
    @assert find_first_winner!(bs, ns) == 4512
    @assert find_last_winner!(bs, ns) == 1924
end

(bs, ns) = read_input("day04.in")
println("Q1: ", find_first_winner!(bs, ns))
println("Q2: ", find_last_winner!(bs, ns))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
using Makie, GeometryTypes
using PortAudio

N = 1024 # size of audio read
N2 = N ÷ 2 + 1 # size of rfft output
D = 200 # number of bins to display
M = 100 # number of lines to draw
S = 0.5 # motion speed of lines
src = PortAudioStream(1, 2)
buf = Array{Float32}(N)
fftbuf = Array{Complex{Float32}}(N2)
magbuf = Array{Float32}(N2)
fftplan = plan_rfft(buf; flags = FFTW.EXHAUSTIVE)

scene = Scene(resolution = (500, 500))
ax = axis(0:0.1:1, 0:0.1:1, 0:0.1:0.5)
center!(scene)

ls = map(1:M) do _
    yoffset = to_node(to_value(scene[:time]))
    offset = lift_node(scene[:time], yoffset) do t, yoff
        Point3f0(0.0f0, (t - yoff) * S, 0.0f0)
    end
    l = lines(
        linspace(0, 1, D),
        0.0f0,
        zeros(Float32, D),
        offset = offset,
        color = (:black, 0.1),
    )
    (yoffset, l)
end

while isopen(scene[:screen])
    for (yoffset, line) in ls
        isopen(scene[:screen]) || break
        read!(src, buf)
        A_mul_B!(fftbuf, fftplan, buf)
        @. magbuf = log(clamp(abs(fftbuf), 0.0001, Inf)) / 10 + 0.5
        line[:z] = magbuf[1:D]
        push!(yoffset, to_value(scene[:time]))
    end
end


text:
struct Arguments
    positional::Vector{Any}
    named::Dict{Symbol, Any}
end

Arguments(v::AbstractVector) = Arguments(v, Dict{Symbol, Any}())

function arguments(args...; kwargs...)
    positional = collect(Any, args)
    named = Dict{Symbol, Any}(kwargs)
    return Arguments(positional, named)
end

Base.get(args::Arguments, i::Int, default) = get(args.positional, i, default)
Base.get(args::Arguments, sym::Symbol, default) = get(args.named, sym, default)

Base.getindex(args::Arguments, i::Int) = args.positional[i]
Base.getindex(args::Arguments, sym::Symbol) = args.named[sym]

Base.setindex!(args::Arguments, val, i::Int) = (args.positional[i] = val)
Base.setindex!(args::Arguments, val, sym::Symbol) = (args.named[sym] = val)

Base.pop!(args::Arguments, i::Int, default) = pop!(args.positional, i, default)
Base.pop!(args::Arguments, sym::Symbol, default) = pop!(args.named, sym, default)

Base.copy(args::Arguments) = Arguments(copy(args.positional), copy(args.named))

function Base.map(f, a::Arguments, as::Arguments...)
    is = eachindex(a.positional)
    ks = keys(a.named)
    function g(i)
        vals = map(t -> t[i], (a, as...))
        return f(vals...)
    end
    positional = collect(Any, Iterators.map(g, is))
    named = Dict{Symbol, Any}(k => g(k) for k in ks)
    return Arguments(positional, named)
end

function Base.mergewith!(op, a::Arguments, b::Arguments)
    la, lb = length(a.positional), length(b.positional)
    for i in 1:lb
        (i ≤ la) ? (a[i] = op(a[i], b[i])) : push!(a.positional, b[i])
    end
    mergewith!(op, a.named, b.named)
    return a
end

latter(_, b) = b

Base.merge!(a::Arguments, b::Arguments) = mergewith!(latter, a, b)

Base.merge(a::Arguments, b::Arguments) = merge!(copy(a), b)

function separate!(continuous::Arguments)
    discrete = Dict{Symbol, Any}()
    for (k, v) in continuous.named
        label, value = getlabel(v), getvalue(v)
        iscontinuous(value) && continue
        discrete[k] = Labeled(label, first(value))
        pop!(continuous.named, k)
    end
    return NamedTuple(discrete) => continuous
end

function recombine!(discrete, continuous)
    for (k, v) in pairs(discrete)
        label, value = getlabel(v), getvalue(v)
        continuous[k] = Labeled(label, fill(value))
    end
    return continuous
end


text:
import  Base: getindex, setindex!, push!, keys

export ConvergenceHistory
export nprods, niters, nrests

########
# Type #
########

"""
Store general and in-depth information about an iterative method.

**Fields**

* `mvps::Int`: number of matrix vector products.
* `mtvps::Int`: number of transposed matrix-vector products
* `iters::Int`: iterations taken by the method.
* `restart::T`: restart relevant information.
    - `T == Int`: iterations per restart.
    - `T == Void`: methods without restarts.

* `isconverged::Bool`: convergence of the method.
* `data::Dict{Symbol,Any}`: Stores all the information stored during the method execution.
It stores tolerances, residuals and other information, e.g. ritz values in [svdl](@ref).

**Constructors**

    ConvergenceHistory()
    ConvergenceHistory(restart)

Create `ConvergenceHistory` with empty fields.

**Arguments**

* `restart`: number of iterations per restart.

**Plots**

Supports plots using the `Plots.jl` package via a type recipe. Vectors are
ploted as series and matrices as scatterplots.

**Implements**

* `Base`: `getindex`, `setindex!`, `push!`

"""
type ConvergenceHistory{T,K}
    mvps::Int
    mtvps::Int
    iters::Int
    restart::K
    isconverged::Bool
    data::Dict{Symbol, Any}
end
function ConvergenceHistory(;restart=nothing, partial=false)
    ConvergenceHistory{!partial,typeof(restart)}(
        0,0,0,restart,false,Dict{Symbol, Any}()
        )
end

"""
Stores information of the current iteration.
"""
typealias PartialHistory ConvergenceHistory{false}

"""
Stores the information of all the iterations.
"""
typealias CompleteHistory ConvergenceHistory{true}

"""
History without resets.
"""
typealias UnrestartedHistory{T} ConvergenceHistory{T, Void}

"""
History with resets.
"""
typealias RestartedHistory{T} ConvergenceHistory{T, Int}

#############
# Functions #
#############

"""
    getindex(ch, s)

Get collection or tolerance associated with key `s` in `ch::ConvergenceHistory`.

    getindex(ch, s, kwargs...)

Access elements of the collection associated with key `s` in `ch::ConvergenceHistory`.
"""
getindex(ch::ConvergenceHistory, s::Symbol) = ch.data[s]
getindex(ch::ConvergenceHistory, s::Symbol, kwargs...) = ch.data[s][kwargs...]

"""
    setindex!(ch, tol, s)

Set tolerance value associated with `s` in `ch::ConvergenceHistory` to `tol`.

    setindex!(ch, val, s, kwargs...)

Set collection element associated with key `s` in `ch::ConvergenceHistory` to val.
"""
setindex!(ch::ConvergenceHistory, val, s::Symbol) = ch.data[s] = val
setindex!(ch::ConvergenceHistory, val, s::Symbol, kwargs...) = ch.data[s][kwargs...] = val

"""
    push!(ch, key, val)

Push contents of `val` to collection associated with `key` in `ch::ConvergenceHistory`.
"""
function push!(ch::ConvergenceHistory, key::Symbol, vec::Union{Vector,Tuple})
    matrix = ch.data[key]
    width = size(matrix,2)
    iter = isa(ch,CompleteHistory) ? ch.iters : 1
    base = (iter-1)*width
    for i in 1:min(width,length(vec))
        matrix[base+i] = vec[i]
    end
end
push!(ch::ConvergenceHistory, key::Symbol, data) = push_custom_data!(ch, key, data)

push_custom_data!(ch::PartialHistory, key::Symbol, data) = ch.data[key] = data
push_custom_data!(ch::CompleteHistory, key::Symbol, data) = ch.data[key][ch.iters] = data

"""
    reserve!(ch, key, maxiter)
    reserve!(typ, key, maxiter)
    reserve!(ch, key, maxiter, size)
    reserve!(typ, ch, key, maxiter, size)

Reserve space for per iteration data in `ch`. If size is provided, intead of a
vector it will reserve matrix of dimensions `(maxiter, size)`.

**Arguments**

* `typ::Type`: Type of the elements to store. Defaults to `Float64` when not given.
* `ch::ConvergenceHistory`: convergence history.
* `key::Union{Symbol,Vector{Symbol}}`: key used to identify the data.
* `maxiter::Int`: number of iterations to save space for.
* `size::Int`: number of elements to store with the `key` identifier.

"""
function reserve!(ch::ConvergenceHistory, keys::Vector{Symbol}, kwargs...)
    for key in keys
        reserve!(ch, key, kwargs...)
    end
end
function reserve!(ch::ConvergenceHistory, key::Symbol, kwargs...)
    reserve!(Float64, ch, key, kwargs...)
end
function reserve!(typ::Type, ch::ConvergenceHistory, key::Symbol, kwargs...)
    _reserve!(typ, ch, key, kwargs...)
end

#If partialhistory, theres no need to store a vector or matrix, instead
#store nothing or store a vector respectively.
_reserve!(typ::Type, ch::PartialHistory, key::Symbol, ::Int) = nothing
function _reserve!(typ::Type, ch::PartialHistory, key::Symbol, ::Int, size::Int)
    ch.data[key] = Vector{typ}(size)
end
function _reserve!(typ::Type, ch::CompleteHistory, key::Symbol, len::Int)
    ch.data[key] = Vector{typ}(len)
end
function _reserve!(typ::Type, ch::CompleteHistory, key::Symbol, len::Int, size::Int)
    ch.data[key] = Matrix{typ}(len, size)
end

"""
    shrink!(ml)
shrinks the reserved space for `ConvergenceHistory` `ch` to the space actually used to log.
"""
shrink!(::PartialHistory) = nothing
function shrink!(ch::CompleteHistory)
    for key in keys(ch)
        elem = ch.data[key]
        if isa(elem, Vector)
            resize!(elem, ch.iters)
        elseif isa(elem, Matrix)
            ch.data[key] = elem[1:ch.iters, :]
        end
    end
end

"""
    nextiter!(ml)

Adds one the the number of iterations in [`ConvergenceHistory`](@ref) `ch`. This is
necessary to avoid overwriting information with `push!(ml)`. It is also able
to update other information of the method.
"""
function nextiter!(ch::ConvergenceHistory; mvps=0,mtvps=0)
    ch.iters+=1
    ch.mvps+=mvps
    ch.mtvps+=mtvps
end

"""
    keys(ch)

Key iterator of the per iteration data logged in `ConvergenceHistory` `ch`.
"""
keys(ch::ConvergenceHistory) = keys(ch.data)

"""
    setconv(ml, val)

Set `val` as convergence status of the method in [`ConvergenceHistory`](@ref) ch.
"""
setconv(ch::ConvergenceHistory, val::Bool) = ch.isconverged=val

"""
    nprods(ch)

Number of matrix-vector products plus transposed matrix-vector products
logged in `ConvergenceHistory` `ch`.
"""
nprods(ch::ConvergenceHistory) = ch.mvps+ch.mtvps

"""
    niters(ch)

Number of iterations logged in `ConvergenceHistory` `ch`.
"""
niters(ch::ConvergenceHistory) = ch.iters

"""
    nrests(ch)

Number of restarts logged in `ConvergenceHistory` `ch`.
"""
nrests(ch::RestartedHistory) = Int(ceil(ch.iters/ch.restart))


text:
# Timestep conversion function
function getindexfromyear_dice_2013(year)
    baseyear = 2010

    if rem(year - baseyear, 5) != 0
        error("Invalid year")
    end

    return div(year - baseyear, 5) + 1
end

# Get parameters from DICE2013 excel sheet
# range is the range of cell values on the excel sheet and must be a string, "B56:B77"
# parameters = :single for just one value, or :all for entire time series
# sheet is the sheet in the excel file to reference (i.e. "Base")
# T is the length of the time period (i.e 60)

# example:   getparams("B15:BI15", :all, "Base",  60)

function getparams(f, range::String, parameters::Symbol, sheet::String, T)

    if parameters == :single
        data = f[sheet][range]
        vals = Float64(data[1])

    elseif parameters == :all
        data = f[sheet][range]
        s = size(data)

        if length(s) == 2 && s[1] == 1
            # convert 2D row vector to 1D col vector
            data = vec(data)
        end

        dims = length(size(data))
        vals = Array{Float64,dims}(data)
    end

    return vals
end


text:
using WAV
using IJulia

function inline_audioplayer(filepath)
  markup = """<audio controls="controls" {autoplay}>
              <source src="$filepath" type="audio/wav" />
              Your browser does not support the audio element.
              </audio>"""
  display(MIME("text/html") ,markup)
end

#function inline_audioplayer(s, fs)
#    buf = IOBuffer()
#    wavwrite(s, buf; Fs=fs)
#    data = base64(unsafe_string(buf))
#    markup = """<audio controls="controls" {autoplay}>
#                <source src="data:audio/wav;base64,$data" type="audio/wav" />
#                Your browser does not support the audio element.
#                </audio>"""
#    display(MIME("text/html") ,markup)
#end


text:
using .AxisArrays: AxisArrays, AxisArray

dimname(a::AxisArray, i) = AxisArrays.axisnames(a)[i]
dimnames(a::AxisArray) = AxisArrays.axisnames(a)
dimvals(a::AxisArray, i) = AxisArrays.axisvalues(a)[i]
iscontdim(a::AxisArray, i) = AxisArrays.axistrait(AxisArrays.axes(a,i)) <: AxisArrays.Dimensional
getdata(a::AxisArray) = parent(a)
function yaxcreate(::Type{<:AxisArray}, data, dnames, dvals, atts)
  d = ntuple(ndims(data)) do i
    dnames[i] => dvals[i]
  end
  AxisArray(data; d...)
end


text:
using SimpleStructs
using Base.Test

# write your own tests here
abstract DummyType

################################################################################
# default values, type assert, value constraints
################################################################################
@defstruct Struct1 <: DummyType (
  i1_field   :: Int = 0,
  (pos_field :: Float64 = 1, pos_field > 0),
  (i2_field  :: Int = 0, i2_field >= i1_field)
)

struct1 = Struct1()
@test isa(struct1, Struct1)
@test isa(struct1, DummyType)

@test struct1.i1_field == 0

@test isa(Struct1(i2_field = 2.0), Struct1)
@test_throws InexactError Struct1(i2_field = 2.5)
@test_throws AssertionError Struct1(i2_field = 2.0, i1_field = 3)

@test_throws MethodError Struct1(i1_field = "foo")

################################################################################
# required fields
################################################################################
@defstruct Struct2 (
  i_field   :: Int,
  (i_field2 :: Int, i_field2 > 0),
  i_field3  :: Int = -1
)

@test_throws AssertionError Struct2()
@test_throws AssertionError Struct2(i_field=0)
@test_throws AssertionError Struct2(i_field=0, i_field2=0)
@test isa(Struct2(i_field=0, i_field2=1), Struct2)

################################################################################
# convert is called on parameters, more user-friendly than type matching
################################################################################
@defstruct Struct3 (
  field1 :: AbstractString = "",
  field2 :: Vector{AbstractString} = []
)

struct3 = Struct3(field1 = "ascii string")
struct3 = Struct3(field1 = "unicode α-β")
struct3 = Struct3(field2 = ["vector", "of", "ASCII", "string"])

################################################################################
# immutable structs
################################################################################
@defimmutable Struct4 (
  field :: Int = 0
)

struct4 = Struct4()
@test_throws ErrorException (struct4.field = 2)

################################################################################
# type parameters
################################################################################
@defimmutable Struct5{R<:Real, T<:Real} (
  field  :: Vector{R},
  field2 :: Vector{T} = zeros(5),
)
@test_throws MethodError struct5 = Struct5()
struct5 = Struct5(field = [1,2,3])
@test typeof(struct5) <: Struct5{Int, Float64}
@test typeof(struct5.field) <: Vector{Int}
@test struct5.field == [1,2,3]
@test typeof(struct5.field2) <: Vector{Float64}
@test struct5.field2 == [0.,0,0,0,0]

################################################################################
# type parameters with base class
################################################################################
@defimmutable Struct6{R<:Real, T<:Real} <: AbstractVector{R} (
  field  :: Vector{R},
  field2 :: Vector{T} = zeros(5),
)
@test_throws MethodError struct6 = Struct6()
struct6 = Struct6(field = [-1,2,3], field2 = [0., 1.])
@test typeof(struct6) <: Struct6{Int, Float64}
@test typeof(struct6) <: AbstractVector{Int}
@test typeof(struct6.field) <: Vector{Int}
@test struct6.field == [-1,2,3]
@test typeof(struct6.field2) <: Vector{Float64}
@test struct6.field2 == [0.,1]

################################################################################
# Base.show
################################################################################
@defimmutable Struct7 (
  f :: Int = 3,
  f2:: AbstractString = "hello"
)
buf = IOBuffer()
show(buf, Struct7())
@test takebuf_string(buf) == """Struct7(f=3, f2="hello")"""

@defimmutable Struct8{T, T2 <: Real} <: AbstractVector{T2} (
  f1 :: T,
  f2 :: T2 = 0.5
)
buf = IOBuffer()
show(buf, Struct8(f1=one(Int64)))
@test takebuf_string(buf) == "Struct8{Int64,Float64}(f1=1, f2=0.5)"


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# High-level utilities

#####################################################################

"""
    tonemapping(input_file::String,
                output_file::String
                ; α::Float32,
                  γ::Float32,
                  luminosity::Union{Float32, Nothing} = nothing,
                  disable_output::Bool = false)

Load a pfm hdr image, apply the tone mapping, and save the generated ldr image.

`input_file` is the path to a valid PFM image. `output_file` is the path in which save the generated image + the name
of the output file. The output format is deduced by the extension of the output file: if the name is `example.jpg`,
the image will be exported in JPEG format. `α` is the normalization coefficient. `γ` is the actor of the γ correction.
Optional parameter `luminosity` can be used for better tuning.

If `disable_output` is `true`, no message is printed.
"""
function tonemapping(input_file::String,
                     output_file::String
                     ; α::Float32 = 0.5f0,
                       γ::Float32 = 1f0,
                       luminosity::Union{Float32, Nothing} = nothing,
                       disable_output::Bool = false)
    # check if valid input_file
    if !(typeof(query(input_file))<:File{format"PFM"})
        error("'$input_file' is not a pfm file!")
    end
    # apply tonemapping
    io = disable_output ? devnull : stdout
    println(io, "\n-> TONE MAPPING PROCESS")
    print(io, "Loading input file '$(input_file)'...")
    image = load(input_file) |> HdrImage
    replace!(image.pixel_matrix, RGB{Float32}(NaN, NaN, NaN) => BLACK)
    print(io, " done!\nApplying tone mapping... ")
    image = (isnothing(luminosity) ? normalize(image, α) : normalize(image, α, luminosity=luminosity)) |> clamp
    image = γ_correction(image, γ)
    # save image
    print(io, " done!\nSaving final image to \"$(output_file)\"...")
    save(output_file, image.pixel_matrix)
    println(io, " done!")
end

#####################################################################

"""
    render(image_tracer::ImageTracer,
           renderer::Renderer
           ; output_file::String = "out.pfm"
             use_threads::Bool = true,
             disable_output::Bool = false)

Render an image given an [`ImageTracer`](@ref) and a [`Renderer`](@ref), and save the generated hdr image in `output_file`.

If `use_threads` is `true`, use macro `Threads.@threads`. If `disable_output` is `true`, no message is printed.
"""
function render(image_tracer::ImageTracer,
                renderer::Renderer
                ; output_file::String = "out.pfm",
                  use_threads::Bool = true,
                  disable_output::Bool = false)
    # check if valid output_file
    if !(typeof(query(output_file))<:File{format"PFM"})
        error("'$output_file' is not a pfm file!")
    end
    # render image
    io = disable_output ? devnull : stdout
    println(io, "\n-> RENDERING")
    println(io, "Rendering image...")
    fire_all_rays!(image_tracer, renderer, use_threads=use_threads, enable_progress_bar=!disable_output)
    print(io, "Saving pfm image...")
    save(output_file, permutedims(image_tracer.image.pixel_matrix))
    println(io, " done!")
end

"""
    render_from_script(input_script::String
                       ; output_file::String = "out.pfm",
                         time::Float32 = 0f0,
                         vars::String = "",
                         use_threads::Bool = true,
                         disable_output::Bool = false)

Render an image given an `input_script` written in SceneLang and save the generated hdr image in `output_file`.

If `use_threads` is `true`, use macro `Threads.@threads`. If `disable_output` is `true`, no message is printed.
"""
function render_from_script(input_script::String
                            ; output_file::String = "out.pfm",
                              time::Float32 = 0f0,
                              vars::String = "",
                              use_threads::Bool = true,
                              disable_output::Bool = false)
    # check if valid output_file
    if !(typeof(query(output_file))<:File{format"PFM"})
        error("'$output_file' is not a pfm file!")
    end
    # parse scene
    scene = Scene(time=time, variables=parse_variables_from_string(vars))
    scene = open_stream(input_script) do stream
        try
            parse_scene(stream, scene)
        catch e
            isa(e, InterpreterException) || rethrow(e)
            print_subsequent_lexer_exceptions(stream, e) |> rethrow
        end
    end
    needs_lights = !isnothing(scene.renderer) && scene.renderer.type ∈ (PointLightRenderer,)
    let are_not_defined = (isempty(scene.world)                  => "No shapes have been spawned.",
                           needs_lights && isempty(scene.lights) => "No lights have been spawned",
                           isnothing(scene.camera)               => "No camera is being used.",
                           isnothing(scene.image)                => "No image is being used",
                           isnothing(scene.renderer)             => "No renderer is being used",
                           isnothing(scene.tracer)               => "No tracer is being used")
        if any(first.(are_not_defined))
            throw(UndefinedSetting(SourceLocation(file_name = input_script), """
            One or more necessary settings have not bene set by a USING or SPAWN command in the given SceneLang script.
            in particular:\n\t$(join(last.(filter(first, are_not_defined)), "\n\t"))
            """))
        end
    end
    # extract render informations
    renderer = if needs_lights
        scene.renderer.type(scene.world, scene.lights; scene.renderer.kwargs...)
    else
        scene.renderer.type(scene.world; scene.renderer.kwargs...)
    end
    image_tracer = ImageTracer(scene.image, scene.camera; scene.tracer.kwargs...)
    # render
    render(image_tracer, renderer, output_file=output_file, use_threads=use_threads, disable_output=disable_output)
end


text:
struct MedianVarianceBinning <: AbstractBinningAlgorithm
    minsplit::Int
    maxbins::Int
end

MedianVarianceBinning() = MedianVarianceBinning(10)
MedianVarianceBinning(minsplit::Int) = MedianVarianceBinning(minsplit, typemax(Int))

function perform(alg::MedianVarianceBinning, predictions::AbstractMatrix{<:Real},
                 labels::AbstractVector{<:Integer})
    @unpack maxbins, minsplit = alg

    # check arguments
    nclasses, nsamples = check_size(predictions, labels)

    # check first if binning is possible
    (nsamples < minsplit || maxbins <= 1) &&
        return Dict{Int,Bin{eltype(predictions)}}(1 => Bin(predictions, labels))

    # create empty priority queue and set of bins
    queue = PriorityQueue{Tuple{typeof(predictions),typeof(labels),Int},
                          Float64}(Base.Order.Reverse)
    bins = Dict{Int,Bin{eltype(predictions)}}()

    # add data set to priority queue
    varfull = vec(var(predictions; dims = 2))
    maxvar, maxvardim = findmax(varfull)
    enqueue!(queue, (predictions, labels, maxvardim), maxvar)
    nbins = 1

    while nbins < maxbins && !isempty(queue)
        # pick set with largest variance
        predictions, labels, maxvardim = dequeue!(queue)

        # check if minimum split size is reached
        nsamples = size(predictions, 2)

        # compute index of median element
        mid = median(view(predictions, maxvardim, :))

        # compute indices of subsets
        idxsbelow = filter(i -> predictions[maxvardim, i] < mid, 1:nsamples)
        idxsabove = filter(i -> predictions[maxvardim, i] >= mid, 1:nsamples)

        # if splitting is not possible fix remaining partition
        if isempty(idxsbelow)
            bins[length(bins)+1] = Bin(predictions[:, idxsabove], labels[idxsabove])

            # there was only one set left, so we are done
            continue
        end
        if isempty(idxsabove)
            bins[length(bins)+1] = Bin(predictions[:, idxsbelow], labels[idxsbelow])

            # there was only one set left, so we are done
            continue
        end

        # otherwise create subsets and update the queue or list of bins
        for idxs in (idxsbelow, idxsabove)
            # create subsets
            newpredictions = predictions[:, idxs]
            newlabels = labels[idxs]

            # try to avoid needless calculations
            if length(newlabels) < minsplit
                bins[length(bins)+1] = Bin(newpredictions, newlabels)
            else
                newvar = vec(var(newpredictions; dims = 2))
                newmaxvar, newmaxvardim = findmax(newvar)
                enqueue!(queue, (newpredictions, newlabels, newmaxvardim), newmaxvar)
            end

            # we created one additional partition in total
            nbins += 1
        end
    end

    # add remaining bins
    while !isempty(queue)
        # pop queue
        predictions, labels = dequeue!(queue)

        # create bin
        bins[length(bins)+1] = Bin(predictions, labels)
    end

    bins
end


text:
# Julia wrapper for header: /usr/include/scip/cons_indicator.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrIndicator(scip)
    ccall((:SCIPincludeConshdlrIndicator, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPcreateConsIndicator(scip, cons, name, binvar, nvars, vars, vals, rhs, initial, separate, enforce, check, propagate, _local, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsIndicator, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{SCIP_VAR}, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Cdouble, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, binvar, nvars, vars, vals, rhs, initial, separate, enforce, check, propagate, _local, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicIndicator(scip, cons, name, binvar, nvars, vars, vals, rhs)
    ccall((:SCIPcreateConsBasicIndicator, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{SCIP_VAR}, Cint, Ptr{Ptr{SCIP_VAR}}, Ptr{Cdouble}, Cdouble), scip, cons, name, binvar, nvars, vars, vals, rhs)
end

function SCIPcreateConsIndicatorLinCons(scip, cons, name, binvar, lincons, slackvar, initial, separate, enforce, check, propagate, _local, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsIndicatorLinCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{SCIP_VAR}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, binvar, lincons, slackvar, initial, separate, enforce, check, propagate, _local, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicIndicatorLinCons(scip, cons, name, binvar, lincons, slackvar)
    ccall((:SCIPcreateConsBasicIndicatorLinCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{SCIP_VAR}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}), scip, cons, name, binvar, lincons, slackvar)
end

function SCIPaddVarIndicator(scip, cons, var, val)
    ccall((:SCIPaddVarIndicator, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, Cdouble), scip, cons, var, val)
end

function SCIPgetLinearConsIndicator(cons)
    ccall((:SCIPgetLinearConsIndicator, libscip), Ptr{SCIP_CONS}, (Ptr{SCIP_CONS},), cons)
end

function SCIPsetLinearConsIndicator(scip, cons, lincons)
    ccall((:SCIPsetLinearConsIndicator, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_CONS}), scip, cons, lincons)
end

function SCIPsetBinaryVarIndicator(scip, cons, binvar)
    ccall((:SCIPsetBinaryVarIndicator, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}), scip, cons, binvar)
end

function SCIPgetBinaryVarIndicator(cons)
    ccall((:SCIPgetBinaryVarIndicator, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_CONS},), cons)
end

function SCIPgetSlackVarIndicator(cons)
    ccall((:SCIPgetSlackVarIndicator, libscip), Ptr{SCIP_VAR}, (Ptr{SCIP_CONS},), cons)
end

function SCIPsetSlackVarUb(scip, cons, ub)
    ccall((:SCIPsetSlackVarUb, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Cdouble), scip, cons, ub)
end

function SCIPisViolatedIndicator(scip, cons, sol)
    ccall((:SCIPisViolatedIndicator, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_SOL}), scip, cons, sol)
end

function SCIPmakeIndicatorFeasible(scip, cons, sol, changed)
    ccall((:SCIPmakeIndicatorFeasible, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_SOL}, Ptr{UInt32}), scip, cons, sol, changed)
end

function SCIPmakeIndicatorsFeasible(scip, conshdlr, sol, changed)
    ccall((:SCIPmakeIndicatorsFeasible, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{SCIP_SOL}, Ptr{UInt32}), scip, conshdlr, sol, changed)
end

function SCIPaddLinearConsIndicator(scip, conshdlr, lincons)
    ccall((:SCIPaddLinearConsIndicator, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{SCIP_CONS}), scip, conshdlr, lincons)
end

function SCIPaddRowIndicator(scip, conshdlr, row)
    ccall((:SCIPaddRowIndicator, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{SCIP_ROW}), scip, conshdlr, row)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule x264_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("x264")
JLLWrappers.@generate_main_file("x264", UUID("1270edf5-f2f9-52d2-97e9-ab00b5d0237a"))
end  # module x264_jll


text:
using Test
using SafeTestsets

@safetestset "Rings Tests" begin include("Rings_test.jl") end
@safetestset "RWGSheet Tests" begin include("RWGSheet_test.jl") end
@safetestset "Elements Tests" begin include("Elements_test.jl") end
@safetestset "RWGData Tests" begin include("RWGData_test.jl") end
@safetestset "PGF Tests" begin include("PGF_test.jl") end
@safetestset "Zint Tests" begin include("Zint_test.jl") end
@safetestset "GSMs Tests" begin include("GSMs_test.jl") end
@safetestset "Modes Tests" begin include("Modes_test.jl") end
@safetestset "FillZY Tests" begin include("FillZY_test.jl") end
@safetestset "SymStripFull Tests" begin include("full_test.jl") end


text:
@testset "MagRef" begin
    @test isequal(MagRef(), Mag())

    x = Arb()
    Arblib.set!(Arblib.radref(x), 2)
    @test Arblib.radref(x) == Mag(2)

    y = Mag(Arblib.radref(x))
    Arblib.set!(y, 3)
    @test Arblib.radref(x) == Mag(2)

    z = Arblib.radref(x)[]
    Arblib.set!(z, 4)
    @test Arblib.radref(x) == Mag(2)

    @test isequal(zero(MagRef), zero(Mag))
    @test isequal(zero(Arblib.radref(x)), zero(Mag))
    @test isequal(one(MagRef), one(Mag))
    @test isequal(one(Arblib.radref(x)), one(Mag))
end

@testset "ArfRef" begin
    @test ArfRef() isa Arf
    @test isequal(ArfRef(), Arf())
    @test precision(ArfRef()) == Arblib.DEFAULT_PRECISION[]
    @test precision(ArfRef(prec = 80)) == 80

    x = Arb()
    ptr = Arblib.midref(x).arf_ptr
    parent = Arblib.parentstruct(x)
    @test ArfRef(ptr, parent) == Arf()
    @test precision(ArfRef(ptr, parent)) == Arblib.DEFAULT_PRECISION[]
    @test precision(ArfRef(ptr, parent, prec = 80)) == 80

    Arblib.set!(Arblib.midref(x), 2)
    @test Arblib.midref(x) == Arf(2)
    y = Arf(Arblib.midref(x))
    Arblib.set!(y, 3)
    @test Arblib.midref(x) == Arf(2)
    z = Arblib.midref(x)[]
    Arblib.set!(z, 4)
    @test Arblib.midref(x) == Arf(2)

    @test isequal(zero(ArfRef), zero(Arf))
    @test isequal(zero(Arblib.midref(x)), zero(Arf))
    @test isequal(one(ArfRef), one(Arf))
    @test isequal(one(Arblib.midref(x)), one(Arf))
end

@testset "ArbRef" begin
    @test ArbRef() isa Arb
    @test isequal(ArbRef(), Arb())
    @test precision(ArbRef()) == Arblib.DEFAULT_PRECISION[]
    @test precision(ArbRef(prec = 80)) == 80

    x = Acb()
    ptr = Arblib.realref(x).arb_ptr
    parent = Arblib.parentstruct(x)
    @test ArbRef(ptr, parent) == Arb()
    @test precision(ArbRef(ptr, parent)) == Arblib.DEFAULT_PRECISION[]
    @test precision(ArbRef(ptr, parent, prec = 80)) == 80

    Arblib.set!(Arblib.realref(x), 2)
    @test Arblib.realref(x) == Arb(2)
    Arblib.set!(Arblib.imagref(x), 3)
    @test Arblib.imagref(x) == Arb(3)
    y = Arb(Arblib.realref(x))
    Arblib.set!(y, 3)
    @test Arblib.realref(x) == Arf(2)
    z = Arblib.realref(x)[]
    Arblib.set!(z, 4)
    @test Arblib.realref(x) == Arf(2)

    @test isequal(zero(ArbRef), zero(Arb))
    @test isequal(zero(Arblib.realref(x)), zero(Arb))
    @test isequal(one(ArbRef), one(Arb))
    @test isequal(one(Arblib.realref(x)), one(Arb))
end

@testset "AcbRef" begin
    @test AcbRef() isa Acb
    @test isequal(AcbRef(), Acb())
    @test precision(AcbRef()) == Arblib.DEFAULT_PRECISION[]
    @test precision(AcbRef(prec = 80)) == 80

    v = AcbRefVector([0])
    x = v[1]
    ptr = x.acb_ptr
    parent = nothing
    @test AcbRef(ptr, parent) == Acb()
    @test precision(AcbRef(ptr, parent)) == Arblib.DEFAULT_PRECISION[]
    @test precision(AcbRef(ptr, parent, prec = 80)) == 80

    @test isequal(zero(AcbRef), zero(Acb))
    @test isequal(zero(x), zero(Acb))
    @test isequal(one(AcbRef), one(Acb))
    @test isequal(one(x), one(Acb))

    Arblib.set!(x, 2)
    @test v[1] == Acb(2)
    y = Acb(x)
    Arblib.set!(y, 3)
    @test v[1] == Acb(2)
    z = x[]
    Arblib.set!(z, 4)
    @test v[1] == Acb(2)
end

@testset "ArfRef from ArbRef" begin
    v = ArbRefVector([1.0, 2.0, 3.0])
    @test v[3] isa ArbRef
    x = Arblib.midref(v[3])
    @test x isa ArfRef
    @test startswith(sprint(show, x), "3")
    y = typeof(x)()
    @test y isa Arf
    y[] = x
    @test startswith(sprint(show, y), "3")
    @test y == x
end

@testset "Refs from AcbRef" begin
    v = AcbRefVector([1, 2, 3])
    @test v[1] isa AcbRef
    @test Arblib.realref(v[1]) isa ArbRef
    @test Arblib.midref(Arblib.realref(v[1])) isa ArfRef
    @test Arblib.imagref(v[1]) isa ArbRef
    @test ComplexF64(v[1]) == 1
end


text:
const lx = SmartGameFormat.Lexer

function read_all_token(str)
    io = lx.TokenStream(IOBuffer(str))
    tkn = lx.Token[]
    while !eof(io)
        push!(tkn, read(io, lx.Token))
    end
    tkn
end

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

@testset "public interface" begin
    # nothing should be exported
    @test_throws UndefVarError Lexer
    @test_throws UndefVarError Token
    @test_throws UndefVarError CharStream
    @test_throws UndefVarError TokenStream
    @test_throws UndefVarError peek
    @test_throws UndefVarError next_token
end

@testset "lexical error" begin
    # all kind of characters that should not be encountered
    # except inside property delimiter "[...]"
    for str in [string.(0:9)...,
                string.('a':'z')...,
                split("'!\"§\$%&/?]\\,<>.:-_#+*","")...]
        io = lx.CharStream(IOBuffer(str))
        @test_throws lx.LexicalError lx.next_token(io)
    end
    # some illegal examples strings
    for str in ("Ab[a]", "ABc[a]", "ab[cd]", "1A[bc]", "AB[]]")
        @testset "$str" begin
            @test_throws lx.LexicalError read_all_token(str)
        end
    end
end

@testset "Token equality" begin
    # make sure comparing Token works properly
    @test lx.Token('S', "test") == lx.Token('S', "test")
    @test lx.Token('S', "test") != lx.Token('S', "test1")
    @test lx.Token('S', "test") != lx.Token('S', "tesT")
    @test lx.Token('S', "test") != lx.Token('T', "test")
    @test hash(lx.Token('S', "test")) == hash(lx.Token('S', "test"))
    @test hash(lx.Token('S', "test")) != hash(lx.Token('S', "test1"))
    @test hash(lx.Token('S', "test")) != hash(lx.Token('S', "tesT"))
    @test hash(lx.Token('S', "test")) != hash(lx.Token('T', "test"))
end

@testset "premature EOF" begin
    # basically when reading an identifier or a property value
    for str in ("AB", "AB[oo", "AB[oo ", "AB[oo\n")
        @testset "$str" begin
            @test_throws EOFError read_all_token(str)
        end
    end
    # however, trailing white spaces will cause different behaviour.
    # this isn't ideal, but the parser should catch it either way
    for str in ("AB\t", "AB ", "AB\n")
        @testset "$str" begin
            @test read_all_token(str) == [lx.Token('I',"AB"), lx.Token('\0')]
        end
    end
end

@testset "normal examples" begin
    for (str, ref) in (
            # minimal legal file
            ("(;)", [lx.Token('('), lx.Token(';'), lx.Token(')')]),
            # will be a parse error, but lexically its fine
            ("()", [lx.Token('('), lx.Token(')')]),
            ("())", [lx.Token('('), lx.Token(')'), lx.Token(')')]),
            ("(()", [lx.Token('('), lx.Token('('), lx.Token(')')]),
            ("(A B)", [lx.Token('('), lx.Token('I',"A"), lx.Token('I',"B"), lx.Token(')')]),
            # typical property
            ("AB[oo]", [lx.Token('I',"AB"), lx.Token('['), lx.Token('S',"oo"), lx.Token(']')]),
            # check that white spaces are tolerated but skipped
            (" AB  \n[oo]",   [lx.Token('I',"AB"), lx.Token('['), lx.Token('S',"oo"), lx.Token(']')]),
            (" AB  \r[oo]",   [lx.Token('I',"AB"), lx.Token('['), lx.Token('S',"oo"), lx.Token(']')]),
            (" AB  \r\n[oo]", [lx.Token('I',"AB"), lx.Token('['), lx.Token('S',"oo"), lx.Token(']')]),
            (" AB  \n\r[oo]", [lx.Token('I',"AB"), lx.Token('['), lx.Token('S',"oo"), lx.Token(']')]),
            # in case of trailing white spaces we need an empty
            # token to avoid throwing EOFError. note that the
            # parser skips the empty token "Token('\0')" entirely.
            (" AB  \r\n[oo]\r\n", [lx.Token('I',"AB"), lx.Token('['), lx.Token('S',"oo"), lx.Token(']'), lx.Token('\0')]),
            # digits in identifiers are allowed for compatibility
            ("A1B[oo]", [lx.Token('I',"A1B"), lx.Token('['), lx.Token('S',"oo"), lx.Token(']')]),
            # a longer more complex example
            (" (; AB[ab] CD[cd])", [lx.Token('('), lx.Token(';'), lx.Token('I',"AB"), lx.Token('['), lx.Token('S',"ab"), lx.Token(']'), lx.Token('I',"CD"), lx.Token('['), lx.Token('S',"cd"), lx.Token(']'), lx.Token(')')]),
           )
        @testset "$str" begin
            @test read_all_token(str) == ref
        end
    end
end

@testset "property values examples" begin
    for (str, ref) in (
            # no property value will lack string token
            ("AB[]", [lx.Token('I',"AB"), lx.Token('['), lx.Token(']')]),
            # lexer does not attempt to interpret property values
            # the following just shows that composite values don't
            # cause any problems
            ("AB[aa:cc]", [lx.Token('I',"AB"), lx.Token('['), lx.Token('S',"aa:cc"), lx.Token(']')]),

            ## ESCAPING
            ("AB[ab(\\)[\\]\\;\\\\cd]", [lx.Token('I',"AB"), lx.Token('['), lx.Token('S',"ab()[];\\cd"), lx.Token(']')]),

            ## SPACES
            # spaces are preserved
            ("AB[ab  cd]", [lx.Token('I',"AB"), lx.Token('['), lx.Token('S',"ab  cd"), lx.Token(']')]),
            # spaces are preserved even when escaped
            ("AB[ab \\ cd]", [lx.Token('I',"AB"), lx.Token('['), lx.Token('S',"ab  cd"), lx.Token(']')]),

            ## TABS
            # tabs are converted to spaces
            ("AB[ab\tcd]", [lx.Token('I',"AB"), lx.Token('['), lx.Token('S',"ab cd"), lx.Token(']')]),
            # tabs are converted to spaces, even when escaped
            ("AB[ab\\\tcd]", [lx.Token('I',"AB"), lx.Token('['), lx.Token('S',"ab cd"), lx.Token(']')]),

            ## NEWLINES
            # new lines are preserved except when escaped by "\"
            ("AB[ab\nc\\\nd]", [lx.Token('I',"AB"), lx.Token('['), lx.Token('S',"ab\ncd"), lx.Token(']')]),
            # supported newlines are CR, LF, CRLF, LFCR
            ("AB[ab\rc\\\rd]", [lx.Token('I',"AB"), lx.Token('['), lx.Token('S',"ab\ncd"), lx.Token(']')]),
            ("AB[ab\n\rc\\\n\rd]", [lx.Token('I',"AB"), lx.Token('['), lx.Token('S',"ab\ncd"), lx.Token(']')]),
            ("AB[ab\r\nc\\\r\nd]", [lx.Token('I',"AB"), lx.Token('['), lx.Token('S',"ab\ncd"), lx.Token(']')]),
           )
        @testset "$str" begin
            @test read_all_token(str) == ref
        end
    end
end