text:
"""
`module  JAC.HighHarmonic
    ... a submodel of JAC that contains all methods to set-up and process high-harmonic computations.
        It defines in particular a number of structs to deal with the observables, approaches, laser pulses and
        time-dependent fields as they frequently occur in the computation of high-harmonic spectra.
"""
module HighHarmonic

    using Printf, ..Basics, ..Defaults, ..Radial, ..ManyElectron, ..Nuclear, ..TableStrings

    """
    `abstract type HighHarmonic.AbstractTargetModel`
        ... defines an abstract and a number of other types to represent an atomic target (cloud) for HHG:

        + struct LocalizedTargetModel     ... to represent a point-like localized target
        + struct GaussianTargetModel      ... to represent a localized, Gaussian-type target (not yet)
        + struct PlanarTargetModel        ... to represent a planar target of finite extent (not yet)
        + struct UserTargetModel          ... to represent a user-defined 3-dimensional target (not yet)
    """
    abstract type  AbstractTargetModel  end

    """
    `struct  HighHarmonic.LocalizedTargetModel  <:  AbstractTargetModel`
        ... defines a struct for point-like localized target cloud.

        + center     ::CartesianVector{Float64}   ... Center of the (point-like) target.
    """
    struct  LocalizedTargetModel  <:  AbstractTargetModel
        center       ::CartesianVector{Float64}
    end

    """
    `struct  HighHarmonic.GaussianTargetModel  <:  AbstractTargetModel`
        ... defines a struct for localized, Gaussian target cloud.

        + center     ::CartesianVector{Float64}   ... Center of the (point-like) target.
        + NoPoints   ::Int64                      ... Number of points to be used for the target.
        + sigma      ::Float64                    ... Gaussian half-widths of the cloud.
    """
    struct  GaussianTargetModel  <:  AbstractTargetModel
        center       ::CartesianVector{Float64}
        NoPoints     ::Int64
        sigma        ::Float64
    end

    """
    `struct  HighHarmonic.TargetCloud`
        ... defines a struct to represent a -- point-like or extended -- target cloud for the computation of high-harmonic spectra

        + model      ::AbstractTargetModel                    ... Model of the target cloud, from which the target is generated.
        + points     ::Array{WeightedCartesian{Float64},1}    ... Cartesian points and weights that define the atomic target cloud.
    """
    struct  TargetCloud
        model        ::AbstractTargetModel
        points       ::Array{WeightedCartesian{Float64},1}
    end

    """
    `abstract type HighHarmonic.AbstractDeterctorModel`
        ... defines an abstract and a number of other types to represent an dectector screen for high-harmonic spectra:

      + struct LocalizedDetectorModel      ... to represent a point-like localized detector
      + struct SphericalDetectorModel      ... to represent a partial-sphere detector screen (not yet)
      + struct PlanarDetectorModel         ... to represent a planar detector screen of finite extent (not yet)
    """
    abstract type  AbstractDetectorModel  end

    """
    `struct  HighHarmonic.LocalizedDetectorModel  <:  AbstractDetectorModel`
        ... defines a struct for point-like localized detector.

        + center     ::CartesianVector{Float64}   ... Center of the (point-like) detector.
    """
    struct  LocalizedDetectorModel  <:  AbstractDetectorModel
        center       ::CartesianVector{Float64}
    end

    """
    `struct  HighHarmonic.SphericalDetectorModel  <:  AbstractDetectorModel`
        ... defines a struct for a partly-spherical detector screen.

        + center     ::CartesianVector{Float64}   ... Center of the detector.
        + NoPoints   ::Int64                      ... Number of points to be used for the detector screen.
        + theta      ::Float64                    ... polar (arc) angle of the screen.
    """
    struct  SphericalDetectorModel  <:  AbstractDetectorModel
        center       ::CartesianVector{Float64}
        NoPoints     ::Int64
        theta        ::Float64
    end

    """
    `struct  HighHarmonic.DetectorScreen`
        ... defines a struct to represent a -- point-like or extended -- detector screen for the computation of high-harmonic spectra.

        + model      ::AbstractDetectorModel                    ... Model of the detector, from which the detection screen is generated.
        + points     ::Array{WeightedCartesian{Float64},1}      ... Cartesian points and weights that define the detection screen.
    """
    struct  DetectorScreen
        model        ::AbstractDetectorModel
        points       ::Array{WeightedCartesian{Float64},1}
    end

    """
    `abstract type HighHarmonic.AbstractHhgObservable`
        ... defines an abstract and a number of singleton types for the observables that can be computed for HHG processes.

      + struct HhgSpectrum
        ... to compute a standard (single-color) high-harmonic spectrum (for a point-like target and point-like detector).

      + struct HhgPolarizedSpectrum
        ... to compute a bi-color high-harmonic spectrum (for a point-like target and point-like detector).

      + struct AttoSecondPulse     (not yet)
      + struct IntensityProfile    (not yet)
      + struct PhaseProfile        (not yet)

      + struct HhgSpectrum
    """
    abstract type  AbstractHhgObservable                              end
    struct         HhgSpectrum           <:  AbstractHhgObservable   end
    struct         HhgPolarizedSpectrum  <:  AbstractHhgObservable   end

    """
    `abstract type HighHarmonic.AbstractHhgApproach`
        ... defines an abstract and a number of singleton types for the computational approach/model that is applied in the
            computation of HHG spectra.

      + struct HhgHydrogenicSaddlePoint
        ... to apply hydrogenic orbitals for the computation of the time-dependent dipole amplitude and the saddle-point
            approximation for the dipole moments.

      + struct HhgHydrogenicFullVolkov
        ... to apply hydrogenic orbitals for the computation of the time-dependent dipole amplitude but the full
            Volkov phase for the dipole moments.
    """
    abstract type  AbstractHhgApproach                                 end
    struct         HhgHydrogenicSaddlePoint  <:  AbstractHhgApproach   end
    struct         HhgHydrogenicFullVolkov   <:  AbstractHhgApproach   end

    """
    `abstract type HighHarmonic.AbstractHhgPulse`
        ... defines an abstract and a number of other types to represent an incident laser pulse/field for HHG:

      + struct HhgLinearPlaneWaveLaser  ... to represent a long, linearly-polarized plane-wave laser field,
                                            propagating into z-direction.
      + struct HhgSineSquaredPulse      ... to represent a sin^2(...) pulse  (not yet)
    """
    abstract type  AbstractHhgPulse  end

    """
    `struct  HighHarmonic.HhgLinearPlaneWaveLaser  <:  AbstractHhgPulse`
        ... defines a struct for a long, plane-wave laser field of given frequency and intensity.

        + omega          ::Float64   ... Frequency of the incident laser field in [a.u.].
        + E0             ::Float64   ... Electric-field amplitude E0 in [a.u.].
        ...
        + intensity      ::Float64   ... Intensity [W/cm^2] of the incident laser field.
        + polarization   ::...       ... Polarization of the incident laser field.
        + beamProfile    ::...
        + envelope       ::...
        + length         ::...
        + centralLength  ::...
    """
    struct  HhgLinearPlaneWaveLaser  <:  AbstractHhgPulse
        omega            ::Float64
        E0               ::Float64
    end

    """
    `HighHarmonic.HhgLinearPlaneWaveLaser()`
        ... constructor for an (empty) instance of a linearly-polarized plane-wave pulse.
    """
    function HhgLinearPlaneWaveLaser()
        HhgLinearPlaneWaveLaser(0., 0.)
    end

    """
    `HighHarmonic.HhgLinearPlaneWaveLaser(omega; intensity::Float64=1.0e14)`
        ... constructor for an linearly-polarized plane-wave pulse of given frequency omega and intensity [W/cm^2].
    """
    function HhgLinearPlaneWaveLaser(omega; intensity::Float64=1.0e14)
        omega_au = Defaults.convertUnits("energy: to atomic", omega);   E0_au = intensity / 3.51e16
        HhgLinearPlaneWaveLaser(omega_au, E0_au)
    end

    # `Base.show(io::IO, pulse::HhgLinearPlaneWaveLaser)`  ... prepares a proper printout of the variable pulse::HhgLinearPlaneWaveLaser.
    function Base.show(io::IO, pulse::HhgLinearPlaneWaveLaser)
        println(io, "omega:              $(pulse.omega)  ")
        println(io, "E0:                 $(pulse.E0)  ")
    end

    """
    `struct  HighHarmonic.Computation`
        ... defines a type for the computation of high-harmonic spectra, atto-second pulses, intensity and phase profiles,
            and several others.

        + observable         ::AbstractHhgObservable    ... Specifies one of the observables to be calculated.
        + approach           ::AbstractHhgApproach      ... Specifies the approach/frame to be applied in the computations.
        + pulse              ::AbstractHhgPulse         ... Specifies the incident laser field.
        + timeMesh           ::Array{Float64,1}         ... time mesh to be used for representing the incident and generated fields.
        + nuclearModel       ::Nuclear.Model            ... Model, charge and parameters of the nucleus.
        + initialOrbital     ::Radial.Orbital           ... (Single-electron) orbital to specify the initial state of atom.
        + target             ::TargetCloud              ... Specifies the target for generating high harmonics.
        + detector           ::DetectorScreen           ... Specifies the detector for observing high harmonics.
    """
    struct  Computation
        observable           ::AbstractHhgObservable
        approach             ::AbstractHhgApproach
        pulse                ::AbstractHhgPulse
        timeMesh             ::Array{Float64,1}
        nuclearModel         ::Nuclear.Model
        initialOrbital       ::Radial.Orbital
        target               ::TargetCloud
        detector             ::DetectorScreen
    end

    # `Base.show(io::IO, computation::HighHarmonic.Computation)`  ... prepares a proper printout of the variable computation::HighHarmonic.Computation.
    function Base.show(io::IO, computation::HighHarmonic.Computation)
        println(io, "observable:           $(computation.observable)  ")
        println(io, "approach:             $(computation.approach)  ")
        println(io, "pulse:                $(computation.pulse)  ")
        println(io, "timeMesh:             $(computation.timeMesh)  ")
        println(io, "nuclearModel:         $(computation.nuclearModel)  ")
        println(io, "initialOrbital:       $(computation.initialOrbital)  ")
        println(io, "target:               $(computation.target)  ")
        println(io, "detector:             $(computation.detector)  ")
    end

    """
    `HighHarmonic.computeTimeMesh(NoCycles::Int64, NoPointsPerCycle::Int64)`
        ... to compute a time mesh for a photon field with given number of cycles (NoCycles) as well as time
            points per cycle; a timeMesh::Array{Float64,1} is returned.
    """
    function computeTimeMesh(NoCycles::Int64, NoPointsPerCycle::Int64)
        dt = 2pi / NoPointsPerCycle;    timeMesh = Float64[]
        for  k = 1:NoCycles*NoPointsPerCycle   push!( timeMesh, k*dt)   end

        println("Time mesh generated in the (time) interval [0., $(timeMesh[end])] with $(length(timeMesh)) points.")
        return( timeMesh )
    end

    """
    `HighHarmonic.computeTimeDipoleMoment(approach::HhgHydrogenicSaddlePoint, pulse::HhgLinearPlaneWaveLaser,
                                          timeMesh::Array{Float64,1}, orbital::Radial.Orbital)`
        ... to compute the time-dependent dipole moment D(t) for the given pulse, approximation and the orbital wave function.
            This dipole moment is calculated for all times as defined by timeMesh. A dipoleMoment::Array{Float64,1} is returned.
    """
    function computeTimeDipoleMoment(approach::HhgHydrogenicSaddlePoint, pulse::HhgLinearPlaneWaveLaser,
                                     timeMesh::Array{Float64,1}, orbital::Radial.Orbital)
        dipoleMoment = CartesianVector{ComplexF64}[]
        for  t in timeMesh
            push!( dipoleMoment, CartesianVector{ComplexF64}(0., 0., 0.) )
        end
        return( dipoleMoment )
    end

    """
    `HighHarmonic.computeElectricField(pulse::HhgLinearPlaneWaveLaser, time::Float64)`
        ... to compute the electric field for the given pulse at time; an eField::CartesianVector{Float64} is returned.
    """
    function computeElectricField(pulse::HhgLinearPlaneWaveLaser, time::Float64)
        return( CartesianVector{Float64}(0., 0., pulse.E0 * cos(pulse.omega*time) ) )
    end

    """
    `HighHarmonic.computeVectorPotential(pulse::HhgLinearPlaneWaveLaser, time::Float64)`
        ... to compute the electric field for the given pulse at time; a aField::CartesianVector{Float64} is returned.
    """
    function computeVectorPotential(pulse::HhgLinearPlaneWaveLaser, time::Float64)
        aField = CartesianVector{Float64}(0., 0., 0.)
        return( aField )
    end

    """
    `HighHarmonic.computeVolkovPhase(approach::HhgHydrogenicSaddlePoint, time::Float64, ip::Float64, aField::CartesianVector{Float64})`
        ... to compute the Volkov phase at time for the given approximation, ionization potential and vector potential (aField);
            a Volkov phase::Float64 is returned.
    """
    function computeVolkovPhase(approach::HhgHydrogenicSaddlePoint, time::Float64, ip::Float64, aField::CartesianVector{Float64})
        phase = 0.
        return( phase )
    end

    """
    `HighHarmonic.computeTimeDipoleAmplitude(approach::HhgHydrogenicSaddlePoint, time::Float64, aField::CartesianVector{Float64},
                                             orbital::Radial.Orbital)`
        ... to compute the dipole amplitude d(t) at time for the given approximation, vector potential (aField) and the
            orbital wave function; a dipole amplitude::CartesianVector{ComplexF64} is returned.
    """
    function computeTimeDipoleAmplitude(approach::HhgHydrogenicSaddlePoint, time::Float64, aField::Array{Float64,1},
                                        orbital::Radial.Orbital)
        dipoleAmplitude = CartesianVector{ComplexF64}(0., 0., 0.)
        return( dipoleAmplitude )
    end

    """
    `HighHarmonic.computeFrequencyDipoleMoment(timeMesh::Array{Float64,1}, timeDipole::Array{CartesianVector{ComplexF64},1})`
        ... to compute the frequency-dependent dipole moment D^~ (omega);
            an omegaMoment::Array{CartesianVector{ComplexF64},1} is returned.
    """
    function computeFrequencyDipoleMoment(timeMesh::Array{Float64,1}, timeDipole::Array{CartesianVector{ComplexF64},1})
        omegaMoment = CartesianVector{ComplexF64}[]
        for  t in timeMesh
           push!( omegaMoment, CartesianVector{ComplexF64}(0., 0., 0.) )
        end
        return( omegaMoment )
    end

    """
    `HighHarmonic.perform(comp::HighHarmonic.Computation; output::Bool=false)`
        ... to perform a computation of (one) selected observable in HHG, such a spectra, attosecond pulses of profiles.
            Different approaches, pulse structures of the incident light as well as different target clouds and detector
            screens can also be specified together with all physical parameters: frequency, intensity and polarization of
            the incident light; size andshape of target clouds; size and shape of detector screens; ....
    """
    function perform(comp::HighHarmonic.Computation; output::Bool=false)
        if  output    results = Dict{String, Any}()    else    results = nothing    end
        nModel = comp.nuclearModel

        if       comp.observable == HighHarmonic.HhgSpectrum()
            timeDipole  = HighHarmonic.computeTimeDipoleMoment(comp.approach, comp.pulse, comp.timeMesh, comp.initialOrbital)
            omegaDipole = HighHarmonic.computeFrequencyDipoleMoment(comp.timeMesh, timeDipole)
        elseif   comp.observable == HighHarmonic.HhgPolarizedSpectrum()
        else     error("Undefined observable for HHG computations.")
        end

        return( 0. )
    end

