text:
using Base.Test
using LinearControl

# Test System
A = [1 2 3; 2 1 3; 3 2 1]
B = [1; 0; 0]
C = [0 0 1]
D = [0]

num = [1]
den = [1 2 10]

# State Space
myss = SS(A,B,C,D)

# Transfer Function
mytf = TF(num,den)

# Linear Model
mylm1 = LinearModel(myss)
mylm2 = LinearModel(mytf)
@assert is(mylm1,mylm2)

@assert is(mylm1.tf, mytf)
@assert is(mylm2.ss, myss)


text:
"""
`ODE`: Ordinary Differential Equation

Defines an initial value problem
```math
\\dot{q} (t) = f(t, q(t)) , \\qquad q(t_{0}) = q_{0} ,
```
with vector field ``f``, initial condition ``q_{0}`` and the solution
``q`` taking values in ``\\mathbb{R}^{d}``.

### Fields

* `d`: dimension of dynamical variable ``q`` and the vector field ``f``
* `f`: function computing the vector field
* `t₀`: initial time
* `q₀`: initial condition

The function `f` providing the vector field must have the interface
```julia
    function f(t, q, fq)
        fq[1] = ...
        fq[2] = ...
        ...
    end
```
where `t` is the current time, `q` is the current solution vector, and
`fq` is the vector which holds the result of evaluating the vector field ``f``
on `t` and `q`.

"""
immutable ODE{T,F} <: Equation{T}
    d::Int
    n::Int
    f::F
    t₀::T
    q₀::Array{T,2}

    function ODE(d, n, f, t₀, q₀)
        @assert d == size(q₀,1)
        @assert n == size(q₀,2)
        @assert T == eltype(q₀)
        @assert ndims(q₀) ∈ (1,2)

        if ndims(q₀) == 1
            q₀ = reshape(q₀, d, n)
        end

        new(d, n, f, t₀, q₀)
    end
end

function ODE{T,F}(f::F, t₀::Real, q₀::DenseArray{T})
    ODE{T,F}(size(q₀, 1), size(q₀, 2), f, t₀, q₀)
end

function ODE(f, q₀)
    ODE(f, 0, q₀)
end

Base.hash(ode::ODE, h::UInt) = hash(ode.d, hash(ode.n, hash(ode.f, hash(ode.t₀, hash(ode.q₀, h)))))
Base.:(==){T1, T2, F1, F2}(ode1::ODE{T1,F1}, ode2::ODE{T2,F2}) = (
                                ode1.d == ode2.d
                             && ode1.n == ode2.n
                             && ode1.f == ode2.f
                             && ode1.t₀ == ode2.t₀
                             && ode1.q₀ == ode2.q₀)

ext: jl


text:
using Test
using StaticArrays
using Nystrom
using LinearAlgebra
using Random
Random.seed!(1)

@testset "Discrete Operator test" begin
    T = SVector{3,ComplexF64}
    T1 = SMatrix{3,3,ComplexF64,9}
    T2 = SMatrix{2,3,ComplexF64,6}

    A = rand(T1, 5, 5)
    B = rand(T1, 5, 5)
    C = rand(T2, 5, 5)
    x = rand(T, 5)
    α = rand(ComplexF64)
    Aop = Nystrom.DiscreteOp(A)
    Bop = Nystrom.DiscreteOp(B)
    Cop = Nystrom.DiscreteOp(C)

    @testset "DiscreteOp test" begin
        @test size(A) == size(Aop)
        @test A*x == Aop*x
        @test A == Nystrom.materialize(Aop)
    end

    @testset "CompositeDiscreteOp test" begin
        @test size(C*A) == size(Cop*Aop)
        @test C*(A*x) == Cop*Aop*x == Cop*(Aop*x)
        @test C*A == Nystrom.materialize(Cop*Aop)

        @test size(C*B*A) == size(Cop*Bop*Aop)
        @test length((Cop*Bop*Aop).maps) == 3
        @test C*(B*(A*x)) == Cop*Bop*Aop*x == Cop*Bop*(Aop*x) == Cop*(Bop*Aop)*x
        @test C*(B*A) == Nystrom.materialize(Cop*Bop*Aop)
    end

    @testset "LinearCombinationDiscreteOp test" begin
        @test size(α*A - B) == size(α*Aop - Bop)
        @test (α*(A*x) - B*x) == (α*Aop*x - Bop*x) == (α*Aop - Bop)*x
        @test (α*A - B) == Nystrom.materialize(α*Aop - Bop)

        @test size(α*A - B + α*B) == size(α*Aop - Bop + α*Bop)
        @test (α*(A*x) - B*x + α*(B*x)) == (α*Aop - Bop + α*Bop)*x
        @test (α*A - B + α*B) == Nystrom.materialize(α*Aop - Bop + α*Bop)
    end

    @testset "Mixed tests" begin
        @test size(C*α*(α*I+B-α*A)+(1-α)*C) == size(Cop*α*(α*I+Bop-α*Aop)+(1-α)*Cop)
        @test C*(α*(α*x+B*x-α*(A*x)))+(1-α)*(C*x) == (Cop*α*(α*I+Bop-α*Aop)+(1-α)*Cop)*x
        @test C*(α*(α*I+B-α*A))+(1-α)*C == Nystrom.materialize(Cop*α*(α*I+Bop-α*Aop)+(1-α)*Cop)
    end

    @testset "GMRES and LU tests" begin
        # mesh
        Geometry.clear_entities!()
        Ω = ParametricSurfaces.Sphere(;radius=1) |> Geometry.Domain
        Γ = boundary(Ω)
        M = ParametricSurfaces.meshgen(Γ,(2,2))
        mesh = NystromMesh(view(M,Γ),order=1)

        mat = α*A-B+α*B
        op = α*Aop-Bop+α*Bop
        Amat = Nystrom.blockmatrix_to_matrix(mat)
        xvec = reinterpret(eltype(eltype(x)), x)
        b = reinterpret(eltype(x), Amat\xvec) |> collect
        σx = Density(x, mesh)
        σb = Density(b, mesh)

        ll = Nystrom.DiscreteOpGMRES(op, σx)
        yy = similar(xvec)
        mul!(yy, ll, xvec)

        @testset "Vector case" begin
            @test σb == (op\σx)
            @test σb ≈ Nystrom.gmres(op, σx)
        end

        op_scalar = Nystrom.DiscreteOp(Amat)
        σx_scalar = Density(collect(xvec), mesh)
        σb_scalar = Density(Amat\xvec, mesh)
        @testset "Scalar case" begin
            @test σb_scalar == (op_scalar\σx_scalar)
            @test σb_scalar ≈ Nystrom.gmres(op_scalar, σx_scalar)
        end
    end
end

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
module OmegaFlux

import ..Omega: ciid, logistic
import Flux
export OmegaDense

# "Dense network layer"
# function Flux.Dense(ω::Omega.Ω, in::Integer, out::Integer, σ = identity;
#                     initW = (ω, dims) -> logistic(ω, 0.0, 0.1, dims),
#                     initb = (ω, dims) -> logistic(ω, 0.0, 0.1, dims))
#   initW_ = initW(ω[@id][1], (out, in))
#   initb_ = initb(ω[@id][2], (out,))
#   Flux.Dense(initW_, initb_, σ)
# end

function OmegaDense(in, out, σ = identity;
  initW = logistic(0.0, 0.1, (in, out,)),
  initb = logistic(0.0, 0.1, (out,)))
  ciid(ω -> Flux.Dense(initW(ω), initb(ω), σ))
end

end


text:
using ReproductionNumbers

# setup
fig_name    =   "nowcasting-Germany"
file_name   =   "../data/raw/nowcasting-data.xlsx"
sheet_name  =   "Nowcast_R"
days_col    =   "Datum des Erkrankungsbeginns"
data_col    =   "Punktschätzer der Anzahl Neuerkrankungen (ohne Glättung)"
ylabel_R    =   "Instantaneous reproduction number R"
ylabel_N    =   "Nowcasted cases"

N, R, df = main_nowcasting(fig_name, file_name, sheet_name, days_col, data_col, ylabel_R, ylabel_N)


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    SDFileBond, SmilesBond, SmartsBond,
    setorder

struct SDFileBond <: Bond
    """Bond
    * Notation
        * Single bond
            * 0: u - v
            * 1: u ◀ v (Up-arrow)
            * 4: u ~ v (Up or down)
            * 6: u ◁ v (Down-arrow)
        * Double bond
            * 0: v ニ u (clockwise, default)
            * 1: u ニ v (counter-clockwise)
            * 2: u ＝ v (equal length, for terminal bond by default)
            * 3: u × v (Cis-Trans Unknown)
    """
    order::Int
    notation::Union{Int,Nothing}
end

SDFileBond() = SDFileBond(1, 0)
SDFileBond(order) = SDFileBond(order, 0)

setorder(edge::SDFileBond, order) = SDFileBond(order, edge.notation)

struct SmilesBond <: Bond
    order::Int
    isaromatic::Union{Bool, Nothing}
    cistrans::Union{Int, Nothing}
end

SmilesBond() = SmilesBond(1, false, nothing)
SmilesBond(order) = SmilesBond(order, false, nothing)

setorder(edge::SmilesBond, order
    ) = SmilesBond(order, edge.isaromatic, edge.cistrans)

struct SmartsBond <: QueryBond
    query::Pair
end

SmartsBond() = SmartsBond(:any => true)


text:
using PooksoftOptionsKit
using DataFrames
using CSV
using Dates

function load_actual_option_price_array(path_to_data_file::String, type::Symbol)::Array{Float64,2}

    # load -
    full_data_frame = CSV.read(path_to_data_file)

    # ok, so we need to extract the reqd cols -
    strike = Symbol("Strike")
    bid = Symbol("Bid")
    ask = Symbol("Ask")

    # find rows of type -
    idx_type = findall(x->x==String(type), full_data_frame[!,:Type])
    number_of_type_rows = length(idx_type)
    data_array = zeros(number_of_type_rows, 3)
    for (index, type_index) in enumerate(idx_type)

        data_array[index,1] = full_data_frame[type_index, strike]
        data_array[index,2] = full_data_frame[type_index, bid]
        data_array[index,3] = full_data_frame[type_index, ask]
    end

    # return -
    return data_array
end

function estimated_options_price_array()

    # setup the call test -
    stock_price = 54.17
    implied_volatility = 0.67
    DTE = (31.0/365.0)
    number_of_levels = 12
    strike_price_array = collect(range(45.0, step=1.0, stop=80.0))
    number_of_strike_prices = length(strike_price_array)
    risk_free_rate = 0.15 # this is in percent
    dividend_rate = 0.0 # this is in percent

    # initialize -
    estimated_price_array = zeros(number_of_strike_prices,2)

    # setup options calculation -
    amd_option_parameters = PSOptionKitPricingParameters(implied_volatility, DTE, number_of_levels, risk_free_rate, dividend_rate)

    # compute -
    for (index,strike_price) in enumerate(strike_price_array)

        # initialize -
        assetSet = Set{PSAbstractAsset}()

        # setup an amd call option -
        callOptionContract = PSCallOptionContract("AMD", Date(2020,6,25), strike_price, 0.0, 1; sense=:buy, contractMultiplier=1.0)
        push!(assetSet, callOptionContract)

        # build the pricing tree -
        tree = build_ternary_price_tree(assetSet, amd_option_parameters, stock_price)

        # compute -
        result = option_contract_price(tree, amd_option_parameters)

        # grab - we should check the return type ... but for now ...
        estimated_price_array[index,1] = strike_price
        estimated_price_array[index,2] = result.value
    end

    return estimated_price_array
end

# compute the estimated prices -
estimated_price_array = estimated_options_price_array()

# load the data -
path_to_data_file = "$(pwd())/tests/data/AMD-31D.csv"
actual_price_array = load_actual_option_price_array(path_to_data_file, :Call)


text:
import Base: +, -, *, /, ÷, sqrt, float
import LinearAlgebra: norm

_genid() = unsafe_trunc(UInt32, uuid4().value)
_indexid(id::UInt32, idx::Integer) = unsafe_trunc(UInt32, uuid5(UUID(id), "[$idx]").value)
_getidstr(x::AbstractBit) = digits(UInt8, x.id)
_getidstr(x::VecOrMat{<:AbstractBit}) = vcat(map(λ -> digits(UInt8, λ.id), x)...)
_getidstr(x) = Vector{UInt8}(string(x))

"""
    SBit

A stochastic bit is a pair of unipolar bits (positive and negative channels).

Fields:
- `bit::Tuple{Bool, Bool}`: a sample of a bitstream
- `value::Float64`: the underlying floating-point number being represented
- `id::UInt32`: a unique identifier for all samples of this bitstream
"""
struct SBit <: AbstractBit
    bit::Tuple{Bool, Bool}
    value::Float64
    id::UInt32

    function SBit(bit::Tuple{Bool, Bool}, value::Float64, id::UInt32)
        if value > 1 || value < -1
            @warn "SBitstream can only be ∈ [-1, 1] (saturation occurring)."
        end

        new(bit, min(max(value, -1), 1), id)
    end
