text:
"""
Trait indicating that the dimension is in the normal forward order.
"""
struct Forward end

"""
Trait indicating that the dimension is in the reverse order.
Selector lookup and plotting will be reverse.
"""
struct Reverse end

"""
An AbstractDimension tags the dimensions in an AbstractArray.

It can also contain spatial coordinates and their metadata. For simplicity,
the same types are used both for storing dimension information and for indexing.
"""
abstract type AbstractDimension{T,M,O} end

"""
AbstractCombined holds mapping that require multiple dimension
when `select()` is used, shuch as for situations where they share an
affine map or similar transformation instead of linear maps. Each dim
that shares the mapping must contain the same (identical) object.

Dimensions holding a DimCombination will work as usual for direct indexing.

All AbstractDimension are assumed to `val` and `metadata` fields.
"""
abstract type AbstractDimCombination end

const AbDim = AbstractDimension
const AbDimType = Type{<:AbDim}
const AbDimTuple = Tuple{Vararg{<:AbDim,N}} where N
const AbDimVector = Vector{<:AbDim}
const DimOrDimType = Union{AbDim,AbDimType}
const AllDimensions = Union{AbDim,AbDimTuple,AbDimType,
                            Tuple{Vararg{AbDimType}},
                            Vector{<:AbDim}}

# Getters
val(dim::AbDim) = dim.val
metadata(dim::AbDim) = dim.metadata
order(dim::AbDim) = dim.order

# DimensionalData interface methods
rebuild(dim::AbDim, val) = basetype(dim)(val, metadata(dim), order(dim))

dims(x::AbDim) = x
dims(x::AbDimTuple) = x
name(dim::AbDim) = name(typeof(dim))
shortname(d::AbDim) = shortname(typeof(d))
shortname(d::Type{<:AbDim}) = name(d)
units(dim::AbDim) = metadata(dim) == nothing ? "" : get(metadata(dim), :units, "")

bounds(a, args...) = bounds(dims(a), args...)
bounds(dims::AbDimTuple, lookupdims::Tuple) = bounds(dims[[dimnum(dims, lookupdims)...]]...)
bounds(dims::AbDimTuple, dim::DimOrDimType) = bounds(dims[dimnum(dims, dim)])
bounds(dims::AbDimTuple) = (bounds(dims[1]), bounds(tail(dims))...)
bounds(dims::Tuple{}) = ()
bounds(dim::AbDim) = bounds(order(dim), dim)
bounds(::Forward, dim::AbDim) = first(val(dim)), last(val(dim))
bounds(::Reverse, dim::AbDim) = last(val(dim)), first(val(dim))

# Base methods
Base.eltype(dim::Type{<:AbDim{T}}) where T = T
Base.length(dim::AbDim) = length(val(dim))
Base.show(io::IO, dim::AbDim) = begin
    printstyled(io, "\n", name(dim), ": "; color=:red)
    show(io, typeof(dim))
    printstyled(io, "\nval: "; color=:green)
    show(io, val(dim))

    # printstyled(io, indent, name(v), color=:green)
    printstyled(io, "\nmetadata: "; color=:blue)
    show(io, metadata(dim))
    print(io, "\n")
end

# AbstractArray methods where dims are the dispatch argument

@inline rebuildsliced(a, data, I) = rebuild(a, data, slicedims(a, I)...)

Base.@propagate_inbounds Base.getindex(a::AbstractArray, dims::Vararg{<:AbDim{<:Number}}) =
    getindex(a, dims2indices(a, dims)...)
Base.@propagate_inbounds Base.getindex(a::AbstractArray, dims::Vararg{<:AbstractDimension}) =
    getindex(a, dims2indices(a, dims)...)

Base.@propagate_inbounds Base.setindex!(a::AbstractArray, x, dims::Vararg{<:AbstractDimension}) =
    setindex!(a, x, dims2indices(a, dims)...)

Base.@propagate_inbounds Base.view(a::AbstractArray, dims::Vararg{<:AbstractDimension}) =
    view(a, dims2indices(a, dims)...)

@inline Base.axes(a::AbstractArray, dims::DimOrDimType) = axes(a, dimnum(a, dims))
@inline Base.size(a::AbstractArray, dims::DimOrDimType) = size(a, dimnum(a, dims))

#= SplitApplyCombine methods?
Should allow groupby using dims lookup to make this worth the dependency
Like group by time/lattitude/height band etc.

SplitApplyCombine.splitdims(a::AbstractArray, dims::AllDimensions) =
    SplitApplyCombine.splitdims(a, dimnum(a, dims))

SplitApplyCombine.splitdimsview(a::AbstractArray, dims::AllDimensions) =
    SplitApplyCombine.splitdimsview(a, dimnum(a, dims))
=#

"""
Dimensions with user-set type paremeters
"""
abstract type AbstractParametricDimension{X,T,M,O} <: AbstractDimension{T,M,O} end

"""
A generic dimension. For use when custom dims are required when loading
data from a file. The sintax is ugly and verbose to use for indexing,
ie `Dim{:lat}(1:9)` rather than `Lat(1:9)`. This is the main reason
they are not the only type of dimension availabile.
"""
struct Dim{X,T,M,O} <: AbstractParametricDimension{X,T,M,O}
    val::T
    metadata::M
    order::O
    Dim{X}(val, metadata, order) where X =
        new{X,typeof(val),typeof(metadata),typeof(order)}(val, metadata, order)
end

@inline Dim{X}(val=:; metadata=nothing, order=Forward()) where X =
    Dim{X}(val, metadata, order)
name(::Type{<:Dim{X}}) where X = "Dim $X"
shortname(::Type{<:Dim{X}}) where X = "$X"
basetype(::Type{<:Dim{X,T,N}}) where {X,T,N} = Dim{X}

"""
    @dim typ [name=string(typ)] [shortname=string(typ)]

Macro to easily define specific dimensions.

Example:
```julia
@dim Lat "Lattitude"
```
"""
macro dim(typ, name=string(typ), shortname=string(typ))
    esc(quote
        struct $typ{T,M,O} <: AbstractDimension{T,M,O}
            val::T
            metadata::M
            order::O
        end
        $typ(val=:; metadata=nothing, order=Forward()) = $typ(val, metadata, order)
        DimensionalData.name(::Type{<:$typ}) = $name
        DimensionalData.shortname(::Type{<:$typ}) = $shortname
    end)
end

# Define some common dimensions
@dim Time
@dim X
@dim Y
@dim Z


text:
#-------------------------------------------------------------------------------
# This file contains utility functions for converting undecorated source
# strings into Julia values.  For example, string->number, string unescaping, etc.

"""
Convert a Julia source code string into a number.
"""
function julia_string_to_number(str::AbstractString, kind)
    str = replace(replace(str, '_'=>""), '−'=>'-')
    if kind == K"Integer"
        x = Base.tryparse(Int, str)
        if Int === Int32 && isnothing(x)
            x = Base.tryparse(Int64, str)
        end
        if isnothing(x)
            x = Base.tryparse(Int128, str)
            if isnothing(x)
                x = Base.parse(BigInt, str)
            end
        end
        return x
    elseif kind == K"Float"
        if !startswith(str,"0x") && 'f' in str
            # This is kind of awful. Should we have a separate Float32 literal
            # type produced by the lexer?  The `f` suffix is nonstandard after all.
            return Base.parse(Float32, replace(str, 'f'=>'e'))
        else
            return Base.parse(Float64, str)
        end
    elseif kind == K"HexInt"
        ndigits = length(str)-2
        return ndigits <= 2  ? Base.parse(UInt8, str)   :
               ndigits <= 4  ? Base.parse(UInt16, str)  :
               ndigits <= 8  ? Base.parse(UInt32, str)  :
               ndigits <= 16 ? Base.parse(UInt64, str)  :
               ndigits <= 32 ? Base.parse(UInt128, str) :
               Base.parse(BigInt, str)
    elseif kind == K"BinInt"
        ndigits = length(str)-2
        return ndigits <= 8   ? Base.parse(UInt8, str)   :
               ndigits <= 16  ? Base.parse(UInt16, str)  :
               ndigits <= 32  ? Base.parse(UInt32, str)  :
               ndigits <= 64  ? Base.parse(UInt64, str)  :
               ndigits <= 128 ? Base.parse(UInt128, str) :
               Base.parse(BigInt, str)
    elseif kind == K"OctInt"
        ndigits = length(str)-2
        x = Base.tryparse(UInt64, str)
        if isnothing(x)
            x = Base.tryparse(UInt128, str)
            if isnothing(x)
                x = Base.parse(BigInt, str)
            elseif ndigits > 43
                x = BigInt(x)
            end
        else
            x = ndigits <= 3  && x <= typemax(UInt8)  ? UInt8(x)   :
                ndigits <= 6  && x <= typemax(UInt16) ? UInt16(x)  :
                ndigits <= 11 && x <= typemax(UInt32) ? UInt32(x)  :
                ndigits <= 22                         ? x          :
                ndigits <= 43                         ? UInt128(x) :
                BigInt(x)
        end
        return x
    end
end

#-------------------------------------------------------------------------------
is_indentation(c) = c == ' ' || c == '\t'

"""
Process Julia source code escape sequences for raw strings
"""
function unescape_raw_string(io::IO, str::AbstractString, is_cmd::Bool)
    delim = is_cmd ? '`' : '"'
    i = firstindex(str)
    lastidx = lastindex(str)
    while i <= lastidx
        c = str[i]
        if c != '\\'
            if c == '\r'
                # convert literal \r and \r\n in strings to \n (issue #11988)
                if i+1 <= lastidx && str[i+1] == '\n'
                    i += 1
                end
                c = '\n'
            end
            write(io, c)
            i = nextind(str, i)
            continue
        end
        # Process \ escape sequences
        j = i
        while j <= lastidx && str[j] == '\\'
            j += 1
        end
        nbackslash = j - i
        if (j <= lastidx && str[j] == delim) || j > lastidx
            # Backslashes before a delimiter must also be escaped
            nbackslash = div(nbackslash,2)
        end
        for k = 1:nbackslash
            write(io, '\\')
        end
        i = j
        if i <= lastidx
            write(io, str[i])
            i = nextind(str, i)
        end
    end
end

