text:
#####
##### Game
#####

export Game, play!, tournament!

mutable struct Game
    table::Table
end

function Base.show(io::IO, game::Game)
    println(io, "\n----------------------- Poker game")
    show(io, game.table, false)
    for player in players_at_table(game)
        show(io, player, false)
    end
    println(io, "-----------------------")
end

function Game(
        players::Tuple;
        deck = ordered_deck(),
        table = nothing,
        dealer_id::Int = default_dealer_id(),
        blinds::Blinds = Blinds(1,2),
    )

    n_player_cards = sum(map(x->cards(x)==nothing ? 0 : length(cards(x)), players))

    @assert 2 ≤ length(players) ≤ 10 "Invalid number of players"

    if length(deck) ≠ 52
        # if the deck isn't full, then players should have been dealt cards.
        @assert n_player_cards > 0

        # If the user is specifying the player cards, then the
        # user should probably also handle the table cards too.
        @assert table ≠ nothing
        @assert cards(table) isa Tuple{Card,Card,Card,Card,Card}
        @assert length(cards(table)) == 5
        @assert length(deck) + n_player_cards + length(cards(table)) == 52
    else # nobody has been dealt yet
        table = Table(;
            deck=deck,
            players=players,
            dealer_id=dealer_id,
            blinds=blinds
        )

        @assert length(deck) == 52
        @assert n_player_cards == 0
        @assert cards(table) == nothing
    end

    return Game(table)
end

players_at_table(game::Game) = players_at_table(game.table)
blinds(game::Game) = blinds(game.table)
any_actions_required(game::Game) = any_actions_required(game.table)
state(game::Game) = state(game.table)
move_buttons!(game) = move_buttons!(game.table)

print_game_state(table, state::PreFlop) = @info "Pre-flop!"
print_game_state(table, state::Flop) =  @info "Flop: $(repeat(" ", 44)) $(table.cards[1:3])"
print_game_state(table, state::Turn) =  @info "Turn: $(repeat(" ", 44)) $(table.cards[4])"
print_game_state(table, state::River) = @info "River: $(repeat(" ", 43)) $(table.cards[5])"

set_preflop_blind_raise!(table::Table, player, ::AbstractGameState, i::Int) = nothing
function set_preflop_blind_raise!(table::Table, player::Player, ::PreFlop, i::Int)
    if 1 ≤ i ≤ length(players_at_table(table))
        if is_first_to_act(table, player)
            # everyone must call big blind to see flop:
            table.current_raise_amt = blinds(table).big
            table.initial_round_raise_amt = blinds(table).big
        end
    end
end
reset_round_bank_rolls!(game::Game, state::PreFlop) = nothing # called separately prior to deal
reset_round_bank_rolls!(game::Game, state::AbstractGameState) = reset_round_bank_rolls!(game.table)

# TODO: compactify. Some of these cases/conditions may be redundant
function end_of_actions(table::Table, player)
    players = players_at_table(table)
    @debug "     last_to_raise.(players) = $(last_to_raise.(players))"
    @debug "     checked.(players) = $(checked.(players))"
    @debug "     folded.(players) = $(folded.(players))"
    @debug "     action_required.(players) = $(action_required.(players))"
    @debug "     !any(action_required.(players)) = $(!any(action_required.(players)))"
    @debug "     all_playing_all_in(table) = $(all_playing_all_in(table))"
    @debug "     all_playing_checked(table) = $(all_playing_checked(table))"
    @debug "     all_all_in_or_checked(table) = $(all_all_in_or_checked(table))"
    @debug "     all_all_in_except_bank_roll_leader(table) = $(all_all_in_except_bank_roll_leader(table))"
    @debug "     all_oppononents_all_in(table, player) = $(all_oppononents_all_in(table, player))"

    case_1 = last_to_raise(player)
    case_2 = all_playing_checked(table)
    case_3 = all_playing_all_in(table)
    case_4 = all_all_in_except_bank_roll_leader(table)
    case_5 = all_all_in_or_checked(table) # TODO: likely replaceable with case_6
    case_6 = !any(action_required.(players))
    case_7 = all_oppononents_all_in(table, player) && paid_current_raise_amount(table, player)
    @debug "     cases = $((case_1, case_2, case_3, case_4, case_5, case_6, case_7))"

    return any((case_1, case_2, case_3, case_4, case_5, case_6, case_7))
end

function last_player_to_raise(table::Table)
    for player in players_at_table(table)
        last_to_raise(player) && return player
    end
    return nothing
end

function all_raises_were_called(table::Table)
    lptr = last_player_to_raise(table)
    lptr===nothing && return true

    players = players_at_table(table)
    ltr = last_to_raise.(players)
    @assert count(ltr) == 1
    players_who_called = []
    @debug "Checking if all raises were called"
    @debug "     table.winners.declared = $(table.winners.declared)"
    arwc = false
    conds = map(players) do player
        cond1 = seat_number(player) == seat_number(lptr)
        cond2 = not_playing(player)
        cond3 = all_in(player)
        cond4 = pot_investment(player) ≈ pot_investment(lptr)
        (cond1, cond2, cond3, cond4)
    end
    @debug begin
        conds_debug = map(players) do player
            sn = seat_number(player)
            cond1 = seat_number(player) == seat_number(lptr)
            cond2 = not_playing(player)
            cond3 = all_in(player)
            cond4 = pot_investment(player) ≈ pot_investment(lptr)
            (sn, cond1, cond2, cond3, cond4)
        end
        for cond in conds_debug
            @debug "sn, cond = $(cond[1]), $(cond[2:end])"
        end
        @debug "snlptr = $(seat_number(lptr))"
    end
    return all(map(cond->any(cond), conds))
end

end_preflop_actions(table, player, ::AbstractGameState) = false
function end_preflop_actions(table::Table, player::Player, ::PreFlop)
    cond1 = is_big_blind(table, player)
    cond2 = checked(player)
    cond3 = still_playing(player)
    return all((cond1, cond2, cond3))
end

function act_generic!(game::Game, state::AbstractGameState)
    table = game.table
    table.winners.declared && return
    set_state!(table, state)
    print_game_state(table, state)
    reset_round_bank_rolls!(game, state)

    any_actions_required(game) || return
    play_out_game(table) && return
    set_play_out_game!(table)
    for (i, player) in enumerate(circle(table, FirstToAct()))
        @debug "Checking to see if it's $(name(player))'s turn to act"
        @debug "     not_playing(player) = $(not_playing(player))"
        @debug "     all_in(player) = $(all_in(player))"
        not_playing(player) && continue # skip players not playing
        set_preflop_blind_raise!(table, player, state, i)
        if end_of_actions(table, player)
            break
        end
        all_in(player) && continue
        @debug "$(name(player))'s turn to act"
        player_option!(game, player)
        table.winners.declared && break
        end_preflop_actions(table, player, state) && break
        if i > n_max_actions(table)
            error("Too many actions have occured, please open an issue.")
        end
    end
    @info "Betting is finished."
    @assert all_raises_were_called(table)
end

function act!(game::Game, state::AbstractGameState)
    act_generic!(game, state)
    reset_round!(game.table)
end

"""
    play!(::Game)

Play a game. Note that this method
expects no cards (players and table)
to be dealt.
"""
function play!(game::Game)
    @info "------ Playing game!"

    table = game.table
    set_active_status!(table)
    players = players_at_table(table)

    initial_brs = deepcopy(bank_roll.(players))
    initial_∑brs = sum(initial_brs)

    @info "Initial bank roll summary: $(bank_roll.(players))"
    did = dealer_id(table)
    sb = seat_number(small_blind(table))
    bb = seat_number(big_blind(table))
    f2a = seat_number(first_to_act(table))
    @info "Buttons (dealer, small, big, 1ˢᵗToAct): ($did, $sb, $bb, $f2a)"

    @assert still_playing(dealer(table)) "The button must be placed on a non-folded player"
    @assert still_playing(small_blind(table)) "The small blind button must be placed on a non-folded player"
    @assert still_playing(big_blind(table)) "The big blind button must be placed on a non-folded player"
    @assert still_playing(first_to_act(table)) "The first-to-act button must be placed on a non-folded player"

    @assert dealer_id(table) ≠ small_blind_id(table) "The button and small blind cannot coincide"
    @assert small_blind_id(table) ≠ big_blind_id(table) "The small and big blinds cannot coincide"
    @assert big_blind_id(table) ≠ first_to_act_id(table) "The big blind and first to act cannot coincide"

    table.transactions = TransactionManager(players)

    @assert all(cards.(players) .== nothing)
    @assert cards(table) == nothing
    reset_round_bank_rolls!(table) # round bank-rolls must account for blinds
    deal!(table, blinds(table))
    @assert cards(table) ≠ nothing

    winners = table.winners

    winners.declared || act!(game, PreFlop())   # Pre-flop bet/check/raise
    winners.declared || act!(game, Flop())      # Deal flop , then bet/check/raise
    winners.declared || act!(game, Turn())      # Deal turn , then bet/check/raise
    winners.declared || act!(game, River())     # Deal river, then bet/check/raise

    winners.players = distribute_winnings!(players, table.transactions, cards(table))
    winners.declared = true

    @debug "amount.(table.transactions.side_pots) = $(amount.(table.transactions.side_pots))"
    @debug "initial_∑brs = $(initial_∑brs)"
    @debug "sum(bank_roll.(players_at_table(table))) = $(sum(bank_roll.(players_at_table(table))))"
    @debug "initial_brs = $(initial_brs)"
    @debug "bank_roll.(players_at_table(table)) = $(bank_roll.(players_at_table(table)))"

    @assert initial_∑brs ≈ sum(bank_roll.(players_at_table(table))) # eventual assertion
    @assert sum(amount.(table.transactions.side_pots)) ≈ 0

    @info "Final bank roll summary: $(bank_roll.(players))"
    @assert winners.declared

    @info "------ Finished game!"
    return winners
end

function set_active_status!(table::Table)
    players = players_at_table(table)
    for player in players
        if zero_bank_roll(player)
            player.active = false
            player.folded = true
            player.action_required = false
        else
            player.active = true
            player.folded = false
            player.action_required = true
        end
    end
end

function reset_game!(game::Game)
    table = game.table
    players = players_at_table(table)

    game.table = Table(;
        deck=ordered_deck(),
        players=players,
        dealer_id=dealer_id(table),
        blinds=table.blinds
    )
    table = game.table
    players = players_at_table(table)
    for player in players
        player.cards = nothing
        player.pot_investment = 0
        player.all_in = false
        player.round_bank_roll = bank_roll(player)
        player.checked = false
        player.last_to_raise = false
        player.round_contribution = 0
        player.active = true
    end
    set_active_status!(table)
    table.play_out_game = false
end

"""
    tournament!(game::Game)

Play until a single player remains!
"""
function tournament!(game::Game)
    @info "********************************"
    @info "******************************** Playing tournament!"
    @info "********************************"
    table = game.table
    players = players_at_table(table)
    while length(players) > 1
        play!(game)
        n_players_remaining = count(map(x->!(bank_roll(x) ≈ 0), players))
        if n_players_remaining ≤ 1
            println("Victor emerges!")
            break
        end
        reset_game!(game)
        move_buttons!(game)
    end
    @info "********************************"
    @info "******************************** Finished tournament!"
    @info "********************************"
    return game.table.winners
end


text:
using CausticsEngineering
using Test

include("triangles.jl")

ext: jl


text:
using Random # , Distributions # not sure if we need that
using Markdown
using LaTeXStrings

export hint, tip, protip, almost, warning_box, danger
export correct, still_missing, still_nothing, wrong_type
export var_not_defined, func_not_defined, keep_working
export not_defined  # deprecated
#export keep_working_if_var_contains_substr
export type_isa, type_eq
export code_for_check_type_funcs
export TODO
export display_msg_if_fail

"Hint box with arguement as text."
hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]));

"Tip box with arguement as text."
tip(text) = Markdown.MD(Markdown.Admonition("tip", "Tip", [text]));

"Tip box with arguement as text."
protip(text, invite = "Curious to learn more?";  boxlabel = "Pro Tip") = Foldable(invite, Markdown.MD(Markdown.Admonition("tip", boxlabel, [text])) );

"Admonition box labeled a warning with arguement as text."
almost(text) = Markdown.MD(Markdown.Admonition("warning", "Almost there!", [text]));

"warning box with arguement as text."
warning_box(text) = Markdown.MD(Markdown.Admonition("warning", "Warning:", [text]));

"Danger box with arguement as text."
danger(text) = Markdown.MD(Markdown.Admonition("danger", "Beware!", [text]));

"Admonition box with reminder to replace missing."
still_missing(text=md"Replace `missing` with your answer.") = Markdown.MD(Markdown.Admonition("warning", "Missing Response", [text]));

"Admonition box with reminder to replace nothing."
still_nothing(text=md"Replace `nothing` with your answer.") = Markdown.MD(Markdown.Admonition("warning", "Here we go!", [text]));

wrong_type() = Markdown.MD(Markdown.Admonition("danger", "Type Error", [md"Check the type of your response."]))

wrong_type(var::Symbol, type::Type; text=md"The type of $var should be $type") = Markdown.MD(Markdown.Admonition("danger", "Type Error", [text]))

"Admonition box with reminder that function name passed is not defined."
func_not_defined(func_name) = Markdown.MD(Markdown.Admonition("danger", "Oopsie!", [md"Make sure that you define a function called **$(Markdown.Code(string(func_name)))**"]));

"Admonition box with reminder that variable name passed is not defined."
var_not_defined(variable_name) = Markdown.MD(Markdown.Admonition("danger", "Oopsie!", [md"Make sure that you define a variable called **$(Markdown.Code(string(variable_name)))**"]));

"Admonition box with reminder that variable name passed is not defined.  (deprecated)"
not_defined(variable_name) = var_not_defined(variable_name)

