text:
export
    LineToCenterlineOverlay,
    LineToFrontOverlay,
    CarFollowingStatsOverlay,
    NeighborsOverlay,
    MOBILOverlay,
    MarkerDistOverlay

mutable struct LineToCenterlineOverlay <: SceneOverlay
    target_id::Int # if -1 does it for all
    line_width::Float64
    color::Colorant

    function LineToCenterlineOverlay(target_id::Int=-1;
        line_width::Float64=0.5, #[m]
        color::Colorant=colorant"blue",
        )

        new(target_id, line_width, color)
    end
end
function render!(rendermodel::RenderModel, overlay::LineToCenterlineOverlay, scene::Scene, roadway::Any)

    if overlay.target_id < 0
        target_inds = 1:length(scene)
    else
        target_inds = overlay.target_id:overlay.target_id
    end

    for ind in target_inds
        veh = scene[ind]
        footpoint = get_footpoint(veh)
        add_instruction!(rendermodel, render_line_segment,
            (veh.state.posG.x, veh.state.posG.y, footpoint.x, footpoint.y, overlay.color, overlay.line_width))
    end

    rendermodel
end

mutable struct LineToFrontOverlay <: SceneOverlay
    target_id::Int # if -1 does it for all
    line_width::Float64
    color::Colorant

    function LineToFrontOverlay(target_id::Int=-1;
        line_width::Float64=0.5, #[m]
        color::Colorant=colorant"blue",
        )

        new(target_id, line_width, color)
    end
end
function render!(rendermodel::RenderModel, overlay::LineToFrontOverlay, scene::Scene, roadway::Roadway)

    if overlay.target_id < 0
        target_inds = 1:length(scene)
    else
        target_inds = overlay.target_id:overlay.target_id
    end

    for ind in target_inds
        veh = scene[ind]
        veh_ind_front = get_neighbor_fore_along_lane(scene, ind, roadway).ind
        if veh_ind_front != 0
            v2 = scene[veh_ind_front]
            add_instruction!(rendermodel, render_line_segment,
                (veh.state.posG.x, veh.state.posG.y, v2.state.posG.x, v2.state.posG.y, overlay.color, overlay.line_width))
        end
    end

    rendermodel
end

mutable struct CarFollowingStatsOverlay <: SceneOverlay
    target_id::Int
    verbosity::Int
    color::Colorant
    font_size::Int

    function CarFollowingStatsOverlay(target_id::Int, verbosity::Int=1;
        color::Colorant=colorant"white",
        font_size::Int=10,
        )

        new(target_id, verbosity, color,font_size)
    end
end
function render!(rendermodel::RenderModel, overlay::CarFollowingStatsOverlay, scene::Scene, roadway::Roadway)

    font_size = overlay.font_size
    text_y = font_size
    text_y_jump = round(Int, font_size*1.2)

    add_instruction!( rendermodel, render_text, (@sprintf("id = %d", overlay.target_id), 10, text_y, font_size, overlay.color), incameraframe=false)
        text_y += text_y_jump

    veh_index = findfirst(scene, overlay.target_id)
    if veh_index != 0
        veh = scene[veh_index]

        if overlay.verbosity ≥ 2
            add_instruction!( rendermodel, render_text, ("posG: " * string(veh.state.posG), 10, text_y, font_size, overlay.color), incameraframe=false)
            text_y += text_y_jump
            add_instruction!( rendermodel, render_text, ("posF: " * string(veh.state.posF), 10, text_y, font_size, overlay.color), incameraframe=false)
            text_y += text_y_jump
        end
        add_instruction!( rendermodel, render_text, (@sprintf("speed: %.3f", veh.state.v), 10, text_y, font_size, overlay.color), incameraframe=false)
        text_y += text_y_jump

        foreinfo = get_neighbor_fore_along_lane(scene, veh_index, roadway; max_distance_fore=Inf)
        if foreinfo.ind != 0
            v2 = scene[foreinfo.ind]
            rel_speed = v2.state.v - veh.state.v
            add_instruction!( rendermodel, render_text, (@sprintf("Δv = %10.3f m/s", rel_speed), 10, text_y, font_size, overlay.color), incameraframe=false)
            text_y += text_y_jump
            add_instruction!( rendermodel, render_text, (@sprintf("Δs = %10.3f m/s", foreinfo.Δs), 10, text_y, font_size, overlay.color), incameraframe=false)
            text_y += text_y_jump

            if overlay.verbosity ≥ 2
                add_instruction!( rendermodel, render_text, ("posG: " * string(v2.state.posG), 10, text_y, font_size, overlay.color), incameraframe=false)
                text_y += text_y_jump
                add_instruction!( rendermodel, render_text, ("posF: " * string(v2.state.posF), 10, text_y, font_size, overlay.color), incameraframe=false)
                text_y += text_y_jump
                add_instruction!( rendermodel, render_text, (@sprintf("speed: %.3f", v2.state.v), 10, text_y, font_size, overlay.color), incameraframe=false)
                text_y += text_y_jump
            end
        else
            add_instruction!( rendermodel, render_text, (@sprintf("no front vehicle"), 10, text_y, font_size, overlay.color), incameraframe=false)
        end
    else
        add_instruction!( rendermodel, render_text, (@sprintf("vehicle %d not found", overlay.target_id), 10, text_y, font_size, overlay.color), incameraframe=false)
    end

    rendermodel
end

mutable struct NeighborsOverlay <: SceneOverlay
    target_id::Int
    color_L::Colorant
    color_M::Colorant
    color_R::Colorant
    line_width::Float64
    textparams::TextParams
    function NeighborsOverlay(target_id::Int;
        color_L::Colorant=colorant"blue",
        color_M::Colorant=colorant"green",
        color_R::Colorant=colorant"red",
        line_width::Float64=0.5, # [m]
        textparams::TextParams=TextParams(),
        )

        new(target_id, color_L, color_M, color_R, line_width, textparams)
    end
end
function render!(rendermodel::RenderModel, overlay::NeighborsOverlay, scene::Scene, roadway::Roadway)

    textparams = overlay.textparams
    yₒ = textparams.y_start
    Δy = textparams.y_jump

    vehicle_index = findfirst(scene, overlay.target_id)
    if vehicle_index != 0

        veh_ego = scene[vehicle_index]
        t = veh_ego.state.posF.t
        ϕ = veh_ego.state.posF.ϕ
        v = veh_ego.state.v
        len_ego = veh_ego.def.length
        # show the ego info
        drawtext(@sprintf("ego t:   %10.3f m", t), yₒ + 0*Δy, rendermodel, textparams)
        drawtext(@sprintf("ego ϕ:   %10.3f rad", ϕ), yₒ + 1*Δy, rendermodel, textparams)
        drawtext(@sprintf("ego v:   %10.3f m/s", v), yₒ + 2*Δy, rendermodel, textparams)

        fore_L = get_neighbor_fore_along_left_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront())
        if fore_L.ind != 0
            veh_oth = scene[fore_L.ind]
            A = get_front(veh_ego)
            B = get_rear(veh_oth)
            v_oth = veh_oth.state.v
            Δv = v_oth - v
            add_instruction!(rendermodel, render_line_segment, (A.x, A.y, B.x, B.y, overlay.color_L, overlay.line_width))
            drawtext(@sprintf("d fore left:   %10.3f m", fore_L.Δs), yₒ + 3*Δy, rendermodel, textparams)
            drawtext(@sprintf("Δv fore left:   %10.3f m/s", Δv), yₒ + 4*Δy, rendermodel, textparams)
        end

        fore_M = get_neighbor_fore_along_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront())
        if fore_M.ind != 0
            veh_oth = scene[fore_M.ind]
            A = get_front(veh_ego)
            B = get_rear(veh_oth)
            v_oth = veh_oth.state.v
            Δv = v_oth - v
            add_instruction!(rendermodel, render_line_segment, (A.x, A.y, B.x, B.y, overlay.color_M, overlay.line_width))
            drawtext(@sprintf("d fore middle: %10.3f m", fore_M.Δs), yₒ + 5*Δy, rendermodel, textparams)
            drawtext(@sprintf("Δv fore middle:   %10.3f m/s", Δv), yₒ + 6*Δy, rendermodel, textparams)
        end

        fore_R = get_neighbor_fore_along_right_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront())
        if fore_R.ind != 0
            veh_oth = scene[fore_R.ind]
            A = get_front(veh_ego)
            B = get_rear(veh_oth)
            v_oth = veh_oth.state.v
            Δv = v_oth - v
            add_instruction!(rendermodel, render_line_segment, (A.x, A.y, B.x, B.y, overlay.color_R, overlay.line_width))
            drawtext(@sprintf("d fore right:  %10.3f m", fore_R.Δs), yₒ + 7*Δy, rendermodel, textparams)
            drawtext(@sprintf("Δv fore right:   %10.3f m/s", Δv), yₒ + 8*Δy, rendermodel, textparams)
        end

        rear_L = get_neighbor_rear_along_left_lane(scene, vehicle_index, roadway, VehicleTargetPointRear(), VehicleTargetPointFront(), VehicleTargetPointRear())
        if rear_L.ind != 0
            veh_oth = scene[rear_L.ind]
            A = get_rear(veh_ego)
            B = get_front(veh_oth)
            v_oth = veh_oth.state.v
            Δv = v_oth - v
            add_instruction!(rendermodel, render_line_segment, (A.x, A.y, B.x, B.y, overlay.color_L, overlay.line_width))
            drawtext(@sprintf("d rear left:   %10.3f m", rear_L.Δs), yₒ + 9*Δy, rendermodel, textparams)
            drawtext(@sprintf("Δv rear left:   %10.3f m/s", Δv), yₒ + 10*Δy, rendermodel, textparams)
        end

        rear_M = get_neighbor_rear_along_lane(scene, vehicle_index, roadway, VehicleTargetPointRear(), VehicleTargetPointFront(), VehicleTargetPointRear())
        if rear_M.ind != 0
            veh_oth = scene[rear_M.ind]
            A = get_rear(veh_ego)
            B = get_front(veh_oth)
            v_oth = veh_oth.state.v
            Δv = v_oth - v
            add_instruction!(rendermodel, render_line_segment, (A.x, A.y, B.x, B.y, overlay.color_M, overlay.line_width))
            drawtext(@sprintf("d rear middle: %10.3f m", rear_M.Δs), yₒ + 11*Δy, rendermodel, textparams)
            drawtext(@sprintf("Δv rear middle:   %10.3f m/s", Δv), yₒ + 12*Δy, rendermodel, textparams)
        end

        rear_R = get_neighbor_rear_along_right_lane(scene, vehicle_index, roadway, VehicleTargetPointRear(), VehicleTargetPointFront(), VehicleTargetPointRear())
        if rear_R.ind != 0
            veh_oth = scene[rear_R.ind]
            A = get_rear(veh_ego)
            B = get_front(veh_oth)
            v_oth = veh_oth.state.v
            Δv = v_oth - v
            add_instruction!(rendermodel, render_line_segment, (A.x, A.y, B.x, B.y, overlay.color_R, overlay.line_width))
            drawtext(@sprintf("d rear right:  %10.3f m", rear_R.Δs), yₒ + 13*Δy, rendermodel, textparams)
            drawtext(@sprintf("Δv rear right:   %10.3f m/s", Δv), yₒ + 14*Δy, rendermodel, textparams)
        end
    end

    rendermodel
end

# type MOBILOverlay <: SceneOverlay
#     egoid::Int
#     mobil::MOBIL
#     rec::SceneRecord
#     textparams::TextParams

#     function MOBILOverlay(
#         egoid::Int,
#         mobil::MOBIL;
#         rec::SceneRecord=SceneRecord(1, 0.1),
#         textparams::TextParams = TextParams(x=275, y_start=120),
#         )

#         retval = new()
#         retval.egoid = egoid
#         retval.mobil = mobil
#         retval.rec = rec
#         retval.textparams = textparams
#         retval
#     end
# end
# function render!(rendermodel::RenderModel, overlay::MOBILOverlay, scene::Scene, roadway::Roadway)

#     rec = overlay.rec
#     update!(rec, scene)

#     mobil = overlay.mobil
#     textparams = overlay.textparams
#     yₒ = textparams.y_start

#     vehicle_index = findfirst(rec, overlay.egoid)
#     veh_ego = scene[vehicle_index]
#     v = veh_ego.state.v
#     egostate_M = veh_ego.state

#     drawtext(@sprintf("speed: %10.3f", v), yₒ, rendermodel, textparams)

#     left_lane_exists = convert(Float64, get(N_LANE_LEFT, rec, roadway, vehicle_index)) > 0
#     right_lane_exists = convert(Float64, get(N_LANE_RIGHT, rec, roadway, vehicle_index)) > 0
#     fore_M = get_neighbor_fore_along_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront())
#     rear_M = get_neighbor_rear_along_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointFront(), VehicleTargetPointRear())

#     drawtext(@sprintf("left lane exists: %s", string(left_lane_exists)), yₒ +   textparams.y_jump, rendermodel, textparams)
#     drawtext(@sprintf("right lane exists: %s", string(right_lane_exists)), yₒ + 2*textparams.y_jump, rendermodel, textparams)

#     # accel if we do not make a lane change
#     accel_M_orig = rand(observe!(reset_hidden_state!(mobil.mlon), scene, roadway, overlay.egoid))
#     dir = DIR_MIDDLE

#     drawtext(@sprintf("accel M orig: %10.3f", accel_M_orig), yₒ + 3*textparams.y_jump, rendermodel, textparams)

#     advantage_threshold = mobil.advantage_threshold
#     if right_lane_exists

#         rear_R = get_neighbor_rear_along_right_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointFront(), VehicleTargetPointRear())

#         # candidate position after lane change is over
#         footpoint = get_footpoint(veh_ego)
#         lane = roadway[veh_ego.state.posF.roadind.tag]
#         lane_R = roadway[LaneTag(lane.tag.segment, lane.tag.lane - 1)]
#         roadproj = proj(footpoint, lane_R, roadway)
#         frenet_R = Frenet(RoadIndex(roadproj), roadway)
#         egostate_R = VehicleState(frenet_R, roadway, veh_ego.state.v)