end # module


text:
using PyCall
import Plots

# Functionality for computing band structures, mostly using pymatgen

function high_symmetry_kpath(model; kline_density=20)
    bandstructure = pyimport("pymatgen.symmetry.bandstructure")
    pystructure = pymatgen_structure(model.lattice, model.atoms)
    symm_kpath = bandstructure.HighSymmKpath(pystructure)

    kcoords, labels = symm_kpath.get_kpoints(kline_density, coords_are_cartesian=false)

    labels_dict = Dict{String, Vector{eltype(kcoords[1])}}()
    for (ik, k) in enumerate(kcoords)
        if length(labels[ik]) > 0
            labels_dict[labels[ik]] = k
        end
    end

    (kcoords=kcoords, klabels=labels_dict, kpath=symm_kpath.kpath["path"])
end

function compute_bands(basis, ρ, kcoords, n_bands;
                       eigensolver=lobpcg_hyper, tol=1e-3, show_progress=true, kwargs...)
    # Create basis with new kpoints, where we cheat by using any symmetry operations.
    ksymops = [[identity_symop()] for _ in 1:length(kcoords)]
    # For some reason rationalize(2//3) isn't supported (julia 1.4)
    myrationalize(x::T) where {T <: AbstractFloat} = rationalize(x, tol=10eps(T))
    myrationalize(x) = x
    bs_basis = PlaneWaveBasis(basis, [myrationalize.(k) for k in kcoords], ksymops)
    ham = Hamiltonian(bs_basis; ρ=ρ)

    band_data = diagonalize_all_kblocks(eigensolver, ham, n_bands + 3;
                                        n_conv_check=n_bands,
                                        tol=tol, show_progress=show_progress, kwargs...)
    if !band_data.converged
        @warn "Eigensolver not converged" iterations=band_data.iterations
    end
    merge((basis=bs_basis, ), select_eigenpairs_all_kblocks(band_data, 1:n_bands))
end

function prepare_band_data(band_data; datakeys=[:λ, :λerror],
                           klabels=Dict{String, Vector{Float64}}())
    # Get pymatgen to compute kpoint distances and get it to split quantities
    # from the band_data object into nicely plottable branches
    # This is a bit of abuse of the routines in pymatgen, but it works ...
    plotter = pyimport("pymatgen.electronic_structure.plotter")

    ret = Dict{Symbol, Any}(:basis => band_data.basis)
    for key in datakeys
        hasproperty(band_data, key) || continue

        # Compute dummy "Fermi level" for pymatgen to be happy
        allfinite = [filter(isfinite, x) for x in band_data[key]]
        eshift = sum(sum, allfinite) / sum(length, allfinite)
        bs = pymatgen_bandstructure(band_data.basis, band_data[key], eshift, klabels)
        data = plotter.BSPlotter(bs).bs_plot_data(zero_to_efermi=false)

        # Check number of k-Points agrees
        @assert length(band_data.basis.kpoints) == sum(length, data["distances"])

        ret[:spins] = [:up]
        spinmap = [("1", :up)]
        if bs.is_spin_polarized
            ret[:spins] = [:up, :down]
            spinmap = [("1", :up), ("-1", :down)]
        end

        ret[:n_branches] = length(data["energy"])
        ret[:n_bands] = size(data["energy"][1]["1"], 1)
        ret[:kdistances] = data["distances"]
        ret[:ticks] = data["ticks"]
        ret[key] = [Dict(spinsym => data["energy"][ibranch][spin]
                         for (spin, spinsym) in spinmap)
                    for ibranch = 1:length(data["energy"])]
    end

    (; ret...)  # Make it a named tuple and return
end

"""
    is_metal(band_data, εF, tol)

Determine whether the provided bands indicate the material is a metal,
i.e. where bands are cut by the Fermi level.
"""
function is_metal(band_data, εF, tol=1e-4)
    # This assumes no spin polarization
    @assert band_data.basis.model.spin_polarization in (:none, :spinless)

    n_bands = length(band_data.λ[1])
    n_kpoints = length(band_data.λ)
    for ib in 1:n_bands
        some_larger = any(band_data.λ[ik][ib] - εF < -tol for ik in 1:n_kpoints)
        some_smaller = any(εF - band_data.λ[ik][ib] < -tol for ik in 1:n_kpoints)
        some_larger && some_smaller && return true
    end
    false
end

function plot_band_data(band_data; εF=nothing,
                        klabels=Dict{String, Vector{Float64}}(), unit=:eV, kwargs...)
    eshift = isnothing(εF) ? 0.0 : εF
    data = prepare_band_data(band_data, klabels=klabels)

    # For each branch, plot all bands, spins and errors
    p = Plots.plot(xlabel="wave vector")
    for ibranch = 1:data.n_branches
        kdistances = data.kdistances[ibranch]
        for spin in data.spins, iband = 1:data.n_bands
            yerror = nothing
            if hasproperty(data, :λerror)
                yerror = data.λerror[ibranch][spin][iband, :] ./ unit_to_au(unit)
            end
            energies = (data.λ[ibranch][spin][iband, :] .- eshift) ./ unit_to_au(unit)

            color = (spin == :up) ? :blue : :red
            Plots.plot!(p, kdistances, energies; color=color, label="", yerror=yerror,
                        kwargs...)
        end
    end

    # X-range: 0 to last kdistance value
    Plots.xlims!(p, (0, data.kdistances[end][end]))
    Plots.xticks!(p, data.ticks["distance"],
                  [replace(l, raw"$\mid$" => " | ") for l in data.ticks["label"]])

    ylims = [-4, 4]
    is_metal(band_data, εF) && (ylims = [-10, 10])
    ylims = round.(ylims * units.eV ./ unit_to_au(unit), sigdigits=2)
    if isnothing(εF)
        Plots.ylabel!(p, "eigenvalues  ($(string(unit))")
    else
        Plots.ylabel!(p, "eigenvalues - ε_f  ($(string(unit)))")
        Plots.ylims!(p, ylims...)
    end

    p
end

function detexify_kpoint(string)
    # For some reason Julia doesn't support this naively: https://github.com/JuliaLang/julia/issues/29849
    replacements = ("\\Gamma" => "Γ",
                    "\\Delta" => "Δ",
                    "\\Sigma" => "Σ")
    for r in replacements
        string = replace(string, r)
    end
    string
end

# TODO This is the top-level function, which should be documented
function plot_bandstructure(basis, ρ, n_bands;
                            εF=nothing, kline_density=20, unit=:eV, kwargs...)
    # Band structure calculation along high-symmetry path
    kcoords, klabels, kpath = high_symmetry_kpath(basis.model; kline_density=kline_density)
    println("Computing bands along kpath:")
    println("       ", join(join.(detexify_kpoint.(kpath), " -> "), "  and  "))
    band_data = compute_bands(basis, ρ, kcoords, n_bands; kwargs...)

    plotargs = ()
    if kline_density ≤ 10
        plotargs = (markersize=2, markershape=:circle)
    end
    plot_band_data(band_data; εF=εF, klabels=klabels, unit=unit, plotargs...)
end
function plot_bandstructure(scfres, n_bands; kwargs...)
    # Convenience wrapper for scfres named tuples
    plot_bandstructure(scfres.ham.basis, scfres.ρ, n_bands; εF=scfres.εF, kwargs...)
end

ext: jl


text:
using Statistics
using DelimitedFiles
using Random
using Printf

struct Lattice
    nsites::Int64
    nup::Int64
    ndn::Int64
    ntot::Int64
    t::Float64
    U::Float64
end

struct QMC
    nblocks::Int64
    samples::Int64
    Δτ::Float64
    L::Int64
    γ::Float64
    Kmatrix::Array{Float64,2}
end

struct Opt  # Optimization
    stab_interval::Int64
    nburn_in::Int64
end

const FloatType = Union{Float64, Complex{Float64}}

include("./MatrixGenerator.jl")
include("./Calculations.jl")

function getinputs()
    lines = readlines("./InputFermiHubbard")

    # basic parameters
    β = parse(Float64, split(lines[1])[2])
    num_sites = parse(Int64, split(lines[2])[2])
    num_spinup = parse(Int64, split(lines[3])[2])
    num_spindn = parse(Int64, split(lines[4])[2])
    num_elec = num_spinup + num_spinup

    Δτ = parse(Float64, split(lines[5])[2])
    t = parse(Float64, split(lines[6])[2])
    U = parse(Float64, split(lines[7])[2])
    num_blocks = convert(Int64, parse(Float64, split(lines[8])[2]))
    num_samples = convert(Int64, parse(Float64, split(lines[9])[2]))

    # optimization parameters
    stab_interval = parse(Int64, split(lines[12])[2])
    num_burnins = convert(Int64, parse(Float64, split(lines[13])[2]))

    # constants for calculations
    L = convert(Int64, round(β / Δτ))
    γ = atanh(sqrt(tanh(Δτ * U / 4)))

    lattice = Lattice(num_sites, num_spinup, num_spindn, num_elec, t, U)

    K_matrix = kinetic_matrix(lattice)
    qmc = QMC(num_blocks, num_samples, Δτ, L, γ, exp(-K_matrix * Δτ / 2))

    opt = Opt(stab_interval, num_burnins)
    return lattice, qmc, opt
end

function mainbody()
    lattice, qmc, opt = getinputs()
    SubEk_array = Vector{Float64}()
    SubEp_array = Vector{Float64}()
    SubE_array = Vector{Float64}()
    SubEk_error_array = Vector{Float64}()
    SubEp_error_array = Vector{Float64}()
    SubE_error_array = Vector{Float64}()
    filename = string("output_Ns", lattice.nsites, "_U", lattice.U, "_beta", qmc.Δτ * qmc.L)
    outputfile = open(filename, "a+")
    @printf(outputfile, "       Canonical Ensemble AFQMC for Hubbard Model     \n\n")
    @printf(outputfile, "Number of Blocks: %d\n", qmc.nblocks)
    @printf(outputfile, "Length of Img Timestep: %.4f\n", qmc.Δτ)
    @printf(outputfile, "Stablization Interval: %d\n\n", opt.stab_interval)
    for i = 1 : qmc.nblocks
        SubE, SubE_error, SubEk, SubEk_error, SubEp, SubEp_error, sgn_mean, NiMax = MCIntegration(lattice, qmc, opt)
        SubE = SubE / sgn_mean
        SubEk = SubEk / sgn_mean
        SubEp = SubEp / sgn_mean

        @printf(outputfile, "Block # %d out of %d\n", i, qmc.nblocks)
        @printf(outputfile, "--------------------\n")
        @printf(outputfile, "Average Sign: %.2f\n", sgn_mean)
        @printf(outputfile, "Maximum Occupancy: %.2f\n", NiMax)
        @printf(outputfile, "Ek            Ep            Etot\n")
        @printf(outputfile, "%.5f     %.5f       %.5f\n", SubEk, SubEp, SubE)
        @printf(outputfile, "Error:\n")
        @printf(outputfile, "%.5f     %.5f       %.5f\n\n", SubEk_error, SubEp_error, SubE_error)

        push!(SubEk_array, SubEk)
        push!(SubEk_error_array, SubEk_error)
        push!(SubEp_array, SubEp)
        push!(SubEp_error_array, SubEp_error)
        push!(SubE_array, SubE)
        push!(SubE_error_array, SubE_error)
    end
    Ek_mean, Ek_error = mean(SubEk_array), mean(SubEk_error_array)
    Ep_mean, Ep_error = mean(SubEp_array), mean(SubEp_error_array)
    E_mean, E_error = mean(SubE_array), mean(SubE_error_array)

    @printf(outputfile, "--------------------\n")
    @printf(outputfile, "Results for this run:\n")
    @printf(outputfile, "Ek            Ep            Etot\n")
    @printf(outputfile, "%.5f     %.5f       %.5f\n", Ek_mean, Ep_mean, E_mean)
    @printf(outputfile, "Error:\n")
    @printf(outputfile, "%.5f     %.5f       %.5f\n\n", Ek_error, Ep_error, E_error)
    close(outputfile)
