text:
@testset "Structural interference" begin
    system = GPSL1()

    sat = ConstantDopplerSatellite(
        system,
        1,
        carrier_doppler = 1000.0Hz,
        carrier_phase = π / 2,
        code_phase = 100.0,
        cn0 = 45dBHz
    )
    si = ConstantDopplerStructuralInterference(
        sat,
        -3dB,
        added_carrier_doppler = 10Hz,
        added_carrier_phase = π / 8,
        added_code_phase = 10.0
    )

    @test @inferred(get_carrier_phase(si)) == π / 2 + π / 8
    @test @inferred(get_code_phase(si)) == 110.0

    rng = MersenneTwister(1234)
    signal = @inferred GNSSSimulator.gen_signal!(
        si,
        2.5e6Hz,
        10.0Hz,
        1/Hz,
        2500,
        rng
    )

    code_doppler = 1000.0 * 1023e3 / 1.57542e9
    reference_signal = get_code.(
        system,
        (0:2499) .* (1023e3 .+ code_doppler) ./ 2.5e6 .+ 100 .+ 10.0,
        1
    ) .* cis.(2π .* (0:2499) .* (1000.0 + 10.0 + 10.0) ./ 2.5e6 .+ π / 2 .+ π / 8) .*
        sqrt(10^(42 / 10))
    @test sqrt(sum(abs2.(signal .- reference_signal)) / 2500) < 1e-2 * sqrt(10^(42 / 10))

    next_si = @inferred GNSSSimulator.propagate(si, 1, 100.0Hz, 2e6Hz, 1.0/Hz, Random.GLOBAL_RNG)
    @test @inferred(get_carrier_doppler(next_si)) == 1010Hz
    @test @inferred(get_code_doppler(next_si)) ≈ 1010Hz / 1540
    @test @inferred(get_carrier_phase(next_si)) ≈ π / 2 + π / 8 + 2π * 1110Hz / 2e6Hz
    @test @inferred(get_code_phase(next_si)) ≈ 110 + (1023e3Hz + 1010Hz / 1540) / 2e6Hz
    @test @inferred(get_existence(next_si)) == true
    @test @inferred(get_carrier_to_noise_density_ratio(next_si)) ≈ 42dBHz
    @test @inferred(get_prn(next_si)) == 1
    @test @inferred(get_gnss_system(next_si)) == system
end


text:
"""
    BiNilplex{T <: Real} = Nilplexification{Nlplex{T}}

A bi-nilplex number is a nilplexification of the nilplex numbers.
"""
const BiNilplex{T <: Real} = Nilplexification{Nlplex{T}}

function BiNilplex(z::Nlplex{T}) where T <: Real
    BiNilplex{T}(z, zero(Nlplex{T}))
end

function BiNilplex(a::T, b::T, c::T, d::T) where T <: Real
    BiNilplex{T}(Nlplex{T}(a, b), Nlplex{T}(c, d))
end

function BiNilplex(a::T, b::T, c::T) where T <: Real
    BiNilplex{T}(Nlplex{T}(a, b), Nlplex{T}(c, zero(T)))
end

function BiNilplex(a::T, b::T) where T <: Real
    BiNilplex{T}(Nlplex{T}(a, b), zero(Nlplex{T}))
end

function BiNilplex(a::T) where T <: Real
    BiNilplex{T}(Nlplex{T}(a, zero(T)), zero(Nlplex{T}))
end

function BiNilplex(a::Real, b::Real)
    BiNilplex(promote(a, b)...)
end

function BiNilplex(a::Real, b::Real, c::Real)
    BiNilplex(promote(a, b, c)...)
end

function BiNilplex(a::Real, b::Real, c::Real, d::Real)
    BiNilplex(promote(a, b, c, d)...)
end

"""
    unreal(z::BiNilplex)

The unreal part of a bi-nilplex number is a 3-vector.
"""
function unreal(z::BiNilplex)
    vcat(unreal(z.l), array(z.r))
end

function show(io::IO, z::BiNilplex)
    print(io, "[1: ")
    print(io, z.l.l)
    print(io, ", A: ")
    print(io, z.l.r)
    print(io, ", B: ")
    print(io, z.r.l)
    print(io, ", AB: ")
    print(io, z.r.r)
    print(io, "]")
end

function random(::Type{BiNilplex{T}}) where T <: Real
    BiNilplex{T}(random(Nlplex{T}), random(Nlplex{T}))
end

function barA(z::BiNilplex{T}) where T <: Real
    BiNilplex{T}(barA(z.l), barA(z.r))
end

function barB(z::BiNilplex{T}) where T <: Real
    BiNilplex{T}(z.l, -(z.r))
end

function quadrance(z::BiNilplex)
    (z * barB(z)).l
end

function norm(z::BiNilplex)
    abs2(quadrance(z))
end

function iszerodivisor(z::BiNilplex)
    iszerodivisor(quadrance(z))
end

function inv(z::BiNilplex)
    if iszerodivisor(z)
        error(ZeroDivisorInverse)
    end

    q = quadrance(z)
    barB(z) * barA(q) / abs2(q)
end

ext: jl


text:
"""

	| X | = n

There is no filtrage on it, only pre-treatment

"""

function Cardinalite(model::Model, X::Variable, n::Int)
	X.cardinalInf = n
	X.cardinalSup = n
end


text:
module TestRmul

include("preamble.jl")
using LinearAlgebra

@testset begin
    @test rmul!!(1, 2) === 2

    A = collect(Float64, reshape(1:4, 2, 2))
    B = UpperTriangular(ones(2, 2))
    AB = A * B
    @test rmul!!(A, B) === A
    @test A == AB
    @test_inferred rmul!!(A, B)

    A = collect(Int, reshape(1:4, 2, 2))
    B = UpperTriangular(ones(2, 2))
    AB = A * B
    @test rmul!!(A, B) :: Matrix{Float64} == AB
    @test_inferred rmul!!(A, B)
end

end  # module


