text:
mutable struct Trace{Tspl <: AbstractSampler, Tvi <: AbstractVarInfo, Tmodel <: Model}
    task  ::  Task
    vi    ::  Tvi
    spl   ::  Tspl
    model ::  Tmodel
    Trace{Tspl, Tvi, Tmodel}() where {Tspl, Tvi, Tmodel} = new()
    function Trace{SampleFromPrior}(m::Model, spl::AbstractSampler, vi::AbstractVarInfo)
        res = new{SampleFromPrior, typeof(vi), typeof(m)}()
        res.vi = vi
        res.model = m
        res.spl = SampleFromPrior()
        return res
    end
    function Trace{T}(m::Model, spl::AbstractSampler, vi::AbstractVarInfo) where T <: Sampler
        res = new{T, typeof(vi), typeof(m)}()
        res.model = m
        res.spl = spl
        res.vi = vi
        return res
    end
end
function Base.copy(trace::Trace)
    newtrace = typeof(trace)()
    newtrace.vi = deepcopy(trace.vi)
    newtrace.task = Base.copy(trace.task)
    newtrace.spl = trace.spl
    newtrace.model = trace.model
    return newtrace
end

# NOTE: this function is called by `forkr`
function Trace(f::Function, m::Model, spl::T, vi::AbstractVarInfo) where {T <: AbstractSampler}
    res = Trace{T}(m, spl, deepcopy(vi));
    # CTask(()->f());
    res.task = CTask( () -> begin res=f(); produce(Val{:done}); res; end )
    if res.task.storage === nothing
        res.task.storage = IdDict()
    end
    res.task.storage[:turing_trace] = res # create a backward reference in task_local_storage
    return res
end
function Trace(m::Model, spl::T, vi::AbstractVarInfo) where {T <: AbstractSampler}
    res = Trace{T}(m, spl, deepcopy(vi));
    # CTask(()->f());
    res.vi.num_produce = 0
    res.task = CTask( () -> begin vi_new=m(vi, spl); produce(Val{:done}); vi_new; end )
    if res.task.storage === nothing
        res.task.storage = IdDict()
    end
    res.task.storage[:turing_trace] = res # create a backward reference in task_local_storage
    return res
end

# step to the next observe statement, return log likelihood
Libtask.consume(t::Trace) = (t.vi.num_produce += 1; consume(t.task))

# Task copying version of fork for Trace.
function fork(trace :: Trace, is_ref :: Bool = false)
    newtrace = copy(trace)
    is_ref && set_retained_vns_del_by_spl!(newtrace.vi, newtrace.spl)
    newtrace.task.storage[:turing_trace] = newtrace
    return newtrace
end

# PG requires keeping all randomness for the reference particle
# Create new task and copy randomness
function forkr(trace :: Trace)
    newtrace = Trace(trace.task.code, trace.model, trace.spl, deepcopy(trace.vi))
    newtrace.spl = trace.spl
    newtrace.vi.num_produce = 0
    return newtrace
end

current_trace() = current_task().storage[:turing_trace]

const Particle = Trace

"""
Data structure for particle filters
- effectiveSampleSize(pc :: ParticleContainer)
- normalise!(pc::ParticleContainer)
- consume(pc::ParticleContainer): return incremental likelihood
"""
mutable struct ParticleContainer{T<:Particle, F, Tvals <: Array{T}, TlogW <: Array{Float64}}
    model :: F
    num_particles :: Int
    vals  :: Tvals
    logWs :: TlogW       # Log weights (Trace) or incremental likelihoods (ParticleContainer)
    logE  :: Float64     # Log model evidence
    # conditional :: Union{Nothing,Conditional} # storing parameters, helpful for implementing rejuvenation steps
    conditional :: Nothing # storing parameters, helpful for implementing rejuvenation steps
    n_consume :: Int # helpful for rejuvenation steps, e.g. in SMC2
end
ParticleContainer{T}(m) where T = ParticleContainer{T}(m, 0)
function ParticleContainer{T}(m, n::Int) where {T}
    ParticleContainer(m, n, Vector{T}(), Vector{Float64}(), 0.0, nothing, 0)
end

Base.collect(pc :: ParticleContainer) = pc.vals # prev: Dict, now: Array
Base.length(pc :: ParticleContainer)  = length(pc.vals)
Base.similar(pc :: ParticleContainer{T}) where T = ParticleContainer{T}(pc.model, 0)
# pc[i] returns the i'th particle
Base.getindex(pc :: ParticleContainer, i :: Real) = pc.vals[i]

# registers a new x-particle in the container
function Base.push!(pc :: ParticleContainer, p :: Particle)
    pc.num_particles += 1
    push!(pc.vals, p)
    push!(pc.logWs, 0)
    pc
end
Base.push!(pc :: ParticleContainer) = Base.push!(pc, eltype(pc.vals)(pc.model))

function Base.push!(pc :: ParticleContainer, n :: Int, spl :: Sampler, varInfo :: VarInfo)
    vals  = Vector{eltype(pc.vals)}(undef,n)
    logWs = zeros(eltype(pc.logWs), n)
    for i=1:n
        vals[i]  = Trace(pc.model, spl, varInfo)
    end
    append!(pc.vals, vals)
    append!(pc.logWs, logWs)
    pc.num_particles += n
    pc
end

# clears the container but keep params, logweight etc.
function Base.empty!(pc :: ParticleContainer)
    pc.num_particles = 0
    pc.vals  = Vector{Particle}()
    pc.logWs = Vector{Float64}()
    pc
end

# clones a theta-particle
function Base.copy(pc :: ParticleContainer)
    particles = collect(pc)
    newpc     = similar(pc)
    for p in particles
        newp = fork(p)
        push!(newpc, newp)
    end
    newpc.logE        = pc.logE
    newpc.logWs       = deepcopy(pc.logWs)
    newpc.conditional = deepcopy(pc.conditional)
    newpc.n_consume   = pc.n_consume
    newpc
end

# run particle filter for one step, return incremental likelihood
function Libtask.consume(pc :: ParticleContainer)
    @assert pc.num_particles == length(pc)
    # normalisation factor: 1/N
    z1 = logZ(pc)
    n = length(pc.vals)

    particles = collect(pc)
    num_done = 0
    for i=1:n
        p = pc.vals[i]
        score = Libtask.consume(p)
        if score isa Real
            score += getlogp(p.vi)
            resetlogp!(p.vi)
            increase_logweight(pc, i, Float64(score))
        elseif score == Val{:done}
            num_done += 1
        else
            error("[consume]: error in running particle filter.")
        end
    end

    if num_done == length(pc)
        res = Val{:done}
    elseif num_done != 0
        error("[consume]: mis-aligned execution traces, num_particles= $(n), num_done=$(num_done).")
    else
        # update incremental likelihoods
        z2 = logZ(pc)
        res = increase_logevidence(pc, z2 - z1)
        pc.n_consume += 1
        # res = increase_loglikelihood(pc, z2 - z1)
    end

    res
end

# compute the normalized weights
weights(pc::ParticleContainer) = softmax(pc.logWs)

# compute the log-likelihood estimate, ignoring constant term ``- \log num_particles``
logZ(pc::ParticleContainer) = logsumexp(pc.logWs)

# compute the effective sample size ``1 / ∑ wᵢ²``, where ``wᵢ```are the normalized weights
function effectiveSampleSize(pc :: ParticleContainer)
    Ws = weights(pc)
    return inv(sum(abs2, Ws))
end

increase_logweight(pc :: ParticleContainer, t :: Int, logw :: Float64) =
    (pc.logWs[t]  += logw)

increase_logevidence(pc :: ParticleContainer, logw :: Float64) =
    (pc.logE += logw)

function resample!(
    pc :: ParticleContainer,
    randcat :: Function = Turing.Inference.resample_systematic,
    ref :: Union{Particle, Nothing} = nothing
)
    n1, particles = pc.num_particles, collect(pc)
    @assert n1 == length(particles)

    # resample
    Ws = weights(pc)

    # check that weights are not NaN
    @assert !any(isnan, Ws)

    n2 = ref === nothing ? n1 : n1 - 1
    indx = randcat(Ws, n2)

    # fork particles
    empty!(pc)
    num_children = zeros(Int,n1)
    for i in indx
        num_children[i] += 1
    end
    for i = 1:n1
        is_ref = particles[i] == ref
        p = is_ref ? fork(particles[i], is_ref) : particles[i]
        num_children[i] > 0 && push!(pc, p)
        for k=1:num_children[i]-1
            newp = fork(p, is_ref)
            push!(pc, newp)
        end
    end

    if isa(ref, Particle)
        # Insert the retained particle. This is based on the replaying trick for efficiency
        #  reasons. If we implement PG using task copying, we need to store Nx * T particles!
        push!(pc, ref)
    end

    pc
end


text:
module IJuliaPortrayals

import Base.writemime
import Base.convert

using Compat

export FromFile
export CSS, IncludeHTML, JavaScript, @JS_str, IFrame
export GIF, JPEG, PNG, SVG
export OGG, MP3, WAV
export YouTube, Vimeo

type FromFile
    path::AbstractString
end

function convert(::Type{AbstractString}, file::FromFile)
    open(file.path) do f
        readall(f)
    end
end

function data_uri(data, mime)
    AbstractString("data:", mime, ";base64,", base64encode(data))
end

embed(path, mime) = open(path) do f
    data_uri(readall(f), mime)
end