#         Δaccel_n = 0.0
#         passes_safety_criterion = true
#         if rear_R.ind != 0
#             id = scene[rear_R.ind].def.id
#             accel_n_orig = rand(observe!(reset_hidden_state!(mobil.mlon), scene, roadway, id))
#             veh_ego.state = egostate_R
#             accel_n_test = rand(observe!(reset_hidden_state!(mobil.mlon), scene, roadway, id))
#             veh_ego.state = egostate_M
#             passes_safety_criterion = accel_n_test ≥ -mobil.safe_decel
#             Δaccel_n = accel_n_test - accel_n_orig
#         end
#         drawtext(@sprintf("Δaccel n: %10.3f", Δaccel_n), yₒ + 5*textparams.y_jump, rendermodel, textparams)

#         if passes_safety_criterion

#             Δaccel_o = 0.0
#             if rear_M.ind != 0
#                 id = scene[rear_M.ind].def.id
#                 accel_o_orig = rand(observe!(reset_hidden_state!(mobil.mlon), scene, roadway, id))
#                 veh_ego.state = egostate_R
#                 accel_o_test = rand(observe!(reset_hidden_state!(mobil.mlon), scene, roadway, id))
#                 veh_ego.state = egostate_M
#                 Δaccel_o = accel_o_test - accel_o_orig
#             end
#             drawtext(@sprintf("Δaccel o: %10.3f", Δaccel_o), yₒ + 6*textparams.y_jump, rendermodel, textparams)

#             veh_ego.state = egostate_R
#             accel_M_test = rand(observe!(reset_hidden_state!(mobil.mlon), scene, roadway, overlay.egoid))
#             veh_ego.state = egostate_M
#             Δaccel_M = accel_M_test - accel_M_orig

#             drawtext(@sprintf("Δaccel M: %10.3f", Δaccel_M), yₒ + 7*textparams.y_jump, rendermodel, textparams)

#             Δaₜₕ = Δaccel_M + mobil.politeness*(Δaccel_n + Δaccel_o)
#             if Δaₜₕ > advantage_threshold
#                 dir = DIR_RIGHT
#                 advantage_threshold = Δaₜₕ
#             end

#             drawtext(@sprintf("Δaₜₕ: %10.3f", Δaₜₕ), yₒ + 8*textparams.y_jump, rendermodel, textparams)
#             drawtext(@sprintf("advantage_threshold: %10.3f", advantage_threshold), yₒ + 9*textparams.y_jump, rendermodel, textparams)
#         end
#     end

#     if left_lane_exists
#         rear_L = get_neighbor_rear_along_left_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointFront(), VehicleTargetPointRear())

#         # candidate position after lane change is over
#         footpoint = get_footpoint(veh_ego)
#         lane = roadway[veh_ego.state.posF.roadind.tag]
#         lane_L = roadway[LaneTag(lane.tag.segment, lane.tag.lane + 1)]
#         roadproj = proj(footpoint, lane_L, roadway)
#         frenet_L = Frenet(RoadIndex(roadproj), roadway)
#         egostate_L = VehicleState(frenet_L, roadway, veh_ego.state.v)

#         Δaccel_n = 0.0
#         passes_safety_criterion = true
#         if rear_L.ind != 0
#             id = scene[rear_L.ind].def.id
#             accel_n_orig = rand(observe!(reset_hidden_state!(mobil.mlon), scene, roadway, id))
#             veh_ego.state = egostate_L
#             # render!(rendermodel, veh_ego, RGBA(0.0,0.0,1.0,0.5))
#             accel_n_test = rand(observe!(reset_hidden_state!(mobil.mlon), scene, roadway, id))
#             # accel_n_test = rand(track_longitudinal!(reset_hidden_state!(mobil.mlon), scene, roadway, id, vehicle_index))
#             drawtext(@sprintf("accel n test: %10.3f", accel_n_test), yₒ + 10*textparams.y_jump, rendermodel, textparams)

#             body = inertial2body(get_rear(scene[rear_L.ind]), get_front(veh_ego)) # project target to be relative to ego
#             s_gap = body.x
#             drawtext(@sprintf("s_gap L: %10.3f", s_gap), yₒ + 9*textparams.y_jump, rendermodel, textparams)

#             veh_ego.state = egostate_M
#             passes_safety_criterion = accel_n_test ≥ -mobil.safe_decel
#             Δaccel_n = accel_n_test - accel_n_orig
#         end
#         drawtext(@sprintf("Δaccel n: %10.3f", Δaccel_n), yₒ + 11*textparams.y_jump, rendermodel, textparams)

#         if passes_safety_criterion

#             Δaccel_o = 0.0
#             if rear_M.ind != 0
#                 id = scene[rear_M.ind].def.id
#                 accel_o_orig = rand(observe!(reset_hidden_state!(mobil.mlon), scene, roadway, id))
#                 veh_ego.state = egostate_L
#                 accel_o_test = rand(observe!(reset_hidden_state!(mobil.mlon), scene, roadway, id))
#                 veh_ego.state = egostate_M
#                 Δaccel_o = accel_o_test - accel_o_orig
#             end
#             drawtext(@sprintf("Δaccel o: %10.3f", Δaccel_o), yₒ + 12*textparams.y_jump, rendermodel, textparams)

#             veh_ego.state = egostate_L
#             accel_M_test = rand(observe!(reset_hidden_state!(mobil.mlon), scene, roadway, overlay.egoid))
#             veh_ego.state = egostate_M
#             Δaccel_M = accel_M_test - accel_M_orig
#             drawtext(@sprintf("Δaccel M: %10.3f", Δaccel_M), yₒ + 13*textparams.y_jump, rendermodel, textparams)

#             Δaₜₕ = Δaccel_M + mobil.politeness*(Δaccel_n + Δaccel_o)
#             if Δaₜₕ > advantage_threshold
#                 dir = DIR_LEFT
#                 advantage_threshold = Δaₜₕ
#             end

#             drawtext(@sprintf("Δaₜₕ: %10.3f", Δaₜₕ), yₒ + 14*textparams.y_jump, rendermodel, textparams)
#             drawtext(@sprintf("advantage_threshold: %10.3f", advantage_threshold), yₒ + 15*textparams.y_jump, rendermodel, textparams)
#         end
#     end

#     drawtext(@sprintf("dir: %10d", dir), yₒ + 17*textparams.y_jump, rendermodel, textparams)

#     rendermodel
# end

# type CollisionOverlay <: SceneOverlay
#     target_id::Int # if -1 does it for all
#     color::Colorant
#     mem::CPAMemory

#     CollisionOverlay(target_id::Int=-1; color::Colorant=RGBA(1.0,0.0,0.0,0.5)) = new(target_id, color, CPAMemory())
# end
# function render!(rendermodel::RenderModel, overlay::CollisionOverlay, scene::Scene, roadway::Roadway)

#     if overlay.target_id < 0
#         target_inds = 1:length(scene)
#     else
#         ind = findfirst(scene, overlay.target_id)
#         target_inds = ind:ind
#     end

#     for ind in target_inds
#         veh = scene[ind]
#         if get_first_collision(scene, ind, overlay.mem).is_colliding
#             render!(rendermodel, scene[ind], overlay.color)
#         end
#     end

#     rendermodel
# end

mutable struct MarkerDistOverlay <: SceneOverlay
    target_id::Int
    textparams::TextParams
    rec::SceneRecord
    function MarkerDistOverlay(target_id::Int;
        textparams::TextParams=TextParams(),
        )

        new(target_id, textparams, SceneRecord(1, 0.1))
    end
end
function render!(rendermodel::RenderModel, overlay::MarkerDistOverlay, scene::Scene, roadway::Roadway)

    textparams = overlay.textparams
    yₒ = textparams.y_start
    Δy = textparams.y_jump

    update!(overlay.rec, scene)

    vehicle_index = findfirst(scene, overlay.target_id)
    if vehicle_index != 0
        veh_ego = scene[vehicle_index]
        drawtext(@sprintf("lane offset:       %10.3f", veh_ego.state.posF.t), yₒ + 0*Δy, rendermodel, textparams)
        drawtext(@sprintf("markerdist left:   %10.3f", convert(Float64, get(MARKERDIST_LEFT, overlay.rec, roadway, vehicle_index))), yₒ + 1*Δy, rendermodel, textparams)
        drawtext(@sprintf("markerdist right:  %10.3f", convert(Float64, get(MARKERDIST_RIGHT, overlay.rec, roadway, vehicle_index))), yₒ + 2*Δy, rendermodel, textparams)
    end

    rendermodel
end

max_stars_repo_name: LongshengJiang/NGSIM-sim

max_issues_repo_name: LongshengJiang/NGSIM-sim

max_forks_repo_name: LongshengJiang/NGSIM-sim


text:
"""
    bocquetMFhy(lnν, z)

Halo multiplicity function for Δm = 200.

*Reference*: Equation (3, 4) and Table 2 of Bocquet et al., MNRAS, 456, 2361 (2016)
- For the parameters of "M200m" and "Hydro" in Table 2
- `bocquetMFhy(lnν, z)` corresponds to f(σ)/2 in this reference, where σ=√(δc/ν).

# Arguments
- `lnν::Real`: natural logarithm of a threshold, ν, i.e., `lnν` = log(ν), defined by ν ≡ [δc/σ(R,z)]^2. Here, δc = 1.6865 and σ(R,z) is the r.m.s. mass fluctuation within a top-hat smoothing of scale R at a redshift `z`.
- `z::Real`: redshift.

Bocquet et al.'s multiplicity function is NOT normalized, i.e,

``∫_-∞^∞ dlnν bocquetMFhy(lnν, z) ≠ 1``

This function is based on [Cosmology Routine Library (CRL)](https://wwwmpa.mpa-garching.mpg.de/~komatsu/crl/).
"""
function bocquetMFhy(lnν::Real, z::Real)
    A0, a0, b0, c0 = 0.228, 2.15, 1.69, 1.30
    Az, az, bz, cz = 0.285, -0.058, -0.366, -0.045
    A = A0 * (1 + z)^Az
    a = a0 * (1 + z)^az
    b = b0 * (1 + z)^bz
    c = c0 * (1 + z)^cz
    ν = exp(lnν)
    σ = 1.6865 / √ν
    mf = 0.5 * A * ((σ / b)^-a + 1) * exp(-c / σ^2)
end

"""
    bocquetMFdm(lnν, z)

Halo multiplicity function for Δm = 200.

*Reference*: Equation (3, 4) and Table 2 of Bocquet et al., MNRAS, 456, 2361 (2016)
- For the parameters of "M200m" and "DMonly" in Table 2
- `bocquetMFdm(lnν, z)` corresponds to f(σ)/2 in this reference, where σ=√(δc/ν).

# Arguments
- `lnν::Real`: natural logarithm of a threshold, ν, i.e., `lnν` = log(ν), defined by ν ≡ [δc/σ(R,z)]^2. Here, δc = 1.6865 and σ(R,z) is the r.m.s. mass fluctuation within a top-hat smoothing of scale R at a redshift `z`.
- `z::Real`: redshift.

Bocquet et al.'s multiplicity function is NOT normalized, i.e,

``∫_-∞^∞ dlnν bocquetMFdm(lnν, z) ≠ 1``

This function is based on [Cosmology Routine Library (CRL)](https://wwwmpa.mpa-garching.mpg.de/~komatsu/crl/).
"""
function bocquetMFdm(lnν::Real, z::Real)
    A0, a0, b0, c0 = 0.175, 1.53, 2.55, 1.19
    Az, az, bz, cz = -0.012, -0.040, -0.194, -0.021
    A = A0 * (1 + z)^Az
    a = a0 * (1 + z)^az
    b = b0 * (1 + z)^bz
    c = c0 * (1 + z)^cz
    ν = exp(lnν)
    σ = 1.6865 / √ν
    mf = 0.5 * A * ((σ / b)^-a + 1) * exp(-c / σ^2)
end

ext: jl


text:
@testset "deprecations" begin
    x = rand(10)
    f = GP(SqExponentialKernel())
    gp = f(x, 0.1)

    plt = @test_deprecated sampleplot(gp, 10)
    @test plt.n == 1

    @test_deprecated sampleplot!(gp, 4)
    @test plt.n == 2

    @test_deprecated sampleplot!(Plots.current(), gp, 3)
    @test plt.n == 3
end

max_stars_repo_name: JuliaGaussianProcesses/AbstractGP

max_issues_repo_name: JuliaGaussianProcesses/AbstractGP

max_forks_repo_name: JuliaGaussianProcesses/AbstractGP


text:
# Automatically-generated headers for cusolver
#
# DO NOT EDIT THIS FILE DIRECTLY.

# Automatically generated using Clang.jl

const CUSOLVER_VER_MAJOR = 10
const CUSOLVER_VER_MINOR = 3
const CUSOLVER_VER_PATCH = 0
const CUSOLVER_VER_BUILD = 89
const CUSOLVER_VERSION = CUSOLVER_VER_MAJOR * 1000 + CUSOLVER_VER_MINOR * 100 + CUSOLVER_VER_PATCH
const cusolverDnContext = Cvoid
const cusolverDnHandle_t = Ptr{cusolverDnContext}
const syevjInfo = Cvoid
const syevjInfo_t = Ptr{syevjInfo}
const gesvdjInfo = Cvoid
const gesvdjInfo_t = Ptr{gesvdjInfo}
const cusolverDnIRSParams = Cvoid
const cusolverDnIRSParams_t = Ptr{cusolverDnIRSParams}
const cusolverDnIRSInfos = Cvoid
const cusolverDnIRSInfos_t = Ptr{cusolverDnIRSInfos}
const cusolver_int_t = Cint

@cenum cusolverStatus_t::UInt32 begin
    CUSOLVER_STATUS_SUCCESS = 0
    CUSOLVER_STATUS_NOT_INITIALIZED = 1
    CUSOLVER_STATUS_ALLOC_FAILED = 2
    CUSOLVER_STATUS_INVALID_VALUE = 3
    CUSOLVER_STATUS_ARCH_MISMATCH = 4
    CUSOLVER_STATUS_MAPPING_ERROR = 5
    CUSOLVER_STATUS_EXECUTION_FAILED = 6
    CUSOLVER_STATUS_INTERNAL_ERROR = 7
    CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED = 8
    CUSOLVER_STATUS_NOT_SUPPORTED = 9
    CUSOLVER_STATUS_ZERO_PIVOT = 10
    CUSOLVER_STATUS_INVALID_LICENSE = 11
    CUSOLVER_STATUS_IRS_PARAMS_NOT_INITIALIZED = 12
    CUSOLVER_STATUS_IRS_PARAMS_INVALID = 13
    CUSOLVER_STATUS_IRS_INTERNAL_ERROR = 14
    CUSOLVER_STATUS_IRS_NOT_SUPPORTED = 15
    CUSOLVER_STATUS_IRS_OUT_OF_RANGE = 16
    CUSOLVER_STATUS_IRS_NRHS_NOT_SUPPORTED_FOR_REFINE_GMRES = 17
    CUSOLVER_STATUS_IRS_INFOS_NOT_INITIALIZED = 18