end

mainbody()

max_stars_repo_name: TongSericus/CanonicalAFQMC

max_issues_repo_name: TongSericus/CanonicalAFQMC

max_forks_repo_name: TongSericus/CanonicalAFQMC


text:
using SailRoute, Test

@testset "Test conversion to polar coordinates" begin
    r_1, theta_1 = SailRoute.calc_polars(1, 1)
    @test r_1 ≈ sqrt(2)
    @test theta_1 ≈ 45.0
    r_2, theta_2 = SailRoute.calc_polars(1, -1)
    @test theta_2 ≈ 135
    r_3, theta_3 = SailRoute.calc_polars(-1, -1)
    @test theta_3 ≈ 225
    r_4, theta_4 = SailRoute.calc_polars(-1, 1)
    @test theta_4 ≈ 315
end

@testset "Generate constant weather sets" begin
    tws_val = 1.0; twa_val = 90.0; cs_val=0.0; cd_val = 90.0;
    wahi_val = 0.0; wadi_val=90.0; n_locs=16;
    t_len = 100
    tws, twa, cs, cd, wahi, wadi = SailRoute.generate_constant_weather(tws_val, twa_val, cs_val, cd_val,
    wahi_val, wadi_val, n_locs, t_len)
    @test tws[1, 1, 1] ≈ tws_val
    @test twa[t_len, n_locs, 2] ≈ twa_val
    @test cs[5, 3, n_locs] ≈ cs_val
    @test cd[1, 1, 1] ≈ cd_val
    @test wahi[1, 1, 1] ≈ wahi_val
    @test wadi[1, 1, 1] ≈ wadi_val
end


text:
"""
```
function ∂evaluate(::Val{track}, f::F, args...) where {track, F}
end
```

To add direct support to ProbabilityModels for differentiating a function, you need to extend ∂evaluate to add an appropriate method of the signature:

`∂evaluate(::Val{track}, f::typeof(extendedfunction), args...)`

returning the tuple:

`f(args...), ∇`

where `∇` is a callable object (eg, an anonymous function or a callable user defined struct) such that
∇(x̄)
returns a gradient (`length(args)==1`) or a tuple of the gradients with respect to the args, where `x̄` contains the gradient with respect to `f(args...)`.

`track` is a N-tuple of booleans, indicating whether the `n`-th argument is being tracked. The gradients where `track==false` are ignored; it may be possible
to use that information to save on computation time, although if `∇` are lazy or everything gets inlined so the compiler can DCE unused results, it may be able
to avoid that computation regardless.
"""
function ∂evaluate end


text:
function main()

	L = [37107287533902102798797998220837590246510135740250
46376937677490009712648124896970078050417018260538
74324986199524741059474233309513058123726617309629
91942213363574161572522430563301811072406154908250
23067588207539346171171980310421047513778063246676
89261670696623633820136378418383684178734361726757
28112879812849979408065481931592621691275889832738
44274228917432520321923589422876796487670272189318
47451445736001306439091167216856844588711603153276
70386486105843025439939619828917593665686757934951
62176457141856560629502157223196586755079324193331
64906352462741904929101432445813822663347944758178
92575867718337217661963751590579239728245598838407
58203565325359399008402633568948830189458628227828
80181199384826282014278194139940567587151170094390
35398664372827112653829987240784473053190104293586
86515506006295864861532075273371959191420517255829
71693888707715466499115593487603532921714970056938
54370070576826684624621495650076471787294438377604
53282654108756828443191190634694037855217779295145
36123272525000296071075082563815656710885258350721
45876576172410976447339110607218265236877223636045
17423706905851860660448207621209813287860733969412
81142660418086830619328460811191061556940512689692
51934325451728388641918047049293215058642563049483
62467221648435076201727918039944693004732956340691
15732444386908125794514089057706229429197107928209
55037687525678773091862540744969844508330393682126
18336384825330154686196124348767681297534375946515
80386287592878490201521685554828717201219257766954
78182833757993103614740356856449095527097864797581
16726320100436897842553539920931837441497806860984
48403098129077791799088218795327364475675590848030
87086987551392711854517078544161852424320693150332
59959406895756536782107074926966537676326235447210
69793950679652694742597709739166693763042633987085
41052684708299085211399427365734116182760315001271
65378607361501080857009149939512557028198746004375
35829035317434717326932123578154982629742552737307
94953759765105305946966067683156574377167401875275
88902802571733229619176668713819931811048770190271
25267680276078003013678680992525463401061632866526
36270218540497705585629946580636237993140746255962
24074486908231174977792365466257246923322810917141
91430288197103288597806669760892938638285025333403
34413065578016127815921815005561868836468420090470
23053081172816430487623791969842487255036638784583
11487696932154902810424020138335124462181441773470
63783299490636259666498587618221225225512486764533
67720186971698544312419572409913959008952310058822
95548255300263520781532296796249481641953868218774
76085327132285723110424803456124867697064507995236
37774242535411291684276865538926205024910326572967
23701913275725675285653248258265463092207058596522
29798860272258331913126375147341994889534765745501
18495701454879288984856827726077713721403798879715
38298203783031473527721580348144513491373226651381
34829543829199918180278916522431027392251122869539
40957953066405232632538044100059654939159879593635
29746152185502371307642255121183693803580388584903
41698116222072977186158236678424689157993532961922
62467957194401269043877107275048102390895523597457
23189706772547915061505504953922979530901129967519
86188088225875314529584099251203829009407770775672
11306739708304724483816533873502340845647058077308
82959174767140363198008187129011875491310547126581
97623331044818386269515456334926366572897563400500
42846280183517070527831839425882145521227251250327
55121603546981200581762165212827652751691296897789
32238195734329339946437501907836945765883352399886
75506164965184775180738168837861091527357929701337
62177842752192623401942399639168044983993173312731
32924185707147349566916674687634660915035914677504
99518671430235219628894890102423325116913619626622
73267460800591547471830798392868535206946944540724
76841822524674417161514036427982273348055556214818
97142617910342598647204516893989422179826088076852
87783646182799346313767754307809363333018982642090
10848802521674670883215120185883543223812876952786
71329612474782464538636993009049310363619763878039
62184073572399794223406235393808339651327408011116
66627891981488087797941876876144230030984490851411
60661826293682836764744779239180335110989069790714
85786944089552990653640447425576083659976645795096
66024396409905389607120198219976047599490197230297
64913982680032973156037120041377903785566085089252
16730939319872750275468906903707539413042652315011
94809377245048795150954100921645863754710598436791
78639167021187492431995700641917969777599028300699
15368713711936614952811305876380278410754449733078
40789923115535562561142322423255033685442488917353
44889911501440648020369068063960672322193204149535
41503128880339536053299340368006977710650566631954
81234880673210146739058568557934581403627822703280
82616570773948327592232845941706525094512325230608
22918802058777319719839450180888072429661980811197
77158542502016545090413245809786882778948721859617
72107838435069186155435662884062257473692284509516
20849603980134001723930671666823555245252804609722
53503534226472524250874054075591789781264330331690]

	s = string(sum(L))

	assert(length(s) >= 10)
	n = s[1:10]
	println(n)
end

main()

max_stars_repo_name: cmey/euler

max_issues_repo_name: cmey/euler

max_forks_repo_name: cmey/euler


text:
# AUTO GENERATED FILE - DO NOT EDIT

export accordion

"""
    accordion(;kwargs...)
    accordion(children::Any;kwargs...)
    accordion(children_maker::Function;kwargs...)

An Accordion component.
Divide content into collapsible sections. For more information, see: https://mantine.dev/core/accordion/
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): <AccordionItem /> components only
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `className` (String; optional): Often used with CSS to style elements with common properties
- `disableIconRotation` (Bool; optional): Should icon rotation be disabled
- `iconPosition` (a value equal to: "right", "left"; optional): Change icon position: left or right
- `multiple` (Bool; optional): Allow multiple items to be opened at the same time
- `offsetIcon` (Bool; optional): Should icon be offset with padding, applicable only when iconPosition is left
- `state` (Dict with Strings as keys and values of type Bool; optional): Controlled state (controls opened state of accordion items)
"""
function accordion(; kwargs...)
        available_props = Symbol[:children, :id, :className, :disableIconRotation, :iconPosition, :multiple, :offsetIcon, :state]
        wild_props = Symbol[]
        return Component("accordion", "Accordion", "dash_mantine_components", available_props, wild_props; kwargs...)
end

accordion(children::Any; kwargs...) = accordion(;kwargs..., children = children)
accordion(children_maker::Function; kwargs...) = accordion(children_maker(); kwargs...)


text:
# File: AUTD3.jl
# Project: src
# Created Date: 11/02/2020
# Author: Shun Suzuki
# -----
# Last Modified: 17/12/2021
# Modified By: Shun Suzuki (suzuki@hapis.k.u-tokyo.ac.jp)
# -----
# Copyright (c) 2020 Hapis Lab. All rights reserved.
#

module AUTD3

using StaticArrays

export Gain, Modulation, PointSequence, GainSequence, Link, AUTD, GainMode, duty_phase_full, phase_full, phase_half
export open_autd, dispose, add_device, synchronize, clear, set_silent_mode, stop, pause, resume, is_open, get_silent_mode, wavelength, set_wavelength
export num_devices
export device_direction_x, device_direction_y, device_direction_z, trans_position
export enumerate_adapters, firmware_info_list
export last_error
export focal_point_gain, focal_point_gain_with_duty, grouped_gain, bessel_beam_gain, bessel_beam_gain_with_duty, plane_wave_gain, plane_wave_gain_with_duty, custom_gain, holo_gain, transducer_test_gain, null_gain
export holo_gain_sdp, holo_gain_evd, holo_gain_gs, holo_gain_gspat, holo_gain_lm, holo_gain_naive
export static_modulation, sine_modulation, sine_pressure_modulation, sine_legacy_modulation, custom_modulation, saw_modulation, square_modulation
export sequence, gain_sequence, add_point, add_gain, set_freq, get_freq, get_sampling_freq, get_sampling_freq_div
export soem_link, twincat_link
export send
export stm, add_stm_gain, start_stm, stop_stm, finish_stm
export NUM_TRANS_IN_UNIT, NUM_TRANS_X, NUM_TRANS_Y, TRANS_SPACING_MM, DEVICE_WIDTH, DEVICE_HEIGHT

include("NativeMethods.jl")

const NUM_TRANS_IN_UNIT = 249
const NUM_TRANS_X = 18
const NUM_TRANS_Y = 14
const TRANS_SPACING_MM = 10.16
const DEVICE_WIDTH = 192.0
const DEVICE_HEIGHT = 151.4

mutable struct Gain
    _ptr::Ptr{Cvoid}
    _disposed::Bool
    function Gain(gain_ptr::Ptr{Cvoid})
        gain = new(gain_ptr, false)
        finalizer(gain -> dispose(gain), gain)
        gain
    end
end

mutable struct Modulation
    _ptr::Ptr{Cvoid}
    _disposed::Bool
    function Modulation(mod_ptr::Ptr{Cvoid})
        modulation = new(mod_ptr, false)
        finalizer(modulation -> dispose(modulation), modulation)
        modulation
    end
end

mutable struct PointSequence
    _ptr::Ptr{Cvoid}
    _disposed::Bool
    function PointSequence(seq_ptr::Ptr{Cvoid})
        seq = new(seq_ptr, false)
        finalizer(seq -> dispose(seq), seq)
        seq
    end
end

mutable struct GainSequence
    _ptr::Ptr{Cvoid}
    _disposed::Bool
    function GainSequence(seq_ptr::Ptr{Cvoid})
        seq = new(seq_ptr, false)
        finalizer(seq -> dispose(seq), seq)
        seq
    end
end

mutable struct Link
    _link_ptr::Ptr{Cvoid}
    function Link(link_ptr::Ptr{Cvoid})
        new(link_ptr)
    end
end

mutable struct AUTD
    _handle::Ptr{Cvoid}
    _disposed::Bool
    function AUTD()
        chandle = Ref(Ptr{Cvoid}(0))
        autd_create_controller(chandle)
        autd = new(chandle[], false)
        finalizer(autd -> dispose(autd), autd)
        autd
    end
end

mutable struct STMController
    _handle::Ptr{Cvoid}
    function STMController(handle::Ptr{Cvoid})
        new(handle)
    end
end

function open_autd(autd::AUTD, link::Link)
    autd_open_controller(autd._handle, link._link_ptr)
end

function dispose(self::AUTD)
    if !self._disposed
        autd_close_controller(self._handle)
        autd_free_controller(self._handle)
        self._disposed = true
    end
