text:
mutable struct CPU
    registers::Array{Int64}
end

struct Instruction
    op::String
    a::Int64
    b::Int64
    c::Int64
end

mutable struct Program
    ip::Int64
    lines::Dict{Int64,Instruction}
end

include("util.jl")

function read_program(path::String)::Program
    lines = getlines(path)
    ip_reg = parse(Int64, split(lines[1], " ")[2])
    proglines = Dict()
    linenum = 0
    for line in lines[2:end]
        parts = split(line, " ")
        op = parts[1]
        a = parse(Int64, parts[2])
        b = parse(Int64, parts[3])
        c = parse(Int64, parts[4])
        proglines[linenum] = Instruction(op, a, b, c)
        linenum += 1
    end
    program = Program(ip_reg, proglines)
    return program
end

function run_instruction(before::CPU, i::Instruction)
    if i.op == "addr"
        before.registers[i.c+1] = before.registers[i.a+1] + before.registers[i.b+1]
    elseif i.op == "addi"
        before.registers[i.c+1] = before.registers[i.a+1] + i.b
    elseif i.op == "mulr"
        before.registers[i.c+1] = before.registers[i.a+1] * before.registers[i.b+1]
    elseif i.op == "muli"
        before.registers[i.c+1] = before.registers[i.a+1] * i.b
    elseif i.op == "banr"
        before.registers[i.c+1] = before.registers[i.a+1] & before.registers[i.b+1]
    elseif i.op == "bani"
        before.registers[i.c+1] = before.registers[i.a+1] & i.b
    elseif i.op == "borr"
        before.registers[i.c+1] = before.registers[i.a+1] | before.registers[i.b+1]
    elseif i.op == "bori"
        before.registers[i.c+1] = before.registers[i.a+1] | i.b
    elseif i.op == "setr"
        before.registers[i.c+1] = before.registers[i.a+1]
    elseif i.op == "seti"
        before.registers[i.c+1] = i.a
    elseif i.op == "gtir"
        if i.a > before.registers[i.b+1]
            before.registers[i.c+1] = 1
        else
            before.registers[i.c+1] = 0
        end
    elseif i.op == "gtri"
        if before.registers[i.a+1] > i.b
            before.registers[i.c+1] = 1
        else
            before.registers[i.c+1] = 0
        end
    elseif i.op == "gtrr"
        if before.registers[i.a+1] > before.registers[i.b+1]
            before.registers[i.c+1] = 1
        else
            before.registers[i.c+1] = 0
        end
    elseif i.op == "eqir"
        if i.a == before.registers[i.b+1]
            before.registers[i.c+1] = 1
        else
            before.registers[i.c+1] = 0
        end
    elseif i.op == "eqri"
        if before.registers[i.a+1] == i.b
            before.registers[i.c+1] = 1
        else
            before.registers[i.c+1] = 0
        end
    elseif i.op == "eqrr"
        if before.registers[i.a+1] == before.registers[i.b+1]
            before.registers[i.c+1] = 1
        else
            before.registers[i.c+1] = 0
        end
    end
end


text:
"
Author: Daniel O'Connor
Email: dt.oconnor@outlook.com
Based off the work published in
https://arxiv.org/abs/1401.7087
S. W. Shin, G. Smith, J. A. Smolin, and U. Vazirani,
How Quantum is the D-Wave Machine?, arXiv:1401.7087[quant-ph] (2014),
"

module SVMC
    include("tools.jl")
    using ..SamplingTools: hj_info, default_temp_range, generate_nbr_dict,
                            generate_nbr, rotor_aggregation, add_noise,
                            add_cross_talk, hj_info!, generate_nbr!
    using Random
    using Distributed
    using Profile
    using BenchmarkTools
    using LinearAlgebra
    using SharedArrays

    function anneal!(spin_vector::Array{Float64, 1},
                    nbs_idxs::Array{Int64, 2},
                    nbs_vals::Array{Float64, 2},
                    a_schedule::Array{Float64, 1},
                    # must be same size or will crash from @inbounds
                    b_schedule::Array{Float64, 1},
                    temp::Float64)
        maxnb = size(nbs_idxs, 1)
        nspins = size(spin_vector, 1)
        perm = collect(1:nspins)
        istep = 1
        z_mag = 0.0
        a_field = 0.0
        b_field = 0.0
        sweep = 1
        ispin = 1
        spin_nb = 1
        deg = 1
        ediff =  0.0
        theta = 0.0
        idx = 1

        @fastmath @inbounds begin
            for istep in eachindex(a_schedule)
                a_field = a_schedule[istep]
                b_field = b_schedule[istep]
                randperm!(perm)
                for idx in eachindex(perm)
                    ispin = perm[idx]
                    ediff = 0.0
                    theta = pi * rand()
                    z_mag = cos(theta) - cos(spin_vector[ispin])
                    # add z components
                    for deg = 1:maxnb
                        spin_nb = nbs_idxs[deg, ispin]
                        if spin_nb === ispin
                            ediff += b_field*nbs_vals[deg, ispin]*z_mag
                        elseif spin_nb !== 0
                            ediff += b_field*z_mag*(nbs_vals[deg, ispin]*cos(spin_vector[spin_nb]))
                        end
                    end
                    # add x components
                    ediff += a_field * (sin(spin_vector[ispin]) - sin(theta))
                    if ediff <= 0.0
                        spin_vector[ispin] = theta
                    elseif exp(-1.0 * ediff / temp) > rand()
                        spin_vector[ispin] = theta
                    end
                end
            end
        end
        # return spin_vector
    end

    function annealiq!(spin_vector::Array{Float64, 1},
                    nbs_idxs::Array{Int64, 2},
                    nbs_vals::Array{Float64, 2},
                    a_schedule::Array{Float64, 2},
                    # must be same size or will crash from @inbounds
                    b_schedule::Array{Float64, 2},
                    temp::Float64)
        maxnb = size(nbs_idxs, 1)
        nspins = size(spin_vector, 1)
        perm = collect(1:nspins)
        num_sweeps = size(a_schedule, 2)
        istep = 1
        z_mag = 0.0
        a_field = 0.0
        b_field = 0.0
        sweep = 1
        ispin = 1
        spin_nb = 1
        deg = 1
        ediff =  0.0
        theta = 0.0
        idx = 1

        @fastmath @inbounds begin
            for istep = 1:num_sweeps
                randperm!(perm)
                for idx in eachindex(perm)
                    ispin = perm[idx]
                    a_field = a_schedule[ispin, istep]
                    b_field = b_schedule[ispin, istep]
                    ediff = 0.0
                    theta = pi * rand()
                    z_mag = cos(theta) - cos(spin_vector[ispin])
                    # add z components
                    for deg = 1:maxnb
                        spin_nb = nbs_idxs[deg, ispin]
                        if spin_nb === ispin
                            ediff += b_field*nbs_vals[deg, ispin]*z_mag
                        elseif spin_nb !== 0
                            ediff += b_field*z_mag*(nbs_vals[deg, ispin]*cos(spin_vector[spin_nb]))
                        end
                    end
                    # add x components
                    ediff += a_field * (sin(spin_vector[ispin]) - sin(theta))
                    if ediff <= 0.0
                        spin_vector[ispin] = theta
                    elseif exp(-1.0 * ediff / temp) > rand()
                        spin_vector[ispin] = theta
                    end
                end
            end
        end
        # return spin_vector
    end

    function anneals!(spin_vector::Array{Float64, 1},
                    nbs_idxs::Array{Int64, 2},
                    nbs_vals::Array{Float64, 2},
                    a_schedule::Array{Float64, 1},
                    # must be same size or will crash from @inbounds
                    b_schedule::Array{Float64, 1},
                    temp::Float64)
        # this svmc is in a sphere rather than a restricted semi-circle
        maxnb = size(nbs_idxs, 1)
        nspins = size(spin_vector, 1)
        perm = collect(1:nspins)
        phi_vector = zeros(Float64, nspins)
        istep = 1
        z_mag = 0.0
        x_mag = 0.0
        a_field = 0.0
        b_field = 0.0
        sweep = 1
        ispin = 1
        spin_nb = 1
        deg = 1
        ediff =  0.0
        theta = 0.0
        phi = 0.0
        idx = 1

        @fastmath @inbounds begin
            for istep in eachindex(a_schedule)
                a_field = a_schedule[istep]
                b_field = b_schedule[istep]
                randperm!(perm)
                for idx in eachindex(perm)
                    ispin = perm[idx]
                    ediff = 0.0
                    theta = pi * rand()
                    phi = 2 * pi * rand()
                    z_mag = cos(theta) - cos(spin_vector[ispin])
                    x_mag = (sin(spin_vector[ispin])*cos(phi_vector[ispin])) - (sin(theta)*cos(phi))
                    # add z components
                    for deg = 1:maxnb
                        spin_nb = nbs_idxs[deg, ispin]
                        if spin_nb === ispin
                            ediff += b_field*nbs_vals[deg, ispin]*z_mag
                        elseif spin_nb !== 0
                            ediff += b_field*z_mag*(nbs_vals[deg, ispin]*cos(spin_vector[spin_nb]))
                        end
                    end
                    # add x components
                    ediff += a_field * x_mag
                    if ediff <= 0.0
                        spin_vector[ispin] = theta
                        phi_vector[ispin] = phi
                    elseif exp(-1.0 * ediff / temp) > rand()
                        spin_vector[ispin] = theta
                        phi_vector[ispin] = phi
                    end
                end
            end
        end
        # return spin_vector
    end

    function annealsiq!(spin_vector::Array{Float64, 1},
                    nbs_idxs::Array{Int64, 2},
                    nbs_vals::Array{Float64, 2},
                    a_schedule::Array{Float64, 2},
                    # must be same size or will crash from @inbounds
                    b_schedule::Array{Float64, 2},
                    temp::Float64)
        # this svmc is in a sphere rather than a restricted semi-circle
        maxnb = size(nbs_idxs, 1)
        nspins = size(spin_vector, 1)
        perm = collect(1:nspins)
        phi_vector = zeros(Float64, nspins)
        num_sweeps = size(a_schedule, 2)
        istep = 1
        z_mag = 0.0
        x_mag = 0.0
        a_field = 0.0
        b_field = 0.0
        sweep = 1
        ispin = 1
        spin_nb = 1
        deg = 1
        ediff =  0.0
        theta = 0.0
        phi = 0.0
        idx = 1

        @fastmath @inbounds begin
            for istep in 1:num_sweeps
                randperm!(perm)
                for idx in eachindex(perm)
                    ispin = perm[idx]
                    a_field = a_schedule[ispin, istep]
                    b_field = b_schedule[ispin, istep]
                    ediff = 0.0
                    theta = pi * rand()
                    phi = 2 * pi * rand()
                    z_mag = cos(theta) - cos(spin_vector[ispin])
                    x_mag = (sin(spin_vector[ispin])*cos(phi_vector[ispin])) - (sin(theta)*cos(phi))
                    # add z components
                    for deg = 1:maxnb
                        spin_nb = nbs_idxs[deg, ispin]
                        if spin_nb === ispin
                            ediff += b_field*nbs_vals[deg, ispin]*z_mag
                        elseif spin_nb !== 0
                            ediff += b_field*z_mag*(nbs_vals[deg, ispin]*cos(spin_vector[spin_nb]))
                        end
                    end
                    # add x components
                    ediff += a_field * x_mag
                    if ediff <= 0.0
                        spin_vector[ispin] = theta
                        phi_vector[ispin] = phi
                    elseif exp(-1.0 * ediff / temp) > rand()
                        spin_vector[ispin] = theta
                        phi_vector[ispin] = phi
                    end
                end
            end
        end
        # return spin_vector
    end

    function sample(h::Dict, J::Dict; num_reads=1000, temp=0.1,
                    # for individual qubit schedules, input a
                    # dictionary of schedules
                    a_sch=1.0 .- collect(LinRange(0, 1, 2000)),
                    b_sch=collect(LinRange(0, 1, 2000)),
                    add_xtalk=false, noisy=false,
                    spherical=false, initial_state=Array{Float64, 1}())

        nspins = length(h)

        # check to see if we have individual qubit schedules
        independent_qubit_sechdules = false
        if isa(a_sch, Dict) | isa(b_sch, Dict)
            # if one input is not a dict but one if, convert to dict of lists
            if !isa(a_sch, Dict)
                schedule = deepcopy(a_sch)
                a_sch = Dict([(i, schedule) for i = 1:nspins])
            elseif !isa(b_sch, Dict)
                schedule = deepcopy(b_sch)
                b_sch = Dict([(i, schedule) for i = 1:nspins])
            end

            # check consistent inputs
            @assert length(a_sch) == nspins && length(b_sch) == nspins
            num_steps = length(a_sch[1])
            @assert all([length(a_sch[i]) == num_steps for i = 1:nspins])
            @assert all([length(b_sch[i]) == num_steps for i = 1:nspins])

            independent_qubit_sechdules = true
            a_sch_all = zeros(Float64, nspins, num_steps)
            b_sch_all = zeros(Float64, nspins, num_steps)
            for i = 1:nspins
                a_sch_all[i, :] = a_sch[i]
                b_sch_all[i, :] = b_sch[i]
            end
        end

        # prepare initial state if defined
        if length(initial_state) == nspins
            spin_vectors = [acos.(deepcopy(initial_state)) for i =1:num_reads]
        else
            # starting in the transverse field direction as default
            spin_vectors = fill(0.5 * pi, nspins, num_reads)
            spin_vectors = [spin_vectors[:, i] for i =1:num_reads]
        end
        spin_biases = Array{Float64, 1}([h[i] for i =1:nspins])
        JM, M, degs = hj_info(h, J)
        Js = copy(JM)
        nbs_idxs, nbs_vals = generate_nbr(nspins, M, maximum(degs) + 1)

        if independent_qubit_sechdules
            if spherical
                fn = annealsiq!
            else
                fn = annealiq!
            end

            if !noisy & !add_xtalk
                @inbounds for i = 1:num_reads
                    fn(spin_vectors[i], nbs_idxs, nbs_vals, a_sch_all, b_sch_all, temp)
                end
            elseif noisy & add_xtalk
                @inbounds for i = 1:num_reads
                    hc, jc = add_cross_talk(add_noise(h, J)...)
                    hj_info!(Js, M, degs, hc, jc)
                    generate_nbr!(nbs_idxs, nbs_vals, nspins, M, maximum(degs) + 1)
                    fn(spin_vectors[i], nbs_idxs, nbs_vals, a_sch_all, b_sch_all, temp)
                end
            elseif add_xtalk & !noisy
                hc, jc = add_cross_talk(h, J)
                hj_info!(Js, M, degs, hc, jc)
                generate_nbr!(nbs_idxs, nbs_vals, nspins, M, maximum(degs) + 1)
                @inbounds for i = 1:num_reads
                    fn(spin_vectors[i], nbs_idxs, nbs_vals, a_sch_all, b_sch_all, temp)
                end
            else
                @inbounds for i = 1:num_reads
                    hc, jc = add_noise(h, J)
                    hj_info!(Js, M, degs, hc, jc)
                    generate_nbr!(nbs_idxs, nbs_vals, nspins, M, maximum(degs) + 1)
                    fn(spin_vectors[i], nbs_idxs, nbs_vals, a_sch_all, b_sch_all, temp)
                end
            end
        else
            if spherical
                fn = anneals!
            else
                fn = anneal!
            end
            if !noisy & !add_xtalk
                @inbounds for i = 1:num_reads
                    fn(spin_vectors[i], nbs_idxs, nbs_vals, a_sch, b_sch, temp)
                end
            elseif noisy & add_xtalk
                @inbounds for i = 1:num_reads
                    hc, jc = add_cross_talk(add_noise(h, J)...)
                    hj_info!(Js, M, degs, hc, jc)
                    generate_nbr!(nbs_idxs, nbs_vals, nspins, M, maximum(degs) + 1)
                    fn(spin_vectors[i], nbs_idxs, nbs_vals, a_sch, b_sch, temp)
                end
            elseif add_xtalk & !noisy
                hc, jc = add_cross_talk(h, J)
                hj_info!(Js, M, degs, hc, jc)
                generate_nbr!(nbs_idxs, nbs_vals, nspins, M, maximum(degs) + 1)
                @inbounds for i = 1:num_reads
                    fn(spin_vectors[i], nbs_idxs, nbs_vals, a_sch, b_sch, temp)
                end
            else
                @inbounds for i = 1:num_reads
                    hc, jc = add_noise(h, J)
                    hj_info!(Js, M, degs, hc, jc)
                    generate_nbr!(nbs_idxs, nbs_vals, nspins, M, maximum(degs) + 1)
                    fn(spin_vectors[i], nbs_idxs, nbs_vals, a_sch, b_sch, temp)
                end
            end
        end
        return rotor_aggregation(spin_vectors, num_reads, collect(UpperTriangular(JM)))
    end

    export sample, anneal!, anneals!, annealiq!, annealsiq!