"Admonition box warning that the answer isn't quite right."
keep_working(text=md"The answer is not quite right.") = Markdown.MD(Markdown.Admonition("danger", "Keep working on it!", [text]));

function keep_working_if_var_contains_substr(var::Symbol,str::Union{String,Markdown.MD},substr::String)
# I had to remove !@isdefined(var) due to how Pluto puts variables into different modules
# not exported, so provide function with same name in notebook
  if     ismissing(var)
         still_missing()
  elseif isnothing(var)
         still_nothing()
  else
      if typeof(str) == Markdown.MD
         str = Markdown.plain(str)
      end
      if occursin(substr,str)
         keep_working(md"Make sure to update the cell setting $var.")
      end
  end
end

type_isa(var, t::Union{Type,Vector{Type},Vector{DataType}}) = (any(typeof(var) .<: t))
type_eq(var, t::Union{Type,Vector{Type},Vector{DataType}}) = (any(typeof(var) .== t))

# To use functions in notebook, run eval(Meta.parse(code_for_check_type_funcs))
code_for_check_type_funcs = """
begin
function check_type_isa(sym::Symbol, var, t::Union{Type,Vector{Type},Vector{DataType}})
   if ismissing(var)
        passed = false
        text = md"The variable \$sym is still set to missing."
        msg = still_missing(text)
   elseif !type_isa(var,t) # (any(typeof(var) .<: t))
        passed = false
        text = "The type of \$sym is not correct.  It should be <: "
        if typeof(t) <: Type
           text = text * "\$t."
        else
           text = text * "one of "
           for tt in t
              if tt == last(t)
                 text = text * " or \$t."
              else
                 text = text * " \$t,"
              end
           end
        end
        msg = Markdown.MD(Markdown.Admonition("danger", "Type Error", [Markdown.MD(text)]))
   else
        passed = true
        msg = Markdown.MD("\$sym has the correct type.")
   end
   return (;passed, msg)
end
function check_type_eq(sym::Symbol, var, t::Union{Type,Vector{Type},Vector{DataType}})
   if ismissing(var)
        passed = false
        text = md"The variable \$sym is still set to missing."
        msg = still_missing(text)
   elseif !type_eq(var, t) # (any(typeof(var) == t))
        passed = false
        if typeof(t) <: Type
           text = "The type of \$sym should be \$t."
        else
           text = "The type of \$sym should be one of "
           for tt in t
              if tt == last(t)
                 text = text * " or \$t."
              else
                 text = text * " \$t,"
              end
           end
        end
        #text = md"The type of \$sym should be \$t."
        msg = Markdown.MD(Markdown.Admonition("danger", "Type Error", [Markdown.MD(text)]))
   else
        passed = true
        msg = Markdown.MD("\$sym has the correct type.")
   end
   return (;passed, msg)
end
end
"""

yays = [md"Great!", md"Yay ❤", md"Great! 🎉", md"Well done!", md"Keep it up!", md"Good job!", md"Awesome!", md"You got the right answer!", md"Let's move on to the next part."];

"Box with random positive message."
correct(text=rand(yays)) = Markdown.MD(Markdown.Admonition("correct", "Got it!", [text]));

TODO_str = html"TODO"
TODO() = TODO_str

function display_msg_if_fail(x; msg_pass = nothing)
   if !x.passed
      x.msg
   else
      msg_pass
   end
end


text:
using EVQueues, Plots

#tiempos de llegada. Debe estar ordenado.
arribos = [1.0;2.0;3.0];
#tiempos de partidas, correlativo al de arribos.
salidas = [12.0;13.0;14.0];
#tiempos de trabajo, correlativo al de arribos.
trabajos = [4.0;5.0;6.0];
#potencias
potencias = [1.0;1.0;1.0];
#cantidad de cargadores.
C=1.0;

#simula usando edf a partir de la traza. Cambiar edf por llf, llr, pf, parallel para las otras politicas.
sim = ev_edf_trace(arribos,trabajos,salidas,potencias,C,snapshots=[4.0])
compute_statistics!(sim)

#TIme plot of occupation and power
p=plot(sim.timetrace)
display(p)

#CDF of departure attained workloads
p=servicecdf(sim.EVs)
display(p)

#State space of the last snapshot
snap = sim.snapshots[end];

w = [ev.currentWorkload for ev in snap.charging];
d = [ev.currentDeadline for ev in snap.charging];
on = [ev.currentPower>0 for ev in snap.charging];

p = plot(   xlabel = "Remaining workload",
            ylabel = "Remaining soj. time",
            title = "State-space snapshot",
            xlims = (0,3*sim.parameters["AvgEnergy"]),
            ylims = (0,3*sim.parameters["AvgDeadline"])
            )

scatter!(p,w[on.==true],d[on.==true],markershape=:circle,markersize=4,color=:blue,label="In service")
scatter!(p,w[on.==false],d[on.==false],markershape=:circle,markersize=4,color=:red,label="Not in service")

display(p)


text:
module Foo
"""
```jldoctest
julia> reshape(Int64[1, 2, 3, 4] * 2, (4,1,1))
4×1×1 Array{Int64,3}:
[:, :, 1] =
 1
 2
 3
 4

julia> reshape(Int64[1, 2, 3, 4], (4,1,1))
4-element Array{Int64,1}:
 1
 2
 3
 4
```
```jldoctest
julia> reshape(Int64[1, 2, 3, 4], (4,1,1))
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> reshape(Int64[1, 2, 3, 4] * 2, (4,1,1))
4-element Array{Int64,1}:
 1
 2
 3
 4
```
```jldoctest
julia> reshape(Int64[1, 2, 3, 4] * 2, (4,1,1))
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> reshape(Int64[1, 2, 3, 4] * 2, (4,1,1))
4-element Array{Int64,1}:
 1
 2
 3
 4
```
```jldoctest
julia> begin
           reshape(Int64[1, 2, 3, 4] * 2, (4,1,1))
       end
4-element Array{Int64,1}:
 1
 2
 3
 4
```
```jldoctest
reshape(Int64[1, 2, 3, 4] * 2, (4,1,1))

# output

4-element Array{Int64,1}:
 1
 2
 3
 4
```
```jldoctest; filter = r"foo"
julia> println("  foobar")
  foobaz
```
```jldoctest
julia> 1 + 2

julia> 3 + 4
```
"""
foo() = 1

    """
    ```jldoctest
    julia> begin
               reshape(Int64[1, 2, 3, 4] * 2, (4,1,1))
           end
    4-element Array{Int64,1}:
     1
     2
     3
     4
    ```
    ```jldoctest
    julia> println(); println("foo")

    bar
    ```
    """
    foo(x) = 1

end # module


text:
"""
    gelombang(f,g,a,b,c,m,n)
berisi program untuk mencari solusi persamaan differensial parsial yaitu persamaan gelombang.
Program ini secara
default berisi 7 masukan, yaitu fungsi nilai batas `f(x)` dan `g(x)`, batas solusi `a` dan `b`,
koefisien persamaan gelombang `c`, dan jumlah grid solusi x dan t yaitu `m` dan `n`.
# Example
```jl
julia> f(x) = sin(pi*x)+sin(2*pi*x);

julia> g(x) = 0;

julia> U = gelombang(f,g,1,0.5,2,11,11)
11×11 Adjoint{Float64,Array{Float64,2}}:
 0.0   0.896802      1.53884     1.76007   …  -0.14204      -0.363271   -0.278768   0.0
 0.0   0.769421      1.32844     1.53884       1.11022e-16  -0.210404   -0.181636   0.0
 0.0   0.431636      0.769421    0.948401      0.360616      0.181636    0.0683644  0.0
 0.0  -2.22045e-16   0.0515989   0.181636      0.769421      0.639384    0.363271   0.0
 ⋮                                         ⋱                                        ⋮
 0.0  -0.363271     -0.639384   -0.769421     -0.181636     -0.0515989   0.0        0.0
 0.0  -0.0683644    -0.181636   -0.360616     -0.948401     -0.769421   -0.431636   0.0
 0.0   0.181636      0.210404    0.0          -1.53884      -1.32844    -0.769421   0.0
 0.0   0.278768      0.363271    0.14204   …  -1.76007      -1.53884    -0.896802   0.0
```
return solusi persamaan gelombang `U`.
"""
function gelombang(f,g,a,b,c,m,n)
  h = a/(m-1);
  x = 0:h:1;
  k = b/(n-1);
  r = c*k/h;
  U = zeros(m,n);
  for i = 2:m-1
    U[i,1] = f(x[i]);
    U[i,2] = (1-r^2)*f(x[i]) + k*g(x[i]) + r^2/2*(f(x[i+1])+f(x[i-1]));
  end
  for j = 3:n
    for i = 2:(m-1)
      U[i,j] = (2-2*r^2)*U[i,j-1]+r^2*(U[i-1,j-1]+U[i+1,j-1])-U[i,j-2];
    end
  end
  U = U';
end


text:
# Visualize output from running multiple experiments on regions of the brain
include("Plot_Helper.jl")
using Plots
using LaTeXStrings

epsi = .1
numregions = 17

for seedtype = 1:4

    prs1, res1, f1s1, Sizes, seednodes =  LoadFlowSeedTrain(epsi,seedtype)

    x_label = "Number of Nodes in Region"
    y_label1 = "Precision"
    y_label2 = "Recall"
    y_label3 = "F1 Score"

    l_place = :bottomright
    s1 = 400
    s2 = 300
    ms = 4
    lw = 2

    # F1 Plots
    plot(f1s1[1,:], labels = "No Penalties",ylim = [0,1],grid = false,size = (s1,s2),
        xlabel = x_label,ylabel = y_label3,legend = l_place,
        xticks = (1:numregions,Sizes), xrotation=50,linewidth = lw,
        color = :blue,markershape = :circle, markersize = ms)

    plot!(f1s1[2,:], labels = "Some Strict",size = (s1,s2),
        color = :red,markershape = :diamond,linewidth = lw,markersize = ms)

    plot!(f1s1[4,:], labels = "Strict+Soft",size = (s1,s2),
        color = :green,markershape = :square,linewidth = lw,markersize = ms)

    savefig("FlowSeed_Train_Seed_"*string(seedtype)*"_eps_"*string(epsi)*"_F1.pdf")

    # Recall Plots
    plot(res1[1,:], labels = "No Penalties",ylim = [0,1],grid = false,size = (s1,s2),
        xlabel = x_label,ylabel = y_label3,legend = l_place,
        xticks = (1:numregions,Sizes), xrotation=50,linewidth = lw,
        color = :blue,markershape = :circle, markersize = ms)

    plot!(res1[2,:], labels = "Some Strict",size = (s1,s2),
        color = :red,markershape = :diamond,linewidth = lw,markersize = ms)

    plot!(res1[4,:], labels = "Strict+Soft",size = (s1,s2),
        color = :green,markershape = :square,linewidth = lw,markersize = ms)

    savefig("FlowSeed_Train_Seed_"*string(seedtype)*"_eps_"*string(epsi)*"_Recall.pdf")

end

max_stars_repo_name: nveldt/FlowSeed

max_issues_repo_name: nveldt/FlowSeed

max_forks_repo_name: nveldt/FlowSeed


text:
module FS

using ..DataMod, ..Helpers
using MultivariateStats, Statistics, StatsBase, StatsModels

include("./manualmodel.jl")
export ManualModel, predict

module Utils
    using DataFrames, LinearAlgebra, StatsBase, StatsModels, ..DataMod, ..Helpers, ..ManualModelMod
    function compose(y::Symbol, terms::Vector{Symbol}; intercept::Bool = true)::FormulaTerm
        return term(y) ~ foldl(+, term.(intercept ? [1; terms] : terms))
    end
    function lt_featurename(f1::Symbol, f2::Symbol)
        s1, s2 = string(f1), string(f2)
        if s1[1] == s2[1]
            return isless(parse(Int, s1[2:end]), parse(Int, s2[2:end]))
        else
            return isless(s2[1], s1[1])
        end
    end
    function get_model(data::Data, sel_xmat::Matrix{Float64}, y_vec::Vector{Float64},
                       desel_i::Vector{Int}, predicted::Symbol, intercept::Bool)::RegressionModel
        # Simply doing `sel_xmat \ y_vec` resulted in a rare error when
        # sel_xmat is square. The \ operator interprets this incorrectly.
        # Therefore, this line was extracted from the \ operator's code:
        weights = qr(sel_xmat, ColumnNorm()) \ y_vec
        weights = copy(weights)
        for i in desel_i
            insert!(weights, intercept ? i + 1 : i, 0.0)
        end
        return ManualModel(weights, intercept, predicted, data)
    end
    function get_models(data::Data, rhs::Vector{Symbol}, intercept::Bool = true)::Vector{RegressionModel}
        desel_i = findall(x -> x ∉ rhs, data.xs)
        n = nrow(data.df)
        sort!(rhs, lt = lt_featurename)
        sel_xmat = isempty(rhs) ? zeros(Float64, n, 0) : (data.xs == rhs ? data.xmat : Matrix(select(data.xdf, rhs)))::Matrix{Float64}
        intercept && (sel_xmat = hcat(ones(Float64, n), sel_xmat))
        y_vecs = Vector{Float64}[data.ydf[!, y] for y in data.ys]
        return map(i -> get_model(data, sel_xmat, y_vecs[i], desel_i, data.ys[i], intercept), 1:length(y_vecs))
    end
    Base.iterate(s::Symbol) = (s, nothing)
    Base.iterate(::Symbol, ::Any) = nothing
    Base.length(s::Symbol) = 1
    export compose, get_models
end
using .Utils
export compose, get_models

include("./fs-stepwise.jl")
include("./fs-lasso.jl")
include("./fs-inipg.jl")
include("./fs-importance.jl")
include("./fs-ga.jl")
include("./fs-l1l2reg.jl")

