text:
function ShowDemoWindowLayout()
    CImGui.CollapsingHeader("Layout") || return

    if CImGui.TreeNode("Child windows")
        ShowHelpMarker("Use child windows to begin into a self-contained independent scrolling/clipping regions within a host window.")
        @cstatic disable_mouse_wheel=false disable_menu=false line=Cint(50) begin
            @c CImGui.Checkbox("Disable Mouse Wheel", &disable_mouse_wheel)
            @c CImGui.Checkbox("Disable Menu", &disable_menu)

            goto_line = CImGui.Button("Goto")
            CImGui.SameLine()
            CImGui.PushItemWidth(100)
            goto_line |= @c CImGui.InputInt("##Line", &line, 0, 0, CImGui.ImGuiInputTextFlags_EnterReturnsTrue)
            CImGui.PopItemWidth()

            # child 1: no border, enable horizontal scrollbar
            begin
                window_flags = CImGui.ImGuiWindowFlags_HorizontalScrollbar | (disable_mouse_wheel ? CImGui.ImGuiWindowFlags_NoScrollWithMouse : 0)
                CImGui.BeginChild("Child1", ImVec2(CImGui.GetWindowContentRegionWidth() * 0.5, 260), false, window_flags)
                for i = 0:100-1
                    CImGui.Text(@sprintf("%04d: scrollable region", i))
                    (goto_line && line == i) && CImGui.SetScrollHereY()
                end
                (goto_line && line >= 100) && CImGui.SetScrollHereY()
                CImGui.EndChild()
            end

            CImGui.SameLine()

            # child 2: rounded border
            begin
                window_flags = (disable_mouse_wheel ? CImGui.ImGuiWindowFlags_NoScrollWithMouse : 0) | (disable_menu ? 0 : CImGui.ImGuiWindowFlags_MenuBar)
                CImGui.PushStyleVar(CImGui.ImGuiStyleVar_ChildRounding, 5.0)
                CImGui.BeginChild("Child2", ImVec2(0, 260), true, window_flags)
                if !disable_menu && CImGui.BeginMenuBar()
                    if CImGui.BeginMenu("Menu")
                        ShowExampleMenuFile()
                        CImGui.EndMenu()
                    end
                    CImGui.EndMenuBar()
                end
                CImGui.Columns(2)
                for i = 0:100-1
                    buf = @sprintf("%03d", i)
                    CImGui.Button(buf, ImVec2(-1.0, 0.0))
                    CImGui.NextColumn()
                end
                CImGui.EndChild()
                CImGui.PopStyleVar()
            end
        end
        CImGui.Separator()

        # Demonstrate a few extra things
        # - Changing ImGuiCol_ChildBg (which is transparent black in default styles)
        # - Using SetCursorPos() to position the child window (because the child window is an item from the POV of the parent window)
        #   You can also call SetNextWindowPos() to position the child window. The parent window will effectively layout from this position.
        # - Using CImGui.GetItemRectMin/Max() to query the "item" state (because the child window is an item from the POV of the parent window)
        #   See "Widgets" -> "Querying Status (Active/Focused/Hovered etc.)" section for more details about this.
        begin
            CImGui.SetCursorPosX(50)
            CImGui.PushStyleColor(CImGui.ImGuiCol_ChildBg, IM_COL32(255, 0, 0, 100))
            CImGui.BeginChild("blah", ImVec2(200, 100), true, CImGui.ImGuiWindowFlags_None)
            for n = 0:50-1
                CImGui.Text("Some test $n")
            end
            CImGui.EndChild()
            child_rect_min = CImGui.GetItemRectMin()
            child_rect_max = CImGui.GetItemRectMax()
            CImGui.PopStyleColor()
            txt = @sprintf("Rect of child window is: (%.0f,%.0f) (%.0f,%.0f)", child_rect_min.x, child_rect_min.y, child_rect_max.x, child_rect_max.y)
            CImGui.Text(txt)
        end

        CImGui.TreePop()
    end

    if CImGui.TreeNode("Widgets Width")
        @cstatic f = Cfloat(0.0) begin
            CImGui.Text("PushItemWidth(100)")
            CImGui.SameLine()
            ShowHelpMarker("Fixed width.")
            CImGui.PushItemWidth(100)
            @c CImGui.DragFloat("float##1", &f)
            CImGui.PopItemWidth()

            CImGui.Text("PushItemWidth(GetWindowWidth() * 0.5)")
            CImGui.SameLine()
            ShowHelpMarker("Half of window width.")
            CImGui.PushItemWidth(CImGui.GetWindowWidth() * 0.5)
            @c CImGui.DragFloat("float##2", &f)
            CImGui.PopItemWidth()

            CImGui.Text("PushItemWidth(GetContentRegionAvail().x * 0.5)")
            CImGui.SameLine()
            ShowHelpMarker("Half of available width.\n(~ right-cursor_pos)\n(works within a column set)")
            CImGui.PushItemWidth(CImGui.GetContentRegionAvail().x * 0.5)
            @c CImGui.DragFloat("float##3", &f)
            CImGui.PopItemWidth()

            CImGui.Text("PushItemWidth(-100)")
            CImGui.SameLine()
            ShowHelpMarker("Align to right edge minus 100")
            CImGui.PushItemWidth(-100)
            @c CImGui.DragFloat("float##4", &f)
            CImGui.PopItemWidth()

            CImGui.Text("PushItemWidth(-1)")
            CImGui.SameLine()
            ShowHelpMarker("Align to right edge")
            CImGui.PushItemWidth(-1)
            @c CImGui.DragFloat("float##5", &f)
            CImGui.PopItemWidth()
        end
        CImGui.TreePop()
    end

    if CImGui.TreeNode("Basic Horizontal Layout")
        CImGui.TextWrapped("(Use CImGui.SameLine() to keep adding items to the right of the preceding item)")

        # text
        CImGui.Text("Two items: Hello"); CImGui.SameLine()
        CImGui.TextColored(ImVec4(1,1,0,1), "Sailor")

        # adjust spacing
        CImGui.Text("More spacing: Hello"); CImGui.SameLine(0, 20)
        CImGui.TextColored(ImVec4(1,1,0,1), "Sailor")

        # button
        CImGui.AlignTextToFramePadding()
        CImGui.Text("Normal buttons"); CImGui.SameLine()
        CImGui.Button("Banana"); CImGui.SameLine()
        CImGui.Button("Apple"); CImGui.SameLine()
        CImGui.Button("Corniflower")

        # button
        CImGui.Text("Small buttons"); CImGui.SameLine()
        CImGui.SmallButton("Like this one"); CImGui.SameLine()
        CImGui.Text("can fit within a text block.")

        # aligned to arbitrary position. Easy/cheap column.
        CImGui.Text("Aligned")
        CImGui.SameLine(150); CImGui.Text("x=150")
        CImGui.SameLine(300); CImGui.Text("x=300")
        CImGui.Text("Aligned")
        CImGui.SameLine(150); CImGui.SmallButton("x=150")
        CImGui.SameLine(300); CImGui.SmallButton("x=300")

        # checkbox
        @cstatic c1=false c2=false c3=false c4=false begin
            @c CImGui.Checkbox("My", &c1)
            CImGui.SameLine()
            @c CImGui.Checkbox("Tailor", &c2)
            CImGui.SameLine()
            @c CImGui.Checkbox("Is", &c3)
            CImGui.SameLine()
            @c CImGui.Checkbox("Rich", &c4)
        end

        # various
        items = ["AAAA", "BBBB", "CCCC", "DDDD"]
        @cstatic f0=Cfloat(1.0) f1=Cfloat(2.0) f2=Cfloat(3.0) item=Cint(-1) begin
            CImGui.PushItemWidth(80)
            @c CImGui.Combo("Combo", &item, items, length(items))
            CImGui.SameLine()
            @c CImGui.SliderFloat("X", &f0, 0.0, 5.0)
            CImGui.SameLine()
            @c CImGui.SliderFloat("Y", &f1, 0.0, 5.0)
            CImGui.SameLine()
            @c CImGui.SliderFloat("Z", &f2, 0.0, 5.0)
            CImGui.PopItemWidth()
        end

        CImGui.PushItemWidth(80)
        CImGui.Text("Lists:")
        @cstatic selection=Cint[0, 1, 2, 3] begin
            for i = 0:4-1
                i > 0 && CImGui.SameLine()
                CImGui.PushID(i)
                CImGui.ListBox("", pointer(selection)+i*sizeof(Cint), items, length(items))
                CImGui.PopID()
                # CImGui.IsItemHovered() && CImGui.SetTooltip("ListBox $i hovered")
            end
            CImGui.PopItemWidth()
        end

        # dummy
        button_sz = ImVec2(40, 40)
        CImGui.Button("A", button_sz)
        CImGui.SameLine()
        CImGui.Dummy(button_sz)
        CImGui.SameLine()
        CImGui.Button("B", button_sz)

        # manually wrapping (we should eventually provide this as an automatic layout feature,
        # but for now you can do it manually)
        CImGui.Text("Manually wrapping:")
        # ImGuiStyle& style = CImGui.GetStyle()
        # buttons_count = 20
        # window_visible_x2 = CImGui.GetWindowPos().x + CImGui.GetWindowContentRegionMax().x
        # for n = 0:buttons_count-1
        #     CImGui.PushID(n)
        #     CImGui.Button("Box", button_sz)
        #     last_button_x2 = CImGui.GetItemRectMax().x
        #     next_button_x2 = last_button_x2 + style.ItemSpacing.x + button_sz.x # expected position if next button was on same line
        #     if n + 1 < buttons_count && next_button_x2 < window_visible_x2
        #         CImGui.SameLine()
        #     end
        #     CImGui.PopID()
        # end

        CImGui.TreePop()
    end

    if CImGui.TreeNode("Tabs")
        if CImGui.TreeNode("Basic")
            tab_bar_flags = CImGui.ImGuiTabBarFlags_None
            if CImGui.BeginTabBar("MyTabBar", tab_bar_flags)
                if CImGui.BeginTabItem("Avocado")
                    CImGui.Text("This is the Avocado tab!\nblah blah blah blah blah")
                    CImGui.EndTabItem()
                end
                if CImGui.BeginTabItem("Broccoli")
                    CImGui.Text("This is the Broccoli tab!\nblah blah blah blah blah")
                    CImGui.EndTabItem()
                end
                if CImGui.BeginTabItem("Cucumber")
                    CImGui.Text("This is the Cucumber tab!\nblah blah blah blah blah")
                    CImGui.EndTabItem()
                end
                CImGui.EndTabBar()
            end
            CImGui.Separator()
            CImGui.TreePop()
        end

        if CImGui.TreeNode("Advanced & Close Button")
            # expose a couple of the available flags. In most cases you may just call BeginTabBar() with no flags (0).
            @cstatic tab_bar_flags=Cint(CImGui.ImGuiTabBarFlags_Reorderable) begin
                @c CImGui.CheckboxFlags("ImGuiTabBarFlags_Reorderable", &tab_bar_flags, CImGui.ImGuiTabBarFlags_Reorderable)
                @c CImGui.CheckboxFlags("ImGuiTabBarFlags_AutoSelectNewTabs", &tab_bar_flags, CImGui.ImGuiTabBarFlags_AutoSelectNewTabs)
                @c CImGui.CheckboxFlags("ImGuiTabBarFlags_TabListPopupButton", &tab_bar_flags, CImGui.ImGuiTabBarFlags_TabListPopupButton)
                @c CImGui.CheckboxFlags("ImGuiTabBarFlags_NoCloseWithMiddleMouseButton", &tab_bar_flags, CImGui.ImGuiTabBarFlags_NoCloseWithMiddleMouseButton)
                if (tab_bar_flags & CImGui.ImGuiTabBarFlags_FittingPolicyMask_) == 0
                    tab_bar_flags |= CImGui.ImGuiTabBarFlags_FittingPolicyDefault_
                end
                result = @c CImGui.CheckboxFlags("ImGuiTabBarFlags_FittingPolicyResizeDown", &tab_bar_flags, CImGui.ImGuiTabBarFlags_FittingPolicyResizeDown)
                if result != 0
                    tab_bar_flags &= ~(CImGui.ImGuiTabBarFlags_FittingPolicyMask_ ⊻ CImGui.ImGuiTabBarFlags_FittingPolicyResizeDown)
                end
                result = @c CImGui.CheckboxFlags("ImGuiTabBarFlags_FittingPolicyScroll", &tab_bar_flags, CImGui.ImGuiTabBarFlags_FittingPolicyScroll)
                if result != 0
                    tab_bar_flags &= ~(CImGui.ImGuiTabBarFlags_FittingPolicyMask_ ⊻ CImGui.ImGuiTabBarFlags_FittingPolicyScroll)
                end

                # Tab Bar
                names = ["Artichoke", "Beetroot", "Celery", "Daikon"]
                @cstatic opened=[true, true, true, true] begin
                    # persistent user state
                    for n = 0:length(opened)-1
                        n > 0 && CImGui.SameLine()
                        CImGui.Checkbox(names[n+1], pointer(opened)+n*sizeof(Bool))
                    end

                    # passing a bool* to BeginTabItem() is similar to passing one to Begin(): the underlying bool will be set to false when the tab is closed.
                    if CImGui.BeginTabBar("MyTabBar", tab_bar_flags)
                        for n = 0:length(opened)-1
                            if opened[n+1] && CImGui.BeginTabItem(names[n+1], pointer(opened)+n*sizeof(Bool))
                                CImGui.Text("This is the $(names[n+1]) tab!")
                                n & 1 != 0 && CImGui.Text("I am an odd tab.")
                                CImGui.EndTabItem()
                            end
                        end
                        CImGui.EndTabBar()
                    end
                end
            end # @cstatic
            CImGui.Separator()
            CImGui.TreePop()
        end
        CImGui.TreePop()
    end

    if CImGui.TreeNode("Groups")
        ShowHelpMarker("Using CImGui.BeginGroup()/EndGroup() to layout items. BeginGroup() basically locks the horizontal position. EndGroup() bundles the whole group so that you can use functions such as IsItemHovered() on it.")
        CImGui.BeginGroup()
        begin
            CImGui.BeginGroup()
            CImGui.Button("AAA")
            CImGui.SameLine()
            CImGui.Button("BBB")
            CImGui.SameLine()
            CImGui.BeginGroup()
            CImGui.Button("CCC")
            CImGui.Button("DDD")
            CImGui.EndGroup()
            CImGui.SameLine()
            CImGui.Button("EEE")
            CImGui.EndGroup()
            CImGui.IsItemHovered() && CImGui.SetTooltip("First group hovered")
        end
        # capture the group size and create widgets using the same size
        size = CImGui.GetItemRectSize()
        values = Cfloat[0.5, 0.20, 0.80, 0.60, 0.25]
        CImGui.PlotHistogram("##values", values, length(values), 0, C_NULL, 0.0, 1.0, size)

        # CImGui.Button("ACTION", ImVec2((size.x - CImGui.GetStyle().ItemSpacing.x)*0.5, size.y))
        CImGui.SameLine()
        # CImGui.Button("REACTION", ImVec2((size.x - CImGui.GetStyle().ItemSpacing.x)*0.5, size.y))
        CImGui.EndGroup()
        CImGui.SameLine()

        CImGui.Button("LEVERAGE\nBUZZWORD", size)
        CImGui.SameLine()

        if CImGui.ListBoxHeader("List", size)
            CImGui.Selectable("Selected", true)
            CImGui.Selectable("Not Selected", false)
            CImGui.ListBoxFooter()
        end

        CImGui.TreePop()
    end
