text:
#####
##### Transactions
#####

"""
    SidePot

A side-pot for player `seat_number`,
who has gone all-in with amount `amt`.
"""
mutable struct SidePot
  seat_number::Int
  amt::Float64
  cap::Float64 # total amount any individual player can contribute to this side-pot
end
seat_number(sp::SidePot) = sp.seat_number
amount(sp::SidePot) = sp.amt
cap(sp::SidePot) = sp.cap

"""
    TransactionManager

Handle pots and side pots
among multiple players.
"""
struct TransactionManager
  sorted_players::Vector{Player}
  initial_brs::Vector{Float64}
  pot_id::Union{Nothing,Vector{Int}}
  side_pots::Vector{SidePot}
  unsorted_to_sorted_map::Vector{Int}
end

function Base.show(io::IO, tm::TransactionManager, include_type = true)
    include_type && println(io, typeof(tm))
    println(io, "Pot(s)           = $(tm.side_pots)")
end

function TransactionManager(players)
    sorted_players = sort(collect(players); by = x->bank_roll(x))

    cap = zeros(length(players))
    for i in 1:length(players)
        if i == 1
            cap[i] = bank_roll(sorted_players[i])
        else
            cap[i] = bank_roll(sorted_players[i]) - sum(cap[1:i-1])
        end
    end

    unsorted_to_sorted_map = collect(map(players) do player
        findfirst(seat_number.(sorted_players) .== seat_number(player))
    end)
    side_pots = [SidePot(sn, 0, cap_i) for (cap_i, sn, amt) in zip(cap, seat_number.(sorted_players), bank_roll.(sorted_players))]

    TransactionManager(
        sorted_players,
        deepcopy(collect(bank_roll.(players))),
        Int[1],
        side_pots,
        unsorted_to_sorted_map,
    )
end
amount(tm::TransactionManager) = amount(tm.side_pots[tm.pot_id[1]])
cap(tm::TransactionManager) = cap(tm.side_pots[tm.pot_id[1]])

function last_action_of_round(table, player, call)
    laor = all_oppononents_all_in(table, player) || (count(action_required.(players_at_table(table))) == 0 && call)
    @debug "last_action_of_round = $laor"
    # @debug "    action_required.(players_at_table(table)) = $(action_required.(players_at_table(table)))"
    # @debug "    all_oppononents_all_in(table, player) = $(all_oppononents_all_in(table, player))"
    # @debug "    call = $call"
    return laor
end

"""
    contribute!(table, player, amt, call=false)

Player `player` contributes amount `amt` to the
appropriate side-pot(s) for calls (for `call = true`)
and raises (for `call = false`). All side-pots are
internally handled.

To properly handle side-pots, each raise and call
is decomposed into side-pots based on the bank roll
of the (sorted) players at the start of the game.

```
                      [1300] __ all-in for sorted_player_6
                    | []
                    | [1100] __ all-in for sorted_player_5
                |   | []
                |   | []
                |   | []
                |   | []
                |   | []
                |   | [500] __ all-in for sorted_player_4
            |   |   | [400] __ all-in for sorted_player_3
        |   |   |   | [300] __ all-in for sorted_player_2
    |   |   |   |   | []
    |   |   |   |   | [100] __ all-in for sorted_player_1
|   |   |   |   |   | []
---------------------
1   2   3   4   5   6
```
"""
function contribute!(table, player, amt, call=false)
    tm = table.transactions
    @debug "$(name(player))'s bank roll (pre-contribute) = \$$(bank_roll(player))"
    if !(0 ≤ amt ≤ bank_roll(player))
        msg1 = "$(name(player)) has insufficient bank"
        msg2 = "roll (\$$(bank_roll(player))) to add \$$amt to pot."
        error(msg1*msg2)
    end
    @assert all_in(player) == false
    @assert !(amt ≈ 0) "Cannot contribute \$$amt to the pot!"

    player.round_contribution += amt
    player.pot_investment += amt
    amt_remaining = amt
    @debug "$(name(player)) contributing $amt to side-pots"
    @debug "caps = $(cap.(tm.side_pots))"
    @debug "$(name(player))'s pot_investment = $(player.pot_investment)"

    for i in 1:length(tm.side_pots)
        @assert 0 ≤ amt_remaining
        cap_i = cap(tm.side_pots[i])
        sp_amt = amount(tm.side_pots[i])
        cond = amt_remaining < cap_i
        amt_contrib = cond ? amt_remaining : cap_i
        contributing = !side_pot_full(tm, i) && !(cap_i ≈ 0)
        # This is a bit noisy:
        # @debug "$(name(player)) potentially contributing $amt_contrib to side-pot $(i) ($cond). cap_i=$cap_i, amt_remaining=$amt_remaining"
        @debug "contributing = $contributing"
        contributing || continue
        @assert !(amt_contrib ≈ 0)
        tm.side_pots[i].amt += amt_contrib
        player.bank_roll -= amt_contrib
        amt_remaining -= amt_contrib
        amt_remaining ≈ 0 && break
    end
    @assert amt_remaining ≈ 0 # pots better be emptied

    if bank_roll(player) ≈ 0 # went all-in, set exactly.
        player.all_in = true
        player.action_required = false
        player.bank_roll = 0
    end

    if is_side_pot_full(tm, table, player, call)
        set_side_pot_full!(tm)
    end
    @debug "$(name(player))'s bank roll (post-contribute) = \$$(bank_roll(player))"
    @debug "all_in($(name(player))) = $(all_in(player))"
end

function is_side_pot_full(tm::TransactionManager, table, player, call)
    players = players_at_table(table)
    # To switch from pot_id = 1 to pot_id = 2, then exactly 1 player  should be all-in:
    # To switch from pot_id = 2 to pot_id = 3, then exactly 2 players should be all-in:
    # ...
    return last_action_of_round(table, player, call) && count(all_in.(players)) == tm.pot_id[1]
end

set_side_pot_full!(tm::TransactionManager) = (tm.pot_id[1]+=1)
side_pot_full(tm::TransactionManager, i) = i < tm.pot_id[1]

sidepot_winnings(tm::TransactionManager, id::Int) = sum(map(x->x.amt, tm.side_pots[1:id]))

function distribute_winnings_1_player_left!(players, tm::TransactionManager, table_cards)
    @assert count(still_playing.(players)) == 1
    n = length(tm.side_pots)
    for (player, initial_br) in zip(players, tm.initial_brs)
        not_playing(player) && continue
        amt_contributed = initial_br - bank_roll(player)
        ∑spw = sidepot_winnings(tm, n)
        prof = ∑spw-amt_contributed
        @info "$(name(player)) wins \$$(∑spw) (\$$(prof) profit) (all opponents folded)"
        player.bank_roll += ∑spw
        for j in 1:n
            tm.side_pots[j].amt = 0 # empty out distributed winnings
        end
        break
    end
    return nothing
end

function distribute_winnings!(players, tm::TransactionManager, table_cards)
    @debug "Distributing winnings..."
    @debug "Pot amounts = $(amount.(tm.side_pots))"
    if count(still_playing.(players)) == 1
        return distribute_winnings_1_player_left!(players, tm, table_cards)
    end
    hand_evals_sorted = map(enumerate(tm.sorted_players)) do (ssn, player)
        fhe = inactive(player) ? nothing : FullHandEval((player.cards..., table_cards...))
        eligible = still_playing(player)
        (; eligible=eligible, player=player, fhe=fhe, ssn=ssn)
    end

    side_pot_winnings = map(players) do player
        zeros(length(tm.side_pots))
    end
    winning_hands = Vector{Symbol}(undef, length(players))
    all_winning_players = Player[]

    for i in 1:length(tm.side_pots)
        sidepot_winnings(tm, length(players)) ≈ 0 && continue # no money left to distribute

        hand_evals_sorted = map(hand_evals_sorted) do (eligible, player, fhe, ssn)
            if ssn ≥ i # only (sorted) players i:end can win the ith side-pot
                (; eligible=eligible, player=player, fhe=fhe, ssn=ssn)
            else
                (; eligible=false, player=player, fhe=fhe, ssn=ssn)
            end
        end
        @debug begin
            s = "Sorted hand evals: \n"
            for hes in hand_evals_sorted
                inactive(hes.player) && continue # (don't have cards to eval)
                s *= "eligible=$(hes.eligible), "
                s *= "sn=$(seat_number(hes.player)), "
                s *= "ssn=$(hes.ssn), "
                s *= "hr=$(hand_rank(hes.fhe)), "
                s *= "ht=$(hand_type(hes.fhe))\n"
            end
            s
        end
        all_valid_min_hrs = min(map(x->hand_rank(x.fhe), filter(x->x.eligible, hand_evals_sorted))...)

        winner_ids = findall(hand_evals_sorted) do x
            x.eligible ? hand_rank(x.fhe)==all_valid_min_hrs : false
        end
        n_winners = length(winner_ids)
        @debug "winner_ids = $(winner_ids)"
        @debug "length(players) = $(length(players))"
        @debug "length(hand_evals_sorted) = $(length(hand_evals_sorted))"
        @debug "length(tm.sorted_players) = $(length(tm.sorted_players))"
        for winner_id in winner_ids
            win_seat = seat_number(tm.sorted_players[winner_id])
            winning_player = players[win_seat]
            not_playing(winning_player) && continue
            push!(all_winning_players, winning_player)
            amt = sidepot_winnings(tm, i) / n_winners
            side_pot_winnings[win_seat][i] = amt
            winning_hands[win_seat] = hand_type(hand_evals_sorted[winner_id].fhe)
        end
        for j in 1:i
            tm.side_pots[j].amt = 0 # empty out distributed winnings
        end
    end
    for (player, initial_br, player_winnings) in zip(players, tm.initial_brs, side_pot_winnings)
        ∑spw = sum(player_winnings)
        ssn = tm.unsorted_to_sorted_map[seat_number(player)]
        if ∑spw ≈ 0
            if active(player)
                hand = hand_evals_sorted[ssn].fhe
                hand_name = hand_type(hand)
                bc = best_cards(hand)
                f_str = folded(player) ? " (folded)" : ""
                @info "$(name(player)) loses$f_str \$$(pot_investment(player)) with $bc ($hand_name)!"
            end
        else
            hand = hand_evals_sorted[ssn].fhe
            hand_name = hand_type(hand)
            bc = best_cards(hand)
            amt_contributed = initial_br - bank_roll(player)
            @debug "$(name(player))'s side-pot wins: \$$(player_winnings)!"
            prof = ∑spw-amt_contributed
            @info "$(name(player)) wins \$$∑spw (\$$prof profit) with $bc ($hand_name)!"
            player.bank_roll += ∑spw
        end
    end

    @debug "Distributed winnings..."
    return Tuple(all_winning_players)
end


text:
using ColorSchemeTools

# code in here is for translating colorschemes from
# https://github.com/matplotlib/matplotlib/blob/master/lib/matplotlib/_cm.py
# into ColorSchemeTools.jl compatible code

# # Gnuplot palette functions
# these are simply translated from Python
_g0(x) = 0
_g1(x) = 0.5
_g2(x) = 1
_g3(x) = x # really
_g4(x) = x^2
_g5(x) = x^3
_g6(x) = x^4
_g7(x) = sqrt(x)
_g8(x) = sqrt(sqrt(x))
_g9(x) = sin(x * π/2)
_g10(x) = cos(x * π/2)
_g11(x) = abs(x - 0.5)
_g12(x) = (2x - 1)^2
_g13(x) = sin(x * π)
_g14(x) = abs(cos(x * π))
_g15(x) = sin(x * 2π)
_g16(x) = cos(x * 2π)
_g17(x) = abs(sin(x * 2π))
_g18(x) = abs(cos(x * 2π))
_g19(x) = abs(sin(x * 4π))
_g20(x) = abs(cos(x * 4π))
_g21(x) = 3x
_g22(x) = 3x - 1
_g23(x) = 3x - 2
_g24(x) = abs(3x - 1)
_g25(x) = abs(3x - 2)
_g26(x) = (3x - 1) / 2
_g27(x) = (3x - 2) / 2
_g28(x) = abs((3x - 1) / 2)
_g29(x) = abs((3x - 2) / 2)
_g30(x) = x / 0.32 - 0.78125
_g31(x) = 2x - 0.84
function _g32(x) # andreas
   ret = zeros(length(x))
       m = (x < 0.25)
       ret[m] = 4 * x[m]
       m = (x >= 0.25) & (x < 0.92)
       ret[m] = -2 * x[m] + 1.84
       m = (x >= 0.92)
       ret[m] = x[m] / 0.08 - 11.5
       return
end
function _g32(x) # kristoffer
    ret = zeros(length(x))
    m = x .< 0.25
    ret[m] = 4 * x[m]
    m = (x .>= 0.25) .& (x .< 0.92)
    ret[m] = -2 .* x[m] .+ 1.84
    m = (x .>= 0.92)
    ret[m] = x[m] ./ 0.08 .- 11.5
    return ret
end

_g33(x) =  abs(2x - 0.5)
_g34(x) =  2x
_g35(x) =  2x - 0.5
_g36(x) =  2x - 1

# Didn't bother with building these functions into gfunc()s
# just pass the functions directly

#
# _gnuplot_data = Dict(
#         :red => gfunc(7),
#         :green => gfunc(5),
#         :blue => gfunc(15))
# =>
gnuplot = make_colorscheme(_g7, _g5, _g15)

#
# _gnuplot2_data = Dict(
#         :red => gfunc(30),
#         :green => gfunc(31),
#         :blue => gfunc(32))
# =>
# this one is the only weird one, so I fudged it thus
wackyblues = _g32(range(0, stop=1, length=100))
gnuplot2 = make_colorscheme(_g30, _g31, (n) -> wackyblues[convert(Int, floor(ColorSchemeTools.lerp(n, 0, 1, 1, 100)))])

