text:
# reduce along dimensions

import ArrayViews: offset
offset(x::Number) = 0

## auxiliary functions

reducindicators(siz::NTuple{1,Int}, d::Integer) = (d == 1 ? (true,) : (false,))

reducindicators(siz::NTuple{2,Int}, d::Integer) = (d == 1 ? (true, false) :
                                                   d == 2 ? (false, true) : (false, false))

reducindicators(siz::NTuple{3,Int}, d::Integer) = (d == 1 ? (true, false, false) :
                                                   d == 2 ? (false, true, false) :
                                                   d == 3 ? (false, false, true) : (false, false, false))

reducindicators{N}(siz::NTuple{N,Int}, d::Integer) = (ds = fill(false,N); ds[d]=true; tuple(ds...))::NTuple{N,Bool}

function reducindicators{N}(siz::NTuple{N,Int}, dims::Union(Dims,Vector))
    ds = fill(false,N)
    for d in dims
        ds[d] = true
    end
    return tuple(ds...)::NTuple{N,Bool}
end

#################################################
#
#   main macro
#
#################################################

macro compose_reducedim(fun, AN)
    # fun: the function name, e.g. sum
    # BT: the base type, e.g. Number, Real

    fun! = symbol(string(fun, '!'))
    _fun! = symbol(string('_', fun!))
    _funimpl! = symbol(string('_', fun, "impl!"))
    _funimpl1d! = symbol(string('_', fun, "impl1d!"))
    _funimpl2d! = symbol(string('_', fun, "impl2d!"))

    # code-gen preparation

    h = codegen_helper(AN)
    he = codegen_helper_ex(AN)
    args = h.args
    aparams = h.dense_aparams

    if AN == 0
        interf_params = [:(op::Functor{2}), aparams...]
    else
        interf_params = [aparams[1], :(op::Functor{2}), aparams[2:end]...]
    end

    saccumf = AN >= 0 ? :saccum : :saccum_fdiff
    paccumf! = AN >= 0 ? :paccum! : :paccum_fdiff!

    quote
        global $(_funimpl1d!)
        function $(_funimpl1d!){Op<:Functor{2}}(op::Op, r1::Bool, n::Int,
                                                dst::ContiguousArray, idst::Int, sdst1::Int,
                                                $(he.imparams1d...))
            if r1
                if $(he.contcol)
                    dst[1] = evaluate(op, dst[1], $(saccumf)(op, n, $(he.kerargs...)))
                else
                    dst[1] = evaluate(op, dst[1], $(saccumf)(op, n, $(he.kerargs1...)))
                end
            else
                if $(he.contcol) && (sdst1 == 1)
                    $(paccumf!)(op, n, dst, idst, $(he.kerargs...))
                else
                    $(paccumf!)(op, n, dst, idst, sdst1, $(he.kerargs1...))
                end
            end
        end

        global $(_funimpl2d!)
        function $(_funimpl2d!){Op<:Functor{2}}(op::Op, r2::Bool, r1::Bool, m::Int, n::Int,
                                                dst::ContiguousArray, idst::Int, sdst1::Int, sdst2::Int,
                                                $(he.imparams2d...))
            if r1
                if r2
                    if $(he.contcol)
                        s = $(saccumf)(op, m, $(he.kerargs...))
                        $(he.nextcol)
                        for j = 2:n
                            s = evaluate(op, s, $(saccumf)(op, m, $(he.kerargs...)))
                            $(he.nextcol)
                        end
                        dst[idst] = evaluate(op, dst[idst], s)
                    else
                        s = $(saccumf)(op, m, $(he.kerargs1...))
                        $(he.nextcol)
                        for j = 2:n
                            s = evaluate(op, s, $(saccumf)(op, m, $(he.kerargs1...)))
                            $(he.nextcol)
                        end
                        dst[idst] = evaluate(op, dst[idst], s)
                    end
                else
                    if $(he.contcol)
                        for j = 1:n
                            dst[idst] = evaluate(op, dst[idst], $(saccumf)(op, m, $(he.kerargs...)))
                            $(he.nextcol)
                            idst += sdst2
                        end
                    else
                        for j = 1:n
                            dst[idst] = evaluate(op, dst[idst], $(saccumf)(op, m, $(he.kerargs1...)))
                            $(he.nextcol)
                            idst += sdst2
                        end
                    end
                end
            else
                if r2
                    if $(he.contcol) && sdst1 == 1
                        for j = 1:n
                            $(paccumf!)(op, m, dst, idst, $(he.kerargs...))
                            $(he.nextcol)
                        end
                    else
                        for j = 1:n
                            $(paccumf!)(op, m, dst, idst, sdst1, $(he.kerargs1...))
                            $(he.nextcol)
                        end
                    end
                else
                    if $(he.contcol) && sdst1 == 1
                        for j = 1:n
                            $(paccumf!)(op, m, dst, idst, $(he.kerargs...))
                            $(he.nextcol)
                            idst += sdst2
                        end
                    else
                        for j = 1:n
                            $(paccumf!)(op, m, dst, idst, sdst1, $(he.kerargs1...))
                            $(he.nextcol)
                            idst += sdst2
                        end
                    end
                end
            end
        end

        global $(_funimpl!)
        function $(_funimpl!){Op<:Functor{2},N}(op::Op, dst::DenseArray, $(aparams...), dim::Int,
                                                insiz::NTuple{N,Int}, rtup::NTuple{N,Bool})

            if dim == 1
                n = insiz[1]::Int
                $(he.getstrides1)
                sdst1 = stride(dst, 1)::Int
                $(he.getoffsets)
                idst = offset(dst) + 1
                $(_funimpl1d!)(op, rtup[1], n, parent(dst), idst, sdst1, $(he.imargs1d...))

            elseif dim == 2
                m = insiz[1]::Int
                n = insiz[2]::Int
                $(he.getstrides2)
                sdst1, sdst2 = strides(dst)::(Int, Int)
                $(he.getoffsets)
                idst = offset(dst) + 1
                $(_funimpl2d!)(op, rtup[2], rtup[1], m, n, parent(dst), idst, sdst1, sdst2, $(he.imargs2d...))

            else
                n = insiz[dim]::Int
                if rtup[dim]
                    _dst = ellipview(dst, 1)
                    for i = 1:n
                        $(_funimpl!)(op, _dst, $(he.eviewargs...), dim-1, insiz, rtup)
                    end
                else
                    for i = 1:n
                        $(_funimpl!)(op, ellipview(dst, i), $(he.eviewargs...), dim-1, insiz, rtup)
                    end
                end
            end
        end

        global $(_fun!)
        function $(_fun!){S<:Number,N}(dst::DenseArray{S,N}, op::Functor{2}, $(aparams...), insiz::NTuple{N,Int}, dims::DimSpec)
            $(_funimpl!)(op, dst, $(args...), length(insiz), insiz, reducindicators(insiz, dims))
            return dst
        end

        global $(fun!)
        function $(fun!)(dst::ContiguousArray, $(interf_params...), dims::DimSpec)
            insiz = $(h.inputsize)
            rsiz = Base.reduced_dims(insiz, dims)
            prod(rsiz) == length(dst) || throw(DimensionMismatch("Incorrect size of dst."))
            $(_fun!)(contiguous_view(dst, rsiz), op, $(args...), insiz, dims)
            return dst
        end

        global $(fun)
        function $(fun)($(interf_params...), dims::DimSpec)
            insiz = $(h.inputsize)
            dst = fill(reduceinit(op, $(h.termtype)), Base.reduced_dims(insiz,dims))
            $(_fun!)(dst, op, $(args...), insiz, dims)
        end
    end
end

@compose_reducedim reducedim 0
@compose_reducedim mapreducedim 1
@compose_reducedim mapreducedim 2
@compose_reducedim mapreducedim 3
@compose_reducedim mapreducedim_fdiff (-2)

#################################################
#
#   basic functions
#
#################################################

sum(a::BitArray, dims::DimSpec) = invoke(sum, (BitArray, Any), a, dims)

macro compose_reducedim_basicfuns(fname, fdiff, OT)
    fname! = symbol(string(fname, '!'))
    fdiff! = symbol(string(fdiff, '!'))

    quote
        global $(fname!)
        $(fname!)(dst::DenseArray, a::NumericArray, dims::DimSpec) = reducedim!(dst, $(OT)(), a, dims)
        $(fname!)(dst::DenseArray, fun::Functor{1}, a::NumericArray, dims::DimSpec) =
            mapreducedim!(dst, fun, $(OT)(), a, dims)
        $(fname!)(dst::DenseArray, fun::Functor{2}, a::DenseArrOrNum, b::DenseArrOrNum, dims::DimSpec) =
            mapreducedim!(dst, fun, $(OT)(), a, b, dims)
        $(fname!)(dst::DenseArray, fun::Functor{3}, a::DenseArrOrNum, b::DenseArrOrNum, c::DenseArrOrNum, dims::DimSpec) =
            mapreducedim!(dst, fun, $(OT)(), a, b, c, dims)

        global $(fdiff!)
        $(fdiff!)(dst::DenseArray, fun::Functor{1}, a::DenseArrOrNum, b::DenseArrOrNum, dims::DimSpec) =
            mapreducedim_fdiff!(dst, fun, $(OT)(), a, b, dims)

        global $(fname)
        $(fname)(a::NumericArray, dims::DimSpec) = reducedim($(OT)(), a, dims)
        $(fname)(fun::Functor{1}, a::NumericArray, dims::DimSpec) = mapreducedim(fun, $(OT)(), a, dims)
        $(fname)(fun::Functor{2}, a::DenseArrOrNum, b::DenseArrOrNum, dims::DimSpec) =
            mapreducedim(fun, $(OT)(), a, b, dims)
        $(fname)(fun::Functor{3}, a::DenseArrOrNum, b::DenseArrOrNum, c::DenseArrOrNum, dims::DimSpec) =
            mapreducedim(fun, $(OT)(), a, b, c, dims)

        global $(fdiff)
        $(fdiff)(fun::Functor{1}, a::DenseArrOrNum, b::DenseArrOrNum, dims::DimSpec) =
            mapreducedim_fdiff(fun, $(OT)(), a, b, dims)
    end
end

@compose_reducedim_basicfuns sum sumfdiff Add
@compose_reducedim_basicfuns maximum maxfdiff _Max
@compose_reducedim_basicfuns minimum minfdiff _Min
@compose_reducedim_basicfuns nonneg_maximum nonneg_maxfdiff NonnegMax

# avoid the ambiguities with BitArrays

_rlen(siz::Dims, d::Int) = siz[d]
_rlen(siz::Dims, reg::(Int,Int)) = siz[reg[1]] * siz[reg[2]]
function _rlen(siz::Dims, reg::Dims)
    p = 1
    for d in reg
        p *= siz[d]
    end
    return p::Int
end

macro compose_meandim(meanf, sumf, AN)
    sumf! = symbol(string(sumf, '!'))
    meanf! = symbol(string(meanf, '!'))
    h = codegen_helper(AN)
    aparams = h.dense_aparams
    args = h.args

    quote
        global $(meanf!)
        function $(meanf!)(r::NumericArray, $(aparams...), dims::DimSpec)
            siz = $(h.inputsize)
            scale!($(sumf!)(fill!(r, 0.0), $(args...), dims), inv(_rlen(siz, dims)))
        end

        global $(meanf)
        function $(meanf)($(aparams...), dims::DimSpec)
            siz = $(h.inputsize)
            scale!($(sumf)($(args...), dims), inv(_rlen(siz, dims)))
        end
    end
end

@compose_meandim mean sum 0
@compose_meandim mean sum 1
@compose_meandim mean sum 2
@compose_meandim mean sum 3
@compose_meandim meanfdiff sumfdiff (-2)

#################################################
#
#   derived functions
#
#################################################

macro mapreducedim_fun1(fname, accum, F, BT)
    fname! = symbol("$(fname)!")
    accum! = symbol("$(accum)!")

    quote
        global $(fname)
        $(fname){T<:$(BT)}(a::DenseArray{T}, dims::DimSpec) = $(accum)(($F)(), a, dims)

        global $(fname!)
        $(fname!){T<:$(BT)}(dst::DenseArray, a::DenseArray{T}, dims::DimSpec) =
            $(accum!)(dst, ($F)(), a, dims)
    end
end

macro mapreducedim_fun2(fname, accum, F, BT)
    fname! = symbol("$(fname)!")
    accum! = symbol("$(accum)!")

    quote
        global $(fname)
        $(fname){TA<:$(BT),TB<:$(BT)}(a::DenseArrOrNum{TA}, b::DenseArrOrNum{TB}, dims::DimSpec) =
            $(accum)(($F)(), a, b, dims)

        global $(fname!)
        $(fname!){TA<:$(BT),TB<:$(BT)}(dst::DenseArray, a::DenseArrOrNum{TA}, b::DenseArrOrNum{TB}, dims::DimSpec) =
            $(accum!)(dst, ($F)(), a, b, dims)
    end
end

# derived functions

@mapreducedim_fun1 sumabs sum AbsFun Number
@mapreducedim_fun1 meanabs mean AbsFun Number
@mapreducedim_fun1 maxabs nonneg_maximum AbsFun Number
@mapreducedim_fun1 minabs minimum AbsFun Number

@mapreducedim_fun1 sumsq sum Abs2Fun Number
@mapreducedim_fun1 meansq mean Abs2Fun Number

@mapreducedim_fun2 sumabsdiff sumfdiff AbsFun Number
@mapreducedim_fun2 meanabsdiff meanfdiff AbsFun Number
@mapreducedim_fun2 maxabsdiff nonneg_maxfdiff AbsFun Number
@mapreducedim_fun2 minabsdiff minfdiff AbsFun Number