text:
# by @walra356
using CairoMakie
CairoMakie.activate!() #HIDE
let
    supertitle = "viewing heatmap detail under conservation of physical scale"
  # define transformation from pixel coordinates to physical coordinates
    edges(p,Δx,x0) = collect(p .* Δx) .-(x0+0.5Δx)
  # parameters camera chip (pixel coordinates)
    dimX = 512; dimY = 512                       # pixel dimensions
    (px, py) = (1:dimX, 1:dimY)                  # pixel positions
    (sx, sy) = (100:150, 380:420)                # idem for submatrix
    (ix, iy) = (125, 400)                        # selected pixel
  # physical coordinates
    (Δx, Δy) = (4.0, 5,0)                       # pixel lattice periods (μm)
    (Ox, Oy) = (0.0, 0.0) # manual offset (μm)
    (x, y) = (edges(px,Δx,Ox), edges(py,Δy,Oy)) # pixel positions (μm)
    (u, v) = (edges(sx,Δx,Ox), edges(sy,Δy,Oy))  # idem for submatrix
  # create exposure matrix σ (model alternative for measurement file)
    valmax = 500.0                              # maximum value exposure matrix
    exposure(i, j , ix, iy, wx, wy) = 0.9valmax *
                        (exp(-(((i-ix)/wx)^2 + ((Δy/Δx)*(j-iy)/wy)^2)) + 0.1)
    σ = round.(Int, [exposure(i, j, ix, iy, 6, 6) for i=1:dimX, j=1:dimY])
    σ1 = σ[px,py] # raw matrix
    σ2 = σ[sx,sy] # submatrix
    footnote = "Origin of coordinate system was shifted by ($Ox μm, $Oy μm)
                w.r.t. left-bottom edge of camera chip
                \nPixel lattice periods: (Δx = $Δx μm, Δy = $Δy μm);
                pixel position of light spot: (ix, y0) = ($(ix), $(iy))"

  # start actual plot
    theme = Theme(fontsize = 12, colormap = :gist_earth, resolution = (750,500))
    set_theme!(theme)
    fig = Figure()
  # collect attributes
    fsize = fig.scene.attributes.fontsize.val
    attr = (xlabelsize = 6fsize/5, ylabelsize = 6fsize/5, titlesize = 7fsize/5,
            xautolimitmargin = (0,0), yautolimitmargin = (0,0), )
    attr1 = (attr..., title = "raw data",
            xticks = [px[1], ix, px[end]] , yticks = [py[1], iy, py[end]] ,
            xlabel = "x (pixel)", ylabel = "y (pixel)",
            aspect = (Δx * size(σ1,1))/(Δy * size(σ1,2)), )
    attr2 = (attr..., title = "physical scale",
            xticks = [x[1]-Δx/2, x[ix], x[end]+Δx/2] ,
            yticks = [y[1]-Δy/2, y[iy], y[end]+Δy/2] ,
            xlabel = "x (μm)", ylabel = "y (μm)",
            aspect = (Δx * size(σ1,1))/(Δy * size(σ1,2)), )
    attr3 = (attr..., title = "physical detail",
            xticks = [u[1]-Δx/2, u[1+ix-sx[1]], u[end]+Δx/2] ,
            yticks = [v[1]-Δy/2, v[1+iy-sy[1]], v[end]+Δy/2] ,
            xlabel = "x (μm)", ylabel = "y (μm)",
            aspect = (Δx * size(σ2,1))/(Δy * size(σ2,2)), )
    attrc = (label = "exposure (counts)", ticksize=6fsize/5, tickalign = 1,
            width = 15, height = Relative(0.54), ) # tweaked
  # create axes and add contents
    ax1 = Axis(fig; attr1...)
    ax2 = Axis(fig; attr2...)
    ax3 = Axis(fig; attr3...)
    hm1 = heatmap!(ax1, px, py, σ1)
    hm2 = heatmap!(ax2, x, y, σ1)
    hm3 = heatmap!(ax3, u, v, σ2)
    axc = Colorbar(fig, hm1; attrc...)
    axn = Label(fig, text = footnote, textsize = 6fsize/5) # footnote
    axt = Label(fig, text = supertitle, textsize = 2fsize) # supertitle
  # create layout and show figure
    fig[1,1] = ax1
    fig[1,2] = ax2
    fig[1,3] = ax3
    fig[1,4] = axc
    fig[2,:] = axn
    fig[0,:] = axt
    fig
    save(joinpath(@__DIR__, "output", "heatmapScaleDetail.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
module TestArrayOfComponents

println("\nTestArrayOfComponents: Demonstrating the handling of arrays of components")

using Modia
using Modia.Electric
using ModiaMath.plot
using Base.Test

@model LPfilter begin
  R=Resistor(R=100.0)
  C=Capacitor(C=2.5E-3, v = Float(start=0.0))
  V=ConstantVoltage(V=10.0)
  ground=Ground()
@equations begin
  connect(V.n, ground.p)
  connect(V.p, R.p)
  connect(R.n, C.p)
  connect(C.n, V.n)
  end
end

# Array of components
@model TwoFilters begin
  F = [LPfilter(); LPfilter(C=Capacitor(C=10*2.5E-3))]
end
checkSimulation(TwoFilters, 1, "F[2].C.v", 3.2967995487381305)

nFilters=10
# Array comprehensions of components
@model ManyFilters begin
  F = [LPfilter() for i in 1:nFilters]
end
checkSimulation(ManyFilters, 1, "F[1].C.v", 9.816405569531037)

nFilters=10
# Array comprehensions of components
@model ManyDifferentFilters begin
  F = [LPfilter(C=Capacitor(C=i*2.5E-3)) for i in 1:nFilters]
end
result = simulate(ManyDifferentFilters, 1)
plot(result, Tuple(["F[$i].C.v" for i in 1:nFilters]), heading="ManyDifferentFilters", figure=5)
@test result["F[1].C.v"][end] == 9.816405569531037
@test result["F[$nFilters].C.v"][end] == 3.2967995487381305

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

# Connecting arrays of components
@model AdvancedLPfilter begin
  M = [Resistor(R=1); Capacitor(C=1E-3)]
  V=ConstantVoltage(V=1)
  ground=Ground()
@equations begin
  connect(V.n, ground.p)
#  connect(V.p, M[1].p)
#  connect(M[1].n, M[2].p)
#  connect(M[2].n, V.n)
  end
end

#=
# Alternative connect constructs:
for i in 1:size(M)
  connect(M[i].p, M[i+1].n)
end
connect(M[1:end-1].p, M[2:end].n)]
[connect(M[i].p, M[i+1].n) for i in 1:size(M)]
connect(M[i].p, M[i+1].n) for i in 1:size(M)
=#

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

# Experimental:
using Modia.@equation
addEquation!(M, e) = begin @show e; push!(M.initializers, Modia.Instantiation.Equations([e])) end

@model LPfilterComponent begin
  R=Resistor(R=100.0)
  C=Capacitor(C=2.5E-3, v = Float(start=0.0))
  ground=Ground()
  in = Pin()
  out = Pin()
@equations begin
  connect(in, R.p)
  connect(R.n, C.p)
  connect(C.n, ground.p)
  connect(C.p, out)
  end
end

@model ManyConnectedDifferentFilters begin
  V=ConstantVoltage(V=10.0)
  ground=Ground()
  F = [LPfilterComponent(C=Capacitor(C=i*2.5E-3)) for i in 1:nFilters]
@equations begin
  connect(V.n, ground.p)
#  connect(V.p, F[1].in)
  end
end

#=
Not working yet
for i in 1:nFilters-1
  addEquation!(ManyConnectedDifferentFilters, @equation(connect(this.F[$i].out, this.F[$(i+1)].in)))
end

simulate(ManyConnectedDifferentFilters, 1)
=#

end


text:
module CreateTablePermissionsRoles

import SearchLight.Migrations: create_table, column, primary_key, add_index, drop_table

function up()
  create_table(:permissionsroles) do
    [
      primary_key()
      column(:permissions_id, :int)
      column(:roles_id, :int)
    ]
  end

  add_index(:permissionsroles, :permissions_id)
end

function down()
  drop_table(:permissionsroles)
end

end


text:
module Indexing

using DocStringExtensions
using TriMatrices: TriLayout, TriLower, TriUpper, TriSymmetric,
                   nelems, hasdiag

include("math.jl")
include("conversion.jl")
include("tri_indices.jl")

end  # module


text:
function make_operation(@nospecialize(op), args)
    if op === (*)
        args = filter(!_isone, args)
        if isempty(args)
            return 1
        end
    elseif op === (+)
        args = filter(!_iszero, args)
        if isempty(args)
            return 0
        end
    end
    return op(args...)
end

function build_expr(head::Symbol, args)
    ex = Expr(head)
    append!(ex.args, args)
    ex
end

# used in parsing
isblock(x) = length(x) == 1 && x[1] isa Expr && x[1].head == :block
function flatten_expr!(x)
    isb = isblock(x)
    if isb
        x = MacroTools.striplines(x[1])
        filter!(z->z isa Symbol || z.head != :line, x.args)
        x = (x.args...,)
    end
    x
end

function detime_dvs(op::Term)
  if op.op isa Sym
      Sym{Real}(nameof(op.op))
  else
    Term(op.op,detime_dvs.(op.args))
  end
end
detime_dvs(op) = op

function retime_dvs(op::Sym,dvs,iv)
    Sym{FnType{Tuple{symtype(iv)}, Real}}(nameof(op))(iv)
end

function retime_dvs(op::Term, dvs, iv)
    similarterm(op, op.op, retime_dvs.(op.args,(dvs,),(iv,)))
end
retime_dvs(op,dvs,iv) = op

is_derivative(O::Term) = isa(O.op, Differential)
is_derivative(::Any) = false

"""
get_variables(O)

Returns the variables in the expression
"""
get_variables(e::Num, varlist=nothing) = get_variables(value(e), varlist)
get_variables!(vars, e, varlist=nothing) = vars

is_singleton(e::Term) = e.op isa Sym
is_singleton(e::Sym) = true
is_singleton(e) = false

get_variables!(vars, e::Number, varlist=nothing) = vars

function get_variables!(vars, e::Symbolic, varlist=nothing)
    if is_singleton(e)
        if isnothing(varlist) || any(isequal(e), varlist)
            push!(vars, e)
        end
    else
        foreach(x -> get_variables!(vars, x, varlist), e.args)
    end
    return (vars isa AbstractVector) ? unique!(vars) : vars
end

function get_variables!(vars, e::Equation, varlist=nothing)
  get_variables!(vars, e.rhs, varlist)
end

get_variables(e, varlist=nothing) = get_variables!([], e, varlist)

modified_states!(mstates, e::Equation, statelist=nothing) = get_variables!(mstates, e.lhs, statelist)

# variable substitution
# Piracy but mild
"""
substitute(expr, s::Pair)
substitute(expr, s::Dict)
substitute(expr, s::Vector)

Performs the substitution `Num => val` on the `expr` Num.
"""
substitute(expr::Num, s::Union{Pair, Vector, Dict}; kw...) = Num(substituter(s)(value(expr); kw...))
# TODO: move this to SymbolicUtils
substitute(expr::Term, s::Pair; kw...) = substituter([s[1] => s[2]])(expr; kw...)
substitute(expr::Term, s::Vector; kw...) = substituter(s)(expr; kw...)

substituter(pair::Pair) = substituter((pair,))
function substituter(pairs)
    dict = Dict(to_symbolic(k) => to_symbolic(v)  for (k, v) in pairs)
    (expr; kw...) -> SymbolicUtils.substitute(expr, dict; kw...)
end

macro showarr(x)
    n = string(x)
    quote
        y = $(esc(x))
        println($n, " = ", summary(y))
        Base.print_array(stdout, y)
        println()
        y
    end
end

@deprecate substitute_expr!(expr,s) substitute(expr,s)

function states_to_sym(states::Set)
    function _states_to_sym(O)
        if O isa Equation
            Expr(:(=), _states_to_sym(O.lhs), _states_to_sym(O.rhs))
        elseif O isa Term
            if isa(O.op, Sym)
                O in states && return tosymbol(O)
                # dependent variables
                return build_expr(:call, Any[O.op.name; _states_to_sym.(O.args)])
            else
                return build_expr(:call, Any[O.op; _states_to_sym.(O.args)])
            end
        elseif O isa Num
            return _states_to_sym(value(O))
        else
            return toexpr(O)
        end
    end
end
states_to_sym(states) = states_to_sym(Set(states))

"""
    toparam(s::Sym) -> Sym{<:Parameter}

Maps the variable to a paramter.
"""
toparam(s::Sym) = Sym{Parameter{symtype(s)}}(s.name)
toparam(s::Sym{<:Parameter}) = s

"""
    tovar(s::Sym) -> Sym{Real}

Maps the variable to a variable (state).
"""
tovar(s::Sym{<:Parameter}) = Sym{symtype(s)}(s.name)
tovar(s::Sym) = s


text:
function GenAmericanBinomialDownOutCall(K, T, S, sig, r, H, N)

        ## as discussed in clewlow and strickland chapter 2, figure 2.12
        ## deviates as different stepping through the matrix requires sometimes
        ## to make use of * instead of /

  dt = T/N ## Time to Maturity/ Number of Steps
  nu = r - 0.5*sig^2
  dxu = sqrt(sig^2*dt + (nu*dt)^2)
  dxd = -dxu
  pu = 1/2 + 1/2*(nu*dt/dxu)
  pd = 1 - pu

  disc = exp(-r*dt)
  dpu = disc*pu
  dpd = disc * pd
  edxud = exp(dxu - dxd)
  edxd = exp(dxd)

  St = zeros(Float64, (N+1, N+1))

        # Initialising asset prices at maturity

  St[1 , N+1] = S / exp(N * dxd)
  for j in 2:(N+1)
    St[j, N+1] = St[j-1, N+1] / edxud
  end

  # create Array for options value tree

  P = zeros(Float64, (N+1, N+1))#linspace(1, N+1, N+1)

  # Initialise Option Values at maturity

for j in 1:(N+1)

  if St[j, N+1] > H ## check if market price is above knock out

    C[j, N+1] = max(0.0, St[j, N+1] - K)

  else

    C[j, N+1] = 0.0 ## if below, assign 0 to the call value

  end

end

  # Stepping back through the tree

for i = linrange(N, 1, N)

  for j in 1:i

    St[j,i] = St[j,i+1] * edxd  ## generating the asset price at node i, j first

    if St[j,i] > H ## checking if asset above knock out

      C[j, i] = dpu * C[j, i+1] + dpd * C[j+1, i+1]
      C[j, i] = max(C[j, i], St[j,i] - K)

    else

      C[j,i] = 0.0  ## apply 0 if not

    end

  end

end

C, St
end


text:
using IHT
using SnpArrays
using DataFrames
using Distributions
using BenchmarkTools
using Random
using LinearAlgebra
using StatsFuns: logistic
using DelimitedFiles

#run n repeats of L0_reg using the same X and b, but different y
function normal_repeat(
    repeats  :: Int,
    x        :: SnpArray,
    xbm      :: SnpBitMatrix,
    true_b   :: Vector{Float64},
    noise    :: Float64,
    position :: Vector{Int} #correct position of the true model
)
    n = size(xbm, 1)
    k = size(position, 1)
    estimated_β = zeros(k, repeats)

    for i in 1:repeats
        # simulate random noise
        ϵ = rand(Normal(0, noise), n)

        #simulate phenotypes (e.g. vector y) via: y = Xb + noise
        y = xbm * true_b + ϵ

        #compute IHT result for less noisy data
        z = ones(n, 1) # intercept
        result = L0_reg(x, z, y, 1, k, debias=false)

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

    return estimated_β
end

#run n repeats of L0_reg using the same X and b, but different y
function logistic_repeat(
    repeats  :: Int,
    x        :: SnpArray,
    xbm      :: SnpBitMatrix,
    true_b   :: Vector{Float64},
    noise    :: Float64,
    position :: Vector{Int} #correct position of the true model
)
    n = size(xbm, 1)
    k = size(position, 1)
    estimated_β = zeros(k, repeats)

    for i in 1:repeats
        #simulate phenotypes
        y_temp = xbm * true_b

        # Apply inverse logit link and sample from the vector of distributions
        prob = logistic.(y_temp) #inverse logit link
        y = [rand(Bernoulli(x)) for x in prob]
        y = Float64.(y)

        #compute IHT result for less noisy data
        z = ones(n, 1) # intercept
        result = L0_logistic_reg(x, z, y, 1, k, glm = "logistic", debias=true, show_info=false)

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

    return estimated_β
end

#run n repeats of L0_reg using the same X and b, but different y
function poisson_repeat(
    repeats  :: Int,
    x        :: SnpArray,
    xbm      :: SnpBitMatrix,
    true_b   :: Vector{Float64},
    noise    :: Float64,
    position :: Vector{Int} #correct position of the true model
)
    n = size(xbm, 1)
    k = size(position, 1)
    estimated_β = zeros(k, repeats)

    for i in 1:repeats
        #simulate phenotypes
        y_temp = xbm * true_b

        # Simulate poisson data
        λ = exp.(y_temp) #inverse log link
        y = [rand(Poisson(x)) for x in λ]
        y = Float64.(y)

        #compute IHT result for less noisy data
        z = ones(n, 1) # intercept
        result = L0_poisson_reg(x, z, y, 1, k, glm = "poisson", debias=false, convg=false, show_info=false)

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

    return estimated_β
end

function run()
    repeats = 10

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

    #simulat data
    n = 2000
    p = 10000
    bernoulli_rates = 0.5rand(p) #minor allele frequencies are drawn from uniform (0, 0.5)
    x = simulate_random_snparray(n, p, bernoulli_rates)
    xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);

    #specify true model size and noise of data
    k = 10   # number of true predictors
    s = 0.1  # noise

    #construct true model b
    true_b = zeros(p)       # model vector
    true_b[1:k] = randn(k)  # Initialize k non-zero entries in the true model
    shuffle!(true_b)        # Shuffle the entries
    correct_position = findall(x -> x != 0, true_b) # keep track of what the true entries are
    true_b[correct_position]

    #run repeats and save to file
    normal_result = normal_repeat(repeats, x, xbm, true_b, s, correct_position)
    logistic_result = logistic_repeat(repeats, x, xbm, true_b, s, correct_position)
    poisson_result = poisson_repeat(repeats, x, xbm, true_b, s, correct_position)

    writedlm("./repeats/normal_$repeats", normal_result)
    writedlm("./repeats/logistic_$repeats", logistic_result)
    writedlm("./repeats/poisson_$repeats", poisson_result)

    println("completed $repeats repeats.")