#
# _ocean_data = Dict(
#         :red => gfunc(23),
#         :green => gfunc(28),
#         :blue => gfunc(3))
# =>
ocean = make_colorscheme(_g23, _g28, _g3)

#
# _afmhot_data = Dict(
#         :red => gfunc(34),
#         :green => gfunc(35),
#         :blue => gfunc(36))
# =>
afmhot = make_colorscheme(_g34, _g35, _g36)

#
# _rainbow_data = Dict(
#         :red => gfunc(33),
#         :green => gfunc(13),
#         :blue => gfunc(10))
# =>
rainbow = make_colorscheme(_g33, _g13, _g10)

#
# _gist_gray_data = {
#         'red': gfunc[3],
#         'green': gfunc[3],
#         'blue': gfunc[3],
# }
# don't really get this next bit, but it's easy to avoid
# def _gist_yarg(x): return 1 - x
#_gist_yarg_data = {'red': _gist_yarg, 'green': _gist_yarg, 'blue': _gist_yarg} ??
# =>

gist_gray = make_colorscheme(identity, identity, identity)
# tada!
gist_yarg = ColorScheme(reverse(gist_gray.colors))

#
# def _gist_heat_red(x): return 1.5 * x
# def _gist_heat_green(x): return 2 * x - 1
# def _gist_heat_blue(x): return 4 * x - 3
# _gist_heat_data = {
#     'red': _gist_heat_red, 'green': _gist_heat_green, 'blue': _gist_heat_blue}
# =>
gist_heat = make_colorscheme((n) -> 1.5n, (n) -> 2n - 1, (n) -> 4n - 3)

#
coolwarm_dict = Dict(
    :red => (
        (0.0, 0.2298057, 0.2298057),
        (0.03125, 0.26623388, 0.26623388),
        (0.0625, 0.30386891, 0.30386891),
        (0.09375, 0.342804478, 0.342804478),
        (0.125, 0.38301334, 0.38301334),
        (0.15625, 0.424369608, 0.424369608),
        (0.1875, 0.46666708, 0.46666708),
        (0.21875, 0.509635204, 0.509635204),
        (0.25, 0.552953156, 0.552953156),
        (0.28125, 0.596262162, 0.596262162),
        (0.3125, 0.639176211, 0.639176211),
        (0.34375, 0.681291281, 0.681291281),
        (0.375, 0.722193294, 0.722193294),
        (0.40625, 0.761464949, 0.761464949),
        (0.4375, 0.798691636, 0.798691636),
        (0.46875, 0.833466556, 0.833466556),
        (0.5, 0.865395197, 0.865395197),
        (0.53125, 0.897787179, 0.897787179),
        (0.5625, 0.924127593, 0.924127593),
        (0.59375, 0.944468518, 0.944468518),
        (0.625, 0.958852946, 0.958852946),
        (0.65625, 0.96732803, 0.96732803),
        (0.6875, 0.969954137, 0.969954137),
        (0.71875, 0.966811177, 0.966811177),
        (0.75, 0.958003065, 0.958003065),
        (0.78125, 0.943660866, 0.943660866),
        (0.8125, 0.923944917, 0.923944917),
        (0.84375, 0.89904617, 0.89904617),
        (0.875, 0.869186849, 0.869186849),
        (0.90625, 0.834620542, 0.834620542),
        (0.9375, 0.795631745, 0.795631745),
        (0.96875, 0.752534934, 0.752534934),
        (1.0, 0.705673158, 0.705673158)),
    :green => (
        (0.0, 0.298717966, 0.298717966),
        (0.03125, 0.353094838, 0.353094838),
        (0.0625, 0.406535296, 0.406535296),
        (0.09375, 0.458757618, 0.458757618),
        (0.125, 0.50941904, 0.50941904),
        (0.15625, 0.558148092, 0.558148092),
        (0.1875, 0.604562568, 0.604562568),
        (0.21875, 0.648280772, 0.648280772),
        (0.25, 0.688929332, 0.688929332),
        (0.28125, 0.726149107, 0.726149107),
        (0.3125, 0.759599947, 0.759599947),
        (0.34375, 0.788964712, 0.788964712),
        (0.375, 0.813952739, 0.813952739),
        (0.40625, 0.834302879, 0.834302879),
        (0.4375, 0.849786142, 0.849786142),
        (0.46875, 0.860207984, 0.860207984),
        (0.5, 0.86541021, 0.86541021),
        (0.53125, 0.848937047, 0.848937047),
        (0.5625, 0.827384882, 0.827384882),
        (0.59375, 0.800927443, 0.800927443),
        (0.625, 0.769767752, 0.769767752),
        (0.65625, 0.734132809, 0.734132809),
        (0.6875, 0.694266682, 0.694266682),
        (0.71875, 0.650421156, 0.650421156),
        (0.75, 0.602842431, 0.602842431),
        (0.78125, 0.551750968, 0.551750968),
        (0.8125, 0.49730856, 0.49730856),
        (0.84375, 0.439559467, 0.439559467),
        (0.875, 0.378313092, 0.378313092),
        (0.90625, 0.312874446, 0.312874446),
        (0.9375, 0.24128379, 0.24128379),
        (0.96875, 0.157246067, 0.157246067),
        (1.0, 0.01555616, 0.01555616)),
    :blue => (
        (0.0, 0.753683153, 0.753683153),
        (0.03125, 0.801466763, 0.801466763),
        (0.0625, 0.84495867, 0.84495867),
        (0.09375, 0.883725899, 0.883725899),
        (0.125, 0.917387822, 0.917387822),
        (0.15625, 0.945619588, 0.945619588),
        (0.1875, 0.968154911, 0.968154911),
        (0.21875, 0.98478814, 0.98478814),
        (0.25, 0.995375608, 0.995375608),
        (0.28125, 0.999836203, 0.999836203),
        (0.3125, 0.998151185, 0.998151185),
        (0.34375, 0.990363227, 0.990363227),
        (0.375, 0.976574709, 0.976574709),
        (0.40625, 0.956945269, 0.956945269),
        (0.4375, 0.931688648, 0.931688648),
        (0.46875, 0.901068838, 0.901068838),
        (0.5, 0.865395561, 0.865395561),
        (0.53125, 0.820880546, 0.820880546),
        (0.5625, 0.774508472, 0.774508472),
        (0.59375, 0.726736146, 0.726736146),
        (0.625, 0.678007945, 0.678007945),
        (0.65625, 0.628751763, 0.628751763),
        (0.6875, 0.579375448, 0.579375448),
        (0.71875, 0.530263762, 0.530263762),
        (0.75, 0.481775914, 0.481775914),
        (0.78125, 0.434243684, 0.434243684),
        (0.8125, 0.387970225, 0.387970225),
        (0.84375, 0.343229596, 0.343229596),
        (0.875, 0.300267182, 0.300267182),
        (0.90625, 0.259301199, 0.259301199),
        (0.9375, 0.220525627, 0.220525627),
        (0.96875, 0.184115123, 0.184115123),
        (1.0, 0.150232812, 0.150232812))
        )

coolwarm = make_colorscheme(coolwarm_dict, notes="blue to white to red")

# _bwr_data = ((0.0, 0.0, 1.0), (1.0, 1.0, 1.0), (1.0, 0.0, 0.0))
# this is another blue to white to red
# _brg_data = ((0.0, 0.0, 1.0), (1.0, 0.0, 0.0), (0.0, 1.0, 0.0))
# this is blue to red to green

bwr_data = (
     (0.0,     (0., 0., 1.)),
     (0.5,     (1., 1., 1.)),
     (1.0,     (1., 0., 0.))
     )

bwr = make_colorscheme(bwr_data)

brg_data  = (
      (0.0,     (0., 0., 1.)),
      (0.5,     (1., 0., 0.)),
      (1.0,     (0., 1., 0.))
      )

brg = make_colorscheme(brg_data)

#
# def _flag_red(x): return 0.75 *  np.sin((x * 31.5 + 0.25) * np.pi) + 0.5
# def _flag_green(x): return       np.sin( x * 31.5         * np.pi)
# def _flag_blue(x): return 0.75 * np.sin((x * 31.5 - 0.25) * np.pi) + 0.5
# _flag_data = {'red': _flag_red, 'green': _flag_green, 'blue': _flag_blue}
# they apparently want red->white->blue->black 16 times
# =>

flag = make_colorscheme(
    (n) -> 0.75 * sin((31.5n + 0.25) * π) + 0.5,
    (n) -> sin(31.5n * π),
    (n) -> 0.75 * sin((31.5n - 0.25) * π) + 0.5,
    # length=300 # looks better with more samples, but probably it will be used for continuous sampling
    )

#
# def _prism_red(x): return 0.75   * np.sin((x * 20.9 + 0.25) * np.pi) + 0.67
# def _prism_green(x): return 0.75 * np.sin((x * 20.9 - 0.25) * np.pi) + 0.33
# def _prism_blue(x): return -1.1  * np.sin((x * 20.9)        * np.pi)
# _prism_data = {'red': _prism_red, 'green': _prism_green, 'blue': _prism_blue}
# The (9) colors of the spectrum repeating 11 times...
# =>

prism = make_colorscheme(
    (n) -> 0.75 * sin((20.9n + 0.25) * π) + 0.67,
    (n) -> 0.75 * sin((20.9n - 0.25) * π) + 0.33,
    (n) -> -1.1 * sin(20.9n * π))

#=
# code to output these colorschemes to files

for cs in (:rainbow, :afmhot, :ocean, :gnuplot, :gnuplot2, :gist_gray, :gist_yarg, :gist_heat, :coolwarm, :bwr, :brg, :flag, :prism)
     colorscheme_to_text(Base.eval(Main, cs), String(cs), "/tmp/$(cs).jl")
 end

# assemble
# who needs Unix tools...
open("/tmp/out.jl", "w") do file
     for f in ("/tmp/rainbow.jl",  "/tmp/prism.jl",  "/tmp/flag.jl", "/tmp/brg.jl",  "/tmp/bwr.jl",  "/tmp/coolwarm.jl",  "/tmp/gist_heat.jl", "/tmp/gist_yarg.jl",  "/tmp/gist_gray.jl",  "/tmp/gnuplot2.jl", "/tmp/gnuplot.jl" , "/tmp/ocean.jl",  "/tmp/afmhot.jl")
        write(file, read(f, String))
        write(file, "\n\n")
    end
end
=#

ext: jl


text:
struct Cache
    words::Set{String}
    synonyms::Dict{String, Set{String}}
    words_by_anagram::Dict{String, Vector{String}}
    abbreviations::Dict{String, Set{String}}
    substrings::PTrie{30}
    prefixes::PTrie{30}
    indicators::Dict{GrammaticalSymbol, Set{String}}
end

function parse_heading(line)
    parts = split(line, '|')
    length(parts) == 2 || @show line parts
    word = parts[1]
    num_entries = parse(Int, parts[2])
    word, num_entries
end

function parse_entry(line)
    normalize.(split(line, '|')[2:end])
end

function add_synonyms!(synonyms, word, list)
    entries = get!(() -> Set{String}(), synonyms, word)
    for entry in list
        push!(entries, entry)
    end
end

function parse_synonyms(fname)
    @info("Parsing synonyms file...")
    synonyms = Dict{String, Set{String}}()
    open(fname) do file
        @assert readline(file) == "ISO8859-1"
        while true
            line = readline(file)
            if isempty(line)
                break
            end
            word, num_entries = parse_heading(line)
            for i in 1:num_entries
                add_synonyms!(synonyms, normalize(word), parse_entry(readline(file)))
            end
        end
    end
    @info("Done!")
    synonyms
end

function make_symmetric!(synonyms)
    for (word, entries) in synonyms
        for entry in entries
            push!(get!(() -> Set{String}(), synonyms, entry), word)
        end
    end
end

function remove_loops!(synonyms)
    for (word, entries) in synonyms
        if word in entries
            delete!(entries, word)
        end
    end
end

"""
Remove synonym entries which contain the original word. For example,
"spin" and "spin around" are listed as synonyms in the dataset, but
crossword rules don't allow the clue to contain part of the answer
like that.
"""
function remove_self_mentions!(synonyms)
    for (word, entries) in synonyms
        disallowed = String[]
        for entry in entries
            if occursin(word, entry) || occursin(entry, word)
                push!(disallowed, entry)
            end
        end
        for entry in disallowed
            delete!(entries, entry)
        end
    end
end

function load_synonyms()
    path = datadep"OpenOffice-MyThes-1.0/openoffice-mythes-mirror-1.0/th_en_US_new.dat"
    synonyms = parse_synonyms(path)
    make_symmetric!(synonyms)
    remove_loops!(synonyms)
    remove_self_mentions!(synonyms)
    synonyms
end

