text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct AdminTeamsSettingsApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_adminTeamsSettingsInfo(_api::AdminTeamsSettingsApi, team_id::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/admin.teams.settings.info", ["slackAuth"])
    Swagger.set_param(_ctx.query, "team_id", team_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Fetch information about settings in a workspace
Param: team_id::String (required)
Return: Dict{String, Any}
"""
function adminTeamsSettingsInfo(_api::AdminTeamsSettingsApi, team_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminTeamsSettingsInfo(_api, team_id; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminTeamsSettingsInfo(_api::AdminTeamsSettingsApi, response_stream::Channel, team_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminTeamsSettingsInfo(_api, team_id; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminTeamsSettingsSetDefaultChannels(_api::AdminTeamsSettingsApi, team_id::String, channel_ids::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.teams.settings.setDefaultChannels", ["slackAuth"])
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.form, "channel_ids", channel_ids)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Set the default channels of a workspace.
Param: team_id::String (required)
Param: channel_ids::String (required)
Return: Dict{String, Any}
"""
function adminTeamsSettingsSetDefaultChannels(_api::AdminTeamsSettingsApi, team_id::String, channel_ids::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminTeamsSettingsSetDefaultChannels(_api, team_id, channel_ids; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminTeamsSettingsSetDefaultChannels(_api::AdminTeamsSettingsApi, response_stream::Channel, team_id::String, channel_ids::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminTeamsSettingsSetDefaultChannels(_api, team_id, channel_ids; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminTeamsSettingsSetDescription(_api::AdminTeamsSettingsApi, team_id::String, description::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.teams.settings.setDescription", ["slackAuth"])
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.form, "description", description)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Set the description of a given workspace.
Param: team_id::String (required)
Param: description::String (required)
Return: Dict{String, Any}
"""
function adminTeamsSettingsSetDescription(_api::AdminTeamsSettingsApi, team_id::String, description::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminTeamsSettingsSetDescription(_api, team_id, description; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminTeamsSettingsSetDescription(_api::AdminTeamsSettingsApi, response_stream::Channel, team_id::String, description::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminTeamsSettingsSetDescription(_api, team_id, description; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminTeamsSettingsSetDiscoverability(_api::AdminTeamsSettingsApi, team_id::String, discoverability::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.teams.settings.setDiscoverability", ["slackAuth"])
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.form, "discoverability", discoverability)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

An API method that allows admins to set the discoverability of a given workspace
Param: team_id::String (required)
Param: discoverability::String (required)
Return: Dict{String, Any}
"""
function adminTeamsSettingsSetDiscoverability(_api::AdminTeamsSettingsApi, team_id::String, discoverability::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminTeamsSettingsSetDiscoverability(_api, team_id, discoverability; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminTeamsSettingsSetDiscoverability(_api::AdminTeamsSettingsApi, response_stream::Channel, team_id::String, discoverability::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminTeamsSettingsSetDiscoverability(_api, team_id, discoverability; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminTeamsSettingsSetIcon(_api::AdminTeamsSettingsApi, image_url::String, team_id::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.teams.settings.setIcon", ["slackAuth"])
    Swagger.set_param(_ctx.form, "image_url", image_url)  # type String
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Sets the icon of a workspace.
Param: image_url::String (required)
Param: team_id::String (required)
Return: Dict{String, Any}
"""
function adminTeamsSettingsSetIcon(_api::AdminTeamsSettingsApi, image_url::String, team_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminTeamsSettingsSetIcon(_api, image_url, team_id; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminTeamsSettingsSetIcon(_api::AdminTeamsSettingsApi, response_stream::Channel, image_url::String, team_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminTeamsSettingsSetIcon(_api, image_url, team_id; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminTeamsSettingsSetName(_api::AdminTeamsSettingsApi, team_id::String, name::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.teams.settings.setName", ["slackAuth"])
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.form, "name", name)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Set the name of a given workspace.
Param: team_id::String (required)
Param: name::String (required)
Return: Dict{String, Any}
"""
function adminTeamsSettingsSetName(_api::AdminTeamsSettingsApi, team_id::String, name::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminTeamsSettingsSetName(_api, team_id, name; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminTeamsSettingsSetName(_api::AdminTeamsSettingsApi, response_stream::Channel, team_id::String, name::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminTeamsSettingsSetName(_api, team_id, name; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export adminTeamsSettingsInfo, adminTeamsSettingsSetDefaultChannels, adminTeamsSettingsSetDescription, adminTeamsSettingsSetDiscoverability, adminTeamsSettingsSetIcon, adminTeamsSettingsSetName


text:
#=
Basic types for strings

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
=#
const STR_KEEP_NUL    = true  # keep nul byte placed by String

# Note: this is still in transition to expressing character set, encoding
# and optional cached info for hashes, UTF-8/UTF-16 encodings, subsets, etc.
# via more type parameters

struct Str{T,SubStr,Cache,Hash} <: AbstractString
    data::String
    substr::SubStr
    cache::Cache
    hash::Hash

    (::Type{Str})(::Type{T}, v::String, s::S, c::C, h::H) where {T<:CSE,S,C,H} =
        new{T,S,C,H}(v,s,c,h)
end

_msk16(v, m) = (v%UInt16) & m
_msk32(v, m) = (v%UInt32) & m

_mskup16(v, m, s) = _msk16(v, m) << s
_mskup32(v, m, s) = _msk32(v, m) << s
_mskdn16(v, m, s) = _msk16(v, m) >>> s
_mskdn32(v, m, s) = _msk32(v, m) >>> s

(::Type{Str})(::Type{C}, v::String) where {C<:CSE} = Str(C, v, nothing, nothing, nothing)
(::Type{Str})(::Type{C}, v::Str) where {C<:CSE} = Str(C, v.data, nothing, nothing, nothing)

# Handle change from endof -> lastindex
@static if !isdefined(Base, :lastindex)
    lastindex(str::AbstractString) = Base.endof(str)
    lastindex(arr::AbstractArray) = Base.endof(arr)
    Base.endof(str::Str) = lastindex(str)
end
@static if !isdefined(Base, :firstindex)
    firstindex(str::AbstractString) = 1
    # AbstractVector might be an OffsetArray
    firstindex(str::Vector) = 1
end

# Definition of built-in Str types

const empty_string = ""

for lst in cse_info
    nam, typ = lst
    str = String(nam)
    sym = symstr(nam, "Str")
    cse = symstr(nam, "CSE")
    @eval const $sym = Str{$cse, Nothing, Nothing, Nothing}
    @eval (::Type{$sym})(v::Vector{UInt8}) = convert($sym, v)
    @eval show(io::IO, ::Type{$sym}) = print(io, $(quotesym(sym)))
    low = lowercase(str)
    if str[1] == '_'
        @eval empty_str(::Type{$cse}) = $(symstr("empty", low))
        @eval @api develop $sym
    else
        emp = symstr("empty_", low)
        @eval const $emp = Str($cse, empty_string)
        @eval empty_str(::Type{$cse}) = $emp
        @eval @api public $sym
    end
    @eval convert(::Type{$sym}, str::$sym) = str
end
empty_str(::Type{<:Str{C}}) where {C<:CSE} = empty_str(C)
empty_str(::Type{String}) = empty_string

typemin(::Type{T}) where {T<:Str} = empty_str(T)
typemin(::T) where {T<:Str} = empty_str(T)

"""Union type for fast dispatching"""
const UniStr = Union{ASCIIStr, _LatinStr, _UCS2Str, _UTF32Str}
show(io::IO, ::Type{UniStr}) = print(io, :UniStr)

# Display BinaryCSE as if String
show(io::IO, str::T) where {T<:Str{BinaryCSE}} = show(io, str.data)
show(io::IO, str::SubString{T}) where {T<:Str{BinaryCSE}} =
    @inbounds show(io, SubString(str.string.data, str.offset+1, str.offset+lastindex(str)))

_allocate(len) = Base._string_n((len+STR_KEEP_NUL-1)%Csize_t)

function _allocate(::Type{T}, len) where {T <: CodeUnitTypes}
    buf = _allocate((len+STR_KEEP_NUL-1) * sizeof(T))
    buf, reinterpret(Ptr{T}, pointer(buf))
end

# Various useful groups of character set types

const MS_UTF8     = MaybeSub{<:Str{UTF8CSE}}
const MS_UTF16    = MaybeSub{<:Str{UTF16CSE}}
const MS_UTF32    = MaybeSub{<:Str{UTF32CSE}}
const MS_SubUTF32 = MaybeSub{<:Str{_UTF32CSE}}
const MS_Latin    = MaybeSub{<:Str{<:Latin_CSEs}}
const MS_ByteStr  = MaybeSub{<:Str{<:Binary_CSEs}}
const MS_RawUTF8  = MaybeSub{<:Union{String,Str{RawUTF8CSE}}}

_wrap_substr(::Type{<:Any}, str) = str
_wrap_substr(::Type{SubString}, str) = SubString(str, 1)
_empty_sub(::Type{T}, ::Type{C}) where {T,C} = _wrap_substr(T, empty_str(C))

## Get the character set / encoding used by a string type
cse(::Type{<:Str{C}}) where {C<:CSE} = C

# Promotion rules for characters

promote_rule(::Type{T}, ::Type{T}) where {T<:Str} = T

promote_rule(::Type{String}, ::Type{<:Str}) = String

promote_rule(::Type{<:Str{S}}, ::Type{<:Str{T}}) where {S,T} =
    (P = promote_rule(S,T)) === Union{} ? Union{} : Str{P}

sizeof(s::Str) = sizeof(s.data) + 1 - STR_KEEP_NUL

"""Codeunits of string as a Vector"""
_data(s::Vector{UInt8}) = s
_data(s::String)        = s
_data(s::Str{<:Byte_CSEs}) = unsafe_wrap(Vector{UInt8}, s.data)

"""Pointer to codeunits of string"""
pointer(s::Str{<:Byte_CSEs}) = pointer(s.data)
pointer(s::Str{<:Word_CSEs}) = reinterpret(Ptr{UInt16}, pointer(s.data))
pointer(s::Str{<:Quad_CSEs}) = reinterpret(Ptr{UInt32}, pointer(s.data))

const CHUNKSZ = sizeof(UInt64) # used for fast processing of strings
const CHUNKMSK = (CHUNKSZ-1)%UInt64

_pnt64(s::Union{String,Vector{UInt8}}) = reinterpret(Ptr{UInt64}, pointer(s))
_pnt64(s::Str) = reinterpret(Ptr{UInt64}, pointer(s.data))

"""Length of string in codeunits"""
ncodeunits(s::Str)              = sizeof(s)
ncodeunits(s::Str{<:Word_CSEs}) = sizeof(s) >>> 1
ncodeunits(s::Str{<:Quad_CSEs}) = sizeof(s) >>> 2

# For convenience
@inline _calcpnt(str, siz) = (pnt = _pnt64(str) - CHUNKSZ;  (pnt, pnt + siz))

@inline _mask_bytes(n) = ((1%UInt) << ((n & CHUNKMSK) << 3)) - 0x1

# Support for SubString of Str

Base.SubString(str::Str{C}) where {C<:SubSet_CSEs} =
    SubString(Str(basecse(C), str))
Base.SubString(str::Str{C}, off::Int) where {C<:SubSet_CSEs} =
    SubString(Str(basecse(C), str), off)
Base.SubString(str::Str{C}, off::Int, fin::Int) where {C<:SubSet_CSEs} =
    SubString(Str(basecse(C), str), off, fin)

# pointer conversions of ASCII/UTF8/UTF16/UTF32 strings:
pointer(str::Str, pos::Integer) = bytoff(pointer(str), pos - 1)

# pointer conversions of SubString of ASCII/UTF8/UTF16/UTF32:
pointer(x::SubString{<:Str}) = bytoff(pointer(x.string), x.offset)
pointer(x::SubString{<:Str}, pos::Integer) = bytoff(pointer(x.string), x.offset + pos - 1)

ext: jl


text:
import Pkg; Pkg.activate(".")

const USE_GPU = false
using Test
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 3);
else
    @init_parallel_stencil(Threads, Float64, 3);
end

_viscosity(x, y, z) = exp(1 - ( (1-x) + y*(1-y) + z*(1-z)))

@parallel function viscosity!(η, x, y, z)
    @all(η) = _viscosity(@idx_x(x), @idx_y(y), @idx_z(z))
    return
end

function testme()

    nx = ny = nz = 8
    x = y = z = LinRange(0.0, 1.0, nx)
    η1 = @zeros(nx, ny, nz)
    @parallel viscosity!(η1, x, y, z)
    η2 = [ _viscosity(xi, yi, zi) for xi in x, yi in y, zi in z]

    return η1 == η2
end

@test  testme() == true

nx, ny, nz = 7, 5, 6
A       =  @rand(nx  , ny  , nz  );
R       = @zeros(nx  , ny  , nz  );
x = LinRange(zero(Data.Number), one(Data.Number), nx)
y = LinRange(zero(Data.Number), one(Data.Number), ny)
z = LinRange(zero(Data.Number), one(Data.Number), nz)
xvec = Data.Array(x)
yvec = Data.Array(y)
zvec = Data.Array(z)

@parallel idx_x!(R, A) = (@all(R) = @idx_x(A); return)
@parallel idx_y!(R, A) = (@all(R) = @idx_y(A); return)
@parallel idx_z!(R, A) = (@all(R) = @idx_z(A); return)

R.=0; @parallel idx_x!(R, x); @test reduce(+, Rcol == x for Rcol in eachcol(collect(eachslice(R, dims=3))[1])) == size(R,2)
R.=0; @parallel idx_y!(R, y); @test reduce(+, Rrow == y for Rrow in eachrow(collect(eachslice(R, dims=3))[1])) == size(R,1)
R.=0; @parallel idx_z!(R, z); @test reduce(+, Rrow == z for Rrow in eachrow(collect(eachslice(R, dims=1))[1])) == size(R,2)

R.=0; @parallel idx_x!(R, xvec); @test reduce(+, Rcol == x for Rcol in eachcol(collect(eachslice(R, dims=3))[1])) == size(R,2)
R.=0; @parallel idx_y!(R, yvec); @test reduce(+, Rrow == y for Rrow in eachrow(collect(eachslice(R, dims=3))[1])) == size(R,1)
R.=0; @parallel idx_z!(R, zvec); @test reduce(+, Rrow == z for Rrow in eachrow(collect(eachslice(R, dims=1))[1])) == size(R,2)


text:
using DataFrames
using Plots

using Futbolista: DataLoaders, Plotter
using Test

@testset "DataLoaders Tests" begin

    sb = DataLoaders.Statsbomb
    ldr = sb.SBLoader()

	comp_df = sb.load_competitions(ldr)
	matches_df = sb.load_matches(ldr, 72, 30)
	events_df = sb.load_events(ldr, 7298)

    @test typeof(comp_df) == DataFrame
    @test typeof(matches_df) == DataFrame
    @test typeof(events_df) == DataFrame
end

@testset "Plotting Tests" begin

    sb = DataLoaders.Statsbomb
    pl = Plotter.StatsbombPitch
    pm = Plotter.Passmap

    ldr = sb.SBLoader()
	events_df = sb.load_events(ldr, 7298)

    pitch = pl.plot_pitch(line_color=:silver, line_width=1)
    passmap = pm.plot_passmap(events_df, "home", "mean")

    @test typeof(pitch) <: Plots.Plot
    @test typeof(passmap) <: Plots.Plot
end


text:
using Pkg, Conda

pypkgs = ["bjoern", "falcon"]

@info "installing Python packages ...: $(join(pypkgs, ' '))"
Conda.add(pypkgs; channel = "conda-forge")


text:
module Signals

import Sparklines: spark

# Functions from packages that Signals extends and re-exports
import Grid: interp # TODO: move to Interpolations.jl
export interp

import SIUnits: Second, ShortUnits.s, ShortUnits.ms, ShortUnits.µs
const μs = µs # \mu != \micro. https://github.com/Keno/SIUnits.jl/issues/23
export Second, s, ms, µs, μs

export Signal,
       RegularSignal,
       SignalVector,
       SignalMatrix,
       regularize,
       isregular,
       samplingfreq,
       samplingrate,
       before,
       after,
       within,
       window

include("signal.jl")
include("show.jl")
include("interpolation.jl")
include("time.jl")

end # module


text:
module GenericSchur
using LinearAlgebra
using LinearAlgebra: Givens, Rotation
using Printf
import LinearAlgebra: lmul!, mul!, checksquare, ldiv!

# This is the public interface of the package.
# Wrappers like `schur` and `eigvals` should just work.
import LinearAlgebra: schur!, eigvals!, eigvecs, eigen!
export triangularize, eigvalscond, subspacesep, balance!

schur!(A::StridedMatrix{T}; kwargs...) where {T} = gschur!(A; kwargs...)

function eigvals!(A::StridedMatrix{T}; kwargs...) where {T}
    S = gschur!(A; wantZ=false, kwargs...)
    S.values
end

# This is probably the best we can do unless LinearAlgebra coöperates
"""
    eigvecs(S::Schur{<:Complex}; left=false) => Matrix

Compute right or left eigenvectors from a Schur decomposition.
Eigenvectors are returned as columns of a matrix.
The returned eigenvectors have unit Euclidean norm, and the largest
elements are real.
"""
function eigvecs(S::Schur{T}; left::Bool=false) where {T <: Complex}
    if left
        v = _gleigvecs!(S.T,S.Z)
    else
        v = _geigvecs!(S.T,S.Z)
    end
    _enormalize!(v)
    v
end
function _enormalize!(v::AbstractMatrix{T}) where T
    n = size(v,1)
    for j=1:n
        s = one(real(T)) / norm(v[:,j],2)
        t = abs2(v[1,j])
        i0 = 1
        for i=2:n
            u = abs2(v[i,j])
            if  u > t
                i0 = i
                t = u
            end
        end
        t = s * conj(v[i0,j]) / sqrt(t)
        for i=1:n
            v[i,j] *=  t
        end
        v[i0,j] = real(v[i0,j])
    end
end

if VERSION < v"1.2-"
    eigsortby = nothing
else
    using LinearAlgebra: eigsortby, sorteig!
end

function eigen!(A::StridedMatrix{T}; permute=true, scale=true,
                sortby::Union{Function,Nothing}=eigsortby) where {T <: Real}
    if permute || scale
        A, B = balance!(A, scale=scale, permute=permute)
    end
    S = triangularize(schur(A))
    v = _geigvecs!(S.T,S.Z)
    if permute || scale
        lmul!(B, v)
    end
    _enormalize!(v)
    if sortby !== nothing
        return LinearAlgebra.Eigen(sorteig!(S.values, v, sortby)...)
    else
        return LinearAlgebra.Eigen(S.values,v)
    end
end

function eigen!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true,
                sortby::Union{Function,Nothing}=eigsortby) where {T <: Complex}
    if permute || scale
        A, B = balance!(A, scale=scale, permute=permute)
    end
    S = schur(A)
    v = _geigvecs!(S.T,S.Z)
    if permute || scale
        lmul!(B, v)
    end
    _enormalize!(v)
    if sortby !== nothing
        return LinearAlgebra.Eigen(sorteig!(S.values, v, sortby)...)
    else
        return LinearAlgebra.Eigen(S.values,v)
    end
end

############################################################################
# Internal implementation follows

macro mydebug(expr); nothing; end

function _fmt_nr(z::Complex)
    s = @sprintf("%10.3e%+10.3eim", reim(z)...)
    return s
end
function _fmt_nr(x::Real)
    s = @sprintf("%10.3e", x)
    return s
end

include("util.jl")
include("lapack_extras.jl")
include("hessenberg.jl")

if VERSION < v"1.3"
    _getdata(H::HessenbergFactorization) = H.data
    const HessenbergArg = HessenbergFactorization
else
    _getdata(H::Hessenberg) = H.H.data
    const HessenbergArg = Hessenberg
end

include("householder.jl")
include("balance.jl")

#
# portions translated from LAPACK::zlahqr
# LAPACK Copyright:
# Univ. of Tennessee
# Univ. of California Berkeley
# Univ. of Colorado Denver
# NAG Ltd.
function _gschur!(H::HessenbergArg{T}, Z=nothing;
                 maxiter = 100*size(H, 1), maxinner = 30*size(H, 1), kwargs...
                 ) where {T <: Complex}
    n = size(H, 1)
    istart = 1
    iend = n
    w = Vector{T}(undef, n)
    HH = _getdata(H)

    RT = real(T)
    ulp = eps(RT)
    smallnum = safemin(RT) * (n / ulp)
    rzero = zero(RT)
    half = 1 / RT(2)
    threeq = 3 / RT(4)

    for j=1:n-1
        HH[j+2:n,j] .= zero(T)
    end

    # iteration count
    it = 0

    @inbounds while iend >= 1
        # key for comparison to zlahqr:
        # I => iend
        # L => istart
        # K => _istart + 1
        istart = 1
        for its=0:maxinner
            it += 1
            if it > maxiter
                throw(ArgumentError("iteration limit $maxiter reached"))
            end

            # Determine if the matrix splits.
            # Find lowest positioned subdiagonal "zero" if any; reset istart if found.
            for _istart in iend-1:-1:istart
                if abs1(HH[_istart+1, _istart]) <= smallnum
                    istart = _istart + 1
                    @mydebug println("Split1! at istart = $istart, subdiag is ",
                                     _fmt_nr(HH[istart, istart-1]))
                    break
                end
                # deflation criterion from Ahues & Tisseur (LAWN 122, 1997)
                tst = abs1(HH[_istart,_istart]) + abs1(HH[_istart+1,_istart+1])
                if tst == 0
                    if (_istart-1 >= 1) tst += abs(real(HH[_istart,_istart-1])) end
                    if (_istart+2 <= n) tst += abs(real(HH[_istart+2,_istart+1])) end
                end
                if abs(real(HH[_istart+1,_istart])) <= ulp*tst
                    ab = max(abs1(HH[_istart+1,_istart]),
                             abs1(HH[_istart,_istart+1]))
                    ba = min(abs1(HH[_istart+1,_istart]),
                             abs1(HH[_istart,_istart+1]))
                    aa = max(abs1(HH[_istart+1,_istart+1]),
                             abs1(HH[_istart,_istart]-HH[_istart+1,_istart+1]))
                    bb = min(abs1(HH[_istart+1,_istart+1]),
                             abs1(HH[_istart,_istart]-HH[_istart+1,_istart+1]))
                    s = aa + ab
                    if ba * (ab / s) <= max(smallnum, ulp * (bb * (aa / s)))
                        istart = _istart + 1
                        @mydebug println("Split2! at istart = $istart, subdiag is ",
                                         _fmt_nr(HH[istart, istart-1]))
                        break
                    end
                end
            end # check for split

            if istart > 1
                # clean up
                HH[istart, istart-1] = zero(T)
            end

            # if block size is one we deflate
            if istart >= iend
                iend -= 1
                @mydebug println("Bottom deflation! Block size is one. New iend is $iend")
                break # from "inner" loop
            end

            # select shift
            # logic adapted from LAPACK zlahqr
            if its % 30 == 10
                s = threeq * abs(real(HH[istart+1,istart]))
                t = s + HH[istart,istart]
            elseif its % 30 == 20
                s = threeq * abs(real(HH[iend,iend-1]))
                t = s + HH[iend,iend]
            else
                # Wilkinson's shift
                t = HH[iend,iend]
                u = sqrt(HH[iend-1,iend]) * sqrt(HH[iend,iend-1])
                s = abs1(u)
                if s ≠ rzero
                    x = half * (HH[iend-1,iend-1] - t)
                    sx = abs1(x)
                    s = max(s, abs1(x))
                    y = s * sqrt( (x/s)^2 + (u/s)^2)
                    if sx > rzero
                        if real(x / sx) * real(y) + imag(x / sx) * imag(y) < rzero
                            y = -y
                        end
                    end
                    t -= u * (u / (x+y))
                end
            end # shift selection

            # run a QR iteration
            @mydebug println("block range $(istart):$(iend) shift ", _fmt_nr(t))

            # following zlahqr, only use single-shift
            singleShiftQR!(HH, Z, t, istart, iend)

        end # inner loop
    end # outer loop
    w = diag(HH)
    return Schur(triu(HH), Z === nothing ? Matrix{T}(undef,0,0) : Z, w)
end

"""
gschur(A::StridedMatrix) -> F::Schur

Computes the Schur factorization of matrix `A` using a generic implementation.
See `LinearAlgebra.schur` for usage.
"""
gschur(A::StridedMatrix{T}; kwargs...) where {T} = gschur!(Matrix(A); kwargs...)

"""
gschur!(A::StridedMatrix) -> F::Schur

Destructive version of `gschur` (q.v.).
"""
function gschur!(A::StridedMatrix{T}; wantZ::Bool=true, scale::Bool=true,
                 kwargs...) where T <: Complex
    n = checksquare(A)
    if scale
        scaleA, cscale, anrm = _scale!(A)
    else
        scaleA = false
    end
    H = _hessenberg!(A)
    if wantZ
        Z = _materializeQ(H)
        S = _gschur!(H, Z; kwargs...)
    else
        S = _gschur!(H; kwargs...)
    end
    if scaleA
        safescale!(S.T, cscale, anrm)
        S.values .= diag(S.T, 0)
    end
    S
end

function singleShiftQR!(HH::StridedMatrix{T}, Z, shift::Number, istart::Integer, iend::Integer) where {T <: Complex}
    n = size(HH, 1)
    ulp = eps(real(eltype(HH)))

    @mydebug Hsave = Z*HH*Z'
    function dcheck(str)
        print(str," decomp err ", norm(Hsave - Z*HH*Z'))
    end

    # key:
    # istart => L
    # iend => I
    # istart1 => M

    # look for two consecutive small subdiagonals
    istart1 = -1
    h11s = zero(eltype(HH))
    h21 = zero(real(eltype(HH)))
    flag = false
    v = zeros(T,2)
    @inbounds for mm = iend-1:-1:istart+1
        # determine the effect of starting the single-shift Francis
        # iteration at row mm: see if this would make HH[mm,mm-1] tiny.
        h11 = HH[mm,mm]
        h22 = HH[mm+1,mm+1]
        h11s = h11 - shift
        h21 = real(HH[mm+1,mm]) # for reflector
        s = abs1(h11s) + abs(h21)
        h11s /= s
        h21 /= s
        v .= (h11s, T(h21))
        h10 = real(HH[mm,mm-1])
        if abs(h10) * abs(h21) <= ulp * (abs1(h11s) * (abs1(h11) + abs1(h22)))
            istart1 = mm
            flag = true
            break
        end
    end
    @inbounds begin
      if !flag
        istart1 = istart
        h11 = HH[istart,istart]
        h22 = HH[istart+1,istart+1]
        h11s = h11 - shift
        h21 = real(HH[istart+1,istart])
        s = abs1(h11s) + abs(h21)
        h11s /= s
        h21 /= s
        v .= (h11s,T(h21))
      end
    end

    @inbounds for k = istart1:iend-1
        if k > istart1
            # prepare to chase bulge
            v .= (HH[k,k-1], HH[k+1,k-1])
        # else
        #    use v from above to create a bulge
        end
        τ1 = LinearAlgebra.reflector!(v)
        if k > istart1
            HH[k,k-1] = v[1]
            HH[k+1,k-1] = zero(T)
        end
        v2 = v[2]
        τ2 = real(τ1*v2)
        # lmul!(R,view(HH,:,k:n))
        for j = k:n
            ss = τ1' * HH[k,j] + τ2 * HH[k+1,j]
            HH[k,j] -= ss
            HH[k+1,j] -= ss * v2
        end
        # rmul!(view(HH,1:min(iend,k+2),:),R)
        for j = 1:min(k+2,iend)
            ss = τ1 * HH[j,k] + τ2 * HH[j,k+1]
            HH[j,k] -= ss
            HH[j,k+1] -= ss * v2'
        end
        if !(Z === nothing)
            # rmul!(Z,R)
            for j = 1:n
                ss = τ1 * Z[j,k] + τ2 * Z[j,k+1]
                Z[j,k] -= ss
                Z[j,k+1] -= ss * v2'
            end
        end
        if (k == istart1) && (istart1 > istart)
            # If the QR step started at a row below istart because two consecutive
            # small subdiagonals were found, extra scaling must be performed
            # to ensure reality of HH[istart1,istart1-1]
            t = 1 - τ1
            t /= abs(t)
            HH[istart1+1,istart1] *= t'
            if istart1 + 2 <= iend
                HH[istart1+2, istart1+1] *= t
            end
            for j = istart1:iend
                if j != istart1 + 1
                    if n > j
                        HH[j,j+1:n] *= t
                    end
                    HH[1:j-1,j] *= t'
                    if !(Z === nothing)
                        Z[1:n, j] *= t'
                    end
                end
            end
        end
        @mydebug dcheck(" QR k=$k")
    end
    # ensure reality of tail
    @inbounds begin
      t = HH[iend,iend-1]
      if imag(t) != 0
        rt = abs(t)
        HH[iend, iend-1] = rt
        t /= rt
        if n > iend
            HH[iend, iend+1:n] *= t'
        end
        HH[1:iend-1, iend] *= t
        if !(Z === nothing)
            Z[1:n, iend] *= t
        end
      end
    end
    @mydebug dcheck(" QR post")
    @mydebug println()
    return HH
end

const _STANDARDIZE_DEFAULT = true

# Mostly copied from GenericLinearAlgebra
function _gschur!(H::HessenbergArg{T}, Z=nothing;
                  tol = eps(real(T)), shiftmethod = :Francis,
                  maxiter = 100*size(H, 1), standardize = _STANDARDIZE_DEFAULT,
                  kwargs...) where {T <: Real}
    n = size(H, 1)
    istart = 1
    iend = n
    HH = _getdata(H)
    τ = Rotation(Givens{T}[])
    w = Vector{Complex{T}}(undef, n)
    iwcur = n
    function putw!(w,z)
        w[iwcur] = z
        iwcur -= 1
    end

    # iteration count
    i = 0

    @inbounds while true
        i += 1
        if i > maxiter
            throw(ArgumentError("iteration limit $maxiter reached"))
        end

        # Determine if the matrix splits. Find lowest positioned subdiagonal "zero"
        for _istart in iend - 1:-1:1
            if abs(HH[_istart + 1, _istart]) < tol*(abs(HH[_istart, _istart]) + abs(HH[_istart + 1, _istart + 1]))
                    istart = _istart + 1
                @mydebug println("Split1! at $istart, subdiagonal element is ",
                                 _fmt_nr(HH[istart, istart-1]))
                standardize && (HH[istart, istart-1] = 0) # clean
                break
            elseif _istart > 1 && abs(HH[_istart, _istart - 1]) < tol*(abs(HH[_istart - 1, _istart - 1]) + abs(HH[_istart, _istart]))
                istart = _istart
                @mydebug println("Split2! at $istart, subdiagonal element is ",
                                 _fmt_nr(HH[_istart, _istart-1]))
                break
            end
            istart = 1
        end

        # if block size is one we deflate
        if istart >= iend

            standardize && putw!(w,HH[iend,iend])
            iend -= 1
            @mydebug println("Bottom deflation! Block size is one. New iend is $iend")
        # and the same for a 2x2 block
        elseif istart + 1 == iend
            @mydebug println("Bottom deflation! Block size is two. New iend is $(iend-2)")

            if standardize
                H2 = HH[iend-1:iend,iend-1:iend]
                G2,w1,w2 = _gs2x2!(H2,iend)
                putw!(w,w2)
                putw!(w,w1)
                lmul!(G2,view(HH,:,istart:n))
                rmul!(view(HH,1:iend,:),G2')
                HH[iend-1:iend,iend-1:iend] .= H2 # clean
                if iend > 2
                    HH[iend-1,iend-2] = 0
#                    HH[iend-2,iend-1] = 0
                end
                Z === nothing || rmul!(Z,G2')
            end
            iend -= 2

        # run a QR iteration
        # shift method is specified with shiftmethod kw argument
        else
            Hmm = HH[iend, iend]
            Hm1m1 = HH[iend - 1, iend - 1]
            if iszero(i % 10)
                # Use Eispack exceptional shift
                β = abs(HH[iend, iend - 1]) + abs(HH[iend - 1, iend - 2])
                d = (Hmm + β)^2 - Hmm*β/2
                t = 2*Hmm + 3*β/2
            else
                d = Hm1m1*Hmm - HH[iend, iend - 1]*HH[iend - 1, iend]
                t = Hm1m1 + Hmm
            end
            # t = iszero(t) ? eps(real(one(t))) : t # introduce a small pertubation for zero shifts
            @mydebug println("block range: $(istart):$(iend) d: ", _fmt_nr(d),
                             " t: ", _fmt_nr(t))

            if shiftmethod == :Francis
                # Run a bulge chase
                if iszero(i % 10)
                    # Vary the shift strategy to avoid dead locks
                    # We use a Wilkinson-like shift as suggested in "Sandia technical report 96-0913J: How the QR algorithm fails to converge and how fix it".

                    @mydebug println("Wilkinson-like shift! Subdiagonal is: ",
                                     _fmt_nr(HH[iend, iend - 1]),
                                     " last subdiagonal is: ",
                                     _fmt_nr(HH[iend-1, iend-2]))
                    _d = t*t - 4d

                    if _d isa Real && _d >= 0
                        # real eigenvalues
                        a = t/2
                        b = sqrt(_d)/2
                        s = a > Hmm ? a - b : a + b
                    else
                        # complex case
                        s = t/2
                    end
                    singleShiftQR!(HH, τ, Z, s, istart, iend)
                else
                    # most of the time use Francis double shifts
                    @mydebug println("Francis double shift! Subdiagonal is: ",
                                     _fmt_nr(HH[iend, iend - 1]),
                                     " last subdiagonal is: ",
                                     _fmt_nr(HH[iend - 1, iend - 2]))
                    doubleShiftQR!(HH, τ, Z, t, d, istart, iend)
                end
            elseif shiftmethod == :Rayleigh
                @mydebug println("Single shift with Rayleigh shift! Subdiagonal is: ",
                                 _fmt_nr(HH[iend, iend - 1]))

                # Run a bulge chase
                singleShiftQR!(HH, τ, Z, Hmm, istart, iend)
            else
                throw(ArgumentError("only support supported shift methods are :Francis (default) and :Rayleigh. You supplied $shiftmethod"))
            end
        end
        if iend <= 2
            if standardize
                if iend == 1
                    putw!(w,HH[iend,iend])
                elseif iend == 2
                    @mydebug println("final std. iend = $iend")
                    H2 = HH[1:2,1:2]
                    G2,w1,w2 = _gs2x2!(H2,2)
                    putw!(w,w2)
                    putw!(w,w1)
                    lmul!(G2,HH)
                    rmul!(view(HH,1:2,:),G2')
                    HH[1:2,1:2] .= H2 # clean
                    Z === nothing || rmul!(Z,G2')
                end
            end
            break
        end
    end

    TT = triu(HH,-1)
    if standardize
        v = w
    else
        v = _geigvals!(TT)
    end
    return Schur{T,typeof(TT)}(TT, Z === nothing ? similar(TT,0,0) : Z, v)
end

# compute Schur decomposition of real 2x2 in standard form
# return corresponding Givens and eigenvalues
# Translated from LAPACK::dlanv2
# Copyright:
# Univ. of Tennessee
# Univ. of California Berkeley
# Univ. of Colorado Denver
# NAG Ltd.
function _gs2x2!(H2::StridedMatrix{T},jj) where {T <: Real}
    a,b,c,d = H2[1,1], H2[1,2], H2[2,1], H2[2,2]
    sgn(x) = (x < 0) ? -one(T) : one(T) # fortran sign differs from Julia
    half = one(T) / 2
    small = 4eps(T) # how big discriminant must be for easy reality check
    if c==0
        cs = one(T)
        sn = zero(T)
    elseif b==0
        # swap rows/cols
        cs = zero(T)
        sn = one(T)
        a,b,c,d = d,-c,zero(T),a
    elseif ((a-d) == 0) && (b*c < 0)
        # nothing to do
        cs = one(T)
        sn = zero(T)
    else
        asubd = a-d
        p = half*asubd
        bcmax = max(abs(b),abs(c))
        bcmis = min(abs(b),abs(c)) * sgn(b) * sgn(c)
        scale = max(abs(p), bcmax)
        z = (p / scale) * p + (bcmax / scale) * bcmis
        # if z is of order machine accuracy: postpone decision
        if z >= small
            # real eigenvalues
            z = p + sqrt(scale) * sqrt(z) * sgn(p)
            a = d + z
            d -= (bcmax / z) * bcmis
            τ = hypot(c,z)
            cs = z / τ
            sn = c / τ
            b -= c
            c = zero(T)
        else
            # complex or almost equal real eigenvalues
            σ = b + c
            τ = hypot(σ, asubd)
            cs = sqrt(half * (one(T) + abs(σ) / τ))
            sn = -(p / (τ * cs)) * sgn(σ)
            # apply rotations
            aa = a*cs + b*sn
            bb = -a*sn + b*cs
            cc = c*cs + d*sn
            dd = -c*sn + d*cs
            a = aa*cs + cc*sn
            b = bb*cs + dd*sn
            c = -aa*sn + cc*cs
            d = -bb*sn + dd*cs
            midad = half * (a+d)
            a = midad
            d = a
            if (c != 0)
                if (b != 0)
                    if b*c >= 0
                        # real eigenvalues
                        sab = sqrt(abs(b))
                        sac = sqrt(abs(c))
                        p = sab*sac*sgn(c)
                        τ = one(T) / sqrt(abs(b+c))
                        a = midad + p
                        d = midad - p
                        b -= c
                        c = 0
                        cs1 = sab*τ
                        sn1 = sac*τ
                        cs, sn = cs*cs1 - sn*sn1, cs*sn1 + sn*cs1
                    end
                else
                    b,c = -c,zero(T)
                    cs,sn = -sn,cs
                end
            end
        end
    end

    if c==0
        w1,w2 = a,d
    else
        rti = sqrt(abs(b))*sqrt(abs(c))
        w1 = a + rti*im
        w2 = d - rti*im
    end
    H2[1,1], H2[1,2], H2[2,1], H2[2,2] = a,b,c,d
    G = Givens(jj-1,jj,cs,sn)
    return G,w1,w2
end

function gschur!(A::StridedMatrix{T}; wantZ::Bool=true, scale::Bool=true,
                 kwargs...) where {T <: Real}
    n = checksquare(A)
    if scale
        scaleA, cscale, anrm = _scale!(A)
    else
        scaleA = false
    end
    H = _hessenberg!(A)
    if wantZ
        Z = _materializeQ(H)
        S = _gschur!(H, Z; kwargs...)
    else
        S = _gschur!(H; kwargs...)
    end
    if scaleA
        safescale!(S.T, cscale, anrm)
        safescale!(S.values, cscale, anrm)
    end
    S
end

function singleShiftQR!(HH::StridedMatrix{T}, τ::Rotation, Z, shift::Number, istart::Integer, iend::Integer) where {T <: Real}
    m = size(HH, 1)
    H11 = HH[istart, istart]
    H21 = HH[istart + 1, istart]
    if m > istart + 1
        Htmp = HH[istart + 2, istart]
        HH[istart + 2, istart] = 0
    end
    G, _ = givens(H11 - shift, H21, istart, istart + 1)
    lmul!(G, view(HH, :, istart:m))
    rmul!(view(HH, 1:min(istart + 2, iend), :), G')
    lmul!(G, τ)
    Z === nothing || rmul!(Z,G')
    for i = istart:iend - 2
        G, _ = givens(HH[i + 1, i], HH[i + 2, i], i + 1, i + 2)
        lmul!(G, view(HH, :, i:m))
        HH[i + 2, i] = Htmp
        if i < iend - 2
            Htmp = HH[i + 3, i + 1]
            HH[i + 3, i + 1] = 0
        end
        rmul!(view(HH, 1:min(i + 3, iend), :), G')
        Z === nothing || rmul!(Z,G')
    end
    return HH
end

function doubleShiftQR!(HH::StridedMatrix{T}, τ::Rotation, Z, shiftTrace::Number, shiftDeterminant::Number, istart::Integer, iend::Integer) where {T <: Real}
    m = size(HH, 1)
    H11 = HH[istart, istart]
    H21 = HH[istart + 1, istart]
    Htmp11 = HH[istart + 2, istart]
    HH[istart + 2, istart] = 0
    if istart + 3 <= m
        Htmp21 = HH[istart + 3, istart]
        HH[istart + 3, istart] = 0
        Htmp22 = HH[istart + 3, istart + 1]
        HH[istart + 3, istart + 1] = 0
    else
        # values doen't matter in this case but variables should be initialized
        Htmp21 = Htmp22 = Htmp11
    end
    G1, r = givens(H11*H11 + HH[istart, istart + 1]*H21 - shiftTrace*H11 + shiftDeterminant, H21*(H11 + HH[istart + 1, istart + 1] - shiftTrace), istart, istart + 1)
    G2, _ = givens(r, H21*HH[istart + 2, istart + 1], istart, istart + 2)
    vHH = view(HH, :, istart:m)
    lmul!(G1, vHH)
    lmul!(G2, vHH)
    vHH = view(HH, 1:min(istart + 3, m), :)
    rmul!(vHH, G1')
    rmul!(vHH, G2')
    lmul!(G1, τ)
    lmul!(G2, τ)
    Z === nothing || rmul!(Z,G1')
    Z === nothing || rmul!(Z,G2')
    for i = istart:iend - 2
        for j = 1:2
            if i + j + 1 > iend break end
            # G, _ = givens(H.H,i+1,i+j+1,i)
            G, _ = givens(HH[i + 1, i], HH[i + j + 1, i], i + 1, i + j + 1)
            lmul!(G, view(HH, :, i:m))
            HH[i + j + 1, i] = Htmp11
            Htmp11 = Htmp21
            # if i + j + 2 <= iend
                # Htmp21 = HH[i + j + 2, i + 1]
                # HH[i + j + 2, i + 1] = 0
            # end
            if i + 4 <= iend
                Htmp22 = HH[i + 4, i + j]
                HH[i + 4, i + j] = 0
            end
            rmul!(view(HH, 1:min(i + j + 2, iend), :), G')
            Z === nothing || rmul!(Z,G')
        end
    end
    return HH
end

# get eigenvalues from a quasitriangular Schur factor
function _geigvals!(HH::StridedMatrix{T}; tol = eps(T)) where {T <: Real}
    # this is used for the non-standard form
    n = size(HH, 1)
    vals = Vector{complex(T)}(undef, n)
    i = 1
    while i < n
        Hii = HH[i, i]
        Hi1i1 = HH[i + 1, i + 1]
        rtest = tol*(abs(Hi1i1) + abs(Hii))
        if abs(HH[i + 1, i]) < rtest
            vals[i] = Hii
            i += 1
        else
            d = Hii*Hi1i1 - HH[i, i + 1]*HH[i + 1, i]
            t = Hii + Hi1i1
            x = 0.5*t
            y = sqrt(complex(x*x - d))
            vals[i] = x + y
            vals[i + 1] = x - y
            i += 2
        end
    end
    if i == n
        vals[i] = HH[n, n]
    end
    return vals
end

include("vectors.jl")
include("triang.jl")

include("norm1est.jl")
include("sylvester.jl")

include("generalized.jl")
include("gsylvester.jl")

include("ordschur.jl")
include("gcondition.jl")

end # module


text:
# first make sure that current working directory is scripts and not Rhone2020
cd(@__DIR__)

using Pkg;
Pkg.activate(joinpath(@__DIR__, ".."))
using Rhone2020;
const R = Rhone2020;

## settings for uncertainty propagation and MCMC:
ind = 1 # set to 1 for faster processing, to 2 for paper
n_iterations = [5*10^4, 10^6][ind] # at least 10^6 for nice distribution in plot
set_n_partcl([1000, 20000][ind]) # fig. 4 used 20000 to make a nice looking histogram.  But that takes time and memory to run!

# Plotting.  Set to `true` for interactive plots
using PyPlot
pygui(false)

# run all scripts for computations
include("read_measurements.jl") # load all data (also runs calibration.jl)
include("derived_quantities.jl") # compute hydraulic gradient, discharge, flow speed, cross-sectional area, friction factor
include("size_evolution_models.jl") # run ct-model and free-gradient model
include("heat_transfer.jl") # compute thermodynamic variables (equilibrium offset-temperature etc.)
include("paper_figures.jl") # plot all figures of the paper and supplements and save them in the products folder
                            # also writes tables and CSV files.
                            # Also runs files: overview-map.jl


text:
module JuCP

using JuMP
using ConstraintProgrammingExtensions

const CP = ConstraintProgrammingExtensions

include("sets.jl")

end


text:
function next_vector(obj::RandomVectorGenerator)
    return jcall(obj, "nextVector", Vector{jdouble}, ())
end


text:
type class_nlo <: internal_AbstractNLPEvaluator
  # linear program with non-linear objective
  # min f(x)
  # A*x = b
  # x >= 0

	_n::Int64 # number of variables
	_m::Int64 # number of constraints
	_A::SparseMatrixCSC{Float64,Int64}
	_b::Array{Float64,1}
  obj::class_nl_function

	function class_nlo(A::SparseMatrixCSC{Float64,Int64},b::Array{Float64,1}, obj::class_nl_function)
		(m, n) = size(A)
		return new(n,m,A,b,obj);
	end
end

################
# METHODS
################

# evaluate objective
function internal_eval_c(nlo::class_nlo, x::Array{Float64,1})
    return nlo.obj.value(x);
end

# evalutate constraints
function internal_eval_a(nlo::class_nlo, x::Array{Float64,1})
    return nlo._A * x - nlo._b;
end

# evaluate gradient of constraints
function internal_eval_jac_a(nlo::class_nlo, x::Array{Float64,1}) # J
    return nlo._A;
end

# hessian of lagrangian
function internal_eval_hesslag_prod(nlo::class_nlo, x::Array{Float64,1}, y::Array{Float64,1})
    return nlo.obj.hessian(x)
end

# gradient of lagrangian
function internal_eval_gradlag(nlo::class_nlo, x::Array{Float64,1}, y::Array{Float64,1})
    return internal_eval_gradc(nlo, x) - nlo._A' * y;
end

# gradient of f
function internal_eval_gradc(nlo::class_nlo, x::Array{Float64,1})
    return nlo.obj.gradient(x);
end

# (\nabla g * x - g)
function internal_eval_b(nlo::class_nlo, x::Array{Float64,1})
    return nlo._b;
end

function n(nlo::class_nlo)
    return nlo._n;
end

function m(nlo::class_nlo)
    return nlo._m;
end

max_forks_repo_name: ohinder/Everest


text:
struct Space{S}
    function Space{S}() where {S}
        new{S::Tuple{Vararg{Union{Int, Symmetry}}}}()
    end
end

@pure Space(dims::Vararg{Union{Int, Symmetry}}) = Space{dims}()
@pure Space(dims::Tuple) = Space{dims}()

_construct(x::Int) = x
@pure _construct(::Type{Symmetry{S}}) where {S} = Symmetry{S}()
@pure function Space(::Type{S}) where {S <: Tuple}
    check_size_parameters(S)
    dims = map(_construct, tuple(S.parameters...))
    Space(dims)
end

_ncomponents(x::Int) = x
_ncomponents(x::Symmetry) = ncomponents(x)
@pure ncomponents(::Space{S}) where {S} = prod(_ncomponents, S)

@pure Base.Dims(::Space{S}) where {S} = flatten_tuple(map(Dims, S))
@pure Base.Tuple(::Space{S}) where {S} = S

@pure Base.length(s::Space) = length(Tuple(s))
Base.getindex(s::Space, i::Int) = Tuple(s)[i]

@pure tensorsize(s::Space) = Dims(s)
@pure tensororder(s::Space) = length(tensorsize(s))
# don't allow to use `size` and `ndims` because their names are confusing.
Base.size(s::Space) = throw(ArgumentError("use `tensorsize` to get size of a tensor instead of `size`"))
Base.ndims(s::Space) = throw(ArgumentError("use `tensororder` to get order of a tensor instead of `ndims`"))

function Base.show(io::IO, ::Space{S}) where {S}
    print(io, "Space", S)
end

# dropfirst/droplast
dropfirst() = error()
dropfirst(x::Int) = ()
@pure dropfirst(::Symmetry{NTuple{2, dim}}) where {dim} = (dim,)
@pure dropfirst(::Symmetry{NTuple{order, dim}}) where {order, dim} = (Symmetry{NTuple{order-1, dim}}(),)
_dropfirst(x, ys...) = (dropfirst(x)..., ys...)
@pure dropfirst(::Space{S}) where {S} = Space(_dropfirst(S...))
@pure droplast(::Space{S}) where {S} = Space(reverse(_dropfirst(reverse(S)...)))
for op in (:dropfirst, :droplast)
    @eval begin
        $op(x::Space, ::Val{0}) = x
        $op(x::Space, ::Val{N}) where {N} = $op($op(x), Val(N-1))
    end
end

# contractions
@pure function contraction(x::Space, y::Space, ::Val{N}) where {N}
    if !(0 ≤ N ≤ tensororder(x) && 0 ≤ N ≤ tensororder(y) && tensorsize(x)[end-N+1:end] === tensorsize(y)[1:N])
        throw(DimensionMismatch("dimensions must match"))
    end
    otimes(droplast(x, Val(N)), dropfirst(y, Val(N)))
end
@pure otimes(x::Space) = x
@pure otimes(x::Space, y::Space) = Space(Tuple(x)..., Tuple(y)...)
@pure otimes(x::Space, y::Space, z::Space...) = otimes(otimes(x, y), z...)
@pure dot(x::Space, y::Space) = contraction(x, y, Val(1))
@pure double_contraction(x::Space, y::Space) = contraction(x, y, Val(2))

# promote_space
promote_space(x::Space) = x
@generated function promote_space(x::Space{S1}, y::Space{S2}) where {S1, S2}
    tensorsize(Space(S1)) == tensorsize(Space(S2)) || return :(throw(DimensionMismatch("dimensions must match")))
    :(Space($(_promote_space(S1, S2, ()))))
end
promote_space(x::Space, y::Space, z::Space...) = promote_space(promote_space(x, y), z...)
## helper functions
_promote_space(x::Tuple{}, y::Tuple{}, promoted::Tuple) = promoted
function _promote_space(x::Tuple, y::Tuple, promoted::Tuple)
    x1 = x[1]
    y1 = y[1]
    if x1 == y1
        _promote_space(Base.tail(x), Base.tail(y), (promoted..., x1))
    else
        x1_len = length(x1)
        y1_len = length(y1)
        if x1_len < y1_len
            common = promote_space(Space(x1),
                                   droplast(Space(y1), Val(y1_len - x1_len))) |> Tuple
            _promote_space(Base.tail(x),
                           Tuple(dropfirst(Space(y), Val(x1_len))),
                           (promoted..., only(common)))
        elseif length(x1) > length(y1)
            common = promote_space(droplast(Space(x1), Val(x1_len - y1_len)),
                                   Space(y1)) |> Tuple
            _promote_space(Tuple(dropfirst(Space(x), Val(y1_len))),
                           Base.tail(y),
                           (promoted..., only(common)))
        else
            error()
        end
    end
end

# tensortype
_typeof(x::Int) = x
_typeof(x::Symmetry) = typeof(x)
@pure function tensortype(x::Space)
    Tensor{Tuple{map(_typeof, Tuple(x))...}, T, tensororder(x), ncomponents(x)} where {T}
end

# LinearIndices/CartesianIndices
for IndicesType in (LinearIndices, CartesianIndices)
    @eval (::Type{$IndicesType})(x::Space) = $IndicesType(tensorsize(x))
end


text:
mutable struct RenderState
  mode::Any
  window::Any
  texture::Any
  sprite::Any
  background::Any
end

function RenderState(; width=32,height=32, background=sfColor_fromRGBA(0,0,0,255))
  mode = sfVideoMode(width, height, 32)

  window = sfRenderWindow_create(mode, "Julia Fractal Renderer", sfResize | sfClose, C_NULL)
  @assert window != C_NULL

  texture = sfTexture_create(width, height)
  @assert texture != C_NULL

  sprite = sfSprite_create()
  sfSprite_setTexture(sprite, texture, sfTrue)
  @assert sprite != C_NULL

  obj = RenderState(mode, window, texture, sprite, background)
  finalizer(destroy!, obj)
  return obj
end

function update!(simState::State, renderState::RenderState)
  sfTexture_updateFromPixels(renderState.texture, map((val, n) -> getColor(val, n, simState.iterations), simState.field, simState.fieldIterations), simState.res..., 0, 0)
end

function render!(simState::State, renderState::RenderState)
  sfRenderWindow_clear(renderState.window, renderState.background)
  sfRenderWindow_drawSprite(renderState.window, renderState.sprite, C_NULL)
  sfRenderWindow_display(renderState.window)
end

function destroy!(renderState::RenderState)
  sfSprite_destroy(renderState.sprite)
  sfTexture_destroy(renderState.texture)
  sfRenderWindow_destroy(renderState.window)
end

max_stars_repo_name: benchislett/Skunkworks

max_issues_repo_name: benchislett/Skunkworks

max_forks_repo_name: benchislett/Skunkworks


text:
"""
    convorder(alg::AbstractIteratedIntegralAlgorithm)

Returns the convergence order of the algorithm w.r.t. the truncation parameter.

See also: [`errcoeff`](@ref)
"""
function convorder end

"""
    errcoeff(dim, stepsize, alg, norm)
    errcoeff(dim, q_12, stepsize, alg, norm)

Returns the coefficient of the truncation parameter in the error estimate.
I.e. the error estimate is of the form

```math
\\lVert I(h)-\\tilde{I}^{(p)}(h) \\rVert_* ≤ \\mathrm{errcoeff}(m,h) \\cdot p^{-γ}
```

where the norm is given by `norm`, the approximation ``\\tilde{I}^{(p)}`` is
calculated using `alg` and ``γ`` is the order of convergence given by [`convorder(alg)`](@ref).
"""
function errcoeff end

"""
    terms_needed(dim, stepsize, eps, alg, norm)

Returns the number of terms in the approximating sum that is needed to ensure an error
in the given norm of at most `eps`.
This depends on the dimension of the Wiener process `dim`, the current stepsize and the chosen algorithm.

See also: [`AbstractIteratedIntegralAlgorithm`](@ref), [`AbstractErrorNorm`](@ref)

# Examples
```jldoctest; setup=:(using LevyArea)
julia> h = 1/128;

julia> terms_needed(10, h, h^(3/2), Milstein(), MaxL2())
7
```

# Implementation
New algorithms should only have to implement [`errcoeff`](@ref) and [`convorder`](@ref).
"""
function terms_needed(dim, stepsize, eps, alg::AbstractIteratedIntegralAlgorithm, norm::AbstractErrorNorm)
    ceil(Int64, (errcoeff(dim, stepsize, alg, norm)/eps)^(1//convorder(alg)) )
end

"""
    terms_needed(dim, q_12, stepsize, eps, alg, norm)

Used for finite-dimensional approximations of a Q-Wiener process with covariance matrix
``Q = Q^\\frac{1}{2}*Q^\\frac{1}{2}``. Here `q_12` is a vector of the eigenvalues of ``Q^\\frac{1}{2}``;
the square root of the covariance matrix. Equivalently these are the square roots of the eigenvalues of ``Q``.

# Examples
```jldoctest; setup=:(using LevyArea)
julia> h = 1/128;

julia> dim = 10;

julia> q = [1/k^2 for k=1:dim];

julia> terms_needed(dim, sqrt.(q), h, h^(3/2), Milstein(), FrobeniusL2())
9
```
"""
function terms_needed(dim, q_12, stepsize, eps, alg::AbstractIteratedIntegralAlgorithm, norm::AbstractErrorNorm)
    length(q_12) == dim || throw(ArgumentError("Length of q_12 must be equal to the dimension."))
    ceil(Int64, (errcoeff(dim, q_12, stepsize, alg, norm)/eps)^(1//convorder(alg)) )
end

"""
    norv(dim, n, alg::AbstractIteratedIntegralAlgorithm)

Returns the number of random numbers needed to simulate the iterated integrals
for a Wiener process of dimension `dim` and with truncation parameter `n`.
"""
norv(dim, n, alg::AbstractIteratedIntegralAlgorithm)

"""
    effective_cost(dim, stepsize, eps, alg, norm)
    effective_cost(dim, q_12, stepsize, eps, alg, norm)

Returns the number of random numbers needed to simulate the iterated integrals
with the given parameters.
"""
function effective_cost(dim, stepsize, eps, alg, norm)
    norv(dim, terms_needed(dim, stepsize, eps, alg, norm), alg)
end
function effective_cost(dim, q_12, stepsize, eps, alg, norm)
    norv(dim, terms_needed(dim, q_12, stepsize, eps, alg, norm), alg)
end

"""
    optimal_algorithm(dim, stepsize, eps=stepsize^(3/2), norm=MaxL2())
    optimal_algorithm(dim, q_12, stepsize, eps, norm=FrobeniusL2())

Returns the optimal algorithm for the given parameters,
i.e. the algorithm that needs to simulate the fewest random numbers
to achieve the desired precision `eps`.

# Examples
```jldoctest; setup=:(using LevyArea)
julia> h = 1/128;

julia> optimal_algorithm(10, h, h^(3/2), MaxL2())
MronRoe()

julia> optimal_algorithm(10, 1.0./(1:10).^2, h, h^(3/2), FrobeniusL2())
Milstein()
```
"""
function optimal_algorithm(dim, stepsize, eps=stepsize^(3/2), norm::AbstractErrorNorm=MaxL2())
    ind = argmin([effective_cost(dim, stepsize, eps, alg, norm) for alg ∈ ITER_INT_ALGS])
    return ITER_INT_ALGS[ind]
end
function optimal_algorithm(dim, q_12, stepsize, eps, norm::AbstractErrorNorm=FrobeniusL2())
    ind = argmin([effective_cost(dim, q_12, stepsize, eps, alg, norm) for alg ∈ ITER_INT_ALGS])
    return ITER_INT_ALGS[ind]
end

max_stars_repo_name: mrbauff/MasterThesis

max_issues_repo_name: mrbauff/MasterThesis


text:
function resume_simulation_from_file(
    input_file::String  ,
    output_file::String ;
    t_end::Float64 = 0.0
    )

    if(t_end == 0)
        println("ERROR: please specify a final time by setting 't_end' argument")
        return
    end

    p0 = Param(

        resume_simulation   = :true,
        input_file_name     = input_file,

        t_sim_final = t_end ,

        folder          = ""   ,
        fname           = output_file
    )

    solve_wave_equation_2D(p0)

end


text:
module LibKLU
const libklu = :libklu

mutable struct klu_symbolic
    symmetry::Cdouble
    est_flops::Cdouble
    lnz::Cdouble
    unz::Cdouble
    Lnz::Ptr{Cdouble}
    n::Cint
    nz::Cint
    P::Ptr{Cint}
    Q::Ptr{Cint}
    R::Ptr{Cint}
    nzoff::Cint
    nblocks::Cint
    maxblock::Cint
    ordering::Cint
    do_btf::Cint
    structural_rank::Cint
    klu_symbolic() = new()
end

mutable struct klu_l_symbolic
    symmetry::Cdouble
    est_flops::Cdouble
    lnz::Cdouble
    unz::Cdouble
    Lnz::Ptr{Cdouble}
    n::Clong
    nz::Clong
    P::Ptr{Clong}
    Q::Ptr{Clong}
    R::Ptr{Clong}
    nzoff::Clong
    nblocks::Clong
    maxblock::Clong
    ordering::Clong
    do_btf::Clong
    structural_rank::Clong
    klu_l_symbolic() = new()
end

mutable struct klu_numeric
    n::Cint
    nblocks::Cint
    lnz::Cint
    unz::Cint
    max_lnz_block::Cint
    max_unz_block::Cint
    Pnum::Ptr{Cint}
    Pinv::Ptr{Cint}
    Lip::Ptr{Cint}
    Uip::Ptr{Cint}
    Llen::Ptr{Cint}
    Ulen::Ptr{Cint}
    LUbx::Ptr{Ptr{Cvoid}}
    LUsize::Ptr{Csize_t}
    Udiag::Ptr{Cvoid}
    Rs::Ptr{Cdouble}
    worksize::Csize_t
    Work::Ptr{Cvoid}
    Xwork::Ptr{Cvoid}
    Iwork::Ptr{Cint}
    Offp::Ptr{Cint}
    Offi::Ptr{Cint}
    Offx::Ptr{Cvoid}
    nzoff::Cint
    klu_numeric() = new()
end

mutable struct klu_l_numeric
    n::Clong
    nblocks::Clong
    lnz::Clong
    unz::Clong
    max_lnz_block::Clong
    max_unz_block::Clong
    Pnum::Ptr{Clong}
    Pinv::Ptr{Clong}
    Lip::Ptr{Clong}
    Uip::Ptr{Clong}
    Llen::Ptr{Clong}
    Ulen::Ptr{Clong}
    LUbx::Ptr{Ptr{Cvoid}}
    LUsize::Ptr{Csize_t}
    Udiag::Ptr{Cvoid}
    Rs::Ptr{Cdouble}
    worksize::Csize_t
    Work::Ptr{Cvoid}
    Xwork::Ptr{Cvoid}
    Iwork::Ptr{Clong}
    Offp::Ptr{Clong}
    Offi::Ptr{Clong}
    Offx::Ptr{Cvoid}
    nzoff::Clong
    klu_l_numeric() = new()
end

mutable struct klu_common_struct
    tol::Cdouble
    memgrow::Cdouble
    initmem_amd::Cdouble
    initmem::Cdouble
    maxwork::Cdouble
    btf::Cint
    ordering::Cint
    scale::Cint
    user_order::Ptr{Cvoid}
    user_data::Ptr{Cvoid}
    halt_if_singular::Cint
    status::Cint
    nrealloc::Cint
    structural_rank::Cint
    numerical_rank::Cint
    singular_col::Cint
    noffdiag::Cint
    flops::Cdouble
    rcond::Cdouble
    condest::Cdouble
    rgrowth::Cdouble
    work::Cdouble
    memusage::Csize_t
    mempeak::Csize_t
    klu_common_struct() = new()
end

const klu_common = klu_common_struct

mutable struct klu_l_common_struct
    tol::Cdouble
    memgrow::Cdouble
    initmem_amd::Cdouble
    initmem::Cdouble
    maxwork::Cdouble
    btf::Clong
    ordering::Clong
    scale::Clong
    user_order::Ptr{Cvoid}
    user_data::Ptr{Cvoid}
    halt_if_singular::Clong
    status::Clong
    nrealloc::Clong
    structural_rank::Clong
    numerical_rank::Clong
    singular_col::Clong
    noffdiag::Clong
    flops::Cdouble
    rcond::Cdouble
    condest::Cdouble
    rgrowth::Cdouble
    work::Cdouble
    memusage::Csize_t
    mempeak::Csize_t
    klu_l_common_struct() = new()
end

const klu_l_common = klu_l_common_struct

function klu_defaults(Common)
    @ccall libklu.klu_defaults(Common::Ptr{klu_common})::Cint
end

function klu_l_defaults(Common)
    @ccall libklu.klu_l_defaults(Common::Ptr{klu_l_common})::Clong
end

function klu_analyze(n, Ap, Ai, Common)
    @ccall libklu.klu_analyze(n::Cint, Ap::Ptr{Cint}, Ai::Ptr{Cint}, Common::Ptr{klu_common})::Ptr{klu_symbolic}
end

function klu_l_analyze(arg1, arg2, arg3, Common)
    @ccall libklu.klu_l_analyze(arg1::Clong, arg2::Ptr{Clong}, arg3::Ptr{Clong}, Common::Ptr{klu_l_common})::Ptr{klu_l_symbolic}
end

function klu_analyze_given(n, Ap, Ai, P, Q, Common)
    @ccall libklu.klu_analyze_given(n::Cint, Ap::Ptr{Cint}, Ai::Ptr{Cint}, P::Ptr{Cint}, Q::Ptr{Cint}, Common::Ptr{klu_common})::Ptr{klu_symbolic}
end

function klu_l_analyze_given(arg1, arg2, arg3, arg4, arg5, arg6)
    @ccall libklu.klu_l_analyze_given(arg1::Clong, arg2::Ptr{Clong}, arg3::Ptr{Clong}, arg4::Ptr{Clong}, arg5::Ptr{Clong}, arg6::Ptr{klu_l_common})::Ptr{klu_l_symbolic}
end

function klu_factor(Ap, Ai, Ax, Symbolic, Common)
    @ccall libklu.klu_factor(Ap::Ptr{Cint}, Ai::Ptr{Cint}, Ax::Ptr{Cdouble}, Symbolic::Ptr{klu_symbolic}, Common::Ptr{klu_common})::Ptr{klu_numeric}
end

function klu_z_factor(Ap, Ai, Ax, Symbolic, Common)
    @ccall libklu.klu_z_factor(Ap::Ptr{Cint}, Ai::Ptr{Cint}, Ax::Ptr{Cdouble}, Symbolic::Ptr{klu_symbolic}, Common::Ptr{klu_common})::Ptr{klu_numeric}
end

function klu_l_factor(arg1, arg2, arg3, arg4, arg5)
    @ccall libklu.klu_l_factor(arg1::Ptr{Clong}, arg2::Ptr{Clong}, arg3::Ptr{Cdouble}, arg4::Ptr{klu_l_symbolic}, arg5::Ptr{klu_l_common})::Ptr{klu_l_numeric}
end

function klu_zl_factor(arg1, arg2, arg3, arg4, arg5)
    @ccall libklu.klu_zl_factor(arg1::Ptr{Clong}, arg2::Ptr{Clong}, arg3::Ptr{Cdouble}, arg4::Ptr{klu_l_symbolic}, arg5::Ptr{klu_l_common})::Ptr{klu_l_numeric}
end

function klu_solve(Symbolic, Numeric, ldim, nrhs, B, Common)
    @ccall libklu.klu_solve(Symbolic::Ptr{klu_symbolic}, Numeric::Ptr{klu_numeric}, ldim::Cint, nrhs::Cint, B::Ptr{Cdouble}, Common::Ptr{klu_common})::Cint
end

function klu_z_solve(Symbolic, Numeric, ldim, nrhs, B, Common)
    @ccall libklu.klu_z_solve(Symbolic::Ptr{klu_symbolic}, Numeric::Ptr{klu_numeric}, ldim::Cint, nrhs::Cint, B::Ptr{Cdouble}, Common::Ptr{klu_common})::Cint
end

function klu_l_solve(arg1, arg2, arg3, arg4, arg5, arg6)
    @ccall libklu.klu_l_solve(arg1::Ptr{klu_l_symbolic}, arg2::Ptr{klu_l_numeric}, arg3::Clong, arg4::Clong, arg5::Ptr{Cdouble}, arg6::Ptr{klu_l_common})::Clong
end

function klu_zl_solve(arg1, arg2, arg3, arg4, arg5, arg6)
    @ccall libklu.klu_zl_solve(arg1::Ptr{klu_l_symbolic}, arg2::Ptr{klu_l_numeric}, arg3::Clong, arg4::Clong, arg5::Ptr{Cdouble}, arg6::Ptr{klu_l_common})::Clong
end

function klu_tsolve(Symbolic, Numeric, ldim, nrhs, B, Common)
    @ccall libklu.klu_tsolve(Symbolic::Ptr{klu_symbolic}, Numeric::Ptr{klu_numeric}, ldim::Cint, nrhs::Cint, B::Ptr{Cdouble}, Common::Ptr{klu_common})::Cint
end

function klu_z_tsolve(Symbolic, Numeric, ldim, nrhs, B, conj_solve, Common)
    @ccall libklu.klu_z_tsolve(Symbolic::Ptr{klu_symbolic}, Numeric::Ptr{klu_numeric}, ldim::Cint, nrhs::Cint, B::Ptr{Cdouble}, conj_solve::Cint, Common::Ptr{klu_common})::Cint
end

function klu_l_tsolve(arg1, arg2, arg3, arg4, arg5, arg6)
    @ccall libklu.klu_l_tsolve(arg1::Ptr{klu_l_symbolic}, arg2::Ptr{klu_l_numeric}, arg3::Clong, arg4::Clong, arg5::Ptr{Cdouble}, arg6::Ptr{klu_l_common})::Clong
end

function klu_zl_tsolve(arg1, arg2, arg3, arg4, arg5, arg6, arg7)
    @ccall libklu.klu_zl_tsolve(arg1::Ptr{klu_l_symbolic}, arg2::Ptr{klu_l_numeric}, arg3::Clong, arg4::Clong, arg5::Ptr{Cdouble}, arg6::Clong, arg7::Ptr{klu_l_common})::Clong
end

function klu_refactor(Ap, Ai, Ax, Symbolic, Numeric, Common)
    @ccall libklu.klu_refactor(Ap::Ptr{Cint}, Ai::Ptr{Cint}, Ax::Ptr{Cdouble}, Symbolic::Ptr{klu_symbolic}, Numeric::Ptr{klu_numeric}, Common::Ptr{klu_common})::Cint
end

function klu_z_refactor(Ap, Ai, Ax, Symbolic, Numeric, Common)
    @ccall libklu.klu_z_refactor(Ap::Ptr{Cint}, Ai::Ptr{Cint}, Ax::Ptr{Cdouble}, Symbolic::Ptr{klu_symbolic}, Numeric::Ptr{klu_numeric}, Common::Ptr{klu_common})::Cint
end

function klu_l_refactor(arg1, arg2, arg3, arg4, arg5, arg6)
    @ccall libklu.klu_l_refactor(arg1::Ptr{Clong}, arg2::Ptr{Clong}, arg3::Ptr{Cdouble}, arg4::Ptr{klu_l_symbolic}, arg5::Ptr{klu_l_numeric}, arg6::Ptr{klu_l_common})::Clong
end

function klu_zl_refactor(arg1, arg2, arg3, arg4, arg5, arg6)
    @ccall libklu.klu_zl_refactor(arg1::Ptr{Clong}, arg2::Ptr{Clong}, arg3::Ptr{Cdouble}, arg4::Ptr{klu_l_symbolic}, arg5::Ptr{klu_l_numeric}, arg6::Ptr{klu_l_common})::Clong
end

function klu_free_symbolic(Symbolic, Common)
    @ccall libklu.klu_free_symbolic(Symbolic::Ptr{Ptr{klu_symbolic}}, Common::Ptr{klu_common})::Cint
end

function klu_l_free_symbolic(arg1, arg2)
    @ccall libklu.klu_l_free_symbolic(arg1::Ptr{Ptr{klu_l_symbolic}}, arg2::Ptr{klu_l_common})::Clong
end

function klu_free_numeric(Numeric, Common)
    @ccall libklu.klu_free_numeric(Numeric::Ptr{Ptr{klu_numeric}}, Common::Ptr{klu_common})::Cint
end

function klu_z_free_numeric(Numeric, Common)
    @ccall libklu.klu_z_free_numeric(Numeric::Ptr{Ptr{klu_numeric}}, Common::Ptr{klu_common})::Cint
end

function klu_l_free_numeric(arg1, arg2)
    @ccall libklu.klu_l_free_numeric(arg1::Ptr{Ptr{klu_l_numeric}}, arg2::Ptr{klu_l_common})::Clong
end

function klu_zl_free_numeric(arg1, arg2)
    @ccall libklu.klu_zl_free_numeric(arg1::Ptr{Ptr{klu_l_numeric}}, arg2::Ptr{klu_l_common})::Clong
end

function klu_sort(Symbolic, Numeric, Common)
    @ccall libklu.klu_sort(Symbolic::Ptr{klu_symbolic}, Numeric::Ptr{klu_numeric}, Common::Ptr{klu_common})::Cint
end

function klu_z_sort(Symbolic, Numeric, Common)
    @ccall libklu.klu_z_sort(Symbolic::Ptr{klu_symbolic}, Numeric::Ptr{klu_numeric}, Common::Ptr{klu_common})::Cint
end

function klu_l_sort(arg1, arg2, arg3)
    @ccall libklu.klu_l_sort(arg1::Ptr{klu_l_symbolic}, arg2::Ptr{klu_l_numeric}, arg3::Ptr{klu_l_common})::Clong
end

function klu_zl_sort(arg1, arg2, arg3)
    @ccall libklu.klu_zl_sort(arg1::Ptr{klu_l_symbolic}, arg2::Ptr{klu_l_numeric}, arg3::Ptr{klu_l_common})::Clong
end

function klu_flops(Symbolic, Numeric, Common)
    @ccall libklu.klu_flops(Symbolic::Ptr{klu_symbolic}, Numeric::Ptr{klu_numeric}, Common::Ptr{klu_common})::Cint
end

function klu_z_flops(Symbolic, Numeric, Common)
    @ccall libklu.klu_z_flops(Symbolic::Ptr{klu_symbolic}, Numeric::Ptr{klu_numeric}, Common::Ptr{klu_common})::Cint
end

function klu_l_flops(arg1, arg2, arg3)
    @ccall libklu.klu_l_flops(arg1::Ptr{klu_l_symbolic}, arg2::Ptr{klu_l_numeric}, arg3::Ptr{klu_l_common})::Clong
end

function klu_zl_flops(arg1, arg2, arg3)
    @ccall libklu.klu_zl_flops(arg1::Ptr{klu_l_symbolic}, arg2::Ptr{klu_l_numeric}, arg3::Ptr{klu_l_common})::Clong
end

function klu_rgrowth(Ap, Ai, Ax, Symbolic, Numeric, Common)
    @ccall libklu.klu_rgrowth(Ap::Ptr{Cint}, Ai::Ptr{Cint}, Ax::Ptr{Cdouble}, Symbolic::Ptr{klu_symbolic}, Numeric::Ptr{klu_numeric}, Common::Ptr{klu_common})::Cint
end

function klu_z_rgrowth(Ap, Ai, Ax, Symbolic, Numeric, Common)
    @ccall libklu.klu_z_rgrowth(Ap::Ptr{Cint}, Ai::Ptr{Cint}, Ax::Ptr{Cdouble}, Symbolic::Ptr{klu_symbolic}, Numeric::Ptr{klu_numeric}, Common::Ptr{klu_common})::Cint
end

function klu_l_rgrowth(arg1, arg2, arg3, arg4, arg5, arg6)
    @ccall libklu.klu_l_rgrowth(arg1::Ptr{Clong}, arg2::Ptr{Clong}, arg3::Ptr{Cdouble}, arg4::Ptr{klu_l_symbolic}, arg5::Ptr{klu_l_numeric}, arg6::Ptr{klu_l_common})::Clong
end

function klu_zl_rgrowth(arg1, arg2, arg3, arg4, arg5, arg6)
    @ccall libklu.klu_zl_rgrowth(arg1::Ptr{Clong}, arg2::Ptr{Clong}, arg3::Ptr{Cdouble}, arg4::Ptr{klu_l_symbolic}, arg5::Ptr{klu_l_numeric}, arg6::Ptr{klu_l_common})::Clong
end

function klu_condest(Ap, Ax, Symbolic, Numeric, Common)
    @ccall libklu.klu_condest(Ap::Ptr{Cint}, Ax::Ptr{Cdouble}, Symbolic::Ptr{klu_symbolic}, Numeric::Ptr{klu_numeric}, Common::Ptr{klu_common})::Cint
end

function klu_z_condest(Ap, Ax, Symbolic, Numeric, Common)
    @ccall libklu.klu_z_condest(Ap::Ptr{Cint}, Ax::Ptr{Cdouble}, Symbolic::Ptr{klu_symbolic}, Numeric::Ptr{klu_numeric}, Common::Ptr{klu_common})::Cint
end

function klu_l_condest(arg1, arg2, arg3, arg4, arg5)
    @ccall libklu.klu_l_condest(arg1::Ptr{Clong}, arg2::Ptr{Cdouble}, arg3::Ptr{klu_l_symbolic}, arg4::Ptr{klu_l_numeric}, arg5::Ptr{klu_l_common})::Clong
end

function klu_zl_condest(arg1, arg2, arg3, arg4, arg5)
    @ccall libklu.klu_zl_condest(arg1::Ptr{Clong}, arg2::Ptr{Cdouble}, arg3::Ptr{klu_l_symbolic}, arg4::Ptr{klu_l_numeric}, arg5::Ptr{klu_l_common})::Clong
end

function klu_rcond(Symbolic, Numeric, Common)
    @ccall libklu.klu_rcond(Symbolic::Ptr{klu_symbolic}, Numeric::Ptr{klu_numeric}, Common::Ptr{klu_common})::Cint
end

function klu_z_rcond(Symbolic, Numeric, Common)
    @ccall libklu.klu_z_rcond(Symbolic::Ptr{klu_symbolic}, Numeric::Ptr{klu_numeric}, Common::Ptr{klu_common})::Cint
end

function klu_l_rcond(arg1, arg2, arg3)
    @ccall libklu.klu_l_rcond(arg1::Ptr{klu_l_symbolic}, arg2::Ptr{klu_l_numeric}, arg3::Ptr{klu_l_common})::Clong
end

function klu_zl_rcond(arg1, arg2, arg3)
    @ccall libklu.klu_zl_rcond(arg1::Ptr{klu_l_symbolic}, arg2::Ptr{klu_l_numeric}, arg3::Ptr{klu_l_common})::Clong
end

function klu_scale(scale, n, Ap, Ai, Ax, Rs, W, Common)
    @ccall libklu.klu_scale(scale::Cint, n::Cint, Ap::Ptr{Cint}, Ai::Ptr{Cint}, Ax::Ptr{Cdouble}, Rs::Ptr{Cdouble}, W::Ptr{Cint}, Common::Ptr{klu_common})::Cint
end

function klu_z_scale(scale, n, Ap, Ai, Ax, Rs, W, Common)
    @ccall libklu.klu_z_scale(scale::Cint, n::Cint, Ap::Ptr{Cint}, Ai::Ptr{Cint}, Ax::Ptr{Cdouble}, Rs::Ptr{Cdouble}, W::Ptr{Cint}, Common::Ptr{klu_common})::Cint
end

function klu_l_scale(arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8)
    @ccall libklu.klu_l_scale(arg1::Clong, arg2::Clong, arg3::Ptr{Clong}, arg4::Ptr{Clong}, arg5::Ptr{Cdouble}, arg6::Ptr{Cdouble}, arg7::Ptr{Clong}, arg8::Ptr{klu_l_common})::Clong
end

function klu_zl_scale(arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8)
    @ccall libklu.klu_zl_scale(arg1::Clong, arg2::Clong, arg3::Ptr{Clong}, arg4::Ptr{Clong}, arg5::Ptr{Cdouble}, arg6::Ptr{Cdouble}, arg7::Ptr{Clong}, arg8::Ptr{klu_l_common})::Clong
end

function klu_extract(Numeric, Symbolic, Lp, Li, Lx, Up, Ui, Ux, Fp, Fi, Fx, P, Q, Rs, R, Common)
    @ccall libklu.klu_extract(Numeric::Ptr{klu_numeric}, Symbolic::Ptr{klu_symbolic}, Lp::Ptr{Cint}, Li::Ptr{Cint}, Lx::Ptr{Cdouble}, Up::Ptr{Cint}, Ui::Ptr{Cint}, Ux::Ptr{Cdouble}, Fp::Ptr{Cint}, Fi::Ptr{Cint}, Fx::Ptr{Cdouble}, P::Ptr{Cint}, Q::Ptr{Cint}, Rs::Ptr{Cdouble}, R::Ptr{Cint}, Common::Ptr{klu_common})::Cint
end

function klu_z_extract(Numeric, Symbolic, Lp, Li, Lx, Lz, Up, Ui, Ux, Uz, Fp, Fi, Fx, Fz, P, Q, Rs, R, Common)
    @ccall libklu.klu_z_extract(Numeric::Ptr{klu_numeric}, Symbolic::Ptr{klu_symbolic}, Lp::Ptr{Cint}, Li::Ptr{Cint}, Lx::Ptr{Cdouble}, Lz::Ptr{Cdouble}, Up::Ptr{Cint}, Ui::Ptr{Cint}, Ux::Ptr{Cdouble}, Uz::Ptr{Cdouble}, Fp::Ptr{Cint}, Fi::Ptr{Cint}, Fx::Ptr{Cdouble}, Fz::Ptr{Cdouble}, P::Ptr{Cint}, Q::Ptr{Cint}, Rs::Ptr{Cdouble}, R::Ptr{Cint}, Common::Ptr{klu_common})::Cint
end

function klu_l_extract(arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9, arg10, arg11, arg12, arg13, arg14, arg15, arg16)
    @ccall libklu.klu_l_extract(arg1::Ptr{klu_l_numeric}, arg2::Ptr{klu_l_symbolic}, arg3::Ptr{Clong}, arg4::Ptr{Clong}, arg5::Ptr{Cdouble}, arg6::Ptr{Clong}, arg7::Ptr{Clong}, arg8::Ptr{Cdouble}, arg9::Ptr{Clong}, arg10::Ptr{Clong}, arg11::Ptr{Cdouble}, arg12::Ptr{Clong}, arg13::Ptr{Clong}, arg14::Ptr{Cdouble}, arg15::Ptr{Clong}, arg16::Ptr{klu_l_common})::Clong
end

function klu_zl_extract(arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9, arg10, arg11, arg12, arg13, arg14, arg15, arg16, arg17, arg18, arg19)
    @ccall libklu.klu_zl_extract(arg1::Ptr{klu_l_numeric}, arg2::Ptr{klu_l_symbolic}, arg3::Ptr{Clong}, arg4::Ptr{Clong}, arg5::Ptr{Cdouble}, arg6::Ptr{Cdouble}, arg7::Ptr{Clong}, arg8::Ptr{Clong}, arg9::Ptr{Cdouble}, arg10::Ptr{Cdouble}, arg11::Ptr{Clong}, arg12::Ptr{Clong}, arg13::Ptr{Cdouble}, arg14::Ptr{Cdouble}, arg15::Ptr{Clong}, arg16::Ptr{Clong}, arg17::Ptr{Cdouble}, arg18::Ptr{Clong}, arg19::Ptr{klu_l_common})::Clong
end

function klu_malloc(n, size, Common)
    @ccall libklu.klu_malloc(n::Csize_t, size::Csize_t, Common::Ptr{klu_common})::Ptr{Cvoid}
end

function klu_free(p, n, size, Common)
    @ccall libklu.klu_free(p::Ptr{Cvoid}, n::Csize_t, size::Csize_t, Common::Ptr{klu_common})::Ptr{Cvoid}
end

function klu_realloc(nnew, nold, size, p, Common)
    @ccall libklu.klu_realloc(nnew::Csize_t, nold::Csize_t, size::Csize_t, p::Ptr{Cvoid}, Common::Ptr{klu_common})::Ptr{Cvoid}
end

function klu_l_malloc(arg1, arg2, arg3)
    @ccall libklu.klu_l_malloc(arg1::Csize_t, arg2::Csize_t, arg3::Ptr{klu_l_common})::Ptr{Cvoid}
end

function klu_l_free(arg1, arg2, arg3, arg4)
    @ccall libklu.klu_l_free(arg1::Ptr{Cvoid}, arg2::Csize_t, arg3::Csize_t, arg4::Ptr{klu_l_common})::Ptr{Cvoid}
end

function klu_l_realloc(arg1, arg2, arg3, arg4, arg5)
    @ccall libklu.klu_l_realloc(arg1::Csize_t, arg2::Csize_t, arg3::Csize_t, arg4::Ptr{Cvoid}, arg5::Ptr{klu_l_common})::Ptr{Cvoid}
end

const KLU_OK = 0

const KLU_SINGULAR = 1

const KLU_OUT_OF_MEMORY = -2

const KLU_INVALID = -3

const KLU_TOO_LARGE = -4

const KLU_DATE = "Mar 12, 2018"

KLU_VERSION_CODE(main, sub) = main * 1000 + sub

const KLU_MAIN_VERSION = 1

const KLU_SUB_VERSION = 3

const KLU_SUBSUB_VERSION = 9

const KLU_VERSION = KLU_VERSION_CODE(KLU_MAIN_VERSION, KLU_SUB_VERSION)

const SuiteSparse_long = Clong

const LONG_MAX = typemax(Clong)

const SuiteSparse_long_max = LONG_MAX

const SuiteSparse_long_idd = "ld"

const SUITESPARSE_DATE = "May 17, 2021"

SUITESPARSE_VER_CODE(main, sub) = main * 1000 + sub

const SUITESPARSE_MAIN_VERSION = 5

const SUITESPARSE_SUB_VERSION = 10

const SUITESPARSE_SUBSUB_VERSION = 1

const SUITESPARSE_VERSION = SUITESPARSE_VER_CODE(SUITESPARSE_MAIN_VERSION, SUITESPARSE_SUB_VERSION)

end


text:
# LEARNING:
# - Julia has "." operators  that perform an operation elementwise on everything in an array
#   https://docs.julialang.org/en/v1/manual/mathematical-operations/#man-dot-operators
#   .. leads to concise syntax like: https://github.com/leesharma/game-of-life-julia/blob/main/rules/standard.jl#L5

using Test;

function run(fname, version=1)
    f = open(fname, "r")
    lines = readlines(f)
    # Each position is either floor (.), an empty seat (L), or an occupied seat (#).

    prev = lines
    while true
        new_lines = take_step(prev, version)
        if prev == new_lines
            break
        end
        prev = new_lines
    end

    occupied_chairs = 0
    for row in prev
        for c in row
            if c == '#'
                occupied_chairs += 1
            end
        end
    end

    result = occupied_chairs
    println("[part ", version, "] Ran file ", fname, " and got result = ", result)
    return result
end

function take_step(lines, version=1)::Array{String}
    # If a seat is empty (L) and there are no occupied seats adjacent to it, the seat becomes occupied.
    # If a seat is occupied (#) and four or more seats adjacent to it are also occupied, the seat becomes empty.
    # Otherwise, the seat's state does not change.
    f = get_new_state
    if version == 2
        f = get_new_state_v2
    end

    out = []
    for row in 1:length(lines)
        new_row = ""
        for col in 1:length(lines[1])
            n = f(lines, row, col)
            new_row *= n
        end
        push!(out, new_row)
    end
    return out
end

row_offsets = [-1, 0, 1]
col_offsets = [-1, 0, 1]
function get_new_state(lines, row, col)::Char
    if lines[row][col] == '.'
        return '.'
    end

    # count occupied neighbors
    occupied = 0
    total_rows = length(lines)
    total_cols = length(lines[1])
    for r in row_offsets
        for c in col_offsets
            # if in bounds
            if !(r == 0 && c == 0) &&
                (row + r >= 1 && row + r <= total_rows) &&
                (col + c >= 1 && col + c <= total_cols)
                if lines[row + r][col + c] == '#'
                    occupied += 1
                end
            end
        end
    end

    seat = lines[row][col]
    if seat == 'L' && occupied == 0
        return '#'
    elseif occupied >= 4
        return 'L'
    else
        return lines[row][col]
    end
end

function get_new_state_v2(lines, row, col)::Char
    if lines[row][col] == '.'
        return '.'
    end

    # count occupied neighbors
    occupied = 0
    total_rows = length(lines)
    total_cols = length(lines[1])
    for r in row_offsets
        for c in col_offsets
            if (r == 0 && c == 0)
                # don't include current seat
                continue
            end

            # while in bounds..
            step_size = 1
            new_r = row + (step_size * r)
            new_c = col + (step_size * c)
            while (new_r >= 1 && new_r <= total_rows) &&
                    (new_c >= 1 && new_c <= total_cols)
                val = lines[new_r][new_c]
                if val == '.'
                    step_size += 1
                    new_r = row + (step_size * r)
                    new_c = col + (step_size * c)
                elseif val == '#'
                    occupied += 1
                    break
                else
                    # val is 'L'
                    break
                end
            end
        end
    end

    seat = lines[row][col]
    if seat == 'L' && occupied == 0
        return '#'
    elseif occupied >= 5
        return 'L'
    else
        return seat
    end
end

ex = readlines(open("11ex.txt", "r"))
ex_1 = readlines(open("11ex.1.txt", "r"))
ex_2 = readlines(open("11ex.2.txt", "r"))
ex_5 = readlines(open("11ex.5.txt", "r"))

@test take_step(ex) == ex_1
# LEARNING: Julia has function composition operator (\circ)
# - https://docs.julialang.org/en/v1/manual/functions/#Function-composition-and-piping
# - https://github.com/leesharma/game-of-life-julia/blob/main/manifolds/grid_torus.jl#L18-L21
take_2_steps = take_step ∘ take_step
@test take_2_steps(ex) == ex_2
@test take_step(take_step(take_step(take_step(take_step(ex))))) == ex_5
# .. stays stable?
@test take_step(take_step(take_step(take_step(take_step(take_step(ex)))))) == ex_5

@test run("11ex.txt") == 37
run("11.txt")

println("VERSION 2")

ex_2_ver2 = readlines(open("11ex.2b.txt", "r"))
@test take_step(ex, 2) == ex_1
@test take_step(take_step(ex, 2), 2) == ex_2_ver2

@test run("11ex.txt", 2) == 26
run("11.txt", 2)


text:
using RobustShortestPath

data = [
 1   4  79   31  66  28;
 1   2  59   97  41  93;
 2   4  31   21  50  40;
 2   3  90   52  95  38;
 2   5   9   23  95  59;
 2   6  32   57  73   7;
 3   9  89  100  38  21;
 3   8  66   13   4  72;
 3   6  68   95  58  58;
 3   7  47   12  56  20;
 4   3  14   19  36  84;
 4   9  95   65  88  42;
 4   8  88   13  62  54;
 5   3  44    8  62  53;
 5   6  83   66  30  19;
 6   7  33    3   7   8;
 6   8  37   99  29  46;
 7  11  79   54  23   3;
 7  12  10   37  35  43;
 8   7  95   71  85  56;
 8  10   0   95  16  64;
 8  12  30   38  16   3;
 9  10   5   69  51  71;
 9  11  44   60  60  17;
10  13  79   78  16  59;
10  14  91   59  64  61;
11  14  53   38  84  77;
11  15  80   85  78   6;
11  13  56   23  26  85;
12  15  75   80  31  38;
12  14   1  100  18  40;
13  14  48   28  45  33;
14  15  25   71  33  56;
]

start_node = data[:,1] #first column of data
end_node = data[:,2] #second column of data
p = data[:,3] #third
q = data[:,4] #fourth
c = data[:,5] #fifth
d = data[:,6] #sixth

# Link length vectors for single-uncertain-coefficient robust shortest paths
cc = p.*c
dd = (p+q).*(c+d) - p.*c

# Setting origin and destination nodes
origin = 1
destination = 15
println("Origin=$origin, Destination=$destination")

println("----------------------------------------------------")
println("Single Coefficient Case")
# For each Gamma from 0 to 6, obtain the robust shortest path
known_cost_solution_one = Dict(0 => 6060, 1=> 15024, 2=> 20864, 3=> 26604, 4=> 31293, 5=> 32291)
for Gamma=0:5
	robust_path, robust_x, worst_case_cost = get_robust_path(start_node, end_node, cc, dd, Gamma, origin, destination)
	println("Gamma=$Gamma: Robust Path is $(robust_path') and the worst-case cost is $worst_case_cost.")

    @assert worst_case_cost == known_cost_solution_one[Gamma]
end

println("----------------------------------------------------")
println("Two Coefficient Case")

for Gamma_u=1:5
    for Gamma_v=1:5
        robust_path, robust_x, worst_case_cost = get_robust_path_two(start_node, end_node, p, q, c, d, Gamma_u, Gamma_v, origin, destination)
        println("(Gamma_u,Gamma_v)=($Gamma_u,$Gamma_v): Robust Path is $(robust_path') and the worst-case cost is $worst_case_cost.")
    end
end

test_Gamma_u = 5
test_Gamma_v = 5
robust_path, robust_x, worst_case_cost = get_robust_path_two(start_node, end_node, p, q, c, d, test_Gamma_u, test_Gamma_v, origin, destination)

println(worst_case_cost)
@assert worst_case_cost==32291.0


text:
function CreateRHD2000Registers(sampleRate)

    r=register(zeros(Int32,40)...,zeros(Int32,1))

    defineSampleRate(sampleRate,r)

    dsp_b = 300.0
    band_a_l = 300.0
    band_a_u = 5000.0

    #Register 0 variables
    r.adcReferenceBw = 3
    setFastSettle(false,r)
    r.ampVrefEnable=1
    r.adcComparatorBias=3
    r.adcComparatorSelect=2

    #Register 1 variables
    r.vddSenseEnable=1

    #Register 3 variables
    r.tempS1=0
    r.tempS2=0
    r.tempEn=0
    setDigOutHiZ(r)

    #Register 4 variables
    r.weakMiso=1
    r.twosComp=0
    r.absMode=0
    enableDsp(true,r)
    setDspCutoffFreq(dsp_b,r)

    #Register 5 variables
    r.zcheckDacPower=0
    r.zcheckLoad=0
    setZcheckScale("ZcheckCs100fF", r)
    r.zcheckConnAll=0
    setZcheckPolarity("ZcheckPositiveInput",r)
    enableZcheck(false,r)

    #Register 7 variable
    setZcheckChannel(0,r)

    #Register 8-13 variables
    r.offChipRH1=0
    r.offChipRH2=0
    r.offChipRL=0
    r.adcAux1En=1
    r.adcAux2En=1
    r.adcAux3En=1

    #these set RHdac and RLdac variables
    setUpperBandwidth(band_a_u,r)
    setLowerBandwidth(band_a_l,r)

    powerUpAllAmps(r)

    r
end

function defineSampleRate(newSampleRate,r)

    r.sampleRate=newSampleRate

    r.muxLoad=0

    if r.sampleRate<3334.0
        r.muxBias=40
        r.adcBufferBias=32
    elseif r.sampleRate<4001.0
        r.muxBias=40
        r.adcBufferBias=16
    elseif r.sampleRate<5001.0
        r.muxBias=40
        r.adcBufferBias=8
    elseif r.sampleRate<6251.0
        r.muxBias=32
        r.adcBufferBias=8
    elseif r.sampleRate<8001.0
        r.muxBias=26
        r.adcBufferBias=8
    elseif r.sampleRate<10001.0
        r.muxBias=18
        r.adcBufferBias=4
    elseif r.sampleRate<12501.0
        r.muxBias=16
        r.adcBufferBias=3
    elseif r.sampleRate<15001.0
        r.muxBias=7
        r.adcBufferBias=3
    else
        r.muxBias=4
        r.adcBufferBias=2
    end

    nothing
end

setFastSettle(enabled, r)=r.ampFastSettle=(enabled ? 1 : 0)

setDigOutLow(r)=(r.digOut=0;r.digOutHiZ=0)

setDigOutHigh(r)=(r.digOut=1;r.digOutHiZ=0)

setDigOutHiZ(r)=(r.digOut=0; r.digOutHiZ=1)

function setDspCutoffFreq(newDspCutoffFreq, r)

    logNewDspCutoffFreq = log10(newDspCutoffFreq)

    fCutoff=zeros(Float64,16)
    logFCutoff=zeros(Float64,16)

    for n=2:16
        x= 2.0 ^ (n-1)
        fCutoff[n]= r.sampleRate * log(x / (x - 1.0)) / (2*pi)
        logFCutoff[n] = log10(fCutoff[n])
    end

    if newDspCutoffFreq > fCutoff[2]
        r.dspCutoffFreq = 1
    else
        temp=zeros(Float64,15)
        for i=1:15
            temp[i]=abs(logNewDspCutoffFreq - logFCutoff[i+1])
        end
        r.dspCutoffFreq=findmin(temp)[2]
    end

    nothing
end

enableAux1(enabled, r)=(r.adcAux1En = (enabled ? 1 : 0))

enableAux2(enabled, r)=(r.adcAux2En = (enabled ? 1 : 0))

enableAux3(enabled, r)=(r.adcAux3En = (enabled ? 1 : 0))

enableDsp(enabled, r)=(r.dspEn = (enabled ? 1 : 0))

enableZcheck(enabled, r)=(r.zcheckEn = (enabled ? 1 : 0))

setZcheckDacPower(enabled, r)=(r.zcheckDacPower = (enabled ? 1 : 0))

function setZcheckScale(scale, r)

    if scale=="ZcheckCs100fF"
        r.zcheckScale = 0x00
    elseif scale=="ZcheckCs1pF"
        r.zcheckScale = 0x01
    elseif scale=="ZcheckCs10pF"
        r.zcheckScale = 0x03
    end
    nothing
end

function setZcheckPolarity(polarity, r)

    if polarity=="ZcheckPositiveInput"
        r.zcheckSelPol=0
    elseif polarity=="ZcheckNegativeInput"
        r.zcheckSelPol=1
    end

    nothing
end

setZcheckChannel(channel, r)=(r.zcheckSelect=channel)

powerUpAllAmps(r)=(r.aPwr=ones(Int32, 64))

powerDownAllAmps(r)=(r.aPwr=zeros(Int32,64))

const RLBase = 3500.0
const RLDac1Unit = 175.0
const RLDac2Unit = 12700.0
const RLDac3Unit = 3000000.0
const RLDac1Steps = 127
const RLDac2Steps = 63

function setLowerBandwidth(lowerBandwidth, r)

    if lowerBandwidth > 1500.0
        lowerBandwidth = 1500.0
    end

    rLTarget = rLFromLowerBandwidth(lowerBandwidth)

    r.rLDac1 = 0
    r.rLDac2 = 0
    r.rLDac3 = 0
    rLActual = RLBase

    if lowerBandwidth < 0.15
        rLActual += RLDac3Unit
        r.rLDac3+=1
    end

    for i=1:RLDac2Steps
        if rLActual < (rLTarget - (RLDac2Unit - RLDac1Unit/2))
            rLActual += RLDac2Unit
            r.rLDac2+=1
        end
    end

    for i=1:RLDac1Steps
        if rLActual < (rLTarget - (RLDac1Unit /2))
            rLActual += RLDac1Unit
            r.rLDac1+=1
        end
    end

    actualLowerBandwidth=lowerBandwidthFromRL(rLActual)

    nothing
end

function setUpperBandwidth(upperBandwidth, r)

    RH1Base = 2200.0
    RH1Dac1Unit = 600.0
    RH1Dac2Unit = 29400.0
    RH1Dac1Steps = 63
    RH1Dac2Steps = 31

    RH2Base = 8700.0
    RH2Dac1Unit = 763.0
    RH2Dac2Unit = 38400.0
    RH2Dac1Steps = 63
    RH2Dac2Steps = 31

    if upperBandwidth > 30000.0
        upperBandwidth = 30000.0
    end

    rH1Target = rH1FromUpperBandwidth(upperBandwidth)

    r.rH1Dac1 = 0
    r.rH1Dac2 = 0
    rH1Actual = RH1Base

    for i=1:RH1Dac2Steps

        if rH1Actual < (rH1Target - (RH1Dac2Unit - RH1Dac1Unit/2))
            rH1Actual += RH1Dac2Unit
            r.rH1Dac2+=1
        end

    end

    for i=1:RH1Dac1Steps

        if rH1Actual < (rH1Target - (RH1Dac1Unit/2))
            rH1Actual += RH1Dac1Unit
            r.rH1Dac1+=1
        end
    end

    rH2Target= rH2FromUpperBandwidth(upperBandwidth)

    r.rH2Dac1=0
    r.rH2Dac2=0
    rH2Actual = RH2Base

    for i=1:RH2Dac2Steps

        if rH2Actual < (rH2Target - (RH2Dac2Unit - RH2Dac1Unit/2))
            rH2Actual += RH2Dac2Unit
            r.rH2Dac2+=1
        end

    end

    for i=1:RH2Dac1Steps

        if rH2Actual < (rH2Target - (RH2Dac1Unit/2))
            rH2Actual += RH2Dac1Unit
            r.rH2Dac1+=1
        end
    end

    actualUpperBandwidth1 = upperBandwidthFromRH1(rH1Actual)
    actualUpperBandwidth2 = upperBandwidthFromRH2(rH2Actual)

    actualUpperBandwidth = sqrt(actualUpperBandwidth1 * actualUpperBandwidth2)

    nothing
end

function rLFromLowerBandwidth(lowerBandwidth)

    log10f = log10(lowerBandwidth)

    if (lowerBandwidth < 4.0)
       return 1.0061 * 10.0 ^ (4.9391 - 1.2088 * log10f + 0.5698 * log10f * log10f + 0.1442 * log10f * log10f * log10f)

    else
        return 1.0061 * 10.0 ^ (4.7351 - 0.5916 * log10f + 0.08482 * log10f * log10f)
    end

end

function rH1FromUpperBandwidth(upperBandwidth)

    log10f = log10(upperBandwidth)

    return 0.9730 * 10.0 ^ (8.0968 - 1.1892 * log10f + 0.04767 * log10f * log10f)

end

function rH2FromUpperBandwidth(upperBandwidth)

    log10f = log10(upperBandwidth)

    return 1.0191 * 10.0 ^ (8.1009 - 1.0821 * log10f + 0.03383 * log10f * log10f)

end

function upperBandwidthFromRH1(rH1)

    a = 0.04767
    b = -1.1892
    c = 8.0968 - log10(rH1/0.9730)

    10.0 ^ ((-b - sqrt(b * b - 4 * a * c))/(2 * a))
end

function upperBandwidthFromRH2(rH2)

    a = 0.03383
    b = -1.0821
    c = 8.1009 - log10(rH2/1.0191)

    10.0 ^ ((-b - sqrt(b * b - 4 * a * c))/(2 * a))
end

function lowerBandwidthFromRL(rL)

    if (rL < 5100.0)
        rL = 5100.0
    end

    if (rL < 30000.0)
        a = 0.08482
        b = -0.5916
        c = 4.7351 - log10(rL/1.0061)
    else
        a = 0.3303
        b = -1.2100
        c = 4.9873 - log10(rL/1.0061)
    end

    10.0 ^ ((-b - sqrt(b * b - 4 * a * c))/(2 * a))
end

function createRhd2000Command(commandType, arg1=0, arg2=0)

    if commandType=="Rhd2000CommandCalibrate"
        return 0x5500
    elseif commandType=="Rhd2000CommandCalClear"
        return 0x6a00
    elseif commandType=="Rhd2000CommandConvert"
        return 0x0000 + (arg1 << 8)
    elseif commandType=="Rhd2000CommandRegRead"
        return 0xc000 + (arg1 << 8)
    elseif commandType=="Rhd2000CommandRegWrite"
        return 0x8000 + (arg1 << 8) + arg2
    end
end

function getRegisterValue(reg, r)

    zcheckDac=128

    aPwr=r.aPwr

    if reg==0
        regout = (r.adcReferenceBw << 6) + (r.ampFastSettle << 5) + (r.ampVrefEnable << 4) +
        (r.adcComparatorBias << 2) + r.adcComparatorSelect
    elseif reg==1
        regout = (r.vddSenseEnable << 6) + r.adcBufferBias
    elseif reg==2
        regout = r.muxBias
    elseif reg==3
        regout = (r.muxLoad << 5) + (r.tempS2 << 4) + (r.tempS1 << 3) + (r.tempEn << 2) +
        (r.digOutHiZ << 1) + r.digOut
    elseif reg==4
        regout = (r.weakMiso << 7) + (r.twosComp << 6) + (r.absMode << 5) + (r.dspEn << 4) +
        r.dspCutoffFreq
    elseif reg==5
        regout = (r.zcheckDacPower << 6) + (r.zcheckLoad << 5) + (r.zcheckScale << 3) +
        (r.zcheckConnAll << 2) + (r.zcheckSelPol << 1) + r.zcheckEn
    elseif reg==6
        regout = zcheckDac
    elseif reg==7
        regout = r.zcheckSelect
    elseif reg==8
        regout = (r.offChipRH1 << 7) + r.rH1Dac1
    elseif reg==9
        regout = (r.adcAux1En << 7) + r.rH1Dac2
    elseif reg==10
        regout = (r.offChipRH2 << 7) + r.rH2Dac1
    elseif reg==11
        regout = (r.adcAux2En << 7) + r.rH2Dac2
    elseif reg==12
        regout = (r.offChipRL << 7) + r.rLDac1
    elseif reg==13
        regout = (r.adcAux3En << 7) + (r.rLDac3 << 6) + r.rLDac2
    elseif reg==14
        regout = (aPwr[8] << 7) + (aPwr[7] << 6) + (aPwr[6] << 5) + (aPwr[5] << 4) +
        (aPwr[4] << 3) + (aPwr[3] << 2) + (aPwr[2] << 1) + aPwr[1]
    elseif reg==15
        regout = (aPwr[16] << 7) + (aPwr[15] << 6) + (aPwr[14] << 5) + (aPwr[13] << 4) +
        (aPwr[12] << 3) + (aPwr[11] << 2) + (aPwr[10] << 1) + aPwr[1]
    elseif reg==16
        regout = (aPwr[24] << 7) + (aPwr[23] << 6) + (aPwr[22] << 5) + (aPwr[21] << 4) +
        (aPwr[20] << 3) + (aPwr[19] << 2) + (aPwr[18] << 1) + aPwr[17]
    elseif reg==17
        regout = (aPwr[32] << 7) + (aPwr[31] << 6) + (aPwr[30] << 5) + (aPwr[29] << 4) +
        (aPwr[28] << 3) + (aPwr[27] << 2) + (aPwr[26] << 1) + aPwr[25]
    elseif reg==18
        regout = (aPwr[40] << 7) + (aPwr[39] << 6) + (aPwr[38] << 5) + (aPwr[37] << 4) +
        (aPwr[36] << 3) + (aPwr[35] << 2) + (aPwr[34] << 1) + aPwr[33]
    elseif reg==19
        regout = (aPwr[48] << 7) + (aPwr[47] << 6) + (aPwr[46] << 5) + (aPwr[45] << 4) +
        (aPwr[44] << 3) + (aPwr[43] << 2) + (aPwr[42] << 1) + aPwr[41]
    elseif reg==20
        regout = (aPwr[56] << 7) + (aPwr[55] << 6) + (aPwr[54] << 5) + (aPwr[53] << 4) +
        (aPwr[52] << 3) + (aPwr[51] << 2) + (aPwr[50] << 1) + aPwr[49]
    elseif reg==21
        regout = (aPwr[64] << 7) + (aPwr[63] << 6) + (aPwr[62] << 5) + (aPwr[61] << 4) +
        (aPwr[60] << 3) + (aPwr[59] << 2) + (aPwr[58] << 1) + aPwr[57]
    end

    return regout
end

#Create a list of 60 commands to program most RAM registers on a RHD2000 chip, read those values back to confirm programming, read ROM registers, and (if calibrate == true) run ADC calibration.

function createCommandListRegisterConfig(commandList, calibrate, r)

    #Start with a few dummy commands in case chip is still powering up
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 63))
    splice!(commandList,1) #remove first dummy command
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 63))

    #Program RAM registers
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 0, getRegisterValue(0, r)))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 1, getRegisterValue(1, r)))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 2, getRegisterValue(2, r)))

    #Don't program Register 3 (MUX Load, Temperature Sensor, and Auxiliary Digital Output)

    #control temperature sensor in another command stream
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 4, getRegisterValue(4, r)))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 5, getRegisterValue(5, r)))

    #Don't program Register 6 (Impedence Check DAC) here; create DAC waveform in another command stream
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 7, getRegisterValue(7, r)))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 8, getRegisterValue(8, r)))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 9, getRegisterValue(9, r)))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 10, getRegisterValue(10, r)))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 11, getRegisterValue(11, r)))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 12, getRegisterValue(12, r)))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 13, getRegisterValue(13, r)))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 14, getRegisterValue(14, r)))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 15, getRegisterValue(15, r)))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 16, getRegisterValue(16, r)))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 17, getRegisterValue(17, r)))

    #Read ROM registers
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 63))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 62))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 61))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 60))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 59))

    #Read chip name from ROM
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 48))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 49))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 50))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 51))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 52))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 53))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 54))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 55))

    #Read Intan name from ROM
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 40))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 41))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 42))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 43))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 44))

    #Read back RAM registers to confirm programming
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 0))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 1))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 2))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 3))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 4))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 5))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 6))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 7))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 8))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 9))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 10))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 11))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 12))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 13))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 14))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 15))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 16))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead", 17))

    #Optionally, run ADC calibration (should only be run once after board is plugged in
    if calibrate
        push!(commandList, createRhd2000Command("Rhd2000CommandCalibrate"))
    else
        push!(commandList, createRhd2000Command("Rhd2000CommandRegRead",63))
    end

    #Program amplifier 31-63 power up/down registers in case RHD2164 is connected
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 18, getRegisterValue(18,r)))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 19, getRegisterValue(19,r)))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 20, getRegisterValue(20,r)))
    push!(commandList, createRhd2000Command("Rhd2000CommandRegWrite", 21, getRegisterValue(21,r)))

    #End with a dummy command
    push!(commandList, createRhd2000Command("Rhd2000CommandRegRead",63))

return commandList

end


text:
function train!(model,X,y,opt=Adam())
    # ps = Flux.params(model)
    loss(X,y) = ELBO(model,X,y)
    # Flux.train!(loss,ps,(X,y),opt)
end