text:
"""
    X = sylvckr(A,B,C)

Solve the continuous Sylvester matrix equation

                AX + XB = C

using the Kronecker product expansion of equations. `A` and `B` are
square matrices, and `A` and `-B` must not have common eigenvalues.
This function is not recommended for large order matrices.
"""
function sylvckr(A, B, C)
    m, n = size(C)
    [m; n] == LinearAlgebra.checksquare(A, B) ||
             throw(DimensionMismatch("A, B and Q have incompatible dimensions"))
    reshape((kron(Array{eltype(A),2}(I, n, n), A) +
            kron(transpose(B), Array{eltype(B),2}(I, m, m))) \ (C[:]),m,n)
end

"""
    X = sylvdkr(A,B,C)

Solve the discrete Sylvester matrix equation

                AXB + X = C

using the Kronecker product expansion of equations. `A` and `B` are
square matrices, and `A` and `-B` must not have common reciprocal eigenvalues.
This function is not recommended for large order matrices.
"""
function sylvdkr(A, B, C)
    m, n = size(C)
    [m; n] == LinearAlgebra.checksquare(A, B) ||
            throw(DimensionMismatch("A, B and C have incompatible dimensions"))
    reshape((kron(transpose(B), A) + I) \ (C[:]), m, n)
end

"""
    X = gsylvkr(A,B,C,D,E)

Solve the generalized Sylvester matrix equation

                AXB + CXD = E

using the Kronecker product expansion of equations. `A`, `B`, `C` and `D` are
square matrices. The pencils `A-λC` and `D+λB` must be regular and
must not have common eigenvalues.
This function is not recommended for large order matrices.
"""
function gsylvkr(A, B, C, D, E)
   m, n = size(E)
   [m; n; m; n] == LinearAlgebra.checksquare(A, B, C, D) ||
                  throw(DimensionMismatch("A, B, C, D and E have incompatible dimensions"))
   reshape((kron(transpose(B), A) + kron(transpose(D), C)) \ (E[:]), m, n)
end

"""
    sylvsyskr(A,B,C,D,E,F) -> (X,Y)

Solve the Sylvester system of matrix equations

                AX + YB = C
                DX + YE = F

using the Kronecker product expansion of equations. `(A,D)`, `(B,E)` are
pairs of square matrices of the same size.
The pencils `A-λD` and `-B+λE` must be regular and must not have common eigenvalues.
This function is not recommended for large order matrices.
"""
function sylvsyskr(A, B, C, D, E, F)
   m, n = size(C)
   (m == size(F, 1) && n == size(F, 2)) ||
         throw(DimensionMismatch("C and F must have the same dimensions"))
   [m; n; m; n] == LinearAlgebra.checksquare(A, B, D, E) ||
                  throw(DimensionMismatch("A, B, C, D, E and F have incompatible dimensions"))
   z = [ kron(Array{eltype(A),2}(I, n, n), A) kron(transpose(B), Array{eltype(B),2}(I, m, m)) ;
         kron(Array{eltype(D),2}(I, n, n), D) kron(transpose(E), Array{eltype(E),2}(I, m, m))] \ [C[:];F[:]]
   (reshape(z[1:m * n], m, n), reshape(z[m * n + 1:end], m, n))
end

"""
    dsylvsyskr(A,B,C,D,E,F) -> (X,Y)

Solve the dual Sylvester system of matrix equations

       AX + DY = C
       XB + YE = F

using the Kronecker product expansion of equations. `(A,D)`, `(B,E)` are
pairs of square matrices of the same size.
The pencils `A-λD` and `-B+λE` must be regular and must not have common eigenvalues.
This function is not recommended for large order matrices.
"""
function dsylvsyskr(A, B, C, D, E, F)
   m, n = size(C)
   (m == size(F, 1) && n == size(F, 2)) ||
         throw(DimensionMismatch("C and F must have the same dimensions"))
   [m; n; m; n] == LinearAlgebra.checksquare(A, B, D, E) ||
                  throw(DimensionMismatch("A, B, C, D, E and F have incompatible dimensions"))
   z = [ kron(Array{eltype(A),2}(I, n, n), A) kron(Array{eltype(D),2}(I, n, n), D);
         kron(transpose(B), Array{eltype(B),2}(I, m, m)) kron(transpose(E), Array{eltype(E),2}(I, m, m))] \ [C[:];F[:]]
   (reshape(z[1:m * n], m, n), reshape(z[m * n + 1:end], m, n))
end

max_stars_repo_name: andreasvarga/MatrixEquations

max_issues_repo_name: andreasvarga/MatrixEquations

max_forks_repo_name: andreasvarga/MatrixEquations


text:
#------------------------------------------------------------------------------

struct ParseException <: Exception
    msg::String
end

Base.showerror(io::IO, e::ParseException) = print(io, e.msg)

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

# returns a model parsed from the given `filename`
# (with matrices in dense format)
# each returned matrix is represented with a dense array
function load_dense_set(filename::String)
    ans = []

    filepath = "$(normpath(joinpath(@__FILE__, "..", "..")))/data/matrices/$(filename)"

    open(filepath, "r") do input_stream
        for arr in JSON.Parser.parse(input_stream)
            push!(ans, arr)
        end
    end

    return ans
end

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