end

run()


text:
using JuExample
using Test

@testset "JuExample.jl" begin
  @test JuExample.test_func(0,1) == 1
  @test JuExample.test_func(1,1) == 2
  @test JuExample.test_func(10,1) == 11
  @test JuExample.test_func(2,3) == 5
  @test JuExample.test_func(2,3) == 5
  @test JuExample.test_func(2,3) == 5
  @test JuExample.test_func(2,3) == 5
end


text:
# This file is based on a part of Julia. License is MIT: https://julialang.org/license

import Base: permutedims, permutedims!
export PermutedDimsQuasiArray

# Some day we will want storage-order-aware iteration, so put perm in the parameters
struct PermutedDimsQuasiArray{T,N,perm,iperm,AA<:AbstractQuasiArray} <: AbstractQuasiArray{T,N}
    parent::AA

    function PermutedDimsQuasiArray{T,N,perm,iperm,AA}(data::AA) where {T,N,perm,iperm,AA<:AbstractQuasiArray}
        (isa(perm, NTuple{N,Int}) && isa(iperm, NTuple{N,Int})) || error("perm and iperm must both be NTuple{$N,Int}")
        isperm(perm) || throw(ArgumentError(string(perm, " is not a valid permutation of dimensions 1:", N)))
        all(map(d->iperm[perm[d]]==d, 1:N)) || throw(ArgumentError(string(perm, " and ", iperm, " must be inverses")))
        new(data)
    end