end
function SBit(bits::VecOrMat{Tuple{Bool, Bool}}, values::VecOrMat{Float64}, id::UInt32)
    arr = similar(bits, SBit)
    for (i, bit) in enumerate(bits)
        arr[i] = SBit(bit, values[i], _indexid(id, i))
    end

    return arr
end

"""
    SBitstream

A stochastic bitstream that represents a real (floating-point) number
between [-1, 1].

Fields:
- `bits::Queue{SBit}`: the underlying bitstream
- `value::Float64`: the underlying floating-point number being represented
- `id::UInt32`: a unique identifier for this bitstream (set automatically)
"""
struct SBitstream <: AbstractBitstream
    bits::Queue{SBit}
    value::Float64
    id::UInt32

    function SBitstream(bits::Queue{SBit}, value::Float64, id::UInt32 = _genid())
        if value > 1 || value < -1
            @warn "SBitstream can only be ∈ [-1, 1] (saturation occurring)."
        end

        new(bits, min(max(value, -1), 1), id)
    end
end
SBitstream(value::Real, id::UInt32 = _genid()) = SBitstream(Queue{SBit}(), float(value), id)
function SBitstream(values::VecOrMat{Float64}, id::UInt32 = _genid())
    arr = similar(values, SBitstream)
    for (i, value) in enumerate(values)
        arr[i] = SBitstream(value, _indexid(id, i))
    end

    return arr
end

"""
    pos(b::SBit)

Return the positive channel bit of a stochastic bit.
"""
pos(b::SBit) = b.bit[1]

"""
    neg(b::SBit)

Return the negative channel bit of a stochastic bit.
"""
neg(b::SBit) = b.bit[2]

"""
    float(b::SBit)
    float(s::SBitstream)

Return the underlying floating-point value of
a stochastic bit or bitstream.
"""
float(b::SBit) = b.value
float(s::SBitstream) = s.value

include("./soperators.jl")
include("./simulatable.jl")

# SBit operator definitions #
@simulatable(SSignedAdder,      +(x::SBit, y::SBit) = x.value + y.value)
@simulatable(SSignedSubtractor, -(x::SBit, y::SBit) = x.value - y.value)
@simulatable(SSignedMultiplier, *(x::SBit, y::SBit) = x.value * y.value)
@simulatable(SSignedDivider,
function /(x::SBit, y::SBit)
    if y.value <= 0
        error("SBit only supports divisors > 0 (y == $y).")
    end

    x.value / y.value
end)
@simulatable(SSignedFixedGainDivider,
function ÷(x::SBit, y::Real)
    if y < 1
        error("SBit only supports fixed-gain divisors >= 1 (y == $y).")
    end

    x.value / y
end)
@simulatable(SSquareRoot,         sqrt(x::SBit) = sqrt(x.value))
@simulatable(SSignedDecorrelator, decorrelate(x::SBit) = x.value)
@simulatable(SSignedMatMultiplier,
    *(x::VecOrMat{SBit}, y::VecOrMat{SBit}) = float.(x) * float.(y), (size(x, 1), size(y, 2)))
@simulatable(SL2Normer,           norm(x::Vector{SBit}) = norm(float.(x)))

# SBitstream operator definitions #
for op in (:+, :-, :*, :/)
    @eval function $(op)(x::SBitstream, y::SBitstream)
        zbit = $(op)(pop!(x), pop!(y))
        z = SBitstream(zbit.value, zbit.id)
        push!(z, zbit)

        return z
    end
end
for op in (:sqrt, :decorrelate)
    @eval function $(op)(x::SBitstream)
        zbit = $(op)(pop!(x))
        z = SBitstream(zbit.value, zbit.id)
        push!(z, zbit)

        return z
    end
end
function ÷(x::SBitstream, y::Real)
    zbit = pop!(x) ÷ y
    z = SBitstream(zbit.value, zbit.id)
    push!(z, zbit)

    return z
end
function *(x::VecOrMat{SBitstream}, y::VecOrMat{SBitstream})
    zbit = pop!.(x) * pop!.(y)
    z = map(λ -> SBitstream(λ.value, λ.id), zbit)
    push!.(z, zbit)

    return z
end
function norm(x::Vector{SBitstream})
    zbit = norm(pop!.(x))
    z = SBitstream(zbit.value, zbit.id)
    push!(z, zbit)

    return z
end

"""
    generate(s::SBitstream, T::Integer = 1)
    generate!(s::SBitstream, T::Integer = 1)

Generate `T` samples of the bitstream.
Add them to its queue for `generate!`.
"""
function generate(s::SBitstream, T::Integer = 1)
    r = rand(T)
    bits = (s.value >= 0) ? zip(r .< abs(s.value), fill(false, T)) :
                            zip(fill(false, T), r .< abs(s.value))
    sbits = map(x -> SBit(x, s.value, s.id), bits)

    return sbits
end
generate!(s::SBitstream, T::Integer = 1) = push!(s, generate(s, T))

push!(s::SBitstream, b::SBit) = enqueue!(s.bits, SBit(b.bit, b.value, s.id))
pop!(s::SBitstream) = isempty(s.bits) ? generate(s)[1] : dequeue!(s.bits)

"""
    estimate!(buffer::AbstractVector, b::SBit)
    estimate!(buffer::AbstractVector, b::VecOrMat{SBit})
    estimate!(buffer::AbstractVector)

Push `b` into the `buffer` and return the current estimate.
"""
function estimate!(buffer::AbstractVector, b::SBit)
    push!(buffer, pos(b) - neg(b))

    return sum(buffer) / length(buffer)
end
function estimate!(buffer::AbstractVector, bs::VecOrMat{SBit})
    push!(buffer, pos.(bs) - neg.(bs))

    return sum(buffer) / length(buffer)
end
estimate!(buffer::AbstractVector) = sum(buffer) / length(buffer)


text:
my_f(x,y) = 2x+3y
derivative_of_my_f(x,y) = ForwardDiff.derivative(x -> my_f(x,y), x)


text:
using InteractiveUtils: versioninfo
versioninfo()

println("-------------------------------------------------------------------------")
println("Julia started from terminal without args: ", length(ARGS) == 0 ? "✔" : "✗")
println("Running from the same folder as this script: ", isfile("./test_setup.jl") ? "✔" : "✗")
println("Running Julia 1.6.0 or above: ", VERSION >= v"1.6.0" ? "✔" : "✗")

name = try
	readchomp(`git config user.name`)
catch
	""
end

mail = try
	readchomp(`git config user.email`)
catch
	""
end

println("Git Config Username: ", length(name) > 0 ? "✔" : "✗")
println("Git Config Email: ", length(mail) > 0 ? "✔" : "✗")

s = try
	using Pkg
	Pkg.activate(".")
	Pkg.add("BenchmarkTools")
	true
catch e
	@warn e
	false
end

println("Installed test environment for later use: ", s ? "✔" : "?")


text:
export DQNLearner

using Random
using Flux

mutable struct DQNLearner{
    Tq<:AbstractApproximator,
    Tt<:AbstractApproximator,
    Tf,
    S<:Union{Int,Nothing},
    R<:AbstractRNG,
} <: AbstractLearner
    approximator::Tq
    target_approximator::Tt
    loss_func::Tf
    stack_size::S
    γ::Float32
    batch_size::Int
    update_horizon::Int
    min_replay_history::Int
    update_freq::Int
    target_update_freq::Int
    update_step::Int
    rng::R
end

function DQNLearner(;
    approximator::Tq,
    target_approximator::Tt,
    loss_func::Tf,
    stack_size::Union{Int,Nothing} = 4,
    γ::Float32 = 0.99f0,
    batch_size::Int = 32,
    update_horizon::Int = 1,
    min_replay_history::Int = 32,
    update_freq::Int = 1,
    target_update_freq::Int = 100,
    update_step::Int = 0,
    seed = nothing,
) where {Tq,Tt,Tf}
    copyto!(approximator, target_approximator)
    rng = MersenneTwister(seed)
    DQNLearner(
        approximator,
        target_approximator,
        loss_func,
        stack_size,
        γ,
        batch_size,
        update_horizon,
        min_replay_history,
        update_freq,
        target_update_freq,
        update_step,
        rng,
    )
end

"""

!!! note
    The state of the observation is assumed to have been stacked,
    if `!isnothing(stack_size)`.
"""
(learner::DQNLearner)(obs) =
    obs |> get_state |>
    x ->
        send_to_device(device(learner.approximator), x) |> learner.approximator |>
        send_to_host

function RLBase.update!(learner::DQNLearner, batch)
    learner.update_step += 1
    learner.update_step % learner.update_freq == 0 || return nothing

    Q, Qₜ, γ, loss_func, update_horizon, batch_size = learner.approximator,
    learner.target_approximator,
    learner.γ,
    learner.loss_func,
    learner.update_horizon,
    learner.batch_size
    states, rewards, terminals, next_states = map(
        x -> send_to_device(device(Q), x),
        (batch.states, batch.rewards, batch.terminals, batch.next_states),
    )
    actions = CartesianIndex.(batch.actions, 1:batch_size)

    gs = gradient(params(Q)) do
        q = batch_estimate(Q, states)[actions]
        q′ = dropdims(maximum(batch_estimate(Qₜ, next_states); dims = 1), dims = 1)
        G = rewards .+ γ^update_horizon .* (1 .- terminals) .* q′
        loss_func(G, q)
    end

    update!(Q, gs)

    if learner.update_step % learner.target_update_freq == 0
        copyto!(Qₜ, Q)
    end
end

function RLBase.extract_experience(t::AbstractTrajectory, learner::DQNLearner)
    s = learner.stack_size
    h = learner.update_horizon
    n = learner.batch_size
    γ = learner.γ
    valid_ind_range = isnothing(s) ? (1:(length(t)-h)) : (s:(1:(length(t)-h)))
    if length(t) > learner.min_replay_history
        inds = rand(learner.rng, valid_ind_range, n)
        states = consecutive_view(get_trace(t, :state), inds; n_stack = s)
        actions = consecutive_view(get_trace(t, :action), inds)
        next_states = consecutive_view(get_trace(t, :state), inds .+ h; n_stack = s)
        consecutive_rewards = consecutive_view(get_trace(t, :reward), inds; n_horizon = h)
        consecutive_terminals =
            consecutive_view(get_trace(t, :terminal), inds; n_horizon = h)
        rewards, terminals = zeros(Float32, n), fill(false, n)

        # make sure that we only consider experiences in current episode
        for i in 1:n
            m = findfirst(view(consecutive_terminals, :, i))
            if isnothing(m)
                terminals[i] = false
                rewards[i] = discount_rewards_reduced(view(consecutive_rewards, :, i), γ)
            else
                terminals[i] = true
                rewards[i] = discount_rewards_reduced(view(consecutive_rewards, 1:m, i), γ)
            end
        end
        (
            states = states,
            actions = actions,
            rewards = rewards,
            terminals = terminals,
            next_states = next_states,
        )
    else
        nothing
    end
end


text:
"""
Produce an elliptic torus with elliptic tube cross
section and overall elliptic contour
"""
module DoubleEllipticTorus

# minimum distance between a point and ellipse (X-Y plane)
# note: the original algorithm only works for ellipse with
# major axis aligned to x-axis
# This is why I need to check whether a < b,
# and make a rotation by -90 degrees.
function min_distance_2d(a, b, coordinates)
    if a < b
        return min_distance_2d(b, a, [coordinates[2], -coordinates[1]])
    end
    x = abs(coordinates[1])
    y = abs(coordinates[2])
    function aux(phi, accum)
        new_phi = atan(((a^2 - b^2)*sin(phi) + y*b)/(x*a))
        if abs(phi - new_phi) < 1e-4
            return new_phi
        else
            return aux(new_phi, accum+1)
        end
    end
    phi = aux(0, 0)
    return sqrt((x - a * cos(phi))^2 + (y - b*sin(phi))^2)
end

function isInside(a, b, xy)
    x, y = xy
    return (x^2/(a^2) + y^2/(b^2) < 1.0)
end

# check whether a point is inside the torus

function get_fn(argd::Dict)
    major_a = float(argd["major"]["a"])
    major_b = float(argd["major"]["b"])
    minor_a = float(argd["minor"]["a"])
    minor_b = float(argd["minor"]["b"])
    center = [float(x) for x in argd["center"]]
    power = convert(Int, argd["power"])
    force_constant = float(argd["k"])

    function fn_compute(raw_coordinates)
        coordinates = raw_coordinates - center
        if norm(coordinates[1:2]) < 1e-6
            d_xy = min(major_a, major_b)
        else
            d_xy = min_distance_2d(major_a, major_b, coordinates[1:2])
        end
        if isInside(major_a, major_b, coordinates[1:2])
            new_coordinates = [-d_xy, coordinates[3]]
        else
            new_coordinates = [d_xy, coordinates[3]]
        end

        # min. distance to the minor ellipse
        if norm(new_coordinates) < 1e-6
            d = min(minor_a, minor_b)
        else
            d = min_distance_2d(minor_a, minor_b, new_coordinates)
        end

        if isInside(minor_a, minor_b, new_coordinates)
            return 0.0
        else
            return (1.0/power) * force_constant * (d^power)
        end
    end
    return fn_compute