end


text:
function play_turn(state::DynamicIndividualState,
    static_state::StaticIndividualState, decision::Tuple{UInt8, UInt8})
    next_state = state

    fm_pending = get_fast_moves_pending(state)
    cmp = get_cmp(state)

    if !iszero(cmp)
        agent = isodd(cmp) ? 0x01 : 0x02
         next_state = evaluate_charged_move(next_state,
            static_state, cmp, decision[agent] == 0x05 ? 0x01 : 0x02,
            0x64, decision[get_other_agent(agent)] == 0x01)
         if !iszero(fm_pending[get_other_agent(agent)])
             next_state = DynamicIndividualState(next_state[0x01],
                next_state[0x02], next_state.data -
                (fm_pending[get_other_agent(agent)] -
                0x0001) * (get_other_agent(agent) == 1 ? UInt32(9) :
                UInt32(63)))
        end
         next_state = DynamicIndividualState(next_state[0x01],
            next_state[0x02], next_state.data - fm_pending[agent] *
            (agent == 1 ? UInt32(9) : UInt32(63)))
    else
         if fm_pending[1] == UInt32(1) || fm_pending[2] == UInt32(1)
            next_state = evaluate_fast_moves(next_state, static_state,
                (fm_pending[1] == UInt32(1) &&
                get_hp(next_state[0x01]) != 0x0000,
                fm_pending[2] == UInt32(1) &&
                get_hp(next_state[0x02]) != 0x0000))
        end

         next_state = step_timers(next_state,
            decision[1] == 0x03 ?
                get_cooldown(static_state[0x01].fast_move) : Int8(0),
            decision[2] == 0x03 ?
                get_cooldown(static_state[0x02].fast_move) : Int8(0))

        if get_hp(next_state[0x01]) != 0x0000 &&
            get_hp(next_state[0x02]) != 0x0000
            if decision[1] == 0x04
                if decision[2] == 0x04
                    atk_cmp = Base.cmp(
                        static_state[0x01].stats.attack,
                        static_state[0x02].stats.attack
                    )
                    if atk_cmp == 1
                        next_state = DynamicIndividualState(next_state[0x01],
                            next_state[0x02], next_state.data + UInt32(1323))
                    elseif atk_cmp == -1
                        next_state = DynamicIndividualState(next_state[0x01],
                            next_state[0x02], next_state.data + UInt32(1764))
                    else
                        next_state = DynamicIndividualState(next_state[0x01],
                            next_state[0x02], next_state.data + UInt32(11025))
                    end
                else
                    next_state = DynamicIndividualState(next_state[0x01],
                        next_state[0x02], next_state.data + UInt32(441))
                end
            elseif decision[2] == 0x04
                next_state = DynamicIndividualState(next_state[0x01],
                    next_state[0x02], next_state.data + UInt32(882))
            end
        end
    end

    return next_state
end

function resolve_chance(state::DynamicIndividualState,
    static_state::StaticIndividualState)
    chance = get_chance(state)
    if chance == UInt32(0)
        return state
    elseif chance == UInt32(5)
        return rand(rb_rng) < 0.5 ?
            # subtract chance, add cmp
            DynamicIndividualState(state[0x01], state[0x02],
                state.data - UInt32(9702)) :
            DynamicIndividualState(state[0x01], state[0x02],
                state.data - UInt32(9261))
    else
        agent = chance < UInt32(3) ? 0x01 : 0x02
        move = isodd(chance) ? static_state[agent].charged_move_1 :
            static_state[agent].charged_move_2
        if buff_applies(move)
            data = apply_buff(state.data, move, agent)
            return DynamicIndividualState(state[0x01], state[0x02],
                data - chance * UInt32(2205))
        else
            return DynamicIndividualState(state[0x01], state[0x02],
                state.data - chance * UInt32(2205))
        end
    end
end

function play_battle(starting_state::DynamicIndividualState,
    static_state::StaticIndividualState; allow_nothing::Bool = false,
    allow_overfarming::Bool = false,)
    state = starting_state
    while true
        state = resolve_chance(state, static_state)
        d1, d2 = get_possible_decisions(state, static_state,
            allow_nothing = allow_nothing,
            allow_overfarming = allow_overfarming)
        (iszero(d1) || iszero(d2)) && return battle_score(
            state, static_state)
        state = play_turn(state, static_state, select_random_decision(d1, d2))
    end
end

function battle_scores(starting_state::DynamicIndividualState,
    static_state::StaticIndividualState, N::Int64)
    return map(x -> play_battle(starting_state, static_state), 1:N)
end

ext: jl


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_forcearrayplot

"""
    ''_forcearrayplot(;kwargs...)

A ForceArrayPlot component.
The ForceArrayPlot component is used to visualize the shapley contributions
to multiple predictions made by a tree-based ML model. This is a wrapper on
top of React implementation published in shapjs package.
Read more about the component here: https://github.com/slundberg/shap
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique to the component.
- `style` (Dict; optional): Inline css of each element
- `title` (String; optional): Plot title
- `className` (String; optional): html class associated with the component, used for styling purposes
- `baseValue` (Real; optional): same as explainer.expected_value
- `plot_cmap` (a value equal to: 'RdBu', 'GnPR', 'CyPU', 'PkYg', 'DrDb', 'LpLb', 'YlDp', 'OrId' | Array of Strings; optional): The colors used for shap contributions that increase/decrease the prediction value.
Should be one of:
-- default colour combinations RdBu, GnPR, CyPU, PkYg, DrDb, LpLb, YlDp, OrId
-- list of two hex codes, e.g., ["#AAAA11", "#6633CC"]
-- list of two rgb values, e.g., ["rgb(255, 13, 87)", "rgb(30, 136, 229)"]
- `link` (a value equal to: 'identity', 'logit'; optional): either 'identity' or 'logit'
- `featureNames` (Dict with Strings as keys and values of type String; optional): Labels corresponding to each feature, should have same set of keys as "features" prop
- `outNames` (Array of Strings; optional): Single element list of prediction variable name.
- `labelMargin` (Real; optional): Margin (in px) for labels on top of the plot
- `ordering_keys` (Array of Bool | Real | String | Dict | Arrays; optional): X-Axis label for each point
- `ordering_keys_time_format` (String; optional): Formatting for temporal axes, one of d3-time-formats
- `explanations` (optional): List of predictions, where each prediction is a dictionary
describing the predicted value, similarity index and shapley
contributions of each feature. explanations has the following type: Array of lists containing elements 'outValue', 'simIndex', 'features'.
Those elements have the following types:
  - `outValue` (Real; required)
  - `simIndex` (Bool | Real | String | Dict | Array; required)
  - `features` (Dict; required)s
- `clickData` (Dict; optional): attribute for attaching callbacks on click events
"""
function ''_forcearrayplot(; kwargs...)
        available_props = Symbol[:id, :style, :title, :className, :baseValue, :plot_cmap, :link, :featureNames, :outNames, :labelMargin, :ordering_keys, :ordering_keys_time_format, :explanations, :clickData]
        wild_props = Symbol[]
        return Component("''_forcearrayplot", "ForceArrayPlot", "dash_shap_components", available_props, wild_props; kwargs...)
end


text:
using ITensors
using Test

@testset "Algorithm" begin
  alg = ITensors.Algorithm("X")

  @test alg isa ITensors.Algorithm"X"
  @test alg == ITensors.Algorithm"X"()

  s = siteinds("S=1/2", 4)
  A = MPO(s, "Id")
  ψ = randomMPS(s)

  @test_throws MethodError contract(alg, A, ψ)
  @test_throws MethodError contract(A, ψ; method="X")
  @test_throws MethodError contract(A, ψ; alg="X")
  @test contract(ITensors.Algorithm("densitymatrix"), A, ψ) ≈ A * ψ
  @test contract(ITensors.Algorithm("naive"), A, ψ) ≈ A * ψ
  @test contract(A, ψ; alg="densitymatrix") ≈ A * ψ
  @test contract(A, ψ; method="densitymatrix") ≈ A * ψ
  @test contract(A, ψ; alg="naive") ≈ A * ψ
  @test contract(A, ψ; method="naive") ≈ A * ψ

  B = copy(A)
  truncate!(ITensors.Algorithm("frobenius"), B)
  @test A ≈ B

  B = copy(A)
  truncate!(B; alg="frobenius")
  @test A ≈ B

  # Custom algorithm
  function ITensors.truncate!(::ITensors.Algorithm"my_new_algorithm", A::MPO; cutoff=1e-15)
    return "my_new_algorithm was called with cutoff $cutoff"
  end
  cutoff = 1e-5
  res = truncate!(A; alg="my_new_algorithm", cutoff=cutoff)
  @test res == "my_new_algorithm was called with cutoff $cutoff"
end


text:
module Analysis

using Oceananigans
using JLD2
using Plots

using Oceananigans.Grids
using Oceananigans.Fields: offset_data
import Oceananigans.Fields: Field

Field(loc::Tuple, grid::AbstractGrid, raw_data::Array) = Field(loc, CPU(), grid, nothing, offset_data(raw_data, grid, loc))
field_from_file(file, loc, name, iter) = Field(loc, file["serialized/grid"], file["timeseries/$name/$iter"])

CellFileField(file, name, iter)  = Field((Cell, Cell, Cell), file["serialized/grid"], file["timeseries/$name/$iter"])
XFaceFileField(file, name, iter) = Field((Face, Cell, Cell), file["serialized/grid"], file["timeseries/$name/$iter"])
YFaceFileField(file, name, iter) = Field((Cell, Face, Cell), file["serialized/grid"], file["timeseries/$name/$iter"])

max_abs_tracer(file, iter) = maximum(abs, interior(CellFileField(file, :c, iter)))

tracer_variance(file, iter) = 1/2 * sum(interior(CellFileField(file, :c, iter)).^2)

kinetic_energy(file, iter) = 1/2 * (sum(interior(XFaceFileField(file, :u, iter)).^2) +
                                    sum(interior(YFaceFileField(file, :v, iter)).^2))

function oceananigans_statistics(filename)
    file = jldopen(filename)

    iterations = parse.(Int, keys(file["timeseries/t"]))

    t  = [ file["timeseries/t/$iter"]  for iter in iterations]
    e  = [ kinetic_energy(file, iter)  for iter in iterations]
    χ  = [ tracer_variance(file, iter) for iter in iterations]
    c★ = [ max_abs_tracer(file, iter)  for iter in iterations]

    close(file)

    Δχ = χ ./ χ[1]
    Δe = e ./ e[1]

    return t, Δe, Δχ, c★
end

function plot_resolutions(t, data, resolutions; kwargs...)
    plt = plot(t[1], data[1]; label = "N² = $(resolutions[1])²", kwargs...)
    for i = 2:length(data)
        plot!(plt, t[i], data[i]; label = "N² = $(resolutions[i])²", kwargs...)
    end
    return plt
end

end # module


text:
function ηs_from_data!(ηs,D,t_star)
    ηs[:,1] = exp.(t_star .* D)
    m = size(ηs,2)
    for i=1:(m-1)
        ηs[:,i+1] = ηs[:,i] .* D ./ i
    end
end
function η_from_data!(η,D,t_star)

    slack = exp.(t_star .* D)
    n = length(slack)
    η[1] = sum(slack)/n
    m = length(η)
    for i=1:(m-1)
        slack = slack .* D ./ i
        η[i+1] = sum(slack)/n
    end
    return η
end
function τ_from_η!(τ,η)
    η[2:end] ./= η[1]
    τ[1] = log(η[1])
    τ[2] = η[2]
    for k in 3:length(η)
        τ[k] = (k-1)*η[k]
        for j in 1:(k-2)
            τ[k] -= τ[j+1] * η[k-j]
        end
    end
    return τ
end

"""
    thorin_moments(D,t_star,m)

D should be a vector of samples from a distribution, t_star is usually taking to be -1, and m is an integer.

This compute the bare thorin moments, defined as a rescaling of the sample biaised `t_star`-shifted cumulants.
"""
function thorin_moments(D::Vector{T},t_star,m) where T
    τ = zeros(T,m+1)
    η = zeros(T,m+1)
    @views η_from_data!(η,D,t_star)
    τ_from_η!(τ,η)
    return τ
end
@inline κ_from_τ(τ,P::PreComp{T}) where T = τ .* [T(1),P.FACTS[1:(length(τ)-1)]...]
function μ_from_κ(κ,P::PreComp{T}) where T
    μ = zero(κ)
    μ[1] = exp(κ[1])
    for k in 2:length(κ)
        for j in 1:k-1
            μ[k] += μ[j] * κ[k-j+1] * P.BINS[j, k-1] #binomial(big(k-2),big(j-1))
        end
    end
    return μ