llsqreg(data::Data)::Vector{RegressionModel} = get_models(data, data.xs)
function ridgereg(data::Data, λs::Vector{Float64} = .5 .^ collect(-5:14))::Vector{RegressionModel}
    ridges = [ridge(data.xmat, data.ymat, λ) for λ in λs]
    c = 1:length(data.ys)
    splits = [[ begin
                    y = r[:, i]
                    [y[end]; y[1:end-1]]
                end for i in c] for r in ridges]
    models = [[ManualModel(split[i], true, data.ys[i], data) for i in c] for split in splits]
    bics = [mean([bic(model[i]) for i in c]) for model in models]
    minind = argmin(bics)
    return models[minind]
end
export llsqreg, ridgereg

end

using .FS

max_stars_repo_name: KasperNooteboom/thesis-rvfl-fs

max_issues_repo_name: KasperNooteboom/thesis-rvfl-fs

max_forks_repo_name: KasperNooteboom/thesis-rvfl-fs


text:
const PWGRID_VERSION = 4

mutable struct GkVectors
    kpts::Array{Float64,2}
    Ngwx::Int
    Ngw::Array{Int,1}
    idx_gk::Array{Array{Int64,1},1}
end

mutable struct GVectors
    Ng::Int
    G::Array{Float64,2}
    G2::Array{Float64}
end

mutable struct PWGrid
    ecutwfc::Float64
    ecutrho::Float64
    Ns::Array{Int64}
    LatVecs::Array{Float64,2}
    RecVecs::Array{Float64,2}
    Ω::Float64
    r::Array{Float64,2}
    gvec::GVectors
    gkvec::GkVectors
end

function PWGrid( ecutwfc::Float64, LatVecs::Array{Float64,2}, kpts_red )
    ecutrho = 4.0*ecutwfc
    #
    println("LatVecs = ", LatVecs)
    RecVecs = 2*pi*inv(LatVecs)'
    println("RecVecs = ", RecVecs)
    Ω = det(LatVecs)
    #
    LatVecsLen = Array{Float64}(undef,3)
    LatVecsLen[1] = norm(LatVecs[1,:])
    LatVecsLen[2] = norm(LatVecs[2,:])
    LatVecsLen[3] = norm(LatVecs[3,:])

    Ns = Array{Int64}(undef,3)
    Ns[1] = 2*round( Int64, sqrt(ecutrho/2)*LatVecsLen[1]/pi ) + 1
    Ns[2] = 2*round( Int64, sqrt(ecutrho/2)*LatVecsLen[2]/pi ) + 1
    Ns[3] = 2*round( Int64, sqrt(ecutrho/2)*LatVecsLen[3]/pi ) + 1

    # Use even sampling numbers
    Ns[1] = Ns[1] % 2 == 1 ? Ns[1] + 1 : Ns[1]
    Ns[2] = Ns[2] % 2 == 1 ? Ns[2] + 1 : Ns[2]
    Ns[3] = Ns[3] % 2 == 1 ? Ns[3] + 1 : Ns[3]

    Npoints = prod(Ns)
    r = init_grid_R( Ns, LatVecs )

    gvec = init_grid_G( Ns, RecVecs )

    #kpts = (kpts_red' * RecVecs)'

    kpts = RecVecs*kpts_red # new following convention of PWDFT.jl

    #kpts = (kpts_red' * RecVecs')'

    #kpts = RecVecs' * kpts_red
    #kpts = RecVecs * kpts_red

    #Nkpts = size(kpts_red)[2]
    #kpts = zeros(Float64, 3, Nkpts)
    #for ik = 1:Nkpts
    #    kpts[:,ik] = RecVecs' * kpts_red[:,ik]
    #end
    #for ik = 1:Nkpts
    #    for i = 1:3
    #        kpts[i,ik] = kpts_red[i,ik] * RecVecs[1,:]
    #    end
    #end
    gkvec = init_gkvec( ecutwfc, gvec.G, kpts )

    return PWGrid( ecutwfc, ecutrho, Ns, LatVecs, RecVecs, Ω, r, gvec, gkvec )
end

function mm_to_nn(mm::Int,S::Int)
    if mm > S/2
        return mm - S
    else
        return mm
    end
end

function init_grid_G( Ns, RecVecs )

    Ng = prod(Ns)

    G  = Array{Float64}(undef,3,Ng)
    G2 = Array{Float64}(undef,Ng)

    ig = 0
    for k in 0:Ns[3]-1
    for j in 0:Ns[2]-1
    for i in 0:Ns[1]-1
        ig = ig + 1
        gi = mm_to_nn( i, Ns[1] )
        gj = mm_to_nn( j, Ns[2] )
        gk = mm_to_nn( k, Ns[3] )
        G[1,ig] = RecVecs[1,1]*gi + RecVecs[1,2]*gj + RecVecs[1,3]*gk
        G[2,ig] = RecVecs[2,1]*gi + RecVecs[2,2]*gj + RecVecs[2,3]*gk
        G[3,ig] = RecVecs[3,1]*gi + RecVecs[3,2]*gj + RecVecs[3,3]*gk
        G2[ig] = G[1,ig]^2 + G[2,ig]^2 + G[3,ig]^2
    end
    end
    end

    return GVectors( Ng, G, G2 )
end

function init_gkvec( ecutwfc, G, kpts )
    Nkpts = size(kpts)[2]
    idx_gk = Array{Array{Int64,1},1}(undef,Nkpts)
    Ng = size(G)[2]
    Gk = zeros( 3, Ng )
    Gk2 = zeros( Ng )
    Ngw = zeros( Int64, Nkpts )
    for ik = 1:Nkpts
        #@printf("\n")
        for ig = 1:Ng
            Gk[:,ig] = G[:,ig] + kpts[:,ik]
            Gk2[ig]  = Gk[1,ig]^2 + Gk[2,ig]^2 + Gk[3,ig]^2
            #@printf("%8d %18.10f %18.10f\n", ig, Gk2[ig], ecutwfc)
        end
        idx_gk[ik] = findall( 0.5*Gk2 .< ecutwfc )
        Ngw[ik] = length(idx_gk[ik])
        #@printf("kpt = %d, Ngw = %d\n", ik, Ngw[ik])
        #@printf("kpts[%d]: %f %f %f\n", ik, kpts[1,ik], kpts[3,ik], kpts[3,ik] )
    end
    Ngwx = maximum(Ngw)
    return GkVectors( kpts, Ngwx, Ngw, idx_gk )
end

function init_grid_R( Ns, LatVecs )
    #
    Npoints = prod(Ns)
    #
    R = Array{Float64}(undef,3,Npoints)
    ip = 0
    for k in 0:Ns[3]-1
    for j in 0:Ns[2]-1
    for i in 0:Ns[1]-1
        ip = ip + 1
        R[1,ip] = LatVecs[1,1]*i/Ns[1] + LatVecs[1,2]*j/Ns[2] + LatVecs[1,3]*k/Ns[3]
        R[2,ip] = LatVecs[2,1]*i/Ns[1] + LatVecs[2,2]*j/Ns[2] + LatVecs[2,3]*k/Ns[3]
        R[3,ip] = LatVecs[3,1]*i/Ns[1] + LatVecs[3,2]*j/Ns[2] + LatVecs[3,3]*k/Ns[3]
    end
    end
    end
    #
    return R
end


text:
# Lotka-Volterra, np = 2
id = :Sahlodin2011c

pL = [2.98; 0.98]
pU = [3.02; 1.02]
Sahlodin2011c_x0(p::Vector{T}) where T = T[1.2*one(T); 1.1*one(T)]
function Sahlodin2011c_f!(du::Vector{T}, u::Vector{T}, p::Vector{T}, t) where T
    du[1] = p[1]*u[1]*(1.0 - u[2])
    du[2] = p[2]*u[2]*(u[1] - 1.0)
    return
end
tspan = (0.0, 2.0)

prob = ODERelaxProb(Sahlodin2011c_f!, tspan, Sahlodin2011c_x0, pL, pU)
url = "https://www.sciencedirect.com/science/article/abs/pii/S0168927411000316"
source = "Sahlodin2011"
desc = "Lotka-Volterra, np = 2"


text:
# check for 45º true-fitted plots for all parameters
# to verify that stats give good identification
using Pkg
Pkg.activate("../../")
using Plots, Econometrics
mcreps = 500 # how many reps?

# load code
include("JDlib.jl")
include("../../src/SNM.jl")
using BSON:@load

function main()
@load "neural_moments.bson" NNmodel transform_stats_info
trueθ = zeros(mcreps,8)
θhat = zeros(mcreps,8)
Threads.@threads for mcrep = 1:mcreps
    println("mc rep: ", mcrep)
    θ = PriorDraw()
    trueθ[mcrep,:] = θ
    θhat[mcrep,:] = NeuralMoments(θ, auxstat, 1, NNmodel, transform_stats_info)
end
return trueθ, θhat
end
trueθ, θhat = main()
savefig(scatter(trueθ[:,1], θhat[:,1]),"mu.svg")
savefig(scatter(trueθ[:,2], θhat[:,2]),"kappa.svg")
savefig(scatter(trueθ[:,3], θhat[:,3]),"alpha.svg")
savefig(scatter(trueθ[:,4], θhat[:,4]),"sigma.svg")
savefig(scatter(trueθ[:,5], θhat[:,5]),"rho.svg")
savefig(scatter(trueθ[:,6], θhat[:,6]),"lambda0.svg")
savefig(scatter(trueθ[:,7], θhat[:,7]),"lambda1.svg")
savefig(scatter(trueθ[:,8], θhat[:,8]),"tau.svg")
for i = 1:8
    ols(θhat[:,i],[ones(mcreps) trueθ[:,i]])
end

max_stars_repo_name: mcreel/ILSNM

max_issues_repo_name: mcreel/ILSNM

max_forks_repo_name: mcreel/ILSNM


text:
module SimulacoesPedro
    # Data
    include("Data.jl")
    export Data

    # Initial point and velocity
    include("initial_point.jl")
    include("initial_velocity.jl")
    export initial_point, initial_velocity

    # Naive potential
    include("rpbc.jl")
    include("upair.jl")
    include("utotal.jl")
    include("linear_index.jl")
    include("utotal_parallel.jl")
    export utotal
    export utotal_parallel1
    export utotal_parallel2
    export utotal_serial2

    # Linked List potential
    include("makell.jl")
    include("wrap_cell.jl")
    include("r2pbc.jl")
    include("upair2.jl")
    include("upair3.jl")
    include("utotal2.jl")
    include("utotal2_generic.jl")
    include("utotal2_parallel.jl")
    include("utotal3.jl")
    include("runsim.jl")
    export utotal2, utotal2_parallel, utotal3!, runsim, runsim!

end


text:
config[@__DIR__] = (
    directory=@__DIR__,
    case="GaAs",
    n_repeat=40,
    kwargs_scf=(randomise_guess=true, random_ampl=1e-3, ),
    firstmixkey="RestaGaAs",
    # firstα=0.1,
    #
    construct_localiser=(lattice, atoms) -> identity,
    mixings=Dict(
        "Kerker"     => (α, localiser) -> KerkerMixing(α=α, kTF=1.0),
        "None"       => (α, localiser) -> SimpleMixing(α=α),
        "RestaGaAs"  => (α, localiser) -> RestaMixing( α=α, kTF=1.0, εr=14),
        "RestaSiO2"  => (α, localiser) -> RestaMixing( α=α, kTF=1.0, εr=2.4),
    ),
)


text:
INTERACTION_DEPVAR_ERROR = "depvar should not be selected for interactions"
INTERACTION_EQUATION_ERROR = "some interaction variable is not in the original equation"
LOG_FUNCTION_ERROR = "log transformation do not allow negative values"
SOME_VARIABLES_NOT_FOUND = "some of the selected variables are not found"


text:
using Sodium
using Sodium.LibSodium
using Test
using Base64

@testset "Sodium" begin
    @test sodium_init() ≥ 0
    @info "libsodium library version : " * unsafe_string(sodium_version_string())
end

@testset "Sealed boxes" begin
    msg = "Message"
    msg_len = length(msg)
    ciphertext_len = crypto_box_SEALBYTES + msg_len
    recipient_pk = Vector{Cuchar}(undef, crypto_box_PUBLICKEYBYTES)
    recipient_sk = Vector{Cuchar}(undef, crypto_box_SECRETKEYBYTES)
    @test crypto_box_keypair(recipient_pk, recipient_sk) == 0
    ciphertext = Vector{Cuchar}(undef, ciphertext_len)
    @test crypto_box_seal(ciphertext, msg, msg_len, recipient_pk) == 0
    decrypted = Vector{Cuchar}(undef, msg_len)
    @test crypto_box_seal_open(decrypted, ciphertext, ciphertext_len, recipient_pk, recipient_sk) == 0
    @test String(decrypted) == msg
end

@testset "Sealed boxes convenience" begin
    msg = "My Message with true Unicode öéâù"
    recipient_pk = Vector{Cuchar}(undef, crypto_box_PUBLICKEYBYTES)
    recipient_sk = Vector{Cuchar}(undef, crypto_box_SECRETKEYBYTES)
    @test crypto_box_keypair(recipient_pk, recipient_sk) == 0
    encrypted = seal(msg, base64encode(recipient_pk))
    ciphertext = base64decode(encrypted)
    decrypted = Vector{Cuchar}(undef, sizeof(ciphertext)-crypto_box_sealbytes())
    @test crypto_box_seal_open(decrypted, ciphertext, sizeof(ciphertext), recipient_pk, recipient_sk) == 0
    @test String(decrypted) == msg
end