end

function dispose(self::Gain)
    if !self._disposed
        autd_delete_gain(self._ptr)
        self._disposed = true
    end
end

function dispose(self::Modulation)
    if !self._disposed
        autd_delete_modulation(self._ptr)
        self._disposed = true
    end
end

function dispose(self::PointSequence)
    if !self._disposed
        autd_delete_sequence(self._ptr)
        self._disposed = true
    end
end

function dispose(self::GainSequence)
    if !self._disposed
        autd_delete_sequence(self._ptr)
        self._disposed = true
    end
end

function enumerate_adapters()
    res = []
    phandle = Ref(Ptr{Cvoid}(0))
    size = autd_get_adapter_pointer(phandle)
    handle::Ptr{Cvoid} = phandle[]

    for i = 0:size-1
        sb_desc = zeros(UInt8, 128)
        sb_name = zeros(UInt8, 128)
        autd_get_adapter(handle, i, sb_desc, sb_name)
        push!(res, [String(strip(String(sb_desc), '\0')), String(strip(String(sb_name), '\0'))])
    end

    autd_free_adapter_pointer(handle)
    res
end

function firmware_info_list(autd::AUTD)
    res = []
    phandle = Ref(Ptr{Cvoid}(0))
    size = autd_get_firmware_info_list_pointer(autd._handle, phandle)
    handle::Ptr{Cvoid} = phandle[]

    for i = 0:size-1
        sb_cpu = zeros(UInt8, 128)
        sb_fpga = zeros(UInt8, 128)
        autd_get_firmware_info(handle, i, sb_cpu, sb_fpga)
        push!(res, [String(strip(String(sb_cpu), '\0')), String(strip(String(sb_fpga), '\0'))])
    end

    autd_free_firmware_info_list_pointer(handle)
    res
end

function add_device(autd::AUTD, pos::SVector{3,Float64}, rot::SVector{3,Float64})
    x, y, z = pos
    az1, ay, az2 = rot
    autd_add_device(autd._handle, x, y, z, az1, ay, az2)
end

function add_device(autd::AUTD, pos::SVector{3,Float64}, qua::SVector{4,Float64})
    x, y, z = pos
    qw, qx, qy, qz = qua
    autd_add_device_quaternion(autd._handle, x, y, z, qw, qx, qy, qz)
end

function clear(autd::AUTD)
    autd_clear(autd._handle)
end

function update_ctrl_flags(autd::AUTD)
    autd_update_ctrl_flags(autd._handle)
end

function stop(autd::AUTD)
    autd_stop(autd._handle)
end

function pause(autd::AUTD)
    autd_pause(autd._handle)
end

function resume(autd::AUTD)
    autd_resume(autd._handle)
end

function is_open(autd::AUTD)
    autd_is_open(autd._handle)
end

function last_error()
    size = autd_get_last_error(Ptr{Cvoid}(0))
    err = zeros(UInt8, size)
    autd_get_last_error(err)
    String(strip(String(err), '\0'))
end

function set_silent_mode(autd::AUTD, silent::Bool)
    autd_set_silent_mode(autd._handle, silent)
end

function get_silent_mode(autd::AUTD)
    autd_get_silent_mode(autd._handle)
end

function set_force_fan(autd::AUTD, force_fan::Bool)
    autd_set_force_fan(autd._handle, force_fan)
end

function get_force_fan(autd::AUTD)
    autd_get_force_fan(autd._handle)
end

function wavelength(autd::AUTD)
    autd_get_wavelength(autd._handle)
end

function set_wavelength(autd::AUTD, wavelength::Float64)
    autd_set_wavelength(autd._handle, wavelength)
end

function attenuation(autd::AUTD)
    autd_get_attenuation(autd._handle)
end

function set_attenuation(autd::AUTD, attenuation::Float64)
    autd_set_attenuation(autd._handle, attenuation)
end

function num_devices(autd::AUTD)
    autd_num_devices(autd._handle)
end

function adjust_amp(amp::Float64)
    d = asin(amp) / pi
    UInt8(510.0 * d)
end

function focal_point_gain_with_duty(position::SVector{3,Float64}; duty::UInt8)
    x, y, z = position
    chandle = Ref(Ptr{Cvoid}(0))
    autd_gain_focal_point(chandle, x, y, z, duty)
    Gain(chandle[])
end

function focal_point_gain(position::SVector{3,Float64}; amp::Float64 = 1.0)
    focal_point_gain_with_duty(position; duty = adjust_amp(amp))
end

function grouped_gain(autd, group_ids::Array{Tuple{Int32,Gain},1})
    chandle = Ref(Ptr{Cvoid}(0))
    autd_gain_grouped(chandle, autd._ptr)
    for (id, gain) in group_ids
        autd_gain_grouped_add(chandle[], id, gain)
    end
    Gain(chandle[])
end

function bessel_beam_gain_with_duty(position::SVector{3,Float64}, direction::SVector{3,Float64}, theta_z::Float64; duty::UInt8)
    x, y, z = position
    nx, ny, nz = direction
    chandle = Ref(Ptr{Cvoid}(0))
    autd_gain_bessel_beam(chandle, x, y, z, nx, ny, nz, theta_z, duty)
    Gain(chandle[])
end

function bessel_beam_gain(position::SVector{3,Float64}, direction::SVector{3,Float64}, theta_z::Float64; amp::Float64 = 1.0)
    bessel_beam_gain_with_duty(position, direction, theta_z; duty = adjust_amp(amp))
end

function plane_wave_gain_with_duty(direction::SVector{3,Float64}; duty::UInt8)
    nx, ny, nz = direction
    chandle = Ref(Ptr{Cvoid}(0))
    autd_gain_plane_wave(chandle, nx, ny, nz, duty)
    Gain(chandle[])
end

function plane_wave_gain(direction::SVector{3,Float64}; amp::Float64 = 1.0)
    plane_wave_gain_with_duty(direction; duty = adjust_amp(amp))
end

function custom_gain(data::Array{UInt16,1})
    len = length(data)
    chandle = Ref(Ptr{Cvoid}(0))
    autd_gain_custom(chandle, data, len)
    Gain(chandle[])
end

function _pack_foci(foci::Array{SVector{3,Float64},1})
    len = length(foci)
    foci_array = zeros(Float64, len * 3)
    for (i, focus) in enumerate(foci)
        foci_array[3*(i-1)+1] = focus[1]
        foci_array[3*(i-1)+2] = focus[2]
        foci_array[3*(i-1)+3] = focus[3]
    end
    foci_array
end

function holo_gain_sdp(foci::Array{SVector{3,Float64},1}, amps::Array{Float64,1}; alpha::Float64 = 1e-3, lambda::Float64 = 0.9, repeat::UInt64 = UInt64(100), normalize::Bool = true)
    len = length(foci)

    foci_array = _pack_foci(foci)
    backend = Ref(Ptr{Cvoid}(0))
    autd_eigen_backend(backend)

    chandle = Ref(Ptr{Cvoid}(0))
    autd_gain_holo_sdp(chandle, backend[], foci_array, amps, len, alpha, lambda, repeat, normalize)
    autd_delete_backend(backend[])
    Gain(chandle[])
end

function holo_gain_evd(foci::Array{SVector{3,Float64},1}, amps::Array{Float64,1}; gamma::Float64 = 1.0, normalize::Bool = true)
    len = length(foci)

    foci_array = _pack_foci(foci)
    backend = Ref(Ptr{Cvoid}(0))
    autd_eigen_backend(backend)

    chandle = Ref(Ptr{Cvoid}(0))
    autd_gain_holo_evd(chandle, backend[], foci_array, amps, len, gamma, normalize)
    autd_delete_backend(backend[])
    Gain(chandle[])
end

function holo_gain_naive(foci::Array{SVector{3,Float64},1}, amps::Array{Float64,1})
    len = length(foci)

    foci_array = _pack_foci(foci)
    backend = Ref(Ptr{Cvoid}(0))
    autd_eigen_backend(backend)

    chandle = Ref(Ptr{Cvoid}(0))
    autd_gain_holo_naive(chandle, backend[], foci_array, amps, len)
    autd_delete_backend(backend[])
    Gain(chandle[])
end

function holo_gain_gs(foci::Array{SVector{3,Float64},1}, amps::Array{Float64,1}; repeat::UInt64 = UInt64(100))
    len = length(foci)

    foci_array = _pack_foci(foci)
    backend = Ref(Ptr{Cvoid}(0))
    autd_eigen_backend(backend)

    chandle = Ref(Ptr{Cvoid}(0))
    autd_gain_holo_gs(chandle, backend[], foci_array, amps, len, repeat)
    autd_delete_backend(backend[])
    Gain(chandle[])
end

function holo_gain_gspat(foci::Array{SVector{3,Float64},1}, amps::Array{Float64,1}; repeat::UInt64 = UInt64(100))
    len = length(foci)

    foci_array = _pack_foci(foci)
    backend = Ref(Ptr{Cvoid}(0))
    autd_eigen_backend(backend)

    chandle = Ref(Ptr{Cvoid}(0))
    autd_gain_holo_gspat(chandle, backend[], foci_array, amps, len, repeat)
    autd_delete_backend(backend[])
    Gain(chandle[])
end

function holo_gain_lm(foci::Array{SVector{3,Float64},1}, amps::Array{Float64,1}; eps_1::Float64 = 1e-8, eps_2::Float64 = 1e-8, tau::Float64 = 1e-3, k_max::UInt64 = UInt64(5), initial::Array{Float64,1} = [])
    len = length(foci)

    foci_array = _pack_foci(foci)
    backend = Ref(Ptr{Cvoid}(0))
    autd_eigen_backend(backend)

    chandle = Ref(Ptr{Cvoid}(0))
    initial_len = length(initial)
    init = initial_len == 0 ? Ptr{Cvoid}(0) : initial
    autd_gain_holo_lm(chandle, backend[], foci_array, amps, len, eps_1, eps_2, tau, k_max, init, initial_len)
    autd_delete_backend(backend[])
    Gain(chandle[])
end

function holo_gain_greedy(foci::Array{SVector{3,Float64},1}, amps::Array{Float64,1}, phase_div::Int32)
    len = length(foci)

    foci_array = _pack_foci(foci)
    backend = Ref(Ptr{Cvoid}(0))
    autd_eigen_backend(backend)

    chandle = Ref(Ptr{Cvoid}(0))
    autd_gain_holo_greedy(chandle, backend, foci_array, amps, len, phase_div)
    autd_delete_backend(backend[])
    Gain(chandle[])
end

function transducer_test_gain(trans_idx::Int32, duty::UInt8, phase::UInt8)
    chandle = Ref(Ptr{Cvoid}(0))
    autd_gain_transducer_test(chandle, trans_idx, duty, phase)
    Gain(chandle[])
end

function null_gain()
    chandle = Ref(Ptr{Cvoid}(0))
    autd_gain_null(chandle)
    Gain(chandle[])
end

function static_modulation(amp::UInt8 = 0xFF)
    chandle = Ref(Ptr{Cvoid}(0))
    autd_modulation_static(chandle, amp)
    Modulation(chandle[])
end

function custom_modulation(data::Array{UInt8,1})
    len = length(data)
    chandle = Ref(Ptr{Cvoid}(0))
    autd_modulation_custom(chandle, data, len)
    Modulation(chandle[])
end

function sine_modulation(freq::Int32, amp::Float64 = 1.0, offset::Float64 = 0.5)
    chandle = Ref(Ptr{Cvoid}(0))
    autd_modulation_sine(chandle, freq, amp, offset)
    Modulation(chandle[])
end

function sine_squared_modulation(freq::Int32, amp::Float64 = 1.0, offset::Float64 = 0.5)
    chandle = Ref(Ptr{Cvoid}(0))
    autd_modulation_sine_squared(chandle, freq, amp, offset)
    Modulation(chandle[])
end

function sine_legacy_modulation(freq::Float64, amp::Float64 = 1.0, offset::Float64 = 0.5)
    chandle = Ref(Ptr{Cvoid}(0))
    autd_modulation_sine_legacy(chandle, freq, amp, offset)
    Modulation(chandle[])
end

function square_modulation(freq::Int32, low::UInt8 = 0x00, high::UInt8 = 0xFF, duty::Float64 = 0.5)
    chandle = Ref(Ptr{Cvoid}(0))
    autd_modulation_square(chandle, freq, low, high, duty)
    Modulation(chandle[])
end

function sequence()
    chandle = Ref(Ptr{Cvoid}(0))
    autd_sequence(chandle)
    PointSequence(chandle[])
end

function gain_sequence(autd, gain_mode::GainMode)
    chandle = Ref(Ptr{Cvoid}(0))
    autd_gain_sequence(chandle, autd._handle, UInt16(gain_mode))
    GainSequence(chandle[])
end

function add_point(seq::PointSequence, point::SVector{3,Float64}, duty::UInt8 = 0xFF)
    x, y, z = point
    autd_sequence_add_point(seq._ptr, x, y, z, duty)
end

function add_gain(seq::GainSequence, gain::Gain)
    autd_sequence_add_gain(seq._ptr, gain._ptr)
end

function set_freq(seq, freq::Float64)
    autd_sequence_set_freq(seq._ptr, freq)
end

function get_freq(seq)
    autd_sequence_freq(seq._ptr)
end

function get_sampling_freq(seq)
    autd_sequence_sampling_freq(seq._ptr)
end