@mapreducedim_fun2 sumsqdiff sumfdiff Abs2Fun Number
@mapreducedim_fun2 meansqdiff meanfdiff Abs2Fun Number

@mapreducedim_fun2 dot sum Multiply Real

@mapreducedim_fun1 sumxlogx sum XlogxFun Real
@mapreducedim_fun2 sumxlogy sum XlogyFun Real

entropy{T<:Real}(a::DenseArray{T}, dims::DimSpec) = negate!(sumxlogx(a, dims))
entropy!{T<:Real}(r::DenseArray{T}, a::ContiguousRealArray, dims::DimSpec) = negate!(sumxlogx!(r, a, dims))


text:
module NeuralNetwork

using LinearAlgebra
using NNlib
using NNlibCUDA
using Zygote
using Zygote: @adjoint

export NetWork, Dense, Denseσ, Conv, Flatten, Dropout, MaxPool, MeanPool, add_layer!, train!, @epochs, @simple_epochs, params, KeepSize

include("utils.jl")
include("activations.jl")
include("progress.jl")
include("Layers/Network.jl")
include("Optimize/train.jl")

end

ext: jl


text:
"""Evaluate `control` at every point of `tlist`.

```julia
values = discretize(control, tlist; via_midpoints=true)
```

discretizes the given `control` to a Vector of values defined on the points of
`tlist`.

If `control` is a function, it will will first be evaluated at the midpoint of
`tlist`, see [`discretize_on_midpoints`](@ref), and then the values on the
midpoints are converted to values on `tlist`. This discretization is more
stable than directly evaluationg the control function at the values of `tlist`,
and ensures that repeated round-trips between [`discretize`](@ref) and
[`discretize_on_midpoints`](@ref) can be done safely, see the note in the
documentation of [`discretize_on_midpoints`](@ref).

The latter can still be achieved by passing `via_midpoints=false`. While such a
direct discretization is suitable e.g. for plotting, but it is unsuitable
for round-trips between [`discretize`](@ref) and
[`discretize_on_midpoints`](@ref)  (constant controls on `tlist` may result in
a zig-zag on the intervals of `tlist`).

If `control` is a vector, it will be returned un-modified if it is of the same
length as `tlist`. Otherwise, `control` must have one less value than `tlist`,
and is assumed to be defined on the midpoins of `tlist`. In that case,
[`discretize`](@ref) acts as the inverse of [`discretize_on_midpoints`](@ref).
See [`discretize_on_midpoints`](@ref) for how control values on `tlist` and
control values on the intervals of `tlist` are related.
"""
function discretize(control::Function, tlist; via_midpoints=true)
    if via_midpoints
        vals_on_midpoints = discretize_on_midpoints(control, tlist)
        return discretize(vals_on_midpoints, tlist)
    else
        return [control(t) for t in tlist]
    end
end

function discretize(control::Vector, tlist)
    if length(control) == length(tlist)
        return control
    elseif length(control) == length(tlist) - 1
        # convert `control` on intervals to values on `tlist`
        # cf. pulse_onto_tlist in Python krotov package
        vals = zeros(eltype(control), length(control) + 1)
        vals[1] = control[1]
        vals[end] = control[end]
        for i in 2:length(vals)-1
            vals[i] = 0.5 * (control[i-1] + control[i])
        end
        return vals
    else
        throw(ArgumentError(
            "control array must be defined on intervals of tlist"
        ))
    end
end

"""Shift time grid values the interval midpoints

```julia
tlist_midpoints = get_tlist_midpoints(tlist)
```

takes a vector `tlist` of length ``n`` and returns a vector of length ``n-1``
containing the midpoint values of each interval. The intervals in `tlist` are
not required to be uniform.
"""
function get_tlist_midpoints(tlist)
    tlist_midpoints = zeros(eltype(tlist), length(tlist) - 1)
    tlist_midpoints[1] = tlist[1]
    tlist_midpoints[end] = tlist[end]
    for i in 2:length(tlist_midpoints) - 1
        dt = tlist[i+1] - tlist[i]
        tlist_midpoints[i] = tlist[i] + 0.5 * dt
    end
    return tlist_midpoints
end

@doc raw"""
Evaluate `control` at the midpoints of `tlist`.

```julia
values = discretize_on_midpoints(control, tlist)
```

discretizes the given `control` to a Vector of values on the midpoints of
`tlist`. Hence, the resulting `values` will contain one less value than
`tlist`.

If `control` is a vector of values defined on `tlist` (i.e., of the same length
as `tlist`), it will be converted to a vector of values on the intervals of
`tlist`. The value for the first and last "midpoint" will remain the original
values at the beginning and end of `tlist`, in order to ensure exact bounary
conditions. For all other midpoints, the value for that midpoint will be
calculated by "un-averaging".

For example, for a `control` and `tlist` of length 5, consider the following
diagram:

~~~
tlist index:       1   2   3   4   5
tlist:             ⋅   ⋅   ⋅   ⋅   ⋅   input values cᵢ (i ∈ 1..5)
                   |̂/ ̄ ̄ ̂\ / ̂\ / ̂ ̄ ̄\|̂
midpoints:         x     x   x     x   output values pᵢ (i ∈ 1..4)
midpoints index:   1     2   3     4
~~~

We will have ``p₁=c₁`` for the first value, ``p₄=c₅`` for the last value. For
all other points, the control values ``cᵢ = \frac{p_{i-1} + p_{i}}{2}`` are the
average of the values on the midpoints. This implies the "un-averaging" for the
midpoint values ``pᵢ = 2 c_{i} - p_{i-1}``.

!!! note

    An arbitrary input `control` array may not be compatible with the above
    averaging formula. In this case, the conversion will be "lossy"
    ([`discretize`](@ref) will not recover the original `control` array; the
    difference should be considered a "discretization error"). However, any
    *further* round-trip conversions between points and intervals are bijective
    and preserve the boundary conditions. In this case, the
    [`discretize_on_midpoints`](@ref) and [`discretize`](@ref) methods are each
    other's inverse. This also implies that for an optimal control procedure,
    it is safe to modify *midpoint* values. Modifying the the values on the
    time grid directly on the other hand may accumulate discretization errors.

If `control` is a vector of one less length than `tlist`, it will be returned
unchanged, under the assumption that the input is already properly discretized.

If `control` is a function, the function will be directly evaluated at the
midpoints marked as `x` in the above diagram..
"""
function discretize_on_midpoints(control::T, tlist) where T<:Function
    return discretize(control, get_tlist_midpoints(tlist); via_midpoints=false)
end

function discretize_on_midpoints(control::Vector, tlist)
    if length(control) == length(tlist) - 1
        return control
    elseif length(control) == length(tlist)
        vals = zeros(eltype(control), length(control) - 1)
        vals[1] = control[1]
        vals[end] = control[end]
        for i in 2:length(vals)-1
            vals[i] = 2 * control[i] - vals[i-1]
        end
        return vals
    else
        throw(ArgumentError(
            "control array must be defined on the points of tlist"
        ))
    end
end

"""Extract a Tuple of controls.

```julia
controls = getcontrols(generator)
```

extracts the controls from a single dynamical generator.

```julia
controls = getcontrols(objectives)
```

extracts the controls from a list of objectives (i.e., from each objective's
`generator`)

In either case, controls that occur multiple times, either in a single
generator, or throughout the different objectives, will occur only once in the
result.

By default, assumes that any `generator` is a nested Tuple, e.g.
`(H0, (H1, ϵ1), (H2, ϵ2), ...)` and extracts (ϵ1, ϵ2)

Each control must be a valid argument for `discretize`.
"""
function getcontrols(generator::Tuple)
    controls = []
    slots_dict = IdDict()  # utilized as Set of controls we've seen
    FuncOrVector = Union{Function, Vector}
    for (i, part) in enumerate(generator)
        if isa(part, Tuple)
            control = part[2]
            control :: FuncOrVector  # assert correct type
            if control in keys(slots_dict)
                # We've seen this control before, so we just record the slot
                # where it is referenced
                push!(slots_dict[control], i)
            else
                push!(controls, control)
                slots_dict[control] = [i]
            end
        end
    end
    return Tuple(controls)
end

function getcontrols(objectives::Vector{T}) where {T<:AbstractControlObjective}
    controls = []
    seen_control = IdDict{Any, Bool}()
    for obj in objectives
        obj_controls = getcontrols(obj.generator)
        for control in obj_controls
            if haskey(seen_control, control)
                # skip: already seen
            else
                push!(controls, control)
                seen_control[control] = true
            end
        end
    end
    return Tuple(controls)
end

"""Replace the controls in `generator` with static values.

```julia
G = evalcontrols(generator, vals_dict)
```

replaces the time-dependent controls in `generator` with the values in
`vals_dict` and returns the static operator `G`.

The `vals_dict` is a dictionary (`IdDict`) mapping controls as returned by
`getcontrols(generator)` to values.
"""
function evalcontrols(generator::Tuple, vals_dict::D) where D<:AbstractDict
    if isa(generator[1], Tuple)
        control = generator[1][2]
        G = vals_dict[control] * generator[1][1]
    else
        G = generator[1]
    end
    for part in generator[2:end]
        if isa(part, Tuple)
            control = part[2]
            G += vals_dict[control] * part[1]
        else
            G += part
        end
    end
    return G
end

"""In-place version of [`evalcontrols`](@ref).

```julia
evalcontrols!(G, generator, vals_dict)
```

acts as [`evalcontrols`](@ref), but modifies `G` in-place.
"""
function evalcontrols!(G, generator::Tuple, vals_dict::D) where D<:AbstractDict
    if isa(generator[1], Tuple)
        control = generator[1][2]
        G .= vals_dict[control] * generator[1][1]
    else
        G .= generator[1]
    end
    for part in generator[2:end]
        if isa(part, Tuple)
            control = part[2]
            G .+= vals_dict[control] * part[1]
        else
            G .+= part
        end
    end
    return G
end

"""Get a vector of the derivatives of `generator` w.r.t. each control.

```julia
getcontrolderivs(generator, controls)
```

return as vector containing the derivative of `generator` with respect to each
control in `controls`. The elements of the vector are either `nothing` if
`generator` does not depend on that particular control, or a function `μ(α)`
that evaluates the derivative for a particular value of the control, see
[`getcontrolderiv`](@ref).
"""
function getcontrolderivs(generator, controls)
    controlderivs = Vector{Union{Function, Nothing}}(nothing, length(controls))
    for (i, control) in enumerate(controls)
        controlderivs[i] = getcontrolderiv(generator, control)
    end
    return controlderivs
end

@doc raw"""
Get the derivative of the generator ``G`` w.r.t. the control ``ϵ(t)``.

```julia
μ  = getcontrolderiv(generator, control)
```

returns `nothing` if the `generator` (Hamiltonian or Liouvillian) does not
depend on `control`, or a function `μ(v)` that evaluates

```math
μ(v) = \left.\frac{∂G}{∂ϵ(t)}\right\vert_{ϵ(t)=v}
```

otherwise. That is, a call `μ(v)` will return the static operator resulting
from evaluating the derivative of the dynamical generator ``G`` with respect to
the control filed ``ϵ(t)`` at a particular point in time where the control
field takes the value ``v``.

Note that for the common case of linear control terms, e.g., ``Ĥ = Ĥ_0 + \sum_l
ϵ_l(t) Ĥ_l``, the derivative ``∂Ĥ/∂ϵ_l(t)`` is simply the control Hamiltonian
``Ĥ_l``. Thus, the resulting function `μ` will simply return ``Ĥ_l``, ignoring
the argument `v`.
"""
function getcontrolderiv(generator::Tuple, control)
    control_generator = nothing
    for part in generator
        if isa(part, Tuple)
            if part[2] === control
                if isnothing(control_generator)
                    control_generator = part[1]
                else
                    control_generator += part[1]
                end
            end
        end
    end
    if isnothing(control_generator)
        return nothing
    else
        return v -> control_generator
    end
end

"""Collect all control parameters from the given objectives.

```julia
get_control_parameters(objectives...; kwargs...)
```

returns a vector of control parameters extracted from the controls in
`objectives`. This first extracts the controls with [`getcontrols`] and then
the "control parameters" from each control. What are the "control parameters"
are depends on the type of the control. For controls that are function ϵ(t),
the control parameters are the values of the control function on intervals of
the time grid. The time grid in the case must be passed as a keyword argument
`tlist`.

The control parameters for multiple controls will be concatenated in the
returned vector.
"""
function get_control_parameters(objectives...; kwargs...)
    controls = getcontrols(objectives)
    parameters = [get_control_parameters(control; kwargs...)
                  for control in controls]
    return vcat(parameters...)
end

"""
```julia
get_control_parameters(func; tlist, on_midpoints=true)
```

returns the control parameters of a control function by discretizing to the
time grid in `tlist`. By default, the discretization is on intervals
(midpoints) of the time grid. With `on_midpoints=false`, the discretization
will be on the actual points of the time grid.
"""
function get_control_parameters(func::Function; tlist, on_midpoints=true,
                                kwargs...)
    if on_midpoints
        return discretize_on_midpoints(func, tlist)
    else
        return discretize(func, tlist)
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/HeatTransfer.jl/blob/master/LICENSE

# # Use of flux boundary condition in 2d heat problem

using FEMBase, Test, HeatTransfer, SparseArrays

X = Dict(
    1 => [0.0,0.0],
    2 => [1.0,0.0])

T = Dict(
    1 => 0.0,
    2 => 0.0)

