text:
#=
  - Julia Version : 1.4.0
  - Author        : Alexander Bills and Abhinav Misalkar
  - Date          : 14th May 2020
  - Title         : Cylindrical Cell Designer
=#

mutable struct struct_electrode
       chem::String                                   # Chemistry type

       bind_wt_fr::Float64                            # binder weight fraction
       cond_wt_fr::Float64                            # conductive weight fraction

       AM_rho::Float64                                # AM density                                      g/cm3
       bind_rho::Float64                              # binder density                                  g/cm3
       cond_rho::Float64                              # conductive density                              g/cm3

       fc_sp_cap::Float64                             # first charge specific capacity                  mAh/g   --> Ah/g
       rev_sp_cap::Float64                            # reversible specific capacity                    mAh/g   --> Ah/g

       th::Float64                                    # total thickness cathode                       microns   --> cm
       por::Float64                                   # porosity

       extra_length::Float64                          # extra length                                     mm     --> cm

       CC_th::Float64                                 # Current Collector Thickness                   microns   --> cm
       CC_rho::Float64                                # Current Collector                               g/cm3

       tab_th::Float64                                # Tab thickness
       tab_width::Float64                             # Tab width
       tab_rho::Float64                               # Tab Density
end

mutable struct struct_cell_general
    form_factor::String              # ["Cyl","Pris","Pou"]
    size::String                     # ["18650", "21700", etc]
    dimensional_delta::Float64
    can_th::Float64
    can_rho::Float64

    sep_th::Float64
    sep_rho::Float64
    sep_por::Float64

    nom_volt::Float64

    electrolyte_type::String
    electrolyte_rho::Float64

    N_P_ratio::Float64               # N to P ratio

    anode::struct_electrode
    cathode::struct_electrode
end

mutable struct struct_cell_design_op
    rev_cap::Float64
    energy::Float64
    mass_cathode_AM::Float64
    mass_anode_AM::Float64
    pos_CC_area::Float64
    neg_CC_area::Float64
    sep_area::Float64
    electrolyte_volm::Float64
    mass_tab::Float64
    mass_of_cannister::Float64
    total_coated_area_cathode::Float64
    pos_AM_wt_fr::Float64
    neg_AM_wt_fr::Float64
end

function electrode_geometry(electrode, volm_jellyroll, jellyroll_length, height_jellyroll)
    # Calculate Eletrode Geometry Characteristics
    # Calculate Weight Fraction

    AM_wt_fr = 1 - electrode.bind_wt_fr - electrode.cond_wt_fr

    dry_electrode_density  = (1 - electrode.por) * electrode.AM_rho                             # g/cm3

    single_side_coating_thickness = (electrode.th - electrode.CC_th) / 2.0                      # microns  --> cm2          - done

    mass_loading  = (dry_electrode_density  *  single_side_coating_thickness) #/ 10.0           # mg/cm2   --> g/cm2        - done

    capacity_loading = (mass_loading  *   electrode.rev_sp_cap) #/ 1000.0                       # mAh/cm2  --> Ah/cm2       - done

    length_electrode = jellyroll_length + electrode.extra_length                                # mm --> cm                 - done

    total_coated_length_electrode = 2 * length_electrode                                        # mm --> cm - done
    total_coated_area   = height_jellyroll * total_coated_length_electrode                      # cm2

    volm_electrode =  (total_coated_length_electrode ) *                                        # cm3
                      (single_side_coating_thickness) *
                      height_jellyroll

    mass_AM = volm_electrode * dry_electrode_density * AM_wt_fr                                 # grams

    mass_electrode    = volm_electrode * dry_electrode_density                                  # grams

    volm_electrode_AM = mass_AM / electrode.AM_rho                                              # cm3

    mass_electrode_bind = (mass_AM / AM_wt_fr) * electrode.bind_wt_fr                           # grams
    mass_electrode_cond = (mass_AM / AM_wt_fr) * electrode.cond_wt_fr                           # grams

    volm_bind  =  mass_electrode_bind * electrode.bind_rho                                      # cm3
    volm_cond  =  mass_electrode_cond * electrode.cond_rho                                      # cm3

    return electrode, mass_AM, capacity_loading, dry_electrode_density,
    single_side_coating_thickness, mass_electrode, total_coated_area,
    mass_electrode_bind, AM_wt_fr, total_coated_length_electrode,
    mass_loading, volm_electrode

end

function CC_geometry(electrode, jellyroll_length, height_jellyroll)
    CC_length = jellyroll_length + electrode.extra_length                                       # mm --> cm - done
    CC_area = (CC_length) * height_jellyroll                                                    # cm2 --> done
    CC_volm = CC_area * (electrode.CC_th)                                                       # cm3 --> done
    CC_mass = CC_volm * electrode.CC_rho                                                        # grams --> done
    return CC_mass, CC_area
end

function np_designer(electrode,capacity_loading)
    mass_loading = (capacity_loading) / (electrode.rev_sp_cap)                                  # mg/cm2 --> g/cm2
    single_side_coating_thickness = (electrode.th - electrode.CC_th) / 2.0                      # microns --> cm
    dry_electrode_density = (mass_loading) / (single_side_coating_thickness)                    # g/cm3
    electrode.por = 1 - (dry_electrode_density / electrode.AM_rho)
    return electrode
end