end

@cenum cusolverEigType_t::UInt32 begin
    CUSOLVER_EIG_TYPE_1 = 1
    CUSOLVER_EIG_TYPE_2 = 2
    CUSOLVER_EIG_TYPE_3 = 3
end

@cenum cusolverEigMode_t::UInt32 begin
    CUSOLVER_EIG_MODE_NOVECTOR = 0
    CUSOLVER_EIG_MODE_VECTOR = 1
end

@cenum cusolverEigRange_t::UInt32 begin
    CUSOLVER_EIG_RANGE_ALL = 1001
    CUSOLVER_EIG_RANGE_I = 1002
    CUSOLVER_EIG_RANGE_V = 1003
end

@cenum cusolverNorm_t::UInt32 begin
    CUSOLVER_INF_NORM = 104
    CUSOLVER_MAX_NORM = 105
    CUSOLVER_ONE_NORM = 106
    CUSOLVER_FRO_NORM = 107
end

@cenum cusolverIRSRefinement_t::UInt32 begin
    CUSOLVER_IRS_REFINE_NOT_SET = 1100
    CUSOLVER_IRS_REFINE_NONE = 1101
    CUSOLVER_IRS_REFINE_CLASSICAL = 1102
    CUSOLVER_IRS_REFINE_CLASSICAL_GMRES = 1103
    CUSOLVER_IRS_REFINE_GMRES = 1104
    CUSOLVER_IRS_REFINE_GMRES_GMRES = 1105
    CUSOLVER_PREC_DD = 1300
    CUSOLVER_PREC_SS = 1301
    CUSOLVER_PREC_SHT = 1302
end

const cusolverSpContext = Cvoid
const cusolverSpHandle_t = Ptr{cusolverSpContext}
const csrqrInfo = Cvoid
const csrqrInfo_t = Ptr{csrqrInfo}


text:
function RHD2164(port::AbstractString)
    if port=="PortA1"
        ports=[PortA1,PortA1Ddr]
    elseif port=="PortA2"
        ports=[PortA2,PortA2Ddr]
    elseif port=="PortB1"
        ports=[PortB1,PortB1Ddr]
    elseif port=="PortB2"
        ports=[PortB2,PortB2Ddr]
    elseif port=="PortC1"
        ports=[PortC1,PortC1Ddr]
    elseif port=="PortC2"
        ports=[PortC2,PortC2Ddr]
    elseif port=="PortD1"
        ports=[PortD1,PortD1Ddr]
    elseif port=="PortD2"
        ports=[PortD2,PortD2Ddr]
    else
        ports=[0,0]
    end
end

function RHD2132(port::AbstractString)
    if port=="PortA1"
        ports=[PortA1]
    elseif port=="PortA2"
        ports=[PortA2]
    elseif port=="PortB1"
        ports=[PortB1]
    elseif port=="PortB2"
        ports=[PortB2]
    elseif port=="PortC1"
        ports=[PortC1]
    elseif port=="PortC2"
        ports=[PortC2]
    elseif port=="PortD1"
        ports=[PortD1]
    elseif port=="PortD2"
        ports=[PortD2]
    else
        ports=[0]
    end
end

getlibrary(fpgas::Array{FPGA,1})=map(getlibrary,fpgas)

function getlibrary(fpga::FPGA)
    fpga.lib=Libdl.dlopen(intan_lib,Libdl.RTLD_NOW)
    fpga.board=ccall(Libdl.dlsym(fpga.lib,:okFrontPanel_Construct), Ptr{Void}, ())
    nothing
end

function init_board!(rhd::RHD2000,fpga::DArray{FPGA,1,Array{FPGA,1}})
    if rhd.debug.state==false
        map(open_board,fpga)
        map(uploadFpgaBitfile,fpga)
    else
        map(getlibrary,fpga)
    end

    @sync for p in procs(fpga)
        @async remotecall_fetch((d,sr,db)->init_board_helper(localpart(d),sr,db),p,fpga,rhd.sr,rhd.debug.state)
    end
    nothing
end

function init_board!(rhd::RHD2000,fpga::Array{FPGA,1})
    if rhd.debug.state==false
        map(open_board,fpga)
        map(uploadFpgaBitfile,fpga)
    else
        for thisfpga in fpga
            getlibrary(thisfpga)
        end
    end

    init_board_helper(fpga,rhd.sr,rhd.debug.state)

    nothing
end

function init_board_helper(fpgas::Array{FPGA,1},sr,mydebug=false)
    for fpga in fpgas
        init_board_helper(fpga,sr,mydebug)
    end
    nothing
end

function init_board_helper(fpga::FPGA,sr,mydebug=false)

    initialize_board(fpga,mydebug)

    #For 64 channel need two data streams, and data will come in
    #on the rising AND falling edges of SCLK
    if fpga.amps[1] == 255
        enableDataStream(fpga,0, false)
    else
        stream=0
        for i in fpga.amps
            enableDataStream(fpga,stream,true)
            if (OPEN_EPHYS)
                setDataSource(fpga,stream,i)
            end
            stream+=1
        end
    end

    #Enable DAC
    enableDac(fpga,0,true)

    calculateDataBlockSizeInWords(fpga)
    calculateDataBlockSizeInBytes(fpga)

    setSampleRate(fpga,sr,mydebug)
    println("Sample Rate set at ",fpga.sampleRate, " on board ", fpga.id)

    ledArray=[1,0,0,0,0,0,0,0]
    if (OPEN_EPHYS)
        setLedDisplay(fpga,ledArray)
    end

    #Set up an RHD2000 register object using this sample rate to optimize MUX-related register settings.
    fpga.r=CreateRHD2000Registers(Float64(fpga.sampleRate))

    if EMG
        setLowerBandwidth(2.0, fpga.r)
        setUpperBandwidth(1000.0, fpga.r)
        setDspCutoffFreq(10.0, fpga.r)
    end
    #Upload version with no ADC calibration to AuxCmd3 RAM Bank 0.
    commandList=createCommandListRegisterConfig(zeros(Int32,1),false,fpga.r)
    uploadCommandList(fpga,commandList, "AuxCmd3", 0)

    #Upload version with ADC calibration to AuxCmd3 RAM Bank 1.
    commandList=createCommandListRegisterConfig(zeros(Int32,1),true,fpga.r)
    uploadCommandList(fpga,commandList, "AuxCmd3", 1)

    selectAuxCommandLength(fpga,"AuxCmd3", 0, length(commandList) - 1)

    if mydebug==false
        for port in ["PortA","PortB","PortC","PortD"]
            if check_port_streams(fpga,port)>0
                determine_delay(fpga,port)
                selectAuxCommandBank(fpga,port, "AuxCmd3", 1)
            end
        end
    end

    setMaxTimeStep(fpga,SAMPLES_PER_DATA_BLOCK)
    setContinuousRunMode(fpga,false)

    if mydebug==false
        runBoard(fpga)
        while (isRunning(fpga))
        end
        flushBoard(fpga)
    end

   for port in ["PortA","PortB","PortC","PortD"]
        if check_port_streams(fpga,port)>0
            selectAuxCommandBank(fpga,port, "AuxCmd3", 0)
        end
    end
    setContinuousRunMode(fpga,true)

    nothing
end

function open_board(fpga::FPGA)

    fpga.lib=Libdl.dlopen(intan_lib,Libdl.RTLD_NOW)
    fpga.board=ccall(Libdl.dlsym(fpga.lib,:okFrontPanel_Construct), Ptr{Void}, ())

    println("Scanning USB for Opal Kelly devices...")
    nDevices=ccall(Libdl.dlsym(fpga.lib,:okFrontPanel_GetDeviceCount), Int, (Ptr{Void},), fpga.board)
    println("Found ", nDevices, " Opal Kelly device(s)")

    #Get Serial Number
    if VERSION > v"0.7-"
        serial=Array{UInt8}(undef,11)
    else
        serial=Array{UInt8}(11)
    end
    ccall(Libdl.dlsym(fpga.lib,:okFrontPanel_GetDeviceListSerial), Int32, (Ptr{Void}, Int, Ptr{UInt8}), fpga.board, fpga.id,serial)
    serial[end]=0
    serialnumber=unsafe_string(pointer(serial))
    println("Serial number of device 0 is ", serialnumber)

    #Open by serial
    if (ccall(Libdl.dlsym(fpga.lib,:okFrontPanel_OpenBySerial), Cint, (Ptr{Void},Ptr{UInt8}),fpga.board,serialnumber)!=0)
        println("Device could not be opened. Is one connected?")
        return -2
    end

    #configure on-board PLL
    ccall(Libdl.dlsym(fpga.lib,:okFrontPanel_LoadDefaultPLLConfiguration), Cint, (Ptr{Void},),fpga.board)

    nothing
end

function open_board(fpgas::Array{FPGA,1})
    for fpga in fpgas
        open_board(fpga)
    end
    nothing
end

function uploadFpgaBitfile(fpgas::Array{FPGA,1})
    for fpga in fpgas
        uploadFpgaBitfile(fpga)
    end
    nothing
end

function uploadFpgaBitfile(rhd::FPGA)

    #upload configuration file
    if rhd.usb3
        if (OPEN_EPHYS)
            errorcode=ccall(Libdl.dlsym(rhd.lib,:okFrontPanel_ConfigureFPGA),Cint,(Ptr{Void},Ptr{UInt8}),rhd.board,usb3bit_open_ephys)
        else
            errorcode=ccall(Libdl.dlsym(rhd.lib,:okFrontPanel_ConfigureFPGA),Cint,(Ptr{Void},Ptr{UInt8}),rhd.board,usb3bit)
        end
    else
        errorcode=ccall(Libdl.dlsym(rhd.lib,:okFrontPanel_ConfigureFPGA),Cint,(Ptr{Void},Ptr{UInt8}),rhd.board,bit)
    end
    if errorcode==0
        println("FPGA configuration loaded.")
    else
        println("FPGA configuration failed.")
    end

    #Check if FrontPanel Support is enabled
    myenable = ccall(Libdl.dlsym(rhd.lib,:okFrontPanel_IsFrontPanelEnabled),Bool,(Ptr{Void},),rhd.board)

    UpdateWireOuts(rhd)
    boardId = GetWireOutValue(rhd,WireOutBoardId)
    boardVersion = GetWireOutValue(rhd,WireOutBoardVersion)

    if rhd.usb3
        if (OPEN_EPHYS)
            target_id = RHYTHM_BOARD_ID_OPENEPHYS
        else
	       target_id = 700
       end
    else
        target_id = RHYTHM_BOARD_ID
    end

    if (boardId != target_id)
        println("FPGA configuration does not support Rythm. Incorrect board ID: ", boardId)
    else
        println("Rhythm configuration file successfully loaded. Rhythm version number: ", boardVersion)
    end

    nothing
end

function initialize_board(fpgas::Array{FPGA,1})
    for fpga in fpgas
        initialize_board(fpga)
    end
    nothing
end

function initialize_board(rhd::FPGA,debug=false)

    resetBoard(rhd)

    if (!OPEN_EPHYS)
        readDigitalInManual(rhd)
    end

    setSampleRate(rhd,30000,debug)
    selectAuxCommandBank(rhd,"PortA", "AuxCmd1", 0)
    selectAuxCommandBank(rhd,"PortB", "AuxCmd1", 0)
    selectAuxCommandBank(rhd,"PortC", "AuxCmd1", 0)
    selectAuxCommandBank(rhd,"PortD", "AuxCmd1", 0)
    selectAuxCommandBank(rhd,"PortA", "AuxCmd2", 0)
    selectAuxCommandBank(rhd,"PortB", "AuxCmd2", 0)
    selectAuxCommandBank(rhd,"PortC", "AuxCmd2", 0)
    selectAuxCommandBank(rhd,"PortD", "AuxCmd2", 0)
    selectAuxCommandBank(rhd,"PortA", "AuxCmd3", 0)
    selectAuxCommandBank(rhd,"PortB", "AuxCmd3", 0)
    selectAuxCommandBank(rhd,"PortC", "AuxCmd3", 0)
    selectAuxCommandBank(rhd,"PortD", "AuxCmd3", 0)
    selectAuxCommandLength(rhd,"AuxCmd1", 0, 0)
    selectAuxCommandLength(rhd,"AuxCmd2", 0, 0)
    selectAuxCommandLength(rhd,"AuxCmd3", 0, 0)

    setContinuousRunMode(rhd,true)

    setMaxTimeStep(rhd,4294967295) #4294967395 == (2^32 - 1)

    setCableLengthFeet(rhd,"PortA", 6.0)  # assume 6 ft cables
    setCableLengthFeet(rhd,"PortB", 6.0)
    setCableLengthFeet(rhd,"PortC", 6.0)
    setCableLengthFeet(rhd,"PortD", 6.0)

    setDspSettle(rhd,false)

    if (OPEN_EPHYS)
        setDataSource(rhd,0, PortA1)
        setDataSource(rhd,1, PortB1)
        setDataSource(rhd,2, PortC1)
        setDataSource(rhd,3, PortD1)
        setDataSource(rhd,4, PortA2)
        setDataSource(rhd,5, PortB2)
        setDataSource(rhd,6, PortC2)
        setDataSource(rhd,7, PortD2)

        if rhd.usb3
            setDataSource(rhd,8, PortA1)
            setDataSource(rhd,9, PortB1)
            setDataSource(rhd,10, PortC1)
            setDataSource(rhd,11, PortD1)
            setDataSource(rhd,12, PortA2)
            setDataSource(rhd,13, PortB2)
            setDataSource(rhd,14, PortC2)
            setDataSource(rhd,15, PortD2)
        end
    else
        SetWireInValue(rhd,WireInDataStreamEn, 0x00000000);
        UpdateWireIns(rhd);
    end

    #remember that julia indexes with 1's instead of 0's to start an array
    enableDataStream(rhd,0, true)
    for i=1:(MAX_NUM_DATA_STREAMS-1)
        enableDataStream(rhd,i,false)
    end

    #clearTtlOut(rhd)

    for i=0:7; enableDac(rhd,i, false); end
    for i=0:7; selectDacDataStream(rhd,i, 0); end
    for i=0:7; selectDacDataChannel(rhd,i, 0); end
    #for i=0:7; selectDacDataStream(rhd,i, 8); end
    #selectDacDataChannel(rhd,0, 0)
    #selectDacDataChannel(rhd,1, 1)
    #for i=2:7; selectDacDataChannel(rhd,i, 0); end

    setDacManual(rhd,32768)    # midrange value = 0 V

    setDacGain(rhd,0)
    setAudioNoiseSuppress(rhd,0)

    setTtlMode(rhd,falses(8))

    #Loop through digital output sequencers
    for i=1:16
        rhd.d[i].channel=i
        update_digital_output(rhd,rhd.d[i])
    end

    for i=0:7; setDacThreshold(rhd,i, 32768, true); end

    enableExternalFastSettle(rhd,false)
    setExternalFastSettleChannel(rhd,15)

    enableExternalDigOut(rhd,"PortA", false)
    enableExternalDigOut(rhd,"PortB", false)
    enableExternalDigOut(rhd,"PortC", false)
    enableExternalDigOut(rhd,"PortD", false)
    setExternalDigOutChannel(rhd,"PortA", 0)
    setExternalDigOutChannel(rhd,"PortB", 0)
    setExternalDigOutChannel(rhd,"PortC", 0)
    setExternalDigOutChannel(rhd,"PortD", 0)

    if (OPEN_EPHYS)
        enableBoardLeds(rhd,true);
    end

    nothing