text:
function uniquearrayrefs_csesummary(ls::LoopSet)
  uniquerefs = ArrayReferenceMeta[]
  # each `Vector{Tuple{Int,Int}}` has the same name
  # 1 to ≥1 maps:
  name_to_array_map = Vector{Int}[] # refname -> uniquerefs
  unique_to_name_and_op_map = Vector{Tuple{Int,Int,Int}}[] # uniquerefs -> (ind in name_to_array, ind in that vector, ind in operations)
  for (i,op) ∈ enumerate(operations(ls))
    arrayref = op.ref
    arrayref === NOTAREFERENCE && continue
    unique = true
    for (j,namev) ∈ enumerate(name_to_array_map)
      vptr(uniquerefs[first(namev)]) === vptr(arrayref) || continue # no unique
      for (k,unique_id) ∈ enumerate(namev)
        if sameref(uniquerefs[unique_id], arrayref)
          # matching name, matching ref
          push!(unique_to_name_and_op_map[unique_id], (j, k, i))
          unique = false
          break
        end
      end
      if unique # matching name, no matching ref
        push!(uniquerefs, arrayref)
        push!(namev, length(uniquerefs))
        push!(unique_to_name_and_op_map, Tuple{Int,Int,Int}[(j,length(namev), i)])
        unique = false
      end
      break
    end
    if unique # no matching name, no matching ref
      push!(uniquerefs, arrayref)
      push!(name_to_array_map, Int[length(uniquerefs)])
      push!(unique_to_name_and_op_map, Tuple{Int,Int,Int}[( length(name_to_array_map), 1, i)])
    end
  end
  uniquerefs, name_to_array_map, unique_to_name_and_op_map
end

function uniquearrayrefs(ls::LoopSet)
    uniquerefs = ArrayReferenceMeta[]
    includeinlet = Bool[]
    # for arrayref ∈ ls.refs_aliasing_syms
    for op ∈ operations(ls)
        arrayref = op.ref
        arrayref === NOTAREFERENCE && continue
        notunique = false
        isonlyname = true
        for ref ∈ uniquerefs
            notunique = sameref(arrayref, ref)
            isonlyname &= vptr(arrayref) !== vptr(ref)
            # if they're not the sameref, they may still have the same name
            # if they have different names, they're definitely not sameref
            notunique && break
        end
        if !notunique
            push!(uniquerefs, arrayref)
            push!(includeinlet, isonlyname)
        end
    end
    uniquerefs, includeinlet
end

otherindexunrolled(loopsym::Symbol, ind::Symbol, loopdeps::Vector{Symbol}) = (loopsym !== ind) && (loopsym ∈ loopdeps)
function otherindexunrolled(ls::LoopSet, ind::Symbol, ref::ArrayReferenceMeta)
    us = ls.unrollspecification
    u₁sym = names(ls)[us.u₁loopnum]
    u₂sym = us.u₂loopnum > 0 ? names(ls)[us.u₂loopnum] : Symbol("##undefined##")
    otherindexunrolled(u₁sym, ind, loopdependencies(ref)) || otherindexunrolled(u₂sym, ind, loopdependencies(ref))
end
function multiple_with_name(n::Symbol, v::Vector{ArrayReferenceMeta})
  found = false
  for ref ∈ v
    tst = vptr(ref) === n
    (found & tst) && return true
    found |= tst
  end
  false
end
# multiple_with_name(n::Symbol, v::Vector{ArrayReferenceMeta}) = sum(ref -> n === vptr(ref), v) > 1
# TODO: DRY between indices_calculated_by_pointer_offsets and use_loop_induct_var
function indices_calculated_by_pointer_offsets(ls::LoopSet, ar::ArrayReferenceMeta)
    indices = getindices(ar)
    ls.isbroadcast && return fill(false, length(indices))
    looporder = names(ls)
    offset = isdiscontiguous(ar)
    gespinds = Expr(:tuple)
    out = Vector{Bool}(undef, length(indices))
    li = ar.loopedindex
    for i ∈ eachindex(li)
        ii = i + offset
        ind = indices[ii]
        if (!li[i]) || (ind === CONSTANTZEROINDEX) || multiple_with_name(vptr(ar), ls.lssm.uniquearrayrefs) ||
            (iszero(ls.vector_width) && isstaticloop(getloop(ls, ind)))# ||
            out[i] = false
        elseif (isone(ii) && (first(looporder) === ind))
            out[i] = otherindexunrolled(ls, ind, ar)
        else
            out[i] = true
        end
    end
    out
end

# @generated function set_first_stride(sptr::StridedPointer{T,N,C,B,R}) where {T,N,C,B,R}
#     minrank = argmin(R)
#     newC = C > 0 ? (C == minrank ? 1 : 0) : C
#     newB = C > 0 ? (C == minrank ? B : 0) : B #TODO: confirm correctness
#     quote
#         $(Expr(:meta,:inline))
#         # VectorizationBase.StridedPointer{$T,1,$newC,$newB,$(R[minrank],)}($(lv(llvmptr))(sptr), (sptr.strd[$minrank],), (Zero(),))
#         VectorizationBase.StridedPointer{$T,1,$newC,$newB,$(R[minrank],)}(VectorizationBase.cpupointer(sptr), (sptr.strd[$minrank],), (Zero(),))
#     end
# end
# set_first_stride(x) = x # cross fingers that this works
# @inline onetozeroindexgephack(sptr::AbstractStridedPointer) = gesp(set_first_stride(sptr), (Static{-1}(),)) # go backwords
# @inline onetozeroindexgephack(sptr::AbstractStridedPointer{T,1}) where {T} = sptr
# @inline onetozeroindexgephack(sptr::StridedPointer{T,1}) where {T} = sptr
# @inline onetozeroindexgephack(x) = x

# # Removes parent/child relationship for all children with ref `ar`
# function freechildren!(op::Operation, ar::ArrayReferenceMeta)
#   newchildren = Operation[]
#   for opc ∈ children(op)
#     if opc.ref === ar
#       for (i,opp) ∈ enumerate(parents(opc))
#         if opp === op
#           deleteat!(parents(opc), i)
#           break
#         end
#       end
#     else
#       push!(newchildren, opc)
#     end
#   end
#   op.children = newchildren
#   nothing
# end
# function replaceparent!(ops::Vector{Operation}, )

# end
function set_ref_loopedindex_and_ind!(ref::ArrayReferenceMeta, i::Int, ii::Int, li::Bool, ind::Symbol)
    ref.loopedindex[i] = li
    getindices(ref)[ii] = ind
end
function set_all_to_constant_index!(
  ls::LoopSet, i::Int, ii::Int, indop::Operation, allarrayrefs::Vector{ArrayReferenceMeta},
  array_refs_with_same_name::Vector{Int}, arrayref_to_name_op_collection::Vector{Vector{Tuple{Int,Int,Int}}}
)
  ops = operations(ls)
  for j ∈ array_refs_with_same_name
    arrayref_to_name_op = arrayref_to_name_op_collection[j]
    set_ref_loopedindex_and_ind!(allarrayrefs[j], i, ii, true, CONSTANTZEROINDEX)
    for (_,__,opid) ∈ arrayref_to_name_op
      op = ops[opid]
      set_ref_loopedindex_and_ind!(op.ref, i, ii, true, CONSTANTZEROINDEX)
      delete_ind = 0
      for (k,opp) ∈ enumerate(parents(op))
        if opp === indop
          delete_ind = k
        end
      end
      # @assert delete_ind ≠ 0
      deleteat!(parents(op), delete_ind)
    end
  end
end
maybeloopvaluename(op::Operation) = isloopvalue(op) ? instruction(op).instr : name(op)
function substitute_ops_all!(
  ls::LoopSet, i::Int, ii::Int, indop::Operation, new_parent::Operation, allarrayrefs::Vector{ArrayReferenceMeta},
  array_refs_with_same_name::Vector{Int}, arrayref_to_name_op_collection::Vector{Vector{Tuple{Int,Int,Int}}}
)
  newindsym = maybeloopvaluename(new_parent)
  isloopval = isloopvalue(new_parent)
  ops = operations(ls)
  for j ∈ array_refs_with_same_name
    arrayref_to_name_op = arrayref_to_name_op_collection[j]
    set_ref_loopedindex_and_ind!(allarrayrefs[j], i, ii, isloopval, newindsym)
    for (_,__,opid) ∈ arrayref_to_name_op
      op = ops[opid]
      set_ref_loopedindex_and_ind!(op.ref, i, ii, isloopval, newindsym)
      sub_ind = 0
      for (k,opp) ∈ enumerate(parents(op))
        if opp === indop
          sub_ind = k
        end
      end
      # @assert sub_ind ≠ 0
      if isloopval
        deleteat!(parents(op), sub_ind)
      else
        parents(op)[sub_ind] = new_parent
      end
    end
  end
end
# function normalize_offsets!(
#   ls::LoopSet, i::Int, allarrayrefs::Vector{ArrayReferenceMeta},
#   array_refs_with_same_name::Vector{Int}, arrayref_to_name_op_collection::Vector{Vector{Tuple{Int,Int,Int}}}
# )
#   ops = operations(ls)
#   length(ops) > 256 && return 0
#   minoffset::Int8 = typemax(Int8)
#   maxoffset::Int8 = typemin(Int8)
#   # we want to store the offsets, because we don't want to require that the `offset` vectors of the variaous `ArrayReferenceMeta`s don't alias
#   # loopsym = Symbol("##DUMMY##NOT#REALLY#A#LOOP##")
#   # stride::Int = typemin(Int)
#   # thereiszerooffset::Bool = false
#   # offsets::Base.RefValue{NTuple{128,Int8}} = Base.RefValue{NTuple{128,Int8}}();
#   # GC.@preserve offsets begin
#   #   poffsets = Base.unsafe_convert(Ptr{Int8}, offsets)
#   for j ∈ array_refs_with_same_name
#     arrayref_to_name_op = arrayref_to_name_op_collection[j]
#     for (_,__,opid) ∈ arrayref_to_name_op
#       op = ops[opid]
#       opref = op.ref
#       off = getoffsets(opref)[i]
#       # thereiszerooffset |= off == zero(Int8)
#       # off == zero(Int8) && return 0
#       minoffset = min(off, minoffset)
#       maxoffset = max(off, maxoffset)
#       # unsafe_store!(poffsets, off, opid)
#       # if loopsym ≢ Symbol("##DUMMY##NOT#REALLY#A#LOOP##")
#       #   stride = Int(getstrides(opref)[i])
#       #   if opref.loopedindex[i]
#       #     loopsym = getindicesonly(op)[i]
#       #   else
#       #     loopsym = Symbol("##NOT#A#LOOP##")
#       #   end
#       # end
#     end
#   end
#   # reaching here means none of the offsets contain `0`
#   # we won't bother if difference between offsets is >127
#   # we don't want `maxoffset` to overflow when subtracting `minoffset`
#   # so we check if it's safe, and give up if it isn't
#   minoffsetint = Int(minoffset)
#   return (((Int(maxoffset) - minoffsetint) > 127)) ? 0 : minoffsetint

#   #   # if loopsym ≢ Symbol("##DUMMY##NOT#REALLY#A#LOOP##")
#   #   #   loop = getloop(ls, loopsym)
#   #   #   if minstride ≠ maxstride
#   #   #     @assert isknown(first(loop)) "Currently, if the same index is used for the same array with different multiples (e.g., `x[i]` and `x[2*i]`), then the start of that loop range must be known at compile time."
#   #   #   end
#   #   #   stride_offset = 1 - gethint(first(loop))
#   #   # else
#   #   #   loop = first(ls.loops)
#   #   # end
#   #   offset_adjust = Int(minoffset)
#   #   if (loopsym ≢ Symbol("##DUMMY##NOT#REALLY#A#LOOP##")) && (loopsym ≢ Symbol("##NOT#A#LOOP##"))
#   #     loopstart = first(getloop(ls, loopsym))
#   #     if isknown(loopstart)
#   #       loopstartval = gethint(loopstart)
#   #       offset_adjust = loopstartval*(stride - 1)
#   #       if offset_adjust + Int(maxoffset) ≤ typemax(Int8)
#   #         minoffset -= Int8(offset_adjust)
#   #         stride = 1
#   #       end
#   #     end
#   #   end
#   #   for j ∈ array_refs_with_same_name
#   #     arrayref_to_name_op = arrayref_to_name_op_collection[j]
#   #     for (_,__,opid) ∈ arrayref_to_name_op
#   #       new_offset = unsafe_load(poffsets, opid) - minoffset
#   #       old_offset = getoffsets(ops[opid].ref)[i]
#   #       @show new_offset, old_offset
#   #       getoffsets(ops[opid].ref)[i] = new_offset
#   #       # getoffsets(ops[opid].ref)[i] = unsafe_load(poffsets, opid) - minoffset
#   #     end
#   #   end
#   # end
#   # return @show offset_adjust, stride
# end
function isloopvalue(ls::LoopSet, ind::Symbol, isrooted::Union{Nothing,Vector{Bool}} = nothing)
  for (i,op) ∈ enumerate(operations(ls))
    if (isrooted ≢ nothing)
      isrooted[i] || continue
    end
    iscompute(op) || continue
    for opp ∈ parents(op)# this is to confirm `ind` still has children
      # (isloopvalue(opp) && instruction(opp).instr === ind) && return true
      if (isloopvalue(opp) && instruction(opp).instr === ind)
        return true
      end
    end
  end
  return false