function cylindrical_cell_designer(cell;verbosity=0)
    cathode = cell.cathode
    anode = cell.anode

                    ###################################### Cell Parameters #####################################

    external_height_of_cell = (parse(Float64, SubString(cell.size, 3:4)))/10.0                              # mm --> cm - done
    internal_height_of_cell = external_height_of_cell * (1 - (cell.dimensional_delta))                      #  cm - done

    external_radius_of_cell = (parse(Float64, SubString(cell.size, 1:2))/2)/10.0                           # mm --> cm - done
    internal_radius_of_cell = (external_radius_of_cell - cell.can_th)                                       # mm --> cm - done

    external_radius_of_jellyroll = internal_radius_of_cell * (1 - (cell.dimensional_delta))                 # mm --> cm - done

    internal_radius_of_jellyroll = external_radius_of_cell/3.0                                              # mm --> cm - done
    height_jellyroll = internal_height_of_cell                                                              # cm

    height_of_tab = height_jellyroll                                                                        # cm

    volm_jellyroll =
        π *(external_radius_of_jellyroll^2 - internal_radius_of_jellyroll^2) *
        (height_jellyroll)                                                                                # mm3 --> cm3 - done

    # print("\nvolm_jellyroll = ", volm_jellyroll, " cm3")

    single_side_anode_coating_thickness = (cell.anode.th - cell.anode.CC_th) / 2.0                          # microns --> cm - done

    single_side_cathode_coating_thickness = (cell.cathode.th - cell.cathode.CC_th) / 2.0                    # microns --> cm - done

    # print("\nsingle_side_anode_coating_thickness = ", single_side_anode_coating_thickness, " cm")
    # print("\nsingle_side_cathode_coating_thickness = ", single_side_cathode_coating_thickness, " cm")
    # print("\nheight_jellyroll = ", height_jellyroll, " cm")

    jellyroll_length =                                                                                # mm --> cm  - done
        (volm_jellyroll) / ((height_jellyroll ) *
               (
                2 * (single_side_anode_coating_thickness ) +
                2 * (single_side_cathode_coating_thickness ) +
                (cell.cathode.CC_th ) +
                (cell.anode.CC_th ) +
                2 * (cell.sep_th )))

    # print("\njellyroll_length = ",jellyroll_length, " cm")

                        ###################################### Volume and Mass #####################################

   pos_electrode, mass_pos_AM, pos_capacity_loading, pos_dry_electrode_density,
   pos_single_side_coating_thickness, mass_pos_electrode, total_coated_area_cathode,
   mass_pos_bind, pos_AM_wt_fr, total_coated_length_cathode, pos_mass_loading,
   volm_pos_electrode = electrode_geometry(cathode, volm_jellyroll, jellyroll_length, height_jellyroll)

   anode_capacity_loading = cell.N_P_ratio * pos_capacity_loading                                             #  mAh/cm2

   anode = np_designer(anode,anode_capacity_loading)

   neg_electrode, mass_neg_AM, neg_capacity_loading, neg_dry_electrode_density,
   neg_single_side_coating_thickness, mass_neg_electrode, total_coated_area_anode,
   mass_neg_bind, neg_AM_wt_fr, total_coated_length_anode, neg_mass_loading,
   volm_neg_electrode  = electrode_geometry(anode,volm_jellyroll, jellyroll_length,height_jellyroll)

   if cell.cathode.por > 1 || cell.cathode.por < 0
        # print(cell.cathode.por)
       @error "Invalid cathode porosity value"
   end

    if cell.anode.por > 1 || cell.anode.por < 0
        print("\n1. Neg por",cell.anode.por)
       @error "Invalid anode porosity value"
   end

    sep_area = 2*(jellyroll_length + max(cathode.extra_length, anode.extra_length)) *                         # cm2
                height_jellyroll

    separator_volm = sep_area * (cell.sep_th)                                                                 # cm3

    separator_mass = separator_volm * cell.sep_rho * (1 - cell.sep_por)                                       # grams

    pos_CC_mass, pos_CC_area = CC_geometry(cathode, jellyroll_length, height_jellyroll)                       # grams, cm2
    neg_CC_mass, neg_CC_area = CC_geometry(anode,   jellyroll_length, height_jellyroll)                       # grams, cm2

    electrolyte_volm =
           (total_coated_area_cathode * ((cell.cathode.th - cell.cathode.CC_th)))* cell.cathode.por +
           (total_coated_area_anode   * ((cell.anode.th - cell.anode.CC_th)/2.0)) * cell.anode.por +
           2* sep_area * (cell.sep_th) * (cell.sep_por)                                                       # cm3

    electrolyte_mass = electrolyte_volm * cell.electrolyte_rho                                        # grams

    volume_canister_material = pi * (
        (((external_radius_of_cell)^2.0) * (external_height_of_cell))   -
        (((internal_radius_of_cell)^2.0) * (internal_height_of_cell)) )                                       # cm3

    mass_of_cannister = volume_canister_material * cell.can_rho                                       # grams

    mass_tab = cathode.tab_rho * (cathode.tab_th) * (height_of_tab) * (cathode.tab_width)                     # grams

    electrode_area = (mass_pos_electrode / pos_dry_electrode_density) / (pos_single_side_coating_thickness)   # cm2

    total_cell_mass =                                                                                # grams
        mass_pos_AM +
        mass_neg_AM +
        separator_mass +
        electrolyte_mass +
        mass_pos_bind +
        mass_neg_bind +
        pos_CC_mass +
        neg_CC_mass +
        2 * mass_tab +
        mass_of_cannister

                        ###################################### Performance #####################################

    first_charge_capacity_cathode = (mass_pos_AM * cell.cathode.fc_sp_cap)                             # Ah
    first_charge_capacity_anode   = (mass_neg_AM * cell.anode.fc_sp_cap)                               # Ah

    reversible_capacity_cathode = (mass_pos_AM * cell.cathode.rev_sp_cap)                              # Ah
    reversible_capacity_anode   = (mass_neg_AM * cell.anode.rev_sp_cap)                                # Ah

    reversible_capacity = min(reversible_capacity_anode, reversible_capacity_cathode)                           # Ah

    energy_cell = reversible_capacity * cell.nom_volt                                                   # Wh

                        ###################################### Verbosity #####################################

    if verbosity == 1
        print("\n")
        print("\n\n\n******************************** \t\t\t\t\t\t\t\tDesign Report \t\t\t\t\t\t\t\t********************************")

        print("\n\n\n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t  Cathode\t\t\t\t\t  Anode ")
        print("\n\n\tTotal Coated Length (cm)                   = \t\t\t",total_coated_length_cathode,"\t\t\t",total_coated_length_anode,)
        print("\n\n\tTotal Coated Area (cm2)                    = \t\t\t",total_coated_area_cathode,"\t\t\t",total_coated_area_anode,)
        print("\n\n\tPorosity                                   = \t\t\t\t",cell.cathode.por,"\t\t\t\t\t\t\t",cell.anode.por,)
        print("\n\n\tCapacity_loading (Ah/cm2)                  = \t\t\t\t",pos_capacity_loading,"\t\t\t\t",neg_capacity_loading,)
        print("\n\n\tMass_loading (g/cm2)                       = \t\t\t",pos_mass_loading,"\t\t\t\t\t",neg_mass_loading,)
        print("\n\n\tDry_electrode_density (g/cm3)              = \t\t\t",pos_dry_electrode_density,"\t\t\t",neg_dry_electrode_density,)
        print("\n\n\tActive Material Used (g)                   = \t\t\t",mass_pos_AM,"\t\t\t",mass_neg_AM,)
        print("\n\n\tMass of Electrode (g)                      = \t\t\t",mass_pos_electrode,"\t\t\t",mass_neg_electrode,)
        print("\n\n\tFirst Charge Capacity (Ah)                 = \t\t\t",first_charge_capacity_cathode,"\t\t\t",first_charge_capacity_anode,)
        print("\n\n\tReversible Capacity (Ah)                   = \t\t\t",reversible_capacity_cathode,"\t\t\t",reversible_capacity_anode,)

        print("\n\n\n******************************** \t\t\t\t\t\t\t\tCell Sizing \t\t\t\t\t\t\t\t********************************")
        print("\n\n\n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\tWeight (g)\t\t\t\t  Volume(mL)")
        print("\n\n\tCathode Active Material                    = \t\t\t",mass_pos_AM)  #,"\t\t\t",volm_cathode_AM,)
        print("\n\n\tAnode Active Material                      = \t\t\t",mass_neg_AM)  #,"\t\t\t\t",volm_anode_AM,)
        print("\n\n\tSeparator                                  = \t\t\t",separator_mass)   #,"\t\t\t",separator_volm,)
        print("\n\n\tElectrolyte                                = \t\t\t",electrolyte_mass) #,"\t\t\t",electrolyte_volm,)
        print("\n\n\tCathode Electrode Binder                   = \t\t\t",mass_pos_bind)    #,"\t\t\t",volm_pos_bind,)
        print("\n\n\tAnode Electrode Binder                     = \t\t\t",mass_neg_bind)    #,"\t\t\t",volm_neg_bind,)
        print("\n\n\tCathode Current Collector                  = \t\t\t",pos_CC_mass)  #,"\t\t\t",pos_CC_volm,)
        print("\n\n\tAnode Current Collector                    = \t\t\t",neg_CC_mass)  #,"\t\t\t",neg_CC_volm,)
        print("\n\n\tcannister                                  = \t\t\t",mass_of_cannister)    #,"\t\t\t",volume_canister_material,)

        print("\n\n\n******************************** \t\t\t\t\t\t\t\tCell Specs \t\t\t\t\t\t\t\t********************************")
        print("\n\n\tCell Mass (g)                              = \t\t\t\t",total_cell_mass,)
        print("\n\n\tCell Energy (Wh)                           = \t\t\t\t",energy_cell,)
        print("\n\n\tSeparator Area (cm2)                       = \t\t\t\t",sep_area,)
        print("\n\n\tSeparator Area (m2)                        = \t\t\t\t",sep_area/10000,)
        print("\n\n\tMass of each Tab (g)                       = \t\t\t\t",mass_tab,)
        print("\n\n\tBATPAC_electrode_area (cm2)                = \t\t\t\t",electrode_area,)

        print("\n\n\n******************************** \t\t\t\t\t\t\t\tOther Details \t\t\t\t\t\t\t\t********************************")
        print("\n\n\n\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t\t  Cathode\t\t\t\t\t\t  Anode ")
        print("\n\n\n\tSingle_side_coating_thickness (microns)  = \t\t\t\t\t\t",single_side_cathode_coating_thickness,"\t\t\t\t\t\t\t",single_side_anode_coating_thickness,)
        print("\n\n\tElectrode Thickness (microns)              = \t\t\t\t\t",cell.cathode.th,"\t\t\t\t\t\t\t",cell.anode.th,)
        print("\n\n\tElectrode Volume (cm3)                     = \t\t\t\t\t",volm_pos_electrode,"\t\t\t",volm_neg_electrode,)
        print("\n\n\tCurrent Collector Area (cm2)               = \t\t\t\t\t",pos_CC_area,"\t\t\t",neg_CC_area,)
        print("\n\n\tCurrent Collector Area (m2)                = \t\t\t\t\t",pos_CC_area/10000,"\t\t",neg_CC_area/10000,)

        print("\n")
    end

    cell_design_op = struct_cell_design_op(reversible_capacity_cathode,
                                             energy_cell,                       # kWh
                                             mass_pos_AM,                       # kg
                                             mass_neg_AM,                       # kg
                                             pos_CC_area,                       # m2
                                             neg_CC_area,                       # m2
                                             sep_area,                          # m2
                                             electrolyte_volm,                  # L
                                             mass_tab,
                                             mass_of_cannister,
                                             total_coated_area_cathode,         # m2
                                             pos_AM_wt_fr,
                                             neg_AM_wt_fr
                                             )

    return cell_design_op
end