element = Element(Seg2, [1, 2])
update!(element, "geometry", X)
update!(element, "temperature", 0.0 => T)
update!(element, "heat flux", 2.0)
problem = Problem(PlaneHeat, "test problem", 1)
add_elements!(problem, [element])
assemble!(problem, 0.0)
f = Vector(sparse(problem.assembly.f)[:])
f_expected = [1.0; 1.0]
@test isapprox(f, f_expected)


text:
add2(x) = x + 2
my_f(x) = 2 + 3x
d_my_f(x) = ForwardDiff.derivative(my_f, x)


text:
###############
###############
###############

T = 100                                 # Number of periods in the data
P = 2                                   # Number of lags in the factor equation
L = 2                                   # Number of AR lags in the observation equation
N = 10                                  # Total number of series
K = 1                                   # Total number of factors (only the global factor)

H = [ones(N) ident(N) zeros(N)]
R = zeros(N, N)
A = zeros(N, 2)
Z = zeros(2, 2)

F = zeros(P + N, P + N)
F[1, 1] = 0.5
F[1, P+N] = 0.1
F[P+N, 1] = 1.0
F[2:(P+N-1), 2:(P+N-1)] = 0.25 * ident(N)

Q = zeros(P + N, P + N)
Q[1, 1] = 1.0
for i in 1:N
    Q[1+i, 1+i] = 2.0
end

μ = zeros(P + N)

num_obs = T

# Gather all SS parameters
ssmodel = SSModel(H, A, F, μ, R, Q, Z)

# Simulate common components model in state space form
data_y, data_z, data_β = simulateSSModel(num_obs, ssmodel)

samplerparams = DFMStruct(P, L, 1000, 50)            # Set model priors

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

results = PCA1LevelEstimator(data_y, samplerparams)

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

medians = Any[]
quant33 = Any[]
quant66 = Any[]
stds = Any[]

j = 1
for i in 1:size(results.F)[1]
    push!(stds, std(results.F[i, j, :]))
    push!(quant33, quantile(results.F[i, :], 0.33))
    push!(quant66, quantile(results.F[i, :], 0.66))
    push!(medians, median(results.F[i, :]))
end

plot(data_β[:,1])
plot!(results.means.F[:, j])
#plot(medians)
plot!(quant33)
plot!(quant66)


text:
mutable struct YAxis <: Axis
    label::Label{String}
    ticks::Vector{Line}
    ticklabels::Vector{Label{String}}
    frame::Arrow
end

YAxis( ; label, ticks, ticklabels, frame) = YAxis(label, ticks, ticklabels, frame)


text:
# %% accessing model
@testset "Accessing model" begin
    model = ABM(Agent0; properties = Dict(:a => 2, :b => "test"))
    # add 5 aents
    for i in 1:5
        add_agent!(model)
    end
    @test model.scheduler == fastest
    @test typeof(model.agents) <: Dict
    @test model.space == nothing
    @test model.properties == Dict(:a => 2, :b => "test")
    a = model[1]
    @test a isa Agent0
    @test a.id == 1
    @test model.a == 2
    @test model.b == "test"

    newa = Agent0(6)
    model[6] = newa
    @test model[6] == newa
    @test_throws ArgumentError model[7] = newa

    prop2 = Agent2(1, 0.5)
    model2 = ABM(Agent0; properties = prop2)
    @test model2.weight == 0.5

end

@testset "model access typestability" begin
    prop1 = Dict(:a => 0.5)
    prop2 = Agent2(1, 0.5)
    model1 = ABM(Agent0; properties = prop1)
    model2 = ABM(Agent0; properties = prop2)

    test1(model) = model.a
    test2(model) = model.weight
    @test_nowarn @inferred test1(model1)
    @test_nowarn @inferred test2(model2)
end


text:
struct LegendreWeight{T} <: AbstractJacobiWeight{T} end
LegendreWeight() = LegendreWeight{Float64}()
legendreweight(d::AbstractInterval{T}) where T = LegendreWeight{float(T)}()[affine(d,ChebyshevInterval{T}())]

function getindex(w::LegendreWeight{T}, x::Number) where T
    x ∈ axes(w,1) || throw(BoundsError())
    one(T)
end

getproperty(w::LegendreWeight{T}, ::Symbol) where T = zero(T)

sum(::LegendreWeight{T}) where T = 2one(T)

broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(*), ::LegendreWeight{T}, ::LegendreWeight{V}) where {T,V} =
    LegendreWeight{promote_type(T,V)}()

broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(sqrt), w::LegendreWeight{T}) where T = w

broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(Base.literal_pow), ::Base.RefValue{typeof(^)}, w::LegendreWeight, ::Base.RefValue{Val{k}}) where k = w

singularitiesbroadcast(_, L::LegendreWeight) = L # Assume we stay smooth
singularitiesbroadcast(::typeof(exp), L::LegendreWeight) = L
singularitiesbroadcast(::typeof(Base.literal_pow), ::typeof(^), L::LegendreWeight, ::Val) = L

for op in (:+, :-, :*)
    @eval begin
        singularitiesbroadcast(::typeof($op), A::LegendreWeight, B::LegendreWeight) = LegendreWeight{promote_type(eltype(A), eltype(B))}()
        singularitiesbroadcast(::typeof($op), L::LegendreWeight, ::NoSingularities) = L
        singularitiesbroadcast(::typeof($op), ::NoSingularities, L::LegendreWeight) = L
    end
end
singularitiesbroadcast(::typeof(^), L::LegendreWeight, ::NoSingularities) = L
singularitiesbroadcast(::typeof(/), ::NoSingularities, L::LegendreWeight) = L # can't find roots

singularities(::AbstractJacobi{T}) where T = LegendreWeight{T}()
singularities(::Inclusion{T,<:AbstractInterval}) where T = LegendreWeight{T}()
singularities(d::Inclusion{T,<:Interval}) where T = LegendreWeight{T}()[affine(d,ChebyshevInterval{T}())]
singularities(::AbstractFillLayout, P) = LegendreWeight{eltype(P)}()

struct Legendre{T} <: AbstractJacobi{T} end
Legendre() = Legendre{Float64}()

weighted(P::Legendre) = P
weighted(P::Normalized{<:Any,<:Legendre}) = P
weighted(P::SubQuasiArray{<:Any,2,<:Legendre}) = P
weighted(P::SubQuasiArray{<:Any,2,<:Normalized{<:Any,<:Legendre}}) = P

legendre() = Legendre()
legendre(d::AbstractInterval{T}) where T = Legendre{float(T)}()[affine(d,ChebyshevInterval{T}()), :]

"""
     legendrep(n, z)

computes the `n`-th Legendre polynomial at `z`.
"""
legendrep(n::Integer, z::Number) = Base.unsafe_getindex(Legendre{typeof(z)}(), z, n+1)

summary(io::IO, ::Legendre) = print(io, "Legendre()")

==(::Legendre, ::Legendre) = true

OrthogonalPolynomial(w::LegendreWeight{T}) where {T} = Legendre{T}()
orthogonalityweight(::Legendre{T}) where T = LegendreWeight{T}()

function qr(P::Legendre)
    Q = Normalized(P)
    QuasiQR(Q, Diagonal(Q.scaling))
end

function ldiv(P::Legendre{V}, f::Inclusion{T}) where {T,V}
    axes(P,1) == f || throw(DimensionMismatch())
    TV = promote_type(T,V)
    Vcat(zero(TV), one(TV), Zeros{TV}(∞))
end

ldiv(P::Legendre{V}, f::AbstractQuasiFill{T,1}) where {T,V} = _op_ldiv(P, f)
function ldiv(::Legendre{V}, f::AbstractQuasiVector) where V
    T = ChebyshevT{V}()
    [cheb2leg(paddeddata(T \ f)); zeros(V,∞)]
end

"""
    legendre_massmatrix

computes the massmatrix by first re-expanding in Legendre
"""
function legendre_massmatrix(Ac, B)
    A = parent(Ac)
    P = Legendre{eltype(B)}()
    (P\A)'*massmatrix(P)*(P\B)
end

@simplify *(Ac::QuasiAdjoint{<:Any,<:Legendre}, B::Legendre) = massmatrix(Legendre{promote_type(eltype(Ac), eltype(B))}())

# massmatrix(P) = Weighted(P)'P
massmatrix(P::Legendre{T}) where T = Diagonal(convert(T,2) ./ (2(0:∞) .+ 1))

########
# Jacobi Matrix
########

jacobimatrix(::Legendre{T}) where T =  Tridiagonal((one(T):∞)./(1:2:∞), Zeros{T}(∞), (one(T):∞)./(3:2:∞))

# These return vectors A[k], B[k], C[k] are from DLMF. Cause of MikaelSlevinsky we need an extra entry in C ... for now.
function recurrencecoefficients(::Legendre{T}) where T
    n = zero(T):∞
    ((2n .+ 1) ./ (n .+ 1), Zeros{T}(∞), n ./ (n .+ 1))
end

# explicit special case for normalized Legendre
# todo: do we want these explicit constructors for normalized Legendre?
# function jacobimatrix(::Normalized{<:Any,<:Legendre{T}}) where T
#     b = (one(T):∞) ./sqrt.(4 .*(one(T):∞).^2 .-1)
#     Symmetric(_BandedMatrix(Vcat(zeros(∞)', (b)'), ∞, 1, 0), :L)
# end
# function recurrencecoefficients(::Normalized{<:Any,<:Legendre{T}}) where T
#     n = zero(T):∞
#     nn = one(T):∞
#     ((2n .+ 1) ./ (n .+ 1) ./ sqrt.(1 .-2 ./(3 .+2n)), Zeros{T}(∞), Vcat(zero(T),nn ./ (nn .+ 1) ./ sqrt.(1 .-4 ./(3 .+2nn))))
# end

###
# Splines
###

function \(A::Legendre, B::HeavisideSpline)
    @assert B.points == -1:2:1
    Vcat(1, Zeros(∞,1))
end

###
# sum
###

function _sum(P::Legendre{T}, dims) where T
    @assert dims == 1
    Hcat(convert(T, 2), Zeros{T}(1,∞))
end

_sum(p::SubQuasiArray{T,1,Legendre{T},<:Tuple{Inclusion,Int}}, ::Colon) where T = parentindices(p)[2] == 1 ? convert(T, 2) : zero(T)

###
# dot
###

_dot(::Inclusion{<:Any,<:ChebyshevInterval}, a, b) = _dot(singularities(a), singularities(b), a, b)
function _dot(::LegendreWeight, ::LegendreWeight, a, b)
    P = Legendre{promote_type(eltype(a),eltype(b))}()
    dot(P\a, (massmatrix(P) * (P\b)))
end


text:
using AxisArrays
using AxisKeys
using Base: ReshapedArray
using DataFrames: DataFrame
using Dates
using Documenter: doctest
using FeatureTransforms
using FeatureTransforms: _periodic
using FeatureTransforms: cardinality, OneToOne, OneToMany, ManyToOne, ManyToMany
using FeatureTransforms.TestUtils
using Tables: columntable, isrowtable, istable, rowtable
using Test
using TimeZones

@testset "FeatureTransforms.jl" begin
    # The doctests fail unless run on 64bit julia 1.6.x, due to printing differences
    Sys.WORD_SIZE == 64 && v"1.6" <= VERSION < v"1.7" && doctest(FeatureTransforms)

    include("linear_combination.jl")
    include("log.jl")
    include("one_hot_encoding.jl")
    include("periodic.jl")
    include("power.jl")
    include("scaling.jl")
    include("temporal.jl")
    include("traits.jl")
    include("test_utils.jl")

    include("types/tables.jl")
    include("types/vector.jl")
    include("types/matrix.jl")
    include("types/cube.jl")
    include("types/xarray.jl")
end


text:
using DataFrames, Random, WeakRefStrings
using JDF

@testset "JDF.jl nrow" begin
    df = DataFrame([collect(1:100) for i =1:3000])
    df[!, :int_missing] =  rand([rand(rand([UInt, Int, Float64, Float32, Bool])), missing], nrow(df))
    df[!, :missing] .= missing
    df[!, :strs] = [randstring(8) for i in 1:nrow(df)]
    df[!, :stringarray] = StringVector([randstring(8) for i in 1:nrow(df)])

    df[!, :strs_missing] = [rand([missing, randstring(8)]) for i in 1:nrow(df)]
    df[!, :stringarray_missing] = StringVector([rand([missing, randstring(8)]) for i in 1:nrow(df)])

    savejdf("a.jdf", df)

    dfjdf = jdf"a.jdf"

    @test ncol(jdf"a.jdf") == 3006
    @test nrow(jdf"a.jdf") == 100
    @test size(jdf"a.jdf") == (100, 3006)
    @test size(jdf"a.jdf", 1) == 100
    @test size(jdf"a.jdf", 2) == 3006

    # clean up
    rm("a.jdf", force=true, recursive=true)
end


text:
function pad!(v::Vector{Int},n::Int64)
    push!(v,repeat([0],n-length(v))...)
end

function pad_left_right(v::Vector{Float64},l::Int64, r::Int64; p::T = 0) where {T}
    vcat(repeat([p],l),v,repeat([p],r))
end

moving_average(time_series::Union{Vector{Float64},Vector{Int64}}, n₋::Int64, n₊::Int64) = [mean(time_series[(i-n₋):(i+n₊)]) for i in (1+n₋):(length(time_series)-n₊)]
moving_average(time_series::Union{Vector{Float64},Vector{Int64}},window::Int64) = moving_average(time_series,0,window-1)


text:
module Expectations

using LinearAlgebra: dot
using QuantEcon: MarkovChain

export expected_value

