text:
module Args
    using ArgParse
    using ..Util

    export read_args, Timespan, VarMapping, SamplingTuple, DomainSize

    VarMapping      = Dict{String, String}
    SamplingTuple   = NTuple{3, Float64}
    DomainSize      = NTuple{2, NTuple{2, Float64}} # ((xmin,xmax), (ymin,ymax))

    struct Timespan{T <: Number}
        t_start :: T
        t_end   :: T
    end

    function ArgParse.parse_item(::Type{Timespan{T}}, x::AbstractString) where T <: Number
        x = strip(x)

        t_start :: T = zero(T)
        t_end   :: T = typemax(T)

        if x == "all"
            return Timespan{T}(t_start, t_end)
        end

        if occursin(',', x)
            timestamps = split(x, ',')

            if length(timestamps) != 2
                throw(ArgumentError(""""$x" does not match the format "start,end"!"""))
            end

            try
                timestamps[1] = strip(timestamps[1])
                if length(timestamps[1]) != 0
                    t_start = parse(T, strip(timestamps[1]))
                end
            catch
                throw(ArgumentError(""""$(timestamps[1])" is not a floating point value!"""))
            end

            try
                timestamps[2] = strip(timestamps[2])
                if length(timestamps[2]) != 0
                    t_end = parse(T, timestamps[2])
                end
            catch
                throw(ArgumentError(""""$(timestamps[2])" is not a floating point value!"""))
            end

            if (t_start < zero(T))
                throw(ArgumentError("Timestamps have to be non-negative!"))
            end

            if t_start > t_end
                throw(ArgumentError("The start timestamp must be less than the end timestamp!"))
            end

            return Timespan{T}(t_start, t_end)
        else
            try
                t = parse(T, x)

                if t < zero(T)
                    throw(ArgumentError(""""$t" must be non-negative!"""))
                end

                return Timespan{T}(t, t)
            catch
                throw(ArgumentError(""""$t" is not a floating point value!"""))
            end
        end
    end

    function ArgParse.parse_item(::Type{VarMapping}, x::AbstractString)
        entries = Set(split(x, ','))
        non_mapping_entries = filter(entry -> !occursin("=>", entry), entries)
        mapping_entries = setdiff(entries, non_mapping_entries)
        mapping_entries = map(str -> split(str, "=>"), collect(mapping_entries))

        entries = merge(Dict(map(str_list   -> String(str_list[1])  => String(str_list[2]), mapping_entries)),
                        Dict(map(entry      -> String(entry)        => String(entry),       collect(non_mapping_entries))))
        return entries
    end

    function ArgParse.parse_item(::Type{SamplingTuple}, x::AbstractString)
        components = split(x, ',')

        if (length(components) ∉ [1, 3])
            throw(ArgumentError(""""$x" has $(length(components)) components but can only have 1 or 3!"""))
        end

        if length(components) == 1
            components = [components[1], components[1], components[1]]
        end

        try
            return tuple(map(x -> parse(Float64, x), components)...)
        catch
            throw(ArgumentError("""The components of "$x" cannot be parsed as numbers!"""))
        end
    end

    function ArgParse.parse_item(::Type{DomainSize}, x::AbstractString)
        floats = split(x, ',')

        if length(floats) != 4
            throw(ArgumentError(""""$x" is not in the format "xmin,xmax,ymin,ymax"!"""))
        end

        try
            floats = tuple(map(x -> parse(Float64, x), floats)...)
            return ((floats[1], floats[2]), (floats[3], floats[4]))
        catch
            throw(ArgumentError("""The components of "$x" cannot be parsed as numbers!"""))
        end
    end

    function read_args()
        parser_settings = ArgParseSettings()
        @add_arg_table! parser_settings begin
            "--output-file", "-o"
                help = "The filename (without extension) to use for the output."
                default = ""
            "--output-time", "-t"
                help = """The time range of the output. Format: "[start,]end" or "all".\n
                                \t start and end are interpreted as timeSTAMPs, not timeSTEP indices!
                                   The end timestep is bounded by the last input timestep.\n
                                \t If only end is given, only that timestep will be rasterized.\n
                                \t If start or end is an empty string, it will be assumed to be 0 or infinity respectively.\n
                                \t By default, every timestep will be output."""
                arg_type = Timespan{Float64}
            "--output-steps", "-s"
                help = """The time range of the output. Format: "[start,]end" or "all".\n
                                \t start and end are interpreted as timeSTEP indices, not timeSTAMPs!
                                   The end timestep is bounded by the last input timestep.\n
                                \t If only end is given, only that timestep will be rasterized.\n
                                \t If start or end is an empty string, it will be assumed to be 0 or infinity respectively.\n
                                \t By default, every timestep will be output."""
                arg_type = Timespan{UInt64}
            "--domain", "-d"
                help = """Specify a custom domain which should be rasterized. Format: "xmin,xmax,ymin,ymax".\n
                                \t By default, the whole domain will be rasterized. The specified domain will be clipped to the actual domain."""
                arg_type = DomainSize
                default = ((-Inf, Inf), (-Inf, Inf))
            "--sampling-rate", "-r"
                help = """The size in meters of one cell edge in the resampled output.\n
                                \t Format: dx[,dy,dz]. If only dx is given, dy and dz will be set equal to dx.\n
                                \t Examples: 100; 50,70,80"""
                arg_type = SamplingTuple
                default = (100., 100., 100.)
            "--water-height"
                help = """Some bathymetry grids place the seafloor at z = 0. Set the water height such that it matches the sea surface's z-coordinate."""
                arg_type = Float64
                default = 0.
            "--memory-limit", "-m"
                help = """The maximum amount of RAM the script should use. This is only a SOFT limit!\n
                                \t [CAUTION] Currently, this limit is quite imprecise, do not specify more than ~half of the available RAM here!\n
                                \t Examples: 8G; 2T; 512M\n
                                \t IEC-Prefixes are used: 1K = 1KiB = 1024B, ..."""
                default = "8G"
            "--tanioka"
                help = """Apply Tanioka's method for converting horizontal displacements to vertical ones during rasterization.\n
                                \t If the variables U or V are not found in --seafloor-vars, they will be added automatically."""
                action = :store_true
            "--seafloor-vars"
                help = """The names of all variables that should be extracted from the 2D grid at the seafloor, separated by commas.\n
                                \t Note that the bathymetry "b" will always be extracted from geometry. Use "b=>NewName" to rename
                                   the output variable.\n
                                \t [CAUTION] You cannot have the same output name for multiple variables,
                                   even across --seafloor-vars, --surface-vars and --volumetric-vars!\n
                                \t Format: [M1[,M2[,...]] where Mi is either a variable name or a mapping like "W=>d".\n
                                \t Examples: U,V,W; U,V,W=>eta; U=>X,V=>Y"""
                arg_type = VarMapping
                default=Dict("W"=>"d", "b"=>"b")
            "--surface-vars"
                help = """The names of all variables that should be extracted from the 2D grid at the sea surface, separated by commas.\n
                                \t The same formatting rules apply as for --seafloor-vars."""
                arg_type = VarMapping
                default=Dict("W"=>"eta")
            "--volumetric-vars"
                help = """The names of all variables that should be extracted from the 3D grid, separated by commas.\n
                                \t The same formatting rules apply as for --seafloor-vars."""
                arg_type = VarMapping
                default=Dict("u"=>"u", "v"=>"v")
            "--seafloor-only"
                help = """Only rasterize the seafloor displacements. The outputs can be processed by kajiura.jl."""
                action = :store_true
            "--load-balancer"
                help = """The static load balancer to choose when assigning simplices to threads. Options are:\n
                                \t naive: Every thread works on a bin that is |D| / n_threads in size\n
                                \t [DO NOT USE] count: Every thread works on a bin that approximately contains n_simplices / n_threads simplices\n
                                \t [DO NOT USE] workload: Every thread has roughly the same ESTIMATED (!) workload (read: processing time) in its bin. Run with --lb-autotune first.
                                """
                default="naive"
            "--lb-autotune"
                help = """[DO NOT USE] Perform a test run with the naive load balancer and collect data to autotune the workload-LB.
                                The found parameters will be saved automatically and the workload-LB can then be used."""
                action = :store_true
            "input-file-2d"
                help = "The SeisSol 2D output in XDMF format. Use the output file with the _surface suffix."
                required = true
            "input-file-3d"
                help = "The SeisSol 3D output in XDMF format. Use the output file without the _surface suffix."
                default = ""
        end

        args = parse_args(parser_settings)
        validate_args!(args)

        return args
    end

    function validate_args!(args::Dict)
        args["memory-limit"] = parse_size(args["memory-limit"])

        if !isnothing(args["output-time"]) && !isnothing(args["output-steps"])
           throw(ArgumentError("Setting both '--output-time' and '--output-steps' simultaneously is not allowed!"))
        end

        if !(args["load-balancer"] in ["naive", "count", "workload"])
            throw(ArgumentError("Load Balancer '$(args["load-balancer"])' does not exist!"))
        end

        return args
    end
end

max_stars_repo_name: SeisSol/SAMPLER

max_issues_repo_name: SeisSol/SAMPLER

max_forks_repo_name: SeisSol/SAMPLER


text:
using .Samsara

#===========================================================================#
# Etabler eit DemoSys : Demonstrasjon om korleis eit Env-system kan lages:
# - subtype AbstractSystem
# - Definer step_system_mechanics!(arg::NY_TYPE)
# Når du seinare skal gjøre det samme, kan dette defineres i ei eiga fil,
# som includeres i prosjektet ditt
#
# Merk at tester for DemoSys er å finne i test_Samsara
#===========================================================================#
mutable struct DemoSys <: AbstractSystem
    _observable_parameters
    _latent_variables
    function DemoSys(observable_param =(1,) , latent_param =nothing)
        latent_param == nothing && (latent_param = 10)
        # Lag eit tomt struct, oppdater med simulation-function (som er overskrevet fra utsida,
        #   og definerer system mechanics), og returner oppdatert struct retval.
        new( observable_param, latent_param )
    end
end

""" Dette er måten man definerer system mechanics: Man definerer step_system_mechanics!(arg)
der arg er av type som definerer kva system det er snakk om. DemoSys har tidsdynamikk som øker
lineært for sys._observable_parameters og minker for sys._latent_variables
ARG: sys: stuct med info nok for å tilfredstille Markov property.
RETURN: sys._observable_param -- observerbar state.
"""
function step_system_mechanics!(arg::DemoSys)
    arg._observable_parameters = arg._observable_parameters .+ arg._latent_variables
    arg._latent_variables -= 1
    return arg._observable_parameters
end

function set_action_in!(sys::DemoSys, action)
    @warn "not implemented in DemoSys"
end

ext: jl

max_stars_repo_name: leikanger/Samsara

max_issues_repo_name: leikanger/Samsara

max_forks_repo_name: leikanger/Samsara


text:
using Statistics
using JLD2
using Printf
using Plots
using Oceananigans
using Oceananigans.Units

using Oceananigans.Fields: interpolate, Field
using Oceananigans.Advection: VelocityStencil
using Oceananigans.Architectures: arch_array
using Oceananigans.Coriolis: HydrostaticSphericalCoriolis
using Oceananigans.BoundaryConditions
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom, inactive_node, peripheral_node
using CUDA: @allowscalar, device!
using Oceananigans.Operators
using Oceananigans.Operators: Δzᵃᵃᶜ
using Oceananigans: prognostic_fields

@inline function visualize(field, lev, dims)
    (dims == 1) && (idx = (lev, :, :))
    (dims == 2) && (idx = (:, lev, :))
    (dims == 3) && (idx = (:, :, lev))

    r = deepcopy(Array(interior(field)))[idx...]
    r[ r.==0 ] .= NaN
    return r
end

#####
##### Grid
#####

arch = GPU()
reference_density = 1029

latitude = (-75, 75)

# 0.25 degree resolution
Nx = 1440
Ny = 600
Nz = 48

const Nyears  = 1
const Nmonths = 12
const thirty_days = 30days

output_prefix = "near_global_lat_lon_$(Nx)_$(Ny)_$(Nz)_fine"
pickup_file   = false