function Cache()
    synonyms = load_synonyms()
    words = Set{String}()
    words_by_anagram = Dict{String, Vector{String}}()
    abbreviations = Dict{String, Set{String}}()
    substrings = PTrie{30}()
    prefixes = PTrie{30}()
    indicators = Dict{GrammaticalSymbol, Set{String}}()

    for word in keys(synonyms)
        push!(words, word)
    end
    for line in eachline(datadep"SCOWL-wordlist-en_US-large/en_US-large.txt")
        push!(words, normalize(line))
    end
    for word in words
        key = join(sort(collect(replace(word, " " => ""))))
        v = get!(Vector{String}, words_by_anagram, key)
        push!(v, word)
    end
    open(joinpath(@__DIR__, "..", "corpora", "mhl-abbreviations", "abbreviations.json")) do file
        for (word, abbrevs) in JSON.parse(file)
            abbreviations[normalize(word)] = Set(normalize.(abbrevs))
        end
    end
    open(joinpath(@__DIR__, "..", "corpora", "abbreviations.json")) do file
        for (word, abbrevs) in JSON.parse(file)
            s = get!(Set{String}, abbreviations, normalize(word))
            for a in abbrevs
                push!(s, a)
            end
        end
    end

    @showprogress "Substrings PTrie" for word in words
        push!(prefixes, word)
        word = replace(word, ' ' => "")
        for i in 1:length(word)
            push!(substrings, word[i:end])
        end
    end

    for (filename, part_of_speech) in [
        ("Anagram", AnagramIndicator()),
        ("Filler", Filler()),
        ("FinalSubstring", FinalSubstringIndicator()),
        ("InitialSubstring", InitialSubstringIndicator()),
        ("InsertAB", InsertABIndicator()),
        ("InsertBA", InsertBAIndicator()),
        ("Initials", InitialsIndicator()),
        ("Reversal", ReversalIndicator())]
        for line in eachline(joinpath(@__DIR__, "..", "corpora", "indicators", filename))
            phrase = normalize(strip(line))
            push!(get!(Set{String}, indicators, part_of_speech), phrase)
        end
    end
    Cache(words, synonyms, words_by_anagram, abbreviations, substrings, prefixes, indicators)
end


text:
# This file is a part of MGVI.jl, licensed under the MIT License (MIT).

using Distributions
using DistributionsAD
using LinearAlgebra
import Zygote

"""
    _fi_cov_only_upper_params(s; offset=0)

Return indices of the covariance part of the FI matrix that
correspond to upper triangular part of the covariance matrix.

Specify offset to get indices shifted to the right and bottom by
the same offset. (useful to remove irrelevant cols from full FI)
"""
function _fi_cov_only_upper_params(s; offset=0)
    slices = []
    for i in 1:s
        start = (i-1)*s+1
        push!(slices, (start:start+i-1))
    end
    cov_slices = union(slices...) .+ offset
    relevant_slices = union(1:offset, cov_slices)
    relevant_slices
end

# fisher information MC implementation

function _grad_logpdf(model::Function, params::AbstractVector, data_point)
    g = Zygote.gradient(dp -> logpdf(model(dp), data_point), params)[1]
    g_flat = collect(g)
    g_flat * g_flat'
end

_cut_params(res, dist::Distribution) = res

"""
    _cut_params(res, dist::Distribution)

noop for most of the distributions.

for MvNormal it removes rows and cols from FI that correspond to lower triangular
part of the covariance matrix. Covariance matrix is symmetric, so lower triangular
parameters are redundand.
"""
function _cut_params(res, dist::MvNormal)
    mean_size = length(dist.μ)
    slices = _fi_cov_only_upper_params(mean_size; offset=mean_size)
    res[slices, slices]
end

function fisher_information_mc(model::Function, params::AbstractVector, n::Integer)
    dist = model(params)
    sample() = _grad_logpdf(model, params, rand(dist))/n

    res = [sample() for _ in 1:Threads.nthreads()]
    Threads.@threads for i in 1:n-Threads.nthreads()
        res[Threads.threadid()] .+= sample()
    end

    _cut_params(sum(res), dist)
end

# end fisher information mc

"""
    explicit_mv_normal_fi(variance::AbstractMatrix)

Compute fisher information of MvNormal based on the wiki article:

https://en.wikipedia.org/wiki/Fisher_information#Multivariate_normal_distribution
"""
function explicit_mv_normal_fi(variance::AbstractMatrix)
    var_dims = length(variance)
    mean_dims = size(variance, 1)
    dims = mean_dims + var_dims

    res = zeros(dims, dims)
    inv_var = inv(variance)
    res[1:mean_dims, 1:mean_dims] .= inv_var
    for m in 1:var_dims
        for n in 1:var_dims
            m_mat = zeros(mean_dims, mean_dims)
            m_mat[m] = 1
            m_mat = Symmetric(m_mat)
            n_mat = zeros(mean_dims, mean_dims)
            n_mat[n] = 1
            n_mat = Symmetric(n_mat)
            res[mean_dims + m, mean_dims + n] = tr(inv_var * m_mat * inv_var * n_mat)/2
        end
    end

    slices = _fi_cov_only_upper_params(mean_dims; offset=mean_dims)
    res[slices, slices]
end


text:
module Runtests

import Test
using  ModiaLang

@time Test.@testset verbose=true "Modia3D" begin
    usePlotPackage("SilentNoPlot")

    Test.@testset "Basic" begin
        include(joinpath("Basic", "AllShapes.jl"))
        include(joinpath("Basic", "PendulumWithBar1.jl"))
        include(joinpath("Basic", "PendulumWithBar2.jl"))
        include(joinpath("Basic", "PendulumWithDamper.jl"))
        include(joinpath("Basic", "PendulumWithFix.jl"))
        include(joinpath("Basic", "PendulumWithParameterizedDamper.jl"))
        include(joinpath("Basic", "PendulumWithSpring.jl"))
        include(joinpath("Basic", "DoublePendulumWithDampers.jl"))
        include(joinpath("Basic", "Mobile.jl"))
        include(joinpath("Basic", "ShaftFreeMotion.jl"))
        Test.@test_throws LoadError include(joinpath("Basic", "Object3DWithoutParentError.jl"))
    end

    Test.@testset "Robot" begin
        include(joinpath("Robot", "ServoWithRamp.jl"))
        include(joinpath("Robot", "ServoWithRampAndPrismatic.jl"))
        include(joinpath("Robot", "ServoWithRampAndRevolute.jl"))
        include(joinpath("Robot", "ServoWithPathAndRevolute.jl"))
        include(joinpath("Robot", "YouBotWithSphere.jl"))
        Test.@test_skip include(joinpath("Robot", "YouBotPingPong.jl")) # too long computation time
        include(joinpath("Robot", "YouBotGripping.jl"))
        Test.@test_skip include(joinpath("Robot", "YouBotsGripping.jl")) # too long computation time
    end

    Test.@testset "Collision" begin
        include(joinpath("Collision", "BouncingSphere.jl"))
        include(joinpath("Collision", "BouncingSphereFreeMotion.jl"))
        include(joinpath("Collision", "BouncingEllipsoid.jl"))
        include(joinpath("Collision", "BouncingCones.jl"))
        include(joinpath("Collision", "BouncingCapsules.jl"))
        include(joinpath("Collision", "BouncingBeams.jl"))
        include(joinpath("Collision", "TwoCollidingBalls.jl"))
        include(joinpath("Collision", "TwoCollidingBoxes.jl"))
        include(joinpath("Collision", "CollidingCylinders.jl"))
        include(joinpath("Collision", "CollidingSphereWithBunnies.jl"))
        include(joinpath("Collision", "NewtonsCradle.jl"))
        include(joinpath("Collision", "Billard4Balls.jl"))
        Test.@test_skip include(joinpath("Collision", "Billard16Balls.jl"))  # too long computation time
    end

    Test.@testset "Tutorial" begin
        include(joinpath("Tutorial", "Pendulum1.jl"))
        include(joinpath("Tutorial", "Pendulum2.jl"))
        include(joinpath("Tutorial", "Pendulum3.jl"))
        include(joinpath("Tutorial", "BouncingSphere.jl"))
    end

    Test.@testset "old" begin
        include(joinpath("old", "Move_Pendulum.jl"))
        include(joinpath("old", "Plot_cor.jl"))
        include(joinpath("old", "Plot_SlidingFriction.jl"))
        include(joinpath("old", "test_Shapes.jl"))
        include(joinpath("old", "Test_PathPlanning.jl"))
        include(joinpath("old", "test_Solids.jl"))
        include(joinpath("old", "Visualize_Beam.jl"))
    end

    usePreviousPlotPackage()
end

end


text:
#=
  Copyright (c) 2021
  https://github.com/Ismael-VC/Varvara.jl/blob/main/CONTRIBUTORS.md

  Permission to use, copy, modify, and distribute this software for any
  purpose with or without fee is hereby granted, provided that the above
  copyright notice and this permission notice appear in all copies.

  THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
  WITH REGARD TO THIS SOFTWARE.
=#

module VarvaraBIOS

end  # module


text:
module relatedness_module

using Distributed
using SharedArrays
using Statistics
using DelimitedFiles
include("sync_processing_module.jl")

#####################
### Sub-functions ###
#####################
### [Weir and Cockerham, 1984 method](https://www.jstor.org/stable/2408641?seq=1)
function Fst_weir1984_func(sync::Array{Any,2}, pool_sizes::Array{Int64,1}, window_size::Int64)
    CHROM=[]; WINDOW=[]; FST=[]
    for chrom in unique(sync[:,1])
        # test
        # chrom = 1
        subsync = sync[sync[:,1] .== chrom, :]
        nWindows = convert(Int, ceil(maximum(subsync[:,2]) / window_size))
        for window in 1:nWindows
            # window = 1
            # println(chrom)
            # println(window)
            # println("#######################")
            _start_ = window_size*(window-1) + 1
            _end_ = window_size*(window)
            windowsync = subsync[(subsync[:,2] .>= _start_) .& (subsync[:,2] .<= _end_), :]
            if (size(windowsync,1) != 0)
                l = size(windowsync,1)     #number of loci in this window
                a = 6                               #number of alleles which is 6 for the sync formal
                r = size(windowsync,2) - 3 #number of pools or sub-populations or populations in the sync file
                p_3d = Array{Float64}(undef, l, a, r)
                for pop in 1:r
                    for locus in 1:l
                        # #test:
                        # pop=1; locus=1
                        COUNTS = parse.(Int64, split(windowsync[locus, 3+pop], ":"))
                        FREQS = COUNTS ./ sum(COUNTS)
                        p_3d[locus, :, pop] = FREQS
                    end
                end
                # # remove fixed loci
                # idx_fixed_loci = Statistics.var(Statistics.maximum(p_3d, dims=2), dims=3) .> 1e-5
                # p_3d = p_3d[idx_fixed_loci[:],:,:]
                # # redefine the number of loci
                # l = size(p_3d)[1]
                # start computations
                n_bar = sum(pool_sizes) / r
                n_c = ( r*n_bar - ((sum(pool_sizes.^2))/(r*n_bar)) ) / (r-1)
                n_3d = reshape(repeat(pool_sizes, inner=l*a), (l,a,r))
                p_bar_2d = sum(n_3d .* p_3d, dims=3) ./ (r*n_bar)
                s_sq_2d = sum(n_3d .* ((p_3d .- p_bar_2d) .^ 2), dims=3) ./ (n_bar*(r-1))
                # building the HWE-based P(heterozygotes)
                h_3d = 2 .* p_3d .* (1 .- p_3d)
                h_bar_2d = sum(n_3d .* h_3d, dims=3) ./ (r * n_bar)
                # calculate the components of Fst
                a_2d = (n_bar/n_c) * ( (s_sq_2d) - ((1/(n_bar-1)) .* ( (p_bar_2d .* (1 .- p_bar_2d)) .- (((r-1)/r) .* s_sq_2d) .- (h_bar_2d ./ 4) )) )
                b_2d = (n_bar/(n_bar-1)) .* ( (p_bar_2d.* (1 .- p_bar_2d)) .- (((r-1)/r) .* s_sq_2d) .- ((((2*n_bar)-1)/(4*n_bar)) .* h_bar_2d) )
                c_2d = h_bar_2d ./ 2
                # calculate Fst and append to output vectors
                if sum(a_2d) <= 0.0 #for Fst values less than or equal to zero are essentially zero, i.e. Fst = (H_between - H_within) / H_between; --> H_between </= H_within
                    Fst = 0.0
                    append!(FST, Fst)
                else
                    Fst = sum(a_2d) / sum( a_2d .+ b_2d .+ c_2d )
                    append!(FST, Fst)
                end
                push!(CHROM, chrom)
                append!(WINDOW, window)
            end #non-zero window size
        end #window
    end #chrom
    OUT_Fst = try
                hcat((chr=convert(Array{String},CHROM), window=convert(Array{Int64},WINDOW), Fst=convert(Array{Float64},FST))...)
            catch
                hcat((chr=convert(Array{String},string.(CHROM)), window=convert(Array{Int64},WINDOW), Fst=convert(Array{Float64},FST))...)
            end
    return(OUT_Fst)
end