expected_value(i_exo, i_endo_next, mc_exo, value) = dot(value[i_endo_next,:], mc_exo.p[i_exo,:])
expected_value(i_exo, mc_exo, value) = value * mc_exo.p[i_exo,:]

using Test
let
  V = [1 2 3;
       4 5 6;
       7 8 9;
       2 5 8]

  mc = MarkovChain([0.1 0.8 0.1;
                    0.8 0.1 0.1;
                    0.1 0.1 0.8])
  @testset "expected value" begin
    @test expected_value(1, 1, mc, V) == dot(V[1,:], mc.p[1,:]) == 2.0
    @test expected_value(2, 3, mc, V) == dot(V[3,:], mc.p[2,:]) ≈ 7.3
    @test expected_value(3, mc, V)[4] == expected_value(3, 4, mc, V) ≈ 7.1
  end

end

end


text:
function available_result(backend)
	pstat = MOI.get(backend, MOI.PrimalStatus())
	dstat = MOI.get(backend, MOI.DualStatus())
	return pstat != MOI.UNKNOWN_RESULT_STATUS &&
	       pstat != MOI.NO_SOLUTION &&
	       dstat != MOI.UNKNOWN_RESULT_STATUS &&
	       dstat != MOI.NO_SOLUTION
end


text:
"""
_ultailL1(C, A, B)

gives L[Block(1,1)] of a block-tridiagonal Toeplitz operator. Based on Delvaux and Dette 2012.
"""

# need to solve [1 -B*inv(L); 0 1] * [C A B; 0 C L] == [C L 0; 0 C L]
# that is
# A - B*inv(L)*C == L
# In the scalar case this is quadratic
# L^2 - A*L + B*C == 0
# so that
# = L = (A ± sqrt(A^2 - 4B*C))/2
# we choose sign(A) to maximise the magnitude as we know inv(T)[1,1] -> 0, hence
# inv(L) -> 0
_ultailL1(c::Number, a::Number, b::Number) = (a + sign(a)*sqrt(a^2-4b*c))/2

# In the matrix case we write inv(L)*C = R (double check) to make it quadratic
# A - B*R == inv(C)inv(R)
# C*A*R - C*B*R^2 == I
# and do eigen decomposition of R to reduce this to a companion matrix problem.

function _ultailL1(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)
    d = size(A,1)
    λs = filter!(λ -> abs(λ) ≤ 1, eigvals([zeros(d,d) -B; -C -A], [B zeros(d,d); zeros(d,d) -B]))
    @assert length(λs) == d
    V = Matrix{eltype(λs)}(undef, d, d)
    j = 1
    for λ in union(λs)
        c = count(==(λ), λs) # multiplicities
        V[:,j:j+c-1] = svd(A-λ*B - C/λ).V[:,end-c+1:end] # nullspace(A-λ*B - C/λ)
        j += c
    end
    C*(V*Diagonal(inv.(λs))/V)
end

function _ul(::TridiagonalToeplitzLayout, J::AbstractMatrix, ::Val{false}; check::Bool = true)
    C = getindex_value(subdiagonaldata(J))
    A = getindex_value(diagonaldata(J))
    B = getindex_value(supdiagonaldata(J))
    L = _ultailL1(C, A, B)
    U = B/L
    UL(Tridiagonal(Fill(convert(typeof(L),C),∞), Fill(L,∞), Fill(U,∞)), OneToInf(), 0)
end

function _ul(::TridiagonalToeplitzLayout, J::AbstractMatrix, ::Val{true}; check::Bool = true)
    C = getindex_value(subdiagonaldata(J))
    A = getindex_value(diagonaldata(J))
    B = getindex_value(supdiagonaldata(J))
    A^2 ≥ 4B*C || error("Pivotting not implemented")
    ul(J, Val(false))
end

function _ul(::BlockTridiagonalToeplitzLayout, J::AbstractMatrix, ::Val{false}; check::Bool = true)
    C = getindex_value(subdiagonaldata(blocks(J)))
    A = getindex_value(diagonaldata(blocks(J)))
    B = getindex_value(supdiagonaldata(blocks(J)))
    L = _ultailL1(C, A, B)
    U = B/L
    UL(mortar(Tridiagonal(Fill(convert(typeof(L),C),∞), Fill(L,∞), Fill(U,∞))), OneToInf(), 0)
end

_inf_getU(::TridiagonalToeplitzLayout, F::UL) = Bidiagonal(Fill(one(eltype(F)),∞),F.factors.du, :U)
_inf_getL(::TridiagonalToeplitzLayout, F::UL) = Bidiagonal(F.factors.d,F.factors.dl, :L)

function _inf_getU(::BlockTridiagonalToeplitzLayout, F::UL)
    d = size(F.factors.blocks.du[1],1)
    II = convert(eltype(F.factors.blocks.du), I(d))
    mortar(Bidiagonal(Fill(II,∞),F.factors.blocks.du, :U))
end

_inf_getL(::BlockTridiagonalToeplitzLayout, F::UL) = mortar(Bidiagonal(F.factors.blocks.d,F.factors.blocks.dl, :L))

getU(F::UL, ::NTuple{2,InfiniteCardinal{0}}) = _inf_getU(MemoryLayout(F.factors), F)
getL(F::UL, ::NTuple{2,InfiniteCardinal{0}}) = _inf_getL(MemoryLayout(F.factors), F)

getU(F::UL{T,<:Tridiagonal}, ::NTuple{2,InfiniteCardinal{0}}) where T = _inf_getU(MemoryLayout(F.factors), F)
getL(F::UL{T,<:Tridiagonal}, ::NTuple{2,InfiniteCardinal{0}}) where T = _inf_getL(MemoryLayout(F.factors), F)


text:
# This file was generated, do not modify it. # hide
predictions = MLJ.predict(mtm, rows=test)
rms_score = round(rms(predictions, targets[test, 1]), sigdigits=4)

@show rms_score

PyPlot.close_figs() # hide

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using CuYao
CUDA_ON = true


text:
"""
    AbstractApprox

Supertype of types specifying how equality or approximate equality should
be computed. See `Equal`, `Approx`, and `EachApprox`.
"""
abstract type AbstractApprox end

"""
    Equal <: AbstractApprox

Approximate equality test that actually demands strict equality.
"""
struct Equal <: AbstractApprox
end

"""
    EachApprox <: AbstractApprox

Demands that each pair of elements are approximately equal.
`kw` are keyword pairs that are forwarded to `isapprox`.
"""
struct EachApprox{T} <: AbstractApprox
    kw::T
end
EachApprox(; kws...) = EachApprox(kws)

"""
   Approx(; kws...) <: AbstractApprox

Specifies using the legacy `isapprox` interface. For example,
for `AbstractMatrix`, matrix norms are used to test closeness.
`kw` are keyword pairs that are forwarded to `isapprox`.
For example, `Approx(atol=1e-9)`.
"""
struct Approx{T} <: AbstractApprox
    kw::T
end
Approx(; kws...) = Approx(kws)

"""
    isapprox(::Equal, x, y)

Return `true` if `x == y`.
"""
Base.isapprox(::Equal, x, y) = (x == y)

"""
    isapprox(a::Union{EachApprox, Approx}, x, y)

Use the definition of approximate equality specified by `a` to determine
if `x` is approximately `y`.
"""
Base.isapprox(a::Union{EachApprox, Approx}, x, y) = isapprox(x, y; pairs(a.kw)...)

"""
    isapprox(a::EachApprox, A::AbstractArray, B::AbstractArray)

Compute element-wise approximate equality.
"""
function Base.isapprox(a::EachApprox, A::AbstractArray, B::AbstractArray)
    n1 = length(A)
    n2 = length(B)
    if n1 != n2
        return false
    end
    for (x, y) in zip(A, B)
        if ! isapprox(a, x, y)
            return false
        end
    end
    return true
end

"""
    UpToPhase <: AbstractApprox

Demands that each pair of elements are approximately equal up to a phase,
that is a number whose absolute value is one.
`kw` are keyword pairs that are forwarded to `isapprox`.
"""
struct UpToPhase{T} <: AbstractApprox
    kw::T
end
UpToPhase(; kws...) = UpToPhase(kws)

function Base.isapprox(a::UpToPhase, x::Number, y::Number)
    aa = Approx(;a.kw...)
    if isapprox(aa, x, zero(x))
        return isapprox(aa, y, zero(y))
    elseif isapprox(aa, y, zero(y))
        return isapprox(aa, x, zero(x))
    end
    return isunitary(x / y, aa)
end

function Base.isapprox(_app::UpToPhase,  A::AbstractArray, B::AbstractArray)
    n1 = length(A)
    n2 = length(B)
    if n1 != n2
        return false
    end
    app = Approx(;_app.kw...)
    seen_non_zero_flag = false
    z = zero(eltype(A)) # TODO use promotion
    for (a, b) in zip(A, B)
        if iszero(a, app)
            !iszero(b, app) && return false
        elseif iszero(b, app)
            !iszero(a, app) && return false
        else
            if ! seen_non_zero_flag
                z = a/b
                isunitary(z, app) || return false
                seen_non_zero_flag = true
            else
                isapprox(app, a/b, z) || return false
            end
        end
    end
    return true
end


text:
using MendelIHT
using SnpArrays
using DataFrames
using Distributions
using BenchmarkTools
using Random
using LinearAlgebra
using GLM
using DelimitedFiles

#run n repeats of L0_reg using the same X and b, but different y
function repeat(
    repeats  :: Int,
    x        :: SnpArray,
    xbm      :: SnpBitMatrix,
    z        :: AbstractMatrix{Float64},
    true_b   :: Vector{Float64},
    cor_pos  :: Vector{Int}, #correct position of the true model
    d        :: UnionAll,
    l        :: Link,
    debias   :: Bool
)
    n = size(xbm, 1)
    k = size(cor_pos, 1)
    estimated_β = zeros(k, repeats)

    for i in 1:repeats
        #simulate phenotypes (e.g. vector y)
        if d == Normal || d == Poisson || d == Bernoulli
            prob = linkinv.(l, xbm * true_b)
            y = [rand(d(i)) for i in prob]
        elseif d == NegativeBinomial
            nn = 10
            μ = linkinv.(l, xbm * true_b)
            clamp!(μ, -20, 20)
            prob = 1 ./ (1 .+ μ ./ nn)
            y = [rand(d(nn, i)) for i in prob] #number of failtures before nn success occurs
        elseif d == Gamma
            μ = linkinv.(l, xbm * true_b)
            β = 1 ./ μ # here β is the rate parameter for gamma distribution
            y = [rand(d(α, i)) for i in β] # α is the shape parameter for gamma
        end
        y = Float64.(y)

        #compute IHT result for less noisy data
        result = L0_reg(x, xbm, z, y, 1, k, d(), l, debias=debias, init=false, show_info=false)

        #store the correct position in estimated model
        estimated_β[:, i] .= result.beta[cor_pos]
    end

    return estimated_β
end

function run()
    #simulat data with k true predictors, from distribution d and with link l.
    repeats = 100 #how many repeats should I run
    n = 5000
    p = 100000
    d = NegativeBinomial
    l = LogLink()
    debias = false

    # set random seed for reproducibility
    Random.seed!(2019)

    # simulat SNP data
    x, = simulate_random_snparray(n, p, undef)
    xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
    z = ones(n, 1) # intercept

    #construct true model b
    true_b = zeros(p)
    true_b[1:4] = [0.1; 0.25; 0.5; 0.8]
    shuffle!(true_b)
    correct_position = findall(x -> x != 0, true_b)

    #run repeats and save to file
    result = repeat(repeats, x, xbm, z, true_b, correct_position, d, l, debias)
    if debias
        writedlm("./repeats/$d" * "_$repeats", result)
    else
        writedlm("./repeats_nodebias/$d" * "_$repeats", result)
    end
end

run()


text:
## PARSE_DSET
#
# Log:
# 11.05.2021 - Created function [1]
#
# Authors:
# [1] Alexander Smolka (alexander.smolka@tum.de)

"""
    parse_dset(raw::CDataRaw)
    parse_dset(path::String)

!!! note "COMSOL Version"
    5.3.0.223

!!! note "Julia Version"
    1.6.0
"""
function parse_dset(raw::CDataRaw; dim="2Dsym")
    # iniate n, u, t, & v
    n = String[]
    u = String[]
    t = Float64[]
    v = Dict{String, Any}()

    # parse dimension argument (get number of spatial dimensions N)
    N, _ = _parsedim(dim)

    # loop through all field variables
    for i in N+1:length(raw.field)
        # get name, unit and time from field string
        name, unit, time = _readfieldall(raw.field[i])

        # add vector from dictionary to correct matrix in dictionary
        v[name] = haskey(v, name) ? hcat(v[name], raw.v[raw.field[i]]) : raw.v[raw.field[i]]

        # always push the time to the vector
        push!(t, time)

        # name MUST be unique
        if name ∈ n; continue; end

        # if no repetition is happening, push variables to vectors
        push!(n, name); push!(u, unit);
    end

    # unique + sort time vector
    t = sort(unique(t))

    # build p (position)
    p = _parseposition(raw.v, dim)

    # return result as CDataParsedDSet object
    return CDataParsedDSet(raw.path, raw.meta, n, u, t, p, v)
end
parse_dset(path::String; dim="2Dsym") = parse_dset(read_spreadsheet(path); dim=dim)

## INTERNAL FUNCTION ======================================================== ##
function _readfieldall(s::String)
    # assumes no parametrix sweep

    # expected format: `<name> (<unit>) @ t=<time>`
    n, u, _, t = split(s, " ")
    return n, u[2:end-1], parse(Float64, t[3:end])
end

