text:
@with_kw struct Parameter

    T=Float32                 # number format

    Tprog=T                   # number format for prognostic variables
    Tcomm=Tprog               # number format for ghost-point copies
    Tini=Tprog                # number format to reduce precision for initial conditions

    # DOMAIN RESOLUTION AND RATIO
    nx::Int=100                         # number of grid cells in x-direction
    Lx::Real=4000e3                     # length of the domain in x-direction [m]
    L_ratio::Real=2                     # Domain aspect ratio of Lx/Ly

    # PHYSICAL CONSTANTS
    g::Real=0.1                         # gravitational acceleration [m/s]
    H::Real=500.                        # layer thickness at rest [m]
    ρ::Real=1e3                         # water density [kg/m^3]
    ϕ::Real=45.                         # central latitude of the domain (for coriolis) [°]
    ω::Real=2π/(24*3600)                # Earth's angular frequency [s^-1]
    R::Real=6.371e6                     # Earth's radius [m]

    # SCALE
    scale::Real=2^6                     # multiplicative scale for the momentum equations u,v
    scale_sst::Real=2^15                # multiplicative scale for sst

    # WIND FORCING OPTIONS
    wind_forcing_x::String="shear"      # "channel", "double_gyre", "shear","constant" or "none"
    wind_forcing_y::String="constant"   # "channel", "double_gyre", "shear","constant" or "none"
    Fx0::Real=0.12                      # wind stress strength [Pa] in x-direction
    Fy0::Real=0.0                       # wind stress strength [Pa] in y-direction
    seasonal_wind_x::Bool=true          # Change the wind stress with a sine of frequency ωFx,ωFy
    seasonal_wind_y::Bool=false         # same for y-component
    ωFx::Real=2                         # frequency [1/year] for x component
    ωFy::Real=2                         # frequency [1/year] for y component

    # BOTTOM TOPOGRAPHY OPTIONS
    topography::String="ridges"         # "ridge", "seamount", "flat", "ridges", "bathtub"
    topo_ridges_positions::Vector = [0.05,0.25,0.45,0.9]
    topo_height::Real=100.               # height of seamount [m]
    topo_width::Real=300e3              # horizontal scale [m] of the seamount

    # SURFACE RELAXATION
    surface_relax::Bool=false           # yes?
    t_relax::Real=100.                  # time scale of the relaxation [days]
    η_refh::Real=5.                     # height difference [m] of the interface relaxation profile
    η_refw::Real=50e3                   # width [m] of the tangent used for the interface relaxation

    # SURFACE FORCING
    surface_forcing::Bool=false         # yes?
    ωFη::Real=1.0                       # frequency [1/year] for surfance forcing
    A::Real=3e-5                        # Amplitude [m/s]
    ϕk::Real=ϕ                          # Central latitude of Kelvin wave pumping
    wk::Real=10e3                       # width [m] in y of Gaussian used for surface forcing

    # TIME STEPPING OPTIONS
    time_scheme::String="RK"            # Runge-Kutta ("RK") or strong-stability preserving RK
                                        # "SSPRK2","SSPRK3","4SSPRK3"
    RKo::Int=4                          # Order of the RK time stepping scheme (2, 3 or 4)
    RKs::Int=3                          # Number of stages for SSPRK2
    RKn::Int=5                          # n^2 = s = Number of stages  for SSPRK3
    cfl::Real=0.9                       # CFL number (1.0 recommended for RK4, 0.6 for RK3)
    Ndays::Real=200.0                   # number of days to integrate for
    nstep_diff::Int=1                   # diffusive part every nstep_diff time steps.
    nstep_advcor::Int=0                 # advection and coriolis update every nstep_advcor time steps.
                                        # 0 means it is included in every RK4 substep
    compensated::Bool=false             # Compensated summation in the time integration?

    # BOUNDARY CONDITION OPTIONS
    bc::String="periodic"               # "periodic" or anything else for nonperiodic
    α::Real=2.                          # lateral boundary condition parameter
                                        # 0 free-slip, 0<α<2 partial-slip, 2 no-slip

    # MOMENTUM ADVECTION OPTIONS
    adv_scheme::String="ArakawaHsu"     # "Sadourny" or "ArakawaHsu"
    dynamics::String="nonlinear"        # "linear" or "nonlinear"

    # BOTTOM FRICTION OPTIONS
    bottom_drag::String="none"          # "linear", "quadratic" or "none"
    cD::Real=1e-5                       # bottom drag coefficient [dimensionless] for quadratic
    τD::Real=300.                       # bottom drag coefficient [days] for linear

    # DIFFUSION OPTIONS
    diffusion::String="constant"        # "Smagorinsky" or "constant", biharmonic in both cases
    νB::Real=500.0                      # [m^2/s] scaling constant for constant biharmonic diffusion
    cSmag::Real=0.15                    # Smagorinsky coefficient [dimensionless]

    # TRACER ADVECTION
    tracer_advection::Bool=true         # yes?
    tracer_relaxation::Bool=true        # yes?
    tracer_consumption::Bool=false      # yes?
    sst_initial::String="waves"         # "west", "south", "linear", "waves","rect", "flat" or "restart"
    sst_rect_coords::Array{Float64,1}=[0.,0.15,0.,1.0]
                                        # (x0,x1,y0,y1) are the size of the rectangle in [0,1]
    Uadv::Real=0.2                      # Velocity scale [m/s] for tracer advection
    SSTmax::Real=1.                     # tracer (sea surface temperature) max for initial conditions
    SSTmin::Real=-1.                    # tracer (sea surface temperature) min for initial conditions
    τSST::Real=100                      # tracer restoring time scale [days]
    jSST::Real=365                      # tracer consumption [days]
    SSTw::Real=5e5                      # width [m] of the tangent used for the IC and interface relaxation
    SSTϕ::Real=0.5                      # latitude/longitude fraction ∈ [0,1] of sst edge
    SSTwaves_ny::Real=4                 # wave crests/troughs in y
    SSTwaves_nx::Real=SSTwaves_ny*L_ratio  # wave crests/troughs in x
    SSTwaves_p::Real=1/2                # power for rectangles (p<1)/smootheness(p>=1) of waves

    # OUTPUT OPTIONS
    output::Bool=false                  # netcdf output?
    output_vars::Array{String,1}=["u","v","η","sst"]  # which variables to output? "du","dv","dη" also allowed
    output_dt::Real=24                  # output time step [hours]
    outpath::String=pwd()               # path to output folder
    compression_level::Int=3            # compression level
    return_time::Bool=false             # return time of simulation of progn vars?

    # INITIAL CONDITIONS
    initial_cond::String="rest"         # "rest" or "ncfile" for restart from file
    initpath::String=outpath            # folder where to pick the restart files from
    init_run_id::Int=0                  # run id for restart from run number
    init_starti::Int=-1                 # timestep to start from (-1 meaning last)
    get_id_mode::String="continue"      # How to determine the run id: "continue" or "fill"
    run_id::Int=-1                      # Output with a specific run id
    init_interpolation::Bool=true       # Interpolate the initial conditions in case grids don't match?

    # ASSERT - CHECK THAT THE INPUT PARAMETERS MAKE SENSE
    @assert all((nx,Lx,L_ratio) .> 0.)  "nx, Lx, L_ratio have to be >0"
    @assert all((g,H,ρ,ω,R) .> 0.)      "g,H,ρ,ω,R have to be >0"
    @assert ϕ <= 90.0 && ϕ >= -90.0     "ϕ has to be in (-90,90), $ϕ given."
    @assert wind_forcing_x in ["channel","double_gyre","shear","constant","none"] "Wind forcing '$wind_forcing_x' unsupported"
    @assert wind_forcing_y in ["channel","double_gyre","shear","constant","none"] "Wind forcing '$wind_forcing_y' unsupported"
    @assert topography in ["ridge","seamount","flat","ridges"] "Topography '$topography' unsupported"
    @assert topo_width > 0.0    "topo_width has to be >0, $topo_width given."
    @assert t_relax > 0.0       "t_relax has to be >0, $t_relax given."
    @assert η_refw > 0.0        "η_refw has to be >0, $η_refw given."
    @assert time_scheme in ["RK","SSPRK2","SSPRK3","4SSPRK3"] "Time scheme $time_scheme unsupported."
    @assert RKo in [2,3,4]        "RKo has to be 2,3 or 4; $RKo given."
    @assert RKs > 1               "RKs has to be >= 2; $RKs given."
    @assert Ndays > 0.0         "Ndays has to be >0, $Ndays given."
    @assert nstep_diff > 0      "nstep_diff has to be >0, $nstep_diff given."
    @assert nstep_advcor >= 0   "nstep_advcor has to be >=0, $nstep_advcor given."
    @assert bc in ["periodic","nonperiodic"]    "boundary condition '$bc' unsupported."
    @assert α >= 0.0 && α <= 2.0    "Tangential boundary condition α has to be in (0,2), $α given."
    @assert adv_scheme in ["Sadourny","ArakawaHsu"] "Advection scheme '$adv_scheme' unsupported"
    @assert dynamics in ["linear","nonlinear"]  "Dynamics '$dynamics' unsupported."
    @assert bottom_drag in ["quadratic","linear","none"] "Bottom drag '$bottom_drag' unsupported."
    @assert cD >= 0.0    "Bottom drag coefficient cD has to be >=0, $cD given."
    @assert τD >= 0.0    "Bottom drag coefficient τD has to be >=0, $τD given."
    @assert diffusion in ["Smagorinsky", "constant"] "Diffusion '$diffusion' unsupported."
    @assert νB > 0.0     "Diffusion scaling constant νB has to be >0, $νB given."
    @assert cSmag > 0.0  "Smagorinsky coefficient cSmag has to be >0, $cSmag given."
    @assert Uadv > 0.0   "Advection velocity scale Uadv has to be >0, $Uadv given."
    @assert output_dt > 0   "Output time step has to be >0, $output_dt given."
    @assert initial_cond in ["rest", "ncfile"] "Initial conditions '$initial_cond' unsupported."
    @assert init_run_id >= 0 "Initial condition run id, init_run_id, has to be >= 0, $init_run_id given."
    @assert init_starti > 0 || init_starti == -1 "Start index, init_starti, has to be >0 || -1, $init_starti given."
    @assert get_id_mode in ["continue","fill","specific"] "get_id_mode $get_id_mode unsupported."
end

"""
Creates a Parameter struct with following options and default values

    T=Float32                 # number format

    Tprog=T                   # number format for prognostic variables
    Tcomm=Tprog               # number format for ghost-point copies
    Tini=Tprog                # number format to reduce precision for initial conditions

    # DOMAIN RESOLUTION AND RATIO
    nx::Int=100                         # number of grid cells in x-direction
    Lx::Real=4000e3                     # length of the domain in x-direction [m]
    L_ratio::Real=2                     # Domain aspect ratio of Lx/Ly

    # PHYSICAL CONSTANTS
    g::Real=0.1                         # gravitational acceleration [m/s]
    H::Real=500.                        # layer thickness at rest [m]
    ρ::Real=1e3                         # water density [kg/m^3]
    ϕ::Real=45.                         # central latitude of the domain (for coriolis) [°]
    ω::Real=2π/(24*3600)                # Earth's angular frequency [s^-1]
    R::Real=6.371e6                     # Earth's radius [m]

    # SCALE
    scale::Real=2^6                     # multiplicative scale for the momentum equations u,v
    scale_sst::Real=2^15                # multiplicative scale for sst

    # WIND FORCING OPTIONS
    wind_forcing_x::String="shear"      # "channel", "double_gyre", "shear","constant" or "none"
    wind_forcing_y::String="constant"   # "channel", "double_gyre", "shear","constant" or "none"
    Fx0::Real=0.12                      # wind stress strength [Pa] in x-direction
    Fy0::Real=0.0                       # wind stress strength [Pa] in y-direction
    seasonal_wind_x::Bool=true          # Change the wind stress with a sine of frequency ωFx,ωFy
    seasonal_wind_y::Bool=false         # same for y-component
    ωFx::Real=2                         # frequency [1/year] for x component
    ωFy::Real=2                         # frequency [1/year] for y component

    # BOTTOM TOPOGRAPHY OPTIONS
    topography::String="ridges"         # "ridge", "seamount", "flat", "ridges", "bathtub"
    topo_ridges_positions::Vector = [0.05,0.25,0.45,0.9]
    topo_height::Real=100.               # height of seamount [m]
    topo_width::Real=300e3              # horizontal scale [m] of the seamount

    # SURFACE RELAXATION
    surface_relax::Bool=false           # yes?
    t_relax::Real=100.                  # time scale of the relaxation [days]
    η_refh::Real=5.                     # height difference [m] of the interface relaxation profile
    η_refw::Real=50e3                   # width [m] of the tangent used for the interface relaxation

    # SURFACE FORCING
    surface_forcing::Bool=false         # yes?
    ωFη::Real=1.0                       # frequency [1/year] for surfance forcing
    A::Real=3e-5                        # Amplitude [m/s]
    ϕk::Real=ϕ                          # Central latitude of Kelvin wave pumping
    wk::Real=10e3                       # width [m] in y of Gaussian used for surface forcing

    # TIME STEPPING OPTIONS
    time_scheme::String="SSPRK3"        # Runge-Kutta ("RK") or strong-stability preserving RK
                                        # "SSPRK2","SSPRK3","4SSPRK3"
    RKo::Int=4                          # Order of the RK time stepping scheme (2, 3 or 4)
    RKs::Int=3                          # Number of stages for SSPRK2
    RKn::Int=5                          # n^2 = s = Number of stages  for SSPRK3
    cfl::Real=4.0                       # CFL number (1.0 recommended for RK4, 0.6 for RK3)
    Ndays::Real=200.0                   # number of days to integrate for
    nstep_diff::Int=1                   # diffusive part every nstep_diff time steps.
    nstep_advcor::Int=0                 # advection and coriolis update every nstep_advcor time steps.
                                        # 0 means it is included in every RK4 substep

    # BOUNDARY CONDITION OPTIONS
    bc::String="periodic"               # "periodic" or anything else for nonperiodic
    α::Real=2.                          # lateral boundary condition parameter
                                        # 0 free-slip, 0<α<2 partial-slip, 2 no-slip

    # MOMENTUM ADVECTION OPTIONS
    adv_scheme::String="ArakawaHsu"     # "Sadourny" or "ArakawaHsu"
    dynamics::String="nonlinear"        # "linear" or "nonlinear"

    # BOTTOM FRICTION OPTIONS
    bottom_drag::String="none"          # "linear", "quadratic" or "none"
    cD::Real=1e-5                       # bottom drag coefficient [dimensionless] for quadratic
    τD::Real=300.                       # bottom drag coefficient [days] for linear

    # DIFFUSION OPTIONS
    diffusion::String="constant"        # "Smagorinsky" or "constant", biharmonic in both cases
    νB::Real=500.0                      # [m^2/s] scaling constant for constant biharmonic diffusion
    cSmag::Real=0.15                    # Smagorinsky coefficient [dimensionless]

    # TRACER ADVECTION
    tracer_advection::Bool=true         # yes?
    tracer_relaxation::Bool=true        # yes?
    tracer_consumption::Bool=false      # yes?
    sst_initial::String="waves"         # "west", "south", "linear", "waves","rect", "flat" or "restart"
    sst_rect_coords::Array{Float64,1}=[0.,0.15,0.,1.0]
                                        # (x0,x1,y0,y1) are the size of the rectangle in [0,1]
    Uadv::Real=0.2                      # Velocity scale [m/s] for tracer advection
    SSTmax::Real=1.                     # tracer (sea surface temperature) max for initial conditions
    SSTmin::Real=-1.                    # tracer (sea surface temperature) min for initial conditions
    τSST::Real=100                      # tracer restoring time scale [days]
    jSST::Real=365                      # tracer consumption [days]
    SSTw::Real=5e5                      # width [m] of the tangent used for the IC and interface relaxation
    SSTϕ::Real=0.5                      # latitude/longitude fraction ∈ [0,1] of sst edge
    SSTwaves_ny::Real=4                 # wave crests/troughs in y
    SSTwaves_nx::Real=SSTwaves_ny*L_ratio  # wave crests/troughs in x
    SSTwaves_p::Real=1/2                # power for rectangles (p<1)/smootheness(p>=1) of waves

    # OUTPUT OPTIONS
    output::Bool=false                  # netcdf output?
    output_vars::Array{String,1}=["u","v","η","sst"]  # which variables to output? "du","dv","dη" also allowed
    output_dt::Real=24                  # output time step [hours]
    outpath::String=pwd()               # path to output folder
    compression_level::Int=3            # compression level

    # INITIAL CONDITIONS
    initial_cond::String="rest"         # "rest" or "ncfile" for restart from file
    initpath::String=outpath            # folder where to pick the restart files from
    init_run_id::Int=0                  # run id for restart from run number
    init_starti::Int=-1                 # timestep to start from (-1 meaning last)
    get_id_mode::String="continue"      # How to determine the run id: "continue" or "fill"
    run_id::Int=-1                      # Output with a specific run id
    init_interpolation::Bool=true       # Interpolate the initial conditions in case grids don't match?
"""
Parameter


