text:
import Base: *, sort
using Random, Statistics, DataFrames, CSV
using BioSequences, FASTX, XLSX, CPUTime
import BioSequences: iscompatible

function *(seq::LongDNASeq, base::DNA)
    seq * LongDNASeq([base])
end

# Overload sort to order DNA codes by degeneracy:
#   A,G,C,T,M,R,W,S,Y,K,V,H,D,B,N
function sort(seq::LongSequence{DNAAlphabet{4}}; rev::Bool=false)
    sort(convert(Vector{DNA}, seq), by=degeneracy, rev=rev)
end

# Determine whether or not a candidate sequence belongs to any member of the list of sites. We only want seq that belong to sites
"""
    isvalid(seq, sites)

Check if any of the sites appear in the sequence.
"""
function isvalid(seq, sites)
    for site in sites
        if occursin(site, seq)
            return false
        end
    end
    return true
end

"""
    find_valid_randomer(prefix, bases, sites)

Randomly select bases to add the prefix, ensuring no site
appears in the randomer.
"""
function find_valid_randomer(prefix, bases, sites)
    n = length(bases)
    pre_n = length(prefix)
    randomer = prefix * (dna"-" ^ n)
    for i = 1:n
        candidates = shuffle(bases[i])
        for cand in candidates
            randomer[pre_n+i] = cand
            if isvalid(randomer[1:pre_n+i], sites)
                break
            else
                randomer[pre_n+i] = DNA_Gap
            end
        end
        if randomer[pre_n+i] == DNA_Gap
            break
        end
    end
    return randomer
end

"""
    degeneracy(base::DNA[, uselog2=true])

Calculate the degeneracy of a single base. The degeneracy
is the number of bases in the code, so
    degeneracy(DNA_N) == 4
and
    degeneracy(DNA_A) == 1

If `uselog2`, the degeneracy is log2 transformed. Because degeneracy(DNA_Gap)=0, we approximate
the log2 transformation as 2^-100.
"""
function degeneracy(base::DNA; uselog2=true)
    deg = 0.0
    if base == DNA_N
        deg = 4.0
    elseif base == DNA_A || base == DNA_C || base == DNA_G || base == DNA_T
        deg = 1.0
    elseif base == DNA_B || base == DNA_D || base == DNA_H || base == DNA_V
        deg = 3.0
    elseif (base == DNA_K || base == DNA_Y || base == DNA_S || base == DNA_W
            || base == DNA_R || base == DNA_M)
        deg = 2.0
    else # DNA_Gap
        deg = 0.0
    end

    if uselog2
        if base == DNA_Gap
            deg = 2^-100  #log2(0) is negative infinity, so we use something close to zero
        end

        return log2(deg)

    else
        return deg
    end
end

"""
    degeneracy(seq::LongSequence{DNAAlphabet{4}}; uselog2=true)

Calculate the degeneracy of a DNA sequence.
"""
function degeneracy(seq::LongSequence{DNAAlphabet{4}}; uselog2=true)
    if uselog2
        return sum(degeneracy.(seq, uselog2=true))
    else
        return prod(degeneracy.(seq, uselog2=false))
    end
end

"""
    pooled_degeneracy(pool)

Calcualte the sum degeneracy of pool.
"""
function pooled_degeneracy(pool)
    deg=0;
    for randomer in pool
        deg +=degeneracy(randomer;uselog2=false);
    end
    return deg
end

"""
    iscompatible(seq_A::LongSequence{DNAAlphabet{4}},seq_B::LongSequence{DNAAlphabet{4}})

Check if seq_A and seq_B contain at least one common sequence. Seq_A and Seq_B must be the same length.
"""
function iscompatible(seq_A::LongSequence{DNAAlphabet{4}},seq_B::LongSequence{DNAAlphabet{4}})
    if length(seq_A) == length(seq_B)
        return prod(iscompatible.(seq_A,seq_B))
    else # Sequences are different lenghts
        return false
    end
end

###############
#OligoRL algorithm
###############

# Run Monte Carlo simulations using a random base policy over a fixed horizon.
function simulate_random(prefix, bases, sites; nsims=1000, horizon=length(bases), reward=degeneracy,kwargs...)
    sim_reward = zeros(nsims)
    true_horizon = min(horizon, length(bases)) # don't go beyond the randomer length
    for i = 1:nsims
        sim_reward[i] = reward(find_valid_randomer(prefix, bases[1:true_horizon], sites))
    end
    return mean(sim_reward)
end

# Run a 1-step greedy lookahead policy.
function simulate_greedy(prefix, bases, sites; horizon=length(bases), reward=degeneracy,kwargs...)
    true_horizon = min(horizon, length(bases)) # don't go beyond the randomer length
    for h = 1:true_horizon
        # Iterate through the bases by decending degeneracy, so the greedy approach
        # is to stop as soon as we've found one.
        found = false # define out here for scoping rules
        for base in sort(bases[h], rev=true)
            found = false
            if isvalid(prefix * base, sites)
                prefix *= base
                found = true
                break
            end
        end
        if ~found
            # None of the bases were valid, so stop and leave the DNA_Gap
            break
        end
    end
    return reward(prefix)
end

function rollout(bases::Array{LongSequence{DNAAlphabet{4}},1}, sites::Array{LongSequence{DNAAlphabet{4}},1}; simulate=simulate_random, kwargs...)
    n = length(bases)
    randomer = dna"-" ^ n

    # If the unwanted site is non-palindromic, we also need to block
    # the reverse complement.
    for i = 1:length(sites)
        if ~ispalindromic(sites[i])
            push!(sites, reverse_complement(sites[i]))
        end
    end

    # Find the longest unwanted site and use this as the horizon.
    max_len = map(length, sites) |> maximum
    horizon = max_len

    for i = 1:n
        best_base = DNA_Gap
        best_score = -1.0
        start = max(i - max_len, 1) # only need to look back max_len when checking the oligos
        for base in bases[i]
            if ~isvalid(randomer[start:i-1] * base, sites)
                # adding this bases creates a restriction site; skip it
                continue
            end
            score = simulate(randomer[start:i-1] * base, bases[i+1:end], sites; horizon=horizon, kwargs...)

            # Update use this base if:
            #   1. the mean randomer score is higher than the previous best
            if score > best_score
                best_score = score
                best_base = base
            end
        end

        if best_base == DNA_Gap
            # The sequence has terminated; there are no bases that can be added without making
            # a restriction site.
            break
        else
            randomer[i] = best_base
        end
    end

    return randomer
end

randomerlen=10
all_ns=dna"AGCTMRWSYKVHDBN"
all_bases=[ all_ns for i =1:randomerlen]
sites=[dna"------"];

randomer=rollout(all_bases,sites;simulate=simulate_random,nsims=100,reward=degeneracy)


text:
using Alpaca
using RDatasets, Test, Distributions
using Random

rng = MersenneTwister(1234)

df = dataset("datasets", "iris")
df.binary = vec(Float64.(rand(rng,0:1,size(df,1),1)))
df[!,:SpeciesDummy] = categorical(df[!,:Species])
idx = rand(rng,1:3,size(df,1),1)
a = ["A","B","C"]
df.Random = vec([a[i] for i in idx])
df.RandomCategorical = categorical(df[!,:Random])

# proper interface

# try it natively in R
# @rput df
# R"mod <- feglm(binary ~ SepalWidth | SpeciesDummy , df, family = binomial())"
# R"summary(mod, \"hessian\")"
# binomial/logit
result = Alpaca.feglm(df, Alpaca.@formula(binary ~ SepalWidth),
    Binomial(),
    fe = :SpeciesDummy
    )
@test StatsBase.coef(result) ≈ [-0.221486] atol = 1e-4
# testing display of RegressionResult
@show result
# with clustering
# R"mod <- feglm(binary ~ SepalWidth | SpeciesDummy | SpeciesDummy, df, family = binomial())"
# R"summary(mod, \"cluster\", cluster = ~ SpeciesDummy)"
result = Alpaca.feglm(df, Alpaca.@formula(binary ~ SepalWidth),
    Binomial(),
    fe = :SpeciesDummy, vcov = :(cluster(SpeciesDummy))
    )