function _parseposition(v::Dict{String, Any}, dim::String)
    if dim == "2Dsym"
        return CPosition2Dsym(v["R"], v["Z"])
    end
end


text:
module AVDevice
  include(joinpath(Pkg.dir("VideoIO"),"src","init.jl"))
  w(f) = joinpath(avdevice_dir, f)

  using AVUtil
  using AVFormat
  using AVCodecs
  include(w("LIBAVDEVICE.jl"))
end


text:
@testset "Chars" begin
    @test FP.EOS == '\0'
    @test all(e -> e in FP.SPACE_CHAR, (' ', '\n', '\t', FP.EOS))
    @test '0' in FP.NUM_CHAR
    @test '5' in FP.NUM_CHAR
    @test '9' in FP.NUM_CHAR
end

@testset "forward_match" begin
    # no post-check
    tf = FP.forward_match("abc")
    @test tf.steps == 2 == length("abc") - 1
    @test tf isa FP.TokenFinder
    @test FP.fixed_lookahead(tf, FP.subs("abc"), false) == (true, 0)
    @test isempty(tf.check.pattern)

    # with post-check
    tf = FP.forward_match("abc", ['α'])
    @test tf.steps == 3 == length("abc") - 1 + 1
    @test FP.fixed_lookahead(tf, FP.subs("abcα"), false) == (true, 1)
    @test FP.fixed_lookahead(tf, FP.subs("abcα"), true) == (false, 0)

    # single char
    tf = FP.forward_match("{")
    @test tf.steps == 0 == length("{") - 1
    @test FP.fixed_lookahead(tf, FP.subs("{"), false) == (true, 0)
    @test FP.fixed_lookahead(tf, FP.subs("{"), true) == (true, 0)
end

@testset "greedy_match" begin
    tf = FP.greedy_match(tail_chars=FP.NUM_CHAR)
    @test tf.steps == -1
    @test FP.greedy_lookahead(tf, 1, '1')
    @test !FP.greedy_lookahead(tf, 1, 'a')
    @test tf isa FP.TokenFinder
    @test isempty(tf.check.pattern)
    # with check
    tf = FP.greedy_match(tail_chars=FP.NUM_CHAR, check=r"^\d{3}$")
    @test FP.greedy_lookahead(tf, 1, '1')
    @test FP.greedy_lookahead(tf, 2, '2')
    @test FP.greedy_lookahead(tf, 3, '3')
    @test FP.check(tf, FP.subs("123"))
    @test !FP.check(tf, FP.subs("1234"))
end

@testset "is_div_open" begin
    tf = FP.F_DIV_OPEN
    @test FP.greedy_lookahead(tf, 1, '@')
    @test !FP.greedy_lookahead(tf, 2, '@')
    for c in ('a', 'Z')
        @test FP.greedy_lookahead(tf, 2, c)
    end
    # not allowed as first char
    for c in ('1', '_', 'α')
        @test !FP.greedy_lookahead(tf, 2, c)
    end
    # allowed as subsequent char
    for c in ('1', '_', 'a', 'Z')
        @test FP.greedy_lookahead(tf, 3, c)
    end
end

@testset "is_latex_command" begin
    tf = FP.F_LX_COMMAND
    @test FP.greedy_lookahead(tf, 1, 'a')
    @test FP.check(tf, FP.subs("\\abcd"))
    @test !FP.check(tf, FP.subs("\\_abcd"))
    @test FP.check(tf, FP.subs("\\abcd1"))
end

@testset "is_lang" begin
    tf = FP.F_LANG_3
    @test FP.greedy_lookahead(tf, 1, '`')
    @test FP.greedy_lookahead(tf, 2, '`')
    @test !FP.greedy_lookahead(tf, 3, '`')
    @test FP.greedy_lookahead(tf, 3, 'j')
    @test FP.check(tf, FP.subs("```julia"))
    @test !FP.check(tf, FP.subs("```013"))

    tf = FP.F_LANG_4
    @test FP.greedy_lookahead(tf, 3, '`')
    @test !FP.greedy_lookahead(tf, 4, '`')
    @test FP.check(tf, FP.subs("````hello"))
    @test !FP.check(tf, FP.subs("```hello"))
    @test !FP.check(tf, FP.subs("`````hello"))

    tf = FP.F_LANG_5
    @test FP.greedy_lookahead(tf, 4, '`')
    @test !FP.greedy_lookahead(tf, 6, '`')
    @test !FP.check(tf, FP.subs("```hello"))
    @test !FP.check(tf, FP.subs("````hello"))
    @test FP.check(tf, FP.subs("`````hello"))
    @test !FP.check(tf, FP.subs("``````hello"))
end

@testset "is_html_entity" begin
    tf = FP.F_HTML_ENTITY
    @test FP.greedy_lookahead(tf, 2, 'a')
    @test !FP.greedy_lookahead(tf, 2, 'α')
    @test FP.check(tf, FP.subs("*"))
end

@testset "is_emoji" begin
    tf = FP.F_EMOJI
    for c in ('a', '+', '_', '-')
        @test FP.greedy_lookahead(tf, 1, c)
    end
    @test FP.check(tf, FP.subs(":smile:"))
end

@testset "is_footnote" begin
    tf = FP.F_FOOTNOTE
    @test FP.check(tf, FP.subs("[^hello]:"))
    @test FP.check(tf, FP.subs("[^h1]"))
end

@testset "regexes" begin
    tf = FP.F_LX_COMMAND
    for c in (FP.subs(raw"\abc"), FP.subs(raw"\abc_def"))
        @test FP.check(tf, c)
    end

    for c in (FP.subs(raw"\abc "), FP.subs(raw"\abc1_def*"), FP.subs(raw"\abc**"),
              FP.subs(raw"\a*bc"), FP.subs(raw"\abc1_"),  FP.subs(raw"\_abc1"))
        @test !FP.check(tf, c)
    end

    tf = FP.F_HTML_ENTITY
    for e in ("⊓", "⊓", "⊓", "⊓",
              "≻", "⅞", "⪌", "⪊")
        @test FP.check(tf, FP.subs(e))
    end
    @test !FP.check(tf, FP.subs("&42;"))
end


text:
function InterSatDirectViewDetector(arg0::OneAxisEllipsoid, arg1::PVCoordinatesProvider)
    return InterSatDirectViewDetector((OneAxisEllipsoid, PVCoordinatesProvider), arg0, arg1)
end