end

"""
    PermutedDimsQuasiArray(A, perm) -> B

Given an AbstractQuasiArray `A`, create a view `B` such that the
dimensions appear to be permuted. Similar to `permutedims`, except
that no copying occurs (`B` shares storage with `A`).

See also: [`permutedims`](@ref).

# Examples
```jldoctest
julia> A = rand(3,5,4);

julia> B = PermutedDimsQuasiArray(A, (3,1,2));

julia> size(B)
(4, 3, 5)

julia> B[3,1,2] == A[1,2,3]
true
```
"""
function PermutedDimsQuasiArray(data::AbstractQuasiArray{T,N}, perm) where {T,N}
    length(perm) == N || throw(ArgumentError(string(perm, " is not a valid permutation of dimensions 1:", N)))
    iperm = invperm(perm)
    PermutedDimsQuasiArray{T,N,(perm...,),(iperm...,),typeof(data)}(data)
end

Base.parent(A::PermutedDimsQuasiArray) = A.parent
Base.size(A::PermutedDimsQuasiArray{T,N,perm}) where {T,N,perm} = genperm(size(parent(A)), perm)
Base.axes(A::PermutedDimsQuasiArray{T,N,perm}) where {T,N,perm} = genperm(axes(parent(A)), perm)

Base.similar(A::PermutedDimsQuasiArray, T::Type, dims::Base.Dims) = similar(parent(A), T, dims)