text:
using Xpress, Test

Xpress.initialize()

@testset "Xpress API" begin

    @test Xpress.getversion() isa VersionNumber

    @test Xpress.getbanner() isa String

    xp = Xpress.XpressProblem()

    @test Xpress.getprobname(xp) == ""
    @test Xpress.setprobname(xp, "xpress-optimization-problem") == nothing
    @test Xpress.getprobname(xp) == "xpress-optimization-problem"

    @test Xpress.getintcontrol(xp, Xpress.Lib.XPRS_DEFAULTALG) == 1
    @test Xpress.getcontrol(xp, Xpress.Lib.XPRS_DEFAULTALG) == 1

    @test Xpress.setcontrol!(xp, Xpress.Lib.XPRS_PRESOLVE, 0) == nothing
    @test Xpress.getcontrol(xp, :XPRS_PRESOLVE) == 0
    @test Xpress.setcontrol!(xp, Xpress.Lib.XPRS_PRESOLVE, 1) == nothing
    @test Xpress.getcontrol(xp, Xpress.Lib.XPRS_PRESOLVE) == 1

    iob = IOBuffer()
    show(iob, xp)
    @test String(take!(iob)) == """
Xpress Problem:
    type   : LP
    sense  : minimize
    number of variables                    = 0
    number of linear constraints           = 0
    number of quadratic constraints        = 0
    number of sos constraints              = 0
    number of non-zero coeffs              = 0
    number of non-zero qp objective terms  = 0
    number of non-zero qp constraint terms = 0
    number of integer entities             = 0
"""

end

tests = ["MOI_Wrapper",
        #  "wordhunt"
         ]

for t in tests
    fp = "$(t).jl"
    println("running $(fp) ...")
    evalfile(fp)
end

ext: jl


text:
### A Pluto.jl notebook ###
# v0.10.0

using Markdown

# ╔═╡ b7a4d6c3-dc52-45c0-bcd6-dc1d32dbd302
md"""
# Factorizaciones y otras diversiones
Autor: Andreas Noack Jensen (MIT) (http://www.econ.ku.dk/phdstudent/noack/)
(con edición de Jane Herriman=

## Sinopsis
 - Factorizaciones
 - Estructuras de matrices especiales
 - Álgebra lineal genérica
"""

# ╔═╡ 43882091-fc1b-409a-92d9-455ca7ad6f64
md"""
Empezamos con un sistema lineal de la forma

`Ax = b`
"""

# ╔═╡ 85ad00b8-4b91-46c9-a130-4430f2986569
A = randn(3,3)

# ╔═╡ 5a8dd4eb-40a6-4750-b2e7-d44b30404363
x = fill(1.0, (3))
b = A*x

# ╔═╡ a98929e4-ce11-49b1-aeca-44ed3680df53
md"""
### Factorización
La función `\` esconde cómo regularmente se resuelve el problema.

Dependiendo de las dimensiones de `A`, distintos métodos son elegidos para resolver el problema.

```
Ax = b
```

Un paso intermedio en la solución es el cálculo de la factorización de la matriz `A`.

Básicamente, una factorización de `A` es una manera de expresar `A` como el producto de matrices triangulares, unitarias, y de permutación.

Julia guarda estas factorizaciones usando un tipo abstracto llamado `Factorization` y varios subtipos.

Un objeto `Factorization` entonces debería ser pensado como una represenatación de la matriz `A`.

#### LU
"""

# ╔═╡ e00b9a3c-3ffc-4aac-b5f7-a88a0daafba4
md"""
Cuando `A` es cuadrada, el sistema linear es resuelto factorizando `A` vía

```
PA=LU
```

donde `P` es una matriz de permutación, `L` es un triangular inferior unitaria y `U` es triangular superior.

Julia permita calcular la facorización LU y define un tipo de factorización compuesta para guardarlo.
"""

# ╔═╡ a6bc3137-f8f6-44ec-bc7e-baa5a2f9c499
md"""
Podemos hacer una factorización LU sobre `A` vía `lu(A)` ó `lufact(A)`.

La función `lu` regresa matrices `l` y `u` y un vector de permutación `p`.
"""

# ╔═╡ 0b6b51ee-ff9f-463e-95e1-2fd1f559bdf1
l,u,p = lu(A)

# ╔═╡ 2aefd8c5-abc3-4d1a-83fb-78f22190e0f3
md"""
El pivoteo está prendido por default, o sea que no podemos assumir que A == LU.
Vamos a comprobar esto viendo la norma de `LU - A`:
"""

# ╔═╡ 94c7fee6-2ada-4223-8fc0-57c7896181aa
norm(l*u - A)

# ╔═╡ edc134d0-de77-429a-8bac-36ef19e80471
md"""
Esto nos muestra que queremos tomar en cuenta el pivoteo!

Podemos pensar en `A[p,:]` como la sintaxis para `PA`, o como el producto de una matriz de permutación con `A`
"""

# ╔═╡ a123c969-8e5c-432b-9bbd-3f1dfa4b8f31
norm(l*u - A[p,:])

# ╔═╡ da792d5e-9ebf-4b93-9d57-61745021ec00
md"""
Por otro lado, podemos apagar el pivoteo usando `Val{false}` en Julia 0.6 ó `Val(false)` en versiones más modernas.
"""

# ╔═╡ 46ab0bb7-ded7-4d98-a452-70242a5bf648
l,u,p = lu(A, Val{false})

# ╔═╡ 8156b8f8-9919-42c7-83cd-49589bfd53e7
md"""
Cuando apagamos el pivoteo, `LU = A`
"""

# ╔═╡ 608bfb30-6530-474b-9374-c5caa3a716f8
norm(l*u - A)

# ╔═╡ cfd8ea62-ff51-4966-ad3e-2d12a937e71e
md"""
Una segunda manera de hacer la factorización es con `lufact`.
"""

# ╔═╡ 289f4361-2962-48f4-8c47-915290a26035
Alu = lufact(A)

# ╔═╡ 039fdbaf-9870-4efc-a54a-f457088c5e88
md"""
Distintas partes de la factorización las puedes accesar con índices especiales
"""

# ╔═╡ d8217c6a-2f5f-4824-99e9-e56bfca07e92
Alu[:P]

# ╔═╡ 7d46f4fb-e092-409f-a0b0-5caad91f4022
Alu[:L]

# ╔═╡ 44aec045-5803-41d2-9574-31a9f14d9cf7
Alu[:U]

# ╔═╡ 466124e7-9978-4c7c-9cf1-6e34c7f55c33
md"""
Podemos calcular la solución de $Ax=b$ del objecto de factorización
"""

# ╔═╡ 9a023eb2-ce3d-413e-a872-231325c556b6
# PA = LU
# A = P'LU
# P'LUx = b
# LUx = Pb
# Ux = L\Pb
# x = U\L\Pb
Alu[:U]\(Alu[:L]\(Alu[:P]b))

# ╔═╡ cf9efbb5-87ad-4795-b132-8877b334f85d
md"""
*Más importantemente,* podemos despachar sobre el tipo de `LU` y simplemente resolver el problema por medio de
"""

# ╔═╡ 45f9ca18-2ba7-4e8c-9087-ebadf2aa7101
Alu\b

# ╔═╡ 424a4d00-0264-45cf-bbba-43fb51ffa394
md"""
Esto puede ser útil si el mismo lado izquierdo es usado para lados derechos.

La factorización también se puede usar para calcular el determinante pues

$\det(A)=\det(PLU)=\det(P)\det(U)=\pm \prod u_{ii}$

porque $U$ es triangular y el signo está definido por $\det(P)$.
"""

# ╔═╡ ba2407dd-ef47-4532-a783-aafb17506e5f
det(A)

# ╔═╡ 299bec2c-626f-40ca-a2f9-4e0450b23f3e
det(Alu)

# ╔═╡ 0fa226f7-8a15-466b-9757-3c019fd4c58e
md"""
#### QR
Cuando `A` es alta,
"""

# ╔═╡ 92962e9b-13c2-465e-ae4d-b1c8a0172063
Atall = randn(3, 2)

# ╔═╡ e802da18-5118-4031-9be8-0fb90a7a0c9d
md"""
Julia calcula la solución de mínimos cuadrados $\hat{x}$ que minimiza $\|Ax-b\|_2$.

Esto se puede hacer factorizando

```
A=QR
```

donde $Q$ es unitaria/ortogonal y

$R=\left(\begin{smallmatrix}R_0\\0\end{smallmatrix}\right)$ y  $R_0$ es triangular superior.

con la factorización QR la norma mínima se puede expresar como

\begin{equation*}
\|Ax-b\|=\|QRx-b\|=\|Q(Rx-Q'b)\|=\|Rx-Q'b\|=\left\|\begin{pmatrix}R_0x-Q_0'b\\Q_1'b\end{pmatrix}\right\|=\|R_0x-Q_0'b\|+\|Q_1'b\|
\end{equation*}

Y entonces el problema se puede reducir a resolver el problema cuadrado $R_0x=Q_0'b$ para $x$.

Podemos factorizar QR sobre `Atall` vía
"""

# ╔═╡ 215ea9e5-1a3a-4750-b249-2ed54f3ea9ac
Aqr = qrfact(Atall)

# ╔═╡ 021e9680-4a33-469e-8014-ae3b831db370
md"""
Otra característica de la factorización QR es que los tipos `Q` para guardar las matrices unitarias $Q$. Se pueden extraer de tipos `QR` con los índices
"""

# ╔═╡ e55a8449-a507-4e90-9f1c-013edb26704d
Aqr[:Q]

# ╔═╡ 60ce24a1-2d71-4efc-95a8-3e210da5c53f
md"""
Similarmente, la matriz superior triangular $R$ se puede extraer con el índice
"""

# ╔═╡ fa8dca75-5aca-454d-b733-8da5b10c6857
Aqr[:R]

# ╔═╡ 60918647-3d3d-4dec-a4e8-e0cf4348d072
md"""
En este caso se guarda R como una matriz 2x2 en vez de 3x2 porque el último renglón de R está lleno de 0's.

Aunque la matriz `Aqr[:Q]` se imprime como $3\times 3$ en el objeto de factorización, en la práctica puede representar la versión delgada también. Así
"""

# ╔═╡ 83e8b0eb-f8c9-4e4d-b319-27134985ffd1
Aqr[:Q]*ones(2)

# ╔═╡ 9aa2af05-bbb0-4570-8958-0ac08c5da7b5
md"""
funciona y representa a $3 x 2$ matrix por un vector de 2-elementos.

Similarmente,
"""

# ╔═╡ 4ce42698-361a-4929-b03e-23b43ad04d33
Aqr[:Q]*ones(3)

# ╔═╡ 32879b97-8789-408f-a992-5f21dcd6e8fe
md"""
funciona representando la matriz $3x3$ y un vector de 3 elementos.

Sin embargo, esto no significa que podemos multiplicar a `Q` por vectores de longitued arbitraria.
"""

# ╔═╡ 4c3f2f14-3bd3-45c1-b36d-1000811fc622
Aqr[:Q]*ones(4)

# ╔═╡ 6a2732c8-241f-426c-8b51-ba130c3a4047
md"""
La matriz tiene representación interna compacta, entonces indexar sólo hace sentido si uno sabe cómo la factorización guarda datos.
"""

# ╔═╡ 45af09e4-8d42-404e-8e69-e0d185f4ea27
Aqr[:Q][1]

# ╔═╡ 8227c126-4417-4d6a-9c40-232ae467b17d
md"""
El objeto QRCompactWY `\` tiene un método para QR y entonces el problema de los mínimos cuadrados es resuelto con
"""

# ╔═╡ cad5d92d-d7cb-444d-9664-e9eccbbb5b06
Aqr\b

# ╔═╡ 6d694031-b7e4-423e-a341-e1e78e26d484
md"""
Y si en vez escribimos
"""

# ╔═╡ ac57e6e1-986c-4efe-ad64-e2bb3a4aeb4c
Atall\b

# ╔═╡ aa4e579d-5b59-494c-be51-dc65e6b57aa6
md"""
En vez de factorizar con QR a `Atall` primero, Julia va a defaultear factorizar *con* pivoteo.
"""

# ╔═╡ 61434524-e698-4a8d-a430-0126ac99157e
md"""
Este default a pivotear la factorización QR le permite a Julia resolver problemas deficientes de rango.

Podemos explícitamente escoger pivotear durante la factorización QR (de una matriz singular, por ejemplo) con `Val{true}`.
"""

# ╔═╡ 04f4492e-5156-4526-b761-ab3c5a209646
v = randn(3)
# Tomar el producto exterio de un vector consigo mismo nos da una matriz singular
singmatriz = v * v'

# ╔═╡ f8bb7598-095d-449a-86e2-f8bd6b9d8b94
Aqrp = qrfact(singmatriz,Val{true})