#####
##### Load forcing files and inital conditions from ECCO version 4
##### https://ecco.jpl.nasa.gov/drive/files
##### Bathymetry is interpolated from ETOPO1 https://www.ngdc.noaa.gov/mgg/global/
#####

using DataDeps

path = "https://github.com/CliMA/OceananigansArtifacts.jl/raw/ss/new_hydrostatic_data_after_cleared_bugs/quarter_degree_near_global_input_data/"

datanames = ["z_faces-50-levels",
             "bathymetry-1440x600",
             "temp-1440x600-latitude-75",
             "salt-1440x600-latitude-75",
             "tau_x-1440x600-latitude-75",
             "tau_y-1440x600-latitude-75",
             "initial_conditions"]

dh = DataDep("quarter_degree_near_global_lat_lon",
    "Forcing data for global latitude longitude simulation",
    [path * data * ".jld2" for data in datanames]
)

DataDeps.register(dh)

datadep"quarter_degree_near_global_lat_lon"

files = [:file_z_faces, :file_bathymetry, :file_temp, :file_salt, :file_tau_x, :file_tau_y, :file_init]
for (data, file) in zip(datanames, files)
    datadep_path = @datadep_str "quarter_degree_near_global_lat_lon/" * data * ".jld2"
    @eval $file = jldopen($datadep_path)
end

bathymetry = file_bathymetry["bathymetry"]

τˣ = zeros(Nx, Ny, Nmonths)
τʸ = zeros(Nx, Ny, Nmonths)
T★ = zeros(Nx, Ny, Nmonths)
S★ = zeros(Nx, Ny, Nmonths)

# Files contain 1 year (1992) of 12 monthly averages
τˣ = file_tau_x["field"] ./ reference_density
τʸ = file_tau_y["field"] ./ reference_density
T★ = file_temp["field"]
S★ = file_salt["field"]

# Remember the convention!! On the surface a negative flux increases a positive decreases
bathymetry = arch_array(arch, bathymetry)

# Stretched faces taken from ECCO Version 4 (50 levels in the vertical)
z_faces = file_z_faces["z_faces"][3:end]

# A spherical domain
@show underlying_grid = LatitudeLongitudeGrid(arch,
                                              size = (Nx, Ny, Nz),
                                              longitude = (-180, 180),
                                              latitude = latitude,
                                              halo = (3, 3, 3),
                                              z = z_faces,
                                              precompute_metrics = true)

grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bathymetry))

τˣ = arch_array(arch, - τˣ)
τʸ = arch_array(arch, - τʸ)

target_sea_surface_temperature = T★ = arch_array(arch, T★)
target_sea_surface_salinity    = S★ = arch_array(arch, S★)

#####
##### Physics and model setup
#####

νz = 5e-3
κz = 1e-4

using Oceananigans.Operators: Δx, Δy
using Oceananigans.TurbulenceClosures

@inline νhb(i, j, k, grid, lx, ly, lz) = (1 / (1 / Δx(i, j, k, grid, lx, ly, lz)^2 + 1 / Δy(i, j, k, grid, lx, ly, lz)^2 ))^2 / 5days

vertical_diffusivity   = VerticalScalarDiffusivity(VerticallyImplicitTimeDiscretization(), ν=νz, κ=κz)
convective_adjustment  = ConvectiveAdjustmentVerticalDiffusivity(VerticallyImplicitTimeDiscretization(), convective_κz = 1.0)
biharmonic_viscosity   = HorizontalDivergenceScalarBiharmonicDiffusivity(ν=νhb, discrete_form=true)

closures = (vertical_diffusivity, convective_adjustment, biharmonic_viscosity)

#####
##### Boundary conditions / time-dependent fluxes
#####

@inline current_time_index(time, tot_months)     = mod(unsafe_trunc(Int32, time / thirty_days), tot_months) + 1
@inline next_time_index(time, tot_months)        = mod(unsafe_trunc(Int32, time / thirty_days) + 1, tot_months) + 1
@inline cyclic_interpolate(u₁::Number, u₂, time) = u₁ + mod(time / thirty_days, 1) * (u₂ - u₁)

Δz_top    = @allowscalar Δzᵃᵃᶜ(1, 1, grid.Nz, grid.underlying_grid)
Δz_bottom = @allowscalar Δzᵃᵃᶜ(1, 1, 1, grid.underlying_grid)

@inline function surface_wind_stress(i, j, grid, clock, fields, τ)
    time = clock.time
    n₁ = current_time_index(time, Nmonths)
    n₂ = next_time_index(time, Nmonths)

    @inbounds begin
        τ₁ = τ[i, j, n₁]
        τ₂ = τ[i, j, n₂]
    end

    return cyclic_interpolate(τ₁, τ₂, time)
end

u_wind_stress_bc = FluxBoundaryCondition(surface_wind_stress, discrete_form = true, parameters = τˣ)
v_wind_stress_bc = FluxBoundaryCondition(surface_wind_stress, discrete_form = true, parameters = τʸ)

# Linear bottom drag:
μ = 0.001 # ms⁻¹

@inline u_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.u[i, j, 1]
@inline v_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.v[i, j, 1]
@inline u_immersed_bottom_drag(i, j, k, grid, clock, fields, μ) = @inbounds - μ * fields.u[i, j, k]
@inline v_immersed_bottom_drag(i, j, k, grid, clock, fields, μ) = @inbounds - μ * fields.v[i, j, k]

u_immersed_bc = ImmersedBoundaryCondition(bottom = FluxBoundaryCondition(u_immersed_bottom_drag, discrete_form = true, parameters = μ))
v_immersed_bc = ImmersedBoundaryCondition(bottom = FluxBoundaryCondition(v_immersed_bottom_drag, discrete_form = true, parameters = μ))

u_bottom_drag_bc = FluxBoundaryCondition(u_bottom_drag, discrete_form = true, parameters = μ)
v_bottom_drag_bc = FluxBoundaryCondition(v_bottom_drag, discrete_form = true, parameters = μ)

@inline function surface_temperature_relaxation(i, j, grid, clock, fields, p)
    time = clock.time

    n₁ = current_time_index(time, Nmonths)
    n₂ = next_time_index(time, Nmonths)

    @inbounds begin
        T★₁ = p.T★[i, j, n₁]
        T★₂ = p.T★[i, j, n₂]
        T_surface = fields.T[i, j, grid.Nz]
    end

    T★ = cyclic_interpolate(T★₁, T★₂, time)

    return p.λ * (T_surface - T★)
end

@inline function surface_salinity_relaxation(i, j, grid, clock, fields, p)
    time = clock.time

    n₁ = current_time_index(time, Nmonths)
    n₂ = next_time_index(time, Nmonths)

    @inbounds begin
        S★₁ = p.S★[i, j, n₁]
        S★₂ = p.S★[i, j, n₂]
        S_surface = fields.S[i, j, grid.Nz]
    end

    S★ = cyclic_interpolate(S★₁, S★₂, time)

    return p.λ * (S_surface - S★)
end

T_surface_relaxation_bc = FluxBoundaryCondition(surface_temperature_relaxation,
                                                discrete_form = true,
                                                parameters = (λ = Δz_top/7days, T★ = target_sea_surface_temperature))

S_surface_relaxation_bc = FluxBoundaryCondition(surface_salinity_relaxation,
                                                discrete_form = true,
                                                parameters = (λ = Δz_top/7days, S★ = target_sea_surface_salinity))

u_bcs = FieldBoundaryConditions(top = u_wind_stress_bc, bottom = u_bottom_drag_bc, immersed = u_immersed_bc)
v_bcs = FieldBoundaryConditions(top = v_wind_stress_bc, bottom = v_bottom_drag_bc, immersed = v_immersed_bc)
T_bcs = FieldBoundaryConditions(top = T_surface_relaxation_bc)
S_bcs = FieldBoundaryConditions(top = S_surface_relaxation_bc)

free_surface = ImplicitFreeSurface(solver_method=:HeptadiagonalIterativeSolver)

buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState())

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    free_surface = free_surface,
                                    momentum_advection = WENO5(vector_invariant = VelocityStencil()),
                                    coriolis = HydrostaticSphericalCoriolis(),
                                    buoyancy = buoyancy,
                                    tracers = (:T, :S),
                                    closure = (horizontal_diffusivity, vertical_diffusivity, convective_adjustment, biharmonic_viscosity),
                                    boundary_conditions = (u=u_bcs, v=v_bcs, T=T_bcs, S=S_bcs),
                                    tracer_advection = WENO5(underlying_grid))

#####
##### Initial condition:
#####

u, v, w = model.velocities
η = model.free_surface.η
T = model.tracers.T
S = model.tracers.S

@info "Reading initial conditions"
T_init = file_init["T"]
S_init = file_init["S"]

set!(model, T=T_init, S=S_init)
fill_halo_regions!(T)
fill_halo_regions!(S)

@info "model initialized"

#####
##### Simulation setup
#####

Δt = 6minutes  # for initialization, then we can go up to 6 minutes?

simulation = Simulation(model, Δt = Δt, stop_time = Nyears*years)

start_time = [time_ns()]

using Oceananigans.Utils

function progress(sim)
    wall_time = (time_ns() - start_time[1]) * 1e-9

    u = sim.model.velocities.u
    η = sim.model.free_surface.η

    @info @sprintf("Time: % 12s, iteration: %d, max(|u|): %.2e ms⁻¹, wall time: %s",
                    prettytime(sim.model.clock.time),
                    sim.model.clock.iteration, maximum(abs, u),
                    prettytime(wall_time))

    start_time[1] = time_ns()

    return nothing
end

simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))

u, v, w = model.velocities
T = model.tracers.T
S = model.tracers.S
η = model.free_surface.η

output_fields = (; u, v, T, S, η)
save_interval = 5days

simulation.output_writers[:surface_fields] = JLD2OutputWriter(model, (; u, v, T, S, η),
                                                              schedule = TimeInterval(save_interval),
                                                              filename = output_prefix * "_surface",
                                                              indices = (:, :, grid.Nz),
                                                              overwrite_existing = true)

simulation.output_writers[:checkpointer] = Checkpointer(model,
                                                        schedule = TimeInterval(1year),
                                                        prefix = output_prefix * "_checkpoint",
                                                        overwrite_existing = true)

# Let's goo!
@info "Running with Δt = $(prettytime(simulation.Δt))"

run!(simulation, pickup = pickup_file)

@info """
    Simulation took $(prettytime(simulation.run_wall_time))
    Free surface: $(typeof(model.free_surface).name.wrapper)
    Time step: $(prettytime(Δt))
"""


text:
#=
Illustrates vectors in the plane.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 07/09/2014
=#

using PyPlot

#=
    Create pyplot figure and axis. Move left and bottom spines to
    intersect at origin. Remove right and top spines. Return the axis
=#

function move_spines()
    fix, ax = subplots()
    for spine in ["left", "bottom"]
        ax[:spines][spine][:set_position]("zero")
    end

    for spine in ["right", "top"]
        ax[:spines][spine][:set_color]("none")
    end
    return ax
end

function plane_fig()  # illustrate vectors in a plane
    ax = move_spines()

    ax[:set_xlim](-5, 5)
    ax[:set_ylim](-5, 5)
    ax[:grid]()
    vecs = {[2, 4], [-3, 3], [-4, -3.5]}
    for v in vecs
        ax[:annotate](" ", xy=v, xytext=[0, 0],
                    arrowprops={"facecolor"=>"blue",
                                "shrink"=>0,
                                "alpha"=>0.7,
                                "width"=>0.5})
        ax[:text](1.1 * v[1], 1.1 * v[2], string(v))
    end
end

function scalar_multiply()  # illustrate scalar multiplication
    ax = move_spines()
    ax[:set_xlim](-5, 5)
    ax[:set_ylim](-5, 5)

    x = [2, 2]
    ax[:annotate](" ", xy=x, xytext=[0, 0],
                  arrowprops={"facecolor"=>"blue",
                              "shrink"=>0,
                              "alpha"=>1,
                              "width"=>0.5})

    ax[:text](x[1] + 0.4, x[2] - 0.2, L"$x$", fontsize="16")

    scalars = [-2, 2]

    for s in scalars
        v = s .* x
        ax[:annotate](" ", xy=v, xytext=[0, 0],
                      arrowprops={"facecolor"=>"red",
                                  "shrink"=>0,
                                  "alpha"=>0.5,
                                  "width"=>0.5})

        ax[:text](v[1] + 0.4, v[2] - 0.2, LaTeXString("\$$s x\$"),
                  fontsize="16")
    end
end

# Plot the first figure --- three vectors in the plane

plane_fig()

max_stars_repo_name: gsulis/structural

max_issues_repo_name: gsulis/structural

max_forks_repo_name: gsulis/structural


text:
using DynamicGrids, Dispersal, Test, Unitful, Dates

using DynamicGrids: applyrule, SimData, extent
using Unitful: d

@testset "exponential growth" begin
    init =  [1.0 4.0 7.0;
             2.0 5.0 8.0;
             3.0 6.0 9.0]

    test1 = [ 1.0  4.0  7.0;
              2.0  5.0  8.0;
              3.0  6.0  9.0]
    test2 = [ 2.0  8.0 14.0;
              4.0 10.0 16.0;
              6.0 12.0 18.0]
    test3 = [ 4.0 16.0 28.0;
              8.0 20.0 32.0;
             12.0 24.0 36.0]

    output = ArrayOutput(init; tspan=1:3)
    rule = Ruleset(ExponentialGrowth(rate=log(2.0), timestep=1))
    sim!(output, rule)

    @test output[1] == test1
    @test output[2] == test2
    @test output[3] == test3

    output = ArrayOutput(init; tspan=1d:1d:3d)
    rule = Ruleset(ExponentialGrowth(rate=log(2.0), timestep=1d); timestep=1d)
    sim!(output, rule)

    @test output[1] == test1
    @test output[2] == test2
    @test output[3] == test3

    output = ArrayOutput(init; tspan=DateTime(2001,1,1):Day(5):DateTime(2001,1,15))
    rule = Ruleset(ExponentialGrowth(rate=log(2.0)/5, timestep=Day(1)); timestep=Day(5))
    sim!(output, rule)

    @test output[1] == test1
    @test output[2] == test2
    @test output[3] == test3
end

@testset "exponential growth with rate from growth map" begin
    init = Float32[1.0 1.0 1.0;
                   1.0 1.0 1.0;
                   1.0 1.0 1.0]

    suit =  log.(Float32[1.0 1.0 2.0;
                         2.0 1.0 0.5;
                         1.0 1.0 0.5])

    output = ArrayOutput(init; tspan=1:3, aux=(suit=suit,))
    output.extent
    rule = Ruleset(ExponentialGrowth(rate=Aux(:suit), timestep=1, nsteps_type=Float32))
    sim!(output, rule)

    @test output[1] == Float32[1.0 1.0 1.0;
                               1.0 1.0 1.0;
                               1.0 1.0 1.0]
    @test output[2] ≈  Float32[1.0 1.0 2.0;
                               2.0 1.0 0.5;
                               1.0 1.0 0.5]
    @test output[3] ≈  Float32[1.0 1.0 4.0;
                               4.0 1.0 0.25;
                               1.0 1.0 0.25]

    @test DynamicGrids.normalise.(output[3], 0.25f0, 4) == Float32[0.2 0.2 1.0;
                                                                   1.0 0.2 0.0;
                                                                   0.2 0.2 0.0]
    @test isinferred(output, rule)
end

@testset "logistic growth" begin

    init =  [1.0 4.0 7.0;
             2.0 5.0 8.0;
             3.0 6.0 9.0]

    test1 = [1.0  4.0  7.0;
             2.0  5.0  8.0;
             3.0  6.0  9.0]

    test2 = [1.81818  5.71429  8.23529;
             3.33333  6.66667  8.88889;
             4.61538  7.5      9.47368]

    test3 = [3.07692  7.27273  9.03226;
             5.0      8.0      9.41176;
             6.31579  8.57143  9.72973]
    rule = Ruleset(LogisticGrowth(rate=log(2.0), carrycap=10, timestep=1))
    output = ArrayOutput(init; tspan=1:3)

    sim!(output, rule)

    @test output[1] == test1
    @test output[2] ≈ test2 atol=1e-4
    @test output[3] ≈ test3 atol=1e-4
end

@testset "logistic growth with rate from a growth map" begin

    init = [1.0 1.0 1.0;
            1.0 1.0 1.0;
            1.0 1.0 1.0]

    test1 = [1.0 1.0 1.0;
             1.0 1.0 1.0;
             1.0 1.0 1.0]

    test2 = [1.0     1.0  1.81818;
             1.81818 1.0  0.5    ;
             1.0     1.0  0.5    ]

    test3 = [1.0      1.0  3.07692;
             3.07692  1.0  0.25   ;
             1.0      1.0  0.25   ]

    suit =  log.([1.0 1.0 2.0;
                  2.0 1.0 0.5;
                  1.0 1.0 0.5])

    output = ArrayOutput(init; tspan=1:3, aux=(suit=suit,))
    rule = Ruleset(LogisticGrowth(rate=Aux(:suit), carrycap=10, timestep=1))
    sim!(output, rule)

    @test output[1] == test1
    @test output[2] ≈ test2 atol=1e-4
    @test output[3] ≈ test3 atol=1e-4

end

@testset "growth map masking" begin
    init =  [1.0 1.0 1.0;
             1.0 1.0 1.0;
             1.0 1.0 1.0]

    suit =  [1.0 1.0 2.0;
             2.0 2.0 0.5;
             1.0 1.0 0.5]

    test1 = [1.0 1.0 1.0;
             1.0 1.0 1.0;
             1.0 1.0 1.0]

    test2 = [0.0 0.0 1.0;
             1.0 1.0 0.0;
             0.0 0.0 0.0]

    output = ArrayOutput(init; tspan=1:3, aux=(suit=suit,))
    maskrule = ThresholdGrowth(rate=Aux(:suit), threshold=1.1)
    ruleset = Ruleset(maskrule)
    data = SimData(extent(output), ruleset)

    @test applyrule(data, maskrule, 5, (1, 1)) === 0
    @test applyrule(data, maskrule, 5, (2, 2)) === 5

    @test applyrule(data, maskrule, 5.0, (1, 1)) === 0.0
    @test applyrule(data, maskrule, 5.0, (2, 2)) === 5.0

    sim!(output, Ruleset(maskrule))

    @test output[1] == test1
    @test output[2] == test2
end


text:
struct ApplyIterable{I, F}
    iter::I
    fun::F
    period::UInt
end

function iterate( iter::ApplyIterable, args... )
    next = iterate( iter.iter, args... )

    if ( next != nothing ) && ( ( next[1][1] % iter.period ) == 0 )
        iter.fun( next[1][2] )
    end

    return next
end

"""
```julia
    apply( iter::I, fun::F, period=1 ) where {I, F}
```

Apply the function `fun` to the iterator `iter` at every `period` samples.
`fun` should accept a single argument, which is the state variable for the iterator
`iter`.
"""
apply( iter::I, fun::F, period=1 ) where {I, F} = ApplyIterable{I, F}( iter, fun, period )


text:
mutable struct EventMetadata
    datetime::Dates.DateTime
    correlation_id::UUID
    causation_id::UUID
    published::Bool
    position::Union{Int,Nothing}
end

struct EventContext
    name::String
end

struct EventAggregate
    name::String
    id::Union{UUID,Nothing}
    revision::Union{Int,Nothing}
end

struct StreamInfo
    type::Symbol
    context_name::String
    aggregate_name::String
    event_name::String
end

function Base.print(io::IO, info::StreamInfo)
    print(io, "$(info.type).$(info.context_name).$(info.aggregate_name).$(info.event_name)")
end

function Base.parse(::Type{StreamInfo}, str::AbstractString)
    matched = match(r"^([^.]+)\.([^.]+)\.([^.]+)\.([^.]+)$", str)
    if isnothing(matched)
        throw(ArgumentError("invalid stream identifier $(repr(str))"))
    end

    StreamInfo(Symbol(matched[1]), matched[2], matched[3], matched[4])
end

const EventData = Dict
const EventCustom = Dict
const EventState = Dict

const EventTypes = Set([
    :command,
    :domain,
    :error
])

struct Event
    id::UUID
    publish_id::Union{String, Nothing}
    name::String
    type::Symbol
    context::EventContext
    aggregate::EventAggregate
    metadata::EventMetadata
    custom::EventCustom
    data::EventData
    state::Union{EventState,Nothing}
end

function Event(
    context_name::String,
    aggregate_name::String,
    event_name::String;
    type::Symbol=:domain,
    id::UUID=UUIDs.uuid4(),
    publish_id::Union{String,Nothing}=nothing,
    aggregate_id::Union{UUID,Nothing}=nothing,
    datetime::Dates.DateTime=Dates.now(Dates.UTC),
    previous::Union{Event,Nothing}=nothing,
    correlation_id::Union{UUID,Nothing}=nothing,
    causation_id::Union{UUID,Nothing}=nothing,
    published::Bool=false,
    custom::EventCustom=EventCustom(),
    data::EventData=EventData(),
    position::Union{Int,Nothing}=nothing,
    revision::Union{Int,Nothing}=nothing,
    state::Union{EventState,Nothing}=nothing
)
    # Validate the event
    if isempty(context_name) || !isnothing(match(r"\.", context_name))
        throw(ArgumentError("invalid context name $(repr(context_name))"))
    end

    if isempty(aggregate_name) || !isnothing(match(r"\.", aggregate_name))
        throw(ArgumentError("invalid aggregate name $(repr(aggregate_name))"))
    end

    if isempty(event_name) || !isnothing(match(r"\.", event_name))
        throw(ArgumentError("invalid event name $(repr(event_name))"))
    end

    if !(type in EventTypes)
        throw(ArgumentError("invalid event type $(repr(type))"))
    end

    if isnothing(correlation_id)
        if isnothing(previous)
            throw(ArgumentError("either a previous event or a correlation id must be provided"))
        end

        correlation_id = previous.metadata.correlation_id
    end

    if isnothing(causation_id)
        if isnothing(previous)
            throw(ArgumentError("either a previous event or a causation id must be provided"))
        end

        causation_id = previous.id
    end

    return Event(
        id,
        publish_id,
        event_name,
        type,
        EventContext(context_name),
        EventAggregate(aggregate_name, aggregate_id, revision),
        EventMetadata(datetime, correlation_id, causation_id, published, position),
        custom,
        data,
        state,
    )
end

function stream_info(event::Event) :: StreamInfo
    StreamInfo(event.type, event.context.name, event.aggregate.name, event.name)
end

function Base.parse(::Type{Event}, doc::Mongoc.BSON) :: Event
    Event(
        doc["id"],
        nothing,
        doc["name"],
        :domain,
        EventContext(
            doc["context"]["name"]
        ),
        EventAggregate(
            doc["aggregate"]["name"],
            doc["aggregate"]["id"],
            doc["aggregate"]["revision"],
        ),
        EventMetadata(
            doc["metadata"]["datetime"],
            doc["metadata"]["correlation_id"],
            doc["metadata"]["causation_id"],
            doc["metadata"]["published"],
            doc["metadata"]["position"],
        ),
        Dict(doc["custom"]),
        Dict(doc["data"]),
        nothing,
    )
end


text:
using Luna

a = 13e-6
gas = :Ar
pres = 5

τfwhm = 30e-15
λ0 = 800e-9
energy = 1e-6

L = 15e-2

grid = Grid.RealGrid(L, λ0, (160e-9, 3000e-9), 1e-12)

a0 = a
aL = 3a/4

afun = let a0=a0, aL=aL, L=L
    afun(z) = a0 + (aL-a0)*z/L
end
# coren, dens = Capillary.gradient(gas, L, pres, pres); # dens is unused
m = Capillary.MarcatiliMode(afun, gas, pres, loss=false, model=:full);

aeff(z) = Modes.Aeff(m, z=z)

densityfun = let dens0=PhysData.density(gas, pres)
    z -> dens0
end

ionpot = PhysData.ionisation_potential(gas)
ionrate = Ionisation.ionrate_fun!_ADK(ionpot)

responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),
             Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot))

linop, βfun! = LinearOps.make_linop(grid, m, λ0)

inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy)

Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff)

statsfun = Stats.default(grid, Eω, m, linop, transform; gas=gas, windows=((150e-9, 300e-9),))
output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)