max_issues_repo_name: battmodels/BatteryCost

max_forks_repo_name: battmodels/BatteryCost


text:
function get_delay_hrs(NBITS, b, db, rs, prec, el=true)
    setprecision(NBITS)
    prec = 1.0e-14

    massfunc(r) = tnfw_m(0.1, 500.0, 20.0, r);
    rhofunc(r) = tnfw_ρ(0.1, 500.0, 20.0, r);
    tg = b/rs
    up_ex, up_ey, up_eu, up_edu, up_eia, up_eoa, up_eda, up_eθ =
    el? numeric_el(massfunc, rhofunc, 0.0, rs, tg, 100000; reltol=prec, abstol=prec) :
        numeric(massfunc, 0.0, rs, tg, 100000; reltol=prec, abstol=prec);
    b += db
    tg = b/rs
    down_ex, down_ey, down_eu, down_edu, down_eia, down_eoa, down_eda, down_eθ =
    el? numeric_el(massfunc, rhofunc, 0.0, rs, tg, 100000; reltol=prec, abstol=prec) :
        numeric(massfunc, 0.0, rs, tg, 100000; reltol=prec, abstol=prec);

    uex,uey,dex,dey = extrapolate_to_meeting(up_ex,up_ey,down_ex,down_ey);

    uet = time(massfunc, uex, uey);
    det = time(massfunc, dex, dey);
    delta_e = abs(uet-det);
    display(string(NBITS)*", "*string(db)*", "*string(delta_e))
    delta_e*1e9*365*24;
end

ext: jl

max_stars_repo_name: skariel/Lensing

max_issues_repo_name: skariel/Lensing

max_forks_repo_name: skariel/Lensing


text:
using Interact, Blink, Colors
export interactive_orbitdiagram, scaleod

"""
    controlwindow(D, n, Ttr, density, i)
Create an Electron control window for the orbit diagram interactive application.

```julia
return n, Ttr, density, i, ▢update, ▢back, ▢reset, α,
       ⬜pmin, ⬜pmax, ⬜umin, ⬜umax
```

All returned values are `Observable`s.
Their value corresponds to the one chosen in the Electron window.
Items with `▢` are buttons and with `⬜` are the boxes with limits.
"""
function controlwindow(D, n, Ttr, density, i)
    n = Interact.textbox(string(n); value = n, label = "n")
    Ttr = Interact.textbox(string(Ttr); value = Ttr, label = "Ttr")
    density = Interact.textbox(string(density); value = density, label = "density")

    i = Interact.dropdown(OrderedDict(string(j) => j for j in 1:D); label = "variable")
    ▢update = Interact.button("update")
    ▢back = Interact.button("← back")
    ▢reset = Interact.button("reset")

    α = Interact.slider(0:0.001:1; value = 0.1, label = "α (transparency)")

    # Limit boxes
    ⬜pmin = Interact.textbox(; value = 0.0, label = "pmin")
    ⬜pmax = Interact.textbox(; value = 1.0, label = "pmax")
    ⬜umin = Interact.textbox(; value = 0.0, label = "umin")
    ⬜umax = Interact.textbox(; value = 1.0, label = "umax")

    w = Window(Dict(:height => 400, :title => "Orbit Diagram controls"))
    s = 5em
    body!(w, Interact.vbox(
        α,
        Interact.hbox(n, Ttr, density),
        Interact.hbox(i, hskip(s), ▢update, hskip(s), ▢back, hskip(s), ▢reset),
        Interact.hline(),
        Interact.hbox(⬜pmin, ⬜pmax),
        Interact.hbox(⬜umin, ⬜umax)
        )
    )

    return n, Ttr, density, i, ▢update, ▢back, ▢reset, observe(α),
           ⬜pmin, ⬜pmax, ⬜umin, ⬜umax
end

"""
    interactive_orbitdiagram(ds::DiscreteDynamicalSystem,
        i::Int, p_index, p_min, p_max;
        density = 500, u0 = get_state(ds), Ttr = 200, n = 500,
        parname = "p"
    )

Open an interactive application for exploring orbit diagrams (ODs) of
discrete systems. The functionality works for _any_ discrete system.

Once initialized it opens a Makie plot window and an Electron control window.

## Interaction
By using the Electron window you are able to update all parameters of the OD
interactively (like e.g. `n` or `Ttr`). You have to press `update` after changing
these parameters. You can even decide which variable to get the OD for,
by choosing one of the variables from the wheel (again, press `update` afterwards).

In the Makie window you can interactively zoom into the OD. Click
and drag with the left mouse button to select a region in the OD. This region is then
**re-computed** at a higher resolution (i.e. we don't "just zoom").

Back in the Electron window, you can press `reset` to bring the OD in the original
state (and variable). Pressing `back` will go back through the history of your exploration
History is stored when any change happens (besides transparency).

## Accessing the data
What is plotted on the application window is a _true_ orbit diagram, not a plotting
shorthand. This means that all data are obtainable and usable directly.
Internally we always scale the orbit diagram to [0,1]² (to allow `Float64` precision
even though plotting is `Float32`-based). This however means that it is
necessary to transform the data in real scale. This is done through the function
[`scaleod`](@ref) which accepts the 5 arguments returned from the current function:
```
od, pmin, pmax, umin, umax = interactive_orbitdiagram(...)
ps, us = scaleod(od, pmin, pmax, umin, umax)
```
"""
function interactive_orbitdiagram(ds::DiscreteDynamicalSystem,
    i::Int, p_index, p_min, p_max;
    density = 1000, u0 = get_state(ds), Ttr = 200, n = 500,
    parname = "p"
    )

    # UI elements
    n, Ttr, density, i, ▢update, ▢back, ▢reset, α, ⬜pmin, ⬜pmax, ⬜umin, ⬜umax =
    controlwindow(dimension(ds), n, Ttr, density, i)

    # Initial Orbit diagram data
    integ = integrator(ds, u0)
    pmin, pmax = p_min, p_max
    odinit, xmin, xmax = minimal_normalized_od(integ, i[], p_index, pmin, pmax, density[], n[], Ttr[], u0)
    od_node = Observable(odinit)
    densityinit = density[]; ninit = n[]; Ttrinit = Ttr[]

    # History stores the variable index and true diagram limits
    history = [(i[], pmin, pmax, xmin, xmax, n[], Ttr[], density[])]
    update_controls!(history[end], i, n, Ttr, density,  ⬜pmin, ⬜pmax, ⬜umin, ⬜umax)

    color = lift(a -> RGBA(0,0,0,a), α)
    scplot = Scene(resolution = (1200, 800))
    scatter!(scplot, od_node, markersize = 0.008, color = color)

    scplot[Axis][:ticks][:ranges] = (collect(0:0.25:1), collect(0:0.25:1))
    scplot[Axis][:ticks][:labels] =
    (["pmin", " ", " ", " ", "pmax"], ["umin", " ", " ", " ", "umax"])
    scplot[Axis][:names][:axisnames] = (parname, "u"*subscript(i[]))

    display(scplot)
    rect = select_rectangle(scplot)

    # Uppon interactively selecting a rectangle, with value `r` (in [0,1]²)
    on(rect) do r
        spmin, sxmin = r.origin
        spmax, sxmax = r.origin + r.widths
        # Convert p,x to true values
        j, ppmin, ppmax, pxmin, pxmax = history[end]
        pdif = ppmax - ppmin; xdif = pxmax - pxmin
        pmin = spmin*pdif + ppmin
        pmax = spmax*pdif + ppmin
        xmin = sxmin*xdif + pxmin
        xmax = sxmax*xdif + pxmin

        od_node[] = minimal_normalized_od(
            integ, j,  p_index, pmin, pmax,
            density[], n[], Ttr[], u0, xmin, xmax
        )
        # update history and controls
        push!(history, (j, pmin, pmax, xmin, xmax, n[], Ttr[], density[]))
        update_controls!(history[end], i, n, Ttr, density,  ⬜pmin, ⬜pmax, ⬜umin, ⬜umax)
    end

    # Upon hitting the update button
    on(▢update) do clicks
        j, pmin, pmax, xmin, xmax, m, T, d = history[end]
        # Check if there was any change:
        if !(⬜pmin[] == pmin && ⬜pmax[] == pmax && ⬜umin[] == xmin && ⬜umax[] == xmax &&
            j == i[] && m == n[] && T == Ttr[] && d == density[])

            pmin, pmax, xmin, xmax = ⬜pmin[], ⬜pmax[], ⬜umin[], ⬜umax[]
            j, m, T, d = i[], n[], Ttr[], density[]

            od_node[] = minimal_normalized_od(
            integ, j, p_index, pmin, pmax, density[], n[], Ttr[], u0, xmin, xmax
            )
            # Update history and controls
            push!(history, (j, pmin, pmax, xmin, xmax, m, T, d))
            update_controls!(history[end], i, n, Ttr, density,  ⬜pmin, ⬜pmax, ⬜umin, ⬜umax)
            scplot[Axis][:names][:axisnames] = (parname, "u"*subscript(i[]))
        end
    end

    # Upon hitting the "reset" button
    on(▢reset) do clicks
        if length(history) > 1
            deleteat!(history, 2:length(history))
            j, pmin, pmax, xmin, xmax, m, T, d = history[end]
            od_node[] = odinit
            # Update controls/labels
            update_controls!(history[end], i, n, Ttr, density,  ⬜pmin, ⬜pmax, ⬜umin, ⬜umax)
            scplot[Axis][:names][:axisnames] = (parname, "u"*subscript(j))
        end
    end

    # Upon hitting the "back" button
    on(▢back) do clicks
        if length(history) > 1
            pop!(history)
            j, pmin, pmax, xmin, xmax, m, T, d = history[end]
            od_node[] = minimal_normalized_od(
                integ, j, p_index, pmin, pmax, d, m, T, u0, xmin, xmax
            )
            # Update limits in textboxes
            update_controls!(history[end], i, n, Ttr, density,  ⬜pmin, ⬜pmax, ⬜umin, ⬜umax)
            scplot[Axis][:names][:axisnames] = (parname, "u"*subscript(j))
        end
    end

    display(scplot)
    return od_node, ⬜pmin, ⬜pmax, ⬜umin, ⬜umax