### [Hivert et al, 2018 method](https://www.biorxiv.org/content/biorxiv/early/2018/03/20/282400.full.pdf)
function Fst_hivert2018_func(sync::Array{Any,2}, pool_sizes::Array{Int64,1}, window_size::Int64)
    CHROM=[]; WINDOW=[]; FST=[]
    for chrom in unique(sync[:,1])
        # # test:
        # chrom = 1
        subsync = sync[sync[:,1] .== chrom, :]
        nWindows = convert(Int, ceil(maximum(subsync[:,2]) / window_size))
        for window in 1:nWindows
            # # test:
            # window = 2
            _start_ = window_size*(window-1) + 1
            _end_ = window_size*(window)
            windowsync = subsync[(subsync[:,2] .>= _start_) .& (subsync[:,2] .<= _end_), :]
            if (size(windowsync,1) != 0)
                l = size(windowsync,1)     #number of loci in this window
                a = 6                               #number of alleles which is 6 for the sync formal
                r = size(windowsync,2) - 3 #number of pools or sub-populations or populations in the sync file
                c_3d = Array{Int64}(undef, l, a, r)
                p_3d = Array{Float64}(undef, l, a, r)
                for pop in 1:r
                    for locus in 1:l
                        COUNTS = parse.(Int64, split(windowsync[locus, 3+pop], ":"))
                        FREQS = COUNTS ./ sum(COUNTS)
                        c_3d[locus, :, pop] = COUNTS
                        p_3d[locus, :, pop] = FREQS
                    end
                end
                C1 = sum(c_3d)
                C1i = reshape(sum(c_3d, dims=[1, 2]), r)
                D2 = sum( (C1i .+ (pool_sizes .- 1)) ./ pool_sizes )
                Pik = reshape(mean(p_3d, dims=1), a, r)'
                Pk = reshape(mean(p_3d, dims=[1,3]), a)'
                # ### filter-out alleles with zero frequencies
                # Pik = Pik[:, (Pk .!= 0.0)[:]]
                # Pk = Pk[:, (Pk .!= 0.0)[:]]
                D2a = sum(C1i .* (C1i .+ (pool_sizes .- 1)) ./ pool_sizes) / C1
                C2 = sum(C1i.^2)
                nc = (C1 - (C2/C1)) / (D2 - 1)
                MSP_k = (C1 - D2) .* sum( C1i .* ((Pik .- Pk).^2), dims=1 )
                MSI_k = (D2 - D2a) .* sum( C1i .* Pik .* (1.0 .- Pik), dims=1 )
                Fst = sum( MSP_k .- MSI_k ) / sum( MSP_k .+ ((nc-1) .* MSI_k) )
                if sum( MSP_k .+ ((nc-1) .* MSI_k) ) == 0.0 ### zero denominator indicating ZERO DIFFERENTIATON!
                    Fst = 0.0
                elseif (Fst < 0.0) ### negative Fst interpreted as zero!
                    Fst = 0.0
                end
                append!(FST, Fst)
                push!(CHROM, chrom)
                append!(WINDOW, window)
            end #non-zero window size
        end #window
    end #chrom
    OUT_Fst = try
                hcat((chr=convert(Array{String},CHROM), window=convert(Array{Int64},WINDOW), Fst=convert(Array{Float64},FST))...)
            catch
                hcat((chr=convert(Array{String},string.(CHROM)), window=convert(Array{Int64},WINDOW), Fst=convert(Array{Float64},FST))...)
            end
    return(OUT_Fst)
end

######################
### Main functions ###
######################
"""
# ___________________________________________________________________
# Fixation index (Fst) estimation from Pool-seq data across all pools

`Fst_across_pools(;filename_sync::String, window_size::Int64, pool_sizes::Array{Int64,1}, METHOD::String)`

Compute **across pools Fst** using [Weir and Cockerham, 1984 method](https://www.jstor.org/stable/2408641?seq=1)
or using [Hivert et al, 2018 method](https://www.biorxiv.org/content/biorxiv/early/2018/03/20/282400.full.pdf)

# Input
1. *filename_sync*: synchronized pileup filename (filtered by minimum allele frequency and depth)
2. *window_size*: size of non-overlapping sliding window in base-pairs
3. *pool-sizes*: array of pool sizes corresponding to each genotype column in the sync file
4. *METHOD*: Fst calculation method to use:
- "WeirCock"
- "Hivert"

# Output
1. Fst estimates per window with the filename: `string(join(split(filename_sync, ".")[1:(end-1)], '.'), "_window_", window_size, "bp_Fst_data.csv")`
2. Mean Fst estimate with the filename: `string(join(split(filename_sync, ".")[1:(end-1)], '.'), "_window_", window_size, "bp_Fst_sumstats.csv")`

# Examples
```
relatedness_module.Fst_across_pools(filename_sync="test/test.sync", window_size=100000, pool_sizes=[20,20,20,20,20], METHOD="WeirCock")
relatedness_module.Fst_across_pools(filename_sync="test/test.sync", window_size=100000, pool_sizes=[20,20,20,20,20], METHOD="Hivert")
```
"""
function Fst_across_pools(;filename_sync::String, window_size::Int64, pool_sizes::Array{Int64,1}, METHOD::String)
    println("#################################################################################################")
    println("Fst estimation across pools:")
    ### load the sync allele frequency data
    sync = DelimitedFiles.readdlm(filename_sync)
    ### input parameters
    nPools = size(sync,2)  - 3
    nLoci = size(sync,1)
    println(string("Synchronized pileup file: ", filename_sync))
    println(string("Non-overlapping window size: ", window_size, " bp"))
    println(string("Number of pools: ", nPools))
    for i in 1:length(pool_sizes)
        println(string("Pool ", i, " size : ", pool_sizes[i]))
    end
    println(string("Number of loci: ", nLoci))
    ### Across pools Fst estimation
    if METHOD == "WeirCock"
        OUT_Fst = Fst_weir1984_func(sync, pool_sizes, window_size)
    elseif METHOD == "Hivert"
        OUT_Fst = Fst_hivert2018_func(sync, pool_sizes, window_size)
    else
        println(string("AwWwwW! SowWwWyyyy! We have not implemented the Pool-seq Fst calculation method: ", METHOD, " yet."))
        println("¯\\_(๑❛ᴗ❛๑)_/¯ ʚ(´◡`)ɞ")
    end
    println("Writing the output into the files:")
    Fst_data_out_fname = string(join(split(filename_sync, ".")[1:(end-1)], '.'), "_window_", window_size, "bp_Fst_data.csv")
    Fst_sumstats_out_fname = string(join(split(filename_sync, ".")[1:(end-1)], '.'), "_window_", window_size, "bp_Fst_sumstats.csv")
    println(string("(1/2) ", Fst_data_out_fname, " for the Fst per chromosome per locus."))
    println(string("(2/2) ", Fst_sumstats_out_fname, " for the summary statistics of Fst across chromosomes and loci (i.e. mean, min, max and var)."))
    println("#################################################################################################")
    ### performing a Shapiro-Wilk's test for normality shows that the Fst across windows is normally distributed!
    SUMSTATS_Fst = (mean_Fst=Statistics.mean(OUT_Fst[.!isnan.(OUT_Fst[:,3]),3]), min_Fst=Statistics.minimum(OUT_Fst[.!isnan.(OUT_Fst[:,3]),3]), max_Fst=Statistics.maximum(OUT_Fst[.!isnan.(OUT_Fst[:,3]),3]), var_Fst=Statistics.var(OUT_Fst[.!isnan.(OUT_Fst[:,3]),3]))
    ### write-out
    writedlm(Fst_data_out_fname, reshape(["CHROM", "POS", "FST"],1,3), ',') ### write header
    io = open(Fst_data_out_fname, "a"); writedlm(io, OUT_Fst, ','); close(io) ### append the column-binded GWAlpha output
    writedlm(Fst_sumstats_out_fname, reshape(vcat(string.(keys(SUMSTATS_Fst))...), 1, length(SUMSTATS_Fst)), ',') ### write header
    io = open(Fst_sumstats_out_fname, "a"); writedlm(io, hcat(SUMSTATS_Fst...), ','); close(io) ### append the column-binded GWAlpha output
    return(0)
end

"""
# _______________________________________________________________
# Pairwise estimation of fixation indices (Fst) from Pool-seq data

`Fst_pairwise(;filename_sync::String, window_size::Int64, pool_sizes::Array{Int64,1}, METHOD::String)`

Compute **pairwise Fst** using [Weir and Cockerham, 1984 method](https://www.jstor.org/stable/2408641?seq=1)
or using [Hivert et al, 2018 method](https://www.biorxiv.org/content/biorxiv/early/2018/03/20/282400.full.pdf)

# Input
1. *filename_sync*: synchronized pileup filename (filtered by minimum allele frequency and depth)
2. *window_size*: size of non-overlapping sliding window in base-pairs
3. *pool-sizes*: array of pool sizes corresponding to each genotype column in the sync file
4. *METHOD*: Fst calculation method to use:
    - "WeirCock"
    - "Hivert"

# Output
Pairwise Fst estimates with the filename: `string(join(split(filename_sync, ".")[1:(end-1)], '.'), "_COVARIATE_FST.csv")`

# Examples
```
relatedness_module.Fst_pairwise(filename_sync="test/test.sync", window_size=100000, pool_sizes=[20,20,20,20,20], METHOD="WeirCock")
relatedness_module.Fst_pairwise(filename_sync="test/test.sync", window_size=100000, pool_sizes=[20,20,20,20,20], METHOD="Hivert")
### Parallel computation:
using Distributed
Distributed.addprocs(length(Sys.cpu_info())-1)
@everywhere include("GWAlpha.jl/src/relatedness_module_parallel.jl")
relatedness_module.Fst_pairwise(filename_sync="test/test.sync", window_size=100000, pool_sizes=[20,20,20,20,20], METHOD="WeirCock")
```
"""
function Fst_pairwise(;filename_sync::String, window_size::Int64, pool_sizes::Array{Int64,1}, METHOD::String=["WeirCock", "Hivert"][1])
    println("#################################################################################################")
    println("Pairwise Fst estimation:")
    ### load the sync allele frequency data
    sync = DelimitedFiles.readdlm(filename_sync)
    ### input parameters
    nPools = size(sync,2)  - 3
    nLoci = size(sync,1)
    println(string("Synchronized pileup file: ", filename_sync))
    println(string("Non-overlapping window size: ", window_size, " bp"))
    println(string("Number of pools: ", nPools))
    for i in 1:length(pool_sizes)
        println(string("Pool ", i, " size : ", pool_sizes[i]))
    end
    println(string("Number of loci: ", nLoci))
    ### Pairwise Fst estimation
    nPairs = convert(Int64, (nPools * (nPools - 1)) / 2)
    PAIRWISE_Fst = SharedArrays.SharedArray{Float64,2}(nPools, nPools)
    if length(Distributed.procs()) > 1
    	println("Performing Pairwise Fst estimation in parallel...")
     	@time x = @sync @distributed for i in 1:(nPools-1)
            for j in (i+1):nPools
                sub = sync[:, [1,2,3,(3+i), (3+j)]]
                if METHOD == "WeirCock"
                    Fst = Fst_weir1984_func(sub, pool_sizes[[i,j]], window_size)[:,3]
                    PAIRWISE_Fst[i, j] = mean(Fst[.!isnan.(Fst)])
                elseif METHOD == "Hivert"
                    Fst = Fst_hivert2018_func(sub, pool_sizes[[i,j]], window_size)[:,3]
                    PAIRWISE_Fst[i, j] = mean(Fst[.!isnan.(Fst)])
                else
                    println(string("AwWwwW! SowWwWyyyy! We have not implemented the Pool-seq Fst calculation method: ", model, " yet."))
                    println("¯\\_(๑❛ᴗ❛๑)_/¯ ʚ(´◡`)ɞ")
                end
                PAIRWISE_Fst[j, i] = PAIRWISE_Fst[i, j]
            end
        end
    else
        @time for i in 1:(nPools-1)
            for j in (i+1):nPools
                sub = sync[:, [1,2,3,(3+i), (3+j)]]
                if METHOD == "WeirCock"
                    Fst = Fst_weir1984_func(sub, pool_sizes[[i,j]], window_size)[:,3]
                    PAIRWISE_Fst[i, j] = mean(Fst[.!isnan.(Fst)])
                elseif METHOD == "Hivert"
                    Fst = Fst_hivert2018_func(sub, pool_sizes[[i,j]], window_size)[:,3]
                    PAIRWISE_Fst[i, j] = mean(Fst[.!isnan.(Fst)])
                else
                    println(string("AwWwwW! SowWwWyyyy! We have not implemented the Pool-seq Fst calculation method: ", model, " yet."))
                    println("¯\\_(๑❛ᴗ❛๑)_/¯ ʚ(´◡`)ɞ")
                end
                PAIRWISE_Fst[j, i] = PAIRWISE_Fst[i, j]
            end
        end
    end
    Fst_data_out_fname = string(join(split(filename_sync, ".")[1:(end-1)], '.'), "_COVARIATE_FST.csv")
    DelimitedFiles.writedlm(Fst_data_out_fname, PAIRWISE_Fst, ',')
    println("#################################################################################################")
    return(0)
end

"""
# _______________________________
# Standardized relatedness matrix

`standardized_relatedness(filename_sync::String)`

Compute a simple standardized relatedness matrix XX'/p, where X is the allele frequency matrix (Pool-seq data) and p is the number of columns of X

# Input
1. *filename_sync*: synchronized pileup filename (filtered by minimum allele frequency and depth)

# Output
Standardized relatedness matrix with the filename: `string(join(split(filename_sync, ".")[1:(end-1)], '.'), "_COVARIATE_RELATEDNESS.csv")`

# Examples
```
relatedness_module.standardized_relatedness("test/test.sync")
```
"""
function standardized_relatedness(filename_sync::String)
    println("#################################################################################################")
    println("Relatedness estimation:")
    ### check if the parsed sync file i.e. the allele frequency csv file exists
    ### if it does then use it, otherwise parse the sync file
    filename_parsed_sync = string(join(split(filename_sync, ".")[1:(end-1)], '.'), "_ALLELEFREQ.csv")
	if isfile(filename_parsed_sync)
		# rm(filename_out)
		println(string("Using the existing parsed sync file in csv format: ", filename_parsed_sync))
    else
        sync_processing_module.sync_parse(filename_sync)
	end
    X_raw = convert(Array{Float64,2},DelimitedFiles.readdlm(filename_parsed_sync, ',')[:,4:end])
    X = X_raw[sum(X_raw, dims=2)[:,1] .!= 0.0, :]'
    K = X*X' ./ size(X,2)
    relatedness_fname = string(join(split(filename_sync, ".")[1:(end-1)], '.'), "_COVARIATE_RELATEDNESS.csv")
    DelimitedFiles.writedlm(relatedness_fname, K, ',')
    println("#################################################################################################")
    return(0)
end

end ### end relatedness_module

max_stars_repo_name: jeffersonfparil/GWAlpha