function g(obj::InterSatDirectViewDetector, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

function get_central_body(obj::InterSatDirectViewDetector)
    return jcall(obj, "getCentralBody", OneAxisEllipsoid, ())
end

function get_slave(obj::InterSatDirectViewDetector)
    return jcall(obj, "getSlave", PVCoordinatesProvider, ())
end


text:
using DataStructures: OrderedDict
#TODO: seed, copy

mutable struct DictSpace <: AbstractSpace
    spaces::OrderedDict{Union{Symbol, AbstractString}, AbstractSpace}
    shape::Tuple
    DictSpace(spaces::OrderedDict{<:Union{Symbol, AbstractString}, <:AbstractSpace}) =  new(spaces, ())
    DictSpace(;space_kwargs...) = new(OrderedDict{Symbol, AbstractSpace}(space_kwargs), ())
end

DictSpace(spaces::Dict{<:Union{Symbol, AbstractString}, <:AbstractSpace})  =
    DictSpace(OrderedDict(sort([(sym, space) for (sym, space) in pairs(spaces)])))

sample(dict_obj::DictSpace) = OrderedDict([(k, sample(space)) for (k, space) in pairs(dict_obj.spaces)])

function contains(x, dict_obj::DictSpace)
    # If x is not a dict or OrderedDict or if x doesn't have the same length as spaces
    if !(isa(x, Dict) || isa(x, OrderedDict)) || length(x) != length(dict_obj.spaces)
        return false
    end

    for (k, space) in pairs(dict_obj.spaces)
        # If k is not in x, or if x[k] ∉ space return false
        (isnothing(get(x, k, nothing)) || !(x[k] ∈ space)) && return false
    end
    return true
end

Base.:(==)(dict_obj::DictSpace, other::DictSpace) = dict_obj.spaces == other.spaces


text:
export get_problem_data

function get_problem_data(T::Symbol, b::Real)
    Lˣ, Lʸ, α, c = problem_specific(Val{T}())
    N = length(α)
    β = b*ones(N)
    ubˣ = min(sqrt(α.*β), Lˣ)
    ubʸ = min(sqrt(α.*β), Lʸ)
    lbˣ = α ./ ubˣ
    lbʸ = α ./ ubʸ

    @assert N == length(α) == size(c,1) == size(c,2) == length(β) ==
            length(lbˣ) == length(lbʸ) == length(ubˣ) == length(ubʸ)
    @assert issym(c)

    Problem(N, Lˣ, Lʸ, lbˣ, lbʸ, ubˣ, ubʸ, α, β, c)
end

function problem_specific{T}(::Val{T}) # fallback to the MCNC benchmarks
    T in [:hp, :apte, :xerox, :ami33, :ami49] || error("Unrecognized benchmark $T")
    (W,H), α, nets = yalparser(joinpath(Pkg.dir("FloorLayout"),"data",string(T)*".yal"))
    N = length(α)
    c = zeros(N,N)
    for net in nets
        k = length(net)
        for i in net, j in net
            if i < j
                c[i,j] += 1/(k-1)
                c[j,i] += 1/(k-1)
            end
        end
    end
    return W, H, α, c
end

function problem_specific(::Val{:small_test})
    N = 2
    α = [.25, .3]#, .35]
    c = rand(N,N)
    c += c'
    return 1, 1, α, c
end

function problem_specific(::Val{:hp4})
    (W,H), α, nets = yalparser(joinpath(Pkg.dir("FloorLayout"),"examples","data","hp.yal"))
    N = 4
    c = zeros(N,N)
    for net in nets
        k = length(net)
        for i in net, j in net
            (i > N || j > N) && continue
            if i < j
                c[i,j] += 1/(k-1)
                c[j,i] += 1/(k-1)
            end
        end
    end
    return 100W, 100H, α[1:N], c
end

function problem_specific(::Val{:hp5})
    (W,H), α, nets = yalparser(joinpath(Pkg.dir("FloorLayout"),"examples","data","hp.yal"))
    N = 5
    c = zeros(N,N)
    for net in nets
        k = length(net)
        for i in net, j in net
            (i > N || j > N) && continue
            if i < j
                c[i,j] += 1/(k-1)
                c[j,i] += 1/(k-1)
            end
        end
    end
    return 100W, 100H, α[1:N], c
end

function problem_specific(::Val{:hp7})
    (W,H), α, nets = yalparser(joinpath(Pkg.dir("FloorLayout"),"examples","data","hp.yal"))
    N = 7
    c = zeros(N,N)
    for net in nets
        k = length(net)
        for i in net, j in net
            (i > N || j > N) && continue
            if i < j
                c[i,j] += 1/(k-1)
                c[j,i] += 1/(k-1)
            end
        end
    end
    return W, H, α[1:N], c
end

function problem_specific(::Val{:Armour62_1})
    Lˣ, Lʸ = 22, 33 # pad a bit
    α = [27, 18 ,27, 18, 18, 18, 9, 9, 9, 24, 60, 42, 18, 24, 27, 75, 64, 41, 27, 45]
    #     1  2  3  4  5  6  7  8  9 10 11 12
    c = [   0   1.8   1.2     0    0    0     0     0    0 1.04  1.12     0    0  1.20    0     0     0     0     0     0
          1.8     0  0.96 24.45 0.78    0 13.95     0 1.20 1.35     0     0    0     0    0     0     0     0  6.90     0
          1.2  0.96     0     0    0 2.21     0     0 3.15 3.90     0     0    0 13.05    0     0     0     0 13.65     0
            0 24.45     0     0 1.08 5.70  7.50     0 2.34    0     0  1.40    0     0    0     0     0  1.50 15.75     0
            0  0.78     0  1.08    0    0  2.25  1.35    0 1.56     0     0    0     0 1.35     0     0     0     0     0
            0     0  2.21  5.70    0    0  6.15     0    0    0     0  0.45    0     0    0     0     0  1.05     0     0
            0 13.95     0  7.50 2.25 6.15     0 24.00    0 1.87     0     0    0  0.96    0     0     0  1.65     0  3.75
            0     0     0     0 1.35    0  24.0     0    0    0     0     0 0.60     0    0     0     0     0  7.50 33.45
            0  1.20  3.15  2.34    0    0     0     0    0    0     0     0    0  7.50    0     0  7.50     0     0     0
         1.04  1.35  3.90     0 1.56    0  1.87     0    0    0  0.36  12.0    0  18.6 1.92     0     0     0  5.25     0
         1.12     0     0     0    0    0     0     0    0 0.36     0  2.25    0  3.00 0.96  2.25     0     0     0     0
            0     0     0  1.40    0 0.45     0     0    0 12.0  2.25     0    0     0 1.65     0 15.00     0  8.40     0
            0     0     0     0    0    0     0  0.60    0    0     0     0    0  8.00 1.04  6.00     0     0     0     0
         1.20     0 13.05     0    0    0  0.96     0 7.50 18.6  3.00     0 8.00     0 9.75     0     0  0.90     0     0
            0     0     0     0 1.35    0     0     0    0 1.92  0.96  1.65 1.04  9.75    0     0  5.25     0     0     0
            0     0     0     0    0    0     0     0    0    0  2.25     0 6.00     0    0     0 12.00     0     0     0
            0     0     0     0    0    0     0     0 7.50    0     0  15.0    0     0 5.25 12.00     0     0  7.50     0
            0     0     0  1.50    0 1.05  1.65     0    0    0     0     0    0  0.90    0     0     0     0  4.65     0
            0  6.90 13.65 15.75    0    0     0  7.50    0 5.25     0  8.40    0     0    0     0  7.50  4.65     0     0
            0     0     0     0    0    0  3.75 33.45    0    0     0     0    0     0    0     0     0     0     0     0]
    return Lˣ, Lʸ, α, c
end

function problem_specific(::Val{:Armour62_2})
    Lˣ, Lʸ = 22, 33 # pad a bit
    #       A    B    C    D    E    F    G    H    J    K    L    M    N    P    R    S    T    U    V    W
    α = [  27,  18,  27,  18,  18,  18,   9,   9,   9,  24,  60,  42,  18,  24,  27,  75,  64,  41,  27,  45]
    c = [   0  120   80    0    0    0    0    0    0   40   80    0    0   80    0    0    0    0    0    0  #  A
            0    0   80 1630   30    0  930    0   80   90    0    0    0    0    0    0    0    0  460    0  #  B
            0    0    0    0    0  130    0    0  210  260    0    0    0  870    0    0    0    0  910    0  #  C
            0    0    0    0   60  380  500    0  130    0    0   70    0    0    0    0    0  100 1050    0  #  D
            0    0    0    0    0    0  150   90    0   60    0    0    0    0   90    0    0    0    0    0  #  E
            0    0    0    0    0    0  410    0    0    0    0   30    0    0    0    0    0   70    0    0  #  F
            0    0    0    0    0    0    0 1600    0  110    0    0    0   60    0    0    0  110    0  250  #  G
            0    0    0    0    0    0    0    0    0    0    0    0   40    0    0    0    0    0  500 2230  #  H
            0    0    0    0    0    0    0    0    0    0    0    0    0  500    0    0  500    0    0    0  #  J
            0    0    0    0    0    0    0    0    0    0   30  800    0 1240  160    0    0    0  350    0  #  K
            0    0    0    0    0    0    0    0    0    0    0  150    0  200   80 1500  350   90    0    0  #  L
            0    0    0    0    0    0    0    0    0    0    0    0    0    0  110    0 1000    0  360    0  #  M
            0    0    0    0    0    0    0    0    0    0    0    0    0  500   40  500    0   40    0    0  #  N
            0    0    0    0    0    0    0    0    0    0    0    0    0    0  650    0    0   60    0    0  #  P
            0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0  350    0    0    0  #  R
            0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0 1000    0    0    0  #  S
            0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0  500    0  #  T
            0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0  320    0  #  U
            0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0  #  V
            0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0] #  W
    c = c + c'

    return Lˣ, Lʸ, α, c
end

function problem_specific(::Val{:Bazaraa75_1})
    Lˣ, Lʸ = 10, 6
    α = [  9,  8, 10,  6,  4,  3,  3,  4,  2,   2,  1,  1,  2]
    c = [  0 288 180  54  72 180  27  72  36    0   0   9   0;
           0   0 240  54  72  24  48 160  16   64   8  16   0;
           0   0   0 120  80   0  60 120  60    0   0  30   0;
           0   0   0   0  72  18  18  48  24   48  12   0   0;
           0   0   0   0   0  12  12  64  16   16   4   8   0;
           0   0   0   0   0   0  18  24   6   12   3   3   0;
           0   0   0   0   0   0   0   0   6    6   3   6   0;
           0   0   0   0   0   0   0   0  16   16  16   4   0;
           0   0   0   0   0   0   0   0   0    4   4   2   0;
           0   0   0   0   0   0   0   0   0    0   2   2   0;
           0   0   0   0   0   0   0   0   0    0   0   2   0;
           0   0   0   0   0   0   0   0   0    0   0   0   0;
           0   0   0   0   0   0   0   0   0    0   0   0   0]
    c = c + c'

    return Lˣ, Lʸ, α, c
end

function problem_specific(::Val{:Bazaraa75_2})
    Lˣ, Lʸ = 9, 7

    α = [  9,  8,  9, 10,  6,  3,  3,  3,  2,   3,  2,  1,  1,  1]
    c = [  0  72 162  90 108  27   0   0  18   27  18   0   0   0;
           0   0  72  80   0  48   0  48  32    0  16   8   0   0;
           0   0   0  45  54  27  27  27   0   27   0   9  18   0;
           0   0   0   0  30   0  30  30  20   0   20  10  10   0;
           0   0   0   0   0  18   0  18  12   18  24   0   0   0;
           0   0   0   0   0   0   9   9   0    0   6   6   6   0;
           0   0   0   0   0   0   0   9  12    9   6   3   0   0;
           0   0   0   0   0   0   0   0   6    9   0   3   0   0;
           0   0   0   0   0   0   0   0   0    6   4   6   2   0;
           0   0   0   0   0   0   0   0   0    0   6   3   6   0;
           0   0   0   0   0   0   0   0   0    0   0   2   0   0;
           0   0   0   0   0   0   0   0   0    0   0   0   4   0;
           0   0   0   0   0   0   0   0   0    0   0   0   0   0;
           0   0   0   0   0   0   0   0   0    0   0   0   0   0]
    c = c + c'

    return Lˣ, Lʸ, α, c
end

function problem_specific(::Val{:Camp91})
    Lˣ, Lʸ = 25, 51
    α = [238, 112, 160, 80, 120, 80, 60, 85, 221, 119]
    N = length(α)
    c = zeros(N,N)
    c[ 1, 6] = 218
    c[ 2, 6] = 148
    c[ 2, 9] = 296
    c[ 3, 4] =  28
    c[ 3, 5] =  70
    c[ 4, 6] =  28
    c[ 4, 7] =  70
    c[ 4, 8] = 140
    c[ 5, 8] = 210
    c[ 7,10] =  28
    c[ 8,10] = 888
    c[ 9,10] =  59.2
    c = c + c'

    return Lˣ, Lʸ, α, c
end

function problem_specific(::Val{:Bozer91})
    Lˣ, Lʸ = 10, 20
    α = [12, 7, 6, 5, 7, 22, 22, 13, 7, 22, 9, 13, 4, 17, 25]
    N = length(α)
    c = zeros(N,N)
    c[ 1,15] =  240
    c[ 2, 1] =  240
    c[ 3, 4] = 1200
    c[ 4,10] = 1200
    c[ 5,14] =  600
    c[ 6, 8] =  480
    c[ 7, 8] =  480
    c[ 8,15] =  120
    c[ 9,10] =  600
    c[10,12] =  600
    c[11, 7] =  480
    c[12,15] =  600
    c[13, 7] =  480
    c[14,12] =  600
    c[15, 2] =   10
    c[15, 3] =   25
    c[15, 5] =   25
    c[15, 6] =   40
    c[15, 9] =   25
    c[15,11] =   40
    c[15,13] =   20
    c = c + c'

    return Lˣ, Lʸ, α, c
end

function problem_specific(::Val{:Bozer97_1})
    Lˣ, Lʸ = 12, 13
    α = [16, 16, 16, 36, 36, 9, 9, 9, 9]
    #    1 2 3 4 5 6 7 8 9
    c = [0 0 0 5 5 0 0 0 1
         0 0 0 3 3 0 0 0 1
         0 0 0 2 2 0 0 0 1
         0 0 0 0 0 4 4 0 0
         0 0 0 0 0 3 0 0 4
         0 0 0 0 0 0 0 0 2
         0 0 0 0 0 0 0 0 1
         0 0 0 0 0 0 0 0 0
         0 0 0 0 0 0 0 0 0]
    c = c + c'

    return Lˣ, Lʸ, α, c
end

function problem_specific(::Val{:Bozer97_2})
    Lˣ, Lʸ = 6, 8
    α = [1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 16, 16]
    #     1  2  3  4  5  6  7  8  9 10 11 12
    c = [ 0  2  0  0 10  0  0  0  9  0  0  0 #  1
          0  0  0  5  0  0  0  0  0  0  7  2 #  2
          0  0  0  0  0  2  9  0  0  0  0  0 #  3
          0  0  0  0  0  0  0  0  0  3  0  3 #  4
          0  0  0  0  0  4  0  1  0  0  0  5 #  5
          0  0  0  0  0  0  0  0  4  0  0  0 #  6
          0  0  0  0  0  0  0  0  0  0  0  0 #  7
          0  0  0  0  0  0  0  0  0  0  5  0 #  8
          0  0  0  0  0  0  0  0  0  0  0  3 #  9
          0  0  0  0  0  0  0  0  0  0  0  0 # 10
          0  0  0  0  0  0  0  0  0  0  0  1 # 11
          0  0  0  0  0  0  0  0  0  0  0  0]# 12
    c = c + c'

    return Lˣ, Lʸ, α, c
end

# α = rand(1:N, N)
# Lˣ = ceil(Int, rand(.8:0.01:1.2) * sqrt(sum(α)))
# Lʸ = ceil(Int, sum(α) / Lˣ)
# c = (tmp = sprand(N,N,.2); tmp += tmp'; for i in 1:N; tmp[i,i] = 0; end; full(tmp))
# @show α;
# @show Lˣ, Lʸ;
# @show c;

function problem_specific(::Val{:Rand9})
    Lˣ, Lʸ = 6, 6
    α = [7,4,4,1,1,1,1,8,8]
    c = [0.0 0.7986094862350401 0.8499233254733563 0.0 0.7157553198075688 0.0 0.7437930539668685 0.0 1.175010054208596
 0.7986094862350401 0.0 0.0 0.0 0.0 1.527404963754088 0.0 0.8993465217540104 0.9770788914543318
 0.8499233254733563 0.0 0.0 0.6928155492041943 0.0 0.327312670502387 0.0 0.01939277023677377 0.9952866480280969
 0.0 0.0 0.6928155492041943 0.0 0.11064546027043032 0.2141258056898243 0.998823070162278 0.06454412519826769 0.4958584789801965
 0.7157553198075688 0.0 0.0 0.11064546027043032 0.0 0.0 0.0 0.16233570993620972 0.0014645204713570337
 0.0 1.527404963754088 0.327312670502387 0.2141258056898243 0.0 0.0 0.0 0.0 0.0
 0.7437930539668685 0.0 0.0 0.998823070162278 0.0 0.0 0.0 0.0 0.0
 0.0 0.8993465217540104 0.01939277023677377 0.06454412519826769 0.16233570993620972 0.0 0.0 0.0 0.4233288909381965
 1.175010054208596 0.9770788914543318 0.9952866480280969 0.4958584789801965 0.0014645204713570337 0.0 0.0 0.4233288909381965 0.0]

    return Lˣ, Lʸ, α, c
end

function problem_specific(::Val{:Rand10})
    Lˣ, Lʸ = 8, 6
    α =  [8,4,1,7,3,7,4,9,3,1]
    #     1  2  3  4  5  6  7  8  9 10
    c =  c = [0.0 0.0 0.9049260691773204 0.0 0.0 0.0 0.0 0.9989655886238922 0.0 0.0
 0.0 0.0 1.5694241789894503 0.00354884895732277 0.3137620799856242 0.0 0.0 0.0 0.0 0.2173791447506559
 0.9049260691773204 1.5694241789894503 0.0 0.9382338874211305 0.0 0.0 0.0 0.6139693075117343 0.0 0.0
 0.0 0.00354884895732277 0.9382338874211305 0.0 0.0 0.0 0.8913110827285136 0.0 0.0 0.3414856679085707
 0.0 0.3137620799856242 0.0 0.0 0.0 0.0 0.0 0.0 0.6328948674240107 0.0
 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.7456193867033032 0.1612888230451428 0.0
 0.0 0.0 0.0 0.8913110827285136 0.0 0.0 0.0 0.1620216504491383 0.0 0.11427952217846271
 0.9989655886238922 0.0 0.6139693075117343 0.0 0.0 0.7456193867033032 0.1620216504491383 0.0 0.0 0.17837661199235666
 0.0 0.0 0.0 0.0 0.6328948674240107 0.1612888230451428 0.0 0.0 0.0 0.9729596485553473
 0.0 0.2173791447506559 0.0 0.3414856679085707 0.0 0.0 0.11427952217846271 0.17837661199235666 0.9729596485553473 0.0]
    return Lˣ, Lʸ, α, c
end

function problem_specific(::Val{:Rand11})
    Lˣ, Lʸ = 9, 8
    α = [8,3,2,8,9,9,5,5,1,10,10]
    c = [0.0 0.0 0.6681298698645621 0.5469619335589575 0.0 0.6656772900595389 0.5422023238674989 0.0 0.7837115852348486 0.0 0.6796644127670064
 0.0 0.0 0.0 0.0 0.9431420891034161 0.0 0.0 1.6639710742325924 0.8441121778543434 0.0 0.0
 0.6681298698645621 0.0 0.0 0.0 0.8257319756858643 0.5105926660160405 0.5592266444786895 0.0 0.4075160359832646 0.0 0.0
 0.5469619335589575 0.0 0.0 0.0 0.0 0.0 0.6821101875680364 0.0 0.0 0.676801127801794 0.1920125959201131
 0.0 0.9431420891034161 0.8257319756858643 0.0 0.0 0.0 0.0 0.0 0.2776509135014058 0.0 0.0
 0.6656772900595389 0.0 0.5105926660160405 0.0 0.0 0.0 0.35272294703294826 0.0 0.0 0.0 0.0
 0.5422023238674989 0.0 0.5592266444786895 0.6821101875680364 0.0 0.35272294703294826 0.0 0.0 0.0 0.0 0.026042138065528686
 0.0 1.6639710742325924 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.7837115852348486 0.8441121778543434 0.4075160359832646 0.0 0.2776509135014058 0.0 0.0 0.0 0.0 0.3023441261401367 0.0
 0.0 0.0 0.0 0.676801127801794 0.0 0.0 0.0 0.0 0.3023441261401367 0.0 0.0
 0.6796644127670064 0.0 0.0 0.1920125959201131 0.0 0.0 0.026042138065528686 0.0 0.0 0.0 0.0]
    Lˣ, Lʸ, α, c
end

function problem_specific(::Val{:Rand12})
    Lˣ, Lʸ = 10, 8
    α = [9,8,3,10,8,8,8,3,1,2,8,4]
    c = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.3461964587900239 0.18380626866158423 0.0 0.37615042725365067
 0.0 0.0 0.0 0.0 0.7630077316908308 0.0 0.3712064996520412 0.39425245684528076 0.15004194657037773 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.7528412087794389 0.0 0.0 0.6581231637329426 0.0 0.0 0.3608842601476572 0.0
 0.0 0.0 0.0 0.0 0.0 0.45365915531690204 0.0 0.0 0.37276051649622444 0.0 0.9459258931718393 0.0
 0.0 0.7630077316908308 0.7528412087794389 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.38318117651430694
 0.0 0.0 0.0 0.45365915531690204 0.0 0.0 0.0 0.2402561854030516 0.0 0.6852890352623946 0.36193208333680826 0.0
 0.0 0.3712064996520412 0.0 0.0 0.0 0.0 0.0 0.0 0.009034050537933602 0.0 0.8895683276627007 0.8646601311994306
 0.0 0.39425245684528076 0.6581231637329426 0.0 0.0 0.2402561854030516 0.0 0.0 0.6149644252149276 0.9786535832072509 0.0 0.0
 0.3461964587900239 0.15004194657037773 0.0 0.37276051649622444 0.0 0.0 0.009034050537933602 0.6149644252149276 0.0 0.0 0.3045850350930539 0.0
 0.18380626866158423 0.0 0.0 0.0 0.0 0.6852890352623946 0.0 0.9786535832072509 0.0 0.0 0.0 0.0
 0.0 0.0 0.3608842601476572 0.9459258931718393 0.0 0.36193208333680826 0.8895683276627007 0.0 0.3045850350930539 0.0 0.0 0.2818960645332804
 0.37615042725365067 0.0 0.0 0.0 0.38318117651430694 0.0 0.8646601311994306 0.0 0.0 0.0 0.2818960645332804 0.0]
    Lˣ, Lʸ, α, c