end

function resetBoard(rhd::FPGA)

    SetWireInValue(rhd,WireInResetRun, 0x01, 0x01)
    UpdateWireIns(rhd)
    SetWireInValue(rhd,WireInResetRun, 0x00, 0x01)
    UpdateWireIns(rhd)

    if rhd.usb3
       SetWireInValue(rhd,WireInMultiUse, div(USB3_BLOCK_SIZE,4))
       UpdateWireIns(rhd)
       if (OPEN_EPHYS)
           ActivateTriggerIn(rhd,TrigInOpenEphys,16)
       else
           ActivateTriggerIn(rhd,TrigInConfig,9)
       end

       SetWireInValue(rhd,WireInMultiUse, DDR_BURST_LENGTH)
       UpdateWireIns(rhd)
       if (OPEN_EPHYS)
           ActivateTriggerIn(rhd,TrigInOpenEphys,17)
       else
           ActivateTriggerIn(rhd,TrigInConfig,10)
       end
    end

    nothing
end

function setSampleRate(rhd::FPGA,newSampleRate::Int64,debug=false)

    if newSampleRate==1000
        M=7
        D=125
    elseif newSampleRate==1250
        M=7
        D=100
    elseif newSampleRate==1500
        M=21
        D=250
    elseif newSampleRate==2000
        M=14
        D=125
    elseif newSampleRate==2500
        M=35
        D=250
    elseif newSampleRate==3000
        M=21
        D=125
    elseif newSampleRate==3333
        M=14
        D=75
    elseif newSampleRate==4000
        M=28
        D=125
    elseif newSampleRate==5000
        M=7
        D=25
    elseif newSampleRate==6250
        M=7
        D=20
    elseif newSampleRate==8000
        M=112
        D=250
    elseif newSampleRate==10000
        M=14
        D=25
    elseif newSampleRate==12500
        M=7
        D=10
    elseif newSampleRate==15000
        M=21
        D=25
    elseif newSampleRate==20000
        M=28
        D=25
    elseif newSampleRate==25000
        M=35
        D=25
    elseif newSampleRate==30000
        M=42
        D=25
    else
    end

    rhd.sampleRate=newSampleRate

    #Wait for DcmProgDone==1 before reprogramming clock synthesizer
    if !debug
        while !isDcmProgDone(rhd)
        end
    end

    #Reprogram clock synthesizer

    SetWireInValue(rhd,WireInDataFreqPll,(256 * convert(Culong,M) + convert(Culong,D)))
    UpdateWireIns(rhd)
    if (OPEN_EPHYS)
        ActivateTriggerIn(rhd,TrigInDcmProg,0)
    else
        ActivateTriggerIn(rhd,TrigInConfig,0)
    end

    #Wait for DataClkLocked = 1 before allowing data acquisition to continue
    if !debug
        while !isDataClockLocked(rhd)
        end
    end

    nothing
end

function isDcmProgDone(rhd::FPGA)

    UpdateWireOuts(rhd)
    value=GetWireOutValue(rhd,WireOutDataClkLocked)
    return ((value & 0x0002) > 1)
end

function isDataClockLocked(rhd::FPGA)

    UpdateWireOuts(rhd)
    value=GetWireOutValue(rhd,WireOutDataClkLocked)
    return ((value & 0x0001) > 0)
end

function setContinuousRunMode(rhd::FPGA,continuousMode)

    if continuousMode
        SetWireInValue(rhd,WireInResetRun,0x02,0x02)
    else
        SetWireInValue(rhd,WireInResetRun,0x00,0x02)
    end

    UpdateWireIns(rhd)

    nothing
end

function setMaxTimeStep(rhd::FPGA,maxTimeStep)

    if (OPEN_EPHYS)
        maxTimeStep=convert(UInt32, maxTimeStep)

        maxTimeStepLsb = maxTimeStep & 0x0000ffff
        maxTimeStepMsb = maxTimeStep & 0xffff0000

        SetWireInValue(rhd,WireInMaxTimeStepLsb,maxTimeStepLsb)
        SetWireInValue(rhd,WireInMaxTimeStepMsb,(maxTimeStepMsb >> 16))
    else
        SetWireInValue(rhd,WireInMaxTimeStep, maxTimeStep);
    end
    UpdateWireIns(rhd)

    nothing
end

setDspSettle(rhd::FPGA,enabled)=(SetWireInValue(rhd,WireInResetRun,(enabled ? 0x04 : 0x00),0x04);UpdateWireIns(rhd))

function setDataSource(rhd::FPGA,stream, dataSource)

    #error checking goes here

    if stream==0
        endPoint=WireInDataStreamSel1234
        bitShift=0
    elseif stream==1
        endPoint=WireInDataStreamSel1234
        bitShift=4
    elseif stream==2
        endPoint=WireInDataStreamSel1234
        bitShift=8
    elseif stream==3
        endPoint=WireInDataStreamSel1234
        bitShift=12
    elseif stream==4
        endPoint=WireInDataStreamSel5678
        bitShift=0
    elseif stream==5
        endPoint=WireInDataStreamSel5678
        bitShift=4
    elseif stream==6
        endPoint=WireInDataStreamSel5678
        bitShift=8
    elseif stream==7
        endPoint=WireInDataStreamSel5678
        bitShift=12
elseif stream==8
	endPoint=WireInDataStreamSel1234
	bitShift=16
elseif stream==9
	endPoint=WireInDataStreamSel1234
	bitShift=20
elseif stream==10
	endPoint=WireInDataStreamSel1234
	bitShift=24
elseif stream==11
	endPoint=WireInDataStreamSel1234
	bitShift=28
elseif stream==12
	endPoint=WireInDataStreamSel5678
	bitShift=16
elseif stream==13
	endPoint=WireInDataStreamSel5678
	bitShift=20
elseif stream==14
	endPoint=WireInDataStreamSel5678
	bitShift=24
elseif stream==15
	endPoint=WireInDataStreamSel5678
	bitShift=28
    end

    bitShift=convert(Int32,bitShift)
    SetWireInValue(rhd,endPoint,(dataSource << bitShift), (0x000f << bitShift))
    UpdateWireIns(rhd)

    nothing
end

function enableDataStream(rhd::FPGA,stream::Int, enabled::Bool)

    #error checking goes here

    stream=convert(Int32,stream)
    if enabled
        if rhd.dataStreamEnabled[stream+1] == 0
            SetWireInValue(rhd,WireInDataStreamEn,0x0001 << stream, 0x0001 << stream)
            UpdateWireIns(rhd)
            rhd.dataStreamEnabled[stream+1] = 1;
            rhd.numDataStreams=rhd.numDataStreams+1;
        end
    else
        if rhd.dataStreamEnabled[stream+1] == 1
            SetWireInValue(rhd,WireInDataStreamEn,0x0000 << stream, 0x0001 << stream)
            UpdateWireIns(rhd)
            rhd.dataStreamEnabled[stream+1] = 0;
            rhd.numDataStreams=rhd.numDataStreams-1;
        end
    end

    nothing
end

function enableExternalDigOut(rhd::FPGA,port, enable)

    SetWireInValue(rhd,WireInMultiUse, (enable ? 1 : 0))
    UpdateWireIns(rhd)

    if port=="PortA"
        if (OPEN_EPHYS)
            ActivateTriggerIn(rhd,TrigInExtDigOut,0)
        else
            ActivateTriggerIn(rhd,TrigInConfig,16)
        end
    elseif port=="PortB"
        if (OPEN_EPHYS)
            ActivateTriggerIn(rhd,TrigInExtDigOut,1)
        else
            ActivateTriggerIn(rhd,TrigInConfig,17)
        end
    elseif port=="PortC"
        if (OPEN_EPHYS)
            ActivateTriggerIn(rhd,TrigInExtDigOut,2)
        else
            ActivateTriggerIn(rhd,TrigInConfig,18)
        end
    elseif port=="PortD"
        if (OPEN_EPHYS)
            ActivateTriggerIn(rhd,TrigInExtDigOut,3)
        else
            ActivateTriggerIn(rhd,TrigInConfig,19)
        end
    end

    nothing
end

function setExternalDigOutChannel(rhd::FPGA,port, channel)

    SetWireInValue(rhd,WireInMultiUse,channel)
    UpdateWireIns(rhd)

    if port=="PortA"
        if (OPEN_EPHYS)
            ActivateTriggerIn(rhd,TrigInExtDigOut,4)
        else
            ActivateTriggerIn(rhd,TrigInConfig,24)
        end
    elseif port=="PortB"
        if (OPEN_EPHYS)
            ActivateTriggerIn(rhd,TrigInExtDigOut,5)
        else
            ActivateTriggerIn(rhd,TrigInConfig,25)
        end
    elseif port=="PortC"
        if (OPEN_EPHYS)
            ActivateTriggerIn(rhd,TrigInExtDigOut,6)
        else
            ActivateTriggerIn(rhd,TrigInConfig,26)
        end
    elseif port=="PortD"
        if (OPEN_EPHYS)
            ActivateTriggerIn(rhd,TrigInExtDigOut,7)
        else
            ActivateTriggerIn(rhd,TrigInConfig,27)
        end
    end

    nothing
end

function setLedDisplay(rhd::FPGA,ledArray)

    ledOut=0
    for i=1:8
        if ledArray[i]>0
            ledOut += 1 << (i-1)
        end
    end

    SetWireInValue(rhd,WireInLedDisplay_openephys,ledOut)
    UpdateWireIns(rhd)

    nothing
end

runBoard(fpgas::Array{FPGA,1})=map(runBoard,fpgas)

runBoard(rhd::FPGA)=ActivateTriggerIn(rhd,TrigInSpiStart,0)

function isRunning(rhd::FPGA)

    UpdateWireOuts(rhd)
    value=GetWireOutValue(rhd,WireOutSpiRunning)

    if ((value & 0x01) == 0)
        return false
    else
        return true
    end
end

function flushBoard(rhd::FPGA)

    if !rhd.usb3
        while (numWordsInFifo(rhd) >= (USB_BUFFER_SIZE/2))
            ReadFromPipeOut(rhd,PipeOutData, USB_BUFFER_SIZE, rhd.usbBuffer)
        end

        while (numWordsInFifo(rhd) > 0)
            ReadFromPipeOut(rhd,PipeOutData, (2 * numWordsInFifo(rhd)), rhd.usbBuffer)
        end
    else
        SetWireInValue(rhd,WireInResetRun, 1<<16, 1<<16)
        UpdateWireIns(rhd)
        if (OPEN_EPHYS)
            while (numWordsInFifo(rhd) > 0)
                ReadFromBlockPipeOut(rhd,PipeOutData, USB3_BLOCK_SIZE, rhd.usbBuffer)
            end
        else
            while (numWordsInFifo(rhd) > USB_BUFFER_SIZE / 2)
                ReadFromBlockPipeOut(rhd,PipeOutData,USB_BUFFER_SIZE,rhd.usbBuffer)
            end
            while (numWordsInFifo(rhd) > 0)
                ReadFromBlockPipeOut(rhd,PipeOutData,USB3_BLOCK_SIZE * max(div(2*numWordsInFifo(rhd) , USB3_BLOCK_SIZE),1),rhd.usbBuffer)
            end
        end
        SetWireInValue(rhd,WireInResetRun,0<<16, 1<<16)
        UpdateWireIns(rhd)
    end
    nothing
end

function numWordsInFifo(rhd::FPGA)

    UpdateWireOuts(rhd)

    if (OPEN_EPHYS)
        GetWireOutValue(rhd,WireOutNumWordsMsb)<<16+GetWireOutValue(rhd,WireOutNumWordsLsb)
    else
        GetWireOutValue(rhd,WireOutNumWords)
    end
end

function compareNumWords(fpgas::Array{FPGA,1})
    out=false
    for fpga in fpgas
        out=out|compareNumWords(fpga)
    end
    out
end

compareNumWords(fpga::FPGA)=numWordsInFifo(fpga) < fpga.numWords

function calibrate_parallel(fpga,s,v,buf,nums,mytime,calnum)

    @sync begin
        for p in procs(fpga)
            @async remotecall_wait((ff,ss)->readDataBlocks_cal(localpart(ff),localpart(ss),v,buf,nums,mytime),p,fpga,s)
        end
    end
    nothing
end