end

ext: jl

max_stars_repo_name: dtoconnor/Annealing-Samplers

max_issues_repo_name: dtoconnor/Annealing-Samplers

max_forks_repo_name: dtoconnor/Annealing-Samplers


text:
function specialize(llvm_ctx::LLVM.Context, ctx::ASTContext, template_decl::ClassTemplateDecl, args...)
    arg_vec = Vector{TemplateArgument}(undef, length(args))
    for (i, arg) in enumerate(args)
        if arg isa Bool
            jlty = typeof(arg)
            clty = jlty_to_clty(jlty, ctx)
            bitwidth_clty = getTypeSize(ctx, clty)
            # this assumes Bool is lowered to int8
            @assert bitwidth_clty == 8
            v = LLVM.GenericValue(jlty_to_llvmty(jlty, llvm_ctx), Int(arg))
            @assert LLVM.intwidth(v) == bitwidth_clty
            arg_vec[i] = TemplateArgument(ctx, v, clty)
            dispose(v)
        elseif arg isa Integer
            int_jlty = typeof(arg)
            int_clty = jlty_to_clty(int_jlty, ctx)
            bitwidth_clty = getTypeSize(ctx, int_clty)
            v = LLVM.GenericValue(jlty_to_llvmty(int_jlty, llvm_ctx), arg)
            @assert LLVM.intwidth(v) == bitwidth_clty
            arg_vec[i] = TemplateArgument(ctx, v, int_clty)
            dispose(v)
        elseif arg isa AbstractQualType
            arg_vec[i] = TemplateArgument(arg)
        else
            error("failed to specialize $arg")
        end
    end
    arg_list = TemplateArgumentList(ctx, arg_vec)
    specialization_decl = findSpecialization(template_decl, arg_list)

    if specialization_decl.ptr == C_NULL
        specialization_decl = ClassTemplateSpecializationDecl(ctx, template_decl, arg_list)
        AddSpecialization(template_decl, specialization_decl)
        if isOutOfLine(template_decl)
            lexical_decl_ctx = getLexicalDeclContext(template_decl)
            setLexicalDeclContext(specialization_decl, lexical_decl_ctx)
        end
    end

    # dispose.(arg_vec)  # FIXME: free NULL pointer

    return specialization_decl
end


text:
"""
    is_discete(k)

Returns, if the `AbstractKoopmanOperator` `k` is discrete in time.
"""
is_discrete(k::AbstractKoopmanOperator) = k.discrete

"""
    is_continouos(k)

Returns, if the `AbstractKoopmanOperator` `k` is continouos in time.
"""
is_continouos(k::AbstractKoopmanOperator) = !k.discrete

"""
    eigen(k)

    Return the eigendecomposition of the `AbstractKoopmanOperator`.
"""
LinearAlgebra.eigen(k::AbstractKoopmanOperator) = eigen(k.operator)

"""
    eigevals(k)

    Return the eigenvalues of the `AbstractKoopmanOperator`.
"""
LinearAlgebra.eigvals(k::AbstractKoopmanOperator) = eigvals(k.operator)

"""
    eigvecs(k)

    Return the eigenvectors of the `AbstractKoopmanOperator`.
"""
LinearAlgebra.eigvecs(k::AbstractKoopmanOperator) = eigvecs(k.operator)

"""
    modes(k)

    Return the eigenvectors of a continuous `AbstractKoopmanOperator`.
"""
modes(k::AbstractKoopmanOperator) = is_continouos(k) ? eigvecs(k) : throw(AssertionError("Koopman is discrete."))

"""
    frequencies(k)

    Return the eigenvalues of a continuous `AbstractKoopmanOperator`.
"""
frequencies(k::AbstractKoopmanOperator) = is_continouos(k) ? eigvals(k) : throw(AssertionError("Koopman is discrete."))

