text:
#####
##### `AbstractLPCMSerializer`
#####

"""
    AbstractLPCMSerializer

A type whose subtypes support:

  - [`deserialize_lpcm`](@ref)
  - [`serialize_lpcm`](@ref)

All definitions of subtypes of the form `S<:AbstractLPCMSerializer` must also support
a constructor of the form `S(::Signal)` and overload `Onda.serializer_constructor_for_file_extension`
with the appropriate file extension.

See also: [`serializer`](@ref), [`LPCM`](@ref), [`LPCMZst`](@ref)
"""
abstract type AbstractLPCMSerializer end

"""
    Onda.serializer_constructor_for_file_extension(::Val{:extension_symbol})

Return a constructor of the form `S(::Signal)::AbstractLPCMSerializer`
corresponding to the provided extension.

This function should be overloaded for new `AbstractLPCMSerializer` subtypes.
"""
function serializer_constructor_for_file_extension(::Val{unknown}) where {unknown}
    return throw(ArgumentError("unknown file extension: $unknown"))
end

function register_file_extension_for_serializer(extension::Symbol,
                                                T::Type{<:AbstractLPCMSerializer})
    return error("""
                 `Onda.register_file_extension_for_serializer(ext, T)` is deprecated; instead, `AbstractLPCMSerializer`
                 authors should define `Onda.serializer_constructor_for_file_extension(::Val{ext}) = T`.
                 """)
end

"""
    serializer(signal::Signal; kwargs...)

Return `S(signal; kwargs...)` where `S` is the `AbstractLPCMSerializer` that
corresponds to `signal.file_extension` (as determined by the serializer author
via `serializer_constructor_for_file_extension`).

See also: [`deserialize_lpcm`](@ref), [`serialize_lpcm`](@ref)
"""
function serializer(signal::Signal; kwargs...)
    T = serializer_constructor_for_file_extension(Val(signal.file_extension))
    return T(signal; kwargs...)
end

"""
    deserialize_lpcm(bytes, serializer::AbstractLPCMSerializer)

Return a channels-by-timesteps `AbstractMatrix` of interleaved LPCM-encoded
sample data by deserializing the provided `bytes` from the given `serializer`.

Note that this operation may be performed in a zero-copy manner such that the
returned sample matrix directly aliases `bytes`.

This function is the inverse of the corresponding [`serialize_lpcm`](@ref)
method, i.e.:

```
serialize_lpcm(deserialize_lpcm(bytes, serializer), serializer) == bytes
```

    deserialize_lpcm(bytes, serializer::AbstractLPCMSerializer, sample_offset, sample_count)

Similar to `deserialize_lpcm(bytes, serializer)`, but deserialize only the segment requested
via `sample_offset` and `sample_count`.

    deserialize_lpcm(io::IO, serializer::AbstractLPCMSerializer[, sample_offset, sample_count])

Similar to the corresponding `deserialize_lpcm(bytes, ...)` methods, but the bytes
to be deserialized are read directly from `io`.

If `sample_offset`/`sample_count` is provided and `io`/`serializer` support
seeking, implementations of this method may read only the bytes required to
extract the requested segment instead of reading the entire stream.
"""
function deserialize_lpcm end

"""
    serialize_lpcm(samples::AbstractMatrix, serializer::AbstractLPCMSerializer)

Return the `AbstractVector{UInt8}` of bytes that results from serializing `samples`
to the given `serializer`, where `samples` is a channels-by-timesteps matrix of
interleaved LPCM-encoded sample data.

Note that this operation may be performed in a zero-copy manner such that the
returned `AbstractVector{UInt8}` directly aliases `samples`.

This function is the inverse of the corresponding [`deserialize_lpcm`](@ref)
method, i.e.:

```
deserialize_lpcm(serialize_lpcm(samples, serializer), serializer) == samples
```

    serialize_lpcm(io::IO, samples::AbstractMatrix, serializer::AbstractLPCMSerializer)

Similar to the corresponding `serialize_lpcm(samples, serializer)` method, but serializes
directly to `io`.
"""
function serialize_lpcm end

#####
##### fallbacks/utilities
#####

function deserialize_lpcm(bytes, serializer::AbstractLPCMSerializer, args...)
    return deserialize_lpcm(IOBuffer(bytes), serializer, args...)
end

function serialize_lpcm(samples::AbstractMatrix, serializer::AbstractLPCMSerializer)
    io = IOBuffer()
    serialize_lpcm(io, samples, serializer)
    return resize!(io.data, io.size)
end

jump(io::IO, n) = (read(io, n); nothing)
jump(io::IOStream, n) = (skip(io, n); nothing)
jump(io::IOBuffer, n) = ((io.seekable ? skip(io, n) : read(io, n)); nothing)

unsafe_vec_uint8(x::AbstractVector{UInt8}) = convert(Vector{UInt8}, x)
function unsafe_vec_uint8(x::Base.ReinterpretArray{UInt8,1})
    return unsafe_wrap(Vector{UInt8}, pointer(x), length(x))
end

#####
##### `LPCM`
#####

const LPCM_SAMPLE_TYPE_UNION = Union{Int8,Int16,Int32,Int64,UInt8,UInt16,UInt32,UInt64}

"""
    LPCM{S}(channel_count)
    LPCM(signal::Signal)

Return a `LPCM<:AbstractLPCMSerializer` instance corresponding to Onda's default
interleaved LPCM format assumed for signal files with the ".lpcm" extension.

`S` corresponds to `signal.sample_type`, while `channel_count` corresponds to
`signal.channel_names`.

Note that bytes (de)serialized via this serializer are little-endian per the
Onda specification.
"""
struct LPCM{S<:LPCM_SAMPLE_TYPE_UNION} <: AbstractLPCMSerializer
    channel_count::Int
end

LPCM(signal::Signal) = LPCM{signal.sample_type}(length(signal.channel_names))

serializer_constructor_for_file_extension(::Val{:lpcm}) = LPCM

function deserialize_lpcm(bytes, serializer::LPCM{S}) where {S}
    sample_count = Int(length(bytes) / sizeof(S) / serializer.channel_count)
    return deserialize_lpcm(bytes, serializer, 0, sample_count)
end

function deserialize_lpcm(bytes, serializer::LPCM{S}, sample_offset, sample_count) where {S}
    i = (serializer.channel_count * sample_offset) + 1
    j = serializer.channel_count * (sample_offset + sample_count)
    return reshape(view(reinterpret(S, bytes), i:j),
                   (serializer.channel_count, sample_count))
end

function deserialize_lpcm(io::IO, serializer::LPCM)
    return deserialize_lpcm(read(io), serializer)
end

function deserialize_lpcm(io::IO, serializer::LPCM{S}, sample_offset,
                          sample_count) where {S}
    bytes_per_sample = sizeof(S) * serializer.channel_count
    jump(io, bytes_per_sample * sample_offset)
    return deserialize_lpcm(read(io, bytes_per_sample * sample_count), serializer)
end

function _validate_lpcm_samples(samples::AbstractMatrix{S}, serializer::LPCM{S}) where {S}
    serializer.channel_count == size(samples, 1) && return nothing
    return throw(ArgumentError("`samples` row count does not match expected channel count"))
end

function serialize_lpcm(io::IO, samples::AbstractMatrix, serializer::LPCM)
    _validate_lpcm_samples(samples, serializer)
    return write(io, samples)
end

function serialize_lpcm(samples::Matrix, serializer::LPCM)
    _validate_lpcm_samples(samples, serializer)
    return reinterpret(UInt8, vec(samples))
end

#####
##### `LPCMZst`
#####

"""
    LPCMZst(lpcm::LPCM; level=3)
    LPCMZst(signal::Signal; level=3)

Return a `LPCMZst<:AbstractLPCMSerializer` instance that corresponds to
Onda's default interleaved LPCM format compressed by `zstd`. This serializer
is assumed for signal files with the ".lpcm.zst" extension.

The `level` keyword argument sets the same compression level parameter as the
corresponding flag documented by the `zstd` command line utility.

See https://facebook.github.io/zstd/ for details about `zstd`.
"""
struct LPCMZst{S} <: AbstractLPCMSerializer
    lpcm::LPCM{S}
    level::Int
    LPCMZst(lpcm::LPCM{S}; level=3) where {S} = new{S}(lpcm, level)
end

LPCMZst(signal::Signal; kwargs...) = LPCMZst(LPCM(signal); kwargs...)

serializer_constructor_for_file_extension(::Val{Symbol("lpcm.zst")}) = LPCMZst

function deserialize_lpcm(bytes, serializer::LPCMZst, args...)
    bytes = unsafe_vec_uint8(bytes)
    return deserialize_lpcm(zstd_decompress(bytes), serializer.lpcm, args...)
end

function deserialize_lpcm(io::IO, serializer::LPCMZst, args...)
    reader = io -> deserialize_lpcm(io, serializer.lpcm, args...)
    return zstd_decompress(reader, io)
end

function serialize_lpcm(samples::AbstractMatrix, serializer::LPCMZst)
    bytes = serialize_lpcm(samples, serializer.lpcm)
    return zstd_compress(unsafe_vec_uint8(bytes), serializer.level)
end

function serialize_lpcm(io::IO, samples::AbstractMatrix, serializer::LPCMZst)
    writer = io -> serialize_lpcm(io, samples, serializer.lpcm)
    return zstd_compress(writer, io, serializer.level)
end


text:
using PyCall

if Sys.iswindows()
	@warn("""
					FEniCS is not known to run on Windows!
					Please try installing from Windows Subsystem for Linux.
		""")
end
	try
		pyimport_conda("fenics", "fenics=2019.1.0", "conda-forge")
		pyimport_conda("mshr", "mshr=2019.1.0", "conda-forge")
	catch ee
		typeof(ee) <: PyCall.PyError || rethrow(ee)
		@warn("""
						Python Dependancies not installed
						Please either:
						 - Rebuild PyCall using the path to FEniCS using
						 - `ENV["PYTHON"]="/path/to/FEniCS"; Pkg.build("PyCall"); Pkg.build("FEniCS")`
						 - Or install the Dependancies using the Docker image provided
				 """
		)
	end

ext: jl


text:
# Normalization factor

@test normalize_factor([1,2]) == [1/3, 2/3]
@test normalize_factor([1 2]) == [1/3 2/3]

# Gaussian Distribution PDF's

@test typeof(normpdf(0,1)) == Function
@test_approx_eq_eps normpdf(0,1)(0) 0.399 0.0001

# Group combinations

groups=Array[["a","b"],["c","d"],["e","f","g"],["z","q"]]
grc=group_combinations(groups)

@test Set({"a","c","e","z"}) in grc
@test length(unique(grc)) == prod(map(length,groups))
@test group_combinations(Array[]) == []


text:
module ExtensionTests

"!!set(name:test)!!get(name)"
set_and_get = ()

"!!var(name:test)"
var = ()

"!!summary(summary)"
summary = ()

"!!longform(...)"
longform = ()

"!!include(includes/file.md)"
includes = ()

"!!set(name:test) !!include(includes/file.md) !!get(name)"
set_includes_get = ()

"!!include(includes/filen_metas.md) Get the inner again: !!get(笔者_inner)"
includes_nested = ()

end


text:
### Domain [mm]
x_max = 1

### Models ###
# The considered PSF for the considered imaging tecnique
	# sigma also defines the resolution of the imaging system.
	sigma = 0.04
	filter = (x, y) -> exp(-(x^2 + y^2)/2/sigma^2)
	filter_dx = (x, y) -> -x/sigma^2*filter(x, y)
	filter_dy = (x, y) -> -y/sigma^2*filter(x, y)
# The number of time samples for reconstruction: 2K+1
	K = 2
# the maximum allowed speed for any considered particle
	v_max = 15 #[mm/s]
        particle_velocity = 1 #[mm/s]
# The sampling rate
	tau = 1/500

### Noise ###
sigma_noise = 0.01

### Medium ###
# total number of the experimets measurements (gives the length of the experience)
	n_im = 500
# Vessels separation
	dx = 0.02
# microbubbles distance to booundary
	bndry_sep = 0.05

### Inverse problem ###
# Considered norm to minimize
	lp_norm = 1
# Threshold on when to consider that there was a difference on the number of particles
# at two consecutive time steps.
	jump_threshold = 0.1

### Sequence ###
# Probability of a particle to activate or deactivate
	activation_probability=0.02

### Function describing a curved domain describing a vessel. At the beggining we impose
### a string line, that afterwards curves

# Elliptic extension of a vessel
using Interpolations
using Roots
using QuadGK

function vessel(t)
    # Straight line parameters and starting point
    line_len = 1/4
    start_pos = 1/4
    # Ellipse curve parameters
    x_semiax = x_max - start_pos*x_max
    y_semiax = 8*(x_max - line_len*x_max)
    max_t = 2
    # when straight line
    if t < line_len*x_max
        return (start_pos*x_max, t)
    end
    # else, the ellipse
    tt = t - line_len*x_max
    ellipse_curve(tt) = (-x_semiax*cos(tt) + x_semiax, y_semiax*sin(tt))
    # We seek to obtain the parametrization with speed one.
    curve_deriv(tt) = (x_semiax*sin(tt), y_semiax*cos(tt))
    deriv_norm(tt) = sqrt(curve_deriv(tt)[1]^2 + curve_deriv(tt)[2]^2)
    arclength(tt) = quadgk(deriv_norm, 0, tt)
    # we invert the arclength
    t_sample = [t for t = linspace(0,max_t,1000)]
    arclength_sample = [arclength(t)[1] for t in t_sample]
    inv_interp = LinearInterpolation(arclength_sample, t_sample)
    natural_curve(tt) = ellipse_curve(inv_interp(tt))
    x0 = start_pos*x_max
    y0 = line_len*x_max
    return (natural_curve(tt)[1] + x0, natural_curve(tt)[2] + y0)
end

target_funct(t) = max(vessel(t)[1] - x_max*(1-bndry_sep), vessel(t)[2] - x_max*(1-bndry_sep))
impossible_time = 1.5
lower_time = 0
if target_funct(lower_time)*target_funct(impossible_time) > 0
    print("This is not going to work")
end
t_max = find_zero(target_funct, (lower_time, impossible_time), Bisection())