end

export get_fn

end

max_stars_repo_name: yuhangwang/SimShape

max_issues_repo_name: yuhangwang/SimShape

max_forks_repo_name: yuhangwang/SimShape


text:
import Base: @deprecate

@deprecate DataFrame!(args...; kwargs...) DataFrame(args...; copycols=false, kwargs...)

# TODO: remove these definitions in year 2021
by(args...; kwargs...) = throw(ArgumentError("by function was removed from DataFrames.jl. " *
                                             "Use the `combine(groupby(...), ...)` or `combine(f, groupby(...))` instead."))

aggregate(args...; kwargs...) = throw(ArgumentError("aggregate function was removed from DataFrames.jl. " *
                                                    "Use the `combine` function instead."))

export categorical, categorical!
function CategoricalArrays.categorical(df::AbstractDataFrame,
                                       cols::Union{ColumnIndex, MultiColumnIndex};
                                       compress::Union{Bool, Nothing}=nothing)
    if compress === nothing
        compress = false
        categoricalstr = "categorical"
    else
        categoricalstr = "(x -> categorical(x, compress=$compress))"
    end
    if cols isa AbstractVector{<:Union{AbstractString, Symbol}}
        Base.depwarn("`categorical(df, cols)` is deprecated. " *
                     "Use `transform(df, cols .=> $categoricalstr, renamecols=false)` instead.",
                     :categorical)
        return transform(df, cols .=> (x -> categorical(x, compress=compress)), renamecols=false)
    elseif cols isa Union{AbstractString, Symbol}
        Base.depwarn("`categorical(df, cols)` is deprecated. " *
                     "Use `transform(df, cols => $categoricalstr, renamecols=false)` instead.",
                     :categorical)
        return transform(df, cols => (x -> categorical(x, compress=compress)), renamecols=false)
    else
        Base.depwarn("`categorical(df, cols)` is deprecated. " *
                     "Use `transform(df, names(df, cols) .=> $categoricalstr, renamecols=false)` instead.",
                     :categorical)
        return transform(df, names(df, cols) .=> (x -> categorical(x, compress=compress)), renamecols=false)
    end
end

function CategoricalArrays.categorical(df::AbstractDataFrame,
                                       cols::Union{Type, Nothing}=nothing;
                                       compress::Bool=false)
    if compress === nothing
        compress = false
        categoricalstr = "categorical"
    else
        categoricalstr = "categorical(x, compress=$compress)"
    end
    if cols === nothing
        cols = Union{AbstractString, Missing}
        Base.depwarn("`categorical(df)` is deprecated. " *
                     "Use `transform(df, names(df, $cols) .=> $categoricalstr, renamecols=false)` instead.",
                     :categorical)
    else
        Base.depwarn("`categorical(df, T)` is deprecated. " *
                     "Use transform(df, names(df, T) .=> $categoricalstr, renamecols=false)` instead.",
                     :categorical)
    end
    return transform(df, names(df, cols) .=> (x -> categorical(x, compress=compress)), renamecols=false)
end

function categorical!(df::DataFrame, cols::Union{ColumnIndex, MultiColumnIndex};
                      compress::Union{Bool, Nothing}=nothing)
    if compress === nothing
        compress = false
        categoricalstr = "categorical"
    else
        categoricalstr = "(x -> categorical(x, compress=$compress))"
    end
    if cols isa AbstractVector{<:Union{AbstractString, Symbol}}
        Base.depwarn("`categorical!(df, cols)` is deprecated. " *
                     "Use `transform!(df, cols .=> $categoricalstr, renamecols=false)` instead.",
                     :categorical!)
        return transform!(df, cols .=> (x -> categorical(x, compress=compress)), renamecols=false)
    elseif cols isa Union{AbstractString, Symbol}
        Base.depwarn("`categorical!(df, cols)` is deprecated. " *
                     "Use `transform!(df, cols => $categoricalstr, renamecols=false)` instead.",
                     :categorical!)
        return transform!(df, cols => (x -> categorical(x, compress=compress)), renamecols=false)
    else
        Base.depwarn("`categorical!(df, cols)` is deprecated. " *
                     "Use `transform!(df, names(df, cols) .=> $categoricalstr, renamecols=false)` instead.",
                     :categorical!)
        return transform!(df, names(df, cols) .=> (x -> categorical(x, compress=compress)), renamecols=false)
    end
end

function categorical!(df::DataFrame, cols::Union{Type, Nothing}=nothing;
                      compress::Bool=false)
    if compress === nothing
        compress = false
        categoricalstr = "categorical"
    else
        categoricalstr = "(x -> categorical(x, compress=$compress))"
    end
    if cols === nothing
        cols = Union{AbstractString, Missing}
        Base.depwarn("`categorical!(df)` is deprecated. " *
                     "Use `transform!(df, names(df, $cols) .=> $categoricalstr, renamecols=false)` instead.",
                     :categorical!)
    else
        Base.depwarn("`categorical!(df, T)` is deprecated. " *
                     "Use `transform!(df, names(df, T) .=> $categoricalstr, renamecols=false)` instead.",
                     :categorical!)
    end
    return transform!(df, names(df, cols) .=> (x -> categorical(x, compress=compress)), renamecols=false)
end


text:
export Element
"""
Subtypes of `Element` `show` FastHenry commands.

Geometric transformations can be preformed on elements with `transform`.
Elements which require a name will automatically generate unique names if
no name is provided.  Groups of elements are elements.
"""
abstract type Element end

# elements
include("node.jl")
include("segment.jl")
include("units.jl")
include("default.jl")
include("external.jl")
include("freq.jl")
include("equiv.jl")
include("end.jl")
include("uniformdiscretizedplanes.jl")
include("comment.jl")
include("group.jl")

include("context.jl")

# functions for elements
include("transform.jl")
include("traversetree.jl")
include("showelement.jl")
#include("mesh.jl")
include("plot.jl")


text:
using ..Continuations: _normalize!
using ..ArrayUtils: _eig

other(::HopfCont) = SaddleNodeCont()
other(::SaddleNodeCont) = HopfCont()

orthogonalize(v, u) = v .- u .* ((u ⋅ v) / (u ⋅ u))

function BifurcationProblem(point::AbstractSpecialPoint,
                            solver::Codim2Solver;
                            t_domain = solver.prob.t_domain,
                            param_axis1 = solver.prob.param_axis1,
                            param_axis2 = solver.prob.param_axis2,
                            kwargs...)

    cd2_prob = solver.prob :: DiffEqCodim2Problem
    de_prob = cd2_prob.de_prob :: AbstractODEProblem

    # Manually cast.  The fact that `resolved.u` is not of type xtype
    # indicates type instability somewhere.  Until it is fixed, let's
    # silently cast always.  After type stability is established, it
    # has to be changed to emit warning before cast.
    xtype = typeof(de_prob.u0)

    resolved = resolve_point(point, solver)
    N = length(de_prob.u0)
    x0 = xtype(resolved.u[1:N])
    t0 = SVector{2}(resolved.u[end-1:end])
    @assert param_axis1 == solver.prob.param_axis1
    @assert param_axis2 == solver.prob.param_axis2

    @assert point.point_type in (PointTypes.bogdanov_takens,
                                 PointTypes.fold_hopf)
    if (contkind(solver) :: ContinuationKind) isa SaddleNodeCont
        next_ckind = HopfCont()
    else
        next_ckind = SaddleNodeCont()
    end

    @assert timekind(point) isa Continuous

    v0, w0 = next_eig(next_ckind, resolved)
    @assert length(v0) == length(x0)
    if next_ckind isa HopfCont
        v0 = v0 .+ 0im
    else
        v0 = real.(v0)
    end
    v0 = cast_container(xtype, v0)

    # Sanity checks:
    if next_ckind isa SaddleNodeCont
        @assert eltype(v0) <: Real
    else
        @assert eltype(v0) <: Complex
    end
    if xtype <: SVector
        v0 :: SVector
    end

    return DiffEqCodim2Problem(
        de_prob,
        param_axis1,
        param_axis2,
        t_domain;
        x0 = x0,
        v0 = v0,
        w0 = w0,
        t0 = t0,
        augmented_system = preferred_augsys(next_ckind),
        kwargs...)
end

next_eig(next_ckind::ContinuationKind, point::SpecialPoint) =
    next_eig(next_ckind, Val(point.point_type), point)

# Switching to Hopf continuation via Bogdanov-Takens from saddle-node
function next_eig(::HopfCont, ::Val{PointTypes.bogdanov_takens},
                  point)
    prev_ckind = SaddleNodeCont()
    v_prev = ds_eigvec(prev_ckind, point.u)
    vals, vecs = _eig(ds_jacobian(prev_ckind, point.J))

    if eltype(vals) <: Complex
        _, i0 = findmin(abs.(real.(vals)))
        v0 = vecs[:, i0]
        w0 = imag(vals[i0])
        if w0 < 0
            w0 = -w0
            v0 = conj(v0)
        end
        return canonicalize(v0), w0
    end

    i1, i2 = sortperm(vals, by=abs)
    v1 = vecs[:, i1]
    v2 = vecs[:, i2]
    if abs(v1 ⋅ v_prev) > abs(v2 ⋅ v_prev)
        v_img = v1
    else
        v_img = v2
    end
    v_next = canonicalize(@. v_prev + v_img * im)
    return v_next, zero(real(eltype(v_prev)))
end

# Switching to saddle-node continuation via Bogdanov-Takens from Hopf
function next_eig(::SaddleNodeCont, ::Val{PointTypes.bogdanov_takens},
                  point)
    prev_ckind = HopfCont()
    v_prev = ds_eigvec(prev_ckind, point.u)
    v_next = _normalize!(real.(v_prev))
    return v_next, zero(eltype(v_next))
end

# Switching to Hopf (saddle-node) continuation via Fold-Hopf from
# saddle-node (Hopf)
function next_eig(next_ckind::ContinuationKind,
                  ::Val{PointTypes.fold_hopf},
                  point)
    prev_ckind = other(next_ckind)
    vals, vecs = _eig(ds_jacobian(prev_ckind, point.J))
    i1, i2, i3 = sortperm(vals, by=abs∘real)
    ir, _, ii = sort((i1, i2, i3), by=i -> abs(imag(vals[i])))
    if next_ckind isa HopfCont
        return _normalize!(vecs[:, ii]), imag(vals[ii])
    else
        return _normalize!(vecs[:, ir]), imag(vals[ir])
    end
end


text:
#! /usr/local/bin/julia
using ArgParse
using FASTX

s = ArgParseSettings()
@add_arg_table s begin
    "--input", "-i"
        arg_type = String
        required = true
        help = "Path to the input file (FASTA or tsv)."
    "--output", "-o"
        arg_type = String
        required = true
        help = "Path to the output file."
end

args = parse_args(ARGS, s)

IS_FASTA = occursin("fasta", args["input"])

if IS_FASTA
    records = [rec for rec in FASTA.Reader(open(args["input"]))]
    fasta_content = ["$(FASTA.identifier(rec))\t$(FASTA.sequence(rec))" for rec in records]
    open(args["output"], "w") do io
        write(io, join(fasta_content, "\n"))
    end
else
    proteins = [split(line, "\t") for line in readlines(args["input"])]
    open(args["output"], "w") do io
        for protein in proteins
            write(io, ">$(protein[1])\n$(protein[2])\n")
        end
    end
end

max_stars_repo_name: GreenCUBIC/PeptideScreening

max_issues_repo_name: GreenCUBIC/PeptideScreening

max_forks_repo_name: GreenCUBIC/PeptideScreening


text:
using PastaQ
using ITensors
using Test
using Random
import Optimisers

@testset "Optimizer update" begin
  θ = rand(10)
  ∇ = rand(10)
  θtest = copy(θ)
  opt = Optimisers.Descent(0.01)
  st  = Optimisers.state(opt, θ)
  st, θ = Optimisers.update(opt, st, θ, ∇)
  @test θ ≈ (θtest - 0.01*∇)
end