"""
    operator(k)

    Return the approximation of the discrete koopman operator stored in `k`.
"""
operator(k::AbstractKoopmanOperator) = is_discrete(k) ? k.operator : throw(AssertionError("Koopman is continouos."))

"""
    generator(k)

    Return the approximation of the continuous koopman generator stored in `k`.
"""
generator(k::AbstractKoopmanOperator) = is_continouos(k) ? k.operator : throw(AssertionError("Koopman is discrete."))

"""
    inputmap(k)

    Return the array `B`, mapping the exogenuos inputs to the koopman space.
"""
inputmap(k::AbstractKoopmanOperator) = k.input

"""
    outputmap(k)

    Return the array `C`, mapping the koopman space back onto the state space.
"""
outputmap(k::AbstractKoopmanOperator) = k.output

"""
    updateable(k)

    Returns `true` if the `AbstractKoopmanOperator` is updateable.
"""
updateable(k::AbstractKoopmanOperator) = !isempty(k.Q) && !isempty(k.P)

"""
    isstable(k)

    Returns `true` if either

    + the koopman operator has just eigenvalues with magnitude less than one or
    + the koopman generator has just eigenvalues with a negative real part
"""
isstable(k::AbstractKoopmanOperator) = is_discrete(k) ? all(abs.(eigvals(k)) .< one(eltype(k.operator))) : all(real.(eigvals(k)) < zero(eltype(k.operator)))


text:
using OrdinaryDiffEq
using Trixi

# flux_central for conservation and lax friedrichs for stability ?
# -> central flux doesnt even work for standard DGSEM

##########  EINSTELLUNGEN     #########
# Entropy conservation test with central surface flux and other initial conditions since it crashes with the blast wave
# initial_condition =  initial_condition_convergence_test
# surface_flux = flux_central

initial_condition = initial_condition_weak_blast_wave
CFL = 0.8
tspan = (0.0, 10)
N = 3
cells_per_dimension = (16, 16)

# CHandrashekar DGSEM Entropy STability
surface_flux = flux_chandrashekar # FluxPlusDissipation(flux_chandrashekar, DissipationLocalLaxFriedrichs(max_abs_speed_naive)) #
volume_flux  = flux_chandrashekar
volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
basis = LobattoLegendreBasis(N)
solver = DGSEM(basis, surface_flux, volume_integral)

# Standard DGSEM Entropy STability
# polydeg = N
# surface_flux = flux_lax_friedrichs
# volume_integral = VolumeIntegralWeakForm()
# solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = volume_integral )

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4)

coordinates_min = (-2.0, -2.0)
coordinates_max = ( 2.0,  2.0)

mesh = CurvedMesh(cells_per_dimension, coordinates_min, coordinates_max)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=CFL)
analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true,
                                                 extra_analysis_errors=(:conservation_error,),
                                                 extra_analysis_integrals=(entropy,)# energy_total,
                                                                           #energy_kinetic, energy_internal)

                                                )