end
@inline a_from_μ(μ,P::PreComp{T}) where T = sqrt(2) .* (P.LAGUERRE[1:length(μ),1:length(μ)]'μ)
function empirical_coefs(D::Vector{T},m::M) where {M<:Int, T}

    P = get_precomp(T,m)::PreComp{T}
    η = zeros(T,m)
    τ = zeros(T,m)
    μ = zeros(T,m)

    # initialize the computations (k=1)
    slack = exp.(-D)
    η[1] = sum(slack)
    μ[1]= η[1]/length(D)

    # Loop (k > 1)
    @inbounds for k in 2:m
        slack .*= D ./ (k-1)
        η[k] = sum(slack)/η[1]
        τ[k] = (k-1)*η[k]
        @inbounds for j in 2:k-1
            τ[k] -= τ[j]*η[k-j+1]
            μ[k] += τ[j]*μ[k-j+1]
        end
        μ[k] += μ[1]*τ[k]
        μ[k] /= (k-1)
    end
    return P.LAGUERRE2[1:m,1:m]'μ
end

function thorin_moments_old(D::Vector{T},t_star,m) where T
    n = length(D)
    τ = zeros(T,m+1)
    η = zeros(T,m+1)
    ηs = zeros(T,(n,m+1))
    @views ηs_from_data!(ηs,D,t_star)
    η = dropdims(Statistics.mean(ηs,dims=1),dims=1)
    τ_from_η!(τ,η)
    return τ
end

"""
resemps_thorin_moments(M,D,t_star,m)

`M` is an integer
`D` is the data, as a (d,n)-shaped matrix, where d is the numebr of dimensions and n the number of samples.
`t_star` is usually taking to be -1
`m` is an integer.

This resamples the bare thorin moments M times, more efficienty than a simple bootstrap of the `thorin_moments(D,t_star,m)` function.
"""
function resemps_thorin_moments(M,D::AbstractMatrix{T},t_star,m) where {T}
    d,n = size(D)
    ηs = zeros(T,(d,n,m+1))
    τ = zeros(T,(M,d,m+1))
    η = zeros(T,(M,d,m+1))
    for dim in 1:d
        @views ηs_from_data!(ηs[dim,:,:],D[dim,:],t_star)
    end

    # Main loop:
    for i in 1:M
        println(i)
        η[i,:,:] = dropdims(Statistics.mean(ηs[:,StatsBase.sample(1:n,n,replace=true),:],dims=2),dims=2)
        for dim in 1:d
           @views τ_from_η!(τ[i,dim,:],η[i,dim,:])
        end
    end
    return τ
end
function resemps_thorin_moments(M,D::Vector{T},t_star,m) where T
    return resemps_thorin_moments(M,reshape(D,(1,length(D))),t_star,m)[:,1,:]
end


text:
nbranch(t::Callsite) = nbranch(flatten(t))
function nbranch(stream)
    count = 0
    for (f, _, _) in stream
        if f == assert
            count+=1
        end
    end
    return count
end

function cut_at_nth_branch(stream, nth_branch)
    cut_stream = Tuple{Any, Any, Tuple}[]
    count = 0
    for (f, ret, args) in stream
        push!(cut_stream, (f, ret, args))
        if f == assert
            count += 1
            if count == nth_branch
                return cut_stream
            end
        end
    end
    return cut_stream
end

function force_and_cut(stream, nth_branch)
    cut_stream = cut_at_nth_branch(stream, nth_branch)
    f, ret, args = pop!(cut_stream)
    @assert f == assert
    @assert anything(ret) isa Bool
    push!(cut_stream, (f, Some(!anything(ret)), args))
end

function checkStream(stream)
    z3model = symbolic(stream)
    sat, model = runZ3(z3model)
    if sat
       inputs = parseZ3(model)
    else
       inputs = ()
    end

    return sat, inputs
end

import InteractiveUtils: subtypes

function concretize(T::Union{DataType, Union})
    ctypes = Any[]
    types = Any[T]

    while !isempty(types)
        T = pop!(types)
        if isconcretetype(T)
            push!(ctypes, T)
        elseif T isa Union
            push!(types, T.a)
            push!(types, T.b)
        else
            append!(types, subtypes(T))
        end
    end
    return ctypes
end

supported(T::DataType)= false
generate(T::DataType) = error("Can't generate values of $T")

supported(::Type{T}) where T<:INTEGERS = true
# generate(::Type{BigInt}) = rand(1:big(typemax(Int))^2)
generate(::Type{T}) where T<:INTEGERS = rand(T)

# UGLY HACK
function enumerateSupportedTypes(i)
    if i == 0
      return Int64
    elseif i == 1
      return UInt64
    elseif i == 2
      return Int128
    elseif i == 3
      return UInt128
    elseif i == 4
      return Int32
    elseif i == 5
      return UInt32
    elseif i == 6
      return Int16
    elseif i == 7
      return UInt16
    elseif i == 8
      return Int8
    elseif i == 9
      return UInt8
    # elseif i == 10
    #   return BigInt
    else
      return Bool
    end
end

function fuzz(f, argtypes...; maxdepth = typemax(Int64))
    worklist = Any[]
    args_ctypes = map(x->Base.filter(supported, concretize(x)), argtypes)
    for types in Iterators.product(args_ctypes...)
        initial_args = map(generate, types)
        push!(worklist, (0, initial_args, nothing))
    end
    fuzz_worklist(f, worklist, maxdepth)
end

function fuzz_wargs(f, initial_args...; maxdepth = typemax(Int64))
    worklist = Any[(0, initial_args, nothing)]
    fuzz_worklist(f, worklist, maxdepth)
end

###
#    Iterative breath first tree search
#      - Invalidated earliest branch and use Z3 to generate an example for the
#        opposite branch
#      - explore all sides and iterate through the program to discover all branches
function fuzz_worklist(f, worklist::Vector{Any}, maxdepth)
    tested = Any[] # did not error
    errored = Any[]

    while !isempty(worklist)
        depth, args, subs = pop!(worklist)
        @debug "Executing item from worklist" depth args subs
        val, _, trace, record = execute(f, args...; subs = subs)

        if val isa Exception
            push!(errored, (val, args, subs))
        else
            push!(tested, (val, args, subs))
        end

        stream = filter(trace)
        branches = Any[]
        for n in (depth+1):nbranch(stream)
            cut_stream = force_and_cut(stream, n)
            push!(branches, (n, cut_stream))
        end

        @info "Found $(length(branches)) new branches to explore"

        for (d, branch) in branches
            if d >= maxdepth
                @info "Terminated fuzzing that went to deep" maxdepth
                continue
            end
            sat, inputs = try
                checkStream(branch)
            catch ex
                @error "Error in Z3 run, skipping" exception=ex branch
                continue
            end
            if sat
                args, subs, _ = inputs
                subs = augment(record, subs)
                # record to subs
                push!(worklist, (d, args, subs))
            end
        end
    end
    tested, errored
end


text:
const testfiles = [
    "ldselect" => "./ldselect.jl",
    "informativeness" => "./informativeness.jl",
]

const args = let x = get(ENV, "TEST_ARGS", "")
    x === "" ? "all" : x
end

const to_test = args === "all" ? [
    (key, value) for (key, value) in testfiles
   ] : [(args, testfiles[findfirst(x -> x[1] === args, testfiles)][2])]

for (key, testfile) in to_test
    @info "Testing $key..."
    @eval module $(gensym())
        include($testfile)
    end
end


text:
@testset "Utilities" begin
  @test_throws ArgumentError readgeotable("doesnotexist.csv")
  @test_throws AssertionError readgeotable(joinpath(datadir,"data3D.tsv"), delim='\t', coordnames=(:a,))

  @testset "spheredir" begin
    @test spheredir(90, 0) ≈ [1,0,0]
    @test spheredir(90,90) ≈ [0,1,0]
    @test spheredir(0,  0) ≈ [0,0,1]
  end
end


text:
using Base.Test
using MPI

MPI.Init()

function scatterv_array(A, counts::Vector{Cint}, root)
    comm = MPI.COMM_WORLD
    T = eltype(A)
    if MPI.Comm_rank(comm) == root
        B = copy(A)
    else
        B = Array(T,1)
    end
    C = MPI.Scatterv(B, counts, root, comm)
end

comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)
root = 0

# Defining this to make ones work for Char
Base.one(::Type{Char}) = '\01'

for typ in MPI.MPIDatatype.types

    A = ones(typ, 3 * div(size,2) + mod(size,2))
    counts = Cint[ mod(i,2) + 1 for i in 0:(size-1)]
    B = scatterv_array(A, counts, root)
    @test B == ones(typ, mod(rank,2) + 1)
end

MPI.Finalize()


text:
#####################
# Vector operations #
#####################

#Reduction operations

import Base: sum, maximum, minimum, any, max, min, countnz

#Sum
sum(a::AFAbstractArray) = af_sum(a)
sum{T}(a::AFAbstractArray{T}, dim::Integer) = AFArray{T}(af_sum(a,dim-1))

#Product
product(a::AFAbstractArray) = af_product(a)
product{T}(a::AFAbstractArray{T}, dim::Integer) = AFArray{T}(af_product(a,dim-1))

#Maximum
maximum(a::AFAbstractArray) = af_max(a)
maximum{T}(a::AFAbstractArray{T}, dim::Integer) = AFArray{T}(af_max_dim(a,dim-1))

#Minimum
minimum(a::AFAbstractArray) = af_min(a)
minimum{T}(a::AFAbstractArray{T}, dim::Integer) = AFArray{T}(af_min_dim(a,dim-1))

#Any
any(a::AFAbstractArray) = af_anyTrue(a)
any(a::AFAbstractArray, dim::Integer) = AFArray{Bool}(af_anyTrue(a, dim-1))

#Alltrue
alltrue(a::AFAbstractArray) = af_allTrue(a)
alltrue(a::AFAbstractArray, dim::Integer) = AFArray{Bool}(icxx"af::allTrue($a, $(dim - 1));")

#Count
countnz(a::AFAbstractArray) = af_count(a)
countnz(a::AFAbstractArray, dim::Integer) = AFArray{UInt32}(af_count(a, dim-1))

#Inclusive Scan Operations

import Base: cumsum, find

#Cumsum
cumsum{T}(a::AFAbstractArray{T}, dim::Integer = 1) = AFArray{T}(af_accum(a,dim-1))

#Find
find(a::AFAbstractArray) = AFArray{UInt32}(af_where(a))

#Numerical differentiation

import Base: diff, gradient

#1D Diff
diff{T}(a::AFAbstractArray{T}, dim::Integer = 1) = AFArray{T}(af_diff(a,dim-1))

#2D Siff
diff2{T}(a::AFAbstractArray{T}, dim::Integer = 1) = AFArray{T}(af_diff2(a,dim-1))

#Gradient
function gradient{T}(a::AFAbstractArray{T})
    dx = AFArray()
    dy = AFArray()
    af_grad(dx, dy, a)
    AFArray{T}(dx), AFArray{T}(dy)
end

#Sorting

import Base: sort
export sortIndex

#Sort
#NOTE: ArrayFire currently restricts dim to the value 0, which means you can sort only columns
function sort{T}(A::AFAbstractArray{T}; rev = false)
    ndims(A) == 1 ||
        error("Must explicitly specify dimension when sorting mutlidimensional array")
    AFArray{T}(af_sort(A,0,!rev))
end
function sort{T}(A::AFAbstractArray{T}, dim; rev = false)
    if dim > 1
        error("ArrayFire currently lets you sort only by columns")
    else
        AFArray{T}(af_sort(A,dim-1,!rev))
    end
end

#sortIndex
function sortIndex(a::AFArray, dim = 1; rev = false)
    (ndims(a) == 2 && dim == 2) && error("ArrayFire currently supports sorting along columns")
    out  = AFArray()
    ind = AFArray()
    af_sort(out, ind, a, dim-1, !rev)
    AFArray{backend_eltype(out)}(out), AFArray{backend_eltype(ind)}(ind)
end

#sortByKey
function sortByKey(keys::AFArray, vals::AFArray, dim = 1; rev = false)
    (ndims(a) == 2 && dim == 2) && error("ArrayFire currently supports sorting along columns")
    out_keys = AFArray()
    out_vals = AFArray()
    af_sort(out_keys, out_vals, keys, vals, dim-1, !rev)
    AFArray{backend_eltype(out_keys)}(out_keys), AFArray{backend_eltype(out_vals)}(out_vals)
end

#Set Operations

import Base: setdiff, union, unique

#Setdiff
setdiff{T}(a::AFAbstractArray{T} , b::AFAbstractArray{T}) = AFArray{T}(af_setIntersect(a,b))

#Union
union{T}(a::AFAbstractArray{T}, b::AFAbstractArray{T}) = AFArray{T}(af_setUnion(a,b))

#Unique
unique{T}(a::AFAbstractArray{T}) = AFArray{T}(af_setUnique(a))

#Reordering functions

import Base:flipdim, vec

flipdim{T}(a::AFAbstractArray{T}, dim::Integer) = AFArray{T}(af_flip(a,dim-1))
vec{T}(a::AFAbstractArray{T}) = AFArray{T}(af_flat(a))
function union{T}(dim::Integer, a::AFAbstractArray{T}, b::AFAbstractArray)
    AFArray{T}(af_join(dim-1, a, b))
end
function union{T}(dim::Integer, a::AFAbstractArray{T}, b::AFAbstractArray{T}, c::AFAbstractArray{T})
    AFArray{T}(af_join(dim-1, a, b, c))
end
function union{T}(dim::Integer, a::AFAbstractArray{T}, b::AFAbstractArray{T},
                    c::AFAbstractArray{T}, d::AFAbstractArray{T})
    AFArray{T}(af_join(dim-1, a, b, c, d))