text:
## IMPORTS ##
using Pkg
Pkg.activate(joinpath(@__DIR__,"../.."))
include(joinpath(@__DIR__,"model.jl"))
include(joinpath(@__DIR__,"utilities.jl"))

## PARAMETERS ##

# Model
N = 8
t = 1
U = 8
ω = 0.2*t
g0 = 0.2*t
g1 = 0.15*t
λ = 0.05*t
doping = 0

# Simulation
T = 40
τ = 0.05
DMRG_numsweeps = 20
DMRG_maxdim = 800
TEBD_maxdim = 800
TEBD_cutoff = 1E-10
DMRG_cutoff = 1E-10

# Number of phonons to try
max_phonons_list = [1, 2, 3, 4, 5, 7]
gs_energy_λ = []
gs_energy_no_λ = []
phonon_density_λ = []
phonon_density_no_λ = []
ΔEs = []

for max_phonons in max_phonons_list

    # Initialize
    println("Initializing with ", max_phonons, " phonons...")
    params = parameters(N=N, t=t, U=U, ω=ω, g0=g0, g1=g1, λ=λ, doping=doping,
                        max_phonons=max_phonons,
                        DMRG_numsweeps=DMRG_numsweeps,
                        DMRG_maxdim=DMRG_maxdim, DMRG_cutoff=DMRG_cutoff,
                        T=T, τ=τ, TEBD_cutoff=TEBD_cutoff)
    hubbholst = HubbardHolsteinModel(params)

    # Run DMRG
    println("Finding ground state for λ=",λ)
    dmrg_results = run_DMRG(hubbholst, params, alg="divide_and_conquer")

    # compensate for offset that comes from having a nonzero λ
    E = dmrg_results.ground_state_energy
    nb = compute_phonon_number(dmrg_results.ground_state)
    ΔE = λ*sum(nb.^2)

    push!(ΔEs, ΔE)
    push!(gs_energy_λ, E)
    push!(phonon_density_λ, nb)

    # Initialize
    println("Initializing with 0 λ")
    params = parameters(N=N, t=t, U=U, ω=ω, g0=g0, g1=g1, λ=0, doping=doping,
                        max_phonons=max_phonons,
                        DMRG_numsweeps=DMRG_numsweeps,
                        DMRG_maxdim=DMRG_maxdim, DMRG_cutoff=DMRG_cutoff,
                        T=T, τ=τ, TEBD_cutoff=TEBD_cutoff)
    hubbholst = HubbardHolsteinModel(params)

    # Run DMRG
    println("Finding ground state for λ=0")
    dmrg_results = run_DMRG(hubbholst, params, alg="divide_and_conquer")

    # compensate for offset that comes from having a nonzero λ
    E = dmrg_results.ground_state_energy
    nb = compute_phonon_number(dmrg_results.ground_state)
    ΔE = λ*sum(nb.^2)

    push!(ΔEs, ΔE)
    push!(gs_energy_no_λ, E)
    push!(phonon_density_no_λ, nb)
end


text:
# This file was generated, do not modify it. # hide
for name in names(datac)
    col = datac[:, name]
    ratio_missing = sum(ismissing.(col)) / nrows(datac) * 100
    println(rpad(name, 30), round(ratio_missing, sigdigits=3))
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/functions/nonlinear/functions.jl
# Defines variable info and function types.
#############################################################################

include(joinpath(@__DIR__, "nonlinear","auxiliary_variables.jl"))

"""
$(TYPEDEF)

An abstract super-type used for representing constraints built by EAGO's backend.
"""
abstract type AbstractEAGOConstraint end

"""

Computes the lower interval bound for `AbstractEAGOConstraint` representing an
inequality constraint.
"""
function lower_interval_bound end

"""

Computes a tuple representing the lower and upper interval bounds for a
`AbstractEAGOConstraint` representing an equality constraint.
"""
function interval_bound end

"""

Eliminate fixed variables by rearrangment or restructuring of `AbstractEAGOConstraint`.
"""
function eliminate_fixed_variables! end

###
###
### Affine function
###
###

"""
$(TYPEDEF)

Current only used for bound tightening. Stores a representation
of an affine inequality.
"""
mutable struct AffineFunctionIneq <: AbstractEAGOConstraint
    terms::Vector{Tuple{Float64,Int}}
    constant::Float64
    len::Int
end
const AFI = AffineFunctionIneq

AffineFunctionIneq() = AffineFunctionIneq(Tuple{Float64,Int}[], 0.0, 0)
function AffineFunctionIneq(f::SAF, s::LT)
    terms = map(x -> (x.coefficient, x.variable.value), f.terms)
    AffineFunctionIneq(terms, f.constant - s.upper, length(f.terms))
end
function AffineFunctionIneq(f::SAF, s::GT)
    terms = map(x -> (-x.coefficient, x.variable.value), f.terms)
    AffineFunctionIneq(terms, s.lower - f.constant, length(f.terms))
end
function AffineFunctionIneq(f::VI; is_max = false)
    if is_max
        return AffineFunctionIneq(Tuple{Float64,Int}[(-1.0, f.value)], 0.0, 1)
    end
    AffineFunctionIneq(Tuple{Float64,Int}[(1.0, f.value)], 0.0, 1)
end

"""
$(TYPEDEF)

Current only used for bound tightening. Stores a representation
of an affine equality.
"""
mutable struct AffineFunctionEq <: AbstractEAGOConstraint
    terms::Vector{Tuple{Float64,Int}}
    constant::Float64
    len::Int
end
const AFE = AffineFunctionEq

AffineFunctionEq() = AffineFunctionEq(Tuple{Float64,Int}[], 0.0, 0)
function AffineFunctionEq(func::SAF, set::ET)
    terms = map(x -> (x.coefficient, x.variable.value), func.terms)
    return AffineFunctionEq(terms, func.constant - set.value, length(func.terms))
end

###
### Parsing definitions
###

function eliminate_fixed_variables!(f::T, v::Vector{VariableInfo}) where T <: Union{AffineFunctionIneq,
                                                                                AffineFunctionEq}
    deleted_count = 0
    i = 1
    while i + deleted_count <= f.len
        coeff, indx = @inbounds f.terms[i]
        variable_info = @inbounds v[indx]
        if variable_info.is_fixed
            f.constant += coeff*variable_info.lower_bound
            deleteat!(f.terms, i)
            deleted_count += 1
        else
            i += 1
        end
    end
    f.len -= deleted_count
    return nothing
end

###
###
### Quadratic function
###
###

"""
$(TYPEDEF)

Stores a general quadratic inequality constraint with a buffer.
"""
mutable struct BufferedQuadraticIneq <: AbstractEAGOConstraint
    func::SQF
    buffer::Dict{Int, Float64}
    saf::SAF
    len::Int
end
const BQI = BufferedQuadraticIneq

"""
$(TYPEDEF)

Stores a general quadratic equality constraint with a buffer.
"""
mutable struct BufferedQuadraticEq <: AbstractEAGOConstraint
    func::SQF
    minus_func::SQF
    buffer::Dict{Int, Float64}
    saf::SAF
    len::Int
end
const BQE = BufferedQuadraticEq

#=
mutable struct BufferedConvexQuadratic <: AbstractEAGOConstraint
    func::SQF
    buffer::Dict{Int, Float64}
    saf::SAF
    len::Int
end
=#

function create_buffer_dict(func::SQF)

    buffer = Dict{Int, Float64}()

    for term in func.quadratic_terms
        buffer[term.variable_1.value] = 0.0
        buffer[term.variable_2.value] = 0.0
    end

    for term in func.affine_terms
        buffer[term.variable.value] = 0.0
    end

    return buffer
end

BufferedQuadraticIneq() = BufferedQuadraticIneq(SQF(SQT[], SAT[], 0.0), Dict{Int, Float64}(), SAF(SAT[], 0.0), 0)

function BufferedQuadraticIneq(func::SQF, set::LT)

    buffer = create_buffer_dict(func)
    saf = SAF([SAT(0.0, VI(k)) for k in keys(buffer)], 0.0)
    len = length(buffer)
    cfunc = copy(func)
    cfunc.constant -= set.upper

    return BufferedQuadraticIneq(cfunc, buffer, saf, len)
end

function BufferedQuadraticIneq(func::SQF, set::GT)

    buffer = create_buffer_dict(func)
    saf = SAF([SAT(0.0, VI(k)) for k in keys(buffer)], 0.0)
    len = length(buffer)
    cfunc = MOIU.operate(-, Float64, func)
    cfunc.constant += set.lower

    return BufferedQuadraticIneq(cfunc, buffer, saf, len)
end

BufferedQuadraticEq() = BufferedQuadraticEq(SQF(SQT[], SAT[], 0.0), SQF(SQT[], SAT[], 0.0), Dict{Int, Float64}(), SAF(SAT[], 0.0), 0)

function BufferedQuadraticEq(func::SQF, set::ET)

    buffer = create_buffer_dict(func)
    saf = SAF([SAT(0.0, VI(k)) for k in keys(buffer)], 0.0)
    len = length(buffer)
    cfunc1 = copy(func)
    cfunc1.constant -= set.value
    cfunc2 = MOIU.operate(-, Float64, func)
    cfunc2.constant += set.value

    return BufferedQuadraticEq(cfunc1, cfunc2, buffer, saf, len)
end

#=
function BufferedConvexQuadratic(f::BufferedQuadraticIneq)
    BufferedConvexQuadratic(copy(f.func), copy(f.buffer), copy(f.saf), f.len)
end
=#

function eliminate_fixed_variables!(f::T, v::Vector{VariableInfo}) where T <: Union{BufferedQuadraticIneq,
                                                                                BufferedQuadraticIneq}
    deleted_count = 0
    i = 1
    while i + deleted_count <= f.len
        term = f.sqf.terms[i]
        variable_info_1 = v[term.variable_1.value]
        variable_info_2 = v[term.variable_2.value]
        if variable_info_1.is_fixed && variable_info_2.is_fixed
            f.sqf.constant += coeff*variable_info_1.lower_bound*variable_info_2.lower_bound
            deleteat!(f.sqf.terms, i)
            deleted_count += 1
        else
            i += 1
        end
    end
    f.len -= deleted_count

    return nothing
end

###
###
### Cone
###
###

"""
$(TYPEDEF)

Stores a second-order cone with a buffer.
"""
mutable struct BufferedSOC <: AbstractEAGOConstraint
    variables::VECOFVAR
    buffer::Dict{Int, Float64}
    saf::SAF
    len::Int
end

function BufferedSOC(func::VECOFVAR, set::SOC)
    len = length(func.variables)
    buffer = Dict{Int, Float64}([(variable.value, 0.0) for variable in func.variables])
    saf = SAF(fill(SAT(0.0, VI(1)), len), 0.0)
    return BufferedSOC(copy(func), buffer, saf, len)
end

include("nonlinear/nonlinear.jl")

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
# types.jl

# Physical Constants

"Planck's constant, (kgm²s⁻¹)."
const hbar = const ħ = 1.054571817e-34
"Electron charge, (kgm²s⁻²)."
const eV = const q = const ElectronVolt = 1.602176634e-19
"Electron mass, (kg)."
const me = MassElectron = 9.1093837015e-31
"Boltzmann's constant, (kgm²K⁻¹)."
const Boltzmann = const kB = 1.380649e-23
"Permittivity of free space, (C²N⁻¹m⁻²)."
const ε_0 = 8.85418682e-12
"Speed of light, (ms⁻¹)."
const c = 299792458

# Structures

"Polaron. Structure to store data of polaron solution and other parameters, for each temperature or frequency."
struct Polaron
    "T, temperature (K)."
    T
    "Kμ, Kadanoff mobility (cm²V⁻¹s⁻¹)."
    Kμ
    "Hμ, Hellwarth mobility (cm²V⁻¹s⁻¹)."
    Hμ
    "FHIPμ, FHIP mobility (cm²V⁻¹s⁻¹)."
    FHIPμ
    "k, spring constant."
    k
    "M, renormalised (phonon-drag) mass (mₑ)."
    M
    "Osaka free energy components (A,B,C) and total (F) (unitless). See Hellwarth et al. 1999 PRB Part IV."
    A; B; C; F
    "Tau, relaxation time from Kadanoff Boltzmann transport equation (s)."
    Tau
    "v and w, raw variational parameters (unitless)."
    v; w
    "βred, reduced thermodynamic beta (unitless)."
    βred
    "rfsi, Feynman polaron radius (Schultz) (m)."
    rfsi
    "rfsmallalpha, small-alpha asymptotic approximation (unitless)."
    rfsmallalpha

    # Setup of simulation. These parameters are sent to the function.

    "α, Fröhlich alpha coupling parameter (unitless)."
    α
    "mb, Band effective mass (mₑ)."
    mb
    "ω, effective dielectric frequency (2π THz)."
    ω
end

# structure initialisation
Polaron()=Polaron([],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[],[])

"NewPolaron. Structure to store data of polaron solution and other parameters, for each temperature or frequency."
struct NewPolaron
    "α, Fröhlich alpha coupling parameter (unitless)"
    α
    "T, temperature (K)."
    T
    "β, reduced Thermodynamic beta (unitless)."
    β
    "ω, phonon frequency (2π THz)."
    ω
    "v, variational parameters (unitless)."
    v
    "w, variational parameters (unitless)."
    w
    "κ, fictitious spring constant (multiples of mₑ) (kgs⁻²)."
    κ
    "M, fictitious particle (multiples of mₑ) (kg)."
    M
    "F, free energy (meV)."
    F
    "Ω, electric field frequencies (multiples of phonon frequency ω) (2π THz)."
    Ω
    "Z, complex impedence (V/A)."
    Z
    "σ, complex conductivity (A/V)."
    σ
    "μ, FHIP mobility (cm²V⁻¹s⁻¹)."
    μ