end
function cse_constant_offsets!(
  ls::LoopSet, allarrayrefs::Vector{ArrayReferenceMeta}, allarrayrefsind::Int, name_to_array_map::Vector{Vector{Int}}, arrayref_to_name_op_collection::Vector{Vector{Tuple{Int,Int,Int}}}
)
  ar = allarrayrefs[allarrayrefsind]
  # vptrar = vptr(ar)
  arrayref_to_name_op = arrayref_to_name_op_collection[allarrayrefsind]
  array_refs_with_same_name = name_to_array_map[first(first(arrayref_to_name_op))]
  us = ls.unrollspecification
  li = ar.loopedindex
  indices = getindices(ar)
  strides = getstrides(ar)
  offset = first(indices) === DISCONTIGUOUS
  # gespindoffsets = fill(Symbol(""), length(li))
  gespindsummary = Vector{Symbol}(undef, length(li))
  for i ∈ eachindex(li)
    gespsymbol::Symbol = Symbol("")
    ii = i + offset
    ind::Symbol = indices[ii]
    licmoffset = !li[i]
    if licmoffset
      # we try to licm and offset so we can set `li[i] = true`
      position_in_array_refs_with_same_name = first(arrayref_to_name_op)[2]
      # if position_in_array_refs_with_same_name ≠ 1, then we already performed the op substitutions
      if length(array_refs_with_same_name) > 1 # if it == 1, then there's only 1
        for j ∈ array_refs_with_same_name
          j == position_in_array_refs_with_same_name && continue
          ref = allarrayrefs[j]
          refinds = getindices(ref)
          # refinds === indices && continue # fast check, should be covered by `j == position_in_array_refs_with_same_name`
          if !((refinds[ii] === ind) & (getstrides(ar)[i] == getstrides(ref)[i]))
            # For now, we'll only bother with `licm` if all share the same indices
            # This is so that we can apply the same `licm` to each and keep the same array name.
            # Otherwise, we'll rely on LLVM to optimize indexing.
            licmoffset = false
            break
          end
        end
      end
      ops = operations(ls)
      # @show ar
      while licmoffset # repeat until we run out
        # ind = indices[ii]
        # indices are all the same across operations, so we look to the first for checking compatibility...
        memop = indop = ops[first(arrayref_to_name_op)[3]] # indop is a dummy placeholder
        for opp ∈ parents(memop)
          if name(opp) === ind
            indop = opp
          end
        end
        # (we found a match) # iscompute(indop) should be guaranteed...
        ((indop ≢ memop) && iscompute(indop)) || break
        instr = instruction(indop).instr
        parents_indop = parents(indop)
        if instr === :(+)
          constants_to_licm = if gespsymbol === Symbol("")
            Expr(:call, GlobalRef(Base,:(+)))
          else
            Expr(:call, GlobalRef(Base,:(+)), gespsymbol)
          end
          new_parent = indop # dummy, for now we will only licm if it lets us remove `indop`
          for opp ∈ parents_indop
            if isconstantop(opp)
              # NOTE: using `name` because this happens at macroexpansion time
              # push!(constants_to_licm.args, constantopname(opp))
              push!(constants_to_licm.args, name(opp))
            elseif new_parent === indop # first
              new_parent = opp
            else # this is the second parent, we give up on `licm`
              licmoffset = false
              break
            end
          end
          licmoffset || continue
          if length(constants_to_licm.args) > 2
            gespsymbol = gensym!(ls, "#gespsym#")
            pushpreamble!(ls, Expr(:(=), gespsymbol, constants_to_licm))
          elseif length(constants_to_licm.args) == 2
            licmoffset = gespsymbol === Symbol("")
            if licmoffset
              gespsymbol = (constants_to_licm.args[2])::Symbol
            end
          else
            licmoffset = false
          end
          if licmoffset # prune
            if new_parent === indop # no parents left
              ind = CONSTANTZEROINDEX
              set_all_to_constant_index!(ls, i, ii, indop, allarrayrefs, array_refs_with_same_name, arrayref_to_name_op_collection)
            else # new_parent is a new parent to replace `indop`
              ind = maybeloopvaluename(new_parent)
              substitute_ops_all!(ls, i, ii, indop, new_parent, allarrayrefs, array_refs_with_same_name, arrayref_to_name_op_collection)
            end
          end
        elseif (instr === :(-)) && length(parents_indop) == 2
          op1 = parents(indop)[1]
          op2 = parents(indop)[2]
          op1const = isconstantop(op1)
          op2const = isconstantop(op2)
          if op1const # op1 - op2
            if op2const
              subexpr = if gespsymbol === Symbol("")
                Expr(:call, GlobalRef(Base,:(-)), name(op1), name(op2))
              else
                Expr(:call, GlobalRef(Base,:(-)), Expr(:call, GlobalRef(Base,:(+)), gespsymbol, name(op1)), name(op2))
              end
              gespsymbol = gensym!(ls, "#gespsym#")
              pushpreamble!(ls, Expr(:(=), gespsymbol, subexpr))
              ind = CONSTANTZEROINDEX
              set_all_to_constant_index!(ls, i, ii, indop, allarrayrefs, array_refs_with_same_name, arrayref_to_name_op_collection)
            else# op1const, op2dynamic
              # won't bother with this for now
              licmoffset = false
            end
          elseif op2const #op1dynamic
            subexpr = if gespsymbol === Symbol("")
              Expr(:call, GlobalRef(Base,:(-)), name(op2))
            else
              Expr(:call, GlobalRef(Base,:(-)), gespsymbol, name(op2))
            end
            gespsymbol = gensym!(ls, "#gespsym#")
            pushpreamble!(ls, Expr(:(=), gespsymbol, subexpr))
            ind = maybeloopvaluename(op1)
            substitute_ops_all!(ls, i, ii, indop, op1, allarrayrefs, array_refs_with_same_name, arrayref_to_name_op_collection)
          else
            licmoffset = false
          end
        else
          licmoffset = false
        end
      end
    end
    # constoffset = normalize_offsets!(ls, i, allarrayrefs, array_refs_with_same_name, arrayref_to_name_op_collection)
    # gespindsummary[i] = (gespsymbol, constoffset)
    gespindsummary[i] = gespsymbol
  end
  return gespindsummary
end
@inline similardims(_, i) = i
@inline similardims(::CartesianIndices{N}, i) where {N} = VectorizationBase.CartesianVIndex(ntuple(_ -> i, Val{N}()))
# @generated function similardims(::CartesianIndices{N}, i) where {N}
#   t = Expr(:tuple)
#   for n ∈ 1:N
#     push!(t.args, :i)
#   end
#   Expr(:block,Expr(:meta,:inline),t)
# end
# function maybebroadcastpush!(q::Expr, loop::Loop, b::Bool, ex)
#   if b
#     push!(q.args, Expr(:call, lv(:similardims), loop.rangesym, ex))
#   else
#     push!(q.args, ex)
#   end
#   return nothing
# end

function adjust_offsets!(
  ls::LoopSet, i::Int,
  array_refs_with_same_name::Vector{Int}, arrayref_to_name_op_collection::Vector{Vector{Tuple{Int,Int,Int}}}
)
  ops = operations(ls)
  if length(ops) ≤ 256
    offsets = Ref{NTuple{256,Int8}}()
    poffsets = Base.unsafe_convert(Ptr{Int8}, offsets)
    GC.@preserve offsets adjust_offsets!(ls, i, poffsets, array_refs_with_same_name, arrayref_to_name_op_collection)
  else
    offsetsv = similar(ops, Int8)
    poffsets = pointer(offsetsv)
    GC.@preserve offsetsv adjust_offsets!(ls, i, poffsets, array_refs_with_same_name, arrayref_to_name_op_collection)
  end
end
function adjust_offsets!(
  ls::LoopSet, i::Int, poffsets::Ptr{Int8},
  array_refs_with_same_name::Vector{Int}, arrayref_to_name_op_collection::Vector{Vector{Tuple{Int,Int,Int}}}
)
  ops = operations(ls)
  minoffset = typemax(Int8)
  maxoffset = typemin(Int8)
  # stridesunequal = false
  for j ∈ array_refs_with_same_name
    arrayref_to_name_op = arrayref_to_name_op_collection[j]
    for (_,__,opid) ∈ arrayref_to_name_op
      opref = ops[opid].ref
      off = getoffsets(opref)[i]
      minoffset = min(off, minoffset)
      maxoffset = max(off, maxoffset)
      unsafe_store!(poffsets, off, opid)
      # stridesunequal |= (stride ≠ getstrides(opref)[i])
    end
  end
  constoffset = Int(minoffset)
  constoffset = Core.ifelse(Int(maxoffset) - constoffset > 127, 0, constoffset)
  if constoffset ≠ 0
    for j ∈ array_refs_with_same_name
      arrayref_to_name_op = arrayref_to_name_op_collection[j]
      for (_,__,opid) ∈ arrayref_to_name_op
        opref = ops[opid].ref
        newoffset = unsafe_load(poffsets, opid) - constoffset
        # if stridesunequal
        #   stride = getstrides(opref)[i]
        #   newoffsetint = Int(newoffset) + (Int(stride) - 1)
        #   # @assert typemin(Int8) ≤ newoffsetint ≤ typemax(Int8)
        #   newoffset = Int8(newoffsetint)
        # end
        getoffsets(ops[opid].ref)[i] = newoffset
      end
    end
  end
  return constoffset#, Core.ifelse(stridesunequal, 1, Int(stride))
end

function calcgespinds(
  ls::LoopSet, ar::ArrayReferenceMeta, gespindsummary::Vector{Symbol}, shouldindbyind::Vector{Bool},
  array_refs_with_same_name::Vector{Int}, arrayref_to_name_op_collection::Vector{Vector{Tuple{Int,Int,Int}}}
)
  gespinds = Expr(:tuple)
  li = ar.loopedindex
  indices = getindicesonly(ar)
  # offsets = getoffsets(ar)
  strides = getstrides(ar)
  for i ∈ eachindex(li)
    ind = indices[i]
    isli = li[i]
    gespsymbol = gespindsummary[i]
    # if isli & (!index_by_index) && (length(operations(ls)) ≤ 256)
    #   ops = operations(ls)
    #   loopfirst = first(getloop(ls, ind))
    #   if isknown(loopfirst)
    #     # copy in case of aliasing
    #   end
    # end
    # constoffset ≠ 0 &&
    constoffset = adjust_offsets!(ls, i, array_refs_with_same_name, arrayref_to_name_op_collection)
    index_by_index = isli ? check_shouldindbyind(ls, ind, shouldindbyind) : true
    # (stridesunequal & isli) && (@assert isknown(first(getloop(ls, ind))))

    # end
    #   stride = strides[i]
    #   if stride ≠ 1
    #     loop = getloop(ls, ind)
    #     if isknown(first(loop))
    #       offsets[i] -= gethint(first(loop))
    #     end
    #   end
    #   # for op ∈ operations(ls)
    #   #   accesses_memory(op) || continue
    #   #   sameref(op.ref, ref) || continue
    #   #   getstrides(op)
    #   # end
    # end
    pushgespind!(gespinds, ls, gespsymbol, constoffset, Int(strides[i]), ind, isli, index_by_index, true)
  end
  gespinds
end

function pushgespind!(
  gespinds::Expr, ls::LoopSet, gespsymbol::Symbol, constoffset::Int, stride::Int, ind::Symbol, isli::Bool, index_by_index::Bool, fromgsp::Bool
)
  if isli
    if ind === CONSTANTZEROINDEX
      if gespsymbol === Symbol("")
        push!(gespinds.args, staticexpr(constoffset))
      elseif constoffset == 0
        push!(gespinds.args, gespsymbol)
      else
        push!(gespinds.args, Expr(:call, GlobalRef(Base,:(+)), gespsymbol, staticexpr(constoffset)))
      end
    else
      if index_by_index
        if gespsymbol === Symbol("")
          if constoffset == 0
            nostep = if fromgsp | (!index_by_index)
              GlobalRef(VectorizationBase, :NullStep)
            else
              GlobalRef(ArrayInterface.Static, :Zero)
            end
            ns = Expr(:call, nostep)
            if fromgsp
              loop = getloop(ls, ind)
              if loop.rangesym ≢ Symbol("")
                ns = Expr(:call, lv(:similardims), loop.rangesym, ns)
              end
            end
            push!(gespinds.args, ns)
          else
            push!(gespinds.args, staticexpr(constoffset))
          end
        elseif constoffset == 0
          push!(gespinds.args, gespsymbol)
        else
          push!(gespinds.args, addexpr(gespsymbol, constoffset))
        end
      else
        loop = getloop(ls, ind)
        if gespsymbol === Symbol("")
          if isknown(first(loop))
            # @show constoffset, gethint(first(loop))
            push!(gespinds.args, staticexpr(constoffset + stride*gethint(first(loop))))
          elseif constoffset == 0
            if stride == 1
              push!(gespinds.args, getsym(first(loop)))
            else
              push!(gespinds.args, mulexpr(getsym(first(loop)), stride))
            end
          elseif stride == 1
            push!(gespinds.args, addexpr(getsym(first(loop)), constoffset))
          else
            push!(gespinds.args, addexpr(mulexpr(getsym(first(loop)), stride), constoffset))
          end
        elseif isknown(first(loop))
          loopfirst = gethint(first(loop))*stride + constoffset
          if loopfirst == 0
            push!(gespinds.args, gespsymbol)
          else
            push!(gespinds.args, Expr(:call, GlobalRef(Base, :(+)), gespsymbol, staticexpr(loopfirst)))
          end
        else
          addedstarts = if stride == 1
            Expr(:call, GlobalRef(Base, :(+)), gespsymbol, getsym(first(loop)))
          else
            Expr(:call, GlobalRef(Base, :(+)), mulexpr(stride, gespsymbol), getsym(first(loop)))
          end
          if constoffset == 0
            push!(gespinds.args, addedstarts)
          else
            push!(gespinds.args, Expr(:call, GlobalRef(Base, :(+)), addedstarts, staticexpr(constoffset)))
          end
        end
      end
    end # if we hit the following elseif/else, not a loopindex
  elseif fromgsp # from gsp means that a loop could be a CartesianIndices, so we may need to expand
    #TODO: broadcast dimensions in case of cartesian indices
    rangesym = ind
    foundind = false
    for op ∈ operations(ls)
      if name(op) === ind
        loopdeps = loopdependencies(op)
        foundind = true
        if length(loopdeps) ≠ 0
          rangesym = getloop(ls, first(loopdeps)).rangesym
        else
          isconstantop(op) || throw(LoopError("Please file an issue with LoopVectorization.jl with a reproducer; tried to eliminate a non-constant operation."))
          rangesym = name(op)
        end
        break
      end
    end
    @assert foundind
    if rangesym === Symbol("") # there is no rangesym, must be statically sized.
      pushgespsym!(gespinds, gespsymbol, constoffset)
    else
      pushsimdims!(gespinds, rangesym, gespsymbol, constoffset)
    end
  else # it is known all inds are 1-dimensional
    pushgespsym!(gespinds, gespsymbol, constoffset)
  end
  return nothing
