text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    ShiftedPath(path, offset)

Produces the same results of `path` but shifted by an `offset`
that can be positive or negative.
"""
struct ShiftedPath{D<:AbstractDomain,P<:AbstractPath{D}} <: AbstractPath{D}
  path::P
  offset::Int
  start::Int
  length::Int
end

function ShiftedPath(path::AbstractPath{D}, offset::Int) where {D<:AbstractDomain}
  _, s = iterate(path)
  start = s - 1
  len = length(path)
  off = offset ≥ 0 ? offset : len + offset

  ShiftedPath{D,typeof(path)}(path, off, start, len)
end

function Base.iterate(path::ShiftedPath, state=1)
  if state > path.length
    nothing
  else
    s = ((state + path.offset - 1) % path.length) + path.start
    loc, _ = iterate(path.path, s)
    loc, state + 1
  end
end

Base.length(path::ShiftedPath) = path.length

# ------------
# IO methods
# ------------
function Base.show(io::IO, path::ShiftedPath)
  print(io, "ShiftedPath")
end

function Base.show(io::IO, ::MIME"text/plain", path::ShiftedPath)
  println(io, path)
end


text:
using TomoForward
using Images
using Plots
using XfromProjections
using StaticArrays
using PyCall
using Logging

include("../tv_primaldual_flow.jl")
include("../optical_flow.jl")
include("../../phantoms/sperm_phantom.jl")

H,W = 609, 609
detmin, detmax = -38.0, 38.0
function radon_operator(height, width, detcount, angles)
    proj_geom = ProjGeom(0.5, detcount, angles)
    A = fp_op_parallel2d_line(proj_geom, height, width, detmin,detmax, detmin,detmax)
    return A
end

images, tracks = get_sperm_phantom(301,grid_size=0.1)
frames = images[:,:,collect(271:10:300)]
r(s) = 1.0
w_tv = 0.01
w_flow  = 0.001
c=10.0

grid = collect(detmin:0.1:detmax)

detcount = Int(floor(H*1.4))

cwd = @__DIR__
path = normpath(joinpath(@__DIR__, "result"))
!isdir(path) && mkdir(path)
cd(path)

angle_numbers = [2,4]
for ang_nr in angle_numbers
    angles = collect(range(π/2,3*π/2,length=ang_nr+1))[1:end-1]

    #match size of input image (generating data)
    As = map(t -> radon_operator(size(frames[:,:,1],1),size(frames[:,:,1],2),detcount, angles),1:size(frames)[3])
    bs = zeros(size(As[1])[1],size(frames)[3])
    map(t -> bs[:,t] = As[t]*vec(frames[:,:,t]), 1:size(frames)[3])

    u0s = zeros(H,W,size(frames)[3])

    As  = map(t -> radon_operator(H,W, detcount, angles),1:size(u0s)[3])

    @info "Reconstruction using tv regularization frame by frame"
    niter=450
    us_tv = zeros(H,W,size(frames)[3])
    for t = 1:size(frames)[3]
        A = As[t]
        p = bs[:,t]
        u0 = u0s[:,:,t]
        us_tv[:,:,t] .= recon2d_tv_primaldual!(us_tv[:,:,t], A, p, niter, w_tv, c)
    end

    @info "Reconstructing using joint motion estimation and reconstruction"
    niter1=50
    niter2=20
    u0s = deepcopy(us_tv)
    us_flow = recon2d_tv_primaldual_flow(As, bs, u0s, niter1, niter2, w_tv, w_flow, c)

    p1 = heatmap(grid, grid, Gray.(us_flow[:,:,1]), yflip=false, aspect_ratio=:equal, framestyle=:none, legend=false)
    savefig(@sprintf "test_%d" length(angles))
    plot!(tracks[271][:,1], tracks[271][:,2], aspect_ratio=:equal, linewidth=5)
    p2 = heatmap(grid, grid, Gray.(us_flow[:,:,2]), yflip=false, aspect_ratio=:equal, framestyle=:none, legend=false)
    plot!(tracks[281][:,1], tracks[281][:,2], aspect_ratio=:equal, linewidth=5)
    p3 = heatmap(grid, grid, Gray.(us_flow[:,:,3]), yflip=false, aspect_ratio=:equal, framestyle=:none, legend=false)
    plot!(tracks[291][:,1], tracks[291][:,2], aspect_ratio=:equal, linewidth=5)
    l = @layout [a b c]
    plot(p1, p2, p3, layout = l, size=(2000,600), linewidth=5)
    savefig(@sprintf "result_all_%d" length(angles))
end

cd(cwd)

ext: jl

max_stars_repo_name: QPImaging/xfromProjections

max_issues_repo_name: QPImaging/xfromProjections

max_forks_repo_name: QPImaging/xfromProjections


text:
struct Tube{T} <: AbstractVolumePrimitive{T, 3} ## Only upright Tubes at the moment
    r_interval::AbstractInterval{T}
    φ_interval::AbstractInterval{T}
    z_interval::AbstractInterval{T}
    translate::Union{CartesianVector{T},Missing}
end

function Tube{T}(dict::Dict{Any, Any}, inputunit_dict::Dict{String,Unitful.Units})::Tube{T} where {T <: SSDFloat}

    z_offset::T = 0.0
    if haskey(dict, "translate")
        translate = CartesianVector{T}( haskey(dict["translate"],"x") ? geom_round(ustrip(uconvert(u"m", T(dict["translate"]["x"]) * inputunit_dict["length"] ))) : T(0),
                                        haskey(dict["translate"],"y") ? geom_round(ustrip(uconvert(u"m", T(dict["translate"]["y"]) * inputunit_dict["length"] ))) : T(0),
                                        haskey(dict["translate"],"z") ? geom_round(ustrip(uconvert(u"m", T(dict["translate"]["z"]) * inputunit_dict["length"] ))) : T(0))
        if translate[1] == T(0.0) && translate[2] == T(0.0) && translate[3] == T(0.0)
            translate = missing
        elseif translate[1] == T(0.0) && translate[2] == T(0.0)
            translate = missing
            z_offset = geom_round(ustrip(uconvert(u"m", T(dict["translate"]["z"]) * inputunit_dict["length"] )))
        end
    else
        translate = missing
    end

    if haskey(dict,"h")
         zStart, zStop = z_offset, geom_round(z_offset + ustrip(uconvert(u"m", T(dict["h"]) * inputunit_dict["length"] )))
    elseif haskey(dict,"z")
         zStart, zStop =  geom_round(z_offset + ustrip(uconvert(u"m", T(dict["z"]["from"]) * inputunit_dict["length"] ))),  geom_round( z_offset + ustrip(uconvert(u"m", T(dict["z"]["to"]) * inputunit_dict["length"])))
    else
        @warn "please specify a height of the Tube 'h'"
    end

    phi_interval =  if haskey(dict, "phi")
        Interval(geom_round(T(deg2rad(dict["phi"]["from"]))), geom_round(T(deg2rad(dict["phi"]["to"]))))
    else
        Interval(geom_round(T(deg2rad(0))), geom_round(T(deg2rad(360))))
    end

    r_interval = if haskey(dict["r"], "from")
        Interval(geom_round(ustrip(uconvert(u"m", T(dict["r"]["from"]) * inputunit_dict["length"] ))), geom_round(ustrip(uconvert(u"m", T(dict["r"]["to"]) * inputunit_dict["length"]))))
    else
        Interval(geom_round(0), geom_round(ustrip(uconvert(u"m", T(dict["r"]) * inputunit_dict["length"]))))
    end

    return Tube{T}(
        r_interval,
        phi_interval,
        Interval(zStart, zStop),
        translate
    )
end

function Geometry(T::DataType, t::Val{:tube}, dict::Dict{Union{Any,String}, Any},inputunit_dict::Dict{String,Unitful.Units})
    return Tube{T}(Dict{Any,Any}(dict), inputunit_dict)
end

function in(point::CartesianPoint{T}, tube::Tube{T}) where T
    ismissing(tube.translate) ? nothing : point -= tube.translate
    point = convert(CylindricalPoint,point)
    if point.r in tube.r_interval && point.φ in tube.φ_interval && point.z in tube.z_interval
        return true
    end
    return false
end

function in(point::CylindricalPoint{T}, tube::Tube{T}) where T
    (ismissing(tube.translate) || tube.translate == CartesianVector{T}(0.0,0.0,0.0)) ? nothing  : point = CylindricalPoint(CartesianPoint(point)-tube.translate)
    if point.r in tube.r_interval && point.φ in tube.φ_interval && point.z in tube.z_interval
        return true
    end
    return false
end

function get_important_points(t::Tube{T}, ::Val{:r})::Vector{T} where {T <: SSDFloat}
    return geom_round.(T[t.r_interval.left, t.r_interval.right])
end

function get_important_points(t::Tube{T}, ::Val{:φ})::Vector{T} where {T <: SSDFloat}
    return geom_round.(T[t.φ_interval.left, t.φ_interval.right])
end

function get_important_points(t::Tube{T}, ::Val{:z})::Vector{T} where {T <: SSDFloat}
    return geom_round.(T[t.z_interval.left, t.z_interval.right])
end

function get_important_points(t::Tube{T}, ::Val{:x})::Vector{T} where {T <: SSDFloat}
    return geom_round.(T[-t.r_interval.left, -t.r_interval.right, t.r_interval.left, t.r_interval.right])
end

function get_important_points(t::Tube{T}, ::Val{:y})::Vector{T} where {T <: SSDFloat}
    return geom_round.(T[-t.r_interval.right, -t.r_interval.left, t.r_interval.left, t.r_interval.right])
end

function sample(c::Tube{T}, stepsize::Vector{T}) where T
    samples = CylindricalPoint[]
    for r in c.r_interval.left:stepsize[1]: c.r_interval.right
        for φ in c.φ_interval.left:stepsize[2]: c.φ_interval.right
            for z in c.z_interval.left:stepsize[3]: c.z_interval.right
                push!(samples, CylindricalPoint{T}(r,φ,z))
            end
        end
    end
    ismissing(c.translate) ? nothing : samples = map(x -> CylindricalPoint(CartesianPoint(x) + c.translate), samples)
    return samples
end

function (+)(t::Tube{T}, translate::Union{CartesianVector{T},Missing})::Tube{T} where {T <: SSDFloat}
    if ismissing(translate)
        return t
    elseif ismissing(t.translate)
        return Tube(t.r_interval, t.φ_interval, t.z_interval, translate)
    else
        return Tube(t.r_interval, t.φ_interval, t.z_interval, t.translate + translate)
    end
 end


text:
module TestSubArray

using Base.Test
using NullableArrays
using Compat.view

nd = rand(3:5)
sz = [ rand(3:10) for i in 1:nd ]
X = NullableArray(rand(sz...), rand(Bool, sz...))

I = [ rand(1:sz[i]) for i in 1:nd ]
H = [ rand(1:sz[i]) for i in 1:nd ]

for i in 1:nd
    J = [ (x->x==i ? I[x] : Colon())(j) for j in 1:nd ]
    S = view(X, J...)
    H = [ (x->x==i ? I[x] : rand(1:sz[x]))(j) for j in 1:nd ]
    _H = H[find(x->x!=i, collect(1:nd))]

    @test values(S, _H...) == X.values[H...]
    @test isnull(S, _H...) == X.isnull[H...]
    @test anynull(S) == anynull(X[J...])
end

end # module TestSubArray


text:
"""
    AbstractGPKernel

Abstract kernel type. User-defined kernels should derive from it.

Implementations have to provide methods for [`get_hyperparameters_size`](@ref)
and [`covariance`](@ref). Methods for [`covariance_training`](@ref), [`covariance_diagonal`](@ref)
and [`covariance_grad`](@ref) are optional.
"""
abstract type AbstractGPKernel end

"""
    covariance(ker::AbstractGPKernel, log_theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2}, z::AbstractArray{Float64, 2})

Return the covariance matrix. Should be overridden by kernel implementations.

# Arguments
- `ker`: The kernel object. Implementations must override with their own subtype.
- `log_theta`: natural logarithm of hyperparameters.
- `x, z`: Input data, reshaped into 2-d arrays.
  `x` must have dimensions ``n \\times d``; `z` must have dimensions ``m \\times d``.

# Return
The covariance matrix, of size ``n \\times m``.
"""
function covariance(ker::AbstractGPKernel, log_theta::AbstractArray{Float64, 1},
        x1::AbstractArray{Float64, 2}, x2::AbstractArray{Float64, 2}) end

"""
    covariance_training(ker::AbstractGPKernel, log_theta::AbstractArray{Float64, 1},
        training_x::AbstractArray{Float64, 2})

This is a speedup version of [`covariance`](@ref), which is only called during
traing sequence. Intermediate matrices computed
in this function for particular hyperparameters can be cached and reused subsequently, either
in this function or in [`covariance_grad`](@ref)

Default method just delegates to [`covariance`](@ref) with `x === z`. Kernel implementations can
optionally override it for betrer performance.

See [`covariance`](@ref) for description of arguments and return values.
"""
function covariance_training(ker::AbstractGPKernel, log_theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2})
    covariance(ker, log_theta, x, x)
end

"""
    covariance_diagonal(ker::AbstractGPKernel, log_theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2})

This is a speedup version of [`covariance`](@ref), which is invoked if the caller
is not interested in the entire covariance matrix, but only needs the variance, i.e. the
diagonal of the covariance matrix.