text:
function main!(r)
    counters = move_West_Nord(r)
    move_forward!(r, Sud)
    putmarker!(r)
    row_length = move_forward!(r, Ost)
    move!(r, Nord)
    snake!(r, row_length)
    move_West_Nord(r)
    go_nachalo!(r, counters)
end

function snake!(r, counter)
    direction = 1
    counter+=1
    counter1 = counter-1
    while counter1!=0
        i = 0
        if direction == 1
            while i != counter - counter1
                if !isborder(r, HorizonSide(1))
                    move!(r, HorizonSide(1))
                    i+=1
                else
                    i += obhod!(r, direction)
                    i+=1
                end
            end
            putmarker!(r)
        else
            putmarker!(r)
        end
        while i != counter-1
            if !isborder(r, HorizonSide(direction))
                move!(r, HorizonSide(direction))
                i += 1
                if direction==1
                    putmarker!(r)
                end
            else
                i += obhod!(r, direction)
                i+=1
                if i<=counter1+1
                    putmarker!(r)
                end
            end
            if i<counter1
                putmarker!(r)
            end
        end
        if !isborder(r,Nord)
            move!(r, Nord)
        else
            return 0
        end
        direction = (direction+2)%4
        counter1 -= 1
    end
end

function move_forward!(r, side)
    c = 0
    while !isborder(r, side)
        move!(r, side)
        c+=1
        if side==Ost
            putmarker!(r)
        end
    end
    return c
end

function obhod!(r, side)
    counter = 0
    counter1 = 0
    while isborder(r, HorizonSide(side))
        move!(r, Nord)
        counter+=1
    end
    move!(r, HorizonSide(side))
    while isborder(r, Sud)
        move!(r, HorizonSide(side))
        counter1 += 1
    end
    for i in 0:counter-1
        move!(r, Sud)
    end
    return(counter1)
end

function move_West_Nord(r)           # идем в левый верхний угол
    counters = []
    while !(isborder(r,Nord) && isborder(r,West))
        counter_N = 0
        while !isborder(r,Nord)
            move!(r, Nord)
            counter_N+=1
        end
        pushfirst!(counters, counter_N)
        counter_W = 0
        while !isborder(r,West)
            move!(r, West)
            counter_W+=1
        end
        pushfirst!(counters, counter_W)
    end
    return counters
end

function go_nachalo!(r, counters)   # возвращаемся в начало
    for i in 1:size(counters, 1)
        if i%2==0
            for j in 0:counters[i]-1
                move!(r, Sud)
            end

        end
        if i%2==1
            for j in 0:counters[i]-1
                move!(r, Ost)
            end
        end
    end
end


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

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

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

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

@parallel function viscosity2!(η, x, y)
    @all(η) = _viscosity(@all(x), @all(y))
    return
end

viscosity(x, y) = [ _viscosity(xi, yi) for xi in x, yi in y]

function testme()

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

    return η1 == η2
end

@test  testme() == true

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

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

R.=0; @parallel idx_x!(R, x); @test reduce(+, Rcol == x for Rcol in eachcol(R)) == size(R,2)
R.=0; @parallel idx_y!(R, y); @test reduce(+, Rrow == y for Rrow in eachrow(R)) == size(R,1)
R.=0; @parallel idx_x!(R, xvec); @test reduce(+, Rcol == xvec for Rcol in eachcol(R)) == size(R,2)
R.=0; @parallel idx_y!(R, yvec); @test reduce(+, Rrow == yvec for Rrow in eachrow(R)) == size(R,1)


text:
"""
    sparsity(A)

Calculates the sparsity of the given array.
The sparsity is defined as number of zero-valued elements divided by
total number of elements.
"""
function sparsity(A::AbstractArray{T}) where T<:Number
    (length(A)-countnz(A))/length(A)
end

"""
    reldiff(A, B)

Relative difference of absolute values of `A` and `B` defined as

``
\\operatorname{reldiff} = 2 \\dfrac{\\operatorname{abs}(A - B)}{\\operatorname{abs}(A+B)}.
``
"""
function reldiff(A::AbstractArray{T}, B::AbstractArray{S}) where T<:Number where S<:Number
    return 2*abs.(A-B)./abs.(A+B)
end

"""
    effreldiff(A, B, threshold=1e-14)

Same as `reldiff(A,B)` but with all elements set to zero where corresponding element of
`absdiff(A,B)` is smaller than `threshold`. This is useful in avoiding artificially large
relative errors.
"""
function effreldiff(
        A::AbstractArray{T}, B::AbstractArray{S}, threshold::Float64=1e-14
    ) where {T <: Number, S <: Number}
    r = reldiff(A,B)
    r[findall(x -> abs.(x)<threshold, absdiff(A,B))] .= 0.
    return r
end

"""
    absdiff(A, B)

Difference of absolute values of `A` and `B`.
"""
function absdiff(A::AbstractArray{T}, B::AbstractArray{S}) where {T<:Number, S<:Number}
    return abs.(A-B)
end

"""
    compare(A, B)

Compares two matrices `A` and `B`, prints out the maximal absolute and relative differences
and returns a boolean indicating wether `isapprox(A,B)`.
"""
function compare(A::AbstractArray{T}, B::AbstractArray{S}) where {T<:Number, S<:Number}
    @printf("max absdiff: %.1e\n", maximum(absdiff(A,B)))
    @printf("mean absdiff: %.1e\n", mean(absdiff(A,B)))
    @printf("max reldiff: %.1e\n", maximum(reldiff(A,B)))
    @printf("mean reldiff: %.1e\n", mean(reldiff(A,B)))

    r = effreldiff(A,B)
    @printf("effective max reldiff: %.1e\n", maximum(r))
    @printf("effective mean reldiff: %.1e\n", mean(r))

    return isapprox(A,B)
end

# NOTE
# Currenlty julia/sparearrays does not implement this function (type signature)
# once it does this can be removed/depracted in favor of mul!
# see also: test/slice_matrices.jl
occursin(
    "SparseArrays",
    string(which(mul!, (Matrix, Matrix, SparseMatrixCSC)).file)
) && @warn(
    "A Method `mul!(::Matrix, ::Matrix, ::SparseMatrixCSC)` now exists in " *
    "`SparseArrays`. The method defined in `helpers.jl` is likely to be  " *
    "unnecessary now."
)

function SparseArrays.mul!(C::StridedMatrix, X::StridedMatrix, A::SparseMatrixCSC)
    mX, nX = size(X)
    nX == A.m || throw(DimensionMismatch())
    fill!(C, zero(eltype(C)))
    rowval = A.rowval
    nzval = A.nzval
    @inbounds for  col = 1:A.n, k=A.colptr[col]:(A.colptr[col+1]-1)
        ki=rowval[k]
        kv=nzval[k]
        for multivec_row=1:mX
            C[multivec_row, col] += X[multivec_row, ki] * kv
        end
    end
    C
end

# Taken from Base
if !isdefined(Base, :splitpath)
    splitpath(p::AbstractString) = splitpath(String(p))

    if Sys.isunix()
        const path_dir_splitter = r"^(.*?)(/+)([^/]*)$"
    elseif Sys.iswindows()
        const path_dir_splitter = r"^(.*?)([/\\]+)([^/\\]*)$"
    else
        error("path primitives for this OS need to be defined")
    end

    _splitdir_nodrive(path::String) = _splitdir_nodrive("", path)
    function _splitdir_nodrive(a::String, b::String)
        m = match(path_dir_splitter,b)
        m === nothing && return (a,b)
        a = string(a, isempty(m.captures[1]) ? m.captures[2][1] : m.captures[1])
        a, String(m.captures[3])
    end

    function splitpath(p::String)
        drive, p = splitdrive(p)
        out = String[]
        isempty(p) && (pushfirst!(out,p))  # "" means the current directory.
        while !isempty(p)
            dir, base = _splitdir_nodrive(p)
            dir == p && (pushfirst!(out, dir); break)  # Reached root node.
            if !isempty(base)  # Skip trailing '/' in basename
                pushfirst!(out, base)
            end
            p = dir
        end
        if !isempty(drive)  # Tack the drive back on to the first element.
            out[1] = drive*out[1]  # Note that length(out) is always >= 1.
        end
        return out
    end
end

"""
    @bm function ... end
    @bm foo(args...) = ...
    @bm "name" begin ... end

Wraps the body of a function with `@timeit_debug <function name> begin ... end`.
This macro can also be used on a code block to generate
`@timeit_debug <name> begin ... end`.

The `@timeit_debug` macro can be disabled per module. Using `@bm` will make
sure that the module is always `MonteCarlo`. One can enable and disable
benchmarking with `enable_benchmarks()` and `disable_benchmarks()`. If they are
disabled they should come with zero overhead. See TimerOutputs.jl for more
details.

Benchmarks/Timings can be retrieved using `print_timer()` and reset with
`reset_timer!()`.
"""
macro bm(args...)
    if length(args) == 1 && args[1].head in (Symbol("="), :function)
        expr = args[1]
        code = TimerOutputs.timer_expr(
            MonteCarlo, true,
            _to_typed_name(expr.args[1]),
            :(begin $(expr.args[2]) end) # inner code block
        )
        Expr(
            expr.head,     # function or =
            esc(expr.args[1]),  # function name w/ args
            code
        )
    else
        # Not a function, just do the same as timeit_debug
        # This is copied from TimerOutputs.jl
        # With __module__ replaced by MonteCarlo because we want to have all
        # timings in the MonteCarlo namespace (otherwise they are not tied to
        # MonteCarlo.timeit_debug_enabled())
        TimerOutputs.timer_expr(MonteCarlo, true, args...)
    end
end

function _to_typed_name(e::Expr)
    if e.head == :where
        _to_typed_name(e.args[1])
    elseif e.head == :call
        string(e.args[1]) * "(" * join(_to_type.(e.args[2:end]), ", ") * ")"
    else
        dump(e)
        "ERROR"
    end
end

_to_type(s::Symbol) = "::Any"
function _to_type(e::Expr)
    if e.head == Symbol("::")
        return "::" * string(e.args[end])
    elseif e.head == Symbol("...")
        return _to_type(e.args[1]) * "..."
    elseif e.head == :parameters # keyword args
        return _to_type(e.args[1])
    elseif e.head == :kw # default value
        return _to_type(e.args[1])
    else
        dump(e)
        return "???"
    end
end

timeit_debug_enabled() = false

"""
    enable_benchmarks()

Enables benchmarking for `MonteCarlo`.

This affects every function with the `MonteCarlo.@bm` macro as well as any
`TimerOutputs.@timeit_debug` blocks. Benchmarks are recorded to the default
TimerOutput `TimerOutputs.DEFAULT_TIMER`. Results can be printed via
`TimerOutputs.print_timer()`.

[`disable_benchmarks`](@ref)
"""
enable_benchmarks() = TimerOutputs.enable_debug_timings(MonteCarlo)

"""
    disable_benchmarks()

Disables benchmarking for `MonteCarlo`.

This affects every function with the `MonteCarlo.@bm` macro as well as any
`TimerOutputs.@timeit_debug` blocks.

[`enable_benchmarks`](@ref)
"""
disable_benchmarks() = TimerOutputs.disable_debug_timings(MonteCarlo)


text:
# Spinnaker.jl: wrapper for FLIR/Point Grey Spinnaker SDK
# Copyright (C) 2019 Samuel Powell

module Spinnaker

using FixedPointNumbers

import Libdl
import Base: unsafe_convert, show, length, getindex, size, convert, range

export System, Camera, CameraList

# Include build configuration
try
  include(joinpath(@__DIR__, "..", "deps", "deps.jl"))
catch
  error("Package configuration file missing, run 'Pkg.build(\"Spinnaker\")' to configure.")
end

const MAX_BUFFER_LEN = Csize_t(1023)

# include API wrapper
include("wrapper/CEnum.jl")
using .CEnum

include("wrapper/spin_common.jl")

function checkerror(err::spinError)
  if err != spinError(0)
    throw(ErrorException("Spinnaker SDK error: $err"))
  end
  return nothing
end

include("wrapper/spin_api.jl")

# export everything spin*
foreach(names(@__MODULE__, all=true)) do s
  if startswith(string(s), "spin")
    @eval export $s
  end
end

# Include interface
include("SpinImage.jl")
include("CameraImage.jl")
include("System.jl")
include("Camera.jl")
include("CameraList.jl")
include("NodeMap.jl")
include("Nodes.jl")

# Create a System object at runtime
function __init__()
  global spinsys = System()
end

end # module