function onlinesort_parallel(fpga,s,v,buf,nums,mytime)

    @sync begin
        for p in procs(fpga)
            @async remotecall_wait((ff,ss)->readDataBlocks_on(localpart(ff),localpart(ss),v,buf,nums,mytime),p,fpga,s)
        end
    end
    nothing
end

function readDataBlocks_cal(fpgas::Array{FPGA,1},s,v,buf,nums,mytime,calnum)

    readDataBlocks_cal(fpgas[1],s,v,buf,nums,mytime,calnum)
    nothing
end

function readDataBlocks_cal(fpga::FPGA,s,v,buf,nums,mytime,calnum)

    #block until there are enough words
    while compareNumWords(fpga)
    end

    if fpga.usb3
       ReadFromBlockPipeOut(fpga,PipeOutData,convert(Clong,fpga.numBytesPerBlock),fpga.usbBuffer)
    else
	   ReadFromPipeOut(fpga,PipeOutData,convert(Clong,fpga.numBytesPerBlock),fpga.usbBuffer)
    end
    fillFromUsbBuffer!(fpga,0,v,mytime)

    #Reference Channels

    cal!(s,v,buf,nums,calnum)

    nothing
end

function readDataBlocks_on(fpgas::Array{FPGA,1},s,v,buf,nums,mytime)

    readDataBlocks_on(fpgas[1],s,v,buf,nums,mytime)
    nothing
end

function readDataBlocks_on(fpga::FPGA,s,v,buf,nums,mytime)

    #block until there are enough words
    while compareNumWords(fpga)
    end

    if fpga.usb3
       ReadFromBlockPipeOut(fpga,PipeOutData,convert(Clong,fpga.numBytesPerBlock),fpga.usbBuffer)
    else
	ReadFromPipeOut(fpga,PipeOutData,convert(Clong,fpga.numBytesPerBlock),fpga.usbBuffer)
    end

    fillFromUsbBuffer!(fpga,0,v,mytime)

    #Reference Channels

    onlinesort!(s,v,buf,nums)

    nothing
end

function cal_update(rhd::RHD2000)

    if rhd.cal==0

        rhd.cal=1

    elseif rhd.cal<3

        if rhd.reads>20
            #rhd.cal=2
            rhd.cal=3

        end

    elseif rhd.cal==3

    end

    rhd.reads+=1

    nothing
end

function readDataBlocks(rhd::RHD2000,numBlocks::Int64,s,myfpga::Array{FPGA,1})

    if compareNumWords(myfpga)
        return false
    end

    #=
    if (numBytesToRead > USB_BUFFER_SIZE)
        println("USB buffer size exceeded")
        return false
    end
    =#

    numRead=0

    numBytesToRead = myfpga[1].numBytesPerBlock * numBlocks
    if length(myfpga)>1
        for fpga in myfpga
            ReadFromPipeOut(fpga,PipeOutData, convert(Clong, fpga.numBytesPerBlock * numBlocks), fpga.usbBuffer)
        end
        numRead=numBytesToRead
    else
        if myfpga[1].usb3
            numRead=ReadFromBlockPipeOut(myfpga[1],PipeOutData, convert(Clong, numBytesToRead), myfpga[1].usbBuffer)
        else
            numRead=ReadFromPipeOut(myfpga[1],PipeOutData, convert(Clong, numBytesToRead), myfpga[1].usbBuffer)
        end
    end

    for i=0:(numBlocks-1)

        #Move data from usbBuffer to v

        fillFromUsbBuffer!(myfpga,i,rhd.v,rhd.time)

        #If analog here, we want to move adc number x into voltage channels
        #
        if myfpga[1].amps[1] == 255
            if (OPEN_EPHYS)
                try
                    rhd.v[:,1] = myfpga[1].adc[:,2] .- 2556 #open ephys special
                catch
                    for jj=1:size(rhd.v,1)
                        if abs(myfpga[1].adc[jj,2] - 2556)>typemax(Int16)
                        else
                            rhd.v[jj,1] = myfpga[1].adc[jj,2] - 2556
                        end
                    end
                end
            else
                try
                    rhd.v[:,1] = myfpga[1].adc[:,2]
                catch
                    for jj=1:size(rhd.v,1)
                        if abs(myfpga[1].adc[jj,2])>typemax(Int16)
                            rhd.v[jj,1] = typemax(Int16) * sign(myfpga[1].adc[jj,2])
                        else
                            rhd.v[jj,1] = myfpga[1].adc[jj,2]
                        end
                    end
                end
            end
        end

        for j=1:size(rhd.v,2)
            if rhd.refs[j]>0
                for k=1:size(rhd.v,1)
                    rhd.v[k,j]=rhd.v[k,j]-rhd.v[k,rhd.refs[j]]
                end
            end
        end

        #Filtering
        for i=1:length(rhd.filts)
            for j=1:length(rhd.filts[i])
                apply_filter(rhd,rhd.filts[i][j],rhd.filts[i][j].chan)
            end
        end

        #Sorting
        applySorting(rhd,s)
    end

    return true
end

function applySorting(rhd::RHD2000,s)

    if rhd.cal==0

        cal!(s,rhd.v,rhd.buf,rhd.nums,rhd.cal)
        rhd.cal=1

    elseif rhd.cal<3

        cal!(s,rhd.v,rhd.buf,rhd.nums,rhd.cal)

        if rhd.reads>20
            #rhd.cal=2
            rhd.cal=3
        end

    elseif rhd.cal==3

        onlinesort!(s,rhd.v,rhd.buf,rhd.nums)
    end

    rhd.reads+=1
    nothing
end

function enableBoardLeds(rhd::FPGA,enable)

    SetWireInValue(rhd,WireInMultiUse, enable ? 1 : 0)
	UpdateWireIns(rhd)
    ActivateTriggerIn(rhd,TrigInOpenEphys, 0)

end

function setClockDivider(rhd,divide_factor)

    SetWireInValue(rhd,WireInMultiUse, divide_factor);
	UpdateWireIns(rhd);
    ActivateTriggerIn(rhd,TrigInOpenEphys, 1);

end

function calculateDataBlockSizeInWords(rhd::FPGA)
    #rhd.numWords = SAMPLES_PER_DATA_BLOCK * (4+2+(rhd.numDataStreams*36)+8+2)

    rhd.numWords = SAMPLES_PER_DATA_BLOCK * (4+2+(rhd.numDataStreams*36) + 8+2)

    nothing
    #4 = magic number; 2 = time stamp; 36 = (32 amp channels + 3 aux commands + 1 filler word); 8 = ADCs; 2 = TTL in/out
end

function calculateDataBlockSizeInBytes(rhd::FPGA)
    rhd.numBytesPerBlock=2 * rhd.numWords
    nothing
end

function checkUsbHeader(usbBuffer,index)

    @inbounds x1 = convert(UInt64,usbBuffer[index])
    @inbounds x2 = convert(UInt64,usbBuffer[index+1])
    @inbounds x3 = convert(UInt64,usbBuffer[index+2])
    @inbounds x4 = convert(UInt64,usbBuffer[index+3])
    @inbounds x5 = convert(UInt64,usbBuffer[index+4])
    @inbounds x6 = convert(UInt64,usbBuffer[index+5])
    @inbounds x7 = convert(UInt64,usbBuffer[index+6])
    @inbounds x8 = convert(UInt64,usbBuffer[index+7])

    header = (x8 << 56) + (x7 << 48) + (x6 << 40) + (x5 << 32) + (x4 << 24) + (x3 << 16) + (x2 << 8) + (x1 << 0)

    return (header == RHD2000_HEADER_MAGIC_NUMBER)|(header==RHD2000_HEADER_MAGIC_NUMBER_OPEN_EPHYS)
end

function fillFromUsbBuffer!(fpgas::Array{FPGA,1},blockIndex::Int64,v,mytime)

    for fpga in fpgas
	fillFromUsbBuffer!(fpga,blockIndex,v,mytime)
    end
    nothing
end

function fillFromUsbBuffer!(fpga::FPGA,blockIndex::Int64,v,mytime)

    index = blockIndex * fpga.numBytesPerBlock + 1

    for t=1:SAMPLES_PER_DATA_BLOCK

	#Header
        if !checkUsbHeader(fpga.usbBuffer,index)
            if t==1 #first header missing, we're fucked. creep ahead and see what you find
                #Add data

                lag=0
                newindex=index

                for i=2:fpga.numBytesPerBlock
                    if checkUsbHeader(fpga.usbBuffer,i)
                        lag=i-index
                        newindex=i
                        break
                    end
                end

            else #something is messed up, probably in the last sample. start incrementing backward until header is found
                lag=0
                newindex=index
                for i=(index-1):-1:1
                    if checkUsbHeader(fpga.usbBuffer,i) #lag specifies number of bytes that are missing
                        lag=index-i
                        newindex=i
                        break
                    end
                end

                if lag>=(fpga.numBytesPerBlock/SAMPLES_PER_DATA_BLOCK) #somehow an extra byte popped up? will want to move forward instead
		lag=0
                newindex=index

                for i=(index+1):fpga.numBytesPerBlock
                    if checkUsbHeader(fpga.usbBuffer,i)
                        lag=i-index
                        newindex=i
                        break
                    end
                end

		end
            end

            #get extra bytes

            while (2*numWordsInFifo(fpga) < lag)
            end

            temp_array=zeros(UInt8,lag)

            ReadFromPipeOut(fpga,PipeOutData, convert(Clong, lag), temp_array)

	    if t==1 #moved forward, so bytes should be added to the end
	       count=1
	       for i=(fpga.numBytesPerBlock+1):(fpga.numBytesPerBlock+lag)
	       	   fpga.usbBuffer[i]=temp_array[count]
		   count+=1
	       end
	    else #moved backward, so skip everything in bad block, then start
	    	 count=1
	       for i=(fpga.numBytesPerBlock+1):(fpga.numBytesPerBlock+lag)
	       	   fpga.usbBuffer[i]=temp_array[count]
		   count+=1
	       end
	       t=t-1 #refill last block
	    end

            #start fresh
            index=newindex
        end

	index+=8
	mytime[t,fpga.id]=convertUsbTimeStamp(fpga.usbBuffer,index)
	index+=4

	#Auxiliary results
	index += (2*3*fpga.numDataStreams)

	#Amplifier
	for i=1:32
	    for j=1:fpga.numDataStreams
		@inbounds v[t,32*(j-1)+i+fpga.shift]=convertUsbWord(fpga.usbBuffer,index)
		index+=2
	    end
	end

	#skip filler word(s)
        if fpga.numDataStreams>0
	       index += 2
       end

	#ADCs
        for i=1:8
            @inbounds fpga.adc[t,i]=convertUsbWord(fpga.usbBuffer,index)
            index+=2
        end

	#TTL in
        @inbounds fpga.ttlin[t]=convertUsbWordu(fpga.usbBuffer,index)
        index += 2

        #TTL out
        @inbounds fpga.ttlout[t]=convertUsbWordu(fpga.usbBuffer,index)
	index += 2
    end

    nothing
end

function queueToFile(rhd::RHD2000,task::Task,fpga)

    if rhd.save.record_mode
        if rhd.save.save_full
            f=open(rhd.save.v, "a+")
            write(f,rhd.v)
            close(f)
        end
        if rhd.save.ttl_s
            #TTL in
            f=open(rhd.save.ttl,"a+")
            write(f,fpga[1].ttlin)
            close(f)

            #TTL out
            f=open(rhd.save.ttl_out,"a+")
            write(f,fpga[1].ttlout)
            close(f)
        end
        if rhd.save.lfp_s
            f=open(rhd.save.lfp,"a+")
            write(f,rhd.lfps)
            close(f)
        end
        if rhd.save.adc_s
            f=open(rhd.save.adc,"a+")
            write(f,fpga[1].adc)
            close(f)
        end

        if rhd.save.ts_s
            writeTimeStamp(rhd)
        end

        save_task(task,rhd)
    end
    #Clear buffers
    @inbounds for i=1:size(rhd.buf,2)
        for j=1:rhd.nums[i]
            rhd.buf[j,i]=Spike()
        end
        rhd.nums[i]=0
    end
    nothing
end

function writeTimeStamp(rhd::RHD2000)

    #write spike times and cluster identity

    f=open(rhd.save.ts, "a+")

    for i=1:size(rhd.time,2)
        write(f,rhd.time[1,i])
    end

    @inbounds for i::UInt16=1:size(rhd.v,2)
        write(f,i) #channel number (UInt16)
        write(f,rhd.nums[i]) #number of spikes coming up (UInt16)
        for j=1:rhd.nums[i]
            write(f,rhd.buf[j,i].inds.start) #Index of start
            write(f,rhd.buf[j,i].inds.stop) #Index of stop
            write(f,rhd.buf[j,i].id) # cluster number (UInt8)
        end
    end

    close(f)

    nothing
end

function convertUsbTimeStamp(usbBuffer, index::Int64)

    @inbounds x1 = convert(UInt32,usbBuffer[index])
    @inbounds x2 = convert(UInt32,usbBuffer[index+1])
    @inbounds x3 = convert(UInt32,usbBuffer[index+2])
    @inbounds x4 = convert(UInt32,usbBuffer[index+3])

    convert(UInt32,((x4<<24) + (x3<<16) + (x2<<8) + (x1<<0)))
end

function convertUsbWord(usbBuffer, index::Int64)

    @inbounds x1 = convert(UInt16,usbBuffer[index])
    @inbounds x2 = convert(UInt16,usbBuffer[index+1])

    convert(Int16,signed((x2<<8)|x1)-typemax(Int16)) #is this right?
end

function convertUsbWordu(usbBuffer,index::Int64)

    @inbounds x1 = convert(UInt16,usbBuffer[index])
    @inbounds x2 = convert(UInt16,usbBuffer[index+1])

    (x2<<8)|x1
end

ReadUsbBuffer(fpgas::Array{FPGA,1})=map(ReadUsbBuffer,fpgas)

function check_port_streams(fpga::FPGA,port)

    if port=="PortA"
        streams=sum((fpga.amps.==0).|(fpga.amps.==1).|(fpga.amps.==8).|(fpga.amps.==9))
    elseif port=="PortB"
        streams=sum((fpga.amps.==2).|(fpga.amps.==3).|(fpga.amps.==10).|(fpga.amps.==11))
    elseif port=="PortC"
        streams=sum((fpga.amps.==4).|(fpga.amps.==5).|(fpga.amps.==12).|(fpga.amps.==13))
    elseif port=="PortD"
        streams=sum((fpga.amps.==6).|(fpga.amps.==7).|(fpga.amps.==14).|(fpga.amps.==15))
    end

    streams