end
moddims{T}(a::AFAbstractArray{T}, dims...) = AFArray{T}(af_moddims(a,dims_to_dim4(dims)))
function reorder{T}(a::AFAbstractArray{T}, x::Integer, y::Integer = 2, z::Integer = 3, w::Integer = 4)
    AFArray{T}(af_reorder(a, x-1, y-1, z-1, w-1))
end
replace{T}(a::AFAbstractArray{T}, cond::AFAbstractArray{Bool}, b::AFAbstractArray{T}) = AFArray{T}(af_replace(a, cond, b))
replace{T}(a::AFAbstractArray{T}, cond::AFAbstractArray{Bool}, b::Real) = AFArray{T}(af_replace(a, cond, b))
select{T}(cond::AFAbstractArray{Bool}, a::AFAbstractArray{T}, b::AFAbstractArray{T}) = AFArray{T}(af_select(cond, a, b))
select{T}(cond::AFAbstractArray{Bool}, a::AFAbstractArray{T}, b::Real) = AFArray{T}(af_select(cond, a, b))
function shift{T}(a::AFAbstractArray{T}, x::Integer, y::Integer = 1, z::Integer = 1, w::Integer = 1)
    AFArray{T}(af_shift(a, x-1, y-1, z-1, w-1))
end

#Some array helper functions
import Base: isnan, isinf
conj{T}(a::AFAbstractArray{T}) = AFArray{T}(af_complex_conj(a))
convert{T,V}(::Type{AFArray{T}}, a::AFAbstractArray{V}) = AFArray{T}(af_as(a, aftype(T)))
isnan(a::AFAbstractArray) = AFArray{Bool}(af_isNaN(a))
iszero(a::AFAbstractArray) = AFArray{Bool}(af_iszero(a))
isinf(a::AFAbstractArray) = AFArray{Bool}(af_isinf(a))


text:
function basis(elem::Hex, N, r, s, t)
    Np = convert(Int, (N+1)^3)
    V, Vr, Vs, Vt = ntuple(x->zeros(length(r), Np), 4)
    tmp_array = zeros(N+1)
    for ii in eachindex(r, s, t)
        sk = 1
        for k in 0:N
            P_k = jacobiP(t[ii], 0, 0, k, tmp_array)
            dP_k = grad_jacobiP(t[ii], 0, 0, k, tuple(tmp_array))
            for i in 0:N
                P_i = jacobiP(r[ii], 0, 0, i, tmp_array)
                dP_i = grad_jacobiP(r[ii], 0, 0, i, tuple(tmp_array))
                for j in 0:N
                    P_j = jacobiP(s[ii], 0, 0, j, tmp_array)
                    dP_j = grad_jacobiP(s[ii], 0, 0, j, tuple(tmp_array))

                    V[ii, sk]  = P_i * P_j * P_k
                    Vr[ii, sk] = dP_i * P_j * P_k
                    Vs[ii, sk] = P_i * dP_j * P_k
                    Vt[ii, sk] = P_i * P_j * dP_k

                    sk += 1
                end
            end
        end
    end

    return V, Vr, Vs, Vt
end

# ===================================================

function nodes(elem::Hex, N)
    r1D, _ = gauss_lobatto_quad(0,0,N)
    return vec.(meshgrid(r1D,r1D,r1D))
end

function equi_nodes(elem::Hex, N)
    r1D = LinRange(-1,1,N+1)
    return vec.(meshgrid(r1D,r1D,r1D))
end

function quad_nodes(elem::Hex, N)
    r1D,w1D = gauss_quad(0,0,N)
    r,s,t = vec.(meshgrid(r1D,r1D,r1D))
    wr,ws,wt = vec.(meshgrid(w1D,w1D,w1D))
    w = @. wr*ws*wt
    return r,s,t,w
end


text:
### Define your own clawpack path ###
if haskey(ENV, "CLAW")
    const CLAW = ENV["CLAW"]
else
    ## CLAW="/path/to/top/level/clawpack"
    println("ENV[\"CLAW\"] is not defined.")
    println("Set the env like the following in your default shell:")
    println("export CLAW = \"/path/to/top/level/clawpack\" ")
    CLAW = ""
    #if !isdir(CLAW); error("CLAW=$CLAW is not correct."); end
end


text:
using PyPlot
using ImplicitDomainQuadrature

xrange = -1:0.2:1
points = ImplicitDomainQuadrature.tensor_product_points(xrange',xrange')

fig,ax = PyPlot.subplots()
ax.scatter(points[1,:],points[2,:])
ax.set_aspect("equal")
ax.grid()
fig.savefig("grid.png")


text:
using ModelingToolkit, OrdinaryDiffEq, Symbolics, Test

@parameters t σ ρ β
@variables x(t) y(t) z(t) F(t) u(t)
D = Differential(t)

eqs = [D(x) ~ σ * (y - x) + F,
    D(y) ~ x * (ρ - z) - y,
    D(z) ~ x * y - β * z]

aliases = [u ~ x + y - z]
lorenz1 = ODESystem(eqs, pins = [F], observed = aliases, name = :lorenz1)
lorenz2 = ODESystem(eqs, pins = [F], observed = aliases, name = :lorenz2)

connections = [lorenz1.F ~ lorenz2.u,
    lorenz2.F ~ lorenz1.u]
connected = ODESystem(Equation[], t, [], [], observed = connections,
                      systems = [lorenz1, lorenz2])

sys = connected

@variables lorenz1₊F lorenz2₊F
@test pins(connected) == Variable[lorenz1₊F, lorenz2₊F]
@test isequal(observed(connected),
              [connections...,
                  lorenz1.u ~ lorenz1.x + lorenz1.y - lorenz1.z,
                  lorenz2.u ~ lorenz2.x + lorenz2.y - lorenz2.z])

collapsed_eqs = [
    D(lorenz1.x) ~ (lorenz1.σ * (lorenz1.y - lorenz1.x) +
                    (lorenz2.x + lorenz2.y - lorenz2.z)),
    D(lorenz1.y) ~ lorenz1.x * (lorenz1.ρ - lorenz1.z) - lorenz1.y,
    D(lorenz1.z) ~ lorenz1.x * lorenz1.y - (lorenz1.β * lorenz1.z),
    D(lorenz2.x) ~ (lorenz2.σ * (lorenz2.y - lorenz2.x) +
                    (lorenz1.x + lorenz1.y - lorenz1.z)),
    D(lorenz2.y) ~ lorenz2.x * (lorenz2.ρ - lorenz2.z) - lorenz2.y,
    D(lorenz2.z) ~ lorenz2.x * lorenz2.y - (lorenz2.β * lorenz2.z)]

simplifyeqs(eqs) = Equation.((x -> x.lhs).(eqs), simplify.((x -> x.rhs).(eqs)))

@test isequal(simplifyeqs(equations(connected)), simplifyeqs(collapsed_eqs))

# Variables indicated to be input/output
@variables x [input = true]
@test hasmetadata(x, Symbolics.option_to_metadata_type(Val(:input)))
@test getmetadata(x, Symbolics.option_to_metadata_type(Val(:input))) == true
@test !hasmetadata(x, Symbolics.option_to_metadata_type(Val(:output)))
@test_throws KeyError getmetadata(x, Symbolics.option_to_metadata_type(Val(:output)))

@variables y [output = true]
@test hasmetadata(y, Symbolics.option_to_metadata_type(Val(:output)))
@test getmetadata(y, Symbolics.option_to_metadata_type(Val(:output))) == true
@test !hasmetadata(y, Symbolics.option_to_metadata_type(Val(:input)))
@test_throws KeyError getmetadata(y, Symbolics.option_to_metadata_type(Val(:input)))


text:
### Begin here
println("""Initialize project:""")
    srcDir  = dirname(@__FILE__)*"\\";
    include(srcDir*"Examples\\Examples.jl")
    include(srcDir*"Examples\\Utils_states.jl")
    include(srcDir*"Moments\\Moments.jl")
    include(srcDir*"Constraints\\C_constraints.jl")
    include(srcDir*"Constraints\\R_constraints.jl")
    include(srcDir*"Constraints\\Utils_cons.jl")
    include(srcDir*"Model\\R_sep_Model.jl")
    include(srcDir*"Model\\C_sep_Model.jl")
    include(srcDir*"Model\\Utils_Model.jl")
    include(srcDir*"sep_Compute.jl")

    using .sep_Compute
    using .Examples
    using .Utils_states
    using .Moments
    using .Utils_cons
    using .C_constraints
    using .R_constraints
    using .Utils_Model
    using .R_sep_Model
    using .C_sep_Model
    using .sep_Compute

    using LinearAlgebra
    using JuMP

## Example
t = (4,4)
d = (2,2)

xx̄yȳMM_blocks = Moments.get_xx̄yȳMM_blocks(d,t)
γγᶥζζᶥ_δδᶥηηᶥ = Moments.get_γγᶥζζᶥ_δδᶥηηᶥ(d,t)
ααᶥββᶥ= B[1,1]

MMexᴿ,MMCoefᴿ = Moments.get_ℜℑααᶥββᶥᴿ(d,xx̄yȳMM_blocks[1,1],γγᶥζζᶥ_δδᶥηηᶥ)
display(MMCoefᴿ[1:5,1:5]) ### Mistake, Zero diagonal ??!

split_expo(ααᶥββᶥ,d) = (ααᶥββᶥ[1:d[1]],
                            ααᶥββᶥ[d[1]+1:2*d[1]],
                            ααᶥββᶥ[1+2*d[1]:2*d[1]+d[2]],
                            ααᶥββᶥ[1+2*d[1]+d[2]:end])
# γ,γᶥ,δ,δᶥ = split_expo([0,0,0, 0,0,0, 0,0,0, 0,0,0],(3,3))
# η,ηᶥ,ζ,ζᶥ = split_expo([0,0,0, 0,0,0, 0,0,0, 0,0,0],(3,3))
γ,γᶥ,ζ,ζᶥ = split_expo(γγᶥζζᶥ,(2,2))
δ,δᶥ,η,ηᶥ = split_expo(δδᶥηηᶥ,(2,2))
# γ,γᶥ,δ,δᶥ = split_expo(γγᶥδδᶥ,d)
# η,ηᶥ,ζ,ζᶥ = split_expo(ζζᶥηηᶥ,d)
ter_i(δᶥ,ηᶥ) = ((-1.0)^sum(δᶥ+ηᶥ))
ter_1(δ,δᶥ,η,ηᶥ) = (1.0im)^sum(δ+δᶥ+η+ηᶥ)
ter_comb(γ,γᶥ,ζ,ζᶥ,δ,δᶥ,η,ηᶥ) = (pf(γ+δ)*pf(γᶥ+δᶥ)*pf(η+ζ)*pf(ηᶥ+ζᶥ))/prod(map(x ->pf(x),[γ,γᶥ,ζ,ζᶥ,δ,δᶥ,η,ηᶥ]))
ter_i(δᶥ,ηᶥ)
ter_1(δ,δᶥ,η,ηᶥ)
ter_comb(γ,γᶥ,ζ,ζᶥ,δ,δᶥ,η,ηᶥ)

ter_i(δᶥ,ηᶥ)*ter_1(δ,δᶥ,η,ηᶥ)*ter_comb(γ,γᶥ,ζ,ζᶥ,δ,δᶥ,η,ηᶥ)
#

pf(x) = prod(factorial.(x))
function get_coef(part1,part2,d)
    γ,γᶥ,ζ,ζᶥ = split_expo(part1,d) ; δ,δᶥ,η,ηᶥ = split_expo(part2,d)
    a = (-1.0^sum(δᶥ+ηᶥ))*(1.0im^sum(δ+δᶥ+η+ηᶥ))*prod(pf.([[γ,γᶥ,ζ,ζᶥ]+[δ,δᶥ,η,ηᶥ]...]))
    return a/prod(pf.([γ,γᶥ,ζ,ζᶥ,δ,δᶥ,η,ηᶥ]))
end
γ,γᶥ,ζ,ζᶥ = split_expo([0,0,0,0,0,0,0,0],d)
δ,δᶥ,η,ηᶥ = split_expo([0,0,0,0,0,0,0,0],d)
a = ((-1.0)^sum(δᶥ+ηᶥ))*(1.0im^sum(δ+δᶥ+η+ηᶥ))*prod(pf.([[γ,γᶥ,ζ,ζᶥ]+[δ,δᶥ,η,ηᶥ]...]))
get_coef(,[0,0,0,0,0,0,0,0],d)

# df = Examples.get_example_meta()
# select(df,df.Example == "Eq14Hor97a")
# ρ  = Examples.get_example(36)
# d  = Examples.get_example_meta(36).Size //Error the names and ρ don't match up
ρ =  [1 0 0 0
      0 0 0 0
      0 0 0 0
      0 0 0 1]
d  = (2,2)
t  = (2,2)
cl = ""
cl = "S∞" # "S₂ S₂₁"
cl = "S₂" # "S₂₁"
cl = "S₂₁sG"

sep_mod_ub = C_sep_Model.Modelξₜˢᵉᵖ(ρ,d,t;con_list = cl,noBlock=true)
sep_mod_b  = C_sep_Model.Modelξₜˢᵉᵖ(ρ,d,t;con_list = cl)
sep_mod_ub = sep_Compute.Computeξₜˢᵉᵖ(sep_mod_ub)
sep_mod_b  = sep_Compute.Computeξₜˢᵉᵖ(sep_mod_b)

# Create variables
##      PSD Moment matrix blocks
PSD_con = C_constraints.make_PSD_con(d,t,Lx;noBlock=false)
PSD_con[0,0][1:3,1:3]

PSD_con[0,1]

display(PSD_con[0,0][5,10])
display(PSD_con[0,0][10,5])
MMexᴿ[0,1]

## f
using JuMP
model = JuMP.Model()
@variable(model, Lx[C_constraints.make_mon_expo_keys(d,t[1])] )## Create variables

d₁,d₂   = d
n       = sum(2 .*d)
MMexᴿ,MMCoefᴿ = Moments.get_ℂ_block_diag(d,t .- 1)

B = MMexᴿ[0,-1]
C = MMCoefᴿ[0,-1]