end

"""
    minimal_normalized_od(integ, i, p_index, pmin, pmax,
                         density, n, Ttr, u0)
    minimal_normalized_od(integ, i, p_index, pmin, pmax,
                         density, n, Ttr, u0, xmin, xmax)

Compute and return a minimal and normalized orbit diagram (OD).

All points are stored in a single vector of `Point2f0` to ensure fastest possible
plotting. In addition all numbers are scaled to [0, 1]. This allows us to have
64-bit precision while display is only 32-bit!

The version with `xmin, xmax` only keeps points with limits between the
real `xmin, xmax` (in the normal units of the dynamical system).
"""
function minimal_normalized_od(integ, i, p_index, pmin, pmax,
                              density::Int, n::Int, Ttr::Int, u0)

    pvalues = range(pmin, stop = pmax, length = density)
    pdif = pmax - pmin
    od = Vector{Point2f0}() # make this pre-allocated
    xmin = eltype(integ.u)(Inf); xmax = eltype(integ.u)(-Inf)
    #= @inbounds =# for (j, p) in enumerate(pvalues)
        pp = (p - pmin)/pdif # p to plot, in [0, 1]
        DynamicalSystemsBase.reinit!(integ, u0)
        integ.p[p_index] = p
        DynamicalSystemsBase.step!(integ, Ttr)
        for z in 1:n
            DynamicalSystemsBase.step!(integ)
            x = integ.u[i]
            push!(od, Point2f0(pp, integ.u[i]))
            # update limits
            if x < xmin
                xmin = x
            elseif x > xmax
                xmax = x
            end
        end
    end
    # normalize x values to [0, 1]
    xdif = xmax - xmin
    #= @inbounds =# for j in eachindex(od)
        x = od[j][2]; p = od[j][1]
        od[j] = Point2f0(p, (x - xmin)/xdif)
    end
    return od, xmin, xmax
end

function minimal_normalized_od(integ, i, p_index, pmin, pmax,
                              density::Int, n::Int, Ttr::Int, u0, xmin, xmax)

    pvalues = range(pmin, stop = pmax, length = density)
    pdif = pmax - pmin; xdif = xmax - xmin
    od = Vector{Point2f0}()
    #= @inbounds =# for p in pvalues
        pp = (p - pmin)/pdif # p to plot, in [0, 1]
        DynamicalSystemsBase.reinit!(integ, u0)
        integ.p[p_index] = p
        DynamicalSystemsBase.step!(integ, Ttr)
        for z in 1:n
            DynamicalSystemsBase.step!(integ)
            x = integ.u[i]
            if xmin ≤ x ≤ xmax
                push!(od, Point2f0(pp, (integ.u[i] - xmin)/xdif))
            end
        end
    end
    return od
end

function  update_controls!(h, i, n, Ttr, density, ⬜pmin, ⬜pmax, ⬜umin, ⬜umax)
    j, pmin, pmax, xmin, xmax, m, T, d = h
    i[] = j; n[] = m; Ttr[] = T; density[] = d
    ⬜pmin[] = pmin; ⬜pmax[] = pmax
    ⬜umin[] = xmin; ⬜umax[] = xmax
    return
end

"""
    scaleod(od, pmin, pmax, umin, umax) -> ps, us
Given the return values of [`interactive_orbitdiagram`](@ref), produce
orbit diagram data scaled correctly in data units. Return the data as a vector of
parameter values and a vector of corresponding variable values.
"""
function scaleod(od, pmin, pmax, umin, umax)
    oddata = od[]; L = length(oddata);
    T = promote_type(typeof(umin[]), Float32)
    ps = zeros(T, L); us = copy(ps)
    udif = umax[] - umin[]; um = umin[]
    pdif = pmax[] - pmin[]; pm = pmin[]
    @inbounds for i ∈ 1:length(oddata)
        p, u = oddata[i]
        ps[i] = pm + pdif*p; us[i] = um + udif*u
    end
    return ps, us
end

# TODO: Use `α` directly after Simon's fix, no need for observe(α)
# TODO: Add exit button that closes both windows
# TODO: Make marker GLMakie.FastPixel()
# TODO: Allow initial state to be a function of paramter (define function `get_u(f, p)`)


text:
"""
Chain all arguments together in a `next` chain. to form a new structure `next` chain.

If `nexts` is empty, `C_NULL` is returned.
"""
function chain(nexts::HighLevelStruct...)
    isempty(nexts) && return C_NULL
    foldr((a, b) -> (@set a.next = b), nexts)
end

"""
Break a `next` chain into its constituents, with all `next` members set to `C_NULL`.
"""
function unchain(x::HighLevelStruct)
    unchained = Any[@set x.next = C_NULL]
    (; next) = x
    while next ≠ C_NULL
        push!(unchained, @set next.next = C_NULL)
        (; next) = next
    end
    unchained
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
using MultivariateStats
using Base.Test

srand(34568)

## data

m = 9
n = 6
n2 = 3

X = randn(m, n)
A = randn(n, n2)
Xt = X'

b = randn(1, n2)

E = randn(m, n2) * 0.1
Y0 = X * A + E
Y1 = X * A .+ b + E

y0 = Y0[:,1]
y1 = Y1[:,1]

## llsq

A = llsq(X, Y0; trans=false, bias=false)
A_r = copy(A)
@test size(A) == (n, n2)
@test_approx_eq X'Y0 X'X * A

a = llsq(X, y0; trans=false, bias=false)
@test size(a) == (n,)
@test_approx_eq a A[:,1]

A = llsq(Xt, Y0; trans=true, bias=false)
@test size(A) == (n, n2)
@test_approx_eq A A_r

a = llsq(Xt, y0; trans=true, bias=false)
@test size(a) == (n,)
@test_approx_eq a A[:,1]

Aa = llsq(X, Y1; trans=false, bias=true)
Aa_r = copy(Aa)
@test size(Aa) == (n+1, n2)
A, b = Aa[1:end-1,:], Aa[end,:]
@test_approx_eq X' * (Y1 .- b) X'X * A

aa = llsq(X, y1; trans=false, bias=true)
@test_approx_eq aa Aa[:,1]

Aa = llsq(Xt, Y1; trans=true, bias=true)
@test_approx_eq Aa Aa_r

aa = llsq(Xt, y1; trans=true, bias=true)
@test_approx_eq aa Aa[:,1]

## ridge (with Real r)

r = 0.1

A = ridge(X, Y0, r; trans=false, bias=false)
A_r = copy(A)
@test size(A) == (n, n2)
@test_approx_eq X'Y0 (X'X + r * eye(n)) * A

a = ridge(X, y0, r; trans=false, bias=false)
@test size(a) == (n,)
@test_approx_eq a A[:,1]

A = ridge(Xt, Y0, r; trans=true, bias=false)
@test size(A) == (n, n2)
@test_approx_eq A A_r

a = ridge(Xt, y0, r; trans=true, bias=false)
@test size(a) == (n,)
@test_approx_eq a A[:,1]

