text:
"""
Object with a *term* and an applicable *BasisFunction* and a eventfield that are later passed to the basisfunction.

$(TYPEDEF)
$(TYPEDSIGNATURES)
$(FIELDS)
# Examples
```julia-repl
julia>  b = TimeExpandedTerm(term,kernel,[:latencyTR,:durationTR])
```
"""
struct TimeExpandedTerm{T<:AbstractTerm} <: AbstractTerm
        "Term that the basis function is applied to. This is regularly called in other functions to get e.g. term-coefnames and timeexpand those"
        term::T
        "Kernel that determines what should happen to the designmatrix of the term"
        basisfunction::BasisFunction
        "Which fields of the event-table should be passed to the basisfunction.Important: The first entry has to be the event-latency in samples!"
        eventfields::Array{Symbol}
end

function TimeExpandedTerm(term::NTuple{N,Union{<:AbstractTerm,<:MixedModels.RandomEffectsTerm}},basisfunction,eventfields::Array{Symbol,1}) where {N}
        # for mixed models, apply it to each Term
        TimeExpandedTerm.(term,Ref(basisfunction),Ref(eventfields))
end
function TimeExpandedTerm(term,basisfunction,eventfields::Nothing)
        # if the eventfield is nothing, call default
        TimeExpandedTerm(term, basisfunction)
end

function TimeExpandedTerm(term,basisfunction,eventfields::Symbol)
        # call with symbol, convert to array of symbol
        TimeExpandedTerm(term, basisfunction,[eventfields])
end

function TimeExpandedTerm(term,basisfunction;eventfields=[:latency])
        # default call, use latency
        TimeExpandedTerm(term, basisfunction,eventfields)
end

function Base.show(io::IO, p::TimeExpandedTerm)
        #print(io, "timeexpand($(p.term), $(p.basisfunction.type),$(p.basisfunction.times))")
        println(io,"$(coefnames(p))")
end

"""
To be explored abstractTerm structure to automatically detect randomeffects. Not in use
"""
struct ZeroCorr2{T<:RandomEffectsTerm} <: AbstractTerm
    term::T
end

"""
Object with a *term* and an applicable *BasisFunction* and a eventfield that are later passed to the basisfunction.

$(TYPEDEF)
$(TYPEDSIGNATURES)
$(FIELDS)
# Examples
```julia-repl
julia>  b = TimeExpandedTerm(term,kernel,[:latencyTR,:durationTR])
```
"""
struct DesignMatrix
        "Array of formulas"
        formulas
        "A concatenated designmatric. In case of Mixed Models an array, where the first one is a FeMat, later ones ReMats. "
        Xs
        "Event table with all events"
        events
end

"""
$(SIGNATURES)
Combine two UnfoldDesignmatrices. This allows combination of multiple events.

This also allows to define events with different lengths.

Not supported for models without timebasis, as it is not needed there (one can simply run multiple models)
# Examples
```julia-repl
julia>  basisfunction1 = firbasis(τ=(0,1),sfreq = 10,name="basis1")
julia>  basisfunction2 = firbasis(τ=(0,0.5),sfreq = 10,name="basis2")
julia>  Xdc1          = designmatrix(UnfoldLinearModel,@formula 0~1,tbl_1,basisfunction1)
julia>  Xdc2          = designmatrix(UnfoldLinearModel,@formula 0~1,tbl_2,basisfunction2)
julia>  combineDesignmatrices(Xdc1,Xdc2)
julia>  Xdc = Xdc1+Xdc2 # equivalently
```

"""
function combineDesignmatrices(X1::DesignMatrix,X2::DesignMatrix)
        X1 = deepcopy(X1)
        X2 = deepcopy(X2)
        Xs1 = get_Xs(X1.Xs)
        Xs2 = get_Xs(X2.Xs)

        sX1 = size(Xs1,1)
        sX2 = size(Xs2,1)

        # append 0 to the shorter designmat
        if sX1 < sX2
                Xs1 = SparseMatrixCSC(sX2, Xs1.n, Xs1.colptr, Xs1.rowval, Xs1.nzval)
        elseif sX2 < sX1
                Xs2 = SparseMatrixCSC(sX1, Xs2.n, Xs2.colptr, Xs2.rowval, Xs2.nzval)
        end
        Xcomb = hcat(Xs1,Xs2)

        #if !(Float64(X1.formulas.rhs.basisfunction.times.step) ≈ Float64(X2.formulas.rhs.basisfunction.times.step))
        #        @warn("Concatenating formulas with different sampling rates. Be sure that this is what you want.")
        #end
        if typeof(X1.Xs) <: Tuple
                # we have random effects
                # combine REMats in single-eventtpe formulas ala y ~ (1|x) + (a|x)
                Xs1 = MixedModels._amalgamate([X1.Xs[2:end]...],Float64)
                Xs2 = MixedModels._amalgamate([X2.Xs[2:end]...],Float64)

                Xcomb = (Xcomb,Xs1...,Xs2...)

                # Next we make the ranefs all equal size
                equalizeReMatLengths!(Xcomb[2:end])

                # check if ranefs can be amalgamated. If this fails, then MixedModels tried to amalgamate over different eventtypes and we should throw the warning
                # if it success, we have to check if the size before and after is identical. If it is not, it tried to amalgamize over different eventtypes which were of the same length

                try
                reterms = MixedModels._amalgamate([Xcomb[2:end]...],Float64)
                if length(reterms) != length(Xcomb[2:end])
                        throw("IncompatibleRandomGroupings")
                end
                catch e
                        @error "Error, you seem to have two different eventtypes with the same random-effect grouping variable. \n
                        This is not allowed, you have to rename one. Example:\n
                        eventA: y~1+(1|item) \n
                        eventB: y~1+(1|item)  \n
                        This leads to this error. Rename the later one\n
                        eventB: y~1+(1|itemB) "
                        throw("IncompatibleRandomGroupings")
                end

        end

        if typeof(X1.formulas) <: FormulaTerm
                DesignMatrix([X1.formulas X2.formulas],Xcomb,[X1.events, X2.events])
        else
                DesignMatrix([X1.formulas... X2.formulas],Xcomb,[X1.events, X2.events])
        end
end

function changeMatSize!(m,fe,remats)
        changeReMatSize!.(remats,Ref(m))
        fe = SparseMatrixCSC(m, fe.n, fe.colptr,fe.rowval, fe.nzval)
        return (fe,remats...,)
end
function changeReMatSize!(remat::MixedModels.AbstractReMat,m ::Integer)

        n = m-length(remat.refs) # missing elements
        if n<0
                deleteat!(remat.refs,range(m+1,stop=length(remat.refs)))
                remat.adjA = remat.adjA[:,1:(m+1)]
                remat.wtz = remat.wtz[:,1:(m+1)]
                remat.z = remat.z[:,1:(m+1)]
        elseif n>0
                append!(remat.refs,repeat([remat.refs[end]],n))
                remat.adjA = [remat.adjA sparse(repeat([0],size(remat.adjA)[1],n))]
                remat.wtz = [remat.wtz zeros((size(remat.wtz)[1],n))]
                remat.z = [remat.z zeros((size(remat.z)[1],n))]
        end

        #ReMat{T,S}(remat.trm, refs, levels, cnames, z, wtz, λ, inds, adjA, scratch)
end
function equalizeReMatLengths!(remats::NTuple{A,MixedModels.AbstractReMat}) where {A}
        # find max length
        m = maximum([x[1] for x in size.(remats)])
        @debug print("combining lengths: $m")
        # for each reMat
        for k = range(1,length=length(remats))
                remat = remats[k]
                if size(remat)[1] == m
                        continue
                end
                # prolong if necessary
                changeReMatSize!(remat,m)

        end
end

Base.:+(X1::DesignMatrix, X2::DesignMatrix) = combineDesignmatrices(X1,X2)

function get_Xs(Xs::Tuple)
        return Xs[1]
end
function get_Xs(Xs::SparseMatrixCSC)
        return Xs
end

"""
$(SIGNATURES)
designmatrix(type, f, tbl; kwargs...)
Return a *DesignMatrix* used to fit the models.
# Arguments
- type::Union{UnfoldLinearMixedModel,UnfoldLinearModel}
- f::FormulaTerm: Formula to be used in this designmatrix
- tbl: Events (usually a data frame) to be modelled
- basisfunction::BasisFunction: basisfunction to be used in modeling (if specified)
- contrasts::Dict: (optional) contrast to be applied to formula
- eventfields::Array: (optional) Array of symbols which are passed to basisfunction event-wise.
First field of array always defines eventonset in samples. Default is [:latency]

# Examples
```julia-repl
julia>  designmatrix(UnfoldLinearModel,f,tbl,basisfunction1)
```

"""
function designmatrix(type,f,tbl,basisfunction;contrasts= Dict{Symbol,Any}(), kwargs...)
        @debug("generating DesignMatrix")
        form = apply_schema(f, schema(f, tbl, contrasts), MixedModels.LinearMixedModel)
        form = apply_basisfunction(form,basisfunction,get(Dict(kwargs),:eventfields,nothing))

        # Evaluate the designmatrix
        if (!isnothing(basisfunction)) & (type<:UnfoldLinearMixedModel)
                X = modelcols.(form.rhs, Ref(tbl))
        else
                X = modelcols(form.rhs, tbl)
        end
        return DesignMatrix(form,X,tbl)
end
function designmatrix(type,f,tbl;kwargs...)
        return designmatrix(type,f,tbl,nothing;kwargs...)
end

"""
$(SIGNATURES)
timeexpand the rhs-term of the formula with the basisfunction

"""

function apply_basisfunction(form,basisfunction::BasisFunction,eventfields)
        @debug("apply_basisfunction")
        return FormulaTerm(form.lhs, TimeExpandedTerm(form.rhs,basisfunction,eventfields))
end

function apply_basisfunction(form,basisfunction::Nothing,eventfields)
        # in case of no basisfunctin, do nothing
        return form
end

"""
$(SIGNATURES)
calculates the actual designmatrix for a timeexpandedterm. Multiple dispatch on StatsModels.modelcols
"""
# Timeexpand the fixed effect part
function StatsModels.modelcols(term::TimeExpandedTerm,tbl)

        X = modelcols(term.term,tbl)
        time_expand(X,term,tbl)
end

mutable struct SparseReMat{T,S} <: MixedModels.AbstractReMat{T}
        trm
        refs::Vector{Int32}
        levels
        cnames::Vector{String}
        z::SparseMatrixCSC{T}
        wtz::SparseMatrixCSC{T}
        λ::LowerTriangular{T,Matrix{T}}
        inds::Vector{Int}
        adjA::SparseMatrixCSC{T,Int32}
        scratch::Matrix{T}
    end
# This function timeexpands the random effects and generates a ReMat object
function StatsModels.modelcols(term::TimeExpandedTerm{<:RandomEffectsTerm},tbl)
#function StatsModels.modelcols(term::TimeExpandedTerm{<:Union{<:RandomEffectsTerm,<:AbstractTerm{<:RandomEffectsTerm}}},tbl)
# exchange this to get ZeroCorr to work

        # get the non-timeexpanded reMat
        reMat = modelcols(term.term,tbl)

        # Timeexpand the designmatrix
        z = transpose(time_expand(transpose(reMat.z),term,tbl))

        # First we check if there is overlap in the timeexpanded term. If so, we cannot continue. Later implementations will remedy that
        #println(dump(term,))
        if hasfield(typeof(term.term),:rhs)
                rhs = term.term.rhs

        elseif hasfield(typeof(term.term.term),:rhs)
                # we probably have something like zerocorr, which does not need to show a .rhs necessarily
                rhs = term.term.term.rhs
        else

                println("term.term: $(dump(term.term))")
                error("unknown RE structure, has no field .rhs:$(typeof(term.term))")

        end

        group = tbl[!,rhs.sym]
        time = tbl[!,term.eventfields[1]]

        # get the from-to onsets of the grouping varibales
        onsets = time_expand_getRandomGrouping(group,time,term.basisfunction)
        #print(size(reMat.z))
        refs = zeros(size(z)[2]).+1
        for (i,o) in enumerate(onsets[2:end])
                # check for overlap
                if (minimum(o) <= maximum(onsets[i+1])) & (maximum(o) <= minimum(onsets[i+1]))
                        error("overlap in random effects structure detected, not currently supported")
                end
        end

        # From now we can assume no overlap
        # We want to fnd which subject is active when
        refs = zeros(size(z)[2]).+1
        uGroup = unique(group)

        for (i,g) = enumerate(uGroup[1:end])

                ix_start = findfirst(g.==group)
                ix_end = findlast(g.==group)
                if i == 1
                        time_start = 1
                else
                        time_start = time[ix_start]
                        time_start = time_start - sum(term.basisfunction.times.<=0)
                end
                if i == length(uGroup)
                        time_stop = size(refs,1)
                else
                        time_stop = time[ix_end]
                        time_stop = time_stop + sum(term.basisfunction.times.>0)
                end
                if time_start < 0
                        time_start = 1
                end

                if time_stop > size(refs,1)
                        time_stop = size(refs,1)
                end

                #println("$g,$time_start,$time_stop")
                refs[Int64(time_start):Int64(time_stop)] .= g
        end

        # Other variables with implementaions taken from the LinerMixedModel function
        wtz = z
        trm = term

        S = size(z, 1)
        T = eltype(z)
        λ  = LowerTriangular(Matrix{T}(I, S, S))

        inds = MixedModels.sizehint!(Int[], (S * (S + 1)) >> 1)
        m = reshape(1:abs2(S), (S, S))
        inds = sizehint!(Int[], (S * (S + 1)) >> 1)
        for j in 1:S
                for i in j:S
                        # We currently restrict to diagonal entries
                        # Once mixedmodels#293 is pushed, we can relax this and use zerocorr()
                        if i == j # for diagonal
                                push!(inds, m[i, j])
                        end
                end
        end

        levels = reMat.levels
        refs =refs

        # reMat.levels doesnt change
        cnames = coefnames(term)
        #print(refs)
        adjA = MixedModels.adjA(refs, z)
        scratch = Matrix{T}(undef, (S, length(uGroup)))

        # Once MixedModels.jl supports it, can be replaced with: SparseReMat
        ReMat{T,S}(rhs,
        refs,
        levels,
        cnames,
        z,
        wtz,
        λ,
        inds,
        adjA,
        scratch)
end

"""
$(SIGNATURES)
Get the timeranges where the random grouping variable was applied
"""
function time_expand_getRandomGrouping(tblGroup,tblLatencies,basisfunction)
        ranges = time_expand_getTimeRange.(tblLatencies,Ref(basisfunction))
end

# helper function to get the ranges from where to where the basisfunction is added
function time_expand_getTimeRange(onset,basisfunction)
        npos = sum(basisfunction.times.>=0)
        nneg = sum(basisfunction.times.<0)

        basis = basisfunction.kernel(onset)

        fromRowIx = floor(onset)-nneg
        toRowIx = floor(onset)+npos

        range(fromRowIx,stop=toRowIx)
end

"""
$(SIGNATURES)
performs the actual time-expansion in a sparse way.

 - Get the non-timeexpanded designmatrix X from StatsModels.
 - evaluate the basisfunction kernel at each event
 - calculate the necessary rows, cols and values for the sparse matrix
 Returns SparseMatrixCSC
"""
function time_expand(X,term,tbl)
        ncolsBasis = size(term.basisfunction.kernel(0),2)
        X = reshape(X,size(X,1),:)

        ncolsX = size(X)[2]
        nrowsX = size(X)[1]
        ncolsXdc = ncolsBasis*ncolsX

        # this is the predefined eventfield, usually "latency"
        onsets = tbl[!,term.eventfields[1]]

        if typeof(term.eventfields) <:Array && length(term.eventfields) == 1
                bases = term.basisfunction.kernel.(tbl[!,term.eventfields[1]])
        else
                bases = term.basisfunction.kernel.(eachrow(tbl[!,term.eventfields]))
        end

        # generate rowindices
        rows =  copy(rowvals.(bases))
        # this shift is necessary as some basisfunction time-points can be negative. But a matrix is always from 1:τ. Thus we have to shift it backwards in time.
        # The onsets are onsets-1 XXX not sure why.
        for r in 1:length(rows)
                rows[r] .+= floor(onsets[r]-1)+term.basisfunction.shiftOnset
        end

        rows = vcat(rows...)
        rows = repeat(rows,ncolsX)

        # generate column indices
        cols = []
        for Xcol in 1:ncolsX
                for b in 1:length(bases)
                        for c in 1:ncolsBasis
                                push!(cols,repeat([c+(Xcol-1)*ncolsBasis],length(nzrange(bases[b],c))))
                        end
                end
        end
        cols = vcat(cols...)

        # generate values
        vals = []
        for Xcol in 1:ncolsX
                push!(vals,vcat(nonzeros.(bases).*X[:,Xcol]...))
        end

        vals = vcat(vals...)
        ix = rows.>0
        A = sparse(rows[ix],cols[ix],vals[ix])

        return A
end

"""
$(SIGNATURES)
coefnames of a TimeExpandedTerm concatenates the basis-function name with the kronecker product of the term name and the basis-function colnames. Separator is ' : '
Some examples for a firbasis:
        basis_313 : (Intercept) : 0.1
        basis_313 : (Intercept) : 0.2
        basis_313 : (Intercept) : 0.3
        ...
"""
function StatsModels.coefnames(term::TimeExpandedTerm)
        terms = coefnames(term.term)
        colnames = term.basisfunction.colnames
        name = term.basisfunction.name
        if typeof(terms) == String
                terms = [terms]
        end
        return name.*" : ".*kron(terms.*" : ",string.(colnames))
end

function StatsModels.coefnames(terms::AbstractArray{<:FormulaTerm})
        return coefnames.(Base.getproperty.(terms,:rhs))
end
function StatsModels.coefnames(term::MixedModels.ZeroCorr)
        coefnames(term.term)
end

function StatsModels.coefnames(term::RandomEffectsTerm)
        coefnames(term.lhs)
end


text:
using CompatHelper
using Pkg
using Test

registries_1 = Pkg.Types.RegistrySpec[Pkg.RegistrySpec(name = "General",
                                                       uuid = "23338594-aafe-5451-b93e-139f81909106",
                                                       url = "https://github.com/JuliaRegistries/General.git")]

registries_2 = Pkg.Types.RegistrySpec[Pkg.RegistrySpec(name = "General",
                                                       uuid = "23338594-aafe-5451-b93e-139f81909106",
                                                       url = "https://github.com/JuliaRegistries/General.git"),
                                      Pkg.RegistrySpec(name = "BioJuliaRegistry",
                                                       uuid = "ccbd2cc2-2954-11e9-1ccf-f3e7900901ca",
                                                       url = "https://github.com/BioJulia/BioJuliaRegistry.git")]

@testset "CompatHelper.jl" begin
    @testset "assert.jl" begin
        @test_nowarn CompatHelper.always_assert(true)
        @test CompatHelper.always_assert(true) isa Nothing
        @test CompatHelper.always_assert(true) == nothing
        @test @test_nowarn CompatHelper.always_assert(true) isa Nothing
        @test @test_nowarn CompatHelper.always_assert(true) == nothing
        @test_throws CompatHelper.AlwaysAssertionError CompatHelper.always_assert(false)
    end
    @testset "ci_service.jl" begin
        withenv("GITHUB_REPOSITORY" => "foo/bar") do
            @test CompatHelper.auto_detect_ci_service() isa CompatHelper.CIService
            @test CompatHelper.auto_detect_ci_service() isa CompatHelper.GitHubActions
            @test_throws ArgumentError CompatHelper.main(; keep_existing_compat = false, drop_existing_compat = false)
        end
        withenv("GITHUB_REPOSITORY" => nothing) do
            @test_throws ErrorException CompatHelper.auto_detect_ci_service()
        end
    end
    @testset "new_versions.jl" begin
        @test_throws ArgumentError CompatHelper.old_compat_to_new_compat("", "", :abc)
    end
end

ext: jl


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
Base.peek(io::Base.AbstractPipe) = Base.peek(Base.pipe_reader(io))

using BinaryBuilder

name = "LibicalBuilder"
version = v"3.0.0"