## Zero propagation: ρ = aaᵀ⊗bbᵀ , ρᵢⱼᵢⱼ = 0 ⟹ aᵢbⱼ = 0 ⟹ L(xᵢyⱼ) = 0
# zp = ccon.zeroprop(ρ,d,t,Lx)
# @constraint(model, zp .== zeros(size(zp)))
##      PSD Moment matrix blocks
PSD_con = C_constraints.make_PSD_con(d,t,Lx)
set_con(PSD_con)

PSD_con[(0,-1)]
M = Array{GenericAffExpr{Float64,VariableRef},2}

L_xx̄ᵀ_tens_yȳᵀ = C_constraints.make_ord4_con(d,Lx)

## Batch
include(srcDir*"batch.jl")
using .batch
t = (2,2)
df = Examples.get_example_meta()
batch.batch_model(t,Examples.get_examples(),df)
batch.batch_Computeξₜˢᵉᵖ("C:\\Users\\andries\\all-my-codes\\ju-sep-rank\\assets\\bounds\\")
nar = batch.unstack_constraints(df)

## Solve the OP.

## The real model.
println("""Specify: ρ, registers: d, constraint list: cl,level of hierarchy: t
ρ =  [1 0 0 0
      0 0 0 0
      0 0 0 0
      0 0 0 1]
d  = (2,2)
cl = "S∞ sG" or "S₂ sG" or "S₂₁ sG"
t  = (2,2)""")
ρ  =  [1. 0. 0. 0.
       0. 0. 0. 0.
       0. 0. 0. 0.
       0. 0. 0. 1.]
d  = (2,2)
cl = "S∞ sG"
t  = (3,3)
println("""Generate the model:
sep_mod = R_sep_Model.Modelξₜˢᵉᵖ(ρ,d,t;con_list = cl)
""")
sep_mod = R_sep_Model.Modelξₜˢᵉᵖ(ρ,d,t;con_list = cl)
display(sep_mod)

## Solve the OP.
println("We initialize the computation module")
sep_mod = sep_Compute.Computeξₜˢᵉᵖ(sep_mod)
Lx = sep_mod[:Lx]
## A closer look at the real constraints
println("""Let us take a deeper look into the constraints:""")
println("These are $(binomial(sum(d)+2*t[1],2*t[1])) the moments present in the hierachy at level t = (2,2)")
mon_expo_keys = R_constraints.make_mon_expo_keys(sum(d),t[1])
display(mon_expo_keys)

println("""These are the blocks of the 'L([x,y]ₜ[x,y]ₜᵀ)'.
 """)

zero_moms = R_constraints.zeroprop(ρ,d,(3,3),Lx)

PSD_con = R_constraints.make_PSD_con(d,t,Lx)
display(PSD_con)
println("""
"ee": moments of the form xᵅyᵝ ⋅ xᵞyᵟ where |α|,|β|,|γ|,|δ| are all even.
"eo": moments of the form xᵅyᵝ ⋅ xᵞyᵟ where |α|,|γ| even and |β|,|δ| odd.
"oe": moments of the form xᵅyᵝ ⋅ xᵞyᵟ where |α|,|γ| odd and |β|,|δ| even.
"oo": moments of the form xᵅyᵝ ⋅ xᵞyᵟ where |α|,|β|,|γ|,|δ| are all odd.
 """)
display(PSD_con["ee"])
display(PSD_con["eo"])
display(PSD_con["oe"])
display(PSD_con["oo"])

println("This is 'L(xxᵀ⊗yyᵀ)'.")
ord4_con = R_constraints.make_ord4_con(d,Lx)
display(ord4_con)

println("This is the 'L(( √ρₘₐₓ - xᵢ²)⋅η), L(( √ρₘₐₓ - yⱼ²)⋅η) ⪰ 0 ∀ i ∈ [d₁], j ∈ [d₂]' constraints.")
cons_S_inf = R_constraints.make_loc_cons_S_inf(ρ,d,t,Lx)
display(cons_S_inf)
display(cons_S_inf[("ee", "x²_1")])
display(cons_S_inf[("eo", "x²_1")])
display(cons_S_inf[("oe", "x²_1")])
# display(cons_S_inf[("oo", "x_1")])

println("This is the L(( (Tr(ρ)) - ∑xᵢ²)⋅η), L(( (Tr(ρ)) - ∑yⱼ²)⋅η) ⪰ 0' constraints.")
cons_S₂    = R_constraints.make_loc_cons_S₂(ρ,d,t,Lx)
display(cons_S₂)
display(cons_S₂[("ee", "x")])
display(cons_S₂[("eo", "x")])
display(cons_S₂[("oe", "x")])

println("This is the ' L(( (Tr(ρ)) - ∑xᵢ²)⋅η) ⪰ 0, L(( 1 - ∑yⱼ²)⋅η) = 0' constraints.")
cons_S₂₁x,cons_S₂₁y   = R_constraints.make_loc_cons_S₂₁(ρ,d,t,Lx)
display(cons_S₂₁x)
display(cons_S₂₁y)

println("This is the 'ρ⊗L(η) - L( (xxᵀ⊗yyᵀ) ⊗ η ) ⪰ 0 ∀ η ∈ even-degree-principle-submatrices of ([x,y]₌ₜ₋₂[x,y]₌ₜ₋₂ᵀ)' constraints.")
G_con = R_constraints.make_G_con(ρ,d,t,Lx)
display(G_con["ee"])
## Inspecting the solution
using JuMP
# println("We can recover the moments as with the following command: Dict(zip(Lx.data,value.(Lx.data)))")
# display(Dict(zip(Lx.data,value.(Lx.data))))
#
# display(value.(PSD_con["ee"]))
#
# println("The fourth order constraint should equal ρ")
# display(value.(ord4_con))
#
# println("One of the S_inf constraints is")
# display(value.(cons_S_inf[("oe", "x²_2")]))
#
# println("The G constriant blocks are")
# display(value.(G_con["ee"]))
## The Real model quick
ex = ["Eq22HHH96a" "RANDa3" "p12HK14b" "T1CD12i"][1]
ρ = Examples.get_example(7)
d = (3,3)#filter(:Example => ==(ex),df).Size[1]
t = (2,2)
cl = "S₂"#"S∞ sG" # "S₂ sG" # "S₂₁ sG"
sepr_mod = R_sep_Model.Modelξₜˢᵉᵖ(ρ,d,t;con_list = cl)
sepr_mod = sep_Compute.Computeξₜˢᵉᵖ(sepr_mod)

display(filter(:Example => ==(ex),df))

#Lrx = sepr_mod[:Lx]
## The complex model quick
include(srcDir*"Model\\C_sep_Model.jl")
using .C_sep_Model
ex = ["RANDℂb2" "Eq2-3BP00b" "RANDa3"][3]
ρ = Examples.get_example(ex)
d = filter(:Example => ==(ex),df).Size[1]
t = (2,2)
cl = "S∞sG" # "S₂ sG" # "S₂₁ sG"
sep_mod = C_sep_Model.Modelξₜˢᵉᵖ(ρ,d,t;con_list = cl)
sep_mod = sep_Compute.Computeξₜˢᵉᵖ(sep_mod)
display(filter(:Example => ==(ex),df))

Lx = sep_mod[:Lx]
## Running things in batches
include(srcDir*"batch.jl")
using .batch
t = (2,2)
batch.batch_model(t,Examples.get_examples(),df)
batch.batch_Computeξₜˢᵉᵖ("C:\\Users\\andries\\all-my-codes\\ju-sep-rank\\assets\\bounds\\")
nar = batch.unstack_constraints(df)

## A closer look at the complex constraints
println("""Let us take a deeper look into the constraints:""")
include(srcDir*"Constraints\\C_constraints.jl")
using .C_constraints
#println("These are $(binomial(sum(d)+2*t[1],2*t[1])) the moments present in the hierachy at level t = (2,2)")
mon_expo_keys = C_constraints.make_mon_expo_keys(d,t[1])
display(mon_expo_keys)

println("""These are the blocks of the 'L([x,y]ₜ[x,y]ₜᵀ)'.
 """)

#zero_moms = C_constraints.zeroprop(ρ,d,(3,3),Lx)
PSD_con = C_constraints.make_PSD_con(d,t,Lx)
display(PSD_con)

println("This is 'L(xx*⊗yy*)'.")
ord4_con = C_constraints.make_ord4_con(d,Lx)
display(ord4_con)

println("This is the 'L(( √ρₘₐₓ - xᵢ²)⋅η), L(( √ρₘₐₓ - yⱼ²)⋅η) ⪰ 0 ∀ i ∈ [d₁], j ∈ [d₂]' constraints.")
cons_S_inf = C_constraints.make_loc_cons_S_inf(ρ,d,t,Lx)
display(cons_S_inf)

println("This is the L(( (Tr(ρ)) - ∑xᵢ²)⋅η), L(( (Tr(ρ)) - ∑yⱼ²)⋅η) ⪰ 0' constraints.")
cons_S₂    = C_constraints.make_loc_cons_S₂(ρ,d,t,Lx)
display(cons_S₂)

println("This is the ' L(( (Tr(ρ)) - ∑xᵢ²)⋅η) ⪰ 0, L(( 1 - ∑yⱼ²)⋅η) = 0' constraints.")
cons_S₂₁x,cons_S₂₁y   = C_constraints.make_loc_cons_S₂₁(ρ,d,t,Lx)
display(cons_S₂₁x)
display(cons_S₂₁y)

println("This is the 'ρ⊗L(η) - L( (xxᵀ⊗yyᵀ) ⊗ η ) ⪰ 0 ∀ η ∈ even-degree-principle-submatrices of ([x,y]₌ₜ₋₂[x,y]₌ₜ₋₂ᵀ)' constraints.")
G_con = C_constraints.make_G_con(ρ,d,t,Lx)
display(G_con)

include(srcDir*"Moments\\Moments.jl")
using .Moments
const mom = Moments

expo_conj(d, ααᶥββᶥ) = vcat(ααᶥββᶥ[d[1]+1:2*d[1]],
                            ααᶥββᶥ[1:d[1]],
                            ααᶥββᶥ[1+2*d[1]+d[2]:end],
                            ααᶥββᶥ[1+2*d[1]:2*d[1]+d[2]])

d = (2,2)
t = (2,2)
ρ  =  [1. 0. 0. 0.
       0. 0. 0. 0.
       0. 0. 0. 0.
       0. 0. 0. 1.]

# include(srcDir*"Constraints\\Utils_cons.jl")
# include(srcDir*"\\Constraints\\C_constraints.jl")
# using .Utils_cons
# using .C_constraints
# using LinearAlgebra
#
# CPSD_con = C_constraints.make_PSD_con(d,t,Lx)
# eigmin(value.(CPSD_con["Default"]))
# eigvals(value.(CPSD_con["Default"]))
# Cord4_con = C_constraints.make_ord4_con(d,Lx)
# value.(Cord4_con["real"])
# value.(Cord4_con["imag"])
#
# Cloc_cons_S₂ = C_constraints.make_loc_cons_S₂(ρ,d,t,Lx)
#
# value.(Cloc_cons_S₂[("Default", "x")])
# Cloc_cons_S₂[("Default", "x")][2,2]
#
# for i in [1,2,3], j in 1:8
#     println(i,j)
# end

#
# include(srcDir*"Examples\\Examples.jl")
# include(srcDir *"Examples\\Utils_states.jl")
# using .Examples
# using .Utils_states
# exd     = Examples.get_examples()
# ex = merge(exd["ent"],exd["sep"])
# keys(ex)
# ρ = ex["RANDa3"]
# Export_example_overview = false
# df           = Examples.get_example_overview(Export_example_overview);

# using LinearAlgebra
# function makediagone(ρ)
#     mpli = 1 ./ map(x -> iszero(x) ? 1 : x,  sqrt.([ρ[i,i] for i in 1:size(ρ)[1]]))
#     return diagm(mpli)*ρ*diagm(mpli)
# end
# ρ = [i+j*im for i in 1:4, j in 1:4]
# makediagone(ρ)

max_stars_repo_name: JAndriesJ/sep-rank-bounding

max_issues_repo_name: JAndriesJ/sep-rank-bounding

max_forks_repo_name: JAndriesJ/sep-rank-bounding


text:
function ExponentialDistribution(arg0::jdouble)
    return ExponentialDistribution((jdouble,), arg0)
end

function cumulative_probability(obj::ExponentialDistribution, arg0::jdouble)
    return jcall(obj, "cumulativeProbability", jdouble, (jdouble,), arg0)
end

function density(obj::ExponentialDistribution, arg0::jdouble)
    return jcall(obj, "density", jdouble, (jdouble,), arg0)
end

function get_mean(obj::ExponentialDistribution)
    return jcall(obj, "getMean", jdouble, ())
end