end

function problem_specific(::Val{:Rand13})
    Lˣ, Lʸ = 12, 9
    α = [8,13,9,12,3,8,2,9,12,10,5,1,11]
    c = [0.0 0.0 0.414286710242429 0.0 0.0 0.0 0.0 0.5332972399929927 0.7234154459659385 0.21173685682394439 0.0 0.0 0.0
 0.0 0.0 0.96343563706096 0.0 0.17619287941756667 0.0 0.0 0.0 0.0 0.0 1.7021339452640891 0.0 0.7813668667684333
 0.414286710242429 0.96343563706096 0.0 0.0 0.2552188035621441 0.0 0.0 0.9001704075319319 0.0 1.1436923143730524 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.19568561775658289 0.048898206368805486 0.0 0.04373639235309823 0.4614440155085817 0.0 0.0 0.0
 0.0 0.17619287941756667 0.2552188035621441 0.0 0.0 0.0 0.0 0.8066651557067757 0.0 0.0 0.0 0.0 0.1546833575219424
 0.0 0.0 0.0 0.19568561775658289 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.8600801439975654 0.0
 0.0 0.0 0.0 0.048898206368805486 0.0 0.0 0.0 0.5899797726504763 0.699413350811348 0.0 0.0 0.8198801301148517 0.0
 0.5332972399929927 0.0 0.9001704075319319 0.0 0.8066651557067757 0.0 0.5899797726504763 0.0 0.0 0.0 0.7072561052323294 0.0 0.8862009474224881
 0.7234154459659385 0.0 0.0 0.04373639235309823 0.0 0.0 0.699413350811348 0.0 0.0 0.0 0.1691967772636247 0.0 0.9937467799481123
 0.21173685682394439 0.0 1.1436923143730524 0.4614440155085817 0.0 0.0 0.0 0.0 0.0 0.0 0.3242846674152646 0.0 0.8821774461706708
 0.0 1.7021339452640891 0.0 0.0 0.0 0.0 0.0 0.7072561052323294 0.1691967772636247 0.3242846674152646 0.0 0.0 0.0
 0.0 0.0 0.0 0.0 0.0 0.8600801439975654 0.8198801301148517 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.7813668667684333 0.0 0.0 0.1546833575219424 0.0 0.0 0.8862009474224881 0.9937467799481123 0.8821774461706708 0.0 0.0 0.0]
    Lˣ, Lʸ, α, c
end


text:
function test_fva_consistency()

    # model
    model = MetNets.ecoli_core_model();
    m, n = size(model)

    # parameters
    niter = 10
    robj = MetNets.ECOLI_MODEL_BIOMASS_IDER
    nths = clamp(nthreads() - 1, 1, 5)
    atol = 1e-5

    # Test fba
    fbaout0 = MetLP.fba(model, robj)
    objval0 = MetLP.av(model, fbaout0, robj)
    @show objval0
    @test objval0 > 0.0

    # Testing fva consistency checking and not checking an obj value
    for (testname, check_obj) in [
            ("No-Checking Obj, FVA  ", nothing),
            ("Checking Obj, FVA  ", robj)
        ]

        # fva
        println("-"^60)
        @info("Testing FVA")
        prog = Progress(niter; desc = testname)
        ref_lbs, ref_ubs = MetLP.fva(model; check_obj, nths, verbose = false)
        for _ in 1:niter
            lbs, ubs = MetLP.fva(model; check_obj, nths, verbose = false)
            @test all(isapprox.(lbs, ref_lbs; atol))
            @test all(isapprox.(ubs, ref_ubs; atol))
            next!(prog)
        end
        finish!(prog)

        @info("Testing FBA")
        for (rxni, rxn) in enumerate(model.rxns)

            println("-"^60)
            @show rxn
            @show model.lb[rxni]
            @show model.ub[rxni]

            # max
            fbaout = MetLP.fba(model, rxn, 1; sense1 = MetLP.MIN_SENSE)
            @test MetLP.av(model, fbaout, rxn) == MetLP.av(fbaout, rxni) # test interface
            val = MetLP.av(model, fbaout, rxn)
            @show val
            @show ref_lbs[rxni]
            @test all(isapprox(ref_lbs[rxni], val; atol))

            # min
            fbaout = MetLP.fba(model, rxn, 1; sense1 = MetLP.MAX_SENSE)
            val = MetLP.av(model, fbaout, rxn)
            @test val == MetLP.av(fbaout, rxni)
            @show val
            @show ref_ubs[rxni]
            @test all(isapprox(ref_ubs[rxni], val; atol))

        end
    end
end
test_fva_consistency()


text:
"""
    Ormsby(; <keyword arguments>)

Create a Ormsby wavelet sampled every dt seconds with corner frequencies
defined by the vector f = [f1, f2, f3, f4]. The final wavelet is multiplied by
a Hamming window.

# Arguments
- `dt::AbstractFloat=0.002`: sampling interval in secs.
- `f::Vector{AbstractFloat}=[2.0, 10.0, 40.0, 60.0]`: corner frequencies in Hz.
      ^
    1 |     ***************
      |    *               *
      |   *                 *
      |  *                   *
      | *                     *
      -----------------------------> f
        f1  f2           f3  f4
# Example
```julia
julia> w = Ormsby(); plot(w);
```
"""
function Ormsby(; dt::Tf=0.002, f::Vector{Tf}=[2.0, 10.0, 40.0, 60.0]) where {Tf<:AbstractFloat}

    f1,f2,f3,f4 = f

    fc = (f2+f3)/2.0
    nw = 2.2/(fc*dt)
    nc = floor(Int, nw/2)
    t = dt*collect(-nc:1:nc)
    nw = 2*nc + 1
    a4 = (pi*f4)^2/(pi*(f4-f3))
    a3 = (pi*f3)^2/(pi*(f4-f3))
    a2 = (pi*f2)^2/(pi*(f2-f1))
    a1 = (pi*f1)^2/(pi*(f2-f1))

    u = [ a4 * (sinc(f4 * t[i] ))^2 - a3 * (sinc(f3 * t[i] ))^2 for i in eachindex(t) ]
    v = [ a2 * (sinc(f2 * t[i] ))^2 - a1 * (sinc(f1 * t[i] ))^2 for i in eachindex(t) ]

    w  = u - v
    Hm = Hamming(nw)
    w  = [w[i] * Hm[i] for i in eachindex(w)]
    w  = w/maximum(w)

end


text:
"""
condprob \\

Calculates the distribution of y conditioned on x \\
Arguments:		dec_x = decimalised version of x \\
							dec_y = decimalised version of y \\
							x_dist_dim = dimension of the probability distribution of x \\
							y_dist_dim = dimension of the probability distribution of y \\
							sampling = the sampling method \\
Returns:			ny, each value of y \\
							pxy, probability of y  given x  \\
"""

function condprob(dec_x::Array{Int,1}, dec_y::Array{Int,1}, x_dist_dim::Int, y_dist_dim::Int, sampling::String)
	ny = zeros(Int, y_dist_dim)
	pxy = zeros(x_dist_dim, y_dist_dim)
	for i in unique(dec_y) # looping through all the potential Y values
		inds = find(i .== dec_y) # the trials where Y = i
		ny[find(i .== unique(dec_y))] = length(inds)
		matching_output = dec_x[inds] # output conditional ensemble
		if length(matching_output) > 0
			pxy[:, i+1] = prob(matching_output, x_dist_dim, sampling)
		else
			println("WARNING: Null output conditional ensemble for output = ", string(i))
		end
	end
	return ny, pxy
end


text:
#
# Least-square function and gradient
#
# f is the inner function and gf is the gradient of the
# function applied to each element of the data vector
#
function sq_residue(p,f,X,Y)
  s = 0.
  for i in 1:length(X)
    s += (Y[i] - f(X[i],p))^2
  end
  s
end
function sq_residue_grad!(p,f,∇f,X,Y,g)
  @. g = 0.
  for i in 1:length(X)
    dsdf = -2*(Y[i] - f(X[i],p))
    grad = ∇f(X[i],p)
    for i in 1:length(g)
      g[i] = g[i] + dsdf*grad[i]
    end
  end
end


text:
function run_script(replace_strings, script_name, script_filepath, module_suffix="")
    file_content = read(script_filepath, String)
    test_script_filepath = script_name * "_test.jl"

    for strs in replace_strings
        new_file_content = replace(file_content, strs[1] => strs[2])

        if new_file_content == file_content
            @warn "$(strs[1]) => $(strs[2]) replacement not found in $script_filepath."
            return false
        else
            file_content = new_file_content
        end

    end

    try
        open(test_script_filepath, "w") do f
            # Wrap test script inside module to avoid polluting namespaces
            write(f, "module _Test_$script_name" * "_$module_suffix\n")
            write(f, file_content)
            write(f, "\nend # module")
        end

        include(test_script_filepath)
    catch err
        @warn "Error while testing script: " * sprint(showerror, err)

        # Print the content of the file to the test log, with line numbers, for debugging
        test_file_content = read(test_script_filepath, String)
        delineated_file_content = split(test_file_content, '\n')
        for (number, line) in enumerate(delineated_file_content)
            @printf("% 3d %s\n", number, line)
        end

        return false
    finally
        # Delete the test script
        rm(test_script_filepath, force=true)
    end

    return true
end


text:
# Various aliases to support the Base.Filesystem API or overwrite base behaviour.

# Fallback non-filepath methods if folks explicitly import these symbols
/(args...) = Base.:/(args...)
join(args...) = Base.join(args...)

# Aliases for Base.Filesystem API
# Filesystem methods that were renamed should still have an alias for better interop.
Base.abspath(fp::AbstractPath) = absolute(fp)
Base.ctime(fp::AbstractPath) = datetime2unix(created(fp))
Base.dirname(fp::AbstractPath) = parent(fp)
Base.filemode(fp::AbstractPath) = mode(fp)
Base.isabspath(fp::AbstractPath) = isabsolute(fp)
Base.ispath(fp::AbstractPath) = exists(fp)
Base.joinpath(root::AbstractPath, pieces::Union{AbstractPath, AbstractString}...) = join(root, pieces...)
Base.mkpath(fp::AbstractPath) = mkdir(fp; recursive=true, exist_ok=true)
Base.mtime(fp::AbstractPath) = datetime2unix(modified(fp))
Base.normpath(fp::AbstractPath) = normalize(fp)
Base.realpath(fp::AbstractPath) = canonicalize(fp)
Base.relpath(fp::AbstractPath) = relative(fp)
Base.relpath(fp::AbstractPath, src::AbstractPath) = relative(fp, src)


text:
using Documenter, SpaceAnalysis, PyPlot

makedocs(;
    modules=[SpaceAnalysis],
    authors="Hongyang Zhou <hyzhou@umich.edu> and contributors",
    sitename="SpaceAnalysis",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://henry2004y.github.io/VisAnaJulia",
        assets=String[],
    ),
    pages = ["Home"          => "index.md",
             "Analysis"      => "man/analysis.md",
             "Example"       => "man/examples.md",
             "Log"           => "man/log.md",
             "Parker Spiral" => "man/parker_spiral.md",
             "Internal"      => "man/internal.md"
    ]
)

deploydocs(
    repo = "github.com/henry2004y/VisAnaJulia.git",
    target = "build",
    branch = "gh-pages"
)