end

function Intan_GUI(myconfig_path=string(base_path,"Intan_config.jl"))

    include(myconfig_path)

end


text:
using IntUtils

function main()
    s = 0
    for i in 1:9999
        if is_amicable(i)
            @printf("%d is amicable\n", i)
            s += i
        end
    end
    @printf("The sum of amicable numbers from 1 to 9999 is %d\n", s)
end

main()


text:
# AUTO GENERATED FILE - DO NOT EDIT

export keywords

"""
    keywords(;kwargs...)

A Keywords component.

Keyword arguments:
- `id` (String; optional)
- `wrapperStyle` (Dict; optional)
- `tagStyle` (Dict; optional)
- `inputStyle` (Dict; optional)
- `tagDeleteStyle` (Dict; optional)
- `placeholder` (String; optional)
- `value` (Array; optional)
"""
function keywords(; kwargs...)
        available_props = Symbol[:id, :wrapperStyle, :tagStyle, :inputStyle, :tagDeleteStyle, :placeholder, :value]
        wild_props = Symbol[]
        return Component("keywords", "Keywords", "dash_cool_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: raulrozza/dash-cool-components

max_issues_repo_name: raulrozza/dash-cool-components

max_forks_repo_name: raulrozza/dash-cool-components


text:
"""
Google Secret Manager API
"""
module _secrets

export secrets

using ..api
using ...root

"""
Google Secret Manager API root.

Secret JSON representation: https://cloud.google.com/secret-manager/docs/reference/rest/v1/projects.secrets#Secret
{
  "name": string,
  "replication": {
    object (Replication)
  },
  "createTime": string,
  "labels": {
    string: string,
    ...
  },

  // Union field expiration can be only one of the following:
  "expireTime": string,
  "ttl": string
  // End of list of possible types for union field expiration.
}

"""
secrets = APIRoot(
    "https://secretmanager.googleapis.com/v1/projects/{project}",
    Dict(
        "cloud-platform" => "Full access to all resources and services in the specified Cloud Platform project.",
    );
    Secrets=APIResource("secrets";  # https://cloud.google.com/secret-manager/docs/reference/rest/v1/projects.secrets/
        list=APIMethod(:GET, "", "Lists matching Secrets"),
        get=APIMethod(:GET, "{secret}", "Gets metadata for a given Secret"),
        create=APIMethod(:POST, "", "Creates a new Secret containing no SecretVersions"),
        addVersion=APIMethod(:POST, "{secret}:addVersion", "Creates a new SecretVersion containing secret data and attaches to an existing Secret"),
        delete=APIMethod(:DELETE, "{secret}", "Deletes a Secret"),
    ),
    Versions=APIResource("secrets/{secret}/versions";  # https://cloud.google.com/secret-manager/docs/reference/rest/v1/projects.secrets.versions/
        list=APIMethod(:GET, "", "Lists SecretVersions.  This call does not return secret data"),
        get=APIMethod(:GET, "{version}", "Gets metadata for a SecretVersion."),
        disable=APIMethod(:POST, "{version}:disable", "Disables a SecretVersion."),
        enable=APIMethod(:POST, "{version}:enable", "Enables a SecretVersion."),
        destroy=APIMethod(:POST, "{version}:destroy", "Destroys a SecretVersion."),
        access=APIMethod(:GET, "{version}:access", "Accesses a SecretVersion. This call returns the secret data."),
    ),
)

end


text:
module MangalTestBackbone
    using Mangal
    using Test
    include("setup.jl")

    # Backbone only
    @test typeof(backbones()) <: Vector{MangalReferenceTaxon}

    @test typeof(backbones("q" => "Salix")) <: Vector{MangalReferenceTaxon}

    @test backbone(1).name == "Abutilon theophrasti"
    @test backbone("Abutilon theophrasti").id == 1

    @test length(backbones(Pair("count", 10))) == 10

end


text:
x=5
y = x+1
println(y)
# if isinteractive()
#     println(y)
# end

# String Operations
# Concatenation: *
# Repetition: ^
s1 = "throat"
s2 = "warbler"
println(s1 * s2)
println(s1^2 * s2^2)


text:
nn = 20
tot = collect(1:nn)
trnidx = view(tot, [1,3,6,9,10,11,12,13,14,15,19])

res


text:
const TrackedComponentArray{V, D, N, DA, A, Ax} = ReverseDiff.TrackedArray{V,D,N,ComponentArray{V,N,A,Ax},DA}

maybe_tracked_array(val::AbstractArray, der, tape, inds, origin) = ReverseDiff.TrackedArray(val, der, tape)
function maybe_tracked_array(val::Real, der, tape, inds, origin::AbstractVector)
    ax = getaxes(ReverseDiff.value(origin))[1]
    i = ax[inds[1]].idx
    return ReverseDiff.TrackedReal(val, der, tape, i, origin)
end

for f in [:getindex, :view]
    @eval function Base.$f(tca::TrackedComponentArray, inds::Union{Symbol, Val}...)
            val = $f(ReverseDiff.value(tca), inds...)
            der = Base.maybeview(ReverseDiff.deriv(tca), inds...)
            t = ReverseDiff.tape(tca)
        return maybe_tracked_array(val, der, t, inds, tca)
    end
end

function Base.getproperty(tca::TrackedComponentArray, s::Symbol)
    if s in (:value, :deriv, :tape)
        return getfield(tca, s)
    else
        val = getproperty(ReverseDiff.value(tca), s)
        der = getproperty(ReverseDiff.deriv(tca), s)
        t = ReverseDiff.tape(tca)
        return maybe_tracked_array(val, der, t, (s,), tca)
    end
end

max_stars_repo_name: jonniedie/ComponentArrays

max_issues_repo_name: jonniedie/ComponentArrays

max_forks_repo_name: jonniedie/ComponentArrays


text:
using InteratomicPotentials

using AtomsBase
using LinearAlgebra
using StaticArrays
using Test
using Unitful
using UnitfulAtomic

include("mocks.jl")

@testset "Potentials.jl" begin
    # The first set of tests are simple unit tests that generally test just one module at a time.
    # These tests collectively have high code coverage and test whether each module is structurally correct.
    # They do not validate the numeric results for scientific accuracy.
    @testset "Unit Tests" begin
        include("unit/constants.jl")
        include("unit/unit_convention.jl")
        include("unit/nnlist.jl")
        @testset "Default Implementations" begin
            include("unit/types/abstract_potential.jl")
            include("unit/types/non_trainable_potential.jl")
            include("unit/types/trainable_potential.jl")
            include("unit/types/empirical_potential.jl")
        end
        @testset "Empirical Potentials" begin
            include("unit/empirical/lj.jl")
            include("unit/empirical/bm.jl")
            include("unit/empirical/coulomb.jl")
            include("unit/empirical/zbl.jl")
            include("unit/empirical/morse.jl")
        end
        @testset "Mixed Potentials" begin
            include("unit/mixed/linear_combination_potential.jl")
        end
    end

    # The second set of tests are integration tests which incorporate multiple modules.
    # These tests validate the numeric results of calcuations to ensure that the modules are logically correct.
    # In other words, they test that the full system works together and produces scientifically accurate results.
    @testset "Integration Tests" begin
        # TODO: Dallas
        include("integration/lj_clusters/lj_150.jl")
        include("integration/lj_clusters/lj_1000.jl")

    end
end


text:
# This file was generated, do not modify it. # hide
y = data.MedV
X = select(data, Not(:MedV));


text:
const TOP_INSTITUTION_FIELDS = [:name, :fisyr, :netassets, :charitytype,
  :categorynamefiltered]

const TOP_INSTITUTION_METHODS = [:rank, :quantile, :cumulative]

const TOP_INSTITUTIONS_YEAR = 2015
const TOP_INSTITUTIONS_QUANTILE = 0.95
const TOP_INSTITUTIONS_CUMULATIVE = 0.95
const TOP_INSTITUTIONS_RANK = 500
const TOP_INSTITUTIONS_PREFIX = "top_"
const WORST_MIN_ASSETS = 100_000_000
const WORST_NUM_YEARS = 5

processmissing(v::AbstractVector) = collect(skipmissing(v))

#writes out the institutions in the subdataframe
function writeTopInstitutions(subdf::SubDataFrame, label::String;
  outputPath::String = TABLE_PATH,
  topInstitutionFields::Vector{Symbol} = TOP_INSTITUTION_FIELDS)::Nothing

  outPath = "$outputPath\\$(label).csv"
  CSV.write(outPath,subdf[topInstitutionFields])

  return nothing
end

getTopInstitutionsRank(subdf::SubDataFrame,
  topInstitutionsRank::Int)::SubDataFrame = view(subdf, 1:topInstitutionsRank, :)

function getTopInstitutionsQuantile(subdf::SubDataFrame,
  topInstitutionsQuantile::Float64)
  v::Vector{Float64} = subdf[:netassets]
  qFunc = StatsBase.ecdf(v)
  dfquantile::Vector{Float64} = (qFunc).(subdf[:netassets])

  return view(subdf, dfquantile .> topInstitutionsQuantile, :)
end

function getTopInstitutionsCumulative(subdf::SubDataFrame,
  topInstitutionsCumulative::Float64)

  v::Vector{Float64} = subdf[:netassets]
  cumv::Vector{Float64} = similar(v) #will hold the running total
  cumv[1] = v[1]
  for i::Int ∈ 2:length(v)
    cumv[i] = cumv[i-1] + v[i]
  end
  cumv .= cumv ./ cumv[end] #normalize

  return view(subdf, cumv .< topInstitutionsCumulative, :)
end

#gets the top institutions given a method and SORTED subdataframe
function getTopInstitutions(subdf::SubDataFrame, topMethod::Symbol;
  topInstitutionsYear::Int = TOP_INSTITUTIONS_YEAR,
  topInstitutionsQuantile::Float64 = TOP_INSTITUTIONS_QUANTILE,
  topInstitutionsRank::Int = TOP_INSTITUTIONS_RANK,
  topInstitutionsCumulative::Float64 = TOP_INSTITUTIONS_CUMULATIVE,
  topInstitutionsPrefix::String = TOP_INSTITUTIONS_PREFIX,
  datasuffix::String = DATA_SUFFIX)

  if topMethod == :rank #this assumes we want the top x ranked asset holders
    subdf = getTopInstitutionsRank(subdf, topInstitutionsRank)

  elseif topMethod == :quantile #this assumes we want the top  x percentile asset holders
    subdf = getTopInstitutionsQuantile(subdf, topInstitutionsQuantile)

  elseif topMethod == :cumulative #this assumes we want the top holders which account for x% of assets
    subdf = getTopInstitutionsCumulative(subdf, topInstitutionsCumulative)
  elseif topMethod ≠ :full
      @warn "Top institutions method not detected, assuming full file"
  end

  label::String =
    "$(topInstitutionsPrefix)_$(topInstitutionsYear)_$(topMethod)_$datasuffix"
  writeTopInstitutions(subdf, label)

  return nothing
end

#takes in an unsorted data object and gets the top institutions by assets
function getTopInstitutions(data::NCCSData, topMethod::Symbol;
    topInstitutionsYear::Int = TOP_INSTITUTIONS_YEAR)
  sort!(data.df, :netassets, rev=true)
  subdf::SubDataFrame =
    view(data.df, (data.df[:fisyr] .== topInstitutionsYear) .& (
    (!ismissing).(data.df[:netassets])), :)

  getTopInstitutions(subdf, topMethod)

  return nothing
end

#convenience method for running multiple methods at once
function getTopInstitutions(data::NCCSData;
  topMethods::Vector{Symbol} = TOP_INSTITUTION_METHODS,
  topInstitutionsYear::Int = TOP_INSTITUTIONS_YEAR,
  sorted::Bool = false)

  (!sorted) && sort!(data.df, [:fisyr, :adjnetassets], (:fisyr, order(:adjnetassets, rev=true)))
  subdf::SubDataFrame =
    view(data.df, (data.df[:fisyr] .== topInstitutionsYear) .& (
    (!ismissing).(data.df[:netassets])),:)

  #Execute the tasks
  for s::Symbol ∈ topMethods
    getTopInstitutions(subdf, s)
  end

  worstPerformers(data)

  return nothing
end

function worstPerformers(data::NCCSData; worstminassets = WORST_MIN_ASSETS,
  numyears = WORST_NUM_YEARS, lastyear = 2015,
  outputpath::String = TABLE_PATH,
  datasuffix::String = DATA_SUFFIX)::Nothing

  firstyear::Int = lastyear - numyears
  outpath = ("$outputpath\\worst_$(firstyear)_$(lastyear)" *
    "_$(worstminassets ÷ 1_000_000)_$datasuffix.csv")

  #partition off the part we want
  validfirmsdf = view(data.df, (data.df[:adjnetassets] .> worstminassets) .&
    (data.df[:fisyr] .≤ lastyear) .& (data.df[:fisyr] .> firstyear),
    [:ein, :name, :excess, :return, :t30ret, :sp500ret, :adjnetassets,
    :monthscovered, :fisyr, :datefiscal])

  #form the consolidated df

  worstdf::DataFrame = DataFrame(ein=validfirmsdf[:ein])
  N::Int = size(worstdf, 1)

  worstdf[:avgwealth] = Vector{MFloat64}(missing, N)
  worstdf[:N] = Vector{MInt}(missing, N)
  worstdf[:monthscovered] = Vector{MInt}(missing, N)
  worstdfindex::Dict = Dict(worstdf[i,:ein] => worstdf[i,:] for i::Int  ∈ 1:N)
  worstdf[:name] = Vector{MString}(missing, N)
  worstdf[:sp500ret] = Vector{MFloat64}(missing, N)
  worstdf[:t30ret] = Vector{MFloat64}(missing, N)
  worstdf[:cagr] = Vector{MFloat64}(missing, N)

  #consolidate the eprformacne data
  for subdf ∈ groupby(validfirmsdf, :ein)
    worstrow = worstdfindex[subdf[1,:ein]]
    worstrow[:name] = subdf[1,:name]

    worstrow[:N] = size(subdf,1)
    worstrow[:monthscovered] = sum(subdf[:monthscovered])

    totalreturn::MFloat64 = prod(1.0 .+ subdf[:return])
    worstrow[:cagr] = totalreturn^(12. / worstrow[:monthscovered]) .- 1

    totalrfr::MFloat64 = prod(1.0 .+ subdf[:t30ret])
    worstrow[:t30ret] = totalrfr^(12. / worstrow[:monthscovered]) .- 1.

    totalsp500::MFloat64 = prod(1.0 .+ subdf[:sp500ret])
    worstrow[:sp500ret] = totalsp500^(12. / worstrow[:monthscovered]) .- 1.

    worstrow[:avgwealth] = mean(subdf[:adjnetassets]) ./ 1_000_000
  end

  #compute excess returns
  worstdf[:excess] = worstdf[:cagr] .- worstdf[:t30ret]
  worstdf = worstdf[completecases(worstdf[[:cagr,:t30ret,:avgwealth,:sp500ret,:excess]]),:]
  worstdf = worstdf[worstdf[:N] .== numyears,:]

  sort!(worstdf, :excess, rev=false)

  worstdf |> CSV.write(outpath)

  return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using LinearAlgebra
using Test

using SetProg
using Polyhedra

using JuMP

@testset "Algebraic" begin
    □ = polyhedron(HalfSpace([1, 0], 1.0) ∩ HalfSpace([-1, 0], 1) ∩ HalfSpace([0, 1], 1) ∩ HalfSpace([0, -1], 1))
    @testset "Ellipsoid" begin
        @testset "John" begin
            mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
            model = JuMP.direct_model(mock);
            # Q = [1 0
            #      0 1]
            # t = √det(Q) = 1                                                                  Q11  Q12  Q22  t
            MOI.Utilities.set_mock_optimize!(mock, mock -> MOI.Utilities.mock_optimize!(mock, [1.0, 0.0, 1.0, 1.0]));

            @variable(model, ◯, Ellipsoid(symmetric=true, dimension=2))
            cref = @constraint(model, ◯ ⊆ □)
            @objective(model, Max, nth_root(volume(◯)))

            SetProg.optimize!(model)
            @test JuMP.termination_status(model) == MOI.OPTIMAL
            @test JuMP.objective_sense(model) == MOI.MAX_SENSE
            @test JuMP.objective_value(model) == 1.0
            @test JuMP.value(◯) isa SetProg.Sets.PolarEllipsoid
            @test JuMP.value(◯).Q == Symmetric([1.0 0.0; 0.0 1.0])
        end
        @testset "Löwner" begin
            mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
            model = JuMP.direct_model(mock);
            # Q = [√2  0
            #       0 √2]
            # t = √det(Q) = 2                                                                  Q11  Q12  Q22  t
            MOI.Utilities.set_mock_optimize!(mock, mock -> MOI.Utilities.mock_optimize!(mock, [0.5, 0.0, 0.5, 0.5]));

            @variable(model, ◯, Ellipsoid(symmetric=true, dimension=2))
            cref = @constraint(model, □ ⊆ ◯)
            @objective(model, Min, nth_root(volume(◯)))

            SetProg.optimize!(model)
            @test JuMP.termination_status(model) == MOI.OPTIMAL
            @test JuMP.objective_sense(model) == MOI.MAX_SENSE
            @test JuMP.objective_value(model) == 0.5
            @test JuMP.value(◯) isa SetProg.Sets.Ellipsoid
            @test JuMP.value(◯).Q == Symmetric([0.5 0.0; 0.0 0.5])
        end
    end
    @testset "Quartic" begin
        @testset "Inner" begin
            mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
            model = JuMP.direct_model(mock);
            # The PSD matrix for the variable  is 3 x 3 so 3 * (3+1) / 2 = 6
            # The PSD matrix for the convexity is 6 x 6 so 6 * (6+1) / 2 = 21
            # entries
            # 1 variable for t
            # hence 28 variables
            MOI.Utilities.set_mock_optimize!(mock, mock -> MOI.Utilities.mock_optimize!(mock, ones(28)))

            @variable(model, ◯, PolySet(symmetric=true, degree=4, dimension=2,
                                        convex=true))
            cref = @constraint(model, ◯ ⊆ □)
            @objective(model, Max, nth_root(volume(◯)))

            SetProg.optimize!(model)
            @test JuMP.termination_status(model) == MOI.OPTIMAL
            @test JuMP.objective_sense(model) == MOI.MAX_SENSE
            @test JuMP.objective_value(model) == 1.0
            @test JuMP.value(◯) isa SetProg.Sets.PolarPolySet{Float64}
            @test JuMP.value(◯).degree == 4
            x, y = SetProg.data(model).polyvars
            @test polynomial(JuMP.value(◯).p) == x^4 + 2x^3*y + 3x^2*y^2 + 2x*y^3 + y^4
            @test JuMP.value(◯).convexity_proof.n == 6
            @test JuMP.value(◯).convexity_proof.Q == ones(21)
            @test JuMP.objective_value(model) == 1.0
        end
        @testset "Outer" begin
            mock = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
            model = JuMP.direct_model(mock);
            MOI.Utilities.set_mock_optimize!(mock, mock -> MOI.Utilities.mock_optimize!(mock, ones(28)))

            @variable(model, ◯, PolySet(symmetric=true, degree=4, dimension=2,
                                        convex=true))
            cref = @constraint(model, □ ⊆ ◯)
            @objective(model, Min, nth_root(volume(◯)))

            SetProg.optimize!(model)
            @test JuMP.termination_status(model) == MOI.OPTIMAL
            @test JuMP.objective_sense(model) == MOI.MAX_SENSE
            @test JuMP.objective_value(model) == 1.0
            @test JuMP.value(◯) isa SetProg.Sets.PolySet{Float64}
            @test JuMP.value(◯).degree == 4
            x, y = SetProg.data(model).polyvars
            @test polynomial(JuMP.value(◯).p) == x^4 + 2x^3*y + 3x^2*y^2 + 2x*y^3 + y^4
            @test JuMP.value(◯).convexity_proof.n == 6
            @test JuMP.value(◯).convexity_proof.Q == ones(21)
        end
    end
end


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.Model{Float64}())
config = MOIT.TestConfig()