@inline function _getindex(::Type{IND}, A::PermutedDimsQuasiArray{T,N,perm,iperm}, I::IND) where {T,N,perm,iperm,IND}
    @boundscheck checkbounds(A, I...)
    @inbounds val = getindex(A.parent, genperm(I, iperm)...)
    val
end
@inline function _setindex!(::Type{IND}, A::PermutedDimsQuasiArray{T,N,perm,iperm}, val, I::IND) where {T,N,perm,iperm,IND}
    @boundscheck checkbounds(A, I...)
    @inbounds setindex!(A.parent, val, genperm(I, iperm)...)
    val
end

@inline genperm(I::NTuple{N,Any}, perm::Dims{N}) where {N} = ntuple(d -> I[perm[d]], Val(N))
@inline genperm(I, perm::AbstractVector{Int}) = genperm(I, (perm...,))

# function permutedims(A::AbstractQuasiArray, perm)
#     dest = similar(A, genperm(axes(A), perm))
#     permutedims!(dest, A, perm)
# end

permutedims(A::AbstractQuasiMatrix) = PermutedDimsQuasiArray(A, (2,1))
# permutedims(v::AbstractQuasiVector) = reshape(v, (1, length(v)))

# function Base.showarg(io::IO, A::PermutedDimsQuasiArray{T,N,perm}, toplevel) where {T,N,perm}
#     print(io, "PermutedDimsQuasiArray(")
#     Base.showarg(io, parent(A), false)
#     print(io, ", ", perm, ')')
#     toplevel && print(io, " with eltype ", eltype(A))
# end

## diagonal special case
permutedims(D::Diagonal{<:Any,<:SubArray{<:Any,1,<:AbstractQuasiArray}}) = D


text:
function invert_rxn!(metnet::MetNet, ider::IDER_TYPE;
        rename::Union{String,Function} = (rxn) -> string(rxn)
    )
    idx = rxnindex(metnet, ider)
    metnet.S[:,idx] .*= -1
    ub_ = -metnet.lb[idx]
    lb_ = -metnet.ub[idx]
    metnet.ub[idx] = ub_
    metnet.lb[idx] = lb_
    metnet.rxns[idx] = rename isa Function ? rename(metnet.rxns[idx]) : rename
    return metnet
end


text:
using ADCPDataProcessing, PIEMetData, TidalDischargeModels, TidalFluxExampleData
TDM = TidalDischargeModels

setADCPdatadir!(Pkg.dir("TidalFluxExampleData","data","adcp"))
setmetdatadir!(Pkg.dir("TidalFluxExampleData","data","met"))

creek = Creek{:sweeney}()

include("deployments.jl")
include("adcpdata.jl")
include("crosssections.jl")
include("dischargedata.jl")

N1 = length(Q1)
N2 = length(Q2)

include("boon.jl")
include("lti.jl")
include("volterra.jl")
include("kmeans.jl")

include("summarystats.jl")


text:
for jy in ("j","y"), nu in (0,1)
    jynu = Expr(:quote, symbol(string(jy,nu)))
    jynuf = Expr(:quote, symbol(string(jy,nu,"f")))
    bjynu = symbol(string("bessel",jy,nu))
    if jy == "y"
        @eval begin
            $bjynu(x::Float64) = nan_dom_err(ccall(($jynu,libm),  Float64, (Float64,), x), x)
            $bjynu(x::Float32) = nan_dom_err(ccall(($jynuf,libm), Float32, (Float32,), x), x)
        end
    else
        @eval begin
            $bjynu(x::Float64) = ccall(($jynu,libm),  Float64, (Float64,), x)
            $bjynu(x::Float32) = ccall(($jynuf,libm), Float32, (Float32,), x)
        end
    end
    @eval begin
        $bjynu(x::Real) = $bjynu(float(x))
        $bjynu(x::Complex) = $(symbol(string("bessel",jy)))($nu,x)
        @vectorize_1arg Number $bjynu
    end
end

type AmosException <: Exception
    info::Int32
end

let
    const ai::Array{Float64,1} = Array(Float64,2)
    const ae::Array{Int32,1} = Array(Int32,2)