# Version 0.4.X seems to have HTML in base/docs/utils.jl
# Version 0.3.X doesnt.
if !isdefined(:HTML)
    export HTML

    type HTML
        doc::AbstractString
    end

    writemime(io::IO, ::MIME"text/html", html::HTML) = print(io, html.doc)
end

type IncludeHTML
    src::AbstractString
end

function writemime(io::IO, ::MIME"text/html", html::IncludeHTML)
    open(html.src) do f
        print(io, readall(f))
    end
end

type JavaScript
    code::AbstractString
end

function writemime(io::IO, ::MIME"text/html", javascript::JavaScript)
    print(io, "")
end

macro JS_str(script)
    script
end

type IFrame
    width::Int
    height::Int
    src::AbstractString
    border::Int
end

IFrame(width, height, src) = IFrame(width, height, src, 0)

function writemime(io::IO, ::MIME"text/html", iframe::IFrame)
    print(io, """

    """)
end

type CSS
    doc::AbstractString
end

function writemime(io::IO, ::MIME"text/html", css::CSS)
    print(io, "")
end

type SVG
    doc::AbstractString
end

writemime(io::IO, ::MIME"image/svg+xml", svg::SVG) = print(io, svg.doc)

type YouTube
    id::AbstractString
    width::Int
    height::Int
end

YouTube(id) = YouTube(id, 420, 315)

function writemime(io::IO, ::MIME"text/html", youtube::YouTube)
    print(io, """

    """)
end

type Vimeo
    id::AbstractString
    width::Int
    height::Int
end

Vimeo(id) = Vimeo(id, 500, 218)

function writemime(io::IO, ::MIME"text/html", vimeo::Vimeo)
    print(io, """

    """)
end

type PNG
    src::AbstractString
end

writemime(io::IO, ::MIME"image/png", img::PNG) = open(img.src) do f
    print(io, readall(f))
end

type GIF
    src::AbstractString
end

function writemime(io::IO, ::MIME"text/html", img::GIF)
    src = embed(img.src, "image/gif")
    print(io, """""")
end

type JPEG
    src::AbstractString
end

writemime(io::IO, ::MIME"image/jpeg", img::JPEG) = open(img.src) do f
    print(io, readall(f))
end

type MP3
    src::AbstractString
    embedded::Bool
end

MP3(src) = MP3(src, false)

function writemime(io::IO, ::MIME"text/html", mp3::MP3)
    src = mp3.embedded ? embed(mp3.src, "audio/mpeg") : mp3.src
    print(io, "<audio controls>")
    print(io, """<source src="$(src)" type="audio/mpeg">""")
    print(io, "</audio>")
end

type OGG
    src::AbstractString
    embedded::Bool
end

OGG(src) = MP3(src, false)

function writemime(io::IO, ::MIME"text/html", ogg::OGG)
    src = ogg.embedded ? embed(ogg.src, "audio/ogg") : ogg.src
    print(io, "<audio controls>")
    print(io, """<source src="$(src)" type="audio/ogg">""")
    print(io, "</audio>")
end

type WAV
    src::AbstractString
    embedded::Bool
end

WAV(src) = WAV(src, false)

function writemime(io::IO, ::MIME"text/html", wav::WAV)
    src = wav.embedded ? embed(wav.src, "audio/mpeg") : wav.src
    print(io, "<audio controls>")
    print(io, """<source src="$(src)" type="audio/wav">""")
    print(io, "</audio>")
end

# The GraphViz portrayal is experimental. There is an existing
# [GraphViz](https://github.com/Keno/GraphViz.jl). It implements a
# full-blown binding, rather than a reading and writing to a subprocess.
# But, I think this portrayal makes sense given the spirit of my package.
export GraphViz

immutable GraphViz
    dot::AbstractString
    engine::AbstractString
    image_type::AbstractString
end

function GraphViz(dot::AbstractString; engine="dot", image_type="svg")
    GraphViz(dot, engine, image_type)
end

const VALID_ENGINES = Set(["neato", "circo", "twopi", "sfdp", "fdp", "dot"])
const VALID_IMAGE_TYPES = Set(["svg", "png"])

function build_command(gv::GraphViz)
    gv.engine ∈ VALID_ENGINES || error("'$(gv.engine)' is not a valid layout")
    gv.image_type ∈ VALID_IMAGE_TYPES || error("'$(gv.image_type)' is not a valid format")

    `dot -K$(gv.engine) -T$(gv.image_type)`
end

function Base.writemime(io::IO, ::MIME"text/html", gv::GraphViz)
    (gv_out, gv_in, gv_process) = readandwrite(build_command(gv))

    write(gv_in, gv.dot)
    close(gv_in)

    emitted = readall(gv_out)
    close(gv_out)

    if gv_process.exitcode == 0
        if(gv.image_type == "svg")
            print(io, emitted)
        else
            data = data_uri(emitted, "image/png")
            print(io, """""")
        end
    else
        error("There was an error in your dot file syntax.")
    end
end

end

ext: jl


text:
"""
    FEMMDeforWinklerModule

Module for operations on boundaries of domains to construct system matrices and
system vectors for linear deformation models with distributed-spring supports
(Winkler foundation model).
"""
module FEMMDeforWinklerModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FENodeSetModule: FENodeSet
import FinEtools.FESetModule: FESet, nodesperelem, manifdim
import FinEtools.IntegDataModule: IntegData, integrationdata, Jacobiansurface
import FinEtools.FieldModule: ndofs, gatherdofnums!
import FinEtools.NodalFieldModule: NodalField
import FinEtools.FEMMBaseModule: FEMMAbstractBase
import FinEtools.AssemblyModule: SysvecAssemblerBase, SysmatAssemblerBase, SysmatAssemblerSparseSymm, startassembly!, assemble!, makematrix!
import FinEtools.MatrixUtilityModule: add_nnt_ut_only!, complete_lt!, locjac!
import FinEtools.SurfaceNormalModule: SurfaceNormal, updatenormal!
import LinearAlgebra: norm, cross

"""
    FEMMDeforWinkler{S<:FESet, F<:Function}

Type for normal spring support  (Winkler).
"""
mutable struct FEMMDeforWinkler{S<:FESet, F<:Function} <: FEMMAbstractBase
    integdata::IntegData{S, F} # geometry data
end

"""
    surfacenormalspringstiffness(self::FEMMDeforWinkler, assembler::A,
                  geom::NodalField{FFlt}, u::NodalField{T},
                  springconstant::FFlt) where {A<:SysmatAssemblerBase, T<:Number}

Compute the stiffness matrix of surface normal spring.
Rationale: consider continuously distributed springs between the surface of the
solid body and the 'ground', in the direction normal to the surface. If the
spring coefficient becomes large, we have an approximate method of enforcing the
normal displacement to the surface.gas
"""
function surfacenormalspringstiffness(self::FEMMDeforWinkler, assembler::A,
    geom::NodalField{FFlt}, u::NodalField{T},
    springconstant::FFlt, surfacenormal::SurfaceNormal) where {A<:SysmatAssemblerBase, T<:Number}
    integdata = self.integdata
    # Constants
    nfes = count(integdata.fes); # number of finite elements in the set
    ndn = ndofs(u); # number of degrees of freedom per node
    nne = nodesperelem(integdata.fes); # number of nodes for element
    sdim = ndofs(geom);            # number of space dimensions
    mdim = manifdim(integdata.fes); # manifold dimension of the element
    Kedim = ndn*nne;             # dimension of the element matrix
    # Precompute basis f. values + basis f. gradients wrt parametric coor
    npts, Ns, gradNparams, w, pc = integrationdata(integdata);
    # Prepare assembler and temporaries
    Ke = zeros(FFlt,Kedim,Kedim);                # element matrix -- used as a buffer
    dofnums = zeros(FInt, Kedim); # degree of freedom array -- used as a buffer
    loc = zeros(FFlt, 1,sdim); # quadrature point location -- used as a buffer
    J = fill(zero(FFlt), sdim,mdim); # Jacobian matrix -- used as a buffer
    Nn = zeros(FFlt, Kedim); # column vector
    startassembly!(assembler, Kedim, Kedim, nfes, u.nfreedofs, u.nfreedofs);
    for i = 1:nfes # Loop over elements
        fill!(Ke, 0.0); # Initialize element matrix
        for j = 1:npts # Loop over quadrature points
            locjac!(loc, J, geom.values, integdata.fes.conn[i], Ns[j], gradNparams[j])
            Jac = Jacobiansurface(integdata, J, loc, integdata.fes.conn[i], Ns[j]);
            n = updatenormal!(surfacenormal, loc, J, integdata.fes.label[i])
            for k= 1:nne
                for r = 1:sdim
                    Nn[(k-1)*sdim+r] = n[r] * Ns[j][k];
                end
            end
            add_nnt_ut_only!(Ke, Nn, springconstant*Jac*w[j])
        end # Loop over quadrature points
        complete_lt!(Ke)
        gatherdofnums!(u, dofnums, integdata.fes.conn[i]);# retrieve degrees of freedom
        assemble!(assembler, Ke, dofnums, dofnums);# assemble symmetric matrix
    end # Loop over elements
    return makematrix!(assembler);
end

function surfacenormalspringstiffness(self::FEMMDeforWinkler,
              geom::NodalField{FFlt}, u::NodalField{T},
              springconstant::FFlt, surfacenormal::SurfaceNormal) where {T<:Number}
    assembler = SysmatAssemblerSparseSymm();
    return surfacenormalspringstiffness(self, assembler, geom, u, springconstant, surfacenormal);