Aa = ridge(X, Y1, r; trans=false, bias=true)
Aa_r = copy(Aa)
@test size(Aa) == (n+1, n2)
A, b = Aa[1:end-1,:], Aa[end,:]
@test_approx_eq X' * (Y1 .- b) (X'X + r * eye(n)) * A

aa = ridge(X, y1, r; trans=false, bias=true)
@test_approx_eq aa Aa[:,1]

Aa = ridge(Xt, Y1, r; trans=true, bias=true)
@test_approx_eq Aa Aa_r

aa = ridge(Xt, y1, r; trans=true, bias=true)
@test_approx_eq aa Aa[:,1]

## ridge (with diagonal r)

r = 0.05 + 0.1 * rand(n)

A = ridge(X, Y0, r; trans=false, bias=false)
A_r = copy(A)
@test size(A) == (n, n2)
@test_approx_eq X'Y0 (X'X + diagm(r)) * A

a = ridge(X, y0, r; trans=false, bias=false)
@test size(a) == (n,)
@test_approx_eq a A[:,1]

A = ridge(Xt, Y0, r; trans=true, bias=false)
@test size(A) == (n, n2)
@test_approx_eq A A_r

a = ridge(Xt, y0, r; trans=true, bias=false)
@test size(a) == (n,)
@test_approx_eq a A[:,1]

Aa = ridge(X, Y1, r; trans=false, bias=true)
Aa_r = copy(Aa)
@test size(Aa) == (n+1, n2)
A, b = Aa[1:end-1,:], Aa[end,:]
@test_approx_eq X' * (Y1 .- b) (X'X + diagm(r)) * A

aa = ridge(X, y1, r; trans=false, bias=true)
@test_approx_eq aa Aa[:,1]

Aa = ridge(Xt, Y1, r; trans=true, bias=true)
@test_approx_eq Aa Aa_r

aa = ridge(Xt, y1, r; trans=true, bias=true)
@test_approx_eq aa Aa[:,1]

## ridge (with qudratic r matrix)

Q = qr(randn(n, n))[1]
r = Q' * diagm(r) * Q

A = ridge(X, Y0, r; trans=false, bias=false)
A_r = copy(A)
@test size(A) == (n, n2)
@test_approx_eq X'Y0 (X'X + r) * A

a = ridge(X, y0, r; trans=false, bias=false)
@test size(a) == (n,)
@test_approx_eq a A[:,1]

A = ridge(Xt, Y0, r; trans=true, bias=false)
@test size(A) == (n, n2)
@test_approx_eq A A_r

a = ridge(Xt, y0, r; trans=true, bias=false)
@test size(a) == (n,)
@test_approx_eq a A[:,1]

Aa = ridge(X, Y1, r; trans=false, bias=true)
Aa_r = copy(Aa)
@test size(Aa) == (n+1, n2)
A, b = Aa[1:end-1,:], Aa[end,:]
@test_approx_eq X' * (Y1 .- b) (X'X + r) * A

aa = ridge(X, y1, r; trans=false, bias=true)
@test_approx_eq aa Aa[:,1]

Aa = ridge(Xt, Y1, r; trans=true, bias=true)
@test_approx_eq Aa Aa_r

aa = ridge(Xt, y1, r; trans=true, bias=true)
@test_approx_eq aa Aa[:,1]


text:
module ImageIO
using ColorTypes, Compat, FileIO
include("types.jl")

export Image
export AbstractImage
using ImageMagick
importall ImageMagick, FileIO

#Define backends
FileIO.read(file::readformats(ImageMagick.BACKEND)) = read(file, ImageMagick.BACKEND)

end


text:
using DistStat, LinearAlgebra,Random,Test, Pkg

type=[Float64,Float32]

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    using CUDA
    A= CuArray
else
    A = Array
end

for T in type
    a=A{T}(reshape(collect(1:81),9,9))
    a_dist=distribute(a)

    @test isapprox(opnorm(a),opnorm(a_dist))

end


text:
const PolicyOrSelector = Union{AbstractPolicy, AbstractActionSelector}

"""
    TDLearner(approximator::Tapp, π::Tp, γ::Float64, α::Float64, n::Int=0) where {Tapp<:AbstractVApproximator, Tp<:PolicyOrSelector} = new{Tapp, Tp, :SRS}(approximator, π, γ, α, n)
    TDLearner(approximator::Tapp, π::Tp, γ::Float64, α::Float64, n::Int=0, method::Symbol=:SARSA) where {Tapp<:AbstractQApproximator, Tp<:PolicyOrSelector}

See more details at Section (7.1) on Page 142 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*
"""
struct TDLearner{Tapp <: AbstractApproximator, Tp <: PolicyOrSelector, method} <: AbstractModelFreeLearner
    approximator::Tapp
    π::Tp
    γ::Float64
    α::Float64
    n::Int
    TDLearner(approximator::Tapp, π::Tp, γ::Float64, α::Float64, n::Int=0) where {Tapp<:AbstractVApproximator, Tp<:PolicyOrSelector} = new{Tapp, Tp, :SRS}(approximator, π, γ, α, n)
    function TDLearner(approximator::Tapp, π::Tp, γ::Float64, α::Float64, n::Int=0, method::Symbol=:SARSA) where {Tapp<:AbstractQApproximator, Tp<:PolicyOrSelector}
        supported_methods = [:SARSA, :ExpectedSARSA]
        !in(method, supported_methods) && throw(ArgumentError("Supported methods are $supported_methods , your input is $method"))
        new{Tapp, Tp, method}(approximator, π, γ, α, n)
    end
end

(learner::TDLearner{<:AbstractQApproximator, <:AbstractActionSelector})(s) = learner.approximator(s) |> learner.π
(learner::TDLearner{<:AbstractApproximator, <:AbstractPolicy})(s) = learner.π(s)
(learner::TDLearner{<:AbstractApproximator, <:AbstractPolicy})(s, ::Val{:dist}) = learner.π(s, Val(:dist))

function update!(learner::TDLearner{<:AbstractQApproximator, <:PolicyOrSelector, :SARSA}, buffer::EpisodeSARDBuffer)
    n = learner.n
    update!(learner,
            @view(buffer.state[max(1, end - n - 1) : end - 1]),
            @view(buffer.action[max(1, end - n - 1) : end - 1]),
            @view(buffer.reward[max(1, end - n) : end]),
            @view(buffer.state[max(1, end - n) : end]),
            @view(buffer.action[max(1, end - n) : end]),
            Val(buffer.isdone[end]))
end

function update!(learner::TDLearner{<:AbstractQApproximator, <:PolicyOrSelector, :ExpectedSARSA}, buffer::EpisodeSARDBuffer)
    n = learner.n
    update!(learner,
            @view(buffer.state[max(1, end - n - 1) : end - 1]),
            @view(buffer.action[max(1, end - n - 1) : end - 1]),
            @view(buffer.reward[max(1, end - n) : end]),
            @view(buffer.state[max(1, end - n) : end]),
            Val(buffer.isdone[end]))
end

function update!(learner::TDLearner{<:AbstractVApproximator, <:PolicyOrSelector, :SRS}, buffer::EpisodeSARDBuffer)
    n = learner.n
    update!(learner,
            @view(buffer.state[max(1, end - n - 1) : end - 1]),
            @view(buffer.reward[max(1, end - n) : end]),
            @view(buffer.state[max(1, end - n) : end]),
            Val(buffer.isdone[end]))
end

function update!(learner::TDLearner{<:AbstractVApproximator,<:PolicyOrSelector, :SRS}, states, rewards, nextstates, ::Val{true})
    n, γ, V, α = learner.n, learner.γ, learner.approximator, learner.α
    # Warning!!! The order of calculation is reversed here for speed. The impact is uncertain!!!
    for (G, s) in zip(reverse_discounted_rewards(rewards, γ), Iterators.reverse(states))
        update!(V, s, α * (G - V(s)))
    end
end

function update!(learner::TDLearner{<:AbstractVApproximator,<:PolicyOrSelector, :SRS}, states, rewards, nextstates, ::Val{false})
    n, γ, V, α = learner.n, learner.γ, learner.approximator, learner.α
    if length(states) ≥ n
        G = discounted_reward(rewards, γ) + γ^n * V(nextstates[end])
        s = states[1]
        update!(V, s, α * (G - V(s)))
    end
end

function update!(learner::TDLearner{<:AbstractQApproximator, <:PolicyOrSelector, :SARSA}, states, actions, rewards, nextstates, nextactions, ::Val{true})
    n, γ, Q, α = learner.n, learner.γ, learner.approximator, learner.α
    # Warning!!! The order of calculation is reversed here for speed. The impact is uncertain!!!
    for (G, s, a) in zip(reverse_discounted_rewards(rewards, γ), Iterators.reverse(states), Iterators.reverse(actions))
        update!(Q, s, a, α * (G - Q(s, a)))
        update!(learner.π, s, Q(s, Val(:argmax)))
    end