global airy
function airy(k::Int, z::Complex128)
    id = int32(k==1 || k==3)
    if k == 0 || k == 1
        ccall((:zairy_,openspecfun), Void,
              (Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
               Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
              &real(z), &imag(z),
              &id, &1,
              pointer(ai,1), pointer(ai,2),
              pointer(ae,1), pointer(ae,2))
        if ae[2] == 0 || ae[2] == 3
            return complex(ai[1],ai[2])
        else
            throw(AmosException(ae[2]))
        end
    elseif k == 2 || k == 3
        ccall((:zbiry_,openspecfun), Void,
              (Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
               Ptr{Float64}, Ptr{Float64}, Ptr{Int32}),
              &real(z), &imag(z),
              &id, &1,
              pointer(ai,1), pointer(ai,2),
              pointer(ae,1))
        if ae[1] == 0 || ae[1] == 3  # ignore underflow
            return complex(ai[1],ai[2])
        else
            throw(AmosException(ae[2]))
        end
    else
        error("invalid argument")
    end
end
end

airy(z) = airy(0,z)
@vectorize_1arg Number airy
airyprime(z) = airy(1,z)
@vectorize_1arg Number airyprime
airyai(z) = airy(0,z)
@vectorize_1arg Number airyai
airyaiprime(z) = airy(1,z)
@vectorize_1arg Number airyaiprime
airybi(z) = airy(2,z)
@vectorize_1arg Number airybi
airybiprime(z) = airy(3,z)
@vectorize_1arg Number airybiprime

airy(k::Number, x::FloatingPoint) = oftype(x, real(airy(k, complex(x))))
airy(k::Number, x::Real) = airy(k, float(x))
airy(k::Number, z::Complex64) = complex64(airy(k, complex128(z)))
airy(k::Number, z::Complex) = airy(convert(Int,k), complex128(z))
@vectorize_2arg Number airy

const cy = Array(Float64,2)
const ae = Array(Int32,2)
const wrk = Array(Float64,2)

function _besselh(nu::Float64, k::Integer, z::Complex128)
    ccall((:zbesh_,openspecfun), Void,
          (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
           Ptr{Int32}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
          &real(z), &imag(z), &nu, &1, &k, &1,
          pointer(cy,1), pointer(cy,2),
          pointer(ae,1), pointer(ae,2))
    if ae[2] == 0 || ae[2] == 3
        return complex(cy[1],cy[2])
    else
        throw(AmosException(ae[2]))
    end
end

function _besseli(nu::Float64, z::Complex128)
    ccall((:zbesi_,openspecfun), Void,
          (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
           Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
          &real(z), &imag(z), &nu, &1, &1,
          pointer(cy,1), pointer(cy,2),
          pointer(ae,1), pointer(ae,2))
    if ae[2] == 0 || ae[2] == 3
        return complex(cy[1],cy[2])
    else
        throw(AmosException(ae[2]))
    end
end

function _besselj(nu::Float64, z::Complex128)
    ccall((:zbesj_,openspecfun), Void,
          (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
           Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
          &real(z), &imag(z), &nu, &1, &1,
          pointer(cy,1), pointer(cy,2),
          pointer(ae,1), pointer(ae,2))
    if ae[2] == 0 || ae[2] == 3
        return complex(cy[1],cy[2])
    else
        throw(AmosException(ae[2]))
    end
end

function _besselk(nu::Float64, z::Complex128)
    ccall((:zbesk_,openspecfun), Void,
          (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32},
           Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),
          &real(z), &imag(z), &nu, &1, &1,
          pointer(cy,1), pointer(cy,2),
          pointer(ae,1), pointer(ae,2))
    if ae[2] == 0 || ae[2] == 3
        return complex(cy[1],cy[2])
    else
        throw(AmosException(ae[2]))
    end
end

function _bessely(nu::Float64, z::Complex128)
    ccall((:zbesy_,openspecfun), Void,
          (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32},
           Ptr{Int32}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32},
           Ptr{Float64}, Ptr{Float64}, Ptr{Int32}),
          &real(z), &imag(z), &nu, &1, &1,
          pointer(cy,1), pointer(cy,2),
          pointer(ae,1), pointer(wrk,1),
          pointer(wrk,2), pointer(ae,2))
    if ae[2] == 0 || ae[2] == 3
        return complex(cy[1],cy[2])
    else
        throw(AmosException(ae[2]))
    end
end

function besselh(nu::Float64, k::Integer, z::Complex128)
    if nu < 0
        s = (k == 1) ? 1 : -1
        return _besselh(-nu, k, z) * complex(cospi(nu),-s*sinpi(nu))
    end
    return _besselh(nu, k, z)
end

function besseli(nu::Float64, z::Complex128)
    if nu < 0
        return _besseli(-nu,z) - 2_besselk(-nu,z)*sinpi(nu)/pi
    else
        return _besseli(nu, z)
    end
end

function besselj(nu::Float64, z::Complex128)
    if nu < 0
        return _besselj(-nu,z)cos(pi*nu) + _bessely(-nu,z)*sinpi(nu)
    else
        return _besselj(nu, z)
    end
end

besselj(nu::Integer, x::FloatingPoint) = typemin(Int32) <= nu <= typemax(Int32) ?
    oftype(x, ccall((:jn, libm), Float64, (Cint, Float64), nu, x)) :
    besselj(float64(nu), x)

besselj(nu::Integer, x::Float32) = typemin(Int32) <= nu <= typemax(Int32) ?
    ccall((:jnf, libm), Float32, (Cint, Float32), nu, x) :
    besselj(float64(nu), x)

besselk(nu::Float64, z::Complex128) = _besselk(abs(nu), z)

function bessely(nu::Float64, z::Complex128)
    if nu < 0
        return _bessely(-nu,z)*cospi(nu) - _besselj(-nu,z)*sinpi(nu)
    else
        return _bessely(nu, z)
    end
end

besselh(nu, z) = besselh(nu, 1, z)
besselh(nu::Real, k::Integer, z::Complex64) = complex64(besselh(float64(nu), k, complex128(z)))
besselh(nu::Real, k::Integer, z::Complex) = besselh(float64(nu), k, complex128(z))
besselh(nu::Real, k::Integer, x::Real) = besselh(float64(nu), k, complex128(x))
@vectorize_2arg Number besselh

besseli(nu::Real, z::Complex64) = complex64(besseli(float64(nu), complex128(z)))
besseli(nu::Real, z::Complex) = besseli(float64(nu), complex128(z))
besseli(nu::Real, x::Integer) = besseli(nu, float64(x))
function besseli(nu::Real, x::FloatingPoint)
    if x < 0 && !isinteger(nu)
        throw(DomainError())
    end
    oftype(x, real(besseli(float64(nu), complex128(x))))
end
@vectorize_2arg Number besseli

function besselj(nu::FloatingPoint, x::FloatingPoint)
    if isinteger(nu)
        if typemin(Int32) <= nu <= typemax(Int32)
            return besselj(int(nu), x)
        end
    elseif x < 0
        throw(DomainError())
    end
    oftype(x, real(besselj(float64(nu), complex128(x))))
end

besselj(nu::Real, z::Complex64) = complex64(besselj(float64(nu), complex128(z)))
besselj(nu::Real, z::Complex) = besselj(float64(nu), complex128(z))
besselj(nu::Real, x::Integer) = besselj(nu, float(x))
@vectorize_2arg Number besselj

besselk(nu::Real, z::Complex64) = complex64(besselk(float64(nu), complex128(z)))
besselk(nu::Real, z::Complex) = besselk(float64(nu), complex128(z))
besselk(nu::Real, x::Integer) = besselk(nu, float64(x))
function besselk(nu::Real, x::FloatingPoint)
    if x < 0
        throw(DomainError())
    end
    if x == 0
        return oftype(x, Inf)
    end
    oftype(x, real(besselk(float64(nu), complex128(x))))
end
@vectorize_2arg Number besselk

bessely(nu::Real, z::Complex64) = complex64(bessely(float64(nu), complex128(z)))
bessely(nu::Real, z::Complex) = bessely(float64(nu), complex128(z))
bessely(nu::Real, x::Integer) = bessely(nu, float64(x))
function bessely(nu::Real, x::FloatingPoint)
    if x < 0
        throw(DomainError())
    end
    if isinteger(nu) && typemin(Int32) <= nu <= typemax(Int32)
        return bessely(int(nu), x)
    end
    oftype(x, real(bessely(float64(nu), complex128(x))))
end
function bessely(nu::Integer, x::FloatingPoint)
    if x < 0
        throw(DomainError())
    end
    return oftype(x, ccall((:yn, libm), Float64, (Cint, Float64), nu, x))
end
function bessely(nu::Integer, x::Float32)
    if x < 0
        throw(DomainError())
    end
    return ccall((:ynf, libm), Float32, (Cint, Float32), nu, x)
end
@vectorize_2arg Number bessely

hankelh1(nu, z) = besselh(nu, 1, z)
@vectorize_2arg Number hankelh1

hankelh2(nu, z) = besselh(nu, 2, z)
@vectorize_2arg Number hankelh2

max_stars_repo_name: vchuravy/julia

max_issues_repo_name: vchuravy/julia

max_forks_repo_name: vchuravy/julia


text:
import FractionalTimeDG.P, FractionalTimeDG.dP

correct_P(x) = [ 1, x, (3x^2-1)/2, (5x^3-3x)/2, (35x^4-30x^2+3)/8 ]

J = 11
τ = range(-1, 1, length=J)

Ψ = zeros(5, J)
FractionalTimeDG.legendre_polys!(Ψ, τ)
err = [ Ψ[i,j] - correct_P(τ[j])[i] for i = 1:5, j = 1:J ]
@test all( abs.(err) .< 1e-10 )

err = [ Ψ[i,j] - P(i-1, τ[j]) for i = 1:5, j = 1:J ]
@test all( abs.(err) .< 1e-10 )

correct_dP(x) = [ 0, 1, 3x, (15x^2-3)/2, (35x^3-15x)/2 ]
dΨ = zeros(5, J)
FractionalTimeDG.deriv_legendre_polys!(dΨ, τ)
err = [ dΨ[i,j] - correct_dP(τ[j])[i] for i = 1:5, j = 1:J ]
@test all( abs.(err) .< 1e-10 )

err = [ dΨ[i,j] - dP(i-1, τ[j]) for i = 1:5, j = 1:J ]
@test all( abs.(err) .< 1e-10 )

u(t) = exp(-t) * sin(2π*t)
α = 0.5
rmax = 6
Mmax = 3rmax
ppImax = 20
N = 15
store = Store(α, rmax, Mmax, ppImax)
q = parse(T, "1.5")
t = OffsetArray(T[ (n/N)^q for n = 0:N ], 0:N)
U = FractionalTimeDG.Fourier_Legendre_coefs(u, rmax, t, store)
τ = range(-1, 1, length=ppImax)
pcwise_t, pcwise_U = FractionalTimeDG.evaluate_pcwise_poly!(U, t, τ, store)
pcwise_u = [ u(t) for t in pcwise_t ]
err = pcwise_U - pcwise_u
@test all( abs.(err) .< 1e-6 )


text:
abstract type AbstractExpectationAlgorithm <: DiffEqBase.DEAlgorithm end
struct Koopman <:AbstractExpectationAlgorithm end
struct MonteCarlo <: AbstractExpectationAlgorithm end

# tuplejoin from https://discourse.julialang.org/t/efficient-tuple-concatenation/5398/8
@inline tuplejoin(x) = x
@inline tuplejoin(x, y) = (x..., y...)
@inline tuplejoin(x, y, z...) = (x..., tuplejoin(y, z...)...)

_rand(x::T) where T <: Sampleable = rand(x)
_rand(x) = x

function __make_map(prob::ODEProblem, args...; kwargs...)
    (u,p) -> solve(remake(prob,u0=u,p=p), args...; kwargs...)
end

function expectation(g::Function, prob::ODEProblem, u0, p, expalg::Koopman, args...;
                        u0_CoV=(u,p)->u, p_CoV=(u,p)->p,
                        maxiters=1000000,
                        batch=0,
                        quadalg=HCubatureJL(),
                        ireltol=1e-2, iabstol=1e-2,
                        nout=1,kwargs...)

    # construct extended state space
    n_states = length(u0)
    n_params = length(p)
    ext_state = [u0; p]
    # ext_state = (u0...,p...)

    # find indices corresponding to distributions, check if sampleable and has non-zero support.
    dist_mask = collect(isa.(ext_state, Sampleable) .& (minimum.(ext_state) .!= maximum.(ext_state)))
    val_mask = .!(dist_mask)

    # get distributions and indx in extended state space
    dists = ext_state[dist_mask]

    # create numerical state space values
    # ext_state_val = minimum.(ext_state)
    T0 = promote_type(eltype.(ext_state)...)
    ext_state_val = [T0(minimum(v)) for v ∈ ext_state] |> collect  #collect needed for zygote for somereason. Otherwise is a tuple

    state_view = @view ext_state_val[dist_mask]
    param_view = @view ext_state_val[val_mask]

    if batch <= 1
        S = __make_map(prob, args...; kwargs...)
        integrand = function (x, p)
            ## Hack to avoid mutating array replacing ext_state_val[ext_state_dist_bitmask] .= x
            x_it = 0
            p_it = 0
            T = promote_type(eltype(x),eltype(p))
            esv = map(1:length(ext_state_val)) do idx
                dist_mask[idx] ? T(x[x_it+=1]) : T(p[p_it+=1])
            end

            _u0 = @view(esv[1:n_states])
            _p = @view(esv[n_states+1:end])

            # Koopman
            w = prod(pdf(a, b) for (a, b) in zip(dists, x))
            Ug = g(S(u0_CoV(_u0,_p), p_CoV(_u0,_p)))

            return Ug*w
        end
    else
        integrand = function (dx, x, p)
            trajectories = size(x, 2)
            T = promote_type(eltype(x),eltype(p))

            prob_func = function (prob, i, repeat)
                x_it = 0
                p_it = 0
                esv = map(1:length(ext_state_val)) do idx
                    dist_mask[idx] ? T(x[x_it+=1,i]) : T(p[p_it+=1])
                end

                u0_view = @view(esv[1:n_states])
                p_view = @view(esv[n_states+1:end])
                remake(prob, u0=u0_CoV(u0_view,p_view),p=p_CoV(u0_view,p_view))
            end

            output_func = function (sol, i)
                w = prod(pdf(a, b) for (a, b) in zip(dists, x[:,i]))
                Ug = g(sol)
                return Ug*w,false #Ug * w, false
            end

            ensembleprob = EnsembleProblem(prob, prob_func=prob_func, output_func=output_func)
            sol = solve(ensembleprob, args...;trajectories=trajectories,kwargs...)
            dx .= hcat(sol.u...) # Why do I need to hcat???
            nothing
        end
    end

    # TODO fix params usage
    lb = minimum.(dists)
    ub = maximum.(dists)

    T = promote_type(eltype(lb),eltype(ub), eltype(ext_state_val))
    intprob = QuadratureProblem(integrand, lb, ub, T.(param_view), batch=batch, nout=nout)
    sol = solve(intprob, quadalg, reltol=ireltol, abstol=iabstol, maxiters=maxiters)
end

function expectation(g::Function, prob::ODEProblem, u0, p, expalg::MonteCarlo, args...;
        trajectories,
        u0_CoV=(u,p)->u, p_CoV=(u,p)->p,
        kwargs...)

    prob_func = function (prob, i, repeat)
        _u0 = _rand.(u0)
        _p = _rand.(p)
        remake(prob, u0=u0_CoV(_u0,_p), p=p_CoV(_u0,_p))
    end

    output_func = (sol, i) -> (g(sol), false)

    monte_prob = EnsembleProblem(prob;
                output_func=output_func,
                prob_func=prob_func)
    sol = solve(monte_prob, args...;trajectories=trajectories,kwargs...)
    mean(sol.u)# , sol
end

"""
    centralmoment(n, g, args...; kwargs) -> [n by 1 Array]

Computes the n central moments of the function g using the Koopman expectation.
The function is a wrapper over expectation, arguments can be piped through with
args and kwargs.

Return: n-length array of the 1 to n central moments

Note: The first central moment is, by definition, always 0

TODO: - add support for vector-valued g functions, currently assumes scalar
      return values.
      - add tests
"""
function centralmoment(n::Int, g::Function, args...; kwargs...)
    if n < 2 return Float64[] end

    # Compute the expectations of g, g^2, ..., g^n
    sol = expectation(x -> [g(x)^i for i in 1:n], args...; nout = n, kwargs...)
    exp_set = sol.u
    mu_g = popfirst!(exp_set)

    # Combine according to binomial expansion
    const_term(n) = (-1)^(n-1) * (n-1) * mu_g^n
    binom_term(n, k, mu, exp_gi) = binomial(n, k) * (-mu)^(n - k) * exp_gi
    binom_sum = function (exp_vals)
        m = length(exp_vals) + 1
        sum([binom_term(m, k + 1, mu_g, v) for (k,v) in enumerate(exp_vals)]) + const_term(m)
    end

    return [zero(exp_set[1]), [binom_sum(exp_set[1:i]) for i in 1:length(exp_set)]...]
end


text:
"""
Embedding operator that randomly samples
between parent's value and the nearest parameter boundary
to get the new valid value if target's parameter is out-of-bounds.
"""
struct RandomBound{S<:RectSearchSpace} <: EmbeddingOperator
    search_space::S

    RandomBound(search_space::S) where {S<:RectSearchSpace} = new{S}(search_space)
end

search_space(rb::RandomBound) = rb.search_space

function apply!(eo::RandomBound, target::AbstractIndividual, ref::AbstractIndividual)
    length(target) == length(ref) == numdims(eo.search_space) ||
        throw(ArgumentError("Dimensions of problem/individuals do not match"))
    ss = search_space(eo)
    ssmins = dimmin(ss)
    ssmaxs = dimmax(ss)

    @inbounds for i in eachindex(target)
        l, u = ssmins[i], ssmaxs[i]

        if target[i] < l
            target[i] = l + rand() * (ref[i]-l)
        elseif target[i] > u
            target[i] = u + rand() * (ref[i]-u)
        else # continuous range doesn't need further checks
            (ss isa MixedPrecisionRectSearchSpace) || continue
        end
        if (ss isa MixedPrecisionRectSearchSpace) && (dimdigits(ss, i) >= 0)
            target[i] = round(target[i], digits=dimdigits(ss, i))
        end
        @assert l <= target[i] <= u "target[$i]=$(target[i]) is out of [$l, $u]"
    end
    return target
end

apply!(eo::RandomBound, target::AbstractIndividual, pop, refIndex::Int) =
    apply!(eo, target, viewer(pop, refIndex))

function apply!(eo::RandomBound, target::AbstractIndividual,
                pop, parentIndices::AbstractVector{Int})
    @assert length(parentIndices) == 1
    apply!(eo, target, pop, parentIndices[1])
end