end

# """
#     surfacenormal(loc::FFltMat,J::FFltMat)

# Compute local normal. This makes sense for bounding surfaces only.
# """
# function  surfacenormal(loc::FFltMat, J::FFltMat)
#     norml= zeros(FFlt, length(loc))
#     # Produce a default normal
#     if (size(J,1)==3) && (size(J,2)==2)# surface in three dimensions
#         norml = cross(vec(J[:,1]),vec(J[:,2]));# outer normal to the surface
#         norml = norml/norm(norml);
#     elseif (size(J,1)==2)  && (size(J,2)==1)# curve in two dimensions
#         norml= [J[2,1];-J[1,1]];# outer normal to the contour
#         norml = norml/norm(norml);
#     else
#         error("No definition of normal vector");
#     end
#     return reshape(norml,length(norml),1) # return a column vector
# end

end


text:
doc"""
    FDist(ν1,ν2)

The *F distribution* has probability density function

$f(x; \nu_1, \nu_2) = \frac{1}{x B(\nu_1/2, \nu_2/2)}
\sqrt{\frac{(\nu_1 x)^{\nu_1} \cdot \nu_2^{\nu_2}}{(\nu_1 x + \nu_2)^{\nu_1 + \nu_2}}},
\quad x>0$

It is related to the [`Chisq`](:func:`Chisq`) distribution via the property that if $X_1
\sim \operatorname{Chisq}(\nu_1)$ and $X_2 \sim \operatorname{Chisq}(\nu_2)$, then
$(X_1/\nu_1) / (X_2 / \nu_2) \sim FDist(\nu_1, \nu_2)`.

```julia
FDist(d1, d2)     # F-Distribution with parameters d1 and d2

params(d)         # Get the parameters, i.e. (d1, d2)
```

External links

* [F distribution on Wikipedia](http://en.wikipedia.org/wiki/F-distribution)
    """
immutable FDist <: ContinuousUnivariateDistribution
    ν1::Float64
    ν2::Float64

    function FDist(ν1::Real, ν2::Real)
        @check_args(FDist, ν1 > zero(ν1) && ν2 > zero(ν2))
        new(ν1, ν2)
    end
end

@distr_support FDist 0.0 Inf

#### Parameters

params(d::FDist) = (d.ν1, d.ν2)

#### Statistics

mean(d::FDist) = (ν2 = d.ν2; ν2 > 2.0 ? ν2 / (ν2 - 2.0) : NaN)

mode(d::FDist) = ((ν1, ν2) = params(d); ν1 > 2.0 ? ((ν1 - 2.0)/ν1) * (ν2 / (ν2 + 2.0)) : 0.0)

function var(d::FDist)
    (ν1, ν2) = params(d)
    ν2 > 4.0 ? 2.0 * ν2^2 * (ν1 + ν2 - 2.0) / (ν1 * (ν2 - 2.0)^2 * (ν2 - 4.0)) : NaN
end

function skewness(d::FDist)
    (ν1, ν2) = params(d)
    if ν2 > 6.0
        return (2.0 * ν1 + ν2 - 2.0) * sqrt(8.0 * (ν2 - 4.0)) / ((ν2 - 6.0) * sqrt(ν1 * (ν1 + ν2 - 2.0)))
    else
        return NaN
    end
end

function kurtosis(d::FDist)
    (ν1, ν2) = params(d)
    if ν2 > 8.0
        a = ν1 * (5. * ν2 - 22.) * (ν1 + ν2 - 2.) + (ν2 - 4.) * (ν2 - 2.)^2
        b = ν1 * (ν2 - 6.) * (ν2 - 8.) * (ν2 - 2.)
        return 12. * a / b
    else
        return NaN
    end
end

function entropy(d::FDist)
    (ν1, ν2) = params(d)
    hν1 = ν1 * 0.5
    hν2 = ν2 * 0.5
    hs = (ν1 + ν2) * 0.5
    return log(ν2 / ν1) + lgamma(hν1) + lgamma(hν2) - lgamma(hs) +
        (1.0 - hν1) * digamma(hν1) + (-1.0 - hν2) * digamma(hν2) +
        hs * digamma(hs)
end

#### Evaluation & Sampling

@_delegate_statsfuns FDist fdist ν1 ν2

rand(d::FDist) = StatsFuns.Rmath.fdistrand(d.ν1, d.ν2)


text:
export SimplexPoint, transferoperator
import StaticArrays: SizedVector, SizedMatrix, MMatrix

include("helper_functions.jl")

"""
    SimplexPoint() <: TriangulationBasedTransferOperator

A transfer operator estimator using a triangulation partition and approximate simplex intersections[^Diego2019].

To use this estimator, the Simplices.jl package must be brought into scope by doing `using Simplices` after
`using PerronFrobenius`.

[^Diego2019]: Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212.
"""
struct SimplexPoint <: TriangulationBasedTransferOperator
    bc::String

    function SimplexPoint(bc::String = "circular")
        isboundarycondition(bc, "triangulation")  || error("Boundary condition '$bc' not valid.")
        new(bc)
    end
end
Base.show(io::IO, se::SimplexPoint) = print(io, "SimplexPoint{$(se.bc)}")

""" Generate a transopergenerator for an exact simplex estimator."""
function transopergenerator(pts, method::SimplexPoint)
    # modified points, where the image of each point is guaranteed to lie within the convex hull of the previous points
    invariant_pts = invariantize(pts)

    # triangulation of the invariant points. The last point is excluded, so that
    # the last vertex also can be mapped forward one step in time.
    triang = DelaunayTriangulation(invariant_pts[1:end-1])

    init = (invariant_pts = invariant_pts, triang = triang)

    TransferOperatorGenerator(method, pts, init)
end