callbacks = CallbackSet(summary_callback, stepsize_callback, analysis_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary

max_stars_repo_path: src/WR2/MS2 Entropy.jl

max_issues_repo_path: src/WR2/MS2 Entropy.jl

max_forks_repo_path: src/WR2/MS2 Entropy.jl


text:
import Base: rand,
             ∈,
             isempty

export Universe

"""
    Universe{N<:Real} <: AbstractPolyhedron{N}

Type that represents the universal set, i.e., the set of all elements.
"""
struct Universe{N<:Real} <: AbstractPolyhedron{N}
    dim::Int
end

isoperationtype(::Type{<:Universe}) = false

# default constructor of type Float64
Universe(dim::Int) = Universe{Float64}(dim)

# --- AbstractPolyhedron interface functions ---

"""
    constraints_list(U::Universe{N}) where {N<:Real}

Return the list of constraints defining a universe.

### Input

- `U` -- universe

### Output

The empty list of constraints, as the universe is unconstrained.
"""
function constraints_list(U::Universe{N}) where {N<:Real}
    return LinearConstraint{N, Vector{N}}[]
end

"""
    constrained_dimensions(U::Universe)::Vector{Int}

Return the indices in which a universe is constrained.

### Input

- `U` -- universe

### Output

The empty vector, as the universe is unconstrained in every dimension.
"""
function constrained_dimensions(U::Universe)::Vector{Int}
    return Int[]
end

# --- LazySet interface functions ---

"""
    dim(U::Universe)

Return the dimension of a universe.

### Input

- `U` -- universe

### Output

The dimension of a universe.
"""
function dim(U::Universe)::Int
    return U.dim
end

"""
    ρ(d::AbstractVector{N}, U::Universe{N}) where {N<:Real}

Return the support function of a universe.

### Input

- `d` -- direction
- `U` -- universe

### Output

The support function in the given direction.

### Algorithm

If the direction is all zero, the result is zero.
Otherwise, the result is `Inf`.
"""
function ρ(d::AbstractVector{N}, U::Universe{N}) where {N<:Real}
    return iszero(d) ? zero(N) : N(Inf)
end

"""
    σ(d::AbstractVector{N}, U::Universe{N}) where {N<:Real}

Return the support vector of a universe.

### Input

- `d` -- direction
- `U` -- universe

### Output

A vector with infinity values, except in dimensions where the direction is zero.
"""
function σ(d::AbstractVector{N}, U::Universe{N}) where {N<:Real}
    return [v == zero(N) ? v : v > zero(N) ? N(Inf) : N(-Inf) for v in d]
end

"""
    ∈(x::AbstractVector{N}, U::Universe{N})::Bool where {N<:Real}

Check whether a given point is contained in a universe.

### Input

- `x` -- point/vector
- `U` -- universe

### Output

The output is always `true`.

### Examples

```jldoctest
julia> [1.0, 0.0] ∈ Universe(2)
true
```
"""
function ∈(x::AbstractVector{N}, U::Universe{N})::Bool where {N<:Real}
    @assert length(x) == dim(U)
    return true
end

"""
    an_element(U::Universe{N}) where {N<:Real}

Return some element of a universe.

### Input

- `U` -- universe

### Output

The origin.
"""
function an_element(U::Universe{N}) where {N<:Real}
    return zeros(N, dim(U))
end

"""
    rand(::Type{Universe}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing
        )::Universe{N}

Create a universe (note that there is nothing to randomize).

### Input

- `Universe` -- type for dispatch
- `N`        -- (optional, default: `Float64`) numeric type
- `dim`      -- (optional, default: 2) dimension
- `rng`      -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`     -- (optional, default: `nothing`) seed for reseeding

### Output

The (only) universe of the given numeric type and dimension.
"""
function rand(::Type{Universe};
              N::Type{<:Real}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing
             )::Universe{N}
    rng = reseed(rng, seed)
    return Universe{N}(dim)
end

"""
    isempty(U::Universe)::Bool

Return if a universe is empty or not.

### Input

- `U` -- universe

### Output

`false`.
"""
function isempty(U::Universe)::Bool
    return false
end

"""
    isbounded(U::Universe)::Bool

Determine whether a universe is bounded.

### Input

- `U` -- universe

### Output

`false` as the universe is unbounded.
"""
function isbounded(U::Universe)::Bool
    return false
end

"""
    isuniversal(U::Universe{N}, [witness]::Bool=false
               )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}

Check whether a universe is universal.

### Input

- `U`       -- universe
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `true`
* If `witness` option is activated: `(true, [])`
"""
function isuniversal(U::Universe{N}, witness::Bool=false
                    )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
    return witness ? (true, N[]) : true
end

"""
    norm(U::Universe, [p]::Real=Inf)

Return the norm of a universe.
It is the norm of the enclosing ball (of the given ``p``-norm) of minimal volume
that is centered in the origin.

### Input

- `U` -- universe
- `p` -- (optional, default: `Inf`) norm

### Output

An error.
"""
function norm(U::Universe, p::Real=Inf)
    error("a universe does not have a norm")
end

"""
    radius(U::Universe, [p]::Real=Inf)

Return the radius of a universe.
It is the radius of the enclosing ball (of the given ``p``-norm) of minimal
volume with the same center.

### Input

- `U` -- universe
- `p` -- (optional, default: `Inf`) norm

### Output

An error.
"""
function radius(U::Universe, p::Real=Inf)
    error("a universe does not have a radius")
end

"""
    diameter(U::Universe, [p]::Real=Inf)

Return the diameter of a universe.
It is the maximum distance between any two elements of the set, or,
equivalently, the diameter of the enclosing ball (of the given ``p``-norm) of
minimal volume with the same center.

### Input

- `U` -- universe
- `p` -- (optional, default: `Inf`) norm

### Output

An error.
"""
function diameter(U::Universe, p::Real=Inf)
    error("a universe does not have a diameter")
end

"""
    translate(U::Universe{N}, v::AbstractVector{N}) where {N<:Real}

Translate (i.e., shift) a universe by a given vector.

### Input

- `U` -- universe
- `v` -- translation vector

### Output

The universe.
"""
function translate(U::Universe{N}, v::AbstractVector{N}) where {N<:Real}
    @assert length(v) == dim(U) "cannot translate a $(dim(U))-dimensional " *
                                "set by a $(length(v))-dimensional vector"
    return U
end


text:
abstract type Storage <: Injection end

include("storage/batteries.jl")


text:
"""
A module which provides wrappers around
`DifferentialEquations` solvers and
`AstrodynamicalModels` for orbit
propagation, iterative Halo orbit solvers,
and manifold calculations.

# Extended Help

**Exports**

$(EXPORTS)

**Imports**

$(IMPORTS)
"""
module Propagation

export Trajectory, ODEProblem, propagate
export Manifold, EnsembleProblem
export initialstate, initialepoch, solution
export isperiodic, halo, monodromy
export manifold, perturb
export stable_eigenvector, unstable_eigenvector

using ..States
using ..Calculations
using ..CoordinateFrames

import Dates: now
import LinearAlgebra: I
import Roots: find_zero

using Unitful
using Requires
using AstroTime
using Distributed
using StaticArrays
using LinearAlgebra
using DocStringExtensions
using AstrodynamicalModels
using DifferentialEquations

@template (FUNCTIONS, METHODS, MACROS) =
    """
    $(SIGNATURES)
    $(DOCSTRING)
    """

@template (TYPES, CONSTANTS) =
    """
    $(TYPEDEF)
    $(DOCSTRING)
    """

include("Trajectories.jl")
include("Propagators.jl")
include("Halos.jl")
include("Manifolds.jl")

function __init__()
    @require DataFrames="a93c6f00-e57d-5684-b7b6-d8193f3e46c0" begin
        DataFrames.DataFrame(traj::Trajectory) = DataFrames.DataFrame(solution(traj))
    end
end

end # module


text:
__precompile__(true)

module StateSpaceReconstruction

using Reexport
using StaticArrays
using LinearAlgebra

include("GroupSlices.jl")

include("Embeddings/Embeddings.jl")
include("Discretization/RectangularPartitions/RectangularPartitions.jl")
include("Discretization/TriangulationPartitions/TriangulationPartitions.jl")

# Rely on multiple dispatch for clashing function definitions in the submodules
dimension(E::AbstractEmbedding) = Embeddings.dimension(E)
dimension(s::T) where {T<:AbstractSimplex} = TriangulationPartitions.Simplices.dimension(s)
dimension(dt::T) where {T<:AbstractDelaunayTriangulation} =
    TriangulationPartitions.Delaunay.dimension(dt)
npoints(s::T) where {T<:AbstractSimplex} = TriangulationPartitions.Simplices.npoints(s)
npoints(E::AbstractEmbedding) = Embeddings.npoints(E)

export dimension, npoints

end # module

# #
# #
# # function estimator_exact(pts::Vector{Vector{T}}) where T
# #
# #     DT = StateSpaceReconstruction.DelaunayTriangulation(pts)
# #
# #     # Tolerance for zero volume (if one of the simplices have nearly zero
# #     # volume, skip computing the intersection volume and just set it to
# #     # zero instead).
# #     n_simplices = length(DT.indices)
# #     ϵ::Float64 = 1e-5/n_simplices
# #
# #     dim = length(DT[1]) - 1
# #
# #     image_simplex = StateSpaceReconstruction.MutableSimplex([zeros(Float64, dim) for i = 1:(dim+1)])
# #     simplex = StateSpaceReconstruction.MutableSimplex([zeros(Float64, dim) for i = 1:(dim+1)])
# #
# #     transfermatrix = zeros(Float64, n_simplices, n_simplices) # intersecting volumes
# #
# #     println("Building transfer matrix")
# #     for i in 1:n_simplices
# #         [image_simplex[k] = pts[DT[i][k]] for k = 1:(dim + 1)]
# #         for j in 1:n_simplices
# #             [simplex[k] = pts[DT[j][k]] for j = 1:(dim + 1)]
# #
#             @show image_simplex
#         end
#     end
#
# end
#
#
#
# @time estimator_exact([rand(3) for i = 1:10])
#
#
#
# TO = zeros(Float64, n_simplices, n_simplices) # intersecting volumes
#
# for i in 1:n_simplices
#     imvol = t.volumes_im[i]
#     for j in 1:n_simplices
#         vol = t.volumes[j]
#         if vol * imvol > 0 && (vol/imvol) > ϵ
#             # Intersecting volume between these two simplices
#             TO[i, j] = simplexintersection(
#                             transpose(t.points[t.simplex_inds[j, :], :]),
#                             transpose(t.impoints[t.simplex_inds[i, :], :])
#                         ) / imvol
#         end
#     end


text:
using DimensionalData, Test
using DimensionalData: X, Y
using BenchmarkTools

for N ∈ (10, 100, 1000)
    x = X(range(1; length = N))
    y = Y(range(1; length = N))
    d = (x, y)
    A = DimArray(rand(length.(d)...), d)

    println("N=$(N). Time to access via X(1:3)")
    @btime $(A)[$(X(1:3))];
end

println("Time to access if X(1:3) is not interpolated constant")
@btime $(A)[X(1:3)];
println("Time to just access A.data")
@btime $(A)[1:3, :]
@btime $(A)[1, 1]

# println("Time to do A+A")
# @btime $(A) .+ $(A);
# println("Time to do A.data + A.data")
# @btime $(A.data) .+ $(A.data);


text:
export dice_image

"""
    dice_image(sheet::Image, rows, cols)
Dice (or split) and image made up of many smaller images into multiple separate images
based on the `rows` and `cols` we have provided as arguments. This is useful because
game assets often come as sprite sheets containing multiple images of the same character
of different tiles.
"""
function dice_image(sheet::Image, rows::Integer, cols::Integer)
    sheet_width, sheet_height = size(sheet)
    w, h = sheet_width/cols, sheet_height/rows

    frames = Image[]
    for row in 1:rows
        y0 = row*h
        for col in 1:cols
            x0 = col*w
            push!(frames, Image[x0:x0+w, y0:y0+h])
        end
    end

    frames
end


text:
include("absorption_line.jl")
include("continuum.jl")

using QuadGK

#export AbstractSpectrum, SimpleSpectrum, SimulatedSpectrum, ConvolvedSpectrum

abstract type AbstractSpectrum end

"Wrapper to turn an arbitrary function of one variable into a sub-type of AbstractSpectrum."
struct SimpleSpectrum <: AbstractSpectrum
    spectrum::Function
    function SimpleSpectrum(f::Function)
        new(f)
    end
end

"Evaluate simple spectrum at wavelength(s) x."
function (s::SimpleSpectrum)(x::Number)
    s.spectrum(x)
end

function (s::SimpleSpectrum)(x::A1) where { T1<:Number, A1<:AbstractArray{T1} }
    s.spectrum(x)
end

"Simulated spectrum consists of stellar lines, telluric lines, and optionally continuum and Doppler shift."
struct SimulatedSpectrum{T1,T2,T3,T4,T5,T6,T7,T8,T9} <: AbstractSpectrum
    star_line_locs::Array{T1,1}
    star_line_widths::Array{T2,1}
    star_line_depths::Array{T3,1}
    telluric_line_locs::Array{T4,1}
    telluric_line_widths::Array{T5,1}
    telluric_line_depths::Array{T6,1}
    continuum_param::Array{T7,1}
    z::T8
    lambda_mid::T9
    limit_line_effect::T9

    function SimulatedSpectrum(star_line_locs::Array{T1,1}, star_line_widths::Array{T2,1}, star_line_depths::Array{T3,1},
        telluric_line_locs::Array{T4,1}, telluric_line_widths::Array{T5,1}, telluric_line_depths::Array{T6,1};
        continuum_param::Array{T7,1} = [1.0], z::T8 = 0.0, lambda_mid::T9 = 0.0, limit_line_effect::T9 = Inf
        ) where {T1<:Number, T2<:Number, T3<:Number, T4<:Number, T5<:Number, T6<:Number, T7<:Number, T8<:Number, T9<:Number}
        @assert length(star_line_locs) == length(star_line_widths) == length(star_line_depths)
        @assert length(telluric_line_locs) == length(telluric_line_widths) == length(telluric_line_depths)
        @assert length(continuum_param) >= 1
        new{T1,T2,T3,T4,T5,T6,T7,T8,T9}(star_line_locs,star_line_widths,star_line_depths,
            telluric_line_locs,telluric_line_widths,telluric_line_depths,
            continuum_param,z,lambda_mid,limit_line_effect)
    end
end

#=
"Simulated spectrum consists of stellar lines, telluric lines, and optionally continuum and Doppler shift."
struct SimulatedSpectrumGpu{T1,T2,T3,T4,T5,T6,T7,T8,T9} <: AbstractSpectrum
    star_line_locs::CuArray{T1,1}
    star_line_widths::CuArray{T2,1}
    star_line_depths::CuArray{T3,1}
    telluric_line_locs::CuArray{T4,1}
    telluric_line_widths::CuArray{T5,1}
    telluric_line_depths::CuArray{T6,1}
    continuum_param::CuArray{T7,1}
    z::T8
    lambda_mid::T9
    limit_line_effect::T9

    function SimulatedSpectrumGpu(star_line_locs::Array{T1,1}, star_line_widths::Array{T2,1}, star_line_depths::Array{T3,1},
        telluric_line_locs::Array{T4,1}, telluric_line_widths::Array{T5,1}, telluric_line_depths::Array{T6,1};
        continuum_param::Array{T7,1} = [1.0], z::T8 = 0.0, lambda_mid::T9 = 0.0, limit_line_effect::T9 = Inf
        ) where {T1<:Number, T2<:Number, T3<:Number, T4<:Number, T5<:Number, T6<:Number, T7<:Number, T8<:Number, T9<:Number}
        @assert length(star_line_locs) == length(star_line_widths) == length(star_line_depths)
        @assert length(telluric_line_locs) == length(telluric_line_widths) == length(telluric_line_depths)
        @assert length(continuum_param) >= 1
        new{T1,T2,T3,T4,T5,T6,T7,T8,T9}(cu(star_line_locs),cu(star_line_widths),cu(star_line_depths),
            cu(telluric_line_locs),cu(telluric_line_widths),cu(telluric_line_depths),
            cu(continuum_param),z,lambda_mid,limit_line_effect)
    end
end
=#

"Evaluate simulated spectrum at wavelength(s) x."
function (spectrum::AbstractSpectrum)(x::Number)
    continuum_model(x*(one(spectrum.z)+spectrum.z),spectrum.continuum_param,center=spectrum.lambda_mid) *
            absorption_lines(x*(one(spectrum.z)+spectrum.z), spectrum.star_line_locs, spectrum.star_line_widths, spectrum.star_line_depths, limit_line_effect=spectrum.limit_line_effect ) *
            absorption_lines(x, spectrum.telluric_line_locs, spectrum.telluric_line_widths, spectrum.telluric_line_depths, limit_line_effect=spectrum.limit_line_effect )
end

function (spectrum::AbstractSpectrum)(x::A) where { T<:Number, A<:AbstractArray{T} }
    #@inbounds spectrum.(x)
    f_cont(y) = continuum_model(y*(one(spectrum.z)+spectrum.z),spectrum.continuum_param,center=spectrum.lambda_mid)
    f_star(y) = absorption_lines(y*(one(spectrum.z)+spectrum.z), spectrum.star_line_locs, spectrum.star_line_widths, spectrum.star_line_depths, limit_line_effect=spectrum.limit_line_effect )
    f_telluric(y) = absorption_lines(y, spectrum.telluric_line_locs, spectrum.telluric_line_widths, spectrum.telluric_line_depths, limit_line_effect=spectrum.limit_line_effect )
    @inbounds f_cont.(x) .* f_star.(x) .* f_telluric.(x)
end

abstract type AbstractConvolutionKernel end

"Spectrum based on conlving another spectrum with a point spread function.  Optionally, limit the width of the convolution."
struct ConvolvedSpectrum{T1,T2} <: AbstractSpectrum
    spectrum::T1
    kernel::T2
    limit_conv_width::Float64
    function ConvolvedSpectrum(spectrum::T1, kernel::T2;
         limit_conv_width = limit_kernel_width_default) where {T1<:AbstractSpectrum, T2<:AbstractConvolutionKernel}
        new{T1,T2}(spectrum,kernel,limit_conv_width)
    end
end

"Evaluate the convolved spectrum a wavelenth(s) x"
function (s::ConvolvedSpectrum)(x::Number)
    quadgk(y -> s.spectrum(x+y) * s.kernel(y),
                -s.limit_conv_width,s.limit_conv_width)[1]
end

# Inefficient since integration allocates arrays rather than scalars
function (s::ConvolvedSpectrum)(x::A) where { T<:Number, A<:AbstractArray{T} }
    @inbounds quadgk(y -> s.spectrum(x.+y) .* s.kernel.(y),
                -s.limit_conv_width,s.limit_conv_width)[1]
end

"Create same simulated spectrum but replacing the Doppler shift parameter (z)."
function doppler_shifted_spectrum(s::AbstractSpectrum, z::Number) end

function doppler_shifted_spectrum(s::SimulatedSpectrum, z::Number)
    SimulatedSpectrum(s.star_line_locs, s.star_line_widths,s.star_line_depths,
            s.telluric_line_locs,s.telluric_line_widths,s.telluric_line_depths,
            continuum_param=s.continuum_param, z=z,
            lambda_mid=s.lambda_mid, limit_line_effect=s.limit_line_effect)
end

function doppler_shifted_spectrum(s::ConvolvedSpectrum{T1,T2}, z::Number) where {
                T1<:AbstractSpectrum, T2<:AbstractConvolutionKernel }
    ConvolvedSpectrum(doppler_shifted_spectrum(s.spectrum,z), s.kernel,
            limit_conv_width=s.limit_conv_width)
end


text:
function setup_routine(p::Param)

    #
    println("\n \n")
    println("___________________________________________________________")
    println(" _  __ _____    _____   ____   _  __      __ ______  _____ ")
    println("| |/ // ____|  / ____| / __ \\ | | \\ \\    / /|  ____||  __ \\ ")
    println("| ' /| |  __  | (___  | |  | || |  \\ \\  / / | |__   | |__) |")
    println("|  < | | |_ |  \\___ \\ | |  | || |   \\ \\/ /  |  __|  |  _  / ")
    println("| . \\| |__| |  ____) || |__| || |____\\  /   | |____ | | \\ \\ ")
    println("|_|\\_\\\\_____| |_____/  \\____/ |______|\\/    |______||_|  \\_\\")
    println("___________________________________________________________")
    println(" ")
    println(" Welcome to the Klein gordon field Solver")
    println("___________________________________________________________")

    #Print Parameters
    print_param_info(p)

    println(" ")
    println("Start up routine:")

    #Folder setup
    if( !set_up_folder(p) )
        return false
    end

    return true

end


text:
# ABIF Reader
# ===========

export get_tags, tagelements

struct AbifDirEntry
    name::String
    number::Int32
    element_type::Int32
    element_size::Int32
    num_elements::Int32
    data_size::Int32
    data_offset::Int32
end

"""
    ABIF.Reader(input::IO)
Create a data reader of the ABIF file format.
# Arguments
* `input`: data source
"""
mutable struct Reader{T<:IO} <: BioCore.IO.AbstractReader
    # input stream
    input::T

    # Tags
    dirs::Vector{AbifDirEntry}
end

function Reader(input::IO)
    header = read_abif_header(input)
    tags   = read_abif_tags(input, header)

    return Reader(input, tags)
end

function BioCore.IO.stream(reader::Reader)
    return reader.input
end

function Base.iterate(p::Reader, k::Int=1)
    if k > length(p)
        return nothing
    else
        return getindex(p, p.dirs[k]), k + 1
    end
end

function Base.length(a::Reader)
    return length(a.dirs)
end

function Base.checkbounds(p::Reader, k::Integer)
    if 1 ≤ k ≤ length(p)
        return true
    end
    throw(BoundsError(p, k))
end

function Base.getindex(a::Reader, k::Integer)
    checkbounds(a, k)
    return Dict([parse_data_tag(a, a.dirs[k])])
end

function Base.getindex(a::Reader, t::AbstractString)
    tag = filter((x) -> x.name == t, a.dirs)

    if isempty(tag)
        throw(KeyError(t))
    end

    if length(tag) > 1
        return Dict([parse_data_tag(a, b) for b in tag])
    end

    return Dict([parse_data_tag(a, first(tag))])
end

function Base.getindex(a::Reader, t::AbifDirEntry)
    return Dict([parse_data_tag(a, t)])
end

function Base.getindex(a::Reader, t::Array{AbifDirEntry})
    return Dict([parse_data_tag(a, k) for k in t])
end

"""
    get_tags(input::Reader)
Returns all existing tags
# Arguments
* `input`: Reader
"""
function get_tags(a::Reader)
    return [tag for tag in a.dirs]
end

"""
    get_tags(input::Reader, tag_name::AbstractString)
Returns all existing tags by name
# Arguments
* `input`: Reader
* `tag_name`: AbstractString
"""
function get_tags(a::Reader, t::AbstractString)
    return [tag for tag in a.dirs if isequal(tag.name, t)]
end

"""
    tagelements(stream::Reader, tag_name::AbstractString)
Returns the number of how many tags exists by the same name
# Arguments
* `stream`: Reader
* `tag_name`: AbstractString
"""
function tagelements(a::Reader, t::AbstractString)
    return length(get_tags(a, t))
end

# extract the header
function read_abif_header(input::IO)
    seekstart(input)
    signature = read(input, 4)

    if is_abif_signature(signature)
        version = ntoh(first(read(input, Int16)))
        return parse_directory(input, position(input))
    else
        error("Invalid File Signature")
    end
end

# extract all tags in file
function read_abif_tags(input::IO, header::AbifDirEntry)
    tags = AbifDirEntry[]
    for index in 0:header.num_elements-1
        start = header.data_offset + index * header.element_size
        tag = parse_directory(input, start)
        push!(tags, tag)
    end
    return tags
end

# extract DirEtry
function parse_directory(input::IO, pos::Int64)
    seek(input, pos)

    name         = String(read(input, 4))
    number       = ntoh(read(input, Int32))
    element_type = ntoh(read(input, UInt16))
    element_size = ntoh(read(input, Int16))
    num_elements = ntoh(read(input, Int32))
    data_size    = ntoh(read(input, Int32))

    if data_size <= 4
        data_offset = position(input)
    else
        data_offset = ntoh(read(input, Int32))
    end

    AbifDirEntry(name, number, element_type, element_size, num_elements, data_size, data_offset)
end

# read bytes according to the element type, other values are unsupported or legacy.
function parse_data_tag(a::Reader, tag::AbifDirEntry)
    if tag.data_offset > 0
        seek(a.input, tag.data_offset)

        data = Tuple{String, Union{AbstractString, Integer}}

        if tag.element_type == 2
            data = String(read(a.input, tag.data_size))

        elseif tag.element_type == 4
            data = [read(a.input, Int16) for _ in 1:tag.num_elements]
            data = convert_to_int(data)

        elseif tag.element_type == 5
            data = [read(a.input, Int32) for _ in 1:tag.num_elements]
            data = convert_to_int(data)

        elseif tag.element_type == 7
            data = [read(a.input, Float32) for _ in 1:tag.num_elements]
            data = convert_to_float(data)

        elseif tag.element_type == 10
            year   = ntoh(read(a.input, Int16))
            month  = Int(read(a.input, UInt8))
            day    = Int(read(a.input, UInt8))
            data   = "$year-$month-$day"

        elseif tag.element_type == 11
            hour    = Int(ntoh(read(a.input, UInt8)))
            minute  = Int(ntoh(read(a.input, UInt8)))
            second  = Int(ntoh(read(a.input, UInt8)))
            hsecond = Int(ntoh(read(a.input, UInt8)))
            data    = "$hour:$minute:$second:$hsecond"

        elseif tag.element_type == 13
            data = Bool(first(read(a.input, tag.data_size)))

        elseif tag.element_type == 18
            data = String(read(a.input, tag.data_size)[2:end])

        elseif tag.element_type == 19
            data = String(read(a.input, tag.data_size)[1:end-1])

        else
            data = read(a.input, tag.data_size)
        end

        return format_data!(data, tag, tagelements(a, tag.name))
    end
end

# if TAG has more than one element, concatenate element number on name
function format_data!(data::Any, tag::AbifDirEntry, elements::Integer)
    if elements > 1
        return ("$(tag.name)$(tag.number)", data)
    end
    return ("$(tag.name)", data)
end

# convert a array of big-endian values
function convert_to_int(data::AbstractArray)
    result = Int[]

    for d in data
        r = Int(ntoh(d))
        push!(result, r)
    end

    if length(result) <= 1
        return first(result)
    end

    return result
end

# convert a array of big-endian values
function convert_to_float(data::AbstractArray)
    result = Float32[]
    for d in data
        r = Float32(ntoh(d))
        push!(result, r)
    end

     if length(result) <= 1
        return first(result)
     end

    return result
end

# Check if the file has a valid Signature
function is_abif_signature(signature::Array{UInt8})
    if signature != b"ABIF"
        return false
    end
    return true
end


text:
# Advent of code 2021, day 1

function day01()
    part = [0, 0]
    input1 = map(s -> parse(Int, s), readlines("AoCdata/AoC_2021_day1.txt"))

    increases(arr) = count(arr[i + 1] > arr[i] for i in 1:(length(arr) - 1))
    part[1] = increases(input1) # 1477

    function slidingincreases(arr)
        return count(
            sum(arr[(i + 1):(i + 3)]) > sum(arr[i:(i + 2)]) for i in 1:(length(arr) - 3)
        )
    end
    part[2] = slidingincreases(input1) # 1523
    return part
end

part = day01()

println("Part 1: ", part[1])
println("Part 2: ", part[2])


text:
using Test
using LogicCircuits
using ProbabilisticCircuits

@testset "pr_constraint Query" begin
    # two nodes
    simplevtree = zoo_vtree_file("simple2.vtree")
    pc, vtree = load_struct_prob_circuit(
                    zoo_psdd_file("simple2.4.psdd"), simplevtree)

    cache = Dict{Tuple{ProbΔNode, Union{ProbΔNode, StructLogicalΔNode}}, Float64}()

    @test abs(pr_constraint(pc[end], pc[end], cache) - 1.0) < 1e-8
    @test abs(pr_constraint(pc[5], pc[3], cache) - 0.2) < 1e-8
    @test abs(pr_constraint(pc[5], pc[4], cache) - 0.8) < 1e-8

    file_circuit = "little_4var.circuit"
    file_vtree = "little_4var.vtree"
    logical_circuit, vtree = load_struct_smooth_logical_circuit(
                                zoo_lc_file(file_circuit), zoo_vtree_file(file_vtree))

    pc = zoo_psdd("little_4var.psdd")

    @test abs(pr_constraint(pc[end], logical_circuit[end - 1], cache) - 1.0) < 1e-8

    # Test with two psdds
    pc1, vtree = load_struct_prob_circuit(zoo_psdd_file("simple2.5.psdd"), simplevtree)
    pc2, vtree = load_struct_prob_circuit(zoo_psdd_file("simple2.6.psdd"), simplevtree)

    pr_constraint_cache = Dict{Tuple{ProbΔNode, Union{ProbΔNode, StructLogicalΔNode}}, Float64}()
    pr_constraint(pc1[end], pc2[end], pr_constraint_cache)
    @test abs(pr_constraint_cache[pc1[1], pc2[1]] - 1.0) < 1e-8
    @test abs(pr_constraint_cache[pc1[1], pc2[2]] - 0.0) < 1e-8
    @test abs(pr_constraint_cache[pc1[3], pc2[4]] - 1.0) < 1e-8
    @test abs(pr_constraint_cache[pc1[3], pc2[5]] - 0.0) < 1e-8
    @test abs(pr_constraint_cache[pc1[9], pc2[8]] - 1.0) < 1e-8
    @test abs(pr_constraint_cache[pc1[5], pc2[4]] - 0.2) < 1e-8
    @test abs(pr_constraint_cache[pc1[5], pc2[5]] - 0.8) < 1e-8
    @test abs(pr_constraint_cache[pc1[2], pc2[3]] - 1.0) < 1e-8
end


text:
using Documenter
using GLM, StatsBase, SparseArrays, LinearAlgebra
using RobustModels

DocMeta.setdocmeta!(RobustModels, :DocTestSetup, :(using RobustModels, StatsBase, GLM); recursive=true)

prettyurls = get(ENV, "CI", "false") == "true"

makedocs(
    modules = [RobustModels, GLM, StatsBase],
    sitename = "RobustModels",
    authors="Bertrand Lacoste <bertrand.lacoste@gmail.com>",
    repo="https://github.com/getzze/RobustModels.jl/blob/{commit}{path}#{line}",
    format = Documenter.HTML(;
        prettyurls=prettyurls,
        canonical="https://getzze.github.io/RobustModels.jl",
        assets = ["assets/favicon.ico"],
    ),
    pages = [
        "Home" => "index.md",
        "manual.md",
        "examples.md",
        "api.md",
    ],
#    doctest = :fix,
    debug = false,
)

deploydocs(
    devbranch = "main",
    repo = "github.com/getzze/RobustModels.jl"
)


text:
dS = 0.01

# b = r gives Black Scholes on Stocks
# b = r - q gives stock option with dividend yield q
# b = 0 give Black76
# b and r = 0 gives margined Future Options
# b = r - rf give Garman Kohlhagen currency options

function Delta(CallPutFlag,  S,  X,  T,  r,  b,  v)

        return (
        (GBlackScholes(CallPutFlag, S + dS, X, T, r, b, v) -
        GBlackScholes(CallPutFlag, S - dS, X, T, r, b, v)) / (2 * dS)
        )

end

function Elasticity(CallPutFlag,  S,  X,  T,  r,  b,  v)

        return ((GBlackScholes(CallPutFlag, S + dS, X, T, r, b, v) -
                GBlackScholes(CallPutFlag, S - dS, X, T, r, b, v)) / (2 * dS)
                * S / GBlackScholes(CallPutFlag, S, X, T, r, b, v)
                )
end

function Gamma(CallPutFlag,  S,  X,  T,  r,  b,  v)

        return ((GBlackScholes(CallPutFlag, S + dS, X, T, r, b, v)
                - 2 * GBlackScholes(CallPutFlag, S, X, T, r, b, v) +
                GBlackScholes(CallPutFlag, S - dS, X, T, r, b, v)) / dS^2
                )

end

function DGammaDVol(CallPutFlag, S, X, T, r, b, v)

        return ((GBlackScholes(CallPutFlag, S + dS, X, T, r, b, v + dS) -
                2 * GBlackScholes(CallPutFlag, S , X, T, r, b, v + dS) +
                GBlackScholes(CallPutFlag, S - dS, X, T, r, b, v + dS) -
                GBlackScholes(CallPutFlag, S + dS, X, T, r, b, v - dS) +
                2 * GBlackScholes(CallPutFlag, S, X, T, r, b, v - dS) -
                GBlackScholes(CallPutFlag, S - dS, X, T, r, b, v - dS)) /
                (2 * dS * dS^2) / 100
                )

end

function GammaP(CallPutFlag, S, X, T, r, b, v)

        return (S / 100 * (GBlackScholes(CallPutFlag, S + dS, X, T, r, b, v) -
                2 * GBlackScholes(CallPutFlag, S , X, T, r, b, v) +
                GBlackScholes(CallPutFlag, S - dS, X, T, r, b, v)) / dS^2
                )

end

function DDeltaDVol(CallPutFlag, S, X, T, r, b, v)

        return (
                1/(4 * dS * dS) * (GBlackScholes(CallPutFlag, S + dS, X, T, r, b, v + dS)
                - GBlackScholes(CallPutFlag, S + dS, X, T, r, b, v - dS)
                - GBlackScholes(CallPutFlag, S - dS, X, T, r, b, v + dS)
                + GBlackScholes(CallPutFlag, S - dS, X, T, r, b, v - dS)) / 100
                )

end

function Vega(CallPutFlag, S, X, T, r, b, v)

        return (
                (GBlackScholes(CallPutFlag, S, X, T, r, b, v + dS) -
                GBlackScholes(CallPutFlag, S, X, T, r, b, v - dS)) / 2
                )

end

function Vomma(CallPutFlag, S, X, T, r, b, v)

        return (
                (GBlackScholes(CallPutFlag, S, X, T, r, b, v + dS)
                - 2 * GBlackScholes(CallPutFlag, S, X, T, r, b, v)
                + GBlackScholes(CallPutFlag, S, X, T, r, b, v - dS)) / dS^2 / 10000
                )
end

function VegaP(CallPutFlag, S, X, T, r, b, v)

        return (
                v / 0.1 * (GBlackScholes(CallPutFlag, S, X, T, r, b, v + dS) -
                GBlackScholes(CallPutFlag, S, X, T, r, b, v - dS)) / 2
                )

end

function DVegaDVol(CallPutFlag, S, X, T, r, b, v)

        return (
                (GBlackScholes(CallPutFlag, S, X, T, r, b, v + dS) -
                2 * GBlackScholes(CallPutFlag, S, X, T, r, b, v) +
                GBlackScholes(CallPutFlag, S, X, T, r, b, v - dS))
                )

end

function Theta(CallPutFlag, S, X, T, r, b, v)

        if T <= (1/365)

                return (
                        GBlackScholes(CallPutFlag, S, X, 1e-05, r, b, v) -
                        GBlackScholes(CallPutFlag, S, X, T, r, b, v)
                        )
        else

                return (
                        GBlackScholes(CallPutFlag, S, X, T - (1/365), r, b, v) -
                        GBlackScholes(CallPutFlag, S, X, T, r, b, v)
                        )
        end

end

function Rho(CallPutFlag, S, X, T, r, b, v)

        return (
                (GBlackScholes(CallPutFlag, S, X, T, r + dS, b + dS, v)
                - GBlackScholes(CallPutFlag, S, X, T, r - dS, b - dS, v)) / 2
                )

end

function RhoFut(CallPutFlag,  S,  X,  T,  r,  b,  v)

        return (
                (GBlackScholes(CallPutFlag, S, X, T, r + dS, b, v) -
                GBlackScholes(CallPutFlag, S, X, T, r - dS, b, v)) / 2
                )

end

function Rho2(CallPutFlag, S, X, T, r, b, v)

        return (
                (GBlackScholes(CallPutFlag, S, X, T, r, b - dS, v) -
                GBlackScholes(CallPutFlag, S, X, T, r, b + dS, v)) / 2
                )

end

function Carry(CallPutFlag, S, X, T, r, b, v)

        return (
                (GBlackScholes(CallPutFlag, S, X, T, r, b + dS, v) -
                GBlackScholes(CallPutFlag, S, X, T, r, b - dS, v)) / 2
                )

end

function Speed(CallPutFlag, S, X, T, r, b, v)

        return (
                1 / dS^3 * (GBlackScholes(CallPutFlag, S + 2 * dS, X, T, r, b, v) -
                3 * GBlackScholes(CallPutFlag, S + dS, X, T, r, b, v) +
                3 * GBlackScholes(CallPutFlag, S, X, T, r, b, v) -
                GBlackScholes(CallPutFlag, S - dS, X, T, r, b, v))
                )

end

function StrikeDelta(CallPutFlag, S, X, T, r, b, v)

        return (
                (GBlackScholes(CallPutFlag, S, X + dS, T, r, b, v) -
                GBlackScholes(CallPutFlag, S, X - dS, T, r, b, v)) / (2 * dS)
                )

end

function RNDensity(CallPutFlag, S, X, T, r, b, v)

        return (
                (GBlackScholes(CallPutFlag, S, X + dS, T, r, b, v) -
                2 * GBlackScholes(CallPutFlag, S, X, T, r, b, v) +
                GBlackScholes(CallPutFlag, S, X - dS, T, r, b, v)) / dS^2
                )

end

# function CrossGamma(CallPutFlag, S1, S2, X, T, r, b1, b2, v1, v2, cor)
#
#         return {
#                 (KirkSpreadOption(CallPutFlag, S1, S2 + dS, X, T, r, b1, b2, v1, v2, cor)
#                 - 2 * KirkSpreadOption(CallPutFlag, S1, S2, X, T, r, b1, b2, v1, v2, cor) +
#                 KirkSpreadOption(CallPutFlag, S1, S2 - dS, X, T, r, b1, b2, v1, v2, cor)) / dS^2
#                 }
#
# end


text:
module LinearSolvers

using ..MPIStateArrays

using StaticArrays, LinearAlgebra

using KernelAbstractions

# just for testing LinearSolvers
LinearAlgebra.norm(A::MVector, p::Real, weighted::Bool) = norm(A, p)
LinearAlgebra.norm(A::MVector, weighted::Bool) = norm(A, 2, weighted)
LinearAlgebra.dot(A::MVector, B::MVector, weighted) = dot(A, B)
LinearAlgebra.norm(A::AbstractVector, p::Real, weighted::Bool) = norm(A, p)
LinearAlgebra.norm(A::AbstractVector, weighted::Bool) = norm(A, 2, weighted)
LinearAlgebra.dot(A::AbstractVector, B::AbstractVector, weighted) = dot(A, B)

export linearsolve!, settolerance!, prefactorize
export AbstractLinearSolver, AbstractIterativeLinearSolver

"""
    AbstractLinearSolver

This is an abstract type representing a generic linear solver.
"""
abstract type AbstractLinearSolver end

"""
    AbstractIterativeLinearSolver

This is an abstract type representing a generic iterative
linear solver.

The available concrete implementations are:

  - [`GeneralizedConjugateResidual`](@ref)
  - [`GeneralizedMinimalResidual`](@ref)
"""
abstract type AbstractIterativeLinearSolver <: AbstractLinearSolver end

"""
    settolerance!(solver::AbstractIterativeLinearSolver, tolerance, relative)

Sets the relative or absolute tolerance of the iterative linear solver
`solver` to `tolerance`.
"""
settolerance!(
    solver::AbstractIterativeLinearSolver,
    tolerance,
    relative = true,
) = (relative ? (solver.rtol = tolerance) : (solver.atol = tolerance))

doiteration!(
    linearoperator!,
    Q,
    Qrhs,
    solver::AbstractIterativeLinearSolver,
    tolerance,
    args...,
) = throw(MethodError(
    doiteration!,
    (linearoperator!, Q, Qrhs, solver, tolerance, args...),
))

initialize!(
    linearoperator!,
    Q,
    Qrhs,
    solver::AbstractIterativeLinearSolver,
    args...,
) = throw(MethodError(initialize!, (linearoperator!, Q, Qrhs, solver, args...)))

"""
    prefactorize(linop!, linearsolver, args...)

Prefactorize the in-place linear operator `linop!` for use with `linearsolver`.
"""
prefactorize(linop!, linearsolver::AbstractIterativeLinearSolver, args...) =
    linop!

"""
    linearsolve!(linearoperator!, solver::AbstractIterativeLinearSolver, Q, Qrhs, args...)

Solves a linear problem defined by the `linearoperator!` function and the state
`Qrhs`, i.e,

```math
L(Q) = Q_{rhs}
```

using the `solver` and the initial guess `Q`. After the call `Q` contains the
solution.  The arguments `args` is passed to `linearoperator!` when it is
called.
"""
function linearsolve!(
    linearoperator!,
    solver::AbstractIterativeLinearSolver,
    Q,
    Qrhs,
    args...;
    max_iters = length(Q),
    cvg = Ref{Bool}(),
)
    converged = false
    iters = 0

    converged, threshold =
        initialize!(linearoperator!, Q, Qrhs, solver, args...)
    converged && return iters

    while !converged && iters < max_iters
        converged, inner_iters, residual_norm =
            doiteration!(linearoperator!, Q, Qrhs, solver, threshold, args...)

        iters += inner_iters

        if !isfinite(residual_norm)
            error("norm of residual is not finite after $iters iterations of `doiteration!`")
        end

        achieved_tolerance = residual_norm / threshold * solver.rtol
    end

    converged || @warn "Solver did not attain convergence after $iters iterations"
    cvg[] = converged

    iters
end

@kernel function linearcombination!(Q, cs, Xs, increment::Bool)
    i = @index(Global, Linear)
    if !increment
        @inbounds Q[i] = -zero(eltype(Q))
    end
    @inbounds for j in 1:length(cs)
        Q[i] += cs[j] * Xs[j][i]
    end
end

end


text:
using ProximalOperators
using ProximalAlgorithms
using BenchmarkTools
using LinearAlgebra
using SparseArrays
using Random

const SUITE = BenchmarkGroup()

for T in [Float64]
    k = "Lasso ($T)"
    SUITE[k] = BenchmarkGroup([k])
    A = T[
        1.0 -2.0 3.0 -4.0 5.0
        2.0 -1.0 0.0 -1.0 3.0
        -1.0 0.0 4.0 -3.0 2.0
        -1.0 -1.0 -1.0 1.0 3.0
    ]
    b = T[1.0, 2.0, 3.0, 4.0]
    m, n = size(A)
    R = real(T)
    lam = R(0.1) * norm(A' * b, Inf)

    SUITE[k]["ForwardBackward"] = @benchmarkable solver(x0, f=f, A=$A, g=g) setup=begin
        solver = ProximalAlgorithms.ForwardBackward(tol=1e-4)
        x0 = zeros($T, size($A, 2))
        f = Translate(SqrNormL2(), -$b)
        g = NormL1($lam)
    end

    SUITE[k]["FastForwardBackward"] = @benchmarkable solver(x0, f=f, A=$A, g=g) setup=begin
        solver = ProximalAlgorithms.ForwardBackward(tol=1e-4, fast=true)
        x0 = zeros($T, size($A, 2))
        f = Translate(SqrNormL2(), -$b)
        g = NormL1($lam)
    end

    SUITE[k]["ZeroFPR"] = @benchmarkable solver(x0, f=f, A=$A, g=g) setup=begin
        solver = ProximalAlgorithms.ZeroFPR(tol=1e-4)
        x0 = zeros($T, size($A, 2))
        f = Translate(SqrNormL2(), -$b)
        g = NormL1($lam)
    end

    SUITE[k]["PANOC"] = @benchmarkable solver(x0, f=f, A=$A, g=g) setup=begin
        solver = ProximalAlgorithms.PANOC(tol=1e-4)
        x0 = zeros($T, size($A, 2))
        f = Translate(SqrNormL2(), -$b)
        g = NormL1($lam)
    end

    SUITE[k]["DouglasRachford"] = @benchmarkable solver(x0, f=f, g=g) setup=begin
        solver = ProximalAlgorithms.DouglasRachford(tol=1e-4, gamma=$R(1))
        x0 = zeros($T, size($A, 2))
        f = LeastSquares($A, $b)
        g = NormL1($lam)
    end

    SUITE[k]["AFBA-1"] = @benchmarkable solver(x0, y0, f=f, g=g, beta_f=beta_f) setup=begin
        beta_f = opnorm($A)^2
        solver = ProximalAlgorithms.AFBA(theta=$R(1), mu=$R(1), tol=$R(1e-6))
        x0 = zeros($T, size($A, 2))
        y0 = zeros($T, size($A, 2))
        f = LeastSquares($A, $b)
        g = NormL1($lam)
    end

    SUITE[k]["AFBA-2"] = @benchmarkable solver(x0, y0, h=h, L=$A, g=g) setup=begin
        beta_f = opnorm($A)^2
        solver = ProximalAlgorithms.AFBA(theta=$R(1), mu=$R(1), tol=$R(1e-6))
        x0 = zeros($T, size($A, 2))
        y0 = zeros($T, size($A, 1))
        h = Translate(SqrNormL2(), -$b)
        g = NormL1($lam)
    end

    SUITE[k]["SFISTA"] = @benchmarkable solver(y0, f=f, Lf=Lf, h=h) setup=begin
        solver = ProximalAlgorithms.SFISTA(tol=$R(1e-3))
        y0 = zeros($T, size($A, 2))
        f = LeastSquares($A, $b)
        h = NormL1($lam)
        Lf = opnorm($A)^2
    end
end


text:
immutable LDLt{T,S<:AbstractMatrix{T}} <: Factorization{T}
    data::S
end

size(S::LDLt) = size(S.data)
size(S::LDLt, i::Integer) = size(S.data, i)

convert{T,S}(::Type{LDLt{T,S}}, F::LDLt) = LDLt{T,S}(convert(S, F.data))
# NOTE: the annotaion <:AbstractMatrix shouldn't be necessary, it is introduced
#       to avoid an ambiguity warning (see issue #6383)
convert{T,S,U<:AbstractMatrix}(::Type{LDLt{T}}, F::LDLt{S,U}) = convert(LDLt{T,U}, F)

convert{T,S,U}(::Type{Factorization{T}}, F::LDLt{S,U}) = convert(LDLt{T,U}, F)

# SymTridiagonal
function ldltfact!{T<:Real}(S::SymTridiagonal{T})
	n = size(S,1)
	d = S.dv
	e = S.ev
	@inbounds for i = 1:n-1
		e[i] /= d[i]
		d[i+1] -= abs2(e[i])*d[i]
		d[i+1] > 0 || throw(PosDefException(i+1))
	end
	return LDLt{T,SymTridiagonal{T}}(S)
end
function ldltfact{T}(M::SymTridiagonal{T})
	S = typeof(one(T)/one(T))
	return S == T ? ldltfact!(copy(M)) : ldltfact!(convert(SymTridiagonal{S}, M))
end
function A_ldiv_B!{T}(S::LDLt{T,SymTridiagonal{T}}, B::AbstractVecOrMat{T})
	n, nrhs = size(B, 1), size(B, 2)
	size(S,1) == n || throw(DimensionMismatch(""))
	d = S.data.dv
	l = S.data.ev
	@inbounds begin
		for i = 2:n
			li1 = l[i-1]
			for j = 1:nrhs
				B[i,j] -= li1*B[i-1,j]
			end
		end
		dn = d[n]
		for j = 1:nrhs
			B[n,j] /= dn
		end
		for i = n-1:-1:1
			di = d[i]
			li = l[i]
			for j = 1:nrhs
				B[i,j] /= di
				B[i,j] -= li*B[i+1,j]
			end
		end
	end
	return B
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
__precompile__()
module CategoricalArrays
    export CategoricalPool, CategoricalValue
    export AbstractCategoricalArray, AbstractCategoricalVector, AbstractCategoricalMatrix,
           CategoricalArray, CategoricalVector, CategoricalMatrix
    export AbstractNullableCategoricalArray, AbstractNullableCategoricalVector,
           AbstractNullableCategoricalMatrix,
           NullableCategoricalArray, NullableCategoricalVector, NullableCategoricalMatrix
    export LevelsException

    export categorical, compress, decompress, droplevels!, levels, levels!, isordered, ordered!
    export cut, recode, recode!

    using Compat

    include("typedefs.jl")

    include("buildfields.jl")

    include("pool.jl")
    include("value.jl")

    include("array.jl")
    include("nullablearray.jl")
    include("subarray.jl")

    include("extras.jl")
    include("recode.jl")

    include("deprecated.jl")
end


text:
"""
    ForwardDiff.derivative(f,z::Complex)

Compute the derivative of function `f` with respect to `z` and `conj(z)`
"""
@inline function derivative(f::F, z::C) where {F,C<:Complex}
    T = typeof(ForwardDiff.Tag(f, C))
    # making f output complex ensures that it gets dispatched to our extract_derivative
    # for complex Duals rather than the native one in ForwardDiff
    return extract_derivative(T,complex.(f(one(ComplexDual{T},z))))
end

"""
    ForwardDiff.extract_derivative(T,d::ComplexDual)

Given a complex dual value `d` and tag `T`, extract the derivatives d/dz
and d/dz* from the partials in `d`.
"""
@inline function extract_derivative(::Type{T},d::ComplexDual{T}) where T
   return dz_partials(d,1)
end

# If the function outputs a tuple of duals, as would be the case for
# nested derivatives, then we deal with each one by one and assemble an overall
# list of derivatives
@inline function extract_derivative(::Type{T},dlist::NTuple{N,ComplexDual{T}}) where {N,T}
   out = ()
   for d in dlist
        out = (out...,extract_derivative(T,d)...)
    end
    out
end


text:
export Degree, Variable, show, iscomplex, isreal, isbool

"""
    Degree(explvar::Int, conjvar::Int)

Define a mathematical degree, which is used to define a variable exponent to
`explvar` and its conjugate exponent to `conjvar`.

### Attributes
- `explvar`: Int
- `conjvar`: Int
"""
mutable struct Degree
    explvar::Int
    conjvar::Int
end

#############################
## Equality, hash
#############################
==(d1::Degree, d2::Degree) = (d1.explvar==d2.explvar) && (d1.conjvar==d2.conjvar)
!=(d1::Degree, d2::Degree) = !(d1 == d2)
hash(deg::Degree, h::UInt) = hash(deg.explvar, hash(deg.conjvar, h))

"""
Variable(varname::String, Complex)

Define a mathematical variable `varname` of a certain mathematical kind,
`Complex` here.

### Attributes
- `name`: String
- `kind`: a type among `Complex`, `Real` and `Bool`
"""
mutable struct Variable <: AbstractPolynomial
    name::String
    kind::Type

    function Variable(name, kind)
        validtypes = Set([Complex, Real, Bool])
        if true ∉ map(x-> kind <: x, validtypes) || typeof(name) ∉ Set([String, SubString{String}])
            error(LOGGER, "Variable() : attempting to define a variable $name of type $kind, supported types are {Complex, Real, Bool}")
        end
        return new(String(name), kind)
    end
end

#############################
## Equality, hash
#############################
==(x::Variable, y::Variable) = (x.name==y.name) && (x.kind==y.kind)
!=(x::Variable, y::Variable) = !(x == y)
hash(var::Variable, h::UInt) = hash(var.name, hash(var.kind, h))

#############################
## Type detection utils
#############################
iscomplex(x::Variable) = x.kind <: Complex
isreal(x::Variable) = x.kind <: Real
isbool(x::Variable) = x.kind <: Bool

#############################
## Show
#############################
function Base.show(io::IO, x::Variable)
    print(io, x.name)
end

function Base.show(io::IO, d::Degree)
    print(io, "(", d.explvar, ",", d.conjvar, ")")
end


text:
module PlutoUtils

using Comonicon
using Pluto
using ExprTools
using Pluto.Configuration: Options, CompilerOptions, ServerOptions, from_flat_kwargs
using Pluto: open_in_default_browser, Notebook

include("patch.jl")

# all the following has a prefix JULIA_PLUTO
const PLUTO_ENV_VARIABLES = [
    "PROJECT",
    "WORKSPACE_USE_DISTRIBUTED",
    "RUN_NOTEBOOK_ON_LOAD",
    "PLUTO_WORKING_DIRECTORY",
]

_get_kw_name(x::Symbol) = x
function _get_kw_name(ex::Expr)
    if ex.head === :kw
        return _get_kw_name(ex.args[1])
    elseif ex.head === :(::)
        return _get_kw_name(ex.args[1])
    else
        error("invalid kwargs expression: $ex")
    end
end

function read_env(env::AbstractDict = ENV)
    for key in PLUTO_ENV_VARIABLES
        if haskey(env, "JULIA_PLUTO_" * key)
            return env[key]
        end
    end
end

# both `run` and `open` can accepet these options
const CLI_OPTIONS = [
    Expr(:kw, :(host::String), "127.0.0.1"),
    Expr(:kw, :(port::Int), 1234),
    Expr(:kw, :(launch_browser::Bool), true),
    Expr(:kw, :(project::String), "@."),
    Expr(:kw, :sysimage, nothing),
]

@static if VERSION > v"1.5.0-"
    push!(CLI_OPTIONS, Expr(:kw, :threads, nothing))
end

const CLI_OPTION_NAMES = map(_get_kw_name, CLI_OPTIONS)

const CLI_OPTION_DOCS = """
- `--host <ip>`: IP address of the host.
- `-p, --port <int>`: port you want to specify.
- `--launch-browser <bool>`: launch browser automatically or not.
- `--project <path>`: specify a project environment to run.
- `--sysimage <path>`: specify the path of a system image to use.
$(if VERSION > v"1.5.0-"
"- `-t, --threads <int>`: specify number of threads."
end
)
"""

defs = Dict{Symbol, Any}()
defs[:name] = :open
defs[:args] = [:(file::String)]
defs[:kwargs] = CLI_OPTIONS
defs[:body] = quote
    isfile(file) || error("file $file does not exist!")
    # from_flat_kwargs(;host=host, port=port, launch_browser=false, ...)
    options = $(Expr(:call, :from_flat_kwargs, Expr(:parameters, [Expr(:kw, each, each) for each in CLI_OPTION_NAMES]...)))
    session = Pluto.ServerSession(;options=options)
    # we overwrite the notebook options temporarily since
    # this session will be used for this notebook only
    nb = Pluto.SessionActions.open(session, file; compiler_options=nothing)
    return run_with_notebook_open(session, nb)
end

@eval begin
    """
    open a Pluto notebook at given path.

    # Args

    - `file`: file path of the Pluto notebook.

    # Options

    $(CLI_OPTION_DOCS)
    """
    @cast $(combinedef(defs))
end

defs = Dict{Symbol, Any}()
defs[:name] = :run
defs[:args] = []
defs[:kwargs] = CLI_OPTIONS
defs[:body] = quote
    options = $(Expr(:call, :from_flat_kwargs, Expr(:parameters, [Expr(:kw, each, each) for each in CLI_OPTION_NAMES]...)))
    s = Pluto.ServerSession(;options=options)
    return Pluto.run(s)
end

@eval begin
    """
    start Pluto notebook server.

    # Options

    $(CLI_OPTION_DOCS)
    """
    @cast $(combinedef(defs))
end

"Pluto CLI - Lightweight reactive notebooks for Julia"
@main

end


text:
using Pkg
Pkg.add("QuadGK")

using QuadGK

func1(x) = exp(-x^2)
res, err = quadgk(func1, -Inf, Inf)

abs(res - sqrt(π)) / sqrt(π)

res, err = quadgk(func1, -Inf, Inf, rtol = 1e-15)

abs(res - sqrt(π)) / sqrt(π)

res, err = quadgk(func1, -Inf, Inf, order = 12)

abs(res - sqrt(π)) / sqrt(π)

sqrt(π)

#%%

res, err = quadgk(func1, BigFloat(-Inf), BigFloat(Inf), rtol = 1e-15, order=12)
abs(res - sqrt(BigFloat(π))) / sqrt(BigFloat(π))
#%%
func2(x, y, z) = x + y^3 + sin(z)

x = 5
z = 3
arg(y) = func2(x, y, z)

quadgk(arg, 1, 3)

quadgk(y -> func2(x, y, z), 1, 3)

res, err = let x=5; z=3
    arg(y) = func2(x, y, z)

    quadgk(arg, 1, 3)
end

#%%
func3(x, y) = x^2 * exp(y)

function test_int(x, ymin, ymax)
    arg(y) = func3(x, y)
    return quadgk(arg, ymin, ymax)[1]
end

test_int(3, 1, 5)

max_stars_repo_name: JeffreySarnoff/techytok-examples

max_issues_repo_name: JeffreySarnoff/techytok-examples

max_forks_repo_name: JeffreySarnoff/techytok-examples


text:
## ------------------------------------------------------------------
function _openbook(dir0::String; reload = false)

    # If cache missed
    book = currbook()
    nipass = 1
    if isnothing(book)

        # bookdir
        bdir = _find_bookdir(dir0)
        isempty(bdir) && error("No '$(_rbook_toml_name())' file found in current or any parent directory!!!")

        # new book
        book = RBook()
        bookdir!(book, bdir)
        currbook!(book)

        nipass = 3
    end

    # Update book
    _include_rbfiles!(book; force = reload, nipass)

    return book
end
openbook(;reload = false) = _openbook(pwd(); reload)

# TODO: reload force option
macro openbook()
    src = srcfile(__source__)
    src = isfile(src) ? src : pwd()
    book = _openbook(dirname(src); reload = false)
    return book
end


text:
@with_kw struct aparameters{R}
  z_ref = 2;
  z_H = 1e-3;
  z_oH = 1e-3;
  z_m = 1e-3;
  z_om = 1e-3;
  d = 0;
end

global atm_parameters=aparameters{Float64}()