# ╔═╡ a8e4b120-54b2-4e16-b993-5a884c29309d
md"""
Notamos que el tipo que resulta del objecto de Factorization es distinto que antes.

`\` también tiene un método de `QRPivoted` y el problema con rango deficiente es entonces calculado como
"""

# ╔═╡ c85c1276-a8d9-4f02-9fd7-0a448da34436
Aqrp\b

# ╔═╡ 9ea6c93a-e36f-4624-9598-7e4741efa50a
md"""
#### Eigendescompisición y los SVDs (Valores de descomposición Singular)
"""

# ╔═╡ 6a578f8f-29fa-4583-a328-e8f33ba80738
md"""
Los resultados de eigendescomposición y de la descomposición singular de valores se guardan en los tipos`Factorization`. Esto también incluye la factorización de Hessenberg y de Schur

La eigendescomposición puede ser calculada como
"""

# ╔═╡ 3e8c3380-7db5-4e35-b6f8-9b94dc417f80
Asym = A + A'
AsymEig = eigfact(Asym)

# ╔═╡ 0c0caea9-3ee9-4e30-8570-4d30ed7cc14a
md"""
Los valores y vectores se pueden recoger del tipo Eigen con un índice especial
"""

# ╔═╡ 3cde0ff8-2325-47b7-a1ff-d15155c2f634
AsymEig[:values]

# ╔═╡ c4d93f45-2191-4624-ade9-277b546ced95
AsymEig[:vectors]

# ╔═╡ b51b8521-6f77-410e-ae49-de138a1aec67
md"""
Una vez más, como la descomposición se guarda en un tipo, podemos despatchar sobre esos tipos y explotar un método especializado para cada factorización, e.g. que $A^{-1}=(V\Lambda V^{-1})^{-1}=V\Lambda^{-1}V^{-1}$.
"""

# ╔═╡ 50f4c300-e49b-4331-b364-582cc7692dd7
inv(AsymEig)*Asym

# ╔═╡ 65386dae-fb08-47d7-a81c-c3133a8403cf
md"""
Julia también tiene una función `eig` que regresa una tupla con los valores y vectores
"""

# ╔═╡ 707539e7-c4b9-46c1-89f1-8d61a19413fa
eig(Asym)

# ╔═╡ 827e2055-5f83-44b1-a1dc-636b807c8b81
md"""
No recomendamos esta versión.

La función `svdfact` calcula la descomposición singular de valores
"""

# ╔═╡ 1afd9e85-1f7d-48aa-baca-8cdf754e4d4b
Asvd = svdfact(A[:,1:2])

# ╔═╡ 78a13f6a-ce68-4060-a1bf-1993ad5e93ad
md"""
y de nuevo `\` tiene un método para el tipo que permite los mínimos cuadrados por SVD
"""

# ╔═╡ a5b49236-f31e-4449-941e-70e1677ae1b2
Asvd\b

# ╔═╡ 9ad38b16-f9b9-4cde-bbfa-204d96ee3a43
md"""
Existen funciones especiales para proporcionar los valores sólamente: `eigvals` and `svdvals`.
"""

# ╔═╡ 0139d8bd-7a60-4cdc-8844-91c0e0eaa9c6
md"""
### Estructuras especiales de matrices

La estructura de matrices es muy importante en el álgebra linear. Ésta estructura se le puede hacer explícita a Julia por medio de los tipos compuestos. Ejemplos: `Diagonal`, `Triangular`, `Symmetric`, `Hermitian`, `Tridiagonal` y `SymTridiagonal`. Se han escrito métodos especializados para cada tipo de matriz especial para aprovechar su estructura. Siguen algunos ejemplos:
"""

# ╔═╡ 7ee1aa1c-7219-46d3-8abd-5a226c0f4abb
A

# ╔═╡ f572bd5d-f1f4-4dce-a5be-293b4404bf71
md"""
Creando una matriz Diagonal
"""

# ╔═╡ 6846c3c4-6f80-43dc-9ece-07efe11133f8
Diagonal(diag(A))

# ╔═╡ 2b5bb260-8997-4b54-ba40-c19c0bd592a5
Diagonal(A)

# ╔═╡ 7f270025-80df-4b91-babc-db876885e479
md"""
Creando una matriz triangular inferior
"""

# ╔═╡ 876fb13e-d540-4931-ad32-f9491fd504c7
LowerTriangular(tril(A))

# ╔═╡ e46ff282-c6a5-4917-8322-9ff2d1df8318
LowerTriangular(A)

# ╔═╡ b1093e1f-65fa-40d1-8b56-a374f15b3f51
md"""
Creando una matriz simétrica
"""

# ╔═╡ 9d354a03-6229-46f1-80eb-cbed5d304d5b
Symmetric(Asym)

# ╔═╡ 6e20dd10-3e4a-411e-bfdf-1534b509e8ce
SymTridiagonal(diag(Asym),diag(Asym,1))

# ╔═╡ fb94650f-06f4-4b65-a598-45c820c27c95
md"""
Cuando se sabe que una matriz es e.g. triangular o simétrica Julia puede que resuelva el problema más rápido convirtiendo a una matriz especial.

Para algunos procedimientos, Julia checa si el input de matriz es triangular o simétrica y lo convierte a tal estructura si lo detecta.

Debería notarse que `Symmetric`, `Hermitian` y `Triangular` no copian la matriz original.

#### Eigenproblema simétrico
La capacidad de Julia para poder detectar si una matriz es simétrica/Hermitian puede influenciar muchísimo sobre qué tan rápido se puede resolver un problema de eigenvalor.
"""

# ╔═╡ 7773e25a-57e8-4917-9629-7329d7ddc60c
n = 1000;
A = randn(n,n);

# ╔═╡ 2f6332ef-caef-423b-bc2d-4299436390f8
md"""
Usamos `A` para genera una matriz simétrica `Asym`
"""

# ╔═╡ b844544e-c513-42af-9b27-b006d50973c0
Asym = A + A';

# ╔═╡ f34e931a-c11b-4776-8d50-90812ef67d5b
md"""
Ahora creemos una matriz Asym para simular una matriz simétrica con errores de punto flotante
"""

# ╔═╡ 09b262f1-6a2d-4696-a194-54b65f129eb3
Asym_noisy = copy(Asym); Asym_noisy[1,2] += 5eps();

# ╔═╡ 44b64772-3daa-4037-9b6f-07d5feba3223
md"""
¿Puede Julia determinar que ambas `Asym` y `Asym_noisy` son matrices simétricas?
"""

# ╔═╡ dcef7621-ef7f-467e-9dec-75fae2bc44f3
println("Is Asym symmetric? ", issymmetric(Asym))
println("Is Asym_noisy symmetric? ", issymmetric(Asym_noisy))

# ╔═╡ e451672e-e44d-45c7-a3c2-3badfcc83d40
md"""
Ahora veamos como el ruido de `Asym_noisy` impacta el tiempo en llevar a cabo una eigendescomposición
"""

# ╔═╡ a994a69a-3879-49ad-a03a-c05d901b4072
@time eigvals(Asym);

# ╔═╡ 63d5a411-bd56-4ff6-bf35-a23995d4c2bb
@time eigvals(Asym_noisy);

# ╔═╡ 9902503c-0448-4bfc-bf6f-7c11321cca78
md"""
Por suerte, le podeemos proveer información explícita sobre la estructura de la matriz a Julia
En este ejemplo, usamos la palabra clave `Symmetric`
"""

# ╔═╡ 9d099581-b364-4ad0-990a-d3151dcdcf96
@time eigvals(Symmetric(Asym_noisy));

# ╔═╡ 8ea1527f-a664-4520-9801-840e77bd6e62
md"""
Y así nuestros cálculos son mucho más eficientes :)
"""

# ╔═╡ 231d2771-137d-4930-80a6-4ca352e88aa5
md"""
### Un gran problema
Usar matrices tridiagonales permite trabajar con problemas potencialmente muy grandes. El siguiente problema no seria posible resolverlo en una laptop si la matriz se tuviera que guardar como tipo `Matrix`.
"""

# ╔═╡ a0afdba1-25a0-409c-b84d-c63e53dbf343
n = 1_000_000;
A = SymTridiagonal(randn(n), randn(n-1));
@time eigmax(A)

# ╔═╡ 49cfb27c-8c00-4458-a77a-38dd5a4c643a
md"""
### Álgebra lineal numérica
La manera usual de agregar soporte para álgebra lineal numérica es haciendo un wrapper para subrutinas de BLAS y LAPACK. Para matrices con elementos `Float32`, `Float64`, `Complex{Float32}` ó `Complex{Float64}` esto es lo que hace Julia. Desde hace rato, Julia ha tnido soport para la multiplicación genérica de tipos. Así, cuando uno multiplica matrices de enteros, obtiene una matriz de enteros.
"""

# ╔═╡ 0f4aec6b-6701-46f9-8a60-1c91f19cbb37
rand(1:10,3,3)*rand(1:10,3,3)

# ╔═╡ 1b26ea8e-d825-48ab-a0a0-a401dd83d42c
md"""
Recientemente, más métodos de álgebra lineal se han añadido a Julia y ahora soporta factorizaciones generales de tipo `LU` y `QR`. Métodos generales para eigenvalores y SVD han sido escritos más recientemente en paquetería externa.

En general, la factorización `LU` puede ser calculada cuando los elementos de la matriz se cierraon sobre los operadores `+`, `-`, `*` y `\`. Por supuesto, la matriz también deben tener rango completo. El método general de factorización `LU` en Julia aplica pivoteo y por lo tanto debe poder soportar `<` y `abs`. Por lo tanto es posible resolver sistemas de ecuaciones de e.g. números racionales como en los ejemplos que siguen.
"""

# ╔═╡ 8b718fad-0ded-42d2-915c-7868a8fc15df
md"""
Para usar números racionales, usamos un doble //:
"""

# ╔═╡ 9858cb5f-ec34-4ce3-aa44-a676bf2d9d35
1//2

# ╔═╡ 4811bad4-fa07-4a22-92e6-c9947bce13d1
md"""
#### Ejemplo 1: Sistemas racionales lineales de ecuaciones
Julia cuenta con números racionales ya instalados. El siguiente ejemplo consta de un sistema lineal de ecuaciones resulto sin promover a tipos de punto flotanto. Puede haber un error de overflow fácilmente al trabajar con racionales, así que usamos `BigInt`s.
"""

# ╔═╡ a756d1d8-57ce-4e17-984f-7fbe1a8117e5
Ar = convert(Matrix{Rational{BigInt}}, rand(1:10,3,3))/10

# ╔═╡ e97ccbb7-fb81-4514-a8b6-0b177ef9620e
x = fill(1, (3))
b = Ar*x

# ╔═╡ a808d3f1-21a0-4d72-8323-9ffdb55084b6
Ar\b

# ╔═╡ d9e7e3c7-cf57-4b86-b2f6-7b24b3ea8bdb
lufact(Ar)

# ╔═╡ 6d5e0b76-5360-4e77-bd5a-93f5aac6011a
md"""
#### Ejemplo 2: Matriz racional de eigenestructura

El siguiente ejemplo muestra como la artimética de matriz racional puede ser usada para calcular una matriz dados los eigenvalores y eigenvectores racionales. Yo he encontrado ésto útil para mostrar ejemplos de sistemas dinámicos lineales.
"""

# ╔═╡ 57d7e0de-a3c9-4fda-b4b1-4e847bd5f8d2
λ1,λ2,λ3 = 1//1,1//2,1//4
v1,v2,v3 = [1,0,0],[1,1,0],[1,1,1]
V,Λ = [v1 v2 v3], Diagonal([λ1,λ2,λ3])
A = V*Λ/V

# ╔═╡ b29effd2-bb03-4f65-bd5c-ba5206378c54
md"""
### Exercises

11.1 ¿Cuáles son los eigenvalores de la Matriz A

```
A =
[
 140   97   74  168  131
  97  106   89  131   36
  74   89  152  144   71
 168  131  144   54  142
 131   36   71  142   36
]
```
"""

# ╔═╡ 9e39d0ed-49b7-48c6-894e-6140c37aebbc

# ╔═╡ b7875dd0-68ff-475f-8c90-e3d33a72b12a
md"""
11.2 Crea una matriz diagonal de los eigenvalores de A
"""

# ╔═╡ 4a6152f4-78be-47e6-906a-25863ba85913

# ╔═╡ b6899c90-d9b9-4476-9971-96fef2aa14d0
md"""
11.3 Realiza un factorización de Hessenberg sobre la matriz A. Verifica que `A = QHQ'`.
"""

# ╔═╡ db8a09c4-2825-47c8-9697-1990d65b0b58