# returns a model parsed from the given `filename`
# (with matrices in coordinate format)
# each returned matrix is represented with an array of
# items (i, j, x) where typeof(x) is Rational
function load_COO_set(filename::String)
    model = []

    filepath = "$(normpath(joinpath(@__FILE__, "..", "..")))/data/matrices/$(filename)"

    str_to_nom = (x -> parse(BigInt, first(split(x, '/'))))
    str_to_den = (x -> (if occursin('/', x) parse(BigInt, last(split(x, '/'))) else 1 end))

    open(filepath, "r") do input_stream
        for arr in JSON.Parser.parse(input_stream)
            normal_repr = [
                arr[1], arr[2],
                map(x -> (x[1], x[2], str_to_nom(x[3]) // str_to_den(x[3])), arr[3])
            ]

            push!(model, normal_repr)
        end
    end

    return model
end

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

# returns an array of all models from data/matrices that satisfy
# the given criterias:
#   `from_dim`  - lower bound for the model's matrices one-side dimension
#   `to_dim`    - upper bound -//-
#   `from_size` - lower bound for the model's total number of matrices
#   `to_size`   - upper bound -//-
#
# each returned model is represented with an array:
#   [filename, dim, size, matrices-data]
# each returned matrix is represented with an array of
# items (i, j, x) of nonzeroes where typeof(x) is Rational
function load_COO_if(;from_dim, to_dim, from_size=0, to_size=Inf)
    models = []

    for filename in readdir("data/matrices/")
        model_data = load_COO_set(filename)
        model_dim = first(first(model_data))
        model_size = length(model_data)

        if from_dim <= model_dim <= to_dim &&
                from_size <= model_size <= to_size
            push!(models, [filename, model_dim, model_size, model_data])
        end
    end

    @info "loaded models: $(length(models))"

    return models
end

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

# returns a model parsed from the given `group\name` file
# (with matrices in MTX format)
# each returned matrix is represented with an array of
#   [filename, dim, size, data]
# where data is stored as items (i, j, x) where typeof(x) is Rational
#
# In case of any parse error occurred, e.g inexact rationals encountered,
# returs an empty array []
function load_MTX(group, name)
    model = []

    # seems to be a default separator
    sep = r"%-+\r?\n"

    filepath = replace(
        "$(normpath(joinpath(@__FILE__, "..", "..")))data/sscollection/$group/$name",
        "\\" => "/"
    )

    open(filepath, "r") do input_stream
        # some models have >3 sections, taking the last
        data = last(split(read(input_stream, String), sep))

        # should be ok
        data = split(data, r"\r?\n")
        meta = map(x -> parse(Int, x), split(first(data), ' '))
        # strange matrices, discard them
        # NEED TO BE HANDLED IN ANOTHER WAY
        if length(meta) != 3
            return []
        end

        m, n, nnz = meta
        # we want only square ones
        if m != n
            return []
        end

        try
            model = [
                group * '/' * name, # filename
                m,                  # dim
                1,                  # size
                collect(            # data
                    map(
                        # `Any` disallows auto-promoting Int to Rational
                        X -> Any[
                            parse(Int, X[1]),
                            parse(Int, X[2]),
                            length(X) > 2 ? Rational(parse(Float64, X[3])) : Rational(1)
                        ] ,
                        filter(!isempty, map(split, view(data, 1:length(data))))
                    )
                )
            ]
        catch e
            if isa(e, InexactError)
                # occurrs if floats are impossibly small/big,
                # discard the model then
                return []
            else
                rethrow(e)
            end
        end
    end

    return model
end

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

# returns an array of all models from data/sscollection that satisfy
# the given criterias:
#   `from_dim`  - lower bound for the model's matrices one-side dimension
#   `to_dim`    - upper bound -//-
#
# each returned model is represented with an array:
#   [filename, dim, size, matrices-data]
# each returned matrix is represented with an array of
# items (i, j, x) of nonzeroes where typeof(x) is Rational
function load_MTX_if(;from_dim=0, to_dim=Inf)
    models = []

    extension = ".mtx"
    dir = "src/data/sscollection/"

    for group in readdir(dir)
        for name in readdir("$dir/$group/")
            if endswith(name, extension)
                model = load_MTX(group, name)
                if ! isempty(model)
                    fn, dim, sz, data = model
                    if from_dim <= model[2] <= to_dim
                        push!(models, [ fn, dim, sz, [ dim, dim, data ] ] )
                    end
                end
            end
        end
    end

    @info "loaded models: $(length(models))"

    return models
end

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

"""
    load_ODEs_recursive_if(path; from_size=-Inf, to_size=Inf)

Finds a basis of the algebra generated with the given `vectors`
using the provived default method `used_algorithm` while reducing input coefficients
modulo `initialprime`. By default, `find_basis_1_β` stands for the algorithm and
`2147483647` for the initial reduction modulo.

Each element in `vector` to be a subtype of `AbstractSparseObject`
"""
function load_ODEs_recursive_if(
                        path;
                        from_size=-Inf, to_size=Inf,
                        exclude=[
                            "curated_erode",
                            "__MACOSX",
                            "non_curated_sbml"
                        ])
    ext = ".ode"

    models = []
    dirs_to_process = [path]

    while !isempty(dirs_to_process)
        path = pop!(dirs_to_process)
        path = normalizepath(path)
        for (pathname, dirs, files) in walkdir(path, onerror=print)
            @info "Walk: " pathname dirs files
            for fname in files
                if endswith(fname, ext)
                    @info "Load: " "$pathname/$fname"
                    try
                        model = load_ODEs("$pathname/$fname")
                        if from_size <= length(model) <= to_size
                            push!(models, [fname, model])
                        end
                    catch ex
                        !isa(ex, ParseException) && rethrow(ex)
                        @warn ex.msg
                    end
                end
            end
            dirs = filter(x -> !(x in exclude), dirs)
            append!(dirs_to_process, dirs)
        end
    end

    @info "loaded models: $(length(models))"

    models
end

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

function load_ODEs_if(;from_size=-Inf, to_size=Inf)
    testspath = "src/data/ODEs/"
    ext = ".ode"

    models = []

    for (path, dirs, files) in walkdir(testspath)
        for fname in files
            if endswith(fname, ext)
                try
                    model = load_ODEs("$path/$fname")
                    if from_size <= length(model) <= to_size
                        push!(models, [fname, model])
                    end
                catch ex
                    #!isa(ex, ParseException) && rethrow(ex)
                    @warn ex.msg
                end
            end
        end
    end

    @info "loaded models: $(length(models))"

    models
end

# filepath - absolute model file path
#
function load_ODEs(filepath)
    # *.ode --> polynomial form of ode system

    lines = []
    open(filepath, "r") do inputs
        lines = map(strip, readlines(inputs))
    end

    # @info "Loading $filepath"

    if isnothing(findlast(startswith("begin"), lines))
        throw(ParseException("bad file encountered at $filepath, skipping it"))
    end

    lines = filter(!isempty, map(strip, lines[
        (findfirst(startswith("begin"), lines) + 1
        :
        findlast(startswith("end"), lines) - 1)
    ]))

    if isnothing(findlast(startswith("begin"), lines))
        throw(ParseException("bad file encountered at $filepath, skipping it"))
    end

    # by default we consider the last begin/end block to yield reactions.
    # If this conjecture is false, we do not work
    reactions = filter(!isempty, map(strip, lines[
        (findlast(startswith("begin"), lines) + 1
        :
        findlast(startswith("end"), lines) - 1)
    ]))

    # removing comments at the end of reactions
    reactions = map(s -> replace(s, r"\[.*\]" => ""), reactions)
    reactions = map(s -> replace(s, "arbitrary" => ""), reactions)

    #=
        Parsing for the two ODE formats is implemented:
            - ERODE "arrows" format for biological processes
            - Ordinary "equations" format

        These are mainly distinct by the presence of the '->' token,
        and we will try to dispatch on it
    =#

    sep = occursin("->", prod(reactions)) ? "->" : "="

    # scan reactions to discover all variables

    sep1 = Regex("($sep)|(,)|([+-/ \t\\*\\(\\)\\[\\]])")
    # assuming there is at least one utf letter in any variable
    strings = map(String, unique(filter!(
        x -> length(x) > 0 && isletter(x[1]),
        split(join(reactions, ' '), sep1)
    )))

    if length(strings) > 1000
        throw(ParseException("too long system encountered at $filepath, skipping it"))
    end

    # we want to agree with CLUE in variable order
    strings = sort(strings)
    #
    @debug "Variables parsed " * join(strings, ", ")
    if isempty(strings)
        throw(ParseException("bad file encountered at $filepath, skipping it"))
    end
    S, xs = Nemo.QQ[strings...]

    # symbol :x to x from QQ[x]
    mapping = Dict{Symbol, fmpq_mpoly}(
        Symbol(x) => gen(S, i)
        for (i, x) in enumerate(strings)
    )

    ODEs = Dict(x => S(0) for x in xs)

    if sep === "->"
        sep2 = r"(->)|(,)"
        for reaction in reactions
            # divide each line into parts,
            # expected to be three of them
            lhs, rhs, speed = map(
                Meta.parse ∘ String ∘ strip,
                split(reaction, sep2)
            )

            reagents, products, speed = map(
                x -> myeval(x, mapping),
                (lhs, rhs, speed)
            )

            concentration = prod(terms(reagents))

            if typeof(concentration * speed) <: AbstractAlgebra.Generic.Frac
                throw(
                    ParseException("We do not work with rational concentrtions yet =(")
                )
            end

            for reagent in monomials(reagents)
                ODEs[reagent] -= concentration * speed
            end
            for product in monomials(products)
                ODEs[product] += concentration * speed
            end
        end
    else
        sep2 = r"(=)|(,)"
        for reaction in reactions
            # divide each line into lhs, rhs
            lhs, rhs = map(
                Meta.parse ∘ String ∘ strip,
                split(reaction, sep2)
            )

            ẋ, rhs = map(
                x -> myeval(x, mapping),
                (lhs, rhs)
            )

            ODEs[ẋ] = rhs
        end
    end

    # @info "loaded a system of $(length(ODEs)) ODEs from $filepath"

    return ODEs
end

ext: jl


text:
export cgls

function cgls{T1,T2}(A::SparseMatrixCSC{T1,Int},b::Array{T2,1}; kwargs...)
    T  = promote_type(T1,T2)
    x1 = zeros(T,size(A,1))
    x2 = zeros(T,size(A,2))

    Af(x,flag) = (flag=='F') ? A_mul_B!(1.0,A,x,0.0,x1) : At_mul_B!(1.0,A,x,0.0,x2)
    return cgls(Af,b;kwargs...)
end

cgls(A,b::Vector;kwargs...) = cgls((x,flag) -> ((flag=='F') ? A*x : A'*x),b::Vector;kwargs...)

"""
x,flag,err,iter,resvec = cg(A,b,tol=1e-2,maxIter=100,x=[],storeInterm=false,out=0)

CGLS Conjugate gradient algorithm applied implicitly to the normal equations

		(A'*A) x = A'*b.

Input:

  A            - function computing A*x = A(x,'F') and A'*x = A(x,'T')
  b            - right hand side vector
  tol          - error tolerance, default 1e-2
  maxIter      - maximum number of iterations, default 100
  x            - starting guess
  storeInterm  - flag for returning intermediate solutions (useful in inverse problems)
  out          - flag for output (-1: no output, 0 : only errors, 1 : final status, 2: error at each iteration)

Output:

  x       - approximate solution (interm==0) or history of approximate solutions (interm==1)
  flag    - exit flag (  0 : desired tolerance achieved,
                        -1 : maxIter reached without converging
                        -2 : Matrix A is not positive definite
                        -9 : right hand side was zero)
  eta     - residual norm: norm(A*x-b)
  rho     - norm of current iterate: norm(x)

"""
function cgls(A::Function,b::Vector; tol::Real=1e-2,maxIter::Int=100,x::Vector=[],storeInterm::Bool=false,out::Int=0)

    m = length(b)
    # Initialization.
    if norm(b)==0; return 0.0,-9,0.0,0.0,[0.0]; end
    if isempty(x) || all(x.==0.0)
        r = copy(b)			# residual r = b - A*x
        s = A(r,'T')		# compute gradient g = A'*(A*x-b)
        n = length(s)
        x = zeros(n)
    else
        r = b-A(x,'F')
        s = A(r,'T')
        n = length(s)
    end

    if storeInterm
        X = zeros(n,maxIter)	# allocate space for intermediates
    end

    p      = copy(s)
    normS0 = BLAS.dot(n,s,1,s,1)       # squared norm of gradients
    normSc = normS0
    Arn    = zeros(maxIter)
    eta    = zeros(maxIter) # norm of residuals
    rho    = zeros(maxIter) # norm or current iterate

    # Iterate.
    flag = -1
    if out==2
        println("=== cgls ===")
        println(@sprintf("%4s\t%8s\t%8s\t%8s","iter","|A'r|","norm(r)","norm(x)"))
    end

    iter = 1 # makes iter available outside the loop
    for iter=1:maxIter
        q     = A(p,'F') # compute A*g
        alpha = normSc/BLAS.dot(m,q,1,q,1)
        BLAS.axpy!(n,alpha,p,1,x,1) # faster than x    += alpha*p

        if storeInterm; X[:,iter] = x; end

        BLAS.axpy!(m,-alpha,q,1,r,1) # faster than r  -= alpha*q
        s   = A(r,'T') # compute gradient, that is A'*r

        normSt = BLAS.dot(n,s,1,s,1)
        if (iter>1) && (normSt <= tol)
            resize!(eta,iter-1)
            resize!(rho,iter-1)
            resize!(Arn,iter-1)
            flag = 0; break
        end

        beta = normSt/normSc

        BLAS.scal!(n,beta,p,1)
        BLAS.axpy!(n,1.0,s,1,p,1)
        # faster than p   = s + beta*p

        # store intermediates and report resuls
        normSc = normSt
        Arn[iter] = normSt
        eta[iter] = BLAS.nrm2(m,r,1) # faster than norm(r)
        rho[iter] = BLAS.nrm2(n,x,1) # faster than norm(x)
        if out==2;  println(@sprintf("%3d\t%1.2e\t%1.2e\t%1.2e",iter,Arn[iter],rho[iter],eta[iter]));end
    end

    if out>=0
        if flag==-1
            println(@sprintf("cgls iterated maxIter (=%d) times witout reaching tolerance. Returned result has residual norm  %1.2e.",
            																										maxIter,Arn[end]))
        elseif flag==0 && out>=1
            println(@sprintf("cgls relative gradient norm below desired tolerance at iteration %d. Returned result has residual norm  %1.2e.",
            																										iter,Arn[end]))
        end
    end

    if storeInterm
        return X[:,1:iter],flag,rho,eta,Arn
    else
        return x,flag,rho,eta,Arn
    end
end


text:
"""
    FieldModule

Module for abstract fields.
"""
module FieldModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import Base.copyto!

"""
    Abstract field.

Expected  attributes:
  + values::FMat{T} # Array of degree of freedom parameters,  indexed by node
  + dofnums::FIntMat # Array of degree of freedom numbers, indexed by node
  + is_fixed::Matrix{Bool} # Array of Boolean flags, indexed by node
  + fixed_values::FMat{T} # Array of fixed values, indexed by node
  + nfreedofs::FInt # Total number of free degrees of freedom
See the macro `add_Field_fields()`.
"""
abstract type Field end

macro add_Field_fields()
    return esc(:(
    values::Array{T, 2};
    dofnums::Array{FInt, 2};
    is_fixed::Matrix{Bool};
    fixed_values::Array{T, 2};;
    nfreedofs::FInt;
    )
    )
end

"""
    ndofs{F<:Field}(self::T)

Dimension of the degree of freedom parameters (i. e.  how many degrees of
freedom per entity).
"""
ndofs(self::Field)  = size(self.values, 2)

"""
    nents{F<:Field}(self::F)

Number of nodes associated with the field.
"""
nents(self::Field)  = size(self.values, 1)

"""
    copyto!(DEST::F,  SRC::F) where {F<:Field}

Copy data from one field to another.
"""
function copyto!(DEST::F,  SRC::F) where {F<:Field}
    copyto!(DEST.values, SRC.values)
    copyto!(DEST.dofnums, SRC.dofnums)
    copyto!(DEST.is_fixed, SRC.is_fixed)
    copyto!(DEST.fixed_values, SRC.fixed_values)
    DEST.nfreedofs = SRC.nfreedofs
    return  DEST
end

"""
    wipe!(self::Field)

Wipe all the data from the field.
"""
function wipe!(self::Field)
    Zer = zero(eltype(self.fixed_values[1]))
    self.nfreedofs = 0
    fill!(self.dofnums, 0)
    fill!(self.is_fixed, false)
    fill!(self.values, Zer)
    fill!(self.fixed_values, Zer)
    return self
end

"""
    gathersysvec{F<:Field}(self::F)

Gather values from the field for the whole system vector.
"""
function gathersysvec(self::F) where {F<:Field}
    nents,dim = size(self.values)
    vec = zeros(typeof(self.values[1,1]),self.nfreedofs,1)
    for i = 1:nents
        for j = 1:dim
            en = self.dofnums[i,j]
            if (en > 0) && (en <= self.nfreedofs)
                vec[en] = self.values[i,j]
            end
        end
    end
    return vec
end

"""
    gathersysvec!(self::F, vec::FVec{T}) where {F<:Field, T}

Gather values from the field for the whole system vector.
"""
function gathersysvec!(self::F, vec::FVec{T}) where {F<:Field, T}
    nents,dim = size(self.values)
    @assert length(vec) == self.nfreedofs
    for i = 1:nents
        for j = 1:dim
            en = self.dofnums[i,j]
            if (en > 0) && (en <= self.nfreedofs)
                vec[en] = self.values[i,j]
            end
        end
    end
    return vec
end

"""
    gathervalues_asvec!(self::Field, dest::AbstractArray{T, 1},
        conn::CC) where {CC, T}

Gather values from the field into a vector.

The order is: for each node  in the connectivity, copy into the buffer all the
degrees of freedom,  then the next node and so on.

`dest` = destination buffer: overwritten  inside,  must be preallocated
in the correct size
"""
function gathervalues_asvec!(self::Field, dest::AbstractArray{T, 1}, conn::CC) where {CC, T}
    en::FInt = 1;
    for i = 1:length(conn)
        for j = 1:size(self.values,2)
            dest[en] = self.values[conn[i],j];
            en = en + 1;
        end
    end
    return dest
end

"""
    gathervalues_asmat!(self::Field, dest::AbstractArray{T, 2},
        conn::CC) where {CC, T}

Gather values from the field into a two-dimensional array.

The order is: for each node  in the connectivity, copy into the corresponding
row of the buffer all the degrees of freedom,  then the next node into the next
row and so on.

`dest` = destination buffer: overwritten  inside,  must be preallocated
in the correct size
"""
function gathervalues_asmat!(self::Field, dest::AbstractArray{T, 2},    conn::CC) where {CC, T}
    for i = 1:length(conn)
        for j = 1:size(self.values,2)
            dest[i, j] = self.values[conn[i], j];
        end
    end
    return dest
end

"""
    gatherfixedvalues_asvec!(self::Field, dest::AbstractArray{T, 1},
        conn::CC) where {CC, T}

Gather FIXED values from the field into a vector.

The order is: for each node  in the connectivity, copy into the buffer all the
fixed degrees of freedom,  then the next node and so on. If a degree of freedom
is NOT fixed, the corresponding entry is  set to zero.

`dest` = destination buffer: overwritten  inside,  must be preallocated
in the correct size
"""
function gatherfixedvalues_asvec!(self::Field, dest::AbstractArray{T, 1},    conn::CC) where {CC, T}
    en::FInt = 1;
    for i = 1:length(conn)
        for j = 1:size(self.fixed_values,2)
            if self.is_fixed[conn[i],j] # free degree of freedom
                dest[en] = self.fixed_values[conn[i], j];
            else
                dest[en] = 0.0
            end
            en = en + 1;
        end
    end
    return dest
end

"""
    gatherfixedvalues_asmat!(self::Field, dest::AbstractArray{T, 2},
        conn::CC) where {CC, T}

Gather FIXED values from the field into a two-dimensional array.

The order is: for each node  in the connectivity, copy into the corresponding
row of the buffer all the degrees of freedom,  then the next node into the next
row and so on.  If a degree of freedom
is NOT fixed, the corresponding entry is  set to zero.

`dest` = destination buffer: overwritten  inside,  must be preallocated
in the correct size
"""
function gatherfixedvalues_asmat!(self::Field, dest::AbstractArray{T, 2},    conn::CC) where {CC, T}
    for i = 1:length(conn)
        for j = 1:size(self.fixed_values,2)
            if self.is_fixed[conn[i],j] # fixed degree of freedom
                dest[i, j] = self.fixed_values[conn[i], j];
            else
                dest[i, j] = 0.0
            end
        end
    end
    return dest
end

function anyfixedvaluenz(self::Field, conn::CC) where {CC}
    for i = 1:length(conn)
        for j = 1:size(self.fixed_values,2)
            if self.is_fixed[conn[i],j] # free degree of freedom
                if  abs(self.fixed_values[conn[i], j]) > 0.0
                    return true
                end
            end
        end
    end
    return false
end

"""
    gatherdofnums!(self::Field, dest::A, conn::CC) where {A, CC}

Gather dofnums from the field.
"""
function gatherdofnums!(self::Field, dest::A, conn::CC) where {A, CC}
    en::FInt = 1;
    for i = 1:length(conn)
        for j = 1:size(self.dofnums,2)
            dest[en] = self.dofnums[conn[i],j];
            en = en+1;
        end
    end
    return dest
end

"""
    numberdofs!(self::Field)

Number the degrees of freedom.

The free components in the field are numbered consecutively. No effort is
made to optimize the numbering in any way. If you'd like to optimize the
numbering of the degrees of freedom, use the above form that sets the
permutation of the degrees of freedom, or the permutation of the nodes.
"""
function numberdofs!(self::Field)
    fixed_dofnum::FInt = 0
    nents,dim = size(self.values)
    self.nfreedofs::FInt =0
    for i=1:nents
        for j=1:dim
            if !self.is_fixed[i,j] # free degree of freedom
                self.nfreedofs = self.nfreedofs + 1
                self.dofnums[i,j] = self.nfreedofs
            else # fixed degree of freedom: no equation
                self.dofnums[i,j] = fixed_dofnum
            end
        end
    end
    return  self
end

"""
    setebc!(self::Field, fenids::FIntVec, is_fixed::Bool, comp::FInt,
      val::FVec{T}) where {T<:Number}

Set the EBCs (essential boundary conditions).

`fenids`         - array of N node identifiers
`is_fixed` = scaler Boolean: are the degrees of freedom being fixed (true)
             or released (false),
`comp` = integer, which  degree of freedom (component),
`val` = array of N values of type T

Note:  Any call to setebc!() potentially changes the current assignment
which degrees of freedom are free and which are fixed
and therefore is presumed to invalidate the
current degree-of-freedom numbering. In such a case this method sets
`nfreedofs = 0`; and  `dofnums=0`.
"""
function setebc!(self::Field, fenids::FIntVec, is_fixed::Bool, comp::FInt,    val::FVec{T}) where {T<:Number}
    @assert comp <= size(self.values,2) "Requested  nonexistent  degree of freedom"
    @assert maximum(fenids) <= size(self.values,1) "Requested nonexistent node"
    @assert size(fenids) == size(val) "Arrays of mismatched sizes"
    for  j = 1:length(fenids)
        self.is_fixed[fenids[j],comp] = is_fixed;
        if self.is_fixed[fenids[j],comp]
            self.fixed_values[fenids[j],comp] = val[j];
        else
            self.fixed_values[fenids[j],comp] = zero(T)
        end
    end
    self.nfreedofs = 0
    fill!(self.dofnums, 0)
    return  self
end

"""
    setebc!(self::Field, fenids::FIntVec, is_fixed::Bool, comp::FInt,
      val::T) where {T<:Number}

Set the EBCs (essential boundary conditions).

`fenids`         - array of N node identifiers
`is_fixed` = scaler Boolean: are the degrees of freedom being fixed (true)
             or released (false),
`comp` = integer, which  degree of freedom (component),
`val` = scalar of type T

Note:  Any call to setebc!() potentially changes the current assignment
which degrees of freedom are free and which are fixed
and therefore is presumed to invalidate the
current degree-of-freedom numbering. In such a case this method sets
`nfreedofs = 0`; and  `dofnums=0`.
"""
function setebc!(self::Field, fenids::FIntVec, is_fixed::Bool, comp::FInt,    val::T) where {T<:Number}
    @assert (comp >= 1 && comp <= size(self.values,2)) "Requested  nonexistent  degree of freedom"
    @assert maximum(fenids) <= size(self.values,1) "Requested nonexistent node"
    @assert minimum(fenids) >= 1 "Requested nonexistent node"
    for  j = 1:length(fenids)
        self.is_fixed[fenids[j],comp] = is_fixed;
        if self.is_fixed[fenids[j],comp]
            self.fixed_values[fenids[j],comp] = val;
        else
            self.fixed_values[fenids[j],comp] = zero(T)
        end
    end
    self.nfreedofs = 0
    fill!(self.dofnums, 0)
    return  self
end

"""
    setebc!(self::Field, fenids::FIntVec, comp::FInt,
      val::FVec{T}) where {T<:Number}

Set the EBCs (essential boundary conditions).

`fenids`         - array of N node identifiers
`comp` = integer, which  degree of freedom (component),
`val` = array of N values of type T

Note:  Any call to setebc!() potentially changes the current assignment
which degrees of freedom are free and which are fixed
and therefore is presumed to invalidate the
current degree-of-freedom numbering. In such a case this method sets
`nfreedofs = 0`; and  `dofnums=0`.
"""
function setebc!(self::Field, fenids::FIntVec, comp::FInt, val::FVec{T}) where {T<:Number}
    return setebc!(self, fenids, true, comp, val)
end

"""
    setebc!(self::Field, fenids::FIntVec, comp::FInt;
      val::T=0.0) where {T<:Number}

Set the EBCs (essential boundary conditions).

`fenids`         - array of N node identifiers
`comp` = integer, which  degree of freedom (component),
`val` = scalar of type T

Note:  Any call to setebc!() potentially changes the current assignment
which degrees of freedom are free and which are fixed
and therefore is presumed to invalidate the
current degree-of-freedom numbering. In such a case this method sets
`nfreedofs = 0`; and  `dofnums=0`.
"""
function setebc!(self::Field, fenids::FIntVec, comp::FInt; val::T=0.0) where {T<:Number}
    return setebc!(self, fenids, true, comp, val)
end

"""
    setebc!(self::Field, fenids::FIntVec, is_fixed::Bool, comp::FInt;
      val::T=0.0) where {T<:Number}

Set the EBCs (essential boundary conditions).

`fenids`         - array of N node identifiers
`comp` = integer, which  degree of freedom (component),
`val` = scalar of type T

Note:  Any call to setebc!() potentially changes the current assignment
which degrees of freedom are free and which are fixed
and therefore is presumed to invalidate the
current degree-of-freedom numbering. In such a case this method sets
`nfreedofs = 0`; and  `dofnums=0`.
"""
function setebc!(self::Field, fenids::FIntVec, is_fixed::Bool, comp::FInt;    val::T=0.0) where {T<:Number}
    j = comp
    @assert (j >= 1) && (j <= ndofs(self))
    setebc!(self, fenids, is_fixed, j, val)
    return self
end

"""
    setebc!(self::Field, fenids::FIntVec)

Set the EBCs (essential boundary conditions).

Suppress all degrees of freedom at the given nodes.

`fenids`         - array of N node identifiers

Note:  Any call to setebc!() potentially changes the current assignment
which degrees of freedom are free and which are fixed
and therefore is presumed to invalidate the
current degree-of-freedom numbering. In such a case this method sets
`nfreedofs = 0`; and  `dofnums=0`.
"""
function setebc!(self::Field, fenids::FIntVec)
    Zer = zero(eltype(self.fixed_values[1]))
    for comp = 1:size(self.values, 2)
        setebc!(self, fenids, true, comp, Zer)
    end
    return self
end

"""
    setebc!(self::Field, fenid::FInt)

Set the EBCs (essential boundary conditions).

Suppress all degrees of freedom at the given node.

`fenid`         - One integer as a node identifier

Note:  Any call to setebc!() potentially changes the current assignment
which degrees of freedom are free and which are fixed
and therefore is presumed to invalidate the
current degree-of-freedom numbering. In such a case this method sets
`nfreedofs = 0`; and  `dofnums=0`.
"""
function setebc!(self::Field, fenid::FInt)
    Zer = zero(eltype(self.fixed_values[1]))
    for comp = 1:size(self.values, 2)
        setebc!(self, [fenid], true, comp, Zer)
    end
    return self
end

"""
    setebc!(self::Field)

Set the EBCs (essential boundary conditions).

All essential boundary conditions are CLEARED.

Note:  Any call to setebc!() potentially changes the current assignment
which degrees of freedom are free and which are fixed
and therefore is presumed to invalidate the
current degree-of-freedom numbering. In such a case this method sets
`nfreedofs = 0`; and  `dofnums=0`.
"""
function setebc!(self::Field)
    self.nfreedofs = 0
    fill!(self.dofnums, 0)
    fill!(self.is_fixed, false)
    fill!(self.fixed_values, zero(eltype(self.fixed_values[1])))
    return  self
end

"""
    applyebc!(self::Field)

Apply EBCs (essential boundary conditions).
"""
function applyebc!(self::Field)
    nents,dim = size(self.values);
    for i = 1:nents
        for j = 1:dim
            if self.is_fixed[i,j]
                self.values[i,j] = self.fixed_values[i,j];
            end
        end
    end
    return  self
end

"""
    scattersysvec!(self::Field, vec::FVec{T}) where {T<:Number}

Scatter values to the field from a system vector.
"""
function scattersysvec!(self::Field, vec::FVec{T}) where {T<:Number}
    nents,dim = size(self.values);
    for i = 1:nents
        for j = 1:dim
            dn = self.dofnums[i,j];
            if (dn > 0) && (dn <= self.nfreedofs)
                self.values[i,j] = vec[dn];
            end
        end
    end
    return  self
end

end


text:
import Base: copy, ctranspose, getindex, showarray, transpose, one, zero
import Base.LinAlg: SingularException

immutable UniformScaling{T<:Number}
    λ::T
end

const I = UniformScaling(1)

eltype{T}(J::UniformScaling{T}) = T
ndims(J::UniformScaling) = 2
getindex(J::UniformScaling, i::Integer,j::Integer) = ifelse(i==j,J.λ,zero(J.λ))

show(io::IO, J::UniformScaling) = print(io, "$(typeof(J))\n$(J.λ)*I")
copy(J::UniformScaling) = UniformScaling(J.λ)

transpose(J::UniformScaling) = J
ctranspose(J::UniformScaling) = UniformScaling(conj(J.λ))

one{T}(::Type{UniformScaling{T}}) = UniformScaling(one(T))
one{T}(J::UniformScaling{T}) = one(UniformScaling{T})
zero{T}(::Type{UniformScaling{T}}) = UniformScaling(zero(T))
zero{T}(J::UniformScaling{T}) = zero(UniformScaling{T})

+(J1::UniformScaling, J2::UniformScaling) = UniformScaling(J1.λ+J2.λ)
+{T}(B::BitArray{2},J::UniformScaling{T}) = bitunpack(B) + J
+(J::UniformScaling, B::BitArray{2})      = J + bitunpack(B)
+(J::UniformScaling, A::AbstractMatrix)   = A + J
+(J::UniformScaling, x::Number)           = J.λ + x
+(x::Number, J::UniformScaling)           = x + J.λ

function +{TA,TJ}(A::AbstractMatrix{TA}, J::UniformScaling{TJ})
    n = chksquare(A)
    B = similar(A, promote_type(TA,TJ))
    copy!(B,A)
    @inbounds for i = 1:n
        B[i,i] += J.λ
    end
    B
end

-(J::UniformScaling)                      = UniformScaling(-J.λ)
-(J1::UniformScaling, J2::UniformScaling) = UniformScaling(J1.λ-J2.λ)
-(B::BitArray{2}, J::UniformScaling)      = bitunpack(B) - J
-(J::UniformScaling, B::BitArray{2})      = J - bitunpack(B)
-(J::UniformScaling, x::Number)           = J.λ - x
-(x::Number, J::UniformScaling)           = x - J.λ

function -{TA,TJ<:Number}(A::AbstractMatrix{TA}, J::UniformScaling{TJ})
    n = chksquare(A)
    B = similar(A, promote_type(TA,TJ))
    copy!(B, A)
    @inbounds for i = 1:n
        B[i,i] -= J.λ
    end
    B
end
function -{TA,TJ<:Number}(J::UniformScaling{TJ}, A::AbstractMatrix{TA})
    n = chksquare(A)
    B = -A
    @inbounds for i = 1:n
        B[i,i] += J.λ
    end
    B
end

*(J1::UniformScaling, J2::UniformScaling) = UniformScaling(J1.λ*J2.λ)
*(B::BitArray{2}, J::UniformScaling) = *(bitunpack(B), J::UniformScaling)
*(J::UniformScaling, B::BitArray{2}) = *(J::UniformScaling, bitunpack(B))
*(S::SparseMatrixCSC, J::UniformScaling) = J.λ == 1 ? S : J.λ*S
*{Tv,Ti}(J::UniformScaling, S::SparseMatrixCSC{Tv,Ti}) = J.λ == 1 ? S : S*J.λ
*(A::AbstractMatrix, J::UniformScaling) = J.λ == 1 ? A : J.λ*A
*(J::UniformScaling, A::AbstractVecOrMat) = J.λ == 1 ? A : J.λ*A

*(x::Number, J::UniformScaling) = UniformScaling(x*J.λ)
*(J::UniformScaling, x::Number) = UniformScaling(J.λ*x)

/(J1::UniformScaling, J2::UniformScaling) = J2.λ == 0 ? throw(SingularException(1)) : UniformScaling(J1.λ/J2.λ)
/(J::UniformScaling, A::AbstractMatrix) = J.λ*inv(A)
/(A::AbstractMatrix, J::UniformScaling) = J.λ == 0 ? throw(SingularException(1)) : A/J.λ

/(J::UniformScaling, x::Number) = UniformScaling(J.λ/x)

\(J1::UniformScaling, J2::UniformScaling) = J1.λ == 0 ? throw(SingularException(1)) : UniformScaling(J1.λ\J2.λ)
\{T<:Number}(A::Union(Bidiagonal{T},Triangular{T}), J::UniformScaling) = inv(A)*J.λ
\(J::UniformScaling, A::AbstractVecOrMat) = J.λ == 0 ? throw(SingularException(1)) : J.λ\A
\(A::AbstractMatrix, J::UniformScaling) = inv(A)*J.λ

\(x::Number, J::UniformScaling) = UniformScaling(x\J.λ)

.*(x::Number,J::UniformScaling) = UniformScaling(x*J.λ)
.*(J::UniformScaling,x::Number) = UniformScaling(J.λ*x)

./(J::UniformScaling,x::Number) = UniformScaling(J.λ/x)

==(J1::UniformScaling,J2::UniformScaling) = (J1.λ == J2.λ)

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2021 Varnerlab
# Robert Frederick Smith School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED AS IS, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #

function Balances(dx,x, parameter_dictionary,t)

    # system dimensions and structural matricies -
    number_of_states = parameter_dictionary["number_of_states"]
    DM = parameter_dictionary["dilution_degradation_matrix"]
    SM = parameter_dictionary["stoichiometric_matrix"]
    half_life_translation = parameter_dictionary["half_life_translation_capacity"]

    # calculate the TX and TL kinetic limit array -
    transcription_kinetic_limit_array = calculate_transcription_kinetic_limit_array(t,x,parameter_dictionary)
    translation_kinetic_limit_array = calculate_translation_kinetic_limit_array(t,x,parameter_dictionary)

    # calculate the TX and TL control array -
    u = calculate_transcription_control_array(t,x,parameter_dictionary)
    w = calculate_translation_control_array(t,x,parameter_dictionary)

    # calculate the rate of transcription and translation -
    rX = transcription_kinetic_limit_array.*u
    rL = translation_kinetic_limit_array.*w
    rV = [rX ; rL]

    # calculate the degradation and dilution rates -
    rd = calculate_dilution_degradation_array(t,x,parameter_dictionary)

    # compute the model equations -
    dxdt = SM*rV + DM*rd

    # package -
    for index = 1:number_of_states
        dx[index] = dxdt[index]
    end

    # extra species: global translation capacity -
	dx[4] = -(log(2)*(half_life_translation^-1))*x[4]
end


text:
############################################################
# Identity Layer
############################################################
@defstruct IdentityLayer Layer (
  name :: AbstractString = "identity",
  (tops :: Vector{Symbol} = Symbol[], length(tops) > 0),
  (bottoms :: Vector{Symbol} = Symbol[], length(bottoms) == length(tops)),
)
@characterize_layer(IdentityLayer,
  can_do_bp => true
)

type IdentityLayerState <: LayerState
  layer      :: IdentityLayer
  blobs      :: Vector{Blob}
  blobs_diff :: Vector{Blob}
end

function setup(backend::Backend, layer::IdentityLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  blobs      = inputs[:] # shallow copy
  blobs_diff = diffs[:] # shallow_copy

  IdentityLayerState(layer, blobs, blobs_diff)
end
function shutdown(backend::Backend, state::IdentityLayerState)
end

function forward(backend::Backend, state::IdentityLayerState, inputs::Vector{Blob})
  # do nothing
end

function backward(backend::Backend, state::IdentityLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})
  # do nothing
end


text:
#
# Basic import for any visualisation that runs the full model.
# Create a base set of results, settings and parameters and
# import everything we might reasonably need.#
#
# fixme extend to multiple systems

function load_system()::TaxBenefitSystem
	sys = load_file( joinpath( Definitions.MODEL_PARAMS_DIR, "sys_2021_22.jl" ))
	#
	# Note that as of Budget21 removing these doesn't actually happen till May 2022.
	#
	load_file!( sys, joinpath( Definitions.MODEL_PARAMS_DIR, "sys_2021-uplift-removed.jl"))
	# uc taper to 55
	load_file!( sys, joinpath( Definitions.MODEL_PARAMS_DIR, "budget_2021_uc_changes.jl"))
	weeklyise!( sys )
	return sys
end

function do_run_a(
	cache_key,
	sys :: TaxBenefitSystem,
	settings :: Settings ) :: AllOutput
	@debug "do_run_a entered"
	obs = Observable(
		Progress(settings.uuid, "",0,0,0,0))
	tot = 0
	of = on(obs) do p
		tot += p.step
		PROGRESS[p.uuid] = (progress=p,total=tot)
	end
	results = do_one_run( settings, [sys], obs )
	settings.poverty_line = make_poverty_line( results.hh[1], settings )

	outf = summarise_frames( results, settings )
	gl = make_gain_lose( BASE_RESULTS.results.hh[1], results.hh[1], settings )
	exres = calc_examples( BASE_PARAMS, sys, settings )
	aout = AllOutput( settings.uuid, cache_key, results, outf, gl, exres )
	return aout;
end

"""
Runner version used in scotbudg without the caching stuff
"""
function do_run( sys :: TaxBenefitSystem, init = false )::NamedTuple
	settings = deepcopy( BASE_SETTINGS )
	settings.uuid = UUIDs.uuid4()
	obs = Observable(Progress(settings.uuid, "",0,0,0,0))
	tot = 0
	of = on(obs) do p
		tot += p.step
		PROGRESS[p.uuid] = (progress=p,total=tot)
	end
    results = do_one_run( settings, [sys], obs )
	outf = summarise_frames( results, settings )
	gl = make_gain_lose( BASE_RESULTS.results.hh[1], results.hh[1], settings )
	delete!( PROGRESS, settings.uuid )
	return (results=results, summary=outf,gain_lose=gl  )
end


text:
"""
    field_default(::Type{T}, name::Symbol)

Return the default value of field `name` of an option type `T`.
"""
function field_default(::Type{T}, name::Symbol) where {T}
    if !isconcretetype(T) && is_option(T)
        error("field_default requires a concrete type, got $T{...}")
    else
        error("field_default is not defined for $T, it may not be an option type")
    end
end

"""
    field_defaults(::Type)

Return default values of given option types.
"""
function field_defaults(::Type{T}) where {T}
    is_option(T) || error("$T is not an option type")
    return Any[field_default(T, each) for each in fieldnames(T)]
end

"""
    type_alias(::Type{OptionType}) -> String

Return the alias name of given `OptionType`.
"""
function type_alias(::Type{T}) where {T}
    is_option(T) || error("$T is not an option type")
    return error("`type_alias` is not defined for $T, it may not be an option type")
end

function set_type_alias(::Type{T}, name::String) where {T}
    isconcretetype(T) || error("cannot set alias for non concrete type")
    type_alias_map = get_type_alias_map(T)::Dict{String, Any}
    type_alias_map[name] = T
    return
end

function get_type_alias_map(::Type{T}) where {T}
    error("`get_type_alias_map` is not defined for `$T`, it may not be an option type")
end


text:
"""
This package provides an interface to the POSIX calls for tty I/O control.

For a complete description of these calls, see the [termios(3)](http://man7.org/linux/man-pages/man3/termios.3.html) manual page.

**Usage**

Example:

    using TERMIOS
    const T = TERMIOS
    term_stdin = T.termios()
    T.tcgetattr(stdin, term_stdin)
    term_stdin.c_iflag |= T.IGNBRK
    T.tcsetattr(stdin, T.TCSANOW, term_stdin)

"""
module TERMIOS

struct TERMIOSError <: Exception
    msg::String
end

Base.showerror(io::IO, e::TERMIOSError) = print(io, "TERMIOSError: $(e.msg)")

const cc_t = Sys.islinux() ? Cuchar : Cvoid
const tcflag_t = Sys.islinux() ? Cuint : Culong
const char = Sys.islinux() ? Cchar : Cchar
const speed_t = Sys.islinux() ? Cuint : Culong

"""End-of-file character ICANON"""
const VEOF = Sys.islinux() ? 4 : 0

"""End-of-line character ICANON"""
const VEOL = Sys.islinux() ? 11 : 1

"""Second EOL character ICANON together with IEXTEN"""
const VEOL2 = Sys.islinux() ? 16 : 2

"""Erase character ICANON"""
const VERASE = Sys.islinux() ? 2 : 3

"""Word-erase character ICANON together with IEXTEN"""
const VWERASE = Sys.islinux() ? 14 : 4

"""Kill-line character ICANON"""
const VKILL = Sys.islinux() ? 3 : 5

"""Reprint-line character ICANON together with IEXTEN"""
const VREPRINT = Sys.islinux() ? 12 : 6

"""Interrupt character ISIG"""
const VINTR = Sys.islinux() ? 0 : 8

"""Quit character ISIG"""
const VQUIT = Sys.islinux() ? 1 : 9

"""Suspend character ISIG"""
const VSUSP = Sys.islinux() ? 10 : 10

"""Delayed suspend character ISIG together with IEXTEN"""
const VDSUSP = Sys.islinux() ? nothing : 11

"""Start (X-ON) character IXON, IXOFF"""
const VSTART = Sys.islinux() ? 8 : 12

"""Stop (X-OFF) character IXON, IXOFF"""
const VSTOP = Sys.islinux() ? 9 : 13

"""Literal-next character IEXTEN"""
const VLNEXT = Sys.islinux() ? 15 : 14

"""Discard character IEXTEN"""
const VDISCARD = Sys.islinux() ? 13 : 15

"""Minimum number of bytes read at once !ICANON"""
const VMIN = Sys.islinux() ? 6 : 16

"""Time-out value (tenths of a second) !ICANON"""
const VTIME = Sys.islinux() ? 5 : 17

"""Status character ICANON together with IEXTEN"""
const VSTATUS = Sys.islinux() ? nothing : 18

const NCCS = Sys.islinux() ? 32 : 20

const VSWTC = Sys.islinux() ? 7 : nothing

#
# Input flags - software input processing
#

"""Ignore BREAK condition"""
const IGNBRK = 0x00000001

"""Map BREAK to SIGINTR"""
const BRKINT = 0x00000002

"""Ignore (discard) parity errors"""
const IGNPAR = 0x00000004

"""Mark parity and framing errors"""
const PARMRK = 0x00000008

"""Enable checking of parity errors"""
const INPCK = 0x00000010

"""Strip 8th bit off chars"""
const ISTRIP = 0x00000020

"""Map NL into CR"""
const INLCR = 0x00000040

"""Ignore CR"""
const IGNCR = 0x00000080

"""Map CR to NL (ala CRMOD)"""
const ICRNL = 0x00000100

"""Enable output flow control"""
const IXON = Sys.islinux() ? 0x00000400 : 0x00000200

"""Enable input flow control"""
const IXOFF = Sys.islinux() ? 0x00001000 : 0x00000400

"""Any char will restart after stop"""
const IXANY = 0x00000800

"""Ring bell on input queue full"""
const IMAXBEL = Sys.islinux() ? nothing : 0x00002000

"""(macos) maintain state for UTF-8 VERASE"""
const IUTF8 = Sys.islinux() ? nothing : 0x00004000

"""(glibc) Translate upper case input to lower case."""
const IUCLC = Sys.islinux() ? (1 << 14) : nothing

#
# Output flags - software output processing
#

"""Enable following output processing """
const OPOST = 0x00000001

"""Map NL to CR-NL (ala CRMOD)"""
const ONLCR = Sys.islinux() ? 0x00000004 : 0x00000002

"""Expand tabs to spaces"""
const OXTABS = 0x00000004

"""Discard EOT's (^D) on output)"""
const ONOEOT = 0x00000008

"""Map CR to NL on output"""
const OCRNL = Sys.islinux() ? 0x00000008 : 0x00000010

"""No CR output at column 0"""
const ONOCR = Sys.islinux() ? 0x00000010 : 0x00000020

"""NL performs CR function"""
const ONLRET = Sys.islinux() ? 0x00000020 : 0x00000040

raw"""\n delay"""
const NLDLY = Sys.islinux() ? 0x00000100 : 0x00000300

"""Horizontal tab delay"""
const TABDLY = Sys.islinux() ? 0x00001800 : 0x00000c04

raw"""\r delay"""
const CRDLY = Sys.islinux() ? 0x00000600 : 0x00003000

"""Form feed delay"""
const FFDLY = Sys.islinux() ? 0x00008000 : 0x00004000

raw"""\b delay"""
const BSDLY = Sys.islinux() ? 0x00002000 : 0x00008000

"""Vertical tab delay"""
const VTDLY = Sys.islinux() ? 0x00004000 : 0x00010000

"""NL type 0."""
const NL0 = 0x00000000

"""NL type 1."""

const NL1 = 0x00000100
const NL2 = Sys.islinux() ? nothing : 0x00000200
const NL3 = Sys.islinux() ? nothing : 0x00000300

"""TAB delay type 0."""
const TAB0 = 0x00000000

"""TAB delay type 1."""
const TAB1 = Sys.islinux() ? 0x00000800 : 0x00000400

"""TAB delay type 2."""
const TAB2 = Sys.islinux() ? 0x00001000 : 0x00000800

"""Expand tabs to spaces."""
const TAB3 = Sys.islinux() ? 0x00001800 : 0x00000004

"""CR delay type 0."""
const CR0 = 0x00000000

"""CR delay type 1."""
const CR1 = Sys.islinux() ? 0x00000200 : 0x00001000

"""CR delay type 2."""
const CR2 = Sys.islinux() ? 0x00000400 : 0x00002000

"""CR delay type 3."""
const CR3 = Sys.islinux() ? 0x00000600 : 0x00003000

"""FF delay type 0."""
const FF0 = 0x00000000

"""FF delay type 1."""
const FF1 = Sys.islinux() ? 0x00008000 : 0x00004000

"""BS delay type 0."""
const BS0 = 0x00000000

"""BS delay type 1."""
const BS1 = Sys.islinux() ? 0x00002000 : 0x00008000

"""VT delay type 0."""
const VT0 = 0x00000000

"""VT delay type 1."""
const VT1 = Sys.islinux() ? 0x00004000 : 0x00010000

"""(glibc) Translate lower case output to upper case."""
const OLCUC = Sys.islinux() ? (1 << 17) : nothing

"""Use fill characters for delay"""
const OFILL = Sys.islinux() ? 0x00000040 : 0x00000080

"""Fill is DEL, else NUL"""
const OFDEL = Sys.islinux() ? nothing : 0x00020000

#
# Control flags - hardware control of terminal
#

"""Ignore control flags"""
const CIGNORE = Sys.islinux() ? nothing : 0x00000001

"""5 bits (pseudo)"""
const CS5 = 0x00000000

"""6 bits"""
const CS6 = Sys.islinux() ? 0x00000010 : 0x00000100

"""7 bits"""
const CS7 = Sys.islinux() ? 0x00000020 : 0x00000200

"""8 bits"""
const CS8 = CS6 | CS7

"""Character size mask"""
const CSIZE = CS5 | CS6 | CS7 | CS8

"""Send 2 stop bits"""
const CSTOPB = Sys.islinux() ? 0x00000040 : 0x00000400

"""Enable receiver"""
const CREAD = Sys.islinux() ? 0x00000080 : 0x00000800

"""Parity enable"""
const PARENB = Sys.islinux() ? 0x00000100 : 0x00001000

"""Odd parity, else even"""
const PARODD = Sys.islinux() ? 0x00000200 : 0x00002000

"""Hang up on last close"""
const HUPCL = Sys.islinux() ? 0x00000400 : 0x00004000

"""Ignore modem status lines"""
const CLOCAL = Sys.islinux() ? 0x00000800 : 0x00008000

"""CTS flow control of output"""
const CCTS_OFLOW = Sys.islinux() ? nothing : 0x00010000

"""RTS flow control of input"""
const CRTS_IFLOW = Sys.islinux() ? nothing : 0x00020000

"""DTR flow control of input"""
const CDTR_IFLOW = Sys.islinux() ? nothing : 0x00040000

"""DSR flow control of output"""
const CDSR_OFLOW = Sys.islinux() ? nothing : 0x00080000

"""DCD flow control of output"""
const CCAR_OFLOW = Sys.islinux() ? nothing : 0x00100000

"""Old name for CCAR_OFLOW"""
const MDMBUF = Sys.islinux() ? nothing : 0x00100000
const CRTSCTS = Sys.islinux() ? nothing : (CCTS_OFLOW | CRTS_IFLOW)

#
# "Local" flags - dumping ground for other state
#
# Warning: some flags in this structure begin with
# the letter "I" and look like they belong in the
# input flag.
#

"""Visual erase for line kill"""
const ECHOKE = Sys.islinux() ? nothing : 0x00000001

"""Visually erase chars"""
const ECHOE = Sys.islinux() ? 0x00000010 : 0x00000002

"""Echo NL after line kill"""
const ECHOK = Sys.islinux() ? 0x00000020 : 0x00000004

"""Enable echoing"""
const ECHO = 0x00000008

"""Echo NL even if ECHO is off"""
const ECHONL = Sys.islinux() ? 0x00000040 : 0x00000010

"""Visual erase mode for hardcopy"""
const ECHOPRT = Sys.islinux() ? nothing : 0x00000020

"""Echo control chars as ^(Char)"""
const ECHOCTL = Sys.islinux() ? nothing : 0x00000040

"""Enable signals INTR, QUIT, [D]SUSP"""
const ISIG = Sys.islinux() ? 0x00000001 : 0x00000080

"""Canonicalize input lines"""
const ICANON = Sys.islinux() ? 0x00000002 : 0x00000100

"""Use alternate WERASE algorithm"""
const ALTWERASE = Sys.islinux() ? nothing : 0x00000200

"""Enable DISCARD and LNEXT"""
const IEXTEN = Sys.islinux() ? 0x00008000 : 0x00000400

"""External processing"""
const EXTPROC = Sys.islinux() ? nothing : 0x00000800

"""Stop background jobs from output"""
const TOSTOP = Sys.islinux() ? 0x00000100 : 0x00400000

"""Output being flushed (state)"""
const FLUSHO = Sys.islinux() ? nothing : 0x00800000

"""No kernel output from VSTATUS"""
const NOKERNINFO = Sys.islinux() ? nothing : 0x02000000

"""XXX retype pending input (state)"""
const PENDIN = Sys.islinux() ? nothing : 0x20000000

"""Don't flush after interrupt"""
const NOFLSH = Sys.islinux() ? 0x00000080 : 0x80000000

#
# Commands passed to tcsetattr() for setting the termios structure.
#

"""Make change immediate"""
const TCSANOW = 0

"""Drain output, then change"""
const TCSADRAIN = 1

"""Drain output, flush input"""
const TCSAFLUSH = 2

"""Flag - don't alter h.w. state"""
const TCSASOFT = Sys.islinux() ? nothing : 0x10

#
# Standard speeds
#

const B0 = Sys.islinux() ? 0 : 0
const B50 = Sys.islinux() ? 1 : 50
const B75 = Sys.islinux() ? 2 : 75
const B110 = Sys.islinux() ? 3 : 110
const B134 = Sys.islinux() ? 4 : 134
const B150 = Sys.islinux() ? 5 : 150
const B200 = Sys.islinux() ? 6 : 200
const B300 = Sys.islinux() ? 7 : 300
const B600 = Sys.islinux() ? 8 : 600
const B1200 = Sys.islinux() ? 9 : 1200
const B1800 = Sys.islinux() ? 10 : 1800
const B2400 = Sys.islinux() ? 11 : 2400
const B4800 = Sys.islinux() ? 12 : 4800
const B7200 = Sys.islinux() ? 13 : 7200
const B19200 = Sys.islinux() ? 14 : 19200
const B38400 = Sys.islinux() ? 15 : 38400
const B9600 = 9600
const B14400 = 14400
const B28800 = 28800
const B57600 = 57600
const B76800 = 76800
const B115200 = 115200
const B230400 = 230400
const EXTA = 19200
const EXTB = 38400

const B460800 = 460800
const B500000 = 500000
const B576000 = 576000
const B921600 = 921600
const B1000000 = 1000000
const B1152000 = 1152000
const B1500000 = 1500000
const B2000000 = 2000000
const B2500000 = 2500000
const B3000000 = 3000000
const B3500000 = 3500000
const B4000000 = 4000000

#
# Values for the QUEUE_SELECTOR argument to `tcflush'.
#
"""Discard data received but not yet read."""
const TCIFLUSH = Sys.islinux() ? 0 : 1

"""Discard data written but not yet sent."""
const TCOFLUSH = Sys.islinux() ? 1 : 2

"""Discard all pending data."""
const TCIOFLUSH = Sys.islinux() ? 2 : 3

#
# Values for the ACTION argument to `tcflow'.
#

"""Suspend output."""
const TCOOFF = Sys.islinux() ? 16 : 1

"""Restart suspended output."""
const TCOON = Sys.islinux() ? 32 : 2

"""Send a STOP character."""
const TCIOFF = Sys.islinux() ? 4 : 3

"""Send a START character."""
const TCION = Sys.islinux() ? 8 : 4

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

# The layout of a termios struct in C must be as follows
#
# ```c
# struct termios {
# 	tcflag_t	c_iflag;
# 	tcflag_t	c_oflag;
# 	tcflag_t	c_cflag;
# 	tcflag_t	c_lflag;
# 	cc_t		c_line;
# 	cc_t		c_cc[NCCS];
# 	speed_t		c_ispeed;
# 	speed_t		c_ospeed;
# };
# ```
#
# We need to create this struct in Julia and set the memory layout in Julia.
# the termios library requires passing in a termios struct that is used to get or set attributes
# There are two ways to create a `struct` in Julia that affect the memory layout.
# `struct termios end` and `mutable struct termios end`
# The first one is a immutable, and cannot be passed as a reference into a library.
# Additionally, because it is a immutable, Julia may make a copy of the struct instead of passing in the reference.
# The second one is what we want.
# Typical workflow in using the termios library involves using `tcgetattr` to initialize a termios struct, changing values appropriately
# and using `tcsetattr` to set the attributes.
# the termios struct has a field `c_cc[NCCS]`.
# This field is laid out sequentially in memory.
# In Julia, there are two ways to lay out memory sequentially. We can either
# 1) use a `NTuple{NCCS, UInt8}`
# 2) lay out the memory manually
# Tuples are immutable, which means if a user wants to change `termios.c_cc`, they would have to create a new tuple with the values needed.
# In both of these approaches, we can use `getproperty` to mimic the interface presented in C

mutable struct termios
    """Input flags"""
    c_iflag::tcflag_t
    """Output flags"""
    c_oflag::tcflag_t
    """Control flags"""
    c_cflag::tcflag_t
    """Local flags"""
    c_lflag::tcflag_t
    @static if Sys.islinux()
        c_line::cc_t
    end
    """Control chars"""
    _c_cc::NTuple{NCCS, UInt8}
    """Input speed"""
    c_ispeed::speed_t
    """Output speed"""
    c_ospeed::speed_t
end

function Base.getproperty(t::termios, name::Symbol)
    if name == :c_cc
        return _C_CC(t)
    else
        return getfield(t, name)
    end
end

function Base.propertynames(t::termios, private = false)
    return @static if Sys.islinux()
        (
            :c_iflag,
            :c_oflag,
            :c_cflag,
            :c_lflag,
            :c_cc,
            :c_ispeed,
            :c_ospeed,
        )
    else
        (
            :c_iflag,
            :c_oflag,
            :c_cflag,
            :c_lflag,
            :c_line,
            :c_cc,
            :c_ispeed,
            :c_ospeed,
        )
    end
end

struct _C_CC
    ref::termios
end

function Base.getindex(c_cc::_C_CC, index)
    return collect(c_cc.ref._c_cc)[index + 1]
end

function Base.setindex!(c_cc::_C_CC, value, index)
    _c_cc = collect(c_cc.ref._c_cc)
    _c_cc[index + 1] = value
    c_cc.ref._c_cc = NTuple{NCCS, UInt8}(_c_cc)
end

function Base.show(io::IO, ::MIME"text/plain", c_cc::_C_CC)
    X = [getfield(c_cc.ref, Symbol("_c_cc$i")) for i in 1:NCCS]
    summary(io, X)
    isempty(X) && return
    print(io, ":")
    if !haskey(io, :compact) && length(axes(X, 2)) > 1
        io = IOContext(io, :compact => true)
    end
    if get(io, :limit, false) && eltype(X) === Method
        # override usual show method for Vector{Method}: don't abbreviate long lists
        io = IOContext(io, :limit => false)
    end
    if get(io, :limit, false) && displaysize(io)[1]-4 <= 0
        return print(io, " …")
    else
        println(io)
    end
    io = IOContext(io, :typeinfo => eltype(X))
    Base.print_array(io, X)
end

function termios()

    term = @static if Sys.islinux()
        termios(
            0,
            0,
            0,
            0,
            0,
            NTuple{NCCS, UInt8}([0 for _ in 1:NCCS]),
            0,
            0,
        )
    else
        termios(
            0,
            0,
            0,
            0,
            NTuple{NCCS, UInt8}([0 for _ in 1:NCCS]),
            0,
            0,
        )
    end

   return term

end

# helper function
_file_handle(s::Base.LibuvStream) = Base._fd(s)

"""
    tcgetattr(fd::RawFD, term::termios)
    tcgetattr(s::Base.LibuvStream, term::termios)
    tcgetattr(f::Int, term::termios)

Get the tty attributes for file descriptor fd
"""
function tcgetattr(fd::RawFD, term::termios)
    r = ccall(:tcgetattr, Cint, (Cint, Ptr{Cvoid}), fd, Ref(term))
    r == -1 ? throw(TERMIOSError("tcgetattr failed: $(Base.Libc.strerror())")) : nothing
end
tcgetattr(s::Base.LibuvStream, term) = tcgetattr(_file_handle(s), term)
tcgetattr(f::Int, term) = tcgetattr(RawFD(f), term)

"""
    tcsetattr(s::Base.LibuvStream, when::Integer, term::termios)

Set the tty attributes for file descriptor fd.
The when argument determines when the attributes are changed:

- `TERMIOS.TCSANOW` to change immediately
- `TERMIOS.TCSADRAIN` to change after transmitting all queued output
- `TERMIOS.TCSAFLUSH` to change after transmitting all queued output and discarding all queued input.
"""
function tcsetattr(fd::RawFD, when::Integer, term::termios)
    r = ccall(:tcsetattr, Cint, (Cint, Cint, Ptr{Cvoid}), fd, when, Ref(term))
    r == -1 ? throw(TERMIOSError("tcsetattr failed: $(Base.Libc.strerror())")) : nothing
end
tcsetattr(s::Base.LibuvStream, when, term) = tcsetattr(_file_handle(s), when, term)
tcsetattr(f::Int, when, term) = tcsetattr(RawFD(f), when, term)

"""
    tcdrain(s::Base.LibuvStream)

Wait until all output written to file descriptor fd has been transmitted.
"""
function tcdrain(fd::RawFD)
    r = ccall(:tcdrain, Cint, (Cint, ), fd)
    r == -1 ? throw(TERMIOSError("tcdrain failed: $(Base.Libc.strerror())")) : nothing
end
tcdrain(s::Base.LibuvStream) = tcdrain(_file_handle(s))
tcdrain(f::Int) = tcdrain(RawFD(f))

"""
    tcflow(s::Base.LibuvStream, action::Integer)

Suspend transmission or reception of data on the object referred to by fd, depending on the value of action:

- `TERMIOS.TCOOFF` to suspend output,
- `TERMIOS.TCOON` to restart output
- `TERMIOS.TCIOFF` to suspend input,
- `TERMIOS.TCION` to restart input.
"""
function tcflow(fd::RawFD, action::Integer)
    r = ccall(:tcflush, Cint, (Cint, Cint), fd, action)
    r == -1 ? throw(TERMIOSError("tcflow failed: $(Base.Libc.strerror())")) : nothing
end
tcflow(s::Base.LibuvStream, action) = tcflow(_file_handle(s), action)
tcflow(fd::Int, action) = tcflow(RawFD(fd), action)

"""
    tcflush(s::Base.LibuvStream, queue::Integer)

Discard data written to the object referred to by fd but not transmitted, or data received but not read, depending on the value of queue_selector:

- `TERMIOS.TCIFLUSH` flushes data received but not read.
- `TERMIOS.TCOFLUSH` flushes data written but not transmitted.
- `TERMIOS.TCIOFLUSH` flushes both data received but not read, and data written but not transmitted.
"""
function tcflush(fd::RawFD, queue::Integer)
    r = ccall(:tcflush, Cint, (Cint, Cint), fd, queue)
    r == -1 ? throw(TERMIOSError("tcflush failed: $(Base.Libc.strerror())")) : nothing
end
tcflush(s::Base.LibuvStream, queue) = tcflush(_file_handle(s), queue)
tcflush(fd::Int, queue) = tcflush(RawFD(fd), queue)

"""
    tcsendbreak(s::Base.LibuvStream, duration::Integer)

Transmit a continuous stream of zero-valued bits for a specific duration, if the terminal is using asynchronous serial data transmission. If duration is zero, it transmits zero-valued bits for at least 0.25 seconds, and not more that 0.5 seconds. If duration is not zero, it sends zero-valued bits for some implementation-defined length of time.

If the terminal is not using asynchronous serial data transmission, tcsendbreak() returns without taking any action.
"""
function tcsendbreak(fd::RawFD, duration::Integer)
    r = ccall(:tcsendbreak, Cint, (Cint, Cint), fd, duration)
    r == -1 ? throw(TERMIOSError("tcsendbreak failed: $(Base.Libc.strerror())")) : nothing
end
tcsendbreak(s::Base.LibuvStream, duration) = tcsendbreak(_file_handle(s), duration)
tcsendbreak(f::Int, duration) = tcsendbreak(RawFD(f), duration)

"""
    cfmakeraw(term::termios)

Set the terminal to something like the "raw" mode of the old Version 7 terminal driver: input is available character by character, echoing is disabled, and all special processing of terminal input and output characters is disabled. The terminal attributes are set as follows:

    term.c_iflag &= ~(IGNBRK | BRKINT | PARMRK | ISTRIP
                 | INLCR | IGNCR | ICRNL | IXON);
    term.c_oflag &= ~OPOST;
    term.c_lflag &= ~(ECHO | ECHONL | ICANON | ISIG | IEXTEN);
    term.c_cflag &= ~(CSIZE | PARENB);
    term.c_cflag |= CS8;
"""
function cfmakeraw(term::termios)
    r = ccall(:cfmakeraw, Cint, (Ref{termios},), Ref(term))
    r == -1 ? throw(TERMIOSError("cfmakeraw failed: $(Base.Libc.strerror())")) : nothing
end

"""
    cfsetspeed(term::termios, speed::Int)

is a 4.4BSD extension. It takes the same arguments as cfsetispeed(), and sets both input and output speed.
"""
function cfsetspeed(term::termios, speed::Integer)
    r = ccall(:cfsetspeed, Cint, (Ref{termios}, speed_t), Ref(term), speed)
    r == -1 ? throw(TERMIOSError("cfsetspeed failed: $(Base.Libc.strerror())")) : nothing
end

"""
    cfgetispeed(term::termios) -> Int

Returns the input baud rate stored in the termios structure.
"""
cfgetispeed(term::termios) = ccall(:cfgetispeed, speed_t, (Ptr{termios}, ), Ref(term))

"""
    cfgetospeed(term::termios) -> Int

Returns the output baud rate stored in the termios structure.
"""
cfgetospeed(term::termios) = ccall(:cfgetospeed, speed_t, (Ptr{termios}, ), Ref(term))

"""
    cfsetispeed(term::termios, speed::Integer)

sets the input baud rate stored in the termios structure to speed, which must be one of these constants:
- `TERMIOS.B0`
- `TERMIOS.B50`
- `TERMIOS.B75`
- `TERMIOS.B110`
- `TERMIOS.B134`
- `TERMIOS.B150`
- `TERMIOS.B200`
- `TERMIOS.B300`
- `TERMIOS.B600`
- `TERMIOS.B1200`
- `TERMIOS.B1800`
- `TERMIOS.B2400`
- `TERMIOS.B4800`
- `TERMIOS.B9600`
- `TERMIOS.B19200`
- `TERMIOS.B38400`
- `TERMIOS.B57600`
- `TERMIOS.B115200`
- `TERMIOS.B230400`

The zero baud rate, B0, is used to terminate the connection. If B0 is specified, the modem control lines shall no longer be asserted. Normally, this will disconnect the line. CBAUDEX is a mask for the speeds beyond those defined in POSIX.1 (57600 and above). Thus, B57600 & CBAUDEX is nonzero.
"""
function cfsetispeed(term::termios, speed::Integer)
    r = ccall(:cfsetispeed, Cint, (Ref{termios}, speed_t), Ref(term), speed)
    r == -1 ? throw(TERMIOSError("cfsetispeed failed: $(Base.Libc.strerror())")) : nothing
end

"""
    cfsetospeed(term::termios, speed::Integer)

sets the output baud rate stored in the termios structure to speed, which must be one of these constants:
- `TERMIOS.B0`
- `TERMIOS.B50`
- `TERMIOS.B75`
- `TERMIOS.B110`
- `TERMIOS.B134`
- `TERMIOS.B150`
- `TERMIOS.B200`
- `TERMIOS.B300`
- `TERMIOS.B600`
- `TERMIOS.B1200`
- `TERMIOS.B1800`
- `TERMIOS.B2400`
- `TERMIOS.B4800`
- `TERMIOS.B9600`
- `TERMIOS.B19200`
- `TERMIOS.B38400`
- `TERMIOS.B57600`
- `TERMIOS.B115200`
- `TERMIOS.B230400`

The zero baud rate, B0, is used to terminate the connection. If B0 is specified, the modem control lines shall no longer be asserted. Normally, this will disconnect the line. CBAUDEX is a mask for the speeds beyond those defined in POSIX.1 (57600 and above). Thus, B57600 & CBAUDEX is nonzero.
"""
function cfsetospeed(term::termios, speed::Integer)
    r = ccall(:cfsetospeed, Cint, (Ref{termios}, speed_t), Ref(term), speed)
    r == -1 ? throw(TERMIOSError("cfsetospeed failed: $(Base.Libc.strerror())")) : nothing
end

end # module


text:
using StackOverflow

function convert_HTTP_Response_To_JSON(response)
    compressed = HTTP.payload(response);

    decompressed = transcode(GzipDecompressor, compressed);

    json = JSON.parse(IOBuffer(decompressed))

    return json
end


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# ╔═╡ 023da279-7a42-4cb1-addb-613d1363a282
md"""
# 2021 project ideas:
Below you can find a list of possible projects. Almost all of these are described in a very generic way and will require some research on your behalf. Some background and/or similar studies are provided for each project, most of these being more complex than what is expected of you.

**Aim of the project**:
1. Describe the setting or problem
2. Describe the desired outcome of the study
3. Describe the model you built, including its limitations and accepted hypotheses.
4. Describe your analysis
5. Describe the outcome or suggestions

Note that for a project with a deterministic outcome, the use of simulation might not be required.
"""

# ╔═╡ 42968b6f-a86f-48dd-b354-d06dfd53624c
md"""
## Swarming (multiple projects)
The use of drone swarms has several applications in a military context. Drones can be equipped with different payloads such as high explosives, sensors (ISR), jammers, etc. Currently, all major military countries are conducting research on both the use of and the defense against drone swarms.

### Possible topics:
- **command and control (C2) of a drone swarm**: for the technology to be scalable, each element of the swarm should be cheap, so you will have a limited amount of control drones that are used for relaying communications with the HQ. The movement patterns of drones in different cases can be studied and optimised.
- **use of a drone swarm for ISR purposes**: drone swarms can be used to [gather intelligence](https://www.timesofisrael.com/in-apparent-world-first-idf-deployed-drone-swarms-in-gaza-fighting/) in addition to more traditional assets. Given a specific scenario and an ISR capacity, determine optimal use cases.
- **defense against a drone swarm**: given an adversarial drone swarm attack, determine optimal techniques of defense (weapon system and/or engagement method).
- **attack using a drone swarm**:  given an adversial defense system, determine optimal techniques/patterns/tactics to maximise the impact of the attack.

### References:
* [General info](https://www.vifindia.org/article/2020/april/27/drone-swarms-bracing-up-with-the-new-threat)
* [More details](http://unsworks.unsw.edu.au/fapi/datastream/unsworks:45320/SOURCE02?view=true)
"""

# ╔═╡ d3195ce3-cc34-458f-a526-f346591b34b2
md"""
## Mine hunting
Mine Counter Measure (MCM) is a general capability that can include several types of missions required to complete a specific operational intent:
- General or detailed survey to prepare Mine Warfare (MW) operations,
- Surveillance or exploration to assess Mine Warfare risk,
- Reconnaissance and Mine avoidance to avoid mine threat,
- Clearance to sanitize,
- Masking and jamming to prevent mine action,
- Provocation and discredit to deceive mine logic.

### Possible projects
- **optimise sweep patterns to maximise detection**: the detection probability of a mine depends on its position with respect to the sonar and the composition of the seabed. Simulate a minefield and try to optimise the sweeping settings to finetune the detection rate.
- **optimise mine field**: different types of (deep) sea mines exist. Given a set of mines, evaluate the composition/layout of a minefield for maximum effect.

### References
- [Synthetic Aperture and 3D Imaging for Mine Hunting Sonar](https://hal.archives-ouvertes.fr/hal-00504862/document)
- [Minhunting sonar perfomance](https://apps.dtic.mil/sti/pdfs/ADA342568.pdf)
- [Modern mines](https://cimsec.org/modern-naval-mines-not-your-grandfathers-weapons-that-wait/)
"""

# ╔═╡ 8bc9e310-3151-4c49-b009-04b4cc1a4a84
md"""
## Airfield destruction
Destroying or rendering unusable an enemy airfield can have a major impact on the course of a conflict (think of the Six Day War). Typically, ballistic missiles and cruise missiles are used, each with a different effect and a different cost. A military airfield usually has a defense system against missiles and rockets where the interception rate is different for each weapon system.
### Possible project
Consider an airport with its different components: the runway, the critical infrastructure (POL + control tower) and the different airplanes. You could look at Kleine Brogel or Florennes air base for inspiration. Determine the effect of using a specific set of weapon systems (ballistic & guided munitions, with submunitions or not with different CEP) and consider the optimal/ideal/minimal configuration that allows you to take out the airfield.

### References
* [Defense systems](https://www.boeing.com/defense/missile-defense/)
"""

# ╔═╡ 279eda92-e1cf-480e-a093-e17625046125
md"""
## Digital communication

### Possible projects
* **TCP/IP**: simulate the TCP/IP model for data transfer and look up the limits with respect to data troughput for given hardware.
* **routing**: analyse how well routing works using different routing protocols on a (dynamic) network.
* **tactical radio network**: compare different scheduling algorithms to maximise data throughput in a tactical radio network.

### References
* [Basic on internet and IP routing, TCP/IP](https://www.khanacademy.org/computing/computers-and-internet/xcae6f4a7ff015e7d:the-internet)
* [Routing protocol](https://en.wikipedia.org/wiki/Routing_protocol)
* [Scheduling](https://en.wikipedia.org/wiki/Scheduling_(computing)#Scheduling_disciplines)
"""

# ╔═╡ cb784ac1-7dde-4163-a325-379f87ccb18b
md"""
## Fine-tuning a grading system
Some universities allow a student to fail a course and still be able to go to the next year or be exempt from retaking the exam. A mishap in an otherwise good course can be forgiven as a result. The "optimal" ground rules of such a system are far from trivial. After all, there are conflicting interests to be realized. One wishes, on the one hand, to limit the number of students who can abuse the system and, on the other hand, to ensure that the majority of students can benefit from the system. In addition, for practical reasons, one may wish to limit the overall number of resits.

### Possible project
By using historical data over several years, it is possible to establish a multinomial distribution that describes the results of students in an academic bachelor. You can use this input to adjust the parameters of your system, taking into account a conflicting desired results (e.g. $<1\%$ “freeloaders”, limiting number of resits, maximal fair pardons).

### References
"""

# ╔═╡ a7659bd3-796b-4155-b01f-1d9ef6b42f12
md"""
## Chemistry
In the various chemistry courses, you have studied the reactions between molecules by means of kinetics. A reaction was modelled using a differential equation. However, it is also possible to model kinetics using stochastic methods.

### Possible projects
Compare a stochastic model  (there are multiple in the reference) with a numerical method that you already know for some known reactions. In addition, use an optimisation method to fine-tune your reaction.

### References
* [Stochastic chemical kinetics](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5425731/)
"""

# ╔═╡ 57e9e568-4bb7-43d1-91cf-0711cc2fa905
md"""
## Finance
Modeling and simulation is widely used in the financial sector. Different techniques can be used.
### Possible projects
* **Portfolio optimisation of crypto currencies in combination with discrete even simulation**: We have already briefly discussed portfolio optimisation in the practical applications. However, the crypto market is much more volatile, so here it may be worthwhile to combine these optimisations with discrete event simulation and real data.
* **Insurance solvency**: Use Monte Carlo simulation or another suited method for a specific insurance to find an optimal balance between insurance premium and solvency.
* **Global trade embargo effects:** In today's world, nations are less likely to declare war on another country, but first try to achieve their goals through diplomatic or economic means. The global trade network has been extensively studied and data is freely available. Study the effect of a trade embargo on a particular country and consider possible unexpected alliances (where you can consider trade as an optimisation problem).
### References
* [Solvency ratio](https://en.wikipedia.org/wiki/Solvency_ratio)
* [Monte Carlo in insurance](https://www.researchgate.net/publication/332637928_Monte_Carlo_Methods_for_Insurance_Risk_Computation)
* [World Trade Network](https://www.researchgate.net/publication/23786718_The_World_Trade_Network)
* [World Trade Network dataset](http://www.cepii.fr/cepii/en/bdd_modele/presentation.asp?id=27), requires free registration, holds all network data.
"""

# ╔═╡ 2e4860bd-b74b-4671-8e57-c6849bc63d88
md"""
## Logistics
Supply chain optimisation and the optimal use of personnel or resources can significantly increase the efficiency of a logistics depot. These problems can be solved with linear programming, but it is often difficult to reach a solution because these problems are typically NP-hard. As a solution, either heuristics (cf. next year) or discrete event simulation is used.

### Possible projects
- **Personnel tasking**: given a specific warehouse, increase efficiency by optimising the tasks each person is executing.
- **Warehouse layout**: given a set of products and historical order data, optimise warehouse or grocery store layout for highest possible output (e.g. something similar to collect and go from the perspective of the store personnel)
- **Seaport operations**: In a cargo port, there are a multitude of actors at work such as pilots, ships, cranes, trains, trucks, etc.  In the context of efficient use of resources, it is interesting to simulate this in order to identify possible bottlenecks in different scenarios.  If necessary, study a part of a Belgian port.
### References
- [complete project personnel optimisation](https://www.mdpi.com/1999-4893/13/12/326)
- [data source](https://www.kaggle.com/msp48731/frequent-itemsets-and-association-rules/data)
- [seaport operations](https://sci-hub.st/https://doi.org/10.1177%2F003754979807000401)
- [layout optimisation](https://sci-hub.st/https://doi.org/10.1111/itor.12852)
"""

# ╔═╡ 4357e4be-1005-470f-978b-90844dd31b1a
md"""
## Healthcare
Within the healthcare sector there are a lot of use cases for discrete event simulation.

### Possible projects
get creative.
### References
- [overview paper](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7491890/)
- [staffing a neonatal ICU](https://journals.sagepub.com/doi/pdf/10.1177/1460458216628314)
- [capacity planning (COVID)](https://arxiv.org/abs/2012.07188)
"""

# ╔═╡ 88b57987-b284-41a7-8c35-3b29a0436d0f
md"""
## Biology
Simulation can also be used in biology.

### References:
- [evolution of a virus](https://academic.oup.com/ve/article/5/1/vez003/5372481)
- [evolution of a virus bis](https://sci-hub.st/https://doi.org/10.1007/s11538-018-00550-4)
"""

# ╔═╡ 72df832c-442f-4693-9953-df873951a1db
md"""
## Energy grid management
Managing the electricity grid is a major challenge.
### Possible projects

### References
- [Introduction to electricity markets](https://www.e-education.psu.edu/ebf483/node/816)
"""

# ╔═╡ 65c30d33-de0e-410d-afa9-9e0e84431012
md"""input ben:
- risicomodellen gecombineerd met discreet event + bayesiaans (covid + splitsen)
-
"""

# ╔═╡ 79db607a-8af7-4f60-a339-0e27cd6edf1f
md"""
## General topic

### Possible projects

### References
"""

# ╔═╡ 3117a9ed-4577-47fb-bc7f-54fb9572e313

# ╔═╡ Cell order:
# ╟─023da279-7a42-4cb1-addb-613d1363a282
# ╟─42968b6f-a86f-48dd-b354-d06dfd53624c
# ╟─d3195ce3-cc34-458f-a526-f346591b34b2
# ╟─8bc9e310-3151-4c49-b009-04b4cc1a4a84
# ╟─279eda92-e1cf-480e-a093-e17625046125
# ╟─cb784ac1-7dde-4163-a325-379f87ccb18b
# ╟─a7659bd3-796b-4155-b01f-1d9ef6b42f12
# ╟─57e9e568-4bb7-43d1-91cf-0711cc2fa905
# ╟─2e4860bd-b74b-4671-8e57-c6849bc63d88
# ╟─4357e4be-1005-470f-978b-90844dd31b1a
# ╟─88b57987-b284-41a7-8c35-3b29a0436d0f
# ╟─72df832c-442f-4693-9953-df873951a1db
# ╠═65c30d33-de0e-410d-afa9-9e0e84431012
# ╠═79db607a-8af7-4f60-a339-0e27cd6edf1f
# ╠═3117a9ed-4577-47fb-bc7f-54fb9572e313

max_stars_repo_path: Exercises/PS09 - Projects.jl

max_issues_repo_path: Exercises/PS09 - Projects.jl

max_forks_repo_path: Exercises/PS09 - Projects.jl


text:
# paste code in the REPL:
using WebSockets, JSON
url = "wss://www.bitmex.com/realtime?subscribe=trade:XBT"

function open_websocket()
  WebSockets.open(url) do ws
    while isopen(ws)                   # 1
      data, success = readguarded(ws)
      !success && break
      data = JSON.parse(String(data))  # 2
      print(data, "\n")
    end

    if !isopen(ws)                     # 3
      @async open_websocket()
    end

  end
end

@async open_websocket()
# Task (runnable) @0x000000000c6eeca0

# Sample output:
# Dict{String, Any}("docs" => "https://ws.bitmex.com/app/wsAPI", "info" => "Welcome to the BitMEX Realtime API.", "version" => "2022-04-29T04:21:29.000Z", "timestamp" => "2022-05-01T11:24:31.350Z", "limit" => Dict{String, Any}("remaining" => 179))
# Dict{String, Any}("success" => true, "request" => Dict{String, Any}("args" => "trade:XBT", "op" => "subscribe"), "subscribe" => "trade:XBTUSD")
# Dict{String, Any}("types" => Dict{String, Any}("grossValue" => "long", "homeNotional" => "float", "price" => "float", "tickDirection" => "symbol", "side" => "symbol", "trdMatchID" => "guid", "foreignNotional" => "float", "symbol" => "symbol", "size" => "long", "timestamp" => "timestamp"), "action" => "partial", "keys" => Any[], "attributes" => Dict{String, Any}("symbol" => "grouped", "timestamp" => "sorted"), "filter" => Dict{String, Any}("symbol" => "XBTUSD"), "data" => Any[Dict{String, Any}("grossValue" => 3682994, "homeNotional" => 0.03682994, "price" => 38012.5, "tickDirection" => "ZeroMinusTick", "side" => "Sell", "trdMatchID" => "a0695101-e4fb-8b27-422c-c168cf25779d", "foreignNotional" => 1400, "symbol" => "XBTUSD", "size" => 1400, "timestamp" => "2022-05-01T11:23:45.185Z")], "table" => "trade", "foreignKeys" => Dict{String, Any}("side" => "side", "symbol" => "instrument"))
# Dict{String, Any}("action" => "insert", "data" => Any[Dict{String, Any}("grossValue" => 263075, "homeNotional" => 0.00263075, "price" => 38012, "tickDirection" => "MinusTick", "side" => "Sell", "trdMatchID"


text:
using CausalInference
using Graphs
using Test
using Random
using Graphs
using Distributions
using CausalInference: disjoint_sorted

@test disjoint_sorted([],[1,2])
@test disjoint_sorted([1,2],[])
@test disjoint_sorted([],[])
@test disjoint_sorted([1,2],[3,4])
@test disjoint_sorted([1,2],[3])
@test disjoint_sorted([1],[3,4])
@test !disjoint_sorted([1,3],[2,3,4])
@test !disjoint_sorted([1,2,3],[3,4,5])

# tests modelled after examples discussed in chapter 2 of
# Pearl, Judea. Causality. Cambridge University Press, 2009.

g = DiGraph(5)
d = nv(g)
for (i,j) in [(1,2), (1,3), (2,4), (3,4), (4,5)]
   add_edge!(g, i, j)
end

h, s = skeleton(d, dseporacle, g)
Z = orientable_unshielded(h, s)
@testset "unshielded" begin
    @test Graph(g) == h
    for z in Z
        u, v, w = z
        n = inneighbors(h, v)
        @test u in n
        @test w in n
        @test !(u in neighbors(h, w))
    end
end

dg = pcalg(d, dseporacle, g)

Random.seed!(123)
N = 1000
p = 0.01
x = randn(N)
v = x + randn(N)*0.25
w = x + randn(N)*0.25
z = v + w + randn(N)*0.25
s = z + randn(N)*0.25

X = [x v w z s]
C = cor(X)
df = (x=x, v=v, w=w, z=z, s=s)

println("Running Gaussian tests")
@time gaussci_g = pcalg(df, p, gausscitest)

println("Running CMI tests")
@time cmi_g = pcalg(df, 0.1, cmitest)

@testset "pcalg_edgde_test" begin
    @test collect(Graphs.edges(cmi_g)) == collect(Graphs.edges(dg))
    @test collect(Graphs.edges(gaussci_g)) == collect(Graphs.edges(dg))
end


text:
#TODO handle identity separately in pointwise operators
abstract type AbstractOperator end

abstract type  LinearOperator <: AbstractOperator end

struct PointwiseOperator <: AbstractOperator
    domain::AbstractDomain
    operations::Array
    gradients::Array
end

struct OperatorChain <: AbstractOperator
    domain::Union{AbstractDomain, Tuple{Vararg{AbstractDomain}}}
    target::Union{AbstractDomain, Tuple{Vararg{AbstractDomain}}}
    operators::Array{O} where O <: AbstractOperator
end

####################################################################################################
####################################################################################################
#Property interfaces
domain(op::AbstractOperator) = op.domain
target(op::AbstractOperator) = op.target

target(op::PointwiseOperator) = op.domain

get_operators(op::AbstractOperator) = identity(op)
get_operators(op::OperatorChain) = op.operators

####################################################################################################
####################################################################################################
#combine_operators Interface
function (op_later::AbstractOperator)(op_first::AbstractOperator)
    return combine_operators(op_later, op_first)
end

function combine_operators(op_later::PointwiseOperator, op_first::PointwiseOperator)
    @assert domain(op_later) == domain(op_first)
    operations = vcat(op_first.operations, op_later.operations)
    gradients = vcat(op_first.gradients, op_later.gradients)
    return PointwiseOperator(op_first.domain, operations, gradients)
end

function combine_operators(op_later, op_first)
    @assert target(op_first) == domain(op_later)
    operators = vcat(get_operators(op_first), get_operators(op_later))
    return OperatorChain(domain(op_first), target(op_later), operators)
end

####################################################################################################
####################################################################################################
#apply Interface

function (op::AbstractOperator)(f::Union{Field, MultiField})
    @assert f.domain == domain(op)
    return apply(op, f.val)
end

#function (op::AbstractOperator)(vals...)
#    return apply(op, vals)
#end

function apply(op::AbstractOperator, vals::NamedTuple)
    vals = [vals[key] for key in keys(domain(op))]
    return apply(op, vals)
end

#TODO generalize to NTuple for val
function apply(op::PointwiseOperator, val::AbstractArray)
    res = @avx @. op.operations[1](val)
    for operation in op.operations[2:end]
        @avx @. res = operation(res)
    end
    return res
end

function apply(op::OperatorChain, val)
    for operator in op.operators
        val = apply(operator, val)
    end
    return val
end

####################################################################################################
####################################################################################################
#apply linearization Interface

#function (op::LinearOperator)(val, jac)
#    res = op(val)
#    jacobian(x) = op(jac(x))
#    return res, jacobian
#end
#
#function (op::PointwiseOperator)(val, jac)
#    gradient = ones(size(val))
#    res = copy(val)
#    for (operation, grad) in zip(op.operations, op.gradients)
#        @avx @. gradient *= grad(res)
#        @avx @. res = operation(res)
#    end
#    jacobian(x) = @avx gradient .* jac(x)
#    return res, jacobian
#end
#
#function (op::OperatorChain)(val, jac)
#    res = copy(val)
#    for operator in op.operators
#        #TODO here, res can be modified by apply, so pointwise operators don't need to make a copy
#        res, jac = operator(res, jac)
#    end
#    return res, jac
#end

####################################################################################################
####################################################################################################
#miscellaneous
adjoint(op::LinearOperator) = LinearOperator(target(op), domain(op), op.adjoints, op.operations)


text:
using GaussianProcesses
using Plots

include("barrier_btg.jl")
include("../../computation/tdist_mle.jl")
include("../../computation/opt_tdist_mle.jl")

btg0 = load_synthetic_btg()

println("lambda nodes")
display(sort(collect((keys(btg0.λbuffer_dict))))')

randseed = getseed()
Random.seed!(randseed);

#btg0 = btg(trainingData0, rangeθ, rangeλ; quadtype = myquadtype)
(pdf0_raw, cdf0_raw, dpdf0_raw, quantInfo0_raw) = solve(btg0);

######
###### Handy test code snippets
######

#v(x) = cdf0_raw(reshape([0.05], 1, 1), reshape([1.0], 1, 1), x)

f = theta -> tdist_mle(btg0, theta, -0.31175)
g(theta, lambda) = tdist_mle(btg0, theta, lambda)
g_vec(arr) = tdist_mle(btg0, abs(arr[1]), arr[2])
h = lambda -> tdist_mle(btg0, 140.0, lambda)
h_boxed = (lambda_boxed) -> tdist_mle(btg0, 140.0, lambda_boxed[1])

MLE_ESTIMATION = false
if MLE_ESTIMATION
    lambda_star = optimize_lambda(h_boxed, [-3, 3])
    println("lambda_star")
    display(lambda_star)

    opt_vec = optimize_theta_lambda_single(g_vec, [1.0, 1000.0], [-5.0, 5.0])
    println("optimal theta-lambda pair")
    display(opt_vec)
end

#barrier(x) = 0.05/(0.001+(x-1)^2)  + 0.05/(0.001+(x-2)^2) - log(1-2(x/2-2)) + 1 + sin(x+exp(x))/10 + exp(-(x-2)^2)
#Plots.plot()
#plt(barrier, 0, 4.999)

if true #regular GP Test
    println("Start Vanilla GP Test")
    x = x_train #reshape(x_train, length(x_train) ,1)
    y = y_train #reshape(y_train, length(y_train), 1)
    # build and fit a GP
    #mymean = MeanLin(zeros(d))
    mymean = MeanZero()
    # mymean = MeanZero()
    kern = SE(0.0, 0.0)
    gp = GP(dropdims(x, dims=2), y, mymean, kern)
    GaussianProcesses.optimize!(gp)
    # predict
    μ, σ² = predict_y(gp, dropdims(x_test, dims=2)); stdv = sqrt.(σ²)
    error_GP = abs.(μ .- y_test)
    error_abs_GP = mean(error_GP)
    error_sq_GP = mean(error_GP.^2)
    CI_test_GP = hcat(-1.96.*stdv .+ μ, 1.96.*stdv .+ μ)
    count_test_GP = sum((y_test .>= CI_test_GP[:, 1]) .* (y_test .<= CI_test_GP[:,2]))/n_test
    nlpd_GP = -mean(log.(pdf.(Normal(), (y_test.-μ)./stdv)./stdv))
    println("MSE:")
    display(error_sq_GP)
    println("Conf interval accuracy:")
    display(count_test_GP)
end
####################################
############### Test ###############
####################################

#TODO:
#scatter training and testing points
#fine mesh for true mean and btg mean
UB = parsed_args["upper_bound_estimate"]
if true
    @info "Start Test"
    before = Dates.now()
    count_test = 0
    error_abs = 0.
    error_sq = 0.
    nlpd = 0.
    id_fail = []
    id_nonproper = []
    CI_set = zeros(n_test, 2)
    median_set = zeros(n_test)
    #summary statistics loop
    for i in 1:n_test
        global error_abs, error_sq, nlpd, count_test
        # mod(i, 20) == 0 ? (@info i) : nothing
        # @info "i" i
        x_test_i = reshape(x_test[i, :], 1, d)
        Fx_test_i = reshape(Fx_test[i, :], 1, p)
        try
            y_test_i = y_test[i]

            median_test_i = fzero(x -> cdf0_raw(x_test_i, reshape([1.0], 1, 1), x) - 0.5, 0, UB)
            median_set[i] = median_test_i
            try #compute confidence interval bounds
                aa = fzero(x -> cdf0_raw(x_test_i, reshape([1.0], 1, 1), x) - 0.975, 0, UB)
                bb = fzero(x -> cdf0_raw(x_test_i, reshape([1.0], 1, 1), x) - 0.025, 0, UB)
                CI_test_i = [bb, aa]
                count_test += (y_test_i >= CI_test_i[1])&&(y_test_i <= CI_test_i[2]) ? 1 : 0
                CI_set[i, :] = CI_test_i
            catch err
                append!(id_fail, i)
            end
            error_abs += abs(y_test_i - median_test_i)
            error_sq += (y_test_i - median_test_i)^2
            nlpd += log(pdf0_raw(y_test_i))
        catch err
            append!(id_nonproper, i)
        end
    end

        #compute summary statitics
        count_test /= n_test
        error_abs  /= n_test
        error_sq   /= n_test
        nlpd       /= -n_test
        after = Dates.now()
        elapsedmin = round(((after - before) / Millisecond(1000))/60, digits=5)
        println("error_sq")
        display(error_sq)
        println("count_test")
        display(count_test)
        println("elapsedmin")
        display(elapsedmin)

        #plotting loop for high res btg mean/conf intervals and true mean
        N_plot = parsed_args["N_plot"]
        plot_CI_set = zeros(N_plot, 2)
        plot_median_set = zeros(N_plot)
        mesh_x = range(xs[1], stop = xs[end], length = N_plot)
        mesh_y = fun.(mesh_x)
        count_test_plot = 0.0
        for i in 1:N_plot
            global error_abs, error_sq, nlpd, count_test_plot
            # mod(i, 20) == 0 ? (@info i) : nothing
            # @info "i" i
            x_test_i = reshape([mesh_x[i]], 1, d)
            Fx_test_i = covariate_fun(x_test_i, 1) #1 indicates constant mean
            try
                y_test_i = fun(x_test_i[1])
                median_test_i = fzero(x -> cdf0_raw(x_test_i, reshape([1.0], 1, 1), x) - 0.5, 0, UB)
                plot_median_set[i] = median_test_i
                try #compute confidence interval bounds
                    aa = fzero(x -> cdf0_raw(x_test_i, reshape([1.0], 1, 1), x) - 0.975, 0, UB)
                    bb = fzero(x -> cdf0_raw(x_test_i, reshape([1.0], 1, 1), x) - 0.025, 0, UB)
                    CI_test_i = [bb, aa]
                    count_test_plot += (y_test_i >= CI_test_i[1])&&(y_test_i <= CI_test_i[2]) ? 1 : 0
                    plot_CI_set[i, :] = CI_test_i
                catch err
                    append!(id_fail, i)
                end
                error_abs += abs(y_test_i - median_test_i)
                error_sq += (y_test_i - median_test_i)^2
                nlpd += log(pdf0_raw(y_test_i))
            catch err
                append!(id_nonproper, i)
            end
        end
        # Plot true mean and btg mean
        PyPlot.close("all") #close existing windows
        if parsed_args["log_scale"]
            PyPlot.yscale("log")
        end
        PyPlot.plot(mesh_x, plot_median_set, label = "BTG median")
        PyPlot.plot(mesh_x, mesh_y, label = "true")
        PyPlot.fill_between(mesh_x, plot_CI_set[:, 1], plot_CI_set[:, 2], alpha = 0.3, label = "95% confidence interval")
        PyPlot.scatter(x_train, y_train, s = 30, color = "blue", marker = "*")
        PyPlot.scatter(x_test, y_test, s = 30, color = "red", marker = "o")
        PyPlot.legend(fontsize=8)
        PyPlot.grid()
        PyPlot.title("BTG $myquadtype", fontsize=10)
        #PyPlot.savefig("exp_synthetic_mle_btg4.pdf")
    end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using LispLikeEval
using Test

@testset "LispLikeEval.jl" begin
    # Write your tests here.
end