@testset "QuadtoSOC" begin
    bridged_mock = MOIB.Constraint.QuadtoSOC{Float64}(mock)
    @testset "Error for non-convex quadratic constraints" begin
        x = MOI.add_variable(bridged_mock)
        @test_throws ErrorException begin
            MOI.add_constraint(bridged_mock,
                               MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Float64}[],
                                                           [MOI.ScalarQuadraticTerm(1.0, x, x)],
                                                           0.0),
                               MOI.GreaterThan(0.0))
        end
        @test_throws ErrorException begin
            MOI.add_constraint(bridged_mock,
                               MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Float64}[],
                                                           [MOI.ScalarQuadraticTerm(-1.0, x, x)],
                                                           0.0),
                               MOI.LessThan(0.0))
        end
    end
    @testset "Quadratic constraints with 2 variables" begin
        MOIU.set_mock_optimize!(mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [0.5, 0.5])
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [0.5, (√13 - 1)/4])
            )
        )
        MOIT.solve_qcp_edge_cases(bridged_mock, config)
        ci = first(MOI.get(mock,
                           MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64},
                                                       MOI.RotatedSecondOrderCone}()))
        x, y = MOI.get(mock, MOI.ListOfVariableIndices())
        # The matrix is
        # 2 1
        # 1 2
        # for which the cholesky factorization is U' * U with U =
        # √2 √2/2
        #  . √3/√2
        expected = MOI.VectorAffineFunction{Float64}(MOI.VectorAffineTerm.([3, 3, 4],
                                                                           MOI.ScalarAffineTerm.([√2, √2/2, √3/√2],
                                                                                [x, y, y])),
                                                     [1.0, 1.0, 0.0, 0.0])
        @test MOI.get(mock, MOI.ConstraintFunction(), ci) ≈ expected
    end
    @testset "QCP tests" begin
        MOIU.set_mock_optimize!(mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1/2, 7/4],
                (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives)  => [zeros(2)],
                (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone)  => [[0.25, 1.0, -1/√2]]))
        MOIT.qcp1test(bridged_mock, config)
        MOIU.set_mock_optimize!(mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [√2],
                (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone)  => [[1/√2, 1/(2 * √2), -1/√2]]))
        MOIT.qcp2test(bridged_mock, config)
        MOIT.qcp3test(bridged_mock, config)
        @testset "Bridge deletion" begin
            ci = first(MOI.get(bridged_mock,
                               MOI.ListOfConstraintIndices{MOI.ScalarQuadraticFunction{Float64},
                                                           MOI.LessThan{Float64}}()))
            test_delete_bridge(bridged_mock, ci, 1, ((MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),))
        end
        MOIU.set_mock_optimize!(mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 1.0],
                (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone)  => [[1.0, 1/3, -1/√2, -1/√6]]))
        MOIT.qcp4test(bridged_mock, config)
        MOIT.qcp5test(bridged_mock, config)
    end
end


text:
all_lines = readlines("day8.txt")
all_lines_vec = hcat(split.(all_lines," | ")...)

outputs = all_lines_vec[2,:]

output_vales = hcat(split.(outputs," ")...)

lenghts = length.(output_vales)

ones = lenghts .== 2
fours = lenghts .== 4
sevens = lenghts .== 3
eights = lenghts .== 7

vals_1478 = ones + fours + sevens + eights
sum(vals_1478)


text:
"""

    SFFMDistribution(
        pm::Array{<:Real},
        distribution::Array{<:Real,3},
        x::Array{<:Real},
        type::String,
    )

- `pm::Array{Float64}`: a vector containing the point masses, the first
    `sum(model.C.<=0)` entries are the left hand point masses and the last
    `sum(model.C.>=0)` are the right-hand point masses.
- `distribution::Array{Float64,3}`: "probability" or "density"`
- `x::Array{Float64,2}`:
    - if `type="probability"` is a `1×NIntervals×NPhases` array
        containing the cell centers.
    - if `type="density"` is a `NBases×NIntervals×NPhases` array
        containing the cell nodes at which the denisty is evaluated.
- `type::String`: either `"probability"` or `"density"`. `"cumulative"` is
    not possible.
"""
struct SFFMDistribution
    pm::Array{<:Real}
    distribution::Array{<:Real,3}
    x::Array{<:Real}
    type::String
end