# ╔═╡ Cell order:
# ╟─b7a4d6c3-dc52-45c0-bcd6-dc1d32dbd302
# ╟─43882091-fc1b-409a-92d9-455ca7ad6f64
# ╠═85ad00b8-4b91-46c9-a130-4430f2986569
# ╠═5a8dd4eb-40a6-4750-b2e7-d44b30404363
# ╟─a98929e4-ce11-49b1-aeca-44ed3680df53
# ╟─e00b9a3c-3ffc-4aac-b5f7-a88a0daafba4
# ╟─a6bc3137-f8f6-44ec-bc7e-baa5a2f9c499
# ╠═0b6b51ee-ff9f-463e-95e1-2fd1f559bdf1
# ╟─2aefd8c5-abc3-4d1a-83fb-78f22190e0f3
# ╠═94c7fee6-2ada-4223-8fc0-57c7896181aa
# ╟─edc134d0-de77-429a-8bac-36ef19e80471
# ╠═a123c969-8e5c-432b-9bbd-3f1dfa4b8f31
# ╟─da792d5e-9ebf-4b93-9d57-61745021ec00
# ╠═46ab0bb7-ded7-4d98-a452-70242a5bf648
# ╟─8156b8f8-9919-42c7-83cd-49589bfd53e7
# ╠═608bfb30-6530-474b-9374-c5caa3a716f8
# ╟─cfd8ea62-ff51-4966-ad3e-2d12a937e71e
# ╠═289f4361-2962-48f4-8c47-915290a26035
# ╟─039fdbaf-9870-4efc-a54a-f457088c5e88
# ╠═d8217c6a-2f5f-4824-99e9-e56bfca07e92
# ╠═7d46f4fb-e092-409f-a0b0-5caad91f4022
# ╠═44aec045-5803-41d2-9574-31a9f14d9cf7
# ╟─466124e7-9978-4c7c-9cf1-6e34c7f55c33
# ╠═9a023eb2-ce3d-413e-a872-231325c556b6
# ╟─cf9efbb5-87ad-4795-b132-8877b334f85d
# ╠═45f9ca18-2ba7-4e8c-9087-ebadf2aa7101
# ╟─424a4d00-0264-45cf-bbba-43fb51ffa394
# ╠═ba2407dd-ef47-4532-a783-aafb17506e5f
# ╠═299bec2c-626f-40ca-a2f9-4e0450b23f3e
# ╟─0fa226f7-8a15-466b-9757-3c019fd4c58e
# ╠═92962e9b-13c2-465e-ae4d-b1c8a0172063
# ╟─e802da18-5118-4031-9be8-0fb90a7a0c9d
# ╠═215ea9e5-1a3a-4750-b249-2ed54f3ea9ac
# ╟─021e9680-4a33-469e-8014-ae3b831db370
# ╠═e55a8449-a507-4e90-9f1c-013edb26704d
# ╟─60ce24a1-2d71-4efc-95a8-3e210da5c53f
# ╠═fa8dca75-5aca-454d-b733-8da5b10c6857
# ╟─60918647-3d3d-4dec-a4e8-e0cf4348d072
# ╠═83e8b0eb-f8c9-4e4d-b319-27134985ffd1
# ╟─9aa2af05-bbb0-4570-8958-0ac08c5da7b5
# ╠═4ce42698-361a-4929-b03e-23b43ad04d33
# ╟─32879b97-8789-408f-a992-5f21dcd6e8fe
# ╠═4c3f2f14-3bd3-45c1-b36d-1000811fc622
# ╟─6a2732c8-241f-426c-8b51-ba130c3a4047
# ╠═45af09e4-8d42-404e-8e69-e0d185f4ea27
# ╟─8227c126-4417-4d6a-9c40-232ae467b17d
# ╠═cad5d92d-d7cb-444d-9664-e9eccbbb5b06
# ╟─6d694031-b7e4-423e-a341-e1e78e26d484
# ╠═ac57e6e1-986c-4efe-ad64-e2bb3a4aeb4c
# ╟─aa4e579d-5b59-494c-be51-dc65e6b57aa6
# ╟─61434524-e698-4a8d-a430-0126ac99157e
# ╠═04f4492e-5156-4526-b761-ab3c5a209646
# ╠═f8bb7598-095d-449a-86e2-f8bd6b9d8b94
# ╟─a8e4b120-54b2-4e16-b993-5a884c29309d
# ╠═c85c1276-a8d9-4f02-9fd7-0a448da34436
# ╟─9ea6c93a-e36f-4624-9598-7e4741efa50a
# ╟─6a578f8f-29fa-4583-a328-e8f33ba80738
# ╠═3e8c3380-7db5-4e35-b6f8-9b94dc417f80
# ╟─0c0caea9-3ee9-4e30-8570-4d30ed7cc14a
# ╠═3cde0ff8-2325-47b7-a1ff-d15155c2f634
# ╠═c4d93f45-2191-4624-ade9-277b546ced95
# ╟─b51b8521-6f77-410e-ae49-de138a1aec67
# ╠═50f4c300-e49b-4331-b364-582cc7692dd7
# ╟─65386dae-fb08-47d7-a81c-c3133a8403cf
# ╠═707539e7-c4b9-46c1-89f1-8d61a19413fa
# ╟─827e2055-5f83-44b1-a1dc-636b807c8b81
# ╠═1afd9e85-1f7d-48aa-baca-8cdf754e4d4b
# ╟─78a13f6a-ce68-4060-a1bf-1993ad5e93ad
# ╠═a5b49236-f31e-4449-941e-70e1677ae1b2
# ╟─9ad38b16-f9b9-4cde-bbfa-204d96ee3a43
# ╟─0139d8bd-7a60-4cdc-8844-91c0e0eaa9c6
# ╠═7ee1aa1c-7219-46d3-8abd-5a226c0f4abb
# ╟─f572bd5d-f1f4-4dce-a5be-293b4404bf71
# ╠═6846c3c4-6f80-43dc-9ece-07efe11133f8
# ╠═2b5bb260-8997-4b54-ba40-c19c0bd592a5
# ╟─7f270025-80df-4b91-babc-db876885e479
# ╠═876fb13e-d540-4931-ad32-f9491fd504c7
# ╠═e46ff282-c6a5-4917-8322-9ff2d1df8318
# ╟─b1093e1f-65fa-40d1-8b56-a374f15b3f51
# ╠═9d354a03-6229-46f1-80eb-cbed5d304d5b
# ╠═6e20dd10-3e4a-411e-bfdf-1534b509e8ce
# ╟─fb94650f-06f4-4b65-a598-45c820c27c95
# ╠═7773e25a-57e8-4917-9629-7329d7ddc60c
# ╟─2f6332ef-caef-423b-bc2d-4299436390f8
# ╠═b844544e-c513-42af-9b27-b006d50973c0
# ╟─f34e931a-c11b-4776-8d50-90812ef67d5b
# ╠═09b262f1-6a2d-4696-a194-54b65f129eb3
# ╟─44b64772-3daa-4037-9b6f-07d5feba3223
# ╠═dcef7621-ef7f-467e-9dec-75fae2bc44f3
# ╟─e451672e-e44d-45c7-a3c2-3badfcc83d40
# ╠═a994a69a-3879-49ad-a03a-c05d901b4072
# ╠═63d5a411-bd56-4ff6-bf35-a23995d4c2bb
# ╟─9902503c-0448-4bfc-bf6f-7c11321cca78
# ╠═9d099581-b364-4ad0-990a-d3151dcdcf96
# ╟─8ea1527f-a664-4520-9801-840e77bd6e62
# ╟─231d2771-137d-4930-80a6-4ca352e88aa5
# ╠═a0afdba1-25a0-409c-b84d-c63e53dbf343
# ╟─49cfb27c-8c00-4458-a77a-38dd5a4c643a
# ╠═0f4aec6b-6701-46f9-8a60-1c91f19cbb37
# ╟─1b26ea8e-d825-48ab-a0a0-a401dd83d42c
# ╟─8b718fad-0ded-42d2-915c-7868a8fc15df
# ╠═9858cb5f-ec34-4ce3-aa44-a676bf2d9d35
# ╟─4811bad4-fa07-4a22-92e6-c9947bce13d1
# ╠═a756d1d8-57ce-4e17-984f-7fbe1a8117e5
# ╠═e97ccbb7-fb81-4514-a8b6-0b177ef9620e
# ╠═a808d3f1-21a0-4d72-8323-9ffdb55084b6
# ╠═d9e7e3c7-cf57-4b86-b2f6-7b24b3ea8bdb
# ╟─6d5e0b76-5360-4e77-bd5a-93f5aac6011a
# ╠═57d7e0de-a3c9-4fda-b4b1-4e847bd5f8d2
# ╟─b29effd2-bb03-4f65-bd5c-ba5206378c54
# ╠═9e39d0ed-49b7-48c6-894e-6140c37aebbc
# ╟─b7875dd0-68ff-475f-8c90-e3d33a72b12a
# ╠═4a6152f4-78be-47e6-906a-25863ba85913
# ╟─b6899c90-d9b9-4476-9971-96fef2aa14d0
# ╠═db8a09c4-2825-47c8-9697-1990d65b0b58

max_stars_repo_name: miguelraz/JuliaTutorials

max_issues_repo_name: miguelraz/JuliaTutorials

max_forks_repo_name: miguelraz/JuliaTutorials


text:
using IntrinsicMetrics
using Base.Test
using StaticArrays
using Suppressor
using LightGraphs, SimpleWeightedGraphs
using NearestNeighbors
using Plots
include("manifolds.jl")

n = 200
repeats_small = 10

# Randomly generate datasets by sampling points from manifolds.
sphere = rand(Sphere(1), n)
atorus = rand(AsymTorus(5, 2, 1), n)
circle = rand(Circle(1), n)
cloud9 = rand(9, n)

# Datasets used in testing.
testsets = [sphere, atorus, circle, cloud9]

# Radius values used in testing.
radii = [0.5, 0.7]

# Metrics used in testing.
metrics = [Euclidean(), Minkowski(3), Cityblock(), Chebyshev()]

function isdistancematrix(m; report = true)
    if m ≉ m'
        report && println("Not symmetric!")
        return false
    elseif any(m .< 0)
        report && println("Not positive!")
        return false
    elseif any(diag(m) .≠ 0)
        report && println("Diagonal not zero!")
        return false
    end
    true
end

@testset "sanity checks" begin
    # Check sampling.
    @test size(testsets[1]) == (3, n)
    @test size(testsets[2]) == (3, n)
    @test size(testsets[3]) == (2, n)
    @test size(testsets[4]) == (9, n)

    # Check isdistancematrix.
    for _ in 1:repeats_small
        rmat = rand(n, n) - 0.5
        @test !isdistancematrix(rmat, report = false)
        symrmat = full(Symmetric(rmat))
        @test !isdistancematrix(symrmat, report = false)
        for i in 1:n
            symrmat[i, i] = 0
        end
        @test !isdistancematrix(symrmat, report = false)
        @test isdistancematrix(abs.(symrmat))
    end
end

@testset "IntrinsicMetric internals" begin
    @testset "ambient metric" begin
        for testset in testsets, r in radii, metric in metrics
            im = IntrinsicMetric(testset, r, metric = metric)
            @test ambientmetric(im) == metric
            @test isa(im, IntrinsicMetric{Float64, size(testset, 1), typeof(metric)})
       end
    end

    @testset "points" begin
        for testset in testsets, r in radii, metric in metrics
            im  = IntrinsicMetric(testset, r, metric = metric)
            dim = size(testset, 1)
            len = size(testset, 2)
            @test points(im) == reinterpret(SVector{dim, eltype(testset)},
                                            testset,
                                            (len, ))
            @test size(testset) == (ambientdim(im), npoints(im))
        end
    end

    @testset "graph properties" begin
        for testset in testsets, r in radii, metric in metrics
            im = IntrinsicMetric(testset, r, metric = metric)
            @test nv(adjgraph(im)) == npoints(im)
            @test isdistancematrix(adjgraph(im).weights)
        end
    end

    @testset "float types" begin
        for T in [Float16, Float32, Float64]
            testset = rand(Sphere(T(1.0)), n)
            im = IntrinsicMetric(testset, 30)
            @test eltype(testset) == T
            @test isa(im, IntrinsicMetric{T, size(testset, 1), Euclidean})
            @test isa(evaluate(im, [1.,0.,0.], [-1.,0.,0.]), T)
        end
    end

    @testset "show" begin
        for testset in testsets, r in radii, metric in metrics
            T = eltype(testset)
            D = size(testset, 1)
            M = typeof(metric)
            n = size(testset, 2)
            N = typeof(BallTree(testset, metric))

            im = IntrinsicMetric(testset, r, metric = metric)
            str = """
                  IntrinsicMetric{$T, $D, $M, SVector{$D,$T}, $N}
                    Number of points: $(length(im.points))
                    Dimensions: $D
                    Radius: $r
                    Metric: $metric
                    Graph: $(im.graph)
                  """
            @test @capture_out(println(im)) == str
            @test @capture_err(println(STDERR, im)) == str
        end
    end
end

@testset "evaluate" begin
    for testset in testsets, r in radii, metric in metrics
        im = IntrinsicMetric(testset, r, metric = metric)
        d1 = Float64[]
        d2 = Float64[]
        d3 = Float64[]
        ds = Float64[]
        # Pick 100 points and check distance properties.
        for _ in 1:100
            p1 = testset[:, rand(1:n)]
            p2 = testset[:, rand(1:n)]
            p3 = testset[:, rand(1:n)]
            push!(d1, evaluate(im, p1, p2))
            push!(ds, evaluate(im, p2, p1))
            push!(d2, evaluate(im, p2, p3))
            push!(d3, evaluate(im, p1, p3))
        end
        # Symmetry.
        @test d1 ≈ ds
        # Non-negativity.
        @test all(d1 .≥ 0)
        # Triangle inequality.
        @test all(d3 .≤ d1 .+ d2 .+ sqrt(eps(Float64)))
    end

    # Going around sphere should be about π.
    im = IntrinsicMetric(hcat(rand(Sphere(1), 1000), [0, 0, 1], [0, 0,-1]), 0.5)
    @test evaluate(im, [0, 0,-1], [0, 0, 1]) ≈ π     atol=0.1
    @test evaluate(im, [0, 0,-2], [0, 0, 1]) ≈ (1+π) atol=0.1
    @test evaluate(im, [0, 0,-1], [0, 0, 2]) ≈ (1+π) atol=0.1
    @test evaluate(im, [0, 0,-2], [0, 0, 2]) ≈ (2+π) atol=0.1

    # Same for circle.
    im = IntrinsicMetric(rand(Circle(1), 1000), 0.1)
    @test evaluate(im, [0,-1], [ 0, 1]) ≈ π atol=0.01
    @test evaluate(im, [1, 0], [-1, 0]) ≈ π atol=0.01
end

@testset "pairwise - size errors" begin
    for testset in testsets, r in radii, metric in metrics
        im = IntrinsicMetric(testset, r, metric = metric)
        d  = size(testset, 1)
        dm = DimensionMismatch

        @test_throws dm evaluate(im, rand(d-1, 1), rand(d,   1))
        @test_throws dm evaluate(im, rand(d,   1), rand(d-1, 1))
        @test_throws dm evaluate(im, rand(d+1, 1), rand(d+1, 1))

        @test_throws dm pairwise(im, rand(d,   10), rand(d+1, 11))
        @test_throws dm pairwise(im, rand(d-1, 10), rand(d-1, 11))

        @test_throws dm pairwise!(zeros(5, 5), im, rand(d-1, 5), rand(d,   5))
        @test_throws dm pairwise!(zeros(5, 5), im, rand(d-1, 5), rand(d-1, 5))
        @test_throws dm pairwise!(zeros(6, 5), im, rand(d,   5), rand(d,   5))
        @test_throws dm pairwise!(zeros(5, 6), im, rand(d,   5), rand(d,   5))

        @test_throws dm pairwise(im, rand(d+1, 10))

        @test_throws dm pairwise!(zeros(10, 11), im, rand(d,   10))
        @test_throws dm pairwise!(zeros(11, 11), im, rand(d+1, 11))

        @test_throws dm pairwise!(zeros(10, 10), im)
    end
end

@testset "pairwise - one matrix" begin
    for testset in testsets, r in radii, metric in metrics
        im = IntrinsicMetric(testset, r, metric = metric)

        # Points from sample.
        pw1 = pairwise(im, testset, testset)
        pw2 = pairwise(im, testset)
        pw3 = pairwise(im)
        @test isdistancematrix(pw1)
        @test isdistancematrix(pw2)
        @test isdistancematrix(pw3)
        @test pw1 ≈ pw2 ≈ pw3
        @test size(pw1) == size(pw2) == size(pw3) == (n, n)

        # Random points.
        pts = (rand(size(testset, 1), 100) - 0.5) * 3
        pw4 = pairwise(im, pts, pts)
        pw5 = pairwise(im, pts)
        @test isdistancematrix(pw4)
        @test isdistancematrix(pw5)
        @test pw4 ≈ pw5
        @test size(pw4) == size(pw5) == (100, 100)

        d1 = Float64[]; d2 = Float64[]; d3 = Float64[]
        for _ in 1:100
            i = rand(1:n)
            j = rand(1:n)
            k = rand(1:n)
            push!(d1, pw2[i, j])
            push!(d2, pw2[j, k])
            push!(d3, pw2[i, k])
        end
        @test all(d3 .≤ d1 .+ d2 .+ sqrt(eps(Float64)))
    end