Default method just returns `diag(covariance(...))`, with `x === z`. Kernel implementations can
optionally override it to achieve betrer performance, by not computing the non diagonal
elements of covariance matrix.

See [`covariance`](@ref) for description of arguments.

# Return
The 1-d array of variances, of size `size(x, 1)`.
"""
function covariance_diagonal(ker::AbstractGPKernel, log_theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2})
    diag(covariance(ker, log_theta, x, x))
end

"""
    covariance_grad(ker::AbstractGPKernel, log_theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2}, R::AbstractArray{Float64, 2})

Return the gradient of the covariance function with respect to logarigthms of hyperparameters,
based on the provided direction matrix.

This function can be optionally overridden by kernel implementations.
If the gradient function is not provided, [`gp_train`](@ref) will fail back to
`NelderMead` algorithm by default.

# Arguments
- `ker`: The kernel object. Implementations must override with their own subtype.
- `log_theta`:  natural logarithm of hyperparameters
- `x`: Training data, reshaped into a 2-d array.
  `x` must have dimensions ``n \\times d``.
- `R` the directional matrix, ``n \\times n``
```math
R = \\frac{1}{\\sigma_n^2}(\\alpha * \\alpha^T - K^{-1}); \\alpha = K^{-1}y
```

# Return
A vector of size `length(log_theta)`, whose ``j``'th element is equal to
```math
tr(R \\frac{\\partial K}{\\partial \\eta_j})
```
"""
function covariance_grad(ker::AbstractGPKernel, theta::AbstractArray{Float64, 1},
        x::AbstractArray{Float64, 2}, R::AbstractArray{Float64, 2})
    :covariance_gradient_not_implemented
end

"""
    get_hyperparameters_size(kernel::AbstractGPKernel, training_data::AbstractArray{Float64, 2})

Return the number of hyperparameters for used by this kernel on this training data set.
Should be overridden by kernel implementations.
"""
function get_hyperparameters_size(kernel::AbstractGPKernel, training_data::AbstractArray{Float64, 2}) end


text:
"""
Structures and functions for the GENERATIVE POMDPs.jl interface of the Pricing MDP
"""

"""
m = PMDPg(edges, products, λ)

PMDP for generative interface
"""
struct PMDPg <: PMDP{State, Action}
    pp::PMDPProblem                # Pricing Problem
    empty_product::Product
    empty_product_id::Int64

    function PMDPg(pp::PMDPProblem)
        new(pp, empty_product(pp), n_products(pp)+1)
    end

    # function PMDPg(E, P, λ, B, A, objective)
    #     # selling_period_ends = get_product_selling_period_ends(E, P)
    #     T = selling_period_ends[1]
    #     empty_product=P[1]
    #     @assert objective in [:revenue, :utilization]
    #     pi = productindices(P)
    #     return new(length(E), T,E,P,λ, selling_period_ends, empty_product,B, A, objective, pi)
    # end
end

"""
Returns next requested product. If in given timestep one of the prodcuts has selling period end, update the product request probs.

TODO: Potential speedup if product_request_probs are not recalculated at every step
"""
function sample_request(m::PMDPg, t::Timestep, rng::AbstractRNG)::Int64
    iₚ = rand(rng, demand(m)[t])
    # prod_index == n_products(pp(m))+1 ? p = empty_product(m) : p = products(m)[prod_index]
    # return p
end

"""
Sample user budget Budget for product requested in state s.
"""
function sample_customer_budget(m::PMDPg, s::State, rng::AbstractRNG)::Action
    # local b::Float64
    if s.iₚ != m.empty_product_id
        budget_distribution = budgets(m)[s.iₚ]
        budget = rand(rng, budget_distribution)
    else
        budget = EMPTY_PRODUCT_USER_BUDGET
    end
    return budget
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export button

"""
    button(;kwargs...)
    button(children::Any;kwargs...)
    button(children_maker::Function;kwargs...)

A Button component.
Render button or link with button styles from mantine theme. For more information, see: https://mantine.dev/core/button/
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): Primary content
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `className` (String; optional): Often used with CSS to style elements with common properties
- `color` (a value equal to: "dark", "gray", "red", "pink", "grape", "violet", "indigo", "blue", "cyan", "teal", "green", "lime", "yellow", "orange"; optional): Button color from theme
- `compact` (Bool; optional): Reduces vertical and horizontal spacing
- `disabled` (Bool; optional): The component can show it is currently unable to be interacted with
- `fullWidth` (Bool; optional): Sets button width to 100% of parent element
- `gradient` (optional): Controls gradient settings in gradient variant only. gradient has the following type: lists containing elements 'from', 'to', 'deg'.
Those elements have the following types:
  - `from` (String; required)
  - `to` (String; required)
  - `deg` (Real; optional)
- `loaderPosition` (a value equal to: "left", "right"; optional): Loader position relative to button label
- `loading` (Bool; optional): Indicate loading state
- `loading_state` (optional): Object that holds the loading state object coming from dash-renderer. loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `n_clicks` (Real; optional): An integer that represents the number of times that this element has been clicked on
- `radius` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional): Button border-radius from theme or number to set border-radius in px
- `size` (a value equal to: "xs", "sm", "md", "lg", "xl"; optional): Predefined button size
- `style` (Dict; optional): Inline style override
- `uppercase` (Bool; optional): Set text-transform to uppercase
- `variant` (a value equal to: "link", "filled", "outline", "light", "gradient", "white", "default"; optional): Controls button appearance
"""
function button(; kwargs...)
        available_props = Symbol[:children, :id, :className, :color, :compact, :disabled, :fullWidth, :gradient, :loaderPosition, :loading, :loading_state, :n_clicks, :radius, :size, :style, :uppercase, :variant]
        wild_props = Symbol[]
        return Component("button", "Button", "dash_mantine_components", available_props, wild_props; kwargs...)
end

button(children::Any; kwargs...) = button(;kwargs..., children = children)
button(children_maker::Function; kwargs...) = button(children_maker(); kwargs...)


text:
using Flux

weights, re = Flux.destructure(Chain(Dense(10, 10, relu), Dense(10, 5)))

# NN = re(fill(1 ./ 1f5, length(weights)))
# NN = re(randn(Float32, length(weights)) ./ 1f5)

NN = re(weights ./ 1f4)
data = [(randn(Float32, 10), randn(Float32, 5)) for i in 1:50]

loss(x, y) = Flux.mse(NN(x), y)

opt = ADAM(1e-4)

@show params(NN)
mean([loss(point[1], point[2]) for point in data])
Flux.train!(loss, Flux.params(NN), data, opt)
mean([loss(point[1], point[2]) for point in data])

weights, _ = Flux.destructure(NN)
W_1, b_1, W_2, b_2 = params(NN)

display(W_1)
display(b_1)
display(W_2)
display(b_2)

NN(rand(10)) .- NN(rand(10))


text:
using PenaltiesLH
using Test

# Scalar penalty functions can accept multiple arguments.
scalar_test_fct(x, y) = sum(x) + sum(y);

function make_test_penalty(pName :: Symbol; isUsed :: Bool = true)
    return Penalty(name = pName,  isUsed = isUsed,
        scalarFct = scalar_test_fct, scalarValue = 0.0);
end

function make_test_penalty_vector(n :: Integer)
    pv = PenaltyVector(n);
    for j = 1 : n
        pv[j] = make_test_penalty(Symbol("p$j"), isUsed = iseven(j));
    end
    return pv
end

function penalty_test()
    @testset "Penalty" begin
        p = make_test_penalty(:p1);
        @test name(p) == :p1
        @test scalar_value(p) == 0.0
        set_scalar_value(p, 2.0)
        @test scalar_value(p) == 2.0
        show(p)

        x = [1.0, 2.0];
        y = [4.3, 2.9];
        # Need to package multiple arguments into a tuple
        sv = compute_penalty(p, (x,y));
        @test scalar_value(p) == sv
    end
end

function pvector_test()
    @testset "PenaltyVector" begin
        n = 7;
        pv = make_test_penalty_vector(n);
        @test length(pv) == n
        @test !isempty(pv)
        for j = 1 : n
            p = pv[j];
            @test isa(p, Penalty)
        end

        j = 0
        for p in pv
            j += 1
            @test isa(p, Penalty)
        end
        @test j == n

        p = retrieve(pv, :p3);
        @test name(p) == :p3
        p = retrieve(pv, :pp3);
        @test isnothing(p)

        pName = :pPushed;
        p7 = make_test_penalty(pName);
        push!(pv, p7);
        @test length(pv) == n+1
        p = retrieve(pv, pName);
        @test name(p) == pName

        d = penalty_dict(pv);
        for p in pv
            if PenaltiesLH.is_used(p)
                @test haskey(d, name(p))
            else
                @test !haskey(d, name(p))
            end
        end

        use(pv, :p3);
        @test is_used(pv, :p3)
        no_use(pv, :p3)
        @test !is_used(pv, :p3)

        # Computing penalties
        x = [3.4, 4.5];
        y = [9.2, 0.8];
        pTotal = total_penalty(pv, data = (x, y));
        compute_all_penalties(pv, (x, y));
        pTotal2 = total_penalty(pv, data = (x, y));
        @test pTotal ≈ pTotal2

        report_penalties(pv)
    end
end

@testset "PenaltiesLH" begin
    penalty_test()
    pvector_test()
end

# ---------


text:
"""```
measurement(m::myModel1010depo{T}, TTT::Matrix{T}, RRR::Matrix{T},
            CCC::Vector{T}) where {T<:AbstractFloat}
```

Assign measurement equation

```
y_t = ZZ*s_t + DD + u_t
```

where