"""
Process Julia source code escape sequences for non-raw strings.
`str` should be passed without delimiting quotes.
"""
function unescape_julia_string(io::IO, str::AbstractString)
    i = firstindex(str)
    lastidx = lastindex(str)
    while i <= lastidx
        c = str[i]
        if c != '\\'
            if c == '\r'
                # convert literal \r and \r\n in strings to \n (issue #11988)
                if i+1 <= lastidx && str[i+1] == '\n'
                    i += 1
                end
                c = '\n'
            end
            write(io, c)
            i = nextind(str, i)
            continue
        end
        # Process \ escape sequences.  See also Base.unescape_string which some
        # of this code derives from (but which disallows \` \' \$)
        i += 1
        if i > lastidx
            break
        end
        c = str[i]
        if c == 'x' || c == 'u' || c == 'U'
            n = k = 0
            m = c == 'x' ? 2 :
                c == 'u' ? 4 : 8
            while (k += 1) <= m && i+1 <= lastidx
                nc = str[i+1]
                n = '0' <= nc <= '9' ? n<<4 + (nc-'0') :
                    'a' <= nc <= 'f' ? n<<4 + (nc-'a'+10) :
                    'A' <= nc <= 'F' ? n<<4 + (nc-'A'+10) : break
                i += 1
            end
            if k == 1 || n > 0x10ffff
                u = m == 4 ? 'u' : 'U'
                throw(ArgumentError("invalid $(m == 2 ? "hex (\\x)" :
                                    "unicode (\\$u)") escape sequence"))
            end
            if m == 2 # \x escape sequence
                write(io, UInt8(n))
            else
                print(io, Char(n))
            end
        elseif '0' <= c <= '7'
            k = 1
            n = c-'0'
            while (k += 1) <= 3 && i+1 <= lastidx
                c = str[i+1]
                n = ('0' <= c <= '7') ? n<<3 + c-'0' : break
                i += 1
            end
            if n > 255
                throw(ArgumentError("octal escape sequence out of range"))
            end
            write(io, UInt8(n))
        else
            u = # C escapes
                c == 'n' ? '\n' :
                c == 't' ? '\t' :
                c == 'r' ? '\r' :
                c == 'e' ? '\e' :
                c == 'b' ? '\b' :
                c == 'f' ? '\f' :
                c == 'v' ? '\v' :
                c == 'a' ? '\a' :
                # Literal escapes allowed in Julia source
                c == '\\' ? '\\' :
                c == '\'' ? '\'' :
                c == '"' ? '"' :
                c == '$' ? '$' :
                c == '`' ? '`' :
                throw(ArgumentError("Invalid escape sequence \\$c"))
            write(io, u)
        end
        i = nextind(str, i)
    end
end

function unescape_julia_string(str::AbstractString, is_cmd::Bool, is_raw::Bool)
    io = IOBuffer()
    if is_raw
        unescape_raw_string(io, str, is_cmd)
    else
        unescape_julia_string(io, str)
    end
    String(take!(io))
end

#-------------------------------------------------------------------------------
# Unicode normalization. As of Julia 1.8, this is part of Base and the Unicode
# stdlib under the name `Unicode.julia_chartransform`. See
# https://github.com/JuliaLang/julia/pull/42561
#
# To allow use on older Julia versions, we reproduce that logic here.

# static wrapper around user callback function
utf8proc_custom_func(codepoint::UInt32, callback::Any) =
    UInt32(callback(codepoint))::UInt32

function utf8proc_decompose(str, options, buffer, nwords, chartransform::T) where T
    ret = ccall(:utf8proc_decompose_custom, Int, (Ptr{UInt8}, Int, Ptr{UInt8}, Int, Cint, Ptr{Cvoid}, Ref{T}),
                str, sizeof(str), buffer, nwords, options,
                @cfunction(utf8proc_custom_func, UInt32, (UInt32, Ref{T})), chartransform)
    ret < 0 && utf8proc_error(ret)
    return ret
end

function utf8proc_map(str::Union{String,SubString{String}}, options::Integer, chartransform=identity)
    nwords = utf8proc_decompose(str, options, C_NULL, 0, chartransform)
    buffer = Base.StringVector(nwords*4)
    nwords = utf8proc_decompose(str, options, buffer, nwords, chartransform)
    nbytes = ccall(:utf8proc_reencode, Int, (Ptr{UInt8}, Int, Cint), buffer, nwords, options)
    nbytes < 0 && utf8proc_error(nbytes)
    return String(resize!(buffer, nbytes))
end

const _julia_charmap = Dict{UInt32,UInt32}(
    0x025B => 0x03B5,
    0x00B5 => 0x03BC,
    0x00B7 => 0x22C5,
    0x0387 => 0x22C5,
    0x2212 => 0x002D,
)

julia_chartransform(codepoint::UInt32) = get(_julia_charmap, codepoint, codepoint)

function normalize_identifier(str)
    flags = Base.Unicode.UTF8PROC_STABLE | Base.Unicode.UTF8PROC_COMPOSE
    utf8proc_map(str, flags, julia_chartransform)
end

ext: jl


text:
@testset "Visualization" begin
    air_passengers = CSV.read(StateSpaceModels.AIR_PASSENGERS, DataFrame)
    log_air_passengers = log.(air_passengers.passengers)
    model = BasicStructural(log_air_passengers, 12)
    fit!(model)
    # forecasting
    forec = forecast(model, 12)
    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), model, forec)
    @test length(rec) == 4
end


text:
using Lerche
using Test

include("trees.jl")
include("test_parser.jl")


text:
# convenience function to calculate the mean-squared error
function mse(arr1::AbstractArray, arr2::AbstractArray)
    if size(arr1) != size(arr2)
        throw(DimensionMismatch("Got $(size(arr1)) and $(size(arr2))"))
    end
    N = length(arr1)
    err = 0.0
    for i in 1:N
        err += (arr2[i] - arr1[i])^2
    end
    err /= N
end


text:
# Fetch wfm
function fetch_wfm(obj::T;ch::String="1") where {T<:NIScope}
	# niScope channels start from 0, so subtract 1
	ch = NIScope_Vertical_Dict[ch]
	y, aux = niScope_Fetch(obj;ch=ch)
	return y, aux.absoluteInitialX, aux.xIncrement
end

# one channel at a time...
function niScope_Fetch(obj::T;ch::String="1") where {T<:NIScope}
#ViStatus niScope_Fetch (ViSession vi, ViConstString channelList, ViReal64 timeout, ViInt32 numSamples, ViReal64* wfm, struct niScope_wfmInfo* wfmInfo);
	chlist  = [UInt8.(collect("$ch")); UInt8(0)] # terminate with NULL char
	timeout = 0.0
	n_samp  = niScope_ActualRecordLength(obj)
	wfmInfo = [niScope_wfmInfo(0.0, 0.0, 0.0, Int32(0), 1.0, 0.0, 0, 0)] # this needs to be an array of structures, as per niScope documentation
	y       = Vector{Cdouble}(undef, n_samp)# arg type is Ref but one needs to pass the actual array
	status  = niScope_Fetch(obj.handle, chlist, timeout, n_samp, y, wfmInfo)
	return y, wfmInfo[1]# return the first wfm in array