end

@testset "pairwise - two matrices" begin
    for testset in testsets, r in radii, metric in metrics
        im = IntrinsicMetric(testset, r, metric = metric)
        pts1 = (rand(size(testset, 1), 100) - 0.5) * 3
        pts2 = (rand(size(testset, 1), 1000) - 0.5) * 3

        pw1 = pairwise(im, pts1, pts2)
        pw2 = pairwise(im, pts2, pts1)

        @test pw1' ≈ pw2
        @test all(pw1 .≥ 0)
        @test size(pw1) == (100, 1000)
    end
end

@testset "pairwise - full graph" begin
    for _ in 1:repeats_small, metric in metrics
        pts = rand(5, n)
        im = IntrinsicMetric(pts, 10, metric = metric)

        @test adjgraph(im).weights ≈ pairwise(im, pts) ≈ pairwise(metric, pts)
    end
end

@testset "plotting" begin
    r = radii[1]
    for testset in testsets[1:3], metric in metrics
        @test plot(IntrinsicMetric(testset, r, metric = metric)) ≠ nothing
    end
    @test_throws ErrorException plot(IntrinsicMetric(testsets[4], r))
end


text:
# code reflection entry-points

using InteractiveUtils

# Return the capability of the current context's device, or a sane fall-back.
function current_capability()
    fallback = minimum(target_support)
    if !initialized[]
        return fallback
    end

    ctx = CuCurrentContext()
    if ctx == nothing
        return fallback
    end

    return capability(device(ctx))
end

#
# code_* replacements
#

"""
    code_llvm([io], f, types; optimize=true, dump_module=false, cap::VersionNumber)

Prints the device LLVM IR generated for the method matching the given generic function and
type signature to `io` which defaults to `stdout`. The IR is optimized according to
`optimize` (defaults to true), and the entire module, including headers and other functions,
is dumped if `dump_module` is set (defaults to false). The device capability `cap` to
generate code for defaults to the current active device's capability, or v"2.0" if there is
no such active context.

See also: [`@device_code_llvm`](@ref), [`Base.code_llvm`](@ref)
"""
function code_llvm(io::IO, @nospecialize(func::Core.Function), @nospecialize(types=Tuple);
                   optimize::Bool=true, dump_module::Bool=false,
                   cap::VersionNumber=current_capability(), kernel::Bool=false, kwargs...)
    tt = Base.to_tuple_type(types)
    ctx = CompilerContext(func, tt, cap, kernel; kwargs...)
    validate_invocation(ctx)

    mod, entry = irgen(ctx)
    if kernel
        entry = promote_kernel!(ctx, mod, entry)
    end
    if optimize
        optimize!(ctx, mod, entry)
    end
    if dump_module
        show(io, mod)
    else
        show(io, entry)
    end
end
code_llvm(@nospecialize(func), @nospecialize(types=Tuple); kwargs...) =
    code_llvm(stdout, func, types; kwargs...)

"""
    code_ptx([io], f, types; cap::VersionNumber, kernel::Bool=false)

Prints the PTX assembly generated for the method matching the given generic function and
type signature to `io` which defaults to `stdout`. The device capability `cap` to generate
code for defaults to the current active device's capability, or v"2.0" if there is no such
active context. The optional `kernel` parameter indicates whether the function in question
is an entry-point function, or a regular device function.

See also: [`@device_code_ptx`](@ref)
"""
function code_ptx(io::IO, @nospecialize(func::Core.Function), @nospecialize(types=Tuple);
                  cap::VersionNumber=current_capability(), kernel::Bool=false, kwargs...)
    tt = Base.to_tuple_type(types)
    ctx = CompilerContext(func, tt, cap, kernel; kwargs...)
    validate_invocation(ctx)

    ptx,_ = compile_function(ctx)
    # TODO: this code contains all the functions in the call chain,
    #       is it possible to implement `dump_module`?
    print(io, ptx)
end
code_ptx(@nospecialize(func), @nospecialize(types=Tuple); kwargs...) =
    code_ptx(stdout, func, types; kwargs...)

"""
    code_sass([io], f, types, cap::VersionNumber)

Prints the SASS code generated for the method matching the given generic function and type
signature to `io` which defaults to `stdout`. The device capability `cap` to generate code
for defaults to the current active device's capability, or v"2.0" if there is no such active
context. The method needs to be a valid entry-point kernel, eg. it should not return any
values.

See also: [`@device_code_sass`](@ref)
"""
function code_sass(io::IO, @nospecialize(func::Core.Function), @nospecialize(types=Tuple);
                   cap::VersionNumber=current_capability(), kernel::Bool=true, kwargs...)
    if !kernel
        error("Can only generate SASS code for kernel functions")
    end
    if ptxas === nothing || cuobjdump === nothing
        error("Your CUDA installation does not provide ptxas or cuobjdump, both of which are required for code_sass")
    end

    tt = Base.to_tuple_type(types)
    ctx = CompilerContext(func, tt, cap, kernel; kwargs...)
    validate_invocation(ctx)

    ptx,_ = compile_function(ctx)

    fn = tempname()
    gpu = "sm_$(cap.major)$(cap.minor)"
    # NOTE: this might not match what is being executed, due to the PTX->SASS conversion
    #       by the driver possibly not matching what `ptxas` (part of the toolkit) does.
    # TODO: see how `nvvp` extracts SASS code when doing PC sampling, and copy that.
    Base.run(`$ptxas --gpu-name $gpu --output-file $fn --input-as-string $ptx`)
    try
        print(io, read(`$cuobjdump --dump-sass $fn`, String))
    finally
        rm(fn)
    end
end
code_sass(@nospecialize(func), @nospecialize(types=Tuple); kwargs...) =
    code_sass(stdout, func, types; kwargs...)

#
# @device_code_* functions
#

export @device_code_lowered, @device_code_typed, @device_code_warntype,
       @device_code_llvm, @device_code_ptx, @device_code_sass,
       @device_code

function emit_hooked_compilation(inner_hook, ex...)
    user_code = ex[end]
    user_kwargs = ex[1:end-1]
    quote
        # wipe the compile cache to force recompilation
        empty!(CUDAnative.compilecache)

        local kernels = 0
        function outer_hook(ctx)
            kernels += 1
            $inner_hook(ctx; $(map(esc, user_kwargs)...))
        end

        if CUDAnative.compile_hook[] != nothing
            error("Chaining multiple @device_code calls is unsupported")
        end
        try
            CUDAnative.compile_hook[] = outer_hook
            $(esc(user_code))
        finally
            CUDAnative.compile_hook[] = nothing
        end

        if kernels == 0
            error("no kernels executed while evaluating the given expression")
        end

        nothing
    end
end

# NOTE: these hooks take both a `f` and an inner `f`, because of how `@cuda`/`_cuda` work:
#       kernels are automatically wrapper in a function returning nothing, for usability.
#
#       Julia-level reflection (lowered/typed/warntype) skips these wrapper, because we
#       can't do call-site inlining and the kernel wrapper would hide any meaningful code.
#
#       at the LLVM level, we inline everything so there's no need to hide the wrapper.

"""
    @device_code_lowered ex

Evaluates the expression `ex` and returns the result of [`Base.code_lowered`](@ref) for
every compiled CUDA kernel.

See also: [`Base.@code_lowered`](@ref)
"""
macro device_code_lowered(ex...)
    quote
        buf = Any[]
        function hook(ctx::CompilerContext)
            append!(buf, code_lowered(something(ctx.inner_f, ctx.f), ctx.tt))
        end
        $(emit_hooked_compilation(:hook, ex...))
        buf
    end
end

"""
    @device_code_typed ex

Evaluates the expression `ex` and returns the result of [`Base.code_typed`](@ref) for
every compiled CUDA kernel.

See also: [`Base.@code_typed`](@ref)
"""
macro device_code_typed(ex...)
    quote
        buf = Any[]
        function hook(ctx::CompilerContext)
            append!(buf, code_typed(something(ctx.inner_f, ctx.f), ctx.tt))
        end
        $(emit_hooked_compilation(:hook, ex...))
        buf
    end
end

"""
    @device_code_warntype [io::IO=stdout] ex

Evaluates the expression `ex` and prints the result of [`Base.code_warntype`](@ref) to `io`
for every compiled CUDA kernel.

See also: [`Base.@code_warntype`](@ref)
"""
macro device_code_warntype(ex...)
    function hook(ctx::CompilerContext; io::IO=stdout, kwargs...)
        code_warntype(io, something(ctx.inner_f, ctx.f), ctx.tt; kwargs...)
    end
    emit_hooked_compilation(hook, ex...)
end

"""
    @device_code_llvm [io::IO=stdout, ...] ex

Evaluates the expression `ex` and prints the result of [`Base.code_llvm`](@ref) to `io` for
every compiled CUDA kernel. For other supported keywords, see
[`CUDAnative.code_llvm`](@ref).

See also: [`Base.@code_llvm`](@ref)
"""
macro device_code_llvm(ex...)
    function hook(ctx::CompilerContext; io::IO=stdout, kwargs...)
        code_llvm(io, ctx.f, ctx.tt; kernel=ctx.kernel, cap=ctx.cap, kwargs...)
    end
    emit_hooked_compilation(hook, ex...)
end

"""
    @device_code_ptx [io::IO=stdout, ...] ex

Evaluates the expression `ex` and prints the result of [`CUDAnative.code_ptx`](@ref) to `io`
for every compiled CUDA kernel. For other supported keywords, see
[`CUDAnative.code_ptx`](@ref).
"""
macro device_code_ptx(ex...)
    function hook(ctx::CompilerContext; io::IO=stdout, kwargs...)
        code_ptx(io, ctx.f, ctx.tt; kernel=ctx.kernel, cap=ctx.cap, kwargs...)
    end
    emit_hooked_compilation(hook, ex...)
end

"""
    @device_code_sass [io::IO=stdout, ...] ex

Evaluates the expression `ex` and prints the result of [`CUDAnative.code_sass`](@ref) to
`io` for every compiled CUDA kernel. For other supported keywords, see
[`CUDAnative.code_sass`](@ref).
"""
macro device_code_sass(ex...)
    function hook(ctx::CompilerContext; io::IO=stdout, kwargs...)
        # we have inlined every function using LLVM, so don't hide the kernel wrapper.
        code_sass(io, ctx.f, ctx.tt; cap=ctx.cap, kwargs...)
    end
    emit_hooked_compilation(hook, ex...)
end

"""
    @device_code dir::AbstractString=... [...] ex

Evaluates the expression `ex` and dumps all intermediate forms of code to the directory
`dir`.
"""
macro device_code(ex...)
    only(xs) = (@assert length(xs) == 1; first(xs))
    function hook(ctx::CompilerContext; dir::AbstractString)
        fn = "$(typeof(something(ctx.inner_f, ctx.f)).name.mt.name)_$(globalUnique+1)"
        mkpath(dir)

        open(joinpath(dir, "$fn.lowered.jl"), "w") do io
            code = only(code_lowered(something(ctx.inner_f, ctx.f), ctx.tt))
            println(io, code)
        end

        open(joinpath(dir, "$fn.typed.jl"), "w") do io
            code = only(code_typed(something(ctx.inner_f, ctx.f), ctx.tt))
            println(io, code)
        end

        open(joinpath(dir, "$fn.ll"), "w") do io
            code_llvm(io, ctx.f, ctx.tt; kernel=ctx.kernel, cap=ctx.cap, dump_module=true)
        end

        open(joinpath(dir, "$fn.ptx"), "w") do io
            code_ptx(io, ctx.f, ctx.tt; kernel=ctx.kernel, cap=ctx.cap)
        end

        open(joinpath(dir, "$fn.sass"), "w") do io
            code_sass(io, ctx.f, ctx.tt; cap=ctx.cap)
        end
    end
    emit_hooked_compilation(hook, ex...)
end


text:
import VectorSphericalWaves.wignerdjmn_recurrence_memoize
using HDF5
using JLD
using Interpolations

create_lookup_table = true
if create_lookup_table
    const wignerdjmn_VS_theta_cache = Dict()
    n_θ_points = Int(1e4);
    θ_array = LinRange(0, π, n_θ_points);
    for s = 1:6
        for m = -s:s
            for n = -s:s
                wignerdjmn_VS_theta_cache[s,m,n] = wignerdjmn_recurrence_memoize.(s, m, n, θ_array)
            end
        end
    end
    save("cach.jld", "wignerdjmn", wignerdjmn_VS_theta_cache, "θ_array", θ_array)

end
"""
    Calculate Wigner-d function, using recurrence and memoize, and saves the cache in a lookup table
"""
const wignerdjmn_VS_theta_cache = load("cach.jld")["wignerdjmn"]
const θ_array_interp = load("cach.jld")["θ_array"]
function wignerdjmn_recurrence_memoize_lookup(s::Int, m::Int, n::Int, θ::R) where R <: Real
    wignerdjmn_array = wignerdjmn_VS_theta_cache[s,m,n]
    # itp = interpolate(θ_array_interp, wignerdjmn_array, BSpline(Linear()))
    itp = LinearInterpolation(θ_array_interp, wignerdjmn_array)
    return itp(θ)
end
s,m,n = 5, 2, 1
t_ = rand(10000);
@time VectorSphericalWaves.wignerdjmn_ELZOUKA.(s,m,n,t_);
@time VectorSphericalWaves.wignerdjmn_recurrence_memoize.(s,m,n,t_);
@time wignerdjmn_recurrence_memoize_lookup.(s,m,n,t_);

"""
Surprisingly, the direct calculation of WignerD is faster than memoizing!
"""


text:
using Test

include("../../src/io.jl")

@testset "ChunkedOuputBuffer Tests" begin

    @testset "test_should_start_empty" begin
        buffer = ChunkedOutputBuffer()
        @test view(buffer) == b""
    end

    @testset "test_should_be_able_to_set_max_chunk_size" begin
        buffer = ChunkedOutputBuffer(max_chunk_size=4)
        @test buffer.max_chunk_size == 4
    end

    @testset "test_small_data_should_be_directly_appended" begin
        buffer = ChunkedOutputBuffer()

        write(buffer, b"hello")

        @test view(buffer) == b"\x00\x05hello"
    end

    @testset "test_overflow_data_should_use_a_new_chunk" begin
        buffer = ChunkedOutputBuffer(max_chunk_size=6)

        write(buffer, b"over")
        write(buffer, b"flow")

        @test view(buffer) == b"\x00\x04over\x00\x04flow"
    end

    @testset "test_big_data_should_be_split_across_chunks" begin
        buffer = ChunkedOutputBuffer(max_chunk_size=2)

        write(buffer, b"octopus")

        @test view(buffer) == b"\x00\x02oc\x00\x02to\x00\x02pu\x00\x01s"
    end

    @testset "test_clear_should_clear_everything" begin
        buffer = ChunkedOutputBuffer()

        write(buffer, b"redacted")
        clear(buffer)

        @test view(buffer) == b""
    end

    @testset "test_cleared_buffer_should_be_reusable" begin
        buffer = ChunkedOutputBuffer()

        write(buffer, b"Windows")
        clear(buffer)
        write(buffer, b"Linux")

        @test view(buffer) == b"\x00\x05Linux"
    end

    @testset "test_should_be_able_to_force_chunks" begin
        buffer = ChunkedOutputBuffer()

        write(buffer, b"hello")
        chunk(buffer)
        write(buffer, b"world")
        chunk(buffer)

        @test view(buffer) == b"\x00\x05hello\x00\x05world"
    end

    @testset "test_should_be_able_to_force_empty_chunk" begin
        buffer = ChunkedOutputBuffer()

        write(buffer, b"hello")
        chunk(buffer)
        chunk(buffer)

        @test view(buffer) == b"\x00\x05hello\x00\x00"
    end