end
function pushgespsym!(gespinds::Expr, gespsymbol::Symbol, constoffset::Int)
  if gespsymbol === Symbol("")
    push!(gespinds.args, staticexpr(constoffset))
  elseif constoffset == 0
    push!(gespinds.args, gespsymbol)
  else
    push!(gespinds.args, addexpr(gespsymbol, constoffset))
  end
  return nothing
end
function pushsimdims!(gespinds::Expr, rangesym::Symbol, gespsymbol::Symbol, constoffset::Int)
  simdimscall = Expr(:call, lv(:similardims), rangesym)
  pushgespsym!(simdimscall, gespsymbol, constoffset)
  push!(gespinds.args, simdimscall)
  return nothing
end

# if using ptr offsets, need to offset by loopstart
# if calculating by index, do not offset by loopstart
"""
Returns a vector of length equal to the number of indices.
A value > 0 indicates which loop number that index corresponds to when incrementing the pointer.
A value < 0 indicates that abs(value) is the corresponding loop, and a `loopvalue` will be used.
"""
function use_loop_induct_var!(
  ls::LoopSet, q::Expr, ar::ArrayReferenceMeta, allarrayrefs::Vector{ArrayReferenceMeta}, allarrayrefsind::Int, includeinlet::Bool
  # array_refs_with_same_name::Vector{Int}, arrayref_to_name_op_collection::Vector{Vector{Tuple{Int,Int,Int}}}
)::Vector{Int}
  us = ls.unrollspecification
  li = ar.loopedindex
  looporder = reversenames(ls)
  uliv = Vector{Int}(undef, length(li))
  indices = getindices(ar)
  offset = first(indices) === DISCONTIGUOUS
  if length(indices) != offset + length(li)
    println(ar)
    throw("Length of indices and length of offset do not match!")
  end
  isbroadcast = ls.isbroadcast
  # no constant offsets when broadcasting
  # if (!isbroadcast) && !(all(li))
  #   checkforgespind = true
  #   gespindoffsets = cse_constant_offsets!(ls, q, ar, allarrayrefs, allarrayrefsind, array_refs_with_same_name, arrayref_to_name_op_collection)
  # else
  #   checkforgespind = false
  #   gespindoffsets = indices#dummy
  # end
  gespinds = Expr(:tuple)
  offsetprecalc_descript = Expr(:tuple)
  use_offsetprecalc = false
  vptrar = vptr(ar)
  Wisz = false#ls.vector_width == 0
  for (i,isli) ∈ enumerate(li)
    ii = i + offset
    ind = indices[ii]
    Wisz && push!(gespinds.args, staticexpr(0)) # wrong for `@_turbo`...
    if !li[i] # if it wasn't set
      uliv[i] = 0
      push!(offsetprecalc_descript.args, 0)
      Wisz || pushgespind!(gespinds, ls, Symbol(""), 0, 1, ind, isli, true, false)
    elseif ind === CONSTANTZEROINDEX
      uliv[i] = 0
      push!(offsetprecalc_descript.args, 0)
      Wisz || pushgespind!(gespinds, ls, Symbol(""), 0, 1, ind, isli, true, false)
    elseif isbroadcast ||
      ((isone(ii) && (last(looporder) === ind)) && !(otherindexunrolled(ls, ind, ar)) ||
      multiple_with_name(vptrar, allarrayrefs)) ||
      (iszero(ls.vector_width) && isstaticloop(getloop(ls, ind)))# ||
      # Not doing normal offset indexing
      uliv[i] = -findfirst(Base.Fix2(===,ind), looporder)::Int
      push!(offsetprecalc_descript.args, 0) # not doing offset indexing, so push 0
      Wisz || pushgespind!(gespinds, ls, Symbol(""), 0, 1, ind, isli, true, false)
    else
      uliv[i] = findfirst(Base.Fix2(===,ind), looporder)::Int
      loop = getloop(ls, ind)
      push!(offsetprecalc_descript.args, max(5,us.u₁+1,us.u₂+1))
      use_offsetprecalc = true
      Wisz || pushgespind!(gespinds, ls, Symbol(""), 0, 1, ind, isli, false, false)
    end
    # cases for pushgespind! and loopval!
    # if !isloopval, same as before
    # if isloopval & `uliv[i] > 0`, then must setup
  end
  # arrayref_to_name_op = arrayref_to_name_op_collection[allarrayrefsind]
  if includeinlet#first(arrayref_to_name_op)[2] == 1 # if it's the first inside `array_refs_with_same_name`, we need to add the expression to the let block
    vpgesped = Expr(:call, lv(:gesp), vptrar, gespinds)
    if use_offsetprecalc
      vpgesped = Expr(:call, lv(:offsetprecalc), vpgesped, Expr(:call, Expr(:curly, :Val, offsetprecalc_descript)))
    end
    push!(q.args, Expr(:(=), vptrar, vpgesped))
    push!(q.args, Expr(:(=), vptr_offset(vptrar), Expr(:call, GlobalRef(VectorizationBase, :increment_ptr), vptrar)))
  end
  uliv
end

# Plan here is that we increment every unique array
function add_loop_start_stop_manager!(ls::LoopSet)
    q = Expr(:block)
    us = ls.unrollspecification
    # Presence of an explicit use of a loopinducation var means we should use that, so we look for one
    # TODO: replace first with only once you add Compat as a dep or drop support for older Julia versions
    loopinductvars = Symbol[]
    for op ∈ operations(ls)
        isloopvalue(op) && push!(loopinductvars, first(loopdependencies(op)))
    end
    # Filtered ArrayReferenceMetas, we must increment each
  # arrayrefs, name_to_array_map, unique_to_name_and_op_map = uniquearrayrefs(ls)
    arrayrefs, includeinlet = uniquearrayrefs(ls)
    use_livs = Vector{Vector{Int}}(undef, length(arrayrefs))
    # for i ∈ eachindex(name_to_array_map)
    for i ∈ eachindex(arrayrefs)
      use_livs[i] = use_loop_induct_var!(ls, q, arrayrefs[i], arrayrefs, i, includeinlet[i])
      #name_to_array_map[first(first(unique_to_name_and_op_map[i]))], unique_to_name_and_op_map)
    end
    # loops, sorted from outer-most to inner-most
    looporder = reversenames(ls)
    # For each loop, we need to choose an induction variable
    nloops = length(looporder)
    # loopstarts = Vector{Vector{ArrayReferenceMeta}}(undef, nloops)
    loopstarts = fill(ArrayReferenceMeta[], nloops)
    terminators = Vector{Int}(undef, nloops) # zero is standard loop induct var
    # loopincrements = Vector{Vector{ArrayReferenceMeta}}(undef, nloops) # Not needed; copy of loopstarts
    # The purpose of the reminds thing here is to pick which of these to use for the terminator
    # We want to pick the one with the fewest outer loops with respect to this one, to minimize
    # the number of redefinitions of max-pointer used for the termination comparison.
    reached_indices = zeros(Int, length(arrayrefs))
    for (i,loopsym) ∈ enumerate(looporder) # iterates from outer to inner
        loopstartᵢ = ArrayReferenceMeta[]; arⱼ = 0; minrem = typemax(Int);
        ric = Tuple{Int,Int}[]
        for j ∈ eachindex(use_livs) # j is array ref number
            for (l,k) ∈ enumerate(use_livs[j])# l is index number, k is loop number
                if k == i
                    push!(loopstartᵢ, arrayrefs[j])
                    push!(ric, ((reached_indices[j] += 1), length(loopstartᵢ)))
                end
            end
        end
        loopstarts[nloops+1-i] = loopstartᵢ
        terminators[nloops+1-i] = if (loopsym ∈ loopinductvars) || (any(r -> any(isequal(-i), r), use_livs)) || iszero(length(loopstartᵢ))
            0
        else
            # @assert !iszero(length(loopstartᵢ))
            last(ric[argmin(first.(ric))]) # index corresponds to array ref's position in loopstart
        end
    end
    ls.lssm = LoopStartStopManager(
        terminators, loopstarts, arrayrefs
    )
    q
end
maxsym(ptr, sub) = Symbol(ptr, "##MAX##", sub, "##")
function pointermax(ls::LoopSet, ar::ArrayReferenceMeta, n::Int, sub::Int, isvectorized::Bool)::Expr
    pointermax(ls, ar, n, sub, isvectorized, getloop(ls, names(ls)[n]))
end
function pointermax(ls::LoopSet, ar::ArrayReferenceMeta, n::Int, sub::Int, isvectorized::Bool, loop::Loop)::Expr
    start = first(loop)
    stop = last(loop)
    incr = step(loop)
    if isknown(start) & isknown(stop)
        return pointermax(ls, ar, n, sub, isvectorized, 1 + gethint(stop) - gethint(start), incr)
    end
    looplensym = isone(start) ? getsym(stop) : loop.lensym
    pointermax(ls, ar, n, sub, isvectorized, looplensym, incr)
end
function pointermax_index(
    ls::LoopSet, ar::ArrayReferenceMeta, n::Int, sub::Int, isvectorized::Bool, stophint::Int, incr::MaybeKnown
)::Tuple{Expr,Int}
    # @unpack u₁loopnum, u₂loopnum, vloopnum, u₁, u₂ = us
    loopsym = names(ls)[n]
    index = Expr(:tuple)
    found_loop_sym = false
    ind = 0
    for (j,i) ∈ enumerate(getindicesonly(ar))
        if i === loopsym
            ind = j
            if iszero(sub)
                push!(index.args, stophint)
            else
                _ind = if isvectorized
                    if isone(sub)
                        Expr(:call, lv(:vsub_nsw), staticexpr(stophint), VECTORWIDTHSYMBOL)
                    else
                        Expr(:call, lv(:vsub_nsw), staticexpr(stophint), mulexpr(VECTORWIDTHSYMBOL, sub))
                    end
                else
                    staticexpr(stophint - sub)
                end
                stride = getstrides(ar)[j]
                if isknown(incr)
                    stride *= gethint(incr)
                else
                    _ind = mulexpr(_ind, getsym(incr))
                end
                if stride ≠ 1
                    @assert stride ≠ 0 "stride shouldn't be 0 if used for determining loop start/stop, but loop $n array $ar was."
                    _ind = lazymulexpr(stride, _ind)
                end
                push!(index.args, _ind)
            end
        else
            push!(index.args, Expr(:call, lv(:Zero)))
        end
    end
    @assert ind ≠ 0 "Failed to find $loopsym"
    index, ind
end
function pointermax_index(ls::LoopSet, ar::ArrayReferenceMeta, n::Int, sub::Int, isvectorized::Bool, stopsym, incr::MaybeKnown)::Tuple{Expr,Int}
    loopsym = names(ls)[n]
    index = Expr(:tuple);
    ind = 0
    for (j,i) ∈ enumerate(getindicesonly(ar))
        if i === loopsym
            ind = j
            if iszero(sub)
                push!(index.args, stopsym)
            else
                _ind = if isvectorized
                    if isone(sub)
                        Expr(:call, lv(:vsub_nsw), stopsym, VECTORWIDTHSYMBOL)
                    else
                        Expr(:call, lv(:vsub_nsw), stopsym, mulexpr(VECTORWIDTHSYMBOL, sub))
                    end
                else
                     subexpr(stopsym, sub)
                end
                stride = getstrides(ar)[j]
                if isknown(incr)
                    stride *= gethint(incr)
                else
                    _ind = mulexpr(_ind, getsym(incr))
                end
                if stride ≠ 1
                    @assert stride ≠ 0 "stride shouldn't be 0 if used for determining loop start/stop, but loop $n array $ar was."
                    _ind = lazymulexpr(stride, _ind)
                end
                push!(index.args, _ind)
            end
        else
            push!(index.args, Expr(:call, lv(:Zero)))
        end
    end
    @assert ind != 0 "Failed to find $loopsym"
    index, ind
end
function pointermax(ls::LoopSet, ar::ArrayReferenceMeta, n::Int, sub::Int, isvectorized::Bool, stopsym, incr::MaybeKnown)::Expr
  index = first(pointermax_index(ls, ar, n, sub, isvectorized, stopsym, incr))
  vptrar = vptr(ar)
  Expr(:call, GlobalRef(VectorizationBase,:increment_ptr), vptrar, vptr_offset(vptrar), index)
end

function defpointermax(ls::LoopSet, ar::ArrayReferenceMeta, n::Int, sub::Int, isvectorized::Bool)::Expr
    Expr(:(=), maxsym(vptr(ar), sub), pointermax(ls, ar, n, sub, isvectorized))
end
function offsetindex(dim::Int, ind::Int, scale::Int, isvectorized::Bool, incr::MaybeKnown)
    index = Expr(:tuple)
    for d ∈ 1:dim
        if d ≠ ind || iszero(scale)
            push!(index.args, Expr(:call, lv(:Zero)))
            continue
        end
        if !isvectorized
            pushmulexpr!(index, scale, incr)
        elseif isone(scale)
            pushmulexpr!(index, VECTORWIDTHSYMBOL, incr)
        else
            push!(index.args, mulexpr(VECTORWIDTHSYMBOL, scale, incr))
        end
    end
    index