end

function update!(learner::TDLearner{<:AbstractQApproximator, <:PolicyOrSelector, :SARSA}, states, actions, rewards, nextstates, nextactions, ::Val{false})
    n, γ, Q, α = learner.n, learner.γ, learner.approximator, learner.α
    if length(states) ≥ n
        G = discounted_reward(rewards, γ) + γ^n * Q(nextstates[end], nextactions[end])
        s, a = states[1], actions[1]
        update!(Q, s, a, α * (G - Q(s, a)))
        update!(learner.π, s, Q(s, Val(:argmax)))
    end
end

function update!(learner::TDLearner{<:AbstractQApproximator, <:PolicyOrSelector, :ExpectedSARSA}, states, actions, rewards, nextstates, ::Val{true})
    n, γ, Q, α = learner.n, learner.γ, learner.approximator, learner.α
    # Warning!!! The order of calculation is reversed here for speed. The impact is uncertain!!!
    for (G, s, a) in zip(reverse_discounted_rewards(rewards, γ), Iterators.reverse(states), Iterators.reverse(actions))
        update!(Q, s, a, α * (G - Q(s, a)))
        update!(learner.π, s, Q(s, Val(:argmax)))
    end
end

function update!(learner::TDLearner{<:AbstractQApproximator, <:PolicyOrSelector, :ExpectedSARSA}, states, actions, rewards, nextstates, ::Val{false})
    n, γ, Q, α, π = learner.n, learner.γ, learner.approximator, learner.α, learner.π
    if length(states) ≥ n
        G = discounted_reward(rewards, γ) + γ^n * sum(Q(nextstates[end], Val(:dist)) .* π(nextstates[end], Val(:dist)))
        s, a = states[1], actions[1]
        update!(Q, s, a, α * (G - Q(s, a)))
        update!(learner.π, s, Q(s, Val(:argmax)))
    end
end

"""

    struct TDλReturnLearner{Tapp <: AbstractApproximator, Tp <: PolicyOrSelector} <: AbstractModelFreeLearner
        approximator::Tapp
        π::Tp
        γ::Float64
        α::Float64
        λ::Float64
    end

See more details at Section (12.2) on Page 292 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*
"""
struct TDλReturnLearner{Tapp <: AbstractApproximator, Tp <: PolicyOrSelector} <: AbstractModelFreeLearner
    approximator::Tapp
    π::Tp
    γ::Float64
    α::Float64
    λ::Float64
end

(learner::TDλReturnLearner)(s) = learner.π(s)
(learner::TDλReturnLearner)(s, ::Val{:dist}) = learner.π(s, Val(:dist))

function update!(learner::TDλReturnLearner{<:AbstractVApproximator}, buffer::EpisodeSARDBuffer)
    λ, γ, V, α, T = learner.λ, learner.γ, learner.approximator, learner.α, length(buffer.state)
    if isfull(buffer)
        for t in 1:(T-1)
            G = 0.
            for n in 1:(T - t - 1)
                G += λ^(n-1) * (discounted_reward(@view(buffer.reward[t:t+n-1]), γ) + γ^n * V(buffer.state[t+n]))
            end
            G *= 1 - λ
            G += λ^(T-t-1) * (discounted_reward(@view(buffer.reward[t:T-1]), γ) + γ^(T-t) * V(buffer.state[T]))
            sₜ = buffer.state[t]
            update!(V, sₜ, α * (G - V(sₜ)))
        end
    end
end

"""
    OffPolicyTDLearner(approximator::Tapp, π_behavior::Tpb, π_target::Tpt, γ::Float64, α::Float64, n::Int=0, method::Symbol=:SARSA_ImportanceSampling) where {Tapp<:AbstractApproximator, Tpb<:AbstractPolicy, Tpt<:AbstractPolicy}

See more details at Section (7.3) on Page 148 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*
"""
struct OffPolicyTDLearner{Tapp <: AbstractApproximator, Tpb <: PolicyOrSelector, Tpt <: PolicyOrSelector, Tα<:Union{Function, Float64}, method} <: AbstractModelFreeLearner
    approximator::Tapp
    π_behavior::Tpb
    π_target::Tpt
    γ::Float64
    α::Tα
    n::Int
    function OffPolicyTDLearner(approximator::Tapp, π_behavior::Tpb, π_target::Tpt, γ::Float64, α::Tα, n::Int=0, method::Symbol=:SARSA_ImportanceSampling) where {Tapp<:AbstractApproximator, Tpb<:PolicyOrSelector, Tpt<:PolicyOrSelector, Tα<:Union{Function, Float64}}
        new{Tapp, Tpb, Tpt, Tα, method}(approximator, π_behavior, π_target, γ, α, n)
    end
end

const QLearner = OffPolicyTDLearner{Tapp, Tp, Tp, Tα, :QLearning} where {Tapp<:AbstractQApproximator, Tp<:PolicyOrSelector, Tα<:Union{Function, Float64}}
QLearner(Q::AbstractQApproximator, π::PolicyOrSelector, γ::Float64, α::Union{Function, Float64}) = OffPolicyTDLearner(Q, π, π, γ, α, 0, :QLearning)

(learner::OffPolicyTDLearner{<:AbstractApproximator, <:AbstractActionSelector})(s) = learner.approximator(s) |> learner.π_behavior
(learner::OffPolicyTDLearner{<:AbstractApproximator, <:AbstractPolicy})(s) = learner.π_behavior(s)
(learner::OffPolicyTDLearner{<:AbstractApproximator, <:AbstractPolicy})(s, ::Val{:dist}) = learner.π_behavior(s, Val(:dist))

priority(learner::QLearner{<:AbstractQApproximator, <:AbstractPolicy, <:Float64}, s, a, r, d, s′) = priority(learner, learner.α, s, a, r, d, s′)
priority(learner::QLearner{<:AbstractQApproximator, <:AbstractPolicy, <:Function}, s, a, r, d, s′) = priority(learner, learner.α((s,a)), s, a, r, d, s′)
function priority(learner::QLearner, α, s, a, r, d, s′)
    γ, Q = learner.γ, learner.approximator
    error = d ? α * (r - Q(s, a)) : α * (r + γ * Q(s′, Val(:max)) - Q(s, a))
    abs(error)
end

function update!(learner::OffPolicyTDLearner{<:AbstractApproximator, <:PolicyOrSelector, <:PolicyOrSelector, <:Union{Function, Float64}, :SARSA_ImportanceSampling}, buffer::EpisodeSARDBuffer)
    n = learner.n
    update!(learner,
            @view(buffer.state[max(1, end - n - 1) : end - 1]),
            @view(buffer.action[max(1, end - n - 1) : end - 1]),
            @view(buffer.reward[max(1, end - n) : end]),
            @view(buffer.state[max(1, end - n) : end]),
            @view(buffer.action[max(1, end - n) : end]),
            Val(buffer.isdone[end]))
end

update!(learner::QLearner{<:AbstractQApproximator, <:PolicyOrSelector, <:Float64}, s, a, r, d, s′) = update!(learner, learner.α, s, a, r, d, s′)
update!(learner::QLearner{<:AbstractQApproximator, <:PolicyOrSelector, <:Function}, s, a, r, d, s′) = update!(learner, learner.α((s, a)), s, a, r, d, s′)
function update!(learner::QLearner, α, s, a, r, d, s′)
    Q, γ, π = learner.approximator, learner.γ, learner.π_target
    error = d ? α * (r - Q(s, a)) : α * (r + γ * Q(s′, Val(:max)) - Q(s, a))
    update!(Q, s, a, error)
    update!(π, s, Q(s, Val(:argmax)))
end

function update!(learner::QLearner, buffer::EpisodeSARDBuffer)
    s, a, r, d, s′ = buffer.state[end-1], buffer.action[end-1], buffer.reward[end], buffer.isdone[end], buffer.state[end]
    update!(learner, s, a, r, d, s′)
end

function update!(learner::QLearner, buffer::EpisodeSARDSBuffer)
    update!(learner, buffer.state[end], buffer.action[end], buffer.reward[end], buffer.isdone[end], buffer.nextstate[end])
end