text:
"""
the eq_const is defined as the ratio of the volume fraction of the solvent in oil to
its volume fraction in water. For instance, For DME, the molar fraction is around 2.0
One criteria is that eq_const should probably be higher than 1.
"""
function water_soluble_solvent_flood(core_props, fluids_ls, fluids_hs, rel_perms_hs,
    rel_perms_ls, core_flood, eq_const)
    # construct the fractional flow curves
    fw_ls, dfw_ls = fractional_flow_function(rel_perms_ls, fluids_ls)
    fw_hs, dfw_hs = fractional_flow_function(rel_perms_hs, fluids_hs)
    # low sal shock (tangent line from (0,0))
    point1 = (-eq_const/(1-eq_const), -eq_const/(1-eq_const))
    sw_shock_ls = tangent_line_saturation(rel_perms_ls, fluids_ls, point1)
    println("low sal sw_shock = $sw_shock_ls")
    t_D_BT_ls = 1/dfw_ls(sw_shock_ls) # breakthrough (BT) time [#PV]
    println("low sal breakthrough time = $t_D_BT_ls")
    # High sal shock (cross point between the ls tangent and the hs fw)
    sw_shock_hs = cross_point_saturation(fw_hs, rel_perms_hs, point1, (sw_shock_ls, fw_ls(sw_shock_ls)))
    # println("high sal sw_shock = $sw_shock_hs")
    sw_init = core_flood.initial_water_saturation
    t_D_BT_hs = (sw_shock_hs-sw_init)/(fw_hs(sw_shock_hs)-fw_hs(sw_init)) # breakthrough (BT) time [#PV]
    # println("high sal breakthrough time = $t_D_BT_hs")

    # tracer
    sw_tracer = tangent_line_saturation(rel_perms_ls, fluids_ls, (0.0, 0.0))
    t_D_tracer = 1/dfw_ls(sw_tracer)
    xt_tracer_shock = dfw_ls(sw_tracer)

    # construct the recovery factor curve versus the # of PV
    # R = zeros(1)
    # pv_R = zeros(1)
    # at breakthrough of the hs brine
    # push!(R, (1-fw_hs(sw_init))*t_D_BT_hs/(1-sw_init)) # recovery at BT
    # push!(pv_R, t_D_BT_hs) # BT time
    # at breakthrough of the ls brine
    # push!(R, R[end]+(1-fw_hs(sw_shock_hs))*(t_D_BT_ls-t_D_BT_hs)/(1-sw_init)) # recovery at BT
    # push!(pv_R, t_D_BT_ls) # BT time

    # after breakthrough
    pv_inj = max(core_flood.injected_pore_volume, 2.0) # at least inject 2 pv
    # f_sw = sw -> (pv_inj-1/dfw_ls(sw)) # find the outlet saturation at the end of injection
    # sw_max = fzero(f_sw, sw_shock_ls)
    # sw_tmp = linspace(sw_shock_ls, sw_max, 100)
    # t_D_tmp = 1.0./dfw_ls.(sw_tmp)

    # s_av_tmp = sw_tmp-(fw_ls.(sw_tmp)-1).*t_D_tmp
    # R_tmp = (s_av_tmp-sw_init)/(1-sw_init)
    # append!(R, R_tmp)
    # append!(pv_R, t_D_tmp)

    # saturation profile
    sor = rel_perms_ls.sor
    sw_inj = core_flood.injected_water_saturation
    phi = core_props.porosity
    ut = core_flood.injection_velocity
    L = core_props.length
    pv_to_t = phi*L/ut
    s1 = collect(linspace(min(sw_inj, 1-sor-eps()), sw_shock_ls, 300))
    xt_s1 = dfw_ls.(s1)
    xt_shock_ls = 1/t_D_BT_ls
    xt_shock_hs = 1/t_D_BT_hs
    xt_prf=[xt_s1; xt_shock_ls+eps(); xt_shock_hs; xt_shock_hs+eps(); 1/0.3]
    sw_prf=[s1; sw_shock_hs; sw_shock_hs; sw_init; sw_init]
    xt_tracer = [0.0, xt_tracer_shock, xt_tracer_shock+eps(), xt_prf[end]]
    c_tracer = [1.0, 1.0, 0.0, 0.0]
    # return pv_R, R, xt_prf, sw_prf # for the time being to test the code

    # extract data for calculating the pressure drop
    k = core_props.permeability
    krw_ls, kro_ls, dkrw_ls, dkro_ls = rel_perm_functions(rel_perms_ls)
    muo_ls, muw_ls = fluids_ls.oil_viscosity, fluids_ls.water_viscosity
    krw_hs, kro_hs, dkrw_hs, dkro_hs = rel_perm_functions(rel_perms_hs)
    muo_hs, muw_hs = fluids_hs.oil_viscosity, fluids_hs.water_viscosity
    λ_ls = s -> krw_ls(s)/muw_ls+kro_ls(s)/muo_ls
    λ_hs = s -> krw_hs(s)/muw_hs+kro_hs(s)/muo_hs
    λ_w_ls = s -> krw_ls(s)/muw_ls
    λ_w_hs = s -> krw_hs(s)/muw_hs
    λ_t = [λ_ls.(s1); λ_hs(sw_shock_hs); λ_hs(sw_shock_hs); λ_hs(sw_init); λ_hs(sw_init)]
    λ_w = [λ_w_ls.(s1); λ_w_hs(sw_shock_hs); λ_w_hs(sw_shock_hs); λ_w_hs(sw_init); λ_w_hs(sw_init)]

    xt = copy(xt_prf)
    sw = copy(sw_prf)
    i=1
    while(true)
        if i==length(xt)
            break
        elseif xt[i]>=xt[i+1]
            deleteat!(xt, i)
            deleteat!(sw, i)
            deleteat!(λ_t, i)
            deleteat!(λ_w, i)
        else
            i+=1
        end
    end

    x = collect(linspace(0,L,200))
    sw_int = Spline1D(xt, sw, k=1, bc="nearest")
    λ_int = Spline1D(xt, λ_t, k=1, bc="nearest")
    λ_w_int = Spline1D(xt, λ_w, k=1, bc="nearest")
    t_inj=pv_inj*pv_to_t
    t = collect([linspace(0,t_D_BT_hs*pv_to_t, 10);
                linspace((t_D_BT_hs+eps())*pv_to_t, t_D_BT_ls*pv_to_t, 40);
                linspace((t_D_BT_ls+eps())*pv_to_t,t_inj, 100)]) # [s] time
    println("co2 breakthrough time is $(t_D_BT_ls*pv_to_t)")
    p_inj = zeros(length(t))
    R_int = zeros(length(t))
    WC_int = zeros(length(t))
    p_inj[1]=L*ut./(k*λ_hs(sw_init))
    WC_int[1]=λ_w_hs(sw_init)/λ_hs(sw_init)
    for i in 2:length(t)
        xt_real = x/t[i]/(ut/phi)
        p_inj[i] = trapz(x, ut./(k*λ_int(xt_real)))
        R_int[i] = (trapz(x, sw_int(xt_real))/L-sw_init)/(1-sw_init)
        WC_int[i] = λ_w_int(xt_real[end])/λ_int(xt_real[end])
    end

    R_int[R_int.<0.0] .= 0.0

    return FracFlowResults([fw_hs, fw_ls], [Line([-eq_const/(1-eq_const), -eq_const/(1-eq_const)],
                            [sw_shock_hs, fw_hs(sw_shock_hs)]),
                            Line([sw_init, fw_hs(sw_init)], [sw_shock_hs, fw_hs(sw_shock_hs)])],
                            [t/pv_to_t R_int], [t R_int], [xt_prf sw_prf], [xt_tracer c_tracer],
                            [t/pv_to_t p_inj], [t p_inj], [t/pv_to_t WC_int], [t WC_int])

end

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
# Autogenerated wrapper script for HGSCVRP_jll for x86_64-linux-gnu-cxx03
export hgs, libhgscvrp

JLLWrappers.@generate_wrapper_header("HGSCVRP")
JLLWrappers.@declare_library_product(libhgscvrp, "libhgscvrp.so")
JLLWrappers.@declare_executable_product(hgs)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libhgscvrp,
        "lib/libhgscvrp.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        hgs,
        "bin/hgs",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "GetData functions tests" begin

datadir = "./sr28/"

@testset "getfooddes Tests" begin
  @test isdefined(FoodRecom, :getfooddes) == true
  fooddes = FoodRecom.getfooddes(datadir)
  @inferred FoodRecom.getfooddes(datadir)
  @test size(fooddes) ==  (8789, 5)
  colnames = [:food, :group, :longdes, :shortdes, :othernames]
  @test names(fooddes) == colnames
end

@testset "getnutdata Tests" begin
  @test isdefined(FoodRecom, :getnutdata) == true
  nutdata = FoodRecom.getnutdata(datadir)
  @inferred FoodRecom.getnutdata(datadir)
  @test size(nutdata) == (100500, 4)
  colnames = [:food, :nut, :amount, :stderr]
  @test names(nutdata) == colnames
end

@testset "getnutrdef Tests" begin
  @test isdefined(FoodRecom, :getnutrdef) == true
  nutrdef = FoodRecom.getnutrdef(datadir)
  @inferred FoodRecom.getnutrdef(datadir)
  @test size(nutrdef) == (150, 4)
  colnames = [:nut, :unit, :nuttag, :nutname]
  @test names(nutrdef) == colnames
  @test nutr_def[31, 2] == "ug"  # successfully converted UTF-8 values
end

@testset "allnutdata Tests" begin
  @test isdefined(FoodRecom, :allnutdata) == true
  nutrdef = FoodRecom.allnutdata(data_folder=datadir)
  @inferred FoodRecom.allnutdata(data_folder=datadir)
  @test size(nutrdef) == (100500, 11)
end

end


text:
module BoundaryValueDiffEq

using Reexport, LinearAlgebra, SparseArrays
@reexport using DiffEqBase

import DiffEqBase: solve
import NLsolve, ForwardDiff, BandedMatrices, FiniteDiff

struct MIRKTableau{T}
    c::Vector{T}
    v::Vector{T}
    b::Vector{T}
    x::Matrix{T}
end

# ODE BVP problem system
mutable struct BVPSystem{T,U<:AbstractArray,P}
    order::Int                  # The order of MIRK method
    M::Int                      # Number of equations in the ODE system
    N::Int                      # Number of nodes in the mesh
    fun!                        # M -> M
    bc!                         # 2 -> 2
    p::P
    x::Vector{T}                # N
    y::Vector{U}                # N{M}
    f::Vector{U}                # N{M}
    residual::Vector{U}         # N{M}
    tmp::Vector{T}
end

include("vector_auxiliary.jl")
include("algorithms.jl")
include("alg_utils.jl")
include("mirk_tableaus.jl")
include("cache.jl")
include("collocation.jl")
include("jacobian.jl")
include("solve.jl")

export Shooting
export MIRK4, GeneralMIRK4

end # module


text:
@inline function C_e(CellData,s,z)
"""
Electrolyte Concentration Transfer Function
# Add Ins and Outs
    # Locations to be computed
    # Sampling Frequency
"""

ζ = (1-CellData.Const.t_plus)/F    #Simplifying Variable
CC_A = CellData.Const.CC_A   # Current-collector area [m^2]
κ_eff_Neg = CellData.Const.κ*CellData.Neg.ϵ_e^CellData.Neg.κ_brug
κ_eff_Pos = CellData.Const.κ*CellData.Pos.ϵ_e^CellData.Pos.κ_brug
σ_eff_Neg = CellData.Neg.σ*CellData.Neg.ϵ_s^CellData.Neg.σ_brug #Effective Conductivity Neg
σ_eff_Pos = CellData.Pos.σ*CellData.Pos.ϵ_s^CellData.Pos.σ_brug #Effective Conductivity Pos

#Defining SOC
θ_neg = CellData.Const.SOC * (CellData.Neg.θ_100-CellData.Neg.θ_0) + CellData.Neg.θ_0
θ_pos = CellData.Const.SOC * (CellData.Pos.θ_100-CellData.Pos.θ_0) + CellData.Pos.θ_0

#Beta's
βn = @. CellData.Neg.Rs*sqrt(s/CellData.Neg.Ds)
βp = @. CellData.Pos.Rs*sqrt(s/CellData.Pos.Ds)

#Prepare for j0
ce0 = CellData.Const.ce0
cs_max_neg = CellData.Neg.cs_max
cs0_neg = cs_max_neg * θ_neg
cs_max_pos = CellData.Pos.cs_max
cs0_pos = cs_max_pos * θ_pos

#Current Flux Density

if CellData.Const.CellTyp == "Doyle_94"
    κ_pos = CellData.Pos.k_norm/(CellData.Pos.cs_max*ce0^(1-CellData.Pos.α))
    κ_neg = CellData.Neg.k_norm/CellData.Neg.cs_max/ce0^(1-CellData.Neg.α)
    j0_neg = κ_neg*(ce0*(cs_max_neg-cs0_neg))^(1-CellData.Neg.α)*cs0_neg^CellData.Neg.α
    j0_pos = κ_pos*(ce0*(cs_max_pos-cs0_pos))^(1-CellData.Pos.α)*cs0_pos^CellData.Pos.α
 else
    j0_neg = CellData.Neg.k_norm*(ce0*(cs0_neg/cs_max_neg*(1-cs0_neg/cs_max_neg)))^(1-CellData.Neg.α)
    j0_pos = CellData.Pos.k_norm*(ce0*(cs0_pos/cs_max_pos*(1-cs0_pos/cs_max_pos)))^(1-CellData.Pos.α)
 end

#Resistances
Rtot_neg = R*CellData.Const.T/(j0_neg*F^2) + CellData.Neg.RFilm
Rtot_pos = R*CellData.Const.T/(j0_pos*F^2) + CellData.Pos.RFilm

#OCP derivative
∂Uocp_pos = CellData.Const.∂Uocp("Pos",θ_pos)/cs_max_pos
∂Uocp_neg = CellData.Const.∂Uocp("Neg",θ_neg)/cs_max_neg

#Condensing Variable
ν_n =  @. CellData.Neg.L*sqrt((CellData.Neg.as/σ_eff_Neg+CellData.Neg.as/κ_eff_Neg)/(Rtot_neg+∂Uocp_neg*(CellData.Neg.Rs/(F*CellData.Neg.Ds))*(tanh(βn)/(tanh(βn)-βn))))
ν_p =  @. CellData.Pos.L*sqrt((CellData.Pos.as/σ_eff_Pos+CellData.Pos.as/κ_eff_Pos)/(Rtot_pos+∂Uocp_pos*(CellData.Pos.Rs/(F*CellData.Pos.Ds))*(tanh(βp)/(tanh(βp)-βp))))

R_ce = roots(CellData,CellData.Const.Ce_M+1)
λ = R_ce[2:end]

#Create all k's
in1 = @. sqrt(λ*CellData.Neg.ϵ_e/CellData.Const.D1)
in2 = @. sqrt(λ*CellData.Sep.ϵ_e/CellData.Const.D2)
in3 = @. sqrt(λ*CellData.Pos.ϵ_e/CellData.Const.D3)

Bound_Neg_1 = CellData.Neg.L * in1
Bound_Sep_0 = CellData.Neg.L * in2
Bound_Sep_1 = (CellData.Neg.L+CellData.Sep.L) * in2
Bound_Pos_0 = (CellData.Neg.L+CellData.Sep.L) * in3
Bound_Pos_1 = (CellData.Neg.L+CellData.Sep.L+CellData.Pos.L) * in3
Bound_Pos_2 = CellData.Pos.L * in3