end

# Broadcast Polaron data.
function Base.show(io::IO, x::NewPolaron)
    flush(stdout)
    print(io, "-------------------------------------------------\n Polaron Information: \n-------------------------------------------------\n", "Fröhlich coupling | α = ", round.(x.α, digits = 3), " | sum(α) = ", round(sum(x.α), digits = 3),"\nTemperatures | T = ", round.(x.T, digits = 3), " K \nReduced thermodynamic | β = ", round.(x.β, digits = 3), "\nPhonon frequencies | ω = ", round.(x.ω, digits = 3), " 2π THz\nVariational parameters | v = ", round.(x.v, digits = 3), " ω | w = ", round.(x.w, digits = 3), " ω\nFictitious spring constant | κ = ", round.(x.κ, digits = 3), " kg/s²\nFictitious mass | M = ", round.(x.M, digits = 3), " kg\nFree energy | F = ", round.(x.F, digits = 3), " meV\nElectric field frequency | Ω = ", round.(Float64.(x.Ω), digits = 3),  " 2π THz\nComplex impedance | Z = ", x.Z .|> y -> round.(ComplexF64.(y), digits = 3), " V/A\nComplex conductivity | σ = ", x.σ .|> y -> round.(ComplexF64.(y), digits = 3), " A/V\nMobility | μ = ", round.(Float64.(x.μ), digits = 3), " cm²/Vs")
end


text:
using JuMP, Gurobi, Combinatorics, Cbc
import XLSX

#### read data ####
xf = XLSX.readxlsx("../data/data1.xlsx")
l = xf["General!B1"]     # work day duration
N = xf["General!B3"] + 1  # number of clients + 1
M = xf["General!B2"]  # number of providers
# + 1 is for including the headquarter (HS14)

# Clients
s = xf["Clients!A2:A$N"]   # client service start times
d = vcat([0], xf["Clients!B2:B$N"])   # durations for clients, set d[1] = 0 for HS14

# Providers
w = xf["Providers!B2:B$(M+1)"]  # start time of providers
f = xf["Providers!A2:A$(M+1)"]  # hiring cost of providers

# make separate s vector for each provider, inserting into the first element the provider's work start time
s_test = ones(M, N)
for k in 1:M
    global s_test[k,:] = vcat([w[k]], s)
end
s = s_test

# Distances
t = xf["Distances"][:]  # 2-d array of travel times
#### ######

L = 10000 # large number

E = collect(powerset(collect(2:N), 2, N-1)) # for subtour 'E'limination constraints

########

model = Model();
set_optimizer(model, Gurobi.Optimizer);
#(13)
@variable(model, x[1:N, 1:N, 1:M], Bin)

######################## routing constraints ######################

# (2) each client served once by one provider
@constraint(model, served[j in 2:N], sum(x[i,j,k] for i in 1:N for k in 1:M) == 1)

# (3) route continuity
@constraint(model, continuity[p in 2:N, k in 1:M], sum(x[i,p,k] for i in 1:N) == sum(x[p,j,k] for j in 1:N) )

# (4) not all providers need to leave 14HS (and start working)
@constraint(model, leave_depot[k in 1:M], sum(x[1,j,k] for j in 2:N) <= 1)

# (5) not all providers need to return to 14HS (i.e. they didn't leave 14HS in the first place)
@constraint(model, return_depot[k in 1:M], sum(x[i,1,k] for i in 2:N) <= 1)

# (6) no self looping
@constraint(model, diag[k in 1:M], sum(x[i,i,k] for i in 1:N) == 0)

# (7) at least 1 provider must leave 14HS and start working to serve clients
@constraint(model, sum(x[1,j,k] for j in 2:N for k in 1:M) >= 1)

# (8) all providers must start their work at 14HS, if they do work
@constraint(model, large[k in 1:M], sum(x[i,j,k] for i in 2:N for j in 2:N) <= L*sum(x[1,j,k] for j in 2:N))

#(9) subtour elimination for each provider
for k in 1:M
    for e in 1:length(E)
        @constraint(model, sum(x[i,j,k] for i in E[e] for j in E[e]) <= length(E[e])-1)
    end
end

######################## end routing constraints ######################

######################## time constraints ######################

# (10) initial client service times
@constraint(model, start[j in 2:N, k in 1:M], x[1,j,k]*(s[k,1] + t[1,j] + d[1]) <= s[k,j])

# (11) in between start and end of work, service times
@constraint(model, bet[i in 2:N, j in 2:N, k in 1:M], x[i,j,k]*(s[k,i] + t[i,j] + d[i]) <= s[k,j])

# (12) come back to 14HS by hour l.
@constraint(model, endofwork[i in 2:N, k in 1:M], x[i,1,k]*(s[k,i] + t[i,1] + d[i]) <= l)

######################## end time constraints ######################

# (1) minimise hiring cost to serve all clients
@objective(model, Min, sum(f[k]*x[1,j,k] for j in 1:N for k in 1:M))

optimize!(model)

(value.(x[:,:,1]))

for k in 1:M
    if sum(value.(x[1,:,k])) == 1
        println("vehicle $k provides service.")
        # come back hr
        for i in 2:N
            if value.(x[i,1,k]) == 1
                println("vehicle $k comes back to 14HS at $(t[i,1] + s[k,i] + d[i])hr. \n")
            end
        end
    end
end

path = [1 2 1]
#path = [1 7 1]
#path = [1 6 1]
#path = [1 3 4 5 1]

worked_hr = 0
#delivered  = 0
for i in 2:length(path)
    global worked_hr += t[path[i-1], path[i]] + d[path[i]]
    #global delivered += d[path[i]]
end
worked_hr


text:
try
    using StaticStoragesTests
    true
catch
    false
end || begin
    push!(LOAD_PATH, joinpath(@__DIR__, "StaticStoragesTests"))
    using StaticStoragesTests
end


text:
# Script to run through the files and fit the zero-inflated negative binomial distrubution

using Plots, Distributions, DelimitedFiles, Base.Printf, Base.Filesystem
include("Utils.jl")
import Main.Utils
gr()

# Set some MCMC parameters
Lchain = 4000
burn = 100
thin = 10
folder = "./"
# restart = "OxyR"

Files = readdir(folder)
filter!(f->occursin(".csv",f), Files)

for (ii,file) in enumerate(Files)
    if occursin(".csv",file)

        name = replace(file, ".csv"=>"")
		println("Running MCMC for "*name)

    	# Load and plot the distribution data
        data = Utils.load_data(file,folder, 151)
        x,y = Utils.genpdf(Integer.(round.(data)))

        # Specify likelihood function and initial guess
        lFunc = p->Utils.log_likelihood_zi_negbinom(data, p)
        guess = [1.0,0.5,0.5]

        # Run the MCMC
        priors = [truncated(Normal(0,20),0,Inf), Uniform(0.0,1.0),Uniform(0.0,1.0)]
		chain = Utils.mcmc_metropolis(guess, lFunc, Lchain; prior=priors, propStd=0.03,scaleProp=false, burn=burn,step=thin);
		chain_red = Utils.rmv_zeros(chain)

        # Chain for K/ν
        tmp = (1.0.-chain_red[:,2])./chain_red[:,2]
        chain_red = [chain_red tmp]

		# Extract parameters and plot
        r = Utils.find_MAP(chain_red,idx=1)
        p = Utils.find_MAP(chain_red,idx=2)
        w = Utils.find_MAP(chain_red,idx=3)
        K = Utils.find_MAP(chain_red,idx=4)

        # Plot and save
        plt = plot(plot(title=name), plot(legend=false), size=(600,800),layout=(2,1))
        bar!(plt[1], x,y, line=0,ylabel="Probability, P(n)", label="Data")
        plot!(plt[2], x,cumsum(y), xlabel="Copy number, n",ylabel="Cumulative probability, C(n)", line=:steppost)
        m = MixtureModel([DiscreteUniform(0,0), NegativeBinomial(r,p)], [w,1-w])
        plot!(plt[1], x, x->pdf(m,x), label="Model")
        plot!(plt[2], 0:maximum(x),n->cdf(m,n))

        # Save data
		try
        	writedlm("Chains/"*name,chain_red)
		catch
			mkdir("Chains")
			writedlm("Chains/"*name,chain_red)
		end
		try
        	Plots.pdf("Distributions/"*name)
		catch
			mkdir("Distributions")
	        Plots.pdf("Distributions/"*name)
		end
    end
end

max_stars_repo_name: rdbrackston/tx-analysis

max_issues_repo_name: rdbrackston/tx-analysis

max_forks_repo_name: rdbrackston/tx-analysis


text:
Base.@deprecate_binding StaticBasket Basket
Base.@deprecate_binding DynamicBasket Basket

@deprecate simplefv(pv, rate, periods) pv * (one(rate) + rate * periods)
@deprecate compoundfv(pv, rate, periods) pv * (one(rate) + rate)^periods


text:
import HTTP
import JSON

include("./request.jl")

function delete_project(server::Bool=true) :: Nothing
    cd("..")
    rm("demo", recursive=true)

    if server
        Dance.CoreEngine.close_server()
    end

    Dance.Router.delete_routes!()
    nothing
end

function dirs_add_multiple() :: Nothing
    mkdir("code")
    mkdir("node_modules")  # test ignores node dir

    cd("code")  # test dir with file & sub-dir
    touch("file1.jl")
    mkdir("sub-dir1")  # test dir with only sub-dir
    cd("sub-dir1")
    mkdir("sub-dir2")  # test dir with file
    cd("sub-dir2")
    touch("file2.jl")

    cd("../../..")
    nothing
end

function dirs_add_single() :: Nothing
    mkdir("code")
    mkdir("node_modules")  # test ignores node dir

    cd("code")  # test dir with file & sub-dir
    touch("file1.jl")

    cd("..")
    nothing
end

function make_and_test_request_get(path::String, status::Int64, headers::Dict{String, String}, content_length::Int64, is_json_body::Bool, body::Any) :: Nothing
    r = HTTP.request("GET", "http://127.0.0.1:8000$path")
    parse_and_test_request(r, status, headers, content_length, is_json_body, body)
end

function make_and_test_request_options(path::String) :: Nothing
    r = HTTP.request("OPTIONS", "http://127.0.0.1:8000$path")
    compare_http_header(r.headers, "Allow", "POST")
    compare_http_header(r.headers, "Access-Control-Allow-Methods", "POST")
    compare_http_header(r.headers, "Access-Control-Allow-Headers", "X-PINGOTHER, Content-Type")
end

function make_and_test_request_post(path::String, payload::Union{Array, Dict}, status::Int64, headers::Dict{String, String}, content_length::Int64, is_json_body::Bool, body::Any) :: Nothing
    r = HTTP.request("POST", "http://127.0.0.1:8000$path", [], JSON.json(payload))
    parse_and_test_request(r, status, headers, content_length, is_json_body, body)
end

function project_settings() :: Nothing
    ## Add `dev.jl` file with 1 overwrite & 1 new entry ##
    cd("demo/settings")

    touch("dev.jl")
    open("dev.jl", "w") do io
        write(io, ":dev = true\n")
        write(io, ":foo = \"bar\"")
    end

    open("Global.jl", "a+") do io
        write(io, "include(\"dev.jl\")")
    end

    cd("..")
end

function routes(file_suffix::String) :: Nothing
    open("routes.jl", "w") do io_routes
        open("../sample/routes/" * file_suffix * ".jl") do io_file
            write(io_routes, io_file)
        end
    end

    nothing
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""NodeStatus is information about the current status of a node.

    IoK8sApiCoreV1NodeStatus(;
        addresses=nothing,
        allocatable=nothing,
        capacity=nothing,
        conditions=nothing,
        config=nothing,
        daemonEndpoints=nothing,
        images=nothing,
        nodeInfo=nothing,
        phase=nothing,
        volumesAttached=nothing,
        volumesInUse=nothing,
    )

    - addresses::Vector{IoK8sApiCoreV1NodeAddress} : List of addresses reachable to the node. Queried from cloud provider, if available. More info: https://kubernetes.io/docs/concepts/nodes/node/#addresses Note: This field is declared as mergeable, but the merge key is not sufficiently unique, which can cause data corruption when it is merged. Callers should instead use a full-replacement patch. See http://pr.k8s.io/79391 for an example.
    - allocatable::Dict{String, IoK8sApimachineryPkgApiResourceQuantity} : Allocatable represents the resources of a node that are available for scheduling. Defaults to Capacity.
    - capacity::Dict{String, IoK8sApimachineryPkgApiResourceQuantity} : Capacity represents the total resources of a node. More info: https://kubernetes.io/docs/concepts/storage/persistent-volumes#capacity
    - conditions::Vector{IoK8sApiCoreV1NodeCondition} : Conditions is an array of current observed node conditions. More info: https://kubernetes.io/docs/concepts/nodes/node/#condition
    - config::IoK8sApiCoreV1NodeConfigStatus : Status of the config assigned to the node via the dynamic Kubelet config feature.
    - daemonEndpoints::IoK8sApiCoreV1NodeDaemonEndpoints : Endpoints of daemons running on the Node.
    - images::Vector{IoK8sApiCoreV1ContainerImage} : List of container images on this node
    - nodeInfo::IoK8sApiCoreV1NodeSystemInfo : Set of ids/uuids to uniquely identify the node. More info: https://kubernetes.io/docs/concepts/nodes/node/#info
    - phase::String : NodePhase is the recently observed lifecycle phase of the node. More info: https://kubernetes.io/docs/concepts/nodes/node/#phase The field is never populated, and now is deprecated.
    - volumesAttached::Vector{IoK8sApiCoreV1AttachedVolume} : List of volumes that are attached to the node.
    - volumesInUse::Vector{String} : List of attachable volumes in use (mounted) by the node.