function get_numerical_mean(obj::ExponentialDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::ExponentialDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

function get_support_lower_bound(obj::ExponentialDistribution)
    return jcall(obj, "getSupportLowerBound", jdouble, ())
end

function get_support_upper_bound(obj::ExponentialDistribution)
    return jcall(obj, "getSupportUpperBound", jdouble, ())
end

function inverse_cumulative_probability(obj::ExponentialDistribution, arg0::jdouble)
    return jcall(obj, "inverseCumulativeProbability", jdouble, (jdouble,), arg0)
end

function is_support_connected(obj::ExponentialDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end

function log_density(obj::ExponentialDistribution, arg0::jdouble)
    return jcall(obj, "logDensity", jdouble, (jdouble,), arg0)
end


text:
# import CSV
import CSVFiles
import DataFrames
import FileIO

# import selected names from PredictMD
import ..convert_value_to_missing!
import ..filename_extension
import ..fix_type
import ..is_nothing

"""
"""
function x_contains_y(
        x::AbstractString,
        y::AbstractVector{<:AbstractString},
        )
    if length(y) == 0
        return false
    end
    for i = 1:length(y)
        if occursin(y[i], x)
            return true
        end
    end
    return false
end

"""
"""
function symbol_begins_with(
    x::Symbol,
    y::AbstractString
    )
    if length(y) <= length(string(x)) && string(x)[1:length(y)] == y
        return true
    else
        return false
    end
    return nothing
end

"""
Given a dataframe, return the column names corresponding to CCS "one-hot"
columns.

# Examples

```julia
import CSVFiles
import FileIO
import PredictMD

df = DataFrames.DataFrame(
    FileIO.load(
        MY_CSV_FILE_NAME;
        type_detect_rows = 30_000,
        )
    )

@info(PredictMD.Cleaning.ccs_onehot_names(df))
@info(PredictMD.Cleaning.ccs_onehot_names(df, "ccs_onehot_"))
```
"""
function ccs_onehot_names(
        df::DataFrames.AbstractDataFrame,
        ccs_onehot_prefix::AbstractString = "ccs_onehot_",
        )
    result = column_names_with_prefix(
        df,
        ccs_onehot_prefix,
        )
    return result
end

"""

"""
function column_names_with_prefix(
        df::DataFrames.AbstractDataFrame,
        prefix::AbstractString,
        )
    all_names = DataFrames.names(df)
    name_begins_with_prefix = Vector{Bool}(length(all_names))
    for j = 1:length(all_names)
        name_begins_with_prefix[j] = symbol_begins_with(
            all_names[j],
            prefix,
            )
    end
    vector_of_matching_names = all_names[name_begins_with_prefix]
    return vector_of_matching_names
end

"""
Given a single ICD 9 code, import the relevant patients from the
Health Care Utilization Project (HCUP) National Inpatient Sample (NIS)
database.

# Examples:

```julia
import CSVFiles
import FileIO
import PredictMD

icd_code_list = ["8841"]
icd_code_type=:procedure
input_file_name_list = [
    "./data/nis_2012_core.csv",
    "./data/nis_2013_core.csv",
    "./data/nis_2014_core.csv",
    ]
output_file_name = "./output/hcup_nis_pr_8841.csv"

PredictMD.Cleaning.clean_hcup_nis_csv_icd9(
    icd_code_list,
    input_file_name_list,
    output_file_name;
    icd_code_type=icd_code_type,
    rows_for_type_detect = 30_000,
    )

df = DataFrames.DataFrame(
    FileIO.load(
        output_file_name;
        type_detect_rows = 30_000,
        )
    )

@info(PredictMD.Cleaning.ccs_onehot_names(df))
```
"""
function clean_hcup_nis_csv_icd9(
        icd_code_list::AbstractVector{<:AbstractString},
        input_file_name_list::AbstractVector{<:AbstractString},
        output_file_name::AbstractString;
        header_row::Bool = true,
        print_every_n_lines::Integer = 1_000_000,
        icd_code_type::Union{Nothing, Symbol} = nothing,
        num_dx_columns::Integer = 25,
        num_pr_columns::Integer = 15,
        ccs_onehot_prefix::AbstractString = "ccs_onehot_",
        rows_for_type_detect::Union{Nothing, Integer} = nothing,
        )
    if is_nothing(rows_for_type_detect)
        error("you need to specify rows_for_type_detect")
    end
    if rows_for_type_detect <= 0
        error("rows_for_type_detect must be > 0")
    end

    if is_nothing(icd_code_type)
        error("you need to specify icd_code_type")
    end
    if icd_code_type==:diagnosis
    elseif icd_code_type==:procedure
    else
        error("\"icd_code_type\" must be one of: :diagnosis, :procedure")
    end

    if length(input_file_name_list) == 0
        error("length(input_file_name_list) == 0")
    end

    input_file_name_list = strip.(input_file_name_list)

    for i = 1:length(input_file_name_list)
        if filename_extension(input_file_name_list[i]) != ".csv"
            error("all input files must be .csv")
        end
    end

    output_file_name = strip.(output_file_name)

    if filename_extension(output_file_name) != ".csv"
        error("output file must be .csv")
    end
    if ispath(output_file_name)
        error(
            string(
                "Output file already exists. ",
                "Rename, move, or delete the file, and then try again.",
                )
        )
    end

    temp_file_name_vector = Vector{String}(length(input_file_name_list))
    for i = 1:length(input_file_name_list)
        temp_file_name_vector[i] = string(tempname(), "_", i, ".csv")
    end

    icd_code_list = strip.(icd_code_list)

    for i = 1:length(input_file_name_list)
        if ispath(temp_file_name_vector[i])
            error("ispath(temp_file_name_vector[i])")
        end
        @info(
            string(
                "Starting to read input file ",
                i,
                " of ",
                length(input_file_name_list),
                ".",
                )
            )
        open(input_file_name_list[i], "r") do f_input
            open(temp_file_name_vector[i], "w") do f_temp_output
                line_number = 1
                for line in eachline(f_input)
                    if line_number == 1 && header_row
                        write(f_temp_output, line)
                        write(f_temp_output, "\n")
                    else
                        if x_contains_y(line, icd_code_list)
                            write(f_temp_output, line)
                            write(f_temp_output, "\n")
                        end
                    end
                    line_number += 1
                    if (print_every_n_lines >= 0) &&
                            (line_number % print_every_n_lines == 0)
                        @info(
                            string(
                                "Input file ",
                                i,
                                " of ",
                                length(input_file_name_list),
                                ". Current line number: ",
                                line_number,
                                )
                            )
                    end
                end
            end
        end
        @info(
            string(
                "Finished reading input file ",
                i,
                " of ",
                length(input_file_name_list),
                ".",
                )
            )
    end

    df_vector = Vector{DataFrames.DataFrame}(length(input_file_name_list))

    for i = 1:length(temp_file_name_vector)
        @info(
            string(
                "Starting to read temporary file ",
                i,
                " of ",
                length(input_file_name_list),
                ".",
                )
            )
        # df_i = DataFrames.readtable(temp_file_name_vector[i])
        # We can't use DataFrames.readtable because it is deprecated.
        df_i = DataFrames.DataFrame(
            FileIO.load(
                temp_file_name_vector[i];
                type_detect_rows = rows_for_type_detect,
                )
            )
        df_vector[i] = df_i
        @info(
            string(
                "Finished reading temporary file ",
                i,
                " of ",
                length(input_file_name_list),
                ".",
                )
            )
    end

    for i = 1:length(temp_file_name_vector)
        Base.Filesystem.rm(
            temp_file_name_vector[i];
            force = true,
            recursive = true,
            )
    end

    all_column_names_vectors = [
        DataFrames.names(df) for df in df_vector
        ]
    shared_column_names = intersect(all_column_names_vectors...)

    for i = 1:length(df_vector)
        extra_column_names = setdiff(
            names(df_vector[i]),
            shared_column_names,
            )
        DataFrames.deletecols!(df_vector[i], extra_column_names,)
    end

    combined_df = vcat(df_vector...)

    for i = 1:length(df_vector)
        df_vector[i] = DataFrames.DataFrame()
    end

    if icd_code_type==:diagnosis
        icd_code_column_names = Symbol[
            Symbol( string("DX", j) ) for j = 1:num_dx_columns
            ]
    elseif icd_code_type==:procedure
        icd_code_column_names = Symbol[
            Symbol( string("PR", j) ) for j = 1:num_pr_columns
            ]
    else
        error("\"icd_code_type\" must be one of: :diagnosis, :procedure")
    end

    row_i_has_kth_icd_code_matrix = Matrix{Bool}(
        size(combined_df, 1),
        length(icd_code_list),
        )
    for k = 1:length(icd_code_list)
        current_icd_code = icd_code_list[k]
        row_i_has_current_icd_code_in_col_j_matrix =
            Matrix{Bool}(
                size(combined_df, 1),
                length(icd_code_column_names),
                )
        for j = 1:length(icd_code_column_names)
            @info(
                string(
                    "icd9 code ",
                    k,
                    " of ",
                    length(icd_code_list),
                    ". DX column ",
                    j,
                    " of ",
                    length(icd_code_column_names),
                    ".",
                    )
                )
            for i = 1:size(combined_df, 1)
                cell_value = combined_df[i, icd_code_column_names[j]]
                if DataFrames.ismissing(cell_value)
                    row_i_has_current_icd_code_in_col_j_matrix[i, j] = false
                else
                    cell_value = strip(string(cell_value))
                    row_i_has_current_icd_code_in_col_j_matrix[i, j] =
                        cell_value == current_icd_code
                end
            end
        end
        row_i_has_current_icd_code_in_any_icdcode_column = vec(
            sum(row_i_has_current_icd_code_in_col_j_matrix, 2) .> 0
            )
        row_i_has_kth_icd_code_matrix[:, k] =
            row_i_has_current_icd_code_in_any_icdcode_column
    end

    matching_rows =
        findall(Bool.(vec(sum(row_i_has_kth_icd_code_matrix, 2).>0)))
    num_rows_before = size(combined_df, 1)
    combined_df = combined_df[matching_rows, :]
    num_rows_after = size(combined_df, 1)
    @info(
        string(
            "I initially identified ",
            num_rows_before,
            " rows that could possibly have matched your ICD code(s).",
            " I checked each row, and ",
            num_rows_after,
            " of those rows actually matched your ICD code(s).",
            "I removed the ",
            num_rows_before - num_rows_after,
            " rows that did not match.",
            )
        )

    dx_column_names = [Symbol(string("DX", i)) for i = 1:num_dx_columns]
    dx_ccs_column_names =
        [Symbol(string("DXCCS", i)) for i = 1:num_dx_columns]

    index_to_ccs = strip.(
        string.(
            unique(
                DataFrames.skipmissing(
                    vcat(
                        [combined_df[:, col] for
                            col in dx_ccs_column_names]...
                        )
                    )
                )
            )
        )
    index_to_ccs = index_to_ccs[findall(index_to_ccs .!= "")]
    index_to_ccs = unique(index_to_ccs)
    index_to_ccs = parse.(Int, index_to_ccs)
    sort!(index_to_ccs)
    index_to_ccs = string.(index_to_ccs)

    ccs_to_index = Dict()
    for k = 1:length(index_to_ccs)
        ccs_to_index[  index_to_ccs[ k ]  ] = k
    end
    ccs_to_index = fix_type(ccs_to_index)

    row_i_has_vcode_dx_in_kth_ccs = Matrix{Bool}(
        size(combined_df, 1),
        length(index_to_ccs),
        )

    for j = 1:length(dx_column_names)
        @info(
            string(
                "Processing DXCCS column ",
                j,
                " of ",
                length(dx_column_names),
                ".",
                )
            )
        jth_dx_col_name = dx_column_names[j]
        jth_dx_ccs_col_name = dx_ccs_column_names[j]
        for i = 1:size(combined_df, 1)
            dx_value = combined_df[i, jth_dx_col_name]
            if DataFrames.ismissing(dx_value)
            else
                dx_value = strip(string(dx_value))
                if length(dx_value) == 0
                elseif dx_value[1] == 'V' || dx_value[1] == "V"
                    ccs_value = combined_df[i, jth_dx_ccs_col_name]
                    if DataFrames.ismissing(ccs_value)
                        error(
                            error(
                                "dx value was not missing but",
                                "ccs value was missing"
                                )
                            )
                    else
                        ccs_value = strip(string(ccs_value))
                        row_i_has_vcode_dx_in_kth_ccs[
                            i,
                            ccs_to_index[ccs_value]
                            ] = true
                    end
                end
            end
        end
    end

    for k = 1:length(index_to_ccs)
        kth_ccs = index_to_ccs[k]
        kth_ccs_onehot_column_name = Symbol(
            string(
                ccs_onehot_prefix,
                kth_ccs,
                )
            )
        temporary_column_ints = Int.(row_i_has_vcode_dx_in_kth_ccs[:, k])
        if sum(temporary_column_ints) > 0
            temporary_column_strings = Vector{String}(size(combined_df, 1))
            for i = 1:size(combined_df, 1)
                if temporary_column_ints[i] > 0
                    temporary_column_strings[i] = "Yes"
                else
                    temporary_column_strings[i] = "No"
                end
            end
            combined_df[kth_ccs_onehot_column_name] =
                temporary_column_strings
        else
        end
    end

    convert_value_to_missing!(
        combined_df,
        "A",
        DataFrames.names(combined_df),
        )
    convert_value_to_missing!(
        combined_df,
        "C",
        DataFrames.names(combined_df),
        )
    convert_value_to_missing!(
        combined_df,
        -99,
        DataFrames.names(combined_df),
        )

    try
        mkpath(dirname(output_file_name))
    catch
    end

    @info(string("Attempting to write output file..."))

    FileIO.save(output_file_name, combined_df)

    @info(
        string(
            "Wrote ",
            size(combined_df, 1),
            " rows to output file: \"",
            output_file_name,
            "\"",
            )
        )

    return output_file_name
end


text:
using DataArrays, DataFrames, MixedModels, Base.Test

include(joinpath(dirname(@__FILE__),"data.jl"))

include("paramlowertriangular.jl")
include("scalarReTerm.jl")
include("vectorReTerm.jl")
include("pls.jl")
include("pirls.jl")
include("throws.jl")
include("inject.jl")
include("cfactor.jl")

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
#
# Function that fills the grid value for each point of the read by reading the
# the file containing solute contributions
#

function gmd3D(mysim :: Simulation, solute :: Vector{Int64}, gmd_solute_file :: String;
               dmax = 7.0, step = 0.5, dmin = 1.0)

  # Compute the grid points and get the closest atom to each grid point

  grid = gmd3D_setgrid(mysim,solute,dmax=dmax,step=step,dmin=dmin)

  # Now will read the solute_contribution file

  gmd_solute = readdlm(gmd_solute_file,comments=true,comment_char='#')

  # We need to get the atom indexes from the header of the file, to see to what column each atom
  # corresponds

  icolumn = Vector{Int64}(undef,length(solute))

  legacy = false
  file = open(gmd_solute_file,"r")
  for line in eachline(file)
    if line[1:1] != "#"
      break
    end
    if findfirst("mass:",line) == nothing
      continue
    end
    data = split(line)
    iat = 0
    try
      iat = parse(Int64,data[3])
    catch
      iat = parse(Int64,data[2])
      legacy = true
    end
    icol = parse(Int64,data[2]) + 2
    icolumn[iat] = icol
  end
  close(file)
  if legacy
    println(" Warning: Reading legacy GMD output (previous to 19.340) ")
    println("          If the atoms of the solute are not the first consecutive atoms, you may have problems.")
  end

  # Now, what we need to do is to, for each point in the grid, find the corresponding
  # atom in the gmd_solute table, and interpolate the gmd provided at the distance dmin

  ngrid = length(grid)

  for igrid in 1:ngrid
    # Find values
    icol = icolumn[grid[igrid].atom]
    i = 1
    rup = gmd_solute[i,1]
    dup = gmd_solute[i,icol]
    rdown = 0.
    ddown = 0.
    while rup < grid[igrid].dmin
      rdown = rup
      ddown = dup
      i = i + 1
      rup = gmd_solute[i,1]
      dup = gmd_solute[i,icol]
    end
    # Interpolating
    grid[igrid].rho = ddown + ((dup-ddown)/(rup-rdown))*(grid[igrid].dmin-rdown)
  end

  # Now, returning the immutable array with the final grid

  gmd3D = Vector{GMD3DGrid}(undef,ngrid)
  for i in 1:ngrid
    gmd3D[i] = GMD3DGrid(grid[i].x,grid[i].atom,grid[i].dmin,grid[i].rho)
  end

  return filter( x -> x.rho > 0. , gmd3D )

end


text:
using Test
using Xpress

@testset "$(folder)" for folder in [
    "MathOptInterface", "xprs_callbacks", "Derivative"
]
    @testset "$(file)" for file in readdir(folder)
        include(joinpath(folder, file))
    end
end

@testset "Xpress tests" begin

    prob = Xpress.XpressProblem()

    @test Xpress.getcontrol(prob, "HEURTHREADS") == 0
    @test Xpress.getcontrol(prob, :HEURTHREADS) == 0

    msg = "Unable to call `Xpress.copyprob`:\n\n91 Error: No problem has been input.\n"
    @test_throws Xpress.XpressError(32,msg) Xpress.copyprob(prob, prob)
end


text:
using PolyesterWeave, Aqua
using Test

@testset "PolyesterWeave.jl" begin

  threads, torelease = PolyesterWeave.request_threads(Threads.nthreads()-1)
  @test threads isa NTuple{Int(PolyesterWeave.worker_mask_count()),PolyesterWeave.UnsignedIteratorEarlyStop{PolyesterWeave.worker_type()}}
  @test sum(map(length, threads)) == (PolyesterWeave.num_threads())-1
  map(PolyesterWeave.free_threads!, torelease)

  r1 = PolyesterWeave.request_threads(Sys.CPU_THREADS, 0x0a)
  @test (r1[2][1] & 0x0a) == r1[2][1]
  r2 = PolyesterWeave.request_threads(Sys.CPU_THREADS, 0xff)
  @test (r2[2][1] & 0xff) == r2[2][1]
  @test count_ones(r1[2][1] ⊻ r2[2][1]) ≤ min(8, Sys.CPU_THREADS)
  @test iszero(r1[2][1] & r2[2][1])
  PolyesterWeave.free_threads!(r1[2][1]); PolyesterWeave.free_threads!(r2[2][1])

  @testset "Valid State" begin
    @test sum(map(count_ones, PolyesterWeave.WORKERS[])) == min(512, PolyesterWeave.dynamic_thread_count() - 1)
  end
end
Aqua.test_all(PolyesterWeave)


text:
# In a recent IEEE Trans. on Evolutionary Computation paper by Gary Yen et al
# it is proposed that a double elimination tournament is held between pareto
# approximation fronts in order to rank the EA's producing these fronts.
# Even though the proposed method is interesting and seems to give robust
# comparison results the paper seems to lack a comparison to some base method
# for ranking. Here I want to study this problem in the abstract using
# simulations and statistics. My base hypothesis is that there are simpler
# and more straightforward ways to accomplish this type ranking.

# The basic setting is that of a pool of performance metrics on which a pool
# of solutions should be compared. Each metric calculates a single, real value
# for each front. Thus we can summarize this information in a matrix where
# each row corresponds to a metric and each column to a front that originates
# from one EA. Lets fix some numbers for our simulations:
Nmetrics = 4
NEAs = 5
Nfronts = 10 # 50 fronts per EA

# Without loss of generality, assume that we normalize all the performance metrics
# to the [0.0, 1.0] range and that higher values indicate better performance.
# Thus we can summarize all the info we have in a randomly sampled matrix where
# each column comes from one of the EA's:
rand_performance() = (rand(Nmetrics, NEAs * Nfronts), [rand(1:5) for i in 1:(NEAs*Nfronts)])

# To hold one tournament among N candidates we randomly select one metric and
# then have N/2 binary tournaments, a pool is an array of column indices into an
# original performance matrix
function hold_tournament(P, pool)
  N = length(pool)
  pool_indices = shuffle(collect(1:N))
  winners = zeros(Int64, int(ceil(N/2)))
  losers = zeros(Int64, int(ceil(N/2)))
  for i in 1:int(floor(N/2))
    i1 = pool[pool_indices[2*i-1]]
    i2 = pool[pool_indices[2*i]]
    metric = rand(1:Nmetrics)
    if P[metric, i1] > P[metric, i2]
      winners[i] = i1
      losers[i] = i2
    else
      losers[i] = i1
      winners[i] = i2
    end
  end

  # Add the odd man out to both winners and losers bracket as stated in their paper.
  if mod(N, 2) == 1
    losers[end] = winners[end] = pool[pool_indices[end]]
  end

  return winners, losers
end

# Lets generate one example:
P, EAs = rand_performance()

# And lets run a double elimination tournament to rank the EAs:
function double_elimination_tournament(P, EAs)
  # Initial round
  winners, losers = hold_tournament(P, collect(1:length(EAs)))
  while length(winners) > 1
    winners, wl = hold_tournament(P, winners, metric)
    lw, eliminated = hold_tournament(P, losers, metric)
    losers, eliminated = hold_tournament(P, [wl, lw], metric)
  end
  winners[1]
end

# Lets run many tournaments and collect stats on who wins the most:
function run_many_tournaments(P, EAs, numtournaments = 1000)
  counts = Dict{Any,Int64}()
  for i in 1:numtournaments
    winner = double_elimination_tournament(P, EAs)
    key = (winner, EAs[winner])
    counts[key] = get(counts, key, 0) + 1
  end
  counts
end

function rank_vector(v, rev=false)
  p = sortperm(v, rev = rev)
  ranks = zeros(Int64, size(v))
  for i in 1:length(v)
    ranks[p[i]] = i
  end
  ranks
end

# We can compare these results to methods based on ranks per metric.
function rank_per_row(matrix, rev = false)
  ranks = zeros(Int64, size(matrix))
  for row in 1:size(matrix, 1)
    ranks[row, :] = rank_vector(matrix[row, :][:], rev)
  end
  ranks
end

ranks = rank_per_row(P, true)
meanranks = mean(ranks, 1)[:]


text:
@testset "Structured Matrices" begin
    @testset "Diagonal" begin
        rng, N = MersenneTwister(123456), 3
        rrule_test(Diagonal, randn(rng, N, N), (randn(rng, N), randn(rng, N)))
        D = Diagonal(randn(rng, N))
        rrule_test(Diagonal, D, (randn(rng, N), randn(rng, N)))
        # Concrete type instead of UnionAll
        rrule_test(typeof(D), D, (randn(rng, N), randn(rng, N)))
    end
    @testset "::Diagonal * ::AbstractVector" begin
        rng, N = MersenneTwister(123456), 3
        rrule_test(
            *,
            randn(rng, N),
            (Diagonal(randn(rng, N)), Diagonal(randn(rng, N))),
            (randn(rng, N), randn(rng, N)),
        )
    end
    @testset "diag" begin
        rng, N = MersenneTwister(123456), 7
        rrule_test(diag, randn(rng, N), (randn(rng, N, N), randn(rng, N, N)))
        rrule_test(diag, randn(rng, N), (Diagonal(randn(rng, N)), randn(rng, N, N)))
        rrule_test(diag, randn(rng, N), (randn(rng, N, N), Diagonal(randn(rng, N))))
        rrule_test(diag, randn(rng, N), (Diagonal(randn(rng, N)), Diagonal(randn(rng, N))))
    end
    @testset "Symmetric" begin
        rng, N = MersenneTwister(123456), 3
        rrule_test(Symmetric, randn(rng, N, N), (randn(rng, N, N), randn(rng, N, N)))
    end
    @testset "$f" for f in (Adjoint, adjoint, Transpose, transpose)
        rng = MersenneTwister(32)
        n = 5
        m = 3
        rrule_test(f, randn(rng, m, n), (randn(rng, n, m), randn(rng, n, m)))
        rrule_test(f, randn(rng, 1, n), (randn(rng, n), randn(rng, n)))
    end
    @testset "$T" for T in (UpperTriangular, LowerTriangular)
        rng = MersenneTwister(33)
        n = 5
        rrule_test(T, T(randn(rng, n, n)), (randn(rng, n, n), randn(rng, n, n)))
    end
end


text:
@testset "Indexes" begin
    @testset "General Interface" begin
        @test isdefined(Deductive, :Index)
        # not sure what else should go here tbh. if anything is wrong the downstream tests will catch it
    end

    @testset "Expressions Index" begin
        @symbols a b c x

        exprs1 = Set{AbstractExpression}([
            a,
            b,
            c
        ])
        idx1 = Deductive.index(exprs1)
        @test idx1.current_id == 4
        @test length(keys(idx1.identifier_map)) == 3
        @test typeof(idx1.table) == Deductive.OperatorIndexTable
        @test length(search(idx1, x)) == 3
        @test search(idx1, x) == Set([a, b, c])
        @test length(search(idx1, x ∨ x)) == 0

        exprs2 = Set{AbstractExpression}([
            a ∧ b,
            b → a,
            c ∧ (a → b),
            a ⟷ b
        ])
        idx2 = Deductive.index(exprs2)
        @test idx2.current_id == 5
        @test length(keys(idx2.identifier_map)) == 4
        @test typeof(idx2.table) == Deductive.OperatorIndexTable
        @test length(search(idx2, x)) == 4
        @test length(search(idx2, x ∨ x)) == 0
        @test search(idx2, x ∧ x) == Set([a ∧ b, c ∧ (a → b)])
        @test search(idx2, x → x) == Set([b → a])
        @test search(idx2, x ⟷ x) == Set([a ⟷ b])
    end
end


text:
"""
barychebinit(n)

Get Chebyshev nodes of first kind and corresponding barycentric Lagrange interpolation weights.
Reference: Berrut, J.P. and Trefethen, L.N., 2004. Barycentric lagrange interpolation. SIAM review, 46(3), pp.501-517.

# Arguments
- `n`: order of the Chebyshev interpolation

# Returns
- Chebyshev nodes
- Barycentric Lagrange interpolation weights
"""
function barychebinit(n)
    x = zeros(Float64, n)
    w = similar(x)
    for i in 1:n
        c = (2i - 1)π / (2n)
        x[n - i + 1] = cos(c)
        w[n - i + 1] = (-1)^(i - 1) * sin(c)
    end
    return x, w
end

function vandermonde(x)
    n = length(x)
    vmat = zeros(Float64, (n, n))
    for i in 1:n
        for j in 1:n
            vmat[i,j] = x[i]^(j-1)
        end
    end
    return vmat
end

function invvandermonde(order)
    n = order
    vmat = zeros(Float64, (n, n))
    for i in 1:n
        for j in 1:n
            c = (2n-2i+1)π/(2n)
            x = cos(c)
            vmat[i,j] = x^(j-1)
        end
    end
    return inv(transpose(vmat))
end

@inline function divfactorial(j, i)
    if i == 0
        return 1
    elseif i == 1
        return 1/j
    elseif i == -1
        return j-1
    else
        return factorial(j-1) / factorial(j+i-1)
    end
end

function weightcoef(a, i::Int, n)
    # integrate when i=1; differentiate when i=-1
    b = zeros(Float64, n)
    for j in 1:n
        if j+i-1 > 0
            b[j] = a^(j+i-1) * divfactorial(j, i)
        elseif j+i-1 == 0
            b[j] = 1
        else
            b[j] = 0
        end
    end
    return b
end

@inline function calcweight(invmat, b)
    return invmat*b
end

"""
function barycheb(n, x, f, wc, xc)

Barycentric Lagrange interpolation at Chebyshev nodes
Reference: Berrut, J.P. and Trefethen, L.N., 2004. Barycentric lagrange interpolation. SIAM review, 46(3), pp.501-517.

# Arguments
- `n`: order of the Chebyshev interpolation
- `x`: coordinate to interpolate
- `f`: array of size n, function at the Chebyshev nodes
- `wc`: array of size n, Barycentric Lagrange interpolation weights
- `xc`: array of size n, coordinates of Chebyshev nodes

# Returns
- Interpolation result
"""
function barycheb(n, x, f, wc, xc)
    for j in 1:n
        if x == xc[j]
            return f[j]
        end
    end

    num, den = 0.0, 0.0
    for j in 1:n
        q = wc[j] / (x - xc[j])
        num += q * f[j]
        den += q
    end
    return num / den
end

function barycheb2(n, x, f, wc, xc)
    for j in 1:n
        if x == xc[j]
            return f[j, :]
        end
    end

    den = 0.0
    num = zeros(eltype(f), size(f)[2])
    for j in 1:n
        q = wc[j] / (x - xc[j])
        num += q .* f[j, :]
        den += q
    end
    return num ./ den
end

function chebint(n, a, b, f, invmat)
    wc = weightcoef(b, 1, n) .- weightcoef(a, 1, n)
    intw = calcweight(invmat, wc)
    return sum(intw .* f)
end

function chebdiff(n, x, f, invmat)
    wc = weightcoef(x, -1, n)
    intw = calcweight(invmat, wc)
    return sum(intw .* f)
end


text:
# # Instructions for Assignment06
#
# ## Introduction
#
# Note: this file is runnable in its current state,
# but is incomplete.
# You can run the file from the command line,
# or use the VS Code julia extension to run individual lines.

## To view additional instructions for this assignment, #jl
## see https://wellesley-bisc195.github.io/BISC195.jl/stable/Assignments/assignment06 #jl

## To view a rendered version of this document, #jl
## see https://wellesley-bisc195.github.io/BISC195.jl/stable/Assignments/assignment06_code #jl

# ## Question 1 - Reading FASTA files
#
# A common way to interact with biological sequence data
# is in a format called [`FASTA`](https://en.wikipedia.org/wiki/FASTA_format).
# This is a very simple format, where each "record" (sequence) contains:
#
# 1. A "header" line, always starting with `>`, that contains a description of the sequence
# 2. The sequence, which can be split across multiple lines
#
# For example, in the `data/` directory of this assignment,
# there is a file called `cov2_genomes.fasta`,
# which contains some genome sequences of Sars-CoV2
# isolated from humans early in the pandemic.
# Use the shell to take a look at the first 10 lines (`head`)
# of the file.
#
# Now, let's work on building a julia function
# to help us use this file.

"""
    count_records(path)

Given a path to a `fasta` file,
counts and returns the number of records it contains.

Example
≡≡≡≡≡≡≡≡≡
    julia> count_records("data/cov2_genomes.fasta")
    10
"""
function count_records(path)
    ## Your code here.
end

# Hint - If you're having trouble, the easiest thing to do
# is probably to count lines that start with `>`.

# ## Question 2 - Split the headers
#
# It's nice to know the number of records the file contains,
# but it's more useful to actually be able to read and store
# those records.
#
# Let's look again at the headers for these sequences.
# This can be accomplished from the terminal using `grep`.
# The simplest way to use `grep` is with 2 arguments -
# the first argument specifies a pattern to look for,
# and the second is the path to a file (or group of files).
#
# `grep` will then print to the screen all of the lines
# from the file that match that pattern.
# To see all of the header lines, for example,
#
# ```sh
# $ grep '>' data/cov2_genomes.fasta
# ```
#
# Each header has 3 components, separated by `|`
# Write a function that divides the header
# into these 3 components, and returns them as a tuple.

"""
    fasta_header(header)

Divides a fasta header into its component parts,
removing any leading or trailing spaces.

Example
≡≡≡≡≡≡≡≡≡
    julia> fasta_header(">M0002 |China|Homo sapiens")
    ("M0002", "China", "Homo sapiens")

    julia> fasta_header("AAATTC")
    Error: Invalid header (headers must start with '>')

    julia> fasta_header(">Another sequence")
    ("Another sequence",)

    julia> fasta_header(">headers| can | have| any number | of | info blocks")
    ("headers", "can", "have", "any number", "of", "info blocks")
"""
function fasta_header(header)
    startswith(header, '>') || error("Invalid header (headers must start with '>')")

    ## Your code here
end

# A couple of things to note about this function if you're strugling:
#
# 1. Notice that, even though each header must start with a `>`,
#    the first item in the returned Tuple should not
# 2. The return item should always be a Tuple, even if there's only 1 item
# 3. The `split()` function can divide a `String` by a delimeter.
# 4. You may need to do some searching on the internet.

# ## Question 3 - parse the file
#
# Alright - let's actually parse the file.
# For now, the easiest thing to do
# is to return two vectors,
# one containing the parsed identifiers as tuples,
# and one containing the sequences as strings.
#
# One thing to keep in mind
# is that the sequences in a fasta file
# can be spread out over multiple lines,
# but they'll be much easier to work with as a single string,
# so your function should have a way to join multiple lines
# (remember the `join()` function?).
#
# I've included a few example files in the "data" directory
# with some example files to test your functions on.
# The examples in the docstring below
# assume that your working directory is the assignment directory.

"""
    function parse_fasta(path)

Reads a fasta-formated file and returns 2 vectors,
one containing parsed headers,
the other containing the entire sequence as a `String`.

Note: function does not validate DNA sequences for correctness.

Example
≡≡≡≡≡≡≡≡≡
    julia> ex1 = parse_fasta("data/ex1.fasta");

    julia> ex1 isa Tuple
    true

    julia> ex1[1]
    2-element Array{Tuple{String,String},1}:
      ("ex1.1", "easy")
      ("ex1.2", "multiline")

    julia> ex1[2]
    2-element Array{String,1}:
    "AATTATAGC"
    "CGCCCCCCAGTCGGATT"

    julia> ex2 = parse_fasta("data/ex2.fasta");

    julia> ex2[1]
    4-element Array{Tuple{String,String},1}:
      ("ex2.1", "oneper")
      ("ex2.2", "wrong")
      ("ex2.3", "wronger")
      ("ex2.4", "wrongest")

    julia> ex2[2]
    4-element Array{String,1}:
      "ATCCGT"
      "ATCGTGGaact"
      "ATCGTGGaact"
      "this isn't a dna string,but parse it anyway"
"""
function parse_fasta(path)
    ## Think through the components you need
    ## Does it make sense to define any containers at the beginning?
    ## How will you loop through the file?
    ## What do you need to get from each line?
end

# Tips:
#
# 1. Build the function in parts
#    - first, see if you can build the vector of headers,
#      ignoring the sequences
#    - next, see if you can collect the sequences,
#      even if they're not joined together properly
# 2. Work on pieces of the function in the REPL.
#    Make sure each expression does what you think it does.
#
#    For example, if you make a `for` loop to go through each line of the file,
#    add `@info` statements to check the value of each variable
#    in each loop.

# ## Question 4 - Sequence Statistics
#
# We now have most of the pieces necessary to get some summary statistics
# of CoV2 genomes.
# Note: You may want to grab some code from previous assignments.
#
# Calculate the mean sequence length and standard deviation
# of the CoV2 geneomes in `data/cov2_genomes.fasta`.

## helper code goes here

mean_cov2_length = 0.0
std_cov2_length = 0.0

# Calculate the mean and standard deviation of GC content
# of the CoV2 genomes in `data/cov2_genomes.fasta`.
# GC content should be a number between 0 and 1.
# Note: if you already parsed the fasta file above
# and assigned the results to variables,
# there's no need to re-parse it.

## helper code goes here

mean_cov2_gc = 0.0
std_cov2_gc = 0.0

# ## Bonus Question
#
# Write a function that will read any proper fasta file
# and returns a tuple with the mean and standard deviation
# of the sequence length and gc content.


text:
"""
    CustomDataset

Minimal Custom Dataset.
`CustomDataset` also support indexing and slicing.

### Fields:

* `length::Int`         - Length of dataset.
* `getdata::Function`   - Function which takes `idx` as input and returns corresponding data

### Available Contructor:

* `CustomDataset(length::Int, getdata::Function)`
* `CustomDataset(;length::Int, getdata::Function)`

### Examples:

```jldoctest
julia> x = rand(10,32)
julia> getdata(idx) = x[idx,:]
julia> dset = CustomDataset(size(x,1), getdata)
julia> [x[1,:], x[2,:]] == dset[1:2]
```
"""
struct CustomDataset <: AbstractDataset
    length::Int
    getdata::Function
end

CustomDataset(; length::Int, getdata::Function) = CustomDataset(length, getdata)

Base.size(d::CustomDataset) = (d.length,)
Base.length(d::CustomDataset) = d.length

Base.firstindex(d::CustomDataset) = 1
Base.lastindex(d::CustomDataset) = length(d)

Base.getindex(d::CustomDataset, idx::Int) = d.getdata(idx)
Base.getindex(d::CustomDataset, r::AbstractArray{<:Any,1}) = [d[ri] for ri in r] #TODO: revise this
Base.getindex(d::CustomDataset, c::Colon) = d[1:length(d)]

function Base.show(io::IO, d::CustomDataset)
    print(io, "Custom Dataset:", "\n    length: $(d.length)", "\n    getdata: $(d.getdata)")
end


text:
function spatial_hessian(poly, reference_x, invjac)
    h = hessian(poly, reference_x)
    H11 = h[1] * invjac[1]^2
    H12 = h[2] * invjac[1] * invjac[2]
    H22 = h[3] * invjac[2]^2
    return [
        H11 H12
        H12 H22
    ]
end

function spatial_gradient(poly, reference_x, invjac)
    return vec(gradient(poly, reference_x)) .* invjac
end

function saye_newton_gradient(xk, lk, xq, pk, ∇pk)
    return vcat(xk - xq + lk * ∇pk, pk)
end

function saye_newton_hessian(lk, ∇pk, ∇2pk)
    H = [
        I+lk*∇2pk ∇pk
        ∇pk' 0.0
    ]
    return H
end

function step_saye_newton_iterate(xk, lk, ∇f, ∇2f, steplimiter)
    dim = length(xk)

    δ = ∇2f \ ∇f
    normδx = norm(δ[1:dim])
    if normδx > steplimiter
        δ *= steplimiter / normδx
    end
    xnext = xk - δ[1:dim]
    lnext = lk - δ[end]
    return xnext, lnext
end

function step_chopp_iterate(xk, xq, pk, ∇pk, steplimiter)
    norm2∇pk = dot(∇pk, ∇pk)

    δ1 = -pk / norm2∇pk * ∇pk
    lnext = (xq - xk)' * ∇pk / norm2∇pk

    δ2 = xq - xk - lnext * ∇pk
    normδ2 = norm(δ2)
    if normδ2 > steplimiter
        δ2 *= steplimiter / normδ2
    end
    xnext = xk + δ1 + δ2

    return xnext, lnext
end

function step_saye_iterate(
    xk,
    Xk,
    lk,
    xq,
    func,
    grad,
    hess,
    condtol,
    steplimiter,
)

    dim = length(xk)
    pk = func(Xk)
    ∇pk = grad(Xk)
    ∇2pk = hess(Xk)

    gf = saye_newton_gradient(xk, lk, xq, pk, ∇pk)
    hf = saye_newton_hessian(lk, ∇pk, ∇2pk)

    if inv(cond(hf)) > condtol
        return step_saye_newton_iterate(xk, lk, gf, hf, steplimiter)
    else
        return step_chopp_iterate(xk, xq, pk, ∇pk, steplimiter / 5)
    end
end

function spatial_closest_point(
    xquery,
    xguess,
    func,
    cellmap,
    tol,
    boundingradius,
    maxiter,
    condtol,
)

    invjac = inverse_jacobian(cellmap)
    grad(X) = spatial_gradient(func, X, invjac)
    hess(X) = spatial_hessian(func, X, invjac)

    x0 = copy(xguess)
    X0 = inverse(cellmap, x0)
    g0 = grad(X0)
    l0 = (xquery - x0)' * g0 / (g0' * g0)

    x1 = copy(x0)
    steplimiter = boundingradius / 2
    for iter = 1:maxiter
        X0 = inverse(cellmap, x0)
        x1, l1 = step_saye_iterate(
            x0,
            X0,
            l0,
            xquery,
            func,
            grad,
            hess,
            condtol,
            steplimiter,
        )

        if norm(x1 - xguess) > boundingradius
            return x0, false
        elseif norm(x1 - x0) < tol
            return x1, true
        else
            x0 = x1
            l0 = l1
        end
    end
    return x1, false
end

function closest_points_on_zero_levelset(
    querypoints,
    seedpoints,
    seedcellids,
    levelset,
    tol,
    boundingradius;
    condtol = 1e4eps(),
    maxiter = 20,
)

    mesh = background_mesh(levelset)
    dim, numquerypoints = size(querypoints)
    closestpoints = zeros(dim, numquerypoints)
    closestcellids = zeros(Int, numquerypoints)
    flags = zeros(Bool, numquerypoints)

    tree = KDTree(seedpoints)
    seedidx, seeddists = nn(tree, querypoints)

    for (idx, sidx) in enumerate(seedidx)
        xguess = seedpoints[:, sidx]
        xquery = querypoints[:, idx]
        guesscellid = seedcellids[sidx]

        cellmap = cell_map(mesh, guesscellid)
        load_coefficients!(levelset, guesscellid)
        func = interpolater(levelset)

        xcp, f = spatial_closest_point(
            xquery,
            xguess,
            func,
            cellmap,
            tol,
            boundingradius,
            maxiter,
            condtol,
        )

        closestpoints[:, idx] = xcp
        closestcellids[idx] = guesscellid
        flags[idx] = f
    end
    return closestpoints, closestcellids, flags
end

function spatial_gradient_at_reference_points(levelset, refpoints, refcellids)
    invjac = inverse_jacobian(background_mesh(levelset))
    dim, numpts = size(refpoints)
    spatialgrads = zeros(dim, numpts)
    for (idx, cellid) in enumerate(refcellids)
        load_coefficients!(levelset, cellid)
        poly = interpolater(levelset)
        point = refpoints[:, idx]
        spatialgrads[:,idx] = spatial_gradient(poly, point, invjac)
    end
    return spatialgrads
end

function signed_distance(
    querypoints,
    spatialclosestpoints,
    refclosestpoints,
    cellids,
    levelset,
)

    spatialgrads = spatial_gradient_at_reference_points(
        levelset,
        refclosestpoints,
        cellids,
    )

    difference = querypoints - spatialclosestpoints
    levelsetsign = sign.(vec(sum(spatialgrads .* difference, dims = 1)))

    signeddistance = levelsetsign .* vec(mapslices(norm, difference, dims = 1))

    return signeddistance
end

function distance_to_zero_levelset(
    querypoints,
    seedpoints,
    seedcellids,
    levelset,
    tol,
    boundingradius,
)
    closestpoints, closestcellids = closest_points_on_zero_levelset(
        querypoints,
        seedpoints,
        seedcellids,
        levelset,
        tol,
        boundingradius,
    )
    mesh = background_mesh(levelset)
    refclosestpoints = map_to_reference(closestpoints, closestcellids, mesh)
    signeddistance = signed_distance(
        querypoints,
        closestpoints,
        refclosestpoints,
        closestcellids,
        levelset,
    )

    return signeddistance
end