function get_sampling_freq_div(seq)
    autd_sequence_sampling_freq_div(seq._ptr)
end

function soem_link(ifname::String, device_num::Int32, cycle_ticks::UInt32 = UInt32(1))
    chandle = Ref(Ptr{Cvoid}(0))
    autd_link_soem(chandle, ifname, device_num, cycle_ticks)
    Link(chandle[])
end

function twincat_link()
    chandle = Ref(Ptr{Cvoid}(0))
    autd_link_twincat(chandle)
    Link(chandle[])
end

function emulator_link(port::UInt16, autd::AUTD)
    chandle = Ref(Ptr{Cvoid}(0))
    autd_link_emulator(chandle, port, autd._handle)
    Link(chandle[])
end

function send(autd::AUTD, header, body)
    if header !== Nothing && body !== Nothing
        autd_send_header_body(autd._handle, header._ptr, body._ptr)
    elseif header !== Nothing
        autd_send_header(autd._handle, header._ptr)
    elseif body !== Nothing
        autd_send_body(autd._handle, body._ptr)
    end
end

function stm(autd::AUTD)
    chandle = Ref(Ptr{Cvoid}(0))
    autd_stm_controller(chandle, autd._handle)
    STMController(chandle[])
end

function add_stm_gain(stm::STMController, gain::Gain)
    autd_add_stm_gain(stm._handle, gain._ptr)
end

function start_stm(stm::STMController, freq::Float64)
    autd_start_stm(stm._handle, freq)
end

function stop_stm(stm::STMController)
    autd_stop_stm(stm._handle)
end

function finish_stm(stm::STMController)
    autd_finish_stm(stm._handle)
end

function trans_position(autd::AUTD, device_idx::Int32, local_trans_idx::Int32)
    x = Ref{Float64}(0)
    y = Ref{Float64}(0)
    z = Ref{Float64}(0)
    autd_trans_position(autd._handle, device_idx, local_trans_idx, x, y, z)
    SVector(x[], y[], z[])
end

function device_direction_x(autd::AUTD, device_idx::Int32)
    x = Ref{Float64}(0)
    y = Ref{Float64}(0)
    z = Ref{Float64}(0)
    autd_device_x_direction(autd._handle, device_idx, x, y, z)
    SVector(x[], y[], z[])
end

function device_direction_y(autd::AUTD, device_idx::Int32)
    x = Ref{Float64}(0)
    y = Ref{Float64}(0)
    z = Ref{Float64}(0)
    autd_device_y_direction(autd._handle, device_idx, x, y, z)
    SVector(x[], y[], z[])
end

function device_direction_z(autd::AUTD, device_idx::Int32)
    x = Ref{Float64}(0)
    y = Ref{Float64}(0)
    z = Ref{Float64}(0)
    autd_device_z_direction(autd._handle, device_idx, x, y, z)
    SVector(x[], y[], z[])
end

end


text:
@testset "395.longest-substring-with-at-least-k-repeating-characters.jl" begin
    @test longest_substring_395("ababbc", 2) == 5
    @test longest_substring_395("aaabb", 3) == 3
end


text:
export addvar, within, allvars, REF, UNKNOWN, PATH_COMPONENT, isemptyid, json

const REF = r"^\?([0-9]+)$"
const LOCAL_ID = r"^@/([0-9]+)$"
const UNKNOWN = "?"
const PATH_COMPONENT = r"^([-[:alnum:]]+|@)/([0-9]+)$"
const VAR_NAME = r"^([0-9]+|\pL\p{Xan}*)(?::((\pL\p{Xan}*)(?:=((?:[^,]|\\,)*))?(?:,(\pL\p{Xan}*)(?:=((?:[^,]|\\,)*))?)*))?$"
const METAPROP = r"(\pL\p{Xan}*)(?:=((?:[^,]|\\,)*))?(,|$)"
const ARRAY_INDEX = r"[1-9][0-9]*"
const JPATH_COMPONENT = r"^([\p{Lu}\p{Ll}\p{Lt}\p{Lm}\p{Lo}\p{Nl}\p{Sc}\p{So}_][\p{Lu}\p{Ll}\p{Lt}\p{Lm}\p{Lo}\p{Nl}\p{Sc}\p{So}_!\p{Nd}\p{No}\p{Mn}\p{Mc}\p{Me}\p{Sk}\p{Pc}]*|[1-9][0-9]*)(?:((?:\.[\p{Lu}\p{Ll}\p{Lt}\p{Lm}\p{Lo}\p{Nl}\p{Sc}\p{So}_][\p{Lu}\p{Ll}\p{Lt}\p{Lm}\p{Lo}\p{Nl}\p{Sc}\p{So}_!\p{Nd}\p{No}\p{Mn}\p{Mc}\p{Me}\p{Sk}\p{Pc}]*)*)(?:\(\))?)?$"

oid(cmd::JusCmd, obj) = oid(connection(cmd), obj)
function oid(con::Connection, obj)
    get!(con.data2oid, obj) do
        oid = con.nextOid += 1
        con.oid2data[oid] = WeakRef(obj)
        oid
    end
end

json(cmd::VarCommand, id::ID) = json(cmd.connection, value)
json(cmd::JusCmd, id::ID) = "$(id.namespace === cmd.namespace ? "@" : id.namespace)/$(id.number)"
json(con::Connection, id::ID) = "$(id.namespace === con.namespace ? "@" : id.namespace)/$(id.number)"
json(id::ID) = "$(id.namespace)/$(id.number)"
json(cmd::JusCmd, data) = json(cmd.config, connection(cmd), data)
json(cmd::VarCommand, value) = json(cmd.config, cmd.connection, value)
function json(cfg::Config, con::Connection, data)
    try
        JSON3.read(JSON3.write(data))
    catch
        if data isa AbstractArray
            map(d-> json(cfg, con, d), data)
        elseif data isa AbstractDict
            for (k, _) in data
                if !(k isa Symbol || k isa AbstractString)
                    throw("Cannot encode object property $(k)")
                end
            end
            Dict([k => json(cfg, con, v) for (k,v) in data]...)
        else
            cfg.verbose ? (; ref = oid(con, data), repr = repr(data)) : (; ref = oid(con, data))
        end
    end
end

Base.haskey(cfg::Config, id::ID) = haskey(cfg.vars, id)
Base.getindex(cfg::Config, id::ID) = cfg.vars[id]
Base.setindex!(cfg::Config, var::Var, id::ID) = cfg.vars[id] = var

Base.haskey(var::Var, name::Symbol) = haskey(var.namedchildren, name)
function Base.getindex(var::Var, name::Symbol)
    get(var.namedchildren, name, ()-> throw("No variable for name $(name)"))
end
Base.setindex!(var::Var, value::Var, name::Symbol) = var.namedchildren[name] = value
Base.haskey(var::Var, index::Integer) = 0 <= index <= length(var.indexedchildren)
Base.length(var::Var) = length(var.indexedchildren)
function Base.getindex(var::Var, index::Integer)
    get(var.indexedchildren, index, ()-> throw("No variable for index $(index)"))
end
function Base.setindex!(var::Var, value::Var, index::Integer)
    if length(var) == index - 1
        push!(var.indexedchildren, value)
    elseif index < 0 || length(var) < index
        throw("Attempt to add variable $(index), when there are only $(length(var))")
    else
        var.indexedchildren[index] = value
    end
end

within(cfg::Config, id::ID, ids::Set{ID}) = within(cfg, get(cfg.vars, id, EMPTYID), ids)

function within(cfg::Config, var::Var, ids::Set{ID})
    (var.id in ids) || ((var.parent != EMPTYID) && within(cfg, cfg[var.parent], ids))
end

function allvars(cfg::Config, vars...)
    results = []
    traverse(id::ID) = traverse(cfg[id])
    function traverse(var::Var)
        push!(results, var)
        for (_, v) in var.namedchildren
            traverse(v)
        end
        for v in var.indexedchildren
            traverse(v)
        end
    end
    for v in vars
        traverse(v)
    end
    results
end

isemptyid(id) = id isa ID && id.number == 0

function parsemetadata(meta::AbstractString, original_meta = nothing)
    metadata = Dict{Symbol, AbstractString}()
    if original_meta !== nothing
        merge!(metadata, original_meta)
    end
    while meta !== ""
        (m = match(METAPROP, meta)) === nothing && throw("Bad metaproperty format: $(meta)")
        metadata[Symbol(m[1])] = m[2] === nothing ? "" : m[2]
        length(meta) == length(m.match) && break
        meta[length(m.match)] != ',' && throw("Bad metaproperty format: $(meta)")
        meta = meta[length(m.match) + 1:end]
    end
    @debug("USING METADATA: $(repr(metadata))")
    metadata
end

function initial_route_meta(cmd::JusCmd, prop::Symbol, var::Var)
    route(var.value, VarCommand(cmd, :metadata, (prop,), var, creating = true))
end

function Var(cmd::JusCmd; id::ID, name = Symbol(""), args...)
    args = (;args...)
    if haskey(args, :value) && !haskey(args, :internal_value)
        args = (; args..., internal_value = args.value)
    end
    if name isa Union{Symbol, AbstractString} && string(name) != ""
        m = match(VAR_NAME, string(name))
        if m !== nothing && m[2] !== nothing
            name = m[1]
            metadata = parsemetadata(m[2], haskey(args, :metadata) ? args.metadata : Dict())
            args = (;args..., metadata, name)
        end
    end
    if name isa AbstractString
        name = Symbol(name)
    end
    ns = cmd.config.namespaces[id.namespace]
    id.number == 0 && (id = ID(id.namespace, ns.curid += 1))
    ns.curid = max(ns.curid, id.number)
    v = cmd.config[id] = Var(; args..., id, name = name isa Number ? Symbol("") : name)
    v.parent !== EMPTYID && name != Symbol("") && (cmd.config[v.parent][name] = v)
    v.parent !== EMPTYID && name != Symbol("") && @debug("ADDED CHILD OF $(v.parent) NAMED $(name) = $(cmd.config[v.parent][name])")
    !isempty(v.metadata) && @debug("VAR $(v.name) metadata: $(v.metadata)")
    # route handle high priority metadata first
    haskey(v.metadata, :create) && initial_route_meta(cmd, :create, v)
    haskey(v.metadata, :access) && initial_route_meta(cmd, :access, v)
    haskey(v.metadata, :path) && initial_route_meta(cmd, :path, v)
    # route remaining metadata
    for (mk, _) in v.metadata
        mk in (:access, :path, :create) && continue
        initial_route_meta(cmd, mk, v)
    end
    route(v.value, VarCommand(cmd, :create, (), v, creating = true))
    v
end

function next_id(config::Config)
    config.namespace.curid += 1
    ID(config.namespace.name, config.namespace.curid)
end

addvar(parent::Var, value) = addvar(parent, "", next_id(parent.config), value)

addvar(parent::Var, name::Union{Integer, Symbol}, id::ID, value, metadata::Dict{Symbol, AbstractString} = Dict{Symbol, AbstractString}()) =
    addvar(parent.config, parent.id, name, id, value, metadata)

addvar(config::Config, value) = addvar(config, EMPTYID, Symbol(""), next_id(config), value)

addvar(config::Config, parent::ID, name::Union{Integer, Symbol}, id::ID, value, metadata::Dict{Symbol, AbstractString} = Dict{Symbol, AbstractString}()) =
    addvar(JusCmd(config, nothing, config.namespace.name, ["set"]), parent, name, id, value, metadata)

function addvar(cmd::JusCmd, parent::ID, name::Union{Integer, Symbol}, id::ID, value, metadata::Dict{Symbol, AbstractString})
    realname = name == Symbol("") && parent != EMPTYID ? length(cmd.config[parent]) + 1 : name
    v = Var(cmd; id, name = realname, value, metadata, parent)
    metadata[:type] = typename(value)
    changed(cmd.config, v)
    changed(cmd.config, v, keys(v.metadata)...)
    cmd.ws !== nothing && addvar(connection(cmd), v)
    v
end

addvar(con::Connection, var::Var) = push!(con.vars, var)

typename(value) = typename(typeof(value))

typename(type::Type) = "$(Base.parentmodule(type)).$type"

function set_type(cmd::VarCommand)
    name = typename(cmd.var.internal_value)
    if name != get(cmd.var.metadata, :type, "")
        set_metadata(cmd, :type, name)
    end
end

qualify(args) = length(args) == 1 ? Symbol(args[1]) : :($(qualify(args[1:end-1])).$(Symbol(args[end])))

function create_from_metadata(cmd::VarCommand)
    type = cmd.var.metadata[:create]
    match(JPATH_COMPONENT, type) === nothing && throw("bad type for create: $type for $(cmd.var)")
    # safe to eval because type is just "ident(.ident)*"
    cmd.var.internal_value = cmd.var.value = Main.eval(:($(qualify(split(type, ".")))()))
    @debug("@@@ CREATED: $(repr(cmd.var.value))")
    VarCommand(cmd, arg = cmd.var.value)
    set_type(cmd)
end

function set_access_from_metadata(var::Var)
    var.action = var.metadata[:access] == "action"
    var.readable = var.metadata[:access] in ["rw", "r"]
    var.writeable = var.metadata[:access] in ["rw", "w", "action"]
end

function set_path_from_metadata(cmd::VarCommand)
    cmd.var.path = []
    path = cmd.var.path
    for el in split(cmd.var.metadata[:path])
        el == ".." && !isempty(filter(x-> x != :.., path)) &&
            throw("'..' not allowed in the middle of a path in $(cmd.var): $(el)")
        if el == ".."
            push!(path, :..)
        else
            m = match(JPATH_COMPONENT, el)
            m === nothing && throw("Bad path component in $(cmd.var): $(el)")
            if endswith(el, "()")
                try
                    # safe to eval because m[1] * m[2] is just "ident(.ident)*"
                    push!(path, Main.eval(Meta.parse(m[1] * m[2])))
                catch err
                    rethrow(CmdException(:path, cmd, "Bad path for variable $(cmd.var.name): '$(m[1] * m[2])'"))
                end
            elseif match(ARRAY_INDEX, m[1]) !== nothing
                push!(path, parse(Int, m[1]))
            else
                push!(path, Symbol(m[1]))
            end
        end
    end
end

dicty(::Union{AbstractDict, NamedTuple}) = true
dicty(x::T) where T = hasmethod(haskey, Tuple{T, Symbol}) && hasmethod(get, Tuple{T, Symbol})

function basic_get_path(cmd::VarCommand, path)
    var = cmd.var
    var.parent == EMPTYID && throw(CmdException(:path, cmd, "no parent"))
    cur = cmd.config[var.parent].internal_value
    for el in path
        if el == :..
            var = cmd.config[var.parent]
            var.parent == EMPTYID && throw(CmdException(:path, cmd, "error going up in path with no parent"))
            cur = cmd.config[var.parent].internal_value
        elseif cur === nothing
            throw(CmdException(:path, cmd, "error getting $(cmd.var) field $(el) in path $path"))
        elseif el isa Symbol
            try
                if dicty(cur) && haskey(cur, el)
                    cur = get(cur, el)
                elseif cur isa AbstractDict && haskey(cur, string(el))
                    cur = get(cur, string(el))
                else
                    cur = getproperty(cur, el)
                end
            catch err
                rethrow(CmdException(:path, cmd, "error getting $(cmd.var) field $(el) in path $path", err))
            end
        elseif el isa Number
            try
                cur = getindex(cur, el)
            catch err
                rethrow(CmdException(:path, cmd, "error getting $(cmd.var) field $(el) in path $path", err))
            end
        elseif hasmethod(el, typeof.((cmd, cur)))
            try
                cur = el(cmd, cur)
            catch err
                rethrow(CmdException(:program, cmd, "error calling $(cmd.var) getter function $(el) in path $path", err))

            end
        elseif hasmethod(el, typeof.((cur,)))
            try
                cur = el(cur)
            catch err
                rethrow(CmdException(:program, cmd, "error calling $(cmd.var) getter function $(el) in path $path", err))
            end
        else
            throw(CmdException(:program, cmd, "No $(cmd.var) getter method $(el) for $(typeof.((cur,))) in path $path"))
        end
    end
    cur
end

function set_path(cmd::VarCommand)
    cmd.creating && (haskey(cmd.var.metadata, :create) || cmd.var.action || !isempty(cmd.var.path)) && return
    !cmd.var.writeable && throw(CmdException(:writeable_error, cmd, "variable $(cmd.var) is not writeable"))
    cur = basic_get_path(cmd, cmd.var.path[1:end - 1])
    el = cmd.var.path[end]
    cmd.arg = cmd.var.value_conversion(cmd.arg)
    if cur === nothing
        throw(CmdException(:path, cmd, "error setting field $(el) in path $(cmd.var.path) for $(cmd.var)"))
    elseif el isa Symbol
        try
            setproperty!(cur, el, cmd.arg)
        catch err
            rethrow(CmdException(:path, cmd, "error setting $(cmd.var) field $(el)", err))
        end
    elseif el isa Number
        if el == length(cur) + 1
            push!(cur, cmd.arg)
        else
            setindex!(cur, el, cmd.arg)
        end
    elseif cmd.var.action
        parent = cmd.config[cmd.var.parent].internal_value
        if :.. in cmd.var.path && hasmethod(el, typeof.((cmd, cur, parent)))
            try
                el(cmd, cur, parent)
            catch err
                rethrow(CmdException(:program, cmd, "error calling $(cmd.var) action function $(el) for $(typeof.((cmd, cur, cmd.var.internal_value))): $(err)", err))
            end
        elseif hasmethod(el, typeof.((cmd, cur)))
            try
                el(cmd, cur)
            catch err
                rethrow(CmdException(:program, cmd, "error calling $(cmd.var) action function $(el) for $(typeof.((cmd, cur,))): $(err)", err))
            end
        elseif :.. in cmd.var.path && hasmethod(el, typeof.((cur, parent)))
            try
                el(cur, parent)
            catch err
                rethrow(CmdException(:program, cmd, "error calling $(cmd.var) action function $(el) for $(typeof.((cur, cmd.var.internal_value))): $(err)", err))
            end
        elseif hasmethod(el, typeof.((cur,)))
            try
                el(cur)
            catch err
                rethrow(CmdException(:program, cmd, "error calling $(cmd.var) action function $(el) for $(typeof.((cur,))): $(err)", err))
            end
        else
            throw(CmdException(:path, cmd, "no $(cmd.var) action function $(el) for $(typeof.((cur,)))"))
        end
    elseif hasmethod(el, typeof.((cmd, cur, cmd.arg)))
        try
            el(cmd, cur, cmd.arg)
        catch err
            rethrow(CmdException(:program, cmd, "error calling $(cmd.var) setter function $(el): $(err)", err))
        end
    elseif hasmethod(el, typeof.((cur, cmd.arg)))
        try
            el(cur, cmd.arg)
        catch err
            rethrow(CmdException(:program, cmd, "error calling $(cmd.var) setter function $(el): $(err)", err))
        end
    else
        throw(CmdException(:path, cmd, "no $(cmd.var) setter function $(el) for $(typeof.((cur, cmd.arg)))"))
    end
end

function get_path(cmd::VarCommand)
    !cmd.var.readable && throw(CmdException(:readable_error, cmd, "variable $(cmd.var) is not readable"))
    cmd.var.value = cmd.var.internal_value = basic_get_path(cmd, cmd.var.path)
    cmd.var.json_value = JSON3.write(json(cmd, cmd.var.value))
    set_type(cmd)
end


text:
# # **7.** Example
#
# This is an example generated with Literate based on this
# source file: [`example.jl`](@__REPO_ROOT_URL__julia_sandbox/literate_test_2.jl).
# You are seeing the
#md # html-output which Documenter have generated based on a markdown
#md # file generated with Literate. The corresponding notebook
#md # can be found here: [`example.ipynb`](@__NBVIEWER_ROOT_URL__generated/example.ipynb),
#nb # generated notebook output. The corresponding markdown (html) output
#nb # can be found here: [`example.html`](https://fredrikekre.github.io/Literate.jl/latest/generated/example.html),
# and the plain script output can be found here: [`example.jl`](./example.jl).

# It is recommended to have the source file
# available when reading this, to better understand how the syntax in the source file
# corresponds to the output you are seeing.

# ### Basic syntax
# The basic syntax for Literate is simple, lines starting with `# ` is interpreted
# as markdown, and all the other lines are interpreted as code. Here is some code:

x = 1//3
y = 2//5

# In markdown sections we can use markdown syntax. For example, we can
# write *text in italic font*, **text in bold font** and use
# [links](https://www.youtube.com/watch?v=dQw4w9WgXcQ).

# It is possible to filter out lines depending on the output using the
# `#md`, `#nb`, `#jl` and `#src` tags (see [Filtering Lines](@ref)):
#md # - This line starts with `#md` and is thus only visible in the markdown output.
#nb # - This line starts with `#nb` and is thus only visible in the notebook output.
#jl # - This line starts with `#jl` and is thus only visible in the notebook output.
#src # - This line starts with `#src` and is thus only visible in the source file.

# The source file is parsed in chunks of markdown and code. Starting a line
# with `#-` manually inserts a chunk break. For example, if we want to
# display the output of the following operations we may insert `#-` in
# between. These two code blocks will now end up in different
# `@example`-blocks in the markdown output, and two different notebook cells
# in the notebook output.

x + y

#-

x * y

# ### Output Capturing
# Code chunks are by default placed in Documenter `@example` blocks in the generated
# markdown. This means that the output will be captured in a block when Documenter is
# building the docs. In notebooks the output is captured in output cells, if the
# `execute` keyword argument is set to true. Output to `stdout`/`stderr` is also
# captured.

function foo()
    println("This string is printed to stdout.")
    return [1, 2, 3, 4]
end

foo()

# Both Documenter's `@example` block and notebooks can display images. Here is an example
# where we generate a simple plot using the
# [Plots.jl](https://github.com/JuliaPlots/Plots.jl) package

using Plots
x = linspace(0, 6π, 1000)
y1 = sin.(x)
y2 = cos.(x)
plot(x, [y1, y2])

# ### Custom processing
#
# It is possible to give Literate custom pre- and post-processing functions.
# For example, here we insert two placeholders, which we will replace with
# something else at time of generation. We have here replaced our placeholders
# with `z` and `1.0 + 2.0im`:

# Generated at DATETIME.


text:
instructions = readlines("2a.txt")

struct Pos
    x
    z
end

function update_position(pos, direction, amount)
    if direction == "forward"
        Pos(pos.x+amount, pos.z)
    elseif direction == "up"
        Pos(pos.x, pos.z - amount)
    elseif direction == "down"
        Pos(pos.x, pos.z + amount)
    end
end

pos = Pos(0,0)
for instruction in instructions
    direction, amount = split(instruction)
    amount = parse(Int, amount)
    pos = update_position(pos, direction, amount)
end

pos.x*pos.z


text:
# Standard LP interface
importall MathProgBase.SolverInterface

###############################################################################
# Solver objects
export IpoptSolver
immutable IpoptSolver <: AbstractMathProgSolver
  options
end
IpoptSolver(;kwargs...) = IpoptSolver(kwargs)

type IpoptMathProgModel <: AbstractMathProgModel
  inner::Any
  options
end
function IpoptMathProgModel(;options...)
  return IpoptMathProgModel(nothing,options)
end
model(s::IpoptSolver) = IpoptMathProgModel(;s.options...)
export model

###############################################################################
# Begin interface implementation
function loadproblem!(model::IpoptMathProgModel, A, l, u, c, lb, ub, sense)
  Asparse = convert(SparseMatrixCSC{Float64,Int32}, A)
  n = int(Asparse.n)
  m = int(Asparse.m)
  nnz = int(length(Asparse.rowval))
  c_correct = float(c)
  if sense == :Max
    c_correct .*= -1.0
  end

  # Objective callback
  function eval_f(x)
    return dot(x,c_correct)
  end

  # Objective gradient callback
  function eval_grad_f(x, grad_f)
    for j = 1:n
      grad_f[j] = c_correct[j]
    end
  end

  # Constraint value callback
  function eval_g(x, g)
    g_val = A*x
    for i = 1:m
      g[i] = g_val[i]
    end
  end

  # Jacobian callback
  function eval_jac_g(x, mode, rows, cols, values)
    if mode == :Structure
      # Convert column wise sparse to triple format
      idx = 1
      for col = 1:n
        for pos = Asparse.colptr[col]:(Asparse.colptr[col+1]-1)
          rows[idx] = Asparse.rowval[pos]
          cols[idx] = col
          idx += 1
        end
      end
    else
      # Values
      idx = 1
      for col = 1:n
        for pos = Asparse.colptr[col]:(Asparse.colptr[col+1]-1)
          values[idx] = Asparse.nzval[pos]
          idx += 1
        end
      end
    end
  end

  x_L = float(l)
  x_U = float(u)
  g_L = float(lb)
  g_U = float(ub)
  model.inner = createProblem(n, x_L, x_U, m, g_L, g_U, nnz, 0,
                              eval_f, eval_g, eval_grad_f, eval_jac_g, nothing)
  model.inner.sense = sense
  addOption(model.inner, "jac_c_constant", "yes")
  addOption(model.inner, "jac_d_constant", "yes")
  addOption(model.inner, "hessian_constant", "yes")
  addOption(model.inner, "hessian_approximation", "limited-memory")
  addOption(model.inner, "mehrotra_algorithm", "yes")
  for (name,value) in model.options
    addOption(model.inner, string(name), value)
  end
end

# generic nonlinear interface
function loadnonlinearproblem!(m::IpoptMathProgModel, numVar::Integer, numConstr::Integer, x_l, x_u, g_lb, g_ub, sense::Symbol, d::AbstractNLPEvaluator)

    initialize(d, [:Grad, :Jac, :Hess])
    Ijac, Jjac = jac_structure(d)
    Ihess, Jhess = hesslag_structure(d)
    @assert length(Ijac) == length(Jjac)
    @assert length(Ihess) == length(Jhess)
    @assert sense == :Min || sense == :Max

    # Objective callback
    if sense == :Min
        eval_f_cb(x) = eval_f(d,x)
    else
        eval_f_cb(x) = -eval_f(d,x)
    end

    # Objective gradient callback
    if sense == :Min
        eval_grad_f_cb(x, grad_f) = eval_grad_f(d, grad_f, x)
    else
        eval_grad_f_cb(x, grad_f) = (eval_grad_f(d, grad_f, x); scale!(grad_f,-1))
    end

    # Constraint value callback
    eval_g_cb(x, g) = eval_g(d, g, x)

    # Jacobian callback
    function eval_jac_g_cb(x, mode, rows, cols, values)
        if mode == :Structure
            for i in 1:length(Ijac)
                rows[i] = Ijac[i]
                cols[i] = Jjac[i]
            end
        else
            eval_jac_g(d, values, x)
        end
    end

    # Hessian callback
    function eval_h_cb(x, mode, rows, cols, obj_factor,
        lambda, values)
        if mode == :Structure
            for i in 1:length(Ihess)
                rows[i] = Ihess[i]
                cols[i] = Jhess[i]
            end
        else
            if sense == :Max
                obj_factor *= -1
            end
            eval_hesslag(d, values, x, obj_factor, lambda)
        end
    end

    m.inner = createProblem(numVar, float(x_l), float(x_u), numConstr,
                                float(g_lb), float(g_ub), length(Ijac), length(Ihess),
                                eval_f_cb, eval_g_cb, eval_grad_f_cb, eval_jac_g_cb,
                                eval_h_cb)
    m.inner.sense = sense

    for (name,value) in m.options
        addOption(m.inner, string(name), value)
    end