"""
mutable struct IoK8sApiCoreV1NodeStatus <: SwaggerModel
    addresses::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1NodeAddress} } # spec name: addresses
    allocatable::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: allocatable
    capacity::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: capacity
    conditions::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1NodeCondition} } # spec name: conditions
    config::Any # spec type: Union{ Nothing, IoK8sApiCoreV1NodeConfigStatus } # spec name: config
    daemonEndpoints::Any # spec type: Union{ Nothing, IoK8sApiCoreV1NodeDaemonEndpoints } # spec name: daemonEndpoints
    images::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1ContainerImage} } # spec name: images
    nodeInfo::Any # spec type: Union{ Nothing, IoK8sApiCoreV1NodeSystemInfo } # spec name: nodeInfo
    phase::Any # spec type: Union{ Nothing, String } # spec name: phase
    volumesAttached::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1AttachedVolume} } # spec name: volumesAttached
    volumesInUse::Any # spec type: Union{ Nothing, Vector{String} } # spec name: volumesInUse

    function IoK8sApiCoreV1NodeStatus(;addresses=nothing, allocatable=nothing, capacity=nothing, conditions=nothing, config=nothing, daemonEndpoints=nothing, images=nothing, nodeInfo=nothing, phase=nothing, volumesAttached=nothing, volumesInUse=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("addresses"), addresses)
        setfield!(o, Symbol("addresses"), addresses)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("allocatable"), allocatable)
        setfield!(o, Symbol("allocatable"), allocatable)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("capacity"), capacity)
        setfield!(o, Symbol("capacity"), capacity)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("conditions"), conditions)
        setfield!(o, Symbol("conditions"), conditions)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("config"), config)
        setfield!(o, Symbol("config"), config)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("daemonEndpoints"), daemonEndpoints)
        setfield!(o, Symbol("daemonEndpoints"), daemonEndpoints)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("images"), images)
        setfield!(o, Symbol("images"), images)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("nodeInfo"), nodeInfo)
        setfield!(o, Symbol("nodeInfo"), nodeInfo)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("phase"), phase)
        setfield!(o, Symbol("phase"), phase)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("volumesAttached"), volumesAttached)
        setfield!(o, Symbol("volumesAttached"), volumesAttached)
        validate_property(IoK8sApiCoreV1NodeStatus, Symbol("volumesInUse"), volumesInUse)
        setfield!(o, Symbol("volumesInUse"), volumesInUse)
        o
    end
end # type IoK8sApiCoreV1NodeStatus

const _property_map_IoK8sApiCoreV1NodeStatus = Dict{Symbol,Symbol}(Symbol("addresses")=>Symbol("addresses"), Symbol("allocatable")=>Symbol("allocatable"), Symbol("capacity")=>Symbol("capacity"), Symbol("conditions")=>Symbol("conditions"), Symbol("config")=>Symbol("config"), Symbol("daemonEndpoints")=>Symbol("daemonEndpoints"), Symbol("images")=>Symbol("images"), Symbol("nodeInfo")=>Symbol("nodeInfo"), Symbol("phase")=>Symbol("phase"), Symbol("volumesAttached")=>Symbol("volumesAttached"), Symbol("volumesInUse")=>Symbol("volumesInUse"))
const _property_types_IoK8sApiCoreV1NodeStatus = Dict{Symbol,String}(Symbol("addresses")=>"Vector{IoK8sApiCoreV1NodeAddress}", Symbol("allocatable")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}", Symbol("capacity")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}", Symbol("conditions")=>"Vector{IoK8sApiCoreV1NodeCondition}", Symbol("config")=>"IoK8sApiCoreV1NodeConfigStatus", Symbol("daemonEndpoints")=>"IoK8sApiCoreV1NodeDaemonEndpoints", Symbol("images")=>"Vector{IoK8sApiCoreV1ContainerImage}", Symbol("nodeInfo")=>"IoK8sApiCoreV1NodeSystemInfo", Symbol("phase")=>"String", Symbol("volumesAttached")=>"Vector{IoK8sApiCoreV1AttachedVolume}", Symbol("volumesInUse")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiCoreV1NodeStatus }) = collect(keys(_property_map_IoK8sApiCoreV1NodeStatus))
Swagger.property_type(::Type{ IoK8sApiCoreV1NodeStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1NodeStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1NodeStatus }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1NodeStatus[property_name]

function check_required(o::IoK8sApiCoreV1NodeStatus)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1NodeStatus }, name::Symbol, val)
end


text:
module LinearFilter

export guassian_kernel, apply_linear_filter, apply_roi_aggregate
# Kernel function
"""
`gaussian_kernel(train::Array{Float64,1}, sigma::Float64)`

The Gaussian Kernel.

```math
w(\\tau) = \\frac{1}{\\sqrt{2 \\pi} \\sigma_w}
\\exp(-\\frac{(\\mathbf{t}-\\tau)^2}{2 * \\sigma ^2})
```
"""
function gaussian_kernel(train::Array{Float64,1}, sigma::Float64=1.0)
    return t->sum(1/sqrt(2*pi)/sigma .* exp.(-broadcast(-, train, t).^2./(2*sigma^2)), 1)
end

function apply_linear_filter(k, tstart, tend; step=100)
    _t = tstart:1/step:tend |> Array{Float64, 1} |> transpose
    k(_t)
end

"""
```apply_roi_aggregate(kernel::LinearFilter.kernel, starts::Array{Float64,1}, roi::Range, backgrounds::Range)```