function update!(learner::OffPolicyTDLearner{<:AbstractQApproximator, <:AbstractPolicy,  <:AbstractPolicy, <:Union{Function, Float64}, :SARSA_ImportanceSampling}, states, actions, rewards, nextstates, nextactions, ::Val{true})
    n, γ, Q, α, π, b = learner.n, learner.γ, learner.approximator, learner.α, learner.π_target, learner.π_behavior
    # Warning!!! The order of calculation is reversed here for speed. The impact is uncertain!!!
    for (G, ρ, s, a) in zip(reverse_discounted_rewards(rewards, γ),
                            reverse_importance_weights(π, b, @view(states[2:end]), @view(actions[2:end])),
                            Iterators.reverse(states),
                            Iterators.reverse(actions))
        update!(Q, s, a, α * ρ * (G - Q(s, a)))
        update!(π, s, Q(s, Val(:argmax)))
    end
end

function update!(learner::OffPolicyTDLearner{<:AbstractVApproximator, <:AbstractPolicy,  <:AbstractPolicy, <:Union{Function, Float64}, :SARSA_ImportanceSampling}, states, actions, rewards, nextstates, nextactions, ::Val{true})
    n, γ, V, α, π, b = learner.n, learner.γ, learner.approximator, learner.α, learner.π_target, learner.π_behavior
    # Warning!!! The order of calculation is reversed here for speed. The impact is uncertain!!!
    for (G, ρ, s, a) in zip(reverse_discounted_rewards(rewards, γ),
                            reverse_importance_weights(π, b, @view(states[2:end]), @view(actions[2:end])),
                            Iterators.reverse(states))
        update!(V, s, α * ρ * (G - V(s)))
    end
end

function update!(learner::OffPolicyTDLearner{<:AbstractQApproximator, <:AbstractPolicy,  <:AbstractPolicy, <:Union{Function, Float64}, :SARSA_ImportanceSampling}, states, actions, rewards, nextstates, nextactions, ::Val{false})
    n, γ, Q, α, π, b = learner.n, learner.γ, learner.approximator, learner.α, learner.π_target, learner.π_behavior
    if length(states) ≥ n
        G = discounted_reward(rewards, γ) + γ^n * Q(nextstates[end], nextactions[end])
        ρ = importance_weight(π, b, @view(states[2:end]), @view(actions[2:end]))
        s, a = states[1], actions[1]
        update!(Q, s, a, α * ρ * (G - Q(s, a)))
        update!(π, s, Q(s, Val(:argmax)))
    end
end

function update!(learner::OffPolicyTDLearner{<:AbstractVApproximator, <:AbstractPolicy,  <:AbstractPolicy, <:Union{Function, Float64}, :SARSA_ImportanceSampling}, states, actions, rewards, nextstates, nextactions, ::Val{false})
    n, γ, V, α, π, b = learner.n, learner.γ, learner.approximator, learner.α, learner.π_target, learner.π_behavior
    if length(states) ≥ n
        G = discounted_reward(rewards, γ) + γ^n * V(nextstates[end])
        ρ = importance_weight(π, b, @view(states[2:end]), @view(actions[2:end]))
        s = states[1]
        update!(V, s, α * ρ * (G - V(s)))
    end
end

"""
    struct DoubleLearner{Tl <: OffPolicyTDLearner, Ts<:AbstractActionSelector} <: AbstractLearner
        Learner1::Tl
        Learner2::Tl
        selector::Ts
    end

See more details at Section (6.7) on Page 126 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*
"""
struct DoubleLearner{Tl <: OffPolicyTDLearner, Ts<:AbstractActionSelector} <: AbstractLearner
    Learner1::Tl
    Learner2::Tl
    selector::Ts
end

(learner::DoubleLearner)(s, ::Val{:dist}) = learner.Learner1(s, Val(:dist)) .+ learner.Learner2(s, Val(:dist))
(learner::DoubleLearner)(s) =  learner.selector(learner(s, Val(:dist)))

function update!(learner::DoubleLearner{<:OffPolicyTDLearner{<:AbstractQApproximator, <:AbstractPolicy, <:AbstractPolicy, <:Union{Function, Float64}, :QLearning}}, buffer::EpisodeSARDBuffer)
    s, a, r, s′ = buffer.state[end-1], buffer.action[end-1], buffer.reward[end], buffer.state[end]

    if rand() < 0.5
        L1, L2 = learner.Learner1, learner.Learner2
    else
        L2, L1 = learner.Learner1, learner.Learner2
    end

    if !buffer.isdone[end]
        update!(L1.approximator, s, a, L1.α * (r + L1.γ * L2.approximator(s′, L1.approximator(s′, Val(:argmax))) - L1.approximator(s, a)))
    else
        update!(L1.approximator, s, a, L1.α * (r - L1.approximator(s, a)))
    end

    update!(L1.π_target, s, L1.approximator(s, Val(:argmax)))
end

"""
    DifferentialTDLearner(approximator::Tapp, π::Tp, α::Float64, β::Float64, R̄::Float64=0., n::Int=0, method::Symbol=:SARSA) where {Tapp<:AbstractApproximator, Tp<:PolicyOrSelector}= new{Tapp, Tp, method}(approximator, π, α, β, R̄, n)

See more details at Section (10.3) on Page 251 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*
"""
mutable struct DifferentialTDLearner{Tapp<:AbstractApproximator, Tp<:PolicyOrSelector, method} <: AbstractLearner
    approximator::Tapp
    π::Tp
    α::Float64
    β::Float64
    R̄::Float64
    n::Int
    DifferentialTDLearner(approximator::Tapp, π::Tp, α::Float64, β::Float64, R̄::Float64=0., n::Int=0, method::Symbol=:SARSA) where {Tapp<:AbstractApproximator, Tp<:PolicyOrSelector}= new{Tapp, Tp, method}(approximator, π, α, β, R̄, n)
end

(learner::DifferentialTDLearner{<:AbstractQApproximator, <:AbstractActionSelector})(s) = learner.approximator(s) |> learner.π
(learner::DifferentialTDLearner{<:AbstractApproximator, <:AbstractPolicy})(s) = learner.π(s)
(learner::DifferentialTDLearner{<:AbstractApproximator, <:AbstractPolicy})(s, ::Val{:dist}) = learner.π(s, Val(:dist))

function update!(learner::DifferentialTDLearner{<:AbstractApproximator, <:PolicyOrSelector, :SARSA}, buffer::Union{EpisodeSARDBuffer, CircularSARDBuffer})
    n = learner.n
    update!(learner,
            @view(buffer.state[max(1, end - n - 1) : end - 1]),
            @view(buffer.action[max(1, end - n - 1) : end - 1]),
            @view(buffer.reward[max(1, end - n) : end]),
            @view(buffer.state[max(1, end - n) : end]),
            @view(buffer.action[max(1, end - n) : end]),
            Val(buffer.isdone[end]))
end

function update!(learner::DifferentialTDLearner{<:AbstractQApproximator, <:PolicyOrSelector, :SARSA}, states, actions, rewards, nextstates, nextactions, ::Val{false})
    n, α, β, Q = learner.n, learner.α, learner.β, learner.approximator
    if length(states) ≥ n
        s, a = states[1], actions[1]
        δ = sum(r -> r - learner.R̄, rewards) + Q(nextstates[end], nextactions[end]) - Q(s, a)
        learner.R̄ += β * δ
        update!(Q, s, a, α * δ)
        update!(learner.π, s, Q(s, Val(:argmax)))
    end
end


text:
# Algorithm 7.2
# Gauss-Seidel Iteration

using LinearAlgebra

function gaussseidel(mat, approx, known, maxIter, tol)
	n, m = size(mat)
	sol = Vector(undef, n)
	a = mat
	for k in 1:maxIter
		for i in 1:n
			sum1 = 0
			for j in 1:(i - 1)
				if j != i
					sum1 += (a[i, j]*x[j])
				end
			end
			sum2 = 0
			for j in (i + 1):n
				if j != i
					sum2 += (a[i, j]*approx[j]) + known[i]
				end
			end
			sol[i] = (-sum1 - sum2)/a[i, i]
		end
		if LinearAlgebra.norm(sol - approx) < tol
			println("Success")
			break
		end
		for i in 1:n
			approx[i] = sol[i]
		end
	end
	return sol
end

A = [[1.0 1.33]; [2.0 -2.55];]
x = [0.0, 0.0]
b = [2.0, 5.0]
println(gaussseidel(A, x, b, 90, 0.001))


text:
module MultistartOptimization

# NOTE: exports in included files

using ArgCheck: @argcheck
using Base.Threads: @spawn, fetch
using DocStringExtensions: FIELDS, FUNCTIONNAME, SIGNATURES, TYPEDEF
using Parameters: @unpack
using Requires: @require
using Sobol: SobolSeq, Sobol

include("generic_api.jl")
include("tiktak.jl")

function __init__()
    @require NLopt="76087f3c-5699-56af-9a33-bf431cd00edd" include("nlopt.jl")