end

max_stars_repo_name: virtualgraham/JuliaBolt

max_issues_repo_name: virtualgraham/JuliaBolt

max_forks_repo_name: virtualgraham/JuliaBolt


text:
# put all using and exports in this package file
# so know where all things coming and leaving are (@31:00 on video)

module aI_Test

using ForwardDiff
using WAV

include("file01.jl")

export my_f, derivative_of_my_f

end # module


text:
using HorizonSideRobots

function draw!(r::Robot, side1::HorizonSide, side2::HorizonSide)
    while (isborder(r, side1) || isborder(r, side2)) == false
        move!(r, side1)
        move!(r, side2)
        putmarker!(r)
    end
end

function back!(r::Robot, side1::HorizonSide, side2::HorizonSide)
    while ismarker(r)
        move!(r, side1)
        move!(r, side2)
    end
end

function diag!(r::Robot)
    vec = West
    for side in [Nord, Sud]
        draw!(r, vec, side)
        back!(r, inverse(vec), inverse(side))
    end

    vec = Ost
    for s in [Nord, Sud]
        draw!(r, vec, side)
        back!(r, inverse(vec), inverse(side))
    end

    putmarker!(r)
end

diag!(r)


text:
## core libc calls ##

@unix_only begin
    _getenv(var::String) = ccall(:getenv, Ptr{Uint8}, (Ptr{Uint8},), var)
    _hasenv(s::String) = _getenv(s) != C_NULL
end
@windows_only begin
_getenvlen(var::String) = ccall(:GetEnvironmentVariableA,stdcall,Uint32,(Ptr{Uint8},Ptr{Uint8},Uint32),var,C_NULL,0)
_hasenv(s::String) = _getenvlen(s)!=0
function _jl_win_getenv(s::String,len::Uint32)
    val=zeros(Uint8,len-1)
    ret=ccall(:GetEnvironmentVariableA,stdcall,Uint32,(Ptr{Uint8},Ptr{Uint8},Uint32),s,val,len)
    if ret==0||ret!=len-1 #Trailing 0 is only included on first call to GetEnvA
        error("unknown system error: ", s, len, ret)
    end
    val
end
end

macro accessEnv(var,errorcase)
@unix_only return quote
     val=_getenv($(esc(var)))
     if val == C_NULL
        $(esc(errorcase))
     end
     bytestring(val)
end
@windows_only return quote
    len=_getenvlen($(esc(var)))
    if len == 0
        $(esc(errorcase))
    end
    bytestring(_jl_win_getenv($(esc(var)),len))
end
end

function _setenv(var::String, val::String, overwrite::Bool)
@unix_only begin
    ret = ccall(:setenv, Int32, (Ptr{Uint8},Ptr{Uint8},Int32), var, val, overwrite)
    systemerror(:setenv, ret != 0)
end
@windows_only begin
    if overwrite||!_hasenv(var)
        ret = ccall(:SetEnvironmentVariableA,stdcall,Int32,(Ptr{Uint8},Ptr{Uint8}),var,val)
        systemerror(:setenv, ret == 0)
    end
end
end

_setenv(var::String, val::String) = _setenv(var, val, true)

function _unsetenv(var::String)
@unix_only begin
    ret = ccall(:unsetenv, Int32, (Ptr{Uint8},), var)
    systemerror(:unsetenv, ret != 0)
end
@windows_only begin
    ret = ccall(:SetEnvironmentVariableA,stdcall,Int32,(Ptr{Uint8},Ptr{Uint8}),var,C_NULL)
    systemerror(:setenv, ret == 0)
end
end

## ENV: hash interface ##

@unix_only type EnvHash <: Associative{ByteString,ByteString}; end
@windows_only type EnvHash <: Associative{ByteString,ByteString}
    block::Ptr{Uint8}
    EnvHash() = new(C_NULL)
end
const ENV = EnvHash()

similar(::EnvHash) = Dict{ByteString,ByteString}()

getindex(::EnvHash, k::String) = @accessEnv k throw(KeyError(k))
get(::EnvHash, k::String, def) = @accessEnv k (return def)
in(k::String, ::KeyIterator{EnvHash}) = _hasenv(k)
pop!(::EnvHash, k::String) = (v = ENV[k]; _unsetenv(k); v)
pop!(::EnvHash, k::String, def) = haskey(ENV,k) ? pop!(ENV,k) : def
function delete!(::EnvHash, k::String)
    warn_once("""
        delete!(ENV,key) now returns the modified environment.
        Use pop!(ENV,key) to retrieve the value instead.
        """)
    _unsetenv(k)
    ENV
end
delete!(::EnvHash, k::String, def) = haskey(ENV,k) ? delete!(ENV,k) : def
setindex!(::EnvHash, v, k::String) = _setenv(k,string(v))
push!(::EnvHash, k::String, v) = setindex!(ENV, v, k)

@unix_only begin
start(::EnvHash) = 0
done(::EnvHash, i) = (ccall(:jl_environ, Any, (Int32,), i) == nothing)

function next(::EnvHash, i)
    env = ccall(:jl_environ, Any, (Int32,), i)
    if env == nothing
        error(BoundsError)
    end
    env::ByteString
    m = match(r"^(.*?)=(.*)$"s, env)
    if m == nothing
        error("malformed environment entry: $env")
    end
    (ByteString[convert(typeof(env),x) for x in m.captures], i+1)
end
end

@windows_only begin
start(hash::EnvHash) = (hash.block = ccall(:GetEnvironmentStringsA,stdcall,Ptr{Uint8},()))
function done(hash::EnvHash, pos::Ptr{Uint8})
    if ccall(:jl_env_done,Any,(Ptr{Uint8},),pos)::Bool
        ccall(:FreeEnvironmentStringsA,stdcall,Int32,(Ptr{Uint8},),hash.block)
        hash.block=C_NULL
        return true
    end
    false
end
function next(hash::EnvHash, pos::Ptr{Uint8})
    len = ccall(:strlen, Uint, (Ptr{Uint8},), pos)
    env=ccall(:jl_pchar_to_string, Any, (Ptr{Uint8},Int), pos, len)::ByteString
    m = match(r"^(.*?)=(.*)$"s, env)
    if m == nothing
        error("malformed environment entry: $env")
    end
    (ByteString[convert(typeof(env),x) for x in m.captures], pos+len+1)
end
end

#TODO: Make these more efficent
function length(::EnvHash)
    i = 0
    for (k,v) in ENV
        i += 1
    end
    return i
end

function show(io::IO, ::EnvHash)
    for (k,v) = ENV
        println(io, "$k=$v")
    end
end

# temporarily set and then restore an environment value
function with_env(f::Function, key::String, val)
    old = get(ENV,key,nothing)
    val != nothing ? (ENV[key]=val) : _unsetenv(key)
    try f()
    finally
        old != nothing ? (ENV[key]=old) : _unsetenv(key)
    catch
        rethrow()
    end
end

## misc environment-related functionality ##

tty_cols() = parseint(Int32, get(ENV,"COLUMNS","80"), 10)
tty_rows() = parseint(Int32, get(ENV,"LINES","25"), 10)

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
# Accessors
getLabelDict(dfg::LightDFG) = copy(dfg.g.labels.sym_int)
getDescription(dfg::LightDFG) = dfg.description
setDescription(dfg::LightDFG, description::String) = dfg.description = description
getAddHistory(dfg::LightDFG) = dfg.addHistory
getSolverParams(dfg::LightDFG) = dfg.solverParams

# setSolverParams(dfg::LightDFG, solverParams) = dfg.solverParams = solverParams
function setSolverParams(dfg::LightDFG, solverParams::P) where P <: AbstractParams
  dfg.solverParams = solverParams
end

# Get user, robot, and session "small" data.
getUserData(dfg::LightDFG)::Dict{Symbol, String} = return dfg.userData
function setUserData(dfg::LightDFG, data::Dict{Symbol, String})::Bool
	dfg.userData = data
	return true
end
getRobotData(dfg::LightDFG)::Dict{Symbol, String} = return dfg.robotData
function setRobotData(dfg::LightDFG, data::Dict{Symbol, String})::Bool
	dfg.robotData = data
	return true
end
getSessionData(dfg::LightDFG)::Dict{Symbol, String} = return dfg.sessionData
function setSessionData(dfg::LightDFG, data::Dict{Symbol, String})::Bool
	dfg.sessionData = data
	return true
end

"""
    $(SIGNATURES)
True if the variable or factor exists in the graph.
"""
function exists(dfg::LightDFG{P,V,F}, node::V) where {P <: AbstractParams, V <: AbstractDFGVariable, F <: AbstractDFGFactor}
	return haskey(dfg.g.variables, node.label)
end

function exists(dfg::LightDFG{P,V,F}, node::F) where {P <: AbstractParams, V <: AbstractDFGVariable, F <: AbstractDFGFactor}
	return haskey(dfg.g.factors, node.label)
end

exists(dfg::LightDFG, nId::Symbol) = haskey(dfg.g.labels, nId)

exists(dfg::LightDFG, node::DFGNode) = exists(dfg, node.label)

"""
    $(SIGNATURES)
Add a DFGVariable to a DFG.
"""
function addVariable!(dfg::LightDFG{<:AbstractParams, V, <:AbstractDFGFactor}, variable::V)::Bool where V <: AbstractDFGVariable
	#TODO should this be an error
	if haskey(dfg.g.variables, variable.label)
		error("Variable '$(variable.label)' already exists in the factor graph")
	end

	#NOTE Internal ID always set to zero as it is not needed?
    # variable._internalId = 0

	FactorGraphs.addVariable!(dfg.g, variable) || return false

    # Track insertion
    push!(dfg.addHistory, variable.label)

    return true
end

"""
    $(SIGNATURES)
Add a DFGFactor to a DFG.
"""
function addFactor!(dfg::LightDFG{<:AbstractParams, V, F}, variables::Vector{V}, factor::F)::Bool where {V <: AbstractDFGVariable, F <: AbstractDFGFactor}
    # if haskey(dfg.g.metaindex[:label], factor.label)
    #     error("Factor '$(factor.label)' already exists in the factor graph")
    # end
	#TODO should this be an error
	if haskey(dfg.g.factors, factor.label)
		error("Factor '$(factor.label)' already exists in the factor graph")
	end
    # for v in variables
    #     if !(v.label in keys(dfg.g.metaindex[:label]))
    #         error("Variable '$(v.label)' not found in graph when creating Factor '$(factor.label)'")
    #     end
    # end

	variableLabels = map(v->v.label, variables)

	resize!(factor._variableOrderSymbols, length(variableLabels))
	factor._variableOrderSymbols .= variableLabels
    # factor._variableOrderSymbols = copy(variableLabels)

    return FactorGraphs.addFactor!(dfg.g, variableLabels, factor)
end

"""
    $(SIGNATURES)
Add a DFGFactor to a DFG.
"""
function addFactor!(dfg::LightDFG{<:AbstractParams, <:AbstractDFGVariable, F}, variableLabels::Vector{Symbol}, factor::F)::Bool where F <: AbstractDFGFactor
	#TODO should this be an error
	if haskey(dfg.g.factors, factor.label)
        error("Factor '$(factor.label)' already exists in the factor graph")
    end

	resize!(factor._variableOrderSymbols, length(variableLabels))
	factor._variableOrderSymbols .= variableLabels
	# factor._variableOrderSymbols = copy(variableLabels)

    return FactorGraphs.addFactor!(dfg.g, variableLabels, factor)
end

function addFactor!(dfg::LightDFG{<:AbstractParams, <:AbstractDFGVariable, F}, factor::F)::Bool where F <: AbstractDFGFactor
	#TODO should this be an error
	if haskey(dfg.g.factors, factor.label)
        error("Factor '$(factor.label)' already exists in the factor graph")
    end

    return FactorGraphs.addFactor!(dfg.g, variableLabels, factor)
end

"""
    $(SIGNATURES)
Get a DFGVariable from a DFG using its label.
"""
function getVariable(dfg::LightDFG, label::Symbol)::AbstractDFGVariable
    return dfg.g.variables[label]
end

"""
    $(SIGNATURES)
Get a DFGFactor from a DFG using its label.
"""
function getFactor(dfg::LightDFG, label::Symbol)::AbstractDFGFactor
    return dfg.g.factors[label]
end

"""
    $(SIGNATURES)
Update a complete DFGVariable in the DFG.
"""
function updateVariable!(dfg::LightDFG, variable::V)::V where V <: AbstractDFGVariable
    if !haskey(dfg.g.variables, variable.label)
        error("Variable label '$(variable.label)' does not exist in the factor graph")
    end
	dfg.g.variables[variable.label] = variable
    return variable
end

"""
    $(SIGNATURES)
Update a complete DFGFactor in the DFG.
"""
function updateFactor!(dfg::LightDFG, factor::F)::F where F <: AbstractDFGFactor
    if !haskey(dfg.g.factors, factor.label)
        error("Factor label '$(factor.label)' does not exist in the factor graph")
    end
	dfg.g.factors[factor.label] = factor
    return factor
end

"""
    $(SIGNATURES)
Delete a DFGVariable from the DFG using its label.
"""
function deleteVariable!(dfg::LightDFG, label::Symbol)::AbstractDFGVariable
    if !haskey(dfg.g.variables, label)
        error("Variable label '$(label)' does not exist in the factor graph")
    end
	variable = dfg.g.variables[label]
	rem_vertex!(dfg.g, dfg.g.labels[label])

    return variable
end

"""
    $(SIGNATURES)
Delete a DFGFactor from the DFG using its label.
"""
function deleteFactor!(dfg::LightDFG, label::Symbol)::AbstractDFGFactor
    if !haskey(dfg.g.factors, label)
        error("Factor label '$(label)' does not exist in the factor graph")
    end
	factor = dfg.g.factors[label]
	variable = rem_vertex!(dfg.g,  dfg.g.labels[label])
    return factor
end