@testset "get/set parameters" begin
  N = 3
  sites = siteinds("Qubit",N)

  # MPS WAVEFUNCTION
  ψ = randomstate(sites; χ = 4)
  θ = PastaQ.getparameters(LPDO(ψ))
  ϕ = LPDO(randomstate(sites; χ = 4))
  @test !(ψ ≈ ϕ.X)
  PastaQ.setparameters!(ϕ, θ)
  @test PastaQ.array(ψ) ≈ PastaQ.array(ϕ.X)

  # LPDO (state)
  ρ = randomstate(sites; χ = 4, ξ = 3)
  θ = PastaQ.getparameters(ρ)
  σ = randomstate(sites; χ = 4, ξ = 3)
  PastaQ.setparameters!(σ, θ)
  @test PastaQ.array(ρ) ≈ PastaQ.array(σ)

  # UNITARY
  U = LPDO(PastaQ.unitary_mpo_to_choi_mps(randomprocess(sites; χ = 4)))
  θ = PastaQ.getparameters(U)
  V = LPDO(PastaQ.unitary_mpo_to_choi_mps(randomprocess(sites; χ = 4)))
  PastaQ.setparameters!(V, θ)
  @test PastaQ.array(U.X) ≈ PastaQ.array(V.X)

  # CHOI
  Λ = randomprocess(sites; χ = 4, ξ = 3)
  θ = PastaQ.getparameters(Λ)
  Γ = randomprocess(sites; χ = 4, ξ = 3)
  PastaQ.setparameters!(Γ, θ)
  @test PastaQ.array(Γ) ≈ PastaQ.array(Λ)
end

@testset "mps-qst: SGD update" begin
  N = 3
  χ = 10
  nsamples = 10
  data = PastaQ.convertdatapoints(randompreparations(N, nsamples))

  ψ = randomstate(N; χ=χ)
  opt = Optimisers.Descent(0.1)
  st = PastaQ.state(opt, ψ)
  ∇, _ = PastaQ.gradients(LPDO(ψ), data)

  ϕ = LPDO(copy(ψ))
  ϕ = PastaQ.update!(ϕ, ∇, (opt,st))
  ψp = copy(ψ)
  for j in 1:N
    ψp[j] = ψp[j] - 0.1*∇[j]
  end
  @test PastaQ.array(ϕ.X) ≈ PastaQ.array(ψp)
end

@testset "lpdo-qst: SGD update" begin
  N = 3
  χ = 10
  nsamples = 10
  data = PastaQ.convertdatapoints(randompreparations(N, nsamples))

  ρ = randomstate(N; χ=χ,ξ = 2)
  opt = Optimisers.Descent(0.1)
  st = PastaQ.state(opt, ρ)
  ∇, _ = PastaQ.gradients(ρ, data)

  γ = copy(ρ)
  γ = PastaQ.update!(γ, ∇, (opt,st))
  ρp = copy(ρ)
  for j in 1:N
    ρp.X[j] = ρp.X[j] - 0.1*∇[j]
  end
  @test PastaQ.array(γ) ≈ PastaQ.array(ρp)
end

@testset "mpo-qpt: SGD update" begin
  N = 3
  χ = 4
  nsamples = 10
  trace_preserving_regularizer = 0.1
  Random.seed!(1234)
  data_in = randompreparations(N, nsamples)
  data_out = PastaQ.convertdatapoints(randompreparations(N, nsamples))
  data = data_in .=> data_out

  U = randomprocess(N; χ=χ)
  Φ = LPDO(PastaQ.unitary_mpo_to_choi_mps(U))
  normalize!(Φ; localnorm=2)

  opt = Optimisers.Descent(0.1)
  st = PastaQ.state(opt, Φ)
  ∇, _ = PastaQ.gradients(Φ, data)

  γ = copy(Φ)
  γ = PastaQ.update!(γ, ∇, (opt,st))
  Φp = copy(Φ)
  for j in 1:N
    Φp.X[j] = Φp.X[j] - 0.1*∇[j]
  end
  @test PastaQ.array(γ.X) ≈ PastaQ.array(Φp.X)
end

@testset "lpdo-qpt: SGD update" begin
  N = 3
  χ = 4
  nsamples = 10
  trace_preserving_regularizer = 0.1
  Random.seed!(1234)
  data_in = randompreparations(N, nsamples)
  data_out = PastaQ.convertdatapoints(randompreparations(N, nsamples))
  data = data_in .=> data_out

  Λ = randomprocess(N; χ=χ, ξ = 3)
  normalize!(Λ; localnorm=2)

  opt = Optimisers.Descent(0.1)
  st = PastaQ.state(opt, Λ)

  ∇, _ = PastaQ.gradients(Λ, data)

  γ = copy(Λ)
  γ = PastaQ.update!(γ, ∇, (opt,st))
  Λp = copy(Λ)
  for j in 1:N
    Λp.X[j] = Λp.X[j] - 0.1*∇[j]
  end
  @test PastaQ.array(γ) ≈ PastaQ.array(Λp)
end


text:
using DataStructures: counter
module PermutationsA

export permlist, permmat, order, idpermlist, matrix, distance, permcycs, list
export cycles, cycletype, cyclelengths
export cprint, lprint, aprint, mprint
export pmap

using  PermPlain
import Base: length, sign, copy, convert, sparse, full, similar
import Base: show,print,writemime
import Base: kron
import PermPlain: permlist, permcycles, permsparse, permmatrix

include("abstractperm.jl")
include("permlist.jl")
include("permcycs.jl")
include("permmat.jl")
include("permsparse.jl")
include("permallrep.jl")
end # module PermutationsA


text:
random_seed = nothing

data_path = "save_load_test/"
data_prefix = "mnm_data"

#Model Parameters
iterations = 39
hard_clustering = false
argmax_sample_stop = 10 #Change to hard assignment from soft at iterations - argmax_sample_stop
split_stop = 10 #Stop split/merge moves at  iterations - split_stop

total_dim = 2

α = Float32(100000.0)

initial_clusters = 1

use_dict_for_global = false

hyper_params = multinomial_hyper(ones(Float32, 100))

#Saving specifics:
enable_saving = true
model_save_interval = 20
save_path = "save_load_test/"
overwrite_prec = false
save_file_prefix = "checkpoint"


text:
struct DimsSelector{N}
    dims::NTuple{N, Int}
end
dims(args...) = DimsSelector(args)

function (d::DimsSelector)(c::CartesianIndex{N}) where N
    t = ntuple(N) do n
        return n in d.dims ? c[n] : 1
    end
    return CartesianIndex(t)
end

getlabeledvector(data, shape, name::StringLike) = string(name), getcolumn(data, Symbol(name))

function getlabeledvector(data, shape, idx::Integer)
    name = columnnames(data)[idx]
    return getlabeledvector(data, shape, name)
end

function getlabeledvector(data, shape::NTuple{N, Any}, d::DimsSelector) where N
    sz = ntuple(N) do n
        return n in d.dims ? length(shape[n]) : 1
    end
    return "", reshape(CartesianIndices(sz), 1, sz...)
end

getlabeledvector(layer::Layer, name::Union{StringLike, Integer, DimsSelector}) =
    getlabeledvector(layer.data, shape(layer), name)

function getlabeledvector(layer::Layer, x::Pair{<:Any, <:StringLike})
    name, label = x
    _, v = getlabeledvector(layer, name)
    return label, v
end

function getlabeledvector(layer::Layer, x::Pair{<:Any, <:Any})
    name, transformation = x
    # consider supporting automated labeling for multiple names here
    label, v = getlabeledvector(layer, name)
    return label, map(transformation, v)
end

function getlabeledvector(layer::Layer, x::Pair{<:Any, <:Pair})
    name, transformation_label = x
    transformation, label = transformation_label
    names = name isa ArrayLike ? name : fill(name)
    vs = map(name -> last(getlabeledvector(layer, name)), names)
    v = map(transformation, vs...)
    return label, v
end

"""
    getlabeledarray(layer::Layer, s)

Return a label and an array from a selector `s`.
"""
getlabeledarray(layer::Layer, s) = getlabeledarray(layer, fill(s))

function getlabeledarray(layer::Layer, selectors::ArrayLike)
    labelsvectors = map(s -> getlabeledvector(layer, s), selectors)
    labels, vectors = map(first, labelsvectors), map(last, labelsvectors)
    arr = unnest(vectors)
    return labels, arr
end


text:
module Whittaker

export  whittaker1,
        whittaker1!,
        whittaker2,
        whittaker2!

include("whittaker1.jl")
include("whittaker2.jl")

end


text:
using sparsePolyChaos
using PolyChaos

###
### Multivariate sparse PCE
###
maxDeg = 10
interactions = 3

model(x,y) = (x^4 - 6 * x^2  + 3) + (x^2 - 1) * (y^2 - 1) + (y^3 - 3 * y) # He_4(x) + He_2(x) * He_2(y) + He_3(y)
# model(x,y) = (x^2 - 1) * (y^2 - 1) # He2(x) * He2(y)
# model(x,y) = (x^2 - 1) # He2(x)

### Setup and compute PCE coefficients of x ###
k = 2
op = GaussOrthoPoly(maxDeg)
mop = MultiOrthoPoly([op for i in 1:k], maxDeg)
L = dim(mop)

pce, Ap, p, R², Q² = sparsePCE(mop, model; pMax = maxDeg, jMax = interactions)

println("Basis polynomials: ", Ap)
println("PCE coefficients: ", pce)
println("Max degree:: ", p)
println("R² error: ", R²)
println("Q² error: ", Q²)

max_stars_repo_name: adriangrupp/sparsePolyChaos

max_issues_repo_name: adriangrupp/sparsePolyChaos

max_forks_repo_name: adriangrupp/sparsePolyChaos


text:
module EMLvol_tank
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...types.EMLtypes
	import EMLtypes.length
	include("vol_tank/vol_tank.jl")
end


text:
function myreverse(phrase::String)
    chars = Char[]
    for c in phrase
        push!(chars, c)
    end
    join(chars[end:-1:1])
end


text:
using Test
using CalibrateEDMF.ReferenceModels
using CalibrateEDMF.TurbulenceConvectionUtils
using Random

pwdir = mktempdir()
@testset "ReferenceModel" begin
    les_dir_test = joinpath(pwdir, "foo/bar/les")
    scm_dir_test = joinpath(pwdir, "foo/bar/scm/Output.DYCOMS_RF01.12345")
    case_name_test = "DYCOMS_RF01"
    y_names = ["thetal_mean", "ql_mean", "qt_mean"]
    ti = 0.0
    tf = 10.0

    ref_model = ReferenceModel(y_names, les_dir_test, scm_dir_test, case_name_test, ti, tf)

    @test y_dir(ref_model) == joinpath(pwdir, "foo/bar/les")
    @test Σ_dir(ref_model) == joinpath(pwdir, "foo/bar/les")
    @test scm_dir(ref_model) == joinpath(pwdir, "foo/bar/scm/Output.DYCOMS_RF01.12345")
    @test num_vars(ref_model) == 3
    @test uuid(ref_model) == "12345"
    @test get_t_start(ref_model) == 0.0
    @test get_t_start_Σ(ref_model) == 0.0
    @test get_t_end(ref_model) == 10.0
    @test get_t_end_Σ(ref_model) == 10.0
    z_obs = get_z_obs(ref_model)
    @test isa(z_obs, Array)
    # Check monotonicity
    for (i, height) in enumerate(z_obs[1:(end - 1)])
        @test height < z_obs[i + 1]
    end

    # Test coarse model
    n_obs = 20
    ref_model_coarse = ReferenceModel(y_names, les_dir_test, scm_dir_test, case_name_test, ti, tf, n_obs = n_obs)
    z_coarse = get_z_obs(ref_model_coarse)
    @test isa(z_coarse, Array)
    @test length(z_coarse) == n_obs

    # Test stretched grid
    namelist_args = [("grid", "stretch", "flag", true)]
    ref_model_stretched =
        ReferenceModel(y_names, les_dir_test, scm_dir_test, case_name_test, ti, tf, namelist_args = namelist_args)
    z_stretch = get_z_obs(ref_model_stretched)
    @test isa(z_stretch, Array)
    # Test stretching
    @test z_stretch[2] - z_stretch[1] < z_stretch[3] - z_stretch[2]
end

@testset "ReferenceModelBatch" begin
    scm_dir_test = joinpath(pwdir, "foo/bar/scm/Output.test_case.12345")
    # Choose same SCM to speed computation
    scm_dirs = repeat([scm_dir_test], 2)
    kwargs_ref_model = Dict(
        :y_names => [["u_mean"], ["v_mean", "thetal_mean"]],
        :y_dir => scm_dirs,
        :scm_dir => scm_dirs,
        :case_name => ["DYCOMS_RF01", "GABLS"],
        :t_start => repeat([0.0], 2),
        :t_end => repeat([10.0], 2),
    )
    ref_model_batch = ReferenceModelBatch(kwargs_ref_model)
    @test length(ref_model_batch.ref_models) == 2
    @test length(ref_model_batch.eval_order) == 2

    # Get a minibatch and check that the eval order shrinks
    ref_models, model_indices = get_minibatch!(ref_model_batch, 1)

    @test isa(ref_models, Vector{ReferenceModel})
    @test isa(model_indices, Vector{Int})
    @test length(ref_model_batch.ref_models) == 2
    @test length(ref_model_batch.eval_order) == 1
end

@testset "ReferenceModel handlers" begin
    # Choose same SCM to speed computation
    data_dir = mktempdir()
    scm_dirs = repeat([joinpath(data_dir, "Output.Bomex.000000")], 2)
    # Reduce resolution and t_max to speed computation as well
    t_max = 4 * 3600.0
    namelist_args = [
        ("time_stepping", "t_max", t_max),
        ("time_stepping", "dt_max", 30.0),
        ("time_stepping", "dt_min", 20.0),
        ("stats_io", "frequency", 720.0),
        ("grid", "dz", 150.0),
        ("grid", "nz", 20),
    ]

    kwargs_ref_model = Dict(
        :y_names => [["u_mean", "v_mean"], ["thetal_mean"]],
        :y_dir => scm_dirs,
        :scm_dir => scm_dirs,
        :case_name => repeat(["Bomex"], 2),
        :t_start => repeat([t_max - 3600], 2),
        :t_end => [t_max, t_max],
        :Σ_t_start => repeat([t_max - 2.0 * 3600], 2),
        :Σ_t_end => repeat([t_max - 0.5 * 3600], 2),
        :n_obs => [nothing, 15],
        :namelist_args => repeat([namelist_args], 2),
    )
    ref_models = construct_reference_models(kwargs_ref_model)
    run_reference_SCM.(ref_models, overwrite = true, run_single_timestep = false)
    Δt = 3.0 * 3600
    ref_model = time_shift_reference_model(ref_models[1], Δt)
    z = get_z_obs(ref_model)

    @test get_t_start(ref_model) == t_max - Δt + (t_max - 3600)
    @test get_t_end(ref_model) == t_max - Δt + t_max
    @test get_t_start_Σ(ref_model) == t_max - Δt + (t_max - 2.0 * 3600)
    @test get_t_end_Σ(ref_model) == t_max - Δt + (t_max - 0.5 * 3600)
    @test length(z) == 20
    @test z[2] - z[1] ≈ 150

    Δt = 2 * t_max
    @test_throws AssertionError time_shift_reference_model(ref_models[1], Δt)

    ref_model = ref_models[2]
    z = get_z_obs(ref_model)

    @test get_t_start(ref_model) == t_max - 3600
    @test get_t_end(ref_model) == t_max
    @test get_t_start_Σ(ref_model) == t_max - 2.0 * 3600
    @test get_t_end_Σ(ref_model) == t_max - 0.5 * 3600
    @test length(z) == 15
    @test z[2] - z[1] ≈ 150.0 * (20 - 1) / (15 - 1)

end


text:
# Science...

# Load modules
using Plots, LaTeXStrings
pyplot()
using ThinFilmsTools

# Substrate index of refraction taken from www.refractiveindex.info
function nacl(x::Float64)
    x /= 1000.0
    n = sqrt(complex(
            1.0 + 0.00055+0.19800/(1.0 - (0.050/x)^2) +
            0.48398/(1.0 - (0.100/x)^2) +
            0.38696/(1.0 - (0.128/x)^2) +
            0.25998/(1.0 - (0.158/x)^2) +
            0.08796/(1.0 - (40.50/x)^2) +
            3.17064/(1.0 - (60.98/x)^2) +
            0.30038/(1.0 - (120.34/x)^2))
    )
    return n
end

function plotSpec(sol, idx)
    plt = plot(Spectrum1D(),
            sol.beam.λ, [sol.Spectra.Rp[:,1,:], sol.Spectra.Rs[:,idx,:]],
            label=["p-wave" "s-wave"],
            yaxis=("Reflectance", (0., 1.)),
            xlims=(sol.beam.λ[1], sol.beam.λ[end]),
            title=("θ = $(sol.beam.θ[idx]) [°]"),
            line=([:solid :dashdot]));
    return plt
end

function main()
    # Define beam
    λ = LinRange(4900,15000,5000) # wavelength range [nm]
    θ = LinRange(0,80,80) # angle of incidence [degrees]
    beam = PlaneWave(λ,θ)
    # Define layers
    l0 = LayerTMMO(RIdb.air(beam.λ))
    l1 = LayerTMMO(RIdb.dummy(beam.λ,4.6,0.0); d=800.)
    l2 = LayerTMMO(RIdb.dummy(beam.λ,1.6,0.0); d=1650.)
    l3 = LayerTMMO(nacl.(beam.λ)) # notice the . in nacl to use the array input
    layers = [l0, l1, l2, l1, l2, l1, l2, l1, l2, l1, l3]
    # Reference wavelength
    λ0 = 12000.0
    # call main script
    sol = tmm_optics(beam, layers; λ0=λ0, pbgflag=true)
    return sol
end

sol = main()

# Plot the R, T and A spectra, at 0, 45 and 80 degress
t45 = Utils.find_closest(sol.beam.θ,45.0)
t45 = Utils.find_closest(sol.beam.θ,80.0)
plot(plotSpec(sol, 1), plotSpec(sol, t45), plotSpec(sol, t80), layout=(3,1))
gui()

# Plot EMF
p1 = plot(Spectrum2D(),
     sol.beam.λ, sol.beam.θ, sol.Spectra.Rs, title=("s-wave Reflectance"),
);
p2 = plot(Spectrum2D(),
     sol.beam.λ, sol.beam.θ, sol.Spectra.Rp, title=("p-wave Reflectance"),
)
plot(p1, p2, layout=(2,1))
gui()

# Plot the refractive index profile
plot(RIprofile(), sol)
gui()

# plot the photonic dispersion with custom function
plot(PBGDispersion2Dalt(), sol.Bloch)
gui()


text:
@testset "ReadStatTable" begin
    tb = ReadStatTable(AbstractVector[], Symbol[])
    @test size(tb) == (0, 0)
    @test length(tb) == 0
    @test isempty(tb)
    @test sprint(show, MIME("text/plain"), tb) == "0×0 ReadStatTable"

    tb = ReadStatTable(AbstractVector[[]], Symbol[:c])
    @test size(tb) == (0, 1)
    @test length(tb) == 1
    @test isempty(tb)
    @test sprint(show, MIME("text/plain"), tb) == "0×1 ReadStatTable"

    c1 = collect(1:10)
    c2 = collect(10.0:-1.0:1.0)
    cols = AbstractVector[c1, c2]
    names = [:c1, :c2]
    tb = ReadStatTable(cols, names)
    @test size(tb) == (10, 2)
    @test size(tb, 1) == 10
    @test_throws ArgumentError size(tb, 3)
    @test length(tb) == 2
    @test !isempty(tb)

    @test Tables.getcolumn(tb, 1) === c1
    @test Tables.getcolumn(tb, :c2) === c2
    @test columnnames(tb) == names
    @test columnnames(tb) !== names

    @test Tables.schema(tb) == Tables.Schema{(:c1, :c2), Tuple{Int, Float64}}()
    @test Tables.columnindex(tb, :c1) == 1
    @test Tables.columntype(tb, :c1) == Int
    @test Tables.rowcount(tb) == 10

    @test values(tb) === cols
    @test haskey(tb, :c1)
    @test haskey(tb, 2)

    @test sprint(show, tb) == "10×2 ReadStatTable"
    @test sprint(show, MIME("text/plain"), tb, context=:displaysize=>(15,80)) == """
        10×2 ReadStatTable:
         Row │    c1       c2
             │ Int64  Float64
        ─────┼────────────────
           1 │     1     10.0
           2 │     2      9.0
           3 │     3      8.0
          ⋮  │   ⋮       ⋮
           8 │     8      3.0
           9 │     9      2.0
          10 │    10      1.0
                4 rows omitted"""

    @test varlabels(tb) === nothing
    @test varformats(tb) === nothing
    @test val_label_keys(tb) === nothing
    @test val_label_dict(tb) === nothing
    @test filelabel(tb) === nothing
    @test filetimestamp(tb) === nothing
    @test fileext(tb) === nothing
end


text:
#patiencesort.jl
#Algorithm 2.5 of Random Eigenvalues by Alan Edelman

function patiencesort(p)
    piles = similar(p, 0)
    for pi in p
        idx = 1
        for pp in piles
            if pi > pp
                idx += 1
            end
        end
        if idx > length(piles)
            d = length(piles)
            resize!(piles, idx)
            piles[d+1:end] = 0
        end
        piles[idx] = pi
    end
    return length(piles)
end


text:
using Caesar

fg = initfg()

cvf = VelPoint2VelPoint2(MvNormal([0.0;0;0;0], Matrix(Diagonal([7.0;7.0;0.1;0.1].^2))))
# cvf = DynPoint2DynPoint2(MvNormal([0.0;0;0;0], Matrix(Diagonal([7.0;7.0;0.1;0.1].^2))))
br = Pose2Point2BearingRange(Normal(pi/2, 0.05), Normal(50.0, 1.0))

addVariable!(fg, :x1, Pose2)
addFactor!(fg, [:x1], PriorPose2(MvNormal([0.0;0;0],Matrix(Diagonal([0.1;0.1;0.05].^2)))))
addVariable!(fg, :l1, DynPoint2(ut=0))
addFactor!(fg, [:x1;:l1], br)

addVariable!(fg, :x2, Pose2)
addFactor!(fg, [:x2], PriorPose2(MvNormal([10.0;0;0],Matrix(Diagonal([0.1;0.1;0.05].^2)))))
addVariable!(fg, :l2, DynPoint2(ut=1_000_000))
addFactor!(fg, [:l1;:l2], cvf)
addFactor!(fg, [:x2;:l2], br)

addVariable!(fg, :x3, Pose2)
addFactor!(fg, [:x3], PriorPose2(MvNormal([20.0;0;0],Matrix(Diagonal([0.1;0.1;0.05].^2)))))
addVariable!(fg, :l3, DynPoint2(ut=2_000_000))
addFactor!(fg, [:l2;:l3], cvf)
addFactor!(fg, [:x3;:l3], br)

addVariable!(fg, :x4, Pose2)
addFactor!(fg, [:x4], PriorPose2(MvNormal([30.0;0;0],Matrix(Diagonal([0.1;0.1;0.05].^2)))))
addVariable!(fg, :l4, DynPoint2(ut=3_000_000))
addFactor!(fg, [:l3;:l4], cvf)
addFactor!(fg, [:x4;:l4], br)

# specific parameters
getSolverParams(fg).drawtree = true

# and solve
tree, smt, hist = solveTree!(fg)

## and plot
using RoMEPlotting
Gadfly.set_default_plot_size(40cm,25cm)

plotKDE(fg, sort(union(ls(fg,r"x"),ls(fg,r"l"))), dims=[1;2], title="Positions")

plotKDE(fg, sort(ls(fg,r"l")), dims=[3;4], title="velocities")

fg = initfg()

cvf = DynPoint2DynPoint2(MvNormal([0.0;0;0;0], Matrix(Diagonal([7.0;7.0;0.1;0.1].^2))))
br = Pose2Point2BearingRange(Normal(-pi/4, 0.05), Normal(20.0, 1.0))

addVariable!(fg, :x1, Pose2)
addFactor!(fg, [:x1], PriorPose2(MvNormal([0.0;0;0],Matrix(Diagonal([0.1;0.1;0.05].^2)))))
addVariable!(fg, :l1, DynPoint2(ut=0))
addFactor!(fg, [:x1;:l1], br)

addVariable!(fg, :x2, Pose2)
addFactor!(fg, [:x2], PriorPose2(MvNormal([10.0;0;0],Matrix(Diagonal([0.1;0.1;0.05].^2)))))
addVariable!(fg, :l2, DynPoint2(ut=1_000_000))
addFactor!(fg, [:l1;:l2], cvf)
addFactor!(fg, [:x2;:l2], br)

addVariable!(fg, :x3, Pose2)
addFactor!(fg, [:x3], PriorPose2(MvNormal([20.0;0;0],Matrix(Diagonal([0.1;0.1;0.05].^2)))))
addVariable!(fg, :l3, DynPoint2(ut=2_000_000))
addFactor!(fg, [:l2;:l3], cvf)
addFactor!(fg, [:x3;:l3], br)

addVariable!(fg, :x4, Pose2)
addFactor!(fg, [:x4], PriorPose2(MvNormal([30.0;0;0],Matrix(Diagonal([0.1;0.1;0.05].^2)))))
addVariable!(fg, :l4, DynPoint2(ut=3_000_000))
addFactor!(fg, [:l3;:l4], cvf)
addFactor!(fg, [:x4;:l4], br)

# specific parameters
getSolverParams(fg).drawtree = true

# and solve
tree, smt, hist = solveTree!(fg)

## and plot
using RoMEPlotting
Gadfly.set_default_plot_size(40cm,25cm)

plotKDE(fg, sort(union(ls(fg,r"x"),ls(fg,r"l"))), dims=[1;2], title="Positions")

plotKDE(fg, sort(ls(fg,r"l")), dims=[3;4], title="velocities")

#
#
# ## Test VelPoint2VelPoint2 factor
#
# res = [0.0;]
#
# userdata = FactorMetadata()
# userdata.variableuserdata = (DynPoint2(ut=0),DynPoint2(ut=1000000))
# idx = 1
# # meas = getSample(cvf)
# meas = (reshape([10;0;0;0.0],4,1),)
#
# Xi = zeros(4,1)
# Xj = zeros(4,1)
# Xj[1,1] = 10.0
#
# cost = (vx,vy) -> cvf(res,userdata,idx,meas,reshape([vy;Xi[2,1];vx;Xi[4,1]],4,1), Xj )
#
#
# cost(0.0,0.0)
# cost(5.0,0.0)
# cost(10.0,0.0)
#
#
#
# ## and look at the cost function
#
# using Gadfly
#
# vx=-20:0.1:20
# vy=-20:0.1:20
#
# Gadfly.plot(x=vx,y=vy,z=cost, Geom.contour)
#
#
#
# 0
#
# # new factor

#
#
# """
#     $(TYPEDEF)
#
# Partial prior belief on Z, Roll, and Pitch of a `Pose3`.
# """
# mutable struct PriorPose3ZRP{T1,T2} <: IncrementalInference.FunctorSingleton where {T1 <: SamplableBelief, T2 <: SamplableBelief}
#   z::T1
#   rp::T2
#   partial::Tuple
#   PriorPose3ZRP{T1,T2}() where {T1, T2} = new{T1,T2}()
#   PriorPose3ZRP{T1,T2}(z::T1,rp::T2) where {T1 <: IIF.SamplableBelief, T2 <: IIF.SamplableBelief} = new{T1,T2}(z, rp, (3,4,5))
# end
# PriorPose3ZRP(z::T1,rp::T2) where {T1 <: SamplableBelief, T2 <: SamplableBelief} = PriorPose3ZRP{T1,T2}(z, rp)
# function getSample(pprz::PriorPose3ZRP, N::Int=1)
#   return ([rand(pprz.z,N)[:]';rand(pprz.rp,N)], )
# end
#
#

#


text:
"""

## Arguments
- **`pt::Vector{Float64}`**: The point
- **`toprow_simplex::Array{Float64, 2}`**: A matrix representation of the simplex where each
    column represents a vertex, but with a row of ones appended on top of it, so that the
    total dimension is `(dim + 1)-by-(dim + 1)`. This is needed when computing determinants.
- **`dim`**: The dimension of the space.
"""
function iscontained_superoptim!(pt, toprow_simplex, dim, signs, startstops)

    # Replace first vertex with the point
    #fill_vert

    found = true
end

"""
Fill the 2nd to last row in the i-th column of `prealloc_matrix` with the
2nd to last row in the `i`-th column of `toprow_simplex` (replacing the `i`th vertex).
"""
function fill_vertex(prealloc_matrix::Array{Float64, 2}, toprow_simplex::Array{Float64, 2},
        startstops::Vector{UnitRange{Int64}}, i::Int)

    prealloc_matrix[startstops[i]] .= toprow_simplex[startstops[i]]
end

function contained_superoptim!(signs, s_arr, sx, point, dim, temp_arr, starts, stops)
    # Redefine the temporary simplex. This is in-place, so we don't allocate
    # memory. We could also have re-initialised `signs`, but since we're never
    # comparing more than two consecutive signs, this is not necessary.
    rezero!(s_arr)
    rezero!(signs)
    fill_into!(s_arr, sx)

    #Replace first vertex with the point
    fill_at_inds!(s_arr, point, 1:dim)

    # Signed volume
    signs[1] = sign(mdet_fm_in_x_view!(temp_arr, s_arr, dim + 1, starts, stops))

    rezero!(s_arr)
    fill_into!(s_arr, sx) #reset

    for κ = 2:dim # Check remaining signs and stop if sign changes
        # Replace the ith vertex with the point we're cheking (leaving the
        # 1 appended to Vi intact.)
        idxs = ((dim + 1)*(κ - 1)+1):((dim + 1)*(κ - 1)+ 1 + dim - 1)
        fill_at_inds!(s_arr, point, idxs) # ith change

        signs[κ] = sign(mdet_fm_in_x_view!(temp_arr, s_arr, dim + 1, starts, stops))

        if !(signs[κ-1] == signs[κ])
            return false
        end

        rezero!(s_arr)
        fill_into!(s_arr, sx)
    end

    # Last the last vertex with the point in question
    idxs = ((dim + 1)*(dim)+1):((dim+1)^2-1)
    fill_at_inds!(s_arr, point, idxs)

    signs[end] = sign(mdet_fm_in_x_view!(temp_arr, s_arr, dim + 1, starts, stops))

    if !(signs[end-1] == signs[end])
       return false
    else
        return true
    end
end

function innerloop_superoptim!(inds::Vector{Int}, signs, s_arr, Sj, pt, dim::Int, M, i::Int, temp_arr, starts, stops)
    for j in 1:length(inds)
        if contained_optim!(signs, s_arr, Sj[j], pt, dim, temp_arr, starts, stops)
            M[inds[j], i] += 1.0
        end
    end
end


text:
# module Bukdu

import Base: pipeline

"""
    pipeline(block::Function, routers...)
"""
function pipeline(block::Function, pipes...)
    for pipe::Symbol in pipes
        pipelines = get(Routing.routing_pipelines, pipe, [])
        push!(pipelines, block)
        Routing.routing_pipelines[pipe] = pipelines
    end
end

# module Bukdu


text:
# Recreate:
# ```R
# angle <- seq(0, 2*pi, length=100)[-11]
# grid.polygon(x=0.25 + 0.15*cos(angle), y=0.5 + 0.3*sin(angle),
#              id=rep(1:2, c(7, 3)), gp=gpar(fill=c("grey", "white")))
# angle[4] <- NA
# grid.polygon(x=0.75 + 0.15*cos(angle), y=0.5 + 0.3*sin(angle),
#              id=rep(1:2, c(7, 3)), gp=gpar(fill=c("grey", "white")))
# ```
# Again this has a lot of behavior I don't like in R, implicit recyclign of the fill in the
# second example. Using id to slice the points. The magical use of NA, to make three
# polygons. And the vectorization of the primatives. But I want to be able to make the same
# kind of figure
using OdinSon
using PyPlot

angle = linspace(0, 2*pi, 11)[1:(end - 1)]
#TODO: I should support giving two arrays instead of having to hcat them
points1 = hcat(0.25 + 0.15*cos(angle), 0.5 + 0.3*sin(angle))
p11 = Polygon(points1[1:(end - 3), :], style=Style(fill=NC"gray"))
p12 = Polygon(points1[(end - 2):end, :], style=Style(fill=NC"white"))
points2 = hcat(0.75 + 0.15*cos(angle), 0.5 + 0.3*sin(angle))
p21 = Polygon(points2[1:3, :], style=Style(fill=NC"gray"))
p22 = Polygon(points2[5:7, :], style=Style(fill=NC"gray"))
p23 = Polygon(points2[8:end, :], style=Style(fill=NC"white"))

render(Canvas([p11, p12, p21, p22, p23]))
#TODO: let this be set in the Canvas
gca()[:set_aspect](0.5)

max_stars_repo_name: gabrielgellner/OdinSon

max_issues_repo_name: gabrielgellner/OdinSon

max_forks_repo_name: gabrielgellner/OdinSon


text:
#=

This file contains all types that are used in the Absorption module:

- `HitranTable` stores HITRAN transition line parameters
- `AbstractBroadeningFunctions` specify the type of line-shape broadening to perform
- `AbstractComplexErrorFunctions` specify the error function to use for Voigt calculations
- `AbstractCrossSectionModel` specify whether to calculate the line-shape from scratch or
  using an interpolator.

=#

"""
    struct HitranTable{FT}

A struct, which provides all HITRAN line parameters needed to compute
absorption cross sections

See https://hitran.org/docs/definitions-and-units/ for details

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct HitranTable{FT<:AbstractFloat}
    "The molecular species identification (ID) number"
    mol::Array{Int,1}
    "The isotopologue ID number"
    iso::Array{Int,1}
    "The wavenumber of the spectral line transition (cm-1) in vacuum"
    νᵢ::Array{FT,1}
    "The spectral line intensity (cm−1/(molecule·cm−2)) at Tref=296K"
    Sᵢ::Array{FT,1}
    "The Einstein-A coefficient (s-1) of a transition"
    Aᵢ::Array{FT,1}
    "The air-broadened half width at half maximum (HWHM) (cm−1/atm) at Tref=296K and reference pressure pref=1atm"
    γ_air::Array{FT,1}
    "The self-broadened half width at half maximum (HWHM) (cm−1/atm) at Tref=296K and reference pressure pref=1atm"
    γ_self::Array{FT,1}
    "The lower-state energy of the transition (cm-1)"
    E″::Array{FT,1}
    "The coefficient of the temperature dependence of the air-broadened half width"
    n_air::Array{FT,1}
    "The pressure shift (cm−1/atm) at Tref=296K and pref=1atm of the line position with respect to the vacuum transition wavenumber νij"
    δ_air::Array{FT,1}
    "The electronic and vibrational quantum numbers and labels for the upper state of a transition"
    global_upper_quanta::Array{String,1}
    "The electronic and vibrational quantum numbers and labels for the lower state of a transition"
    global_lower_quanta::Array{String,1}
    "Rotational, hyperfine and other quantum numbers and labels for the upper state of a transition"
    local_upper_quanta::Array{String,1}
    "Rotational, hyperfine and other quantum numbers and labels for the lower state of a transition"
    local_lower_quanta::Array{String,1}
    "Ordered list of indices corresponding to uncertainty estimates of transition parameters"
    ierr::Array{String,1}
    "Ordered list of reference identifiers for transition parameters"
    iref::Array{String,1}
    "A flag indicating the presence of additional data and code relating to line-mixing"
    line_mixing_flag::Array{String,1}
    "The upper state degeneracy"
    g′::Array{FT,1}
    "The lower state degeneracy"
    g″::Array{FT,1}
end

"""
    struct AbscoTable{FT}

A struct, which provides all ABSCO tabulated cross section parameters

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct AbscoTable{FT<:AbstractFloat}
    "The molecular species identification (ID) number"
    mol::Int
    "The isotopologue ID number"
    iso::Int
    "The wavenumber array (cm⁻¹)"
    ν::Array{Float64,1}
    "The 4D cross section table"
    σ::Array{FT,4}
    "The pressure coordinates (hPa)"
    p::Array{FT,1}
    "The temperature coordinates (K)"
    T::Array{FT,2}
end

#=

Types of broadening functions (currently Doppler, Lorentz, and Voigt)

=#

"""
    type AbstractBroadeningFunction
Abstract Broadening Function type for generic line broadening function
"""
abstract type AbstractBroadeningFunction end

"Doppler line broadening"
struct Doppler <: AbstractBroadeningFunction end

"Lorentz line broadening"
struct Lorentz <: AbstractBroadeningFunction end

"Voigt line broadening"
struct Voigt <: AbstractBroadeningFunction end

#=

Types of complex error functions for use in Voigt computations

=#

"""
    type AbstractComplexErrorFunction
Abstract Complex Error Function type for generic complex error functions
"""
abstract type AbstractComplexErrorFunction end

"Humlicek only formulation for Complex Error Function"
struct HumlicekErrorFunction <: AbstractComplexErrorFunction end

"Mix of Humlicek and Weidemann (N=32) Error Function, suggested for Voigt function"
struct HumlicekWeidemann32VoigtErrorFunction <: AbstractComplexErrorFunction end

"Mix of Humlicek and Weidemann (N=32) Error Function, suggested for Speed Dependent Voigt function"
struct HumlicekWeidemann32SDErrorFunction <: AbstractComplexErrorFunction end

"Humlicek with a single rational approximation."
struct CPF12ErrorFunction <: AbstractComplexErrorFunction end

"Mix of Humplicek and erfc Special Function for Voigt"
struct ErfcHumliErrorFunctionVoigt  <: AbstractComplexErrorFunction end

"Mix of Humplicek and erfc Special Function for Voigt"
struct ErfcHumliErrorFunctionSD  <: AbstractComplexErrorFunction end

"erfc Special Function for Voigt"
struct ErfcErrorFunction  <: AbstractComplexErrorFunction end

#=

Types of models that can be used to calculate an absorption cross-section:

- HitranModel
- InterpolationModel

=#

"""
    type AbstractCrossSectionModel
Abstract Cross Section Model type for generic cross section calculations
Currently either a HitranModel or an InterpolationModel
"""
abstract type AbstractCrossSectionModel end

"""
    struct HitranModel{FT}

A struct which provides all model parameters needed for cross-section
calculations using HITRAN data

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct HitranModel <: AbstractCrossSectionModel

    "Struct with hitran data"
    hitran::HitranTable
    "Broadening function (Doppler/Lorentz/Voigt)"
    broadening::AbstractBroadeningFunction
    "Wing cutoff [cm-1]"
    wing_cutoff::Real
    "VMR of gas itself [0-1]"
    vmr::Union{Real, Vector}
    "Complex Error Function to Use"
    CEF::AbstractComplexErrorFunction
    "Computer `Architecture` on which `Model` is run"
    architecture::AbstractArchitecture
end

"""
    struct InterpolationModel{FT}

A struct which provides all model parameters needed for cross-section
calculations using an Interpolator

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct InterpolationModel <: AbstractCrossSectionModel

    "The interpolator"
    itp

    # Everything below is metadata associated with the interpolator

    "The molecular species identification (ID) number"
    mol::Int
    "The isotopologue ID number"
    iso::Int
    "Wavelength grids"
    ν_grid::AbstractRange{<:Real}
    "Wavelength grids"
    p_grid::AbstractRange{<:Real}
    "Wavelength grids"
    t_grid::AbstractRange{<:Real}

end

# Types of Errors that may be thrown

struct HitranEmptyError <: Exception end
Base.showerror(io::IO, e::HitranEmptyError) = print(io, e, "No HITRAN records match the parameters")


text:
# http://www.neos-server.org/neos/NEOS-API.html
module NEOS

warn("All models submitted to NEOS become part of the public domain. For more see http://www.neos-server.org")

using LightXML
using Requests
using Codecs
using Libz

importall MathProgBase.SolverInterface

include("NEOSServer.jl")
include("NEOSSolverInterface.jl")
include("writer.jl")

include("solvers/CPLEX.jl")
include("solvers/MOSEK.jl")
include("solvers/SYMPHONY.jl")
include("solvers/XpressMP.jl")

export NEOSServer,
	NEOSCPLEXSolver, NEOSMOSEKSolver, NEOSSYMPHONYSolver, NEOSXpressMPSolver,

	addparameter!, addemail!,

 	# NEOS API functions
	neosHelp, emailHelp, welcome, version, ping, printQueue,
	listAllSolvers, listCategories,
	getSolverTemplate,
	listSolversInCategory,
	submitJob,
	getJobStatus, killJob, getFinalResults, getFinalResultsNonBlocking,
	getJobInfo,
	getIntermediateResults, getIntermediateResultsNonBlocking
end


text:
# Generated library core
"""
	WaylandInterface

The supertype of all Wayland interfaces. Subtypes must either have `id` property which holds their object id or custom [`send_request`](@ref) and [`receive_event`](@ref).

Direct analogy to Wayland's object model. The subtypes represent different interfaces and their instances represent instances of objects implementing these interfaces.
"""
abstract type WaylandInterface end
"""
	WaylandRequestMeta(target, name, opcode, args)

Request meta object. Create an instance only if you want to describe a request.

`target`: the interface the request will be targeted for
`name`: the name of the request
`opcode`: opcode of the request. Depends on the order the request is listed in, counting from 0. For example, wl_display::sync is listed first (and therefore is in the wl_display's requests array at index 1), so its opcode is 0 (index - 1).
`args`: A Vector of name=>type pairs, describing the arguments in order.

Objects of this type are created by the parser from [`SRequest`](@ref)s obtained from the xml file. Opcode is the index of the request in the [`SInterface`](@ref), starting from 0.
"""
struct WaylandRequestMeta
	target::Type{<: WaylandInterface}
	name::Symbol
	opcode::UInt16
	args::Union{Vector{Pair{Symbol, TypeofAbstractWlMsgType}}, Nothing}
end
"""
	WaylandEventMeta

Describes an event. Create an instance only if you want to describe an event - for operating on events use [`WaylandEvent`](@ref)
"""
struct WaylandEventMeta
	source::Type{<: WaylandInterface}
	name::Symbol
	opcode::UInt16
	args::Union{Vector{Pair{Symbol, TypeofAbstractWlMsgType}}, Nothing}
end
"""
	WaylandEvent{object, ename}(opcode, args)

Decoded event message. Parametrised so you can match an event "name" from wayland object "obj" with WaylandEvent{obj, :name}
"""
struct WaylandEvent{object, ename}
	source::WaylandInterface
	name::Symbol
	opcode::UInt16
	args::Vector{WlMsgType}
	WaylandEvent{object, ename}(opcode, args) where ename where object = new(object, ename, opcode, args)
end

struct WaylandEnum
	owner # WaylandInterfaceMeta, but circular dependency
	name::Symbol
	entries::AbstractDict{Symbol, WlUInt}
end

const RequestDict = Dict{Symbol, WaylandRequestMeta}
const EventDict = Dict{Symbol, WaylandEventMeta}
const EnumDict = Dict{Symbol, WaylandEnum}

struct WaylandInterfaceMeta
	requests::RequestDict
	events::EventDict
	enums::EnumDict
end
const InterfaceDict = Dict{Symbol, WaylandInterfaceMeta}

let current_newid = 0 # Because we return the value of `+= 1`, the first ID returned is 1.
	"""
		get_newid()

	Gets the next free ID in client range, which is [1, 0xff000000)
	"""
	global get_newid() = current_newid < 0xff000000 ? current_newid += 1 : current_newid = 1 # TODO: Check if it's not still used, even after a full loop. (Not necessary for testing, as a testing client will never go through 16 777 215 IDs)
end
#="""
	used_ids

A Set of IDs that are in use, updated every time an object gets created by us, and when it gets destroyed. NYI
"""
used_ids = Set{WlID}=#
"""
	send_request(io::IO, from::WlID, opcode, args::WlMsgType...)

Create the request message and send it to the IO.

This is part of the low level interface used by the generated library. You probably want the generated functions instead.
"""
function send_request(io::IO, from::WlID, opcode, args::WlMsgType...)
	size = 8 # just the header
	if isempty(args)
		payload = nothing
	else
		payload = IOBuffer()
		for arg in args
			size += write(payload, arg)
		end
	end
	send(io, GenericMessage(from, size, opcode, payload))
end
"""
	send_request(to::WaylandInterface, from::WlID, opcode, args::WlMsgType...)

[`send_request`](@ref) to the given object's ID.
"""
send_request(to::WaylandInterface, from::WlID, opcode, args::WlMsgType...) = send_request(to.id, from, opcode, args...)
"""
    receive_event(io::IO; msg_type=GenericMessage)

Receive an event message and decode it into a `WaylandEvent` object.

This is part of the low level interface used by the generated library. You probably want the generated high-level event listening utilities instead.
"""
function receive_event(io::IO, msg_type=GenericMessage)
	receive(io, msg_type)
end

# Library generation
"""
	genrequest(meta)

Generate a request function from meta.

meta must contain fields:
.target: a `WaylandInterface` subtype - type of the object the request will be sent to
.name: a `Symbol` which will become the function name
.opcode: a `UInt16` - the request opcode
.args: an ordered iterable collection of `name=>Type` `Pair{Symbol,TypeofWlMsgType}`s describing the (`name::Type`) arguments.

new_id arguments aren't presented to the user, but are acquired by the get_newid call in send_request arguments.
"""
function genrequest(meta)
	if meta.args == nothing
		:( $(meta.name)(target::$(meta.target)) = send_request(target, $(meta.opcode)) )
	else
		head_args = Vector{Expr}()
		tail_args = Vector{Expr}()
		for (name, type) in meta.args
			if type == WlNewID
				push!(tail_args, :(get_newid()))
			else
				push!(head_args, :($name::$type))
				push!(tail_args, name)
			end
		end
		:( $(meta.name)(target::$(meta.target), $(head_args...)) = send_request(target, $(meta.opcode), $(tail_args...)) )
	end
end
"""
	ListenerRegistry

As the name says. Dictionary mapping event types to `Set`s of their listeners.
"""
ListenerRegistry = Dict{WaylandEvent, Set{Function}}
"""
	GlobalListenerRegistry

Global, default instance of [`ListenerRegistry`](@ref).
"""
global_listener_registry = ListenerRegistry()
"""
	addlistener(f, event::Type{<: WaylandEvent}, data = nothing)

Register `f` to be called every time an event of `event` type gets fired.

The default implementations only pass one argument to `f` – the event object itself – if `data` is set to nothing, and two arguments – the event object as first and `data` as second – otherwise.
"""
addlistener(f, event::Type{WaylandEvent}) = error("Generic `addlistener` called for $(event)!")
"""
	remlistener(f, event::Type{WaylandEvent}, throws = true)

Remove `f` from registered listeners on event type `event`, returning `f`. If `throws` is true, throws ArgumentError if `f` is not registered as listener on `event`, otherwise returns false in that case.
"""
remlistener(f, event::Type{WaylandEvent}, throws = true) = error("Generic `remlistener` called for $(event)!")
"""
    genlibclient(protocols::Set{SProtocol})

Generate the client library. Entry point of the generation step.

The generated library consists of a internally used interface dictionary of type InterfaceDict holding metadata, and structs and methods meant for direct use.

The interface dictionary is a Dict which maps a Symbol interface name to a simple WaylandInterfaceMeta object, consisting only of three dictionaries: `requests::RequestDict`, `events::EventDict` and `enums::EnumDict`. It's the easiest way to check enum values and events to register listeners to.

The generated structs are subtypes of [`WaylandInterface`](@ref) and are described there, while methods are requests – which have the names of the requests they represent – and `addlistener` + `remlistener` duo for listening to events. For example, parsing `wl_display` interface will generate a `WlDisplay` struct for representing such object, a [`get_registry(target::WlDisplay)`](@ref) method for sending its `get_registry` request and [`addlistener(f, event::Type{WaylandEvent{ID, :error}}) where ID`](@ref), [`remlistener(f, event::Type{WaylandEvent{ID, :error}}) where ID`](@ref) duo for interacting with its "error" events.

By default there is only one InterfaceDict holding all known interfaces, returned from one call to genlibclient. However, if the separation of loaded protocols is desired, it is entirely possible to call this function multiple times, resulting in multiple working dictionaries and multiple but indistinguishable sets of generated methods. (Do note that in the event of repeating interface name there will be a namespace clash of generated globally-visible structs and, if they have the same named requests/events, the methods)

To be implemented:

- Treating destructor requests differently so that they get processed automatically without the need for the user to manually call them.
"""
function genlibclient(protocols::Set{SProtocol})
	ifs = Dict{Symbol, WaylandInterfaceMeta}()
	for prot in protocols
		for interface in prot.interfaces
			# Generate the interface type.
			@eval begin
				interface = $interface
				iftypename = tostructid(interface.name)
				"""
					$iftypename

				Represents a $(interface.name) object.
				"""
				struct $iftypename <: WaylandInterface # Visible in global scope since that's where [`eval`](@ref) executes.
					id::WlID
				end
				iftype = $iftypename
			end
			reqs = RequestDict()
			for (index, request) in enumerate(interface.requests)
				if request.args == nothing
					args = nothing
				else
					args = Vector{Pair{Symbol, TypeofAbstractWlMsgType}}()
					for arg in request.args
						push!(args, Symbol(arg.name) => arg.type)
					end
				end
				request_meta = WaylandRequestMeta(iftype, Symbol(request.name), index - 1, args)
				push!(reqs, request_meta.name => request_meta)
				# [`eval`](@ref) executes in module's global scope, so the request methods along with their docs will too be visible.
				docargstring = "target::$(request_meta.target)"
				if args != nothing
					for (name, type) in args
						docargstring *= ", $name::$type"
					end
				end
				"""
					$(request_meta.name)($docargstring)

				$(repr(request.description))
				"""
				eval(genrequest(request_meta))
			end
			evs = EventDict()
			for (index, event) in enumerate(interface.events)
				if event.args == nothing
					args = nothing
				else
					args = Vector{Pair{Symbol, TypeofAbstractWlMsgType}}()
					for arg in event.args
						push!(args, Symbol(arg.name) => arg.type)
					end
				end
				event_meta = WaylandEventMeta(iftype, Symbol(event.name), index - 1, args)
				push!(evs, event_meta.name => event_meta)
				# Add the event to global listener registry
				event_type = WaylandEvent{iftype, event_meta.name}
				push!(global_listener_registry, event_type => Set{Function}())
				# not sure if eval is needed here
				@eval begin
					function addlistener(f::Function, ::Type{WaylandEvent{$iftype, $(event_meta.name)}})
						push(global_listener_registry[$event_type], f)
					end
					function remlistener(f::Function, ::Type{WaylandEvent{$iftype, $(event_meta.name)}}, throws = true)
						if f in global_listener_registry[$event_type]
							delete!(global_listener_registry, f)
							return f
						else
							if throws
								throw(ArgumentError("$f wasn't registered in global_listener_registry."))
							else
								return false
							end
						end
					end
				end
			end
			ens = EnumDict()
			ifmeta = WaylandInterfaceMeta(reqs, evs, ens) # needed here and not in the push!() below because it's needed for WaylandEnum construction
			for enum in interface.enums
				entries = Dict{Symbol, WlUInt}()
				for (name, entry) in enum.entries
					push!(entries, name => entry.value)
				end
				push!(ens, Symbol(enum.name) => WaylandEnum(ifmeta, Symbol(enum.name), entries))
			end
			push!(ifs, Symbol(interface.name) => ifmeta)
		end
	end
	return ifs
end
"""
	genlibclient(path::AbstractString)

Generate the client library from a protocol XML file denoted by `path`.

Equivalent to calling `genlibclient(wlparse(path))`.
"""
genlibclient(path::AbstractString) = genlibclient(wlparse(path))
"""
	genlibclient(paths::Set{AbstractString})

Generate the client library from a Set of protocol XML files.

Equivalent to calling `genlibclient(reduce(append!, wlparse.(paths)))`.
"""
genlibclient(paths::Set{AbstractString}) = genlibclient(reduce(append!, wlparse.(paths)))


text:
"""
Type: NLPAtX

Methods: update!, reinit!

NLPAtX contains the information concerning a nonlinear problem at
the iteration x.

min_{x ∈ ℜⁿ} f(x) subject to lcon <= c(x) <= ucon, lvar <= x <= uvar.

Basic information is:
 - x : the current candidate for solution to our original problem
 - fx : which is the funciton evaluation at x
 - gx : which is the gradient evaluation at x
 - Hx : which is the hessian representation at x

 - mu : Lagrange multiplier of the bounds constraints

 - cx : evaluation of the constraint function at x
 - Jx : jacobian matrix of the constraint function at x
 - lambda : Lagrange multiplier of the constraints

 - current_time  : time
 - current_score : score
 - evals : number of evaluations of the function (import the type NLPModels.Counters)

Note: * by default, unknown entries are set to nothing (except evals).
      * All these information (except for x and lambda) are optionnal and need to be update when
        required. The update is done trhough the update! function.
      * x and lambda are mandatory entries. If no constraints lambda = [].
      * The constructor check the size of the entries.
"""
mutable struct 	NLPAtX <: AbstractState

#Unconstrained State
    x            :: AbstractVector     # current point
    fx           :: FloatVoid   # objective function
    gx           :: Iterate     # gradient size: x
    Hx           :: MatrixType  # hessian size: |x| x |x|

#Bounds State
    mu           :: Iterate     # Lagrange multipliers with bounds size of |x|

#Constrained State
    cx           :: Iterate     # vector of constraints lc <= c(x) <= uc
    Jx           :: MatrixType  # jacobian matrix, size: |lambda| x |x|
    lambda       :: AbstractVector    # Lagrange multipliers

 #Resources State
    current_time   :: FloatVoid
    current_score  :: FloatVoid
    evals          :: Counters

 function NLPAtX(x             :: AbstractVector,
                 lambda        :: AbstractVector;
                 fx            :: FloatVoid    = nothing,
                 gx            :: Iterate      = nothing,
                 Hx            :: MatrixType   = nothing,
                 mu            :: Iterate      = nothing,
                 cx            :: Iterate      = nothing,
                 Jx            :: MatrixType   = nothing,
                 current_time  :: FloatVoid    = nothing,
                 current_score :: FloatVoid    = nothing,
                 evals         :: Counters     = Counters())

  _size_check(x, lambda, fx, gx, Hx, mu, cx, Jx)

  return new(x, fx, gx, Hx, mu, cx, Jx, lambda, current_time, current_score, evals)
 end
end

"""
An additional constructor for unconstrained problems
"""
function NLPAtX(x             :: AbstractVector;
                fx            :: FloatVoid    = nothing,
                gx            :: Iterate      = nothing,
                Hx            :: MatrixType   = nothing,
                mu            :: Iterate      = nothing,
                current_time  :: FloatVoid    = nothing,
                current_score :: FloatVoid    = nothing,
                evals         :: Counters     = Counters())

    _size_check(x, zeros(0), fx, gx, Hx, mu, nothing, nothing)

	return NLPAtX(x, zeros(0), fx = fx, gx = gx,
                  Hx = Hx, mu = mu, current_time = current_time,
                  current_score = current_score, evals = evals)
end

"""
reinit!: function that set all the entries at void except the mandatory x

Warning: if x, lambda or evals are given as a keyword argument they will be
prioritized over the existing x, lambda and the default Counters.
"""
function reinit!(stateatx :: NLPAtX, x :: AbstractVector, l :: AbstractVector; kwargs...)

 for k ∈ fieldnames(typeof(stateatx))
   if !(k ∈ [:x,:lambda,:evals]) setfield!(stateatx, k, nothing) end
 end

 return update!(stateatx; x=x, lambda = l, evals = Counters(), kwargs...)
end
"""
reinit!: short version of reinit! reusing the x in the state

Warning: if x, lambda or evals are given as a keyword argument they will be
prioritized over the existing x, lambda and the default Counters.
"""
function reinit!(stateatx :: NLPAtX; kwargs...)
 return reinit!(stateatx, stateatx.x, stateatx.lambda; kwargs...)
end

"""
_size_check!: check the size of the entries in the State
"""
function _size_check(x, lambda, fx, gx, Hx, mu, cx, Jx)

    if gx != nothing && length(gx) != length(x)
     throw(error("Wrong size of gx in the NLPAtX."))
    end
    if Hx != nothing && size(Hx) != (length(x), length(x))
     throw(error("Wrong size of Hx in the NLPAtX."))
    end
    if mu != nothing && length(mu) != length(x)
     throw(error("Wrong size of mu in the NLPAtX."))
    end

    if lambda != zeros(0)
        if cx != nothing && length(cx) != length(lambda)
         throw(error("Wrong size of cx in the NLPAtX."))
        end
        if Jx != nothing && size(Jx) != (length(lambda), length(x))
         throw(error("Wrong size of Jx in the NLPAtX."))
        end
    end

end


text:
#Extract variable names from a formula before creating schema
function varnames(f::FormulaTerm)
  rhs = varnames(f.rhs)
  if rhs isa Symbol rhs = [rhs] end
  return (lhs=varnames(f.lhs), rhs=rhs)
end

function varnames(f::Term)
  return Symbol(f)
end

function varnames(f::StatsModels.TupleTerm)
  raw = map(x -> varnames(x), f)
  return get_unique(raw)
end

function varnames(f::ConstantTerm)
  return nothing
end

function varnames(f::InteractionTerm)
  return map(x -> varnames(x), f.terms)
end

function get_unique(f::Tuple)
  out = Symbol[]
  fl(v) = for x in v
    if x isa Tuple fl(x) else push!(out, x) end
  end
  fl(f)
  return unique(out)
end

function combine_formulas(formulas::Union{Array,FormulaTerm}...)
  formulas = map(x -> x isa FormulaTerm ? [x] : x, formulas)
  form_combs = collect(Base.product(formulas...))
  return reshape(form_combs, length(form_combs))
end

"""

    allsub(f::FormulaTerm)

Takes a formula `f` and returns an array containing formulas for all
subsets of the formula terms.
"""
function allsub(f::FormulaTerm)
  if f.rhs isa Union{ConstantTerm,Term,Tuple{Term}} return [f] end
  combs = Tuple.(collect(combinations(Symbol.(f.rhs))))
  terms = map(x -> term.(x), combs)
  forms = map(x -> term(f.lhs.sym) ~ sum(x), terms)
  return [term(f.lhs.sym) ~ term(1); forms]
end

max_stars_repo_name: kenkellner/julia-hme

max_issues_repo_name: kenkellner/julia-hme

max_forks_repo_name: kenkellner/julia-hme


text:
module TestIncMat
include("../src/lina.jl")
include("square.jl")

using Test
import SparseArrays: nnz
import .LinA: IncMat

@testset "Create D: Square               " begin
    edges, n = square_edges()
    m = length(edges)
    lambda = fill(0.5, length(edges))
    D = IncMat(edges, lambda, n)
    D_expected =
        [0.5  -0.5   0.0   0.0
         0.5   0.0  -0.5   0.0
         0.0   0.5   0.0  -0.5
         0.0   0.0   0.5  -0.5]

    @test size(D) == (m, n)
    @test nnz(D) == 2m
    @test D_expected ≈ Array(D)
end

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
# LibTOAST.jl: interface to the TOAST++ library
# Copyright (C) 2019 Samuel Powell

# Construct a one-based Julia CSC matrix from the row pointer, column index
# and values of an m x n zero-based CSR matrix.
function _CSR0_to_CSC1(m,n,rowptr,colidx,values)

  _CSR0_to_CSC1(m,n,rowptr,colidx,pointer(values))

end

# Construct a one-based Julia CSC matrix from the row pointer, column index
# and values of an m x n zero-based CSR matrix.
function _CSR0_to_CSC1(m,n,rowptr,colidx,values::Ptr)

  # From toastmatlab.cc
  # void CopyMatrix (mxArray **array, const RCompRowMatrix &mat)

  # Array of pointers to row index, column pointer, real pointer
  cscval = Array{Float64}(undef, length(colidx))
  rowidx = Array{Cint}(undef, length(colidx))
  colptr = Array{Cint}(undef, n+1)

  icxx"""
    int i, j, k;
    int m = $m;
    int n = $n;
    int *rowptr = $(pointer(rowptr));
    int *colidx = $(pointer(colidx));

    const double *pval = $values;
    int nz = rowptr[m];

    // This should maybe the larger than int!
    int *rcount = new int[n];
    int idx;

    //mxArray *tmp = mxCreateSparse (m, n, nz, mxREAL);
    for (i = 0; i < n; i++) rcount[i] = 0;
    for (i = 0; i < nz; i++) rcount[colidx[i]]++;

    double *pr = $(pointer(cscval));
    int *ir = $(pointer(rowidx));
    int *jc = $(pointer(colptr));

    jc[0] = 0;
    for (i = 0; i < n; i++) jc[i+1] = jc[i]+rcount[i];
    for (i = 0; i < n; i++) rcount[i] = 0;
    for (i = 0; i < m; i++)
    for (k = rowptr[i]; k < rowptr[i+1]; k++) {
      j = colidx[k];
      idx = jc[j]+rcount[j];
      ir[idx] = i;
      pr[idx] = pval[k];
      rcount[j]++;
    }

    delete []rcount;
  """

  # Build sparse with 1-based indexing and Int size pointers
  return SparseMatrixCSC(n,m,Vector{Int}(colptr.+1),Vector{Int}(rowidx.+1),cscval)

end


text:
require(Pkg.dir("MethylUtils", "test","mlml-test-2.jl"))
include(Pkg.dir("MethylUtils", "src","Features.jl"))

#d = Dzeros()
#mlml_parallel_performance()

## TIMING
TIMING_PARALLEL = false
TIMING = false
PROFILING = false

PARALLEL_CYTOSINE_COUNTS = true

println("Get random counts")
@everywhere meth_counts = rand_bs_ox_counts(10,100,5)
println("Finished generating random counts")
#println(meth_counts)

if TIMING_PARALLEL
  d = dzeros(10)
  println("Start timing (parallel) EM")
  cgs = @time mlml_parallel_performance(d,meth_counts)
  type_cgs = typeof(cgs)
  println("Got $type_cgs")
  println(cgs)
end

if TIMING
  println("Start timing EM")
  cgs = @time mlml_performance(meth_counts)
  println("EM Done in time: $elapsed")
end

if PROFILING
    println("Start profiling EM")
    @profile mlml_performance(meth_counts)
    println ("Done profiling")
end

println("Get random CytosineCounts")
@everywhere bs_counts = rand_CytosineCounts(10,100,5)
@everywhere ox_counts = rand_CytosineCounts(10,100,5)

println("Finished generating random counts")
println(bs_counts)
println(ox_counts)

if PARALLEL_CYTOSINE_COUNTS
  d = dzeros(10)
  println("Start timing (parallel) EM")
  cgs = @time mlml_parallel_performance(d,bs_counts, ox_counts)
  type_cgs = typeof(cgs)
  println("Got $type_cgs")
  println(cgs)
end


text:
module Knockoffs

using LinearAlgebra
using Statistics
using Convex
using SCS
using SnpArrays
using DelimitedFiles
using ProgressMeter
using Distributions
using CSV
using DataFrames
using fastPHASE

import Base: eltype, getindex, size

export fixed_knockoffs, modelX_gaussian_knockoffs, normalize_col!,
    coefficient_diff, threshold, extract_beta, extract_combine_beta,
    partition, rapid, snpknock2, decorrelate_knockoffs,
    process_fastphase_output, fastphase,
    # functions for hmm
    get_haplotype_transition_matrix, get_genotype_transition_matrix,
    get_initial_probabilities,
    GenotypeState, MarkovChainTable, pair_to_index, index_to_pair,
    hmm_knockoff, forward_backward_sampling, forward_backward_sampling!,
    form_emission_prob_matrix, get_genotype_emission_probabilities,
    get_haplotype_emission_probabilities, markov_knockoffs, markov_knockoffs!,
    sample_markov_chain, sample_markov_chain!,
    update_normalizing_constants!, single_state_dmc_knockoff!

include("struct.jl")
include("fixed.jl")
include("modelX.jl")
include("dmc.jl")
include("feature_stats.jl")
include("threshold.jl")
include("hmm_wrapper.jl")
include("hmm.jl")

end