```
Var(ϵ_t) = QQ
Var(u_t) = EE
Cov(ϵ_t, u_t) = 0
```
"""
function measurement(m::myModel1010depo{T},
                     TTT::Matrix{T},
                     RRR::Matrix{T},
                     CCC::Vector{T}) where {T<:AbstractFloat}
    endo     = m.endogenous_states
    endo_new = m.endogenous_states_augmented
    exo      = m.exogenous_shocks
    obs      = m.observables

    # If shocks = true, then return measurement equation matrices with rows and columns for
    # anticipated policy shocks
    _n_observables = n_observables(m)
    _n_states = n_states_augmented(m)
    _n_shocks_exogenous = n_shocks_exogenous(m)

    ZZ = zeros(_n_observables, _n_states)
    DD = zeros(_n_observables)
    EE = zeros(_n_observables, _n_observables)
    QQ = zeros(_n_shocks_exogenous, _n_shocks_exogenous)

    ## GDP growth - Quarterly!
    ZZ[obs[:obs_gdp], endo[:y_t]]          = 1.0
    ZZ[obs[:obs_gdp], endo_new[:y_t1]]     = -1.0
    ZZ[obs[:obs_gdp], endo[:z_t]]          = 1.0
    ZZ[obs[:obs_gdp], endo_new[:e_gdp_t]]  = 1.0
    ZZ[obs[:obs_gdp], endo_new[:e_gdp_t1]] = -m[:me_level]
    DD[obs[:obs_gdp]]                      = 100*(exp(m[:z_star])-1)

    ## GDI growth- Quarterly!
    ZZ[obs[:obs_gdi], endo[:y_t]]          = m[:γ_gdi]
    ZZ[obs[:obs_gdi], endo_new[:y_t1]]     = -m[:γ_gdi]
    ZZ[obs[:obs_gdi], endo[:z_t]]          = m[:γ_gdi]
    ZZ[obs[:obs_gdi], endo_new[:e_gdi_t]]  = 1.0
    ZZ[obs[:obs_gdi], endo_new[:e_gdi_t1]] = -m[:me_level]
    DD[obs[:obs_gdi]]                      = 100*(exp(m[:z_star])-1) + m[:δ_gdi]

    ## Hours growth
    ZZ[obs[:obs_hours], endo[:L_t]] = 1.0
    DD[obs[:obs_hours]]             = m[:Lmean]

    ## Labor Share/real wage growth
    ZZ[obs[:obs_wages], endo[:w_t]]      = 1.0
    ZZ[obs[:obs_wages], endo_new[:w_t1]] = -1.0
    ZZ[obs[:obs_wages], endo[:z_t]]      = 1.0
    DD[obs[:obs_wages]]                  = 100*(exp(m[:z_star])-1)

    ## Inflation (GDP Deflator)
    ZZ[obs[:obs_gdpdeflator], endo[:π_t]]            = m[:Γ_gdpdef]
    ZZ[obs[:obs_gdpdeflator], endo_new[:e_gdpdef_t]] = 1.0
    DD[obs[:obs_gdpdeflator]]                        = 100*(m[:π_star]-1) + m[:δ_gdpdef]

    ## Inflation (Core PCE)
    ZZ[obs[:obs_corepce], endo[:π_t]]             = 1.0
    ZZ[obs[:obs_corepce], endo_new[:e_corepce_t]] = 1.0
    DD[obs[:obs_corepce]]                         = 100*(m[:π_star]-1)

    ## Nominal interest rate
    ZZ[obs[:obs_nominalrate], endo[:R_t]] = 1.0
    DD[obs[:obs_nominalrate]]             = m[:Rstarn]

    ## Nominal deposit rate
    if get_setting(m, :add_deposits)
        ZZ[obs[:obs_depositrate], endo[:Rd_t]] = 1.0
        DD[obs[:obs_depositrate]]              = m[:Rstard]
    end

    ## Consumption Growth
    ZZ[obs[:obs_consumption], endo[:c_t]]      = 1.0
    ZZ[obs[:obs_consumption], endo_new[:c_t1]] = -1.0
    ZZ[obs[:obs_consumption], endo[:z_t]]      = 1.0
    DD[obs[:obs_consumption]]                  = 100*(exp(m[:z_star])-1)

    ## Investment Growth
    ZZ[obs[:obs_investment], endo[:i_t]]      = 1.0
    ZZ[obs[:obs_investment], endo_new[:i_t1]] = -1.0
    ZZ[obs[:obs_investment], endo[:z_t]]      = 1.0
    DD[obs[:obs_investment]]                  = 100*(exp(m[:z_star])-1)

    ## 20 yrs forward transition matrix
    TTT20 = if subspec(m) == "ss16"
        eye(size(TTT,1))
    else
        (1/80)*((UniformScaling(1.) - TTT) \ (UniformScaling(1.) - TTT^80))
    end
    ## BAA or BBB Spread
    ZZ[obs[:obs_BBBspread], endo[:ERtil_k_t]]    = 1.0
    ZZ[obs[:obs_BBBspread], endo[:R_t]]          = -1.0
    # the previous 2 lines create a row of ZZ which is a vector of
    # zeros except for the coefficients on the relevant states. We
    # multiply this by TTT20 to adjust for the 20-year maturity of the
    # bond, and reassign the result to the BBB spread row of the ZZ
    # matrix, overwriting the original coefficients. Then we add a
    # coefficient for the measurement error.
    ZZ[obs[:obs_BBBspread], :]                   = ZZ[obs[:obs_BBBspread], :]' * TTT20
    ZZ[obs[:obs_BBBspread], endo_new[:e_BBB_t]]  = 1.0
    DD[obs[:obs_BBBspread]]                      = 100*log(m[:spr]*m[:lnb_safe]*m[:lnb_liq])

    ## AAA Spread
    ZZ[obs[:obs_AAAspread], endo[:b_liq_t]]       = -(1.0-m[:λ_AAA])*(m[:σ_c]*(1 + m[:h]*exp(-m[:z_star])))/(1 - m[:h]*exp(-m[:z_star]))
    ZZ[obs[:obs_AAAspread], endo[:ERtil_k_t]]     = m[:λ_AAA]
    ZZ[obs[:obs_AAAspread], endo[:R_t]]           = -m[:λ_AAA]
    # the previous 3 lines create a row of ZZ which is a vector of
    # zeros except for the coefficients on the relevant states. We
    # multiply this by TTT20 to adjust for the 20-year maturity of the
    # bond, and reassign the result to the AAA spread row of the ZZ
    # matrix, overwriting the original coefficients. Then we add a
    # coefficient for the measurement error.
    # λ_AAA is fixed to 0 in this model.
    ZZ[obs[:obs_AAAspread], :]                    = ZZ[obs[:obs_AAAspread], :]' * TTT20
    ZZ[obs[:obs_AAAspread], endo_new[:e_AAA_t]]   = 1.0
    DD[obs[:obs_AAAspread]]                       = 100*log(m[:lnb_liq]) + 100*m[:λ_AAA]*log(m[:spr]*m[:lnb_safe])

    ## 10 yrs infl exp
    TTT10                          = (1/40)*((UniformScaling(1.) - TTT)\(UniformScaling(1.) - TTT^40))
    ZZ[obs[:obs_longinflation], :] = TTT10[endo[:π_t], :]
    DD[obs[:obs_longinflation]]    = 100*(m[:π_star]-1)

    ## Long Rate
    ZZ[obs[:obs_longrate], :]               = ZZ[obs[:obs_nominalrate], :]' * TTT10
    ZZ[obs[:obs_longrate], endo_new[:lr_t]] = 1.0
    DD[obs[:obs_longrate]]                  = m[:Rstarn]

    ## TFP
    ZZ[obs[:obs_tfp], endo[:z_t]]       = (1-m[:α])*m[:Iendoα] + 1*(1-m[:Iendoα])
    ZZ[obs[:obs_tfp], endo_new[:tfp_t]] = 1.0
    ZZ[obs[:obs_tfp], endo[:u_t]]       = m[:α]/( (1-m[:α])*(1-m[:Iendoα]) + 1*m[:Iendoα] )
    ZZ[obs[:obs_tfp], endo_new[:u_t1]]  = -(m[:α]/( (1-m[:α])*(1-m[:Iendoα]) + 1*m[:Iendoα]) )

    QQ[exo[:g_sh], exo[:g_sh]]                  = m[:σ_g]^2
    QQ[exo[:b_liqtil_sh], exo[:b_liqtil_sh]]    = m[:σ_b_liqtil]^2
    QQ[exo[:b_liqp_sh], exo[:b_liqp_sh]]        = m[:σ_b_liqp]^2
    QQ[exo[:b_safetil_sh], exo[:b_safetil_sh]]  = m[:σ_b_safetil]^2
    QQ[exo[:b_safep_sh], exo[:b_safep_sh]]      = m[:σ_b_safep]^2
    QQ[exo[:μ_sh], exo[:μ_sh]]                  = m[:σ_μ]^2
    QQ[exo[:z_sh], exo[:z_sh]]                  = m[:σ_z]^2
    QQ[exo[:λ_f_sh], exo[:λ_f_sh]]              = m[:σ_λ_f]^2
    QQ[exo[:λ_w_sh], exo[:λ_w_sh]]              = m[:σ_λ_w]^2
    QQ[exo[:rm_sh], exo[:rm_sh]]                = m[:σ_r_m]^2
    QQ[exo[:σ_ω_sh], exo[:σ_ω_sh]]              = m[:σ_σ_ω]^2
    QQ[exo[:μ_e_sh], exo[:μ_e_sh]]              = m[:σ_μ_e]^2
    QQ[exo[:γ_sh], exo[:γ_sh]]                  = m[:σ_γ]^2
    QQ[exo[:π_star_sh], exo[:π_star_sh]]        = m[:σ_π_star]^2
    QQ[exo[:lr_sh], exo[:lr_sh]]                = m[:σ_lr]^2
    QQ[exo[:zp_sh], exo[:zp_sh]]                = m[:σ_z_p]^2
    QQ[exo[:tfp_sh], exo[:tfp_sh]]              = m[:σ_tfp]^2
    QQ[exo[:gdpdef_sh], exo[:gdpdef_sh]]        = m[:σ_gdpdef]^2
    QQ[exo[:corepce_sh], exo[:corepce_sh]]      = m[:σ_corepce]^2
    QQ[exo[:gdp_sh], exo[:gdp_sh]]              = m[:σ_gdp]^2
    QQ[exo[:gdi_sh], exo[:gdi_sh]]              = m[:σ_gdi]^2
    QQ[exo[:AAA_sh], exo[:AAA_sh]]              = m[:σ_AAA]^2
    QQ[exo[:BBB_sh], exo[:BBB_sh]]              = m[:σ_BBB]^2

    # These lines set the standard deviations for the anticipated shocks. They
    # are here no longer calibrated to the std dev of contemporaneous shocks,
    # as we had in 904
    for i = 1:n_anticipated_shocks(m)
        ZZ[obs[Symbol("obs_nominalrate$i")], :]              = ZZ[obs[:obs_nominalrate], :]' * (TTT^i)
        DD[obs[Symbol("obs_nominalrate$i")]]                 = m[:Rstarn]

        if subspec(m) == "ss21"
            QQ[exo[Symbol("rm_shl$i")], exo[Symbol("rm_shl$i")]] = m[Symbol("σ_r_m")]^2 /
            n_anticipated_shocks(m)
        else
            QQ[exo[Symbol("rm_shl$i")], exo[Symbol("rm_shl$i")]] = m[Symbol("σ_r_m$i")]^2
        end
    end

    # Adjustment to DD because measurement equation assumes CCC is the zero vector
    if any(CCC .!= 0)
        DD += ZZ*((UniformScaling(1) - TTT)\CCC)
    end

    return Measurement(ZZ, DD, QQ, EE)
end


text:
import Base: mod

function mod{D, T <: Integer}(i :: NTuple{D, T}, n ::NTuple{D, T}, offset ::Integer = 0)
  (mod(x - offset, y) + offset for (x, y) in zip(i, n))
end

function mod{T1 <: Integer, T2 <: Integer, T3 <: Integer}(i ::T1, n ::T2, offset ::T3 = 0)
  mod(i - offset, n) + offset
end

function within{D, I1 <: Integer, I2 <: Integer}(idx  :: NTuple{D, I1},
                                                 size :: NTuple{D, I2})
  for (i, n) in zip(idx, size)
    if !(1 <= i <= n)
      return false
    end
  end
  return true
end


text:
# Here we take the data from the photodynamic chains
# and write it to a text file so it can be read in
# Python:

using JLD2

@load "../../../data/T1_pd_MCMC_run_000_001.jld2"
chain001 = chain
flatchain001 = flatchain
@load "../../../data/T1_pd_MCMC_run_000_002.jld2"
chain002 = chain
flatchain002 = flatchain
@load "../../../data/T1_pd_MCMC_run_000_003.jld2"
chain003 = chain
flatchain003 = flatchain
# (This doesn't have variable contamination factor).

#istart = 15000; iend = 500000
istart = 1001; iend = 200000
flatchain = hcat(flatchain001[:,istart:iend],flatchain002[:,istart:iend],flatchain003[:,istart:iend])
nsamples = size(flatchain)[2]

using DelimitedFiles
writedlm("../../../data/T1_photdyn_chain_noprior.txt",flatchain')


text:
export quickstart, size, drop_monomorphic, drop_monomorphic!

## experimental and not exported or documented!
function adjacency_matrix(data::PopData)
    data_loci = groupby(data.loci, :locus)
    out_vec = Vector{Array{Int8,2}}(undef, length(data_loci))
    for (j,i) in enumerate(data_loci)
        uniq = unique(skipmissing(i.genotype))
        adj_mat = fill(Int8(0), length(samples(data)), length(uniq))
        for (j,k) in zip(i.genotype, eachrow(adj_mat))
            k .= Ref(j) .=== uniq
        end
        out_vec[j] = adj_mat
    end
    return out_vec
end

#TODO change location in API docs and rename allele_pool?
#TODO replace alleles with universal Symbol type?
"""
    alleles(locus::T) where T<:GenoArray
Return an array of all the non-missing alleles of a locus.
"""
@inline function alleles(locus::T) where T<:GenoArray
    if all(ismissing.(locus))
        return Vector{Union{Missing, eltype(locus).b.types[1]}}(undef, length(locus))
    end
    alle_out = Base.Iterators.flatten(skipmissing(locus)) |> collect
end

"""
    alleles(locus::T, miss::Bool = false) where T<:GenoArray
Return an array of all the non-missing alleles of a locus. Use the second positional
argument as `true` to include missing values.
"""
@inline function alleles(locus::T, miss::Bool) where T<:GenoArray
    int_type = eltype(locus).b.types[1]
    if all(ismissing.(locus))
        return Vector{Union{Missing, int_type}}(undef, length(locus))
    end
    alle_out = Vector{Union{Missing, int_type}}(Base.Iterators.flatten(skipmissing(locus)) |> collect)
    if miss == true
        append!(alle_out, locus[locus .=== missing])
    end
    return alle_out
end

# TODO add to docs: Utils.jl API
function Base.size(data::PopData)
    return (samples = size(data.meta)[1], loci = length(loci(data)))
end

"""
    unique_alleles(locus::T) where T<:GenoArray
Return an array of all the unique non-missing alleles of a locus.
"""
@inline function unique_alleles(locus::GenoArray)
    unique(alleles(locus))
end

"""
    convert_coord(coordinate::String)