"""
    $(SIGNATURES)
List the DFGVariables in the DFG.
Optionally specify a label regular expression to retrieves a subset of the variables.
"""
function getVariables(dfg::LightDFG, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[])::Vector{AbstractDFGVariable}

	# variables = map(v -> v.dfgNode, filter(n -> n.dfgNode isa DFGVariable, vertices(dfg.g)))
	variables = collect(values(dfg.g.variables))
    if regexFilter != nothing
        variables = filter(v -> occursin(regexFilter, String(v.label)), variables)
    end
	if length(tags) > 0
        mask = map(v -> length(intersect(v.tags, tags)) > 0, variables )
        return variables[mask]
    end
	return variables
end

"""
    $(SIGNATURES)
Get a list of IDs of the DFGVariables in the DFG.
Optionally specify a label regular expression to retrieves a subset of the variables.

Example
```julia
getVariableIds(dfg, r"l", tags=[:APRILTAG;])
```

Related

ls
"""

function getVariableIds(dfg::LightDFG, regexFilter::Union{Nothing, Regex}=nothing; tags::Vector{Symbol}=Symbol[])::Vector{Symbol}

	# variables = map(v -> v.dfgNode, filter(n -> n.dfgNode isa DFGVariable, vertices(dfg.g)))
	if length(tags) > 0
		return map(v -> v.label, getVariables(dfg, regexFilter, tags=tags))
	else
		variables = collect(keys(dfg.g.variables))
		regexFilter != nothing && (variables = filter(v -> occursin(regexFilter, String(v)), variables))
		return variables
    end
end

"""
    $(SIGNATURES)
List the DFGFactors in the DFG.
Optionally specify a label regular expression to retrieves a subset of the factors.
"""
function getFactors(dfg::LightDFG, regexFilter::Union{Nothing, Regex}=nothing)::Vector{AbstractDFGFactor}
	# factors = map(v -> v.dfgNode, filter(n -> n.dfgNode isa DFGFactor, vertices(dfg.g)))
	factors = collect(values(dfg.g.factors))
	if regexFilter != nothing
		factors = filter(f -> occursin(regexFilter, String(f.label)), factors)
	end
	return factors
end

"""
    $(SIGNATURES)
Get a list of the IDs of the DFGFactors in the DFG.
Optionally specify a label regular expression to retrieves a subset of the factors.
"""
function getFactorIds(dfg::LightDFG, regexFilter::Union{Nothing, Regex}=nothing)::Vector{Symbol}
	# factors = map(v -> v.dfgNode, filter(n -> n.dfgNode isa DFGFactor, vertices(dfg.g)))
	factors = collect(keys(dfg.g.factors))
	if regexFilter != nothing
		factors = filter(f -> occursin(regexFilter, String(f)), factors)
	end
	return factors
end

"""
    $(SIGNATURES)
Checks if the graph is fully connected, returns true if so.
"""
function isFullyConnected(dfg::LightDFG)::Bool
    return length(LightGraphs.connected_components(dfg.g)) == 1
end

function _isready(dfg::LightDFG, label::Symbol, ready::Int)::Bool

	haskey(dfg.g.variables, label) && (return dfg.g.variables[label].ready == ready)
	haskey(dfg.g.factors, label) && (return dfg.g.factors[label].ready == ready)

	#TODO should this be a breaking error?
	@error "Node not in factor or variable"
	return false
end

function _isbackendset(dfg::LightDFG, label::Symbol, backendset::Int)::Bool
	haskey(dfg.g.variables, label) && (return dfg.g.variables[label].backendset == backendset)
	haskey(dfg.g.factors, label) && (return dfg.g.factors[label].backendset == backendset)

	#TODO should this be a breaking error?
	@error "Node not a factor or variable"
	return false
end

"""
    $(SIGNATURES)
Retrieve a list of labels of the immediate neighbors around a given variable or factor.
"""
function getNeighbors(dfg::LightDFG, node::DFGNode; ready::Union{Nothing, Int}=nothing, backendset::Union{Nothing, Int}=nothing)::Vector{Symbol}
	label = node.label
    if !exists(dfg, label)
        error("Variable/factor with label '$(node.label)' does not exist in the factor graph")
    end

	neighbors_il =  FactorGraphs.outneighbors(dfg.g, dfg.g.labels[label])
	neighbors_ll = [dfg.g.labels[i] for i in neighbors_il]
    # Additional filtering
    ready != nothing && filter!(lbl -> _isready(dfg, lbl, ready), neighbors_ll)
	backendset != nothing && filter!(lbl -> _isbackendset(dfg, lbl, backendset), neighbors_ll)

    # Variable sorting (order is important)
    if typeof(node) <: AbstractDFGFactor
        order = intersect(node._variableOrderSymbols, neighbors_ll)#map(v->v.dfgNode.label, neighbors))
        return order
    end

    return neighbors_ll
end

"""
    $(SIGNATURES)
Retrieve a list of labels of the immediate neighbors around a given variable or factor specified by its label.
"""
function getNeighbors(dfg::LightDFG, label::Symbol; ready::Union{Nothing, Int}=nothing, backendset::Union{Nothing, Int}=nothing)::Vector{Symbol}
	if !exists(dfg, label)
        error("Variable/factor with label '$(label)' does not exist in the factor graph")
    end

	neighbors_il =  FactorGraphs.outneighbors(dfg.g, dfg.g.labels[label])
	neighbors_ll = [dfg.g.labels[i] for i in neighbors_il]
    # Additional filtering
    ready != nothing && filter!(lbl -> _isready(dfg, lbl, ready), neighbors_ll)
	backendset != nothing && filter!(lbl -> _isbackendset(dfg, lbl, backendset), neighbors_ll)

    # Variable sorting (order is important)
    if haskey(dfg.g.factors, label)
        order = intersect(dfg.g.factors[label]._variableOrderSymbols, neighbors_ll)#map(v->v.dfgNode.label, neighbors))
        return order
    end

    return neighbors_ll

end

function _copyIntoGraph!(sourceDFG::LightDFG{<:AbstractParams, V, F}, destDFG::LightDFG{<:AbstractParams, V, F}, ns::Vector{Int}, includeOrphanFactors::Bool=false)::Nothing where {V <: AbstractDFGVariable, F <: AbstractDFGFactor}
	#kan ek die in bulk copy, soos graph en dan nuwe map maak
	# Add all variables first,
	labels = [sourceDFG.g.labels[i] for i in ns]

	for v in values(sourceDFG.g.variables)
	    if v.label in labels
			exists(destDFG, v) ? (@warn "_copyIntoGraph $(v.label) exists, ignoring") : addVariable!(destDFG, deepcopy(v))
	    end
	end

    # And then all factors to the destDFG.
    for f in values(sourceDFG.g.factors)
        if f.label in labels
            # Get the original factor variables (we need them to create it)
			neigh_ints = FactorGraphs.outneighbors(sourceDFG.g, sourceDFG.g.labels[f.label])

			neigh_labels = [sourceDFG.g.labels[i] for i in neigh_ints]
            # Find the labels and associated variables in our new subgraph
            factVariables = V[]
            for v_lab in neigh_labels
                if haskey(destDFG.g.variables, v_lab)
                    push!(factVariables, getVariable(destDFG, v_lab))
                    #otherwise ignore
                end
            end

            # Only if we have all of them should we add it (otherwise strange things may happen on evaluation)
            if includeOrphanFactors || length(factVariables) == length(neigh_labels)
                exists(destDFG, f.label) ? (@warn "_copyIntoGraph $(f.label) exists, ignoring") : addFactor!(destDFG, factVariables, deepcopy(f))
            end
        end
    end
    return nothing
end

"""
    $(SIGNATURES)
Retrieve a deep subgraph copy around a given variable or factor.
Optionally provide a distance to specify the number of edges should be followed.
Optionally provide an existing subgraph addToDFG, the extracted nodes will be copied into this graph. By default a new subgraph will be created.
Note: By default orphaned factors (where the subgraph does not contain all the related variables) are not returned. Set includeOrphanFactors to return the orphans irrespective of whether the subgraph contains all the variables.
"""
function getSubgraphAroundNode(dfg::LightDFG{P,V,F}, node::DFGNode, distance::Int64=1, includeOrphanFactors::Bool=false, addToDFG::LightDFG=LightDFG{P,V,F}())::LightDFG where {P <: AbstractParams, V <: AbstractDFGVariable, F <: AbstractDFGFactor}
    if !exists(dfg,node.label)
        error("Variable/factor with label '$(node.label)' does not exist in the factor graph")
    end

	# Get a list of all unique neighbors inside 'distance'
	ns = neighborhood(dfg.g, dfg.g.labels[node.label], distance)

	# Copy the section of graph we want
	_copyIntoGraph!(dfg, addToDFG, ns, includeOrphanFactors)
	return addToDFG

end
# dfg::LightDFG{P,V,F}
# where {P <: AbstractParams, V <: AbstractDFGVariable, F <: AbstractDFGFactor}

"""
    $(SIGNATURES)
Get a deep subgraph copy from the DFG given a list of variables and factors.
Optionally provide an existing subgraph addToDFG, the extracted nodes will be copied into this graph. By default a new subgraph will be created.
Note: By default orphaned factors (where the subgraph does not contain all the related variables) are not returned. Set includeOrphanFactors to return the orphans irrespective of whether the subgraph contains all the variables.
"""
function getSubgraph(dfg::LightDFG{P,V,F}, variableFactorLabels::Vector{Symbol}, includeOrphanFactors::Bool=false, addToDFG::LightDFG=LightDFG{P,V,F}())::LightDFG where {P <: AbstractParams, V <: AbstractDFGVariable, F <: AbstractDFGFactor}
	for label in variableFactorLabels
		if !exists(dfg, label)
            error("Variable/factor with label '$(label)' does not exist in the factor graph")
        end
    end

	variableFactorInts = [dfg.g.labels[s] for s in variableFactorLabels]

    _copyIntoGraph!(dfg, addToDFG, variableFactorInts, includeOrphanFactors)
    return addToDFG
end

"""
    $(SIGNATURES)
Get an adjacency matrix for the DFG, returned as a Matrix{Union{Nothing, Symbol}}.
Rows are all factors, columns are all variables, and each cell contains either nothing or the symbol of the relating factor.
The first row and first column are factor and variable headings respectively.
"""
function getAdjacencyMatrix(dfg::LightDFG)::Matrix{Union{Nothing, Symbol}}
	#TODO Why does it need to be sorted?
	varLabels = sort(collect(keys(dfg.g.variables)))#ort(map(v->v.label, getVariables(dfg)))
    factLabels = sort(collect(keys(dfg.g.factors)))#sort(map(f->f.label, getFactors(dfg)))
    vDict = Dict(varLabels .=> [1:length(varLabels)...].+1)

    adjMat = Matrix{Union{Nothing, Symbol}}(nothing, length(factLabels)+1, length(varLabels)+1)
    # Set row/col headings
    adjMat[2:end, 1] = factLabels
    adjMat[1, 2:end] = varLabels
    for (fIndex, factLabel) in enumerate(factLabels)
        factVars = getNeighbors(dfg, getFactor(dfg, factLabel))
        map(vLabel -> adjMat[fIndex+1,vDict[vLabel]] = factLabel, factVars)
    end
    return adjMat
end

function getAdjacencyMatrixSparse(dfg::LightDFG)::Tuple{LightGraphs.SparseMatrixCSC, Vector{Symbol}, Vector{Symbol}}
	varLabels = collect(keys(dfg.g.variables))
    factLabels = collect(keys(dfg.g.factors))
	varIndex = [dfg.g.labels[s] for s in varLabels]
	factIndex = [dfg.g.labels[s] for s in factLabels]

	adj = adjacency_matrix(dfg.g)

	adjvf = adj[factIndex, varIndex]
	return adjvf, varLabels, factLabels
end

"""
    $(SIGNATURES)
Gets an empty and unique CloudGraphsDFG derived from an existing DFG.
"""
function _getDuplicatedEmptyDFG(dfg::LightDFG{P,V,F})::LightDFG where {P <: AbstractParams, V <: AbstractDFGVariable, F <: AbstractDFGFactor}
    newDfg = LightDFG{P,V,F}(;
		userId=dfg.userId, robotId=dfg.robotId, sessionId=dfg.sessionId,
		params=deepcopy(dfg.solverParams))
	newDfg.description ="(Copy of) $(dfg.description)"
	return newDfg
end


text:
const PUDL_DESCRIPTION = """
	This archive contains raw input data for the Public Utility Data Liberation (PUDL) software developed by [Catalyst Cooperative](https://catalyst.coop/).
	It is organized into [Frictionless Data Packages](https://specs.frictionlessdata.io/data-package/).
	For additional information about this data and PUDL, see the following resources:
	* [The PUDL Repository on GitHub](https://github.com/catalyst-cooperative/pudl)
	* [PUDL Documentation](https://readthedocs.org/projects/catalystcoop-pudl/)
	* [Other Catalyst Cooperative data archives](https://zenodo.org/communities/catalyst-cooperative/)
	"""

function unpack_if_zipped(filename)
    ext = last(splitext(filename))
    f = ext == ".tgz" || ext == ".zip" ? unpack : identity
    return f(str)
end

"""
    pudl_datasets() -> Vector{String}

List the names of the datasets available as `DataDep`s.

Each dataset is available via `dataset"name"`. See DataDeps.jl for details.
"""
pudl_datasets() = sort(collect(filter(startswith("pudl"), keys(DataDeps.registry))))


text:
#===============================================================================

    RawData.jl

    struct containing 3D array of raw data
    parameterized on unit type

    Dimensions:
    1. Unit: `cycle` for amplification data (type Int_T)
       `temperature` for melting curve data (type Float_T).
    2. Well: enumerated type, maximum 16.
    3. Channel: 1 or 2.

    Author: Tom Price
    Date:   July 2019

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

abstract type AbstractRaw end

struct RawData{F <: Union{Int_T,Float_T}} <: AbstractRaw
    data ::Array{F,3}
end

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
abstract type AbstractProcess <: AbstractEvent end
abstract type DiscreteProcess <: AbstractProcess end

struct InterruptException <: Exception
  by :: AbstractProcess
  cause :: Any
end

struct EmptySchedule <: Exception end

struct EventKey
  time :: Float64
  priority :: Int
  id :: UInt
end

function isless(a::EventKey, b::EventKey) :: Bool
  (a.time < b.time) || (a.time == b.time && a.priority > b.priority) || (a.time == b.time && a.priority == b.priority && a.id < b.id)
end

mutable struct Simulation <: Environment
  time :: Float64
  heap :: DataStructures.PriorityQueue{BaseEvent, EventKey}
  eid :: UInt
  sid :: UInt
  active_proc :: Union{AbstractProcess, Nothing}
  function Simulation(initial_time::Number=zero(Float64))
    new(initial_time, DataStructures.PriorityQueue{BaseEvent, EventKey}(), zero(UInt), zero(UInt), nothing)
  end
end

function step(sim::Simulation)
  isempty(sim.heap) && throw(EmptySchedule())
  (bev, key) = DataStructures.peek(sim.heap)
  DataStructures.dequeue!(sim.heap)
  sim.time = key.time
  bev.state = processed
  for callback in bev.callbacks
    callback()
  end
end

function now(sim::Simulation)
  sim.time
end