# Collection of sources required to build LibicalBuilder
sources = [
    "https://github.com/libical/libical/archive/3.0.zip" =>
    "c23d077479635b95e3819838d7b00e7a6c273a029f4ed497a90a575eed3f8c0c",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libical-3.0/
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DICAL_GLIB=OFF -DWITH_CXX_BINDINGS=OFF ..
make
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, :glibc),
    MacOS(:x86_64),
    FreeBSD(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libical", :libical)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

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

max_stars_repo_name: jgoldfar/LibicalBuilder

max_issues_repo_name: jgoldfar/LibicalBuilder

max_forks_repo_name: jgoldfar/LibicalBuilder


text:
using Interpolations

""" Model for one order of an echelle spectrum
Arrays organized into (basis_component,pixel).
"""
struct SpectralOrderModel{T1<:Real, T2<:Union{Real,Base.TwicePrecision{T1}}, T3<:Union{Real,Base.TwicePrecision{T1}},
      T4<:Real, T5<:Real, BlazeT<:AbstractBlaze1DType} <: AbstractSpectralOrderModel
   log_lambda::StepRangeLen{T1,T2,T3}
   star_basis::AbstractArray{T4,2}
   telluric_basis::AbstractArray{T5,2}
   blaze::BlazeT
   psf::AbstractPSF
   #=
   function SpectralOrderModel{T1,T2,T3,T4,T5}(log_lambda::LambdaT,
      star_basis::StarBasisT,
      telluric_basis::TelluricBasisT,
      blaze::BlazeT{T4,T5}
      )  where {    T1<:Real, T2<:Real, T3<:Real, T4<:Real, T5<:Real,
         LambdaT <: AbstractRange{T1}, StarBasisT <: AbstractArray{T2,2},
         TelluricBasisT <: AbstractArray{T3,2},
         #PixelRangeT <: AbstractRange{T4},
         BlazeT <: AbstractBlazeType
         }
      @assert length(log_lambda)==size(star_basis,2)==size(telluric_basis,2)
      @assert 1 <= size(star_basis,1) <= max_num_basis_vectors
      @assert 0 <= size(telluric_basis,1) <= max_num_basis_vectors
      new(log_lambda,star_basis,telluric_basis,blaze) # pixel_range,blaze)
   end
   =#
end

function SpectralOrderModel(data::ObsT, order::Integer;
   log_lambda::LambdaT, # TODO FIX = make_model_wavelength_range(data,order),
   star_basis::StarBasisT = default_stellar_template(length(log_lambda)), #num_pixels_per_order(data)),
   telluric_basis::TelluricBasisT = default_telluric_template(length(log_lambda)), #num_pixels_per_order(data)),
   pixel_range::PixelRangeT = 1:num_pixels_per_order(data),
   degree::Integer = default_blaze_degree,
   blaze::BlazeT = make_flat_blaze_polynomial(data,order,degree=degree,pixel_range=pixel_range),
   psf::PSFT = PSFDeltaFunction()
   )  where {    T1<:Real, T2<:Real, T3<:Real, T4<:Real,
      ObsT <: AbstractEchelleObservationOrSet,
      LambdaT <: AbstractRange{T1},
      StarBasisT <: AbstractArray{T2,2}, TelluricBasisT <: AbstractArray{T3,2},
      PixelRangeT <: AbstractRange{T4}, BlazeT <: AbstractBlazeType,
      PSFT <: AbstractPSF
      }
   # TODO: move asserts to inner constructor?
   @assert 1 <= order <= num_orders(data)
   @assert 1 <= length(log_lambda) <= max_template_length
   @assert size(star_basis,2) == size(telluric_basis,2)
   @assert 1 <= size(star_basis,2)
   @assert min_num_star_basis_vectors <= size(star_basis,1) <= max_num_star_basis_vectors
   @assert min_num_telluric_basis_vectors <= size(telluric_basis,1) <= max_num_telluric_basis_vectors
   @assert blaze.order == order

   SpectralOrderModel(log_lambda,star_basis,telluric_basis,blaze,psf)
end

function get_order(som::SOMT) where { SOMT <: AbstractSpectralOrderModel }
   som.blaze.order
end

""" Model for one order of an echelle spectrum
Arrays organized into (basis_component,pixel).
"""
struct SpectralModel{T1<:Real, T2<:Real, T3<:Real, BlazeT<: AbstractBlaze2DType } <: AbstractSpectralModel
   log_lambda::StepRangeLen{T1,Base.TwicePrecision{T1},Base.TwicePrecision{T1} }
   star_basis::AbstractArray{T2,2}
   telluric_basis::AbstractArray{T3,2}
   blaze_list::BlazeT
   psf::AbstractPSF
end

function SpectralModel(data::ObsT;
   log_lambda::LambdaT = make_model_wavelength_range(data),
   star_basis::StarBasisT = default_stellar_template(length(log_lambda)), # num_orders(data)*num_pixels_per_order(data)),
   telluric_basis::TelluricBasisT = default_telluric_template(length(log_lambda)), #num_orders(data)*num_pixels_per_order(data)),
   pixel_range::PixelRangeT = 1:num_pixels_per_order(data),
   degree::Integer = default_blaze_degree,
   blaze_list::BlazeT = make_flat_blaze_polynomial_2d(data,pixel_rage=pixel_range,degree=degree),
      #map(order->make_flat_blaze(data,order,degree=degree,pixel_range=pixel_range),1:num_orders(data)),
   psf::PSFT = PSFDeltaFunction()
   )  where {    T1<:Real, T2<:Real, T3<:Real, T4<:Real,
      ObsT <: AbstractEchelleObservationOrSet,
      LambdaT <: AbstractRange{T1},
      StarBasisT <: AbstractArray{T2,2}, TelluricBasisT <: AbstractArray{T3,2},
      PixelRangeT <: AbstractRange{T4}, BlazeT <: AbstractBlaze2DType, # BT<: AbstractBlazeType, BlazeT <: AbstractArray{BT,1},
      PSFT <: AbstractPSF
      }
   # TODO: Move asserts to inner constructor?
   @assert 1 <= length(log_lambda) <= max_template_length
   @assert length(log_lambda)==size(star_basis,2) == size(telluric_basis,2)
   @assert 1 <= size(star_basis,1)
   @assert min_num_star_basis_vectors <= size(star_basis,1) <= max_num_star_basis_vectors
   @assert min_num_telluric_basis_vectors <= size(telluric_basis,1) <= max_num_telluric_basis_vectors
   @assert all( 1 .<= map(b->b.order,blaze_list.blaze_list) .<= num_orders(data) )
   @assert length(blaze_list.blaze_list) == num_orders(data)  # TODO: Could relax, but would need to remove assumption index=blaze.order

   SpectralModel(log_lambda,star_basis,telluric_basis,blaze_list,psf)
end

function SpectralOrderModel(model::SMT, order::Integer; max_abs_z::Real = 1e-4 ) where { SMT <: AbstractSpectralModel}
   @assert 1 <= order <= length(model.blaze_list.blaze_list)
   log_lambda_lo = min_log_lambda(model.blaze_list.blaze_list[order]) - log1p(max_abs_z)
   log_lambda_hi = max_log_lambda(model.blaze_list.blaze_list[order]) + log1p(max_abs_z)
   # TODO OPT: Above use model.log_lambda's step to figure out how much to move idx without taking log's
   idx_lo = log_lambda_lo > model.log_lambda[1] ? findfirst(l->l>=log_lambda_lo, model.log_lambda) : 1
   idx_hi = log_lambda_hi < model.log_lambda[end] ? findlast(l->l<=log_lambda_hi, model.log_lambda) : length(model.log_lambda)
   len = idx_hi-idx_lo+1
   #log_lambda_narrow = StepRangeLen{Float64,Float64,Base.TwicePrecision{Float64}}(model.log_lambda[1],model.log_lambda.step,len,idx_lo-1)
   log_lambda_narrow = StepRangeLen{Float64,Float64,Base.TwicePrecision{Float64}}(model.log_lambda[idx_lo],model.log_lambda.step,len)
   #log_lambda_narrow = range(model.log_lambda[idx_lo],step=model.log_lambda.step,length=len)
   #log_lambda_narrow = range(model.log_lambda[idx_lo],model.log_lambda.step,len)
   #log_lambda_narrow = StepRangeLen(model.log_lambda[1],model.log_lambda.step,len,idx_lo-1)
   SpectralOrderModel(log_lambda_narrow,view(model.star_basis,:,idx_lo:idx_hi),view(model.telluric_basis,:,idx_lo:idx_hi),model.blaze_list.blaze_list[order],model.psf)
end

# Default templates for star and tellurics
""" default_telluric_template(n; num_basis_vectors = 0)
   Returns array with (num_basis_vectors, n) zeros. """
function default_telluric_template(n::Integer; num_basis_vectors::Integer = 0)
   @assert 0 <= min_num_telluric_basis_vectors <= num_basis_vectors <= max_num_telluric_basis_vectors
   @assert 1 <= n <= max_template_length
   zeros(num_basis_vectors,n)
end

""" default_stellar_template(n; num_basis_vectors = 0)
   Returns array with one row of n 1's followed by num_basis_vectors rows of zeros. """
function default_stellar_template(n::Integer; num_basis_vectors::Integer = 0)
   @assert 0 <= min_num_star_basis_vectors-1 <= num_basis_vectors <= max_num_star_basis_vectors
   @assert 1 <= n <= max_template_length
   result = vcat(ones(1,n),zeros(num_basis_vectors,n))
   return result
end

function make_model_wavelength_range(data::EOT, order::OrdersType;
            oversample_factor::Real = 1.0, max_abs_z::Real = 1e-4,
            num_knots::Integer=round(Int64,oversample_factor*(1+max_abs_z)^2*num_pixels_per_order(data))
            ) where { EOT <: AbstractEchelleObservation }
      min_log_lambda = log(minimum(data.lambda[1,order])  / (1+max_abs_z) )
      max_log_lambda = log(maximum(data.lambda[end,order]) * (1+max_abs_z) )
      range(min_log_lambda, stop=max_log_lambda, length=num_knots )
end

function make_model_wavelength_range(data::EOST, order::OrdersType;
            oversample_factor::Real = 1.0, max_abs_z::Real = 1e-4,
            num_knots::Integer=round(Int64,oversample_factor*(1+max_abs_z)^2*num_pixels_per_order(data))
            ) where { EOST <: AbstractEchelleObservationSet }
      min_log_lambda = log(minimum(data.lambda[1,order])  / (1+max_abs_z) )
      max_log_lambda = log(maximum(data.lambda[end,order]) * (1+max_abs_z) )
      range(min_log_lambda, stop=max_log_lambda, length=num_knots )
end

function make_model_wavelength_range(data::EOST;
            oversample_factor::Real = 1.0, max_abs_z::Real = 1e-4,
            num_knots::Integer=round(Int64,oversample_factor*(1+max_abs_z)^2*num_orders(data)*num_pixels_per_order(data))
            ) where { EOST <: AbstractEchelleObservationSet }
      min_log_lambda = log(minimum(data.lambda[1,1,:])      / (1+max_abs_z) )
      max_log_lambda = log(maximum(data.lambda[end,end,:])  * (1+max_abs_z) )
      range(min_log_lambda, stop = max_log_lambda, length = num_knots )
end

""" Model for one echelle observation (potentially multiple orders)
Arrays over (basis_component).
"""
struct ObservationModel{T1<:Real, T2<:Real, T3<:Real, T4<:Real } <: AbstractObservationModel
   star_score::Array{T1,1}
   telluric_score::Array{T2,1}
   z_bc::T3
   normalization::T4

   function ObservationModel{T1,T2,T3,T4}( star_score::AbstractArray{T1,1},
      telluric_score::AbstractArray{T2,1}, z_bc::T3, normalization::T4
      )  where {T1<:Real, T2<:Real, T3<:Real, T4<:Real }
      @assert min_num_star_basis_vectors <= size(star_score,1) <= max_num_star_basis_vectors
      @assert min_num_telluric_basis_vectors <= size(telluric_score,1) <= max_num_telluric_basis_vectors
      new(star_score,telluric_score,z_bc,normalization)
   end

end
export ObservationModel

function ObservationModel(;
   num_star_basis_vectors::Integer = 1,
   num_telluric_basis_vector::Integer = 0,
   star_score::AbstractArray{T1,1} = zeros(num_star_basis_vectors),
   telluric_score::AbstractArray{T2,1} = zeros(num_telluric_basis_vector),
   z_bc::T3 = 0.0,
   normalization::T4 = 1.0
   )  where {T1<:Real, T2<:Real, T3<:Real, T4<:Real }
   ObservationModel{eltype(star_score),eltype(telluric_score),typeof(z_bc),typeof(normalization)}(star_score,telluric_score,z_bc,normalization)
end

""" Model for a time series of echelle observations (potentially multiple orders)
Arrays of scores over (basis_component,observation).
Arrays of z and normalization over (observation).
"""
struct ObsSeriesModel{T1<:Real, T2<:Real, T3<:Real, T4<:Real } <: AbstractObsSeriesModel
   star_score::Array{T1,2}
   telluric_score::Array{T2,2}
   z_bc::Array{T3,1}
   normalization::Array{T4,1}
   function ObsSeriesModel{T1,T2,T3,T4}( star_score::AbstractArray{T1,2},
      telluric_score::AbstractArray{T2,2}, z_bc::AbstractArray{T3,1}, normalization::AbstractArray{T4,1}
      )  where {T1<:Real, T2<:Real, T3<:Real, T4<:Real }
      @assert min_num_star_basis_vectors <= size(star_score,1) <= max_num_star_basis_vectors
      @assert min_num_telluric_basis_vectors <= size(telluric_score,1) <= max_num_telluric_basis_vectors
      @assert size(star_score,2) == size(telluric_score,2) == length(z_bc) == length(normalization)
      @assert length(z_bc) <= max_num_obs
      new(star_score,telluric_score,z_bc,normalization)
   end
end
export ObsSeriesModel

function ObsSeriesModel(num_obs::Integer;
   num_star_basis_vectors::Integer = 1,
   num_telluric_basis_vector::Integer = 0,
   star_score::AbstractArray{T1,2} = vcat(ones(1,num_obs),zeros(num_star_basis_vectors-1,num_obs)),
   telluric_score::AbstractArray{T2,2} = zeros(num_telluric_basis_vector,num_obs),
   z_bc::AbstractArray{T3,1} = zeros(num_obs),
   normalization::AbstractArray{T4,1} = ones(num_obs)
   )  where {T1<:Real, T2<:Real, T3<:Real, T4<:Real }
   @assert min_num_star_basis_vectors <= size(star_score,1) <= max_num_star_basis_vectors
   @assert min_num_telluric_basis_vectors <= size(telluric_score,1) <= max_num_telluric_basis_vectors
   @assert size(star_score,2) == size(telluric_score,2) == length(z_bc) == length(normalization)

   ObsSeriesModel{eltype(star_score),eltype(telluric_score),eltype(z_bc),eltype(normalization)}(star_score,telluric_score,z_bc,normalization)
end

function ObservationModel(obs_models::ObsSeriesModel, i::Integer) where {ObsSeriesModel <: AbstractObsSeriesModel}
   @assert 1 <= i <= length(obs_models.normalization)
   ObservationModel( star_score=view(obs_models.star_score,:,i), telluric_score=view(obs_models.telluric_score,:,i), z_bc=obs_models.z_bc[i], normalization=obs_models.normalization[i] )
end

function num_obs(osm::OSMT) where { OSMT <: AbstractObsSeriesModel }
  return length(osm.normalization)
end

#=
function eval_model_baseline(log_lambda_obs::LambdaT,
      spec_model::SOMT;
      pixel_range_obs::AbstractRange = 1:length(log_lambda_obs),
      output::OutputT = zeros(eltype(spec_model.star_basis),length(log_lambda_obs))
      ) where {  T1 <: Real, T2 <: Real,
      LambdaT <: AbstractArray{T1,1} #= Wavelength1DType{T1} =#, SOMT <: AbstractSpectralOrderModel,
      OMT <: AbstractObservationModel, OutputT <: AbstractVector{T2} }
   @assert length(log_lambda_obs) == length(pixel_range_obs)
   @assert length(log_lambda_obs) == length(output)
   # Compute each term on a grid
   blaze = calc_blaze(spec_model.blaze,pixel_range_obs) # TODO: Use sub-pixels once integrating over PSF
   star_on_grid = spec_model.star_basis[1,:] # * obs.star_score
   interp_star = scale(interpolate( star_on_grid, BSpline(Cubic(Line(OnGrid())))), spec_model.log_lambda)
   # Interpolate each term to observed wavelength
   flux_star = interp_star.(log_lambda_obs)
   telluric_factor = 1.0
   # TODO: Replace product of interpolations with integration over PSF
   output .= obs.normalization .* blaze .* telluric_factor .* flux_star
end
=#

function eval_model_star(spec_model::SOMT, obs_model::OMT;
      log_lambda_obs::LambdaT = spec_model.log_lambda, # view(spec_model.log_lambda,2:(length(spec_model.log_lambda)-1)),
      output::OutputT = zeros(eltype(spec_model.star_basis),length(log_lambda_obs))
      ) where {  T1 <: Real, T2 <: Real,
      SOMT <: AbstractSpectralOrderModel, OMT <: AbstractObservationModel,
      LambdaT <: AbstractArray{T1,1} #=Wavelength1DType{T1}=#, #= PixelRangeT <: AbstractRange, =#
      OutputT <: AbstractVector{T2} }
   @assert length(log_lambda_obs) == length(output)
   @assert size(spec_model.star_basis,1) == size(obs_model.star_score,1)
   @assert size(spec_model.telluric_basis,1) == size(obs_model.telluric_score,1)
   # Compute each term on a grid
   star_on_grid = spec_model.star_basis' * obs_model.star_score
   # TODO: Replace with my interpolator to avoid having to allocate memory to store values at knots
   interp_star = scale(extrapolate(interpolate( star_on_grid, BSpline(Cubic(Line(OnGrid())))),Line()), spec_model.log_lambda)
   # Interpolate each term to observed wavelength
   log_lambda_source = log_lambda_ssb_from_obs.(log_lambda_obs, -obs_model.z_bc)
   flux_star = interp_star.(log_lambda_source)
end

function eval_model_telluric(spec_model::SOMT, obs_model::OMT;
      log_lambda_obs::LambdaT = spec_model.log_lambda,
      output::OutputT = zeros(eltype(spec_model.star_basis),length(log_lambda_obs))
      ) where {  T1 <: Real, T2 <: Real,
      SOMT <: AbstractSpectralOrderModel, OMT <: AbstractObservationModel,
      LambdaT <: AbstractArray{T1,1} #=Wavelength1DType{T1}=#, #= PixelRangeT <: AbstractRange, =#
      OutputT <: AbstractVector{T2} }
   @assert length(log_lambda_obs) == length(output)
   @assert size(spec_model.star_basis,1) == size(obs_model.star_score,1)
   # Compute each term on a grid
   telluric_on_grid = 1.0 .- spec_model.telluric_basis' * obs_model.telluric_score
   # TODO: Replace with my interpolator to avoid having to allocate memory to store values at knots
   interp_telluric = scale(interpolate( telluric_on_grid, BSpline(Cubic(Line(OnGrid())))), spec_model.log_lambda)
   # Interpolate each term to observed wavelength
   telluric_factor = interp_telluric.(log_lambda_obs)
end

function eval_model_star_telluric(spec_model::SOMT, obs_model::OMT;
      log_lambda_obs::LambdaT = spec_model.log_lambda,
      output::OutputT = zeros(eltype(spec_model.star_basis),length(log_lambda_obs))
      ) where {  T1 <: Real, T2 <: Real,
      SOMT <: AbstractSpectralOrderModel, OMT <: AbstractObservationModel,
      LambdaT <: AbstractArray{T1,1} #=Wavelength1DType{T1}=#, OutputT <: AbstractVector{T2} }
   @assert length(log_lambda_obs) == length(output)
   @assert size(spec_model.telluric_basis,1) == size(obs_model.telluric_score,1)
   output .= eval_model_star(spec_model,obs_model,log_lambda_obs=log_lambda_obs) .*
      eval_model_telluric(spec_model,obs_model,log_lambda_obs=log_lambda_obs) .*
      obs_model.normalization
end

function eval_model_blaze(spec_model::SOMT #=, obs_model::OMT=#;
      pixel_range_obs::PixelRangeT = spec_model.blaze.pixel_range,
      output::OutputT = zeros(eltype(spec_model.star_basis),length(pixel_range_obs))
      ) where {  #= T1 <: Real, =# T2 <: Real,
      SOMT <: AbstractSpectralOrderModel, OMT <: AbstractObservationModel,
      PixelRangeT <: AbstractRange, OutputT <: AbstractVector{T2} }
   @assert length(pixel_range_obs) == length(output)
   output .= calc_blaze(spec_model.blaze,pixel_range_obs)
end

function eval_model(spec_model::SOMT, obs_model::OMT;
      #pixel_range_obs::PixelRangeT = spec_model.blaze.pixel_range,
      log_lambda_obs::LambdaT = spec_model.log_lambda,
      output::OutputT = zeros(eltype(spec_model.star_basis),length(log_lambda_obs))
      ) where { T1 <: Real, T2 <: Real,
      SOMT <: AbstractSpectralOrderModel, OMT <: AbstractObservationModel,
      LambdaT <: AbstractArray{T1,1} #=Wavelength1DType{T1}=# ,  #=PixelRangeT <: AbstractRange, =#
      OutputT <: AbstractVector{T2} }
   @assert length(output) == length(log_lambda_obs)
   # @assert length(log_lambda_obs) == length(spec_model.blaze.pixel_range) # ERROR TODO: UNDERSTAND WHY?
   #blaze = eval_model_blaze(spec_model)
   blaze = calc_blaze_at_log_lambda(spec_model.blaze, log_lambda_obs)
   incoming_flux = eval_model_star_telluric(spec_model,obs_model,log_lambda_obs=log_lambda_obs)
   output .= blaze .* incoming_flux
end

function eval_model(spec_model::SMT, obs_model::OMT, log_lambda_obs::LambdaT;
      #= pixel_range_obs::PixelRangeTA = fill(1:size(log_lambda_obs,1),size(log_lambda_obs,2)), =#
      output::OutputT = zeros(eltype(spec_model.star_basis),size(log_lambda_obs))
      ) where {  T1 <: Real, #= T2 <: Real, =# T3 <: Real,
      LambdaT <: AbstractArray{T1,2}, SMT <: AbstractSpectralModel,
      OMT <: AbstractObservationModel, #= PixelRangeT <: AbstractRange{T2},
      PixelRangeTA <: AbstractArray{PixelRangeT,1}, =#
      OutputT <: AbstractArray{T3,2} }
   @assert size(output,2) == size(log_lambda_obs,2)
   map(i->eval_model(SpectralOrderModel(spec_model,i),obs_model,
               log_lambda_obs=view(log_lambda_obs,:,i),
               #= pixel_range_obs=pixel_range_obs[i], =#
               output=view(output,:,i)), 1:size(log_lambda_obs,2) )
   return output
end

function eval_model(spec_model::SMT, obs_series_model::OSMT, log_lambda_obs::LambdaT;
      pixel_range_obs::PixelRangeTA = fill(1:size(log_lambda_obs,1),size(log_lambda_obs,2)),
      output::OutputT = zeros(eltype(spec_model.star_basis),size(log_lambda_obs))
      ) where {  T1 <: Real, T2 <: Real, T3 <: Real,
      LambdaT <: AbstractArray{T1,3}, SMT <: AbstractSpectralModel,
      OSMT <: AbstractObsSeriesModel, PixelRangeT <: AbstractRange{T2},
      PixelRangeTA <: AbstractArray{PixelRangeT,1},
      OutputT <: AbstractArray{T3,3} }
   @assert size(output,3) == size(log_lambda_obs,3)
   map(obsid -> eval_model(spec_model,ObservationModel(obs_series_model,obsid),
                           view(log_lambda_obs,:,:,obsid),
                           #pixel_range_obs=pixel_range_obs,
                           output=view(output,:,:,obsid)),
               1:num_obs(obs_series_model) )

   return output
end

#=
function eval_model(log_lambda_obs::LambdaT,
      spec_model::AbstractArray{SOMT,1}, obs_model::OMT;
      pixel_range_obs::AbstractArray{AbstractRange,1} = fill(1:size(log_lambda_obs,1),size(log_lambda_obs,2)),
      output::OutputT = zeros(eltype(model.star_basis),size(log_lambda_obs,1),size(log_lambda_obs,2))
      ) where {  T1 <: Real, T2 <: Real,
      LambdaT <: AbstractArray{AbstractRange{T1},1}, SOMT <: AbstractSpectralOrderModel,
      OMT <: AbstractObservationModel, OutputT <: AbstractArray{T2,2} }
   map(i->eval_model(log_lambda_obs[i],spec_model[i],obs_model,
         pixel_range=pixel_range_obs[i],output=view(output,:,i)), 1:size(log_lambda_obs,2) )
   return output
end
=#

function make_stellar_template(data::EOST, sm::SMT, obsm::OSMT
   ) where { EOST <: AbstractEchelleObservationSet, SMT <: AbstractSpectralModel, OSMT <: AbstractObsSeriesModel }
   model = eval_model(log.(data.lambda), sm, obsm)
   template = reshape(mean(data.flux ./ model, dims=3), num_pixels_per_order(data), num_orders(data) )
end

function eval_model_fit(data::EOST, spec_model::SMT, obs_model::OSMT
      #;output::OutputT = zeros(eltype(spec_model.star_basis),size(data.lambda))
      ) where {  T2<:Real, T3 <: Real,
      EOST <: AbstractEchelleObservationSet, SMT <: AbstractSpectralModel,
      OSMT <: AbstractObsSeriesModel, PixelRangeT <: AbstractRange{T2},
      PixelRangeTA <: AbstractArray{PixelRangeT,1},
      OutputT <: AbstractArray{T3,2} }
   sum_sq = 0.
   for obsid in 1:num_obs(data)
      obsm = ObservationModel(obs_model,obsid)
      for i in 1:num_orders(data)
         pixel_range = spec_model.blaze_list.blaze_list[i].pixel_range
         sum_sq += sum( ( data.flux[pixel_range] .-
                  eval_model(log.(view(data.lambda,pixel_range,i,obsid)),
                     SpectralOrderModel(spec_model,i),obsm,pixel_range=pixel_range ) ).^2 ./
                        data.var[pixel_range,i,obsid] )
      end
   end
   return sum_sq
end

#=
function eval_model_fit(data::EOST, sm::SMT, obsm::OSMT
   ) where { EOST <: AbstractEchelleObservationSet, SMT <: AbstractSpectralModel, OSMT <: AbstractObsSeriesModel }
   model = eval_model(log.(data.lambda), sm, obsm)
   data.flux .- model
   template = reshape(mean(data.flux ./ model, dims=3), num_pixels_per_order(data), num_orders(data) )
end
=#


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""
    Generator

struct with inner constructor:
```julia
function Generator(;
    existing_kw::Real=0,
    min_kw::Real=0,
    max_kw::Real=1.0e6,
    installed_cost_per_kw::Real=500.0,
    om_cost_per_kw::Real=10.0,
    om_cost_per_kwh::Float64=0.0,
    fuel_cost_per_gallon::Float64 = 3.0,
    fuel_slope_gal_per_kwh::Float64 = 0.076,
    fuel_intercept_gal_per_hr::Float64 = 0.0,
    fuel_avail_gal::Float64 = 660.0,
    min_turn_down_pct::Float64 = 0.0,
    only_runs_during_grid_outage::Bool = true,
    sells_energy_back_to_grid::Bool = false,
    can_net_meter::Bool = false,
    can_wholesale::Bool = false,
    can_export_beyond_nem_limit = false,
    macrs_option_years::Int = 0,
    macrs_bonus_pct::Float64 = 1.0,
    macrs_itc_reduction::Float64 = 0.0,
    federal_itc_pct::Float64 = 0.0,
    federal_rebate_per_kw::Float64 = 0.0,
    state_ibi_pct::Float64 = 0.0,
    state_ibi_max::Float64 = 1.0e10,
    state_rebate_per_kw::Float64 = 0.0,
    state_rebate_max::Float64 = 1.0e10,
    utility_ibi_pct::Float64 = 0.0,
    utility_ibi_max::Float64 = 1.0e10,
    utility_rebate_per_kw::Float64 = 0.0,
    utility_rebate_max::Float64 = 1.0e10,
    production_incentive_per_kwh::Float64 = 0.0,
    production_incentive_max_benefit::Float64 = 1.0e9,
    production_incentive_years::Int = 0,
    production_incentive_max_kw::Float64 = 1.0e9,
)
```
"""
struct Generator <: AbstractGenerator
    existing_kw
    min_kw
    max_kw
    installed_cost_per_kw
    om_cost_per_kw
    om_cost_per_kwh
    fuel_cost_per_gallon
    fuel_slope_gal_per_kwh
    fuel_intercept_gal_per_hr
    fuel_avail_gal
    min_turn_down_pct
    only_runs_during_grid_outage
    sells_energy_back_to_grid
    can_net_meter
    can_wholesale
    can_export_beyond_nem_limit
    can_curtail
    macrs_option_years
    macrs_bonus_pct
    macrs_itc_reduction
    federal_itc_pct
    federal_rebate_per_kw
    state_ibi_pct
    state_ibi_max
    state_rebate_per_kw
    state_rebate_max
    utility_ibi_pct
    utility_ibi_max
    utility_rebate_per_kw
    utility_rebate_max
    production_incentive_per_kwh
    production_incentive_max_benefit
    production_incentive_years
    production_incentive_max_kw

    function Generator(;
        existing_kw::Real=0,
        min_kw::Real=0,
        max_kw::Real=1.0e6,
        installed_cost_per_kw::Real=500.0,
        om_cost_per_kw::Real=10.0,
        om_cost_per_kwh::Float64=0.0,
        fuel_cost_per_gallon::Float64 = 3.0,
        fuel_slope_gal_per_kwh::Float64 = 0.076,
        fuel_intercept_gal_per_hr::Float64 = 0.0,
        fuel_avail_gal::Float64 = 660.0,
        min_turn_down_pct::Float64 = 0.0,
        only_runs_during_grid_outage::Bool = true,
        sells_energy_back_to_grid::Bool = false,
        can_net_meter::Bool = false,
        can_wholesale::Bool = false,
        can_export_beyond_nem_limit = false,
        can_curtail::Bool = false,
        macrs_option_years::Int = 0,
        macrs_bonus_pct::Float64 = 1.0,
        macrs_itc_reduction::Float64 = 0.0,
        federal_itc_pct::Float64 = 0.0,
        federal_rebate_per_kw::Float64 = 0.0,
        state_ibi_pct::Float64 = 0.0,
        state_ibi_max::Float64 = 1.0e10,
        state_rebate_per_kw::Float64 = 0.0,
        state_rebate_max::Float64 = 1.0e10,
        utility_ibi_pct::Float64 = 0.0,
        utility_ibi_max::Float64 = 1.0e10,
        utility_rebate_per_kw::Float64 = 0.0,
        utility_rebate_max::Float64 = 1.0e10,
        production_incentive_per_kwh::Float64 = 0.0,
        production_incentive_max_benefit::Float64 = 1.0e9,
        production_incentive_years::Int = 0,
        production_incentive_max_kw::Float64 = 1.0e9,
        )

        new(
            existing_kw,
            min_kw,
            max_kw,
            installed_cost_per_kw,
            om_cost_per_kw,
            om_cost_per_kwh,
            fuel_cost_per_gallon,
            fuel_slope_gal_per_kwh,
            fuel_intercept_gal_per_hr,
            fuel_avail_gal,
            min_turn_down_pct,
            only_runs_during_grid_outage,
            sells_energy_back_to_grid,
            can_net_meter,
            can_wholesale,
            can_export_beyond_nem_limit,
            can_curtail,
            macrs_option_years,
            macrs_bonus_pct,
            macrs_itc_reduction,
            federal_itc_pct,
            federal_rebate_per_kw,
            state_ibi_pct,
            state_ibi_max,
            state_rebate_per_kw,
            state_rebate_max,
            utility_ibi_pct,
            utility_ibi_max,
            utility_rebate_per_kw,
            utility_rebate_max,
            production_incentive_per_kwh,
            production_incentive_max_benefit,
            production_incentive_years,
            production_incentive_max_kw
        )
    end
end


text:
struct NoDefault end

struct Verifier
    type
    f
    default
    optional::Bool
    forfield::String
end
struct VerifierTree
    dict::Dict
    optional::Bool
    forfield::String
end
VerifierTree(dict::Dict; optional=false, forfield="") = VerifierTree(dict, optional, forfield)
VerifierTree(pairs::Pair...; optional=false, forfield="") = VerifierTree(Dict(pairs...), optional, forfield)

struct VerificationError <: Exception
    msg::String
end
function Verifier(::Type{T}, f=@λ x->true; default=NoDefault(), optional=false, forfield="") where T
    if default !== NoDefault()
        if !(default isa T)
            error("The type of default value `$default` does not match the required type `$T`!")
        elseif !f(default)
            error("The default value does not pass the checker!")
        end
    end
    Verifier(T, f, default, optional, forfield)
end

# Better printing
using AbstractTrees

function AbstractTrees.children(ne::VerifierTree)
    [k=>v for (k, v) in ne.dict]
end
AbstractTrees.children(ne::Pair{String, <:VerifierTree}) = children(ne.second)
AbstractTrees.children(ne::Pair{String, <:Verifier}) = []
function AbstractTrees.printnode(io::IO, x::VerifierTree)
    print(io, "ROOT")
    x.optional && print(io, " optional")
end
function AbstractTrees.printnode(io::IO, x::Pair{String, <:VerifierTree})
    print(io, "$(x.first)")
    x.second.optional && print(io, ", optional")
end
function AbstractTrees.printnode(io::IO, x::Pair{String, <:Verifier})
    e = x.second
    print(io, "$(x.first)::$(e.type)")
    print(io, ", checker = ")
    show(io, e.f)
    e.optional && print(io, ", optional")
    e.default isa NoDefault || print(io, ", default = $(e.default)")
end
function Base.show(io::IO, e::VerifierTree)
    print_tree(io, e)
end
Base.show(io::IO, ::MIME"text/plain", e::VerifierTree) = Base.show(io, e)

function verify(tree::VerifierTree, dict)
    updated = Dict()
    # sort the fields for that the dependent fields are updated the last
    kvpairs = [(k,v) for (k, v) in tree.dict]
    sort!(kvpairs, by=x->!isempty(x[2].forfield))
    for (k, v) in kvpairs
        if !isempty(v.forfield)  # check target field
            method = updated[v.forfield]
            if "$(v.forfield).$(method)" != k   # we do not need to extract this conditional field!
                continue
            end
        end
        if v isa VerifierTree
            # recurse
            if !haskey(dict, k)
                if !v.optional
                    throw(VerificationError("Required argument absent: $k"))
                else
                    updated[k] = verify(v, Dict())
                end
            else
                updated[k] = verify(v, dict[k])
            end
        else
            # verify
            if !haskey(dict, k)
                if v.default === NoDefault()
                    if !v.optional
                        throw(VerificationError("Required argument absent: $k"))
                    else
                        # leave it empty
                    end
                else
                    updated[k] = v.default
                end
            else
                val = dict[k]
                # type check
                if !(val isa v.type)
                    try
                        val = v.type(val)   # type cast
                    catch
                        throw(VerificationError("Input argument type should be $(v.type), got: $(typeof(val))"))
                    end
                end
                # f check
                if !v.f(val)
                    throw(VerificationError("Input argument does not pass verifier $(v.f), got: $val"))
                end
                updated[k] = val
            end
        end
    end
    return updated
end

using MLStyle

macro terms(args...)
    terms_impl(args...)
end

function terms_impl(args...)
    x, verifier = verifier_impl(:(x = @terms $(args...)))
    return verifier
end

function verifier_impl(expr)
    @match expr begin
        :($x = @terms $line $body $(kwi...)) => begin
            sym, forfield = @match x begin
                :($b.$a) => ("$b.$a", String(b))
                :($a) => ("$a", "")
            end
            optional = false
            for arg in kwi
                @match arg begin
                    :(optional) => (optional = true)
                end
            end
            @match body begin
                :(begin $(args...) end) => begin
                    keys = String[]
                    vals = []
                    for arg in args
                        arg isa LineNumberNode && continue
                        k, v = verifier_impl(arg)
                        push!(keys, String(k))
                        push!(vals, v)
                    end
                    ex = :($VerifierTree($([:($k=>$v) for (k, v) in zip(keys, vals)]...); optional=$optional, forfield=$forfield))
                    sym, ex
                end
                Expr(:$, subexpr) => begin
                    sym, :($VerifierTree($subexpr.dict; optional=$optional, forfield=$forfield))
                end
            end
        end
        :($x = @check $line $(args...)) => begin
            sym, type, forfield = @match x begin
                :($b.$a::$T) => ("$b.$a", T, String(b))
                :($a::$T) => ("$a", T, "")
                :($b.$a) => ("$b.$a", :Any, String(b))
                :($a) => ("$a", :Any, "")
            end
            f, optional, default = :(x->true), false, NoDefault()
            for arg in args
                @match arg begin
                    Expr(:(->), a, b) => (f = arg)
                    :(optional) => (optional = true)
                    :(default = $val) => (default = val)
                end
            end
            lambda = LegibleLambdas.parse_lambda(f)
            ex = :($Verifier($type, $lambda; default=$default, optional=$optional, forfield=$forfield))
            sym, ex
        end
    end
end

const graph_verifier = @terms begin
    nv::Int=@check x->1<=x
    edges::Vector{Vector{Int}}=@check x->all(e->length(e)==2, x)
end

const optimizer_verifier = @terms begin
    method::String = @check x->x ∈ ["TreeSA", "GreedyMethod"] default="GreedyMethod"
    method.TreeSA = @terms begin
        sc_target::Float64 = @check default=20.0
        sc_weight::Float64 = @check default=1.0
        rw_weight::Float64 = @check default=1.0
        ntrials::Int = @check x -> x > 0 default=3
        niters::Int = @check x -> x > 0 default=10
        nslices::Int = @check x-> x>= 0 default=0
        betas::Vector{Float64} = @check x->length(x)>0 default=collect(0.01:0.05:30.0)
    end optional
    method.GreedyMethod = @terms begin
        nrepeat::Int = @check x->length(x)>0 default=10
    end optional
end

const application_verifier = @terms begin
    api::String = @check x->x∈["solve", "opteinsum", "graph", "help"]
    api.solve = @terms begin
        property::String = @check x -> x ∈ ["SizeMax", "SizeMin",
            "SizeMax3", "SizeMin3", "CountingAll",
            "CountingMax", "CountingMin", "CountingMax3", "CountingMin3", "GraphPolynomial",
            "SingleConfigMax", "SingleConfigMin", "SingleConfigMax3", "SingleConfigMin3",
            "ConfigsMaxTree", "ConfigsMin",
            "ConfigsMaxTree3", "ConfigsMin3",
            "ConfigsAll", "ConfigsAllTree"
            ]
        problem::String = @check x -> x ∈ ["MaximalIS", "IndependentSet", "DominatingSet",
                "MaxCut", "Coloring{3}", "Matching"]
        weights::Vector = @check optional
        openvertices::Vector = @check default=[]
        fixedvertices::Dict = @check default=Dict()
        cudadevice::Int = @check default=-1
        optimizer = @terms $(optimizer_verifier) optional
        graph = @terms $(graph_verifier)
    end
    api.opteinsum = @terms begin
        inputs::Vector{Vector{Int}}=@check
        output::Vector{Int}=@check
        sizes::Dict{Int,Int}=@check
        optimizer = @terms $(optimizer_verifier) optional
        simplifier::String = @check x->x ∈ ["MergeGreedy", "MergeVectors"] optional
    end
    api.graph = @terms begin
        type::String = @check x->x ∈["kings", "square", "regular", "smallgraph"]
        type.kings = @terms begin
            m::Int = @check x->x>0
            n::Int = @check x->x>0
            filling::Float64= @check x->x>0 default=1.0
            seed::Int = @check default=42
        end
        type.square = @terms begin
            m::Int = @check x->x>0
            n::Int = @check x->x>0
            filling::Float64 = @check x->x>0 default=1.0
            seed::Int = @check default=42
        end
        type.smallgraph = @terms begin
            name::String = @check
        end
        type.regular = @terms begin
            d::Int = @check x->x>0
            n::Int = @check x->x>0
            seed::Int = @check default=42
        end
    end
end


text:
using Markdown
using REPL
using Documenter
using Documenter: Utilities, Expanders, Documents
using Documenter.Utilities: Selectors
using Documenter.Expanders: ExpanderPipeline, iscode, _any_color_fmt, droplines, prepend_prompt, remove_sandbox_from_output
using UUIDs

abstract type CastBlocks <: ExpanderPipeline end

Selectors.order(::Type{CastBlocks}) = 12.0
Selectors.matcher(::Type{CastBlocks}, node, page, doc) = iscode(node, r"^@cast")

# Slightly modified from:
# https://github.com/JuliaDocs/Documenter.jl/blob/68dbd53d4ff6b339e795a4a3328955ad5c689e0e/src/Expanders.jl#L638-L696
function Selectors.runner(::Type{CastBlocks}, x, page, doc)
    matched = match(r"^@cast(?:\s+([^\s;]+))?\s*(;.*)?$", x.language)
    matched === nothing && error("invalid '@cast' syntax: $(x.language)")
    name, kwargs = matched.captures
    # The sandboxed module -- either a new one or a cached one from this page.
    mod = Utilities.get_sandbox_module!(page.globals.meta, "atexample", name)

    # "parse" keyword arguments to repl
    ansicolor = _any_color_fmt(doc)
    delay = 0.5
    if kwargs !== nothing
        matched = match(r"\bansicolor\s*=\s*(true|false)\b", kwargs)
        if matched !== nothing
            ansicolor = matched[1] == "true"
        end

        # match integer delay
        matched = match(r"\bdelay\s*=\s*([0-9]+)", kwargs)
        if matched !== nothing
            delay = convert(Float64, parse(Int, matched[1]))
        else
            # match float delay
            matched = match(r"\bdelay\s*=\s*((?:[0-9]*[.])?[0-9]+)", kwargs)
            if matched !== nothing
                delay = parse(Float64, matched[1])
            end
        end
    end

    multicodeblock = Markdown.Code[]
    cast = Cast(IOBuffer(); delay=delay)
    cast_from_string!(x.code, cast; doc=doc, page=page, ansicolor=ansicolor, mod=mod, multicodeblock=multicodeblock)

    name = "$(uuid4()).cast"
    raw_html = Documents.RawHTML("""<asciinema-player src="./assets/casts/$(name)" idle-time-limit="1" autoplay="true" start-at="0.3"></asciinema-player >""")

    path = joinpath(page.workdir, "assets", "casts", name)
    mkpath(dirname(path))
    write(path, take!(cast.write_handle))

    page.mapping[x] = Documents.MultiOutput([Documents.MultiCodeBlock("julia-repl", multicodeblock), raw_html])
end

Base.@kwdef struct FakeDoc
    internal = (; errors = [])
end

Base.@kwdef struct FakePage
    workdir = pwd()
end
Utilities.locrepr(::FakePage) = nothing

function cast_from_string!(code_string::AbstractString, cast::Cast; doc=FakeDoc(), page=FakePage(), ansicolor=true, mod = Module(), multicodeblock = Markdown.Code[])
    linenumbernode = LineNumberNode(0, "REPL") # line unused, set to 0
    @debug "Evaluating @cast:\n$(x.code)"
    for (ex, str) in Utilities.parseblock(code_string, doc, page; keywords = false,
                                          linenumbernode = linenumbernode)
        buffer = IOBuffer()
        input  = droplines(str)
        if !isempty(input)
            push!(multicodeblock, Markdown.Code("julia-repl", prepend_prompt(input)))
            write_event!(cast, InputEvent, input)
            write_event!(cast, OutputEvent, prepend_prompt(input) * "\n")
        end

        if VERSION >= v"1.5.0-DEV.178"
            # Use the REPL softscope for REPLBlocks,
            # see https://github.com/JuliaLang/julia/pull/33864
            ex = REPL.softscope!(ex)
        end
        c = capture(cast; rethrow = InterruptException, color = ansicolor,
        process = str -> remove_sandbox_from_output(str, mod)) do
            cd(page.workdir) do
                Core.eval(mod, ex)
            end
        end
        Core.eval(mod, Expr(:global, Expr(:(=), :ans, QuoteNode(c.value))))
        output = if !c.error
            hide = REPL.ends_with_semicolon(input)
            Documenter.DocTests.result_to_string(buffer, hide ? nothing : c.value)
        else
            Documenter.DocTests.error_to_string(buffer, c.value, [])
        end

        out = IOBuffer()
        if isempty(input) || isempty(output)
            println(out)
        else
            println(out, output, "\n")
        end
        outstr = String(take!(out))
        # Replace references to gensym'd module with Main
        outstr = remove_sandbox_from_output(outstr, mod)
        write_event!(cast, OutputEvent, outstr)
    end
end

"""
    @cast_str(code_string, delay=0) -> Cast

Creates a [`Cast`](@ref) object by executing the code in `code_string` in a REPL-like environment.

## Example

```julia
using Asciicast

c = cast"x=1"0.5 # note we set a delay of 0.5 here

Asciicast.save("test.cast", c)
```

"""
macro cast_str(code_string, delay=0)
    cast = Cast(IOBuffer(); delay=delay)
    cast_from_string!(code_string, cast)
    return cast
end


text:
module SubstitutionTilings

export substitute, dilate, empirical_frequency, draw, embed_aff, Tiling, SubSystem, id, vertices, in_border, transition_matrix, DGroupElem
export transition_matrix
export collar_in, is_interior, frequency, total_collaring, total_collaring, UnrecognizedCollar
export embed_center
export autocorrelation

include("CoreDefs.jl")
using .CoreDefs
include("NumFields.jl")
using .NumFields
#include("Collaring.jl")
#using .Collaring
#include("systems/heisenberg.jl")
#using .Heisenberg
include("systems/penrose.jl")
using .Penrose
include("systems/pinwheel.jl")
using .Pinwheel
#include("systems/heisenberg-improper.jl")
#using .HeisenbergImproper
include("systems/nilpotent.jl")
using .Nilpotent
include("systems/chair.jl")
using .Chair

end


text:
using Minerva, Flux
using Flux: mse

# TODO: BatchNorm?
# TODO: Dropout?
# TODO: Regularization?

# agent settings
max_memory = 1000
discount_rate = 1.00
update_episodes = 5

# network settings
batch_size = 32
learning_rate = 1e-2

Q = Chain(
    Dense(4, 64),
    Dense(64, 64),
    Dense(64, 2)
)
agent = DeepQAgent(Q, learning_rate)
env = CartPole()
ϵ_schedule = StepDecay(1.0, 0.05, 0.1, 1)
recorder = train!(agent, env, max_episodes = 500, ϵ_schedule = ϵ_schedule)
evaluate(agent, env, 1, ϵ = 0.05, rendered = true)


text:
using InvariantSets
using Test

@testset "InvariantSets.jl" begin
    @testset "minimal_invariant_set" begin
        include(raw"minimal_invariant_set/rakovic_algorithm.jl")
        include(raw"minimal_invariant_set/schulze_algorithm.jl")
    end
    @testset "maximal_invariant_set" begin
        include(raw"maximal_invariant_set/maximal_RPI_set.jl")
        include(raw"maximal_invariant_set/maximal_CPI_set.jl")
    end
    @testset "concrete_operations" begin
        include(raw"concrete_operator.jl")
    end
    @testset "Helper Functionalities" begin
        include(raw"helper/control_systems.jl")
        include(raw"helper/jump.jl")
    end
end

max_stars_repo_name: mforets/InvariantSets

max_issues_repo_name: mforets/InvariantSets


text:
using LinearAlgebra, Distributions, Plots
include("GaussianProcesses.jl")

# distance function (squared)
# for euclidean, just use norm(v, p)
d_euclid(a::Array, b::Array) = norm(a-b, 2)

# example function
f_target(x) = sqrt(x) * sin(x/10)
noise_distro = Normal(0.8)

x_train = collect(range(0, 100, length = 20));
#= x_train = [collect(range(0, 49, length = 50))..., collect(range(50, 100, length = 10))...]; =#
y_train = f_target.(x_train) .+ rand.(noise_distro);
plot(x_train, y_train)

# plotting grid
samplex = collect(range(0, 100, length = 10));

# instantiate a GP
myGP = GaussianProcesses.GP_Regressor(x_train, y_train, samplex; σ=0.8, h=20.0, d=d_euclid);

# draw one sample function from prior
plot(myGP.x_grid, rand(myGP.prior))

# check posterior mean
plot(myGP.x_grid, mean(myGP.posterior))

# check posterior standard deviation
plot(myGP.x_grid, sqrt.(diag(cov(myGP.posterior))))

# plot equivalent kernel weights of train data
plot(myGP.x_train, myGP.equiv_weights[10,:])

# multidimensional example

using DataFrames

dt = DataFrame(
	x1=[-5,-5,-5,0,0,0,5,5,5],
	x2=[-5,0,5,-5,0,5,-5,0,5]
);

f(x1,x2) = -x1^2 - x2^2

dt[!, :y] .= f.(dt.x1, dt.x2);
dt

max_stars_repo_name: tobanw/gaussian-process-regression

max_issues_repo_name: tobanw/gaussian-process-regression

max_forks_repo_name: tobanw/gaussian-process-regression


text:
cd(@__DIR__)
KitML.load_data("dataset.csv", dlm = ",")

nn1 = KitML.load_model("jmodel.jld2")
nn2 = KitML.load_model("tfmodel.h5")

KitML.save_model(nn1)
KitML.save_model(nn2; mode = :tf)


text:
function test_confrac_iterator()
    cf = confrac(Int[])
    @test collect(cf) == []
    @test eltype(cf) == Int
    @test length(cf) == 0

    cf = confrac([1])
    @test collect(cf) == [1]
    @test eltype(cf) == Int
    @test length(cf) == 1

    cf = confrac(BigInt[1])
    @test collect(cf) == [1]
    @test eltype(cf) == BigInt
    @test length(cf) == 1

    cf = confrac([1, 2, 3, 4])
    @test collect(cf) == [1, 2, 3, 4]
    @test eltype(cf) == Int
    @test length(cf) == 4

    cf = confrac([1], 1)
    @test collect(take(cf, 3)) == [1, 1, 1]
    @test eltype(cf) == Int

    cf = confrac([1, 2, 3], 2)
    @test collect(take(cf, 5)) == [1, 2, 3, 2, 3]
    @test eltype(cf) == Int
end

function test_confrac_rationalize()
    cf = confrac(Int[])
    @test rationalize(cf) == 1 // 0
    @test rationalize(cf, 0) == 1 // 0

    cf = confrac([1])
    @test rationalize(cf) == 1
    @test rationalize(cf, 0) == 1 // 0
    @test rationalize(cf, 1) == 1

    cf = confrac([1, 2, 3, 4])
    @test rationalize(cf) == 43 // 30
    @test rationalize(cf, 0) == 1 // 0
    @test rationalize(cf, 1) == 1
    @test rationalize(cf, 2) == 3 // 2
    @test rationalize(cf, 3) == 10 // 7
    @test rationalize(cf, 4) == 43 // 30

    cf = confrac([1], 1)
    @test rationalize(cf) == 1
    @test rationalize(cf, 0) == 1 // 0
    @test rationalize(cf, 1) == 1
    @test rationalize(cf, 2) == 2
    @test rationalize(cf, 3) == 3 // 2

    cf = confrac([1, 2, 3], 2)
    @test rationalize(cf) == 10 // 7
    @test rationalize(cf, 0) == 1 // 0
    @test rationalize(cf, 1) == 1
    @test rationalize(cf, 2) == 3 // 2
    @test rationalize(cf, 3) == 10 // 7
    @test rationalize(cf, 4) == 23 // 16
    @test rationalize(cf, 5) == 79 // 55
end

function test_confrac_confrac()
    @test confrac(0 // 1) == confrac([0])

    @test confrac(3 // 7) == confrac([0, 2, 3])
    @test confrac(10 // 7) == confrac([1, 2, 3])

    @test confrac(-3 // 7) == confrac([-1, 1, 1, 3])
    @test confrac(-10 // 7) == confrac([-2, 1, 1, 3])

    @test confrac(3, 0, 7) == confrac([0, 2, 3])
    @test confrac(-3, 0, 7) == confrac([-1, 1, 1, 3])
    @test confrac(1, 4, 7) == confrac([0, 2, 3])

    @test confrac(0, 2, 1) == confrac([1, 2], 2)
    @test confrac(0, 14, 1) == confrac([3, 1, 2, 1, 6], 2)
    @test confrac(0, 19, 1) == confrac([4, 2, 1, 3, 1, 2, 8], 2)
    @test confrac(0, 42, 1) == confrac([6, 2, 12], 2)
    @test confrac(0, 48, 4) == confrac([1, 1, 2], 2)

    @test confrac(-41, 313, 36) == confrac([-1, 2, 1, 5, 8, 1, 2, 17, 2, 1, 8, 5, 1, 3, 1, 1, 2, 2, 1, 1, 3], 3)
    @test confrac(41, 313, -36) == confrac([-2, 2, 1, 2, 2, 1, 1, 3, 1, 5, 8, 1, 2, 17, 2, 1, 8, 5, 1, 3, 1], 3)
    @test confrac(-725, 313, 608) == confrac([-2, 1, 5, 8, 5, 1, 3, 1, 1, 2, 2, 1, 1, 3, 1, 5, 8, 1, 2, 17, 2, 1], 4)
    @test confrac(725, 313, -608) == confrac([-2, 1, 3, 1, 1, 17, 2, 1, 8, 5, 1, 3, 1, 1, 2, 2, 1, 1, 3, 1, 5, 8, 1, 2], 6)
end

function test_confrac_convergents()
    cf = confrac(Int[])
    @test collect(convergents(cf)) == []
    @test eltype(convergents(cf)) == Tuple{Int,Int}
    @test length(convergents(cf)) == 0

    cf = confrac([1])
    @test collect(convergents(cf)) == [(1, 1)]
    @test eltype(convergents(cf)) == Tuple{Int,Int}
    @test length(convergents(cf)) == 1

    cf = confrac(BigInt[1])
    @test collect(convergents(cf)) == [(1, 1)]
    @test eltype(convergents(cf)) == Tuple{BigInt,BigInt}
    @test length(convergents(cf)) == 1

    cf = confrac([1, 2, 3, 4])
    @test collect(convergents(cf)) == [(1, 1), (3, 2), (10, 7), (43, 30)]
    @test eltype(convergents(cf)) == Tuple{Int,Int}
    @test length(convergents(cf)) == 4

    cf = confrac([1], 1)
    @test collect(take(convergents(cf), 3)) == [(1, 1), (2, 1), (3, 2)]
    @test eltype(convergents(cf)) == Tuple{Int,Int}

    cf = confrac([1, 2, 3], 2)
    @test collect(take(convergents(cf), 5)) == [(1, 1), (3, 2), (10, 7), (23, 16), (79, 55)]
    @test eltype(convergents(cf)) == Tuple{Int,Int}
end

function test_confrac_all()
    print(rpad("Math.NumberTheory.ConFrac...", 50, ' '))

    test_confrac_iterator()
    test_confrac_rationalize()
    test_confrac_confrac()
    test_confrac_convergents()

    println("PASS")
end


text:
# Adding Semester to Date.DatePeriod
struct Semester <: Dates.DatePeriod
    value::Int64
    Semester(v::Number) = new(v)
end

# Adding Undated to Dates.DatePeriod
struct Undated <: Dates.DatePeriod
    value::Int64
    Undated(v::Number) = new(v)
end
value(d::Semester) = d.value
value(d::Undated) = d.value

abstract type SimpleDate <: Dates.TimeType end

SIMPLEDATES = (:Year, :Semester, :Quarter, :Month, :Week, :Day, :Undated)

# weeks computations
function week_per_year()
    weeks = Dict()
    weeks[0] = Dict([("max_weeks", 52), ("cum_weeks", 0)])
    cummulativeweeks = 0
    for y = 1:2500
        w = Dates.week(Dates.Date(y, 12, 28))
        cummulativeweeks += w
        weeks[y] = Dict([("max_weeks", w), ("cum_weeks", cummulativeweeks)])
    end
    return weeks
end

if !@isdefined WEEKTABLE
    const WEEKTABLE = week_per_year()
end

"""
    UTInstant{T}

The `UTInstant` represents a machine timeline based on UT time (1 day = one revolution of
the earth). The `T` is a `Period` parameter that indicates the resolution or precision of
the instant.
"""
struct UTInstant{P<:Dates.DatePeriod} <: Dates.Instant
    periods::P
end

# Adding new date types to SimpleDate

for (T1, T2) in (
    (:Year, :YearDate),
    (:Quarter, :QuarterDate),
    (:Month, :MonthDate),
    (:Week, :WeekDate),
    (:Day, :DayDate),
)
    @eval begin
        struct $T2 <: SimpleDate
            instant::Dates.UTInstant{Dates.$T1}
            $T2(instant::Dates.UTInstant{Dates.$T1}) = new(instant)
        end
    end
end

# Adding new periods to TimeType
for (T1, T2) in ((:Semester, :SemesterDate), (:Undated, :UndatedDate))
    @eval begin
        struct $T2 <: SimpleDate
            instant::UTInstant{$T1}
            $T2(instant::UTInstant{$T1}) = new(instant)
        end
    end
end

# frequency traits
abstract type Frequency end
struct FYear <: Frequency; end
struct FSemester <: Frequency; end
struct FQuarter <: Frequency; end
struct FMonth <: Frequency; end
struct FWeek <: Frequency; end
struct FDay <: Frequency; end
struct FUndated <: Frequency; end

simpleperiod(::Type{FYear}, y::Int64) = YearDate(y)
simpleperiod(::Type{FSemester}, y::Int64, s::Int64) = SemesterDate(y, s)
simpleperiod(::Type{FQuarter}, y::Int64, q::Int64) = QuarterDate(y, q)
simpleperiod(::Type{FMonth}, y::Int64, m::Int64) = MonthDate(y, m)
simpleperiod(::Type{FWeek}, y::Int64, w::Int64) = WeekDate(y, w)
simpleperiod(::Type{FDay}, y::Int64, m::Int64, d::Int64) = DayDate(y, m, d)
simpleperiod(::Type{FUndated}, u::Int64) = UndatedDate(u)

DAYMULTIPLEDATE = Union{DayDate,WeekDate}
DAYMULTIPLEPERIOD = Union{Day,Week}
MONTHMULTIPLEDATE = Union{MonthDate,QuarterDate,SemesterDate,YearDate}
MONTHMULTIPLEPERIOD = Union{Month,Quarter,Semester,Year}

# Convenience default constructors
UTY(x) = Dates.UTInstant(Dates.Year(x))
UTS(x) = UTInstant(Semester(x))
UTQ(x) = Dates.UTInstant(Dates.Quarter(x))
UTM(x) = Dates.UTInstant(Dates.Month(x))
UTW(x) = Dates.UTInstant(Dates.Week(x))
UTD(x) = Dates.UTInstant(Dates.Day(x))
UTU(x) = UTInstant(Undated(x))

"""
    argerror([msg]) -> Union{ArgumentError, Nothing}

Return an `ArgumentError` object with the given message,
or [`nothing`](@ref) if no message is provided. For use by `validargs`.
"""
argerror(msg::String) = ArgumentError(msg)
argerror() = nothing

### CONSTRUCTORS ###
# Core constructors
"""
    YearDate(y) -> Year

Construct a `YearDate` type. Argument must be convertible to [`Int64`](@ref).
"""
function YearDate(y::Int64)
    return YearDate(UTY(yearfromepoch(y)))
end

function yearfromepoch(y)
    return y - 1
end

"""
    SemesterDate(y, [q]) -> SemesterDate

Construct a `SemesterDate` type by parts. Arguments must be convertible to [`Int64`](@ref).
"""
function SemesterDate(y::Int64, s::Int64 = 1)
    err = validargs(Semester, y, s)
    err === nothing || throw(err)
    return SemesterDate(UTS(totalsemesters(y, s)))
end

function validargs(::Type{Semester}, y::Int64, s::Int64)
    0 < s < 3 || return argerror("Semester: $s out of range (1:2)")
    return argerror()
end

function totalsemesters(y, s)
    return 2 * yearfromepoch(y) + s
end

"""
    QuarterDate(y, [q]) -> QuarterDate

Construct a `QuarterDate` type by parts. Arguments must be convertible to [`Int64`](@ref).
"""
function QuarterDate(y::Int64, q::Int64 = 1)
    err = validargs(Quarter, y, q)
    err === nothing || throw(err)
    return QuarterDate(UTQ(totalquarters(y, q)))
end

function validargs(::Type{Quarter}, y::Int64, q::Int64)
    0 < q < 5 || return argerror("Quarter: $q out of range (1:4)")
    return argerror()
end

function totalquarters(y, q)
    return 4 * yearfromepoch(y) + q
end

"""
    MonthDate(y, [m]) -> MonthDate

Construct a `MonthDate` type by parts. Arguments must be convertible to [`Int64`](@ref).
"""
function MonthDate(y::Int64, m::Int64 = 1)
    err = validargs(Month, y, m)
    err === nothing || throw(err)
    return MonthDate(UTM(totalmonths(y, m)))
end

function validargs(::Type{Month}, y::Int64, m::Int64)
    0 < m < 13 || return argerror("Month: $m out of range (1:12)")
    return argerror()
end

function totalmonths(y, m)
    return 12 * yearfromepoch(y) + m
end

"""
    WeekDate(y, [w]) -> WeekDate

Construct a `WeekDate` type by parts. Arguments must be convertible to [`Int64`](@ref).
"""
function WeekDate(y::Int64, w::Int64 = 1)
    err = validargs(Week, y, w)
    err === nothing || throw(err)
    return WeekDate(UTW(totalweeks(y, w)))
end

function validargs(::Type{Week}, y::Int64, w::Int64)
    0 < y < 2500 || return argerror("Week: $y out of range (0:2500)")
    0 < w < 54 || return argerror("Week: $w out of range (1:53)")
    return argerror()
end

function totalweeks(y, w)
    return WEEKTABLE[y-1]["cum_weeks"] + w
end

"""
    DayDate(y, [m, d]) -> DayDate

Construct a `DayDate` type by parts. Arguments must be convertible to [`Int64`](@ref).
A `DayDate``is identical to  a `Dates.Date`
"""
function DayDate(y::Int64, m::Int64 = 1, d::Int64 = 1)
    err = Dates.validargs(Dates.Date, y, m, d)
    err === nothing || throw(err)
    return DayDate(Dates.UTD(Dates.totaldays(y, m, d)))
end

"""
    UndatedDate(u) -> UndatedDate

Construct a `UndatedDate` type. Argument must be convertible to [`Int64`](@ref).
"""
function UndatedDate(u::Int64)
    return UndatedDate(UTU(u))
end

# Using Dates.Periods
DayDate(y::Dates.Year, m::Dates.Month = Dates.Month(1), d::Dates.Day = Dates.Day(1)) =
    DayDate(Dates.value(y), Dates.value(m), Dates.value(d))
WeekDate(y::Dates.Year, w::Dates.Week = Dates.Week(1)) =
    WeekDate(Dates.value(y), Dates.value(w))
MonthDate(y::Dates.Year, m::Dates.Month = Dates.Month(1)) =
    MonthDate.(Dates.value(y), Dates.value(m))
QuarterDate(y::Dates.Year, q::Dates.Quarter = Dates.Quarter(1)) =
    QuarterDate(Dates.value(y), Dates.value(q))
SemesterDate(y::Dates.Year, s::Semester = Semester(1)) =
    SemesterDate(Dates.value(y), Dates.value(s))
YearDate(y::Dates.Year) = YearDate(Dates.value(y))

# To allow any order/combination of Dates.Periods
"""
    Semester(period::Dates.Period...) -> Semester

Construct a `Semester` type by `Period` type parts. Arguments may be in any order. `Semester` parts
not provided will default to the value of `Dates.default(period)`.
"""
function SemesterDate(period::Dates.Period, periods::Dates.Period...)
    y = Dates.Year(1)
    s = Semester(1)
    for p in (period, periods...)
        isa(p, Dates.Year) && (y = p::Dates.Year)
        isa(p, Semester) && (s = p::Semester)
    end
    return SemesterDate(y, s)
end

function QuarterDate(period::Dates.Period, periods::Dates.Period...)
    y = Dates.Year(1)
    q = Dates.Quarter(1)
    for p in (period, periods...)
        isa(p, Dates.Year) && (y = p::Dates.Year)
        isa(p, Dates.Quarter) && (q = p::Dates.Quarter)
    end
    return QuarterDate(y, q)
end

function MonthDate(period::Dates.Period, periods::Dates.Period...)
    y = Dates.Year(1)
    m = Dates.Month(1)
    for p in (period, periods...)
        isa(p, Dates.Year) && (y = p::Dates.Year)
        isa(p, Dates.Month) && (m = p::Dates.Month)
    end
    return MonthDate(y, m)
end

function WeekDate(period::Dates.Period, periods::Dates.Period...)
    y = Dates.Year(1)
    w = Dates.Week(1)
    for p in (period, periods...)
        isa(p, Dates.Year) && (y = p::Dates.Year)
        isa(p, Dates.Week) && (w = p::Dates.Week)
    end
    return WeekDate(y, w)
end

function DayDate(period::Dates.Period, periods::Dates.Period...)
    y = Dates.Year(1)
    m = Dates.Month(1)
    d = Dates.Day(1)
    for p in (period, periods...)
        isa(p, Dates.Year) && (y = p::Dates.Year)
        isa(p, Dates.Month) && (m = p::Dates.Month)
        isa(p, Dates.Day) && (d = p::Dates.Day)
    end
    return DayDate(y, m, d)
end

# Fallback constructors
YearDate(y) = YearDate(Int64(y))
SemesterDate(y, s = 1) = SemesterDate(Int64(y), Int64(s))
QuarterDate(y, q = 1) = QuarterDate(Int64(y), Int64(q))
MonthDate(y, m = 1) = MonthDate(Int64(y), Int64(m))
WeekDate(y, w = 1) = WeekDate(Int64(y), Int64(w))
DayDate(y, m = 1, d = 1) = DayDate(Int64(y), Int64(m), Int64(d))
UndatedDate(y) = UndatedDate(Int64(y))

# Traits, Equality
Base.isfinite(::Union{Type{T},T}) where {T<:SimpleDate} = true
calendar(dt::T) where {T<:SimpleDate} = Dates.ISOCalendar

"""
    eps(::Type{YearDate}) -> Dates.Year
    eps(::Type{Semester}) -> Semester
    eps(::Type{Quarter}) -> Dates.Quarter
    eps(::Type{Month}) -> Dates.Month
    eps(::Type{Week}) -> Dates.Week
    eps(::Type{Day}) -> Dates.Day
    eps(::Type{Undated}) -> 1
    eps(::SimpleDate) -> Dates.Period

Return the smallest unit value supported by the `SimpleDate`.

# Examples
```jldoctest
julia> eps(Year)
1 year

julia> eps(Day)
1 day
"""
Base.eps(
    ::Union{
        Type{YearDate},
        Type{SemesterDate},
        Type{QuarterDate},
        Type{MonthDate},
        Type{WeekDate},
        Type{DayDate},
        Type{UndatedDate},
        SimpleDate,
    },
)

Base.eps(::Type{YearDate}) = Dates.Year(1)
Base.eps(::Type{SemesterDate}) = Semester(1)
Base.eps(::Type{QuarterDate}) = Dates.Quarter(1)
Base.eps(::Type{MonthDate}) = Dates.Month(1)
Base.eps(::Type{WeekDate}) = Dates.Week(1)
Base.eps(::Type{DayDate}) = Dates.Day(1)
Base.eps(::Type{UndatedDate}) = Int64(1)
Base.eps(::T) where {T<:SimpleDate} = eps(T)::Dates.Period

# zero returns dt::T - dt::T
Base.zero(::Type{YearDate}) = Dates.Year(0)
Base.zero(::Type{SemesterDate}) = Semester(0)
Base.zero(::Type{QuarterDate}) = Dates.Quarter(0)
Base.zero(::Type{MonthDate}) = Dates.Month(0)
Base.zero(::Type{WeekDate}) = Dates.Week(0)
Base.zero(::Type{DayDate}) = Dates.Day(0)
Base.zero(::Type{UndatedDate}) = Int64(0)
Base.zero(::T) where {T<:SimpleDate} = zero(T)::Dates.Period

Base.typemax(::Union{YearDate,Type{YearDate}}) = YearDate(252522163911149)
Base.typemin(::Union{YearDate,Type{YearDate}}) = YearDate(-252522163911150)
Base.typemax(::Union{SemesterDate,Type{SemesterDate}}) = SemesterDate(252522163911149, 2)
Base.typemin(::Union{SemesterDate,Type{SemesterDate}}) = SemesterDate(-252522163911150, 1)
Base.typemax(::Union{QuarterDate,Type{QuarterDate}}) = QuarterDate(252522163911149, 4)
Base.typemin(::Union{QuarterDate,Type{QuarterDate}}) = QuarterDate(-252522163911150, 1)
Base.typemax(::Union{MonthDate,Type{MonthDate}}) = MonthDate(252522163911149, 12)
Base.typemin(::Union{MonthDate,Type{MonthDate}}) = MonthDate(-252522163911150, 1)
Base.typemax(::Union{WeekDate,Type{WeekDate}}) = WeekDate(2499, 52)
Base.typemin(::Union{WeekDate,Type{WeekDate}}) = WeekDate(1, 1)
Base.typemax(::Union{DayDate,Type{DayDate}}) = DayDate(252522163911149, 12, 31)
Base.typemin(::Union{DayDate,Type{DayDate}}) = DayDate(-252522163911150, 1, 1)
Base.typemax(::Union{UndatedDate,Type{UndatedDate}}) = UndatedDate(typemax(Int64))
Base.typemin(::Union{UndatedDate,Type{UndatedDate}}) = UndatedDate(typemin(Int64))

# Periods promotion, isless, ==
Base.promote_rule(::Type{Dates.Year}, ::Type{Semester}) = Semester
Base.promote_rule(::Type{Semester}, ::Type{Dates.Quarter}) = Dates.Quarter
Base.promote_rule(::Type{Semester}, ::Type{Dates.Month}) = Dates.Month
Base.promote_rule(::Type{Semester}, ::Type{Dates.Week}) = Dates.Week
Base.promote_rule(::Type{Semester}, ::Type{Dates.Day}) = Dates.Day

# other periods with fixed conversions but which aren't fixed time periods
const OtherPeriod = Union{Dates.Month,Dates.Quarter,Semester,Dates.Year}

let vmax = typemax(Int64) ÷ 2, vmin = typemin(Int64) ÷ 2
    @eval function Base.convert(::Type{Semester}, x::Dates.Year)
        $vmin ≤ value(x) ≤ $vmax || throw(InexactError(:convert, Semester, x))
        Semester(value(x) * 2)
    end
end
Base.convert(::Type{Dates.Year}, x::Semester) = Dates.Year(divexact(value(x), 2))
#Base.promote_rule(::Type{Dates.Year}, ::Type{Semester}) = Semester

let vmax = typemax(Int64) ÷ 6, vmin = typemin(Int64) ÷ 6
    @eval function Base.convert(::Type{Dates.Month}, x::Semester)
        $vmin ≤ value(x) ≤ $vmax || throw(InexactError(:convert, Dates.Month, x))
        Dates.Month(value(x) * 6)
    end
end
Base.convert(::Type{Semester}, x::Dates.Month) = Semester(divexact(value(x), 6))
#Base.promote_rule(::Type{Semester}, ::Type{Dates.Month}) = Dates.Month

let vmax = typemax(Int64) ÷ 2, vmin = typemin(Int64) ÷ 2
    @eval function Base.convert(::Type{Dates.Quarter}, x::Semester)
        $vmin ≤ value(x) ≤ $vmax || throw(InexactError(:convert, Dates.Quarter, x))
        Dates.Quarter(value(x) * 3)
    end
end
Base.convert(::Type{Semester}, x::Dates.Quarter) = Dates.Quarter(divexact(value(x), 2))
#Base.promote_rule(::Type{Semester}, ::Type{Dates.Quarter}) = Dates.Quarter

# truncating conversions to milliseconds, nanoseconds and days:
# overflow can happen for periods longer than ~300,000 years
days(c::Year) = Dates.days(c)
days(c::Semester) = 182.62125 * value(c)
days(c::Quarter) = Dates.days(c)
days(c::Month) = Dates.days(c)
days(c::Week) = Dates.days(c)
days(c::Day) = Dates.days(c)

"""
    islongyear(y)

Return `true` if and only if year `y` has 53 ISO-8601 weeks.

# Examples
```julia-repl
julia> islongyear(2020)
true

julia> islongyear(2021)
false
```
"""
function islongyear(y::Integer)
    mod(trunc(Int, y) + trunc(Int, y / 4) - trunc(Int, y / 100) + trunc(Int, y / 400), 7) ==
    4 ||
        mod(
            trunc(Int, y - 1) + trunc(Int, (y - 1) / 4) - trunc(Int, (y - 1) / 100) +
            trunc(Int, (y - 1) / 400),
            7,
        ) == 3
end

"""
    weeksinyear(y)

Return the number of ISO-8601 weeks in year `y`.

# Examples
```julia-repl
julia> islongyear(2020)
53

julia> islongyear(2021)
52
```
"""
function weeksinyear(y::Integer)
    islongyear(y) ? 53 : 52
end

function weeksinyears(y1::Integer, y2::Integer)
    return sum(weeksinyear, y1:y2)
end

"""
    numberofweeks(y, w)

Return the number of ISO-8601 weeks since Base Year = 1 in week `w` of year `y`.

# Examples
```julia-repl
julia> islongyear(2020)
53

julia> islongyear(2021)
52
```
"""
function numberofweeks(y::Integer, w::Integer)
    sum(weeksinyear.(1:(y-1))) + w
end

function get_thursday(w::WeekDate)
    return 7 * (value(w) - 1) + 4
end


text:
"""
Write mdp transitions in a ".tra" file to give as an input to storm.
The transition probabilities are written according to storm explicit format as follows:
```
mdp
0 0 1 0.3
0 0 4 0.7
0 1 0 0.5
0 1 1 0.5
1 0 1 1.0
...
```
In our usual MDP formulation a row is: s, a, s', T(s,a,s')

    write_mdp_transition(mdp::MDP, filename::String="mdp.tra")
"""
function write_mdp_transition(mdp::MDP, filename::String="mdp.tra")
    states = ordered_states(mdp)
    actions = ordered_actions(mdp)
    open(filename, "w") do f
        write(f, "mdp \n")
        for s in states
            si = stateindex(mdp, s)
            si -= 1 # 0-indexed
            for a in actions
                ai = actionindex(mdp, a) - 1
                d = transition(mdp, s, a)
                for (sp, p) in weighted_iterator(d)
                    spi = stateindex(mdp, sp) - 1
                    p == 0.0 ? continue : nothing
                    line = string(si, " ", ai, " ", spi, " ", p, "\n")
                    write(f, line)
                end
            end
        end
    end
end

"""
Write state labels in a ".lab" file to give as input to storm

File format:
Define labels and write the mapping from states to labels in a file using the following syntax:
```
#DECLARATION
label1, label2, label3
#END
0 label1 label2
15 label2
23 label3
```

All the possible labels must be declared first. A state can have several labels. In the file the states must be ordered.

    write_mdp_labels{S, A}(mdp::MDP{S, A}, labeling::Dict{S, Vector{String}}, filename::String="mdp.lab")
"""
function write_mdp_labels{S, A}(mdp::MDP{S, A}, labeling::Dict{S, Vector{String}}, filename::String="mdp.lab")
    labels = Set{String}()
    for labv in Set(values(labeling))
        for lab in labv
            push!(labels, lab)
        end
    end
    open(filename, "w") do f
        write(f, "#DECLARATION\n")
        for label in labels
            write(f, "$label ")
        end
        write(f, "\n")
        write(f, "#END\n")
        # states must be sorted by index
        for (i, s) in enumerate(ordered_states(mdp))
            if haskey(labeling, s)
                si = stateindex(mdp, s) - 1 #0 indexed
                write(f, string(si, " "))
                for lab in labeling[s]
                    write(f, lab, " ")
                end
                write(f, "\n")
            end
        end
    end
end


text:
@testset "version" begin

@test isa(CUDAdrv.version(), VersionNumber)

end