end

end # module


text:
module JavaCall
export JavaObject, JavaMetaClass,
       jint, jlong, jbyte, jboolean, jchar, jshort, jfloat, jdouble,
       JObject, JClass, JMethod, JString,
       @jimport, jcall, jfield, isnull,
       getname, listmethods, getreturntype, getparametertypes

using Base.Dates
using WinReg
using Compat
import Compat.String

import Base.convert, Base.isnull, Base.unsafe_convert

# if VERSION < v"0.5.0-dev+4612"
# 	const unsafe_string = Base.bytestring
# else
# 	const unsafe_string = Base.unsafe_string
# end

if VERSION < v"0.5.0-dev+4612"
    import Compat.unsafe_string
else
    import Base.unsafe_string
end

include("jvm.jl")
include("jnienv.jl")
include("core.jl")
include("convert.jl")
include("reflect.jl")

function __init__()
	findjvm()
	global create = Libdl.dlsym(libjvm, :JNI_CreateJavaVM)
end

end # module


text:
struct TensorKeyIterator{G<:Sector,F₁<:FusionTree{G},F₂<:FusionTree{G}}
    rowr::SectorDict{G, FusionTreeDict{F₁, UnitRange{Int}}}
    colr::SectorDict{G, FusionTreeDict{F₂, UnitRange{Int}}}
end
struct TensorPairIterator{G<:Sector, F₁<:FusionTree{G}, F₂<:FusionTree{G}, A<:DenseMatrix}
    rowr::SectorDict{G, FusionTreeDict{F₁, UnitRange{Int}}}
    colr::SectorDict{G, FusionTreeDict{F₂, UnitRange{Int}}}
    data::SectorDict{G, A}
end

const TensorIterator{G<:Sector,F₁<:FusionTree{G},F₂<:FusionTree{G}} = Union{TensorKeyIterator{G,F₁,F₂},TensorPairIterator{G,F₁,F₂}}

Base.IteratorSize(::Type{<:TensorIterator}) = Base.HasLength()
Base.IteratorEltype(::Type{<:TensorIterator}) = Base.HasEltype()
Base.eltype(T::Type{TensorKeyIterator{G,F₁,F₂}}) where {G,F₁,F₂} = Tuple{F₁,F₂}

function Base.length(t::TensorKeyIterator)
    l = 0
    for (rowdict, coldict) in zip(values(t.rowr), values(t.colr))
        l += length(rowdict) * length(coldict)
    end
    return l
end
function Base.iterate(it::TensorKeyIterator)
    i = 1
    i > length(it.rowr) && return nothing
    rowit, colit = it.rowr.values[i], it.colr.values[i]

    rownext = iterate(rowit)
    colnext = iterate(colit)
    # while rownext === nothing || colnext === nothing: Julia did not infer that after while loop, both were not nothing
    while true
        if rownext === nothing
            i += 1
        elseif colnext === nothing
            i += 1
        else
            break
        end
        i > length(it.rowr) && return nothing
        rowit, colit = it.rowr.values[i], it.colr.values[i]
        rownext = iterate(rowit)
        colnext = iterate(colit)
    end
    (f1, r1), rowstate = rownext
    (f2, r2), colstate = colnext

    return (f1, f2), (f2, i, rowstate, colstate)
end
function Base.iterate(it::TensorKeyIterator, state)
    (f2, i, rowstate, colstate) = state
    rowit, colit = it.rowr.values[i], it.colr.values[i]
    rownext = iterate(rowit, rowstate)
    if rownext !== nothing
        (f1, r1), rowstate = rownext
        return (f1, f2), (f2, i, rowstate, colstate)
    end
    colnext = iterate(colit, colstate)
    if colnext !== nothing
        rownext = iterate(rowit) # should not be nothing
        @assert rownext !== nothing
        (f1, r1), rowstate = rownext
        (f2, r2), colstate = colnext
        return (f1, f2), (f2, i, rowstate, colstate)
    end
    while true
        if rownext === nothing
            i += 1
        elseif colnext === nothing
            i += 1
        else
            break
        end
        i > length(it.rowr) && return nothing
        rowit, colit = it.rowr.values[i], it.colr.values[i]
        rownext = iterate(rowit)
        colnext = iterate(colit)
    end
    (f1, r1), rowstate = rownext
    (f2, r2), colstate = colnext

    return (f1, f2), (f2, i, rowstate, colstate)
end

# WARNING: This only works if both SectorDict and FusionTreeDict are VectorDict
# function Base.iterate(it::TensorKeyIterator, s = (1,1,1))
#     i,j,k = s
#     length(it.rowr) < i && return nothing
#     @inbounds begin
#         f1 = it.rowr.values[i].keys[j]
#         f2 = it.colr.values[i].keys[k]
#     end
#
#     if j < length(it.rowr.values[i])
#         j += 1
#     elseif k < length(it.colr.values[i])
#         j = 1
#         k += 1
#     else
#         j = 1
#         k = 1
#         i += 1
#     end
#
#     return (f1,f2), (i,j,k)
# end


text:
function benchmark_tracking_loop()
    #init signals
    carrier_doppler = 0Hz
    start_code_phase = 0
    code_frequency = 1023kHz
    sampling_frequency = 2.5MHz
    prn = 1
    start_carrier_phase = 0
    cpustate = Tracking.TrackingState(GPSL1, carrier_doppler, start_code_phase)
    gpustate = GNSSBenchmarks.gpuTrackingState(GPSL1, carrier_doppler, start_code_phase)
    sgpustate = GNSSBenchmarks.sgpuTrackingState(GPSL1, carrier_doppler, start_code_phase)
    #init data frame
    results = DataFrame(
        Samples = SAMPLES,
        sCPU_median = zeros(Float32,length(SAMPLES)),
        GPU_median = zeros(Float32,length(SAMPLES)),
        sGPU_median = zeros(Float32,length(SAMPLES))
    )
    counter = Int32(1)
    for N in SAMPLES
        cpusignal = cis.(
            2π .* carrier_doppler .* (1:N) ./ sampling_frequency .+ start_carrier_phase
        ) .*
        GNSSSignals.get_code.(
            GPSL1,
            code_frequency .* (1:N) ./ sampling_frequency .+ start_code_phase,
            prn
        )
        println("Benchmarking the tracking loop on CPU: StructArray{ComplexF32}(Array, Array) ", N," samples...")
        result = median(@benchmark Tracking.track(
            $cpusignal,
            $cpustate,
            $prn,
            $sampling_frequency,
        )).time
        println(result)
        results.sCPU_median[counter] = result
        println("Benchmarking the tracking loop on GPU: CuArray{ComplexF32} ", N, " samples...")
        result = median(@benchmark gpu_track(
            CuArray{ComplexF32}(cpusignal),
            $gpustate,
            $prn,
            $sampling_frequency,
        )).time
        println(result)
        results.GPU_median[counter] = result
        println("Benchmarking the tracking loop GPU: StructArray{ComplexF32}(CuArray,CuArray) ", N, " samples...")
        result = median(@benchmark Tracking.track(
            StructArray{ComplexF32}(
                (real(CuArray{ComplexF32}(cpusignal)),
                imag(CuArray{ComplexF32}(cpusignal)))),
            $sgpustate,
            $prn,
            $sampling_frequency,
        )).time
        println(result)
        results.sGPU_median[counter] = result
        counter += 1
    end
    CSV.write("data/tracking_loop.csv", results)
end

function benchmark_tracking_loop_feedback(
    signal::Union{CuArray, StructArray},
    state::Union{gpuTrackingState, sgpuTrackingState},
    prn::Integer,
    sampling_frequency::typeof(1.0Hz)
)
    counter = 1
    results = gpu_track(signal, state, prn, sampling_frequency)
    while counter < 9999
        # tracking loop that feeds the state back
        results = gpu_track(signal, get_state(results), prn, sampling_frequency)
        counter += 1
    end
end


text:
struct TransDateTimeBoolVariate <: AbstractBoolVariate
    name::String
    basevariate::AbstractDateTimeVariate
    transform::Function
end

Base.length(var::TransDateTimeBoolVariate) = length(var.basevariate)

function slice(variate::TransDateTimeBoolVariate, fromobs::Integer, toobs::Integer, slicelength::Integer)
    basevar = variate.basevariate
    f = variate.transform
    g = (ms::Int64) -> begin
        ms == zero(Int64) ? false : begin
            dt = Dates.epochms2datetime(ms)
            f(dt)
        end
    end
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slices = slice(basevar, fromobs, toobs, slicelength)
    mapslice(g, slices, slicelength, Bool)
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML