text:
# __precompile__(true)  #  SEMBLE N'AVOIR AUCUN EFFET !
module Args

# import Base: tac

export get_syntaxe, get_usage, parse_commandline
export to_s_dict, get_args, show_args, get, set
export actions_desc, action, APPDIR
export actions, x_solvers, t_algos, presorts, informats, outformats
export reset_seed

global const APPDIR = dirname(dirname(@__FILE__()))
global args=nothing
global settings=nothing

include("utils/abbrev_dict.jl")

using ArgParse
using DataStructures # risque de pollution du scope
using Random
using Statistics

# import DataStructures: OrderedDict # résoudrait le warning de conflit potentiel

# using Globals

"""
Accesseur à un paramètre de la ligne de commande.
Sans argument, cette méthode retourne le dict qui peut être modifié par
l'appelant
"""
function get(k="")
    global args
    if args === nothing
        error("Dans args.jl/get : args n'est pas encore défini !")
    end
    # # On s'assure que k est une String, même si un symbole était fourni.
    # k = String(k)
    # On s'assure que k est un Symbol, même si une String était fournie.
    k = Symbol(k)
    if k == ""
        return args
    else
        return args[k]
    end
end

"""
Permet de modifier dynamiquement un paramètre du programme global.
e.g
```
old_level = Args.get(:level)
Args.set(:level, 4) # on passe en très verbose provisoirement
my_function_to_debug()
Args.set(:level, old_level) # on revient à la verbosité iniiale
```
"""
function set(k, v)
    global args
    if args === nothing
        error("Dans args.jl/set : args n'est pas encore défini !")
    end
    k = Symbol(k)
    args[k] = v
end

function actions_desc()
    OrderedDict{String, String}(
        "annealing" => "résolution par recuit simulé (sim. annealing)",
        "recuit"    => "résolution par recuit simulé (sim. annealing)",

        "explore" =>
        "résolution par exploration pure",

        "carlo" =>
        "résolution par méthode de Monte-Carlo",

        "descent" =>
        "résolution par steepest descent",

        "greedy" =>
        "résolution par algorithme glouton",

        "mip" =>
        "résolution par solver MIP (CPLEX, GLPK, CLP, ...)",

        "dmip" =>
        "résolution par solver DiscretMip (CPLEX, GLPK, CLP, ...)
         (générique mais lent, avec une variable binaire par date d'atterrissage)",

        "timing" =>
        "résoud le timing d'une solution à partir d'une instance et d'une liste
        de noms.
        --infile : nom du fichier d'instance traitée
        --planes : liste de noms d'avions imposant l'ordre d'une solution.
        exemple :  \"[3,4,5,6,8,9,7,1,10,2]\"
        ou bien :  \" 3 4 5 6 8 9 7 1 10 2 \"",

        "validate" =>
        "vérifie une solution à partir d'une instance et d'un solution complète.
        --infile : nom du fichier d'instance
        --solfile : le fichier solution à valider (au format alp).
        ",

        "stats" =>
        "statistiques sur l'instance (dont test de l'inégalité triangulaire, ...)",

        "test" =>
        "exécute la procédure de test (pour mise au point d'une nouvelle action) ",

        "help" =>
        "Syntaxe de ce programme et description des actions disponibles"
   )
end

# function actions()
#     actions = collect(keys(actions_desc()))
# end
actions() = collect(keys(actions_desc()))
x_solvers() = ["cplex", "cbc", "clp", "glpk"]
t_algos() = ["dp", "lp", "earliest"]
presorts() = ["none", "target", "rtarget", "lb", "rlb", "ub", "rub", "mean", "bary", "shuffle"]
# informats() = ["AUTO", "alp", "alpx", "orlib", "ampl"]
# outformats() = ["AUTO", "alp", "alpx", "ampl"]
# reductors() = ["mean", "median"]

function get_syntaxe(error_msg="")
    exe = basename(Base.PROGRAM_FILE) # myapp.jl
    # exe = Base.PROGRAM_FILE # /homa/uma/diam/.../myapp.jl

    if error_msg != ""
        error_msg = "\n    $(error_msg)\n"
    end
    txt = "$(error_msg)
    Syntaxe : $(exe) [action] [--key1 val1] [--flag3] [--key3 val3] ...]\n
    Actions autorisées : $(join(actions(), ","))
    Aide en ligne : $(exe) --help    (ou -h)
    "
end

function get_actions_dest()
    io = IOBuffer()
    println(io, """
    Une action est toute abréviation non ambiguë de la liste suivante :
    """)
    for (key,val) in actions_desc()
        # println(io, string(key, "\n  $val\n"))
        println(io, key, "\n  ", val)
    end
    String(take!(io))
end

function get_usage(error_msg="")
    global settings

    io = IOBuffer()

    println(io,"="^70)
    println(io, "### Liste des actions disponibles")
    println(io, get_actions_dest())
    println(io,"="^70)
    println(io, "### Résumé de la syntaxe")
    println(io, usage_string(settings))
    println(io, """
    Utiliser -h ou --help pour la liste détaillée des options
    """)

    String(take!(io))
end

function to_s_dict(args)
    io = IOBuffer()
    for (k,v) in args
        # @show typeof(k)
        if typeof(v) <: Dict
            println(io, " $(rpad(k, 10))  => Dict:")
            # for (local,local) in variable_parent => FONCTIONNE !
            for (k2,v2) in v
                # println(io, "     $k2 => $v2")
                println(io, "     $(rpad(k2, 10))  =>  $v2")
            end
        else
            # println("$k => $v")
            # @show k, v
            if v === nothing
                println(io, " $(rpad(k, 10))  =>  nothing")
            else
                println(io, " $(rpad(k, 10))  =>  $v")
            end
        end
    end
    String(take!(io))
end
function show_dict(args)
    println(to_s_dict(args))
end

function show_args()
    println("Args : état des paramètres")
    global args
    show_dict(args)
end

"""
réinitialise la seed du générateur aléatoire avec la valeur
passée en paramètre.
Cette méthode peut (doit) être appelée par exemple avant chaque lancement
d'un test de performance multiple pour rendre le résultat déterministe.

ATTENTION : si ne paramètre vaut 0, le germe n'est jamais le même.
Cette méthode est donc surtout utile si l'on vaut un comportement déterministe
à différent endroit du programme.

- si seed == -1 => on utilise la valeur de Args.get[:seed]
- si seed == 0 => aléatoire (idem si Args.get[:seed] == 0 )
- si seed > 0 : impose ce germe au générateur aléatoire

Exemple :
```reset_seed()```
=> réinitialise le germe avec la valeur passé en paramètre de l'application

```reset_seed(seed=Args.get(:seed)+i)```
=> initialise le germe avec une valeur dépendante de l'itération i
"""
function reset_seed(; seed=-1, level=2)
    if seed == -1
        seed = Args.get(:seed)
    end
    if seed == 0
        used_seed = time_ns()
    else
        used_seed = seed
    end
    if level >= 2
        println("Reset random generator with seed=$used_seed")
    end
    Random.seed!(used_seed)
end

function parse_commandline(inargs= ["test"])
    global args, settings

    # Si l'argument est une String : on le casse en tableau de String
    # Pratique surtout pour les tests mais inutile en fonctionnement normal)
    # e.g parse_commandline("val -p p3,p1,p6,...,p10")
    if isa(inargs, String)
        inargs = split(inargs)
    end

    settings = ArgParseSettings()
    # settings.prog = "seqata.jl"
    settings.description = "Validateur/Solver pour le projet Sequata"
    settings.version = "Seqata-0.0.7"
    settings.add_version = true
    settings.allow_ambiguous_opts = true
    settings.autofix_names = true
    @add_arg_table settings begin
        "action"
        help = "action à exécuter.
                Actions autorisées : [$(join(actions(), ", "))]"
        # help = "action à exécuter INCOMPLET....."
        # required = true
        default = "none"

        "--infile", "-i"
        help = "Nom du fichier d'instance."
        # required = true  # 20/06/2018 inhibé pour permettre action help
        default = "NO_INFILE"

        "--outdir", "-d"
        help = "Répertoire d'enregistrement des fichiers de sortie
                (pour les solutions ou la génération d'instance).
                Le répertoire doit exister, sinon \".\"est utilisé."
        default = "_tmp"
        # range_tester = (v->isdir(v))

        "--solfile", "-s" #, "--sol", "--solution", "-S"
        help = "Non du fichier d'une solution initiale ou à valider."
        default = "NO_SOLFILE"

        "--planes", "-p"
        help = "Liste de noms d'avions imposant l'ordre d'une solution. " *
               "e.g \"[3,1,4,2...]\" ou \"3,1,4,2...\" ou \"3 1 4 2...\""

        "--force", "-f"
        help = "Force l'écrasement d'un fichier existant"
        action = :store_true

        "--itermax", "-n"
        help = "Nombre d'itération (sa fonction dépend de l'action). "
        arg_type = Int
        default = 0    # 0 pour AUTOMATIC

        "--nb-no-improv-max", "--nim"
        help = "Nombre de non-amélioration maxi (pour arrêt du tabou)."
        arg_type = Int
        default = 0    # 0 pour AUTOMATIC

        "--presort"
        help = "tri initial des avions selon la date d'atterrissage cible. "*
               "Tri possible :  [$(join(presorts(), ", "))] "
        default = "target"
        range_tester = (v->v in presorts())

        "--external-mip-solver", "-x"
        help = "Solver externe : [$(join(x_solvers(), ", "))]"
        default = "glpk"
        range_tester = (v->v in x_solvers())

        "--timing-algo-solver", "-t"
        help = "Algo de résolution pour le sous-problème de timing."*
               "(possible : [$(join(t_algos(), ", "))])"
        # default = "earliest"
        default = "lp"
        range_tester = (v->v in t_algos())

        "--cost-precision", "--prc"
        help = "nombre de décimales pour l'affichage des coûts." *
               "Utilisé aussi en interne pour nomaliser le coûts de chaque " *
               "avions en chaque date possible."
        arg_type = Int
        range_tester = (v->v>=0) # peut être nul (mais non testé !)
        # default = 5
        default = 7 # modif le 30/10/2018

        "--level", "--loglevel", "-L"
        help = "Niveau de verbosité (pour debug)"
        arg_type = Int
        default = 0

        "--test", "-T"
        help = "Mode test (pour mise au point d'une variante du code)"
        action = :store_true

        "--loop"
        help = "Loop over instances"
        action = :store_true

        "--plot"
        help = "Plot current solution cost over time"
        action = :store_true

        "--save-costs"
        help = "Saves solution cost over time"
        action = :store_true

        "--seed"
        help = "Graine pour le générateur aléatoire"
        arg_type = Int
        default = 0
        range_tester = v->begin
            if v < 0
                println("\n\nLe germe aléatoire doit être positif ou nul\n")
                return false
            end
        return true
        end

        "files"
        nargs = '*'
        help = "arguments supplémentaires (e.g liste de fichiers...).
                L'utilisation est réservée à certaines actions"
        required = false
        # default = Any["popo"] # Vector vide de Any
        default = Any[] # Vector vide de Any

    end

    if false
        for arg in ARGS
            println("ARGS AVANT extraction action arg=$(arg)")
        end
    end

    global args
    # args = ArgParse.parse_args(settings)
    args = ArgParse.parse_args(inargs, settings, as_symbols=true)
    args[:external_mip_solver] = Symbol(args[:external_mip_solver])
    args[:timing_algo_solver]  = Symbol(args[:timing_algo_solver])
    args[:presort]             = Symbol(args[:presort])
    if false
        show_args()
    end

    #
    # ON FAIT QUELQUES VERIFICATIONS ET POST-TRAITEMENTS
    #

    # on extrait la vraie action à partir de son abréviation
    #
    if args[:action] != "none"
        abbrev = args[:action]
        abbrevs = abbrevsdict(actions()) # voir fichier abbrev.jl
        # @show abbrev
        if !haskey(abbrevs, abbrev)
            msg =  get_syntaxe() *
                "\nErreur : action \"$(abbrev)\" inconnue (dans args.jl)"
            throw(ArgumentError(msg))
        end
        # On récupère le nom complet de l'action
        # (si on a passé l'action "val" on remplace par "validate")
        args[:action] = abbrevs[abbrev]
    end
    # e.g. On transforme la chaine "validate" en symbole :validate
    args[:action] = Symbol(args[:action])

    # si action Help : on affiche la liste des actions dispos, puis on quitte
    #
    if args[:action] == :help
        println("="^70)
        println("### Valeur des paramètres après analyse :")
        show_args()
        println(get_usage())
        exit(0)
    end

    # Si on a pas passé --infile, on récupère le premier argument sans clé
    # et on le supprime de la liste
    #
    if args[:infile] == "NO_INFILE" && length(args[:files]) >= 1
        args[:infile] = popfirst!(args[:files])
    end

    # Vérification du fichier d'instance à traiter
    #
    if args[:infile] != "NO_INFILE" && !isfile(args[:infile])
        # println("ERREUR : fichier inexistant ou illisible : $(args[:infile])")
        # exit(1)
        msg = "ERREUR : fichier inexistant ou illisible : $(args[:infile])"
        throw(ArgumentError(msg))
    end

    # Si on a pas passé --solfile, on récupère le premier argument sans clé
    # et on le supprime de la liste
    #
    if args[:solfile] == "NO_SOLFILE" && length(args[:files]) >= 1
        args[:solfile] = popfirst!(args[:files])
    end

    # si action est timing : on doit extraire la liste des avions à valider
    #
    if args[:action] == :timing
        # planes = matchall(r"[\w]+", args[:planes])
        planes = collect((m.match for m = eachmatch(r"[\w]+", args[:planes]) ))
        if length(planes) < 3
            msg = "ERREUR : pas assez d'avions : $(ength(planes)) (au moins trois !)"
            throw(ArgumentError(msg))
        end
        args[:planes] = planes
        println("planes=",join(args[:planes]))
    end

    # si action est validate : on extrait le nom du fichier solution à valider
    #
    if args[:action] == :validate
        # Le paramètre solfile doit exister (via option --solfile ou sans clé)
        if args[:solfile] != "NO_SOLFILE" && !isfile(args[:solfile])
            msg = "ERREUR : fichier inexistant ou illisible : $(args[:solfile])"
            throw(ArgumentError(msg))
        end
    end

    # Vérification du répertoire de sortie
    #
    # @show Args.args
    # show_args()
    if !isdir(args[:outdir])
        # println("ERREUR : répertoire inexistant ou illisible : $(args[:outdir])")
        # exit(1)
        if Args.args[:level] >= 3
            println("L3 WARNING : répertoire inexistant ou illisible : $(args[:outdir])")
            println("L3 => on utilisera le répertoire courant : (\".\")")
        end
        args[:outdir] = "."
    end

    # L'action :recuit n'est qu'un alias pour l'action :annealing
    if args[:action] == :recuit
        args[:action] = :annealing
    end

    Args.reset_seed(level=1)

    return args
end

# Args.parse_commandline()

end # module


text:
#Demo 13: A simple scatter plot
#-------------------------------------------------------------------------------
using InspectDR
using Colors
using Random

red = RGB24(1, 0, 0)
green = RGB24(0, 1, 0)
blue = RGB24(0, 0, 1)

#==Input
===============================================================================#
x = collect(0:100) #Must vectorize using collect - ranges not yet supported
Random.seed!(11) #Reseed
y = randn(length(x))

#==Generate plot object
===============================================================================#
plot = InspectDR.Plot2D(:lin, :lin,
	title = "Very Minimal Example: Scatter Plot",
	xlabel = "X-Values",
	ylabels = ["Y-Values"]
)
plot.layout[:enable_legend] = true
plot.layout[:halloc_legend] = 150 #Default: 100

#Set grid:
graph = plot.strips[1]
graph.grid = InspectDR.GridRect(vmajor=true, vminor=true, hmajor=true)

wfrm = add(plot, x, y, id="Random Sample")
	#No line on a scatter plot; width used for glyph:
	wfrm.line = line(width=3, style=:none)
	wfrm.glyph = glyph(shape=:o, size=10, color=red, fillcolor=blue)

#==Render plot
===============================================================================#
#Show Gtk GUI
gplot = display(InspectDR.GtkDisplay(), plot)

#Don't need to show GUI if simple plot image is desired:
InspectDR.write_png("export_simpleplot.png", plot)
#InspectDR.write_svg("export_simpleplot.svg", plot)
#InspectDR.write_eps("export_simpleplot.eps", plot)
#InspectDR.write_pdf("export_simpleplot.pdf", plot)

:DONE

ext: jl


text:
# ### m8.4t

using Pkg, DrWatson

@quickactivate "StatisticalRethinkingTuring"
using Turing
using StatisticalRethinking
Turing.setprogress!(false)

# Can't really set a Uniform[-Inf,Inf] on σ

# Turing model
@model ppl8_4(y) = begin
    α₁ ~ Uniform(-3000, 1000)
    α₂ ~ Uniform(-1000, 3000)
    σ ~ truncated(Cauchy(0,1), 0, Inf)

    y .~ Normal(α₁ + α₂, σ)
end

# Observations

y = rand(Normal(0,1), 100);

# Sample

m8_4t = ppl8_4(y)
nchains = 4; sampler = NUTS(0.65); nsamples=2000
chns8_4t = mapreduce(c -> sample(m8_4t, sampler, nsamples), chainscat, 1:nchains)

# Results rethinking

m8_4rethinking = "
         mean      sd     5.5%   94.5% n_eff Rhat
 a1    -861.15 558.17 -1841.89  -31.04     7 1.43
 a2     861.26 558.17    31.31 1842.00     7 1.43
 sigma    0.97   0.07     0.89    1.09     9 1.17
";

# Describe the posterior samples

chns |> display

# End of `08/m8.4t.jl`


text:
struct Int64Ctx
    value::Int64
end
@context Int64Ctx

struct Int64TypeVarCtx{V} end
@context Int64TypeVarCtx

@testset "outline contextual" begin

    f(r) = r+3

    @contextualized function f(r) @with {offset::Int64Ctx}
        x = r + offset.value
        return x
    end

    @test f(2) == 5
    @test (@withctx Int64Ctx(5) f(2)) == 7
end

@testset "outline contextual (inline)" begin

    f(r) = r+3

    @inline @contextualized function f(r) @with {offset::Int64Ctx}
        x = r + offset.value
        return x
    end

    @test f(2) == 5
    @test (@withctx Int64Ctx(5) f(2)) == 7
end

@testset "outline contextual (short)" begin

    f(r) = r+3

    @contextualized f(r) = (@with {offset::Int64Ctx}; r + offset.value)

    @test f(2) == 5
    @test (@withctx Int64Ctx(5) f(2)) == 7
end

@testset "outline contextual (short, inline)" begin

    f(r) = r+3

    @inline @contextualized f(r) = (@with {offset::Int64Ctx}; r + offset.value)

    @test f(2) == 5
    @test (@withctx Int64Ctx(5) f(2)) == 7
end

@testset "outline contextual without context assignment" begin
    g() = 0

    @contextualized function g() @with {_::Int64Ctx}
        1
    end

    @test g() == 0
    @test (@withctx Int64Ctx(5) g()) == 1
end

@testset "outline contextual without context assignment (short)" begin
    g() = 0

    @contextualized g() = (@with {_::Int64Ctx}; 1)

    @test g() == 0
    @test (@withctx Int64Ctx(5) g()) == 1
end

@testset "outline contextual with access to type" begin
    g() = 0
    @contextualized function g() where {V}
        @with {_ :: Int64TypeVarCtx{V}}
        V
    end

    @test g() == 0
    @test (@withctx Int64TypeVarCtx{3}() g()) == 3
end

@testset "outline contextual with access to type (short)" begin
    g() = 0
    @contextualized g() where {V} = (@with {_ :: Int64TypeVarCtx{V}}; V)

    @test g() == 0
    @test (@withctx Int64TypeVarCtx{3}() g()) == 3
end
abstract type OutlineTypeA end
struct OutlineTypeB <: OutlineTypeA end

@testset "outline contextual subtype" begin
    g() = 0

    @contextualized function g() @with {ctx<:OutlineTypeA}
        1
    end

    @test g() == 0
    @test (@withctx OutlineTypeB() g()) == 1
end

@testset "outline contextual subtype (short)" begin
    g() = 0

    @contextualized g() = (@with {ctx<:OutlineTypeA}; 1)

    @test g() == 0
    @test (@withctx OutlineTypeB() g()) == 1
end


text:
using Test

ENV["TESTING"] = "TRUE"
const GROUP = get(ENV, "GROUP", "test")

OUTER_GROUP = first(split(GROUP, ' '))

# Run the tests.
if OUTER_GROUP == "test" || OUTER_GROUP == "all"

    # Determines which group of tests should be run.
    group_info = split(GROUP, ' ')
    TEST_GROUP = length(group_info) == 1 ? "all" : group_info[2]

    using AbstractGPs
    using BlockDiagonals
    using ChainRulesCore
    using FillArrays
    using FiniteDifferences
    using LinearAlgebra
    using KernelFunctions
    using Random
    using StaticArrays
    using StructArrays
    using TemporalGPs

    using Zygote

    using FiniteDifferences: rand_tangent
    using AbstractGPs: var
    using TemporalGPs: AbstractLGSSM, _filter, NoContext
    using Zygote: Context, _pullback

    include("test_util.jl")

    @show TEST_GROUP GROUP

    @testset "TemporalGPs.jl" begin

        if TEST_GROUP == "util" || GROUP == "all"
            println("util:")
            @testset "util" begin
                include(joinpath("util", "harmonise.jl"))
                include(joinpath("util", "scan.jl"))
                include(joinpath("util", "zygote_friendly_map.jl"))
                include(joinpath("util", "zygote_rules.jl"))
                include(joinpath("util", "gaussian.jl"))
                include(joinpath("util", "mul.jl"))
                include(joinpath("util", "regular_data.jl"))
            end
        end

        if TEST_GROUP ∈ ["models", "models-lgssm", "gp", "space_time"] || GROUP == "all"
            include(joinpath("models", "model_test_utils.jl"))
            include(joinpath("models", "test_model_test_utils.jl"))
        end

        if TEST_GROUP == "models" || GROUP == "all"
            println("models:")
            @testset "models" begin
                include(joinpath("models", "linear_gaussian_conditionals.jl"))
                include(joinpath("models", "gauss_markov_model.jl"))
                include(joinpath("models", "missings.jl"))
            end
        end

        if TEST_GROUP == "models-lgssm" || GROUP == "all"
            println("models (lgssm):")
            @testset "models (lgssm)" begin
                include(joinpath("models", "lgssm.jl"))
            end
        end

        if TEST_GROUP == "gp" || GROUP == "all"
            println("gp:")
            @testset "gp" begin
                include(joinpath("gp", "lti_sde.jl"))
                include(joinpath("gp", "posterior_lti_sde.jl"))
            end
        end

        if TEST_GROUP == "space_time" || GROUP == "all"
            println("space_time:")
            @testset "space_time" begin
                include(joinpath("space_time", "rectilinear_grid.jl"))
                include(joinpath("space_time", "regular_in_time.jl"))
                include(joinpath("space_time", "separable_kernel.jl"))
                include(joinpath("space_time", "to_gauss_markov.jl"))
                include(joinpath("space_time", "pseudo_point.jl"))
            end
        end
    end
end

# Run the examples.
if GROUP == "examples"

    using Pkg

    Pkg.activate(joinpath("..", "examples"))
    Pkg.develop(path="..")
    Pkg.resolve()
    Pkg.instantiate()

    include(joinpath("..", "examples", "exact_time_inference.jl"))
    include(joinpath("..", "examples", "exact_time_learning.jl"))
    include(joinpath("..", "examples", "exact_space_time_inference.jl"))
    include(joinpath("..", "examples", "exact_space_time_learning.jl"))
    include(joinpath("..", "examples", "approx_space_time_inference.jl"))
    include(joinpath("..", "examples", "approx_space_time_learning.jl"))
    include(joinpath("..", "examples", "augmented_inference.jl"))
end


text:
endpoint(::GitHubAPI, ::typeof(is_member), org::AStr, user::AStr) =
    Endpoint(:GET, "/orgs/$org/members/$user"; allow_404=true)
postprocessor(::GitHubAPI, ::typeof(is_member)) = DoSomething(ismemberorcollaborator)
into(::GitHubAPI, ::typeof(is_member)) = Bool


text:
using Lumira
using Gadfly
using Base.Test

plot = plot([sin, cos], 0, 25)

setToWriteToFile(false)

exportToSVG(plot)

setId("trigo")
setName("HelloWorld")

@test getId() == "trigo"
@test getName() == "HelloWorld"
@test createTemplate("dummyPath") == true
@test createChartCode("dummyPath") == true
@test createLumiraExtension() == true


text:
using FileIO
using LasIO
using Mmap
using BenchmarkTools

workdir = dirname(@__FILE__)
# source: http://www.liblas.org/samples/
filename = "libLAS_1.2.las" # point format 0
testfile = joinpath(workdir, filename)
writefile = joinpath(workdir, "libLAS_1.2-out.las")

function read_original(io::IO, ::Type{LasPoint0})
    x = read(io, Int32)
    y = read(io, Int32)
    z = read(io, Int32)
    intensity = read(io, UInt16)
    flag_byte = read(io, UInt8)
    raw_classification = read(io, UInt8)
    scan_angle = read(io, Int8)
    user_data = read(io, UInt8)
    pt_src_id = read(io, UInt16)
    LasPoint0(
        x,
        y,
        z,
        intensity,
        flag_byte,
        raw_classification,
        scan_angle,
        user_data,
        pt_src_id
    )
end

function test_orignal()
    open(testfile) do s
        LasIO.skiplasf(s)
        header = read(s, LasHeader)

        n = header.records_count
        pointtype = pointformat(header)
        pointdata = Vector{pointtype}(undef, n)
        for i=1:n
            pointdata[i] = read_original(s, pointtype)
        end
        header, pointdata
    end
end

function test_new()
    open(testfile) do s
        LasIO.skiplasf(s)
        header = read(s, LasHeader)

        n = header.records_count
        pointtype = pointformat(header)
        pointdata = Vector{pointtype}(undef, n)
        for i=1:n
            pointdata[i] = read(s, pointtype)
        end
        header, pointdata
    end
end

function test_stream()
    open(testfile) do s
        LasIO.skiplasf(s)
        header = read(s, LasHeader)

        n = header.records_count
        pointtype = pointformat(header)

        pointsize = Int(header.data_record_length)
        pointbytes = Mmap.mmap(s, Vector{UInt8}, n*pointsize, position(s))
        pointdata = PointVector{pointtype}(pointbytes, pointsize)

        for i=1:n
            pointdata[i]
        end

        header, pointdata
    end
end

oh, op = test_orignal()
nh, np = test_new()
sh, sp = test_stream()

println("Original function")
@btime test_orignal()
println("New function using generated read")
@btime test_new()
println("New function using streaming")
@btime test_stream()

@show op, np, sp
@assert op[5] == np[5] == sp[5]


text:
const ResultsByTime = SortedDict{Dates.DateTime, DataFrames.DataFrame}
const FieldResultsByTime = Dict{OptimizationContainerKey, ResultsByTime}

abstract type OperationModelSimulationResults end
# Subtypes need to implement the following methods for SimulationProblemResults{T}
# - read_aux_variables_with_keys
# - read_variables_with_keys
# - read_parameters_with_keys
# - read_duals_with_keys
# - read_expressions_with_keys
# - load_results!

"""
Holds the results of a simulation problem for plotting or exporting.
"""
mutable struct SimulationProblemResults{T} <:
               IS.Results where {T <: OperationModelSimulationResults}
    problem::String
    base_power::Float64
    execution_path::String
    results_output_folder::String
    timestamps::StepRange{Dates.DateTime, Dates.Millisecond}
    results_timestamps::Vector{Dates.DateTime}
    values::T
    system::Union{Nothing, PSY.System}
    system_uuid::Base.UUID
    resolution::Dates.TimePeriod
    store::Union{Nothing, SimulationStore}
end

function SimulationProblemResults{T}(
    store::SimulationStore,
    model_name::AbstractString,
    problem_params::ModelStoreParams,
    sim_params::SimulationStoreParams,
    path,
    vals::T;
    results_output_path=nothing,
    system=nothing,
) where {T <: OperationModelSimulationResults}
    if results_output_path === nothing
        results_output_path = joinpath(path, "results")
    end

    time_steps = range(
        sim_params.initial_time,
        length=problem_params.num_executions * sim_params.num_steps,
        step=problem_params.interval,
    )
    return SimulationProblemResults{T}(
        model_name,
        problem_params.base_power,
        path,
        results_output_path,
        time_steps,
        Vector{Dates.DateTime}(),
        vals,
        system,
        problem_params.system_uuid,
        get_resolution(problem_params),
        store isa HdfSimulationStore ? nothing : store,
    )
end

get_model_name(res::SimulationProblemResults) = res.problem
get_system(res::SimulationProblemResults) = res.system
get_resolution(res::SimulationProblemResults) = res.resolution
get_execution_path(res::SimulationProblemResults) = res.execution_path
get_model_base_power(res::SimulationProblemResults) = res.base_power
IS.get_timestamp(result::SimulationProblemResults) = result.results_timestamps
get_interval(res::SimulationProblemResults) = res.timestamps.step
IS.get_base_power(result::SimulationProblemResults) = result.base_power

# Subtypes can override if they don't have these fields.
get_aux_variables(res::SimulationProblemResults) = res.values.aux_variables
get_duals(res::SimulationProblemResults) = res.values.duals
get_expressions(res::SimulationProblemResults) = res.values.expressions
get_parameters(res::SimulationProblemResults) = res.values.parameters
get_variables(res::SimulationProblemResults) = res.values.variables

"""
Return an array of variable names (strings) that are available for reads.
"""
list_variable_names(res::SimulationProblemResults) =
    encode_keys_as_strings(keys(get_variables(res)))

"""
Return an array of dual names (strings) that are available for reads.
"""
list_dual_names(res::SimulationProblemResults) =
    encode_keys_as_strings(keys(get_duals(res)))

"""
Return an array of parmater names (strings) that are available for reads.
"""
list_parameter_names(res::SimulationProblemResults) =
    encode_keys_as_strings(keys(get_parameters(res)))

"""
Return an array of auxillary variable names (strings) that are available for reads.
"""
list_aux_variable_names(res::SimulationProblemResults) =
    encode_keys_as_strings(keys(get_aux_variables(res)))

"""
Return an array of expression names (strings) that are available for reads.
"""
list_expression_names(res::SimulationProblemResults) =
    encode_keys_as_strings(keys(get_expressions(res)))

"""
Return an array of VariableKeys that are available for reads.
"""
list_variable_keys(res::SimulationProblemResults) = collect(keys(get_variables(res)))

"""
Return an array of ConstraintKeys that are available for reading duals.
"""
list_dual_keys(res::SimulationProblemResults) = collect(keys(get_duals(res)))

"""
Return an array of ParameterKeys that are available for reads.
"""
list_parameter_keys(res::SimulationProblemResults) = collect(keys(get_parameters(res)))

"""
Return an array of AuxVarKeys that are available for reads.
"""
list_aux_variable_keys(res::SimulationProblemResults) =
    collect(keys(get_aux_variables(res)))

"""
Return an array of ExpressionKeys that are available for reads.
"""
list_expression_keys(res::SimulationProblemResults) = collect(keys(get_expressions(res)))

"""
Return a reference to a StepRange of available timestamps.
"""
get_timestamps(result::SimulationProblemResults) = result.timestamps

"""
Return the system used for the problem. If the system hasn't already been deserialized or
set with [`set_system!`](@ref) then deserialize and store it.
"""
function get_system!(results::SimulationProblemResults)
    file = joinpath(
        results.execution_path,
        "problems",
        results.problem,
        make_system_filename(results.system_uuid),
    )
    results.system = PSY.System(file, time_series_read_only=true)
    return results.system
end

"""
Set the system in the results instance.

Throws InvalidValue if the system UUID is incorrect.

# Arguments

  - `results::SimulationProblemResults`: Results object
  - `system::AbstractString`: Path to the system json file

# Examples

```julia
julia > set_system!(res, "my_path/system_data.json")
```
"""
function set_system!(results::SimulationProblemResults, system::AbstractString)
    set_system!(results, System(system))
end

function set_system!(results::SimulationProblemResults, system::PSY.System)
    sys_uuid = IS.get_uuid(system)
    if sys_uuid != results.system_uuid
        throw(
            IS.InvalidValue(
                "System mismatch. $sys_uuid does not match the stored value of $(results.system_uuid)",
            ),
        )
    end

    results.system = system
    return
end

function _deserialize_key(
    ::Type{<:OptimizationContainerKey},
    results::SimulationProblemResults,
    name::AbstractString,
)
    !haskey(results.values.container_key_lookup, name) && error("$name is not stored")
    return results.values.container_key_lookup[name]
end

function _deserialize_key(
    ::Type{T},
    results::SimulationProblemResults,
    args...,
) where {T <: OptimizationContainerKey}
    return make_key(T, args...)
end

function _get_containers(x::SimulationProblemResults)
    return (
        get_aux_variables(x),
        get_duals(x),
        get_expressions(x),
        get_parameters(x),
        get_variables(x),
    )
end

function _validate_keys(existing_keys, container_keys)
    existing = Set(existing_keys)
    for key in container_keys
        if key ∉ existing
            @error "$key is not stored", existing_keys
            throw(IS.InvalidValue("$key is not stored"))
        end
    end
    return
end

"""
Return the final values for the requested variables for each time step for a problem.

Decision problem results are returned in a Dict{String, Dict{DateTime, DataFrame}}.

Emulation problem results are returned in a Dict{String, DataFrame}.

Limit the data sizes returned by specifying `initial_time` and `count` for decision problems
or `start_time` and `len` for emulation problems.

See also [load_results!](@ref) to preload data into memory.

# Arguments

  - `variables::Vector{Union{String, Tuple}}`: Variable name as a string or a Tuple with
    variable type and device type. If not provided then return all variables.
  - `initial_time::Dates.DateTime`: Initial time of the requested results. Decision problems
    only.
  - `count::Int`: Number of results. Decision problems only.
  - `start_time::Dates.DateTime`: Start time of the requested results. Emulation problems
    only.
  - `len::Int`: Number of rows in each DataFrame. Emulation problems only.

# Examples

```julia-repl
julia > variables_as_strings =
    ["ActivePowerVariable__ThermalStandard", "ActivePowerVariable__RenewableDispatch"]
julia > variables_as_types =
    [(ActivePowerVariable, ThermalStandard), (ActivePowerVariable, RenewableDispatch)]
julia > read_realized_variables(results, variables_as_strings)
julia > read_realized_variables(results, variables_as_types)
```
"""
function read_realized_variables(res::SimulationProblemResults; kwargs...)
    return read_realized_variables(res, collect(keys(get_variables(res))); kwargs...)
end

function read_realized_variables(res::SimulationProblemResults, variables; kwargs...)
    return read_realized_variables(res, [VariableKey(x...) for x in variables]; kwargs...)
end

function read_realized_variables(
    res::SimulationProblemResults,
    variables::Vector{<:AbstractString};
    kwargs...,
)
    return read_realized_variables(
        res,
        [_deserialize_key(VariableKey, res, x) for x in variables];
        kwargs...,
    )
end

function read_realized_variables(
    res::SimulationProblemResults,
    variables::Vector{<:OptimizationContainerKey};
    kwargs...,
)
    result_values = read_variables_with_keys(res, variables; kwargs...)
    return Dict(encode_key_as_string(k) => v for (k, v) in result_values)
end

"""
Return the final values for the requested variable for each time step for a problem.

Decision problem results are returned in a Dict{DateTime, DataFrame}.

Emulation problem results are returned in a DataFrame.

Limit the data sizes returned by specifying `initial_time` and `count` for decision problems
or `start_time` and `len` for emulation problems.

See also [load_results!](@ref) to preload data into memory.

# Arguments

  - `variable::Union{String, Tuple}`: Variable name as a string or a Tuple with
    variable type and device type.
  - `initial_time::Dates.DateTime`: Initial time of the requested results. Decision problems
    only.
  - `count::Int`: Number of results. Decision problems only.
  - `start_time::Dates.DateTime`: Start time of the requested results. Emulation problems
    only.
  - `len::Int`: Number of rows in each DataFrame. Emulation problems only.

# Examples

```julia-repl
julia > read_realized_variable(results, "ActivePowerVariable__ThermalStandard")
julia > read_realized_variable(results, (ActivePowerVariable, ThermalStandard))
```
"""
function read_realized_variable(
    res::SimulationProblemResults,
    variable::AbstractString;
    kwargs...,
)
    return first(
        values(
            read_realized_variables(
                res,
                [_deserialize_key(VariableKey, res, variable)];
                kwargs...,
            ),
        ),
    )
end

function read_realized_variable(res::SimulationProblemResults, variable...; kwargs...)
    return first(
        values(read_realized_variables(res, [VariableKey(variable...)]; kwargs...)),
    )
end

"""
Return the final values for the requested auxiliary variables for each time step for a problem.

Refer to [`read_realized_variables`](@ref) for help and examples.
"""
function read_realized_aux_variables(res::SimulationProblemResults; kwargs...)
    return read_realized_aux_variables(
        res,
        collect(keys(get_aux_variables(res)));
        kwargs...,
    )
end

function read_realized_aux_variables(
    res::SimulationProblemResults,
    aux_variables;
    kwargs...,
)
    return read_realized_aux_variables(
        res,
        [AuxVarKey(x...) for x in aux_variables];
        kwargs...,
    )
end

function read_realized_aux_variables(
    res::SimulationProblemResults,
    aux_variables::Vector{<:AbstractString};
    kwargs...,
)
    return read_realized_aux_variables(
        res,
        [_deserialize_key(AuxVarKey, res, x) for x in aux_variables];
        kwargs...,
    )
end

function read_realized_aux_variables(
    res::SimulationProblemResults,
    aux_variables::Vector{<:OptimizationContainerKey};
    kwargs...,
)
    result_values = read_aux_variables_with_keys(res, aux_variables; kwargs...)
    return Dict(encode_key_as_string(k) => v for (k, v) in result_values)
end

"""
Return the final values for the requested auxiliary variable for each time step for a problem.

Refer to [`read_realized_variable`](@ref) for help and examples.
"""
function read_realized_aux_variable(
    res::SimulationProblemResults,
    aux_variable::AbstractString;
    kwargs...,
)
    return first(
        values(
            read_realized_aux_variables(
                res,
                [_deserialize_key(AuxVarKey, res, aux_variable)];
                kwargs...,
            ),
        ),
    )
end

function read_realized_aux_variable(
    res::SimulationProblemResults,
    aux_variable...;
    kwargs...,
)
    return first(
        values(
            read_realized_aux_variables(res, [aux_variableKey(aux_variable...)]; kwargs...),
        ),
    )
end

"""
Return the final values for the requested parameters for each time step for a problem.

Refer to [`read_realized_variables`](@ref) for help and examples.
"""
function read_realized_parameters(res::SimulationProblemResults; kwargs...)
    return read_realized_parameters(res, collect(keys(get_parameters(res))); kwargs...)
end

function read_realized_parameters(res::SimulationProblemResults, parameters; kwargs...)
    return read_realized_parameters(
        res,
        [ParameterKey(x...) for x in parameters];
        kwargs...,
    )
end

function read_realized_parameters(
    res::SimulationProblemResults,
    parameters::Vector{<:AbstractString};
    kwargs...,
)
    return read_realized_parameters(
        res,
        [_deserialize_key(ParameterKey, res, x) for x in parameters];
        kwargs...,
    )
end

function read_realized_parameters(
    res::SimulationProblemResults,
    parameters::Vector{<:OptimizationContainerKey};
    kwargs...,
)
    result_values = read_parameters_with_keys(res, parameters; kwargs...)
    return Dict(encode_key_as_string(k) => v for (k, v) in result_values)
end

"""
Return the final values for the requested parameter for each time step for a problem.

Refer to [`read_realized_variable`](@ref) for help and examples.
"""
function read_realized_parameter(
    res::SimulationProblemResults,
    parameter::AbstractString;
    kwargs...,
)
    return first(
        values(
            read_realized_parameters(
                res,
                [_deserialize_key(ParameterKey, res, parameter)];
                kwargs...,
            ),
        ),
    )
end

function read_realized_parameter(res::SimulationProblemResults, parameter...; kwargs...)
    return first(
        values(read_realized_parameters(res, [ParameterKey(parameter...)]; kwargs...)),
    )
end

"""
Return the final values for the requested duals for each time step for a problem.

Refer to [`read_realized_variables`](@ref) for help and examples.
"""
function read_realized_duals(res::SimulationProblemResults; kwargs...)
    return read_realized_duals(res, collect(keys(get_duals(res))); kwargs...)
end

function read_realized_duals(res::SimulationProblemResults, duals; kwargs...)
    return read_realized_duals(res, [ConstraintKey(x...) for x in duals]; kwargs...)
end

function read_realized_duals(
    res::SimulationProblemResults,
    duals::Vector{<:AbstractString};
    kwargs...,
)
    return read_realized_duals(
        res,
        [_deserialize_key(ConstraintKey, res, x) for x in duals];
        kwargs...,
    )
end

function read_realized_duals(
    res::SimulationProblemResults,
    duals::Vector{<:OptimizationContainerKey};
    kwargs...,
)
    result_values = read_duals_with_keys(res, duals; kwargs...)
    return Dict(encode_key_as_string(k) => v for (k, v) in result_values)
end

"""
Return the final values for the requested dual for each time step for a problem.

Refer to [`read_realized_variable`](@ref) for help and examples.
"""
function read_realized_dual(res::SimulationProblemResults, dual::AbstractString; kwargs...)
    return first(
        values(
            read_realized_duals(
                res,
                [_deserialize_key(ConstraintKey, res, dual)];
                kwargs...,
            ),
        ),
    )
end

function read_realized_dual(res::SimulationProblemResults, dual...; kwargs...)
    return first(values(read_realized_duals(res, [ConstraintKey(dual...)]; kwargs...)))
end

"""
Return the final values for the requested expressions for each time step for a problem.

Refer to [`read_realized_variables`](@ref) for help and examples.
"""
function read_realized_expressions(res::SimulationProblemResults; kwargs...)
    return read_realized_expressions(res, collect(keys(get_expressions(res))); kwargs...)
end

function read_realized_expressions(res::SimulationProblemResults, expressions; kwargs...)
    return read_realized_expressions(
        res,
        [ExpressionKey(x...) for x in expressions];
        kwargs...,
    )
end

function read_realized_expressions(
    res::SimulationProblemResults,
    expressions::Vector{<:AbstractString};
    kwargs...,
)
    return read_realized_expressions(
        res,
        [_deserialize_key(ExpressionKey, res, x) for x in expressions];
        kwargs...,
    )
end

function read_realized_expressions(
    res::SimulationProblemResults,
    expressions::Vector{<:OptimizationContainerKey};
    kwargs...,
)
    result_values = read_expressions_with_keys(res, expressions; kwargs...)
    return Dict(encode_key_as_string(k) => v for (k, v) in result_values)
end

"""
Return the final values for the requested expression for each time step for a problem.

Refer to [`read_realized_variable`](@ref) for help and examples.
"""
function read_realized_expression(
    res::SimulationProblemResults,
    expression::AbstractString;
    kwargs...,
)
    return first(
        values(
            read_realized_expressions(
                res,
                [_deserialize_key(ExpressionKey, res, expression)];
                kwargs...,
            ),
        ),
    )
end

function read_realized_expression(res::SimulationProblemResults, expression...; kwargs...)
    return first(
        values(read_realized_expressions(res, [ExpressionKey(expression...)]; kwargs...)),
    )
end

"""
Return the optimizer stats for the problem as a DataFrame.

# Accepted keywords

  - `store::SimulationStore`: a store that has been opened for reading
"""
function read_optimizer_stats(res::SimulationProblemResults; store=nothing)
    if store === nothing && res.store !== nothing
        # In this case we have an InMemorySimulationStore.
        store = res.store
    end
    return _read_optimizer_stats(res, store)
end

function _read_optimizer_stats(res::SimulationProblemResults, ::Nothing)
    open_store(
        HdfSimulationStore,
        joinpath(get_execution_path(res), "data_store"),
        "r",
    ) do store
        _read_optimizer_stats(res, store)
    end
end

"""
Save the realized results to CSV files for all variables, paramaters, duals, auxiliary variables,
expressions, and optimizer statistics.

# Arguments

  - `res::Union{ProblemResults, SimulationProblmeResults`: Results
  - `save_path::AbstractString` : path to save results (defaults to simulation path)
"""
function export_realized_results(res::SimulationProblemResults)
    save_path = mkpath(joinpath(res.results_output_folder, "export"))
    return export_realized_results(res, save_path)
end

function export_realized_results(
    res::Union{ProblemResults, SimulationProblemResults},
    save_path::AbstractString,
)
    if !isdir(save_path)
        throw(IS.ConflictingInputsError("Specified path is not valid."))
    end
    write_data(read_variables_with_keys(res, list_variable_keys(res)), save_path)
    !isempty(list_dual_keys(res)) &&
        write_data(read_duals_with_keys(res, list_dual_keys(res)), save_path; name="dual")
    !isempty(list_parameter_keys(res)) && write_data(
        read_parameters_with_keys(res, list_parameter_keys(res)),
        save_path;
        name="parameter",
    )
    !isempty(list_aux_variable_keys(res)) && write_data(
        read_aux_variables_with_keys(res, list_aux_variable_keys(res)),
        save_path;
        name="aux_variable",
    )
    !isempty(list_expression_keys(res)) && write_data(
        read_expressions_with_keys(res, list_expression_keys(res)),
        save_path;
        name="expression",
    )
    export_optimizer_stats(res, save_path)
    files = readdir(save_path)
    compute_file_hash(save_path, files)
    @info("Files written to $save_path folder.")
    return save_path
end

"""
Save the optimizer statistics to CSV or JSON

# Arguments

  - `res::Union{ProblemResults, SimulationProblmeResults`: Results
  - `directory::AbstractString` : target directory
  - `format = "CSV"` : can be "csv" or "json
"""
function export_optimizer_stats(
    res::Union{ProblemResults, SimulationProblemResults},
    directory::AbstractString;
    format="csv",
)
    data = read_optimizer_stats(res)
    isnothing(data) && return
    if uppercase(format) == "CSV"
        CSV.write(joinpath(directory, "optimizer_stats.csv"), data)
    elseif uppercase(format) == "JSON"
        JSON.write(joinpath(directory, "optimizer_stats.json"), JSON.json(to_dict(data)))
    else
        throw(error("writing optimizer stats only supports csv or json formats"))
    end
end


text:
#=
rect_sino.jl
2019-07-12, Helena H
2019-07-13, Jeff Fessler, refactor to use sg.grid
=#

export rect_sino

using Plots: plot
#using MIRT: sino_geom, MIRT_sino_geom, downsample2

"""
`sino = rect_sino(sg, rects ; oversample=1, xscale=1, yscale=1)`

Create sinogram projection of one or more rectangles.
Works for any sinogram geometry.

in
- `sg::MIRT_sino_geom,`		sinogram geometry object from `sino_geom()`
- `rects::Matrix`			`[ne 6]` rectangle parameters
							`[centx centy widthx widthy angle_degrees value]`

options
- `oversample::Int`	oversampling factor for emulating "strips"
						default 1: just 1 ray per detector element
- `xscale::Int`		use -1 to flip in x (not recommended); default 1
- `yscale::Int`		use -1 to flip in y (not recommended); default 1

out
- `sino`		`[nb na]` sinogram

To get the sample locations, use `(pos,ang) = sg.grid`
"""
function rect_sino(sg::MIRT_sino_geom, rects::AbstractMatrix{<:Real} ;
		oversample::Int = 1, kwargs...)

	sg = sg.over(oversample)
	(rg, ϕg) = sg.grid
	sino = rect_sino(rg, ϕg, rects; kwargs...)
	if oversample > 1
		sino = downsample2(sino, [oversample, 1])
	end
	return sino
end

"""
`sino = rect_sino(rg::AbstractArray{<:Real}, ϕg::AbstractArray{<:Real},
		rects ; oversample=1, xscale=1, yscale=1)`

Create sinogram projection of one or more rectangles,
for arbitrary radial/angular sampling grid locations `(rg, ϕg)`

in
- `rg::AbstractArray{<:Real}`	radial sampling locations
- `ϕg::AbstractArray{<:Real}`	angular sampling locations (radians)
- `rects::Matrix`			`[n 6]` rectangle parameters
							`[centx centy widthx widthy angle_degrees value]`

options
- `xscale::Int`		use -1 to flip in x (not recommended); default 1
- `yscale::Int`		use -1 to flip in y (not recommended); default 1

out
- `sino::AbstractArray{Float32}` same size as `rg` and `ϕg`
"""
function rect_sino(rg::AbstractArray{<:Real}, ϕg::AbstractArray{<:Real},
		rects::AbstractMatrix{<:Real} ;
		xscale::Int = 1,
		yscale::Int = 1,
	)

	size(rects,2) != 6 && throw("6 parameters per rectangle")
	size(rg) != size(ϕg) && throw("rg and ϕg size mismatch")

	sino = zeros(Float32, size(rg))

	cangs = cos.(ϕg)
	sangs = sin.(ϕg)

	#loop over rects
	#ticker reset
	ne = size(rects, 1)
	for ie in 1:ne
		#ticker(mfilename, ie, ne)
		rect = rects[ie, :]

		cx = rect[1] * xscale
		cy = rect[2] * yscale
		wx = rect[3]
		wy = rect[4]
		(wx <= 0) || (wy <= 0) && throw("need positive rectangle sizes")
		eang = deg2rad(rect[5])
		val = rect[6]

		if yscale == -1
			eang = -eang
		end
		if xscale == -1
			eang = pi - eang
		end

		cos_rot = cangs .* cos(eang) + sangs .* sin(eang)
		sin_rot = sangs * cos(eang) - cangs * sin(eang)
		rp = sqrt.((wx * cos_rot).^2 + (wy * sin_rot).^2) # projected radius

		sp = cx .* cangs + cy .* sangs # radial shift
		dis = (rg - sp) ./ rp # scaled distance from center

		# projection angle after affine scaling and rotate
		abs_cos_ang_pi = wx * abs.(cos_rot) ./ rp
		abs_sin_ang_pi = wy * abs.(sin_rot) ./ rp

		# break points of the trapezoid
		len = 1 ./ max.(abs_cos_ang_pi, abs_sin_ang_pi)
		dmax = (abs_cos_ang_pi + abs_sin_ang_pi) / 2
		dbreak = abs.(abs_cos_ang_pi - abs_sin_ang_pi) / 2
		dmax_break = dmax - dbreak
		scale = val * wx * wy ./ rp .* len #lmax

		sino += scale .* trapezoid.(dis, -dmax, -dbreak, dbreak, dmax)
	end

	return sino
end

"""
`trapezoid(t::Real, t1, t2, t3, t4)`
"""
function trapezoid(t::Real, t1::Real, t2::Real, t3::Real, t4::Real)
	if t1 < t < t2
		return (t - t1)/(t2 - t1)
	elseif t2 <= t <= t3
		return 1
	elseif t3 < t < t4
		return (t4 - t)/(t4 - t3)
	else
		return 0
	end
end

"""
`rect_sino()`

show doc strings
"""
function rect_sino()
	@doc rect_sino
end

"""
`rect_sino_test()`

internal test routine: standard sampling
"""
function rect_sino_test()
	rect = [60 10.5 38 27 0 1]
	sg = sino_geom(:ge1, down=8)
	rect_sino(sg, rect; xscale=-1, yscale=-1) # test scale
	t = LinRange(-2, 5, 101)
	f = trapezoid.(t, 1, 2, 3, 4.)
#	plot(t,f)
	true
end

"""
`rect_sino_show()`
show examples
"""
function rect_sino_show( ;
		down::Int = 4,
		rect::Matrix = [70 20 70 30 0 1],
		orbit::Real = 360,
		na::Int = 400,
		oversample::Int = 2,
	)

#=
	ig = image_geom(nx=512, ny=504, fov=500)
	ig = ig.down(down)
	xtrue = rect_im(ig, rect; oversample=4)
=#

	arg = (na=na, down=down, orbit=orbit, offset=0.25)
	geoms = (
		sino_geom(:par, nb = 444, d = 1 ; arg...),
		sino_geom(:fan, nb = 888, d = 1 ; arg..., dsd = 949, dod = 408),
		sino_geom(:fan, nb = 888, d = 1 ; arg...,
			dsd = 949, dod = 408, dfs = Inf), # flat fan
		sino_geom(:moj, nb = 666, d = 1 ; arg...),
	)

	ngeom = length(geoms)
	pl = Array{Plot}(undef, ngeom)

	for ii=1:ngeom
		sg = geoms[ii]
		sino = rect_sino(sg, rect; oversample=oversample)
		dfs = sg.how === :fan ? " dfs=$(sg.dfs)" : ""
		pl[ii] = jim(sg.r, sg.ad, sino, title="$(sg.how)$dfs")
	end
	plot(pl...)
end

"""
`rect_sino(:test)`
self test
"""
function rect_sino(test::Symbol)
	test === :show && return rect_sino_show()
	test != :test && throw(ArgumentError("test $test"))
	rect_sino() # doc
	rect_sino_test()
	rect_sino(:show)
	true
end


text:
Base.@kwdef mutable struct Config
    meta :: String = ""
    auto_sys_last_traceback :: Bool = true
    auto_fix_qt_plugin_path :: Bool = true
    auto_ipython_display :: Bool = true
end

const CONFIG = Config()

# TODO: load_config(), save_config()


text:
# quick script i use during development …
using StaticArrays
using RadialBasisFunctionModels
import Zygote as Zyg

f = x -> [ 1 + x[1]; sum(x.^2) ]
X = [ -3 .+ 6 * rand(2) for i = 1 : 5]
Y = f.(X)

#%%
# Initialize the `RadialFunction` to use for the RBF model:
φ = Gaussian()

# Construct an interpolating model with linear polynomial tail:
rbf = RBFModel( X, Y, φ, 1; centers = C, interpolation_indices = I )

#=
using StaticArrays
ξ = @SVector(rand(2))
jac( rbf, ξ )
auto_jac( rbf, ξ)
=#
#%%

l = function( x )
    global X,Y,φ
    rbf_m = RBFModel( X, Y, φ, 1; centers = C, interpolation_indices = I )
    rbf_m(x, 1)
end

Zyg.gradient( l, rand(2))


text:
# Comparison

function (==)(a::MultiValue,b::MultiValue)
  a.array == b.array
end

function (≈)(a::MultiValue,b::MultiValue)
  a.array ≈ b.array
end

function (≈)(
  a::AbstractArray{<:MultiValue}, b::AbstractArray{<:MultiValue})
  if size(a) != size(b); return false; end
  for (ai,bi) in zip(a,b)
    if !(ai≈bi); return false; end
  end
  true
end

# Addition / subtraction

for op in (:+,:-)
  @eval begin

    function ($op)(a::MultiValue{S}) where S
      r = $op(a.array)
      MultiValue(r)
    end

    function ($op)(a::MultiValue{S},b::MultiValue{S}) where S
      r = $op(a.array, b.array)
      MultiValue(r)
    end

  end
end

# Matrix Division

function (\)(a::TensorValue, b::MultiValue)
  r = a.array \ b.array
  MultiValue(r)
end

# Operations with other numbers

for op in (:+,:-,:*)
  @eval begin
    ($op)(a::MultiValue,b::Number) = MultiValue($op(a.array,b))
    ($op)(a::Number,b::MultiValue) = MultiValue($op(a,b.array))
  end
end

(/)(a::MultiValue,b::Number) = MultiValue(a.array/b)

# Dot product (simple contraction)

(*)(a::VectorValue{D}, b::VectorValue{D}) where D = inner(a,b)

function (*)(a::MultiValue,b::MultiValue)
  r = a.array * b.array
  MultiValue(r)
end

@inline dot(u::VectorValue,v::VectorValue) = inner(u,v)

# Inner product (full contraction)

inner(a::Real,b::Real) = a*b

"""
"""
@generated function inner(a::MultiValue{S,T,N,L}, b::MultiValue{S,W,N,L}) where {S,T,N,L,W}
  str = join([" a.array.data[$i]*b.array.data[$i] +" for i in 1:L ])
  Meta.parse(str[1:(end-1)])
end

# Reductions

for op in (:sum,:maximum,:minimum)
  @eval begin
    $op(a::MultiValue) = $op(a.array)
  end
end

# Outer product (aka dyadic product)

"""
"""
outer(a::Real,b::Real) = a*b

outer(a::MultiValue,b::Real) = a*b

outer(a::Real,b::MultiValue) = a*b

@generated function outer(a::VectorValue{D},b::VectorValue{Z}) where {D,Z}
  str = join(["a.array[$i]*b.array[$j], " for j in 1:Z for i in 1:D])
  Meta.parse("MultiValue(SMatrix{$D,$Z}($str))")
end

@generated function outer(a::VectorValue{D},b::MultiValue{Tuple{A,B}}) where {D,A,B}
  str = join(["a.array[$i]*b.array[$j,$k], "  for k in 1:B for j in 1:A for i in 1:D])
  Meta.parse("MultiValue(SArray{Tuple{$D,$A,$B}}($str))")
end

# Linear Algebra

det(a::TensorValue) = det(a.array)

inv(a::TensorValue) = MultiValue(inv(a.array))

# Measure

"""
"""
meas(a::VectorValue) = sqrt(inner(a,a))

meas(a::TensorValue) = abs(det(a))

@inline norm(u::VectorValue) = sqrt(inner(u,u))

@inline norm(u::VectorValue{0,T}) where T = sqrt(zero(T))

# conj

conj(a::MultiValue) = MultiValue(conj(a.array))

# Trace

@generated function tr(v::TensorValue{D}) where D
  str = join([" v.array.data[$i+$((i-1)*D)] +" for i in 1:D ])
  Meta.parse(str[1:(end-1)])
end

@generated function tr(v::MultiValue{Tuple{A,A,B}}) where {A,B}
  str = ""
  for k in 1:B
    for i in 1:A
      if i !=1
        str *= " + "
      end
      str *= " v.array[$i,$i,$k]"
    end
    str *= ", "
  end
  Meta.parse("VectorValue($str)")
end

# Adjoint and transpose

function adjoint(v::TensorValue)
  t = adjoint(v.array)
  TensorValue(t)
end

function transpose(v::TensorValue)
  t = transpose(v.array)
  TensorValue(t)
end

# Symmetric part

"""
"""
@generated function symmetic_part(v::TensorValue{D}) where D
  str = "("
  for j in 1:D
    for i in 1:D
      str *= "0.5*v.array.data[$i+$((j-1)*D)] + 0.5*v.array.data[$j+$((i-1)*D)], "
    end
  end
  str *= ")"
  Meta.parse("TensorValue($str)")
end


text:
# ---
# title: 1465. Maximum Area of a Piece of Cake After Horizontal and Vertical Cuts
# id: problem1465
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/maximum-area-of-a-piece-of-cake-after-horizontal-and-vertical-cuts/description/>
# hidden: true
# ---
#
# Given a rectangular cake with height `h` and width `w`, and two arrays of
# integers `horizontalCuts` and `verticalCuts` where `horizontalCuts[i]` is the
# distance from the top of the rectangular cake to the `ith` horizontal cut and
# similarly, `verticalCuts[j]` is the distance from the left of the rectangular
# cake to the `jth` vertical cut.
#
# _Return the maximum area of a piece of cake after you cut at each horizontal
# and vertical position provided in the arrays`horizontalCuts` and
# `verticalCuts`. _Since the answer can be a huge number, return this modulo
# 10^9 + 7.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/05/14/leetcode_max_area_2.png)
#
#
#
#     Input: h = 5, w = 4, horizontalCuts = [1,2,4], verticalCuts = [1,3]
#     Output: 4
#     Explanation: The figure above represents the given rectangular cake. Red lines are the horizontal and vertical cuts. After you cut the cake, the green piece of cake has the maximum area.
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/05/14/leetcode_max_area_3.png)**
#
#
#
#     Input: h = 5, w = 4, horizontalCuts = [3,1], verticalCuts = [1]
#     Output: 6
#     Explanation: The figure above represents the given rectangular cake. Red lines are the horizontal and vertical cuts. After you cut the cake, the green and yellow pieces of cake have the maximum area.
#
#
# **Example 3:**
#
#
#
#     Input: h = 5, w = 4, horizontalCuts = [3], verticalCuts = [3]
#     Output: 9
#
#
#
#
# **Constraints:**
#
#   * `2 <= h, w <= 10^9`
#   * `1 <= horizontalCuts.length < min(h, 10^5)`
#   * `1 <= verticalCuts.length < min(w, 10^5)`
#   * `1 <= horizontalCuts[i] < h`
#   * `1 <= verticalCuts[i] < w`
#   * It is guaranteed that all elements in `horizontalCuts` are distinct.
#   * It is guaranteed that all elements in `verticalCuts` are distinct.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# Group calls