#Scaled coefficients
k3_s = @. cos(Bound_Neg_1)*cos(Bound_Sep_0) + CellData.Const.D1*in1*sin(Bound_Neg_1)*sin(Bound_Sep_0)/(CellData.Const.D2*in2)
k4_s = @. cos(Bound_Neg_1)*sin(Bound_Sep_0) - CellData.Const.D1*in1*cos(Bound_Sep_0)*sin(Bound_Neg_1)/(CellData.Const.D2*in2)
k5_s = @. k3_s*(cos(Bound_Sep_1)*cos(Bound_Pos_0)+CellData.Const.D2*in2*sin(Bound_Sep_1)*sin(Bound_Pos_0)/(CellData.Const.D3*in3))+k4_s*(sin(Bound_Sep_1)*cos(Bound_Pos_0)-CellData.Const.D2*in2*cos(Bound_Sep_1)*sin(Bound_Pos_0)/(CellData.Const.D3*in3))
k6_s = @. k3_s*(cos(Bound_Sep_1)*sin(Bound_Pos_0)-CellData.Const.D2*in2*sin(Bound_Sep_1)*cos(Bound_Pos_0)/(CellData.Const.D3*in3))+k4_s*(sin(Bound_Sep_1)*sin(Bound_Pos_0)+CellData.Const.D2*in2*cos(Bound_Sep_1)*cos(Bound_Pos_0)/(CellData.Const.D3*in3))

#Solving for k1:
Int_ψ1 = @. CellData.Neg.ϵ_e*(2*Bound_Neg_1+sin(2*Bound_Neg_1))/(4*in1)
Int_ψ2 = @. CellData.Sep.ϵ_e/(4*in2)*(2*(k3_s^2+k4_s^2)*CellData.Sep.L*in2+2*k3_s*k4_s*cos(2*Bound_Sep_0)-2*k3_s*k4_s*cos(2*Bound_Sep_1)-(k3_s-k4_s)*(k3_s+k4_s)*(sin(2*Bound_Sep_0)-sin(2*Bound_Sep_1)))
Int_ψ3 = @. CellData.Pos.ϵ_e/(4*in3)*(2*(k5_s^2+k6_s^2)*CellData.Pos.L*in3+2*k5_s*k6_s*cos(2*Bound_Pos_0)-2*k5_s*k6_s*cos(2*Bound_Pos_1)-(k5_s-k6_s)*(k5_s+k6_s)*(sin(2*Bound_Pos_0)-sin(2*Bound_Pos_1)))

k1 = @. 1/(sqrt(Int_ψ1+Int_ψ2+Int_ψ3))
k3 = @. k1*k3_s
k4 = @. k1*k4_s
k5 = @. k1*k5_s
k6 = @. k1*k6_s

j_Neg = @. k1*ζ*ν_n*(Bound_Neg_1*(κ_eff_Neg+σ_eff_Neg*cosh(ν_n))*sin(Bound_Neg_1)+(κ_eff_Neg+σ_eff_Neg*cos(Bound_Neg_1))*sinh(ν_n)*ν_n)/(CC_A*(κ_eff_Neg+σ_eff_Neg)*(Bound_Neg_1^2+ν_n^2)*sinh(ν_n))
zero_tf_neg = @. k1*ζ*sin(Bound_Neg_1)/(CC_A*Bound_Neg_1)
j_Neg[:,findall(s.==0)] .= zero_tf_neg[:,findall(s.==0)]

j_Pos = @. -ζ*ν_p/(CC_A*(κ_eff_Pos+σ_eff_Pos)*(Bound_Pos_2^2+ν_p^2)*sinh(ν_p))*(-k6*Bound_Pos_2*cos(Bound_Pos_1)*(σ_eff_Pos+κ_eff_Pos*cosh(ν_p))+Bound_Pos_2*(κ_eff_Pos+σ_eff_Pos*cosh(ν_p))*(k6*cos(Bound_Pos_0)-k5*sin(Bound_Pos_0))+k5*Bound_Pos_2*(σ_eff_Pos+ κ_eff_Pos*cosh(ν_p))*sin(Bound_Pos_1)+sinh(ν_p)*(k5*σ_eff_Pos*cos(Bound_Pos_0)+k5*κ_eff_Pos*cos(Bound_Pos_1)+k6*σ_eff_Pos*sin(Bound_Pos_0)+k6*κ_eff_Pos*sin(Bound_Pos_1))*ν_p)
zero_tf = @. -ζ*(k6*(cos(Bound_Pos_0)-cos(Bound_Pos_1))+k5*(sin(Bound_Pos_1)-sin(Bound_Pos_0)))/(CC_A*Bound_Pos_2)
j_Pos[:,findall(s.==0)] .= zero_tf[:,findall(s.==0)]

C_e =  @. ((j_Neg + j_Pos)/(s+λ))
i=Int64(1)
ψ = fill(0.0,length(z),length(λ))
for loop in 1:length(λ)
    i=1
    for x in z #Eigen Weighting
        if x < CellData.Neg.L+eps()
            ψ[i,loop] = k1[loop]*cos(in1[loop]*x) #negative electrode
        elseif x > (CellData.Neg.L+CellData.Sep.L)-eps()
            ψ[i,loop] = k5[loop]*cos(in3[loop]*x)+k6[loop]*sin(in3[loop]*x) # postive electrode
        else
            ψ[i,loop] = k3[loop]*cos(in2[loop]*x)+k4[loop]*sin(in2[loop]*x) # separator
        end
    i = i+1
    end
end

C_e = ψ*C_e
D =  zeros(length(z))
D_term = "zeros(length($z))"
res0 = zeros(length(z))

return C_e, D, res0, D_term
end

function roots(CellData,roots_n) #Change these functions to be inclusive of CellData struct -> one input for roots
    root = Float64[0]
    i = 1e-5
    ∇ = 1e-5
    if(roots_n > 1)
        while length(root) <= roots_n-1
            if flambda(CellData,(i-∇))*flambda(CellData,i)<0
            push!(root, _bisection(flambda,CellData,(i-∇),i,(i-∇),i,√eps(),0,100)[1])
            end
        i = i+∇
        end
    end
    return root
end

function flambda(CellData,λ)
    ϵ1 = CellData.Neg.ϵ_e
    ϵ2 = CellData.Sep.ϵ_e
    ϵ3 = CellData.Pos.ϵ_e
    D1 = CellData.Const.D1
    D2 = CellData.Const.D2
    D3 = CellData.Const.D3
    Lneg = CellData.Neg.L
    Lnegsep = CellData.Const.Lnegsep
    Ltot = CellData.Neg.L+CellData.Sep.L+CellData.Pos.L
    k1 = Float64(1)
    sle1 = sqrt(λ*ϵ1/D1)
    sle2 = sqrt(λ*ϵ2/D2)
    sle3 = sqrt(λ*ϵ3/D3)
    k3 = k1*(cos(sle1*Lneg).*cos(sle2*Lneg) + D1*sle1.*sin(sle1*Lneg).*sin(sle2*Lneg)./(D2*sle2))
    k4 = k1*(cos(sle1*Lneg).*sin(sle2*Lneg) - D1*sle1.*cos(sle2*Lneg).*sin(sle1*Lneg)./(D2*sle2))
    k5 = k3*(cos(sle2*Lnegsep).*cos(sle3*Lnegsep) + D2*sle2.*sin(sle2*Lnegsep).*sin(sle3*Lnegsep)./(D3*sle3))+k4*(sin(sle2*Lnegsep).*cos(sle3*Lnegsep) - D2*sle2.*cos(sle2*Lnegsep).*sin(sle3*Lnegsep)./(D3*sle3))
    k6 = k3*(cos(sle2*Lnegsep).*sin(sle3*Lnegsep) - D2*sle2.*sin(sle2*Lnegsep).*cos(sle3*Lnegsep)./(D3*sle3))+k4*(sin(sle2*Lnegsep).*sin(sle3*Lnegsep) + D2*sle2.*cos(sle2*Lnegsep).*cos(sle3*Lnegsep)./(D3*sle3))
    Psiprime = -k5.*sle3.*sin(sle3*Ltot) + k6.*sle3.*cos(sle3*Ltot)
    return Psiprime
end

max_stars_repo_name: BradyPlanden/RealisationAlgorithms

max_issues_repo_name: BradyPlanden/RealisationAlgorithms

max_forks_repo_name: BradyPlanden/RealisationAlgorithms


text:
using BinaryBuilder

name = "SLATE"
version = v"2022.05.00"

# Collection of sources required to build PETSc. Avoid using the git repository, it will
# require building SOWING which fails in all non-linux platforms.
sources = [
    GitSource("https://bitbucket.org/icl/slate.git", "5fb57877effa06d2ef090402a39341ebeb44f180")
]

# Bash recipe for building across all platforms

# Needs to add -Dcapi eventually once it's added to the cmake build system. Note yet available under CMAKAE toolchain.
script = raw"""
cd slate
git submodule update --init
mkdir build && cd build
cmake \
  -DCMAKE_INSTALL_PREFIX=${prefix} \
  -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
  -DCMAKE_BUILD_TYPE="Release" \
  -Dblas=openblas \
  -Dbuild_tests=no \
  -DMPI_RUN_RESULT_CXX_libver_mpi_normal="0" \
  -DMPI_RUN_RESULT_CXX_libver_mpi_normal__TRYRUN_OUTPUT="" \
  -Drun_result="0" \
  -Drun_result__TRYRUN_OUTPUT="ok" \
  -Dc_api=yes \
  ..
make -j${nproc}
make install
"""

# We attempt to build for all defined platforms
platforms = expand_gfortran_versions(expand_cxxstring_abis(supported_platforms(; exclude=!Sys.islinux)))
platforms = filter(p -> libgfortran_version(p) ≠ v"3", platforms)
products = [
    LibraryProduct("libslate", :libslate),
    LibraryProduct("libslate_lapack_api", :libslate_lapack_api)
    # LibraryProduct("libslate_scalapack_api, :libslate_scalapack_api) ** Not yet available under CMAKE toolchain.
]

dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("OpenBLAS32_jll"),
    Dependency("MPICH_jll")
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"7")


text:
## ExoplanetsSysSim/src/transit_detection_model.jl
## (c) 2015 Eric B. Ford

# Several functions below based on https://github.com/christopherburke/KeplerPORTs/blob/master/KeplerPORTs_utils.py
# That follows the procedure outlined in Burke et al.(2015).
# However we don't currently interpolate the CDDP or mesthreshold to the relevant duration

function real_log_choose(m::Float64, n::Float64)::Float64
  lgamma(m+1)-lgamma(n+1)-lgamma(m-n+1.0)
end

function real_binom(k::Float64, BigM::Float64, f::Float64)::Float64
    F1 = real_log_choose(BigM,k)
    F2 = k*log(f)
    F3 = (BigM-k)*log(1.0-f)
    x = exp(F1+F2+F3)
    return x
end

function kepler_window_function_binomial_model(exp_num_transits_no_gaps::Float64, duty_cycle::Float64; min_transits::Float64 = 3.0)::Float64
  if exp_num_transits_no_gaps < min_transits
     return 0.0
  else
     return max(1.0 - real_binom(min_transits,exp_num_transits_no_gaps,duty_cycle), 0.0)
  end
end

function kepler_window_function_binomial_model(t::KeplerTarget, exp_num_transits_no_gaps::Float64, period::Float64, duration::Float64; min_transits::Float64 = 3.0)::Float64
   kepler_window_function_binomial_model(exp_num_transits_no_gaps, t.duty_cycle, min_transits=min_transits)
end

function kepler_window_function_dr25_model(t::KeplerTarget, exp_num_transits_no_gaps::Float64, period::Float64, duration::Float64)::Float64
   ExoplanetsSysSim.WindowFunction.eval_window_function(t.window_function_id, Duration=duration, Period=period)
end

#kepler_window_function = kepler_window_function_binomial_model
kepler_window_function = kepler_window_function_dr25_model

# See Christiansen et al. (2015)  This assumes a linear limbdarkening coefficient of 0.6
function frac_depth_to_tps_depth(frac_depth::Float64)
    const alp = 1.0874
    const bet = 1.0187
    const REALDEPTH2TPSSQUARE = 1.0  # WARNING: Waiting for this to be confirmed
    k = sqrt(frac_depth)
    tps_depth = min( (alp-bet*k) * frac_depth* REALDEPTH2TPSSQUARE, 1.0)   # NOTE: I added the max based on common sense
    return tps_depth::Float64
end

function detection_efficiency_theory(mes::Float64, expected_num_transits::Float64; min_pdet_nonzero::Float64 = 0.0)
   const muoffset =  0.0
   const sig =  1.0
   const mesthresh = 7.1
   mes *= 1.003
   if mes > (9.0 - mesthresh - muoffset)
      return 1.0
   else
      pdet = 0.5 + 0.5*erf((mes - mesthresh - muoffset) / sqrt(2.0*sig*sig))
      pdet = pdet >= min_pdet_nonzero ? pdet : 0.0
      return pdet
   end
end

function detection_efficiency_fressin2013(mes::Float64, expected_num_transits::Float64)
   const mesmin =  6.0
   const mesmax =  16.0
   mes *= 1.003
   if mes <= mesmin
      return 0.0
   elseif mes >= mesmax
      return 1.0
   else
     return (mes - mesmin) / (mesmax - mesmin)
  end
end

function detection_efficiency_christiansen2015(mes::Float64, expected_num_transits::Float64; mes_threshold::Float64 = 7.1, min_pdet_nonzero::Float64 = 0.0)
   const a =  4.65  # from code for detection_efficiency(...) at https://github.com/christopherburke/KeplerPORTs/blob/master/KeplerPORTs_utils.py
   const b =  0.98
   #const a =  4.35 # from arxiv abstract. Informal testing showed it didn't matter
   #const b =  1.05
   mes *= 1.003
   usemes = max(0.0,mes - 7.1 - (mes_threshold - 7.1))
   pdet = cdf(Gamma(a,b), usemes)
   pdet = pdet >= min_pdet_nonzero ? pdet : 0.0
   return pdet
end

#global const detection_efficiency_dr25_simple_a =  30.87  # from pg 16 of https://exoplanetarchive.ipac.caltech.edu/docs/KSCI-19110-001.pdf
#global const detection_efficiency_dr25_simple_b =  0.271
#global const detection_efficiency_dr25_simple_dist = Gamma(detection_efficiency_dr25_simple_a,detection_efficiency_dr25_simple_b)
function detection_efficiency_dr25_simple(mes::Float64, expected_num_transits::Float64; min_pdet_nonzero::Float64 = 0.0)::Float64
   const a = 30.87  # from pg 16 of https://exoplanetarchive.ipac.caltech.edu/docs/KSCI-19110-001.pdf
   const b = 0.271
   const c = 0.940
   const dist = Gamma(a,b)
   mes *= 1.003
   #pdet = c*cdf(Gamma(a,b), mes)
   pdet::Float64 = c*cdf(dist, mes)::Float64
   #pdet::Float64 = c*cdf(detection_efficiency_dr25_simple_dist, mes)::Float64
   pdet = pdet >= min_pdet_nonzero ? pdet : 0.0
   return pdet
end

function get_param_for_fgk_detection_and_vetting_efficiency_depending_on_num_transits(num_tr::Integer)
    if num_tr <= 3
        return (33.3884, 0.264472, 0.699093)
    elseif num_tr <= 4
        return  (32.886, 0.269577, 0.768366)
    elseif num_tr <= 5
        return (31.5196, 0.282741, 0.833673)
    elseif num_tr <= 6
        return (30.9919, 0.286979, 0.859865)
    elseif num_tr <= 9
        return (30.1906, 0.294688, 0.875042)
    elseif num_tr <= 18
        return (31.6342, 0.279425, 0.886144)
    elseif num_tr <= 36
        return (32.6448, 0.268898, 0.889724)
    else
        return (27.8185, 0.32432, 0.945075)
    end
end

function get_param_for_m_detection_and_vetting_efficiency_depending_on_num_transits(num_tr::Integer)
    if num_tr <= 3
        return (33.3884, 0.264472, 0.699093)
    elseif num_tr <= 4
        return  (32.886, 0.269577, 0.768366)
    elseif num_tr <= 5
        return (31.5196, 0.282741, 0.833673)
    elseif num_tr <= 6
        return (30.9919, 0.286979, 0.859865)
    elseif num_tr <= 9
        return (30.1906, 0.294688, 0.875042)
    elseif num_tr <= 18
        return (30.39071943, 0.2933492,  0.89018648)
    elseif num_tr <= 36
        return (29.59808918, 0.29182242, 0.90915166)
    else
        return (31.25980585, 0.30160317, 0.86708737)
    end
end

# WARNING: Combined detection and vetting efficiency model - do NOT include additional vetting efficiency
function detection_and_vetting_efficiency_model_v1(mes::Float64, expected_num_transits::Float64; min_pdet_nonzero::Float64 = 0.0)::Float64
    mes *= 1.003
    num_transit_int = convert(Int64,floor(expected_num_transits))
    num_transit_int += rand() < expected_num_transits-num_transit_int ? 1 : 0
    a, b, c = get_param_for_fgk_detection_and_vetting_efficiency_depending_on_num_transits(num_transit_int)
    dist = Gamma(a,b)
    pdet::Float64 = c*cdf(dist, mes)::Float64
    pdet = pdet >= min_pdet_nonzero ? pdet : 0.0
    return pdet
end

# WARNING: Hardcoded choice of transit detection efficiency here for speed and so as to not have it hardcoded in multiple places
#detection_efficiency_model = detection_efficiency_christiansen2015
# detection_efficiency_model = detection_efficiency_dr25_simple
detection_efficiency_model = detection_and_vetting_efficiency_model_v1

function vetting_efficiency_none(R_p::Real, P::Real)
    return 1.0
end

# From Mulders e-mail (Gaia, no score)
# WARNING: Assumes inputs are in R_sol and days
function vetting_efficiency_dr25_mulders(R_p::Real, P::Real)
    const c = 0.93
    const a_R = -0.03
    const P_break = 205.
    const a_P = 0.00
    const b_P = -0.24
    pvet = c*(R_p/earth_radius)^a_R
    if P < P_break
        pvet *= (P/P_break)^a_P
    else
        pvet *= (P/P_break)^b_P
    end
    return pvet
end

# From Mulders et al. 2018 (arXiv 1805.08211)
# WARNING: Assumes inputs are in R_sol and days
function vetting_efficiency_dr25_mulders_vetcut(R_p::Real, P::Real)
    const c = 0.63
    const a_R = 0.19
    const P_break = 53.
    const a_P = -0.07
    const b_P = -0.39
    pvet = c*(R_p/earth_radius)^a_R
    if P < P_break
        pvet *= (P/P_break)^a_P
    else
        pvet *= (P/P_break)^b_P
    end
    return pvet
end

# Resume code original to SysSim

function interpolate_cdpp_to_duration(t::KeplerTarget, duration::Float64)::Float64
   duration_in_hours = duration *24.0
   dur_idx = searchsortedlast(cdpp_durations,duration_in_hours)   # cdpp_durations is defined in constants.jl
   if dur_idx <= 0
      cdpp = t.cdpp[1,1]
   elseif dur_idx==length(cdpp_durations) && (duration_in_hours >= cdpp_durations[end]) # Should be 15 cdpp_durations.
      cdpp = t.cdpp[length(cdpp_durations),1]
   else
      w = ((duration_in_hours)-cdpp_durations[dur_idx]) / (cdpp_durations[dur_idx+1]-cdpp_durations[dur_idx])
      cdpp = w*t.cdpp[dur_idx+1,1] + (1-w)*t.cdpp[dur_idx,1]
   end
   return cdpp
end

#this function calculates snr using cdpp instead of osd
function calc_snr_if_transit_cdpp(t::KeplerTarget, depth::Float64, duration::Float64, cdpp::Float64, sim_param::SimParam; num_transit::Float64 = 1)
  # depth_tps = frac_depth_to_tps_depth(depth)                  # TODO SCI:  WARNING: Hardcoded this conversion.  Remove once depth calculated using limb darkening model
  # snr = depth_tps*sqrt(num_transit*duration*LC_rate)/cdpp     # WARNING: Assumes measurement uncertainties are uncorrelated & CDPP based on LC
  snr = depth*sqrt(num_transit*duration*LC_rate)/cdpp     # WARNING: Assumes measurement uncertainties are uncorrelated & CDPP based on LC
end

function calc_snr_if_transit(t::KeplerTarget, depth::Real, duration::Real, osd::Real, sim_param::SimParam; num_transit::Real = 1)
  # depth_tps = frac_depth_to_tps_depth(depth)                  # WARNING: Hardcoded this conversion
  # snr = depth_tps/osd*1.0e6 # osd is in ppm
  snr = depth/osd*1.0e6 # osd is in ppm
end

function calc_snr_if_transit_central(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  depth = calc_transit_depth(t,s,p)
  duration_central = calc_transit_duration_central(t,s,p)
  kepid = StellarTable.star_table(t.sys[s].star.id, :kepid)
  osd_duration_central = get_legal_durations(period,duration_central)	#tests if durations are included in Kepler's observations for a certain planet period. If not, returns nearest possible duration
  osd_central = WindowFunction.interp_OSD_from_table(kepid, period, osd_duration_central)
  if osd_duration_central > duration_central				#use a correcting factor if this duration is lower than the minimum searched for this period.
      osd_central = osd_central*osd_duration_central/duration_central
  end
  num_transit = calc_expected_num_transits(t,s,p,sim_param)
  calc_snr_if_transit(t,depth,duration_central,osd_central, sim_param,num_transit=num_transit)
end

function calc_snr_if_transit_central_cdpp(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  depth = calc_transit_depth(t,s,p)
  duration_central = calc_transit_duration_central(t,s,p)
  cdpp = interpolate_cdpp_to_duration(t, duration_central)
  num_transit = calc_expected_num_transits(t,s,p,sim_param)
  calc_snr_if_transit_cdpp(t,depth,duration_central,cdpp, sim_param,num_transit=num_transit)
end

function calc_prob_detect_if_transit(t::KeplerTarget, snr::Float64, period::Float64, duration::Float64, sim_param::SimParam; num_transit::Float64 = 1)
  const min_pdet_nonzero = 1.0e-4                                                # TODO OPT: Consider raising threshold to prevent a plethora of planets that are very unlikely to be detected due to using 0.0 or other small value here
  wf = kepler_window_function(t, num_transit, period, duration)
  return wf*detection_efficiency_model(snr, num_transit, min_pdet_nonzero=min_pdet_nonzero)
end

function calc_prob_detect_if_transit(t::KeplerTarget, depth::Float64, period::Float64, duration::Float64, osd::Float64, sim_param::SimParam; num_transit::Float64 = 1)
  snr = calc_snr_if_transit(t,depth,duration,osd, sim_param, num_transit=num_transit)
  return calc_prob_detect_if_transit(t, snr, period, duration, sim_param, num_transit=num_transit)
end

function calc_prob_detect_if_transit_cdpp(t::KeplerTarget, depth::Float64, period::Float64, duration::Float64, cdpp::Float64, sim_param::SimParam; num_transit::Float64 = 1)
  snr = calc_snr_if_transit_cdpp(t,depth,duration,cdpp, sim_param, num_transit=num_transit)
  return calc_prob_detect_if_transit(t, snr, period, duration, sim_param, num_transit=num_transit)
end

function calc_prob_detect_if_transit_central(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  period = t.sys[s].orbit[p].P
  depth = calc_transit_depth(t,s,p)
  duration_central = calc_transit_duration_central(t,s,p)
  kepid = StellarTable.star_table(t.sys[s].star.id, :kepid)
  osd_duration_central = get_legal_durations(period,duration_central)	#tests if durations are included in Kepler's observations for a certain planet period. If not, returns nearest possible duration
  osd_central = WindowFunction.interp_OSD_from_table(kepid, period, osd_duration_central)
  if osd_duration_central > duration_central				#use a correcting factor if this duration is lower than the minimum searched for this period.
      osd_central = osd_central*osd_duration_central/duration_central
  end
  ntr = calc_expected_num_transits(t,s,p,sim_param)
  calc_prob_detect_if_transit(t,depth,period,duration_central,osd_central, sim_param, num_transit=ntr)
end

function calc_prob_detect_if_transit_central_cdpp(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  period = t.sys[s].orbit[p].P
  depth = calc_transit_depth(t,s,p)
  duration_central = calc_transit_duration_central(t,s,p)
  cdpp = interpolate_cdpp_to_duration(t, duration_central)
  ntr = calc_expected_num_transits(t,s,p,sim_param)
  calc_prob_detect_if_transit_cdpp(t,depth,period,duration_central,cdpp, sim_param, num_transit=ntr)
end

function calc_prob_detect_if_transit_with_actual_b(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  period = t.sys[s].orbit[p].P
  size_ratio = t.sys[s].planet[p].radius/t.sys[s].star.radius
  depth = calc_transit_depth(t,s,p)
  duration = calc_transit_duration(t,s,p)
  b = calc_impact_parameter(t.sys[s],p)
  snr_correction = calc_depth_correction_for_grazing_transit(b,size_ratio)
  depth *= snr_correction
  kepid = StellarTable.star_table(t.sys[s].star.id, :kepid)
  osd_duration = get_legal_durations(period,duration)	#tests if durations are included in Kepler's observations for a certain planet period. If not, returns nearest possible duration
  osd = WindowFunction.interp_OSD_from_table(kepid, period, osd_duration)
  if osd_duration > duration				#use a correcting factor if this duration is lower than the minimum searched for this period.
      osd = osd*osd_duration/duration
  end
  ntr = calc_expected_num_transits(t,s,p,sim_param)
  calc_prob_detect_if_transit(t,depth,period,duration,osd, sim_param, num_transit=ntr)
end

function calc_prob_detect_if_transit_with_actual_b_cdpp(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  period = t.sys[s].orbit[p].P
  size_ratio = t.sys[s].planet[p].radius/t.sys[s].star.radius
  depth = calc_transit_depth(t,s,p)
  duration = calc_transit_duration(t,s,p)
  b = calc_impact_parameter(t.sys[s],p)
  snr_correction = calc_depth_correction_for_grazing_transit(b,size_ratio)
  depth *= snr_correction
  cdpp = interpolate_cdpp_to_duration(t, duration)
  ntr = calc_expected_num_transits(t,s,p,sim_param)
  calc_prob_detect_if_transit_cdpp(t,depth,period,duration,cdpp, sim_param, num_transit=ntr)
end

# Compute probability of detection if we average over impact parameters b~U[0,1)
function calc_ave_prob_detect_if_transit_from_snr(t::KeplerTarget, snr_central::Float64, period::Float64, duration_central::Float64, size_ratio::Float64, osd_central::Float64, sim_param::SimParam; num_transit::Float64 = 1)
  const min_pdet_nonzero = 1.0e-4
  wf = kepler_window_function(t, num_transit, period, duration_central)

  detection_efficiency_central = detection_efficiency_model(snr_central, num_transit, min_pdet_nonzero=min_pdet_nonzero)
  if wf*detection_efficiency_central <= min_pdet_nonzero
     return 0.
  end

  # Breaking integral into two sections [0,1-b_boundary) and [1-b_boundary,1], so need at least 5 points to evaluate integral via trapezoid rule
  const num_impact_param_low_b =  7                            # Number of points to evaluate integral over [0,1-b_boundary) via trapezoid rule
  const num_impact_param_high_b = 5 # (size_ratio<=0.05) ? 5 : 11  # Number of points to evaluate integral over [1-b_boudnary,1) via trapezoid rule.  If using 2*size_ratio for bondary for small planets, then keep this odd, so one point lands on 1-size_ratio.
  @assert(num_impact_param_low_b >= 5)
  @assert(num_impact_param_high_b >= 3)
  const num_impact_param = num_impact_param_low_b+num_impact_param_high_b-1 # One point is shared
  const b_boundary = (size_ratio <= 0.15) ? 2*size_ratio : min(max(0.3,size_ratio),0.5)
  b = Array{Float64}(num_impact_param)
  weight = Array{Float64}(num_impact_param)
  b[1:num_impact_param_low_b] = linspace(0.0,1-b_boundary,num_impact_param_low_b)
  b[num_impact_param_low_b:num_impact_param] = linspace(1-b_boundary,1.0,num_impact_param_high_b)
  weight[1:num_impact_param_low_b] = (1-b_boundary)/(num_impact_param_low_b-1)  # Points for first integral
  weight[1] *= 0.5                        # Lower endpoint of first integral
  weight[num_impact_param_low_b] *= 0.5   # Upper endpoint of first integral
  weight[num_impact_param_low_b] += 0.5*(b_boundary)/(num_impact_param_high_b-1) # Also lower endpoint of second integral
  weight[(num_impact_param_low_b+1):num_impact_param] = b_boundary/(num_impact_param_high_b-1)
  weight[num_impact_param] *= 0.5         # Upper endpoint of second integral
  @assert isapprox(sum(weight),1.0)

  function integrand(b::Float64)::Float64
     depth_factor = calc_depth_correction_for_grazing_transit(b,size_ratio)
     duration_factor = calc_transit_duration_factor_for_impact_parameter_b(b,size_ratio)
     kepid = StellarTable.star_table(t.sys[1].star.id, :kepid)
     osd_duration = get_legal_durations(period,duration_central*duration_factor)	#tests if durations are included in Kepler's observations for a certain planet period. If not, returns nearest possible duration
     osd = WindowFunction.interp_OSD_from_table(kepid, period, osd_duration)
     if osd_duration > duration_central*duration_factor				#use a correcting factor if this duration is lower than the minimum searched for this period.
	osd = osd*osd_duration/(duration_central*duration_factor)
     end
     snr_factor = depth_factor*(osd_central/osd)
     detection_efficiency_model(snr_central*snr_factor, num_transit, min_pdet_nonzero=min_pdet_nonzero)
  end

  ave_detection_efficiency = sum(weight .* map(integrand,b)::Vector{Float64} )

  return wf*ave_detection_efficiency
end

# Compute probability of detection if we average over impact parameters b~U[0,1)
function calc_ave_prob_detect_if_transit_from_snr_cdpp(t::KeplerTarget, snr_central::Float64, period::Float64, duration_central::Float64, size_ratio::Float64, cdpp_central::Float64, sim_param::SimParam; num_transit::Float64 = 1)
  const min_pdet_nonzero = 1.0e-4
  wf = kepler_window_function(t, num_transit, period, duration_central)

  detection_efficiency_central = detection_efficiency_model(snr_central, num_transit, min_pdet_nonzero=min_pdet_nonzero)
  if wf*detection_efficiency_central <= min_pdet_nonzero
     return 0.
  end

  # Breaking integral into two sections [0,1-b_boundary) and [1-b_boundary,1], so need at least 5 points to evaluate integral via trapezoid rule
  const num_impact_param_low_b =  7                            # Number of points to evaluate integral over [0,1-b_boundary) via trapezoid rule
  const num_impact_param_high_b = 5 # (size_ratio<=0.05) ? 5 : 11  # Number of points to evaluate integral over [1-b_boudnary,1) via trapezoid rule.  If using 2*size_ratio for bondary for small planets, then keep this odd, so one point lands on 1-size_ratio.
  @assert(num_impact_param_low_b >= 5)
  @assert(num_impact_param_high_b >= 3)
  const num_impact_param = num_impact_param_low_b+num_impact_param_high_b-1 # One point is shared
  const b_boundary = (size_ratio <= 0.15) ? 2*size_ratio : min(max(0.3,size_ratio),0.5)
  b = Array{Float64}(num_impact_param)
  weight = Array{Float64}(num_impact_param)
  b[1:num_impact_param_low_b] = linspace(0.0,1-b_boundary,num_impact_param_low_b)
  b[num_impact_param_low_b:num_impact_param] = linspace(1-b_boundary,1.0,num_impact_param_high_b)
  weight[1:num_impact_param_low_b] = (1-b_boundary)/(num_impact_param_low_b-1)  # Points for first integral
  weight[1] *= 0.5                        # Lower endpoint of first integral
  weight[num_impact_param_low_b] *= 0.5   # Upper endpoint of first integral
  weight[num_impact_param_low_b] += 0.5*(b_boundary)/(num_impact_param_high_b-1) # Also lower endpoint of second integral
  weight[(num_impact_param_low_b+1):num_impact_param] = b_boundary/(num_impact_param_high_b-1)
  weight[num_impact_param] *= 0.5         # Upper endpoint of second integral
  @assert isapprox(sum(weight),1.0)

  function integrand(b::Float64)::Float64
     depth_factor = calc_depth_correction_for_grazing_transit(b,size_ratio)
     duration_factor = calc_transit_duration_factor_for_impact_parameter_b(b,size_ratio)
     cdpp = interpolate_cdpp_to_duration(t,duration_central*duration_factor)
     snr_factor = depth_factor*sqrt(duration_factor)*(cdpp_central/cdpp)
     detection_efficiency_model(snr_central*snr_factor, num_transit, min_pdet_nonzero=min_pdet_nonzero)
  end

  ave_detection_efficiency = sum(weight .* map(integrand,b)::Vector{Float64} )

  return wf*ave_detection_efficiency
end

function calc_ave_prob_detect_if_transit_cdpp(t::KeplerTarget, depth::Float64, period::Float64, duration_central::Float64, size_ratio::Float64, sim_param::SimParam; num_transit::Float64 = 1)
  cdpp_central = interpolate_cdpp_to_duration(t, duration_central)
  snr_central = calc_snr_if_transit(t,depth,duration_central,cdpp_central, sim_param, num_transit=num_transit)
  return calc_ave_prob_detect_if_transit_from_snr_cdpp(t, snr_central, period, duration_central, size_ratio, cdpp_central, sim_param, num_transit=num_transit)
end

function calc_ave_prob_detect_if_transit_cdpp(t::KeplerTarget, s::Integer, p::Integer, sim_param::SimParam)
  size_ratio = t.sys[s].planet[p].radius/t.sys[s].star.radius
  depth = calc_transit_depth(t,s,p)
  period = t.sys[s].orbit[p].P
  duration_central = calc_transit_duration_central(t,s,p)
  ntr = calc_expected_num_transits(t,s,p,sim_param)
  calc_ave_prob_detect_if_transit_cdpp(t,depth,period,duration_central, size_ratio, sim_param, num_transit=ntr)
end


text:
"""
    Triangulate the polytope faces.

`βs` is an npoints-by-(2*dim + 2) dimensional array of convex expansion coefficients of the
polytope generators (vertices) in terms of the original simplices s₁ and s₂. The first
(dim + 1) columns correspond to s₁, while the remaining (dim + 1) columns correspond to s₂.

`n` is the the dimension of the space, and `dim` number of generators of the polytope.
"""
function TriangulationPolytopeFaces(βs::AbstractArray{Float64, 2}, n_generators::Int, dim::Int)
    # We adopt the notations
    # NSF ≡ Nonsimplicial faces
    # APF ≡ All polytope faces.
    triang_NSF = zeros(Int, 0, dim + 1)
    triang_APF = zeros(Int, 0, dim + 1)

    faces_containing_intvers = βs[βs .== 0]

    # Find the faces containing intersecting vertices.
    faces_containing_intvers = zeros(Int, size(βs))
    faces_containing_intvers[βs .== 0] .= 1 # set zero convex parameters to 1
    faces_containing_intvers[βs .> 0] .= 0 # set positive convex parameters to 0

    # What is the number of intersecting points in each of the faces?
    n_intpoints_ineachface = ones(Int, 1, n_generators) * faces_containing_intvers
    n_intpoints_ineachface = sum(faces_containing_intvers, dims=1)

    # Indices of the potential faces of the polytope. Different indices might correspond
    # to the same face.
    faceindices = transpose(heaviside0(n_intpoints_ineachface .- dim)) .* (1:2*dim+2)
    faceindices = (LinearIndices(faceindices))[findall(x->x!=0, faceindices)]

    # The number of faces of the intersecting volume polytope
    numofPolFaces = length(faceindices)

    # The number vertices furnishing each polytope face.
    n_intpoints_ineachface = round.(Int64, n_intpoints_ineachface[faceindices])

    # The indices of the intersecting points (in IntVert) furnishing each polytope face.
    # Array of dimension dim x numofPolFaces. Each column represents a potential polytope face.
    npts_eachface = faces_containing_intvers[:, faceindices] .*
                            repeat(collect(1:n_generators), 1, numofPolFaces)

    # Go through each faces and decide whether it is a true face or a boundary of a face.
    # Disregard stuff if it is not a true face. 'NSF_inds' starts out with only
    # zeros. If the potential polytope face is a true face, then set the value for that
    # face to 1.
    NSF_inds = classify_faces(numofPolFaces, βs, dim, npts_eachface, faceindices)
    n_NSF = length(NSF_inds) # The number of nonsingular faces

    if n_NSF >= dim + 1
        # all these arrays contain the corresponding information but only for the
        # non singular faces. And therefore, the second dimension goes over 1:n_NSF
        NonSingularPolytopeFaces = faceindices[NSF_inds]
        NumOfVertInEachFace = n_intpoints_ineachface[NSF_inds]
        NonSingularPointsInEachFace = npts_eachface[1:end, vec(NSF_inds)]
        # THe number of polytope faces that are actually simplices
        Simplicial = vec(findall(x->x!=0, heaviside0(dim .- NumOfVertInEachFace) .* collect(1:n_NSF)))
        NonSimplicialFaces = 0

        if length(Simplicial) == 0
            #println("No polytope faces are simplices")
            NonSimplicialFaces = NonSingularPolytopeFaces
            NonSimplicial = 1:n_NSF
            VerticesSFaces = 0
        elseif length(Simplicial) == n_NSF
            #println("All polytope faces are simplices")

            VerticesSFaces = NonSingularPointsInEachFace
            inner = reshape(VerticesSFaces, size(Simplicial, 1) * n_generators, 1)
            inner_nonzeros = inner[findall(x->x!=0, inner)]
            inner_reshaped_transposed = copy(transpose(reshape(inner_nonzeros, dim, size(Simplicial, 1))))
            VerticesSFaces = inner_reshaped_transposed
        else
            #println("Some polytope faces are simplices")

            VerticesSFaces = NonSingularPointsInEachFace[:, Simplicial]
            inner = reshape(VerticesSFaces, size(Simplicial, 1)*n_generators, 1)

            inner_nonzeros = inner[(LinearIndices(inner))[findall(x->x!=0, inner)]]
            inner_reshaped_transposed = transpose(reshape(inner_nonzeros, dim, size(Simplicial, 1)))
            VerticesSFaces = inner_reshaped_transposed
            NonSimplicial = complementary(Simplicial, n_NSF)
            NonSimplicialFaces = NonSingularPolytopeFaces[NonSimplicial]
        end

        if NonSimplicialFaces[1] > 0
            #println("One or more nonsimplical faces.")
            VerticesNSFaces = NonSingularPointsInEachFace[:, NonSimplicial]
            SimplexIndexNS = round.(Int64, ceil.(NonSimplicialFaces / (dim + 1)))
            SimplexFaceIndexNS = round.(Int64, NonSimplicialFaces .- (SimplexIndexNS .- 1) * (dim + 1))
            #t1 = time_ns()
            triang_NSF = TriangulationNonSimplicialFaces(VerticesNSFaces, SimplexIndexNS, SimplexFaceIndexNS, βs, dim)
            #t2 = time_ns()
            #println("Triangulation of nonsimplicial faces took ", (t2 - t1)/10^6, " ms")
        end
        if !isempty(triang_NSF)
            triang_APF = Update(VerticesSFaces, triang_NSF)
        end
    else
        #println("# Non-simplicial faces < n + 1\n")
    end
    return triang_APF, n_NSF
end

"""
Go through each faces and decide whether it is a true face or a boundary of a face.
Returns a vector of length equal to the number of faces in the intersecting polytope.
The ith entry of this vector is _i_ if the corresponding face is a true face, 0 otherwise.
"""
function classify_faces(numofPolFaces::Int,
                        βs::AbstractArray{Float64, 2},
                        dim::Int,
                        npts_eachface::AbstractArray{Int, 2},
                        faceindices::AbstractArray{Int, 1})

    NSF_inds = zeros(Int, 0)

    # Loop over faces.
    for faceindex = 1:numofPolFaces
        inds = findall(x->x!=0, npts_eachface[:, faceindex])
        Aux = heaviside0(-sum(βs[inds, :], dims=1))
        multiplicity = [Aux * ones(2*dim + 2, 1);
                        Aux * [ones(dim + 1, 1);
                        2 * ones(dim + 1, 1)]]
        #multiplicity(1): number of times that the face with index faceindices(a) appears
        #multiplicity(2): (number of faces of simplex1 containing the face) +
        #                   2*(number of faces of simplex2 containing the face)
        if multiplicity[1] == 1 || (multiplicity[1] == 2 && multiplicity[2] == 3 &&
                                    faceindices[faceindex] <= dim + 1)
            push!(NSF_inds, faceindex)
        end
    end

    return NSF_inds
end


text:
"""
Clamped cylinder, transient vibration.

The cylinder is banged at the initial time (i.e. it is given an initial velocity).
"""
module clamp_cyl_forc_expl_examples

using LinearAlgebra
using SparseArrays
using Arpack
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: default_layout_3d, plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
using PlotlyJS
using Gnuplot
using FinEtools.MeshExportModule.VTKWrite: vtkwrite

function loop!(M, K, U0, V0, tend, dt, force!, peek)
    U1 = deepcopy(U0)
    V1 = deepcopy(U0)
    A0 = deepcopy(U0)
    A1 = deepcopy(U0)
    F = deepcopy(U0)
    E = deepcopy(U0)
    invM = deepcopy(U0)
    invM .= 1.0 ./ vec(diag(M))
    nsteps = Int64(round(tend/dt))
    if nsteps*dt < tend
        dt = tend / (nsteps+1)
    end
    t = 0.0
    A0 .= M\force!(F, t);
    peek(0, U0, t)
    @time for step in 1:nsteps
        # External loading
        F = force!(F, t);
        # Displacement update
        @. U1 = U0 + dt*V0 + ((dt^2)/2)*A0;
        #
        F .-=  mul!(E, K, U1)
        # Compute the new acceleration.
        A1 .= invM .* F
        # # Update the velocity
        @. V1 = V0 + (dt/2)* (A0+A1);
        t = t + dt
        U0, U1 = U1, U0
        V0, V1 = V1, V0
        A0, A1 = A1, A0
        peek(step, U0, t)
    end
end

E = 200e9;
nu = 0.3;
rho = 7800.0
R = 0.1
L = 0.8
omegaf = 2*pi*10^4
tend = 15 * 2*pi/omegaf

cylindrical!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) = begin
    r = vec(XYZ); r[2] = 0.0;
    csmatout[:, 3] .= vec(r)/norm(vec(r))
    csmatout[:, 2] .= (0.0, 1.0, 0.0) #  this is along the axis
    cross3!(view(csmatout, :, 1), view(csmatout, :, 2), view(csmatout, :, 3))
    return csmatout
end

function _execute(n = 8, thickness = 0.001, visualize = true)

    @info "Mesh: $n elements per side"
    # Mesh
    tolerance = min(R, L) / n  / 100

    tolerance = R/n/1000
    fens, fes = T3block(360.0, L, n, 4*n);
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        a=fens.xyz[i, 1]; y=fens.xyz[i, 2];
        fens.xyz[i, :] .= (R*sin(a/180*pi), y-L/2, R*cos(a/180*pi))
    end
    fens, fes = mergenodes(fens, fes, tolerance)
    bfes = meshboundary(fes)

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)
    ocsys = CSys(3, 3, cylindrical!)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = FEMMShellT3FFModule.make(IntegDomain(fes, TriRule(1), thickness), ocsys, mater)
    # Set up
    femm.drilling_stiffness_scale = 1.0

# Construct the requisite fields, geometry and displacement
# Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

# Apply EBC's
    l1 = connectednodes(meshboundary(fes))
    for i in 1:6
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

# Assemble the system matrix
    FEMMShellT3FFModule.associategeometry!(femm, geom0)
    K = FEMMShellT3FFModule.stiffness(femm, geom0, u0, Rfield0, dchi);
    M = FEMMShellT3FFModule.mass(femm, geom0, dchi);
    # Check that the mass matrix is diagonal
    # Make sure the matrix is truly diagonal: delete all tiny off-diagonals
    I, J, V = findnz(M)
    for i in 1:length(I)
        if I[i] != J[i]
            V[i] = 0.0
        end
    end
    M = sparse(I, J, V, dchi.nfreedofs, dchi.nfreedofs)

# Solve
    function pwr(A, B)
        fa = cholesky(B)
        v = rand(size(B, 1))
        w = fill(0.0, size(B, 1))
        for i in 1:30
            w = (A * v)
            w = w / norm(w)
            v .= fa \ w
            v = v / norm(v)
        end
        sqrt((v' * A * v) / (v' * B * v))
    end
    @time omega_max = pwr(K, M)
    @show omega_max

    # @time evals, evecs, nconv = eigs(Symmetric(K), Symmetric(M); nev=1, which=:LM, explicittransform=:none)
    # @show sqrt(evals)
    # @show omega_max = sqrt(evals[1])

    @show dt = Float64(0.9* 2/omega_max)

    U0 = gathersysvec(dchi)
    V0 = deepcopy(U0)

    qpoint = selectnode(fens; nearestto=[0 -L/4 R])[1]
    cpoint = selectnode(fens; nearestto=[0 0 -R])[1]
    applyebc!(dchi)
    numberdofs!(dchi);
    qpointdof = dchi.dofnums[qpoint, 3]
    cpointdof = dchi.dofnums[cpoint, 3]
    cpointdof6 = dchi.dofnums[cpoint, 6]

    function computetrac!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        dx = XYZ[1] - fens.xyz[qpoint, 1]
        dy = XYZ[2] - fens.xyz[qpoint, 2]
        dz = XYZ[3] - fens.xyz[qpoint, 3]
        mag = 10^5*exp(-20/R*sqrt(dx^2+dy^2+dz^2))
        r = vec(XYZ); r[2] = 0.0
        r .= vec(r)/norm(vec(r))
        theta = atan(r[3], r[1])
        forceout[1:3] .= r*mag
        forceout[4:6] .= 0.0
        return forceout
    end

    lfemm = FEMMBase(IntegDomain(fes, TriRule(3)))
    fi = ForceIntensity(FFlt, 6, computetrac!);
    Fmag = distribloads(lfemm, geom0, dchi, fi, 2);

    function force!(F, t)
        F .= sin(omegaf*t) .* Fmag
        # @show norm(F)
        return F
    end

    nsteps = Int(round(tend/dt))
    cdeflections = fill(0.0, nsteps+1)
    displacements = []
    nbtw = Int(round(nsteps/100))

    peek(step, U, t) = begin
        cdeflections[step+1] = U[cpointdof]
        if rem(step+1, nbtw) == 0
            push!(displacements, deepcopy(U))
        end
        nothing
    end
    loop!(M, K, U0, V0, nsteps*dt, dt, force!, peek)

    # @gp  "set terminal windows 0 "  :-

    color = "red"
    @gp  :- collect(0.0:dt:(nsteps*dt)) cdeflections " lw 2 lc rgb '$color' with lines title 'Deflection at the center' "  :-

    @gp  :- "set xlabel 'Time'" :-
    @gp  :- "set ylabel 'Deflection'" :-
    @gp  :- "set title 'Free-floating plate'"

# Visualization
    @show length(displacements)
    box = boundingbox(fens.xyz)
    box[1] -= 0.05*L
    box[2] += 0.05*L
    box[3] -= 0.05*L
    box[4] += 0.05*L
    box[5] -= 0.5*L
    box[6] += 0.5*L

    tbox = plot_space_box(reshape(box, 2, 3))
    tenv0 = plot_midsurface(fens, fes; x=geom0.values, u=0.0 .* dchi.values[:, 1:3], R=Rfield0.values, facecolor="rgb(155, 155, 155)", opacity=0.3);

    plots = cat(tbox, tenv0; dims=1)
    layout = default_layout_3d(;width=600, height=600)
    layout[:scene][:aspectmode] = "data"
    pl = render(plots; layout=layout, title = "Step ")
    sleep(2.5)

    scale = 10^6;
    for i in 1:length(displacements)
        scattersysvec!(dchi, scale .* displacements[i])
        Rfield1 = deepcopy(Rfield0)
        update_rotation_field!(Rfield1, dchi)
        tenv1 = plot_midsurface(fens, fes; x=geom0.values, u=dchi.values[:, 1:3], R=Rfield1.values, facecolor="rgb(50, 125, 125)",  lighting = attr(diffuse=0.3, fresnel=4, specular=0.5), lightposition = attr(x = 400, y = -100, z = -5));
        plots = cat(tbox, tenv0, tenv1; dims=1)
        pl.plot.layout[:title] = "Step $(i)"
        react!(pl, plots, pl.plot.layout)
        sleep(0.115)
    end

end

function test_convergence()
    @info "Clamped cylinder vibration"
    for n in [64]
        _execute(n, 0.01, !false)
    end
    return true
end

function allrun()
    println("#####################################################")
    println("# test_convergence ")
    test_convergence()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
"""
# module Matdyn

# Examples

```jldoctest
julia>
```
"""
module Matdyn

end


text:
#@testset "recompilation" begin
#    @eval f() = 1
#
#    @test f() == 1
#    @test TinyCassette.execute(nothing, f) == 1
#
#    @eval f() = 2
#    @test f() == 2
#    @test TinyCassette.execute(nothing, f) == 2
#end


text:
function computeRelativityCorrection(x::Float64,y::Float64,z::Float64,
    vx::Float64,vy::Float64,vz::Float64,)
    relcorr = -2.0 * (x * vx + y * vy + z * vz) / C_MPS^2
    return relcorr
end

function getSidTime(dt::DateTime)
    jd = datetime2julian(DateTime(year(dt),
                            month(dt), day(dt)))
    jc = (jd - 2451545.0) / 36525.0
    sid = 24110.54841 + jc * (8640184.812866 +
                            jc * (0.093104 -
                                jc * 0.0000062 ) )
    sid /= 3600.0
    sid = sid % 24.0
    if sid < 0.0 sid += 24.0 end
    return sid
end

function secondofday(dt::DateTime)
    return hour(dt) * 3600.0 + minute(dt) * 60.0 +
                second(dt) + millisecond(dt) / 100000.0
end

function derivative(inState::Array{Float64,1},
                    accel::Array{Float64,1},
                    j20::Float64=PZ90Ellipsoid.j20,
                    mu::Float64=PZ90Ellipsoid.gm_km,
                    ae::Float64=PZ90Ellipsoid.a_km,)
    x  = inState[1]
    y  = inState[3]
    z  = inState[5]
    r² = x^2 + y^2 + z^2
    r  = sqrt(r²)
    xmu = mu / r²
    rho = ae / r
    xr  = x / r
    yr  = y / r
    zr  = z / r
    zr2 = zr^2
    k1  = j20 * xmu * 1.5 * rho^2
    cm  = k1 * (1.0 - 5.0 * zr2)
    cmz = k1 * (3.0 - 5.0 * zr2)
    k2  = cm - xmu

    gloAx = k2 * xr + accel[1]
    gloAy = k2 * yr + accel[2]
    gloAz = (cmz - xmu) * zr + accel[3]

    dxt = zeros(6)

    dxt[1] = inState[2]
    dxt[2] = gloAx
    dxt[3] = inState[4]
    dxt[4] = gloAy
    dxt[5] = inState[6]
    dxt[6] = gloAz
    dxt
end

function calculate_Fit_pos(epoch::DateTime,
        ephTime::DateTime,
        clkbias::Float64, # tau
        clkdrift::Float64, # gamma
        ix::Float64,iy::Float64,iz::Float64,
        ivx::Float64,ivy::Float64,ivz::Float64,
        iax::Float64,iay::Float64,iaz::Float64,
        health::Float64=-1.0,
        step::Float64=1.0,
        we::Float64=PZ90Ellipsoid.ge.angVelocity,
)
    elapte = (epoch - ephTime).value / 1000.0 # TODO  second?
    # Values to be returned will be stored here
    # If the exact epoch is found, let's return the values
    if epoch == ephTime        # exact match for epoch
        oxx = ix * 1.e3;   # m
        oxy = iy * 1.e3;   # m
        oxz = iz * 1.e3;   # m
        ovx = ivx * 1.e3;  # m/sec
        ovy = ivy * 1.e3;  # m/sec
        ovz = ivz * 1.e3;  # m/sec
    # In the GLONASS system, 'clkbias' already includes the
    # relativistic correction, therefore we must substract the late
    # from the former.
        relcorr = computeRelativityCorrection(oxx, oxy, oxz, ovx, ovy, ovz)
        clkbias = clkbias + clkdrift * elapte - relcorr;
        clkdrift = clkdrift;
        frame = fmPZ90;

    # We are done, let's return
        return Xvt(ECEF(oxx, oxy, oxz), (ovx, ovy, ovz),
    relcorr,clkbias,clkdrift,frame,
    health == 1.0 ? hsHealthy : hsUnhealthy)
    end

    # Get the data out of the GloRecord structure
    px = ix   # X coordinate (km)
    vx = ivx   # X velocity   (km/s)
    ax = iax   # X acceleration (km/s^2)

    py = iy   # Y coordinate
    vy = ivy   # Y velocity
    ay = iay   # Y acceleration

    pz = iz   # Z coordinate
    vz = ivz   # Z velocity
    az = iaz   # Z acceleration

    # We will need some PZ-90 ellipsoid parameters

    # Get sidereal time at Greenwich at 0 hours UT

    gst = getSidTime(ephTime)
    s0 =  gst * π / 12.0
    numSeconds = secondofday(ephTime)

    s =  s0 + we * numSeconds
    cs =  cos(s)
    ss =  sin(s)

    # Initial state matrix
    initialState = zeros(6)
    accel = zeros(3)
    # dxt1  = zeros(6)
    # dxt2  = zeros(6)
    # dxt3  = zeros(6)
    # dxt4  = zeros(6)
    tempRes = zeros(6)

    # Get the reference state out of GloEphemeris object data. Values
    # must be rotated from PZ-90 to an absolute coordinate system
    # Initial x coordinate (m)
    initialState[1]  = (px * cs - py * ss);
    # Initial y coordinate
    initialState[3]  = (px * ss + py * cs);
    # Initial z coordinate
    initialState[5]  = pz;

    # Initial x velocity   (m/s)
    initialState[2]  = (vx * cs - vy * ss - we * initialState[3] );
    # Initial y velocity
    initialState[4]  = (vx * ss + vy * cs + we * initialState[1] );
    # Initial z velocity
    initialState[6]  = vz;

    # Integrate satellite state to desired epoch using the given step
    rkStep = step

    if epoch < ephTime
        rkStep = step * (-1.0)
    end

    rkStepNs = Nanosecond(rkStep * 1e9)
    workEpoch = TimeDate(ephTime)
    tolerance = Nanosecond(1) # 1e-9

    while true
    # If we are about to overstep, change the stepsize appropriately
    # to hit our target final time.
        if rkStep > 0.0
            if (workEpoch + rkStepNs) > epoch
                rkStepNs = epoch - workEpoch
                rkStep = rkStepNs.value / 1e9
            end
        else
            if (workEpoch + rkStepNs) < epoch
                rkStepNs = (epoch - workEpoch)
                rkStep = rkStepNs.value / 1e9
            end
        end

        numSeconds += rkStep;
        s = s0 + we * ( numSeconds );
        cs = cos(s);
        ss = sin(s);

        # Accelerations are computed once per iteration
        accel[1] = ax * cs - ay * ss
        accel[2] = ax * ss + ay * cs
        accel[3] = az
        dxt1 = derivative(initialState, accel)
        for j = 1:6
            tempRes[j] = initialState[j] + rkStep * dxt1[j] / 2.0
        end
        dxt2 = derivative(tempRes, accel)
        for j = 1:6
            tempRes[j] = initialState[j] + rkStep * dxt2[j] / 2.0
        end
        dxt3 = derivative(tempRes, accel)
        for j = 1:6
            tempRes[j] = initialState[j] + rkStep * dxt3[j];
        end
        dxt4 = derivative(tempRes, accel)
        for j = 1:6
            initialState[j] = initialState[j] + rkStep *
                            ( dxt1[j] +
                                2.0 * ( dxt2[j] + dxt3[j] ) +
                                dxt4[j] ) / 6.0;
        end
        # @bp
        # If we are within tolerance of the target time, we are done.
        workEpoch += rkStepNs
        if ( abs(epoch - workEpoch) < tolerance )
            break
        end
    end
    # End of 'while (!done)...'
    px = initialState[1]
    py = initialState[3]
    pz = initialState[5]
    vx = initialState[2]
    vy = initialState[4]
    vz = initialState[6]
    # @bp

    oxx = 1000.0 * ( px * cs + py * ss );         # X coordinate
    oxy = 1000.0 * (-px * ss + py * cs);          # Y coordinate
    oxz = 1000.0 * pz;                        # Z coordinate

    ovx = 1000.0 * ( vx * cs + vy * ss + we * (oxy / 1000.0) ) # X velocity
    ovy = 1000.0 * (-vx * ss + vy * cs - we * (oxx / 1000.0) ) # Y velocity
    ovz = 1000.0 * vz;                        # Z velocity

    # In the GLONASS system, 'clkbias' already includes the relativistic
    # correction, therefore we must substract the late from the former.
    ox = ECEF(oxx, oxy, oxz)

    relcorr = computeRelativityCorrection(oxx, oxy, oxz, ovx, ovy, ovz)
    clkbias = clkbias + clkdrift * elapte - relcorr
    clkdrift = clkdrift;
    frame = fmPZ90;

    # We are done, let's return
    return Xvt(ECEF(oxx, oxy, oxz), (ovx, ovy, ovz),
                    relcorr,clkbias,clkdrift,frame,
                    health == 1.0 ? hsHealthy : hsUnhealthy)

end

function calculate_Fit_pos(ees::EphEpochStore, dt::DateTime,
                            vaild_d::Minute=Minute(15),)
    k = ees.keys
    if length(k) == 0
        return ECEF{Float64}(NaN, NaN, NaN)
    end
    d = k .- dt
    i = argmin(abs.(d))

    md = d[i]
    if md > vaild_d || md < -vaild_d
        return ECEF{Float64}(NaN, NaN, NaN)
    end

    mk = k[i]
    me = ees.eph[mk].d

    xvt = calculate_Fit_pos(dt, mk,
                        me[1], me[2],

                        me[4], me[8],  me[12],
                        me[5], me[9],  me[13],
                        me[6], me[10], me[14],

                        me[7],)
    return xvt.x
end


text:
using JuMP
using Cbc

LP=Model(Cbc.Optimizer)
@variable(LP, x11>=0)
@variable(LP, x12>=0)
@variable(LP, x21>=0)
@variable(LP, x22>=0)
@variable(LP, x31>=0)
@variable(LP, x32>=0)
@variable(LP, z>=0)
@constraint(LP, x11+x21+x31<=z)
@constraint(LP, x12+x22+x32<=24)
@constraint(LP, x11+2*x12<=15)
@constraint(LP, 1.2*x21+3x22<=24)
@constraint(LP, 2*x31+3x32<=18)
@objective(LP, Max, x11+2*x12+1.2*x21+3*x22+2*x31+3*x32)
status= JuMP.optimize!(LP)

LP=Model(Cbc.Optimizer)
@variable(LP, x11>=0)
@variable(LP, x12>=0)
@variable(LP, x21>=0)
@variable(LP, x22>=0)
@variable(LP, x31>=0)
@variable(LP, x32>=0)
@variable(LP, z>=0)
@constraint(LP, x11+x21+x31<=z)
@constraint(LP, x12+x22+x32<=z)
@constraint(LP, x11+2*x12>=15)
@constraint(LP, 1.2*x21+3*x22>=24)
@constraint(LP, 2*x31+3*x32>=18)
@objective(LP, Min, z)
status= JuMP.optimize!(LP)

LP=Model(Cbc.Optimizer)
@variable(LP, x>=0)
@variable(LP, y>=0)
@constraint(LP, 4x+6y>=30)
@constraint(LP, 6x+2y>=20)
@constraint(LP, x+2y>=12)
@objective(LP, Min, 2x+3y)
status= JuMP.optimize!(LP)
println("X: $(JuMP.value(x))")
println("Y : $(JuMP.value(y))")

max_stars_repo_name: valenmgama/optimization-tec

max_issues_repo_name: valenmgama/optimization-tec

max_forks_repo_name: valenmgama/optimization-tec


text:
n = 79 # rand(40:100)
A = matrixdepot("gilbert", n)
@test issymmetric(A)

B = matrixdepot("gilbert", n, 0.2)
C = matrixdepot("gilbert", Int, n, 0.5)
@test nnz(B) <= nnz(C)

@test matrixdepot("gilbert", Int, 1) !== nothing
println("'gilbert' passed test ...")


text:
module RethinkDB

import JSON

include("query.jl")

type RethinkDBConnection
  socket :: Base.TCPSocket
end

# TODO: handle error is not connected or incorrect handshake
function connect(server::AbstractString = "localhost", port::Int = 28015)
  c = RethinkDBConnection(Base.connect(server, port))
  handshake(c)
  c
end

function handshake(conn::RethinkDBConnection)
  # Version.V0_4
  version = UInt32(0x400c2d20)

  # Key Size
  key_size = UInt32(0)

  # Protocol.JSON
  protocol = UInt32(0x7e6970c7)

  handshake = pack_command([version, key_size, protocol])
  write(conn.socket, handshake)
  is_valid_handshake(conn)
end

function is_valid_handshake(conn::RethinkDBConnection)
  readstring(conn.socket) == "SUCCESS"
end

function readstring(sock::TCPSocket, msg = "")
  c = read(sock, UInt8)
  s = convert(Char, c)
  msg = string(msg, s)
  if (s == '\0')
    return chop(msg)
  else
    readstring(sock, msg)
  end
end

function pack_command(args...)
  o = Base.IOBuffer()
  for enc_val in args
    write(o, enc_val)
  end
  o.data
end

function disconnect(conn::RethinkDBConnection)
  close(conn.socket)
end

function exec(conn::RethinkDBConnection, q)
  j = JSON.json([1 ; Array[q.value]])
  send_command(conn, j)
end

function run(conn::RethinkDBConnection, q)
  o = exec(conn, q)
  response_type = o["t"]

  if response_type == 1
    return o["r"][1]
  elseif response_type == 2
    return o["r"]
  end
  o["r"]
end

function token()
  t = Array{UInt64}(1)
  t[1] = object_id(t)
  return t[1]
end

function send_command(conn::RethinkDBConnection, json)
  t = token()
  q = pack_command([ t, convert(UInt32, length(json)), json ])

  write(conn.socket, q)
  read_response(conn, t)
end

function read_response(conn::RethinkDBConnection, token)
  remote_token = read(conn.socket, UInt64)
  if remote_token != token
    return "Error"
  end

  len = read(conn.socket, UInt32)
  res = read(conn.socket, len)

  output = convert(UTF8String, res)
  JSON.parse(output)
end

include("sample.jl")

end

max_stars_repo_name: mrcsparker/julia-rethinkdb

max_issues_repo_name: mrcsparker/julia-rethinkdb

max_forks_repo_name: mrcsparker/julia-rethinkdb


text:
export MidiPerformanceEncoder, TextMelodyEncoder
export encode_notesequence, decode_to_notesequence

using NoteSequences
using NoteSequences.PerformanceRepr, NoteSequences.MelodyRepr
import NoteSequences: encode_event, decode_event

const PAD_TOKEN = 1
const EOS_TOKEN = 2

"""
    MidiPerformanceEncoder <: Any

Encoder to convert between performance event indices and `NoteSequences`.

`MidiPerformanceEncoder` accounts for padding and EOS tokens, while
`PerformanceOneHotEncoding` does not.

## Fields
* `encoding::PerformanceOneHotEncoding`: Performance encoder object that is used for encoding.
* `steps_per_second::Int`: Amount of quantization steps per second.
* `num_velocitybins::Int`: Number of velocity bins.
* `add_eos::Bool`: If true, add the EOS token at the end while encoding.
"""
struct MidiPerformanceEncoder
    encoding::PerformanceOneHotEncoding
    steps_per_second::Int
    num_velocitybins::Int
    add_eos::Bool

    function MidiPerformanceEncoder(steps_per_second::Int,
                                    num_velocitybins::Int,
                                    minpitch::Int,
                                    maxpitch::Int,
                                    add_eos::Bool)

        encoding = PerformanceOneHotEncoding(num_velocitybins=num_velocitybins,
                                             max_shift_steps=steps_per_second,
                                             minpitch=minpitch,
                                             maxpitch=maxpitch)

        new(encoding, steps_per_second, num_velocitybins, add_eos)
    end
end

function Base.getproperty(encoder::MidiPerformanceEncoder, sym::Symbol)
    if sym === :num_reserved_tokens
        return 2 # PAD, EOS
    elseif sym === :vocab_size
        return encoder.encoding.num_classes + encoder.num_reserved_tokens
    elseif sym === :default_event
        return encoder.encoding.defaultevent
    else
        getfield(encoder, sym)
    end
end

"""
    encode_notesequence(ns::NoteSequence, encoder::MidiPerformanceEncoder)

Convert the given `NoteSequence` into performance event indices.
If `encoder.add_eos` is true, `EOS_TOKEN` will be added at the end.
"""
function encode_notesequence(ns::NoteSequence, encoder::MidiPerformanceEncoder)
    NoteSequences.absolutequantize!(ns, encoder.steps_per_second)
    performance = Performance(ns, velocity_bins=encoder.num_velocitybins)
    indices = [encode_event(event, encoder) for event in performance]
    encoder.add_eos && push!(indices, EOS_TOKEN)

    indices
end

"""
    decode_to_notesequence(indices::Vector{Int}, encoder::MidiPerformanceEncoder)

Convert the given performance event indices into a `NoteSequence`.
"""
function decode_to_notesequence(indices::Vector{Int}, encoder::MidiPerformanceEncoder)
    performance = Performance(encoder.steps_per_second, velocity_bins=encoder.num_velocitybins)
    performance_events = [decode_event(index, encoder) for index in indices]
    append!(performance, performance_events)

    getnotesequence(performance)
end

"""
    encode_event(event::PerformanceEvent, encoder::MidiPerformanceEncoder)

Encode a `PerformanceEvent` into it's one-hot index using a `MidiPerformanceEncoder`.
"""
function encode_event(event::PerformanceEvent, encoder::MidiPerformanceEncoder)
    index = encode_event(event, encoder.encoding)
    index + encoder.num_reserved_tokens
end

"""
    decode_event(index::Int, encoder::MidiPerformanceEncoder)

Decode a one-hot index into it's `PerformanceEvent` using a `MidiPerformanceEncoder`.
"""
function decode_event(index::Int, encoder::MidiPerformanceEncoder)
    index -= encoder.num_reserved_tokens
    decode_event(index, encoder.encoding)
end

"""
    TextMelodyEncoder <: Any

Encoder to convert `Melody` events and `NoteSequences`
to `Melody` one-hot indices.

`TextMelodyEncoder` accounts for padding and EOS tokens, while
`MelodyOneHotEncoding` does not.

## Fields
* `encoding::MelodyOneHotEncoding`: The `Melody` encoder object used for encoding.
* `steps_per_second::Int`: Quantization steps per second.
"""
struct TextMelodyEncoder
    encoding::MelodyOneHotEncoding
    steps_per_second::Int

    function TextMelodyEncoder(minpitch::Int, maxpitch::Int, steps_per_second::Int)
        encoding = MelodyOneHotEncoding(minpitch, maxpitch)
        new(encoding, steps_per_second)
    end
end

function Base.getproperty(encoder::TextMelodyEncoder, sym::Symbol)
    if sym === :num_reserved_tokens
        return 2 # PAD, EOS
    elseif sym === :vocab_size
        return encoder.encoding.num_classes + encoder.num_reserved_tokens
    elseif sym === :default_event
        return encoder.encoding.defaultevent
    else
        getfield(encoder, sym)
    end
end

"""
    encode_notesequence(ns::NoteSequence, encoder::TextMelodyEncoder)

Convert the given `NoteSequence` into melody event indices.
"""
function encode_notesequence(ns::NoteSequence, encoder::TextMelodyEncoder)
    !ns.isquantized && NoteSequences.absolutequantize!(ns, encoder.steps_per_second)
    melody = fill(NoteSequences.MelodyRepr.MELODY_NO_EVENT, ns.total_time)
    sort!(ns.notes, by=note->note.start_time)

    for note in ns.notes
        melody[note.start_time + 1] = note.pitch
        if note.end_time < ns.total_time
            melody[note.end_time] = NoteSequences.MelodyRepr.MELODY_NOTE_OFF
        end
    end

    [encode_event(event, encoder) for event in melody]
end

"""
    encode_event(event::Int, encoder::TextMelodyEncoder)

Encode a `Melody` event into it's one-hot index using a `TextMelodyEncoder`.
"""
function encode_event(event::Int, encoder::TextMelodyEncoder)
    index = encode_event(event, encoder.encoding)
    index + encoder.num_reserved_tokens
end


text:
using StateSpaceModels, Distributions, LinearAlgebra, Plots, Random

# Fix seed
Random.seed!(1)

# Number of observations
n = 400
# State dimension (2d position + 2d speed)
m = 4
# Measurements dimension (2d noisy position)
p = 2
# Control dimension (2d acceleration)
q = 2

# Damping ratio
ρ = 0.05
# Time delta
Δ = 1.0

# State transition matrix
T = kron(Matrix{Float64}(I, p, p), [1.0 (1.0 - ρ * Δ / 2.0) * Δ; 0.0 (1.0 - ρ * Δ)])
# Input matrix
R = kron(Matrix{Float64}(I, p, p), [0.5 * Δ^2; Δ])
# Output (measurement) matrix
Z = kron(Matrix{Float64}(I, p, p), [1.0 0.0])

# Generate random actuators
Q = 0.5 * Matrix{Float64}(I, q, q)
η = MvNormal(Q)

# Generate random measurement noise
H = 2.0 * Matrix{Float64}(I, p, p)
ε = MvNormal(H)

# Simulate vehicle trajectory
α = zeros(n + 1, m)
y = zeros(n, p)
for t in 1:n
    y[t, :] = Z * α[t, :] + rand(ε)
    α[t + 1, :] = T * α[t, :] + R * rand(η)
end
α = α[1:n, :]

# User defined model
model = StateSpaceModel(y, Z, T, R)
# Estimate vehicle speed and position
ss = statespace(model)

ss.model.H
ss.model.Q

anim = @animate for i in 1:n
    plot(y[1:i, 1], y[1:i, 2], label="Measured position", line=:scatter, lw=2, markeralpha=0.2, color=:black, title="Vehicle tracking")
    plot!(α[1:i, 1], α[1:i, 3], label="True position", lw=3, color=:indianred)
    plot!(ss.filter.a[2:i+1, 1], ss.filter.a[2:i+1, 3], label="Estimated position", lw=2, color=:forestgreen)
end

gif(anim, "vehicle_tracking.gif", fps = 15)