end

getsense(m::IpoptMathProgModel) = m.inner.sense
numvar(m::IpoptMathProgModel) = m.inner.n
numconstr(m::IpoptMathProgModel) = m.inner.m
optimize!(m::IpoptMathProgModel) = solveProblem(m.inner)
function status(m::IpoptMathProgModel)
  if m.inner.status == 0 || m.inner.status == 1
    return :Optimal
  end
  return :Infeasible
end
getobjval(m::IpoptMathProgModel) = m.inner.obj_val * (m.inner.sense == :Max ? -1 : +1)
getsolution(m::IpoptMathProgModel) = m.inner.x
getconstrsolution(m::IpoptMathProgModel) = m.inner.g
getreducedcosts(m::IpoptMathProgModel) = zeros(m.inner.n)
getconstrduals(m::IpoptMathProgModel) = zeros(m.inner.m)
getrawsolver(m::IpoptMathProgModel) = m.inner
setwarmstart!(m::IpoptMathProgModel, x) = copy!(m.inner.x, x) # starting point


text:
function addparameter! end

"""
$(TYPEDEF)

Represents a 'placeholder' for data; a value that may be filled in later.

`Parameter`s can be evaluated by simply calling them with no arguments.

`Parameter`s keep track of whether they've already been evaluated using `dirty` flag.
To reevalute a parameter, the dirty flag must first be set using `setdirty!(parameter)`.
The update function will then be called when the parameter itself is called.

# Examples

```jldoctest; setup = :(using Parametron; model = Parametron.mock_model())
julia> value = Ref(1)
Base.RefValue{Int64}(1)

julia> p = Parameter{Int}(() -> value[], model)
Parameter{Int64, …}(…)

julia> p()
1

julia> value[] = 2
2

julia> p()
1

julia> Parametron.setdirty!(p); p()
2
```
"""
struct Parameter{T, F, InPlace}
    dirty::Base.RefValue{Bool}
    f::F
    val::Base.RefValue{T}

    # out-of-place
    @doc """
    $(SIGNATURES)

    Create a new 'out-of-place' `Parameter` with an update function `f` that takes
    no arguments and returns a value of type `T`.
    """ ->
    Parameter{T}(f::F, model) where {T, F} = addparameter!(model, new{T, F, false}(Ref(true), f, Base.RefValue{T}()))

    # in-place
    @doc """
    $(SIGNATURES)

    Create a new 'in-place' `Parameter` with an update function `f` that takes
    `val` as its argument and updates it in place.
    """ ->
    Parameter(f::F, val::T, model) where {T, F} = addparameter!(model, new{T, F, true}(Ref(true), f, Base.RefValue(val)))
end

"""
$(SIGNATURES)

Create a new 'out-of-place' `Parameter` with an update function `f` that takes
no arguments. The type of the output is determined upon construction by calling
the update function.

!!! warning

    Explicitly specifying the return value type using the `Parameter{T}(f, model)`
    constructor is preferred, as using this constructor can lead to type inference
    issues.
"""
Parameter(f, model) = Parameter{typeof(f())}(f, model)

"""
$(SIGNATURES)

Create a new 'in-place' `Parameter` that always returns `val`. This constructor
may be used to create `Parameter`s that use `val` as a work buffer that is
manually/externally updated.

!!! warning

    By using this constructor, the automated mechanism for lazily updating a `Parameter`'s
    value when necessary is essentially circumvented, and the user is responsible for
    ensuring that the `Parameter`'s value is updated at the appropriate time.
"""
Parameter(model; val::T) where {T} = Parameter(identity, val, model)

isinplace(::Type{Parameter{T, F, InPlace}}) where {T, F, InPlace} = InPlace
Base.show(io::IO, param::Parameter{T, F, InPlace}) where {T, F, InPlace} = print(io, "Parameter{$T, …}(…)")

@inline function (parameter::Parameter{T})() where {T}
    if parameter.dirty[]
        update!(parameter)
        parameter.dirty[] = false
    end
    parameter.val[]::T
end

update!(parameter::Parameter{T, F, true}) where {T, F} = (parameter.f(parameter.val[]); nothing)
update!(parameter::Parameter{T, F, false}) where {T, F} = (parameter.val[] = parameter.f()::T; nothing)

setdirty!(parameter::Parameter) = (parameter.dirty[] = true; nothing)


text:
using PowerVega
using Test

using PowerModels
# import Ipopt

data = PowerModels.parse_file("$(joinpath(dirname(pathof(PowerModels)), ".."))/test/data/matpower/case5.m")

@testset "PowerVega.jl" begin
    # # Write your own tests here.
    #  graph = build_graph_network(data; exclude_sources=true)
    #  plot_network(data)
    #  plot_network!(data)

    #  plot_network_status!(graph)
    #  plot_system_voltage(data)

    #  result = run_opf(data, DCPPowerModel, Ipopt.Optimizer)
    #  PowerModels.update_data!(data, result["solution"])
    #  plot_power_flow(data)

    #  # these functions work if we got this far
    #  @test true

     ## Test vegalite plotting
     case = PowerModels.parse_file("$(joinpath(dirname(pathof(PowerModels)), ".."))/test/data/matpower/case5.m")
    #  plot_vega(case)
     @test true

end


text:
module ModelVisualisations

    using Plots#, HeavyChain

    include("src/invpendulum.jl")
    include("src/ballbot.jl")
    #include("src/heavychain.jl")
    include("src/segwayrobot.jl")

    """
    *function visualise(modelname::Symbol, X::Array{Float64, 2}, t::AbstractVector, filename::String,; options::Dict{Symbol})*\\
    Visualise simulationdata and save .mp4 file.
    Valid modelnames are :InvPendulum, :BallBot, :HeavyChain, :SegwayRobot .
    Valid options are:\\
    - fps: the fps of the video, defaults to 25
    - U: provide input data, which will be plotted below the main window #TODO
    - title: provide a title for the window, simulation time will be added
    - size: resolution of the whole plot in pixels
    """

    function visualise(modelname::Symbol, X::Array{Float64, 2}, t::AbstractVector, filename::String, des_traj::Array{Float64,2}; options::Dict{Symbol}=Dict{Symbol, Any}())
        @assert(issorted(t))
        @assert(length(t)==size(X, 2))

        if modelname == :InvPendulum
            @assert(size(X, 1)==4)
            visualise_invpendulum(X, t, filename; options = options)

        elseif modelname == :BallBot
            @assert(size(X, 1)==12)
            visualise_ballbot(X, t, filename, des_traj; options = options)

        elseif modelname == :HeavyChain
            @assert(size(X, 1)==21)
            visualise_heavychain(X, t, filename; options = options)

        elseif modelname == :SegwayRobot
            @assert(size(X, 1)==7)
            visualise_segwayrobot(X, t, filename, des_traj; options = options)

        end
        return nothing
    end

    export visualise, plot_all_segways, plot_all_ballbots,  visualise_segwayrobot_compare, visualise_ballbot_compare

end # module

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
module FGPEexamples

using Pkg, Weave

cssfile = joinpath(@__DIR__,"templates","skeleton_css.css")
latexfile = joinpath(@__DIR__, "../docs", "templates", "julia_tex.tpl")

source_path = @__DIR__
html_path = joinpath(source_path,"../docs","html")
pdf_path = joinpath(source_path,"../docs","pdf")
notebook_path = joinpath(source_path,"../docs","notebooks")

# flist = ["1dharmonic" "1dbrightsoliton" "2dharmonic" "3dquenchiso" "3dquenchtrap" "3dquenchslab" "3dquenchtube" "jrsoliton" "1dtemporalsoliton"]
flist = ["jrsoliton"] # to only build this example

function weavedocs(file)
    tmp = joinpath(source_path,file)*".jmd"
    weave(tmp,out_path=html_path,doctype = "md2html")
    # Weave.convert_doc(tmp,joinpath(notebook_path,file)*".ipynb")
    # weave(tmp,out_path=pdf_path,doctype="md2pdf"; template=latexfile)
end

function weave_some(flist)
    for file in flist
        weavedocs(file)
    end
end

weave_all() = weave_some(flist)
end


text:
@testset "Small graphs" begin
    g = smallgraph("diamondgraph")
    @test nv(g) == 4 && ne(g) == 5

    g = smallgraph("diamond")
    @test nv(g) == 4 && ne(g) == 5

    g = smallgraph(:diamond)
    @test nv(g) == 4 && ne(g) == 5

    g = smallgraph(:bull)
    @test nv(g) == 5 && ne(g) == 5

    g = smallgraph(:chvatal)
    @test nv(g) == 12 && ne(g) == 24

    g = smallgraph(:cubical)
    @test nv(g) == 8 && ne(g) == 12

    g = smallgraph(:desargues)
    @test nv(g) == 20 && ne(g) == 30

    g = smallgraph(:dodecahedral)
    @test nv(g) == 20 && ne(g) == 30

    g = smallgraph(:frucht)
    @test nv(g) == 12 && ne(g) == 18

    g = smallgraph(:heawood)
    @test nv(g) == 14 && ne(g) == 21

    g = smallgraph(:house)
    @test nv(g) == 5 && ne(g) == 6

    g = smallgraph(:housex)
    @test nv(g) == 5 && ne(g) == 8

    g = smallgraph(:icosahedral)
    @test nv(g) == 12 && ne(g) == 30

    g = smallgraph(:krackhardtkite)
    @test nv(g) == 10 && ne(g) == 18

    g = smallgraph(:moebiuskantor)
    @test nv(g) == 16 && ne(g) == 24

    g = smallgraph(:octahedral)
    @test nv(g) == 6 && ne(g) == 12

    g = smallgraph(:pappus)
    @test nv(g) == 18 && ne(g) == 27

    g = smallgraph(:petersen)
    @test nv(g) == 10 && ne(g) == 15

    g = smallgraph(:sedgewickmaze)
    @test nv(g) == 8 && ne(g) == 10

    g = smallgraph(:tetrahedral)
    @test nv(g) == 4 && ne(g) == 6

    g = smallgraph(:truncatedcube)
    @test nv(g) == 24 && ne(g) == 36

    g = smallgraph(:truncatedtetrahedron)
    @test nv(g) == 12 && ne(g) == 18 && !is_directed(g)

    g = smallgraph(:truncatedtetrahedron_dir)
    @test nv(g) == 12 && ne(g) == 18 && is_directed(g)

    g = smallgraph(:tutte)
    @test nv(g) == 46 && ne(g) == 69

    @test_throws ArgumentError g = smallgraph(:nonexistent)
end


text:
module _K2En
export K2En
using Statistics: std
    """
        K2, Ci = K2En(Sig)

    Returns the Kolmogorov entropy estimates `K2` and the correlation
    integrals `Ci` for `m` = [1,2] estimated from the data sequence `Sig`
    using the default parameters: embedding dimension = 2, time delay = 1,
    r = 0.2*SD(`Sig`), logarithm = natural

        K2, Ci = K2En(Sig::AbstractArray{T,1} where T<:Real; m::Int=2, tau::Int=1, r::Real=0.2*std(Sig,corrected=false), Logx::Real=exp(1))

    Returns the Kolmogorov entropy estimates `K2` for dimensions = [1,...,`m`]
    estimated from the data sequence `Sig` using the 'keyword' arguments:

    # Arguments:
    `m`     - Embedding Dimension, a positive integer\n
    `tau`   - Time Delay, a positive integer\n
    `r`     - Radius, a positive scalar \n
    `Logx`  - Logarithm base, a positive scalar\n

    # See also `DistEn`, `XK2En`, `MSEn`

    # References:
        [1] Peter Grassberger and Itamar Procaccia,
            "Estimation of the Kolmogorov entropy from a chaotic signal."
            Physical review A 28.4 (1983): 2591.

        [2] Lin Gao, Jue Wang  and Longwei Chen
            "Event-related desynchronization and synchronization
            quantification in motor-related EEG by Kolmogorov entropy"
            J Neural Eng. 2013 Jun;10(3):03602

    """
    function K2En(Sig::AbstractArray{T,1} where T<:Real; m::Int=2, tau::Int=1, r::Real=0.2*std(Sig,corrected=false), Logx::Real=exp(1))

    N = size(Sig)[1]
    (N>10) ? nothing : error("Sig:   must be a numeric vector")
    (m > 0) ? nothing : error("m:     must be an integer > 0")
    (tau > 0) ? nothing : error("tau:   must be an integer > 0")
    (r>0) ? nothing :  error("r:     must be a positive scalar > 0")
    (Logx>0) ? nothing :  error("Logx:     must be a positive number > 0")

    m = m+1
    Zm = zeros(N,m)
    Ci = zeros(m)
    for n = 1:m
        N2 = N-(n-1)*tau
        Zm[1:N2,n] = Sig[(n-1)*tau + 1:N]
        Norm = Inf*ones(N2-1,N2-1)
        for k = 1:N2-1
            Temp = repeat(transpose(Zm[k,1:n]),N2-k,1) - Zm[k+1:N2,1:n]
            Norm[k,k:N2-1] = sqrt.(sum(Temp.^2, dims=2))
        end
        Ci[n] = 2*sum(Norm .< r)/(N2*(N2-1))
    end

    K2 = log.(Logx, Ci[1:m-1]./Ci[2:m])/tau
    K2[isinf.(K2)] .= NaN

    return K2, Ci
    end