"""
Convert from a vector of coefficients for the DG system to a distribution.

    Coeffs2Dist(
        model::SFFM.Model,
        mesh::SFFM.Mesh,
        Coeffs;
        type::String = "probability",
        probTransform::Bool = true,
    )

# Arguments
- `model`: a Model object
- `mesh`: a Mesh object as output from MakeMesh
- `Coeffs::Array`: a vector of coefficients from the DG method
- `type::String`: an (optional) declaration of what type of distribution you
    want to convert to. Options are `"probability"` to return the probabilities
    ``P(X(t)∈ D_k, φ(t) = i)`` where ``D_k``is the kth cell, `"cumulative"` to
    return the CDF evaluated at cell edges, or `"density"` to return an
    approximation to the density ar at the SFFM.CellNodes(mesh).
- `probTransform::Bool` a boolean value specifying whether to transform to a
    probabilistic interpretation or not. Valid only for lagrange basis.

# Output
- a tuple with keys
(pm=pm, distribution=yvals, x=xvals, type=type)
    - `pm::Array{Float64}`: a vector containing the point masses, the first
        `sum(model.C.<=0)` entries are the left hand point masses and the last
        `sum(model.C.>=0)` are the right-hand point masses.
    - `distribution::Array{Float64,3}`:
        - if `type="cumulative"` returns a `2×NIntervals×NPhases` array
            containing the CDF evaluated at the cell edges as contained in
            `x` below. i.e. `distribution[1,:,i]` returns the cdf at the
            left-hand edges of the cells in phase `i` and `distribution[2,:,i]`
            at the right hand edges.
        - if `type="probability"` returns a `1×NIntervals×NPhases` array
            containing the probabilities ``P(X(t)∈ D_k, φ(t) = i)`` where ``D_k``
            is the kth cell.
        - if `type="density"` returns a `NBases×NIntervals×NPhases` array
            containing the density function evaluated at the cell nodes as
            contained in `x` below.
    - `x::Array{Float64,2}`:
        - if `type="cumulative"` returns a `2×NIntervals×NPhases` array
            containing the cell edges as contained. i.e. `x[1,:]`
            returns the left-hand edges of the cells and `x[2,:]` at the
            right-hand edges.
        - if `type="probability"` returns a `1×NIntervals×NPhases` array
            containing the cell centers.
        - if `type="density"` returns a `NBases×NIntervals×NPhases` array
            containing the cell nodes.
    - `type`: as input in arguments.
"""
function Coeffs2Dist(
    model::SFFM.Model,
    mesh::SFFM.DGMesh,
    Coeffs::AbstractArray;
    type::String = "probability",
    probTransform::Bool = true,
)
    V = SFFM.vandermonde(NBases(mesh))
    N₋ = sum(model.C .<= 0)
    N₊ = sum(model.C .>= 0)
    if !probTransform
        temp = reshape(Coeffs[N₋+1:end-N₊], NBases(mesh), NIntervals(mesh), NPhases(model))
        temp = V.w.*temp.*(Δ(mesh)./2.0)'
        Coeffs = [Coeffs[1:N₋]; temp[:]; Coeffs[end-N₊+1:end]]
    end
    if type == "density"
        xvals = SFFM.CellNodes(mesh)
        if Basis(mesh) == "legendre"
            yvals = reshape(Coeffs[N₋+1:end-N₊], NBases(mesh), NIntervals(mesh), NPhases(model))
            for i in 1:NPhases(model)
                yvals[:,:,i] = V.V * yvals[:,:,i]
            end
            pm = [Coeffs[1:N₋]; Coeffs[end-N₊+1:end]]
        elseif Basis(mesh) == "lagrange"
            yvals =
                Coeffs[N₋+1:end-N₊] .* repeat(1.0 ./ V.w, NIntervals(mesh) * NPhases(model)) .*
                (repeat(2.0 ./ Δ(mesh), 1, NBases(mesh) * NPhases(model))'[:])
            yvals = reshape(yvals, NBases(mesh), NIntervals(mesh), NPhases(model))
            pm = [Coeffs[1:N₋]; Coeffs[end-N₊+1:end]]
        end
        if NBases(mesh) == 1
            yvals = [1;1].*yvals
            xvals = [SFFM.CellNodes(mesh)-Δ(mesh)'/2;SFFM.CellNodes(mesh)+Δ(mesh)'/2]
        end
    elseif type == "probability"
        if NBases(mesh) > 1
            xvals = SFFM.CellNodes(mesh)[1, :] + (Δ(mesh) ./ 2)
        else
            xvals = SFFM.CellNodes(mesh)
        end
        if Basis(mesh) == "legendre"
            yvals = (reshape(Coeffs[N₋+1:NBases(mesh):end-N₊], 1, NIntervals(mesh), NPhases(model)).*Δ(mesh)')./sqrt(2)
            pm = [Coeffs[1:N₋]; Coeffs[end-N₊+1:end]]
        elseif Basis(mesh) == "lagrange"
            yvals = sum(
                reshape(Coeffs[N₋+1:end-N₊], NBases(mesh), NIntervals(mesh), NPhases(model)),
                dims = 1,
            )
            pm = [Coeffs[1:N₋]; Coeffs[end-N₊+1:end]]
        end
    elseif type == "cumulative"
        if NBases(mesh) > 1
            xvals = SFFM.CellNodes(mesh)[[1;end], :]
        else
            xvals = [SFFM.CellNodes(mesh)-Δ(mesh)'/2;SFFM.CellNodes(mesh)+Δ(mesh)'/2]
        end
        if Basis(mesh) == "legendre"
            tempDist = (reshape(Coeffs[N₋+1:NBases(mesh):end-N₊], 1, NIntervals(mesh), NPhases(model)).*Δ(mesh)')./sqrt(2)
            pm = [Coeffs[1:N₋]; Coeffs[end-N₊+1:end]]
        elseif Basis(mesh) == "lagrange"
            tempDist = sum(
                reshape(Coeffs[N₋+1:end-N₊], NBases(mesh), NIntervals(mesh), NPhases(model)),
                dims = 1,
            )
            pm = [Coeffs[1:N₋]; Coeffs[end-N₊+1:end]]
        end
        tempDist = cumsum(tempDist,dims=2)
        temppm = zeros(Float64,1,2,NPhases(model))
        temppm[:,1,model.C.<=0] = pm[1:N₋]
        temppm[:,2,model.C.>=0] = pm[N₊+1:end]
        yvals = zeros(Float64,2,NIntervals(mesh),NPhases(model))
        yvals[1,2:end,:] = tempDist[1,1:end-1,:]
        yvals[2,:,:] = tempDist
        yvals = yvals .+ reshape(temppm[1,1,:],1,1,NPhases(model))
        pm[N₋+1:end] = pm[N₋+1:end] + yvals[end,end,model.C.>=0]
    end

    out = SFFMDistribution(pm, yvals, xvals, type)
    println("UPDATE: distribution object created with keys ", fieldnames(SFFMDistribution))
    return out
end
function Coeffs2Dist(
    model::SFFM.Model,
    mesh::Union{FRAPMesh, FVMesh},
    Coeffs::AbstractArray;
    type::String = "probability",
)

    if type != "probability"
        args = [
            model;
            mesh;
            Coeffs;
            type;
        ]
        error("Input Error: no functionality other than 'probability' implemented")
    end

    N₋ = sum(model.C .<= 0)
    N₊ = sum(model.C .>= 0)

    xvals = SFFM.CellNodes(mesh)

    yvals = sum(
        reshape(Coeffs[N₋+1:end-N₊], NBases(mesh), NIntervals(mesh), NPhases(model)),
        dims = 1,
    )
    pm = [Coeffs[1:N₋]; Coeffs[end-N₊+1:end]]

    out = SFFMDistribution(pm, yvals, xvals, type)
    println("UPDATE: distribution object created with keys ", fieldnames(SFFMDistribution))
    return out
end

"""
Converts a distribution as output from `Coeffs2Dist()` to a vector of DG
coefficients.

    Dist2Coeffs(
        model::SFFM.Model,
        mesh::SFFM.Mesh,
        Distn::SFFMDistribution;
        probTransform::Bool = true,
    )

# Arguments
- `model`: a Model object
- `mesh`: a Mesh object as output from MakeMesh
- `Distn::SFFMDistribution
    - if `type="probability"` is a `1×NIntervals×NPhases` array containing
        the probabilities ``P(X(t)∈ D_k, φ(t) = i)`` where ``D_k``
        is the kth cell.
    - if `type="density"` is a `NBases×NIntervals×NPhases` array containing
        either the density function evaluated at the cell nodes which are in
        `x` below, or, the inner product of the density function against the
        lagrange polynomials.
- `probTransform::Bool` a boolean value specifying whether to transform to a
    probabilistic interpretation or not. Valid only for lagrange basis.

# Output
- `coeffs` a row vector of coefficient values of length
    `TotalNBases*NPhases + N₋ + N₊` ordered according to LH point masses, RH
    point masses, interior basis functions according to basis function, cell,
    phase. Used to premultiply operators such as B from `MakeB()`
"""
function Dist2Coeffs(
    model::SFFM.Model,
    mesh::SFFM.DGMesh,
    Distn::SFFMDistribution;
    probTransform::Bool = true,
)
    V = SFFM.vandermonde(NBases(mesh))
    theDistribution =
        zeros(Float64, NBases(mesh), NIntervals(mesh), NPhases(model))
    if Basis(mesh) == "legendre"
        if Distn.type == "probability"
            # for the legendre basis the first basis function is ϕ(x)=Δ√2 and
            # all other basis functions are orthogonal to this. Hence, we map
            # the cell probabilities to the first basis function only.
            theDistribution[1, :, :] = Distn.distribution./Δ(mesh)'.*sqrt(2)
        elseif Distn.type == "density"
            # if given density coefficients in lagrange form
            theDistribution = Distn.distribution
            for i = 1:NPhases(model)
                theDistribution[:, :, i] = V.inv * theDistribution[:, :, i]
            end
        end
        # also put the point masses on the ends
        coeffs = [
            Distn.pm[1:sum(model.C .<= 0)]
            theDistribution[:]
            Distn.pm[sum(model.C .<= 0)+1:end]
        ]
    elseif Basis(mesh) == "lagrange"
        theDistribution .= Distn.distribution
        if !probTransform
            theDistribution = (1.0./V.w) .* theDistribution .* (2.0./Δ(mesh)')
        end
        if Distn.type == "probability"
            # convert to probability coefficients by multiplying by the
            # weights in V.w/2
            theDistribution = (V.w .* theDistribution / 2)[:]
        elseif Distn.type == "density"
            # convert to probability coefficients by multiplying by the
            # weights in V.w/2 and cell widths Δ
            theDistribution = ((V.w .* theDistribution).*(Δ(mesh) / 2)')[:]
        end
        # also put the point masses on the ends
        coeffs = [
            Distn.pm[1:sum(model.C .<= 0)]
            theDistribution
            Distn.pm[sum(model.C .<= 0)+1:end]
        ]
    end
    coeffs = Matrix(coeffs[:]')
    return coeffs
end
function Dist2Coeffs(
    model::SFFM.Model,
    mesh::Union{SFFM.FRAPMesh,SFFM.FVMesh},
    Distn::SFFMDistribution;
    probTransform::Bool = true,
)

    # also put the point masses on the ends
    coeffs = [
        Distn.pm[1:sum(model.C .<= 0)]
        Distn.distribution[:]
        Distn.pm[sum(model.C .<= 0)+1:end]
    ]

    coeffs = Matrix(coeffs[:]')
    return coeffs
end

"""
Computes the error between distributions.

    starSeminorm(
        d1::SFFMDistribution,
        d2::SFFMDistribution,
        )

# Arguments
- `d1`: a distribution object as output from `Coeffs2Dist` with
    `type="probability"``
- `d2`: a distribution object as output from `Coeffs2Dist` with
    `type="probability"``
"""
function starSeminorm(
    d1::SFFMDistribution,
    d2::SFFMDistribution,
    )
    if ((d1.type!="probability") || (d2.type!="probability"))
        throw(ArgumentError("distributions need to be of type probability"))
    end
    return sum(abs.(d1.pm-d2.pm)) + sum(abs.(d1.distribution-d2.distribution))
end

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
using DashBootstrapComponents, DashHtmlComponents

buttons = html_div([
    dbc_button("Large button", size = "lg", className = "me-1"),
    dbc_button("Regular button", className = "me-1"),
    dbc_button("Small button", size = "sm"),
]);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
## PART 1

println("Part 1");

# Define functions for this part
function f1(v̄)
    return v̄[:,1] ./ v̄[:,2];
end

function f2(v̄)
    return exp( v̄[:,3] - v̄[:,4] );
end

function f3(v̄)
    logs = Float64[];
    for vv in v̄[:,5]
        push!(logs,log(vv));
    end
    return ( v̄[:,1]./v̄[:,2] + v̄[:,3]./v̄[:,4] ) .* logs;
end

# Calculate f1, f2, and f3 for these means
println("RESUMING QUESTION 2");
println("\n\tCalculating functions of means f1v1, f2v1, ...");
f1N1 = f1(v̄1);
f2N1 = f2(v̄1);
f3N1 = f3(v̄1);

f1N2 = f1(v̄2);
f2N2 = f2(v̄2);
f3N2 = f3(v̄2);

# Find standard deviations
println("\n\tCalculating standard deviations of the functions");
σ̂truef1N1 = std(f1N1);
σ̂truef2N1 = std(f2N1);
σ̂truef3N1 = std(f3N1);

σ̂truef1N2 = std(f1N2);
σ̂truef2N2 = std(f2N2);
σ̂truef3N2 = std(f3N2);


text:
"""
Kernel
"""

struct LegendreBasis <: Basis
    beta::Float64
    o::PyObject
    u::PiecewiseLegendrePoly
    uhat::PiecewiseLegendreFT
    v::Nothing
    statistics::Statistics
    size::Int64
end

function LegendreBasis(o::PyObject)
    return LegendreBasis(o.beta, o, o.u, o.uhat, o.v,
                         o.statistics == "F" ? fermion : boson,
                         o.size)
end

function LegendreBasis(statistics::Statistics, beta::Real, size::Int64;
                       cl::Vector{Float64}=ones(Float64, size))
    return LegendreBasis(pyaugment.LegendreBasis(statistics, beta, size, cl))
end

function default_tau_sampling_points(basis::LegendreBasis)
    return basis.o.default_tau_sampling_points()
end

function default_matsubara_sampling_points(basis::LegendreBasis; mitigate=true)
    return basis.o.default_matsubara_sampling_points()
end


text:
abstract type Model end

mutable struct SnakeModel <: Model
    status::Symbol
    size_properties::Tuple{Int, Int, Int, Int} # node_h, node_w, parent_h, parent_w
    direction::Symbol
    node_anchors::Vector{Tuple{Int, Int}}
    global_event_queue::Channel
end

function SnakeModel(parent_h::Int, parent_w::Int, event_queue::Channel; init_pos=(4, 50))
    # ┌┐
    # └┘
    node_h = node_w = 2
    anchors = [init_pos]
    for i=1:3
        push!(anchors, (init_pos[1], init_pos[2]-node_w*i))
    end

    return SnakeModel(:alive, (node_h, node_w, parent_h, parent_w), :right, anchors, event_queue)
end,

struct UpdateEvent <: T.Event end

struct LossEvent <: T.Event end

function bump(snake_model::SnakeModel)
    put!(snake_model.global_event_queue, LossEvent())
    snake_model.status = :dead
    return false
end

function move(snake_model::SnakeModel, direction::Symbol)
    (snake_model.status === :dead) && return false
    (snake_model.direction===:up && direction===:down) && return true
    (snake_model.direction===:down && direction===:up) && return true
    (snake_model.direction===:right && direction===:left) && return true
    (snake_model.direction===:left && direction===:right) && return true

    h, w, ph, pw = snake_model.size_properties
    if direction === :up
        new_anchor = snake_model.node_anchors[1].-(h, 0)
        (new_anchor[1] < 1) && return bump(snake_model)
    elseif direction === :down
        new_anchor = snake_model.node_anchors[1].+(h, 0)
        (new_anchor[1]+h > ph) && return bump(snake_model)
    elseif direction === :right
        new_anchor = snake_model.node_anchors[1].+(0, w)
        (new_anchor[2]+w > pw) && return bump(snake_model)
    elseif direction === :left
        new_anchor = snake_model.node_anchors[1].-(0, w)
        (new_anchor[2] < 1) && return bump(snake_model)
    end

    pop!(snake_model.node_anchors)
    pushfirst!(snake_model.node_anchors, new_anchor)
    snake_model.direction = direction

    put!(snake_model.global_event_queue, UpdateEvent())

    return true
end

function auto_move(snake_model::SnakeModel)
    Base.Threads.@spawn begin
        while move(snake_model, snake_model.direction)
            sleep(0.5)
        end
    end
end


text:
using StanSample, Test

@testset "StanRun cmdstan run and results" begin
    # setup environment
    MODELDIR = mktempdir()
    SAMPLEDIR = mktempdir()
    SRC = joinpath(MODELDIR, "test.stan")
    cp(joinpath(@__DIR__, "test.stan"), SRC)

    # run a model
    model = StanModel(SRC)
    @test repr(model) ==
        "Stan model at $(SRC)\n    (CmdStan home: $(StanRun.get_cmdstan_home()))"
    exec_path = StanRun.ensure_executable(model)
    @test isfile(exec_path)
    n_chains = 5
    OUTPUT_BASE = joinpath(SAMPLEDIR, "test")
    time_before = time()
    @test stan_compile(model) ≡ nothing
    chains = stan_sample(model, (N = 100, x = randn(100)), 5; output_base = OUTPUT_BASE)
    for (sample, logfile) in chains
        @test ctime(sample) ≥ time_before
        @test ctime(sample) ≥ time_before
    end
    @test first.(chains) == sort(StanRun.find_samples(OUTPUT_BASE)) ==
        [joinpath(SAMPLEDIR, "test_chain_$(i).csv") for i in 1:n_chains]
end