Returns:
- _result :: Array{Float64, 2}(signal, trial)
"""
function apply_roi_aggregate(kernel, starts, roi, backgrounds)
    _result = zeros(length(starts),length(roi))
    for idx = 1:length(starts)
        _mark = starts[idx]
        _x = (_mark + roi) |> Array{Float64, 1} |> transpose
        _backgrounds = (_mark + backgrounds) |> Array{Float64, 1} |> transpose
        _target = kernel(_x)
        _backgrounds = kernel(_backgrounds) |> mean

        _result[idx,:] = _target - _backgrounds
    end
    _result'
end

end


text:
"""
A `Solution` comprises a single solution vector and its parent
"""
type Solution
	data::Vector{Float64}
	parent # ::Branch # circular type dependency...
	Solution(d::Vector{Float64}) = new(d,nothing)
	Solution(d::Vector{Float64}, p) = new(d,p)
	Solution(s::Solution, p) = new(deepcopy(s.data), p)
end
@relay(Solution, :data)

Base.show(io::Base.IO, P::Solution) = write(io, "***SOLUTION***")


text:
module Scenes

using CthulhuVision.Light
using CthulhuVision.Math
using CthulhuVision.Triangles
using CthulhuVision.Worlds
using CthulhuVision.Materials
using CthulhuVision.BVH

export Scene, SceneSettings, Mesh
export addmesh!, addinstance!

struct SceneSettings
    ambientemission::RGB
end

struct MeshLocation
    first::UInt32
    last::UInt32
end

struct Scene
    vertexes::Vector{Vector3}
    meshtriangles::Vector{MeshTriangle}
    meshes::Vector{MeshLocation}
    instances::Vector{MeshInstance}
    hitables::Vector{Hitable}
    settings::SceneSettings

    Scene(settings::SceneSettings) = new([], [], [], [], [], settings)
end

struct Mesh
    vertexes::Vector{Vector3}
    triangles::Vector{MeshTriangle}
end

function offsetmeshtriangle(mt::MeshTriangle, vertexoffset::UInt32) :: MeshTriangle
    MeshTriangle(
        mt.vertex1 + vertexoffset,
        mt.vertex2 + vertexoffset,
        mt.vertex3 + vertexoffset
    )
end

function addmesh!(scene::Scene, mesh::Mesh) :: UInt32
    # MeshTriangle indexes are one-indexed when used as relative indexes inside a mesh.
    # Therefore, we offset all indexes by the existing length of the vertexes vector.
    vertexoffset = UInt32(length(scene.vertexes))
    append!(scene.vertexes, mesh.vertexes)

    meshtriangles = [offsetmeshtriangle(mt, vertexoffset) for mt in mesh.triangles]
    firstmeshtriangle = length(scene.meshtriangles) + 1
    lastmeshtriangle = firstmeshtriangle + length(meshtriangles) - 1
    append!(scene.meshtriangles, meshtriangles)

    location = MeshLocation(firstmeshtriangle, lastmeshtriangle)
    push!(scene.meshes, location)
    meshindex = length(scene.meshes)

    meshindex
end

function addinstance!(scene::Scene, instance::MeshInstance)
    push!(scene.instances, instance)
    instanceindex = length(scene.instances)

    # This tells us where all the triangle indexes are.
    locations = scene.meshes[instance.meshindex]

    for index = locations.first:locations.last
        meshtriangle = scene.meshtriangles[index]
        v1 = instance.tform * scene.vertexes[meshtriangle.vertex1]
        v2 = instance.tform * scene.vertexes[meshtriangle.vertex2]
        v3 = instance.tform * scene.vertexes[meshtriangle.vertex3]

        box = AABB(v1, v2, v3)

        hitable = Hitable(box, index, instanceindex)
        push!(scene.hitables, hitable)
    end
end

end


text:
module Scallops

using RayTraceEllipsoids, HCubature, Statistics, Optim, StaticArrays, DataFrames, CSV, Interpolations
export createobj, raytrace_sensitivity, main_figure
RETINAS = (:dis_ret, :pro_ret)
mutable struct Membrane
    passthrough::Float64
    positive::Bool
    rxy::Float64
    rz::Float64
    cz::Float64
    θ0::Float64
    Membrane(passthrough::Float64, positive::Bool, rxy::Real, rz::Real) = new(passthrough, positive, rxy, rz, 0.0, 0.0)
end
Membrane(p::Float64, d::Bool, r::Real) = Membrane(p, d, r, r)
Membrane() = Membrane(1.0, false, 0)
mutable struct Tissue
    ri::Float64
    thick::Float64
    absorption_coefficient::Float64
end
Tissue() = Tissue(0.0, 0.0, 0.0)
struct Layer
    dis_membrane::Membrane
    tissue::Tissue
    pro_membrane::Membrane
end
Layer() = Layer(Membrane(), Tissue(), Membrane())
function getlayers(absorption_coefficient, water_ac, reflectance)
    water_ri = 1.334
    sky              = Membrane(1., true, Inf, Inf)
    water            = Tissue(water_ri, 0, water_ac)
    dis_cornea       = Membrane(1., true, 244, 227)
    cornea           = Tissue(1.37, 23, water_ac)
    dis_lens         = Membrane(1., true, 225, 213)
    lens             = Tissue(1.42, 215, water_ac)
    dis_dis_ret      = Membrane(1., false, 337 - 6)
    dis_ret          = Tissue(1.35, 12, absorption_coefficient[:dis_ret])
    pro_dis_ret      = Membrane(1., false, 337 + 6)
    dis_ret2prox_ret = Tissue(1.34, 102-12-30, water_ac)
    dis_pro_ret      = Membrane(1., false, 337 - 15)
    pro_ret          = Tissue(1.35, 30, absorption_coefficient[:pro_ret])
    pro_pro_ret      = Membrane(1., false, 337 + 15)
    pro_ret2mirror   = Tissue(1.34, 136-102, water_ac)
    mirror           = Membrane(reflectance, false, 417)
    behind_mirror    = Tissue(0, Inf, Inf)
    earth            = Membrane(1., false, Inf)
    (water = Layer(sky, water, dis_cornea),
     cornea = Layer(dis_cornea, cornea, dis_lens),
     lens = Layer(dis_lens, lens, dis_dis_ret),
     dis_ret = Layer(dis_dis_ret, dis_ret, pro_dis_ret),
     dis_ret2prox_ret = Layer(pro_dis_ret, dis_ret2prox_ret, dis_pro_ret),
     pro_ret = Layer(dis_pro_ret, pro_ret, pro_pro_ret),
     pro_ret2mirror = Layer(pro_pro_ret, pro_ret2mirror, mirror),
     mirror = Layer(mirror, behind_mirror, earth))
end
function norefraction!(layers)
    ks = collect(keys(layers))
    todo = setdiff(ks, [:water, :mirror])
    ri = layers.water.tissue.ri
    for i in todo
        layers[i].tissue.ri = ri
    end
    layers
end
function morphlayer!(l, morphz, morphxy)
    l.dis_membrane.rxy *= morphxy
    l.dis_membrane.rz *= morphz
    l.tissue.thick *= morphz
end
function morph2layers(morphz, absorption_coefficient, water_ac, reflectance)
    morphxy = sqrt(1/morphz)
    layers = getlayers(absorption_coefficient, water_ac, reflectance)
    for l in layers
        morphlayer!(l, morphz, morphxy)
    end
    z = sum(l.tissue.thick for l in layers if !isinf(l.tissue.thick))/2
    for l in layers
        l.dis_membrane.cz =  z - signit(l.dis_membrane.positive)*l.dis_membrane.rz
        z -= l.tissue.thick
    end
    layers
end
signit(positive::Bool) = -(-1)^positive
x2h(x, a, b, positive) = signit(positive)*b/a*sqrt(a^2 - x^2)
getzs(l::Layer, photoreceptor_radius) = map([:pro_membrane, :dis_membrane]) do membrane
    m = getfield(l, membrane)
    h = x2h(photoreceptor_radius, m.rxy, m.rz, m.positive)
    m.cz + h
end
function _getmedium(l, photoreceptor_radius)
    z1, z2 = getzs(l, photoreceptor_radius)
    Retina(Cylinder(photoreceptor_radius, z1, z2), l.tissue.absorption_coefficient)
end
getmedium(::Val{:dis_ret}, l, photoreceptor_radius) = _getmedium(l, photoreceptor_radius)
getmedium(::Val{:pro_ret}, l, photoreceptor_radius) = _getmedium(l, photoreceptor_radius)
getmedium(_, l, _) = NonRetina(l.tissue.absorption_coefficient)
function getopticunits(layers, photoreceptor_radius, ignore_incoming)
    nlayers = length(layers)
    mediums1 = [getmedium(Val(k), v, photoreceptor_radius) for (k,v) in pairs(layers)]
    iₙ = findfirst(isequal(:dis_ret), keys(layers))
    medium_inds = [1:nlayers - 1; nlayers - 1:-1:iₙ]
    ns_ind1 = [1:nlayers - 1; nlayers - 1:-1:iₙ]
    ns_ind2 = [2:nlayers; nlayers - 2:-1:iₙ - 1]
    pointin_inds = [2:nlayers; nlayers-1:-1:iₙ]
    funs = [fill(!, nlayers-1); fill(identity, nlayers - iₙ)]
    surface_inds = [2:nlayers; nlayers - 1:-1:iₙ]
    mediums = mediums1[medium_inds]
    if ignore_incoming
        for retina in RETINAS
            i = findfirst(isequal(retina), keys(layers))
            mediums[i] = NonRetina(layers[1].tissue.absorption_coefficient)
        end
    end
    nss = [layers[i1].tissue.ri/layers[i2].tissue.ri for (i1, i2) in zip(ns_ind1, ns_ind2)]
    pointins = [fun(layers[i].dis_membrane.positive) for (i, fun) in zip(pointin_inds, funs)]
    ellipsoids = [cspheroid(layers[i].dis_membrane.cz, layers[i].dis_membrane.rxy, layers[i].dis_membrane.rz, layers[i].dis_membrane.positive) for i in surface_inds]
    ks = keys(layers)
    names = [i < nlayers ? Symbol(ks[j], :_in) : ks[j] for (i,j) in enumerate(medium_inds)]
    passes = [layers[i].dis_membrane.passthrough for i in surface_inds]
    NamedTuple{Tuple(names)}(OpticUnit.(ellipsoids, pointins, nss, mediums, passes))
end
angulardiff(r::TraceRay, θ) = abs(acos(r.dir[3]) - abs(θ))
angulardiff(r::DeadRay, θ) = abs(θ)
function findnodalz(lb, ub, ous, rc)
    θs = range(-0.02, stop=0.1, length = 4)
    function fun(nodalz, θ)
        l = Light(10^4, rc.rxy, rc, θ, nodalz)
        r = l(TraceRay, 0.5)
        r = raytrace!(ous, r)
        angulardiff(r, θ)
    end
    nodals = map(θs) do θ
        res = optimize(nodalz -> fun(nodalz, θ), lb, ub)
        @assert Optim.minimum(res) < 1e-4 "failed to find optimum for the nodal point"
        # Optim.minimum(res) > 1e-4 && println("failed nodalz")
        Optim.minimizer(res)
    end
    itp = interpolate(nodals, BSpline(Cubic(Line(OnGrid()))))
    sitp = scale(itp, θs)
    sitp(0)
end
function createobj(morphz, photoreceptor_radius, absorption_coefficient; ignore_incoming=false, water_ac = 0.03e-6, reflectance = 0.9)
    layers = morph2layers(morphz, absorption_coefficient, water_ac, reflectance)
    # norefraction!(layers)
    mem = layers[1].pro_membrane
    rc = (rxy = mem.rxy, rz = mem.rz, cz = mem.cz)
    ous = getopticunits(layers, photoreceptor_radius, ignore_incoming)
    nodalz = findnodalz(layers.mirror.dis_membrane.cz - layers.mirror.dis_membrane.rz, layers.mirror.dis_membrane.cz + layers.mirror.dis_membrane.rz, ous, rc)
    fls = [abs(mean((layers[retina].dis_membrane.cz - layers[retina].dis_membrane.rz, layers[retina].pro_membrane.cz - layers[retina].pro_membrane.rz)) - nodalz) for retina in RETINAS]
    focallengths = NamedTuple{RETINAS}(fls)
    (layers = layers, rc = rc, ous = ous, nodalz = nodalz, focallengths = focallengths)
end
function raytrace_fraction(distance, aperture, rc, ous, nodalz)
    l = Light(distance, aperture, rc, 0, nodalz)
    function fun(b)
        for retina in RETINAS
            ous[retina].medium.signal.photoreceptor = 0.0
            ous[retina].medium.signal.retina = 0.0
        end
        r = l(Ray, b[1], b[2])
        raytrace!(ous, r)
        SVector(ous[RETINAS[1]].medium.signal.photoreceptor, ous[RETINAS[1]].medium.signal.retina, ous[RETINAS[2]].medium.signal.photoreceptor, ous[RETINAS[2]].medium.signal.retina)
    end
    s, _ = hcubature(fun, SVector(0., 0.), SVector(1., 0.5), initdiv=20)#, maxevals=10^4)
    NamedTuple{RETINAS}(((photoreceptor = 2s[1], retina = 2s[2]), (photoreceptor = 2s[3], retina = 2s[4])))
end
function height2sigma(signal)
    h = signal.photoreceptor/signal.retina
    1/(sqrt(2π)*h)
end
sigma2fwhm(σ) = 2sqrt(2log(2))*σ
spatial2angular(fwhm, focallength) = 2atan(fwhm/2/focallength)
sigma2angular(σ, focallength) = spatial2angular(sigma2fwhm(σ), focallength)
height2angular(signal, focallength) = sigma2angular(height2sigma(signal), focallength)
function fwhm_height(distance, aperture, rc, ous, nodalz, focallengths)
    ss = raytrace_fraction(distance, aperture, rc, ous, nodalz)
    NamedTuple{RETINAS}(([height2angular(ss[retina], focallengths[retina]) for retina in RETINAS]))
end
function optimize_morth4fwhm(retina, distance, aperture, photoreceptor_radius, absorption_coefficient)
    function fun(morphz)
        layers, rc, ous, nodalz, focallengths = createobj(morphz, photoreceptor_radius, absorption_coefficient)
        fwhm = fwhm_height(distance, aperture, rc, ous, nodalz, focallengths)
        fwhm[retina]
    end
    res = optimize(fun, 1, 1.5)#, Brent(), abs_tol = 1e-2)
    Optim.minimizer(res)#, Optim.minimum(res)
end
function optimize_morth(retina, distance, aperture, photoreceptor_radius, absorption_coefficient)
    function fun(morphz)
        layers, rc, ous, nodalz, focallengths = createobj(morphz, photoreceptor_radius, absorption_coefficient)
        Δθ = getΔΘ(focallengths[retina], photoreceptor_radius)
        cosa = cos(Δθ/2)
        s = raytrace_fraction(distance, aperture, rc, ous, nodalz, cosa, retina)
        1/s
    end
    res = optimize(fun, 1, 1.5, Brent(), abs_tol = 1e-2)
    Optim.minimizer(res)#, Optim.minimum(res)
end
getΔΘ(focallength, photoreceptor_radius) = 2atan(photoreceptor_radius/focallength)
modified_photons(D, Δθ) = (π/4)^2*D^2*Δθ^2
modified_land(D, Δθ, k, l) = modified_photons(D, Δθ)*RayTraceEllipsoids.absorption(l, k)
function main_figure(; min_aperture = 200, max_aperture = 400, n_apertures = 8, distance = 10^5, photoreceptor_radius = 5, absorption_coefficient = (dis_ret = 0.00667, pro_ret = 0.00667), filename = joinpath(@__DIR__(), "data"))
    apertures = range(min_aperture, stop = max_aperture, length = n_apertures)
    morphz = 1.1
    layers, rc, ous, nodalz, focallengths = createobj(morphz, photoreceptor_radius, absorption_coefficient)
    df = DataFrame(retina = Symbol[], morpah_factor = Float64[], aperture = Float64[], focallength = Float64[], Δθ = Float64[], dis_ret_ac = Float64[], pro_ret_ac = Float64[], outer_segment_length = Float64[], theoretical_sensitivity = Float64[], raytrace_sensitivity = Float64[], fwhm = Float64[])
    allowmissing!(df)
    for (i, aperture) in enumerate(apertures)
        fwhm = fwhm_height(distance, aperture, rc, ous, nodalz, focallengths)
        for (j, retina) in enumerate(RETINAS)
            Δθ = getΔΘ(focallengths[retina], photoreceptor_radius)
            l = layers[retina].tissue.thick
            theo_sen = modified_land(aperture, δθ, absorption_coefficient[retina], l)
            raytrace_sen = raytrace_sensitivity(distance, aperture, rc, nodalz, ous, retina)
            push!(df, (retina, morphz, aperture, focallengths[retina], Δθ, absorption_coefficient[:dis_ret], absorption_coefficient[:pro_ret], l, theo_sen, raytrace_sen, rad2deg(fwhm[retina])))
        end
    end
    CSV.write("$(filename)dis.csv", df[df[:retina] .== :dis_ret,2:end])
    CSV.write("$(filename)pro.csv", df[df[:retina] .== :pro_ret,2:end])
    CSV.write("$filename.csv", df)
end
function burn!(x)
    for i in length(x)-1:-1:1
        if x[i] < x[i+1]
            x[i] = x[i+1]
        end
    end
end
function getΔΘ(distance, aperture, rc, nodalz, ous, retina)
    function θ2signal(θ)
        l = Light(distance, aperture, rc, θ, nodalz)
        ous[retina].medium.signal.photoreceptor = 0.0
        count = raytrace!(ous, l)
        ous[retina].medium.signal.photoreceptor/count
    end
    maxθ = retina == :dis_ret ? 0.2 : 0.6
    x = range(0, maxθ, length=25)
    y = θ2signal.(x)
    burn!(y)
    itp = interpolate(y, BSpline(Quadratic(Reflect(OnCell()))))
    sitp = scale(itp, x)
    x2 = range(0, maxθ, length = 1000)
    hm = sitp(0)/2
    for x in x2
        if sitp(x) < hm
            return 2x
        end
    end
end
function raytrace_fraction(distance, aperture, rc, ous, nodalz, cosa, retina)
    count = 0
    ous[retina].medium.signal.photoreceptor = 0.0
    n = 1000
    for b in range(0, 1, length=n)
        θ = π - acos(b*(cosa - 1) - cosa)
        l = Light(distance, aperture, rc, θ, nodalz)
        count += raytrace!(ous, l, n = n)
    end
    ous[retina].medium.signal.photoreceptor/count
end
function unmodified_photons(distance, cosa, aperture, rc, nodalz)
    Sₑ = 2π*distance^2*(1 - cosa)
    Sₐ = π*(aperture/2)^2
    h = rc.rz*sqrt(1 - (aperture/2/rc.rxy)^2) + rc.cz
    D = distance + nodalz - h
    Sₑ*Sₐ/D^2
end
function raytrace_sensitivity(distance, aperture, rc, nodalz, ous, retina)
    Δθ = getΔΘ(distance, aperture, rc, nodalz, ous, retina)
    cosa = cos(Δθ/2)
    s = raytrace_fraction(distance, aperture, rc, ous, nodalz, cosa, retina)
    Fₐ = unmodified_photons(distance, cosa, aperture, rc, nodalz)
    s*Fₐ
end

end


text:
# Testing continuous univariate distributions

using Distributions, Random
using Test

using Calculus: derivative

n_tsamples = 100

# additional distributions that have no direct counterparts in R references
@testset "Testing $(distr)" for distr in [Biweight(),
                                          Biweight(1,3),
                                          Epanechnikov(),
                                          Epanechnikov(1,3),
                                          Triweight(),
                                          Triweight(2),
                                          Triweight(1, 3),
                                          Triweight(1)]

    test_distr(distr, n_tsamples; testquan=false)
end

# Test for non-Float64 input
using ForwardDiff
@test string(logpdf(Normal(0,1),big(1))) == "-1.418938533204672741780329736405617639861397473637783412817151540482765695927251"
@test derivative(t -> logpdf(Normal(1.0, 0.15), t), 2.5) ≈ -66.66666666666667
@test derivative(t -> pdf(Normal(t, 1.0), 0.0), 0.0) == 0.0

# Test for numerical problems
@test pdf(Logistic(6,0.01),-2) == 0

@testset "Normal with std=0" begin
    d = Normal(0.5,0.0)
    @test pdf(d, 0.49) == 0.0
    @test pdf(d, 0.5) == Inf
    @test pdf(d, 0.51) == 0.0

    @test cdf(d, 0.49) == 0.0
    @test cdf(d, 0.5) == 1.0
    @test cdf(d, 0.51) == 1.0

    @test ccdf(d, 0.49) == 1.0
    @test ccdf(d, 0.5) == 0.0
    @test ccdf(d, 0.51) == 0.0

    @test quantile(d, 0.0) == -Inf
    @test quantile(d, 0.49) == 0.5
    @test quantile(d, 0.5) == 0.5
    @test quantile(d, 0.51) == 0.5
    @test quantile(d, 1.0) == +Inf

    @test rand(d) == 0.5
    @test rand(MersenneTwister(123), d) == 0.5
end

# Test for parameters beyond those supported in R references
@testset "VonMises with large kappa" begin
    d = VonMises(1.1, 1000)
    @test var(d) ≈ 0.0005001251251957198
    @test entropy(d) ≈ -2.034688918525470
    @test cf(d, 2.5) ≈ -0.921417 + 0.38047im atol=1e-6
    @test pdf(d, 0.5) ≈ 1.758235814051e-75
    @test logpdf(d, 0.5) ≈ -172.1295710466005
    @test cdf(d, 1.0) ≈ 0.000787319 atol=1e-9
end