end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
using Printf
using LinearAlgebra
BLAS.set_num_threads(1)  # 1 thread only

import Tulip

function run_rmp_tulip(finst)
    tlp = Tulip.Model{Float64}()
    Tulip.load_problem!(tlp, finst)

    # Set parameters
    tlp.params.OutputLevel = 1
    tlp.params.Presolve = 0

    # Solve model
    t = @elapsed Tulip.optimize!(tlp)

    println("Solve time: $t")
    println("Barrier iterations: ", Tulip.get_attribute(tlp, Tulip.BarrierIterations()))
    return nothing
end

if abspath(PROGRAM_FILE) == @__FILE__
    run_rmp_tulip(ARGS[1])
end

max_stars_repo_name: mtanneau/LPBenchmarks

max_issues_repo_name: mtanneau/LPBenchmarks

max_forks_repo_name: mtanneau/LPBenchmarks


text:
## There are many missing in the data. The following functions help dealing with missings.
function check_missing(x)
    for i in x
        if ismissing(i) ==true
            return true
        end
    end
    return false
end

function count_missing(x)
    """Counts the number of missing values in an array"""
        count = 0
        for y in x
            if ismissing(y)==true
            count = count+1
            end
        end
        return count
    end

function max_missing(x)
"""Finds the maximum value of a array where missing are present."""
    y = replace(copy(x), missing => 0)
    return(maximum(y))
end

function mean_missing(x)
"""Finds the mean value of a array where missing are present."""
    flat = vcat(x...)
    missings = count_missing(flat)
    y = replace(x, missing => 0)
    return sum(y)/(length(flat)- missings)
end

function replace_missing(datacube, replacement = 0)
"""Replaces all the missing values in a 3D array with replacement value provided as a parameter."""
    for i in 1:size(datacube)[1]
        for j in 1:size(datacube)[2]
            for k in 1:size(datacube)[3]
                if ismissing(datacube[i,j,k])
                    datacube[i,j,k] = replacement
                end
            end
        end
    end
    return datacube
end


text:
###################################################################################################
struct CHV{Tode} <: AbstractCHVIterator
	ode::Tode	# ODE solver to use for the flow in between jumps
end

function (chv::CHV{Tode})(xdot, x, caract::PDMPCaracteristics, t) where {Tode}
	tau = x[end]
	rate = get_rate(caract.ratecache, x)
	sr = caract.R(rate, x, caract.xd, caract.parms, tau, true)[1]
	caract.F(xdot, x, caract.xd, caract.parms, tau)
	xdot[end] = 1.0
	@inbounds for i in eachindex(xdot)
		xdot[i] = xdot[i] / sr
	end
	return nothing
end

###################################################################################################
### implementation of the CHV algo using DiffEq

# The following function is a callback to discrete jump. Its role is to perform the jump on the solution given by the ODE solver
# callable struct
function chvjump(integrator, prob::PDMPProblem, save_pre_jump, save_rate, verbose)
	# we declare the characteristics for convenience
	caract = prob.caract
	ratecache = caract.ratecache
	simjptimes = prob.simjptimes

	# final simulation time
	tf = prob.tspan[2]

	# find the next jump time
	t = integrator.u[end]
	simjptimes.lastjumptime = t

	verbose && printstyled(color=:green, "--> Jump detected at t = $t !!\n")
	verbose && printstyled(color=:green, "--> jump not yet performed, xd = ", caract.xd,"\n")

	if (save_pre_jump) && (t <= tf)
		verbose && printstyled(color=:green, "----> saving pre-jump\n")
		pushXc!(prob, (integrator.u[1:end-1]))
		pushXd!(prob, copy(caract.xd))
		pushTime!(prob, t)
		#save rates for debugging
		save_rate && push!(prob.rate_hist, sum(ratecache.rate))
	end

	# execute the jump
	caract.R(get_rate(ratecache, integrator.u), integrator.u, caract.xd, caract.parms, t, false)
	if (t < tf)
		#save rates for debugging
		save_rate && push!(prob.rate_hist, sum(ratecache.rate))

		# Update event
		ev = pfsample(ratecache.rate)

		# we perform the jump
		affect!(caract.pdmpjump, ev, integrator.u, caract.xd, caract.parms, t)

		u_modified!(integrator, true)

		@inbounds for ii in eachindex(caract.xc)
			caract.xc[ii] = integrator.u[ii]
		end
	end
	verbose && printstyled(color=:green,"--> jump computed, xd = ",caract.xd,"\n")
	# we register the next time interval to solve the extended ode
	simjptimes.njumps += 1
	simjptimes.tstop_extended += -log(rand())
	add_tstop!(integrator, simjptimes.tstop_extended)
	verbose && printstyled(color=:green,"--> End jump\n\n")
end

function chv_diffeq!(problem::PDMPProblem,
			ti::Tc, tf::Tc, X_extended::vece,
			verbose = false; ode = Tsit5(), save_positions = (false, true), n_jumps::Td = Inf64, reltol=1e-7, abstol=1e-9, save_rate = false, finalizer = finalizer) where {Tc, Td, vece}
	verbose && println("#"^30)
	verbose && printstyled(color=:red,"Entry in chv_diffeq\n")

	ti, tf = problem.tspan
	algopdmp = CHV(ode)

	# initialise the problem. If I call twice this solve function, it should give the same result...
	init!(problem)

	# we declare the characteristics for convenience
	caract = problem.caract
	ratecache = caract.ratecache
	simjptimes = problem.simjptimes

#ISSUE HERE, IF USING A PROBLEM p MAKE SURE THE TIMES in p.sim ARE WELL SET
	# set up the current time as the initial time
	t = ti
	# previous jump time, needed because problem.simjptimes.lastjumptime contains next jump time even if above tf
	tprev = t

	# vector to hold the state space for the extended system
	# X_extended = similar(problem.xc, length(problem.xc) + 1)
	# @show typeof(X_extended) vece

	for ii in eachindex(caract.xc)
		X_extended[ii] = caract.xc[ii]
	end
	X_extended[end] = ti

	# definition of the callback structure passed to DiffEq
	cb = DiscreteCallback(problem, integrator -> chvjump(integrator, problem, save_positions[1], save_rate, verbose), save_positions = (false, false))

	# define the ODE flow, this leads to big memory saving
	# prob_CHV = ODEProblem((xdot,x,data,tt) -> problem(xdot, x, data, tt), X_extended, (0.0, 1e9))
	prob_CHV = ODEProblem((xdot, x, data, tt) -> algopdmp(xdot, x, caract, tt), X_extended, (0.0, 1e9))
	integrator = init(prob_CHV, ode, tstops = simjptimes.tstop_extended, callback = cb, save_everystep = false, reltol = reltol, abstol = abstol, advance_to_tstop = true)

	# current jump number
	njumps = 0
	simjptimes.njumps = 1

	while (t < tf) && simjptimes.njumps < n_jumps
		verbose && println("--> n = $(problem.simjptimes.njumps), t = $t, δt = ", simjptimes.tstop_extended)
		step!(integrator)

		@assert( t < simjptimes.lastjumptime, "Could not compute next jump time $(simjptimes.njumps).\nReturn code = $(integrator.sol.retcode)\n $t < $(simjptimes.lastjumptime),\n solver = $ode. dt = $(t - simjptimes.lastjumptime)")
		t, tprev = simjptimes.lastjumptime, t

		# the previous step was a jump! should we save it?
		if njumps < simjptimes.njumps && save_positions[2] && (t <= tf)
			verbose && println("----> save post-jump, xd = ",problem.Xd)
			pushXc!(problem, copy(caract.xc))
			pushXd!(problem, copy(caract.xd))
			pushTime!(problem, t)
			njumps +=1
			verbose && println("----> end save post-jump, ")
		end
		finalizer(ratecache.rate, caract.xc, caract.xd, caract.parms, t)
	end
	# we check that the last bit [t_last_jump, tf] is not missing
	if t>tf
		verbose && println("----> LAST BIT!!, xc = ", caract.xc[end], ", xd = ", caract.xd, ", t = ", problem.time[end])
		prob_last_bit = ODEProblem((xdot,x,data,tt) -> caract.F(xdot, x, caract.xd, caract.parms, tt), copy(caract.xc), (tprev, tf))
		sol = DiffEqBase.solve(prob_last_bit, ode)
		verbose && println("-------> xc[end] = ",sol.u[end])
		pushXc!(problem, sol.u[end])
		pushXd!(problem, copy(caract.xd))
		pushTime!(problem, sol.t[end])
	end
	return PDMPResult(problem, save_positions)
end

function solve(problem::PDMPProblem{Tc, Td, vectype_xc, vectype_xd, Tcar}, algo::CHV{Tode}, X_extended; verbose = false, n_jumps = Inf64, save_positions = (false, true), reltol = 1e-7, abstol = 1e-9, save_rate = false, finalizer = finalize_dummy) where {Tc, Td, vectype_xc, vectype_xd, vectype_rate, Tnu, Tp, TF, TR, Tcar, Tode <: DiffEqBase.DEAlgorithm}

	return chv_diffeq!(problem, problem.tspan[1], problem.tspan[2], X_extended, verbose; ode = algo.ode, save_positions = save_positions, n_jumps = n_jumps, reltol = reltol, abstol = abstol, save_rate = save_rate, finalizer = finalizer)
end

"""
	solve(problem::PDMPProblem, algo; verbose = false, n_jumps = Inf64, save_positions = (false, true), reltol = 1e-7, abstol = 1e-9, save_rate = false, finalizer = finalize_dummy)

Simulate the PDMP `problem` using the CHV algorithm.

# Arguments
- `problem::PDMPProblem`
- `alg` can be `CHV(ode)` (for the [CHV algorithm](https://arxiv.org/abs/1504.06873)), `Rejection(ode)` for the Rejection algorithm and `RejectionExact()` for the rejection algorithm in case the flow in between jumps is known analytically. In this latter case, `prob.F` is used for the specification of the Flow. The ODE solver `ode` can be any solver of [DifferentialEquations.jl](https://github.com/JuliaDiffEq/DifferentialEquations.jl) like `Tsit5()` for example or anyone of the list `[:cvode, :lsoda, :adams, :BDF, :euler]`. Indeed, the package implement an iterator interface which does not work yet with `ode = LSODA()`. In order to have access to the ODE solver `LSODA()`, one should use `ode = :lsoda`.
- `verbose` display information during simulation
- `n_jumps` maximum number of jumps to be computed
- `save_positions` which jump position to record, pre-jump (save_positions[1] = true) and/or post-jump (save_positions[2] = true).
- `reltol`: relative tolerance used in the ODE solver
- `abstol`: absolute tolerance used in the ODE solver
- `ind_save_c`: which indices of `xc` should be saved
- `ind_save_d`: which indices of `xd` should be saved
- `save_rate = true`: requires the solver to save the total rate. Can be useful when estimating the rate bounds in order to use the Rejection algorithm as a second try.
-  `X_extended = zeros(Tc, 1 + 1)`: (advanced use) options used to provide the shape of the extended array in the [CHV algorithm](https://arxiv.org/abs/1504.06873). Can be useful in order to use `StaticArrays.jl` for example.
-  `finalizer = finalize_dummy`: allows the user to pass a function `finalizer(rate, xc, xd, p, t)` which is called after each jump. Can be used to overload / add saving / plotting mechanisms.

!!! note "Solvers for the `DiffEqJump` wrapper"
    We provide a basic wrapper that should work for `VariableJumps` (the other types of jumps have not been thoroughly tested). You can use `CHV` for this type of problems. The `Rejection` solver is not functional yet.

"""
function solve(problem::PDMPProblem{Tc, Td, vectype_xc, vectype_xd, Tcar}, algo::CHV{Tode}; verbose = false, n_jumps = Inf64, save_positions = (false, true), reltol = 1e-7, abstol = 1e-9, save_rate = false, finalizer = finalize_dummy) where {Tc, Td, vectype_xc, vectype_xd, vectype_rate, Tnu, Tp, TF, TR, Tcar, Tode <: DiffEqBase.DEAlgorithm}

	# resize the extended vector to the proper dimension
	X_extended = zeros(Tc, length(problem.caract.xc) + 1)

	return chv_diffeq!(problem, problem.tspan[1], problem.tspan[2], X_extended, verbose; ode = algo.ode, save_positions = save_positions, n_jumps = n_jumps, reltol = reltol, abstol = abstol, save_rate = save_rate, finalizer = finalizer )
end


text:
using Revise
using gtpmip

grounded_plan = solve_hpd("experiments/real/domain.hpd",
                        "experiments/real/serve_problem.hpd"; max_levels=400, has_placement_constraint=false)
grounded_plan.plan