end
function append_pointer_maxes!(
    loopstart::Expr, ls::LoopSet, ar::ArrayReferenceMeta, n::Int, submax::Int, isvectorized::Bool, stopindicator, incr::MaybeKnown
)
  vptr_ar = vptr(ar)
  if submax < 2
    for sub ∈ 0:submax
      push!(loopstart.args, Expr(:(=), maxsym(vptr_ar, sub), pointermax(ls, ar, n, sub, isvectorized, stopindicator, incr)))
    end
  else
    index, ind = pointermax_index(ls, ar, n, submax, isvectorized, stopindicator, incr)
    pointercompbase = maxsym(vptr_ar, submax)
    ip = GlobalRef(VectorizationBase, :increment_ptr)
    push!(loopstart.args, Expr(:(=), pointercompbase, Expr(:call, ip, vptr_ar, vptr_offset(vptr_ar), index)))
    dim = length(getindicesonly(ar))
    # OFFSETPRECALCDEF = true
    # if OFFSETPRECALCDEF
    strd = getstrides(ar)[ind]
    for sub ∈ 0:submax-1
      ptrcmp = Expr(:call, ip, vptr_ar, pointercompbase, offsetindex(dim, ind, (submax - sub)*strd, isvectorized, incr))
      push!(loopstart.args, Expr(:(=), maxsym(vptr_ar, sub), ptrcmp))
    end
  end
end
function append_pointer_maxes!(loopstart::Expr, ls::LoopSet, ar::ArrayReferenceMeta, n::Int, submax::Int, isvectorized::Bool)
  loop = getloop(ls, n)
  @assert loop.itersymbol == names(ls)[n]
  start = first(loop)
  stop = last(loop)
  incr = step(loop)
  if isknown(start) & isknown(stop)
    return append_pointer_maxes!(loopstart, ls, ar, n, submax, isvectorized, startstopΔ(loop)+1, incr)
  end
  looplensym = isone(start) ? getsym(stop) : loop.lensym
  append_pointer_maxes!(loopstart, ls, ar, n, submax, isvectorized, looplensym, incr)
end

function maxunroll(us::UnrollSpecification, n)
  @unpack u₁loopnum, u₂loopnum, u₁, u₂ = us
  if n == u₁loopnum# && u₁ > 1
    u₁
  elseif n == u₂loopnum# && u₂ > 1
    u₂
  else
    1
  end
end

function startloop(ls::LoopSet, us::UnrollSpecification, n::Int, staticinit::Bool = false)
  @unpack u₁loopnum, u₂loopnum, vloopnum, u₁, u₂ = us
  lssm = ls.lssm
  termind = lssm.terminators[n]
  ptrdefs = lssm.incrementedptrs[n]
  loopstart = Expr(:block)
  firstloop = n == num_loops(ls)
  for ar ∈ ptrdefs
    ptr_offset = vptr_offset(ar)
    push!(loopstart.args, Expr(:(=), ptr_offset, ptr_offset))
  end
  if iszero(termind)
    loopsym = names(ls)[n]
    push!(loopstart.args, startloop(getloop(ls, loopsym), loopsym, staticinit))
  else
    isvectorized = n == vloopnum
    submax = maxunroll(us, n)
    append_pointer_maxes!(loopstart, ls, ptrdefs[termind], n, submax, isvectorized)
  end
  loopstart
end
function offset_ptr(
    ar::ArrayReferenceMeta, us::UnrollSpecification, loopsym::Symbol, n::Int, UF::Int, offsetinds::Vector{Bool}, loop::Loop
)
  indices = getindices(ar)
  strides = getstrides(ar)
  offset = first(indices) === DISCONTIGUOUS
  gespinds = Expr(:tuple)
  li = ar.loopedindex
  for i ∈ eachindex(li)
    ii = i + offset
    ind = indices[ii]
    if !offsetinds[i] || ind !== loopsym
      push!(gespinds.args, Expr(:call, lv(:Zero)))
    else
      incrementloopcounter!(gespinds, us, n, UF * strides[i], loop)
    end
    # ind == loopsym && break
  end
  vpoff = vptr_offset(ar)
  call = Expr(:call, GlobalRef(VectorizationBase, :increment_ptr), vptr(ar), vpoff, gespinds)
  Expr(:(=), vpoff, call)
end
function incrementloopcounter!(q::Expr, ls::LoopSet, us::UnrollSpecification, n::Int, UF::Int)
    @unpack u₁loopnum, u₂loopnum, vloopnum, u₁, u₂ = us
    lssm = ls.lssm
    ptrdefs = lssm.incrementedptrs[n]
    looporder = names(ls)
    loopsym = looporder[n]
    termind = lssm.terminators[n]
    loop = getloop(ls, n)
    if iszero(termind) # increment liv
        push!(q.args, incrementloopcounter(us, n, loopsym, UF, loop))
    end
    for (j,ar) ∈ enumerate(ptrdefs)
        offsetinds = indices_calculated_by_pointer_offsets(ls, ar)
        push!(q.args, offset_ptr(ar, us, loopsym, n, UF, offsetinds, loop))
    end
    nothing
end
function terminatecondition(ls::LoopSet, us::UnrollSpecification, n::Int, inclmask::Bool, UF::Int)
  lssm = ls.lssm
  termind = lssm.terminators[n]
  if iszero(termind)
    loop = getloop(ls, n)
    return terminatecondition(loop, us, n, loop.itersymbol, inclmask, UF)
  end

  termar = lssm.incrementedptrs[n][termind]
  ptr = vptr(termar)
  optr = vptr_offset(ptr)
  if inclmask && isvectorized(us, n)
    Expr(:call, GlobalRef(VectorizationBase, :vlt), optr, maxsym(ptr, 0), ptr)
  else
    Expr(:call, GlobalRef(VectorizationBase, :vle), optr, maxsym(ptr, UF), ptr)
  end
end


text:
module EXT_TEST

export EXT_TEST_hmmalign
export EXT_TEST_hmmsearch
export EXT_TEST_cmsearch

######################################################################
# FUNCTION:              EXT_TEST_hmmalign
######################################################################
function EXT_TEST_hmmalign()
	try
	    run(pipeline(`hmmalign -h`, DevNull))
	catch
	    return false
	end
    return true
end

######################################################################
# FUNCTION:              EXT_TEST_hmmsearch
######################################################################
function EXT_TEST_hmmsearch()
	try
	    run(pipeline(`hmmsearch -h`, DevNull))
	catch
	    return false
	end
    return true
end

######################################################################
# FUNCTION:              EXT_TEST_cmsearch
######################################################################
function EXT_TEST_cmsearch()
	try
	    run(pipeline(`cmsearch -h`, DevNull))
	catch
	    return false
	end
    return true
end

#/module
end


text:
const ScalarLinearSet{T} = Union{MOI.EqualTo{T}, MOI.LessThan{T}, MOI.GreaterThan{T}}

vector_set(::MOI.EqualTo) = MOI.Zeros(1)
vector_set(::MOI.LessThan) = MOI.Nonpositives(1)
vector_set(::MOI.GreaterThan) = MOI.Nonnegatives(1)

vector_set_type(::Type{<:MOI.EqualTo}) = MOI.Zeros
vector_set_type(::Type{<:MOI.LessThan}) = MOI.Nonpositives
vector_set_type(::Type{<:MOI.GreaterThan}) = MOI.Nonnegatives

"""
    VectorizeBridge{T, F, S, G}

Transforms a constraint `G`-in-`scalar_set_type(S, T)` where
`S <: VectorLinearSet` to `F`-in-`S`.
"""
mutable struct VectorizeBridge{T, F, S, G} <: AbstractBridge
    vector_constraint::CI{F, S}
    set_constant::T # constant in scalar set
end
function bridge_constraint(::Type{VectorizeBridge{T, F, S, G}},
                           model::MOI.ModelLike, g::G,
                           set::MOI.AbstractScalarSet) where {T, F, S, G}
    set_constant = MOI.constant(set)
    h = MOIU.operate(-, T, g, set_constant)
    if -set_constant != MOI.constant(h)[1]
        # This means the constant in `f` was not zero
        constant = MOI.constant(h)[1] + set_constant
        throw(MOI.ScalarFunctionConstantNotZero{typeof(constant), G,
                                                typeof(set)}(constant))
    end
    f = MOIU.operate(vcat, T, h)
    vector_constraint = MOI.add_constraint(model, f, vector_set(set))
    return VectorizeBridge{T, F, S, G}(vector_constraint, set_constant)
end

function MOI.supports_constraint(::Type{VectorizeBridge{T}},
                                ::Type{<:MOI.AbstractScalarFunction},
                                ::Type{<:ScalarLinearSet{T}}) where T
    return true
end
function added_constraint_types(::Type{<:VectorizeBridge{T, F, S}}) where {T, F, S}
    return [(F, S)]
end
function concrete_bridge_type(::Type{<:VectorizeBridge{T}},
                              G::Type{<:MOI.AbstractScalarFunction},
                              S::Type{<:ScalarLinearSet{T}}) where T
    H = MOIU.promote_operation(-, T, G, T)
    F = MOIU.promote_operation(vcat, T, H)
    return VectorizeBridge{T, F, vector_set_type(S), G}
end

# Attributes, Bridge acting as an model
function MOI.get(::VectorizeBridge{T, F, S},
                 ::MOI.NumberOfConstraints{F, S}) where {T, F, S}
    return 1
end
function MOI.get(bridge::VectorizeBridge{T, F, S},
                 ::MOI.ListOfConstraintIndices{F, S}) where {T, F, S}
    return [bridge.vector_constraint]
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::VectorizeBridge)
    MOI.delete(model, bridge.vector_constraint)
end

# Attributes, Bridge acting as a constraint

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintPrimal,
                 bridge::VectorizeBridge)
    x = MOI.get(model, attr, bridge.vector_constraint)
    @assert length(x) == 1
    y = x[1]
    status = MOI.get(model, MOI.PrimalStatus(attr.N))
    if status != MOI.INFEASIBILITY_CERTIFICATE &&
       status != MOI.NEARLY_INFEASIBILITY_CERTIFICATE
       # If it is an infeasibility certificate, it is a ray and satisfies the
       # homogenized problem, see https://github.com/JuliaOpt/MathOptInterface.jl/issues/433
       # Otherwise, we need to add the set constant since the ConstraintPrimal
       # is defined as the value of the function and the set_constant was
       # removed from the original function
       y += bridge.set_constant
   end
   return y
end
function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintDual,
                 bridge::VectorizeBridge)
    x = MOI.get(model, attr, bridge.vector_constraint)
    @assert length(x) == 1
    return x[1]
end
function MOI.modify(model::MOI.ModelLike, bridge::VectorizeBridge,
                    change::MOI.ScalarCoefficientChange)
    MOI.modify(model, bridge.vector_constraint,
               MOI.MultirowChange(change.variable,
                                  [(1, change.new_coefficient)]))
end
function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintSet,
                 bridge::VectorizeBridge, new_set::ScalarLinearSet)
    bridge.set_constant = MOI.constant(new_set)
    MOI.modify(model, bridge.vector_constraint,
               MOI.VectorConstantChange([-bridge.set_constant]))
end
function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintFunction,
                 bridge::VectorizeBridge{T, F, S, G}) where {T, F, S, G}
    f = MOIU.scalarize(MOI.get(model, attr, bridge.vector_constraint), true)
    @assert isone(length(f))
    # If `G` is `MOI.SingleVariable`, `f` will be `MOI.ScalarAffineFunction`.
    return convert(G, f[1])
end
function MOI.get(model::MOI.ModelLike, ::MOI.ConstraintSet,
                 bridge::VectorizeBridge{T, F, S}) where {T, F, S}
    return scalar_set_type(S, T)(bridge.set_constant)
end


text:
@testset "_parsefeterm" begin
    @test _parsefeterm(term(:male)) === nothing
    @test _parsefeterm(fe(:hhidpn)) == ([:hhidpn]=>Symbol[])
    @test _parsefeterm(fe(:hhidpn)&fe(:wave)&term(:male)) == ([:hhidpn, :wave]=>[:male])
end

@testset "getfename" begin
    @test getfename([:hhidpn]=>Symbol[]) == "fe_hhidpn"
    @test getfename([:hhidpn, :wave]=>[:male]) == "fe_hhidpn&fe_wave&male"
end


text:
using KernelDensityEstimate
using Base.Test
# include("BallTree01.jl")
# include("BallTreeDensity01.jl")
# include("KDE01.jl")
# include("DualTree01.jl")

# type BallTree
#   dims::Int64                     # dimension of data
#   num_points::Int64               # of points
#   centers::Array{Float64,1}       # ball centers, dims numbers per ball
#   ranges::Array{Float64,1}        # bounding box ranges, dims per ball, dist from center to one side
#   weights::Array{Float64,1}       # total weight in each ball
#
#   left_child::Array{Int64,1}
#   right_child::Array{Int64,1}     # left, right children; no parent indices
#   lowest_leaf::Array{Int64,1}
#   highest_leaf::Array{Int64,1}    # lower & upper leaf indices for each ball
#   permutation::Array{Int64,1}     # point's position in the original data
#
#   next::Int64                     # internal var for placing the non-leaf nodes
#
#   swapHandle::Function
#   calcStatsHandle::Function
#   data
# end
# type BallTreeDensity
#   bt::BallTree
#
#   KernelType
#   multibandwidth::Int64               # flag: is bandwidth uniform?
#
#   means::Array{Float64,1}                  # Weighted mean of points from this level down
#   bandwidth::Array{Float64,1}              # Variance or other multiscale bandwidth
#   bandwidthMin::Array{Float64,1}
#   bandwidthMax::Array{Float64,1}           # Bounds on BW in non-uniform case
#
#   calcStatsHandle::Function
#   swapHandle::Function
# end

# parse the output from matlab process
function parseMatPrintKDE(filename::ASCIIString)
  fid = open(filename,"r")
  dict = Dict{ASCIIString, Array{Float64,1}}()
  for line in eachline(fid)
    twos = split(split(line,"\n")[1],"=")
    s = (split(split(twos[2],"[")[2],"]")[1])
    dict[twos[1]] = vec(readdlm(IOBuffer(s),','))
  end
  close(fid)
  return dict
end