Luna.run(Eω, grid, linop, transform, FT, output)

##
Plotting.pygui(true)
Plotting.stats(output)
Plotting.prop_2D(output)
Plotting.time_1D(output, [5e-2, 10e-2, 11e-2])
Plotting.spec_1D(output, [5e-2, 10e-2, 11e-2])


text:
module PrimaryRecursive

# Write your package code here.
include("PRBase.jl")
include("arithmetic.jl")

export p1, p2, p3, PrimRec, succ, zro, Comb, Proj,
    add, mult, pred, cosub, dff, sgn, nsgn, const1, remainder, ge,
    ifel, square, floor_sqrt, traingle, pair, fib, fst, scd
end


text:
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{T, D},
    ::Type{S},
) where {
    T <: PSY.HybridSystem,
    D <: AbstractHybridFormulation,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(T, sys)

    # Variables
    add_variables!(container, ActivePowerVariable, devices, D())
    add_variables!(container, ReactivePowerVariable, devices, D())
    add_variables!(container, ActivePowerInVariable, devices, D())
    add_variables!(container, ActivePowerOutVariable, devices, D())
    add_variables!(container, ComponentActivePowerVariable, devices, D())
    add_variables!(container, EnergyVariable, devices, D())
    if get_attribute(model, "reservation")
        add_variables!(container, ReservationVariable, devices, D())
    end

    # Parameters
    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)

    # Initial Conditions
    initial_conditions!(container, devices, D())

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ComponentActivePowerRangeExpressionLB,
        ComponentActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ComponentActivePowerRangeExpressionUB,
        ComponentActivePowerVariable,
        devices,
        model,
        S,
    )
    add_feedforward_arguments!(container, model, devices)

    if has_service_model(model)
        add_variables!(container, ComponentActivePowerReserveUpVariable, devices, D())
        add_variables!(container, ComponentActivePowerReserveDownVariable, devices, D())

        add_to_expression!(
            container,
            ComponentActivePowerRangeExpressionLB,
            ComponentActivePowerReserveDownVariable,
            devices,
            model,
            S,
        )
        add_to_expression!(
            container,
            ComponentActivePowerRangeExpressionUB,
            ComponentActivePowerReserveUpVariable,
            devices,
            model,
            S,
        )
        add_expressions!(container, ComponentReserveDownBalanceExpression, devices, model)
        add_expressions!(container, ComponentReserveUpBalanceExpression, devices, model)
    end
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{T, D},
    ::Type{S},
) where {
    T <: PSY.HybridSystem,
    D <: AbstractHybridFormulation,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(T, sys)

    # Constraints
    #TODO: Change to Expressions for Services support
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        OutputActivePowerVariableLimitsConstraint,
        ActivePowerOutVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        InputActivePowerVariableLimitsConstraint,
        ActivePowerInVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ComponentActivePowerVariableLimitsConstraint,
        ComponentActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ComponentActivePowerVariableLimitsConstraint,
        ComponentActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )

    add_constraints!(container, EnergyBalanceConstraint, devices, model, S)
    add_constraints!(container, PowerOutputRangeConstraint, devices, model, S)

    if has_service_model(model)
        add_constraints!(container, ReserveEnergyConstraint, devices, model, S)
        add_constraints!(container, RangeLimitConstraint, devices, model, S)
        add_constraints!(container, ComponentReserveUpBalance, devices, model, S)
        add_constraints!(container, ComponentReserveDownBalance, devices, model, S)
    end

    add_feedforward_constraints!(container, model, devices)

    # Cost Function
    objective_function!(container, devices, model, S)

    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{T, D},
    ::Type{S},
) where {T <: PSY.HybridSystem, D <: AbstractHybridFormulation, S <: PM.AbstractPowerModel}
    devices = get_available_components(T, sys)

    # Variables
    add_variables!(container, ActivePowerVariable, devices, D())
    add_variables!(container, ActivePowerInVariable, devices, D())
    add_variables!(container, ActivePowerOutVariable, devices, D())
    add_variables!(container, ComponentActivePowerVariable, devices, D())
    add_variables!(container, EnergyVariable, devices, D())

    add_variables!(container, ReactivePowerVariable, devices, D())
    add_variables!(container, ComponentReactivePowerVariable, devices, D())
    if get_attribute(model, "reservation")
        add_variables!(container, ReservationVariable, devices, D())
    end

    # Parameters
    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)

    # Initial Conditions
    initial_conditions!(container, devices, D())

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ReactivePowerBalance,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ComponentActivePowerRangeExpressionLB,
        ComponentActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ComponentActivePowerRangeExpressionUB,
        ComponentActivePowerVariable,
        devices,
        model,
        S,
    )

    add_feedforward_arguments!(container, model, devices)

    if has_service_model(model)
        add_variables!(container, ComponentActivePowerReserveUpVariable, devices, D())
        add_variables!(container, ComponentActivePowerReserveDownVariable, devices, D())
        add_to_expression!(
            container,
            ComponentActivePowerRangeExpressionLB,
            ComponentActivePowerReserveDownVariable,
            devices,
            model,
            S,
        )
        add_to_expression!(
            container,
            ComponentActivePowerRangeExpressionUB,
            ComponentActivePowerReserveUpVariable,
            devices,
            model,
            S,
        )
        add_expressions!(container, ComponentReserveDownBalanceExpression, devices, model)
        add_expressions!(container, ComponentReserveUpBalanceExpression, devices, model)
    end
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{T, D},
    ::Type{S},
) where {T <: PSY.HybridSystem, D <: AbstractHybridFormulation, S <: PM.AbstractPowerModel}
    devices = get_available_components(T, sys)

    # Constraints
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ComponentActivePowerVariableLimitsConstraint,
        ComponentActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ComponentActivePowerVariableLimitsConstraint,
        ComponentActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        InputActivePowerVariableLimitsConstraint,
        ActivePowerInVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        OutputActivePowerVariableLimitsConstraint,
        ActivePowerOutVariable,
        devices,
        model,
        S,
    )

    # Reactive power Constraints
    add_constraints!(
        container,
        ReactivePowerVariableLimitsConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ComponentReactivePowerVariableLimitsConstraint,
        ComponentReactivePowerVariable,
        devices,
        model,
        S,
    )

    # Constraints
    add_constraints!(container, EnergyBalanceConstraint, devices, model, S)
    add_constraints!(container, PowerOutputRangeConstraint, devices, model, S)
    add_constraints!(container, ReactivePowerConstraint, devices, model, S)
    add_constraints!(container, InterConnectionLimitConstraint, devices, model, S)
    if has_service_model(model)
        add_constraints!(container, ReserveEnergyConstraint, devices, model, S)
        add_constraints!(container, RangeLimitConstraint, devices, model, S)
        add_constraints!(container, ComponentReserveUpBalance, devices, model, S)
        add_constraints!(container, ComponentReserveDownBalance, devices, model, S)
    end
    add_feedforward_constraints!(container, model, devices)

    # Cost Function
    objective_function!(container, devices, model, S)

    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{T, D},
    ::Type{S},
) where {
    T <: PSY.HybridSystem,
    D <: StandardHybridDispatch,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(T, sys)

    # Variables
    add_variables!(container, ActivePowerVariable, devices, D())
    add_variables!(container, ReactivePowerVariable, devices, D())
    add_variables!(container, ActivePowerInVariable, devices, D())
    add_variables!(container, ActivePowerOutVariable, devices, D())
    add_variables!(container, ComponentActivePowerVariable, devices, D())
    add_variables!(container, EnergyVariable, devices, D())
    if get_attribute(model, "reservation")
        add_variables!(container, ReservationVariable, devices, D())
    end

    # Parameters
    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)

    # Initial Conditions
    initial_conditions!(container, devices, D())

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ComponentActivePowerRangeExpressionLB,
        ComponentActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ComponentActivePowerRangeExpressionUB,
        ComponentActivePowerVariable,
        devices,
        model,
        S,
    )

    add_feedforward_arguments!(container, model, devices)

    if has_service_model(model)
        add_variables!(container, ComponentActivePowerReserveUpVariable, devices, D())
        add_variables!(container, ComponentActivePowerReserveDownVariable, devices, D())
        add_to_expression!(
            container,
            ComponentActivePowerRangeExpressionLB,
            ComponentActivePowerReserveDownVariable,
            devices,
            model,
            S,
        )
        add_to_expression!(
            container,
            ComponentActivePowerRangeExpressionUB,
            ComponentActivePowerReserveUpVariable,
            devices,
            model,
            S,
        )
        add_expressions!(container, ComponentReserveDownBalanceExpression, devices, model)
        add_expressions!(container, ComponentReserveUpBalanceExpression, devices, model)
    end
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{T, D},
    ::Type{S},
) where {
    T <: PSY.HybridSystem,
    D <: StandardHybridDispatch,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(T, sys)

    # Constraints
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ComponentActivePowerVariableLimitsConstraint,
        ComponentActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ComponentActivePowerVariableLimitsConstraint,
        ComponentActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        InputActivePowerVariableLimitsConstraint,
        ActivePowerInVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        OutputActivePowerVariableLimitsConstraint,
        ActivePowerOutVariable,
        devices,
        model,
        S,
    )

    add_constraints!(container, EnergyBalanceConstraint, devices, model, S)
    add_constraints!(container, PowerOutputRangeConstraint, devices, model, S)

    if has_service_model(model)
        add_constraints!(container, ReserveEnergyConstraint, devices, model, S)
        add_constraints!(container, RangeLimitConstraint, devices, model, S)
        add_constraints!(container, ComponentReserveUpBalance, devices, model, S)
        add_constraints!(container, ComponentReserveDownBalance, devices, model, S)
    end

    add_feedforward_constraints!(container, model, devices)

    # Cost Function
    objective_function!(container, devices, model, S)

    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{T, D},
    ::Type{S},
) where {T <: PSY.HybridSystem, D <: StandardHybridDispatch, S <: PM.AbstractPowerModel}
    devices = get_available_components(T, sys)

    # Variables
    add_variables!(container, ActivePowerVariable, devices, D())
    add_variables!(container, ActivePowerInVariable, devices, D())
    add_variables!(container, ActivePowerOutVariable, devices, D())
    add_variables!(container, ComponentActivePowerVariable, devices, D())
    add_variables!(container, EnergyVariable, devices, D())

    add_variables!(container, ReactivePowerVariable, devices, D())
    add_variables!(container, ComponentReactivePowerVariable, devices, D())
    if get_attribute(model, "reservation")
        add_variables!(container, ReservationVariable, devices, D())
    end

    # Parameters
    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)

    # Initial Conditions
    initial_conditions!(container, devices, D())

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ReactivePowerBalance,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ComponentActivePowerRangeExpressionLB,
        ComponentActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ComponentActivePowerRangeExpressionUB,
        ComponentActivePowerVariable,
        devices,
        model,
        S,
    )

    add_feedforward_arguments!(container, model, devices)

    if has_service_model(model)
        add_variables!(container, ComponentActivePowerReserveUpVariable, devices, D())
        add_variables!(container, ComponentActivePowerReserveDownVariable, devices, D())
        add_to_expression!(
            container,
            ComponentActivePowerRangeExpressionLB,
            ComponentActivePowerReserveDownVariable,
            devices,
            model,
            S,
        )
        add_to_expression!(
            container,
            ComponentActivePowerRangeExpressionUB,
            ComponentActivePowerReserveUpVariable,
            devices,
            model,
            S,
        )
        add_expressions!(container, ComponentReserveDownBalanceExpression, devices, model)
        add_expressions!(container, ComponentReserveUpBalanceExpression, devices, model)
    end
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{T, D},
    ::Type{S},
) where {T <: PSY.HybridSystem, D <: StandardHybridDispatch, S <: PM.AbstractPowerModel}
    devices = get_available_components(T, sys)
    # Constraints
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ComponentActivePowerVariableLimitsConstraint,
        ComponentActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ComponentActivePowerVariableLimitsConstraint,
        ComponentActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        InputActivePowerVariableLimitsConstraint,
        ActivePowerInVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        OutputActivePowerVariableLimitsConstraint,
        ActivePowerOutVariable,
        devices,
        model,
        S,
    )

    # Reactive power Constraints
    add_constraints!(
        container,
        ReactivePowerVariableLimitsConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ComponentReactivePowerVariableLimitsConstraint,
        ComponentReactivePowerVariable,
        devices,
        model,
        S,
    )

    # Constraints
    add_constraints!(container, EnergyBalanceConstraint, devices, model, S)
    add_constraints!(container, PowerOutputRangeConstraint, devices, model, S)
    add_constraints!(container, ReactivePowerConstraint, devices, model, S)
    add_constraints!(container, InterConnectionLimitConstraint, devices, model, S)

    if has_service_model(model)
        add_constraints!(container, ReserveEnergyConstraint, devices, model, S)
        add_constraints!(container, RangeLimitConstraint, devices, model, S)
        add_constraints!(container, ComponentReserveUpBalance, devices, model, S)
        add_constraints!(container, ComponentReserveDownBalance, devices, model, S)
    end

    add_feedforward_constraints!(container, model, devices)

    # Cost Function
    objective_function!(container, devices, model, S)

    return