function now(ev::AbstractEvent)
  return now(environment(ev))
end

function active_process(sim::Simulation) :: AbstractProcess
  sim.active_proc
end

function reset_active_process(sim::Simulation)
  sim.active_proc = nothing
end

function set_active_process(sim::Simulation, proc::AbstractProcess)
  sim.active_proc = proc
end


text:
@testset "Generic elliptic curve" begin

  @testset "Constructors" begin
    @test_throws ErrorException EllipticCurve([1])
    @test_throws ErrorException EllipticCurve([1, 2, 3])
    @test_throws ErrorException EllipticCurve([1, 2, 3, 4, 5, 6])
    @test_throws ErrorException EllipticCurve([0, 0])
    @test_throws ErrorException EllipticCurve([0, 0, 0, 0, 0])

    E = @inferred EllipticCurve([1, 2], check = false)
    @test typeof(E) == EllCrv{fmpq}
    @test a_invars(E) == (0, 0, 0, 1, 2)
    @test coefficients(E) == (0, 0, 0, 1, 2)

    E = @inferred EllipticCurve([1, 2, 3, 4, 5])
    @test typeof(E) == EllCrv{fmpq}
    @test a_invars(E) == (1, 2, 3, 4, 5)

    # this is Cremona: 11a2, lmfdb: 11.a1
    E11_a1 = @inferred EllipticCurve([0, -1, 1, -7820, -263580], check = false)

    # this is Cremona: 41a1, lmfdb: 41.a1
    E43_a1 = @inferred EllipticCurve([0, 1, 1, 0, 0])

    E = @inferred EllipticCurve([0, 0, 0, 1, 2])
    @test is_short_weierstrass_model(E)

    E = @inferred EllipticCurve([0, 0, 0, 1, 2], check = false)
    @test is_short_weierstrass_model(E)

    Qx, x = PolynomialRing(FlintQQ, "x")

    f1 = x^3+3*x+5
    g1 = x+2
    E = EllipticCurve(f1, g1)
    f2, g2 = hyperelliptic_polynomials(E)
    @test f1 == f2 && g1 == g2

    E = @inferred EllipticCurve(f1)
    @test E isa EllCrv{fmpq}
    E = @inferred EllipticCurve(f1, check = false)
    @test E isa EllCrv{fmpq}

    E = EllipticCurve(f1, 1)
    f2, g2 = hyperelliptic_polynomials(E)
    @test f1 == f2 && 1 == g2

    E = EllipticCurve(x^3 + 1, zero(Qx))
    @test E isa EllCrv{fmpq}

    @test_throws ArgumentError EllipticCurve(x^10-21, x^3+5)
    @test_throws ArgumentError EllipticCurve(x^3+3, x^3+5)
    @test_throws ArgumentError EllipticCurve(3*x^3 + 1)

    K, a = NumberField(x^2 - x - 1, "a")
    OK = maximal_order(K)

    E31_1_a1 = @inferred EllipticCurve([K(1), a + 1, a, a, K(0)])
    @test typeof(E31_1_a1) == EllCrv{nf_elem}

    # lmfdb: 116.1-a1
    E116_1_a1 = @inferred EllipticCurve([K(1), K(-1), a, -a, K(0)] )
    @test typeof(E31_1_a1) == EllCrv{nf_elem}

    # short example
    Eshort = @inferred EllipticCurve([4, 0])
    @test typeof(Eshort) == EllCrv{fmpq}
    @test a_invars(Eshort) == (0, 0, 0, 4, 0)
  end

  QQx, = PolynomialRing(QQ, 5)
  F = FractionField(QQx)
  a1, a2, a3, a4, a6 = F.(gens(QQx))
  E = EllipticCurve([a1, a2, a3, a4, a6])
  @test b_invars(E) == (a1^2 + 4*a2, a1*a3 + 2*a4, a3^2 + 4*a6, a1^2*a6 - a1*a3*a4 + a2*a3^2 + 4*a2*a6 - a4^2)
  b2,b4,b6,b8 = b_invars(E)
  c4, c6 = c_invars(E)
  @test (c4, c6) == c_invars(E) # to test caching
  @test (c4, c6) == (a1^4 + 8*a1^2*a2 - 24*a1*a3 + 16*a2^2 - 48*a4, -a1^6 - 12*a1^4*a2 + 36*a1^3*a3 - 48*a1^2*a2^2 + 72*a1^2*a4 + 144*a1*a2*a3 - 64*a2^3 + 288*a2*a4 - 216*a3^2 - 864*a6)
  @test c4 == b2^2 - 24*b4
  @test c6 == -b2^3 + 36*b2*b4 - 216*b6
  d = discriminant(E)
  @test d == -b2^2*b8 - 8*b4^3 - 27*b6^2 + 9*b2*b4*b6
  @test 1728*d == c4^3 - c6^2

  # base change
  E = EllipticCurve([1, 0, 0, 0, 1])
  K, = quadratic_field(-1)
  EK = @inferred base_change(K, E)
  @test a_invars(EK) == (1, 0, 0, 0, 1)
  @test base_field(EK) === K

  F = GF(2)
  f = x -> divexact(F(numerator(x)), F(denominator(x)))
  EF = base_change(f, E)
  @test a_invars(EF) == (1, 0, 0, 0, 1)
  @test base_field(EF) === F

  # The following curves will be used in later tests
  # Creation of these was tested in previous testset
  E11_a1 = EllipticCurve([0, -1, 1, -7820, -263580], check =false)

  E43_a1 = EllipticCurve([0, 1, 1, 0, 0])

  Qx, x = PolynomialRing(FlintQQ, "x")
  K, a = NumberField(x^2 - x - 1, "a")
  OK = maximal_order(K)

  E31_1_a1 = EllipticCurve([K(1), a + 1, a, a, K(0)])

  E116_1_a1 =EllipticCurve([K(1), K(-1), a, -a, K(0)] )

  Eshort = EllipticCurve([4, 0])

  @testset "Field access" begin
    @test base_field(E11_a1) == FlintQQ
    @test base_field(E43_a1) == FlintQQ
    @test base_field(E31_1_a1) == K
    @test base_field(E116_1_a1) == K
    @test base_field(Eshort) == FlintQQ
  end

  @testset "Point construction" begin
    P = @inferred E43_a1([FlintQQ(-1), FlintQQ(0)])
    @test typeof(P) == EllCrvPt{fmpq}
    @test parent(P) == E43_a1
    @test @inferred is_finite(P)
    @test @inferred !is_infinite(P)

    P = @inferred E43_a1([-1, 0], check = false)
    @test typeof(P) == EllCrvPt{fmpq}
    @test parent(P) == E43_a1
    @test @inferred is_finite(P)
    @test @inferred !is_infinite(P)

    P = @inferred E43_a1([fmpz(-1), fmpz(0)])
    @test typeof(P) == EllCrvPt{fmpq}
    @test parent(P) == E43_a1
    @test @inferred is_finite(P)
    @test @inferred !is_infinite(P)

# the error is/was from doing QQ(K(0)) - which is possible now
#    @test_throws MethodError E43_a1([gen(K), gen(K)])

    @test_throws ErrorException E43_a1([2, 2])

    P = @inferred infinity(E43_a1)
    @test @inferred !is_finite(P)
    @test @inferred is_infinite(P)

    P = @inferred E43_a1([FlintQQ(-1), FlintQQ(0)])

    P = @inferred E116_1_a1([K(0), -K(a)])
    @test typeof(P) == EllCrvPt{nf_elem}
    @test parent(P) == E116_1_a1
    @test @inferred is_finite(P)
    @test @inferred !is_infinite(P)

    P = @inferred infinity(E116_1_a1)
    @test !is_finite(P)
    @test is_infinite(P)

    @test_throws ErrorException E116_1_a1([1, 1])

    P = @inferred Eshort([2, 4], check = false)
    @test @inferred is_finite(P)
    @test typeof(P) == EllCrvPt{fmpq}
    @test parent(P) == Eshort

    E = EllipticCurve(GF(7,2),[1,2,3,4,5])
    L = @inferred points_with_x(E,0)
    @test E([0,5]) in L && E([0, 6]) in L

  end

  @testset "Equation" begin
    E = EllipticCurve( [1, 2, 3, 4, 5])
    Kxy, (x,y) = PolynomialRing(base_field(E), ["x","y"])
    @test y^2 + x*y + 3*y - x^3 - 2*x^2 - 4*x - 5 == @inferred Kxy(equation(E))
  end

  @testset "Discriminant" begin
    @test (2*a + 10)*OK == @inferred (discriminant(E116_1_a1)*OK)
    @test -43 == @inferred discriminant(E43_a1)
    @test -4096 == @inferred discriminant(Eshort)
    E43_a1 = EllipticCurve([0, 1, 1, 0, 0], check = false)
    @test discriminant(E43_a1) == -43
    Eshort = @inferred EllipticCurve([4, 0], check = false)
    @test discriminant(Eshort) == -4096
  end

  @testset "j-invariant" begin
    b = (fmpq(-215055, 58) * a - fmpq(65799, 29))
    @test  b == @inferred j_invariant(E116_1_a1)
    @test fmpq(-4096, 43) == @inferred j_invariant(E43_a1)
    @test 1728 == @inferred j_invariant(Eshort)

    E = @inferred elliptic_curve_from_j_invariant(23)
    @test j_invariant(E) == 23

    E = @inferred elliptic_curve_from_j_invariant(1728)
    @test j_invariant(E) == 1728

    E = @inferred elliptic_curve_from_j_invariant(GF(3,2)(0))
    @test j_invariant(E) == 0

    E = @inferred elliptic_curve_from_j_invariant(GF(2,4)(0))
    @test j_invariant(E) == 0

  end

  @testset "Point aritmetic" begin
    #addition
    P = @inferred E43_a1([FlintQQ(-1), FlintQQ(0)])
    O = infinity(E43_a1)

    @test P == @inferred P + O
    @test E43_a1([2, -4]) == P + P
    @test E43_a1([fmpq(-2, 9), fmpq(1, 27)]) == P + P + P

    P = @inferred E116_1_a1([K(0), -K(a)])
    @test E116_1_a1([1, -1]) == @inferred P + P
    @test E116_1_a1([K(1), -a]) == @inferred P + P + P
    @test infinity(E116_1_a1) == @inferred P + P + P + P + P
    @test P == @inferred P + infinity(E116_1_a1)

    P = @inferred Eshort([2, 4])
    O = @inferred infinity(Eshort)
    @test_throws ErrorException Eshort([-1, -1])

    @test Eshort([0, 0]) == @inferred P + P
    @test P == @inferred O + P

    P = Eshort([2, 4])
    @test Eshort([2, -4]) == @inferred -P
    P = infinity(Eshort)
    @test P == @inferred -P

    P = E43_a1([-1, 0])
    @test E43_a1([-1, -1]) == @inferred -P
    P = infinity(E43_a1)
    @test P == @inferred -P

    # inversion
    P = @inferred E116_1_a1([K(0), -K(a)])
    @test E116_1_a1([0, 0]) == @inferred -P
    P = infinity(E116_1_a1)
    @test P == @inferred -P

    # equality
    P1 = Eshort([2, 4])
    @test @inferred ==(P1, P1)
    P2 = infinity(Eshort)
    @test @inferred ==(P2, P2)
    @test @inferred !==(P2, P1)

    P1 = E43_a1([-1, 0])
    @test @inferred ==(P1, P1)
    P2 = infinity(E43_a1)
    @test @inferred ==(P2, P2)
    @test @inferred !==(P2, P1)

    P1 = E116_1_a1([K(0), -K(a)])
    @test @inferred ==(P1, P1)
    P2 = infinity(E116_1_a1)
    @test @inferred ==(P2, P2)
    @test @inferred !==(P2, P1)

    # scalar multiplication
    P1 = Eshort([2, 4])
    @test Eshort([0, 0]) == @inferred 2*P1
    @test infinity(Eshort) == @inferred 4*P1

    P1 = E43_a1([-1, 0])
    @test E43_a1([fmpq(11, 49), fmpq(-363, 343)]) == @inferred 4*P1

    P1 = E116_1_a1([K(0), -K(a)])
    @test E116_1_a1([K(0), K(0)]) == @inferred 4*P1
    @test infinity(E116_1_a1) == @inferred 5*P1

    #division

    P1 = Eshort([2, 4])
    Q = (2*P1)//2
    @test Q == P1 || Q == -P1

    P1 = E116_1_a1([K(0), -K(a)])
    @test (3*P1)//3 == P1
    @test @inferred 3*(3*P1)//(-3) == -3*P1
    @test_throws ErrorException P1//5
  end
end


text:
###############################################################################
#
# Photosynthesis model parameter set -- temperature dependencies and etc
#
###############################################################################
#= AbstractPhotoModelParaSet type tree
---> C3ParaSet
---> C4ParaSet
=#
"""
    abstract type AbstractPhotoModelParaSet{FT}

Hierarchy of the `AbstractPhotoModelParaSet`:
- [`C3ParaSet`](@ref)
- [`C4ParaSet`](@ref)
"""
abstract type AbstractPhotoModelParaSet{FT} end

"""
    mutable struct C3Paraset{FT}

Parameter sets for C3 photosynthesis.

# Fields
$(DocStringExtensions.FIELDS)
"""
mutable struct C3ParaSet{FT<:AbstractFloat} <: AbstractPhotoModelParaSet{FT}
    "Jmax temperature dependency"
    JT ::AbstractTDParameterSet{FT}
    "Kc temperature dependency"
    KcT::AbstractTDParameterSet{FT}
    "Ko temperature dependency"
    KoT::AbstractTDParameterSet{FT}
    "Respiration temperature dependency"
    ReT::AbstractTDParameterSet{FT}
    "Vcmax temperature dependency"
    VcT::AbstractTDParameterSet{FT}
    "Γ_star temperature dependency"
    ΓsT::AbstractTDParameterSet{FT}
    "Fluorescence model"
    Flu::AbstractFluoModelParaSet{FT}
    "Vcmax25 and respiration correlation"
    VR   ::FT
    "Coefficient 4.0/4.5 for NADPH/ATP requirement stochiometry, respectively"
    Eff_1::FT
    "Coefficient 8.0/10.5 for NADPH/ATP requirement stochiometry, respectively"
    Eff_2::FT
end

"""
    mutable struct C4ParaSet{FT}

Parameter sets for C3 photosynthesis.

# Fields
$(DocStringExtensions.FIELDS)
"""
mutable struct C4ParaSet{FT<:AbstractFloat} <: AbstractPhotoModelParaSet{FT}
    "Kpep temperature dependency"
    KpT::AbstractTDParameterSet{FT}
    "Respiration temperature dependency"
    ReT::AbstractTDParameterSet{FT}
    "Vcmax temperature dependency"
    VcT::AbstractTDParameterSet{FT}
    "Vpmax temperature dependency"
    VpT::AbstractTDParameterSet{FT}
    "Fluorescence model"
    Flu::AbstractFluoModelParaSet{FT}
    "Vcmax25 and respiration correlation"
    VR::FT
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land