function constructBTD(dict::Dict{ASCIIString, Array{Float64,1}})
  refbt = BallTree(dict["dims"][1],dict["num_points"][1],
  dict["centers"], #
  dict["ranges"], # =
  dict["weights"], #
  round(Int64,dict["left_child"]), # =
  round(Int64,dict["right_child"]), # =
  round(Int64,dict["lowest_leaf"]), # =
  round(Int64,dict["highest_leaf"]), # =
  round(Int64,dict["permutation"]), # =
  0, +, +, Union{} )

  refbtd = BallTreeDensity(refbt, Union{}, 0,
  dict["means"], # =
  dict["bandwidth"], # =
  dict["bwMin"],
  dict["bwMax"],
  +, +)
  return refbtd
end

function testSubtract(b1::BallTreeDensity, b2::BallTreeDensity, tol=1e-10)
    #@show b1.bt.dims, b2.bt.dims, b1.bt.num_points, b2.bt.num_points
    if (b1.bt.dims - b2.bt.dims) != 0 || (b1.bt.num_points - b2.bt.num_points) != 0
        error("testSubtract -- b1 and b2 are not the same size.")
    end
    function testDiffs(a1::Array{Float64},a2::Array{Float64}; tol=1e-10, name="")
        n = norm(a1 - a2)
        if n > tol
          @show name, n
          error("testSubtract -- b1 and b2 are not the same.")
        end
        nothing
    end
    function testInds(a1::Array{Int64},a2::Array{Int64})
        for i in 1:length(a1)
            res = 0
            if (a1[i] < 0)
                res = (a2[i]+1)
            else
                res = (a1[i]+1-a2[i])
            end
            if res != 0
              error("testSubtract -- b1 and b2 have different indices.")
            end
        end
        nothing
    end
    testDiffs(b1.bt.centers, b2.bt.centers, name="centers", tol=tol)
    testDiffs(b1.bt.ranges, b2.bt.ranges, name="ranges", tol=tol)
    testDiffs(b1.bt.weights, b2.bt.weights, name="weights", tol=tol)
    testInds(b1.bt.left_child, b2.bt.left_child)
    testInds(b1.bt.right_child, b2.bt.right_child)
    testInds(b1.bt.lowest_leaf, b2.bt.lowest_leaf)
    testInds(b1.bt.highest_leaf, b2.bt.highest_leaf)
    testInds(b1.bt.permutation[(b1.bt.num_points+1):end], b2.bt.permutation[(b2.bt.num_points+1):end])
    testDiffs(b1.means, b2.means, name="means", tol=tol)
    testDiffs(b1.bandwidth, b2.bandwidth, tol=tol)
    testDiffs(b1.bandwidthMin, b2.bandwidthMin, tol=tol)
    testDiffs(b1.bandwidthMax, b2.bandwidthMax, tol=tol)
    println("Success")
    true
end

function UnitTest1D01()
  print("Running UnitTest1D01...")
  p = kde!([.1,.45,.55,3.8],[0.08])

  d = parseMatPrintKDE("testdata/test1DResult.txt")
  refbtd = constructBTD(d)

  # printBallTree(p)
  testSubtract(refbtd, p, 1e-5)
end

function UnitTest1Dlcv01()
  print("Running UnitTest1Dlcv01...")
  x = vec(readdlm("testdata/test1Dlcv100.txt")')
  p = kde!(x,"lcv")

  d = parseMatPrintKDE("testdata/test1Dlcv100Result.txt")
  refbtd = constructBTD(d)

  # printBallTree(p)
  testSubtract(refbtd, p, 1e-4)
end

function UnitTest2D01()
  print("Running UnitTest2D01...")
  pts = [[0.5172, 0.7169, 0.4049]';
         [0.0312, 1.0094, 2.0204]']
  p = kde!(pts,[0.1])

  d = parseMatPrintKDE("testdata/test2DResult.txt")
  refbtd = constructBTD(d)

  # printBallTree(p)
  testSubtract(refbtd, p, 1e-5)
end

function UnitTest2Dlcv01()
  print("Running UnitTest2Dlcv01...")
  x = readdlm("testdata/test2Dlcv100.txt")'
  p = kde!(x,"lcv")

  d = parseMatPrintKDE("testdata/test2Dlcv100Result.txt")
  refbtd = constructBTD(d)

  # printBallTree(p)
  testSubtract(refbtd, p, 1e-4)
end

function UnitTest2Dvar01()
  print("Running UnitTest2Dvar01...")
  pts = [[0.5172, 7.169, 4.049]';
         [0.0312, 10.0094, -2.0204]']
  p=kde!(pts,[0.1; 1.0]);
  d = parseMatPrintKDE("testdata/test2DvarResult.txt")
  refbtd = constructBTD(d)

  # printBallTree(p)
  testSubtract(refbtd, p, 1e-4)
end

function UnitTest2Dvarlcv01()
  print("Running UnitTest2Dvarlcv01...")
  x = readdlm("testdata/test2Dvarlcv100.txt")'
  p = kde!(x,"lcv")

  d = parseMatPrintKDE("testdata/test2Dvarlcv100Result.txt")
  refbtd = constructBTD(d)

  # printBallTree(p)
  testSubtract(refbtd, p, 2e-3)
end

function testProds(;D=3,M=6,N=100,n=100, dev=1.0, MCMC=5)
  P = BallTreeDensity[];
  [push!(P, kde!(dev*randn(D,N))) for i in 1:M];
  dummy = kde!(randn(D,n),[1.0]);
  pGM, = prodAppxMSGibbsS(dummy, P, Union{}, Union{}, MCMC);
  sum(abs(pGM))<1e-14 ? error("testProds -- prodAppxMSGibbsS, nothing in pGM, len $(length(P))") : nothing
  prodDev = sqrt(dev^(2*M)/(M*(dev^2)))
  T1 = norm(Base.mean(pGM,2)) < 1.0*prodDev
  T2 = true
  for i in 1:D
    tv = Base.std(pGM[i,:])
    lv = (0.66*prodDev < tv < 1.33*prodDev)
    T2 = T2 && lv
  end
  T1 && T2
end

function rangeTestProds(;D=3,M=6,N=100,n=100, dev=1.0, MCMC=5)
  @show v = [testProds(D=D,M=M,N=N,n=n, dev=dev, MCMC=MCMC) for i in 1:10]
  sum(map(Int,v)) >= 5
end

function rangeUnitTests()
  passt = true
  @show passt = passt && rangeTestProds(D=2,M=2);
  @show passt = passt && rangeTestProds(D=2,M=4);
  @show passt = passt && rangeTestProds(D=2,M=6);
  @show passt = passt && rangeTestProds(D=3,M=6, MCMC=10);
  @show passt = passt && rangeTestProds(D=4,M=6, MCMC=10);
  # @show passt = passt && rangeTestProds(D=3,M=10, MCMC=10);
  @show passt = passt && rangeTestProds(D=3,M=5,N=300);
  @show passt = passt && rangeTestProds(D=2,M=7,n=300);
  @show passt = passt && rangeTestProds(D=3,M=2,MCMC=100);
  return passt
end

function intgAppxGaussianOffs(;offs::Float64=0.0, N::Int=201, dim::Int=1)
  p = kde!(randn(dim,100));
  pts = randn(dim,150)
  pts[1,:] += offs
  q = kde!(pts);
  return intersIntgAppxIS(p,q, N=N)
end

function integralAppxUnitTests()
  testflag = true
  @show a = intgAppxGaussianOffs(offs=0.0,dim=1)
  @show testflag = testflag && 0.2 < a < 0.35
  @show a = intgAppxGaussianOffs(offs=1.0,dim=1,N=1000)
  @show testflag = testflag && 0.1 < a < 0.3
  @show a = intgAppxGaussianOffs(offs=-2.0,dim=1,N=1000)
  @show testflag = testflag && 0.01 < a < 0.15
  @show a = intgAppxGaussianOffs(offs=0.0,dim=2)
  @show testflag = testflag && 0.05 < a < 0.12

  return testflag
end

function testRand()
  println("testing rand functionality")
  p = kde!(rand(2,100))
  pts = rand(p,100)
  return true
end

@test UnitTest1D01()
@test UnitTest1Dlcv01()
@test UnitTest2D01()
# UnitTest2Dlcv01()
@test UnitTest2Dvar01()
#UnitTest2Dvarlcv01()
@test rangeUnitTests()
@test integralAppxUnitTests()

@test testRand()

  # check for errors on plotting code
p = kde!(rand(100));
plotKDE([p],c=["red"]);
@test true


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

# Collection of sources required to build DSFMTBuilder
sources = [
    "http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/dSFMT-src-2.2.3.tar.gz" =>
    "82344874522f363bf93c960044b0a6b87b651c9565b6312cf8719bb8e4c26a0e",

    # dSFMT patches to create a shared library
    "patches",
]

# Bash recipe for building across all platforms
script = raw"""
cd dSFMT-src-2.2.3

# Apply all our patches
for f in $WORKSPACE/srcdir/*.patch; do
    patch -p1 < ${f}
done

DSFMT_CFLAGS="-DNDEBUG -DDSFMT_MEXP=19937 -fPIC -DDSFMT_DO_NOT_USE_OLD_NAMES  -O3 -finline-functions -fomit-frame-pointer -fno-strict-aliasing --param max-inline-insns-single=1800 -Wmissing-prototypes -Wall  -std=c99 -shared"

if [[ $target == x86_64-*  ]]; then
DSFMT_CFLAGS="$DSFMT_CFLAGS -msse2 -DHAVE_SSE2"
fi

mkdir -p $prefix/lib
$CC $DSFMT_CFLAGS dSFMT.c -o $prefix/lib/libdSFMT.$dlext

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    BinaryProvider.Linux(:i686, :glibc, :blank_abi),
    BinaryProvider.Linux(:x86_64, :glibc, :blank_abi),
    BinaryProvider.Windows(:x86_64, :blank_libc, :blank_abi),
    BinaryProvider.Windows(:i686, :blank_libc, :blank_abi),
    BinaryProvider.MacOS(:x86_64, :blank_libc, :blank_abi),
    BinaryProvider.Linux(:aarch64, :glibc, :blank_abi),
    BinaryProvider.Linux(:armv7l, :glibc, :eabihf),
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libdSFMT", :libdSFMT)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "DSFMTBuilder", sources, script, platforms, products, dependencies)


text:
abstract Sensor

type Electrode <: Sensor
    label::AbstractString
    coordinate::Coordinate
    info::Dict
end

import Base.show
function show{S <: Sensor}(s::S)
    println("Sensor: $(s.label) $(typeof(s)) - ($(s.coordinate.x), $(s.coordinate.y), $(s.coordinate.z)) ($(typeof(s.coordinate)))")
end

function show{S <: Sensor}(s::Array{S})
    println("$(length(s)) sensors: $(typeof(s[1])) ($(typeof(s[1].coordinate)))")
end

label{S <: Sensor}(s::S) = s.label
label{S <: Sensor}(s::Array{S, 1}) = AbstractString[si.label for si in s]
labels{S <: Sensor}(s::S) = label(s)
labels{S <: Sensor}(s::Array{S}) = label(s)

x{S <: Sensor}(s::S) = s.coordinate.x
y{S <: Sensor}(s::S) = s.coordinate.y
z{S <: Sensor}(s::S) = s.coordinate.z
x{S <: Sensor}(s::Array{S}) = AbstractFloat[si.coordinate.x for si in s]
y{S <: Sensor}(s::Array{S}) = AbstractFloat[si.coordinate.y for si in s]
z{S <: Sensor}(s::Array{S}) = AbstractFloat[si.coordinate.z for si in s]


text:
using BenchmarkTools
using Configurations

struct Asset
    package::Union{Nothing,String}
    path::String
end

function Asset(s::String)
    parts = strip.(split(s, ":"))
    if length(parts) == 1
        Asset(nothing, parts[1])
    elseif length(parts) == 2
        Asset(parts[1], parts[2])
    else
        error("invalid syntax for asset string: $s")
    end
end

macro asset_str(s::String)
    return Asset(s)
end

function Base.show(io::IO, x::Asset)
    print(io, "asset\"")
    if x.package !== nothing
        printstyled(io, x.package, ": "; color=:green)
    end
    printstyled(io, x.path, "\""; color=:cyan)
end

Base.convert(::Type{Asset}, s::String) = Asset(s)

@option struct Install
    path::String = "~/.julia"
    completion::Bool = true
    quiet::Bool = false
    compile::String = "yes"
    optimize::Int = 2
end

@option struct Precompile
    execution_file::Vector{String} = String[]
    statements_file::Vector{String} = String[]
end

@option struct SysImg
    path::String = "deps"
    incremental::Bool = true
    filter_stdlibs::Bool = false
    cpu_target::String = "pentium4;sandybridge,-xsaveopt,clone_all"
    precompile::Precompile = Precompile()
end

@option struct Download
    host::String = "github.com"
    user::String
    repo::String
end

@option struct Application
    path::String = "build"
    assets::Vector{Asset} = Asset[]
    incremental::Bool = false
    filter_stdlibs::Bool = true
    cpu_target::String = "pentium4;sandybridge,-xsaveopt,clone_all"
    precompile::Precompile = Precompile()
    c_driver_program::Union{String,Nothing} = nothing
end

@option struct Comonicon
    name::String

    install::Install = Install()
    sysimg::Maybe{SysImg} = nothing
    download::Union{Download,Nothing} = nothing
    application::Union{Application,Nothing} = nothing
end

option = Comonicon("foo",
    Install("~/.julia", true, false, "min", 2),
    SysImg(
        "deps", true, false, "native",
        Precompile(["deps/precopmile.jl"], String[])
    ),
    Download("github.com", "Roger-luo", "Foo.jl"),
    Application(
        "build",
        Asset[asset"PkgTemplate: templates", asset"assets/images"],
        true, false,
        "generic;sandybridge,-xsaveopt,clone_all;haswell,-rdrnd,base(1)",
        Precompile(String[], String[]),
        nothing
    )
)

d = to_dict(option)
@btime from_dict(Comonicon, d)

using Configurations: from_dict_dynamic, from_dict_specialize, from_dict_generated

@btime from_dict_dynamic(Comonicon, d)