"""
n: The number of points sampled inside each simplex.
randomsampling: Sample randomly (`random == true`), or using an even simplex splitting routine (`random == false`)
"""
function (tog::TransferOperatorGenerator{<:SimplexPoint})(; tol = 1e-8, randomsampling::Bool = false, n::Int = 100)
    pts, inds = getfield.(Ref(tog.init), (:invariant_pts, :triang))
    D = length(pts[1])
    N = length(inds)
    ϵ = tol / N

    # Prepare memory-efficient representations of the simplices
    simplices = reshape_simplices(pts[1:(end- 1)], inds)
    image_simplices = reshape_simplices(pts[2:end], inds)

    #=
    Compute the convex expansions coefficients needed to generate points
    within a simplex. Also, update number of points if a shape-preserving
    simplex subdivision algorithm was used (in that case, we get a
    *minimum* of `n_sample_pts` set of coefficients.
    =#
    convex_coeffs = SimplexTypes.subsample_coeffs(D, n, randomsampling)
    n_coeffs::Int = size(convex_coeffs, 2)
    convex_coeffs = [convex_coeffs[:, i] for i = 1:n_coeffs]

    # Pre-allocated arrays (SizedArrays, for efficiency)
    L = (D+1)^2
    NVERTICES = D+1
    pt          = SizedMatrix{1, D}(zeros(1, D))
    s_arr       = SizedVector{L}(zeros(L))
    signs       = SizedVector{NVERTICES}(zeros(NVERTICES))

    # Re-arrange simplices so that look-up is a bit more efficient
    simplex_arrs = Vector{Array{Float64, 2}}(undef, N)
    imsimplex_arrs = Vector{Array{Float64, 2}}(undef, N)
    for i in 1:N
        simplex_arrs[i] = hcat(pts[inds[i]]...,)
        imsimplex_arrs[i] = hcat(pts[inds[i] .+ 1]...,)
    end

    # Pre-allocate some information to avoid memory allocation when
    # computing the fractional overlaps. These are a temporary matrix
    # which we use to compute determinants (to check if a point lies
    # inside a simplex), and the indices needed to access its values
    # at different stages of the checking stage (of whether a point
    # lies inside a simplex). Pre-allocating this information and
    # passing it on reduces runtime and memory allocation
    # by half an order of magnitude.
    temp_arr = rand(MMatrix{D + 1, D + 1})
    starts = [((D + 1)*i)+1 for i = 0:D]
    stops = starts .+ D

    # The Markov matrix
    M = zeros(Float64, N, N)

    for i in 1:N
        # Find the indices of the simplices who potentially intersect
        # with the image of simplex #i
        idxs = idxs_potentially_intersecting_simplices(pts, inds, i)
        Sj = Vector{AbstractArray}(undef, length(idxs))
        get_simplices_at_inds!(Sj, idxs, simplices)

        # Generate points within the image of simplex #i and check
        # what fraction of those points fall into the simplices.
        # We approximate the intersecting volume between the
        # image simplex and the individual simplices by those
        # fractions.
        @views is = imsimplex_arrs[i]
        for k in 1:n_coeffs
            pt = transpose(convex_coeffs[k]) * transpose(is)
            innerloop_optim!(idxs, signs, s_arr, Sj, pt, D, M, i, temp_arr, starts, stops)
        end
    end
    # Need to normalise, because all we have up until now is counts
    # of how many points inside the image simplex falls into
    # the simplices.
    params = (tol = tol, randomsampling = randomsampling, n = n)
    M_normalized = transpose(M) ./ n_coeffs
    return TransferOperatorApproximation(tog, M_normalized, params)

end

function transferoperator(pts, method::SimplexPoint; tol = 1e-8, randomsampling::Bool = false, n::Int = 100)
    tog = transopergenerator(pts, method)
    tog(tol = tol, randomsampling = randomsampling, n = n)
end


text:
using NCDatasets
using Test

fname = "test_cont_ragged_array.nc"
ds = NCDataset(fname,"c")
# Dimensions

ds.dim["obs"] = 7
ds.dim["profile"] = 3

# Declare variables

ncprofile = defVar(ds,"profile", Int32, ("profile",))
ncprofile.attrib["cf_role"] = "profile_id"

nctime = defVar(ds,"time", Float64, ("profile",))
nctime.attrib["standard_name"] = "time"
nctime.attrib["long_name"] = "time"
nctime.attrib["units"] = "days since 1970-01-01 00:00:00"

nclon = defVar(ds,"lon", Float32, ("profile",))
nclon.attrib["standard_name"] = "longitude"
nclon.attrib["long_name"] = "longitude"
nclon.attrib["units"] = "degrees_east"

nclat = defVar(ds,"lat", Float32, ("profile",))
nclat.attrib["standard_name"] = "latitude"
nclat.attrib["long_name"] = "latitude"
nclat.attrib["units"] = "degrees_north"

ncrowSize = defVar(ds,"rowSize", Int32, ("profile",))
ncrowSize.attrib["long_name"] = "number of obs for this profile "
ncrowSize.attrib["sample_dimension"] = "obs"

ncz = defVar(ds,"z", Float32, ("obs",))
ncz.attrib["standard_name"] = "altitude"
ncz.attrib["long_name"] = "height above mean sea level"
ncz.attrib["units"] = "km"
ncz.attrib["positive"] = "up"
ncz.attrib["axis"] = "Z"

ncpressure = defVar(ds,"pressure", Float32, ("obs",))
ncpressure.attrib["standard_name"] = "air_pressure"
ncpressure.attrib["long_name"] = "pressure level"
ncpressure.attrib["units"] = "hPa"
ncpressure.attrib["coordinates"] = "time lon lat z"

nctemperature = defVar(ds,"temperature", Float32, ("obs",))
nctemperature.attrib["standard_name"] = "surface_temperature"
nctemperature.attrib["long_name"] = "skin temperature"
nctemperature.attrib["units"] = "Celsius"
nctemperature.attrib["coordinates"] = "time lon lat z"

nchumidity = defVar(ds,"humidity", Float32, ("obs",))
nchumidity.attrib["standard_name"] = "relative_humidity"
nchumidity.attrib["long_name"] = "relative humidity"
nchumidity.attrib["units"] = "%"
nchumidity.attrib["coordinates"] = "time lon lat z"

# Global attributes

ds.attrib["featureType"] = "profile"

# Define variables

ncprofile[:] = [1,2,3]
nctime[:] = [1,1,1]
nclon[:] = [0,0,0]
nclat[:] = [0,0,0]
ncrowSize[:] = [3,2,2]
ncz[:] = [1.,2.,3.,  10.,20.,  100., 200]
ncpressure[:] = [1.,2.,3.,  10.,20.,  100., 200]
nctemperature[:] = [1.,2.,3.,  10.,20.,  100., 200]
nchumidity[:] = [1.,2.,3.,  10.,20.,  100., 200]

close(ds)

ds = NCDataset(fname);
ncvar = ds["z"]
data = loadragged(ncvar,:)
@test data == [[1.,2.,3.],  [10.,20.],  [100., 200.]]
close(ds)

rm(fname)


text:
"""
  scalarDirichlet!(dirichletNodes,dBCarr,Stiff,F)

Applies Dirichlet boundary conditions for Scalar equation
"""
function scalarDirichlet!(dirichletNodes,dBCarr,Stiff,F)
  Nnodes  = length(F)
  dLength = length(dirichletNodes)

  for i=1:dLength
    d = dirichletNodes[i]

    for j=1:Nnodes
      #for j in rvals, adjust RHS with dirichlet BC
      F[j] -= Stiff[j,d]*dBCarr[i]

      # zero out (dirichlet) BC locations
      Stiff[j,d] = 0.0
      Stiff[d,j] = 0.0
    end

    # identity in correct location
    Stiff[d,d] = 1.0

    # set dirichlet nodes to exact boundary condition
    F[d] = dBCarr[i]
  end
end

"""
  scalarNeumann!(xy,cm,order,neumannNodes,nBCarr,F)

Applies Neumann boundary conditions for arbitrary scalar equation
"""
function scalarNeumann!(xy,cm,order,neumannNodes,nBCarr,F)
  Nel = size(cm,1)
  gaussorder=0

  # generate basis values at gauss points
  if order == :Linear
    order=1
    gaussorder = 2
  elseif order == :Quadratic
    order=2
    gaussorder = 3
  end

  if length(neumannNodes)>0
    for el=1:Nel
      xNodes = [xy[cm[el].NodeList[i]].x for i in 1:Int(gaussorder^2)]
      yNodes = [xy[cm[el].NodeList[i]].y for i in 1:Int(gaussorder^2)]
      # loop over ONLY QUAD NODES of each element
      for i=1:4
        # check that nodes i is a neumann Boundary node
        if cm[el].NodeList[i] in neumannNodes
          # check that next node i+1 is a neumann Boundary node
          j = mod(i,4)+1
          if cm[el].NodeList[j] in neumannNodes
            ci = cm[el].NodeList[i]
            cj = cm[el].NodeList[j]

            w,s,t = GaussEdge(i,gaussorder)

            jac = sqrt((xNodes[i]-xNodes[j])^2+(yNodes[i]-yNodes[j])^2)

            ki = findfirst(s->s==ci,neumannNodes)
            kj = findfirst(s->s==cj,neumannNodes)

            fluxi = nBCarr[ki]
            fluxj = nBCarr[kj]

            for gpt=1:(order+1)
              phi,_,_,_ = derivShape2D(s[gpt],t[gpt],xNodes,yNodes,order)
              fluxVal = shapeEval([fluxi,fluxj],[phi[i],phi[j]])

              F[ci] += w[gpt]*fluxVal*phi[i]*jac
              F[cj] += w[gpt]*fluxVal*phi[j]*jac
            end
          end
        end
      end
    end
  end
end

"""
  scalarBC!(xy,cm,order,dirichletNodes,neumannNodes,
                   dBCarr,nBCarr,Stiff,F)

Applies Dirichlet and Neumann boundary conditions to arbitrary scalar equation
"""
function scalarBC!(xy,cm,order,dirichletNodes,neumannNodes,
                   dBCarr,nBCarr,Stiff,F)
  scalarDirichlet!(dirichletNodes,dBCarr,Stiff,F)
  scalarNeumann!(xy,cm,order,neumannNodes,nBCarr,F)
end

"""
  fluidUDirichlet!(xy,xyp,dUNodes,dUBCarr,Stiff,F)

Applies Dirichlet boundary conditions for horizontal velocity components for Stokes's equation.

Updates the Stiff matrix and F vector
"""
function fluidUDirichlet!(xy,xyp,dUNodes,dUBCarr,Stiff,F)
  nu = length(xy)
  np = length(xyp)
  Nnodes  = length(F)
  uLength = length(dUNodes)

  d = dUNodes
  s1 = 1:nu
  s2 = (nu+1):2*nu
  s3 = (2*nu+1):(2*nu+np)

  s12 = 1:2*nu

  F[s1] -= Stiff[s1,d]*dUBCarr
  F[s2] -= Stiff[s2,d]*dUBCarr
  F[s3] -= Stiff[s3,d]*dUBCarr

  for i=1:uLength
    for j in s12
      Stiff[j,d[i]] = 0.0
      Stiff[d[i],j] = 0.0
    end

    # identity in correct location
    Stiff[d[i],d[i]] = 1.0

    # set dirichlet nodes to exact boundary condition
    F[d[i]] = dUBCarr[i]

    for j in s3
      Stiff[d[i],j] = 0.0
      Stiff[j,d[i]] = 0.0
    end
  end
end

"""
  fluidVDirichlet!(xy,xyp,dVNodes,dVBCarr,Stiff,F)

Applies Dirichlet boundary conditions for vertical velocity components for Stokes's equation.

Updates the Stiff matrix and F vector
"""
function fluidVDirichlet!(xy,xyp,dVNodes,dVBCarr,Stiff,F)
  nu = length(xy)
  np = length(xyp)
  vLength = length(dVNodes)

  d  = nu .+ dVNodes
  s1 = 1:nu
  s2 = (nu+1):2*nu
  s3 = (2*nu+1):(2*nu+np)

  s12 = 1:2*nu

  F[s1] -= Stiff[s1,d]*dVBCarr
  F[s2] -= Stiff[s2,d]*dVBCarr
  F[s3] -= Stiff[s3,d]*dVBCarr

  for i=1:vLength
    for j in s12
      Stiff[j,d[i]] = 0.0
      Stiff[d[i],j] = 0.0
    end

    # identity in correct location
    Stiff[d[i],d[i]] = 1.0

    # set dirichlet nodes to exact boundary condition
    F[d[i]] = dVBCarr[i]

    for j in s3
      Stiff[d[i],j] = 0.0
      Stiff[j,d[i]] = 0.0
    end
  end
end

"""
  fluidDirichlet!(xy,xyp,dUNodes,dVNodes,dUBCarr,dVBCarr,Stiff,F)

computes u and v dirichlet conditions for fluid equation
"""
function fluidDirichlet!(mesh::FluidMesh,prob::Problem,LinOp::LinearOperator)
  # get xy into correct form
  dUNodes = prob.bcid[:dUNodes]
  dUBCarr = prob.bcval[:dUBC]
  dVNodes = prob.bcid[:dVNodes]
  dVBCarr = prob.bcval[:dVBC]

  fluidUDirichlet!(mesh.xy,mesh.xyp,dUNodes,dUBCarr,LinOp.Op,LinOp.rhs)
  fluidVDirichlet!(mesh.xy,mesh.xyp,dVNodes,dVBCarr,LinOp.Op,LinOp.rhs)
end

"""
  fluidNeumann(xy,nNodes,nBCarr,Stiff,F) -> Stiff, F

Because we have assumed nothing else than zero neumann boundary conditions, i.e. 'outflow' conditions, this just checks whether there are no outflow conditions, in which case it adds an additional equation to normalize the pressure about 0.
"""
function fluidNeumann(mesh::FluidMesh,prob::Problem,LinOp::LinearOperator)
  if !(:nNodes in bc(prob))
    NumUnodes = length(mesh.xy)
    NumPnodes = length(mesh.xyp)
    ucol = spzeros(2*NumUnodes,1)
    pcol = ones(NumPnodes,1)/NumPnodes
    urow = ucol'
    prow = pcol'

    s1 = vcat(ucol,pcol)
    s2 = hcat(urow,prow,1)

    LinOp.Op = hcat(LinOp.Op,s1)
    LinOp.Op = vcat(LinOp.Op,s2)

    push!(LinOp.rhs,0)
  end

  return LinOp.Op,LinOp.rhs
end

"""
  fluidBC!(xy,xyp,dUNodes,dVNodes,dUBCarr,dVBCarr,nNodes,nBCarr,Stiff,F)
    -> Stiff, F

Computes Dirichlet and Neumann boundary conditions
"""
function fluidBC!(mesh::FluidMesh,prob::Problem,
                  LinOp::AbstractLinearOperator)

  fluidDirichlet!(mesh::FluidMesh,prob::Problem,LinOp::LinearOperator)
  LinOp.Op, LinOp.rhs = fluidNeumann(mesh::FluidMesh,prob::Problem,
                                   LinOp::LinearOperator)
end


text:
function testLinear(f::F, m::Int64, p::Float64, C::Float64) where F<:Function
  ϵ::Float64 = (1.0 - C*p)/2.0
  testFunction(f, m, C*p, makef(p), C, ϵ, GLOBAL_RNG)
end

Random.seed!(192837465)

testLinear(BernoulliFactory.linear, ntrials, 0.1, 2.0)
testLinear(BernoulliFactory.linear, ntrials, 0.4, 2.0)
testLinear(BernoulliFactory.linear, ntrials, 0.1, 3.0)
testLinear(BernoulliFactory.linear, ntrials, 0.6, 1.5)
testLinear(BernoulliFactory.linear, ntrials, 0.6, 0.4)

include("huber2017_test.jl")


text:
"""
    first_po(n::Int64, k::Int64) -> po_idx::Vector{Int64}

Notation-Generator: Calculate the first generator matrix in PO notation (padding ones indices) by n and k.
"""
function first_po(n::Int64, k::Int64)::Vector{Int64}
    return collect(1:(n-k))
end

"""
    next_po(current_po::Vector{Int64}, n::Int64, i=0) -> next::Vector{Int64}

Notation-Generator: Calculates the next generator matrix in PO notation (padding ones indices) in relation to a given generator matrix in PO notation.
"""
function next_po(current_po::Vector{Int64}, n::Int64, i=0)::Vector{Int64}
    if i > size(current_po,1)
        println("Error: i=$i")
        return nothing
    end
    next = copy(current_po)
    if current_po[end-i] + 1 > n - i
        return next_po(current_po, n, i+1)
    else
        next[end-i] += 1
        for j in (size(next,1) - i + 1):size(next,1)
           next[j] = next[j-1] +1
        end
        return next
    end
end

"""
    construct_generators(n::Int64, k::Int64, prnt = false) -> Ω::Int64, all_pos, all_generating_matrices

Matrix- and Notation-Generator: Calculate all generating matrices (including PO notations) and the number of generating matrices Ω.
"""
function construct_generators(n::Int64, k::Int64, prnt = false)
    # calculate number of possible generating matrices
    Ω = binomial(n,k)
    println("There are Ω = ($n choose $k) = $Ω  possible basic permutations.")

    # Generate all possible padding ones positions first
    # number of padding ones
    num_pos = n - k
    # save first PO
    all_pos = []
    push!(all_pos, first_po(n,k))
    # save first PO as matrix
    all_generating_matrices = SparseMatrixCSC{ℤ₂, Int64}[]
    push!(all_generating_matrices, po2matrix(all_pos[1],n))

    if prnt
        println(
            "1:\t ", all_pos[1],
            "\t↔ AI:\t ", po2ai(all_pos[1],n),
            "\tA² Test: ", multiplyZ2(all_generating_matrices[1],all_generating_matrices[1]) == id(n),
            ": ", all_generating_matrices[1], "\n\n")
    end

    # construct all other generating matrices
    for i in 2:Ω
        next = next_po(all_pos[end],n)
        push!(all_pos, next)
        push!(all_generating_matrices, po2matrix(all_pos[end],n))
        if prnt
            println(
                "$i:\t", next,
                "\t↔ AI:\t ", po2ai(all_pos[end],n),
                "\tA² Test: ", multiplyZ2(all_generating_matrices[end],all_generating_matrices[end]) == id(n),
                ": ", all_generating_matrices[end], "\n\n")
        end
    end

    # return all necessary variables
    return Ω, all_pos, all_generating_matrices
end


text:
using Compat

"""
Any entity that is tracked through the pooled experiment. For CRISPR screens,
this is equivalent to the sgRNA. This object stores properties relating to the
performance of this entity in the screen.

$(FIELDS)
"""
type Barcode
    "The target gene id"
    gene::Int

    "The knockdown efficiency"
    knockdown::Float64

    "The theoretical phenotype exerted"
    theo_phenotype::Float64

    "The observed phenotype (only relevant for FACS screens)"
    obs_phenotype::Float64

    "The gene behavior, either sigmoidal or linear"
    behavior::Symbol

    "The gene activity class, either inactive, up or down"
    class::Symbol

    "Initial frequency after transfection"
    initial_freq::Float64

    function Barcode(gene, knockdown, phenotype, behavior, class)
        if (0 <= knockdown <= 1)
            new(gene, knockdown, phenotype, NaN, behavior, class, -Inf)
        else
            error("Knockdown must be between 0 and 1, inclusive")
        end
    end
end

as_array(bc::Barcode) = [getfield(bc, fld) for fld in fieldnames(bc) if fld != :obs_phenotype]

"""
A description of screen parameters
"""
@compat abstract type ScreenSetup end

as_array(ss::ScreenSetup) = [getfield(ss, fld) for fld in fieldnames(ss)]

"""
A type representing the parameters used in a typical FACS screen.

$(FIELDS)
"""
type FacsScreen <: ScreenSetup
    "Number of genes targeted by the screen"
    num_genes::Int

    "Number of guides per gene"
    coverage::Int

    """Number of cells with each guide. `representation=10` means that there are
    10 times as many cells as guides during transfection. The actual number of
    cells per guide post-transfection will be less depending on the MOI"""
    representation::Int

    """The multiplicity of infection, ``\\lambda``, of the screen. We model this as a Poisson
    process during transfection (see [`Simulation.transfect`](@ref)).

    !!! note

        We **do not** model multiple infections. We assume that the MOI is properly
        selected and less than half the cells are transfected by any virus, i.e.
        ``\\lambda \\lt 0.5`` and then select only the cells that have a single
        transfection occurrence:

        ```math
        P(x = 1; Poisson(λ))
        ```

        For ``\\lambda = 0.25`` this works out to being ``\\approx 19.5\\%`` of the
        number of cells (`num_genes` * `coverage` * `representation`).
    """
    moi::Float64

    """The standard deviation expected for cells during FACS sorting. This should
    be set according to the biological variance experimentally observed, e.g. in
    the fluorescence intensity of isogenic cells"""
    σ::Float64

    """Range of guide phenotypes to collect in each bin

    # Example

    In the following example

    ```julia
    p = 0.05
    bin_info = Dict(:bin1 => (0.0, p), :bin2 => (1.0-p, 1.0))
    ```

    The 5th percentile of cells sorted according to their phenotype (fluorescence,
    size, etc) will be compared to the 95th percentile.
    """
    bin_info::Associative{Symbol, Tuple{Float64, Float64}}

    "Number of cells sorted expressed as an integer multiple of the number of guides"
    bottleneck_representation::Int

    """Sequencing depth as a integer multiple of the number of guides, i.e.
    `seq_depth=10` is equivalent to `10 * num_genes * coverage` reads.
    """
    seq_depth::Int

    function FacsScreen()
        new(500, 5, 100, 0.25, 1.0, OrderedDict(:bin1 => (0.0, 1/3), :bin2 => (2/3, 1.0)), 1000, 1000)
    end
end

"""
A type representing the parameters used in a typical growth-based screen.

$(FIELDS)
"""
type GrowthScreen <: ScreenSetup
    "Number of genes targeted by the screen"
    num_genes::Int

    "Number of guides per gene"
    coverage::Int

    """Number of cells with each guide. `representation=10` means that there are
    10 times as many cells as guides during transfection. The actual number of
    cells per guide post-transfection will be less depending on the MOI"""
    representation::Int

    """The multiplicity of infection, ``\\lambda``, of the screen. We model this as a Poisson
    process during transfection (see [`Simulation.transfect`](@ref)).

    !!! note

        We **do not** model multiple infections. We assume that the MOI is properly
        selected and less than half the cells are transfected by any virus, i.e.
        ``\\lambda \\lt 0.5`` and then select only the cells that have a single
        transfection occurrence:

        ```math
        P(x = 1; Poisson(λ))
        ```

        For ``\\lambda = 0.25`` this works out to being ``\\approx 19.5\\%`` of the
        number of cells (`num_genes` * `coverage` * `representation`).
    """
    moi::Float64

    """Sequencing depth as a integer multiple of the number of guides, i.e.
    `seq_depth=10` is equivalent to `10 * num_genes * coverage` reads.
    """
    seq_depth::Int

    """For growth screens, how much of a bottleneck is applied. This is minimum
    number of cells that is kept when passaging the pool of cells. This is
    expressed as an integer multiple of the number of guides."""
    bottleneck_representation::Int

    """For growth screens, how many bottlenecks are applied. This is the integer
    number of passages during the growth screen."""
    num_bottlenecks::Int

    """Before each passage, the theoretical phenotype of each cell is convolved
    with a normal noise distribution with a standard deviation, σ, dictated by
    this parameter. This should be set based on an expectation of noisiness of
    the subsampling"""
    noise::Float64

    function GrowthScreen()
        new(500, 5, 100, 0.25, 1000, 1000, 10, 0.01)
    end
end

function _print(io::IO, s::ScreenSetup)
    print(io, typeof(s), " <: ScreenSetup \n")

    info = (
        "Number of genes: $(s.num_genes)",
        "Number of guides per gene: $(s.coverage) ($(s.num_genes*s.coverage) guides total)",
        "Cells per guide at transfection: $(s.representation) ($(s.num_genes*s.coverage*s.representation) cells total)",
        "Multiplicity of infection (M.O.I.): $(s.moi)",
        "Number of cells per guide at sequencing: $(s.seq_depth) ($(s.num_genes*s.coverage*s.seq_depth) cells total)",
    )

    for line in info
        print(io, "\n    ", line)
    end
end

function Base.show(io::IO, s::FacsScreen)
    _print(io::IO, s)

    info = (
        "",
        "Number of cells per guide sorted: $(s.bottleneck_representation) ($(s.num_genes*s.coverage*s.bottleneck_representation) cells total)",
        "Gaussian standard deviation of sorting: $(s.σ) (phenotype units)",
        "Bins:"
    )

    for line in info
        print(io, "\n    ", line)
    end

    for (bin, range) in s.bin_info
        print(io, "\n        $bin: $(round(Int, range[1]*100))-$(round(Int, range[2]*100)) percentile of cells",)
    end
end

function Base.show(io::IO, s::GrowthScreen)
    _print(io::IO, s)

    info = (
        "",
        "Minimum number of cells kept per subsampling: $(s.bottleneck_representation) ($(s.num_genes*s.coverage*s.bottleneck_representation) cells total)",
        "Number of rounds of subsampling: $(s.num_bottlenecks)",
        "Gaussian standard deviation of subsampling: $(s.noise) (phenotype units)",
    )

    for line in info
        print(io, "\n    ", line)
    end
end


text:
function meanstd( img )

	T  = ( eltype(img) <: Integer ) ? Int64 : Float32;

	I  = T(0)
	I2 = T(0)
	N  = length(img);
	@inbounds @simd for idx in 1:N
		I  += img[idx]
		I2 += img[idx]*img[idx]
	end
	# mean = sumI / N
	# std  : sum( ( I - mean )^2 ) = sum( I^2 + mean^2 - 2*I*mean ) = sumI2 + Nmean^2 - 2mean*sumI

	mean  = I/N
	mean2 = N*mean*mean
	std   = sqrt( ( I2 + mean2 - I*2*mean )/(N-1) )

	return mean*N, std
end

function linearBinSearch( array::Array{T1,1}, query::T2 ) where {T1<:Number,T2<:Number}

    idx = 0;
    @inbounds for e in 1:length(array)
        query <= array[ e ] && break
        idx = e;
    end

    return idx
end

function binaryBinSearch( array::Array{T1,1}, query::T2 ) where {T1<:Number,T2<:Number}

    hi, lo, probe = length( array ) + 1, 0, 0;

    @inbounds while 1 < hi - lo
        probe = ( hi + lo ) >>> 1
        if array[probe] < query
            lo = probe
        else
            hi = probe
        end
    end

    return lo
end

function downsample( image, step )
    h, w = size(image)
    rowindices = 1:step:h
    colindices = 1:step:w

    subsampled_img = zeros( eltype(image), length(rowindices), length(colindices) )

    ccont = 1
    for col in colindices
        rcont = 1
        for row in rowindices

            subsampled_img[rcont,ccont] = image[ row, col ]
            rcont += 1
        end
        ccont += 1
    end

    return subsampled_img
end

function meanDownsample( image::Array{}, scale )

	h, w   = size( image )
	nh, nw = length( scale:scale:h ), length( scale:scale:w )
	N      = scale*scale
	imgd   = zeros( eltype(image), nh, nw )

	ccont = 1;
	for col in scale:scale:w
		rcont = 1;
		for row in scale:scale:h
			mean = sum( image[ row-scale+1:row, col-scale+1:col ] )/N
			imgd[ rcont, ccont ] = mean;

		rcont += 1;
		end
	ccont += 1;
	end

	return imgd
end

function reverse!( image::AbstractArray{T,2} ) where {T<:Real }

	h, w = size( image );
	halfh, halfw = div.( size( image ), 2 );

	for col in 1:halfw
		t = image[ :, col ]
		image[ :, col ] .= image[ :, w-col+1 ]
		image[ :, w-col+1 ] .= t
	end

	for row in 1:halfh
		t = image[ row, : ]
		image[ row, : ] .= image[ h-row+1, : ]
		image[  h-row+1, : ] .= t
	end
end


text:
module IntensityScans

using TcpInstruments
using ThorlabsLTStage
using HydrophoneCalibrations
using IntensityMetrics

using Base: @kwdef
using RecipesBase
using Dates
using ProgressMeter
using Unitful
using Unitful: mm, m, μs, ms, s
using UnitfulRecipes

export IntensityScan
export scan_x, scan_y, scan_z
export scan_xy, scan_xz, scan_yz
export scan_xyz
export save, load
export mm, m, μs, ms, s

# Metrics
export ScanParameters
export ScanMetric
export compute_metrics, get_metrics
export intensity, intensity_sppa, intensity_spta, mechanical_index

# IntensityMetrics
export Medium, Excitation

include("util.jl")

include("models.jl")
include("plots.jl")

include("single_axis_scan.jl")
include("double_axis_scan.jl")
include("scan_xyz.jl")

include("metrics/metrics.jl")

end


text:
module YahooFinanceAPI

#using Dates

# QuoteProperties submodule
# Used for the Quotes API
# see online documentation of Yahoo API for full list
# https://code.google.com/p/yahoo-finance-managed/wiki/enumQuoteProperty
module QuoteProperties
include("yahoo_finance_api/QuoteProperties.jl");
end

# MarketQuoteProperties submodule
# used for Industry and Sector APIs
# see online documentation of Yahoo API for full list
# https://code.google.com/p/yahoo-finance-managed/wiki/enumMarketQuoteProperties
module MarketQuoteProperties
include("yahoo_finance_api/MarketQuoteProperties.jl");
end

# IndustryCodes submodule
# Used for using the "Companies" API
# See full list source at
# http://pastie.org/6419646#205
module IndustryCodes
include("yahoo_finance_api/IndustryCodes.jl")
end

# Sectors submodule
# Used Industries API. Must select the sector number
module Sectors
include("yahoo_finance_api/Sectors.jl")
end

# Localization submodule
# used to localize RSS feed
module Localization
include("yahoo_finance_api/Localization.jl")
end

#include library files
include("yahoo_finance_api/types.jl")
include("yahoo_finance_api/apiCalls.jl")

# export
export HistoricalData, fetchHistoricalData
export QuoteProperties, fetchQuotes
export IndustryCodes, MarketQuoteProperties
export fetchSectors
export Sectors, fetchIndustries, fetchCompanies
export Localization, fetchRSS

end


text:
"""
    TimeVaryingCLF

Time varying Control Lyapunov function V for a control affine system
# Fields
- `V`: function V(x, t) that represents the CLF
- `∇V`: function ∇V(x, t) defining the spatial gradient of the CLF
- `∂ₜV`: function ∂ₜV(x, t) defining the time derivative of the CLF
- `α`: class K function α(V(x, t))
"""
struct TimeVaryingCLF <: LyapunovFunction
    V
    ∇V
    ∂ₜV
    α
end

"""
    TimeVaryingCLF(V, α)

Construct a CLF given V and class K function α.
"""
function TimeVaryingCLF(V, α)
    ∇V(x, t) = gradient((r, s) -> V(r, s), x, t)[1]'
    ∂ₜV(x, t) = gradient((r, s) -> V(r, s), x, t)[2]
    return TimeVaryingCLF(V, ∇V, ∂ₜV, α)
end

"""
	(V::TimeVaryingCLF)(x, t)

Evaluate CLF at state x and time t.
"""
function (V::TimeVaryingCLF)(x, t)
    return V.V(x, t)
end


text:
using ConicNonlinearBridge
using Convex
using Ipopt
using FactCheck

tests = ["test_utilities.jl",
         "test_affine.jl",
         "test_lp.jl"]
tests_socp = ["test_socp.jl"]
tests_exp = ["test_exp.jl"]

println("Running tests:")

# The following syntax can be used to solve it using other solvers
set_default_solver(ConicNLPWrapper(nlp_solver=IpoptSolver(print_level=0)))

for curtest in tests
    info(" Test: $(curtest)")
    include(Pkg.dir("Convex","test",curtest))
end

for curtest in tests_socp
    info(" Test: $(curtest)")
    include(Pkg.dir("Convex","test",curtest))
end

for curtest in tests_exp
    info(" Test: $(curtest)")
    include(Pkg.dir("Convex","test",curtest))
end


text:
@doc """
    pairwise(dist::StringDistance, xs::AbstractVector, ys::AbstractVector = xs; preprocess = nothing)

Compute distances between all pairs of elements in `xs`  and `ys` according to the
`StringDistance` `dist`. Returns a matrix R such that `R[i, j]` corrresponds to the distance between `xs[i]` and `ys[j]`.

For AbstractQGramDistances preprocessing will be used either if `preprocess` is set
to true or if there are more than 5 elements in `xs`. Set `preprocess` to
false if no preprocessing should be used, regardless of length.

Both symmetric and asymmetric versions are available.

### Examples
```julia-repl
julia> using StringDistances
julia> iter = ["New York", "Princeton"]
julia> pairwise(Levenshtein(), iter)
2×2 Array{Float64,2}:
 0.0  9.0
 9.0  0.0
julia> iter2 = ["San Francisco"]
julia> pairwise(Levenshtein(), iter, iter2)
2×1 Array{Float64,2}:
 12.0
 10.0
```
"""
function pairwise(dist::StringDistance, xs::AbstractVector, ys::AbstractVector = xs; preprocess = nothing)
    T = result_type(dist, eltype(xs), eltype(ys))
    if Missing <: Union{eltype(xs), eltype(ys)}
        T = Union{T, Missing}
    end
    R = Matrix{T}(undef, length(xs), length(ys))
    pairwise!(R, dist, xs, ys; preprocess = preprocess)
end

@doc """
    pairwise!(R::AbstractMatrix, dist::StringDistance, xs::AbstractVector, ys::AbstractVector = xs; preprocess = nothing)

Compute distances between all pairs of elements in `xs` and `ys` according to the
`StringDistance` `dist` and write the result in `R`. `R[i, j]` corresponds to the distance between `xs[i]` and `ys[j]`.

For AbstractQGramDistances preprocessing will be used either if `preprocess` is set
to true or if there are more than 5 elements in `xs`. Set `preprocess` to
false if no preprocessing should be used, regardless of length.
"""
function pairwise!(R::AbstractMatrix, dist::StringDistance, xs::AbstractVector, ys::AbstractVector = xs; preprocess = nothing)
    length(xs) == size(R, 1) || throw(DimensionMismatch("inconsistent length"))
    length(ys) == size(R, 2) || throw(DimensionMismatch("inconsistent length"))
    ((xs === ys) & (dist isa SemiMetric)) ?
        _symmetric_pairwise!(R, dist, xs; preprocess = preprocess) :
        _asymmetric_pairwise!(R, dist, xs, ys; preprocess = preprocess)
end

function _symmetric_pairwise!(R::AbstractMatrix, dist::StringDistance, xs::AbstractVector; preprocess = nothing)
    objs = _preprocess(xs, dist, preprocess)
    for i in 1:length(objs)
        # handle missing
        R[i, i] = objs[i] != objs[i]
        Threads.@threads for j in (i+1):length(objs)
            R[i, j] = R[j, i] = evaluate(dist, objs[i], objs[j])
        end
    end
    return R
end

function _asymmetric_pairwise!(R::AbstractMatrix, dist::StringDistance, xs::AbstractVector, ys::AbstractVector; preprocess = nothing)
    objsxs = _preprocess(xs, dist, preprocess)
    objsys = xs === ys ? objsxs : _preprocess(ys, dist, preprocess)
    for i in 1:length(objsxs)
        Threads.@threads for j in 1:length(objsys)
            R[i, j] = evaluate(dist, objsxs[i], objsys[j])
        end
    end
    return R
end

function _preprocess(xs, dist::StringDistance, preprocess)
    if preprocess === nothing
        preprocess = length(xs) >= 5
    end
    if (dist isa AbstractQGramDistance) && preprocess
        return fetch.(map(x -> (Threads.@spawn x === missing ? x : QGramSortedVector(x, dist.q)), xs))
    else
        return xs
    end
end


text:
# defines variables found in INDEXP.INC for use in Julia

const KPVACX =  1
const KPVACY =  2
const KPVACZ =  3
const KPRACX =  4
const KPRACY =  5
const KPRACZ =  6
const KPVELX =  7
const KPVELY =  8
const KPVELZ =  9
const KPROTX = 10
const KPROTY = 11
const KPROTZ = 12
const KPPOSX = 13
const KPPOSY = 14
const KPPOSZ = 15
const KPBANK = 16
const KPELEV = 17
const KPHEAD = 18
const KPFLP1 = 19
const KPFLP2 = 20
const KPFLP3 = 21
const KPFLP4 = 22
const KPFLP5 = 23
const KPFLP6 = 24
const KPFLP7 = 25
const KPFLP8 = 26
const KPFLP9 = 27
const KPFLP10= 28
const KPFLP11= 29
const KPFLP12= 30
const KPFLP13= 31
const KPFLP14= 32
const KPFLP15= 33
const KPFLP16= 34
const KPFLP17= 35
const KPFLP18= 36
const KPFLP19= 37
const KPFLP20= 38
const KPENG1 = 39
const KPENG2 = 40
const KPENG3 = 41
const KPENG4 = 42
const KPENG5 = 43
const KPENG6 = 44
const KPENG7 = 45
const KPENG8 = 46
const KPENG9 = 47
const KPENG10= 48
const KPENG11= 49
const KPENG12= 50
const KPVIDT = 51
const KPBEDT = 52
const KPALDT = 53
const KPBADT = 54
const KPELDT = 55
const KPHEDT = 56
const KPWDTX = 57
const KPWDTY = 58
const KPWDTZ = 59
const KPADTX = 60
const KPADTY = 61
const KPADTZ = 62
const KPRCGX = 63
const KPRCGY = 64
const KPRCGZ = 65
const KPVINF = 66
const KPBETA = 67
const KPALFA = 68
const KPVIAS = 69
const KPBREF = 70
const KPAREF = 71
const KPDRAG = 72
const KPSIDE = 73
const KPLIFT = 74
const KPCDIN = 75
const KPCDRG = 76
const KPCSID = 77
const KPCLIF = 78
const KPCROL = 79
const KPCMOM = 80
const KPCYAW = 81
const KPSEFF = 82
const KPAFOX = 83
const KPAFOY = 84
const KPAFOZ = 85
const KPAMOX = 86
const KPAMOY = 87
const KPAMOZ = 88
const KPRFOX = 89
const KPRFOY = 90
const KPRFOZ = 91
const KPRMOX = 92
const KPRMOY = 93
const KPRMOZ = 94
const KPGFOX = 95
const KPGFOY = 96
const KPGFOZ = 97
const KPGMOX = 98
const KPGMOY = 99
const KPGMOZ =100
const KPEFOX =101
const KPEFOY =102
const KPEFOZ =103
const KPEMOX =104
const KPEMOY =105
const KPEMOZ =106
const KPTFOX =107
const KPTFOY =108
const KPTFOZ =109
const KPTMOX =110
const KPTMOY =111
const KPTMOZ =112
const KPMACH =113
const KPALTK =114
const KPDENS =115
const KPVSOU =116
const KPVISC =117
const KPGRAV =118
const KPLOAD =119
const KPTOT  =119
const KPFREE = KPADTZ
const KPPLOT = KPSEFF

const IPVACX =  1
const IPVACY =  2
const IPVACZ =  3
const IPRACX =  4
const IPRACY =  5
const IPRACZ =  6
const IPVIAS =  7
const IPBETA =  8
const IPALFA =  9
const IPROTX = 10
const IPROTY = 11
const IPROTZ = 12
const IPPOSX = 13
const IPPOSY = 14
const IPPOSZ = 15
const IPBANK = 16
const IPELEV = 17
const IPHEAD = 18
const IPFLP1 = 19
const IPFLP2 = 20
const IPFLP3 = 21
const IPFLP4 = 22
const IPFLP5 = 23
const IPFLP6 = 24
const IPFLP7 = 25
const IPFLP8 = 26
const IPFLP9 = 27
const IPFLP10= 28
const IPFLP11= 29
const IPFLP12= 30
const IPFLP13= 31
const IPFLP14= 32
const IPFLP15= 33
const IPFLP16= 34
const IPFLP17= 35
const IPFLP18= 36
const IPFLP19= 37
const IPFLP20= 38
const IPENG1 = 39
const IPENG2 = 40
const IPENG3 = 41
const IPENG4 = 42
const IPENG5 = 43
const IPENG6 = 44
const IPENG7 = 45
const IPENG8 = 46
const IPENG9 = 47
const IPENG10= 48
const IPENG11= 49
const IPENG12= 50
const IPFORX = 51
const IPFORY = 52
const IPFORZ = 53
const IPMOMX = 54
const IPMOMY = 55
const IPMOMZ = 56
const IPLIFT = 57
const IPGAMM = 58
const IPRACC = 59
const IPUSR1 = 60
const IPUSR2 = 61
const IPLSQD = 62
const IPTOT  = 62

const IQEPS = 19
const IQTAU = 20
const IQCH  = 21
const IQCL  = 22
const IQCM  = 23
const IQFL  = 24
const IQFR  = 25
const IQDA  = 26
const IQAF  = 27
const IQTOT = 27

const KSVACX =  1
const KSVACY =  2
const KSVACZ =  3
const KSRACX =  4
const KSRACY =  5
const KSRACZ =  6
const KSVELX =  7
const KSVELY =  8
const KSVELZ =  9
const KSROTX = 10
const KSROTY = 11
const KSROTZ = 12
const KSPOSX = 13
const KSPOSY = 14
const KSPOSZ = 15
const KSBANK = 16
const KSELEV = 17
const KSHEAD = 18
const KSVINF = 19
const KSBETA = 20
const KSALFA = 21
const KSFORX = 22
const KSFORY = 23
const KSFORZ = 24
const KSMOMX = 25
const KSMOMY = 26
const KSMOMZ = 27
const KSXYZX = 28
const KSXYZY = 29
const KSXYZZ = 30
const KSUVWX = 31
const KSUVWY = 32
const KSUVWZ = 33
const KSWDTX = 34
const KSWDTY = 35
const KSWDTZ = 36
const KSADTX = 37
const KSADTY = 38
const KSADTZ = 39
const KSVIDT = 40
const KSBEDT = 41
const KSALDT = 42
const KSBADT = 43
const KSELDT = 44
const KSHEDT = 45
const KSGAMA = 46
const KSEPSB = 47
const KSTAUT = 48
const KSTOT  = 48
const KSPLOT = KSTOT

const IRVACX =  1
const IRVACY =  2
const IRVACZ =  3
const IRRACX =  4
const IRRACY =  5
const IRRACZ =  6
const IRVIAS =  7
const IRBETA =  8
const IRALFA =  9
const IRROTX = 10
const IRROTY = 11
const IRROTZ = 12
const IRPOSX = 13
const IRPOSY = 14
const IRPOSZ = 15
const IRVELX = 16
const IRVELY = 17
const IRVELZ = 18
const IRBANK = 19
const IRELEV = 20
const IRHEAD = 21
const IRFORX = 22
const IRFORY = 23
const IRFORZ = 24
const IRMOMX = 25
const IRMOMY = 26
const IRMOMZ = 27
const IRGAMA = 28
const IREPSB = 29
const IRTAUT = 30
const IRTOT  = 30


text:
# Copyright (c) 2015, Intel Corporation
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
#     * Redistributions of source code must retain the above copyright notice,
#       this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above copyright
#       notice, this list of conditions and the following disclaimer in the
#       documentation and/or other materials provided with the distribution.
#     * Neither the name of Intel Corporation nor the names of its contributors
#       may be used to endorse or promote products derived from this software
#       without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.

using FactCheck
using Latte

net = Net(8)
data, data_value = MemoryDataLayer(net, :data, (10,))
rand!(data_value)
fc1 = InnerProductLayer(:fc1, net, data, 20)
_reshape = ReshapeLayer(:reshape, net, fc1, (10, 2))

init(net)
facts("Testing Reshape Layer") do
    @fact size(get_buffer(net, :reshapevalue)) --> (10, 2, 8)
    @fact size(get_buffer(net, :reshape∇)) --> (10, 2, 8)

    forward(net; phase=Latte.Test)
    @fact get_buffer(net, :reshapevalue)[:] --> get_buffer(net, :fc1value)[:]

    rand!(get_buffer(net, :reshape∇))
    backward(net)
    @fact get_buffer(net, :reshape∇)[:] --> get_buffer(net, :fc1∇)[:]
end

FactCheck.exitstatus()


text:
@testset "approximators" begin
    include("v_approximators.jl")
    include("q_approximators/q_approximators.jl")
end


text:
function test_shared_parameters_layers(backend::Backend, layer_type, T, eps)
  println("-- Testing $layer_type layer with shared param on $(typeof(backend)){$T}...")
  registry_reset(backend)

  w,h,c,n = 2,3,4,5
  input = rand(T, w,h,c,n)

  layer_data = MemoryDataLayer(tops=[:data], batch_size=n, data=Array[input])
  layer_split = SplitLayer(tops=[:data1,:data2], bottoms=[:data])

  key="key"
  if layer_type == "convolution"
    l1 = ConvolutionLayer(name="conv1", param_key=key, neuron=Neurons.Sigmoid(), tops=[:out1], bottoms=[:data1])
    l2 = ConvolutionLayer(name="conv2", param_key=key, neuron=Neurons.Sigmoid(), tops=[:out2], bottoms=[:data2])
  elseif layer_type == "inner-product"
    out_dim = 5
    l1 = InnerProductLayer(name="ip1", output_dim=out_dim, neuron=Neurons.ReLU(), param_key=key, tops=[:out1], bottoms=[:data1])
    l2 = InnerProductLayer(name="ip2", output_dim=out_dim, neuron=Neurons.ReLU(), param_key=key, tops=[:out2], bottoms=[:data2])
  else
    error("Unknown layer_type $layer_type")
  end

  layer_sub = ElementWiseLayer(operation=ElementWiseFunctors.Subtract(),bottoms=[:out1,:out2],tops=[:diff])

  net = Net("test-shared-params", backend, [layer_data, layer_split, l1, l2, layer_sub])
  init(net)
  forward(net)

  @test net.layers[end] == layer_sub
  output = zeros(T, size(net.states[end].blobs[1]))
  copy!(output, net.states[end].blobs[1])
  @test all(abs(output) .< eps)

  destroy(net)
end

function test_shared_parameters_layers(backend::Backend, T, eps)
  test_shared_parameters_layers(backend, "convolution", T, eps)
  test_shared_parameters_layers(backend, "inner-product", T, eps)
end
function test_shared_parameters_layers(backend::Backend)
  test_shared_parameters_layers(backend, Float64, 1e-9)
  test_shared_parameters_layers(backend, Float32, 1e-4)
end

if test_cpu
  test_shared_parameters_layers(backend_cpu)
end
if test_gpu
  test_shared_parameters_layers(backend_gpu)
end


text:
module DifferentiableObjects

using   ForwardDiff,
        DiffResults,
        # StaticArrays, # Explicitly support both
        NaNMath,
        SIMDPirates,
        PaddedMatrices,
        Parameters,
        LinearAlgebra,
        VectorizationBase

using SIMDPirates: vadd, vmul, evmul, vbroadcast, Vec
using PaddedMatrices:   AbstractFixedSizeVector, AbstractFixedSizeMatrix,
                        AbstractMutableFixedSizeVector, AbstractMutableFixedSizeMatrix,
                        ConstantFixedSizeVector, ConstantFixedSizeMatrix,
                        FixedSizeVector, FixedSizeMatrix,
                        AbstractFixedSizeArray, PtrVector, PtrMatrix

# import  NLSolversBase: AbstractObjective,
#         value!!,
#         value!,
#         value,
#         gradient,
#         gradient!,
#         value_gradient!,
#         value_gradient!!,
#         hessian!,
#         hessian!!

export  AbstractDifferentiableObject,
        GradientConfiguration,
        HessianConfiguration,
        OnceDifferentiable,
        TwiceDifferentiable,
        optimize_light!,
        optimize_scale!,
        value,
        value!,
        gradient,
        gradient!,
        hessian,
        hessian!,
        value_gradient!

abstract type AbstractDifferentiableObject{P,T} end #<: AbstractObjective end

# const SizedVector{P,T} = Union{AbstractFixedSizeVector{P,T},MVector{P,T}}
# const SizedSquareMatrix{P,T} = Union{AbstractFixedSizeMatrix{P,P,T},MMatrix{P,P,T}}
const SizedVector{P,T} = AbstractFixedSizeVector{P,T}
const SizedSquareMatrix{P,T} = AbstractFixedSizeMatrix{P,P,T}
const SymmetricMatrix{P,T} = Symmetric{T, <: AbstractFixedSizeMatrix{P,P,T}}
# debug() = true

###
### This replaces the default fallback, which uses unsafe_load.
###
@inline function VectorizationBase.load(ptr::Ptr{D}) where {T, V, N, D <: ForwardDiff.Dual{T,V,N}}
    vptr = Base.unsafe_convert(Ptr{V}, ptr)
    D(
        VectorizationBase.load(vptr),
        ForwardDiff.Partials{N,V}(
            ntuple(n -> VectorizationBase.load(vptr + n*sizeof(V)), Val(N))
        )
    )
end
### This replaces the default fallback, which uses unsafe_store!
@inline function VectorizationBase.store!(ptr::Ptr{D}, v::D) where {T,V,N, D <: ForwardDiff.Dual{T,V,N}}
    vptr = Base.unsafe_convert(Ptr{V}, ptr)
    val = v.value
    partials = v.partials
    VectorizationBase.store!(vptr, val)
    # Compiler should unroll
    @inbounds @simd for n ∈ 1:N
        VectorizationBase.store!(vptr + n*sizeof(V), partials[n])
    end
    v
end

include("initial_linesearch_guess.jl")
include("backtracking.jl")

include("forward_diff_differentiable.jl")
include("scaled_gradients.jl")
include("bfgs.jl")
# include("optimize.jl")
include("optimize_light.jl")
# include("miscellaneous.jl")

end # module