end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_bbc_m4_opt_rho0_69(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_bbc_m4_opt_rho0_69!(A_copy)
end

@inline function exp_bbc_m4_opt_rho0_69!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=7
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 T2k7
    # Computing Bb2 = x*I+x*A
    coeff1=0.08830117968963386
    coeff2=1.6185302637301089
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.08830117968963386
    coeff2=1.6185302637301089
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=-0.012459653416356822
    coeff2=0.10918913720213257
    coeff3=1.0031491738964604
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=0.1634190055195114
    coeff2=0.9966322449217647
    coeff3=0.02185426146644069
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=0.07611513553783929
    coeff2=-0.1833134793841606
    coeff3=-0.005277078018897659
    coeff4=-0.0019621792320780053
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=0.07611513553783929
    coeff2=-0.1833134793841606
    coeff3=-0.005277078018897659
    coeff4=-0.0019621792320780053
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=5.151267490244333
    coeff2=0.3931553522425758
    coeff3=-0.0022244347351211098
    coeff4=0.000393833045141096
    coeff5=0.2966584239338796
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.14452058311618007
    coeff2=0.15667930242227487
    coeff3=0.20054402254365047
    coeff4=0.015735447774198674
    coeff5=1.4175497249260038
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing T2k7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-0.11353832118393399
    coeff2=0.06878048308436636
    coeff3=-1.4102780317055799
    coeff4=-0.19117716256815653
    coeff5=0.39486010387414244
    coeff6=1.4114856620223937
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    return memslots1 # Returning T2k7
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
using MCMCChains: promote_eltype_tuple_type
using Dates,JLD2,Statistics,Optim,Parameters,Distributions,DataFrames,CSV,MCMCChains
import KenyaCoVSD
include("fitting_methods.jl");
@load("forecasts/condensed_county_forecasts.jld2")
## Data S1
## Get all parameter names
modelfitfilenames = readdir("modelfits",join = true)
modeldict = load(modelfitfilenames[1])
model = modeldict[first(keys(modeldict))]
chn = model.MCMC_results.chain
parameternames = names(chn)

##Fill in a matrix of strings with the posterior mean and 95% CIs for each parameter
q = quantile(chn)
m = mean(chn)

parameterfits = fill("",47,length(parameternames))
mean_parameterfits = fill(0.0,47,length(parameternames))

for (i,filename) in enumerate(modelfitfilenames)
    modeldict = load(filename)
    model = modeldict[first(keys(modeldict))]
    chn = model.MCMC_results.chain
    parameternames = names(chn)
    q = quantile(chn)
    m = mean(chn)
    for j = 1:length(parameternames)
        str = "$(round(m.nt.mean[j],digits = 3)) ($(round(q.nt[2][j],digits = 3)),$(round(q.nt[end][j],digits = 3)))"
        parameterfits[i,j] = str
        mean_parameterfits[i,j] = m.nt.mean[j]
    end
end

##Fill in a matrix of IFR AND detection %s from the condensed fitting information
fatality_detection_percent = zeros(47,2)
for i = 1:47,j = 1:2
    fit = condensed_county_forecasts[i]
    fatality_detection_percent[i,j] = round(fit.pred.μ_fit[j]*100,sigdigits = 3)
end

##Get the population exposure estimates
@load("forecasts/county_pop_exposure_estimates_firstjune2021.jld2")
population_exposure_fits = fill("",47)
for i = 1:47
    fit = population_exposure_estimates[i]
    mean_est = round(fit.pred_pop_exposure*100,digits = 1)
    lest = round((fit.pred_pop_exposure - fit.lerr_pop_exposure)*100,digits = 1)
    uest = round((fit.pred_pop_exposure + fit.uerr_pop_exposure)*100,digits = 1)
    str = "$(mean_est) ($(lest),$(uest))"
    population_exposure_fits[i] = str
end

##Combine fits into one document
countynames = [fit.name for fit in condensed_county_forecasts]
data = [countynames parameterfits string.(fatality_detection_percent) population_exposure_fits]
_data = [countynames mean_parameterfits]
colnames = [:county_name;parameternames;[:fatality_detection_lower_SES,:fatality_detection_upper_SES,:population_exposure]]
_colnames = [:county_name;parameternames]
dataS1 = DataFrame()
parameter_post_means = DataFrame()
for (j,colname) in enumerate(colnames)
    dataS1[!,string(colname)] = data[:,j]
end

for (j,colname) in enumerate(_colnames)
    parameter_post_means[!,string(colname)] = _data[:,j]
end

CSV.write("opendatacsvs/dataS1.csv",dataS1)
CSV.write("forecasts/parameter_posterior_means.csv",parameter_post_means)

## Data S2, S3 and S4

dataS2 = DataFrame()
dataS3 = DataFrame()
dataS4 = DataFrame()

for (i,filename) in enumerate(modelfitfilenames)
    modeldict = load(filename)
    model = modeldict[first(keys(modeldict))]
    if i == 1
        dataS2[!,"Date"] = [Date(2020,2,20) + Day(k) for k = 1:size(model.PCR_cases,1)]
        dataS3[!,"Date"] = [Date(2020,2,20) + Day(k) for k = 1:size(model.sero_cases,1)]
        dataS4[!,"Date"] = [Date(2020,2,20) + Day(k) for k = 1:size(model.deaths,1)]
    end
    dataS2[!,model.areaname*"_PCR_pos"] = model.PCR_cases[:,1]
    dataS2[!,model.areaname*"_PCR_neg"] = model.PCR_cases[:,2]
    dataS3[!,model.areaname*"_sero_pos"] = model.sero_cases[:,1,1]
    dataS3[!,model.areaname*"_sero_neg"] = model.sero_cases[:,2,1]
    dataS4[!,model.areaname*"_deaths_pos"] = model.deaths
end
CSV.write("opendatacsvs/dataS2.csv",dataS2)
CSV.write("opendatacsvs/dataS3.csv",dataS3)
CSV.write("opendatacsvs/dataS4.csv",dataS4)


text:
# ---
# title: 502. IPO
# id: problem502
# author: Indigo
# date: 2021-09-08
# difficulty: Hard
# categories: Heap, Greedy
# link: <https://leetcode.com/problems/ipo/description/>
# hidden: true
# ---
#
# Suppose LeetCode will start its IPO soon. In order to sell a good price of its
# shares to Venture Capital, LeetCode would like to work on some projects to
# increase its capital before the IPO. Since it has limited resources, it can
# only finish at most **k** distinct projects before the IPO. Help LeetCode
# design the best way to maximize its total capital after finishing at most
# **k** distinct projects.
#
# You are given several projects. For each project **i** , it has a pure profit
# **P i** and a minimum capital of **C i** is needed to start the corresponding
# project. Initially, you have **W** capital. When you finish a project, you
# will obtain its pure profit and the profit will be added to your total
# capital.
#
# To sum up, pick a list of at most **k** distinct projects from given projects
# to maximize your final capital, and output your final maximized capital.
#
# **Example 1:**
#
#
#
#     Input: k=2, W=0, Profits=[1,2,3], Capital=[0,1,1].
#
#     Output: 4
#
#     Explanation: Since your initial capital is 0, you can only start the project indexed 0.
#                  After finishing it you will obtain profit 1 and your capital becomes 1.
#                  With capital 1, you can either start the project indexed 1 or the project indexed 2.
#                  Since you can choose at most 2 projects, you need to finish the project indexed 2 to get the maximum capital.
#                  Therefore, output the final maximized capital, which is 0 + 1 + 3 = 4.
#
#
# **Note:**
#
#   1. You may assume all numbers in the input are non-negative integers.
#   2. The length of Profits array and Capital array will not exceed 50,000.
#   3. The answer is guaranteed to fit in a 32-bit signed integer.
#
#
## @lc code=start
using LeetCode

function find_maximized_capital(k::Int, w::Int, profits::Vector{Int}, capital::Vector{Int})
    arr = @view(collect(zip(capital, profits))[:])
    sort!(arr)
    hp = Int[]
    for _ in 1:k
        idx = searchsortedlast(arr, (w + 1, -1))
        for item in @view(arr[1:idx])
            heappush!(hp, -item[2])
        end
        arr = @view(arr[idx+1:end])
        w -= isempty(hp) ? 0 : hp[1]
    end
    w
end
## @lc code=end


text:
"""
FaultTree
"""

export AbstractFaultTreeNode, AbstractFaultTreeOperation, FaultTreeOperation, FaultTreeEvent, FaultTreeKoutofN
export ftevent

abstract type AbstractFaultTreeNode end
abstract type AbstractFaultTreeOperation <: AbstractFaultTreeNode end

struct FaultTreeOperation <: AbstractFaultTreeOperation
    params::Set{Symbol}
    op::Symbol
    args::Vector{AbstractFaultTreeNode}
end

struct FaultTreeKoutofN <: AbstractFaultTreeOperation
    params::Set{Symbol}
    op::Symbol
    k::Int
    args::Vector{AbstractFaultTreeNode}
end

struct FaultTreeEvent <: AbstractFaultTreeNode
    params::Set{Symbol}
    var::Symbol
end

"""
ftree
"""

function ftevent(var::Symbol)
    FaultTreeEvent(Set([var]), var)
end

function ftevent(var::Vararg)
    s = Symbol(var...)
    FaultTreeEvent(Set([s]), s)
end


text:
# ------------------------------------------------------------------------------------------
# # Entropies & Dimensions
#
#
# Topics:
# * Entropies from **DynamicalSystems.jl**
# * Generalized Dimension
# * Automated dimension estimation!
# * Other related concepts
# * Docstrings
#
# ---
#
# # Generalized Entropy
#
# * In the study of dynamical systems there are many quantities that identify as "entropy".
# * These quantities are not the more commonly known [thermodynamic
# ones](https://en.wikipedia.org/wiki/Entropy), used in Statistical Physics.
# * Rather, they are more like the entropies of [information
# theory](https://en.wikipedia.org/wiki/Entropy_(information_theory), which represent
# information contained within a dataset.
# * In general, the more "uncertain" or "random" the dataset is, the larger its entropy will
# be. On the other hand, the lower the entropy, the more "predictable" the dataset becomes.
#
#
# Let $p$ be an array of probabilities (such that it sums to 1). Then the generalized
# entropy is defined as
#
# $$
# R_\alpha(p) = \frac{1}{1-\alpha}\log\left(\sum_i p[i]^\alpha\right)
# $$
#
# and is also called [Rényi entropy](https://en.wikipedia.org/wiki/R%C3%A9nyi_entropy).
# Other entropies, like e.g. the [Shannon
# entropy](https://en.wikipedia.org/wiki/Entropy_(information_theory) are generalized by it,
# since for $\alpha = 1$, the Rényi entropy becomes the Shannon entropy,
#
# $$
# R_1(p) = -\left(\sum_i p[i] \log (p[i]) \right)
# $$
#
# The Rényi entropy can be computed for a specific dataset, given $p$. But how does one get
# $p$?
# 1. $p$ represents the probability that a point of a dataset falls into a specific "bin".
# 2. It is nothing more than the (normalized) histogram of the dataset!
# ------------------------------------------------------------------------------------------

using DynamicalSystems

# ------------------------------------------------------------------------------------------
# Let's generate a dataset so that we can practice calculating entropies.
# ------------------------------------------------------------------------------------------

N = 100000
randomdata = Dataset(rand(N,3))

# ------------------------------------------------------------------------------------------
# ---
#
# ```julia
# genentropy(α, ε, dataset::AbstractDataset; base = e)
# ```
# * This function calculates the generalized entropy of order `α`.
# * It first calculates the probability array $p$.
# * The "histogram" is created by partitioning the `dataset` into boxes of size `ε`.
#
# ------------------------------------------------------------------------------------------

genentropy(2, 0.1, randomdata)

genentropy(2, 0.01, randomdata)

genentropy(2, 0.001, randomdata)

genentropy(2, 0.0001, randomdata)

# ------------------------------------------------------------------------------------------
# Note that the output of `genentropy` changed with changing $\varepsilon$ until we hit
# $\varepsilon = 0.001$.
#
# At this point the value for the entropy has already saturated. There's no use in
# partitioning the dataset in smaller boxes.
#
# ---
#
# Now let's calculate the entropy of a coin toss!
#
# First, let's create an array, `y`, that stores the results of coin tosses as `0`s or `1`s
# for `1000000` tosses.
# ------------------------------------------------------------------------------------------

y = Float64.(rand(Bool, 1000000))

sh = genentropy(1, 0.1, y) # Renyi entropy with α = 1 is the Shannon entropy

# ------------------------------------------------------------------------------------------
# The above number should be log(2) [by
# definition](https://en.wikipedia.org/wiki/Shannon_(unit)
# ------------------------------------------------------------------------------------------

isapprox(sh, log(2), rtol = 1e-4)

# ------------------------------------------------------------------------------------------
# `genentropy` is conveniently used with `trajectory` outputs.
#
# Here we create a trajectory for a towel map,
# ------------------------------------------------------------------------------------------

towel = Systems.towel()
tr = trajectory(towel, N-1);
summary(tr)

# ------------------------------------------------------------------------------------------
# and calculate its entropy:
# ------------------------------------------------------------------------------------------

genentropy(1, 0.01, tr) # The result is with log base-e !

# ------------------------------------------------------------------------------------------
# Let's also compare the entropy of the above dataset (a trajectory of the towel map) with
# that of a random dataset:
# ------------------------------------------------------------------------------------------

genentropy(1, 0.01, randomdata)

# ------------------------------------------------------------------------------------------
# * As expected, the entropy of the random dataset is higher.
#
# ---
#
# How much time does the computation take?
# ------------------------------------------------------------------------------------------

using BenchmarkTools
@btime genentropy(1, 0.01, $tr);

# ------------------------------------------------------------------------------------------
# ## Specialized histogram
# * Partitioning the dataset (i.e. generating a "histogram") is in general a costly
# operation that depends exponentially on the number of dimensions.
# * In this specific application however, we can tremendously reduce the memory allocation
# and time spent!
#
# Let's get the array of probabilities $p$ for size ε where `tr` is the trajectory of the
# towel map
# ------------------------------------------------------------------------------------------

ε = 0.01
p = non0hist(ε, tr)

# ------------------------------------------------------------------------------------------
# Here's a sanity check, showing our probabilities should sum to roughly `1`.
# ------------------------------------------------------------------------------------------

sum(p)

# ------------------------------------------------------------------------------------------
# How long does generating the histogram take?
# ------------------------------------------------------------------------------------------

@btime non0hist($ε, $tr);

# ------------------------------------------------------------------------------------------
# How long does this take if we create 9-dimensional data and compare again?
# ------------------------------------------------------------------------------------------

nine = Dataset(rand(N, 9))
@btime non0hist($ε, $nine);

# ------------------------------------------------------------------------------------------
# * We went from dimension 3 to dimension 9 but the time roughly only tripled
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Attractor Dimension
# 1. There are numerous methods that one can use to calculate a so-called "dimension" of a
# dataset, like for example the [Fractal
# dimension](https://en.wikipedia.org/wiki/Fractal_dimension).
#
# 2. Most of the time these dimensions indicate some kind of scaling behavior.
#
# 3. For example, the scaling of `genentropy` with decreasing `ε` gives the so-called
# "generalized dimension".
#
#
# $ E \approx -D\log(\varepsilon)$ with $E$ the entropy and $D$ the "dimension".
#
# ---
# I want to know dimension of attractor of the Towel Map!
#
#
# ------------------------------------------------------------------------------------------

towel = Systems.towel()
towel_tr = trajectory(towel, 1000000);
summary(towel_tr)

# ------------------------------------------------------------------------------------------
# Note that more points = more precision = more computations = more time!
#
# Now I want to compute `genentropy` for different ε.
#
# Which ε should we use...?
#
# Let's do a "random" logspace based guess...
# ------------------------------------------------------------------------------------------

ες = logspace(-4, 1, 12)

Es = zeros(ες)
for (i, ε) ∈ enumerate(ες)
    Es[i] = genentropy(1, ε, towel_tr)
end
Es

# ------------------------------------------------------------------------------------------
# **Shorter version:**
# ------------------------------------------------------------------------------------------

Es = genentropy.(1, ες, towel_tr)

# ------------------------------------------------------------------------------------------
# Alright. Remember that it should be that $E \approx -D\log(\varepsilon)$
#  with $E$ the entropy and $D$ the "dimension".
#
# Let's plot and see:
# ------------------------------------------------------------------------------------------

using PyPlot; figure(figsize = (6,4))
x = -log.(ες)
plot(x, Es); xlabel("-log(ε)"); ylabel("Entropy");
plot([x[4], x[4]], [0, 15], color = "C1", alpha = 0.5)
plot([x[end-3], x[end-3]], [0, 15], color = "C1", alpha = 0.5);

# ------------------------------------------------------------------------------------------
# What typically happens is that there is some region where this scaling behavior holds, but
# then it stops holding due to the finite amount of data points.
#
# Above, the expected scaling behavior holds between the orange vertical lines.
#
# Let's choose the curve points that do fall in the linear regime of the above plot,
# ------------------------------------------------------------------------------------------

x, y = -log.(ες)[4:end-2], Es[4:end-2]

# ------------------------------------------------------------------------------------------
# and find the slope of the curve there, to calculate the dimension, D.
# ------------------------------------------------------------------------------------------

offset, slope = linreg(x, y)
D = slope

# ------------------------------------------------------------------------------------------
# This is actually a correct result, the information dimension of the attractor of the towel
# map is around 2.
#
# ---
#
# * Are the values of `ες` we used good?
# * For a general dataset, how can we determine them?
#
# the function `estimate_boxsizes(dataset; kwargs...)` can help with that!
# ------------------------------------------------------------------------------------------

ες = estimate_boxsizes(towel_tr)

# ------------------------------------------------------------------------------------------
# Let's plot $E$ vs. $-\log \epsilon$ again
# ------------------------------------------------------------------------------------------

Es = genentropy.(1, ες, towel_tr)
figure(figsize = (6,4))
plot(-log.(ες), Es); xlabel("-log(ε)"); ylabel("E");

# ------------------------------------------------------------------------------------------
# ---
# # Automated Dimension Estimation
#
# Given some arbitrary plot like the one above, is there any algorithm to deduce a scaling
# region???
#
# The function `linear_regions(x, y; kwargs...)` decomposes the function `y(x)` into regions
# where  the function is linear.
#
# It returns the indices of `x` that correspond to linear regions and the approximated
# tangents at each region!
# ------------------------------------------------------------------------------------------

x = -log.(ες)
lrs, slopes = linear_regions(x, Es)

for i in 1:length(slopes)
    println("linear region $(i) starts from index $(lrs[i]) and ends at index $(lrs[i+1])")
    println("with corresponding slope $(slopes[i])")
    println()
end

# ------------------------------------------------------------------------------------------
# The linear region which is biggest is "probably correct one".
# Here the last linear region is the largest; thus the slope is
# ------------------------------------------------------------------------------------------

slopes[end]

# ------------------------------------------------------------------------------------------
# This `linear_regions` function seems very shady... Is there any "easy" way to visualize
# what it does? *Say no more!*
#
# In this next example, we'll use `ChaosTools`, which belows to `DynamicalSystems`.
# ------------------------------------------------------------------------------------------

using PyPlot, ChaosTools

# ------------------------------------------------------------------------------------------
# `ChaosTools` provides a function, `plot_linear_regions` !
# ------------------------------------------------------------------------------------------

figure(figsize=(6,4))
plot_linear_regions(x, Es)
xlabel("-log(ε)"); ylabel("E");

# ------------------------------------------------------------------------------------------
# Adjust the tolerance of `linear_regions` using keyword argument `tol`
# ------------------------------------------------------------------------------------------

figure(figsize=(6,4))
plot_linear_regions(x, Es; tol = 0.8)
xlabel("-log(ε)"); ylabel("E");

# ------------------------------------------------------------------------------------------
# Notice how the color schemes allow us to visualize that adjusting the tolerance reduces
# our number of linear regimes!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## `generalized_dim` function
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Let's summarize what we just did to estimate the dimension of an attractor.
#
# 1. We decided on some partition sizes `ες` to use (the function `estimate_boxsizes` can
# give an estimate for that).
# 2. For each `ε` in `ες` we calculated the entropy via `genentropy`. We stored these
# entropies in an array `Es`.
# 3. We tried to find a "linear scaling region" of the curve `Es` vs. `-log.(ες)`.
# 4. The slope of this "linear scaling region" is the dimension we estimated.
#
# Wouldn't it be **cool** if all of this process could happen with one function call?
#
# This is *exactly* what the following function does:
# ```julia
# generalized_dim(α, dataset, ες = estimate_boxsizes(tr))
# ```
# which computes the `α`-order generalized dimension.
# ------------------------------------------------------------------------------------------

generalized_dim(2.0, tr)

generalized_dim(1.0, tr)

# ------------------------------------------------------------------------------------------
# The first input to `generalized_dim` is our value for $\alpha$!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Similarly, let's calculate the dimension of the Henon map that we have seen in previous
# tutorials,
# ------------------------------------------------------------------------------------------

hen = Systems.henon()
tr = trajectory(hen, 200000)
generalized_dim(0, tr)

# ------------------------------------------------------------------------------------------
# and of the Lorenz system:
# ------------------------------------------------------------------------------------------

lor = Systems.lorenz()
tr_lor = trajectory(lor, 1000.0; dt = 0.05);

generalized_dim(2.0, tr_lor)

# ------------------------------------------------------------------------------------------
# The correlation dimension of the Lorenz attractor (for default parameters) is reported
# *somewhere* around `2.0` (Grassberger and Procaccia, 1983).
#
# ## `generalized_dim` is but a crude estimate!
#
# **It is important to understand that `generalized_dim` is only a crude estimate! You must
# check and double-check and triple-check if you want more accuracy!**
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Other related concepts
# ## Kaplan-Yorke dimension
# The Kaplan-Yorke dimension is defined as simply the (interpolated) number where the sum of
# Lyapunov exponents crosses zero.
#
# This simple interpolation is done by the function `kaplanyorke_dim(ls)`.
# ------------------------------------------------------------------------------------------

lor = Systems.lorenz()

ls = lyapunovs(lor, 4000.0; Ttr = 100.0)

kaplanyorke_dim(ls)

# ------------------------------------------------------------------------------------------
# ## Permutation Entropy
#
# Another entropy-like quantity that you can compute with **DynamicalSystems.jl** is the
# [permutation entropy](https://juliadynamics.github.io/DynamicalSystems.jl/latest/chaos/ent
# ropies/#permutation-entropy).
#
# This is done with the function `permentropy(s, order; ...)`. The permutation entropy is
# different because it requires a *timeseries* `s` as the input.
# ------------------------------------------------------------------------------------------

# create timeseries:
s = (ds = Systems.towel(); trajectory(towel, 10000)[:, 1])
# order `o` permutation entropy:
o = 6
permentropy(s, o)

# ------------------------------------------------------------------------------------------
# # Docstrings
# ------------------------------------------------------------------------------------------

?genentropy

?generalized_dim

?non0hist

?estimate_boxsizes

?linear_regions

?kaplanyorke_dim

?permentropy

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
"""
    pspole(psys::PeriodicStateSpace{PM}[,K]; kwargs...) -> val

Return for the periodic system `psys = (A(t),B(t),C(t),D(t))` the complex vector `val` containing
the characteristic exponents of the periodic matrix `A(t)` (also called the _poles_ of the system `psys`).

Depending on the underlying periodic matrix type `PM`, the optional argument `K` and keyword arguments `kwargs` may have the following values:

- if `PM = PeriodicFunctionMatrix`, or `PM = PeriodicSymbolicMatrix`, or `PM = PeriodicTimeSeriesMatrix`, then `K` is the number of factors used to express the monodromy matrix of `A(t)` (default: `K = 1`)  and `kwargs` are the keyword arguments of  [`pseig(::PeriodicFunctionMatrix)`](@ref);

- if `PM = HarmonicArray`, then `K` is the number of harmonic components used to represent the Fourier series of `A(t)` (default: `K = max(10,nh-1)`, `nh` is the number of harmonics terms of `A(t)`)  and `kwargs` are the keyword arguments of  [`psceig(::HarmonicArray)`](@ref);

- if `PM = FourierFunctionMatrix`, then `K` is the number of harmonic components used to represent the Fourier series of `A(t)` (default: `K = max(10,nh-1)`, `nh` is the number of harmonics terms of `A(t)`)  and `kwargs` are the keyword arguments of  [`psceig(::FourierFunctionMatrix)`](@ref);

- if `PM = PeriodicMatrix` or `PM = PeriodicArray`, then `K` is the starting sample time (default: `K = 1`)  and `kwargs` are the keyword arguments of  [`psceig(::PeriodicMatrix)`](@ref);
"""
function pspole(psys::PeriodicStateSpace, K::Int; kwargs...)
    return psceig(psys.A,K; kwargs...)
end
pspole(psys::PeriodicStateSpace; kwargs...) = psceig(psys.A; kwargs...)
"""
    pszero(psys::PeriodicStateSpace{HarmonicArray}[, N]; P, atol, rtol, fast) -> val
    pszero(psys::PeriodicStateSpace{PeriodicFunctionMatrix}[, N]; P, atol, rtol, fast) -> val
    pszero(psys::PeriodicStateSpace{PeriodicSymbolicMatrix}[, N]; P, atol, rtol, fast) -> val

Compute the finite and infinite zeros of a continuous-time periodic system `psys = (A(t), B(t), C(t), D(t))` in `val`,
where the periodic system matrices `A(t)`, `B(t)`, `C(t)`, and `D(t)` are in a _Harmonic Array_, or
_Periodic Function Matrix_, or _Periodic Symbolic Matrix_ representation
(the last two representation are automatically converted to a _Harmonic Array_ representation).
`N` is the number of selected harmonic components in the Fourier series of the system matrices (default: `N = max(20,nh-1)`,
where `nh` is the maximum number of harmonics terms) and the keyword parameter `P` is the number of full periods
to be considered (default: `P = 1`) to build
a frequency-lifted LTI representation based on truncated block Toeplitz matrices (see [`ps2fls`](@ref)).

The computation of the zeros of the _complex_ lifted system is performed by reducing the corresponding system pencil
to an appropriate Kronecker-like form which exhibits the finite and infinite eigenvalues.
The reduction is performed using orthonal similarity transformations and involves rank decisions based on rank revealing QR-decompositions with column pivoting,
if `fast = true`, or, the more reliable, SVD-decompositions, if `fast = false`. For a system `psys` of period `T`,
the finite zeros are determined as those eigenvalues which have imaginary parts in the interval `[-ω/2, ω/2]`, where `ω = 2π/(P*T)`.
To eliminate possible spurious finite eigenvalues, the intersection of two finite eigenvalue sets is computed
for two lifted systems obtained for `N` and `N+2` harmonic components.
The infinite zeros are determined as the infinite zeros of the LTI system `(A(ti), B(ti), C(ti), D(ti))`
resulting for a random time value `ti`. _Warning:_ While this evaluation of the number of infinite zeros mostly
provides the correct result, there is no theoretical assessment of this approach (counterexamples are welcome!).

The keyword arguments `atol`  and `rtol` specify the absolute and relative tolerances for the nonzero
elements of the underlying lifted system pencil, respectively.
The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of the pencil, and `ϵ` is the
working machine epsilon.
"""
function pszero(psys::PeriodicStateSpace{<: HarmonicArray}, N::Union{Int,Missing} = missing; P::Int= 1, fast::Bool = true, atol::Real = 0, rtol::Real = 0)
    ismissing(N) && (N = max(20, max(size(psys.A.values,3),size(psys.B.values,3),size(psys.C.values,3),size(psys.D.values,3))-1))
    (N == 0 || islti(psys) ) && (return spzeros(dssdata(psaverage(psys))...; fast, atol1 = atol, atol2 = atol, rtol)[1])
    period = psys.A.period
    ωhp2 = pi/P/period
    n = size(psys.A,1)
    T = promote_type(Float64, eltype(psys.A))
    # employ heuristics to determine fix finite zeros by comparing two sets of computed zeros
    z = spzeros(dssdata(ps2fls(psys, N; P))...; fast, atol1 = atol, atol2 = atol, rtol)[1]
    zf = z[isfinite.(z)]
    ind = sortperm(imag(zf),by=abs);
    nf = count(abs.(imag(zf[ind[1:min(4*n,length(ind))]])) .<=  ωhp2*(1+sqrt(eps(T))))
    zf = zf[ind[1:nf]]
    z2 = spzeros(dssdata(ps2fls(psys, N+2; P))...; fast, atol1 = atol, atol2 = atol, rtol)[1]
    zf2 = z2[isfinite.(z2)]
    ind = sortperm(imag(zf2),by=abs);
    nf2 = count(abs.(imag(zf2[ind[1:min(4*n,length(ind))]])) .<=  ωhp2*(1+sqrt(eps(T))))
    zf2 = zf2[ind[1:nf2]]
    σf = Complex{T}[]
    nf < nf2 || ((zf, zf2) = (zf2, zf))
    atol > 0 || (norms = max(norm(psys.A.values,Inf),norm(psys.B.values,Inf),norm(psys.C.values,Inf),norm(psys.D.values,Inf)))
    tol = atol > 0 ? atol : (rtol > 0 ? rtol*norms : sqrt(eps(T))*norms)
    for i = 1:min(nf,nf2)
        minimum(abs.(zf2 .- zf[i])) < tol  && push!(σf,zf[i])
    end
    isreal(σf) && (σf = real(σf))
    if any(isinf.(z))
       # Conjecture: The number of infinite zeros is the same as that of the time-evaluated system!
       zm = spzeros(dssdata(psteval(psys, period*rand()))...; fast, atol1 = atol, atol2 = atol, rtol)[1]
       return [σf; zm[isinf.(zm)]]
    else
       return σf
    end
end
"""
    pszero(psys::PeriodicStateSpece{FourierFunctionMatrix}[, N]; P, atol, rtol, fast) -> val

Compute the finite and infinite zeros of a continuous-time periodic system `psys = (Af(t), Bf(t), Cf(t), Df(t))` in `val`,
where the periodic system matrices `Af(t)`, `Bf(t)`, `Cf(t)`, and `Df(t)` are in a _Fourier Function Matrix_ representation.
`N` is the number of selected harmonic components in the Fourier series of the system matrices (default: `N = max(20,nh-1)`,
where `nh` is the maximum number of harmonics terms) and the keyword parameter `P` is the number of full periods
to be considered (default: `P = 1`) to build
a frequency-lifted LTI representation based on truncated block Toeplitz matrices (see [`ps2frls`](@ref)).

The computation of the zeros of the _real_ lifted system is performed by reducing the corresponding system pencil
to an appropriate Kronecker-like form which exhibits the finite and infinite eigenvalues.
The reduction is performed using orthonal similarity transformations and involves rank decisions based on rank revealing QR-decompositions with column pivoting,
if `fast = true`, or, the more reliable, SVD-decompositions, if `fast = false`. For a system `psys` of period `T`,
the finite zeros are determined as those eigenvalues which have imaginary parts in the interval `[-ω/2, ω/2]`, where `ω = 2π/(P*T)`.
To eliminate possible spurious finite eigenvalues, the intersection of two finite eigenvalue sets is computed
for two lifted systems obtained for `N` and `N+2` harmonic components.
The infinite zeros are determined as the infinite zeros of the LTI system `(Af(ti), Bf(ti), Cf(ti), Df(ti))`
resulting for a random time value `ti`. _Warning:_ While this evaluation of the number of infinite zeros mostly
provides the correct result, there is no theoretical assessment of this approach (counterexamples are welcome!).

The keyword arguments `atol`  and `rtol` specify the absolute and relative tolerances for the nonzero
elements of the underlying lifted system pencil, respectively.
The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of the pencil, and `ϵ` is the
working machine epsilon.
"""
function pszero(psys::PeriodicStateSpace{<: FourierFunctionMatrix}, N::Union{Int,Missing} = missing; P::Int= 1, fast::Bool = true, atol::Real = 0, rtol::Real = 0)
    ismissing(N) && (N = max(20, maximum(ncoefficients.(Matrix(psys.A.M))), maximum(ncoefficients.(Matrix(psys.B.M))),
                                   maximum(ncoefficients.(Matrix(psys.C.M))), maximum(ncoefficients.(Matrix(psys.A.M)))))
    (N == 0 || islti(psys) ) && (return spzeros(dssdata(psaverage(psys))...; fast, atol1 = atol, atol2 = atol, rtol)[1])

    # employ heuristics to determine fix finite zeros by comparing two sets of computed zeros
    z = spzeros(dssdata(ps2frls(psys, N; P))...; fast, atol1 = atol, atol2 = atol, rtol)[1]

    period = psys.A.period
    ωhp2 = pi/P/period
    n = size(psys.A,1)
    T = promote_type(Float64, eltype(psys.A))
    zf = z[isfinite.(z)]
    ind = sortperm(imag(zf),by=abs);
    nf = count(abs.(imag(zf[ind[1:min(4*n,length(ind))]])) .<=  ωhp2*(1+sqrt(eps(T))))
    zf = zf[ind[1:nf]]

    z2 = spzeros(dssdata(ps2frls(psys, N+2; P))...; fast, atol1 = atol, atol2 = atol, rtol)[1]
    zf2 = z2[isfinite.(z2)]
    ind = sortperm(imag(zf2),by=abs);
    nf2 = count(abs.(imag(zf2[ind[1:min(4*n,length(ind))]])) .<=  ωhp2*(1+sqrt(eps(T))))
    zf2 = zf2[ind[1:nf2]]
    σf = Complex{T}[]
    nf < nf2 || ((zf, zf2) = (zf2, zf))
    atol > 0 || (norms = max(norm(coefficients(psys.A.M),Inf),norm(coefficients(psys.B.M),Inf),norm(coefficients(psys.C.M),Inf),norm(coefficients(psys.D.M),Inf)))
    tol = atol > 0 ? atol : (rtol > 0 ? rtol*norms : sqrt(eps(T))*norms)
    for i = 1:min(nf,nf2)
        minimum(abs.(zf2 .- zf[i])) < tol  && push!(σf,zf[i])
    end
    isreal(σf) && (σf = real(σf))

    if any(isinf.(z))
       # Conjecture: The number of infinite zeros is the same as that of the time-evaluated system!
       zm = spzeros(dssdata(psteval(psys, period*rand()))...; fast, atol1 = atol, atol2 = atol, rtol)[1]
       zf = [σf; zm[isinf.(zm)]]
    end
    nz = length(zf)
    nz > n && (@warn "$(nz-n) spurious finite zero(s) present")
    return zf
end
"""
     pszero(psys::PeriodicStateSpace{PeriodicMatrix}[, K]; atol, rtol, fast) -> val
     pszero(psys::PeriodicStateSpace{PeriodicArray}[, K]; atol, rtol, fast) -> val

Compute the finite and infinite zeros of a discrete-time periodic system `psys = (A(t), B(t), C(t), D(t))` in `val`,
where the periodic system matrices `A(t)`, `B(t)`, `C(t)`, and `D(t)` are in either _Periodic Matrix_ or _Periodic Array_ representation.
The optional argument `K` specifies a desired time to start the sequence of periodic matrices (default: `K = 1`).

The computation of zeros relies on the _fast_ structure exploiting reduction of singular periodic pencils as described in [1],
which separates a matrix pencil `M-λN` which contains the infinite and finite zeros
and the left and right Kronecker indices.
The reduction is performed using orthonal similarity transformations and involves rank decisions based
on rank revealing QR-decompositions with column pivoting,
if `fast = true`, or, the more reliable, SVD-decompositions, if `fast = false`.

The keyword arguments `atol` and `rtol` specify the absolute and relative tolerances for the nonzero
elements of the periodic system matrices, respectively.
The default relative tolerance is `n*ϵ`, where `n` is proportional with the largest dimension of the periodic matrices,
and `ϵ` is the working machine epsilon.

_References_

[1] A. Varga and P. Van Dooren. Computing the zeros of periodic descriptor systems.
    Systems & Control Letters 50:371–381, 2003.
"""
function pszero(psys::PeriodicStateSpace{<: AbstractPeriodicArray{:d,T}}, K::Int = 1; fast::Bool = true, atol::Real = 0, rtol::Real = 0) where {T}
    islti(psys)  && (return spzeros(psys.A.M[1], I, psys.B.M[1], psys.C.M[1], psys.D.M[1]; fast, atol1 = atol, atol2 = atol, rtol)[1])

    (na, nb, nc, nd) = (psys.A.dperiod, psys.B.dperiod, psys.C.dperiod, psys.D.dperiod)
    N = na*psys.A.nperiod
    p, m = size(psys)
    ndx, nx = size(psys.A)
    patype = length(nx) == 1
    si = [getpm(psys.A,K,na) getpm(psys.B,K,nb); getpm(psys.C,K,nc) getpm(psys.D,K,nd)]
    ndxi = ndx[patype ? 1 : mod(K-1,na)+1]
    nxi1 = nx[patype ? 1 : mod(K,na)+1]
    ti = [ -I zeros(T,ndxi,m); zeros(T,p,nxi1+m)]
    tolr = atol
    atol == 0 && (sc = sum(nx.+p)* eps(float(T)))
    n1 = size(si,2)
    for i = K:K+N-3
        m1 = size(si,1)
        si1 = [getpm(psys.A,i+1,na) getpm(psys.B,i+1,nb); getpm(psys.C,i+1,nc) getpm(psys.D,i+1,nd)]
        mi1 = size(si1,1)
        ndxi1 = mi1-p
        nxi2 = nx[patype ? 1 : mod(i+1,na)+1]
        ti1 = [ -I zeros(T,ndxi1,m); zeros(T,p,nxi2+m)]

        F = qr([ ti; si1 ], ColumnNorm())
        nr = minimum(size(F.R))
        # compute rank of r
        ss = abs.(diag(F.R[1:nr,1:nr]))
        atol == 0 && ( tolr = sc * maximum(ss))
        rankr = count(ss .> tolr)
        si = F.Q'*[si; zeros(T,mi1,n1)]; si=si[rankr+1:end,:]
        ti = F.Q'*[ zeros(T,m1,nxi2+m); ti1]; ti = ti[rankr+1:end,:]
    end
    sn = [getpm(psys.A,K+N-1,na) getpm(psys.B,K+N-1,nb); getpm(psys.C,K+N-1,nc) getpm(psys.D,K+N-1,nd)]
    ndxi = ndx[patype ? 1 : mod(K+N-2,na)+1]
    nxi1 = nx[patype ? 1 : mod(K+N-1,na)+1]
    tn = [ I zeros(T,ndxi,m); zeros(T,p,nxi1+m)]
    a = [ zeros(T,size(tn)...) sn; si ti]
    e = [ tn zeros(T,size(sn)...); zeros(T,size(si)...) zeros(T,size(ti)...)]
    return complex(pzeros(a, e; fast, atol1 = atol, atol2 = atol, rtol)[1]).^(1/N)
end
function pszero(psys::PeriodicStateSpace{<: Union{PeriodicFunctionMatrix, PeriodicSymbolicMatrix, PeriodicTimeSeriesMatrix}}, K::Union{Int,Missing} = missing; kwargs...)
    return pszero(convert(PeriodicStateSpace{HarmonicArray},psys), K; kwargs...)
end
"""
     isstable(psys[, K = 1]; smarg = 1, fast = false, offset = sqrt(ϵ), kwargs...) -> Bool

Return `true` if the continuous-time periodic system `psys` has only stable characteristic multipliers and `false` otherwise.

To assess the stability, the absolute values of the characteristic multipliers (i.e., the eigenvalues of the monodromy matrix)
must be less than `smarg-β`, where `smarg` is the discrete-time stability margin (default: `smarg = 1`)  and
`β` is an offset specified via the keyword parameter `offset = β` to be used to numerically assess the stability
of eigenvalues. The default value used for `β` is `sqrt(ϵ)`, where `ϵ` is the working machine precision.

The monodromy matrix is determined as a product `K` state transition matrices (default: `K = 1`)
computed by integrating numerically a homogeneous linear ODE with periodic coefficients
(see function [`monodromy`](@ref) for options which can be specified via the keyword arguments `kwargs`).
If `fast = false` (default) then the characteristic multipliers are computed using an approach
based on the periodic Schur decomposition [1], while if `fast = true`
the structure exploiting reduction [2] of an appropriate lifted pencil is employed.
This later option may occasionally lead to inaccurate results for large number of matrices.

_References_

[1] A. Bojanczyk, G. Golub, and P. Van Dooren,
    The periodic Schur decomposition. Algorithms and applications, Proc. SPIE 1996.

[2] A. Varga & P. Van Dooren. Computing the zeros of periodic descriptor systems.
    Systems and Control Letters, 50:371-381, 2003.
"""
function isstable(psys::PeriodicStateSpace{<: AbstractPeriodicArray{:c,T}}, K::Int = 1; smarg::Real = 1,
                  fast::Bool = false, offset::Real = sqrt(eps(float(real(T)))), kwargs...) where T
    ev = pseig(monodromy(convert(PeriodicFunctionMatrix,psys.A), K; kwargs...); fast)
    return all(abs.(ev) .< smarg-abs(offset))
end
"""
    isstable(psys; smarg = 1, fast = false, offset = sqrt(ϵ)) -> Bool

Return `true` if the discrete-time periodic system `psys` has only stable characteristic multipliers and `false` otherwise.

To assess the stability, the absolute values of the characteristic multipliers (i.e., the eigenvalues of the monodromy matrix)
must be less than `smarg-β`, where `smarg` is the discrete-time stability margin (default: `smarg = 1`)  and
`β` is an offset specified via the keyword parameter `offset = β` to be used to numerically assess the stability
of eigenvalues. The default value used for `β` is `sqrt(ϵ)`, where `ϵ` is the working machine precision.

If `fast = false` (default) then the characteristic multipliers are computed using an approach
based on the periodic Schur decomposition [1], while if `fast = true`
the structure exploiting reduction [2] of an appropriate lifted pencil is employed.
This later option may occasionally lead to inaccurate results for large number of matrices.

_References_

[1] A. Bojanczyk, G. Golub, and P. Van Dooren,
    The periodic Schur decomposition. Algorithms and applications, Proc. SPIE 1996.

[2] A. Varga & P. Van Dooren. Computing the zeros of periodic descriptor systems.
    Systems and Control Letters, 50:371-381, 2003.
"""
function isstable(psys::PeriodicStateSpace{<: AbstractPeriodicArray{:d,T}}; smarg::Real = 1, fast::Bool = false,
                  offset::Real = sqrt(eps(float(real(T))))) where T
    ev = pseig(psys.A; fast)
    return all(abs.(ev) .< smarg-abs(offset))
end


text:
module Spline

## spline construction and evaluation using "dierckx" FORTRAN library via "Dierckx" julia package
## http://www.netlib.org/dierckx/
## https://github.com/kbarbary/Dierckx.jl
using Dierckx
using QuadGK
##using Gadfly

export spline, hazard, cumhaz, eval_basis, plot_spline, pad_knots, const_spline

## FIXME: consistent use of Q vs k

## construct spline of order k+1 on interval (outer[1], outer[2])
## with given internal knots and weights theta
## knots = vector of interior knots
## w = vector of weights (one for each basis spline)
function spline(outer::NTuple{2,Float64}, knots::Vector{Float64}, w::Vector{Float64}, k::Int=3)
    ## number of basis splines
    m = length(knots) + k + 1

    @assert (outer[1] >= 0.0) & (outer[2] >= 0.0)

    if m != length(w)
        error("number of coefficients must match number of basis splines")
        return NaN
    end

    ## if !issorted(knots)
    ##   knots = sort(knots)
    ## end

    ## total number of knots
    n = m + k + 1

    ## pad knots by repeating outer knot (k+1)-times on each side
    t = Vector{Float64}(undef, n) ##Vector{Float64}(undef,n)  # All knots
    t[(k+2):(end-k-1)] = knots
    t[1:(k+1)] .= outer[1]
    t[(end-k):end] .= outer[2]

    Spline1D(t, w, k, 3, 0.0)
end

"""
    const_spline(val::Float64, outer::NTuple{2,Float64}; Q::Int=4)

Create spline of order `Q` with outer knots at `outer` and all weights equal to `val`. Resulting spline
is a constant function with value `val` everywhere.
"""
function const_spline(val::Float64, outer::NTuple{2,Float64}; Q::Int=4)
    spline(outer, [outer[1] + (outer[2]-outer[1])/2], val .+ zeros(Float64, Q+1), Q-1)
end

get_inner_knots(s::Dierckx.Spline1D) = s.t[(s.k+2):(end-s.k-1)]
get_outer_knots(s::Dierckx.Spline1D) = (s.t[1], s.t[end])

function pad_knots(knots::Vector{Float64}, outer::NTuple{2,Float64}, p::Int)
    pk = Vector{Float64}(undef, length(knots) + 2*(p+1))
    pk[(p+2):(end-p-1)] = knots
    pk[1:(p+1)] .= outer[1]
    pk[(end-p):end] .= outer[2]
    pk
end

hazard(x, s::Dierckx.Spline1D) = evaluate(s, x)

## evaluate cumulative hazard from 0.0 to x
function cumhaz(x::Vector{T}, s::Dierckx.Spline1D) where T <: Real
    res = zeros(eltype(x), length(x))
    for i in eachindex(x)
        if (s.t[1] <= 0.0) & (x[i] <= s.t[end])
            res[i] = integrate(s, 0.0, x[i])
        else
            res[i] = quadgk(t -> hazard(t,s), 0.0, x[i])[1]
        end
    end
    res
end

cumhaz(x::T, s::Dierckx.Spline1D) where T <: Real = cumhaz([x], s)[1]

function plot_spline(s::Dierckx.Spline1D)
    a = s.t[1]
    b = s.t[end] + 10.0

    x = Vector(a:((b-a)/1000):b)
    y = evaluate(s, x)

    plot(x=x, y=y, Geom.line)
end

"""
    eval_basis(s::Dierckx.Spline1D, x::Vector{Float64}, y::Vector{Float64})

Evaluate all basis functions of 1D spline `s` at points `x`, and all integrals of spline functions at `y`.
"""
function eval_basis(s::Dierckx.Spline1D, x::Vector{T}, y::Vector{T}=x) where T <: Real
    eval_basis(x, y, get_inner_knots(s), get_outer_knots(s), s.k+1)
end

function eval_basis(x::Vector{T}, y::Vector{T}, knots::Vector{T}, outer=NTuple{2,T}, Q::Int=4) where T <: Real
    K = length(knots) + Q
    A = Matrix{T}(undef, length(x), K)
    B = Matrix{T}(undef, length(y), K)

    for k in 1:K
        w = zeros(T, K)
        w[k] = 1.0
        s = spline(outer, knots, w, Q-1)
        A[:,k] = hazard(x, s)
        B[:,k] = cumhaz(y, s)
    end

    A, B
end

end

function test_spline()
    outer = (0.0, 1.0)
    knots = [0.2, 0.4, 0.6, 0.8]

    x = rand(100)
    theta = rand(Normal(0, 1), 8)
    A = Array{Float64,2}(undef, 100, 8)

    for k in 1:8
        w = zeros(Float64, 8)
        w[k] = 1.0
        s = spline(outer, knots, w, 3)
        A[:,k] = hazard(x, s)
    end

    a = A*exp.(theta)
    b = hazard(x, spline(outer, knots, exp.(theta), 3))

    display(hcat(a, b))
    isapprox(a, b)
end

function test_splineA(M=10, n=100)
    outer = (0.0, 1.0)
    knots = [0.2, 0.4, 0.6, 0.8]

    x = rand(n)
    h = zeros(Float64, n)
    ch = zeros(Float64, n)
    A = Array{Float64,2}(undef, n, 8)
    B = Array{Float64,2}(undef, n, 8)

    for k in 1:8
        w = zeros(Float64, 8)
        w[k] = 1.0
        s = spline(outer, knots, w, 3)
        A[:,k] = hazard(x, s)
        B[:,k] = cumhaz(x, s)
    end

    for m in 1:M
        theta = rand(Normal(0, 1), 8)
        et = exp.(theta)
        h += A*et
        ch += B*et
    end

    h ./ M, ch ./ M
end

function test_splineB(M=10, n=100)
    outer = (0.0, 1.0)
    knots = [0.2, 0.4, 0.6, 0.8]

    x = rand(n)
    h = zeros(Float64, n)
    ch = zeros(Float64, n)

    for m in 1:M
        theta = rand(Normal(0, 1), 8)
        s = spline(outer, knots, exp.(theta), 3)
        h += hazard(x, s)
        ch += cumhaz(x, s)
    end

    h ./ M, ch ./ M
end


text:
export get_steady_states
export get_single_solution

"""
$(TYPEDSIGNATURES)
Return a dictionary specifying all variables and parameters of the solution
in `result` on `branch` at the position `index`.
"""
function get_single_solution(res::Result; branch::Int64, index)

    # check if the dimensionality of index matches the solutions
    if length(size(res.solutions)) !== length(index)
        # if index is a number, use linear indexing
        index = length(index) == 1 ? CartesianIndices(res.solutions)[index] : error("Index ", index, " undefined for a solution of size ", size(res.solutions))
    else
        index = CartesianIndex(index)
    end

    vars = Dict(zip(res.problem.variables, res.solutions[index][branch]))

    # collect the swept parameters required for this call
    swept_params = Dict( key => res.swept_parameters[key][index[i]] for (i,key) in enumerate(keys(res.swept_parameters)))
    full_solution = merge(vars, res.fixed_parameters, swept_params)

    return Dict(zip(keys(full_solution), ComplexF64.(values(full_solution))))
end

get_single_solution(res::Result, index) = [get_single_solution(res, index=index, branch = b) for b in 1:length(res.solutions[1])]

"""
    get_steady_states(prob::Problem,
                        swept_parameters::ParameterRange,
                        fixed_parameters::ParameterList;
                        random_warmup=false,
                            threading=false,
                            sorting="hilbert")

Solves `prob` over the ranges specified by `swept_parameters`, keeping `fixed_parameters` constant.
`swept_parameters` accepts pairs mapping symbolic variables to arrays or `LinRange`.
`fixed_parameters` accepts pairs mapping symbolic variables to numbers.

Keyword arguments
- `random_warmup`: If `true`, a problem similar to `prob` but with random complex parameters is first solved to find all non-singular paths. The subsequent tracking to find results for all swept_parameters is then much faster than the initial solving. If `random_warmup=false`, each parameter point is solved separately by tracking the maximum number of paths (employs a total degree homotopy).
This takes far longer but can be more reliable.
- `threading`: If `true`, multithreaded support is activated. The number of available threads is set by the environment variable `JULIA_NUM_THREADS`.
- `sorting`: the method used by `sort_solutions` to get continuous solutions branches.  The current options are `"hilbert"` (1D sorting along a Hilbert curve), `"nearest"` (nearest-neighbor sorting) and `"none"`.

Example: solving a simple harmonic oscillator ``m \\ddot{x} + γ \\dot{x} + ω_0^2 x = F \\cos(ωt)``
to obtain the response as a function of ``ω``
```julia-repl
# having obtained a Problem object, let's find steady states
julia> range = ParameterRange(ω => LinRange(0.8,1.2,100) ) # 100 parameter sets to solve
julia> fixed = ParameterList(m => 1, γ => 0.01, F => 0.5, ω_0 => 1)
julia> get_steady_states(problem, range, fixed)

A steady state result for 100 parameter points

    Solution branches:   1
       of which real:    1
       of which stable:  1

    Classes: stable, physical, Hopf, binary_labels

```

It is also possible to create multi-dimensional solutions plots.
```julia-repl
# The swept parameters take precedence over fixed -> use the same fixed
julia> range = ParameterRange(ω => LinRange(0.8,1.2,100), F => LinRange(0.1,1.0,10) ) # 100x10 parameter sets

# The swept parameters take precedence over fixed -> the F in fixed is now ignored
julia> get_steady_states(problem, range, fixed)

A steady state result for 1000 parameter points

    Solution branches:   1
       of which real:    1
       of which stable:  1

    Classes: stable, physical, Hopf, binary_labels
```

"""
function get_steady_states(prob::Problem, swept_parameters::ParameterRange, fixed_parameters::ParameterList; random_warmup=false, threading=false, sorting="nearest")
    # make sure the variables are in our namespace to make them accessible later
    declare_variable.(string.(cat(prob.parameters, prob.variables, dims=1)))

    # prepare an array of vectors, each representing one set of input parameters
    # an n-dimensional sweep uses an n-dimensional array
    unique_fixed = filter_duplicate_parameters(swept_parameters, fixed_parameters)
    input_array = _prepare_input_params(prob, swept_parameters, unique_fixed)
    # feed the array into HomotopyContinuation, get back an similar array of solutions
    raw = _get_raw_solution(prob, input_array, sweep=swept_parameters, random_warmup=random_warmup, threading=threading)

    # extract all the information we need from results
    #rounded_solutions = unique_points.(HomotopyContinuation.solutions.(getindex.(raw, 1)); metric = EuclideanNorm(), atol=1E-14, rtol=1E-8)
    rounded_solutions = HomotopyContinuation.solutions.(getindex.(raw, 1));
    solutions = pad_solutions(rounded_solutions)

    compiled_J = _compile_Jacobian(prob, swept_parameters, unique_fixed)

    result = Result(solutions, swept_parameters, unique_fixed, prob, Dict(), compiled_J) # a "raw" solution struct

    sort_solutions!(result, sorting=sorting) # sort into branches
    classify_solutions!(result, is_physical, "physical")
    classify_solutions!(result, is_stable, "stable")
    classify_solutions!(result, is_Hopf_unstable, "Hopf")
    order_branches!(result, ["physical", "stable", "Hopf"]) # shuffle the branches to have relevant ones first
    classify_binaries!(result) # assign binaries to solutions depending on which branches are stable

    return result

end

get_steady_states(p::Problem, swept, fixed; kwargs...) = get_steady_states(p, ParameterRange(swept), ParameterList(fixed); kwargs...)
get_steady_states(eom::HarmonicEquation, swept, fixed; kwargs...) = get_steady_states(Problem(eom), swept, fixed; kwargs...)

""" Compile the Jacobian from `prob`, inserting `fixed_parameters`.
    Returns a function that takes a dictionary of variables and `swept_parameters` to give the Jacobian."""
function _compile_Jacobian(prob::Problem, swept_parameters::ParameterRange, fixed_parameters::ParameterList)
    J_variables = cat(prob.variables, collect(keys(swept_parameters)), dims=1)

    if prob.jacobian != false
        compiled_J = compile_matrix(prob.jacobian, J_variables, fixed_parameters)
    else
        compiled_J = LinearResponse.get_implicit_Jacobian(prob.eom)
    end
    compiled_J
end

"""
Take a matrix containing symbolic variables `variables` and keys of `fixed_parameters`.
Substitute the values according to `fixed_parameters` and compile into a function that takes numerical arguments
    in the order set in `variables`.
"""
function compile_matrix(matrix, variables, fixed_parameters)
    J = substitute_all(matrix, fixed_parameters)
    matrix = [build_function(el, variables) for el in J]
    matrix = eval.(matrix)
    function m(s::Dict)
        vals = [s[var] for var in variables]
        return [ComplexF64(Base.invokelatest(el, vals)) for el in matrix]
    end
    return m
end

"Order the solution branches in `res` such that close classified positively by `class` are first."
function order_branches!(res::Result, class::String)
    indices = findall( x-> x==1, any.(classify_branch(res, class)))
    order = cat(indices, setdiff(1:length(res.solutions[1]), indices), dims=1) # permutation of indices which puts stable branches first
    reorder_solutions!(res, order)
end

"Order the solution branches in `res` such that close classified positively by `classes` are first.
The order of classes has descending precedence."
function order_branches!(s::Result, classes::Vector{String})
    for class in reverse(classes)
        order_branches!(s, class)
    end
end

"Reorder the solutions in `res` to match the index permutation `order`."
function reorder_solutions!(res::Result, order::Vector{Int64})
    res.solutions = reorder_array(res.solutions, order)
    for key in keys(res.classes)
        res.classes[key] = reorder_array(res.classes[key], order)
    end
end

"Reorder the elements of `a` to match the index permutation `order`."
function reorder_array(a::Array, order::Vector{Int64})
    #typeof(a[1]) ∉ [Vector, BitVector] && error("cannot unambiguously reorder a non-vector")
    new_array = similar(a)
    for (i,el) in enumerate(a)
        new_array[i] = el[order]
    end
    return new_array
end

"prepares an input vector to be parsed to the 2D phase diagram with parameters to sweep and kwargs"
function _prepare_input_params(prob::Problem, sweeps::ParameterRange, fixed_parameters::ParameterList)
    # sweeping takes precedence over fixed_parameters
    fixed_parameters = filter_duplicate_parameters(sweeps, fixed_parameters)

    # fixed order of parameters
    all_keys = cat(collect(keys(sweeps)), collect(keys(fixed_parameters)), dims=1)
    # the order of parameters we have now does not correspond to that in prob!
    # get the order from prob and construct a permutation to rearrange our parameters
    error = ArgumentError("some input parameters are missing or do not appear in the equation")
    permutation = try
         p = [findall(x->isequal(x, par), all_keys) for par in prob.parameters] # find the matching position of each parameter
         all((length.(p)) .== 1) || throw(error) # some parameter exists more than twice!
         p = getindex.(p, 1) # all exist once -> flatten
         isequal(all_keys[getindex.(p, 1)], prob.parameters) || throw(error) # parameters sorted wrong!
         #isempty(setdiff(all_keys, prob.parameters)) || throw(error) # extra parameters present!
         p
    catch
        throw(error)
    end

    param_ranges = collect(values(sweeps)) # Vector of the sweep LinRanges
    input_array = collect(Iterators.product(param_ranges..., values(fixed_parameters)...)) # array of all permutations (fixed_params do not change)

    # order each parameter vector to match the order in prob
    input_array = getindex.(input_array, [permutation])
    # HC wants arrays, not tuples
    return tuple_to_vector.(input_array)
end

"Remove occurrences of `sweeps` elements from `fixed_parameters`."
function filter_duplicate_parameters(sweeps, fixed_parameters)
    new_params = copy(fixed_parameters)
    for par in keys(sweeps)
        delete!(new_params, par)
    end
    return new_params
end

"Uses HomotopyContinuation to solve `problem` at specified `parameter_values`."
function _get_raw_solution(problem::Problem, parameter_values::Array{ParameterVector}; sweep=[],random_warmup=false, threading=false)
    # HomotopyContinuation accepts 1D arrays of parameter sets
    params_1D = reshape(parameter_values, :, 1)

    if random_warmup
        complex_pert = [1E-2 * issubset(p, keys(sweep))*randn(ComplexF64) for p in problem.parameters] # complex perturbation of the warmup parameters
        warmup_parameters = params_1D[Int(round(length(params_1D)/2))] .* (ones(length(params_1D[1])) + complex_pert)
        warmup_solution = HomotopyContinuation.solve(problem.system,  target_parameters=warmup_parameters, threading=threading)
        result_full = HomotopyContinuation.solve(problem.system, HomotopyContinuation.solutions(warmup_solution), start_parameters=warmup_parameters, target_parameters=parameter_values, threading=threading)
    else
        result_full = Array{Vector{Any}, 1}(undef, length(parameter_values))
        bar = Progress(length(parameter_values), 1, "Solving via total degree homotopy ...", 50)
        for i in 1:length(parameter_values) # do NOT thread this
            p = parameter_values[i]
            next!(bar)
            result_full[i] = [HomotopyContinuation.solve(problem.system, start_system=:total_degree, target_parameters=p, threading=threading, show_progress=false), p]
        end
    end

    # reshape back to the original shape
    return reshape(result_full, size(parameter_values)...)
end

"Add `padding_value` to `solutions` in case their number changes in parameter space."
function pad_solutions(solutions::Array{Vector{Vector{ComplexF64}}}; padding_value=NaN)
    Ls    = length.(solutions)
    nvars = length(solutions[1][1]) # number of variables
    max_N = maximum(Ls) # length to be fixed
    padded_solutions = deepcopy(solutions)
    for (i,s) in enumerate(solutions)
        if Ls[i]<max_N
            padded_solutions[i] = vcat(solutions[i],[padding_value*ones(nvars) for k in 1:(max_N-length(s))])
        end
    end
    return padded_solutions
end

tuple_to_vector(t::Tuple) = [i for i in t]

###
# DEPRECATED
###

#=

"Returns the parameters of res which are fixed throughout all solution points."
function get_fixed_parameters(res::Result)
    swept_parameters = keys(res.swept_parameters)
    all_parameters = res.problem.parameters
    return setdiff(all_parameters, swept_parameters)
end

function is_Hopf_unstable(solns::Result; index::Int64=0, branch::Int64=0, im_tol=im_tol)
    unstable(b, i) = is_Hopf_unstable(solns.solutions[i][b], solns.parameters[i], solns.eom, im_tol=im_tol)
    (index != 0 && branch != 0) && return unstable(branch, index)
    (index == 0 && branch != 0) && return BitVector([unstable(branch, i) for i in 1:length(solns.solutions)])
    (index == 0 && branch == 0) && return [BitVector([unstable(b,i) for b in 1:length(solns.solutions[1])]) for i in 1:length(solns.solutions)]
    (index != 0 && branch == 0) && error("you must specify an index and a branch, only a branch or neither")
end =#

#is_stable(solutions::Vector{SteadyState}, parameters::ParameterVector, eom::Problem; im_tol::Float64) = BitVector([is_stable(single, parameters, eom, im_tol=im_tol) for single in solutions])
#is_Hopf_unstable(solutions::Vector{SteadyState}, parameters::ParameterVector, eom::Problem; im_tol::Float64) = BitVector([is_Hopf_unstable(single, parameters, eom, im_tol=im_tol) for single in solutions])


text:
"""
We distinguish between models objects and fitted model objects generally,
for OLSModel specifically:

We separate OLSModel {LinearModel} from its fitted counterpart FittedOLSModel
 {LinearModelFit}  - OLSModel and FittedOLSModel are immutable structs. However,
 FittedOLSModel, whilst immutable, contains a mutable attribute
 modelfit::FitOutput. fit!(model::OLSModel) returns an object of type FitOutput
 in place but fit(model::OLSModel) returns an object of type FittedOLSModel with
 updated modelfit::FitOutput.

 This means in simulations we can save memory allocation by just calling fit!
 without having to pass around y's, X's etc. It also means struct attributes
 are immutable so we can't accidentally overwrite data whilst also leaving
 FitOutput accessible to manipulation after estimating models etc.
"""

"""
    OLSModel <: LinearModel

Object describing data and variance covariance matrix to be computed using an
Ordinary Least Squares Model. Q, R are the QR decomposition pre-calculated and
carried around by OLSModel.
"""
struct OLSModel <: LinearModel
    y::Vector
    X::Matrix
    vcov::vcov
    Q::Matrix # TODO: Create a lighweight version that doesn't carry data
    R::Matrix
    X_names::Vector
    N::Int64
    K::Int64
    # Default standard errors homoscedastic
    function OLSModel(y::Vector,
                      X::Matrix;
                      vcov::vcov = vcovIID(),
                      X_names = nothing)
        N = size(X, 1)
        K = size(X, 2)
        length(y) == N || error("y and x have differing numbers
        of observations.")
        Q, R = qr(X)

        if isnothing(X_names)
            X_names = "x_" .* string.(1:size(X, 2))
        end
        new(y, X, vcov, Q, R, X_names, N, K)
    end
end

# Adding DataFrames functionality
function OLSModel(y::Vector, X::DataFrame; vcov::vcov = vcovIID())
    X_names = names(X)
    X = Matrix(X) # Backcompability warning requires Juila > 0.7
    return OLSModel(y, X, vcov = vcov, X_names = X_names)
end

mutable struct FitOutput
    β::Vector
    se_β::Vector
    pval::Vector
    resid::Vector
    σ_sq::Float64
    vcov_matrix::Matrix
end

struct FittedOLSModel <: LinearModelFit
    y::Vector
    X::Matrix
    vcov::vcov
    Q::Matrix
    R::Matrix
    X_names::Vector
    N::Int64
    K::Int64
    modelfit::FitOutput
end

"""
    fit!(model::OLSModel)

Use the QR decomposition to estimate y = X β and return FitOutput.
"""
function fit!(model::OLSModel)
    y = model.y
    X = model.X
    # QR decomposition for numerical stability
    Q = model.Q
    R = model.R
    # Find beta
    β = inv(R) * Q' * y
    # SEs
    resid = y - X*β
    se_β, pval, σ_sq, vcov_matrix = inference(resid, β, model, model.vcov)
    return FitOutput(β, se_β, pval, resid, σ_sq, vcov_matrix)
end

"""
    fit(model::OLSModel)

Instead of returning FitOutput, returns a FittedOLSModel.
"""
function fit(model::OLSModel)
    FittedOLSModel(
        model.y,
        model.X,
        model.vcov,
        model.Q,
        model.R,
        model.X_names,
        model.N,
        model.K,
        fit!(model))
end

max_stars_repo_name: EdJeeOnGitHub/jeeves

max_issues_repo_name: EdJeeOnGitHub/jeeves

max_forks_repo_name: EdJeeOnGitHub/jeeves


text:
# %% imports

include("create_mat.jl")
include("clusterize_mat.jl")
include("hclustrecipe.jl")

# %% code to classify and visualize the clusters

# Load the data
dic = load("Julia/MatrixCluster/mc.jld")
mc = dic["mc"]

# run a param research select most appropriate one manually
# (best compromise between lowest order, higher n_clust, smaller size_bigger_clust)
param_search = find_top_param(mc)

# %% Label the matrices section

params = param_search[26]       ## "cosine" :: 19, 23, 36, 39 ## "euclidean" ::15, 26, 35
labels = label_mat(mc, params)
uni = unique_label(labels)

# %% plot the dendrogram
params = param_search[26]
plot_dend(mc, metric = params[2])

# %% plot the matrice distribution

#from uni extract a vector 1D
di = []
sizes = []
for idx in 1:length(uni)-1
    push!(di, uni[idx]...)
    push!(sizes, length(uni[idx]))
end

histogram(di, bins = 1:1:23,
            xlabel = "Clusters",
            ylabel = "Number of M_alpha using each cluster",
            legend = false,
            )

m = mean(sizes)

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
using ADCME
using AdFem

include("chip_unstructured_solver.jl")
include("chip_unstructured_geometry.jl")
include("plot_inverse_one_iter.jl")

trialnum = 1

k_mold = 0.014531
k_chip_ref = 2.60475
k_air = 0.64357

function k_exact(x, y)
    k_mold + 1000 * k_chip_ref * (x-0.49)^2 / (1 + x^2)
end

nu = 0.47893  # equal to 1/Re
power_source = 82.46295  #82.46295 = 1.0e6 divide by air rho cp   #0.0619 = 1.0e6 divide by chip die rho cp
buoyance_coef = 299102.83

u_std = 0.001
p_std = 0.000001225
T_infty = 300

NT = 14    # number of iterations for Newton's method, 8 is good for m=400

k_chip = Variable(ones(length(chip_fem_idx))) .* k_mold

heat_source_fem = zeros(ndof)
heat_source_fem[chip_fem_top_idx] .= power_source
heat_source_gauss = dof_to_gauss_points(heat_source_fem, mesh)
heat_source = compute_fem_source_term1(constant(heat_source_gauss), mesh)

B = constant(compute_interaction_matrix(mesh)) # function not exist
Laplace = constant(compute_fem_laplace_matrix1(nu * constant(ones(ngauss)), mesh))

# apply Dirichlet to velocity and temperature; set left bottom two points to zero to fix rank deficient problem for pressure
bd = [bd; bd .+ ndof;
      fvm_bd .+ 2*ndof;
     bd .+ (2*ndof+nelem)]

# add solid region into boundary condition for u, v, p, i.e. exclude solid when solving Navier Stokes
bd = [bd; solid_fem_idx; solid_fem_idx .+ ndof; solid_fvm_idx .+ 2*ndof]

S0 = constant(zeros(nelem+3*ndof))
S = solve_navier_stokes(S0, NT, k_chip)
S_computed = S[end, :]
S_data = matread("fn$trialnum/xz_chip_unstructured_data.mat")["V"]

loss =  mean((S_computed .- S_data)^2)
loss = loss * 1e10
# g = gradients(loss, θ)

# ---------------------------------------------------
_loss = Float64[]
cb = (vs, iter, loss)->begin
    global _loss
    push!(_loss, loss)
    printstyled("[#iter $iter] loss=$loss\n", color=:green)
    if mod(iter, 10)==1
        close("all")
        pix_plot_velo_pres_temp_cond(iter)
        matwrite("fn$trialnum/pix_loss$iter.mat", Dict("iter"=>iter,"loss"=>_loss, "k_chip"=>vs[1]))
    end
end

sess = Session(); init(sess)
BFGS!(sess, loss, vars = [k_chip], callback = cb)


text:
function locinfo(frame::Frame)
    if frame.framecode.scope isa Method
        meth = frame.framecode.scope
        ret = JuliaInterpreter.whereis(meth)
        ret === nothing && return nothing
        deffile, _ = ret
        ret = JuliaInterpreter.whereis(frame)
        ret === nothing && return nothing
        current_file, current_line = ret
        local body, defline
        try # https://github.com/timholy/CodeTracking.jl/issues/31
            body, defline = CodeTracking.definition(String, meth)
        catch
            return nothing
        end
        if deffile != current_file || defline > current_line
            isfile(current_file) || return nothing
            body = read(current_file, String)
            defline = 0 # We are not sure where the context start in cases like these, could be improved?
        end
        return defline, current_line, body
    else
        println("not yet implemented")
    end
end

# Used for the tests
const _print_full_path = Ref(true)

function locdesc(frame::Frame; current_line=false)
    sprint() do io
        if frame.framecode.scope isa Method
            locdesc(io, frame; current_line=current_line)
        else
            println(io, "not yet implemented")
        end
    end
end

function locdesc(io, frame::Frame; current_line=false)
    framecode = frame.framecode
    meth = framecode.scope
    @assert meth isa Method

    argnames = framecode.src.slotnames[2:meth.nargs]
    spectypes = Any[Any for i=1:length(argnames)]

    is_kw = false
    if frame.caller !== nothing
        is_kw = occursin("#kw##", string(frame.caller.framecode.scope))
    end
    if is_kw
        i = 0
        for arg in argnames
            if arg == Symbol("")
                break
            end
            i += 1
        end
        kw_indices = 1:i
        positional_indices = i+2:length(argnames)
    else
        kw_indices = 1:0
        positional_indices = 1:length(argnames)
    end

    methname = string(meth.name)
    if is_kw
        m = match(r"#(.*?)#(?:[0-9]*)$", methname)
        m === nothing || (methname = m.captures[1])
    end
    print(io, methname,'(')
    function print_indices(indices)
        first = true
        for (argname, argT) in zip(argnames[indices], spectypes[indices])
            first || print(io, ", ")
            first = false
            print(io, argname)
            argT === Any || print(io, "::", argT)
        end
    end
    print_indices(positional_indices)
    if !isempty(kw_indices)
        print(io, "; ")
        print_indices(kw_indices)
    end
    line = current_line ? JuliaInterpreter.linenumber(frame) : meth.line
    path = string(_print_full_path[] ? meth.file : basename(String(meth.file)), ":", line)
    path = CodeTracking.replace_buildbot_stdlibpath(String(path))
    print(io, ") at ", path)
end


text:
module MDP

    export Model,
           ProbModel,
           DynamicModel,
           valueIteration,
           finiteHorizon

    using Printf

    abstract type Model end

    struct ProbModel <: Model
        bellmanUpdate! :: Function
        objective  :: Function
        contractionFactor :: Float64
        stateSize  :: Int
        actionSize :: Int

        P :: Array{ Matrix{Float64}, 1}
        cost   :: Matrix{Float64}
        reward :: Matrix{Float64}

        function ProbModel(c, P; objective=:Max)
            (n, m) = size(c)

            if length(P) != m
                error("Number of transition matrices does not match number of actions.")
            end

            P_concatenated = vcat(P...)
            if size(P_concatenated) != (n*m, n)
                error("Size of transition and reward matrices are inconsistent")
            end

            is_square(Pi)          = size(Pi) == (n,n)
            is_row_stochastic(Pi)  = isapprox(sum(Pi; dims=2), ones(n); atol = 100*eps(Float64))
            is_stopping_action(Pi) = Pi == zero(Pi)
            for Pi in P
                if !is_square(Pi)
                    error("Transition matrix is not a square matrix")
                elseif !(is_row_stochastic(Pi) || is_stopping_action(Pi))
                    error("Transition matrix is not row stochastic.")
                end
            end

            if objective != :Max && objective != :Min
                error("Model objective must be :Max or :Min")
            else
                obj, compare = (objective == :Max) ? (maximum, >) : (minimum, <)

                function bellman!(vUpdated, gUpdated, v; discount=1)
                    Q = c + discount * reshape(P_concatenated * v, n, m);

                    # vUpdated = m.objective(Q, (), 2)
                    withIndex!(vUpdated, gUpdated, compare, Q)
                end

                # See Puterman Thm 6.6.6
                contractionFactor = 1 - sum(minimum(P_concatenated; dims=1))

                new(bellman!, obj, contractionFactor, n, m, P, c, c)
            end
        end
    end

    struct DynamicModel <: Model
        bellmanUpdate! :: Function # (valueFunction; discount=1.0) -> (valueFunction, policy)
        objective  :: Function
        contractionFactor :: Float64

        function DynamicModel(bellmanUpdate!; contractionFactor=1, objective=:Max)
            if objective != :Max && objective != :Min
                error("Model objective must be :Max or :Min")
            else
                obj = (objective == :Max) ? maximum : minimum
                new(bellmanUpdate!, obj, contractionFactor)
            end
        end

    end

    function valueIteration(m::Model, initial_v;
                    discount   = 0.95,
                    iterations = 1_000,
                    tolerance  = 1e-4)

        scaledDiscount  = (discount < 1) ? (1-discount)/discount : 1
        scaledTolerance = scaledDiscount * tolerance / 2

        update!(vUpdated, gUpdated, v) = m.bellmanUpdate!(vUpdated, gUpdated, v; discount=discount)

        v_previous = copy(initial_v)
        v          = copy(initial_v)
        g          = zeros(Int, size(v))

        update!(v, g, initial_v)
        copy!(v_previous, v)

        v_precision = spanNorm(v, initial_v)

        if abs( 1 - m.contractionFactor * discount ) < 4*eps(Float64)
            warn("Contraction factor too small to guarantee convergece. Value iteration may not converge.")
        else
            # See Puterman Prop 6.6.5
            # We compare with zero to allow overflow errors when v_precision is 0.
            if abs(v_precision) < 4*eps(Float64)
                iteration_bound = 1
            else
                iteration_bound = floor(Int, log( scaledTolerance/v_precision ) / log( m.contractionFactor*discount ))
            end

            @info "value iteration will converge in at most $iteration_bound iterations"
            if (iterations <= iteration_bound)
                @warn "Value iteration may not converge. Iterations $iterations less than estimated bound $iteration_bound"
            end
        end

        iterationCount  = 1;
        while (v_precision > scaledTolerance && iterationCount < iterations)
            iterationCount += 1

            copy!(v_previous, v)
            update!(v, g, v_previous)

            v_precision = spanNorm(v, v_previous)
        end

        if (v_precision > scaledTolerance)
            warn(@sprintf("Value iteration did not converge.
                 Reached precision %e at iteration %d", 2*v_precision/scaledDiscount, iterationCount))
        else
            @info @sprintf("Reached precision %e at iteration %d", 2*v_precision/scaledDiscount, iterationCount)
        end

        # Renormalize v -- See Puterman 6.6.12 for details
        v .+= m.objective(v - v_previous)/scaledDiscount

        return (v, g)
    end

    # For the probability model, it is possible to automatically initialize
    # the initial vector
    function valueIteration(m::ProbModel;
                    initial_v  = zeros(m.stateSize),
                    discount   = 0.95,
                    iterations = 1_000,
                    tolerance  = 1e-4)

        return valueIteration(m, initial_v;
                    discount   = discount,
                    iterations = iterations,
                    tolerance  = tolerance)
    end

    function finiteHorizon(m::Model, final_v;
                           discount = 1.0,
                           horizon :: Int = 10)

        update!(vUpdated, gUpdated, v) = m.bellmanUpdate!(vUpdated, gUpdated, v; discount=discount)

        v = [ zero(final_v) for stage = 1 : horizon ]
        g = [ zeros(Int,   size(final_v)) for stage = 1 : horizon ]

        update!(v[horizon], g[horizon], final_v)
        for stage in horizon-1: -1 : 1
          update!(v[stage], g[stage], v[stage+1])
        end
        return (v,g)
    end

    function spanNorm(x, y)
        # z = x - y
        # return max(z) - min(z)
        # Optimized code. This could have been done using @devec, but
        # I don't want to add a dependency just for one function.
        z = vec(x - y)
        max_z, min_z = -Inf, Inf
        for elem in z
            if max_z < elem
                max_z = elem
            end
            if min_z > elem
                min_z = elem
            end
        end
        return max_z - min_z
    end

    # Julia does not have an in-built function that returns the minimum and the
    # arg min.
    # function withIndex{T}(f::Function, x::Matrix{T})
    #     idx = vec(mapslices(f,x,2))
    #     val :: Vector{T} = [x[i,idx[i]] for i = 1:size(x,1) ]

    #     return val,idx
    # end

    # A more direct implementation
    function withIndex!(val::AbstractArray{T,1}, idx::AbstractArray{Int,1}, compare::Function, x::AbstractArray{T,2}) where T
        (n, m) = size(x)

        for i=1:n
            idx[i], val[i] = 1, x[i,1]
            for j=2:m
                if compare(x[i,j], val[i])
                    idx[i], val[i] = j, x[i,j]
                end
            end
        end
    end
end