end


text:
module TestUtils
using Folds.Implementations: infer_eltype
using Test

@testset "infer_eltype" begin
    @test infer_eltype(1:10) == Int
    @test infer_eltype(string(x) for x in 1:10) == String
    @test infer_eltype(isodd(x) ? x : missing for x in 1:10) == Union{Int,Missing}
    @test infer_eltype(isodd(x) ? x : missing for x in (0, true)) == Union{Int,Bool,Missing}
end

end  # module


text:
using Documenter, SWIP

makedocs(;
    modules=[SWIP],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/zenna/SWIP.jl/blob/{commit}{path}#L{line}",
    sitename="SWIP.jl",
    authors="Zenna Tavares",
    assets=String[],
)

deploydocs(;
    repo="github.com/zenna/SWIP.jl",
)


text:
function hessian{T<:FloatingPoint}(f::Function,x::Array{T,1})

    m = length(x)

	e = eps(T)^(1/4)*Base.maxabs([x,one(T)])
	dh = eye(m)*e
    hess = Array(T,m,m)

    for i = 1:m
	    f1 = f(x+2*dh[:,i])
		f2 = f(x+dh[:,i])
		f3 = f(x)
		f4 = f(x-dh[:,i])
		f5 = f(x-2*dh[:,i])
		hess[i:i,i:i] = (-f1+16*f2-30*f3+16*f4-f5)/(12*e^2)
    end

    for i = 1:m-1
        for j = i+1:m
		    f1 = f(x+dh[:,i]+dh[:,j])
			f2 = f(x-dh[:,i]-dh[:,j])
			f3 = f(x-dh[:,i]+dh[:,j])
			f4 = f(x+dh[:,i]-dh[:,j])
			hess[i:i,j:j] = (f1+f2-f3-f4)/(4*e^2)
            hess[j,i] = hess[i,j]
        end
    end

    return hess

end

max_stars_repo_name: dingliumath/SolveDSGE

max_issues_repo_name: dingliumath/SolveDSGE

max_forks_repo_name: dingliumath/SolveDSGE


text:
i,j,m,d,seed = 55,24,10^9,34,292929     # parameters
s = Array{Int32}(undef,i); r = similar(s)
s[1:2] = [seed,1]                       # table initialization
for n = 3:i; (s[n] = s[n-2]-s[n-1]) < 0 && (s[n] += m) end
t = 1; for u=1:i; (global t+=d)>i && (t-=i); r[u]=s[t] end # permutation, r = s[(d*(1:i) .% i).+1]

u,v,n = i,i-j,i-1
while (n += 1) > 0
    (global u += 1) > i && (u = 1)      # circular indexing: u,v = ((n,n-j) .% i).+1
    (global v += 1) > i && (v = 1)
    (r[u] -= r[v]) < 0 && (r[u] += m)   # table update
    n < 220 && continue                 # 165 silent values
    print((n,r[u]))                     # show (index,value) of next pseudorandom number
    x = readline(stdin)                 # wait until the ENTER key is pressed
    length(x) > 0 && break              # any other key before ENTER => exit
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CUDNN_CUDA102_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("CUDNN_CUDA102")
JLLWrappers.@generate_main_file("CUDNN_CUDA102", UUID("8ba91e8f-b0e5-5ca0-b631-aadbb9431ebf"))
end  # module CUDNN_CUDA102_jll


text:
################################################################################
# MIT License
#
# Copyright (c) 2021 Hajime Nakagami<nakagami@gmail.com>
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
################################################################################
using Decimals

function dpb_bit_to_int64(dpd::UInt64, mask::UInt64)::UInt64
    (dpd & mask) != 0 ? Int64(1) : Int64(0)
end

function dpd_to_int64(dpd::UInt)::UInt64
    # Convert DPD encodined value to int (0-999)
    # dpd: DPD encoded value. 10bit unsigned int

    b = zeros(UInt64, 10)
    b[10] = dpdBitToInt(dpd, 0x0200)
    b[9] = dpdBitToInt(dpd, 0x0100)
    b[8] = dpdBitToInt(dpd, 0x0080)
    b[7] = dpdBitToInt(dpd, 0x0040)
    b[6] = dpdBitToInt(dpd, 0x0020)
    b[5] = dpdBitToInt(dpd, 0x0010)
    b[4] = dpdBitToInt(dpd, 0x0008)
    b[3] = dpdBitToInt(dpd, 0x0004)
    b[2] = dpdBitToInt(dpd, 0x0002)
    b[1] = dpdBitToInt(dpd, 0x0001)

    b = zeros(UInt64, 3)
    if b[4] == 0
        d[3] = b[10]*4 + b[9]*2 + b[8]
        d[2] = b[7]*4 + b[6]*2 + b[5]
        d[1] = b[3]*4 + b[2]*2 + b[1]
    elseif b[4] == 1 && b[3] == 0 && b[2] == 0
        d[3] = b[10]*4 + b[9]*2 + b[8]
        d[2] = b[7]*4 + b[6]*2 + b[5]
        d[1] = 8 + b[1]
    elseif b[4] == 1 && b[3] == 0 && b[2] == 1
        d[3] = b[10]*4 + b[9]*2 + b[8]
        d[2] = 8 + b[5]
        d[1] = b[7]*4 + b[6]*2 + b[1]
    elseif b[4] == 1 && b[3] == 1 && b[2] == 0
        d[3] = 8 + b[8]
        d[2] = b[7]*4 + b[6]*2 + b[5]
        d[1] = b[10]*4 + b[9]*2 + b[1]
    elseif b[7] == 0 && b[6] == 0 && b[4] == 1 && b[3] == 1 && b[2] == 1
        d[3] = 8 + b[8]
        d[2] = 8 + b[5]
        d[1] = b[10]*4 + b[9]*2 + b[1]
    elseif b[7] == 0 && b[6] == 1 && b[4] == 1 && b[3] == 1 && b[2] == 1
        d[3] = 8 + b[8]
        d[2] = b[10]*4 + b[9]*2 + b[5]
        d[1] = 8 + b[1]
    elseif b[7] == 1 && b[6] == 0 && b[4] == 1 && b[3] == 1 && b[2] == 1
        d[3] = b[10]*4 + b[9]*2 + b[8]
        d[2] = 8 + b[5]
        d[1] = 8 + b[1]
    elseif b[7] == 1 && b[6] == 1 && b[4] == 1 && b[3] == 1 && b[2] == 1
        d[3] = 8 + b[8]
        d[2] = 8 + b[5]
        d[1] = 8 + b[1]
    else
        throw(DomainError(plugin_name, "Invalid DPD encoding"))
    end

    int64(d[3])*100 + int64(d[2])*10 + int64(d[1])
end

function calc_significand(prefix::Int64, dpd_bits::BigInt, numBits::Int)::BigInt
    # prefix: High bits integer value
    # dpdBits: dpd encoded bits
    # numBits: bit length of dpd_bits
    # https://en.wikipedia.org/wiki/Decimal128_floating-point_format#Densely_packed_decimal_significand_field
    num_segments = div(numBits, 10)
    segments::Vector{BigInt} = []

    for _ in 1:num_segments
        segments.append(dpd_bits & 0b1111111111)
        dpd_bits >>= 10
    end

    v = prefix
    for dpd in segments
        v = v + 1000 + dpd_to_int64(dpd)
    end
    v
end

function decimal128_to_sign_digits_exponent(b::Vector{UInt8})::Union{Decimal, Tuple{Int, BigInt, Int32}}
    # https://en.wikipedia.org/wiki/Decimal128_floating-point_format
    sign::Int = 0
    digits::BigInt = 0
    exponent::Int32 = 0

    prefix::Int64 = 0

    if (b[0] & 0x80) == 0x80
        sign = 1
    end
    cf = (UInt32(b[0]&0x7f) << 10) + UInt32(b[1]<<2) + UInt32(b[2]>>6)
    if (cf & 0x1F000) == 0x1F000
        if sign == 1
            return Deicmal(-NaN)
        else
            return Decimal(Nan)
        end
    elseif (cf & 0x1F000) == 0x1E000
        if sign == 1
            return Decimal(Inf)
        else
            return Decimal(-Inf)
        end
    elseif (cf & 0x18000) == 0x00000
        exponent = int32(0x0000 + (cf & 0x00fff))
        prefix = int64((cf >> 12) & 0x07)
    elseif (cf & 0x18000) == 0x08000
        exponent = int32(0x1000 + (cf & 0x00fff))
        prefix = Int64((cf >> 12) & 0x07)
    elseif (cf & 0x18000) == 0x10000
        exponent = Int32(0x2000 + (cf & 0x00fff))
        prefix = Int64((cf >> 12) & 0x07)
    elseif (cf & 0x1e000) == 0x18000
        exponent = int32(0x0000 + (cf & 0x00fff))
        prefix = Int64(8 + (cf>>12)&0x01)
    elseif (cf & 0x1e000) == 0x1a000
        exponent = Int32(0x1000 + (cf & 0x00fff))
        prefix = Int64(8 + (cf>>12)&0x01)
    elseif (cf & 0x1e000) == 0x1c000
        exponent = int32(0x2000 + (cf & 0x00fff))
        prefix = int64(8 + (cf>>12)&0x01)
    else
        throw(DomainError(b, "decimal128 value error"))
    end
    exponent -= 6176

    dbd_bits &= big"0x3fffffffffffffffffffffffffff"
    digits = calc_significand(prefix, dpd_bits, 110)

    (sign, digits, exponent)
end

function decimal_fixed_to_decimal(b::Vector{UInt8}, scale::Int32)::Decimal
    v = decimal128_to_sign_digits_exponent(b)
    if isa(v, Decimal)
        return v
    end
    sign, digits, _ = v
    if sign != 0
        digits *= -1
    end

    Decimal(digits, scale)
end

function decimal64_to_decimal(b::Vector{UInt8})::Decimal
    # https://en.wikipedia.org/wiki/Decimal64_floating-point_format
    prefix::Int64 = 0
    sign::Int = 0
    if (b[0] & 0x80) == 0x80
        sign = 1
    end
    cf = (UInt32(b[0]) >> 2) & 0x1f
    exponent = ((Int32(b[0]) & 3) << 6) + ((Int32(b[1]) >> 2) & 0x3f)

    dpd_bits = reinterpret(BigInt, b)
    dpd_bits &= big"0x3ffffffffffff"

    if cf == 0x1f
        if sign == 1
            return Deicmal(-NaN)
        else
            return Deicmal(NaN)
        end
    elseif cf == 0x1e
        if sign == 1
            return Decimal(-Inf)
        else
            return Decimal(Inf)
        end
    elseif (cf & 0x18) == 0x00
        exponent = 0x000 + exponent
        prefix = Int64(cf & 0x07)
    elseif (cf & 0x18) == 0x08
        exponent = 0x100 + exponent
        prefix = Int64(cf & 0x07)
    elseif (cf & 0x18) == 0x10
        exponent = 0x200 + exponent
        prefix = Int64(cf & 0x07)
    elseif (cf & 0x1e) == 0x18
        exponent = 0x000 + exponent
        prefix = Int64(8 + cf&1)
    elseif (cf & 0x1e) == 0x1a
        exponent = 0x100 + exponent
        prefix = Int64(8 + cf&1)
    elseif (cf & 0x1e) == 0x1c
        exponent = 0x200 + exponent
        prefix = Int64(8 + cf&1)
    else
        throw(DomainError(b, "decimal64 value error"))
    end
    digits = calc_significand(prefix, dpd_bits, 50)
    exponent -= 398

    if sign != 0
        digits *= big"-1"
    end

    Decimal(digits, eponent)
end

function decimal128_to_decimal(b::Vector{UInt8})::Decimal
    # https://en.wikipedia.org/wiki/Decimal64_floating-point_format
    v = decimal128_to_sign_digits_exponent(b)
    if isa(v, Decimal)
        return v
    end
    sign, _, exponent = v
    if sign != 0
        digits *= -1
    end

    Decimal(digits, exponent)
end


text:
using VectorMathOptInterface
using Test

@testset "VectorMathOptInterface.jl" begin
    # Write your tests here.
end


text:
module ReinforcementLearningBase

const RLBase = ReinforcementLearningBase
export RLBase

using Random

include("inline_export.jl")
include("interface.jl")
include("CommonRLInterface.jl")
include("base.jl")

end # module


text:
## Create an extractor that extracts the modified original matrix, e.g.
## LD for BunchKaufman, UL for CholeskyDense, LU for LUDense and
## define size methods for Factorization types using it.

immutable BunchKaufman{T} <: Factorization{T}
    LD::Matrix{T}
    ipiv::Vector{BlasInt}
    uplo::Char
    symmetric::Bool
end

function bkfact!{T<:BlasReal}(A::StridedMatrix{T}, uplo::Symbol=:U, symmetric::Bool=issym(A))
    symmetric || error("The Bunch-Kaufman decomposition is only valid for symmetric matrices")
    LD, ipiv = LAPACK.sytrf!(string(uplo)[1] , A)
    BunchKaufman(LD, ipiv, string(uplo)[1], symmetric)
end
function bkfact!{T<:BlasComplex}(A::StridedMatrix{T}, uplo::Symbol=:U, symmetric::Bool=issym(A))
    LD, ipiv = (symmetric ? LAPACK.sytrf! : LAPACK.hetrf!)(string(uplo)[1] , A)
    BunchKaufman(LD, ipiv, string(uplo)[1], symmetric)
end
bkfact{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U, symmetric::Bool=issym(A)) = bkfact!(copy(A), uplo, symmetric)
bkfact{T}(A::StridedMatrix{T}, uplo::Symbol=:U, symmetric::Bool=issym(A)) = bkfact!(convert(Matrix{promote_type(Float32,typeof(sqrt(one(T))))},A),uplo,symmetric)

convert{T}(::Type{BunchKaufman{T}},B::BunchKaufman) = BunchKaufman(convert(Matrix{T},B.LD),B.ipiv,B.uplo,B.symmetric)
convert{T}(::Type{Factorization{T}}, B::BunchKaufman) = convert(BunchKaufman{T}, B)

size(B::BunchKaufman) = size(B.LD)
size(B::BunchKaufman,d::Integer) = size(B.LD,d)
issym(B::BunchKaufman) = B.symmetric
ishermitian(B::BunchKaufman) = !B.symmetric

inv{T<:BlasReal}(B::BunchKaufman{T})=copytri!(LAPACK.sytri!(B.uplo, copy(B.LD), B.ipiv), B.uplo, true)

function inv{T<:BlasComplex}(B::BunchKaufman{T})
    if issym(B)
        copytri!(LAPACK.sytri!(B.uplo, copy(B.LD), B.ipiv), B.uplo)
    else
        copytri!(LAPACK.hetri!(B.uplo, copy(B.LD), B.ipiv), B.uplo, true)
    end
end

A_ldiv_B!{T<:BlasReal}(B::BunchKaufman{T}, R::StridedVecOrMat{T}) = LAPACK.sytrs!(B.uplo, B.LD, B.ipiv, R)
function A_ldiv_B!{T<:BlasComplex}(B::BunchKaufman{T}, R::StridedVecOrMat{T})
    (issym(B) ? LAPACK.sytrs! : LAPACK.hetrs!)(B.uplo, B.LD, B.ipiv, R)
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
module VerFedLogistic

using LinearAlgebra
using Printf
using SparseArrays
using Random
using SharedArrays
using Distributed
using Combinatorics
using StatsBase
using Flux
using Zygote

export Client, AsynClient
export Server, AsynServer
export connect
export update_batch
export send_embedding, update_embedding
export update_model, update_grads, compute_mini_batch_gradient, send_gradient
export eval
export softmax, neg_log_loss
export load_data, split_data, generate_batches
export vertical_lr_train, evaluation
export read_libsvm

include("./utils.jl")
include("./client.jl")
include("./server.jl")
include("./training.jl")

end # module


text:
const _DEFAULT_BEFORE_SCRIPT = """
    import Pkg

    Pkg.Registry.add("General")
"""

const _DEFAULT_SCRIPT = """
    import Pkg

    Pkg.instantiate()
    Pkg.update()
"""

const _DEFAULT_AFTER_SCRIPT = """
    import Pkg
"""


text:
using OmegaGrav
using Dierckx
using Tables, CSV
# %% Compute the matter power spectrum using CLASS
include("compute_pk_class.jl")
# %% Compute Ωtherm at seven redshifts
redshift = [0, 0.3, 0.5, 0.7, 1, 1.3, 1.5]
nred = length(redshift)
Ωth = zeros(nred)
Ωthu = zeros(nred)
Ωthl = zeros(nred)
ΩthB1 = zeros(nred)
for ired = 1:nred
   z = redshift[ired]
   # Define functions to return linear baryon+CDM power spectrum
   # Note: The CLASS code takes wavenumbers in units of 1/Mpc (no h) and
   # return power spectra in units of Mpc^3 (no 1/h^3).
   pkcb_class(kovh) = cosmo.pk_cb_lin(kovh * h0, z) * h0^3
   # Spline interpolate in log(k)
   lnk = log(1e-4):0.05:log(100)
   pkcb = Spline1D(lnk, pkcb_class.(exp.(lnk)))
   # %% Compute Ωtherm from Halos, excluding the neutrino contribution
   Ωth[ired] = otherm_upp(x -> pkcb(log(x)), z, Ωm, h0, Ωcb)
   # %% Compute Ωtherm for upper and lower 68% confidence level in B
   Ωthl[ired] = otherm_upp(x -> pkcb(log(x)), z, Ωm, h0, Ωcb, massbias = 1.315)
   Ωthu[ired] = otherm_upp(x -> pkcb(log(x)), z, Ωm, h0, Ωcb, massbias = 1.221)
   # %% Compute Ωtherm for no mass bias case for comparison
   ΩthB1[ired] = otherm_upp(x -> pkcb(log(x)), z, Ωm, h0, Ωcb, massbias = 1)
end

#%% Save table data to otherm_upp.csv
fb = params["omega_b"] / (params["omega_cdm"] + params["omega_b"])
t = Tables.table([redshift Ωth Ωth - Ωthl Ωthu - Ωth ΩthB1])
header = [
   "z",
   "Omega_th_B1p266",
   "Omega_th_err_low",
   "Omega_th_err_high",
   "Omega_th_B1",
]
CSV.write("otherm_upp.csv", t, header = header)

# %% Clean CLASS (the equivalent of the struct_free() in the `main`
# of CLASS. This step is primordial when running in a loop over different
# cosmologies, as you will saturate your memory very fast if you ommit
# it.
cosmo.struct_cleanup()
# If you want to change completely the cosmology, you should also
# clean the arguments, otherwise, if you are simply running on a loop
# of different values for the same parameters, this step is not needed
cosmo.empty()


text:
module GRMakie

using AbstractPlotting
using Colors, Printf
import GR

import AbstractPlotting: backend_showable, backend_show

using AbstractPlotting.IntervalSets: ClosedInterval

include("utils.jl")
include("primitives.jl")
include("saving.jl")

function activate!()
    AbstractPlotting.current_backend[] = GRBackend()
end

function __init__()
    activate!()
    ENV["GKS_ENCODING"] = haskey(ENV, "GKS_ENCODING") ? ENV["GKS_ENCODING"] : "utf-8"
    GR.selntran(0)
end

end # module


text:
# # Resource management

# Let's see how resources get freed automatically, and when they aren't. First let's set the preference `"LOG_DESTRUCTION"` to `"true"` to see what's going on:

using SwiftShader_jll # hide
using Vulkan
@set_driver :SwiftShader # hide

Vulkan.set_preferences!("LOG_DESTRUCTION" => "true")

# Now let's create a bunch of handles and see what happens when the GC runs:

function do_something()
    instance = Instance([], [])
    physical_devices = unwrap(enumerate_physical_devices(instance))
    physical_device = first(physical_devices)
    device = Device(physical_device, [DeviceQueueCreateInfo(0, [1.0])], [], [])

    command_pool = CommandPool(device, 0)
    buffers = unwrap(
        allocate_command_buffers(
            device,
            CommandBufferAllocateInfo(command_pool, COMMAND_BUFFER_LEVEL_PRIMARY, 10),
        ),
    )
    ## they won't be automatically freed individually
    free_command_buffers(device, command_pool, buffers)
    nothing
end

do_something()

## to force some finalizers to run
GC.gc()

#=

Not all handles were destroyed upon finalization. In particular, the physical device and the buffers were not destroyed. That's because a physical device is not owned by the application, so you can't destroy it, and buffers must be freed in batches with `free_command_buffers` (as done above). See more information in [Automatic finalization](@ref).

!!! note
    Not all finalizers have run. In some cases (e.g. when a finalizer must be run to release resources), it may be preferable to run them directly.
    You can do this by calling `finalize` (exported from `Base`):

    ```@example resource_management
    instance = Instance([], [])
    finalize(instance)
    ```
=#

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
sample_dict = Dict(
    "Key" => "Value",
    "Some other key" => 2323,
)

sample_dict["Key"]

sample_dict["New Key"] = "New Value"
sample_dict["Key"] = "Override"
sample_dict

pop!(sample_dict, "Some other key")
sample_dict

# tuples, looks like python
animals = ("penguins", "cats", "dogs")

# Julia is 1 indexed...
animals[1]

# Tuples are immutable, cant do this:
# animals[1] = "otters"

# Arrays/Vectors
people = ["Tom", "Scott", "Dave", "Mike", "John"]
fibo = [1, 1, 2, 3, 5, 8, 13]

# 4-element Vector{Any}:
mix = [ people[1], fibo[2], people[3], fibo[4] ]

# Push/Pop
push!(fibo, 21)
pop!(fibo)
push!(fibo, 21, 34, 55)

array_of_arrays = [
    people,
    fibo,
    mix,
]
array_of_new_arrays = [
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9],
    [Nothing, 0, "."]
]

# 1D, 2D, 3D arrays of random numbers.
rand(3)
rand(4, 5)
rand(2, 2, 3)

max_stars_repo_name: AustinTSchaffer/DailyProgrammer

max_issues_repo_name: AustinTSchaffer/DailyProgrammer

max_forks_repo_name: AustinTSchaffer/DailyProgrammer


text:
using SCS
const SCSOPT = with_optimizer(SCS.Optimizer, verbose = false)

#TODO
# @testset "SCS Exponential Cone Problems" begin
#     list_of_exp_problems = [
#         exp1_test,
#         exp2_test
#     ]
#     test_strong_duality(list_of_exp_problems, SCSOPT)
# end

@testset "SCS Power Cone Problems" begin
    list_of_pow_problems = [
        pow1_test,
        pow2_test
    ]
    test_strong_duality(list_of_pow_problems, SCSOPT; atol = 1e-3, rtol = 1e-3)
end


text:
using Test, Ariadne

@testset "ariadne" begin
    # Sadly just a stub
end

max_stars_repo_name: basp/ariadne

max_issues_repo_name: basp/ariadne

max_forks_repo_name: basp/ariadne


text:
## Indexing

function range_indecies(idx::Tuple)
    x = Vector{Any}(undef, length(idx))
    for (i,n) in enumerate(idx)
        x[i] = 1:n
    end
    x
end

replace_last_index!(idx::Vector, x) = (idx[end] = x; idx)