@test StatsBase.stderror(result) ≈ [0.317681] atol = 1e-4
# two-way clustering
# R"mod <- feglm(binary ~ SepalWidth | SpeciesDummy | SpeciesDummy + RandomCategorical, df, family = binomial())"
# R"summary(mod, \"cluster\", cluster = ~ SpeciesDummy + RandomCategorical)"
result = Alpaca.feglm(df, Alpaca.@formula(binary ~ SepalWidth),
    Binomial(),
    fe = :SpeciesDummy, vcov = :(cluster(SpeciesDummy + RandomCategorical))
    )
@test StatsBase.stderror(result) ≈ [0.228229] atol = 1e-4
# poisson
# R"mod <- feglm(binary ~ SepalWidth | SpeciesDummy, df, family = poisson())"
# R"summary(mod, \"hessian\")"
result = Alpaca.feglm(df, Alpaca.@formula(binary ~ SepalWidth),
    Poisson(),
    fe = :SpeciesDummy
    )
@test StatsBase.coef(result) ≈ [-0.103008] atol = 1e-4
# two fe's
# R"mod <- feglm(binary ~ SepalWidth | SpeciesDummy + RandomCategorical, df, family = poisson())"
# R"summary(mod, \"hessian\")"
result = Alpaca.feglm(df, Alpaca.@formula(binary ~ SepalWidth),
    Poisson(),
    fe =:(SpeciesDummy + RandomCategorical)
    )
@test StatsBase.coef(result) ≈ [-0.153311] atol = 1e-4

# rng = MersenneTwister(1234)
# N = 1_000_000
# K = 100
# id1 = rand(rng, 1:(round(Int64,N/K)), N)
# id2 = rand(rng, 1:K, N)
# x1 =  randn(rng, N) ./ 10.0
# x2 =  randn(rng, N) ./ 10.0
# y= exp.(3.0 .* x1 .+ 2.0 .* x2 .+ sin.(id1) .+ cos.(id2).^2 .+ randn(rng, N))
# df = DataFrame(id1 = categorical(id1), id2 = categorical(id2), x1 = x1, x2 = x2, y = y)
# df = DataFrame(id1_noncat = id1, id1 = categorical(id1), id2 = categorical(id2), x1 = x1, x2 = x2, y = y)
#
# f = Alpaca.@formula(y ~ x1 + x2)
# result = Alpaca.feglm(df, Alpaca.@formula(y ~ x1 + x2),
#     Poisson(),
#     fe = :id2,
#     start = [0.2;0.2], trace = 2
#     )

# string interface (not essential)
# rr1 = Alpaca.feglm(df, "binary ~ SepalWidth | SpeciesDummy", "Binomial")
# @test StatsBase.coef(rr1) ≈ [-0.221486] atol = 1e-4
# rr1 = Alpaca.feglm(df, "binary ~ SepalWidth | SpeciesDummy", "Binomial", vcov = :(cluster(SpeciesDummy)))
# @test StatsBase.stderror(rr1) ≈ [0.317681] atol = 1e-4

ext: jl


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdnotification

"""
    'feffery'_antdnotification(;kwargs...)

An AntdNotification component.

Keyword arguments:
- `id` (String; optional)
- `bottom` (Real; optional)
- `className` (String; optional)
- `description` (String; optional)
- `duration` (Real; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `message` (String; optional)
- `placement` (a value equal to: 'topLeft', 'topRight', 'bottomLeft', 'bottomRight'; optional)
- `style` (Dict; optional)
- `top` (Real; optional)
- `type` (optional)
"""
function 'feffery'_antdnotification(; kwargs...)
        available_props = Symbol[:id, :bottom, :className, :description, :duration, :loading_state, :message, :placement, :style, :top, :type]
        wild_props = Symbol[]
        return Component("'feffery'_antdnotification", "AntdNotification", "feffery_antd_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: limmyforget/feffery-antd-components

max_issues_repo_name: limmyforget/feffery-antd-components

max_forks_repo_name: limmyforget/feffery-antd-components


text:
if basename(pwd()) == "aoc"
    cd("2016/2")
end

offset(c::Char) =
    if c == 'U'
        CartesianIndex(-1, 0)
    elseif c == 'D'
        CartesianIndex(1, 0)
    elseif c == 'L'
        CartesianIndex(0, -1)
    elseif c == 'R'
        CartesianIndex(0, 1)
    end

function part1(filename::AbstractString)
    keypad = [1 2 3; 4 5 6; 7 8 9]
    position = CartesianIndex(2, 2)
    result = Int[]
    for line in eachline(filename)
        for c in line
            newposition = position + offset(c)
            if checkbounds(Bool, keypad, newposition)
                position = newposition
            end
        end
        push!(result, keypad[position])
    end
    foldl((a, b) -> 10a + b, result, init = 0)
end

@assert part1("example.txt") == 1985
part1("input.txt")

function part2(filename::AbstractString)
    keypad = [
        '0' '0' '1' '0' '0'
        '0' '2' '3' '4' '0'
        '5' '6' '7' '8' '9'
        '0' 'A' 'B' 'C' '0'
        '0' '0' 'D' '0' '0'
    ]
    position = CartesianIndex(3, 1)
    result = ""
    for line in eachline(filename)
        for c in line
            newposition = position + offset(c)
            if checkbounds(Bool, keypad, newposition) && keypad[newposition] != '0'
                position = newposition
            end
        end
        result *= keypad[position]
    end
    result
end

@assert part2("example.txt") == "5DB3"
part2("input.txt")

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
mutable struct QROperator{CO,MT,T} <: Operator{T}
    R_cache::CO
    H::MT # Contains the Householder reflections
    ncols::Int   # number of cols already upper triangularized
end

QROperator(R::CachedOperator,H::AbstractArray,ncs::Int) =
    QROperator{typeof(R),typeof(H),eltype(H)}(R,H,ncs)

convert(::Type{Operator{T}},QR::QROperator) where {T} =
    QROperator(convert(Operator{T},QR.R_cache), convert(AbstractArray{T}, QR.H),QR.ncols)

qr(QR::QROperator) = QR
factorize(QR::QROperator) = QR

for OP in (:domainspace,:rangespace)
    @eval $OP(QR::QROperator) = $OP(QR.R_cache)
end

getindex(QR::QROperator, k::Integer, j::Integer) = QR.R_cache.op[k,j]

struct QROperatorR{QRT,T} <: Operator{T}
    QR::QRT
end

QROperatorR(QR) = QROperatorR{typeof(QR),eltype(QR)}(QR)
domainspace(R::QROperatorR) = domainspace(R.QR)
rangespace(R::QROperatorR) = ℓ⁰

function getindex(R::QROperatorR,k::Integer,j::Integer)
    if j < k
        zero(eltype(R))
    else
        resizedata!(R.QR,:,j)
        R.QR.R_cache[k,j]
    end
end

bandwidths(R::QROperatorR) = 0,bandwidth(R.QR.R_cache,2)

struct QROperatorQ{QRT,T} <: Operator{T}
    QR::QRT
end

QROperatorQ(QR) = QROperatorQ{typeof(QR),eltype(QR)}(QR)

domainspace(Q::QROperatorQ) = ℓ⁰
rangespace(Q::QROperatorQ) = rangespace(Q.QR)

getindex(Q::QROperatorQ, k::Integer, j::Integer) = (mul_coefficients(Q',eltype(Q)[zeros(k-1);1]))[j]

function getproperty(F::QROperator, d::Symbol)
    if d == :R
        return QROperatorR(F)
    elseif d == :Q
        return QROperatorQ(F)
    else
        getfield(F, d)
    end
end

# iteration for destructuring into components
Base.iterate(S::QROperator) = (S.Q, Val(:R))
Base.iterate(S::QROperator, ::Val{:R}) = (S.R, Val(:done))
Base.iterate(S::QROperator, ::Val{:done}) = nothing

# override for custom data types
QROperator(R::CachedOperator{T,AM}) where {T,AM<:AbstractMatrix} =
    error("Cannot create a QR factorization for $(typeof(R))")

adjoint(Q::QROperatorQ) = Adjoint(Q)
size(Q::Adjoint{<:Any,<:QROperatorQ}) = (size(parent(Q),2), size(parent(Q),1))

function qr!(A::CachedOperator; cached::Int=0)
    QR = QROperator(A)
    if cached ≠ 0
        resizedata!(QR,:,cached)
    end
    QR
end

"""
    qr(A::Operator)

returns a cached QR factorization of the Operator `A`.  The result `QR`
enables solving of linear equations: if `u=QR\b`, then `u`
approximately satisfies `A*u = b`.
"""
function qr(A::Operator; cached::Int=0)
    if isambiguous(domainspace(A)) || isambiguous(rangespace(A))
        throw(ArgumentError("Only non-ambiguous operators can be factorized."))
    end
    qr!(cache(A;padding=true);cached=cached)
end

factorize(A::Operator) = qr(A)

for OP in (:qr, :factorize)
    @eval begin
        $OP(A::AbstractVector{<:Operator}) = $OP(interlace(A))
        $OP(A::AbstractMatrix{<:Operator}) = $OP(interlace(A))
        $OP(A::AbstractArray{<:Operator}) = $OP(interlace(A))
    end
end

function det(R::QROperatorR;maxiterations::Int=10_000)
    QR = R.QR
    RD = R.QR.R_cache
    resizedata!(QR,:,1)
    ret = -RD[1,1]
    k = 2
    while abs(abs(RD[k-1,k-1])-1) > eps(eltype(R))
        resizedata!(QR,:,k)
        ret *= -RD[k,k]
        k+=1
        k > maxiterations && error("Determinant unlikely to converge after 10_000 iterations.")
    end

    ret
end
det(QR::QROperatorQ) = 1

det(QR::QROperator) = det(QR.R)
det(A::Operator) = det(qr(A))

## Multiplication routines

# Q

mul_coefficients(At::Transpose{T,<:QROperatorQ{T}},B::AbstractVector{T}) where {T<:Real} = parent(At)'*B
mul_coefficients(At::Transpose{T,<:QROperatorQ{T}},B::AbstractMatrix{T}) where {T<:Real} = parent(At)'*B

mul_coefficients(Ac::Adjoint{T,<:QROperatorQ{QR,T}},B::AbstractVector{T};tolerance=eps(eltype(Ac))/10,maxlength=1000000) where {QR,T} =
        mulpars(Ac,B,tolerance,maxlength)

mul_coefficients(Ac::Adjoint{T,<:QROperatorQ{QR,T}},B::AbstractVector{V};opts...) where {QR,T,V} =
    mul_coefficients(Ac,AbstractVector{T}(B); opts...)

function *(Ac::Adjoint{<:Any,<:QROperatorQ}, b::AbstractVector; kwds...)
    A = parent(Ac)
    Fun(domainspace(A),mul_coefficients(Ac,coefficients(b,rangespace(A));kwds...))
end

function *(Ac::Adjoint{<:Any,<:QROperatorQ}, b; kwds...)
    A = parent(Ac)
    Fun(domainspace(A),mul_coefficients(Ac,coefficients(b,rangespace(A));kwds...))
end

ldiv_coefficients(A::QROperatorQ, B; opts...) = mul_coefficients(A', B; opts...)
\(A::QROperatorQ, B::Fun; opts...) = *(A', B; opts...)

# R
function ldiv_coefficients(R::QROperatorR, b::AbstractVector)
    if length(b) > R.QR.ncols
        # upper triangularize columns
        resizedata!(R.QR, :, length(b))
    end
    UpperTriangular(view(R.QR.R_cache.data, 1:length(b), 1:length(b))) \ b
end

\(R::QROperatorR,b::Fun{SequenceSpace};kwds...) =
    Fun(domainspace(R),ldiv_coefficients(R,b.coefficients;kwds...))
\(A::QROperatorR,b::Fun;kwds...) = error("\\ not implement for $(typeof(b)) right-hand sides")

# QR

for TYP in (:Real,:Complex,:Number)
    @eval ldiv_coefficients(QR::QROperator{CO,MT,T},b::AbstractVector{T}; kwds...) where {CO,MT,T<:$TYP} =
        ldiv_coefficients(QR.R, mul_coefficients(QR.Q',b;kwds...))
end

function ldiv_coefficients(QR::QROperator{CO,MT,T},b::AbstractVector{V};kwds...) where {CO,MT,T,V<:Number}
    TV = promote_type(T,V)
    ldiv_coefficients(convert(Operator{TV}, QR),convert(Vector{TV}, b);kwds...)
end

function ldiv_coefficients(QR::QROperator{CO,MT,T},b::AbstractVector{V};kwds...) where {CO,MT,T<:Real,V<:Complex}
    a=ldiv_coefficients(QR,real(b);kwds...)
    b=im*ldiv_coefficients(QR,imag(b);kwds...)
    n=max(length(a),length(b))
    pad!(a,n)+pad!(b,n)
end
ldiv_coefficients(QR::QROperator{CO,MT,T},b::AbstractVector{V};kwds...) where {CO,MT,T<:Complex,V<:Real} =
    ldiv_coefficients(QR,Vector{T}(b);kwds...)

\(A::QROperator,b::Fun;kwds...) =
    Fun(domainspace(A),ldiv_coefficients(A,coefficients(b,rangespace(A));kwds...))

\(A::QROperator,B::MatrixFun;kwds...) = \(A,Array(B);kwds...)


text:
push!(LOAD_PATH,"../src/")

using Pkg
Pkg.activate("../")

using Documenter
using sparse

makedocs(
    sitename = "sparse",
    format = Documenter.HTML(),
    modules = [sparse]
)


text:
using StochasticDiffEq,  DiffEqNoiseProcess, Test

f1(u, p, t) = 1.01u
g1(u, p, t) = 1.01u
dt = 1//2^(4)
prob1 = SDEProblem(f1,g1,1.0,(0.0,1.0))
sol1 = solve(prob1,EM(),dt=dt,save_noise = true)

W2 = NoiseWrapper(sol1.W)
prob2 = SDEProblem(f1,g1,1.0,(0.0,1.0),noise=W2)
dt = 1//2^(5)
sol2 = solve(prob2,EM(),dt=dt)

W3 = NoiseWrapper(sol1.W)
prob3 = SDEProblem(f1,g1,1.0,(0.0,1.0))
dt = 1//2^(5)
sol3 = solve(prob3,EM(),dt=dt,noise=W3)

using Plots
plotly()
plot(sol1, label="reference")
plot!(sol2, label="good coarse reconstruction")
plot!(sol3, label="bad coarse reconstruction")


text:
abstract type RenderState end
abstract type AlphaBlendState <: RenderState end
abstract type StencilState <: RenderState end
abstract type DepthState <: RenderState end
abstract type CullState <: RenderState end

mutable struct AlphaBlendDisabled <: AlphaBlendState
end

function setstate(::AlphaBlendDisabled)
	glDisable(GL_BLEND)
end

mutable struct AlphaBlendSrcAlpha <: AlphaBlendState
end

function setstate(::AlphaBlendSrcAlpha)
	glEnable(GL_BLEND)
	glBlendEquationSeparate(GL_FUNC_ADD, GL_FUNC_ADD)
	glBlendFuncSeparate(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA, GL_SRC_ALPHA, GL_ZERO)
end

mutable struct AlphaBlendAdditive <: AlphaBlendState
end

function setstate(::AlphaBlendAdditive)
	glEnable(GL_BLEND)
	glBlendEquationSeparate(GL_FUNC_ADD, GL_FUNC_ADD)
	glBlendFuncSeparate(GL_SRC_ALPHA, GL_ONE, GL_SRC_ALPHA, GL_ZERO)
end

mutable struct AlphaBlendConstant <: AlphaBlendState
	color::Tuple{Float32, Float32, Float32, Float32}
end

function setstate(state::AlphaBlendConstant)
	glEnable(GL_BLEND)
	glBlendColor(state.color...)
	glBlendEquationSeparate(GL_FUNC_ADD, GL_FUNC_ADD)
	@assert glGetError() == GL_NO_ERROR
	glBlendFuncSeparate(GL_CONSTANT_COLOR, GL_ONE_MINUS_CONSTANT_COLOR, GL_CONSTANT_ALPHA, GL_ZERO)
	@assert glGetError() == GL_NO_ERROR
end

mutable struct StencilStateDisabled <: StencilState
end

function setstate(::StencilStateDisabled)
	glDisable(GL_STENCIL_TEST)
end

mutable struct DepthStateDisabled <: DepthState
end

function setstate(::DepthStateDisabled)
	glDisable(GL_DEPTH_TEST)
end

mutable struct DepthStateLess <: DepthState
end

function setstate(::DepthStateLess)
	glEnable(GL_DEPTH_TEST)
	glDepthFunc(GL_LESS)
end

mutable struct CullStateDisabled <: CullState
end

function setstate(::CullStateDisabled)
	glDisable(GL_CULL_FACE)
end

mutable struct CullStateCCW <: CullState
end

function setstate(::CullStateCCW)
	glEnable(GL_CULL_FACE)
	# Poly orientation is CCW and we cull the back faces
	glFrontFace(GL_CCW)
	glCullFace(GL_BACK)
end

mutable struct RenderStateHolder
	# we store the states with their type's supertype as key, i.e. all the AlphaBlendStates will go to the same position
	states::Dict{DataType, RenderState}

	RenderStateHolder() = new(Dict{DataType, RenderState}())
end

function resetstate(holder::RenderStateHolder, state::RenderState)
	@assert supertype(T) == RenderState
	delete!(holder.states, state)
end

function setstate(holder::RenderStateHolder, state::RenderState)
	holder.states[supertype(typeof(state))] = state
end

function getstate(holder::RenderStateHolder, ::Type{T}, default) where T <: RenderState
	@assert supertype(T) == RenderState
	get(holder.states, T, default)
end

function set_and_apply(holder::RenderStateHolder, state::RenderState)
	key = supertype(typeof(state))
	if !haskey(holder.states, key) || holder.states[key] != state
		holder.states[key] = state
		setstate(state)
	end
end


text:
# Unit test for the lin-solve mehtods (in src/LinSolver.jl)
# The "Gun" problem form gun_native.jl

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra

@testset "linsolvers" begin
    TOL = 1e-10;

    nep = nep_gallery("nlevp_native_gun")

    n = size(nep,1);
    λ = 250^2+1im;

    x=ones(n);

    M = compute_Mder(nep, λ)

    @bench @testset "basic functionality" begin
        linsolver1=create_linsolver(BackslashLinSolverCreator(),nep,λ);
        y1=lin_solve(linsolver1,x)

        linsolver2=create_linsolver(FactorizeLinSolverCreator(),nep,λ);
        y2=lin_solve(linsolver1,x)

        @test y1 ≈ y2

        y3=M\x;
        @test y1 ≈ y3
        @test y2 ≈ y3
        F=factorize(M);
        y4=F\x;
        @test y1 ≈ y4
        @test y2 ≈ y4
        @test y3 ≈ y4

        normy1=norm(y1);
        normy2=norm(y2);
        normy3=norm(y3);
        normy4=norm(y4);

        @info "norm(y1,1)=$normy1"
        # To display further digits
        @info "normy1-round(128*normy1)/128=$(normy1-round(128*normy1)/128) (for further digits)"
        @info "norm(y2,1)=$(norm(y2))"
        @info "normy2-round(128*normy2)/128=$(normy2-round(128*normy2)/128)"
        @info "norm(y3,1)=$(norm(y3))"
        @info "normy3-round(128*normy3)/128=$(normy3-round(128*normy3)/128)"
        @info "norm(y4,1)=$(norm(y4))"
        @info "normy4-round(128*normy4)/128=$(normy4-round(128*normy4)/128)"

        M1norm=norm(M,1);
        r1=M*y1-x;
        r2=M*y2-x;
        r3=M*y3-x;
        r4=M*y4-x;
        @test norm(r1)/M1norm < eps()

        # By first multiplying with M
        z=M*x;
        normz=norm(z)
        @info "z=M*x; norm(z)=$normz x'*z=$(x'*z)"
        @info "normz-floor(normz*10)/10=$(normz-Int(floor(normz*10))/10)"
        t=M\z;
        @info "norm(M\\z-x)=$(norm(t-x))"
        @test norm(t-x)/M1norm < eps()
    end

    @bench @testset "gmres" begin

    end

    @bench @testset "umfpack iterative refinements" begin
        # Test using max 0 vs max 2 iterative refinements. There is little point in testing
        # max 1 refinement, since in practice only a single refinement is used anyway.
        c0 = create_linsolver(FactorizeLinSolverCreator(umfpack_refinements=0),nep, λ)
        c2 = create_linsolver(FactorizeLinSolverCreator(),nep, λ)

        # do a warmup-solve to force JIT compilation to get accurate memory measurements
        lin_solve(c0, x)

        r0 = @timed lin_solve(c0, x)
        r2 = @timed lin_solve(c2, x)

        # test that the norm is at least 10 times bigger with no iterative refinements
        norm0 = norm(M*r0[1] - x)
        norm2 = norm(M*r2[1] - x)
        @test norm0 > 10 * norm2

        # test that much less memory is allocated with no iterative refinements
        memory0 = r0[3]
        memory2 = r2[3]
        @test memory0 < memory2 * 0.75
    end

    @testset "eigsolvers" begin

        @bench @testset "DefaultEigSolver - full" begin
            Random.seed!(0)
            n = 20
            A = rand(ComplexF64, n, n)
            B = rand(ComplexF64, n, n)

            eigsolver1 = DefaultEigSolver(A,B)
            eigsolver2 = DefaultEigSolver(B\A)

            λ1,v1 = eig_solve(eigsolver1, nev=3)
            λ2,v2 = eig_solve(eigsolver2, nev=3)

            for i = 1:3
                @test λ1[i] ≈ λ2[i]
            end
        end

        @bench @testset "DefaultEigSolver - sparse" begin
            Random.seed!(0)
            n = 20
            A = sprand(ComplexF64, n, n, 0.25) + I
            B = sprand(ComplexF64, n, n, 0.25) + I

            eigsolver1 = DefaultEigSolver(A,B)
            eigsolver2 = DefaultEigSolver(sparse(Matrix(B)\Matrix(A)))

            λ1,v1 = eig_solve(eigsolver1, nev=3)
            λ2,v2 = eig_solve(eigsolver2, nev=3)

            for i = 1:3
                @test λ1[i] ≈ λ2[i]
            end
        end

        @bench @testset "EigenEigSolver" begin
            Random.seed!(0)
            n = 20
            A = rand(ComplexF64, n, n)
            B = rand(ComplexF64, n, n)

            eigsolver1 = EigenEigSolver(A,B)
            eigsolver2 = EigenEigSolver(B\A)

            λ1,v1 = eig_solve(eigsolver1, nev=3)
            λ2,v2 = eig_solve(eigsolver2, nev=3)

            for i = 1:3
                @test λ1[i] ≈ λ2[i]
            end
        end

        @bench @testset "ArnoldiEigSolver" begin
            Random.seed!(0)
            n = 20
            A = sprand(ComplexF64, n, n, 0.25) + I
            B = sprand(ComplexF64, n, n, 0.25) + I

            eigsolver1 = ArnoldiEigSolver(A,B)
            eigsolver2 = ArnoldiEigSolver(sparse(Matrix(B)\Matrix(A)))

            λ1,v1 = eig_solve(eigsolver1, nev=3)
            λ2,v2 = eig_solve(eigsolver2, nev=3)

            for i = 1:3
                @test λ1[i] ≈ λ2[i]
            end
        end

    end
end


text:
# This file contains exceptions for methods
# related to the time series generator

@exception(
    "Exception thrown when the file is not found.",
    "Gen",
    "ScatsGenNotAFile",
    "The file is not found (\\\"\", e.file, \"\\\").",
)

@exception(
    "Exception thrown when the path is a directory.",
    "Gen",
    "ScatsGenIsADir",
    "Specified path is a directory (\\\"\", e.file, \"\\\").",
)

@exception(
    "Exception thrown when an end of file occurred.",
    "Gen",
    "ScatsGenEOF",
    "Unexpected end of file (\\\"\", e.file, \"\\\").",
)

@exception(
    "Exception thrown when an input for `N` is wrong.",
    "Gen",
    "ScatsGenWR_N",
    "Wrong input: N (\\\"\", e.file, \"\\\").",
)

@exception(
    "Exception thrown when an input for `Δt` is wrong.",
    "Gen",
    "ScatsGenWR_Δt",
    "Wrong input: Δt (\\\"\", e.file, \"\\\").",
)

@exception(
    "Exception thrown when an input for `q` is wrong.",
    "Gen",
    "ScatsGenWR_q",
    "Wrong input: q (\\\"\", e.file, \"\\\").",
)

@exception(
    "Exception thrown when an input for `α` is wrong.",
    "Gen",
    "ScatsGenWR_α",
    "Wrong input: α (\\\"\", e.file, \"\\\").",
)

@exception(
    "Exception thrown when an input for `β` is wrong.",
    "Gen",
    "ScatsGenWR_β",
    "Wrong input: β (\\\"\", e.file, \"\\\").",
)

@exception(
    "Exception thrown when an input for `r` is wrong.",
    "Gen",
    "ScatsGenWR_r",
    "Wrong input: r (\\\"\", e.file, \"\\\").",
)

@exception(
    "Exception thrown when an input for `A` is wrong.",
    "Gen",
    "ScatsGenWR_A",
    "Wrong input: A (\\\"\", e.file, \"\\\").",
)

@exception(
    "Exception thrown when an input for `ν` is wrong.",
    "Gen",
    "ScatsGenWR_ν",
    "Wrong input: ν (\\\"\", e.file, \"\\\").",
)

@exception(
    "Exception thrown when an input for `ϕ` is wrong.",
    "Gen",
    "ScatsGenWR_ϕ",
    "Wrong input: ϕ (\\\"\", e.file, \"\\\").",
)

@exception(
    "Exception thrown when an input for `γ` is wrong.",
    "Gen",
    "ScatsGenWR_γ",
    "Wrong input: γ (\\\"\", e.file, \"\\\").",
)


text:
# include("kernel.jl")
# include("DCC.jl")
# From this line
# This function is: util(c,η,s) = (1+η) c^(1-σ)/(1-σ)
mutable struct Utility
    σ::Float64
    Utility(σ::Real) = new(Float64(σ))
    function (self::Utility)(c::Union{Real,RealVector})
        if (self.σ == 1)
            return(log.(c))
        else
            return(float.(c).^(1 - self.σ)./(1 - self.σ));
        end
    end

    function (self::Utility)(c::Union{Real,RealVector},η::Union{Real,RealVector})
        if (self.σ == 1)
            return(exp.(η).*log.(c))
        else
            return(exp.(η).* (float.(c).^(1 - self.σ))./(1 - self.σ));
        end
    end

end

mutable struct Transition
    α::Real
    A::Real
    Transition(α,A)=new(α,A);
    function (self::Transition)(c::Union{Real,RealVector},s::Union{Real,RealVector})
    # return(1.1 * s.^(self.α) .- c)
    # return(α*(s-c.+ 2.0))
    return (self.A * (s.^(self.α)) .- c)
    end
end

# function Transition(s::Real,c::Real)
    # return(1.05 * ( s - c ) + 3 )
# end

# function Transition(s::RealVector,c::RealVector)
    # return(1.05 * ( s - c ) .+ 3 )
# end

# How to use sub types?
function find_optim(xin::Real,ηin::Real, lb::Real,ub::Real, β::Float64, trans::Transition, util::Utility,ValueFn::ApproxFn)
    @suppress begin
        # η = 0.01;
        ff = c ->  - (util(c' * [1],ηin) + β * ValueFn(trans(c' * [1],xin)));
        # ub = trans(0,xin[1])-lb;
        @show init_val = (ub + lb) / 2;
        f_opt = optimize(ff,[lb],[ub],[lb],SAMIN(),Optim.Options(g_tol = 1e-12,
                         iterations = 1000,
                     store_trace = false,
                     show_trace = false));
        return(f_opt.minimizer[1])
    end
end

mutable struct State
    s::Union{Real,RealVector} #Observed
    η::Real #Unobserved
    State()=new(0,0);
    function State(s::Real, η::Real)
        new([s],η);
    end

    function State(s::Array, η::Real)
        new(s,η);
    end
end

function vectorize(s::State)
    return(vcat(s.s,s.η))
end

mutable struct DynamicDecisionProcess
    σ::Float64
    util::Utility
    trans::Transition
    PolicyFn::ApproxFn
    ValueFn::ApproxFn
    β::Float64
    dtrans::Function
    dutil::Function
    nSolve::Int
    function DynamicDecisionProcess(σ::Real,β::Float64,α::Real,A::Real)
        nSolve = 400;
        v_tol = 0.001;
        ϵ     = 0.05;
        # step_length= 2/nSolve;
        step_length=ϵ;
        util  = Utility(σ);
        trans = Transition(α,A);
        # x = hcat(range(2*ϵ,step= ϵ,length=nSolve)); #Convert it into 2dimension
        s = convert(Array{Float64,1},range(ϵ,step=step_length,length=nSolve));
        # TN=TruncatedNormal(0, 1, 0, 10);
        # s = rand(TN,nSolve);
        #Convert it into 2dimension
        η = randn(nSolve)/100; #Normal distributed error
        # η=zeros(nSolve);
        # v = (util(s)./(1 -β));
        sdat = hcat(s,η);

        s = vcat([1e-10*ϵ],s);
        c_opt = zeros(nSolve);
        v_first=zeros(nSolve+1);
        PolicyFn = ApproxFn(sdat,c_opt,:gaussian,2);
        ValueFn = ApproxFn(s,v_first,:gaussian,2);
        # PolicyFn.h =  0.5 * PolicyFn.h;
        # ValueFn.h =  0.5 * ValueFn.h ;
        # ValueFn.h = 1; #Try this
        dtrans = (c,s) -> Tracker.gradient(trans,c,s);
        dutil  = (c,η) -> Tracker.gradient(util,c,η)[1].data;
        ddc = new(float(σ),util,trans,PolicyFn,ValueFn,β,dtrans,dutil,nSolve);
        return(ddc);
    end
end

# dtrans = (c,s) -> Tracker.gradient(ddc.trans,c,s) #Transition derivatives, can be broadcasted
# Compute equilibrium(): similar to Paul Schrimpf's dynamic choice problem

# While normV + normS2 + normV2 > tol (What is tole)
    # UpateValue, updateStrategy     : Update the solved value for policy and value function
    # gridUpdate && abs2(strategy.val,stratGrid.val)>strategy.val.size()*0.2
    # # UpdateSolvedGrid: To make the distribution stationary
    # stratGrid = strategy

function computeDistance(ddc::DynamicDecisionProcess,old_policy::ApproxFn,old_value::ApproxFn)
     normPolicy = abs2(old_policy.y - ddc.PolicyFn.y)./ddc.PolicyFn.n;
     normValue  = abs2(old_value.y  - ddc.ValueFn.y)./ddc.ValueFn.n;
     normData   = mean(abs2(old_value.xdata - ddc.ValueFn.xdata));
     return(normPolicy+normValue+normData);
end

function UpdateSolvedGrid!(ddc::DynamicDecisionProcess,T)
    # Update the grade for solving the points
    # ddc.policy
    η_new = randn(ddc.nSolve)/100; #Generate η
    # η_new = zeros(ddc.nSolve); #Generate η
    s_new = ddc.ValueFn.xdata;
    x_new = hcat(s_new,η_new);
    s_old = s_new;
    val_new = val_old = ddc.ValueFn.y;
    t = 0;
    while t < T
        # @show t;
        c_new = ddc.PolicyFn(x_new);
        s_new = ddc.trans(c_new,s_old);
        UpdateData(ddc.ValueFn,s_old,val_old);

        val_new = ddc.util(c_new,x_new[:,2]) + ddc.β * ddc.ValueFn(s_new);
        val_old = deepcopy(val_new);
        s_new[s_new.<0].=minimum(s_old);
        η_new = randn(ddc.nSolve); #Generate η
        x_new = hcat(s_new,η_new);
        s_diff= abs2(s_old - s_new)/ddc.nSolve;
        s_old = deepcopy(s_new);
        t += 1;
    end

end

function UpdateVal!(ddc::DynamicDecisionProcess,T::Int)
    c_opt = ddc.PolicyFn.y;
    s = ddc.PolicyFn.xdata[:,1];
    η = ddc.PolicyFn.xdata[:,2];
    iter = 0
    tol = 1e-8
    v_diff = Inf
    while (iter < T && v_diff > tol)
        ϵ = minimum(s)*0.9;
        lb = ϵ;
        for n = 1:ddc.nSolve
            ub = max(0.9 * ddc.trans(0,s[n]),ddc.trans(0,s[n])-2*ϵ);
            c_opt[n] = find_optim(s[n],η[n],lb,ub,ddc.β,ddc.trans,ddc.util,ddc.ValueFn);
        end

        y = ddc.util(c_opt) + ddc.β * ddc.ValueFn(ddc.trans(c_opt,s));
        # println("Iteration:",iter);
        # @show iter;
        # @show v_diff = mean(abs.(y - ddc.ValueFn.y));
        UpdateVal(ddc.ValueFn,vcat([ddc.ValueFn.y[1]],y));
        @show iter += 1;
    end
    UpdateVal(ddc.PolicyFn,c_opt);
end

function computeEquilibrium(ddc::DynamicDecisionProcess)
    i = 0;
    diff_v=Inf;
    diff_v_old = 0;
    diff_r=Inf;
    # tol = 1/ddc.nSolve;
    tol_r = 0.1;
    tol_v = 1e-5;
    # while (i < 30) & (diff_v > tol)
    while (i < 30) & (diff_r > tol_r) & (diff_v > tol_v)
        old_value=deepcopy(ddc.ValueFn); old_policy=deepcopy(ddc.PolicyFn);
        UpdateVal!(ddc,10);
        @show diff_v=computeDistance(ddc,old_policy,old_value);
        @show diff_r = abs(diff_v-diff_v_old) / diff_v;
        diff_v_old = diff_v;
        # UpdateSolvedGrid!(ddc::DynamicDecisionProcess,2);
        i+=1;
    end
end

# simulate dynamic discrete choice problem from the solved problem
function simulate_ddc(nM,nT,ddc::DynamicDecisionProcess)
    # lb = minimum(ddc.PolicyFn.xdata,dims=1);
    # ub = maximum(ddc.PolicyFn.xdata,dims=1);
    # if lb[2]==ub[2]
    #     ub[2] = lb[2]+1e-6
    # end
    # x0 = zeros(nM,ddc.PolicyFn.q); #Initial states
    # for i = 1:ddc.PolicyFn.q
    #     du = Uniform(lb[i],ub[i]);
    #     x0[:,i] = rand(du,nM);
    # end
    # du=DiscreteUniform(1,ddc.nSolve); #Draw from the nSolved states
    # s = zeros(nM,nT+1);
    # a = zeros(nM,nT);
    # s[:,1] = x0[:,1];
    # du = Uniform(lb[2],ub[2]);

    # du=DiscreteUniform(1,ddc.nSolve); #Draw from the nSolved states
    du=Uniform(2,20);
    s = [ State() for i=zeros(nM,nT)]; #nM x nT state
    a = [ 0.0 for i=zeros(nM,nT)]; #nM x nT action
    for n = 1:nM
        ind_1=rand(du);
        # s[n,1] = State(ddc.PolicyFn.xdata[ind_1,1:end-1],randn()/100);
        s[n,1] = State(rand(du),randn()/100);
        for t = 1:(nT-1)
        # x = hcat(s[:,t],randn(nM)/3);
            x = vectorize(s[n,t]); #This is a bit problematic
            a[n,t] = ddc.PolicyFn(x);
            s_1 = ddc.trans(a[n,t],s[n,t].s);
            if s_1[1] < 0
                @show n,t
                break
            end
            s[n,t+1] = State(s_1,randn()/100);

        end
        x = vectorize(s[n,nT]); #T
        a[n,nT] = ddc.PolicyFn(x);
    end

    return(state=s,action=a);
end

function check_ee(ddc::DynamicDecisionProcess)
    s0=rand(Uniform(3,15));
    ϵ0=0;
    a0=ddc.PolicyFn([s0,ϵ0])[1];
    s1=ddc.trans(a0,s0);
    a1=ddc.PolicyFn([s1,0])[1];
    ds=ddc.dtrans(a1,s1)[2].data;
    return(ddc.dutil(a0,0)- ds *ddc.β  * ddc.dutil(a1,0));
end

mutable struct EulerEquation
    σ::Float64
    util::Utility
    trans::Transition
    PolicyFn::ApproxFn
    β::Float64
    dtrans::Function
    dutil::Function
    nSolve::Int
    function EulerEquation(σ::Real,β::Float64,α::Real)
        v_tol = 0.001;
        ϵ     = 0.005;
        util  = Utility(σ);
        trans = Transition(α);
        nSolve = 500;
        # x = hcat(range(2*ϵ,step= ϵ,length=nSolve)); #Convert it into 2dimension
        s = convert(Array{Float64,1},range(2*ϵ+1,step=ϵ,length=nSolve)); #Convert it into 2dimension
        η = randn(nSolve)/100; #Normal distributed error
        # η=zeros(nSolve);
        # v = (util(s)./(1 -β));
        sdat = hcat(s,η);
        c_opt = zeros(nSolve).+ϵ;
        tol=1e-3; diff=Inf;max_iter=10;iter=0;
        dtrans = (c,s) -> Tracker.gradient(trans,c,s);
        dutil  = (c) -> Tracker.gradient(util,c,0)[1].data;

        PolicyFn = ApproxFn(sdat,c_opt,:gaussian,2);
        dtrans_s = (c,s) -> dtrans(c,s)[2].data;

        function f_obj(c)
            diff_val =
            dutil.(c) - β.* dtrans_s.(c,PolicyFn.xdata[:,1]) .* PolicyFn( hcat(trans.(c,PolicyFn.xdata[:,1]),zeros(nSolve)) );
            return(sum(diff_val.^2))
        end

        while (diff > tol) & (iter<max_iter)

            result = optimize(f_obj,c_opt,PolicyFn.xdata[:,1],zeros(nSolve).+ϵ,SAMIN());
            # c_opt=result.minimizer;;
            @show diff = sum( (result.minimizer-PolicyFn.y).^2);
            UpdateVal(PolicyFn,result.minimizer);
            iter+=1;
        end
        ee = new(float(σ),util,trans,PolicyFn,β,dtrans,dutil,nSolve);
        return(ee);
    end
end

function check_ee(ee::EulerEquation)
    s0=rand();
    ϵ0=0;
    a0=ee.PolicyFn([s0,ϵ0])[1];
    s1=ee.trans(a0,s0);
    a1=ee.PolicyFn([s1,0])[1];
    return(ee.dutil(a0)- 1.05*ee.β * ee.dutil(a1));
end

function simulate_ee(nM,nT,ee::EulerEquation)
    lb = minimum(ee.PolicyFn.xdata,dims=1);
    ub = maximum(ee.PolicyFn.xdata,dims=1);
    if lb[2]==ub[2]
        ub[2] = lb[2]+1e-6
    end
    x0 = zeros(nM,ee.PolicyFn.q); #Initial states
    for i = 1:ee.PolicyFn.q
        du = Uniform(lb[i],ub[i]);
        x0[:,i] = rand(du,nM);
    end
    s = zeros(nM,nT+1);
    a = zeros(nM,nT);
    s[:,1] = x0[:,1];
    for t = 1:nT
        # x = hcat(s[:,t],randn(nM)/3);
        x = hcat(s[:,t],zeros(nM));
        a[:,t] = ee.PolicyFn(x);
        s[:,t+1] = ee.trans(a[:,t],s[:,t]);
    end
    return(state=s,action=a);
end


text:
module DummyModule

using JuMP

export solveOptProb2, MyResults

function solveOptProb2(p)
	a = p[1]
	b = p[2]
	m = Model()
	@variable(m, 0 <= x <= 2 )
	@variable(m, 0 <= y <= 30 )

	@objective(m, Max, a*x + 3*y )
	@constraint(m, 1x + 5y <= b )

	status = solve(m)

	return MyResults(getobjectivevalue(m), getvalue(x), getvalue(y), status)

end

type MyResults
	obj::Float64
	x::Float64
	y::Float64
	status::Symbol
end

function print(r::MyResults)

	println(" ")
	println("objective = ",r.obj)
	println("x = ",r.x,"\t y = ",r.y)
	println(r.status)

	return nothing

end

end


text:
#tes Alt parsing and repr update saving writing to fasta
using st16SseqJuliaTools
using CSV

@testset "AlignmentMatrix" begin
    aln = Alignment("test/playdata/sample.fasta")
    @test aln.repr_n == [3196, 2454, 2280, 1676, 1425, 1340, 1183, 1018, 696, 663, 280, 251, 208, 203, 76, 66, 65, 43, 32, 32, 28, 28, 20, 20, 15, 12, 10, 8, 8, 8, 7, 5, 4, 4, 4, 4, 4, 4]
    write_to_fasta(aln,"test/playdata/out1.fasta")
    aln2 = Alignment("test/playdata/out1.fasta")
    @test aln.repr_n == aln2.repr_n
    write_to_fasta(aln,2:20,"test/playdata/out2.fasta")
    @test success(`cmp --quiet test/playdata/out2.fasta test/playdata/out2_ref.fasta`)
    aln3 = sub_alignment(aln,2:20)
    write_to_fasta(aln3,"test/playdata/out3.fasta")
    @test success(`cmp --quiet test/playdata/out3.fasta test/playdata/out2_ref.fasta`)

    # test one fragment clustering
    cluster1 = Alignment("test/playdata/cluster1.fasta")
    cluster1_ref = Alignment("test/playdata/cluster1_ref.fasta")
    vsearh_cluster!(cluster1,11:110,0.97)
    @test cluster1.M == cluster1_ref.M

    cluster1 = Alignment("test/playdata/cluster1.fasta")
    cluster1_ref = Alignment("test/playdata/cluster1_ref.fasta")
    vsearh_cluster!(cluster1,100,1,0.97)
    @test cluster1.M == cluster1_ref.M

    cluster1 = Alignment("test/playdata/clusteringtest_mafft.fasta")
    cluster1_ref = Alignment("test/playdata/clusteringtest_mafft_ref.fasta")
    vsearh_cluster!(cluster1,100,5,0.97)
    @test cluster1.M == cluster1_ref.M

    #test conservation score

    @test (2.0, 0.0, 'C')  == char_array_entropy(['A','T','G','C'],[0,0,0,4])
    @test (0.0, 0.0, 'A')  == char_array_entropy(['A','T','G','C'],[4,4,4,4])

    #test data collection for chimeric detection
    aln = Alignment("test/playdata/short_chim_4_test.fasta")
    recomb_data = eval_recomb(aln)
    CSV.write("test/playdata/short_chim_4_test_out.csv",recomb_data)
    @test success(`cmp --quiet test/playdata/short_chim_4_test_out.csv test/playdata/short_chim_4_test_ref.csv `)
end


text:
module AmbienceParamTrigger

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Ambience Param" => Ahorn.EntityPlacement(
        Maple.AmbienceParamTrigger,
        "rectangle"
    )
)

function Ahorn.editingOptions(trigger::Maple.AmbienceParamTrigger)
    return Dict{String, Any}(
        "direction" => Maple.trigger_position_modes
    )
end

end


text:
age                                            :: Integer
sex                                            :: Sex
hid                                            :: Int128
pid                                            :: Int128
age                                            :: Integer

ethnic_group                                   :: Ethnic_Group_Type
marital_status                                 :: Marital_Status_Type
highest_qualification                          :: Qualification_Type
industrial_classification                      :: Standard_Industrial_Classification_2007
occupational_classification                    :: Standard_Occupational_Classification
# FIXME needs work on the mapping - leeve out for now
socio_economic_grouping                        :: Socio_Economic_Grouping_Type
age_completed_full_time_education              :: Integer tea
years_in_full_time_work                        :: Integer  ftwk
employment_status                              :: Employment_Status_ILO_Definition
hours_worked                                   :: Natural "tothours" from frs-yr or "injwk" from adult

income                                         :: Incomes_Dict
assets                                         :: Asset_Dict
pension_contributions                          :: Real  job.udeduc1  job.udeduc2
contracted_out_of_serps                        :: Bool ???

registered_blind                               :: Bool spcreg1 adult|child
registered_partially_sighted                   :: Bool spcreg2 adult|child
registered_deaf                                :: Bool spcreg3 adult|child
has_learning_difficulty                        :: Bool disd05 | cdisd05
disabilities                                   :: Dict{Disability_Type,Bool} dis01 .. dis10
health_status                                  :: Health_Status
wealth                                         :: Wealth_List
is_informal_carer                              :: Bool carefl adult|child
receives_informal_care_from_non_householder    :: Bool ??
hours_of_care_received                         :: Real  hourcare
hours_of_care_given                            :: Real  hourtot
hid                           :: BigInt
year                          :: Unsigned
current_simulated_date        :: Date
tenure                        :: Tenure_Type
region                        :: Region_Type

ct_band                       :: Council_Tax_Band_Type
council_tax                   :: "ctamt"

sewerage                      :: Real
water                         :: Real
mortgage_payment              :: Real
mortgage_interest             :: Real
years_outstanding_on_mortgage :: Integer
mortgage_outstanding          :: Real
year_house_bought             :: Unsigned
gross_rent                    :: Real # rentg Gross rent including Housing Benefit  or rent Net amount of last rent payment
rent_includes_water_and_sewerage :: Bool
other_housing_charges         :: Real # rent Net amount of last rent payment
gross_housing_costs           :: Real
total_income                  :: Real
total_wealth                  :: Real
# wealth                        :: Incomes_Dict
house_value                   :: Real
is_care_home                  :: Bool


text:
using Test
using LinearAlgebra
using IntervalArithmetic
using PolynomialBasis
using ImplicitDomainQuadrature
using CartesianMesh
using HDGElasticity

HDGE = HDGElasticity

function allequal(v1,v2)
    return all(v1 .== v2)
end

function allapprox(v1,v2)
    return all(v1 .≈ v2)
end

@test HDGE.default_integer_type() == typeof(1)
@test HDGE.default_float_type() == typeof(1.0)

@test HDGE.reference_element_length() == 2.0

@test HDGElasticity.in_reference_interval(-0.25)
@test !HDGElasticity.in_reference_interval(+1.2)

@test HDGElasticity.reference_interval_1d() == (-1.0,+1.0)

@test_throws ArgumentError HDGE.reference_cell(3)
box = HDGE.reference_cell(1)
testbox = IntervalBox([-1.0],[1.0])
@test allequal(box,testbox)

box = HDGE.reference_cell(2)
testbox = IntervalBox([-1.,-1.],[1.,1.])
@test allequal(box,testbox)

@test_throws ArgumentError HDGElasticity.number_of_faces(1)
@test_throws ArgumentError HDGElasticity.number_of_faces(4)
@test HDGElasticity.number_of_faces(2) == 4
@test HDGElasticity.number_of_faces(3) == 6

normals = HDGE.reference_normals()
@test allequal(normals[1],[0.0,-1.0])
@test allequal(normals[2],[1.0,0.0])
@test allequal(normals[3],[0.0,1.0])
@test allequal(normals[4],[-1.0,0.0])

@test HDGE.reference_cell_volume(1) == 2.0
@test HDGE.reference_cell_volume(2) == 4.0
@test HDGE.reference_cell_volume(3) == 8.0

map = InterpolatingPolynomial(2,1,1)
coeffs = [0.0  0.0
          0.0  1.0]
update!(map,coeffs)
@test HDGE.determinant_jacobian(map,0.0) ≈ 0.5

map = InterpolatingPolynomial(2,1,2)
coeffs = [-1. 0. 1.
           0. -1 0]
update!(map,coeffs)
@test HDGE.determinant_jacobian(map,-0.5) ≈ sqrt(2)

mesh = UniformMesh([0.0,0.0],[1.,1.],[1,1])
basis = TensorProductBasis(2,2)
coords = HDGElasticity.nodal_coordinates(mesh,basis)
testcoords = [0.0 0.0 0.0 0.5 0.5 0.5 1.0 1.0 1.0
              0.0 0.5 1.0 0.0 0.5 1.0 0.0 0.5 1.0]
@test allequal(coords,testcoords)

@test HDGElasticity.neighbor_faceid(1) == 3
@test HDGElasticity.neighbor_faceid(2) == 4
@test HDGElasticity.neighbor_faceid(3) == 1
@test HDGElasticity.neighbor_faceid(4) == 2
@test_throws ArgumentError HDGElasticity.neighbor_faceid(5)
@test_throws ArgumentError HDGElasticity.neighbor_faceid(0)

f(x) = 2x[1] + 3x[2]
box = HDGE.reference_cell(2)
restricted_funcs = HDGE.restrict_on_faces(f,box)
@test allequal(restricted_funcs[1].([-0.5,0.5]),[-4.,-2.0])
@test allequal(restricted_funcs[2].([-0.5,0.5]),[0.5,3.5])
@test allequal(restricted_funcs[3].([-0.5,0.5]),[2.0,4.0])
@test allequal(restricted_funcs[4].([-0.5,0.5]),[-3.5,-0.5])

basis = TensorProductBasis(2,4)
@test HDGElasticity.number_of_basis_functions(basis) == 25
@test HDGElasticity.dimension(basis) == 2
@test HDGElasticity.number_of_basis_functions(basis.basis) == 5

poly = InterpolatingPolynomial(1,1,3)
@test HDGE.dimension(poly) == 1
poly = InterpolatingPolynomial(1,2,3)
@test HDGElasticity.dimension(poly) == 2

quad = tensor_product_quadrature(1,5)
@test HDGElasticity.number_of_quadrature_points(quad) == 5
@test HDGElasticity.dimension(quad) == 1
quad = tensor_product_quadrature(2,6)
@test HDGElasticity.dimension(quad) == 2
@test HDGElasticity.number_of_quadrature_points(quad) == 36

function plane_distance_function(coords,n,x0)
    return [n'*(coords[:,idx]-x0) for idx in 1:size(coords)[2]]
end
poly = InterpolatingPolynomial(1,2,1)
coords = [0. 0. 2.  2.
          0. 1. 0.  1.]
testn = [1.,1.]/sqrt(2.)
coeffs = plane_distance_function(coords,testn,[0.5,0.])
update!(poly,coeffs)
cellmap = HDGElasticity.CellMap([0.,0.],[2.,1.])
p = [-0.25,-0.5]
n = HDGElasticity.levelset_normal(poly,p,cellmap)
@test allapprox(testn,n)

levelsetcoeffs = HDGElasticity.levelset_coefficients(
    x->plane_distance_function(x,[1.,0.],[0.5,0.]),mesh,poly.basis)
testlevelsetcoeffs = [-0.5,-0.5,0.5,0.5]
@test allapprox(levelsetcoeffs,testlevelsetcoeffs)

points = [-0.4  -0.9
          -0.8   -0.2]
n = HDGElasticity.levelset_normal(poly,points,cellmap)
@test allapprox(n,hcat(testn,testn))

normal = [1,2]
@test HDGElasticity.tangents(normal) == [2,-1]

normals = [1 2 3 4
           5 6 7 8]
testt = [5 6 7 8
         -1 -2 -3 -4]
@test allapprox(HDGElasticity.tangents(normals),testt)

cellmap = HDGElasticity.CellMap([0.,0.],[2.,1.])
@test HDGElasticity.scale_area(cellmap,[1.,0.]) ≈ 0.5
@test HDGElasticity.scale_area(cellmap,[0.,1.]) ≈ 1.0
n = [1. 0.
     0. 1.]
@test allapprox(HDGElasticity.scale_area(cellmap,n),[0.5,1.])

facequads = [[],[1,2,3],[1,2,3,4,5],[]]
isactiveface = HDGElasticity.active_faces(facequads)
@test allequal(isactiveface,[0,1,1,0])

cellmap = HDGElasticity.CellMap([3.,2.],[5.,7.])
midpoints = HDGElasticity.face_midpoints(cellmap)
@test length(midpoints) == 4
@test allapprox(midpoints[1],[4.,2.])
@test allapprox(midpoints[2],[5.,4.5])
@test allapprox(midpoints[3],[4.,7.])
@test allapprox(midpoints[4],[3.,4.5])

imap = InterpolatingPolynomial(2,1,2)
coeffs = [-1.,0.,0.,1.,1.,0.]
update!(imap,coeffs)
p = [-0.5,0.5]
normals = HDGElasticity.curve_normals(imap,p)
testn = [-1. 1.
          1. 1.]
testn .*= 1.0/sqrt(2.)
@test allapprox(testn,normals)

mesh = UniformMesh([0.,0.],[1.,1.],[1,1])
levelsetcoeffs = reshape([-1.,-1.,1.,1.],4,1)
levelset = InterpolatingPolynomial(1,2,1)
update!(levelset,levelsetcoeffs[:,1])
dgmesh = HDGElasticity.DGMesh(mesh,levelsetcoeffs,levelset)
ufs = HDGElasticity.UniformFunctionSpace(dgmesh,1,2,levelsetcoeffs,levelset)
hybrid_element_numbering = HDGElasticity.HybridElementNumbering(dgmesh,ufs)
number_of_hybrid_elements = hybrid_element_numbering.number_of_hybrid_elements
ft = typeof(1.0)
coordinates = Vector{Any}(undef,number_of_hybrid_elements)
visited = zeros(Bool,number_of_hybrid_elements)
facetohelid = hybrid_element_numbering.facetohelid
interfacehelid = hybrid_element_numbering.interfacehelid
refcoordinates = [[fm(p) for p in ufs.sbasis.points'] for fm in dgmesh.facemaps]
HDGElasticity.assign_uniform_cell_hybrid_coordinates!(coordinates,visited,1,1,
    dgmesh,facetohelid,interfacehelid,refcoordinates,ufs.icoeffs)
testvisited = zeros(Bool,number_of_hybrid_elements)
testvisited[1:4] .= true
@test allequal(testvisited,visited)
@test allapprox(coordinates[1][1],[0.,0.])
@test allapprox(coordinates[1][2],[1.,0.])
@test allapprox(coordinates[2][1],[1.,0.])
@test allapprox(coordinates[2][2],[1.,1.])
@test allapprox(coordinates[3][1],[0.,1.])
@test allapprox(coordinates[3][2],[1.,1.])
@test allapprox(coordinates[4][1],[0.5,0.])
@test allapprox(coordinates[4][2],[0.5,1.0])

coordinates = HDGElasticity.hybrid_element_coordinates(dgmesh,ufs,hybrid_element_numbering)


text:
#=
  Calculating π

  from Foundations of Probabilistic Programming,
  page 4

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

# Well, using PPL for this is probably not the best approach. :-)
@model function pi_calc()
    limit = 10^6
    i = 0 # ~ DiscreteUniform(1,limit)
    n = 1 # ~ DiscreteUniform(1,limit)

    while i < limit
        x ~ Uniform()
        y ~ Uniform()
        if (x*x+y*y) < 1
            n = n + 1
        end
        i = i + 1
    end
    ppi ~ Uniform(3,4)
    ppi = 4*n/limit
end

model = pi_calc()

num_chns = 4
# chns = sample(model, MH(), MCMCThreads(), 1000, num_chns)
chns = sample(model, MH(), 10)

display(chns)
# show_var_dist_pct(chns, :x)
# show_var_dist_pct(chns, :y)
# show_var_dist_pct(chns, :i)
# show_var_dist_pct(chns, :n)
show_var_dist_pct(chns, :ppi)

#=
# Plain simulation
ppi:3.141694144
π  :π
diff: -0.00010149041020701688
  7.368581 seconds (1.03 M allocations: 49.375 MiB)
=#
function pi_calc2(limit=10^9)

    i = 0
    n = 1

    while i < limit
        x = rand()
        y = rand()
        if (x*x+y*y) < 1
            n = n + 1
        end
        i = i + 1
    end

    return 4*n/limit
end

#=
@time ppi= pi_calc2()
println("ppi:$ppi")
println("π  :$(π)")
println("diff: ", π-ppi)
=#

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
"""
    StraightRoadway
A simple type representing a one lane, one dimensional straight roadway
# Fields
- `length::Float64`
"""
struct StraightRoadway
    length::Float64
end

"""
    mod_position_to_roadway(s::Float64, roadway::StraightRoadway)
performs a modulo of the position `s` with the length of `roadway`
"""
function mod_position_to_roadway(s::Float64, roadway::StraightRoadway)
    while s > roadway.length
        s -= roadway.length
    end
    while s < 0.0
        s += roadway.length
    end
    return s
end

"""
    get_headway(s_rear::Float64, s_fore::Float64, roadway::StraightRoadway)
returns a positive distance between s_rear and s_fore.
"""
function get_headway(s_rear::Float64, s_fore::Float64, roadway::StraightRoadway)
    while s_fore < s_rear
        s_fore += roadway.length
    end
    while s_fore > s_rear + roadway.length
        s_fore -= roadway.length
    end
    return s_fore - s_rear # positive distance
end


text:
using ArrayIterators, Test

mysum(x::EachRow) = vec(sum(parent(x),dims=1))
mysum(x::EachCol) = vec(sum(parent(x),dims=2))

A = reshape(1:20,4,5)
@test mysum(eachrow(A)) == sum(eachrow(A))
@test mysum(eachslice(A,dims=1)) == sum(eachrow(A))

@test mysum(eachcol(A)) == sum(eachcol(A))
@test mysum(eachslice(A,dims=2)) == sum(eachcol(A))