export groupStart, groupEnd

groupStart() = ncclGroupStart()
groupEnd()   = ncclGroupEnd()


text:
using ApproximateDPs
const ADP = ApproximateDPs
using FlightSims
const FS = FlightSims
using Transducers
using Plots
using Random, LinearAlgebra, ComponentArrays
using DynamicPolynomials, UnPack
using DifferentialEquations, DataFrames
using NumericalIntegration: integrate

function initialise()
    # setting
    A, B = [-10 1; -0.002 -2], [0; 2]
    Q, R = I, I
    env = LinearSystem(A, B)
    __t = 0.0
    irl = CTLinearValueIterationIRL(Q, R)
    env, irl
end

function train!(env, irl; Δt=0.01, tf=3.0, w_tol=1e-3)
    @unpack A, B = env
    args_linearsystem = (A, B)
    linearsystem, integ = FSimZoo.LinearSystem_SingleIntegrator(args_linearsystem)  # integrated system with scalar integrator ∫r
    x0 = State(linearsystem, integ)([0.4, 4.0])
    irl.V̂.param = zeros(size(irl.V̂.param))  # zero initialisation
    û = ADP.ApproximateOptimalInput(irl, B)
    _û = (X, p, t) -> û(X.x, p, t)  # for integrated system

    cb_train = ADP.update_params_callback(irl, w_tol)
    cb = CallbackSet(cb_train)
    running_cost = ADP.RunningCost(irl)
    prob, df = sim(
                   x0,
                   apply_inputs(Dynamics!(linearsystem, integ, running_cost); u=_û);
                   tf=tf,
                   callback=cb,
                   savestep=Δt
                  )
    ts = df.time
    xs = df.sol |> Map(datum -> datum.linearsystem.state) |> collect
    plot(ts, hcat(xs...)')
    # ∫rs = df.sol |> Map(datum -> datum.integ.integral) |> collect
    # plot(hcat(∫rs...)')
end

function main(; seed=1)
    Random.seed!(seed)
    env, irl = initialise()
    train!(env, irl; w_tol = 1e-3)
end


text:
#Find Volumes of Various Shapes.
#Wikipedia reference: https://en.wikipedia.org/wiki/Volume
# contributed By:- [Qingyu Qu](https://github.com/ErikQQY)

"""
    vol_cube()

Compute the volume of a cube.

# Examples

```julia-repl
julia> vol_cube(1)
1
julia> vol_cube(3)
27
julia> vol_cube(-1)
DomainError
```
"""
function vol_cube(side)
    if side<0
        throw(DomainError("vol_cube() only works for non negative values."))
    end
    return side^3
end

"""
    vol_cuboid(width, height, length)

Compute the volume of a vol_cuboid

# Examples

```julia-repl
julia> vol_cuboid(1, 1, 1)
1
julia> vol_cuboid(1, 2, 3)
6
```
"""
function vol_cuboid(width, height, length)
    if width < 0 || height <0 || length < 0
        throw(DomainError("vol_cuboid() only works for non negative values."))
    end
    return width*height*length
end

"""
    vol_cone(area_of_base, height)

Compute the Volume of a Cone

# Examples

```julia-repl
julia> vol_cone(10, 3)
10.0
julia> vol_cone(1, 1)
0.3333333333333333
```
"""
function vol_cone(area_of_base, height)
    if area_of_base < 0 || height <0
        throw(DomainError("vol_cone() only works for non negative values."))
    end
    return area_of_base*height/3.0
end

"""
    vol_right_circ_cone(radius, height)

Compute the Volume of a Right Circular Cone.

# Examples

```julia-repl
julia> vol_right_circ_cone(2, 3)
12.566370614359172
```
"""
function vol_right_circ_cone(radius, height)
    if radius < 0 || height < 0
        throw(DomainError("vol_right_circ_cone() only works for non negative values."))
    end
    return π*radius^2*height/3.0
end

"""
    vol_prism(area_of_base, height)

Compute the Volume of a Prism.

# Examples

```julia-repl
julia> vol_prism(10, 2)
20.0
julia> vol_prism(11, 1)
11.0
```
"""
function vol_prism(area_of_base, height)
    if area_of_base < 0 || height < 0
        throw(DomainError("vol_prism() only works for non negative values."))
    end
    return area_of_base*height
end

"""
    vol_pyramid(area_of_base, height)

Compute the volume of a Pyramid.

# Examples

```julia-repl
julia> vol_pyramid(10, 3)
10.0
julia> vol_pyramid(1.5, 3)
1.5
```
"""
function vol_pyramid(area_of_base, height)
    if area_of_base < 0 || height < 0
        throw(DomainError("vol_pyramid() only works for non negative values."))
    end
    return area_of_base*height/3.0
end

"""
    vol_sphere(radius)

Compute the volume of a sphere.

# Examples

```julia
vol_sphere(5)  # returns 523.5987755982989
vol_sphere(1)  # returns 4.1887902047863905
vol_sphere(-1) # returns DomainError
```
"""
function vol_sphere(radius)
    if radius < 0
        throw(DomainError("vol_sphere() only works for non negative values."))
    end
    return 4/3*pi*(radius)^3
end

"""
    vol_circular_cylinder(area_of_, height)

Compute the Volume of a Circular Cylinder.

# Examples

```julia-repl
julia> vol_circular_cylinder(1, 1)
3.141592653589793
julia> vol_circular_cylinder(4, 3)
150.79644737231007
"""
function vol_circular_cylinder(radius, height)
    if radius < 0 || height < 0
        throw(DomainError("vol_circular_cylinder() only works for non negative values."))
    end
    return pi*radius^2*height
end

max_stars_repo_name: KohRongSoon/Julia

max_issues_repo_name: KohRongSoon/Julia

max_forks_repo_name: KohRongSoon/Julia


text:
#---------------------------------------------------------------------#
#This function performs the time-integration
#Written by F.X. Giraldo on April 19, 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
include("ti_SSP.jl")
include("ti_LSRK.jl")

function time_integration!(q0,u,coord,Dhat,Fhat,intma,periodicity,time,ntime,dt,stages,DFloat,ti_type,space_method,plot_movie)

    if (ti_type == "SSP")
        q0 = ti_SSP!(q0,u,coord,Dhat,Fhat,intma,periodicity,time,ntime,dt,stages,space_method,plot_movie,DFloat)
    elseif (ti_type == "LSRK")
        q0 = ti_LSRK!(q0,u,coord,Dhat,Fhat,intma,periodicity,time,ntime,dt,space_method,plot_movie,DFloat)
    end
    return (q0)
end


text:
# Sundials.jl algorithms

# Abstract Types
abstract type SundialsODEAlgorithm{Method, LinearSolver} <: DiffEqBase.AbstractODEAlgorithm end
abstract type SundialsDAEAlgorithm{LinearSolver} <: DiffEqBase.AbstractDAEAlgorithm end
abstract type SundialsNonlinearSolveAlgorithm{LinearSolver} end

# ODE Algorithms
"""
CVODE_BDF: CVode Backward Differentiation Formula (BDF) solver.

method - This is the method for solving the implicit equation. For BDF this defaults to :Newton while for Adams this defaults to :Functional. These choices match the recommended pairing in the Sundials.jl manual. However, note that using the :Newton method may take less iterations but requires more memory than the :Function iteration approach.
linear_solver - This is the linear solver which is used in the :Newton method.
The choices for the linear solver are:

:Dense - A dense linear solver.
:Band - A solver specialized for banded Jacobians. If used, you must set the position of the upper and lower non-zero diagonals via jac_upper and jac_lower.
:LapackDense - A version of the dense linear solver that uses the Julia-provided OpenBLAS-linked LAPACK for multithreaded operations. This will be faster than :Dense on larger systems but has noticable overhead on smaller (<100 ODE) systems.
:LapackBand - A version of the banded linear solver that uses the Julia-provided OpenBLAS-linked LAPACK for multithreaded operations. This will be faster than :Band on larger systems but has noticable overhead on smaller (<100 ODE) systems.
:Diagonal - This method is specialized for diagonal Jacobians.
:GMRES - A GMRES method. Recommended first choice Krylov method
:BCG - A Biconjugate gradient method.
:PCG - A preconditioned conjugate gradient method. Only for symmetric linear systems.
:TFQMR - A TFQMR method.
:KLU - A sparse factorization method. Requires that the user specifies a Jacobian. The Jacobian must be set as a sparse matrix in the ODEProblem type.

Example:

CVODE_BDF() # BDF method using Newton + Dense solver
CVODE_BDF(method=:Functional) # BDF method using Functional iterations
CVODE_BDF(linear_solver=:Band,jac_upper=3,jac_lower=3) # Banded solver with nonzero diagonals 3 up and 3 down
CVODE_BDF(linear_solver=:BCG) # Biconjugate gradient method

All of the additional options are available. The full constructor is:

CVODE_BDF(;method=:Newton,linear_solver=:Dense,
          jac_upper=0,jac_lower=0,
          stored_upper = jac_upper + jac_lower,
          non_zero=0,krylov_dim=0,
          stability_limit_detect=false,
          max_hnil_warns = 10,
          max_order = 5,
          max_error_test_failures = 7,
          max_nonlinear_iters = 3,
          max_convergence_failures = 10,
          prec = nothing, prec_side = 0)
"""
struct CVODE_BDF{Method, LinearSolver, P, PS} <: SundialsODEAlgorithm{Method, LinearSolver}
    jac_upper::Int
    jac_lower::Int
    krylov_dim::Int
    stability_limit_detect::Bool
    max_hnil_warns::Int
    max_order::Int
    max_error_test_failures::Int
    max_nonlinear_iters::Int
    max_convergence_failures::Int
    prec::P
    psetup::PS
    prec_side::Int
end
Base.@pure function CVODE_BDF(;
    method = :Newton,
    linear_solver = :Dense,
    jac_upper = 0,
    jac_lower = 0,
    non_zero = 0,
    krylov_dim = 0,
    stability_limit_detect = false,
    max_hnil_warns = 10,
    max_order = 5,
    max_error_test_failures = 7,
    max_nonlinear_iters = 3,
    max_convergence_failures = 10,
    prec = nothing,
    psetup = nothing,
    prec_side = 0,
)
    if linear_solver == :Band && (jac_upper == 0 || jac_lower == 0)
        error("Banded solver must set the jac_upper and jac_lower")
    end
    if !(
        linear_solver in (
            :None,
            :Diagonal,
            :Dense,
            :LapackDense,
            :Band,
            :LapackBand,
            :BCG,
            :GMRES,
            :FGMRES,
            :PCG,
            :TFQMR,
            :KLU,
        )
    )
        error("Linear solver choice not accepted.")
    end
    CVODE_BDF{method, linear_solver, typeof(prec), typeof(psetup)}(
        jac_upper,
        jac_lower,
        krylov_dim,
        stability_limit_detect,
        max_hnil_warns,
        max_order,
        max_error_test_failures,
        max_nonlinear_iters,
        max_convergence_failures,
        prec,
        psetup,
        prec_side,
    )
end
"""
CVODE_Adams: CVode Adams-Moulton solver.

method - This is the method for solving the implicit equation. For BDF this defaults to :Newton while for Adams this defaults to :Functional. These choices match the recommended pairing in the Sundials.jl manual. However, note that using the :Newton method may take less iterations but requires more memory than the :Function iteration approach.
linear_solver - This is the linear solver which is used in the :Newton method.
The choices for the linear solver are:

:Dense - A dense linear solver.
:Band - A solver specialized for banded Jacobians. If used, you must set the position of the upper and lower non-zero diagonals via jac_upper and jac_lower.
:LapackDense - A version of the dense linear solver that uses the Julia-provided OpenBLAS-linked LAPACK for multithreaded operations. This will be faster than :Dense on larger systems but has noticable overhead on smaller (<100 ODE) systems.
:LapackBand - A version of the banded linear solver that uses the Julia-provided OpenBLAS-linked LAPACK for multithreaded operations. This will be faster than :Band on larger systems but has noticable overhead on smaller (<100 ODE) systems.
:Diagonal - This method is specialized for diagonal Jacobians.
:GMRES - A GMRES method. Recommended first choice Krylov method
:BCG - A Biconjugate gradient method.
:PCG - A preconditioned conjugate gradient method. Only for symmetric linear systems.
:TFQMR - A TFQMR method.
:KLU - A sparse factorization method. Requires that the user specifies a Jacobian. The Jacobian must be set as a sparse matrix in the ODEProblem type.

All of the additional options are available. The full constructor is:

CVODE_Adams(;method=:Functional,linear_solver=:None,
            jac_upper=0,jac_lower=0,
            stored_upper = jac_upper + jac_lower,
            krylov_dim=0,
            stability_limit_detect=false,
            max_hnil_warns = 10,
            max_order = 12,
            max_error_test_failures = 7,
            max_nonlinear_iters = 3,
            max_convergence_failures = 10,
            prec = nothing, psetup = nothing, prec_side = 0)
"""
struct CVODE_Adams{Method, LinearSolver, P, PS} <:
       SundialsODEAlgorithm{Method, LinearSolver}
    jac_upper::Int
    jac_lower::Int
    krylov_dim::Int
    stability_limit_detect::Bool
    max_hnil_warns::Int
    max_order::Int
    max_error_test_failures::Int
    max_nonlinear_iters::Int
    max_convergence_failures::Int
    prec::P
    psetup::PS
    prec_side::Int
end
Base.@pure function CVODE_Adams(;
    method = :Functional,
    linear_solver = :None,
    jac_upper = 0,
    jac_lower = 0,
    krylov_dim = 0,
    stability_limit_detect = false,
    max_hnil_warns = 10,
    max_order = 12,
    max_error_test_failures = 7,
    max_nonlinear_iters = 3,
    max_convergence_failures = 10,
    prec = nothing,
    psetup = nothing,
    prec_side = 0,
)
    if linear_solver == :Band && (jac_upper == 0 || jac_lower == 0)
        error("Banded solver must set the jac_upper and jac_lower")
    end
    if !(
        linear_solver in (
            :None,
            :Diagonal,
            :Dense,
            :LapackDense,
            :Band,
            :LapackBand,
            :BCG,
            :GMRES,
            :FGMRES,
            :PCG,
            :TFQMR,
            :KLU,
        )
    )
        error("Linear solver choice not accepted.")
    end
    CVODE_Adams{method, linear_solver, typeof(prec), typeof(psetup)}(
        jac_upper,
        jac_lower,
        krylov_dim,
        stability_limit_detect,
        max_hnil_warns,
        max_order,
        max_error_test_failures,
        max_nonlinear_iters,
        max_convergence_failures,
        prec,
        psetup,
        prec_side,
    )
end
"""
ARKODE: Explicit and ESDIRK Runge-Kutta methods of orders 2-8 depending on choice of options.

The main options for ARKODE are the choice between explicit and implicit and the method order, given via:

ARKODE(Sundials.Explicit()) # Solve with explicit tableau of default order 4
ARKODE(Sundials.Implicit(),order = 3) # Solve with explicit tableau of order 3
The order choices for explicit are 2 through 8 and for implicit 3 through 5. Specific methods can also be set through the etable and itable options for explicit and implicit tableaus respectively. The available tableaus are:

etable:
    HEUN_EULER_2_1_2: 2nd order Heun's method
    BOGACKI_SHAMPINE_4_2_3:
    ARK324L2SA_ERK_4_2_3: explicit portion of Kennedy and Carpenter's 3rd order method
    ZONNEVELD_5_3_4: 4th order explicit method
    ARK436L2SA_ERK_6_3_4: explicit portion of Kennedy and Carpenter's 4th order method
    SAYFY_ABURUB_6_3_4: 4th order explicit method
    CASH_KARP_6_4_5: 5th order explicit method
    FEHLBERG_6_4_5: Fehlberg's classic 5th order method
    DORMAND_PRINCE_7_4_5: the classic 5th order Dormand-Prince method
    ARK548L2SA_ERK_8_4_5: explicit portion of Kennedy and Carpenter's 5th order method
    VERNER_8_5_6: Verner's classic 5th order method
    FEHLBERG_13_7_8: Fehlberg's 8th order method

itable:
    SDIRK_2_1_2: An A-B-stable 2nd order SDIRK method
    BILLINGTON_3_3_2: A second order method with a 3rd order error predictor of less stability
    TRBDF2_3_3_2: The classic TR-BDF2 method
    KVAERNO_4_2_3: an L-stable 3rd order ESDIRK method
    ARK324L2SA_DIRK_4_2_3: implicit portion of Kennedy and Carpenter's 3th order method
    CASH_5_2_4: Cash's 4th order L-stable SDIRK method
    CASH_5_3_4: Cash's 2nd 4th order L-stable SDIRK method
    SDIRK_5_3_4: Hairer's 4th order SDIRK method
    KVAERNO_5_3_4: Kvaerno's 4th order ESDIRK method
    ARK436L2SA_DIRK_6_3_4: implicit portion of Kennedy and Carpenter's 4th order method
    KVAERNO_7_4_5: Kvaerno's 5th order ESDIRK method
    ARK548L2SA_DIRK_8_4_5: implicit portion of Kennedy and Carpenter's 5th order method

These can be set for example via:

ARKODE(Sundials.Explicit(),etable = Sundials.DORMAND_PRINCE_7_4_5)
ARKODE(Sundials.Implicit(),itable = Sundials.KVAERNO_4_2_3)

All of the additional options are available. The full constructor is:

ARKODE(stiffness=Sundials.Implicit();
      method=:Newton,linear_solver=:Dense,
      jac_upper=0,jac_lower=0,stored_upper = jac_upper+jac_lower,
      non_zero=0,krylov_dim=0,
      max_hnil_warns = 10,
      max_error_test_failures = 7,
      max_nonlinear_iters = 3,
      max_convergence_failures = 10,
      predictor_method = 0,
      nonlinear_convergence_coefficient = 0.1,
      dense_order = 3,
      order = 4,
      set_optimal_params = false,
      crdown = 0.3,
      dgmax = 0.2,
      rdiv = 2.3,
      msbp = 20,
      adaptivity_method = 0,
      prec = nothing, psetup = nothing, prec_side = 0
      )
"""
struct ARKODE{Method, LinearSolver, MassLinearSolver, T, T1, T2, P, PS} <:
       SundialsODEAlgorithm{Method, LinearSolver}
    stiffness::T
    jac_upper::Int
    jac_lower::Int
    mass_upper::Int
    mass_lower::Int
    krylov_dim::Int
    mass_krylov_dim::Int
    max_hnil_warns::Int
    max_error_test_failures::Int
    max_nonlinear_iters::Int
    max_convergence_failures::Int
    predictor_method::Int
    nonlinear_convergence_coefficient::Float64
    dense_order::Int
    order::Int
    set_optimal_params::Bool
    crdown::Float64
    dgmax::Float64
    rdiv::Float64
    msbp::Int
    itable::T1
    etable::T2
    prec::P
    psetup::PS
    prec_side::Int
end

Base.@pure function ARKODE(
    stiffness = Implicit();
    method = :Newton,
    linear_solver = :Dense,
    mass_linear_solver = :Dense,
    jac_upper = 0,
    jac_lower = 0,
    mass_upper = 0,
    mass_lower = 0,
    non_zero = 0,
    krylov_dim = 0,
    mass_krylov_dim = 0,
    max_hnil_warns = 10,
    max_error_test_failures = 7,
    max_nonlinear_iters = 3,
    max_convergence_failures = 10,
    predictor_method = 0,
    nonlinear_convergence_coefficient = 0.1,
    dense_order = 3,
    order = 4,
    set_optimal_params = false,
    crdown = 0.3,
    dgmax = 0.2,
    rdiv = 2.3,
    msbp = 20,
    adaptivity_method = 0,
    itable = nothing,
    etable = nothing,
    prec = nothing,
    psetup = nothing,
    prec_side = 0,
)
    if linear_solver == :Band && (jac_upper == 0 || jac_lower == 0)
        error("Banded solver must set the jac_upper and jac_lower")
    end
    if !(
        linear_solver in (
            :None,
            :Diagonal,
            :Dense,
            :LapackDense,
            :Band,
            :LapackBand,
            :BCG,
            :GMRES,
            :FGMRES,
            :PCG,
            :TFQMR,
            :KLU,
        )
    )
        error("Linear solver choice not accepted.")
    end
    if !(
        mass_linear_solver in (
            :None,
            :Diagonal,
            :Dense,
            :LapackDense,
            :Band,
            :LapackBand,
            :BCG,
            :GMRES,
            :FGMRES,
            :PCG,
            :TFQMR,
            :KLU,
        )
    )
        error("Mass Matrix Linear solver choice not accepted.")
    end
    ARKODE{
        method,
        linear_solver,
        mass_linear_solver,
        typeof(stiffness),
        typeof(itable),
        typeof(etable),
        typeof(prec),
        typeof(psetup),
    }(
        stiffness,
        jac_upper,
        jac_lower,
        mass_upper,
        mass_lower,
        krylov_dim,
        mass_krylov_dim,
        max_hnil_warns,
        max_error_test_failures,
        max_nonlinear_iters,
        max_convergence_failures,
        predictor_method,
        nonlinear_convergence_coefficient,
        dense_order,
        order,
        set_optimal_params,
        crdown,
        dgmax,
        rdiv,
        msbp,
        itable,
        etable,
        prec,
        psetup,
        prec_side,
    )
end

# DAE Algorithms
"""
IDA: This is the IDA method from the Sundials.jl package.

Note that the constructors for the Sundials algorithms take a main argument:
linearsolver - This is the linear solver which is used in the Newton iterations. The choices are:

:Dense - A dense linear solver.
:Band - A solver specialized for banded Jacobians. If used, you must set the position of the upper and lower non-zero diagonals via jac_upper and jac_lower.
:LapackDense - A version of the dense linear solver that uses the Julia-provided OpenBLAS-linked LAPACK for multithreaded operations. This will be faster than :Dense on larger systems but has noticable overhead on smaller (<100 ODE) systems.
:LapackBand - A version of the banded linear solver that uses the Julia-provided OpenBLAS-linked LAPACK for multithreaded operations. This will be faster than :Band on larger systems but has noticable overhead on smaller (<100 ODE) systems.
:GMRES - A GMRES method. Recommended first choice Krylov method
:BCG - A Biconjugate gradient method.
:PCG - A preconditioned conjugate gradient method. Only for symmetric linear systems.
:TFQMR - A TFQMR method.
:KLU - A sparse factorization method. Requires that the user specifies a Jacobian. The Jacobian must be set as a sparse matrix in the ODEProblem type.

Note that the preconditioner for iterative linear solvers (if supplied) should be a left preconditioner.

Example:

IDA() # Newton + Dense solver
IDA(linear_solver=:Band,jac_upper=3,jac_lower=3) # Banded solver with nonzero diagonals 3 up and 3 down
IDA(linear_solver=:BCG) # Biconjugate gradient method

All of the additional options are available. The constructor is:

IDA(;linear_solver=:Dense,jac_upper=0,jac_lower=0,krylov_dim=0,
    max_order = 5,
    max_error_test_failures = 7,
    max_nonlinear_iters = 3,
    nonlinear_convergence_coefficient = 0.33,
    nonlinear_convergence_coefficient_ic = 0.0033,
    max_num_steps_ic = 5,
    max_num_jacs_ic = 4,
    max_num_iters_ic = 10,
    max_num_backs_ic = 100,
    use_linesearch_ic = true,
    max_convergence_failures = 10,
    init_all = false,
    prec = nothing, psetup = nothing)
"""
struct IDA{LinearSolver, P, PS} <: SundialsDAEAlgorithm{LinearSolver}
    jac_upper::Int
    jac_lower::Int
    krylov_dim::Int
    max_order::Int
    max_error_test_failures::Int
    nonlinear_convergence_coefficient::Float64
    max_nonlinear_iters::Int
    max_convergence_failures::Int
    nonlinear_convergence_coefficient_ic::Float64
    max_num_steps_ic::Int
    max_num_jacs_ic::Int
    max_num_iters_ic::Int
    max_num_backs_ic::Int
    use_linesearch_ic::Bool
    init_all::Bool
    prec::P
    psetup::PS
end
Base.@pure function IDA(;
    linear_solver = :Dense,
    jac_upper = 0,
    jac_lower = 0,
    krylov_dim = 0,
    max_order = 5,
    max_error_test_failures = 7,
    max_nonlinear_iters = 3,
    nonlinear_convergence_coefficient = 0.33,
    nonlinear_convergence_coefficient_ic = 0.0033,
    max_num_steps_ic = 5,
    max_num_jacs_ic = 4,
    max_num_iters_ic = 10,
    max_num_backs_ic = 100,
    use_linesearch_ic = true,
    init_all = false,
    max_convergence_failures = 10,
    prec = nothing,
    psetup = nothing,
)
    if linear_solver == :Band && (jac_upper == 0 || jac_lower == 0)
        error("Banded solver must set the jac_upper and jac_lower")
    end
    if !(
        linear_solver in (
            :None,
            :Diagonal,
            :Dense,
            :LapackDense,
            :Band,
            :LapackBand,
            :BCG,
            :GMRES,
            :FGMRES,
            :PCG,
            :TFQMR,
            :KLU,
        )
    )
        error("Linear solver choice not accepted.")
    end
    IDA{linear_solver, typeof(prec), typeof(psetup)}(
        jac_upper,
        jac_lower,
        krylov_dim,
        max_order,
        max_error_test_failures,
        nonlinear_convergence_coefficient,
        max_nonlinear_iters,
        max_convergence_failures,
        nonlinear_convergence_coefficient_ic,
        max_num_steps_ic,
        max_num_jacs_ic,
        max_num_iters_ic,
        max_num_backs_ic,
        use_linesearch_ic,
        init_all,
        prec,
        psetup,
    )
end
"""
KINSOL: Newton-Krylov technique solver
"""
struct KINSOL{LinearSolver} <: SundialsNonlinearSolveAlgorithm{LinearSolver}
    jac_upper::Int
    jac_lower::Int
    userdata::Any
end
Base.@pure function KINSOL(;
    linear_solver = :Dense,
    jac_upper = 0,
    jac_lower = 0,
    userdata = nothing,
)
    if !(
        linear_solver in (
            :None,
            :Diagonal,
            :Dense,
            :LapackDense,
            :Band,
            :LapackBand,
            :BCG,
            :GMRES,
            :FGMRES,
            :PCG,
            :TFQMR,
            :KLU,
        )
    )
        error("Linear solver choice not accepted.")
    end
    KINSOL{linear_solver}(
        jac_upper,
        jac_lower,
        userdata,
    )
end

method_choice(alg::SundialsODEAlgorithm{Method}) where {Method} = Method
method_choice(alg::SundialsDAEAlgorithm) = :Newton
linear_solver(
    alg::SundialsODEAlgorithm{Method, LinearSolver},
) where {Method, LinearSolver} = LinearSolver
linear_solver(alg::SundialsDAEAlgorithm{LinearSolver}) where {LinearSolver} = LinearSolver
linear_solver(alg::SundialsNonlinearSolveAlgorithm{LinearSolver}) where {LinearSolver} = LinearSolver


text:
"""
    gillespie(em::AbstractEpidemicModel, state_0; kwargs)

Compute an exact trajectory of an `AbstractEpidemicModel` with initial state `state_0` using a generic gillespie algorithm.

`em` must implement the `AbstractEpidemicModel` interface.
Keyword arguments:
* `tmax=100.0`: maximum time of the simulation
* `nmax=1000`: maximum number of iterations of the simulation
* `sampling_method=:tree`: Sampling method for choosing what reaction happens at each iteration. The default is `:tree`, which takes ``O(\\log N)`` per iteration, where ``N`` is the length of the reaction rate vector. Set to `:array` for a method with ``O(N)`` sampling time, but faster for small numbers of reactions (≈ under 1000).
"""
function gillespie(em::AbstractEpidemicModel, state_0; tmax=100.0, nmax=1000, sampling_method=:tree, should_stop=default_stop, process_output=true)
    t = 0.0
    n = 1
    state = init_state(em, state_0)
    output = init_output(em, state, nmax)
    ts = Vector{Float64}(undef, nmax)
    ts[1] = t
    if sampling_method == :tree
        a = CategoricalTree(init_rates(em, state))
    else
        a = init_rates(m, state)
    end
    a0 = sum(a)
    while t < tmax && n < nmax && a0 > 0 && !should_stop(em, state)
        τ = log(1/rand())/a0
        t += τ
        n += 1
        k = rand_categorical(a, a0)
        update_state!(state, a, k, em)
        a0 = sum(a)
        update_output!(output, state, n, k, em)
        ts[n] = t
    end
    if process_output
        return ts[1:n], finalize_output(em, output[1:n])
    else
        return ts[1:n], output[1:n]
    end
end


text:
export SymmorphicSymmetryEmbedding
export embed
export get_irrep_components
export iscompatible
export fractional_momentum

"""
    SymmorphicSymmetryEmbedding{S1<:AbstractSymmetry, S2<:AbstractSymmetry}<:AbstractSymmetryEmbedding

"""
struct SymmorphicSymmetryEmbedding{
    S1<:AbstractSymmetry,
    S2<:AbstractSymmetry
}<:AbstractSymmetryEmbedding
    lattice::Lattice
    normal::SymmetryEmbedding{S1}
    rest::SymmetryEmbedding{S2}

    function SymmorphicSymmetryEmbedding(
        lattice::Lattice,
        symmetry::SymmorphicSymmetry{S1, S2, E},
    ) where {S1, S2, E}
        if !iscompatible(lattice, symmetry)
            throw(ArgumentError("lattice and symmetry are not compatible"))
        end
        return new{S1, S2}(
            lattice,
            embed(lattice, symmetry.normal),
            embed(lattice, symmetry.rest),
        )
    end

    function SymmorphicSymmetryEmbedding(normal::SymmetryEmbedding, rest::SymmetryEmbedding)
        lattice = normal.lattice
        if rest.lattice != lattice
            throw(ArgumentError("two symmetry embeddings should have the same lattice"))
        end
        symmetry = SymmorphicSymmetry(normal.symmetry, rest.symmetry)
        return SymmorphicSymmetryEmbedding(lattice, symmetry)
    end
end

"""
    embed(lattice::Lattice, ssym::SymmorphicSymmetry)
"""
function embed(lattice::Lattice, ssym::SymmorphicSymmetry)
    return SymmorphicSymmetryEmbedding(lattice, ssym)
end

# ⋊(normal::SymmetryEmbedding, rest::SymmetryEmbedding) = SymmorphicSymmetryEmbedding(normal, rest)
# ⋊ˢ(normal::SymmetryEmbedding, rest::SymmetryEmbedding) = SymmorphicSymmetryEmbedding(normal, rest)

# ⋉(rest::SymmetryEmbedding, normal::SymmetryEmbedding) = SymmorphicSymmetryEmbedding(normal, rest)
# ⋉ˢ(rest::SymmetryEmbedding, normal::SymmetryEmbedding) = SymmorphicSymmetryEmbedding(normal, rest)

Base.eltype(::SymmorphicSymmetryEmbedding) = SitePermutation
Base.valtype(::SymmorphicSymmetryEmbedding) = SitePermutation
Base.eltype(::Type{<:SymmorphicSymmetryEmbedding}) = SitePermutation
Base.valtype(::Type{<:SymmorphicSymmetryEmbedding}) = SitePermutation

function Base.:(==)(lhs::SSE, rhs::SSE) where {SSE<:SymmorphicSymmetryEmbedding}
    return lhs.lattice == rhs.lattice && lhs.normal == rhs.normal && lhs.rest == rhs.rest
end

GroupTools.group_order(sym::SymmorphicSymmetryEmbedding) = group_order(sym.normal) * group_order(sym.rest)

function fractional_momentum(sym::SymmorphicSymmetryEmbedding{<:FiniteTranslationSymmetry, S2}, args...) where S2
    return fractional_momentum(sym.normal, args...)
end

"""
    iscompatible(lattice::Lattice, ssym::SymmorphicSymmetry)
"""
function iscompatible(lattice::Lattice, ssym::SymmorphicSymmetry)
    return iscompatible(lattice, ssym.normal) && iscompatible(lattice, ssym.rest)
end


text:
function vec_isless(x::AbstractArray, y::AbstractArray)
    for (xi,yi) in zip(x,y)
        if xi >= yi return false end
    end
    return true
end

mutable struct SurvivorApproximator
    sample::MatF64
    α::Float64
end

function (surv::SurvivorApproximator)(x::AbstractVector{T}) where T<:Real
    count=0
    N = size(surv.sample, 2)
    # count = @parallel (+) for i in 1:N
    #     convert(Int, vec_isless(x, view(surv.sample, :, i)))
    # end
    for i in 1:N
        count += convert(Int, vec_isless(x, view(surv.sample, :, i)))
    end
    prob=count/N
    se=sqrt((prob-prob^2)/N)
    err=norminvcdf((1 + surv.α)/2)*se
    return prob,err
end

function below_nonrisk_frontier(x::AbstractVector{T}, frontier::LinkedList{VecF64}) where T <: Real
    for p in frontier
        if vec_isless(x,p) return true end
    end
    return false
end

function above_risk_frontier(x::AbstractVector{T}, frontier::LinkedList{VecF64}) where T <: Real
    for p in frontier
        if vec_isless(p, x) return true end
    end
    return false
end

function decapitate!(l::Cons)
    l.head=l.tail.head
    l.tail=l.tail.tail
    return l
end


text: include("testbase.jl")