function assign!(A::AbstractArray, B::AbstractArray{T,N}; last_dim=1:size(B,N)) where {T,N}
    A_dims, B_dims = size(A), size(B)
    @assert A_dims[1:end-1] == B_dims[1:end-1] "Inconsistent dimensions with $(A_dims[1:end-1]) and $(B_dims[1:end-1])"
    A_dims = replace_last_index!(range_indecies(A_dims), last_dim)
    B_dims = range_indecies(B_dims)
    A_idxs = CartesianIndices(Tuple(A_dims))
    B_idxs = CartesianIndices(Tuple(B_dims))
    for (Aidx, Bidx) = zip(A_idxs, B_idxs)
        A[Aidx] = B[Bidx]
    end
    A
end

## Top-k pooling

function topk_index(y::AbstractVector, k::Integer)
    v = nlargest(k, y)
    return collect(1:length(y))[y .>= v[end]]
end

topk_index(y::Adjoint, k::Integer) = topk_index(y', k)

## Get feature with defaults

get_feature(::Nothing, i::Integer) = zeros(0)
get_feature(A::AbstractMatrix, i::Integer) = (i ≤ size(A,2)) ? view(A, :, i) : zeros(0)

function bypass_graph(nf_func=identity, ef_func=identity, gf_func=identity)
    return function (fg::FeaturedGraph)
        FeaturedGraph(graph(fg), nf_func(node_feature(fg)), ef_func(edge_feature(fg)),
                      gf_func(global_feature(fg)))
    end
end


text:
import HTTP
import JSON

const _httpjson_github_owner = "bcbi"
const _httpjson_github_repo = "ModelSanitizer.jl"

function _httpjson__get_github_pull_request_title_without_number_unauthenticated(travis_pull_request::AbstractString)::String
    _travis_pull_request::String = lowercase(strip(travis_pull_request))
    if length(_travis_pull_request) == 0 || _travis_pull_request == "false"
        return ""
    else
        pull_request_json = HTTP.request(
            "GET",
            "https://api.github.com/repos/$(github_owner)/$(github_repo)/pulls/$(_travis_pull_request)",
            [
                "User-Agent" => "bcbi-bot",
                ],
            )
        pull_request_dict = JSON.parse(String(deepcopy(pull_request_json.body)))
        pull_request_title::String = strip(pull_request_dict["title"])
        return pull_request_title
    end
end

function _httpjson_get_github_pull_request_title_without_number_unauthenticated(d::AbstractDict)
    travis_pull_request::String = lowercase(strip(get(d, "TRAVIS_PULL_REQUEST", "false")))
    result = _httpjson__get_github_pull_request_title_without_number_unauthenticated(travis_pull_request)
    return result
end

function _httpjson_get_github_pull_request_title_without_number_unauthenticated()::String
    result::String = _httpjson_get_github_pull_request_title_without_number_unauthenticated(ENV)
    return result
end

function _httpjson_get_github_pull_request_title_unauthenticated(d::AbstractDict)::String
    travis_pull_request::String = lowercase(strip(get(d, "TRAVIS_PULL_REQUEST", "false")))
    if length(travis_pull_request) == 0 || travis_pull_request == "false"
        return ""
    else
        result_without_number::String = _httpjson_get_github_pull_request_title_without_number_unauthenticated(d)
        result_with_number::String = string(travis_pull_request, ": ", result_without_number)
    end
end

function _httpjson_get_github_pull_request_title_unauthenticated()::String
    result::String = _httpjson_get_github_pull_request_title_unauthenticated(ENV)
    return result
end


text:
import ForwardDiff

struct ForwardDiffGradAlg <: GradAlg end
const ForwardDiffGrad = ForwardDiffGradAlg()

function lineargradient(rv, ω::Ω, ::ForwardDiffGradAlg)
  rv(ω) # Init

  # Unlinearizes xs into ω::Ω and applies rv(ω)
  function unlinearizeapl(xs)
    # Replace the tag tag here
    # @assert false
    # zazoom(1)
    # ω_ = tag(ω, dsofttrue(ForwardDiff.Dual))
    # T = ForwardDiff.Dual{Any,Float64,0}
    # rv2 = ciid(ω -> apl(rv, @show(Omega.Soft.tagerror(ω, dsofttrue(T)))))
    # ω_ = Omega.tagerror(ω, dsofttrue(T))
    # then we dont want to override the tag
    # because rv is going to add a tagz
    ωu = unlinearize(xs, ω)
    res = apl(rv, ωu)
    # logerr(res)
    # println(res)
    # res
    # @assert false
  end
  xs = linearize(ω)
  ForwardDiff.gradient(unlinearizeapl, xs)
end

value(x::ForwardDiff.Dual) = x.value


text:
using StanQuap
using Test

include("./../scripts/howell1.jl")


text:
abstract type AbstractIntervalEigenSolver end

struct Hertz <: AbstractIntervalEigenSolver end
struct Rohn <: AbstractIntervalEigenSolver end

"""
    eigenbox(A[, method=Rohn()])

Returns an enclosure of all the eigenvalues of `A`. If `A` is symmetric, than the
output is a real interval, otherwise it is a complex interval.

### Input

- `A` -- square interval matrix
- `method` -- method used to solve the symmetric interval eigenvalue problem (bounding
    eigenvalues of general matrices is also reduced to the symmetric case).
    Possible values are

      - `Rohn` -- (default) fast method to compute an enclosure of the eigenvalues of
            a symmetric interval matrix
      - `Hertz` -- finds the exact hull of the eigenvalues of a symmetric interval
            matrix, but has exponential complexity.

### Algorithm

The algorithms used by the function are described in [[HLA13]](@ref).

### Notes

The enclosure is not rigorous, meaning that the real eigenvalue problems solved internally
utilize normal floating point computations.

### Examples

```jldoctest
julia> A = [0 -1 -1; 2 -1.399.. -0.001 0; 1 0.5 -1]
3×3 Matrix{Interval{Float64}}:
 [0, 0]  [-1, -1]                       [-1, -1]
 [2, 2]       [-1.39901, -0.000999999]    [0, 0]
 [1, 1]        [0.5, 0.5]               [-1, -1]

julia> eigenbox(A)
[-1.90679, 0.970154] + [-2.51903, 2.51903]im

julia> eigenbox(A, Hertz())
[-1.64732, 0.520456] + [-2.1112, 2.1112]im
```
"""
function eigenbox(A::Symmetric{Interval{T}, Matrix{Interval{T}}}, ::Rohn) where {T}

    AΔ = Symmetric(radius.(A))
    Ac = Symmetric(mid.(A))

    ρ = eigmax(AΔ)
    λmax = eigmax(Ac)
    λmin = eigmin(Ac)
    return Interval(λmin - ρ, λmax + ρ)

end

function eigenbox(A::Symmetric{Interval{T}, Matrix{Interval{T}}}, ::Hertz) where {T}

    n = checksquare(A)
    Amax = Matrix{T}(undef, n, n)
    Amin = Matrix{T}(undef, n, n)

    λmin = Inf
    λmax = -Inf
    @inbounds for z in Orthants(n)
        first(z) < 0 && continue
        for j in 1:n
            for i in 1:j
                if z[i] == z[j]
                    Amax[i, j] = sup(A[i, j])
                    Amin[i, j] = inf(A[i, j])
                else
                    Amax[i, j] = inf(A[i, j])
                    Amin[i, j] = sup(A[i, j])
                end
            end
        end

        candmax = eigmax(Symmetric(Amax))
        candmin = eigmin(Symmetric(Amin))
        λmin = min(λmin, candmin)
        λmax = max(λmax, candmax)
    end
    return IA.Interval(λmin, λmax)
end

function eigenbox(A::AbstractMatrix{Interval{T}},
                  method::AbstractIntervalEigenSolver) where {T}

    λ = eigenbox(Symmetric(0.5*(A + A')), method)

    n = size(A, 1)
    μ = eigenbox(Symmetric([zeros(n, n) 0.5*(A - A');
                            0.5*(A' - A) zeros(n, n)]), method)

    return λ + μ*im
end

function eigenbox(M::AbstractMatrix{Complex{Interval{T}}},
                  method::AbstractIntervalEigenSolver) where {T}
    A = real.(M)
    B = imag.(M)
    λ = eigenbox(Symmetric(0.5*[A+A' B'-B;
                                B-B' A+A']), method)

    μ = eigenbox(Symmetric(0.5*[B+B' A-A';
                                A'-A B+B']), method)

    return λ + μ*im
end

function eigenbox(M::Hermitian{Complex{Interval{T}}, Matrix{Complex{Interval{T}}}},
                  method::AbstractIntervalEigenSolver) where {T}
    A = real(M)
    B = imag(M)
    return eigenbox(Symmetric([A B';B A]), method)
end

# default
eigenbox(A) = eigenbox(A, Rohn())


text:
"""
Functor to update JS part when an observable changes.
We make this a Functor, so we can clearly identify it and don't sent
any updates, if the JS side requires to update an Observable
(so we don't get an endless update cycle)
"""
struct JSUpdateObservable
    session::Session
    id::String
end

function (x::JSUpdateObservable)(value)
    # Sent an update event
    send(x.session, payload = value, id = x.id, msg_type = UpdateObservable)
end

"""
Update the value of an observable, without sending changes to the JS frontend.
This will be used to update updates from the forntend.
"""
function update_nocycle!(obs::Observable, value)
    setindex!(obs, value, notify = (f-> !(f isa JSUpdateObservable)))
end

function jsrender(session::Session, obs::Observable)
    html = map(obs) do data
        repr_richest(jsrender(session, data))
    end
    dom = DOM.m_unesc("span", html[])
    onjs(session, html, js"""
        function (html){
            var dom = $(dom);
            if(dom){
                dom.innerHTML = html;
                return true;
            }else{
                //deregister the callback if the observable dom is gone
                return false;
            }
        }
    """)
    return dom
end


text:
export get_cstind

"""
    test_cstind(
        rom_name::String,
        ind::NSGA2ECInd,
        ecfg::NamedTuple,
        ccfg::NamedTuple,
        test_action::Int32
    )

Test that this individual only outputs this action.
"""
function test_cstind(
    rom_name::String,
    e_chromosome::Vector{Float64},
    c_chromosome::Vector{Float64},
    mcfg::NamedTuple,
    ecfg::NamedTuple,
    ccfg::NamedTuple,
    test_action::Int32,
    reducer::Reducer
)
    game = Game(rom_name, 0)
    enco = IPCGPInd(ecfg, e_chromosome)
    cont = CGPInd(ccfg, c_chromosome)
    for i in 1:100
        s = get_state(game, mcfg.grayscale, mcfg.downscale)
        output = IICGP.process(enco, reducer, cont, ccfg, s)
        action = game.actions[argmax(output)]
        @assert action == test_action
        act(game.ale, action)
        if game_over(game.ale)
            break
        end
    end
    close!(game)
end

"""
    get_cstind(
        indtype::Type,
        rom_name::String,
        mcfg::NamedTuple,
        ecfg::NamedTuple,
        ccfg::NamedTuple,
        reducer::Reducer
    )

Build individuals constantly outputing the same actions.
"""
function get_cstind(
    indtype::Type,
    rom_name::String,
    mcfg::NamedTuple,
    ecfg::NamedTuple,
    ccfg::NamedTuple,
    reducer::Reducer
)
    game = Game(rom_name, 0)
    actions = game.actions
    close!(game)
    @assert ccfg.n_cst_inputs > 2
    cstinds = Vector{indtype}()
    R = ccfg.rows
    C = ccfg.columns
    P = ccfg.n_parameters
    M = R * C + ccfg.n_in
    for i in eachindex(actions)
        enco = IPCGPInd(ecfg)
        cont = CGPInd(ccfg)
        for j in eachindex(cont.outputs)
            if j != i
                #cont.outputs[j] = ccfg.n_in - ccfg.n_cst_inputs + 1
                cont.chromosome[R*C*(3+P)+j] = (ccfg.n_in - ccfg.n_cst_inputs + 1) / M
            else
                #cont.outputs[j] = ccfg.n_in
                cont.chromosome[R*C*(3+P)+j] = (ccfg.n_in) / M
            end
        end
        new_ind = indtype(mcfg, enco.chromosome, cont.chromosome)
        test_cstind(rom_name, new_ind.e_chromosome, new_ind.c_chromosome, mcfg,
            ecfg, ccfg, actions[i], reducer)
        push!(cstinds, new_ind)
    end
    cstinds
end


text:
module SparseSymmetricTensors
#=------------------------------------------------------------------------------
   Main File of the COOTensor class, Only type definitions and Constructor
   functions should be placed here.
------------------------------------------------------------------------------=#

using Base.Cartesian
using Printf
using SparseArrays
using DataStructures
#using StaticArrays

import Combinatorics.permutations, Combinatorics.multinomial
import LinearAlgebra.eigen, LinearAlgebra.norm, LinearAlgebra.dot
import Arpack.eigs, Arpack.svds
import Base: ==

export COOTen
export SSSTensor

#TODO: SSSTensor([i,j,k]) returns abnormal result, this should be altered to set
# default hyperedge weight to 1.

#TODO: should raise an error message for an empty iterator in the contstuctor
#TODO: possible to create a zero tensor from Dense Array, but not list constructor
#      need to standardize this.
#TODO: possible to condense constructor definitions with default params
#TODO: Remove triangles_iterator constructor

abstract type AbstractSSTen end

#TODO: Determine if we should make mutable in a later version
#TODO: Can make this parametric, but only supporting real floats currently.
struct COOTen <: AbstractSSTen

	cubical_dimension::Int
	order::Int
	unique_nnz::Int
	indices::Array{Int,2} #must be sorted lexigraphically for efficient Inner Products
	vals::Array{Float64,1}

	#edges::Vector{Tuple{Array{Int,1},Float64}}

	function COOTen(indices::Array{Int,2};nocheck::Bool=false)
		if nocheck
			unique_nnz, order = size(indices)
			return new(maximum(indices),order,unique_nnz,indices, ones(unique_nnz))
	    else
			cubical_dimension,unique_nnz,order,edges = COOTenVerifier(indices)
		    return new(cubical_dimension,order,unique_nnz,sort(indices,dims=1),ones(unique_nnz))
		end
	end

	function COOTen(indices::Array{Int,2},values::Array{N,1};
		            nocheck::Bool=false) where N <: AbstractFloat
		if nocheck
			unique_nnz, order = size(indices)
			return new(maximum(indices),order,unique_nnz,indices, values)
	    else
			cubical_dimension,unique_nnz,order,indices,values = COOTenVerifier(indices,values)
			perm = sort(1:unique_nnz, by=i->indices[i,:])
	   	    return new(cubical_dimension,order,unique_nnz,indices[perm,:],values[perm])
		end
	end

	function COOTen(indices::Array{Int,2},values::Array{N,1},n::Int;
		            nocheck::Bool=false) where N <: AbstractFloat
		unique_nnz, order = size(indices)
		if nocheck
			return new(n,order,unique_nnz,indices,values)
	    else
    		perm = sort(1:unique_nnz, by=i->indices[i,:])
			_,unique_nnz,order,indices,values = COOTenVerifier(indices,values,n)
			return new(n,order,unique_nnz,indices[perm,:],values[perm])
		end
	end

	COOTen() = error("empty hyper-edge list.")
end

#iterators for COOTen
Base.iterate(A::COOTen, state=1) =
    state > length(A) ? nothing : ((A.indices[state,:],A.vals[state]),state +1)
Base.length(A::COOTen) = A.unique_nnz

#TODO: rename to DICTen
mutable struct SSSTensor <: AbstractSSTen
    edges::Dict{Array{Int,1},AbstractFloat}
    cubical_dimension::Int

    function SSSTensor(e::Union{Array{Tuple{Array{Int,1},N},1},
							    Dict{Array{Int64,1},N}},
					   n::Int;nocheck::Bool=false) where {N <: AbstractFloat}
	     if nocheck
    		 new(reduce_edges(e),n)
		 else

			 if SSSTensor_verifier(e,n)
				 new(reduce_edges(e),n)
			 else
				 error("Cubical dimension $(n) too small.")
			 end
		 end
       # Need to adjust the error message on this constructor
	end

	#Edge List constructor
	function SSSTensor(e::Union{Array{Tuple{Array{Int,1},N},1},
							    Dict{Array{Int64,1},N}};
					   nocheck::Bool=false) where {N <: AbstractFloat}
		n = SSSTensor_verifier(e,nocheck)
		new(reduce_edges(e),n)
	end

	SSSTensor() = error("empty hyper-edge list.")

	#Dense tensor constructor
	function SSSTensor(A::Array{N,k},n::Int=typemax(Int)) where {N <: AbstractFloat,k}

		if any(size(A) .> n)
		    error("input Tensor size ",size(A),
			      " has dimension too large for input cubical dimension ",n)
		else
			edge_dict, n = SSSTensor_from_Array(A)
			new(edge_dict,n)
		end
	end

	#function SSSTensor(path_file::String)
	#TODO: load .ssten file
end

#write iterators for COOTen
Base.iterate(A::SSSTensor, state=1) = Base.iterate(A.edges,state)
Base.length(A::SSSTensor) = Base.length(A.edges)

# Include the files from the other functions
include("Helpers.jl") # makes helpers accessible to other files
include("CoreOperations.jl")
include("Contraction.jl")
include("HyperGraphAlgos.jl")
include("NumericalRoutines.jl")

@generated function SSSTensor_from_Array(A::Array{N,k}) where {N<:AbstractFloat,k,p}
    quote

        shape = size(A)
        n = shape[1]
        for i in 2:length(shape)
            if shape[i] != n
                error(string("input tensor of shape ",shape," is not cubical"))
            end
        end

        y = Dict{Array{Int64,1},N}()
        string_as_varname_function(@sprintf("i_%d",$k+1),n)
        @nloops $k i d-> 1:i_{d+1} begin
            indices = @ntuple $k i
            add_perm!(y,A,$k,indices)
        end
        return y,n
    end
end

#=------------------------------------------------------------------------------
						          Input Verifiers
--------------------------------------------------------------------------------
						          SSTen Verifier
------------------------------------------------------------------------------=#
"""-----------------------------------------------------------------------------
    SSSTensor_verifier(edges,n)

  This function takes a list of edges and a cubical dimension and checks whether
or not the edges are appropriate for a a super symmetric cubical tensor of n
dimensions.

Inputs:
-------
* edges - (Array{Tuple{Array{Int,1},Float}}):

  An array which contains Tuples of index arrays and Tupleed edge values
  associated. The indices must be sorted.
* n  - (Int):

  An iteger indicating the desired dimension of the cubical tensor.

Outputs:
--------
* is_valid - (Bool)

    An integer indicating whether or not the edges are appropriate for the
    tensor specified.
-----------------------------------------------------------------------------"""
function SSSTensor_verifier(edges::Union{Array{Tuple{Array{Int,1},N},1},
                                         Dict{Array{Int64,1},N}},
										 n::Int) where N <: AbstractFloat
    max_index = SSSTensor_verifier(edges,false)
    return max_index <= n
end

#UNTESTED
"""-----------------------------------------------------------------------------
    SSSTensor_verifier(edges,nocheck)

  This function takes in a list of edges and checks whether or not the edges are
appropriate for a super symmetric tensor. Finds the largest index and over all
the edges and sets that as the cubical dimension of the tensor. Used as a helper
function for the SSSTensor constructors.

Input:
------
* edges - (Dict{Array{Int64,1},AbstractFloat}):

    An array which contains Tuples of index arrays and paired edge values
    associated. The indices must be sorted.

* nocheck - (optional Bool):

Output:
-------
* max_index - (Int):

    An integer indicating the maximum index, returns 0 if an edge is found not
    to be sorted.
-----------------------------------------------------------------------------"""
function SSSTensor_verifier(edges::Dict{Array{Int64,1},N},
                            nocheck::Bool=false) where N <: AbstractFloat
    max_index = -Inf
    order = -1
    for (indices,_) in edges
        if order == -1
            order = length(indices)
        end

		if !nocheck  #assure input meets invariant standards
			if length(indices) != order
                error("hyperedge $(indices) must be order $(order)")
            end
			if !issorted(indices)
			  error("$(indices) must be sorted ascendingly")
			end
			if any(x -> x < 1,indices)
			  error("$(indices) has an index < 1")
			end
		end

		if indices[end] > max_index
			max_index = indices[end]
		end
    end

    return max_index
end
#=
"""-----------------------------------------------------------------------------
    SSSTensor_verifier(edges)

  This function takes in a dictionary of edges and a cubical dimension and checks
whether or not the edges are appropriate for a super symmetric tensor with
dimension n. Used as a helper function for the SSSTensor constructors.

Input:
------
* edges - (Array{Tuple{Array{Int,1},Float}}):

    An array which contains Tuples of index arrays and paired edge values
    associated. The indices must be sorted.

Output:
-------
* max_index - (Int):

    An integer indicating the maximum index, returns 0 if an edge is found not
    to be sorted.
-----------------------------------------------------------------------------"""
=#
function SSSTensor_verifier(edges::Array{Tuple{Array{Int,1},N},1},
                            nocheck::Bool=false) where N <: AbstractFloat
  max_index = -Inf
  order = -1
  for (indices,_) in edges
    if order == -1
	  order = length(indices)
    else
	  if length(indices) != order
	    error(string("hyperedge ",indices," must be order ",order))
	  end
    end

    if !issorted(indices)
	  error(string(indices, " must be sorted ascendingly"))
    end
    if any(x -> x < 1,indices)
	  error(string(indices," has an index < 1"))
    end
    if indices[end] > max_index
	  max_index = indices[end]
    end
  end

  return max_index
end

#=------------------------------------------------------------------------------
						          COOTen Verifiers
------------------------------------------------------------------------------=#

"""-----------------------------------------------------------------------------
    COOTen_verifier(indices, values)

  This function takes in a 2-D array and checks to see that the indices in each
row are sorted, that the sizes of the passed in values are correct and if the
invariants are preserved, returns the order, cubical_dimension, and unique
non-zero count back for the constructor.

Input:
------
* indices - (Array{Int64,2}):

    Array storing all the indices of the tensor to be formed. Rows must be
    sorted in increasing order.

* values - (optional Array{Float64,1}):

    Array storing the non-zero values of the tensor to be formed. Must have the
    same number of rows as indices as each row of each correspond to the same
    hyperedge. If none are passed in, one is constructed containing all ones.

* n - (optional Int)

    The desired cubical dimension, allows for the construction of tensors with
    zero subtensors. Throws an error if not the largest
Output:
-------

* cubical_dimension

  The dimension of each mode for the cubical tensor formed.

* unique_nnz

  The number of unique non-zeros in the tensor up to symmetry.

* order

  The number of modes of the tensor.

* indices - (Array{Int,2})

  A static matrix made from the indices passed in.

* values - (Array{AbstractFloat,1})

  A static vector made from the values passed in.
-----------------------------------------------------------------------------"""
function COOTenVerifier(indices::Array{Int,2},
    					values::Union{Array{Float64,1},Nothing}=nothing,
                        n::Union{Int,Nothing}=nothing)

    rows, order = size(indices)
	if values === nothing
		values = ones(rows)
	end

	if n === nothing
		n = Inf
	end
	unique_nnz = length(values)

	if rows != unique_nnz
		error("Indices and values have mis-matched sizes. indices has $(rows)",
		      " rows and values has $(unique_nnz) entries. These must be the ",
			  "same.")
	end

	cubical_dimension = 0

    for i in 1:rows
		#check for sorted rows
		if !issorted(indices[i,:])
			error("row $(i) has unsorted indices: $(indices[i,:]). Each row ",
			      "must be sorted in increasing order.")
		end

		#check for valid indices
		if indices[i,1] < 1
			error("row $(i) has index less than 1 in entry $(indices[i,:]), ",
			      "Julia is indexed by 1.")
		elseif indices[i,order] > cubical_dimension

			if cubical_dimension > n
				error("indices contain vertex index,$(indices[i,order]), ",
				      "greater than specified cubical dimension, n:$(n), ",
					  "specified.")
			end
			cubical_dimension = indices[i,order]
		end
	end

	return cubical_dimension,unique_nnz,order,indices, values
	#       [((indices[i,:]),values[i]) for i in 1:unique_nnz]

end

end #module end