Takes non-decimal-degree format as a `String` and returns it as a decimal degree
`Float32`. Can be broadcasted over an array of coordinate strings to convert them.
## Formatting requirements
- Coordinates as a `String` separated by spaces (`"11 43 41"`) or colons (`"11:43:41"`)
- Must use negative sign (`"-11 43.52"`) or single-letter cardinal direction (`"11 43.52W"`)
- Missing data should be coded as the string `"missing"` (can be accomplished with `replace!()`)
- Can mix colons and spaces (although it's bad practice)
### Example
```
julia> convert_coord("-41 31.52")
-41.5253f0

julia> convert_coord.(["-41 31.52", "25 11:54S"])
2-element Array{Float32,1}:
-41.5253
-25.1983
```
"""
function convert_coord(coordinate::String)
    lowercase(coordinate) == "missing" && return missing
    coord_strip = replace(uppercase(coordinate), r"[NSEW]" => "")
    split_coord = parse.(Float32, split(coord_strip, r"\s|:"))
    split_coord[2] /=60.0
    if length(split_coord) == 3
        split_coord[3] /=3600.0
    end
    conv = sum(abs.(split_coord))
    # N + E are positive | S + W are negative
    if split_coord[1] < 0 || occursin(r"[SW]", uppercase(coordinate))
        # negative
        return round(conv * -1, digits = 4)
    else
        # positive
        return round(conv, digits = 4)
    end
end

function Base.copy(data::PopData)
    PopData(copy.([data.meta,data.loci])...)
end

"""
    drop_monomorphic(data::PopData)
Return a `PopData` object omitting any monomorphic loci. Will inform you which
loci were removed.
"""
function drop_monomorphic(data::PopData)
    rm_loci = Vector{String}()
    for (loc, loc_sdf) in pairs(groupby(data.loci, :locus))
        length(unique(skipmissing(loc_sdf[:, :genotype]))) == 1 && push!(rm_loci, loc.locus)
    end

    if length(rm_loci) == 0
        return data
    elseif length(rm_loci) == 1
        @info "Dropping monomorphic locus " * rm_loci[1]
    else
        @info "Dropping $(length(rm_loci)) monomorphic loci" * "\n $rm_loci"
    end
    exclude(data, locus = rm_loci)
end

"""
    drop_monomorphic!(data::PopData)
Edit a `PopData` object in place by omitting any monomorphic loci. Will inform you which
loci were removed.
"""
function drop_monomorphic!(data::PopData)
    rm_loci = Vector{String}()
    for (loc, loc_sdf) in pairs(groupby(data.loci, :locus))
        length(unique(skipmissing(loc_sdf[:, :genotype]))) == 1 && push!(rm_loci, loc.locus)
    end
    if length(rm_loci) == 0
        return data
    elseif length(rm_loci) == 1
        @info "Dropping monomorphic locus " * rm_loci[1]
    else
        @info "Dropping $(length(rm_loci)) monomorphic loci" * "\n $rm_loci"
    end
    exclude!(data, locus = rm_loci)
end

##TODO add to docs API
"""
    loci_dataframe(data::PopData)
Return a wide `DataFrame` of samples as columns, ommitting population information.

**Example**
```
julia> loci_dataframe(@nancycats)
9×237 DataFrame. Omitted printing of 232 columns
│ Row │ N215       │ N216       │ N217       │ N218       │ N219       │
│     │ Tuple…?    │ Tuple…?    │ Tuple…?    │ Tuple…?    │ Tuple…?    │
├─────┼────────────┼────────────┼────────────┼────────────┼────────────┤
│ 1   │ missing    │ missing    │ (135, 143) │ (133, 135) │ (133, 135) │
│ 2   │ (136, 146) │ (146, 146) │ (136, 146) │ (138, 138) │ (140, 146) │
│ 3   │ (139, 139) │ (139, 145) │ (141, 141) │ (139, 141) │ (141, 145) │
│ 4   │ (116, 120) │ (120, 126) │ (116, 116) │ (116, 126) │ (126, 126) │
│ 5   │ (156, 156) │ (156, 156) │ (152, 156) │ (150, 150) │ (152, 152) │
│ 6   │ (142, 148) │ (142, 148) │ (142, 142) │ (142, 148) │ (142, 148) │
│ 7   │ (199, 199) │ (185, 199) │ (197, 197) │ (199, 199) │ (193, 199) │
│ 8   │ (113, 113) │ (113, 113) │ (113, 113) │ (91, 105)  │ (113, 113) │
│ 9   │ (208, 208) │ (208, 208) │ (210, 210) │ (208, 208) │ (208, 208) │
```
"""
function loci_dataframe(data::PopData)
    unstack(select(data.loci, Not(:population)), :name, :genotype)[:, Not(:locus)]
end

#TODO add to docs API
#TODO make a SMatrix instead?
"""
    loci_matrix(data::PopData)
Return a matrix of genotypes with dimensions `samples × loci`.
Rows are samples and columns are loci. Will return an error if ploidy varies between samples.

**Example**
```
julia> loci_matrix(@nancycats)
237×9 Array{Union{Missing, Tuple{Int16,Int16}},2}:
 missing     (136, 146)  (139, 139)  …  (199, 199)  (113, 113)  (208, 208)
 missing     (146, 146)  (139, 145)     (185, 199)  (113, 113)  (208, 208)
 (135, 143)  (136, 146)  (141, 141)     (197, 197)  (113, 113)  (210, 210)
 (133, 135)  (138, 138)  (139, 141)     (199, 199)  (91, 105)   (208, 208)
 (133, 135)  (140, 146)  (141, 145)     (193, 199)  (113, 113)  (208, 208)
 (135, 143)  (136, 146)  (145, 149)  …  (193, 195)  (91, 113)   (208, 208)
 (135, 135)  (136, 146)  (139, 145)     (199, 199)  (105, 113)  (208, 208)
 (135, 143)  (136, 146)  (135, 149)     (193, 197)  (91, 91)    (208, 212)
 (137, 143)  (136, 146)  (139, 139)     (197, 197)  (105, 113)  (208, 212)
 (135, 135)  (132, 132)  (141, 145)     (197, 197)  (91, 105)   (208, 208)
 (137, 141)  (130, 136)  (137, 145)  …  (193, 199)  (91, 91)    (182, 182)
 (129, 133)  (130, 136)  (135, 145)     (193, 199)  (91, 113)   (182, 208)
 ⋮                                   ⋱
 (133, 135)  (136, 136)  (135, 139)  …  (199, 199)  (113, 113)  (182, 182)
 (133, 141)  (136, 136)  (135, 139)     (197, 197)  (113, 113)  (182, 208)
 (133, 141)  (130, 146)  (141, 141)     (191, 199)  missing     (208, 208)
 (123, 133)  (138, 138)  (141, 145)     (191, 197)  missing     (208, 208)
 (123, 133)  (138, 138)  (139, 139)     (197, 199)  missing     (208, 208)
 (133, 141)  (136, 146)  (139, 139)  …  (197, 197)  missing     (208, 208)
 (133, 141)  (130, 136)  (139, 145)     (191, 199)  missing     (208, 208)
 (133, 141)  (136, 146)  (139, 145)     (199, 199)  missing     (208, 220)
 (133, 143)  (130, 130)  (135, 145)     (197, 197)  missing     (208, 208)
 (135, 141)  (136, 144)  (143, 143)     (191, 197)  (113, 117)  (208, 208)
 (137, 143)  (130, 136)  (135, 145)  …  (193, 199)  (113, 117)  (208, 208)
 (135, 141)  (130, 146)  (135, 139)     (197, 197)  missing     (208, 208)
 ```
"""
function loci_matrix(data::PopData)
    dims = size(data)
    sort_df = issorted(data.loci, [:name, :locus]) ? sort(data.loci, [:name, :locus]) : data.loci
    reshape(sort_df.genotype, (dims.samples, dims.loci)) |> collect
end

"""
    multitest_missing(pvals::Vector{T}, method::String) where T <: Union{Missing, <:AbstractFloat}
Modification to `MultipleTesting.adjust` to include `missing` values in the
returned array. See MultipleTesting.jl docs for full more detailed information.

**Example**
```
julia> multitest_missing([0.1, 0.01, 0.005, 0.3], "bh")

````

### `correction` methods (case insensitive)
- `"bonferroni"` : Bonferroni adjustment
- `"holm"` : Holm adjustment
- `"hochberg"` : Hochberg adjustment
- `"bh"` : Benjamini-Hochberg adjustment
- `"by"` : Benjamini-Yekutieli adjustment
- `"bl"` : Benjamini-Liu adjustment
- `"hommel"` : Hommel adjustment
- `"sidak"` : Šidák adjustment
- `"forwardstop"` or `"fs"` : Forward-Stop adjustment
- `"bc"` : Barber-Candès adjustment
"""
@inline function multitest_missing(pvals::Vector{T}, method::String) where T <: Union{Missing, <:AbstractFloat}
    # make a dict of all possible tests and their respective functions
    d = Dict(
        "bonferroni" => Bonferroni(),
        "holm" => Holm(),
        "hochberg" => Hochberg(),
        "bh" => BenjaminiHochberg(),
        "by" => BenjaminiYekutieli(),
        "bl" => BenjaminiLiu(),
        "hommel" => Hommel(),
        "sidak" => Sidak(),
        "forwardstop" => ForwardStop(),
        "fs" => ForwardStop(),
        "bc" => BarberCandes(),
    )
    p_copy = copy(pvals)
    p_copy[.!ismissing.(p_copy)] .= adjust(p_copy[.!ismissing.(p_copy)] |> Vector{Float64}, d[lowercase(method)])
    return p_copy
end

"""
    nonmissing(vec::T) where T<:AbstractArray
Convenience function to count the number of non-`missing` values
in a vector.
"""
@inline function nonmissing(vec::T) where T<:AbstractArray
    count(!ismissing, vec)
end

#TODO add to API docs
"""
    nonmissing(data::PopData, locus::String)
Convenience function to count the number of non-`missing` samples
at a locus.
"""
@inline function nonmissing(data::PopData, locus::String)
    data.loci[data.loci.locus .== locus, :genotype] |> nonmissing
end

"""
    nonmissings(vec1::AbstractVector, vec2::AbstractVector)
Return a vector of indices where neither input vectors have a `missing` value, i.e. an
intersection of the indices of their non-missing elements.
"""
@inline function nonmissings(vec1::T, vec2::T) where T <: AbstractVector
    intersect(map(i -> findall(!ismissing, i), (vec1, vec2))...)
end

#TODO add to docs API
"""
    pairwise_pairs(smp_names::Vector{String})
Given a vector of sample names, returns a vector of tuples of unique all x
all combinations of sample pairs, excluding self-comparisons.

**Example**
```
julia> samps = ["red_1", "red_2", "blue_1", "blue_2"] ;

julia> pairwise_pairs(samps)
6-element Array{Tuple{String,String},1}:
 ("red_1", "red_2")
 ("red_1", "blue_1")
 ("red_1", "blue_2")
 ("red_2", "blue_1")
 ("red_2", "blue_2")
 ("blue_1", "blue_2")
```
"""
@inline function pairwise_pairs(smp_names::AbstractVector{String})
    [tuple(smp_names[i], smp_names[j]) for i in 1:length(smp_names)-1 for j in i+1:length(smp_names)]
end

#TODO add to docs API
"""
    phase(data::PopData)
Return a `Vector` of length `ploidy` composed of allele matrices with dimensions `samples × loci`.
Rows are samples and columns are loci. Will return an error if ploidy varies between samples.

**Example**
```
julia> mtx = phase(@nancycats)
2-element Array{Array{Union{Missing, Int16},2},1}:
 [missing 136 … 113 208; missing 146 … 113 208; … ; 137 130 … 113 208; 135 130 … missing 208]
 [missing 146 … 113 208; missing 146 … 113 208; … ; 143 136 … 117 208; 141 146 … missing 208]

julia> mtx[1]
237×9 Array{Union{Missing, Int16},2}:
    missing  136  139  116         156  142  199  113         208
    missing  146  139  120         156  142  185  113         208
 135         136  141  116         152  142  197  113         210
 133         138  139  116         150  142  199   91         208
 133         140  141  126         152  142  193  113         208
 135         136  145  120         150  148  193   91         208
 135         136  139  116         152  142  199  105         208
 135         136  135  120         154  142  193   91         208
 137         136  139  116         150  142  197  105         208
 135         132  141  120         150  148  197   91         208
 137         130  137  128         152  142  193   91         182
 129         130  135  126         144  140  193   91         182
   ⋮                                      ⋮
 133         136  135     missing  146  142  199  113         182
 133         136  135     missing  150  142  197  113         182
 133         130  141     missing  148  142  191     missing  208
 123         138  141     missing  148  142  191     missing  208
 123         138  139     missing  150  142  197     missing  208
 133         136  139     missing  150  142  197     missing  208
 133         130  139     missing  152  142  191     missing  208
 133         136  139     missing  150  142  199     missing  208
 133         130  135     missing  148  142  197     missing  208
 135         136  143     missing  144  142  191  113         208
 137         130  135     missing  150  142  193  113         208
 135         130  135     missing  150  142  197     missing  208
```
"""
function phase(data::PopData)
    dims = size(data)
    ploidy = unique(data.meta.ploidy)
    ploidy = length(ploidy) != 1 ? error("Phasing will not work on mixed-ploidy samples") : ploidy[1]
    sort_df = issorted(data.loci, [:name, :locus]) ? sort(data.loci, [:name, :locus]) : data.loci
    matrices = map(j -> map(i -> ismissing(i) ? missing : i[j] , sort_df.genotype), 1:ploidy)
    map(i -> collect(reshape(i, (dims.samples, dims.loci))), matrices)
end

"""
    quickstart()
Prints helpful text of how to get started using PopGen.
"""
function quickstart()
    printstyled("\n        Quickstart for PopGen\n\n", bold = true)
    println("Documentation: https://pdimens.github.io/PopGen.jl/")
    println("Motivational(?) quote: ", motivational_quote())
    println("\nA few things things you can do to get started:")

    printstyled("\nLoad in data\n\n", color = :magenta)
    println("- read_from(filename; kwargs...)")
    println("- genepop(infile; kwargs...)  or similar file-specific importer")
    println("- use available @gulfsharks or @nancycats datasets")

    printstyled("\nExplore PopData\n\n", color = :blue)
    println("- populations(PopData) to view population information")
    println("- loci(PopData) to view locus names")
    println("- samples(PopData) to view sample names")
    println("- missing(PopData, by = ...) to view missing information")

    printstyled("\nManipulate PopData\n\n", color = :light_red)
    println("- populations!(PopData, ...) to rename populations")
    println("- locations!(PopData, ...) to add geographical coordinates")
    println("- exclude!(PopData, kwargs...) to selectively remove data")

    printstyled("\nAnalyses\n\n", color = :green)
    println("- richness(PopData) to calculate allelic richness")
    println("- allele_avg(PopData) to calculate average # of alleles")
    println("- summary(PopData) to calculate F-statistics, heterozygosity, etc.")
    println("- hwe_test(PopData) to test for Hardy-Weinberg Equilibrium")

    return
end

"""
    reciprocal(num::T) where T <: Signed
Returns the reciprocal (1/number) of a number. Will return `0` when
the number is `0` instead of returning `Inf`.
"""
function reciprocal(num::T) where T <: Real
    !iszero(num) ? 1.0/float(num) : 0.0
end

"""
    sim_pairs(data::Vector{String})
Takes a Vector of sample names and returns a Tuple of sample pairs, grouped by simulation
number. This is an internal function used for isolating sibship pairs from simulated shipship
pairs (via `PopGenSims.jl`) to perform `relatedness` estimates only on those pairs.

**Example**
julia> a = ["sim1_off1", "sim1_off2", "sim2_off1", "sim2_off2"] ;

julia> sim_pairs(a)
("sim1_off1", "sim1_off2")
("sim2_off1", "sim2_off2")
"""
function sim_pairs(data::Vector{String})
    n = length(data)
    isodd(n) && error("Expected an even number of samples, but got $n")
    Tuple.(Base.Iterators.partition(sort(data), 2))
end

"""
    Base.sort(x::NTuple{N,T}) where N where T <: Signed
Sort the integers within a Tuple and return the sorted Tuple.
"""
function Base.sort(x::NTuple{N,T}) where N where T <: Signed
    Tuple(sort(SVector(x)))
end

#TODO add to API docs
"""
    strict_shuffle(x::T) where T <: AbstractArray
Shuffle only the non-missing values of a Vector and return a copy of the vector,
keeping the `missing` values at their original locations.
Use `strict_shuffle!` to edit in-place instead of returning a copy.
"""
@inline function strict_shuffle(x::T) where T <: AbstractArray
    # get indices of where original missing are
    miss_idx = findall(i -> i === missing, x)
    out_vec = shuffle(Xoroshiro128Star(), x[.!ismissing.(x)])

    insert!.(Ref(out_vec), miss_idx, missing)
    return out_vec
end

"""
    strict_shuffle!(x::T)! where T <: AbstractArray
Shuffle only the non-missing values of a Vector, keeping the
`missing` values at their original locations. Use `strict_shuffle`
to return a copy instead of editing in-place.
"""
function strict_shuffle!(x::T) where T <: AbstractArray
    @inbounds shuffle!(Xoroshiro128Star(), @view x[.!ismissing.(x)])
    return x
end

"""
    motivational_quote()
Randomly returns one of ~50 _somewhat_ motivational quotes.
"""
function motivational_quote()
# quotes retrieved on 06/11/2020 from https://www.coburgbanks.co.uk/blog/friday-funnies/50-funny-motivational-quotes/
    quotes =  [
    "\"The elevator to success is out of order. You’ll have to use the stairs, one step at a time.\" Joe Girard",
    "\"People often say that motivation doesn’t last. Well, neither does bathing – that’s why we recommend it daily.\" Zig Ziglar",
    "\"I always wanted to be somebody, but now I realise I should have been more specific.\" Lily Tomlin",
    "\"I am so clever that sometimes I don’t understand a single word of what I am saying.\" Oscar Wilde",
    "\"People say nothing is impossible, but I do nothing every day.\" Winnie the Pooh",
    "\"Life is like a sewer… what you get out of it depends on what you put into it.\" Tom Lehrer",
    "\"You can’t have everything. Where would you put it?\" Steven Wright",
    "\"Work until your bank account looks like a phone number.\" Unknown ",
    "\"Change is not a four letter word… but often your reaction to it is!\" Jeffrey Gitomer",
    "\"If you think you are too small to make a difference, try sleeping with a mosquito.\" Dalai Lama",
    "\"Bad decisions make good stories.\" Ellis Vidler",
    "\"I’ll probably never fully become what I wanted to be when I grew up, but that’s probably because I wanted to be a ninja princess.\" Cassandra Duffy",
    "\"When life gives you lemons, squirt someone in the eye.\" Cathy Guisewite",
    "\"A clear conscience is a sure sign of a bad memory.\" Mark Twain",
    "\"Well-behaved women seldom make history.\" Laurel Thatcher Ulrich",
    "\"I didn’t fail the test. I just found 100 ways to do it wrong.\" Benjamin Franklin",
    "\"I used to think I was indecisive, but now I’m not so sure.\" Unknown",
    "\"Don’t worry about the world coming to an end today. It’s already tomorrow in Australia.\" Charles Schulz",
    "\"Think like a proton. Always positive.\" Unknown",
    "\"Be happy – it drives people crazy.\" Unknown",
    "\"Optimist: someone who figures that taking a step backward after taking a step forward is not a disaster, it’s more like a cha-cha.\" Robert Brault",
    "\"The question isn’t who is going to let me, it’s who is going to stop me.\" Ayn Rand.",
    "\"You’re only given a little spark of madness. You mustn’t lose it.\" Robin Williams",
    "\"I am an early bird and a night owl… so I am wise and I have worms\" Michael Scott",
    "\"If you let your head get too big, it’ll break your neck.\" Elvis Presley",
    "\"The road to success is dotted with many tempting parking spaces.\" Will Rogers",
    "\"Live each day like it’s your second to the last. That way you can fall asleep at night.\" Jason Love",
    "\"Even a stopped clock is right twice every day. After some years, it can boast of a long series of successes.\" Marie Von Ebner-Eschenbach",
    "\"Honest criticism is hard to take, particularly from a relative, a friend, an acquaintance, or a stranger.\" Franklin P. Jones",
    "\"Opportunity is missed by most people because it is dressed in overalls and looks like work.\" Thomas Eddison",
    "\"A diamond is merely a lump of coal that did well under pressure.\" Unknown",
    "\"Nothing is impossible, the word itself says “I’m possible!\" Audrey Hepburn",
    "\"Friendship is like peeing on yourself: everyone can see it, but only you get the warm feeling that it brings.\" Robert Bloch",
    "\"Women who seek to be equal with men lack ambition.\" Marilyn Monroe",
    "\"By working faithfully eight hours a day you may eventually get to be boss and work twelve hours a day.\" Robert Frost",
    "\"The trouble with having an open mind, of course, is that people will insist on coming along and trying to put things in it.\" Terry Pratchett",
    "\"Age is of no importance unless you’re a cheese.\" Billie Burke",
    "\"When tempted to fight fire with fire, remember that the Fire Department usually uses water.\" Unknown",
    "\"Trying is the first step toward failure.\" Homer Simpson",
    "\"Happiness is just sadness that hasn’t happened yet.\" Unknown",
    "\"The best things in life are actually really expensive.\" Unknown",
    "\"A few harmless flakes working together can unleash an avalanche of destruction.\" Justin Sewell",
    "\"It could be that your purpose in life is to serve as a warning to others.\"  Ashleigh Brilliant",
    "\"If the world didn’t suck we’d all fly into space.\" Unknown",
    "\"Always remember that you are unique – just like everybody else.\" Unknown"
    ]
    return quotes[rand(1:length(quotes))]
end

#TODO add to docs API
"""
    generate_meta(data::DataFrame)
Given a genotype DataFrame formatted like `PopData.loci`, generates a corresponding
`meta` DataFrame. In other words, it creates the `.meta` part of `PopData` from the `.loci` part.

**Example**:
```
julia> cats = @nancycats ;

julia> cats_nometa = cats.loci ;

julia> cats_meta = generate_meta(cats_nometa)
237×5 DataFrame
 Row │ name    population  ploidy  longitude  latitude
     │ String  String      Int8    Float32?   Float32?
─────┼─────────────────────────────────────────────────
   1 │ N215    1                2   missing   missing
   2 │ N216    1                2   missing   missing
   3 │ N217    1                2   missing   missing
   4 │ N218    1                2   missing   missing
   5 │ N219    1                2   missing   missing
   6 │ N220    1                2   missing   missing
   7 │ N221    1                2   missing   missing
  ⋮  │   ⋮         ⋮         ⋮         ⋮         ⋮
 232 │ N295    17               2   missing   missing
 233 │ N296    17               2   missing   missing
 234 │ N297    17               2   missing   missing
 235 │ N281    17               2   missing   missing
 236 │ N289    17               2   missing   missing
 237 │ N290    17               2   missing   missing
                                       224 rows omitted
```
"""
function generate_meta(data::DataFrame)
    grp = groupby(data, :name)
    nms = map(z -> z.name, keys(grp))
    pops = map(z -> first(z.population), grp)
    ploids = map(z -> find_ploidy(z.genotype), grp)
    DataFrame(
        :name => nms,
        :population => pops,
        :ploidy => ploids,
        :longitude => Vector{Union{Missing, Float32}}(undef, (length(nms))),
        :latitude => Vector{Union{Missing, Float32}}(undef, (length(nms)))
    )
end


text:
module TestPolyAAnalsis

using BioSequences
using Base.Test
NEW_PATH=join(split(Base.source_path(),"/")[1:end-3],"/")
push!(LOAD_PATH, NEW_PATH)
using DataFrames
using DataStructures
import GenomicFeatures: GFF3
import GenomicFeatures: Interval
import GenomicFeatures: Strand
import GenomicFeatures: IntervalCollection
using PolyAAnalysis

@testset "PolyAtrimming" begin
include("trim_polyA_3end.jl")
include("trim_polyA_from_fastq_record.jl")
include("detect_polyA_in_a_string.jl")
include("get_polyA_prefixes.jl")
end

@testset "AnnotatePolyA" begin
include("AnnotatePolyA.jl")
end

@testset "ParseGFF" begin
include("ParseGFF.jl")
end
end #module TestPolyAanalsis

max_stars_repo_name: laurabiggins/polyAterminus

max_issues_repo_name: laurabiggins/polyAterminus

max_forks_repo_name: laurabiggins/polyAterminus


text:
# AES256CBC.jl

VERSION >= v"0.4.0-dev+6521" && __precompile__()
module AES256CBC

import OpenSSL
export UBytes
export string2bytes, genRandUBytes, genKey32Iv16, genPadding
export encryptAES256CBC, decryptAES256CBC

#  (key32, iv16) = genKey32Iv16(passwd, salt)
#    input:
#      passwd (any bytes) (from command line)
#      salt (8 bytes) (from Salted__????????)
#    process:
#      s1 = md5(passwd + salt)      (16 bytes)
#      s2 = md5(s1 + passwd + salt) (16 bytes)
#      s3 = md5(s2 + passwd + salt) (16 bytes)
#    output:
#      key32 = s1 + s2 (32 bytes)
#      iv16 = s3 (16 bytes)

typealias UBytes Array{UInt8, 1}

function string2bytes(s::AbstractString)
  return read(IOBuffer(s), UInt8, length(s)) # UBytes
end

function genRandUBytes(n::Int)
  # srand()
  return rand(UInt8(0):UInt8(255), n)
end

function genKey32Iv16(passwd::UBytes, salt::UBytes)
  OpenSSL.init()
  s1 = OpenSSL.Digest.digest("MD5", [passwd; salt])
  s2 = OpenSSL.Digest.digest("MD5", [s1; passwd; salt])
  s3 = OpenSSL.Digest.digest("MD5", [s2; passwd; salt])
  OpenSSL.cleanup()
  return ([s1; s2], s3) # (key32::UBytes, iv16::UBytes)
end

function genPadding(n::Int)
  return UBytes(map(b -> n, 1:n))
end

function encryptAES256CBC(key32::UBytes, iv16::UBytes, plain::UBytes)
  OpenSSL.init()
  e = OpenSSL.Cipher.encrypt("aes_256_cbc", key32, iv16, plain) # auto padding
  OpenSSL.cleanup()
  return e
end

function decryptAES256CBC(key32::UBytes, iv16::UBytes, cipher::UBytes)
  OpenSSL.init()
  d = OpenSSL.Cipher.decrypt("aes_256_cbc", key32, iv16, cipher) # auto padding
  OpenSSL.cleanup()
  return d
end

end


text:
struct PrisonersDilemma end

n_agents(simpleGame::PrisonersDilemma) = 2

ordered_actions(simpleGame::PrisonersDilemma, i::Int) = [:cooperate, :defect]
ordered_joint_actions(simpleGame::PrisonersDilemma) = vec(collect(Iterators.product([ordered_actions(simpleGame, i) for i in 1:n_agents(simpleGame)]...)))

n_joint_actions(simpleGame::PrisonersDilemma) = length(ordered_joint_actions(simpleGame))
n_actions(simpleGame::PrisonersDilemma, i::Int) = length(ordered_actions(simpleGame, i))

function reward(simpleGame::PrisonersDilemma, i::Int, a)
    if i == 1
        noti = 2
    else
        noti = 1
    end

    if a[i] == :cooperate && a[noti] == :cooperate
        return -1.0
    elseif a[i] == :defect && a[noti] == :cooperate
        return 0.0
    elseif a[i] == :cooperate && a[noti] == :defect
        return -4.0
    elseif a[i] == :defect && a[noti] == :defect
        return -3.0
    end
end

function joint_reward(simpleGame::PrisonersDilemma, a)
    return [reward(simpleGame, i, a) for i in 1:n_agents(simpleGame)]
end

function SimpleGame(simpleGame::PrisonersDilemma)
    return SimpleGame(
        0.9,
        vec(collect(1:n_agents(simpleGame))),
        [ordered_actions(simpleGame, i) for i in 1:n_agents(simpleGame)],
        (a) -> joint_reward(simpleGame, a)
    )
end


text:
abstract type DistributionFunc
end

struct BoltzmanDist <: DistributionFunc
end

struct FermiDist <: DistributionFunc
end

struct BlakemoreDist <: DistributionFunc
    b::Number
end
BlakemoreDist() = BlakemoreDist(0.27)

function (b::BoltzmanDist)(E,Ef,T)
    return b(E-Ef,T)
end

function (b::BoltzmanDist)( ΔE,T)
    return b((ΔE)/(kb*T))
end

function (b::BoltzmanDist)( ΔE)
    return exp(ΔE)
end

function (b::FermiDist)(E,Ef,T)
    return b(E-Ef,T)
end

function (b::BlakemoreDist)(E,Ef,T)
    return b(E-Ef,T)
end

"""
# Arguments
    eta
Calculation One Half Order Fermi Integral of eta, using trapezoidal integral to approximat it
the accuracy was compared to Rational Chebychev approximations, and Gaussian Quadrature Calculations
with maximum absolute error smaller than 1e-12.
The defition of the Fermi Integral is such it approaches exp(-η) without normalization constants
"""
function trapFermi(eta)
    aj = 0.5
    range = 8.0
    if eta > 0.0
        range = sqrt(eta + 64.)
    end
    h = 0.5
    nmax = range / h
    sum = 0.
    for i = 1:nmax
        u = i * h
        ff = 2 * (u^(2 * aj + 1)) / (1 + exp(u * u - eta))
        sum = sum + ff
    end
    pol = 0.
    npole = 0
    bk1 = 0
    while bk1 <= 14 * pi
        npole = npole + 1
        bk = (2 * npole - 1) * pi
        rho = sqrt(eta * eta + bk * bk)
        t1 = 1
        t2 = 0
        if eta < 0
            tk = - aj * (atan(-bk / eta) + pi)
        elseif eta == 0
        tk = 0.5 * pi * aj
        else
            tk = aj * atan(bk / eta)
        end
        rk = - (rho^aj)
        tk = tk + 0.5 * atan(t2 / t1)
        if eta < 0
            rk = -rk
        end
        ak = (2. * pi / h) * sqrt(0.5 * (rho + eta))
        bk1 = (2. * pi / h) * sqrt(0.5 * (rho - eta))
        if bk1 <= (14. * pi)
            gama = exp(bk1)
            t1 = gama * sin(ak + tk) - sin(tk)
            t2 = 1.0 - 2.0 * gama * cos(ak) + gama * gama
            pol = pol + 4.0 * pi * rk * t1 / t2
        end
    end
    npole = npole - 1
    fdp = (sum * h + pol)*2/sqrt(pi)
end

"""
# Arguments
    eta
Calculation One Minus Half Order Fermi Integral of eta, using trapezoidal integral to approximat it
the accuracy was compared to Rational Chebychev approximations, and Gaussian Quadrature Calculations
with maximum absolute error smaller than 1e-12.
The defition of the Fermi Integral is such it approaches exp(-η) without normalization constants
"""
function DtrapFermi(eta)
    aj = -0.5
    range = 8.0
    if eta > 0.0
        range = sqrt(eta + 64.)
    end
    h = 0.5
    nmax = range / h
    sum = 0.
    for i = 1:nmax
        u = i * h
        ff = 2 * (u^(2 * aj + 1)) / (1 + exp(u * u - eta))
        sum = sum + ff
    end
    pol = 0.
    npole = 0
    bk1 = 0
    while bk1 <= 14 * pi
        npole = npole + 1
        bk = (2 * npole - 1) * pi
        rho = sqrt(eta * eta + bk * bk)
        t1 = 1
        t2 = 0
        if eta < 0
            tk = - aj * (atan(-bk / eta) + pi)
        elseif eta == 0
        tk = 0.5 * pi * aj
        else
            tk = aj * atan(bk / eta)
        end
        rk = - (rho^aj)
        tk = tk + 0.5 * atan(t2 / t1)
        if eta < 0
            rk = -rk
        end
        ak = (2. * pi / h) * sqrt(0.5 * (rho + eta))
        bk1 = (2. * pi / h) * sqrt(0.5 * (rho - eta))
        if bk1 <= (14. * pi)
            gama = exp(bk1)
            t1 = gama * sin(ak + tk) - sin(tk)
            t2 = 1.0 - 2.0 * gama * cos(ak) + gama * gama
            pol = pol + 4.0 * pi * rk * t1 / t2
        end
    end
    npole = npole - 1
    fdp = (sum * h + pol)*2/sqrt(pi)
end

"""
# Arguments
    eta
Calculation Three Half Order Fermi Integral of eta, using trapezoidal integral to approximat it
the accuracy was compared to Rational Chebychev approximations, and Gaussian Quadrature Calculations
with maximum absolute error smaller than 1e-12.
The defition of the Fermi Integral is such it approaches exp(-η) without normalization constants
"""
function trapFermi32(eta)
    aj = 1.5
    range = 8.0
    if eta > 0.0
        range = sqrt(eta + 64.)
    end
    h = 0.5
    nmax = range / h
    sum = 0.
    for i = 1:nmax
        u = i * h
        ff = 2 * (u^(2 * aj + 1)) / (1 + exp(u * u - eta))
        sum = sum + ff
    end
    pol = 0.
    npole = 0
    bk1 = 0
    while bk1 <= 14 * pi
        npole = npole + 1
        bk = (2 * npole - 1) * pi
        rho = sqrt(eta * eta + bk * bk)
        t1 = 1
        t2 = 0
        if eta < 0
            tk = - aj * (atan(-bk / eta) + pi)
        elseif eta == 0
            tk = 0.5 * pi * aj
        else
            tk = aj * atan(bk / eta)
        end
        rk = - (rho^aj)
        tk = tk + 0.5 * atan(t2 / t1)
        if eta < 0
            rk = -rk
        end
        ak = (2. * pi / h) * sqrt(0.5 * (rho + eta))
        bk1 = (2. * pi / h) * sqrt(0.5 * (rho - eta))
        if bk1 <= (14. * pi)
            gama = exp(bk1)
            t1 = gama * sin(ak + tk) - sin(tk)
            t2 = 1.0 - 2.0 * gama * cos(ak) + gama * gama
            pol = pol + 4.0 * pi * rk * t1 / t2
        end
    end
    npole = npole - 1
    fdp = (sum * h + pol)*2/sqrt(pi)
end

function Fermi12(x)
    Integrand(t,x) =  sqrt(t)/(exp(t-x)+1)
    t = Fun(t-> Integrand(t,x),0..10)
    return sum(t)
end

function (b::FermiDist)(ΔE,T)
    return b(ΔE/kb/T)
end
function (b::FermiDist)(ΔE)
    return trapFermi(ΔE)
end

function (b::BlakemoreDist)(ΔE,T)
    return b(ΔE/kb/T)
end
function (b::BlakemoreDist)(ΔE)
    return 1/(exp(-ΔE)+b.b)
end

function Primitive(b::BlakemoreDist,x)
    return log(b.b*exp(x)+1)/(b.b)
end

function Primitive(b::BoltzmanDist,x)
    return b(x)
end

function Primitive(b::FermiDist,x)
    return trapFermi32(x)
end


text:
"""
    extract_fault_property(fault::LoE)

Extract LoE's fault property
"""
function extract_fault_property(fault::LoE)
    @unpack time, index, level = fault
    (; time=time, index=index, level=level,)
end

"""
faults can be a vector of `AbstractFault`s, e.g. FaultSet
"""
function extract_fault_property(faults::Vector{AbstractFault})
    faults |> Map(extract_fault_property) |> collect |> unique |> sort
end

"""
    get_fault_kind(faults::Vector{AbstractFault})

Get fault kind from `faults::Vector{AbstractFault}`.

# Notes
- failure: fully faulted
- fault: partially faulted
"""
function get_fault_kind(faults::Vector{AbstractFault})
    fault_kind = nothing
    # single fault or failure
    if length(faults) == 1
        fault = faults[1]
        @unpack index, level = fault
        # single failure
        if level == 0.0
            fault_kind = "single_failure_$(string(faults[1].index))"
        # single fault
        else
            fault_kind = "single_fault_$(string(faults[1].index))"
        end
    # double fault or failure
    elseif length(faults) == 2
        fault1, fault2 = faults
        # double fault or failure
        if fault1.level == fault2.level
            # double failure
            if fault1.level == 0.0
                fault_kind = "double_failure"
            # double fault
            else
                fault_kind = "double_fault"
            end
        # single fauliure and single fault
        else
            fault_kind = "single_failure_single_fault"
        end
    else
        error("Undefined fault kind")
    end
    fault_kind
end


text:
mutable struct MultiplicativeScalePrior <: MyXMLElement
    mults::Parameter
    precs::MultiplicativeParameter
    mult_prior::MultivariateGammaLikelihood
    scale_prior::MyXMLElement
end

function make_xml(mps::MultiplicativeScalePrior)
    make_xml(mps.mults)
    make_xml(mps.precs)
    make_xml(mps.mult_prior)
    make_xml(mps.scale_prior)
    return nothing
end

function get_loggables(mps::MultiplicativeScalePrior)
    make_xml(mps)
    return [mps.mults.el]
end

function get_loadings_prior(msp::MultiplicativeScalePrior)
    return msp.scale_prior
end

function NormalGammaPrecisionOperatorXMLElement(
    multiplicativeLikelihood::MultiplicativeScalePrior)

    ggp = MultipilcativeGammaGibbsProvider(multiplicativeLikelihood.mults,
                multiplicativeLikelihood.scale_prior)

    return NormalGammaPrecisionOperatorXMLElement(nothing,
                ggp,
                multiplicativeLikelihood.mult_prior,
                Int[],
                1.0)
end

function set_shrinkage_mults!(msp::MultiplicativeScalePrior;
                        shapes::Vector{Float64} = Float64[],
                        scales::Vector{Float64} = Float64[],
                        set_mults::Bool = true)
    μ = shapes .* scales
    set_shrinkage_mults!(msp.mult_prior, shapes=shapes, scales=scales)
    if set_mults
        set_value(msp.mults, μ)
    end
    scale_mode = [prod(μ[1:i]) for i = 1:length(μ)]
    scale_mode = scale_mode.^(-0.5)
    return scale_mode
end

function get_priors(msp::MultiplicativeScalePrior)
    make_xml(msp)
    return [msp.mult_prior.el, msp.scale_prior.el]
end

function get_normal_prior(msp::MultiplicativeScalePrior)
    return get_normal_prior(msp.scale_prior)
end


text:
abstract type Ensemble end

struct CanonicalEnsemble <: Ensemble end
struct MicrocanonicalEnsemble <: Ensemble end
struct PotentialEnergyLandscapeEnsemble <: Ensemble
    E_L::Float64
end


text:
using Printf
using KiteModels, KitePodModels, KiteUtils

# change this to KPS3 or KPS4
const Model = KPS4

if ! @isdefined kcu;  const kcu = KCU(se());   end
if ! @isdefined kps4; const kps4 = Model(kcu); end

# the following values can be changed to match your interest
dt = 0.05
STEPS = 600
PLOT = true
FRONT_VIEW = false
ZOOM = false
PRINT = false
STATISTIC = false
# end of user parameter section #

if PLOT
    using Plots
    include("plot2d.jl")
end

v_time = zeros(STEPS)
v_speed = zeros(STEPS)
v_force = zeros(STEPS)

function simulate(integrator, steps, plot=false)
    start = integrator.p.iter
    for i in 1:steps
        if PRINT
            lift, drag = KiteModels.lift_drag(kps4)
            @printf "%.2f: " round(integrator.t, digits=2)
            println("lift, drag  [N]: $(round(lift, digits=2)), $(round(drag, digits=2))")
        end
        acc = 0.0
        if kps4.t_0 > 15.0
            acc = 0.1
        end
        v_ro = kps4.sync_speed+acc*dt
        v_time[i] = kps4.t_0
        v_speed[i] = kps4.v_reel_out
        v_force[i] = winch_force(kps4)
        KiteModels.next_step!(kps4, integrator, v_ro = v_ro, dt=dt)

        if plot
            reltime = i*dt
            if mod(i, 5) == 0
                p = plot2d(kps4.pos, reltime; zoom=ZOOM, front=FRONT_VIEW, segments=se().segments)
                display(p)
            end
        end
    end
    (integrator.p.iter - start) / steps
end

integrator = KiteModels.init_sim!(kps4, stiffness_factor=0.04, prn=STATISTIC)
kps4.sync_speed = 0.0

if PLOT
    av_steps = simulate(integrator, STEPS, true)
else
    println("\nStarting simulation...")
    simulate(integrator, 100)
    runtime = @elapsed av_steps = simulate(integrator, STEPS-100)
    println("\nTotal simulation time: $(round(runtime, digits=3)) s")
    speed = (STEPS-100) / runtime * dt
    println("Simulation speed: $(round(speed, digits=2)) times realtime.")
end
lift, drag = KiteModels.lift_drag(kps4)
println("lift, drag  [N]: $(round(lift, digits=2)), $(round(drag, digits=2))")
println("Average number of callbacks per time step: $av_steps")

p1 = plot(v_time, v_speed, ylabel="v_reelout  [m/s]", legend=false)
p2 = plot(v_time, v_force, ylabel="tether_force [N]", legend=false)
plot(p1, p2, layout = (2, 1), legend = false)


text:
include("RCF/Types.jl")
include("RCF/cyclo.jl")
include("RCF/rcf.jl")
include("RCF/conductor.jl")
include("RCF/conductor_sieve.jl")


text:
@testset "Felzenszwalb" begin
    if Base.VERSION < v"1.4.0"
        _only(c) = (@assert length(c) == 1; c[1])
    else
        _only = only
    end

    img = zeros(Gray{N0f8}, 10, 10)
    img[2:3, 2:3] .= 0.2
    img[2, 8] = 0.5
    img[5:6, 5:6] .= 0.8
    T = Gray{N0f8}

    result = felzenszwalb(img, 1)

    @test length(result.segment_labels) == 4
    @test result.segment_labels == collect(1:4)

    @test all(label->(label==result.image_indexmap[2,2]), result.image_indexmap[2:3, 2:3])
    @test all(label->(label==result.image_indexmap[5,5]), result.image_indexmap[5:6, 5:6])

    TM = ImageSegmentation.meantype(T)
    @test result.segment_means[result.image_indexmap[1,1]] == zero(TM)
    @test result.segment_means[result.image_indexmap[2,2]] == TM(img[2,2])
    @test result.segment_means[result.image_indexmap[2,8]] == TM(img[2,8])
    @test result.segment_means[result.image_indexmap[5,5]] == TM(img[5,5])

    @test result.segment_pixel_count[result.image_indexmap[1,1]] == 91
    @test result.segment_pixel_count[result.image_indexmap[2,2]] == 4
    @test result.segment_pixel_count[result.image_indexmap[2,8]] == 1
    @test result.segment_pixel_count[result.image_indexmap[5,5]] == 4

    img = zeros(RGB{Float64}, 10, 10)
    img[2:3, 2:3] .= RGB(0.2,0.2,0.2)
    img[2, 8] = RGB(0.5,0.5,0.5)
    img[5:6, 5:6] .= RGB(0.8,0.8,0.8)
    T = RGB{Float64}
    TM = ImageSegmentation.meantype(T)

    result = felzenszwalb(img, 1, 2)

    @test length(result.segment_labels) == 3
    @test result.segment_labels == collect(1:3)

    @test result.image_indexmap[1, 1] == result.image_indexmap[2, 8]
    @test all(label->(label==result.image_indexmap[2,2]), result.image_indexmap[2:3, 2:3])
    @test all(label->(label==result.image_indexmap[5,5]), result.image_indexmap[5:6, 5:6])

    @test result.segment_means[result.image_indexmap[1,1]] ≈ RGB{Float64}(0.5/92, 0.5/92, 0.5/92)
    @test result.segment_means[result.image_indexmap[2,2]] == TM(img[2,2])
    @test result.segment_means[result.image_indexmap[5,5]] == TM(img[5,5])

    @test result.segment_pixel_count[result.image_indexmap[1,1]] == 92
    @test result.segment_pixel_count[result.image_indexmap[2,2]] == 4
    @test result.segment_pixel_count[result.image_indexmap[5,5]] == 4

    # issue 20
    img = falses(10, 10)
    img[2:3, 2:3] .= true
    img[2, 8] = true
    img[5:6, 5:6] .= true
    T = Bool
    TM = ImageSegmentation.meantype(T)

    result = felzenszwalb(img, 1, 2)

    @test length(result.segment_labels) == 3
    @test result.segment_labels == collect(1:3)

    @test result.image_indexmap[1, 1] == result.image_indexmap[2, 8]
    @test all(label->(label==result.image_indexmap[2,2]), result.image_indexmap[2:3, 2:3])
    @test all(label->(label==result.image_indexmap[5,5]), result.image_indexmap[5:6, 5:6])

    @test result.segment_means[result.image_indexmap[1,1]] ≈ 1/92
    @test result.segment_means[result.image_indexmap[2,2]] == TM(img[2,2])
    @test result.segment_means[result.image_indexmap[5,5]] == TM(img[5,5])

    @test result.segment_pixel_count[result.image_indexmap[1,1]] == 92
    @test result.segment_pixel_count[result.image_indexmap[2,2]] == 4
    @test result.segment_pixel_count[result.image_indexmap[5,5]] == 4

    img = falses(10, 10, 10)
    img[2:3, 2:3, 2:3] .= true
    img[5:6, 5:6, 5:6] .= true
    result = felzenszwalb(img, 1, 2)
    @test result.segment_labels == collect(1:3)
    ibkg = result.image_indexmap[1, 1, 1]
    i = _only(unique(result.image_indexmap[2:3, 2:3, 2:3]))
    @test i != ibkg
    result.image_indexmap[2:3, 2:3, 2:3] .= ibkg
    j = _only(unique(result.image_indexmap[5:6, 5:6, 5:6]))
    @test j != ibkg && j != i
    result.image_indexmap[5:6, 5:6, 5:6] .= ibkg
    @test all(==(ibkg), result.image_indexmap)
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Helpers.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export @help, @stdout_to_pager, @out2pr

"""
    @help(f)

Open the documentation of the function `f` in pager.

# Examples

```julia-repl
julia> @help write
```
"""
macro help(f)
    ex_out = quote
        @doc($(esc(f))) |> pager
    end

    return ex_out
end

"""
    @stdout_to_pager(ex_in)

Capture the `stdout` generated by `ex_in` and show inside a pager.

!!! note
    The command **must** write to `stdout` explicitly. For example,
    `@stdout_to_pager 1` shows a blank screen since `1` does not write to
    `stdout`, but returns `1`. `@stdout_to_pager show(1)`, on the other hand,
    shows the number `1` inside the pager.

!!! note
    This macro can also be called using the shorter name `@out2pr`.
"""
macro stdout_to_pager(ex_in)
    ex_out = quote
        hascolor = get(stdout, :color, true)
        old_stdout = stdout
        buf = IOBuffer()
        io = IOContext(buf, :color => hascolor, :limit => false)
        try
            Base.eval(:(stdout = $io))
            $(esc(ex_in))
            Base.eval(:(stdout = $old_stdout))
            String(take!(buf)) |> pager
            close(io)
        finally
            Base.eval(:(stdout = $old_stdout))
        end
    end

    return ex_out
end

macro out2pr(ex)
    return :(@stdout_to_pager $(esc(ex)))
end


text:
# ----------------------------------------------------
# --- Finding USRP
# ----------------------------------------------------
mutable struct uhd_string_vector
end

"""
Find all connected USRP devices
out = uhd_find_devices(args)
- args : String with uhd_find_devices argument (such as "addr=192.168.10.16"
"""
function uhd_find_devices(stringArgs)
	args	=  Base.unsafe_convert(Cstring,stringArgs);
    # --- Create an handler to get the string
    global stringHandler = uhd_string_vector_make()

    # --- Call UHD Find device
    ccall((:uhd_usrp_find,libUHD),uhd_error,(Cstring,Ptr{Ptr{uhd_string_vector}}),args,stringHandler)

    vectorSize = uhd_string_vector_size(stringHandler)
    allStr = String[]
    if vectorSize == 0
        @info "No UHD devices found. Try with \"addr=xxx.xxx.x.x\" to specify the USRP IP address"
    else
        @info "Found $(vectorSize) USRP devices"
        for n = 0 : 1 : vectorSize - 1
            # --- Recreate the string for a fancy output
            stringOut = uhd_string_vector_at(stringHandler,0)
            global allS = String.(split(stringOut,","))
            ee = "UHD Device $n \n"*join([" "*allS[i]*"\n" for i ∈ 1:length(allS)])
            customPrint(ee,"UHD",bold=true,color=:magenta)
            # --- FIXME Populate a dict to get USRP IP addr
            push!(allStr,ee)
        end
    end
    # --- Free the structure
    uhd_string_vector_free(stringHandler)
    # Return number of UHD device found
    return allStr
end

# In case of no input parameters, equivalent to call with ""
uhd_find_devices() = uhd_find_devices("")

"""
Create a string handler usefull to get/take messages from the USRP
"""
function uhd_string_vector_make()
    stringHandler = Ref{Ptr{uhd_string_vector}}()
    ccall((:uhd_string_vector_make,libUHD),uhd_error,(Ptr{Ptr{uhd_string_vector}},),stringHandler)
    return stringHandler
end

"""
Free a string handler craeted by uhd_string_vector_make()
"""
function uhd_string_vector_free(stringHandler)
    # --- Free the structure
    ccall((:uhd_string_vector_free,libUHD),uhd_error,(Ptr{Ptr{uhd_string_vector}},),stringHandler);
end

"""
Get the UHD internal string size. Takes a stringHandler as input parameter and returns a Int (Number of strings)
The string handler can be obtained with uhd_string_vector_make
"""
function uhd_string_vector_size(stringHandler)
    # --- Pointer for answer
    vectorSize = Ref{Csize_t}(0)
    # --- CCall
    errCall = ccall((:uhd_string_vector_size,libUHD),uhd_error,(Ptr{uhd_string_vector},Ref{Csize_t}),stringHandler[],vectorSize)
    return vectorSize[]
end

"""
Get the string at position index in the stringHandler container
The string handler can be obtained with uhd_string_vector_make
"""
function uhd_string_vector_at(stringHandler,index)
    # --- Init container
    v = Vector{UInt8}(undef,1024)
    # --- Get the data
    errCall = ccall((:uhd_string_vector_at,libUHD),uhd_error,(Ptr{uhd_string_vector},Csize_t,Ptr{Cchar},Csize_t),stringHandler[],index,v,1024)
    # --- Convert it into String
    strOut = String(v)
    # --- Find end of string
    indexM = findfirst('\0',strOut)
    if indexM > 1
        # Non null
        stringOut = strOut[1:indexM-1]
    else
        stringOut = ""
    end
    return stringOut
end


text:
mutable struct MarkovChain{S <: DiseaseStateSequence, M <: ILM}
  iterations::Int64
  events::Vector{Events{S}}
  transmission_network::Vector{TransmissionNetwork}
  risk_parameters::Vector{RiskParameters{S}}
  substitution_model::Union{Nothing, Vector{NucleicAcidSubstitutionModel}}
  log_posterior::Vector{Float64}
  Σrp::OnlineStats.CovMatrix
  Σsm::Union{Nothing, OnlineStats.CovMatrix}

  function MarkovChain{S, M}(e::Events{S},
                             n::TransmissionNetwork,
                             rp::RiskParameters{S},
                             lp::Float64) where {
                             S <: DiseaseStateSequence,
                             M <: TNILM}
    return new{S, M}(0, [e], [n], [rp], nothing, [lp], OnlineStats.CovMatrix(), nothing)
  end

  function MarkovChain{S, M}(e::Events{S},
                             n::TransmissionNetwork,
                             rp::RiskParameters{S},
                             sm::NucleicAcidSubstitutionModel,
                             lp::Float64) where {
                             S <: DiseaseStateSequence,
                             M <: PhyloILM}
    return new{S, M}(0, [e], [n], [rp], [sm], [lp], OnlineStats.CovMatrix(), OnlineStats.CovMatrix())
  end
end

function Base.length(x::MarkovChain)
  return x.iterations
end

function Base.show(io::IO, x::MarkovChain{S, M}) where {
                   S <: DiseaseStateSequence,
                   M <: ILM}
  return print(io, "$S $M Markov chain (iterations = $(length(x)))")
end

function TransmissionNetworkDistribution(x::MarkovChain; burnin::Int64=0, thin::Int64=1)
  return TNDistribution(x.transmission_network[1+burnin:thin:end])
end


text:
@testset "spectral operators" begin
    function P5(r::AbstractVector{T}) where {T}
        return T(1) / T(8) * (T(15) * r - T(70) * r .^ 3 + T(63) * r .^ 5)
    end
    function P6(r::AbstractVector{T}) where {T}
        return T(1) / T(16) * (-T(5) .+ T(105) * r .^ 2 - T(315) * r .^ 4 + T(231) * r .^ 6)
    end
    function DP6(r::AbstractVector{T}) where {T}
        return T(1) / T(16) * (T(2 * 105) * r - T(4 * 315) * r .^ 3 + T(6 * 231) * r .^ 5)
    end

    IPN(::Type{T}, N) where {T} = T(2) / T(2 * N + 1)

    for T in (Float32, Float64, BigFloat)
        r, w = legendregausslobatto(T, 7)
        x = LinRange{T}(-1, 1, 101)
        D = spectralderivative(r)
        P = spectralinterpolation(r, x)

        @test sum(P5(r) .^ 2 .* w) ≈ IPN(T, 5)
        @test D * P6(r) ≈ DP6(r)
        @test P * P6(r) ≈ P6(x)
    end

    for T in (Float32, Float64, BigFloat)
        r, w = legendregauss(T, 7)
        D = spectralderivative(r)

        @test sum(P5(r) .^ 2 .* w) ≈ IPN(T, 5)
        @test sum(P6(r) .^ 2 .* w) ≈ IPN(T, 6)
        @test D * P6(r) ≈ DP6(r)
    end
end


text:
"""
    BIC(sem_fit::SemFit)

Return the bayesian information criterion.
"""
BIC(sem_fit::SemFit) = minus2ll(sem_fit) + log(n_obs(sem_fit))*n_par(sem_fit)


text:
using EngEconomics, Roots, Plots

# Given
aPW = 56740
aLife = 5
aInit = -180000
MARR = 0.10
# Equal cash flows at the end of each year

# Find
aAnnual = (aPW - aInit) / seriesPresentAmountFactor(MARR, aLife)

IRR_eq(x) = aInit + aAnnual * seriesPresentAmountFactor(x, aLife)

xVec = collect(0:0.01:1.0)
yVec = zeros(size(xVec))
plot(xVec, IRR_eq.(xVec))
plot!(xVec, yVec)
irr = find_zero(IRR_eq, 0.12)


text:
##############################################################
#                   ALGORITMO MEMÉTICO
##############################################################

using Random

include("../algoritmos-geneticos/generar-cromosoma.jl")
include("../algoritmos-geneticos/torneo-binario.jl")
include("../utils/distancias.jl")
include("../utils/funcion-objetivo.jl")
include("busqueda-local.jl")
"""
    AuxCrucePorIndices(i,j,g, FunciónCruce)
Función auxiliar para ahorrarse el cruce si i==j
"""
function AuxCrucePorIndices(i,j,g, FunciónCruce)
    if i == j
        return g[i], g[i]
    else
        return FunciónCruce(g[i], g[j])
    end
end
"""
    Memetico(     evaluaciones_máximas_FE,
                    numero_cromosomas_por_generación,
                    probabilidad_cruce,
                    probabilidad_mutación,
                    tamaño_cromosoma,
                    función_evaluación,
                    función_cruce,
                    función_mutación,
                    porcentaje_cromosomas_busqueda_local,
                    funcion_selección_índices_para_busqueda_local,
                    numero_generaciones_aplicar_bl = 10,
    )
    Algoritmo genético generacional, tras las iteraciones convenientes devuelve el cromosoma mejor y su función de valuación
"""
function Memetico(   evaluaciones_máximas_FE,
                numero_cromosomas_por_generación,
                probabilidad_cruce,
                probabilidad_mutación,
                tamaño_cromosoma,
                función_evaluación,
                función_cruce,
                función_mutación,
                porcentaje_cromosomas_busqueda_local,
                funcion_selección_índices_para_busqueda_local,
                numero_generaciones_aplicar_bl = 10,
                )
    # Generamos la primera generación
    generación = [GeneraCromosoma(tamaño_cromosoma, 0.1) for i in 1:numero_cromosomas_por_generación]
    evaluaciones_función_evaluación = map(x -> función_evaluación(x), generación)

    # --- Cálculo de valores auxiliares ---
    # Mitad población
    M = floor(Int, numero_cromosomas_por_generación/2)
    ### Cruce
    # Número de parejas a cruzar
    numero_cruces = round(Int, probabilidad_cruce* numero_cromosomas_por_generación/2)
    permutaciones = randperm(M)
    parejas_a_cruzar = permutaciones[1:numero_cruces]
    parejas_que_no_se_cruzan = permutaciones[(numero_cruces+1):M]
    # Índices de las parejas que se van a cruzar
    índices_cruce = [
        [2*i-1, 2*i] for i in parejas_a_cruzar
    ]
    # Indice de los elementos que no se cruzan
    índices_no_cruce = vcat([
        2*i-1 for i in parejas_que_no_se_cruzan
    ],[
        2*i for i in parejas_que_no_se_cruzan
    ])
    ### Mutación
    cantidad_cromosomas_a_mutar = round(Int, probabilidad_mutación* numero_cromosomas_por_generación)
    índices_a_mutar = randperm(M)[1:cantidad_cromosomas_a_mutar]

    # Variable que acumula el número de evaluaciones de la función de evaluación
    evaluaciones = numero_cromosomas_por_generación
    # Array que contendrá a los elementos seleccionados
    Seleccionados = generación
    numero_generacion = 1
    while evaluaciones < evaluaciones_máximas_FE
        # Paso 1: Selección (torneo binario)
        índices_seleccionados = TorneoBinario(evaluaciones_función_evaluación, numero_cromosomas_por_generación)

        # Paso 2: Cruce
        Seleccionados = [
            x
            for v in índices_cruce
                # Devuelve dos elementos del cruce
                for x in
                AuxCrucePorIndices(
                    índices_seleccionados[v[1]],
                    índices_seleccionados[v[2]],
                    generación,
                    función_cruce
                )
        ]
        Seleccionados = vcat(Seleccionados, [generación[i] for i in índices_no_cruce])

        # Paso 3: Mutación
        for i in índices_a_mutar
            Seleccionados[i] = función_mutación(Seleccionados[i])
        end
        # Paso 4: Reemplazo
        # Calculamos mejor cromosoma generación anterior
        índice_mejor = argmax(evaluaciones_función_evaluación)
        mejor_cromosa = generación[índice_mejor]
        evaluación_mejor_cromosoma = evaluaciones_función_evaluación[índice_mejor]

        generación = copy(Seleccionados) # Por tratarse de algoritmo generacional se reemplaza totalmente la población
        # Calculas la evaluación de la generación nueva
        evaluaciones_función_evaluación = map(x -> función_evaluación(x), generación)
        evaluaciones += numero_cromosomas_por_generación

        # para mantener el elitismo si no está el mejor de la generación lo añadimos
        if !(mejor_cromosa in generación)
            índice_peor = argmin(evaluaciones_función_evaluación)
            generación[índice_peor] = mejor_cromosa
            evaluaciones_función_evaluación[índice_peor] = evaluación_mejor_cromosoma
        end
        numero_generacion += 1

        # Paso 5: se comprueba si hay que aplicar la búsqueda local
        if(numero_generacion % numero_generaciones_aplicar_bl == 0)
            for i in funcion_selección_índices_para_busqueda_local(numero_cromosomas_por_generación,
                porcentaje_cromosomas_busqueda_local,
                evaluaciones_función_evaluación)
                generación[i], evaluaciones_función_evaluación[i] = PrimeroMejor(generación[i],
                                2*tamaño_cromosoma,
                                función_evaluación)
                evaluaciones += 2*tamaño_cromosoma
            end
        end
    end

    # Devolvemos el mejor cromosoma encontrado
    índice_mejor = argmax(evaluaciones_función_evaluación)
    evaluación_mejor = evaluaciones_función_evaluación[índice_mejor]
    mejor_cromosa = generación[índice_mejor]
    return mejor_cromosa, evaluación_mejor
end

function Memetico_LearnerOneNN(data::Matrix{<:Real},labels,
    evaluaciones_máximas_FE,
    numero_cromosomas_por_generación,
    probabilidad_cruce,
    probabilidad_mutación,
    tamaño_cromosoma,
    función_cruce,
    función_mutación,
    porcentaje_cromosomas_busqueda_local,
    funcion_selección_índices_para_busqueda_local,
    numero_generaciones_aplicar_bl = 10,
    )
     # Creamos función objetivo
    a = 0.5
    umbral_tasa_reduccion = 0.1
    función_evaluación = CrearFuncionObjetivo(data, labels, a, umbral_tasa_reduccion )

    # Buscamos pesos de acorde al algoritmo Memetico
    w, f_w = Memetico(evaluaciones_máximas_FE,
                        numero_cromosomas_por_generación,
                        probabilidad_cruce,
                        probabilidad_mutación,
                        tamaño_cromosoma,
                        función_evaluación,
                        función_cruce,
                        función_mutación,
                        porcentaje_cromosomas_busqueda_local,
                        funcion_selección_índices_para_busqueda_local,
                        numero_generaciones_aplicar_bl,
                        )

    return WeightedLearnerEuclideanOneNN(w, data, labels), f_w, w

end


text:
module Simrel

include("utilities.jl")

end # module


text:
@testset "magic.jl" begin
    @printf("> magic\n")

    [@test length(Bb.getOccupancyString(i)) == 2^i for i = 1:6]
    @test Bb.getOccupancyString(1) == ("0", "1")

    [@test Bb.getNumberOfSquares(p2u["a1"], true, d) == 6 for d in [1, 3]]
    [@test Bb.getNumberOfSquares(p2u["a1"], true, d) == 0 for d in [2, 4]]
    @test Bb.getNumberOfSquares(p2u["a1"], false, 2) == 6
    [@test Bb.getNumberOfSquares(p2u["a1"], false, d) == 0 for d in [1, 3, 4]]

    [@test Bb.getNumberOfBits(p2u[s], true) == 12 for s in ["a1", "a8", "h1", "h8"]]
    [@test Bb.getNumberOfBits(p2u[s], true) == 11 for s in ["a4", "f8", "e1", "h2"]]
    [@test Bb.getNumberOfBits(p2u[s], true) == 10 for s in ["b4", "d3", "f7", "g2"]]
    [@test Bb.getNumberOfBits(p2u[s], false) == 6 for s in ["a1", "a8", "h1", "h8"]]
    [@test Bb.getNumberOfBits(p2u[s], false) == 5 for s in ["a3", "g5", "e2", "f8"]]
    [@test Bb.getNumberOfBits(p2u[s], false) == 7 for s in ["c3", "f5", "e3"]]
    [@test Bb.getNumberOfBits(p2u[s], false) == 9 for s in ["d4", "e5"]]
    [@test Bb.ROOK_BITS[p2u[s]] == 12 for s in ["a1", "a8", "h1", "h8"]]
    [@test Bb.ROOK_BITS[p2u[s]] == 11 for s in ["a4", "f8", "e1", "h2"]]
    [@test Bb.ROOK_BITS[p2u[s]] == 10 for s in ["b4", "d3", "f7", "g2"]]
    [@test Bb.BISHOP_BITS[p2u[s]] == 6 for s in ["a1", "a8", "h1", "h8"]]
    [@test Bb.BISHOP_BITS[p2u[s]] == 5 for s in ["a3", "g5", "e2", "f8"]]
    [@test Bb.BISHOP_BITS[p2u[s]] == 7 for s in ["c3", "f5", "e3"]]
    [@test Bb.BISHOP_BITS[p2u[s]] == 9 for s in ["d4", "e5"]]

    @test Base.summarysize(Bb.DIAGO_OCCS) == 49080
    @test Base.summarysize(Bb.ORTHO_OCCS) == 826296
    @test (Bb.MASK_RANKS[1] & ~p2u["h1"] ⊻ Bb.MASK_FILES[1] & ~p2u["a8"]) in Bb.ORTHO_OCCS[p2u["a1"]]

    @test Bb.randomMagic() != Bb.randomMagic()

    blockers = reduce(|, p2u[s] for s in ["e2", "e5", "b3", "d4", "h3", "h6"])
    rook_attack = reduce(|, p2u[s] for s in ["e2", "e4", "e5", "b3", "c3", "d3", "f3", "g3", "h3"])
    bishop_attack = reduce(|, p2u[s] for s in ["d4", "f2", "g1", "f4", "g5", "h6", "d2", "c1"])
    @test Bb.getMagicAttack(p2u["e3"], blockers, true) == rook_attack
    @test Bb.getMagicAttack(p2u["e3"], blockers, false) == bishop_attack
end


text:
#=
An optimization oracle is a linear optimization oracle for the set S, i.e., it solves
    z^\ast(c) :=    min c^Tw
                    s.t. w \in S
The solution of the above problem returned by the oracle is referred to as w^\ast(\cdot) in
the paper

A specific implementation of an optimization oracle should take a single input argument c,
and should return the pair (z^\ast(c), w^\ast(c))

This file contains the code for the portfolio problem optimization oracle, which
relies on the JuMP and Gurobi packages.
=#

using JuMP, Gurobi, LinearAlgebra

"""
    portfolio_simplex_oracle_jump_basic(c, Sigma, gamma)

Optimization oracle in JuMP for the portfolio problem:
    min  c^T w
    s.t. w^T Sigma w <= gamma
         e^T w <= 1
         w >= 0
Here c is equal to minus the vector of returns, Sigma is the covariance matirx of the returns,
and gamma is the desired risk level in terms of the overall variance of the portfolio. Short
selling is not allowed.
"""
function portfolio_simplex_oracle_jump_basic(c::Vector{Float64}, Sigma::Matrix{Float64}, gamma::Float64)
    d = length(c)

    mod = Model(with_optimizer(Gurobi.Optimizer, OutputFlag = 0))
    @variable(mod, w[1:d] >= 0)
    @constraint(mod, sum(w[i] for i = 1:d) <= 1)
    @constraint(mod, w'*Sigma*w <= gamma)

    @objective(mod, Min, dot(c, w))

    optimize!(mod)
    z_ast = objective_value(mod)
    w_ast = value.(w)

    return (z_ast, w_ast)
end

"""
Smarter way of implementing the oracle in JuMP. Constructs the feasible region from
(Sigma, gamma) ahead of time and then returns an oracle that is just a function of c.
"""
function portfolio_simplex_jump_setup(Sigma::Matrix{Float64}, gamma::Float64; gurobiEnv = Gurobi.Env())
    (d, d2) = size(Sigma)

    if d != d2
        error("Sigma dimensions don't match")
    end

    mod = Model(with_optimizer(Gurobi.Optimizer, gurobiEnv))
    @variable(mod, w[1:d] >= 0)
    @constraint(mod, sum(w[i] for i = 1:d) <= 1)
    @constraint(mod, w'*Sigma*w <= gamma)

    function local_portfolio_oracle(c::Vector{Float64})
        @objective(mod, Min, dot(c, w))

        optimize!(mod)
        z_ast = objective_value(mod)
        w_ast = value.(w)

        return (z_ast, w_ast)
    end

    return c -> local_portfolio_oracle(c)
end

max_stars_repo_name: paulgrigas/SmartPredictThenOptimize

max_issues_repo_name: paulgrigas/SmartPredictThenOptimize

max_forks_repo_name: paulgrigas/SmartPredictThenOptimize