text:
#=
    The File contains `CSR5_Sparsik<T>` and related funcs implementation

    `CSR5_Sparsik` implements `AbstractSparseMatrix` interface
=#

# I HAVE BROKE SOMETHING IN MATRIX-VECTOR MULTIPLICATION,
# SO `apply_vector` DOES NOT WORK NOW

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

import Base: ==, !=, +, -, *

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

include("csr_sparsik.jl")

LI = LinearIndices

#using Base

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

# It is convenient to have 0-based indexing in this file

struct Tile
    bit_flag::OffsetMatrix{Bool}
    y_offset::OffsetVector{Int}
    seg_offset::OffsetVector{Int}
    empty_offset::OffsetVector{Int}
end

Base.iterate(t::Tile, s::Int64) = iterate(map(i -> getfield(t, i), 1:nfields(t)), s)

mutable struct CSR5_Sparsik{T<:Field} <: AbstractSparseMatrix{T}
    # the struct uses the matrix storing strategy described in the article
    # "CSR5: An Efficient Storage Format for Cross-Platform
    #        Sparse Matrix-Vector Multiplication" (‡)
    #
    # The paper (‡) is available at
    # https://arxiv.org/pdf/1503.05032.pdf
    #
    # Моя отрада – мысленный полет
    # По книгам, со страницы на страницу.
    # Зимой за чтеньем быстро ночь пройдет,
    # Тепло по телу весело струится,
    # А если попадется редкий том,
    # От радости я на небе седьмом.
    #

    m::Int
    n::Int

    field::T

    # inherited from CSR_Sparsik
    row_ptr::OffsetVector{Int}
    col_idx::OffsetVector{Int}
    val::OffsetVector{<:FieldElem}

    # newcomers
    ω::Int
    σ::Int
    tile_ptr::OffsetVector{Int}
    tiles::OffsetVector{Tile}
end

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

Base.size(A::CSR5_Sparsik) = (A.n, A.m)

nnz(A::CSR5_Sparsik) = length(A.val)

issquare(A::CSR5_Sparsik) = A.n == A.m

order(A::CSR5_Sparsik) = A.n

base_ring(v::CSR5_Sparsik) = v.field

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

# let typemax(Int) == (2^63 - 1) be negative
isneg(x::Int) = x < 0 || x == typemax(x)
rmsign(x::Int) = (x == typemax(x) ? 0 : abs(x))
norm(x::Int) = (x == typemax(x) ? 0 : x)

# Note that a default matrix is linearly iterated column-wise,
# that is, linear iteration over
#   |1  2|
#   |4  8|
# yields [1, 4, 2, 8]
#
# Arrays are stored in a column major order
#

function CSR_to_CSR5(A::CSR_Sparsik)

    m, n = size(A)
    field = base_ring(A)
    row_ptr = OffsetVector(A.row_ptr, Origin(0))
    col_idx = OffsetVector(A.col_idx, Origin(0))
    val = OffsetVector(A.val, Origin(0))

    ω = 4
    σ = 4
    p = ceil(Int, nnz(A) / (ω * σ))

    tile_ptr = OffsetVector(zeros(Int, p + 1), Origin(0))
    for tid in 0 : p
        bnd = tid * ω * σ
        tile_ptr[tid] = searchsortedlast(row_ptr, bnd)
    end
    for tid in 0 : p - 1
        for rid in tile_ptr[tid] : (tile_ptr[tid + 1] - 1)
            if row_ptr[rid] == row_ptr[rid + 1]
                tile_ptr[tid] = - tile_ptr[tid]
                (tile_ptr[tid] == 0) && (tile_ptr[tid] = typemax(Int))
                break
            end
        end
    end

    tiles = OffsetVector(Tile[], Origin(0))

    pos = 1
    for tid in 0 : (p - 2)

        # generating bit_flag
        bit_flag = OffsetArray(fill(false, (ω, σ)), Origin(0, 0))
        bit_flag[1] = true
        for i in pos : length(row_ptr)
            # iterating thorough columns
            idx = row_ptr[i - 1] - tid * (ω * σ)
            if idx > ω * σ
                pos = i
                break
            end
            bit_flag[idx + 1] = true
        end

        # generating y_offset and seg_offset
        y_offset = OffsetVector(zeros(Int, ω), Origin(0))
        seg_offset = OffsetVector(zeros(Int, ω), Origin(0))
        tmp_bit = OffsetVector(fill(false, ω), Origin(0))

        for i in 0 : ω - 1
            for j in 0 : σ - 1
                y_offset[i] = y_offset[i] + bit_flag[j, i]
                tmp_bit[i] = tmp_bit[i] || bit_flag[j, i]
            end
            seg_offset[i] = 1 - tmp_bit[i]
        end
        exclusive_prefix_sum_scan!(y_offset)
        segmented_sum!(seg_offset, tmp_bit)

        # generating empty_offset
        empty_offset = OffsetVector(Int[], Origin(0))
        if isneg(tile_ptr[tid])
            l = reduce(+, bit_flag)
            fill!(resize!(empty_offset, l), 0)
            eid = 0
            for i in 1 : ω * σ
                if bit_flag[i]
                    ptr = tid * ω * σ + i - 1
                    idx = searchsortedlast(row_ptr, ptr)
                    idx -= rmsign(tile_ptr[tid])
                    empty_offset[eid] = idx
                    eid += 1
                end
            end
        end

        push!(tiles, Tile(bit_flag, y_offset, seg_offset, empty_offset))

    end

    CSR5_Sparsik(m, n, field, row_ptr, col_idx, val, ω, σ, tile_ptr, tiles)
end

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

function Base.prod(A::CSR5_Sparsik, x)
    m, n = size(A)
    ω, σ = A.ω, A.σ
    field = base_ring(A)

    val = A.val
    col_idx = A.col_idx
    tile_ptr = A.tile_ptr

    y = OffsetVector(fill(zero(field), n), Origin(0))

    for (tid, tile) in zip(Iterators.countfrom(0), A.tiles)
        bit_flag, y_offset, seg_offset, empty_offset = tile

        tmp = ZeroBased(fill(zero(field), ω))
        last_tmp = ZeroBased(fill(zero(field), ω))
        RED = false

        for i in 0 : ω - 1
            offset = y_offset[i]
            offset = isneg(tile_ptr[tid]) ? empty_offset[offset] : offset

            sum = zero(field)

            for j in 0 : σ - 1

                ptr = tid * ω * σ + j * ω + i
                sum += val[ptr] * x[col_idx[ptr] + 1]

                # the junction
                if bit_flag[i, j]
                    if RED
                        tmp[i - 1] = sum
                        sum = zero(field)
                        RED = false
                    else
                        y[norm(tile_ptr[tid]) + offset] = sum
                        if ! isneg(tile_ptr[tid])
                            y_offset[i] += 1
                        else
                            empty_offset[y_offset[i]] += 1
                        end

                        sum = zero(field)
                    end
                end

            end

            if i + 1 < ω && !bit_flag[i + 1, σ - 1]
                RED = true
            end
            last_tmp[i] = sum
        end

        fast_segmented_sum!(tmp, seg_offset)

        for i in 0 : ω - 1
            last_tmp[i] += tmp[i]
            if ! isneg(tile_ptr[tid])
                y[min(norm(tile_ptr[tid]) + y_offset[i], 7)] = last_tmp[i]
            else
                y[norm(tile_ptr[tid]) + empty_offset[y_offset[i]]] = last_tmp[i]
            end
        end

    end

    return from_dense(y, field)
end

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

function exclusive_prefix_sum_scan!(data)
    # the function uses linear indexing,
    # constructing the result inplace

    # strange to see an OffsetVector failing this test
    # hmm..
    Base.has_fast_linear_indexing(data) || throw(DomainError("too slow"))

    tmp = zero(first(data))
    for i in LI(data)
        data[i], tmp = tmp, tmp + data[i]
    end
    data
end

function inclusive_prefix_sum_scan!(data)
    # the function uses linear indexing,
    # constructing the result inplace

    Base.has_fast_linear_indexing(data) || throw(DomainError("too slow"))

    for i in LI(data)
        i == firstindex(data) && continue;
        data[i] += data[i - 1]
    end
    data
end

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

function fast_segmented_sum!(data, seg_offset)
    # the function implements Algorithm 6 presented in the article (‡)

    LI(data) == LI(seg_offset) || throw(DimensionMismatch("indexing must match"))

    tmp = copy(data)
    inclusive_prefix_sum_scan!(data)

    for i in LI(data)
        data[i] = data[i + seg_offset[i]] - data[i] + tmp[i]
    end

    data
end

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

function segmented_sum!(data, flag)
    # the function implements Algorithm 2 presented in the article (‡)

    LI(data) == LI(flag) || throw(DimensionMismatch("indexing must match"))

    for i in LI(data)
        if flag[i]
            j = i + 1
            while j <= lastindex(flag) && ! flag[j]
                data[i] += data[j]
                j += 1
            end
        else
            data[i] = zero(data[i])
        end
    end

    data
end

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

A = [1 0 2 3 0 0 4 5;
     0 1 0 2 0 0 0 0;
     0 0 0 0 0 0 0 0;
     1 2 3 4 5 0 6 7;
     0 1 0 2 0 3 0 0;
     1 2 0 0 0 0 0 0;
     0 1 2 3 4 5 6 7;
     1 2 3 4 5 6 7 8;]

A_csr = dense_to_CSR(A, QQ);

A_csr5 = CSR_to_CSR5(A_csr)


text:
using Revise

using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using SparseArrays
using Random
# Random.seed!(233)

# TODO: specify your input parameters
A = sprand(10,5,0.3)
f = rand(10)
sol = A\f
u = constant(A)\f

sess = Session()
init(sess)
@show run(sess, u)-sol
# error()

# TODO: change your test parameter to `m`
# gradient check -- v
function scalar_function(m)
    return sum((constant(A)\m)^2)
    B = SparseTensor(ii, jj, m, size(A)...)
    return sum((B\Array([f f]'))^2)
end

ii, jj, vv = find(constant(A))

# TODO: change `m_` and `v_` to appropriate values
# m_ = constant(rand(length(vv)))
# v_ = rand(length(vv))
m_ = constant(rand(5,10))
v_ = rand(5,10)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session()
init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
using Distributions, Plots, LaTeXStrings; pyplot()

dist = TriangularDist(0,2,1)
xGrid = 0:0.01:2
uGrid = 0:0.01:1

p1 = plot( xGrid, pdf.(dist,xGrid), c=:blue,
		xlims=(0,2), ylims=(0,1.1),
		xlabel="x", ylabel="f(x)")
p2 = plot( xGrid, cdf.(dist,xGrid), c=:blue,
		xlims=(0,2), ylims=(0,1),
		xlabel="x", ylabel="F(x)")
p3 = plot( uGrid,quantile.(dist,uGrid), c=:blue,
		xlims=(0,1), ylims=(0,2),
		xlabel="u", ylabel=(L"F^{-1}(u)"))

plot(p1, p2, p3, legend=false, layout=(1,3), size=(1200, 400))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
function SwingOptionPricerSD_advanced(K, T, SMA, SDA, days, Strike, sigMA, sigDA, alpha1, alpha2, Vbar1, Vbar2, xi1, xi2, rhoz, MR_alphaMA, MR_alphaDA, MeanLNMA, MeanLNDA, r, N, M)

## setting parameters

dt = T/N
alpha1dt = alpha1 * dt
alpha2dt = alpha2 * dt

xi1sdt = xi1 * sqrt(dt)
xi2sdt = xi2 * sqrt(dt)
sdt = sqrt(dt)

lnSMA = log(SMA)
lnSDA = log(SDA)

z = zeros(Float64, 4)
pathMA = zeros(Float64, N, M)
pathDA = zeros(Float64, N, M)

## outer loop is the number of simulations bigger is better

for j = 1:M

	lnStMA = lnSMA
	lnStDA = lnSDA
	Vt1 = Vbar1
	Vt2 = Vbar2

	## inner loop is simulating each price path of length N
	for i = 1:N

		z = generateCorrRandom(rhoz)

		pathMA[i,j] = lnStMA
		pathDA[i,j] = lnStDA

		###
		Vt1 = Vt1 + alpha1dt * (Vbar1 - Vt1) + xi1sdt * sqrt(Vt1) * z[3]
		Vt2 = Vt2 + alpha2dt * (Vbar2 - Vt2) + xi2sdt * sqrt(Vt2) * z[4]

		lnStMA = lnStMA + (MR_alphaMA * (MeanLNMA - lnStMA)-0.5*Vt1) * dt + sqrt(Vt1) * sdt * z[1]
		lnStDA = lnStDA + (MR_alphaDA * (MeanLNDA - lnStDA)-0.5*Vt2) * dt + sqrt(Vt2) * sdt * z[2]
	end

end

means = zeros(Float64, 12, M)

for j in 1:M

	for i in 1:12

		means[i, j] = mean(exp(pathMA[21*(i-1)+(1:21), j]))
	end
end

MAAverage = zeros(Float64, (N-21), M)

for j in 1:M

	for i in 1:12

		MAAverage[21*(i-1)+(1:21), j] = fill(means[i,j], 21)
	end

end

diff = exp(pathDA[22:N,1:M]) - (MAAverage + Strike)

resultVector = zeros(Float64, M)

for i in 1:M

	resultVector[i] = mean(sort(diff[1:(N-21),i], rev=true)[1:days])

end

swingValue = mean(resultVector) * exp(-r*T)

SD = std(resultVector)
SE = SD / sqrt(M)

println("The Call Value is $swingValue")
println("The SD is $SD")
println("The SE is $SE")

MC = [swingValue , SD, SE]

end