text:
module UPOs

export PO,search_in_interval,approximate_period,
        get_period,follow_PO_family_from_period,
        find_p_zero,integrate,monodromy_method_constant_period,
        monodromy_method_constant_energy,po_coordinates,QP,qp,QPp,QPq,mirror_Qq,
        mirror_Pp,scarring_measure,lyapunov,follow_PO_family_from_energy,energy,
        family_A,family_B
    using ..ClassicalSystems
    using ..ClassicalDicke
    using LinearAlgebra
    using DiffEqCallbacks,OrdinaryDiffEq, DiffEqBase
    using ProgressMeter
    using ..PhaseSpaces
    using ..DickeBCE
    using ..EnergyShellProjections
    Id=Matrix{Float64}(I, 4, 4);
    """
    ```julia
    struct PO
    ```
    This struct represents a periodic orbit in the classical Dicke model. To generate,
    use
    ```julia
    po = PO(system, u, T)
    ```
    where `system` is an instance of [`ClassicalDicke.DickeSystem`](@ref),
    `u` is a real vector of an initial condition in the form `[Q,q,P,p]`, and `T`
    is a real number representing the period. You can retrieve this values using
    `po.system`, `po.u`, and `po.T`.
    """
    struct PO{System <: DickeSystem}
        system::System
        u::Vector{Float64}
        T::Real
        PO{System}(system::System, u::AbstractVector{<:Real},
            T::Real) where {System<:DickeSystem} = new(system,Float64.(u),T)
    end
    PO(system::System, u::AbstractVector{<:Real},
        T::Real) where {System<:DickeSystem} = PO{System}(system,u,T)

    function Base.show(io::IO, m::PO)
         print(io,string("PO @ u=",m.u,", T=", m.T))
    end
    """
    ```julia
    Base.:(==)(po1::PO,po2::PO;Ttol::Real=1e-5,tol::Real=1e-6)
    ```
    The comparison `po1 == po2` returns `true` if `po1` and `po2` represent the same periodic orbit, regardless
    of which initial condition they have.
    This is done by evaluating `∈(po1.u,  po2, tol = tol) && abs(po1.T-po2.T) < Ttol`
    """
    function Base.:(==)(po1::PO,po2::PO;Ttol::Real=1e-5,tol::Real=1e-6)
        if abs(po1.T-po2.T)>= Ttol
            return false
        end
        return ∈(po1.u,  po2,tol = tol)
    end
    """
    ```julia
    Base.in(u::AbstractVector{<:Number}, po::PO;tol::Real=1e-6)
    Base.∈(u::AbstractVector{<:Number}, po::PO;tol::Real=1e-6)
    ```

    You may write `u ∈ po`.

    Returns `true` if the point `u` is part of the periodic orbit `po`. This
    is done by integrating `po.u` and seeing if it comes within `tol` of `u`.
    """
    function Base.in(u::AbstractArray{<:Number,1}, po::PO;tol::Real=1e-6)
        if sum(abs2,u-po.u)<tol
            return true
        end

        result=false

        function distance(uvar,t,integrator)
            m=sum(abs2,uvar-u)
            if m<tol
                result=true
                terminate!(integrator)
            end
            return m*sign(uvar[4]-u[4])
        end
        cb=ContinuousCallback(distance,integrator->distance(integrator.u,integrator.t,integrator),save_positions=(false,false),rootfind=true,interp_points=10,reltol=10^-8,abstol=10^-8)
        ClassicalSystems.integrate(po;tol=1e-8,save_everystep=false,callback=cb)
        return result
    end
    """
    ```julia
    function PO(system::DickeSystem, u::AbstractVector{<:Real})
    ```
    Generates a [`PO`](@ref), given a periodic condition `u = [Q,q,P,p]`, where
    the period is calculated using [`approximate_period(system,u,bound=1e-2)`](@ref).
    """
    PO(system::DickeSystem,u::AbstractVector{<:Real})=PO(system,u,approximate_period(system,u,bound=1e-2))

    """
    ```julia
    function integrate(po::PO;tol::Real=1e-16, kargs...)
    ```
    Returns an instance of [`OrdinaryDiffEq.ODESolution`](https://diffeq.sciml.ai/dev/basics/solution/), resulting
    from integrating `po.u` from `t = 0` to `t = po.T`.

    # Arguments
    - `po` should be an instance of [`PO`](@ref).
    # Keyword arguments
    - `tol` is a real number indicating the precision for the integrator. Defaults to `1e-16`.
    - `kargs` are redirected to [`ClassicalSystems.integrate`](@ref ClassicalSystems.integrate(::ClassicalSystems.ClassicalSystem,::AbstractVector{<:Real},::Real)).
    """
    ClassicalSystems.integrate(po::PO;tol::Real=1e-16, kargs...)=ClassicalSystems.integrate(po.system,po.u,po.T;tol=tol, kargs...);

    """
    ```julia
    function action(po::PO;kargs...)
    ```
    Numerically computes the [action](https://en.wikipedia.org/wiki/Action_(physics)#Abbreviated_action_(functional))
    ```math
        S = \\int p\\,\\text{d} q + P \\, \\text{d} Q
    ```
    along the periodic orbit `po`.
    # Arguments
    - `po` should be an instance of [`PO`](@ref)
    # Keyword arguments
    - `kargs` are redirected to [`integrate`](@ref integrate(::ClassicalSystem,::AbstractVector{<:Real},::Real))
    """
    function action(po::PO;kargs...)
        us=integrate(po;kargs...).u
        return sum(us[i][3]*(us[i][1]-us[i-1][1]) + us[i][4]*(us[i][2]-us[i-1][2]) for i in 2:length(us))
    end
    """
    ```julia
    function average_over_PO(po::Union{OrdinaryDiffEq.ODESolution,PO},
        f::Function;
        kargs...)
    ```
    Computes the average of a function ``f`` over the periodic orbit `po`
    ```math
        \\frac{1}{T} \\int_0^T f(u(t)) \\, \\text{d} t,
    ```
    where ``T=`` `po.T` and ``u=`` `po.u`.
    # Arguments
    - `po` should be an instance of [`PO`](@ref) or an instance of [`OrdinaryDiffEq.ODESolution`](https://diffeq.sciml.ai/dev/basics/solution/),
      in which case the integration is ommited.
    - `f` should be a function with a method `f([Q,q,P,p])`, which returns values that
      may be added together and multiplied by scalars (e.g. numbers or arrays).
    # Keyword arguments
    - `kargs` are redirected to  [`ClassicalSystems.integrate`](@ref ClassicalSystems.integrate(::ClassicalSystems.ClassicalSystem,::AbstractVector{<:Real},::Real)).
    """
    function average_over_PO(po::Union{OrdinaryDiffEq.ODESolution,PO},
        f::Function;
        tol::Real=1e-12,
        kargs...)

        tot=nothing

        if isa(po,PO)
            function guardar(u,t,integrator)
                v=f(u)*(t-integrator.tprev)/po.T
                if tot===nothing
                    tot=v
                else
                    tot+=v
                end
            end
            cb=FunctionCallingCallback(guardar;func_start = false,func_everystep =true)
            ClassicalSystems.integrate(po;kargs...,tol=tol,save_everystep=false,callback=cb)
        else
            T=po.t[end]-po.t[1]
            for i in 2:length(po.t)
                v=f(po.u[i])*(po.t[i]-po.t[i-1])/T
                if tot===nothing
                    tot=v
                else
                    tot+=v
                end
            end
        end
        return tot
    end
    """
    ```julia
    function average_over_PO_QP(po::Union{OrdinaryDiffEq.ODESolution,PO},
        f::Function;
        kargs...)
    ```
    Same as [`average_over_PO`](@ref), but `f` only depends on `Q,P`, that is, the signature
    of `f` is `f(Q,P)`.
    """
    average_over_PO_QP(PO::PO,f::Function;kargs...)=average_over_PO(PO,x->f(x[1],x[3]);kargs...)
    """
    ```julia
    function average_over_PO_qp(po::Union{OrdinaryDiffEq.ODESolution,PO},
        f::Function;
        kargs...)
    ```
    Same as [`average_over_PO`](@ref), but `f` only depends on `q,p`, that is, the signature
    of `f` is `f(q,p)`.
    """
    average_over_PO_qp(PO::PO,f::Function;kargs...)=average_over_PO(PO,x->f(x[2],x[4]);kargs...)

    """
    ```julia
    function jz_PO_average(po::Union{OrdinaryDiffEq.ODESolution,PO},
        f::Function;
        kargs...)
    ```
    Returns `average_over_PO_QP(PO,PhaseSpaces.jz;kargs...)`.
    """
    jz_PO_average(PO::PO;kargs...)=average_over_PO_QP(PO,PhaseSpaces.jz;kargs...)

    """
    ```julia
    function jy_PO_average(po::Union{OrdinaryDiffEq.ODESolution,PO},
        f::Function;
        kargs...)
    ```
    Returns `average_over_PO_QP(PO,PhaseSpaces.jy;kargs...)`.
    """
    jy_PO_average(PO::PO;kargs...)=average_over_PO_QP(PO,PhaseSpaces.jy;kargs...)

    """
    ```julia
    function jx_PO_average(po::Union{OrdinaryDiffEq.ODESolution,PO},
        f::Function;
        kargs...)
    ```
    Returns `average_over_PO_QP(PO,PhaseSpaces.jx;kargs...)`.
    """
    jx_PO_average(PO::PO;kargs...)=average_over_PO_QP(PO,PhaseSpaces.jx;kargs...)

    """
    ```julia
    function mirror_Qq(po::PO)
    ```
    Returns a periodic orbit `po1`, which results from changing the sign
    of the ``Q`` and ``q`` coordinates of the periodic orbit `po`.
    """
    mirror_Qq(po::PO)=PO(po.system,po.u.*[-1,-1,1,1],po.T)

    """
    ```julia
    function mirror_Pp(po::PO)
    ```
    Returns a periodic orbit `po1`, which results from changing the sign
    of the ``P`` and ``p`` coordinates of the periodic orbit `po`.
    """
    mirror_Pp(po::PO)=PO(po.system,po.u.*[1,1,-1,-1],po.T)

    """
    ```julia
    function approximate_period(system::DickeSystem,
        u₀::AbstractVector{<:Real};
        bound::Real = 0.1,
        tol::Real = 1e-8,
        max_p_crossings::Real = Inf,
        verbose::Bool = true)
    ```
    Computes the time necessary for the initial condition `u₀` to come back to the
    same `p`-plane within a neighborhood of radius `bound` of `u₀`. This is useful
    to approximate the period of an *almost* periodic condition `u₀`. The maximum time
    of integration is `10000`.
    # Arguments
    - `system` should be an instance of [`ClassicalDicke.DickeSystem`](@ref).
    - `u₀` is a point `[Q,q,P,p]` in the phase space of the Dicke model.
    # Keyword arguments
    - `bound` is a positive real number indicating how close does the evolution has
      to come back to `u₀`. Defaults to `0.1`.
    - `tol` is the integration tolerance. Not to be confused with `bound`. Defaults to  1e-8.
    - `max_p_crossings` is the maximum number of times the condition may cross the `p` plane
      before aborting. Defaults to `Inf` (no maximum).
    - `verbose` is a boolean indicating whether to log information messages. Defaults to `true`.
    """
    function approximate_period(system::DickeSystem,
        u₀::AbstractVector{<:Real};
        bound::Real = 0.1,
        tol::Real = 1e-8,
        max_p_crossings::Real = Inf,
        verbose::Bool = true)

        notreturnedcount=0
        lastsavet=0
        function guardar(integrator)
            c=integrator.u[:]

            returned=bound > norm(c-u₀)
            if !returned
                notreturnedcount+=1
                if notreturnedcount>=max_p_crossings
                    terminate!(integrator)
                    if verbose
                        @info "Exceeded max_p_crossings = $max_p_crossings"
                    end
                    return nothing
                end
            else
                lastsavet = integrator.t
                terminate!(integrator)
                return nothing

            end
            u_modified!(integrator,false)
            return nothing
        end
        cb=ContinuousCallback((uvar,t,integrator)-> (uvar[4]-u₀[4]),nothing,guardar,save_positions=(false,false),rootfind=true,interp_points=3,reltol=10^-8,abstol=10^-8)
        t=ClassicalSystems.integrate(system,u₀,10000;callback=cb,save_everystep=false,tol=tol).t[end]
        return lastsavet

    end

    function monodromy_method_step_constant_period(system::DickeSystem,
        u₀::AbstractVector{<:Real},
        T::Real;
        tol::Real=1e-12)

        u₁,M=ClassicalSystems.integrate(system,u₀,T;save_everystep=false,get_fundamental_matrix=true,tol=tol).u[end].x
        return u₀-(M-Id)^-1*(u₁-u₀),norm(u₁-u₀)
    end
    function hamiltonian_gradient(system::DickeSystem,u::AbstractVector{<:Real})
        Fu₁=Float64[0,0,0,0]
        ClassicalSystems.step(system)(Fu₁,u,[ClassicalSystems.parameters(system);1.0],1.0)
        gradH=[-Fu₁[3],-Fu₁[4],Fu₁[1],Fu₁[2]]
        return gradH
    end

    function monodromy_method_step_constant_energy(system::DickeSystem,u₀::AbstractVector{<:Real},T::Real;tol::Real=1e-12)
        u₁,M=ClassicalSystems.integrate(system,u₀,T;save_everystep=false,get_fundamental_matrix=true,tol=tol).u[end].x
        gradH=hamiltonian_gradient(system,u₁)
        ξ=Float64[0,0,0,1]
        Fu₁=-[-gradH[3],-gradH[4],gradH[1],gradH[2]]
        Mat=[M-Id Fu₁;
            transpose(gradH) 0;
            transpose(ξ) 0]
        res=pinv(Mat)*[u₀-u₁;0;0]
        du=res[1:4]
        dT=res[5]
        return u₀+du,T+dT,norm(u₁-u₀)
    end

    """
    ```julia
    function monodromy_method_constant_period(system::DickeSystem,
        u₀::AbstractVector{<:Real},
        T::Real;
        maxiters::Integer=100,
        tol::Real=1e-8,
        inttol::Real=tol/10)
    ```
    Returns a [`PO`](@ref) with a given period `T` that is found by iteratively perturbating the initial condition `u₀` in
    the direction that minimizes the distance between ``u_0`` and ``u_0(T)``. The perturbation
    is obtained from Eq. (31) of Ref. [Simonovi1999](@cite), taking ``\\Lambda`` as the fundamental
    matrix evaluated at time `T`. Note that the energy of `u₀` will change. To conserve the energy and
    vary `T`, use [`monodromy_method_constant_energy`](@ref) instead.

    # Arguments
    - `system` should be an instance of [`ClassicalDicke.DickeSystem`](@ref).
    - `u₀` is a point `[Q,q,P,p]` in the phase space of the Dicke model, which is used
      as the starting point to find an orbit.
    - `T` is a positive real number, indicating the desired period.
    # Keyword arguments
    - `maxiters` is a integer which sets the maximum number of iterations. Defaults
      to `100`.
    - `tol` is a tolerance. If  ``\\|u-u(T)\\|<`` `tol`,
      ``u`` is considered a periodic condition with period `T`. The smaller the tolerance, the
      more iterations are needed to converge. Default is `1e-8`
    - `inttol` is the tolerance to be passed to [`ClassicalSystems.integrate`](@ref ClassicalSystems.integrate(::ClassicalSystems.ClassicalSystem,::AbstractVector{<:Real},::Real)). Default is
      `tol/10`.
    """
    function monodromy_method_constant_period(system::DickeSystem,
        u₀::AbstractVector{<:Real},
        T::Real;
        maxiters::Integer=100,
        tol::Real=1e-8,
        inttol::Real=tol/10)
    u₁=u₀
    for j in 1:maxiters
        u₁,Δ=monodromy_method_step_constant_period(system,u₁,T,tol=inttol)
        if Δ<tol
            break
        end
        if  j==maxiters
            @error("Maximum iterations reached")
        end
    end

    PO(system,u₁,T)
    end
    """
    ```julia
    function monodromy_method_constant_energy(system::DickeSystem,
        u₀::AbstractVector{<:Real},
        T::Real;
        maxiters::Integer=100,
        tol::Real=1e-8,
        inttol::Real=tol/10,
        correct_energy::Bool=true)
    ```
    Returns a [`PO`](@ref) with energy approximatelly (see below) the same energy as ``u₀`` obtained
    by iteratively perturbating period `T` and the initial condition `u₀`  in
    the direction that minimizes the distance between ``u_0`` and ``u_0(T)`` and is perpendicular to the
    Hamiltonian gradient. The algorithm is described in App. A, section A.1. of Ref. [Pilatowsky2021](@cite).

    # Arguments
    - `system` should be an instance of [`ClassicalDicke.DickeSystem`](@ref).
    - `u₀` is a point `[Q,q,P,p]` in the phase space of the Dicke model, which is used
      as the starting point to find an orbit.
    - `T` is a positive real number, which is used as the starting point to find an orbit.
    # Keyword arguments
    - `maxiters` is a integer which sets the maximum number of iterations. Defaults
      to `100`.
    - `tol` is a tolerance. If  ``\\|u-u(T)\\|<`` `tol`,
      ``u`` is considered a periodic condition with period `T`. The smaller the tolerance, the
      more iterations are needed to converge. Default is `1e-8`
    - `inttol` is the tolerance to be passed to [`ClassicalSystems.integrate`](@ref ClassicalSystems.integrate(::ClassicalSystems.ClassicalSystem,::AbstractVector{<:Real},::Real)). Default is
      `tol/10`.
    - `correct_energy` is a Boolean. As described in  Ref. [Pilatowsky2021](@cite), this algorithm
      only approximately conserves energy. If `correct_energy == true`, the initial condition
      is projected back to the energy shell between each iteration. This allows energy to be truly
      conserved, however, it makes the algorithm more unstable. The default is `true`.
    """
    function monodromy_method_constant_energy(system::DickeSystem,
        u₀::AbstractVector{<:Real},
        T::Real;
        maxiters::Integer=100,
        tol::Real=1e-8,
        inttol::Real=tol/10,
        correct_energy::Bool=true)
    u₁=u₀
    T₁=T
    ϵ = ClassicalDicke.hamiltonian(system)(u₀)
    for j in 1:maxiters
        if correct_energy
            sgn = ClassicalDicke.q_sign(system,u₁)
            u₁=Point(system,Q=u₁[1],P=u₁[3],p=u₁[4],sgn=sgn,ϵ=ϵ)
        end
        u₁,T₁,Δ=monodromy_method_step_constant_energy(system,u₁,T₁;tol=inttol)
        if T₁<0
            error("Convergence error")
        end
        if Δ<tol
            break
        end
        if j==maxiters
            error("Maximum iterations")
        end
    end

    PO(system,u₁,T₁)
    end

    """
    ```julia
    function find_p_zero(system::DickeSystem,
        u₀::AbstractVector{<:Real};
        tol::Real=1e-6,
        negative::Bool=false)
    ```
    Evolves `u₀` until it crosses the ``p=0`` plane from negative to positive and returns the result.
    # Arguments
    - `system` should be an instance of [`ClassicalDicke.DickeSystem`](@ref).
    - `u₀` is a point `[Q,q,P,p]` in the phase space of the Dicke model.
    # Keyword arguments
    - `tol` is the numerical tolerance.
    - If `negative` is `true`, then the crossing is from positive to negative. Default is `false`.
    """
    function find_p_zero(system::DickeSystem,
        u₀::AbstractVector{<:Real};
        tol::Real=1e-6,
        negative::Bool=false)
        nu=Float64[0,0,0,0]

        function guardar(integrator)

            nu=integrator.u

            terminate!(integrator)
            return nothing
        end
        try

            negguardar=nothing
            if negative
                negguardar,guardar=guardar,negguardar
            end
            cb=ContinuousCallback((uvar,t,integrator)-> (uvar[4]-0),guardar,negguardar,save_positions=(false,false),rootfind=true,interp_points=3,abstol=tol)
            ClassicalSystems.integrate(system,u₀,10000;callback=cb,save_everystep=false,tol=tol)
            return nu
        catch
            if abs(u₀[4])< tol  #hay un error cuando el punto ya satisface la condicion por juliaDiff
                return u₀
            end
            return NaN
        end
    end
    """
    ```julia
    function lyapunov(po::PO)
    ```
    Returns the Lyapunov exponent of the periodic orbit `po`, as given by Eq. (B3) of
    Ref. [Pilatowsky2021](@cite).
    """
     function lyapunov(po::PO)
         system=po.system
         T=po.T
         maximum(log.(abs.(
            LinearAlgebra.eigvals(ClassicalSystems.integrate(po;
                    get_fundamental_matrix=true,save_everystep=false)[end].x[2])
            ))/T)
    end
    """
    ```julia
    energy(po::PO) = ClassicalDicke.hamiltonian(po.system)(po.u)
    ```
    Returns the energy of the periodic orbit `po`.

    """
    energy(po::PO) =ClassicalDicke.hamiltonian(po.system)(po.u)
    """
    ```julia
    function follow_PO_family_from_period(po::PO;
        step::Real=0.1,
        tol::Real=1e-5,
        initalperturbation::AbstractVector{<:Real}=Float64[0,1,0,0],
        verbose::Bool =true)
    ```
    Returns a function `T -> po1` that returns a PO from the same family as `po`
    but with period `T`. This algorithm applies the function [`monodromy_method_constant_period`](@ref)
    repeatedly, increasing or decreasing the period in small perturbations to reach the target.

    Note that results are cached, so when `po` is returned by calling the function
    at some period `T`, a second call with the same period will not perform any calculation and will return
    `po` again.
    # Arguments
    - `po` should be an instance of [`PO`](@ref)
    # Keyword arguments
    - `step` is the initial size of the perturbations in period, this is decreased and increased
      dynamically (default is `0.1`).
    - `tol` is the tolerance to pass to [`monodromy_method_constant_period`](@ref) (default is `1e-5`).
    - `verbose` is a Boolean indicating whether to print the progress. Default is `true`.

    The function returned accepts the following keyword arguments:
    - `tol` overrides `tol` above.
    - `minstep` is the minimum step (which is varied dynamically). Default is `step/1000`.
    - `maxstep` is the maximum step. Default is `step`.
    - `step` overrides `step` above. Default is `maxstep`.
    - `maxiters` is the maximum number of iterations. Default is `1000`.
    """
    function follow_PO_family_from_period(po::PO;
        step::Real=0.1,
        tol::Real=1e-5,
        initalperturbation::AbstractVector{<:Real}=Float64[0,1,0,0],
        verbose::Bool =true)

        system=po.system
        u₀=po.u
        period = po.T
        us=[u₀]

        periods=Float64[period]

        function get(finalperiod::Real;tol::Real=tol, minstep::Real=step/1000,maxstep::Real=step,
            step::Real=maxstep,maxiters::Integer=1000)
            index =NaN
            if verbose
                prog=ProgressUnknown("Remaining T")
            end
            _it=0
            while true
                _it+=1
                if _it>maxiters
                    error("Maximum iterations")
                end

                index=findmin(abs.(finalperiod.-periods))[2]
                last_p=periods[index]
                falta=abs(last_p-finalperiod)

                if falta==0
                    break
                end
                step=min(step,falta)
                direction=sign(finalperiod-last_p)
                newu=us[index]
                if verbose
                    prog.desc="Current period = $(round(last_p,digits=3)), step = 10^$(round(log10(step),digits=1)), iterations ="
                    next!(prog)

                end
               try
                    ΔT=direction*step
                    o= monodromy_method_constant_period(system,newu,last_p+ΔT,inttol=min(1e-12,tol/100),tol=tol)

                    newu = o.u
                    period = o.T
                catch e
                    if step/1.2<minstep# (!converror && H(newu) in energies) ||

                            error("The algorithm is not converging. Got to ",PO(system,us[index],periods[index]))

                    end
                    step/=1.2
                    continue
                end

                step=min(maxstep,step*1.1)

                nindex=searchsortedlast(periods,period)
                insert!(us,nindex + 1,newu)
                insert!(periods,nindex + 1,period)
                index=nindex+1

            end
            if verbose
                finish!(prog)
            end
            return PO(system,us[index],periods[index])
        end
        return get
    end

    """
    ```julia
    function follow_PO_family_from_energy(po::PO;
        step::Real=0.1,
        tol::Real=1e-5,
        energytol::Real=1e-3,
        initalperturbation::AbstractVector{<:Real}=Float64[0,1,0,0],
        verbose::Bool =true,
        dontallowlessenergy::Bool=false)
    ```
    Returns a function `ϵ -> po1` that returns a PO from the same family as `po`
    but with energy `ϵ`. This algorithm is detailed in App. A.2. of Ref. [Pilatowsky2021](@cite)

    Note that results are cached, so when `po` is returned by calling the function
    at some energy `ϵ`, a second call with the same energy will not perform any calculation and will return
    `po` again.
    # Arguments
    - `po` should be an instance of [`PO`](@ref).
    # Keyword arguments
    - `step` is the initial size of the perturbations in energy, this is decreased and increased
      dynamically (default is `0.1`). If the orbits seem to suddenly jump, try to decrease this number.
    - `tol` is the tolerance to pass to [`monodromy_method_constant_energy`](@ref) (default is `1e-5`).
    - `energytol` is the precision in energy below which a PO is returned   (default is `1e-3`).
    - `initalperturbation` is de direction of the first perturbation in the form `[Q,q,P,p]` (default is `[0,1,0,0]`).
    - `verbose` is a Boolean indicating whether to print the progress. Default is `true`.
    - `dontallowlessenergy` forbids evaluation below the energy of `po`. Default is `false`.

    The function returned accepts the following keyword arguments:
    - `tol` overrides  `tol` above.
    - `energytol` overrides  `energytol` above.
    - `maxstep` is the maximum step. Default is `step`.
    - `step` overrides `step` above. Default is `maxstep`.
    - `minstep` is the minimum step (which is varied dynamically). Default is `min(step,energytol/10)`.
    - `correct_energy` is passed to [`monodromy_method_constant_energy`](@ref). Default is `true`.
    - `maxiters` is the maximum number of iterations. Default is `1000`.
    - `energywiggletolerance` is the tolerance for oscillations. Oscillating behaviour indicates
      instability, but sometimes it happens and there must be some tolerance. The default is `1e-2`.
    """
    function follow_PO_family_from_energy(po::PO;
        step::Real=0.1,
        tol::Real=1e-5,
        energytol::Real=1e-3,
        initalperturbation::AbstractVector{<:Real}=Float64[0,1,0,0],
        verbose::Bool =true,
        correct_energy = true,
        dontallowlessenergy::Bool=false)
        system=po.system
        u₀=po.u
        period = po.T
        H = ClassicalDicke.hamiltonian(system)
        us=[u₀]

        E₀=H(u₀)

        energies=Float64[H(u₀)]
        periods=Float64[period]

        cacheJ=ClassicalSystems.eye(4)
        Λ=[0 0 -1 0
           0 0 0 -1
           1 0 0 0
           0 1 0 0]

        function get(finalEnergy::Real;tol::Real=tol,
            energytol::Real=energytol,
            maxstep=step::Real,step::Real=step, minstep::Real=min(step,energytol/10),
            maxiters::Integer=1000,energywiggletolerance::Real=1e-2)
            index =NaN
            _it=0
            if verbose
                prog=ProgressUnknown("Current Energy")
            end
            if dontallowlessenergy
                finalEnergy=max(finalEnergy,E₀)
            end
            while true #abs(energies[index]-finalEnergy)>energytol
                _it+=1
                if _it>maxiters
                    error("Maximum iterations")
                end
                #index=findmin(abs.(energies.-finalEnergy))[2]
                index=max(searchsortedlast(energies,finalEnergy),1)
                last_e=energies[index]
                falta=abs(last_e-finalEnergy)
                if falta<energytol
                    break
                end
                step=min(step,falta)
                direction=sign(finalEnergy-last_e)

                newu=us[index]
                if verbose
                    prog.desc="Current energy = $(round(last_e,digits=3)), step = 10^$(round(log10(step),digits=1)), iterations ="
                    next!(prog)
                end
                converror=false
               try
                    ClassicalSystems.step(system).jac(cacheJ,newu,[ClassicalSystems.parameters(system);1.0],1.0)
                    grad=hamiltonian_gradient(system,newu)
                    hessian=Λ*cacheJ
                    if index==1
                        Δu=initalperturbation
                    else

                       Δu=hamiltonian_gradient(system,newu)
                    end

                    Δu=Δu/norm(Δu)
                    ΔE=direction*step

                    a=0.5*transpose(Δu)*hessian*Δu
                    b=dot(grad,Δu)
                    c=- ΔE
                    #Taylor series: -c=Δϵ=b*s + a s^2
                    s=(-b+sqrt(b^2-4*a*c))/(2*a)
                    if abs(H(newu + s*Δu) - (last_e + ΔE)) >abs(H(newu - s*Δu) - (last_e + ΔE))
                        s=-s
                    end
                    lastu=newu
                    E=H(newu+s*Δu)
                    periodGuess=periods[index]
                    if index>1
                        periodGuess+= (E-energies[index-1])*(periods[index]-periods[index-1])/(energies[index]-energies[index-1])
                    end

                    o= monodromy_method_constant_energy(system,newu+s*Δu,periodGuess,inttol=min(1e-12,tol/100),tol=tol,correct_energy=correct_energy)
                    newu = o.u
                    period = o.T
                    if(E-H(newu))>energywiggletolerance
                        if index>1
                            deleteat!(energies,index)
                            deleteat!(us,index)
                            deleteat!(periods,index)
                        end
                        converror=true
                        index-=1
                    end

                catch e
                    converror=true
                end

                notbetter=abs(finalEnergy-H(newu))> falta
                if  converror || H(newu) in energies || notbetter

                    if step/1.2<minstep# (!converror && H(newu) in energies) ||

                        error("The algorithm is not converging. Got to ",PO(system,us[index],periods[index]), " with energy ",H(us[index]) )

                    end
                    step/=1.2
                    continue

                end
                step=min(maxstep,step*1.1)

                nindex=searchsortedlast(energies,finalEnergy)
                insert!(us,nindex+1,newu)
                insert!(energies,nindex+1,H(newu))
                insert!(periods,nindex+1,period)
                index=nindex+1

            end
            if verbose
                finish!(prog)
            end
            return PO(system,us[index],periods[index])
        end
        return get
    end
    """
    ```julia
    function family_A(system::DickeSystem;kargs...)
    ```
    Returns a function `ϵ -> po` that returns a PO from family ``\\mathcal{A}`` of Ref. [Pilatowsky2021](@cite). This is done
    by passing the ground state fixed point with the positive normal frequency as period to [`follow_PO_family_from_energy`](@ref).

    # Arguments
    - `system` should be an instance of [`ClassicalDicke.DickeSystem`](@ref ClassicalDicke.DickeSystem).
      The system must be in the supperradiant regime.
    # Keyword arguments
    - `kargs...` are redirected to [`follow_PO_family_from_energy`](@ref)
    """
    family_A(system::DickeSystem;kargs...)=follow_PO_family_from_energy(PO(system,
        ClassicalDicke.minimum_energy_point(system,+),2*pi/ClassicalDicke.normal_frequency(system,+));
        tol=1e-8,
        correct_energy =false, #for some reason, the algorithm with correct_energy=true crashes at energy -1.4 for family_A...
        energytol=1e-3,
        kargs...)
    """
    ```julia
    function family_B(system::DickeSystem;kargs...)
    ```
    Returns a function `ϵ -> po` that returns a PO from family ``\\mathcal{B}`` of Ref. [Pilatowsky2021](@cite). This is done
    by passing the ground state fixed point with the negative normal frequency as period to [`follow_PO_family_from_energy`](@ref).

    # Arguments
    - `system` should be an instance of [`ClassicalDicke.DickeSystem`](@ref ClassicalDicke.DickeSystem).
      The system must be in the supperradiant regime.
    # Keyword arguments
    - `kargs...` are redirected to [`follow_PO_family_from_energy`](@ref)
    """
    family_B(system::DickeSystem;kargs...)=follow_PO_family_from_energy(PO(system,ClassicalDicke.minimum_energy_point(system,+),2*pi/ClassicalDicke.normal_frequency(system,-));
        tol=1e-8,energytol=1e-3,kargs...)

    """
    ```julia
    function po_coordinates(coords::AbstractArray{<:Integer},
        po::PO;tol::Real=1e-12)
    ```
    Returns an array `[Tuple(x1[coords]), Tuple(x2[coords]), ...]`, where is an array containing up
    to 4 integers from 1 to 4, and `xi` are the points in `po` obtained with numerical
    integration with tolerance `tol`.
    """
    function po_coordinates(coords::AbstractArray{<:Integer},
        po::PO;tol::Real=1e-12)
        u=integrate(po)
        us=u.u
        return [Tuple(u[i] for i in coords) for u in us];
    end
    """
    ```julia
    QP(po::PO; tol::Real=1e-12) = po_coordinates([1,3], po, tol=tol)
    ```
    Returns an array `[(Q1, P1), (Q2, P2), ...]` with the coordinates of `po`, integrated with tolerance `tol`.
    This can be directly passed to [`Plots.plot`](https://docs.juliaplots.org/latest/tutorial/)
    """
    QP(po::PO;tol::Real=1e-12)=po_coordinates([1,3],po,tol=tol)
    """
    ```julia
    qp(po::PO; tol::Real=1e-12) = po_coordinates([2,4], po, tol=tol)
    ```
    Returns an array `[(q1, p1), (q2, p2), ...]` with the coordinates of `po`, integrated with tolerance `tol`.
    This can be directly passed to [`Plots.plot`](https://docs.juliaplots.org/latest/tutorial/)
    """
    qp(po::PO;tol::Real=1e-12) = po_coordinates([2,4],po,tol=tol)
    """
    ```julia
    QPp(po::PO; tol::Real=1e-12) = po_coordinates([1,3,4], po, tol=tol)
    ```
    Returns an array `[(Q1, P1, p1), (Q2, P2, p2), ...]` with the coordinates of `po`, integrated with tolerance `tol`.
    This can be directly passed to [`Plots.plot`](https://docs.juliaplots.org/latest/tutorial/) to generate a 3D plot.
    """
    QPp(po::PO; tol::Real=1e-12)=po_coordinates([1,3,4],po,tol=tol)
    """
    ```julia
    QPq(po::PO; tol::Real=1e-12) = po_coordinates([1,3,2], po, tol=tol)
    ```
    Returns an array `[(Q1, P1, q1), (Q2, P2, q2), ...]` with the coordinates of `po`, integrated with tolerance `tol`.
    This can be directly passed to [`Plots.plot`](https://docs.juliaplots.org/latest/tutorial/) to generate a 3D plot.
    """
    QPq(po::PO; tol::Real=1e-12)=po_coordinates([1,3,2], po, tol=tol)

    """
    ```julia
    function overlap_of_tube_with_homogenous_state(po::PO{DickeBCE.QuantumDickeSystem};
                                            time_integral_tolerance::Real=1e-7,
                                            phase_space_integral_resolution::Real=0.1)
    ```
    Returns the overlap ``\\text{tr}(\\hat{\\rho}_\\epsilon \\hat{\\rho}_{\\mathcal{O}} )``
    of a tubular state ``\\hat{\\rho}_{\\mathcal{O}}`` around the periodic orbit ``\\mathcal{O}=`` `po`, (Eq. (15) of Ref. [Pilatowsky2021](@cite))
    with a totally delocalized state ``\\hat{\\rho}_\\epsilon`` (Eq. (16) of Ref. [Pilatowsky2021](@cite)).
    # Arguments
    - `po` should be an instance of [`PO`](@ref). The system passed to create `po`
      should have been a [`QuantumDickeSystem`](@ref DickeBCE.QuantumDickeSystem).
    # Keyword arguments
    - `time_integral_tolerance` is the numerical tolerance for the integral in Eq. (15) of Ref. [Pilatowsky2021](@cite). Default is `1e-7`.
    - `phase_space_integral_resolution` is the phase space resolution for the integral in Eq. (16) of Ref. [Pilatowsky2021](@cite), that is, `res` in
      [`EnergyShellProjections.energy_shell_average`](@ref). Default is `0.1`.
    """
    function overlap_of_tube_with_homogenous_state(po::PO{DickeBCE.QuantumDickeSystem};
        time_integral_tolerance::Real=1e-7,
        phase_space_integral_resolution::Real=0.1)

        system=po.system
        orbit=integrate(po,tol=time_integral_tolerance)
        ∫dtHtx(x)=average_over_PO(orbit,u-> DickeBCE.husimi_of_coherent(system,u,x))
        res=phase_space_integral_resolution
        ϵ=energy(po)
        symmP = (po == mirror_Pp(po))
        symmQP = (po == mirror_Qq(po)) && symmP
        return EnergyShellProjections.energy_shell_average(po.system,
                ϵ=ϵ,
                f=∫dtHtx,
                res=res,
                symmetricP =symmP,
                symmetricQP=symmQP)
    end
    """
    ```julia
    function scarring_measure(
        po::PO{DickeBCE.QuantumDickeSystem},
        quantum_state::AbstractVector{<:Number};
        chop::Real=1e-3,
        kargs...)
    ```
    Returns the scarring measure ``\\mathcal{P}(\\mathcal{O},\\hat{\\rho})`` as defined in Eq. (17) of Ref. [Pilatowsky2021](@cite).
    # Arguments
    - `po` should be an instance of [`PO`](@ref) representing ``\\mathcal{O}`` above.
      The system passed to create `po` should have been a [`QuantumDickeSystem`](@ref DickeBCE.QuantumDickeSystem).
    - `quantum_state` should be a vector representing the quantum state  ``\\hat{\\rho}`` in the coherent efficient basis.
    # Keyword arguments
    - `chop` is the tolerance to be passed to [`DickeBCE.husimi`](@ref) (see the documentation of [`DickeBCE.coherent_overlap`](@ref))
    - `kargs` are redirected to [`overlap_of_tube_with_homogenous_state`](@ref), which gives the denominator in Eq. (17) of Ref. [Pilatowsky2021](@cite).
    """
    function scarring_measure(
        po::PO{DickeBCE.QuantumDickeSystem},
        quantum_state::AbstractVector{<:Number};
        chop::Real=1e-3,
        kargs...)
        system=po.system
        po_and_state = average_over_PO(po,x -> DickeBCE.husimi(system,x,quantum_state;chop=chop))
        po_and_hom_state = overlap_of_tube_with_homogenous_state(po;kargs...)
        return po_and_state/po_and_hom_state
    end

end


text:
# Normed{T,f} maps UInts from 0 to 2^f-1 to the range [0.0, 1.0]
# For example, Normed{UInt8,8} == N0f8 maps 0x00 to 0.0 and 0xff to 1.0

struct Normed{T<:Unsigned,f} <: FixedPoint{T,f}
    i::T

    Normed{T, f}(i::Integer,_) where {T,f} = new{T, f}(i%T)   # for setting by raw representation
    Normed{T, f}(x) where {T,f} = convert(Normed{T,f}, x)
    Normed{T, f}(x::Normed{T,f}) where {T,f} = x
    Normed{T, f}(x::Char) where {T,f} = throw(ArgumentError("Normed cannot be constructed from a Char"))
    Normed{T, f}(x::Complex) where {T,f} = Normed{T, f}(convert(real(typeof(x)), x))
    Normed{T, f}(x::Base.TwicePrecision) where {T,f} = Normed{T, f}(convert(Float64, x))
end

typechar(::Type{X}) where {X <: Normed} = 'N'
signbits(::Type{X}) where {X <: Normed} = 0

for T in (UInt8, UInt16, UInt32, UInt64)
    for f in 0:sizeof(T)*8
        sym = Symbol(String(take!(showtype(_iotypealias, Normed{T,f}))))
        @eval begin
            const $sym = Normed{$T,$f}
            export $sym
        end
    end
end

reinterpret(::Type{Normed{T,f}}, x::T) where {T <: Unsigned,f} = Normed{T,f}(x, 0)

zero(::Type{Normed{T,f}}) where {T,f} = Normed{T,f}(zero(T),0)
function oneunit(::Type{T}) where {T <: Normed}
    T(typemax(rawtype(T)) >> (8*sizeof(T)-nbitsfrac(T)), 0)
end
one(::Type{T}) where {T <: Normed} = oneunit(T)
zero(x::Normed) = zero(typeof(x))
oneunit(x::Normed) =  one(typeof(x))
one(x::Normed) = oneunit(x)
rawone(v) = reinterpret(one(v))

# Conversions
convert(::Type{U}, x::U) where {U <: Normed} = x
convert(::Type{Normed{T1,f}}, x::Normed{T2,f}) where {T1 <: Unsigned,T2 <: Unsigned,f} = Normed{T1,f}(convert(T1, x.i), 0)
function convert(::Type{Normed{T,f}}, x::Normed{T2}) where {T <: Unsigned,T2 <: Unsigned,f}
    U = Normed{T,f}
    y = round((rawone(U)/rawone(x))*reinterpret(x))
    (0 <= y) & (y <= typemax(T)) || throw_converterror(U, x)
    reinterpret(U, _unsafe_trunc(T, y))
end
convert(::Type{U}, x::Real) where {U <: Normed} = _convert(U, rawtype(U), x)

convert(::Type{N0f16}, x::N0f8) = reinterpret(N0f16, convert(UInt16, 0x0101*reinterpret(x)))
function _convert(::Type{U}, ::Type{T}, x) where {U <: Normed,T}
    y = round(widen1(rawone(U))*x)
    (0 <= y) & (y <= typemax(T)) || throw_converterror(U, x)
    U(_unsafe_trunc(T, y), 0)
end
# Prevent overflow (https://discourse.julialang.org/t/saving-greater-than-8-bit-images/6057)
_convert(::Type{U}, ::Type{T}, x::Float16) where {U <: Normed,T} =
    _convert(U, T, Float32(x))
_convert(::Type{U}, ::Type{UInt128}, x::Float16) where {U <: Normed} =
    _convert(U, UInt128, Float32(x))
function _convert(::Type{U}, ::Type{UInt128}, x) where {U <: Normed}
    y = round(rawone(U)*x)   # for UInt128, we can't widen
    (0 <= y) & (y <= typemax(UInt128)) & (x <= Float64(typemax(U))) || throw_converterror(U, x)
    U(_unsafe_trunc(UInt128, y), 0)
end

rem(x::T, ::Type{T}) where {T <: Normed} = x
rem(x::Normed, ::Type{T}) where {T <: Normed} = reinterpret(T, _unsafe_trunc(rawtype(T), round((rawone(T)/rawone(x))*reinterpret(x))))
rem(x::Real, ::Type{T}) where {T <: Normed} = reinterpret(T, _unsafe_trunc(rawtype(T), round(rawone(T)*x)))
rem(x::Float16, ::Type{T}) where {T <: Normed} = rem(Float32(x), T)  # avoid overflow

# convert(::Type{AbstractFloat}, x::Normed) = convert(floattype(x), x)
float(x::Normed) = convert(floattype(x), x)

convert(::Type{BigFloat}, x::Normed) = reinterpret(x)*(1/BigFloat(rawone(x)))
function convert(::Type{T}, x::Normed) where {T <: AbstractFloat}
    y = reinterpret(x)*(one(rawtype(x))/convert(T, rawone(x)))
    convert(T, y)  # needed for types like Float16 which promote arithmetic to Float32
end
convert(::Type{Bool}, x::Normed) = x == zero(x) ? false : true
convert(::Type{Integer}, x::Normed) = convert(Integer, x*1.0)
convert(::Type{T}, x::Normed) where {T <: Integer} = convert(T, x*(1/oneunit(T)))
convert(::Type{Rational{Ti}}, x::Normed) where {Ti <: Integer} = convert(Ti, reinterpret(x))//convert(Ti, rawone(x))
convert(::Type{Rational}, x::Normed) = reinterpret(x)//rawone(x)

# Traits
abs(x::Normed) = x

(-)(x::T) where {T <: Normed} = T(-reinterpret(x), 0)
(~)(x::T) where {T <: Normed} = T(~reinterpret(x), 0)

+(x::Normed{T,f}, y::Normed{T,f}) where {T,f} = Normed{T,f}(convert(T, x.i+y.i),0)
-(x::Normed{T,f}, y::Normed{T,f}) where {T,f} = Normed{T,f}(convert(T, x.i-y.i),0)
*(x::T, y::T) where {T <: Normed} = convert(T,convert(floattype(T), x)*convert(floattype(T), y))
/(x::T, y::T) where {T <: Normed} = convert(T,convert(floattype(T), x)/convert(floattype(T), y))

# Comparisons
 <(x::T, y::T) where {T <: Normed} = reinterpret(x) < reinterpret(y)
<=(x::T, y::T) where {T <: Normed} = reinterpret(x) <= reinterpret(y)

# Functions
trunc(x::T) where {T <: Normed} = T(div(reinterpret(x), rawone(T))*rawone(T),0)
floor(x::T) where {T <: Normed} = trunc(x)
function round(x::Normed{T,f}) where {T,f}
    mask = convert(T, 1<<(f-1))
    y = trunc(x)
    return convert(T, reinterpret(x)-reinterpret(y)) & mask>0 ?
            Normed{T,f}(y+oneunit(Normed{T,f})) : y
end
function ceil(x::Normed{T,f}) where {T,f}
    k = 8*sizeof(T)-f
    mask = (typemax(T)<<k)>>k
    y = trunc(x)
    return convert(T, reinterpret(x)-reinterpret(y)) & (mask)>0 ?
            Normed{T,f}(y+oneunit(Normed{T,f})) : y
end

trunc(::Type{T}, x::Normed) where {T <: Integer} = convert(T, div(reinterpret(x), rawone(x)))
round(::Type{T}, x::Normed) where {T <: Integer} = round(T, reinterpret(x)/rawone(x))
floor(::Type{T}, x::Normed) where {T <: Integer} = trunc(T, x)
 ceil(::Type{T}, x::Normed) where {T <: Integer} =  ceil(T, reinterpret(x)/rawone(x))

isfinite(x::Normed) = true
isnan(x::Normed) = false
isinf(x::Normed) = false

bswap(x::Normed{UInt8,f}) where {f} = x
bswap(x::Normed)  = typeof(x)(bswap(reinterpret(x)),0)

function minmax(x::T, y::T) where {T <: Normed}
    a, b = minmax(reinterpret(x), reinterpret(y))
    T(a,0), T(b,0)
end

# Iteration
# The main subtlety here is that iterating over 0x00uf8:0xffuf8 will wrap around
# unless we iterate using a wider type
@inline start(r::StepRange{T}) where {T <: Normed} = widen1(reinterpret(r.start))
@inline next(r::StepRange{T}, i::Integer) where {T <: Normed} = (T(i,0), i+reinterpret(r.step))
@inline function done(r::StepRange{T}, i::Integer) where {T <: Normed}
    i1, i2 = reinterpret(r.start), reinterpret(r.stop)
    isempty(r) | (i < min(i1, i2)) | (i > max(i1, i2))
end

function decompose(x::Normed)
    g = gcd(reinterpret(x), rawone(x))
    div(reinterpret(x),g), 0, div(rawone(x),g)
end

# Promotions
promote_rule(::Type{T}, ::Type{Tf}) where {T <: Normed,Tf <: AbstractFloat} = promote_type(floattype(T), Tf)
promote_rule(::Type{T}, ::Type{R}) where {T <: Normed,R <: Rational} = R
function promote_rule(::Type{T}, ::Type{Ti}) where {T <: Normed,Ti <: Union{Signed, Unsigned}}
    floattype(T)
end
@generated function promote_rule(::Type{Normed{T1,f1}}, ::Type{Normed{T2,f2}}) where {T1,T2,f1,f2}
    f = max(f1, f2)  # ensure we have enough precision
    T = promote_type(T1, T2)
    # make sure we have enough integer bits
    i1, i2 = 8*sizeof(T1)-f1, 8*sizeof(T2)-f2  # number of integer bits for each
    i = 8*sizeof(T)-f
    while i < max(i1, i2)
        T = widen1(T)
        i = 8*sizeof(T)-f
    end
    :(Normed{$T,$f})
end

_unsafe_trunc(::Type{T}, x::Integer) where {T} = x % T
_unsafe_trunc(::Type{T}, x) where {T}          = unsafe_trunc(T, x)

ext: jl


text:
function _fix_compound_string(chars::Array{Char,1})

    # init tmp
    tmp_array = Array{Char,1}()

    if (isempty(chars) == true)
        return tmp_array
    end

    # get next char -
    next_char = pop!(chars)
    if (isnumeric(next_char) == true)

        push!(tmp_array, next_char)

        values = _fix_compound_string(chars)
        for value in values
            push!(tmp_array, value)
        end

    elseif (isnumeric(next_char) == false)

        # ok, so pop another char -
        next_next_char = pop!(chars)
        if (isnumeric(next_next_char) == false)

            push!(tmp_array, next_char)
            push!(tmp_array, '1')

            # push the test char back onto the stack -
            push!(chars, next_next_char)
        else
            push!(tmp_array, next_char)
            push!(tmp_array, next_next_char)
        end

        values = _fix_compound_string(chars)
        for value in values
            push!(tmp_array, value)
        end
    end

    # return -
    return tmp_array
end

function check(result::Some)::(Union{Nothing,T} where {T<:Any})

    # ok, so check, do we have an error object?
    # Yes: log the error if we have a logger, then throw the error.
    # No: return the result.value

    # Error case -
    if (isa(something(result), Exception) == true)

        # get the error object -
        error_object = result.value

        # get the error message as a String -
        error_message = sprint(showerror, error_object, catch_backtrace())
        @error(error_message)

        # throw -
        throw(result.value)
    end

    # default -
    return result.value
end

function read_file_from_path(path_to_file::String)::Array{String,1}

    # initialize -
    buffer = String[]

    # Read in the file -
    open("$(path_to_file)", "r") do file
        for line in eachline(file)
            +(buffer,line)
        end
    end

    # return -
    return buffer
end

function +(buffer::Array{String,1}, content_array::Array{String,1})
    for line in content_array
        push!(buffer, line)
    end
end

function +(buffer::Array{String,1}, line::String)
    push!(buffer, line)
end

function extract_db_file_section(file_buffer_array::Array{String,1}, single_line_section_marker::String)::Union{Nothing,String}

    # find the SECTION START on a single line -
    section_line_start = 1
    for (index, line) in enumerate(file_buffer_array)
        if (occursin(single_line_section_marker, line) == true)
            section_line_start = index
        end
    end

    # return -
    return section_line_start == 1 ? nothing : file_buffer_array[section_line_start]
end

function extract_db_file_section(file_buffer_array::Array{String,1}, start_section_marker::String,
    end_section_marker::String)::Array{String,1}

    # initialize -
    section_buffer = String[]

    # find the SECTION START AND END -
    section_line_start = 1
    section_line_end = 1
    for (index, line) in enumerate(file_buffer_array)

        if (occursin(start_section_marker, line) == true)
            section_line_start = index
        elseif (occursin(end_section_marker, line) == true || length(line) == index)
            section_line_end = index
        end
    end

    for line_index = (section_line_start+1):(section_line_end-1)
        line_item = file_buffer_array[line_index]
        push!(section_buffer, line_item)
    end

    # return -
    return section_buffer
end

function Base.:(==)(c1::KEGGCompound, c2::KEGGCompound)
    return ((c1.compound_name == c2.compound_name) &&
            (c1.compound_id == c2.compound_id) &&
            (c1.compound_formula == c2.compound_formula) &&
            (c1.compound_mw == c2.compound_mw))
end

function Base.isless(c1::KEGGCompound, c2::KEGGCompound)
    compound_id_c1 = c1.compound_id
    compound_id_c2 = c2.compound_id
    return (isless(compound_id_c1,compound_id_c2))
end

function extract_metabolite_symbols(reaction_phrase::String)

    # initialize -
    metabolite_symbol_array = Array{String,1}()

    # ok, so we need to check - do we have a <=> in this reaction?
    if (contains(reaction_phrase, "<=>") == true)

        # ok, so if we get here, then we have a top level reaction phrase -
        # split this, and go down a level -
        reaction_components = string.(split(reaction_phrase, "<=>") .|> lstrip .|> rstrip)
        for reaction_component in reaction_components

            # get the symbols -
            values = extract_metabolite_symbols(reaction_component)
            for value in values
                push!(metabolite_symbol_array, value)
            end
        end
    elseif (contains(reaction_phrase, "+") == true)

        # ok, so if we get here, then we have a left or right reaction phrase
        # split this, and go down a level -
        reaction_components = string.(split(reaction_phrase, "+") .|> lstrip .|> rstrip)
        for reaction_component in reaction_components

            # get the symbols -
            values = extract_metabolite_symbols(reaction_component)
            for value in values
                push!(metabolite_symbol_array, value)
            end
        end
    elseif (isnumeric(reaction_phrase[1]) == true &&
            isequal(reaction_phrase[2], ' ') == true)

        # split -
        reaction_component = string.(split(reaction_phrase, ' ')[2])
        values = extract_metabolite_symbols(reaction_component)
        for value in values
            push!(metabolite_symbol_array, value)
        end
    else

        # finally -> we just have a symbol -
        return [reaction_phrase |> lstrip |> rstrip]
    end

    return metabolite_symbol_array
end

function extract_stoichiometric_dictionary(reaction_phrase::String,
    direction::Symbol = :reactants)

    # initialize -
    stoichiometric_dictionary = Dict{String,Any}()

    # ok, so we need to check - do we have a <=> in this reaction?
    if (contains(reaction_phrase, "<=>") == true)

        # ok, so if we get here, then we have a top level reaction phrase -
        # split this, and go down a level -
        reaction_components = string.(split(reaction_phrase, "<=>") .|> lstrip .|> rstrip)
        for (index, reaction_component) in enumerate(reaction_components)

            # setup the direction -
            index == 1 ? direction = :reactants : direction = :products

            # get the symbols -
            value_dictionary = extract_stoichiometric_dictionary(reaction_component, direction)
            for (key, value) in value_dictionary
                stoichiometric_dictionary[key] = value
            end
        end
    elseif (contains(reaction_phrase, "+") == true)

        # ok, so if we get here, then we have a left or right reaction phrase
        # split this, and go down a level -
        reaction_components = string.(split(reaction_phrase, "+") .|> lstrip .|> rstrip)
        for reaction_component in reaction_components

            # get the symbols -
            value_dictionary = extract_stoichiometric_dictionary(reaction_component, direction)
            for (key, value) in value_dictionary
                stoichiometric_dictionary[key] = value
            end
        end
    elseif (length(reaction_phrase) >= 2 &&
            isnumeric(reaction_phrase[1]) == true &&
            isequal(reaction_phrase[2], ' ') == true)

        # what is my multiple?
        direction == :reactants ? α = -1 : α = 1

        # split -
        st_value = string.(split(reaction_phrase, ' ')[1])
        metabolite_key = string.(split(reaction_phrase, ' ')[2])

        # if we get here, then we have a coefficient -
        stoichiometric_coefficient = α * parse(Float64, st_value)

        # add to dictionary -
        stoichiometric_dictionary[metabolite_key] = stoichiometric_coefficient

    else

        # when we get here, we just have a single bare metabolite -

        # what is my multiple?
        direction == :reactants ? α = -1 : α = 1

        # if we get here, then we have a coefficient -
        stoichiometric_coefficient = α * 1.0
        metabolite_key = reaction_phrase

        # finally -> we just have a symbol -
        stoichiometric_dictionary[metabolite_key] = stoichiometric_coefficient
    end

    # remove spaces ... not sure where this is coming from ... but this hack should fix ...
    delete!(stoichiometric_dictionary, "")

    # return -
    return stoichiometric_dictionary
end

function extract_atom_dictionary(formula::String)

    # test -
    atom_dictionary = Dict()
    local_array = Array{Char,1}()
    element_key_array = Array{Char,1}()

    # turn string into char array -
    formula_char_array = collect(formula)

    # add extra 1 if last char is a letter -
    if (isnumeric(last(formula_char_array)) == false)
        push!(formula_char_array, '1')
    end

    # clean the array -
    C_tmp = _fix_compound_string(reverse(formula_char_array))

    @show C_tmp

    # read from the bottom -
    reverse!(C_tmp)

    # how many chars do we have?
    while (isempty(C_tmp) == false)

        # clean out the array from the last pass -
        empty!(local_array)
        empty!(element_key_array)

        # grab the next value -
        next_value = pop!(C_tmp)
        if (isnumeric(next_value) == false)

            # we have an element -> read until I hit another element -
            is_ok_to_loop = true
            while (is_ok_to_loop)

                if (isempty(C_tmp) == true)
                    break
                end

                read_one_ahead = pop!(C_tmp)
                if (isnumeric(read_one_ahead) == true)
                    push!(local_array, read_one_ahead)
                else

                    # ok: so if we get here - then we read the next char, but it was a
                    # letter (element) - so we need to push it back on the stack ...
                    push!(C_tmp, read_one_ahead)

                    # shutodown -
                    is_ok_to_loop = false
                end
            end
        end

        # we need to turn local array into a string -
        buffer = ""
        [buffer *= string(x) for x in local_array]
        if (isempty(buffer) == false)
            atom_dictionary[string(next_value)] = parse(Int64, buffer)
        end
    end

    # return -
    return atom_dictionary
end


text:
int [] doubleArray (int [] a){
  int [] res = new int [a . length];
  int i = 0 ;
  for (int n : a){
    res [i] = 2 * n ;
    i ++ ;
  }
  return res ;
}

void shiftLeft (int [] a){
  int x = a [0];
  int i = 0 ;
  while (i < a.length - 1){
    a [i] = a [i + 1];
    i ++ ;
  }
  a[a.length - 1]= x ;
  return;
 }

int scalProd(int[] a, int[] b) {
  int res = 0;
  int i = 0;
  while (i < a.length) {
    res = res + a[i] * b[i];
    i++;
  }
  return res;
}

int main () {
  int [] a = new int [5];
  int i = 0 ;
  while (i < a.length){
    a [i]= i ;
    i ++ ;
    }
  shiftLeft (a);
  int [] b = doubleArray (a);
  for (int x : a)printInt (x);
  for (int x : b)printInt (x);
  printInt(scalProd(a,b));
  return 0 ;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
### struct for defining a mixed-model potential
mutable struct MixedModel <: AbstractMixedModel
  dim::Int64
  potentials::Vector{T} where T <: AbstractPotential
  function MixedModel(dim::Int64,potentials::Vector{T}) where T <: AbstractPotential
    new(dim,potentials)
  end
end

function getdimensionality(mm::MixedModel)
  return mm.dim
end

"""
Force functions with return value
"""
function force(system::AbstractSystem, mm::MixedModel)
  fout::Array{Float64} = zeros(Float64,size(system.x))
  for pot in mm.potentials
    fout .+= force(system,pot)
  end
  return fout
end

"""
Force functions with in-place changing
"""
function force!(system::AbstractSystem, mm::MixedModel)
  for pot in mm.potentials
    force!(system,pot)
  end
end

"""
Force function with in-place changing of a particular model
"""
function force!(ind::Int64,system::AbstractSystem,mm::MixedModel)
  force!(system,mm.potentials[ind])
end

function addforce!(dt::Float64,system::AbstractSystem,mm::MixedModel)
  for pot in mm.potentials
    addforce!(dt,system,pot)
  end
end

function addforce!(dt::Float64,sys::ActiveBrownianSystem,mm::MixedModel)
  for pot in mm.potentials
    if pot.dim == sys.dim
      sys.x += dt .* pot.f
    else
      @. @views sys.x[1:sys.dim-1] .+= dt .* pot.f[1:sys.dim-1]
    end
  end
  @. @views sys.x[1:sys.dim-1] .+= (dt*sys.v0) .* [cos(sys.x[end]),sin(sys.x[end])]
end

"""
Compute the gradient of the parameters for a particular model
"""
function gradient!(ind::Int64,grad::Array{Float64},system::AbstractSystem,mm::MixedModel)
  gradient!(grad,system,mm.potentials[ind])
end


text:
# --------------------------------------------------------------------------
# ACE1.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

using BenchmarkTools
globalsuite = BenchmarkGroup()

include("profile_ylm.jl")

# judge(minimum(results), minimum(results1)) |> display


text:
cd(@__DIR__)
using DocStringExtensions
using Statistics
using StatsBase
using LinearAlgebra
using DataFrames
using FastTransforms
using Distributions
using Random
# Base.:|>(args::Tuple, f) = f(args...)

C(r, a = 2.0) = (1 + r^2)^(-a / 2)

periodic_dist(x, L) = x ≤ L / 2 ? x : L - x

"""
    generate_corr_matrix(dims::Tuple, C::Function; T::Type = Float64)

Calculate the function `C` for each distance `r` and returns an array where each C(r) is stored.
The distnace `r` is measured from upper left corner and implies periodic boundary conditions.

# Arguments:
- `dims::Tuple`: dimensions of the desired lattice
- `C::Function`: desired correlation function `C(r)`
- `T::Type = Float64`: type of the retured array
"""
function generate_corr_matrix(dims::Tuple, C::Function; T::Type = Float64)
	Λ = Array{T}(undef, dims)
	fill_periodic_distances!(Λ, C)
	return Λ
end

"""
    fill_periodic_distances!(Λ::AbstractArray, C::Function)

Helper function for the `generate_corr_matrix` function.
The type of `Λ` is known on call and therefore can this inner function works faster.

# Arguments:
- `Λ::AbstractArray`: array to save the `C(r)` to
- `C::Function`: desired correlation function `C(r)`
"""
function fill_periodic_distances!(Λ::AbstractArray, C::Function)
	for i ∈ CartesianIndices(Λ)
		r = norm(periodic_dist.(Tuple(i) .- 1, size(Λ)))
		Λ[i] = C(r)
	end
end

"""
    generate_continuous_disorder(dims::Tuple, S::AbstractArray)

Produces an array with dimensions `dims` with correlated disorder where the correlation follows the function `C(r)`

# Arguments:
- `dims::Tuple`: dimensions of the desired lattice
- `C::Function`: desired correlation function `C(r)`
"""
function generate_continuous_disorder(dims::Tuple, S::AbstractArray)
	φq = sqrt.(max.(S, 0)) .* rand.(Normal(0, √(2 * prod(dims)))) # 2V = 2prod(dims) is used because the random numbers are produced only in real space
	φ = real.(ifft(φq))
	return φ
end

"""
    generate_S(dims::Tuple, C::Function)

Produces a spectral density array with dimensions `dims` of the array with `C(r)` entries

# Arguments:
- `dims::Tuple`: dimensions of the desired lattice
- `C::Function`: desired correlation function `C(r)`
"""
generate_S(dims::Tuple, C::Function) = real.(fft(generate_corr_matrix(dims, C)))

"""
    transform_to_discrete_disorder(Λ::AbstractArray, p::Number; T = Int8(1), F = Int8(0))

Produces a truncated array similar to `Λ` but only with `T` and `F` values.
The concentration of `T` is given through `p`.

# Arguments:
- `Λ::AbstractArray`: expect an array with continuous correlated numbers
- `p::Number`: concentration of `T` values
- `T = Int8(1)`: values with concentration `p` to fill into the returned array
- `F = Int8(0)`: values with concentration `1 - p` to fill into the returned array
"""
transform_to_discrete_disorder(Λ::AbstractArray, p::Number; T = Int8(1), F = Int8(0)) = [i ≤ quantile(Normal(0, 1), p) ? T : F for i ∈ Λ]

######################### TEST #########################
Random.seed!(1)
L1 = 1000
L2 = 2000
a = Inf
p = 0.8

S = generate_S((L1, L2), r -> C(r, a))
Λc = generate_continuous_disorder((L1, L2), S)
Λd = transform_to_discrete_disorder(Λc, p)
@show mean(Λc)
@show std(Λc)
@show check_p = counts(Λd, 1:1)[1] / prod(size(Λd))

using PyPlot
clf()
imshow(Λc; interpolation = "none")
savefig("$(L1)x$(L2)_a$(a)_p$(p)_continuous_disorder.pdf")

clf()
imshow(Λd, interpolation = "none")
savefig("$(L1)x$(L2)_a$(a)_p$(p)_discrete_disorder.pdf")

# ######################### TEST SPEED #########################
# using BenchmarkTools

# @btime generate_continuous_disorder((L1, L2), r -> C(r, a))
# @btime transform_to_discrete_disorder(Λc, p)

A = -1.38

Int(round(abs(A))) # reads: perform Int of round of abs of A
round(abs(A)) |> Int
abs(A) |> round |> Int
A |> abs |> round |> Int # reads: take A, do abs, round it, make it an Int

### dot version

A = [0.1, -1.3, 2.5]

Int.(round.(abs.(A)))
round.(abs.(A)) .|> Int
abs.(A) .|> round .|> Int
A .|> abs .|> round .|> Int

### absurd

(1 + 3)^(1 / 2)
(1 + 3) |> x -> ^(x, 1 / 2)

############### in-place operations and |>

A = rand(10000)

B = round.(abs.(A))
A .= round.(abs.(A)) # in-place
A .= A .|> abs .|> round # in-place

max_stars_repo_name: zierenberg/CorrelatedNoiseCollection

max_issues_repo_name: zierenberg/CorrelatedNoiseCollection

max_forks_repo_name: zierenberg/CorrelatedNoiseCollection


text:
using Dates

mutable struct FittedBondDiscountCurve{C <: BusinessCalendar, B <: BondHelper, DC <: DayCount, F <: FittingMethod} <: Curve
  lazyMixin::LazyMixin
  settlementDays::Int
  referenceDate::Date
  calendar::C
  bonds::Vector{B}
  dc::DC
  fittingMethod::F
  accuracy::Float64
  maxEvaluations::Int
  simplexLambda::Float64

  FittedBondDiscountCurve{C, B, DC, F}(settlementDays::Int,
                          referenceDate::Date,
                          calendar::C,
                          bonds::Vector{B},
                          dc::DC,
                          fittingMethod::F,
                          accuracy::Float64,
                          maxEvaluations::Int,
                          simplexLambda::Float64) where {C, B, DC, F} =

                          (x = new(LazyMixin(), settlementDays, referenceDate, calendar, bonds, dc, fittingMethod, accuracy, maxEvaluations, simplexLambda);
                          x.fittingMethod.commons.costFunction.curve = x)

  # FittedBondDiscountCurve(settlementDays::Integer, referenceDate::Date, calendar::BusinessCalendar, bonds::Vector{B}, dc::DayCount, fittingMethod::FittingMethod, accuracy::Float64=1e-10,
  #                                     maxEvaluations::Integer=10000, simplexLambda::Float64=1.0) =
  #                         new(settlementDays, referenceDate, calendar, bonds, dc, fittingMethod, accuracy, maxEvaluations, simplexLambda)
    #n = length(bonds)
  #   println("hi")
  #   this = new(settlementDays, referenceDate, calendar, bonds, dc, fittingMethod, accuracy, maxEvaluations, simplexLambda)
  #   println("also hi")
  #   this.fittingMethod.costFunction.curve = this
  #
  #   return this
  # end
end

FittedBondDiscountCurve(settlementDays::Int,
                        referenceDate::Date,
                        calendar::C,
                        bonds::Vector{B},
                        dc::DC,
                        fittingMethod::F,
                        accuracy::Float64=1e-10,
                        maxEvaluations::Int=10000,
                        simplexLambda::Float64=1.0) where {C <: BusinessCalendar, B <: BondHelper, DC <: DayCount, F <: FittingMethod} =
                            FittedBondDiscountCurve{C, B, DC, F}(settlementDays, referenceDate, calendar, bonds, dc, fittingMethod, accuracy, maxEvaluations, simplexLambda)

### Fitted curve methods ###
discount_impl(curve::FittedBondDiscountCurve, t::Float64) = discount_function(curve.fittingMethod, curve.fittingMethod.commons.solution, t)

function initialize!(curve::FittedBondDiscountCurve)
 # yield conventions
 dc = curve.dc
 yield_comp = CompoundedCompounding()
 freq = Annual()

 n = length(curve.bonds)
 cost_f = curve.fittingMethod.commons.costFunction

 squared_sum = 0.0
 @inbounds for i in eachindex(curve.bonds)
   bond = curve.bonds[i].bond
   leg = bond.cashflows
   clean_price = bond.faceAmount
   bond_settlement = get_settlement_date(bond)

   # get the ytm of the bond
   ytm = yield(bond, clean_price, dc, yield_comp, freq, bond_settlement)
   dur = duration(bond, ytm, dc, yield_comp, freq, ModifiedDuration(), bond_settlement)

   curve.fittingMethod.commons.weights[i] = 1.0 / dur
   squared_sum += curve.fittingMethod.commons.weights[i] * curve.fittingMethod.commons.weights[i]

   cf = bond.cashflows
   for k = 1:length(cf.coupons) # for redemption
     if !has_occurred(cf.coupons[k], bond_settlement)
       cost_f.firstCashFlow[i] = k
       break
     end
   end
 end

 divide_array_by_self!(curve.fittingMethod.commons.weights, sqrt(squared_sum))

 return curve
end

function perform_calculations!(curve::FittedBondDiscountCurve)
 cost_f = curve.fittingMethod.commons.costFunction
 constraint = NoConstraint()

 x = zeros(curve.fittingMethod.size)
 # x = Vector{DD}(curve.fittingMethod.size)

 if length(curve.fittingMethod.commons.guessSolution) > 0
   x = curve.fittingMethod.commons.guessSolution
 end

 simplex = Simplex(curve.simplexLambda)
 problem = Problem(cost_f, constraint, x)

 max_stationary_state_iterations = 100
 root_epsilon = curve.accuracy
 function_epsilon = curve.accuracy
 gradient_norm_epsilon = curve.accuracy

 end_criteria = EndCriteria(curve.maxEvaluations, max_stationary_state_iterations, root_epsilon, function_epsilon, gradient_norm_epsilon)

 minimize!(simplex, problem, end_criteria)
 curve.fittingMethod.commons.solution = problem.currentValue

 number_of_iterations = problem.functionEvaluation
 cost_value = problem.functionValue

 curve.fittingMethod.commons.guessSolution = curve.fittingMethod.commons.solution
 curve.fittingMethod.commons.numberOfIterations = number_of_iterations
 curve.fittingMethod.commons.minimumCostValue = cost_value

 # we have calculated

 return curve
end


text:
# Test qubit Clifford algebra relations
@testset "Qubit Clifford relations" begin
    @test comm(q('x'),q('y')) == 2*im*q('z')
    @test comm(q('y'),q('z')) == 2*im*q('x')
    @test comm(q('z'),q('x')) == 2*im*q('y')
    @test acomm(q('x'),q('y')) == spzeros(QComp,q.dim,q.dim)
    @test acomm(q('y'),q('z')) == spzeros(QComp,q.dim,q.dim)
    @test acomm(q('z'),q('x')) == spzeros(QComp,q.dim,q.dim)
    @test acomm(q('x'),q('x')) == 2*q('i')
    @test acomm(q('y'),q('y')) == 2*q('i')
    @test acomm(q('z'),q('z')) == 2*q('i')
    @test acomm(q('x'),q('y')) == spzeros(QComp,q.dim,q.dim)
    @test acomm(q('y'),q('z')) == spzeros(QComp,q.dim,q.dim)
    @test acomm(q('z'),q('x')) == spzeros(QComp,q.dim,q.dim)
end

# Test superpositions and bra
xplus = (q[1] + q[2])/sqrt(2)
xplus2 = (q(0) + q(1))/sqrt(2)
@test xplus == xplus2
@test bra(xplus)(xplus) ≈ dot(xplus,xplus)
yplus = (q[1] + im*q[2])/sqrt(2)

# Test transition operator
xtoy = transition(yplus, xplus)
@test tr(xtoy) ≈ dot(xplus, yplus)

# Test expectations and weak values
@test weakvaluevec(xplus, yplus, q('z')) == QComp(-im)
xplusp = projector(xplus)
yplusp = projector(yplus)
@test weakvalue(xplusp, yplusp, q('z')) == QComp(-im)


text:
# readme_example.jl
using Kyulacs: Observable, QuantumCircuit, QuantumState
using Kyulacs.Gate: CNOT, Y, merge

state = QuantumState(3)
seed = 0  # set random seed
state.set_Haar_random_state(seed)

circuit = QuantumCircuit(3)
circuit.add_X_gate(0)
merged_gate = merge(CNOT(0, 1), Y(1))
circuit.add_gate(merged_gate)
circuit.add_RX_gate(1, 0.5)
circuit.update_quantum_state(state)

observable = Observable(3)
observable.add_operator(2.0, "X 2 Y 1 Z 0")
observable.add_operator(-3.0, "Z 2")
value = observable.get_expectation_value(state)
# 0.2835596510287872
println(value)


text:
module OperatorFlux

using Flux
using FFTW
using Tullio
using ChainRulesCore
using GaussQuadrature
using LinearAlgebra

import FFTW: r2r

export FourierTransform,
    ChebyshevTransform,
    LegendreTransform,
    SpectralConv,
    SpectralKernelOperator,
    SpectralCovariance,
    forward,
    inverse,
    truncate_modes,
    pad_modes

"""
    Supertype for all spectral transforms.
"""
abstract type AbstractTransform end

"""
    Supertype for all operators.
"""
abstract type AbstractOperator end

"""
    Supertype for all spectral operators.
"""
abstract type AbstractSpectralOperator end

"""
    forward(trafo, x)

Forward transform of `trafo` applied to `x`.
"""
forward(trafo::AbstractTransform, _...) =
    error("forward not implemented for given trafo")

"""
    inverse(trafo, x)

Inverse transform of `trafo` applied to `x`.
"""
inverse(trafo::AbstractTransform, _...) =
    error("inverse not implemented for given trafo")

"""
    truncate_modes(trafo, c)

Spectral mode truncation for transform of `trafo` applied to `c`.
"""
truncate_modes(trafo::AbstractTransform, _...) =
    error("truncate_modes not implemented for given trafo")

"""
    pad_modes(trafo, c, size_pad)

Spectral mode padding of size `size_pad` for transform of `trafo` applied to `c`.
"""
pad_modes(trafo::AbstractTransform, _...) =
    error("pad_modes not implemented for given trafo")

include("utils.jl")
include("transform_fourier.jl")
include("transform_chebyshev.jl")
include("transform_legendre.jl")
include("chainrules.jl")
include("operators.jl")

end # module


text:
module NamasteJL

#
# TODO:
# + _get_namaste, _make_namaste, _set_namaste might be interesting as macros
# + Get should support taking an option arg to return specific value (e.g. type, who, what, where when
# + _make_namaste should probably be encode namaste, and we need an decode namaste
#

#
# Namaste.jl is a port of github.com/mjgiarlo/namaste from Python
# to Julia.
#

export Encode, Decode, GetType, GetTypes, DirType, Who, What, When, Where

function normalizeTagName(tag::AbstractString)
    names = Dict(
        "type" => "0",
        "who" => "1",
        "what" => "2",
        "when" => "3",
        "where" => "4",
    )
    default = tag[0:end]
    return get(names(lowercase(tag)), default)
end

function charEncode(value::AbstractString)
    #FIXME: Decide on encoding/decoding scheme, probably use the one
    # from Pairtree spec.
    return value
end

function charDecode(value::AbstractString)
    #FIXME: Decide on encoding/decoding scheme, probably use the one
    # from Pairtree spec.
    return value
end

function Encode(tag::AbstractString, value::AbstractString)
    if tag == ""
        return f"$(charEncode(value))"
    end
    f"$(normalizeTagName(tag))=$(charEncode(value))"
end

function Decode(value::AbstractString)
    #FIXME: Need to add support for x_ fields
    if startswith(value, "0=") || startswith(value, "1=") || startswith(value, "2=") || startswith(value, "3=") || startswith(value, "4=")
        return f"$(value[1:2])$(charDecode(value[3:end]))"
    end
    f"$(charDecode(value))"
end

function _make_namaste(tag::Int, value::AbstractString)
    Encode("$tag", value)
end

function _get_namaste(d::AbstractString, tag::Int)
    tags = String[]
	d_info = readdir(d)
	namaste = String[]
	for f in filter(x -> startswith(x, "$tag="), readdir(d))
		push!(namaste, f)
	end
    if length(namaste) == 0
        return "", "namaste missing"
    end
    "namaste: $(join(namaste, ", "))", ""
end

function _set_namaste(d::AbstractString, tag::Int, value::AbstractString)
	if value == ""
		return "", "$value is empty"
    end
    if isdir(d) == false
		return "", "$d is not a directory"
	end
    namaste = joinpath(d, _make_namaste(tag, value))
    open(namaste, "w") do f
        write(f, "$value\n")
    end
    namaste, ""
end

function DirType(d::AbstractString, value::AbstractString; verbose::Bool = true)
    (namaste, err) = _set_namaste(d, 0, value)
    if err != ""
        return "", err
    end
    if verbose
        println(namaste)
    end
    namaste, ""
end

function Who(d::AbstractString, value::AbstractString; verbose::Bool = true)
    (namaste, err) = _set_namaste(d, 1, value)
    if err != ""
        return "", err
    end
    if verbose
        println(namaste)
    end
    namater, ""
end

function What(d::AbstractString, value::AbstractString; verbose::Bool = true)
    (namaste, err) = _set_namaste(d, 2, value)
    if err != ""
        return "", err
    end
    if verbose
        println(namaste)
    end
    namaste, ""
end

function When(d::AbstractString, value::AbstractString; verbose::Bool = true)
    (namaste, err) = _set_namaste(d, 3, value)
    if err != ""
        return "", err
    end
    if verbose
        println(namaste)
    end
    namaste, ""
end

function Where(d::AbstractString, value::AbstractString; verbose::Bool = true)
    (namaste, err) = _set_namaste(d, 4, value)
    if err != ""
        return "", err
    end
    if verbose
        println(namaste)
    end
    namaste, ""
end

function GetType(d::AbstractString; verbose::Bool = true)
    tags = String[]
    if d == ""
        d = "."
    end
    if isdir(d) == false
        return "", "$d is not a directory"
    end
    for t in [0,1,2,3,4]
        tags = vcat(tags, _get_namaste(d, t))
    end
    if verbose
        println("namastes: $(join(tags, ", "))")
    end
    tags, ""
end

function GetTypes(d::AbstractString; verbose::Bool = true)
    if d == ""
        d = "."
    end
    if isdir(d) == false
        return "", "$d is not a directory"
    end
    types = Dict()
    (tags, err) = _get_namaste(d, 0)
    if err != ""
        return "", err
    end
    if length(tags) > 0
        for tag in tags
            (name, major, minor) = ("", "", "")
            if contains(tag, "_")
                (name, version) = split(tag, "_", limit = 2)
                (major, minor) = split(version, ".", limit = 2)
                types[name] = Dict("name" => name, "major" => major, "minor" => minor)
                if verbose
                    println("namaste - directory type $name - version $major $minor")
                end
            end
        end
    end
    types, ""
end

end # module

max_stars_repo_name: rsdoiel/NamasteJL

max_issues_repo_name: rsdoiel/NamasteJL

max_forks_repo_name: rsdoiel/NamasteJL


text:
using BioMedQuery.Entrez
using SQLite
using BioMedQuery.DBUtils

db_path = "./test_db.sqlite"

config = Dict(:db_path=>db_path, :overwrite=>true)
db = Entrez.DB.init_pubmed_db_sqlite(config)

#check collection of tables
tables_query = BioMedQuery.DBUtils.select_all_tables(db)
tables = ["article","author","author2article","mesh_descriptor","mesh_heading","mesh_qualifier", "sqlite_sequence"]
@test sort(tables) == sort(tables_query)

#check minimum insert
BioMedQuery.DBUtils.insert_row!(db, "article", Dict(:pmid => 1234,
:title=>"Test Article",
:pubYear=>nothing))

sel = BioMedQuery.DBUtils.db_select(db, ["pmid"], "article", Dict(:title=>"Test Article", :pmid=>1234))

@test length(sel[1]) == 1
@test sel[1][1].value == 1234

if isfile(db_path)
    rm(db_path)
end


text:
using Documenter
using PlotlyJS
using PlotlyBase

PlotlyJS.set_default_renderer(PlotlyJS.DOCS)

const THIS_DIR = dirname(@__FILE__)

# used in examples
using Distributions, HTTP, DataFrames, RDatasets, Colors, CSV, JSON
using Random, Dates, LinearAlgebra, DelimitedFiles

# to override display_dict below
# import Documenter.Utilities: display_dict, limitstringmime
# using Base64:stringmime

# function display_dict(p::PlotlyBase.Plot)
#     out = Dict{MIME,Any}()
#     # Always generate text/plain
#     out[MIME"text/plain"()] = limitstringmime(MIME"text/plain"(), p)
#     svg_m = MIME"image/svg+xml"()
#     out[svg_m] = stringmime(svg_m, p)
#     out
# end

# display_dict(p::PlotlyJS.SyncPlot) = display_dict(p.plot)

## handle examples
# Walk through each example in a file and get the markdown from `single_example`
function single_example_file(filename::String)
    base_fn = split(filename, ".")[1]
    start_example = "```@example $(base_fn)"
    end_example = "```"
    # Open a file to write to
    open(joinpath(THIS_DIR, "src", "examples", "$(base_fn).md"), "w") do outfile

        write_example(ex) = println(outfile, start_example, "\n", ex, "\n", end_example, "\n")

        fn_h1 = titlecase(replace(base_fn, "_" => " "))
        println(outfile, "# $(fn_h1)\n")

        # Read lines from a files
        fulltext = open(
            f -> read(f, String),
            joinpath(THIS_DIR, "..", "examples", filename),
            "r"
        )
        all_lines = split(fulltext, "\n")
        l = 1
        regex = r"^function ([^_].+?)\("
        regex_end = r"^end$"

        # find preamble
        if base_fn == "subplots"  # special case
            preamble = "using PlotlyJS, Dates\ninclude(\"../../../examples/line_scatter.jl\")"
            write_example(preamble)
        else
            first_line = findfirst(x -> match(regex, x) !== nothing, all_lines)
            if first_line !== nothing
                preamble = strip(join(all_lines[1:first_line - 1], "\n"))
                write_example(preamble)
            end
        end

        while true
            # Find next function name (break if none)
            l = findnext(x -> match(regex, x) !== nothing, all_lines, l + 1)
            if l == 0 || l === nothing
                break
            end
            # find corresponding end for this function
            end_l = findnext(x -> match(regex_end, x) !== nothing, all_lines, l + 1)

            # Pull out function text
            func_block = join(all_lines[l:end_l], "\n")
            fun_name = match(regex, all_lines[l])[1]
            # println("adding $fun_name")
            an_ex = string(func_block, "\n", fun_name, "()")
            write_example(an_ex)
            l = end_l
        end
    end  # do outfile

    return nothing
end

function handle_examples()
    @info "Updating examples!"
    # Read all file names in
    if length(ARGS) == 0
        all_file_names = readdir(joinpath(THIS_DIR, "..", "examples"))
    else
        all_file_names = [endswith(i, ".jl") ? i : "$(i).jl" for i in ARGS]
    end
    all_julia_files = filter(x -> endswith(x, ".jl"), all_file_names)

    foreach(single_example_file, all_julia_files)
end

handle_examples()

makedocs(
    sitename="PlotlyJS",
    format=Documenter.HTML(
        assets=[
           "include_plotlyjs.js"
        ]
    ),
    modules=[PlotlyJS, PlotlyBase],
    linkcheck=true,
    pages=[
        "Home" => "index.md",
        "User Guide" => [
            "Preliminaries" => "basics.md",
            "Building Blocks" => "building_traces_layouts.md",
            "Putting it together" => "syncplots.md",
            "Working with plots" => "manipulating_plots.md",
            "Contributing" => "contributing.md",
        ],
        "Examples" => [
            "3d" => "examples/3d.md",
            "Area" => "examples/area.md",
            "Bar" => "examples/bar.md",
            "Box Plots" => "examples/box_plots.md",
            "Contour" => "examples/contour.md",
            "Financial Charts" => "examples/finance.md",
            "Heatmaps" => "examples/heatmaps.md",
            "Histograms" => "examples/histograms.md",
            "Line and Scatter" => "examples/line_scatter.md",
            "Maps" => "examples/maps.md",
            "Shapes" => "examples/shapes.md",
            "Subplots" => "examples/subplots.md",
            "Tables" => "examples/tables.md",
            "Ternary" => "examples/ternary.md",
            "Time Series" => "examples/time_series.md",
            "Violin" => "examples/violin.md",
        ],
        "API Docs" => "api.md"
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#= deploydocs(
    repo = "<repository url>"
) =#

deploydocs(
    repo="github.com/JuliaPlots/PlotlyJS.jl.git",
)


text:
#=
primal dual point
=#

mutable struct Point{T <: Real}
    vec::Vector{T}

    x::SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int}}, true}
    y::SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int}}, true}
    z::SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int}}, true}
    tau::SubArray{T, 0, Vector{T}, Tuple{Int}, true}
    s::SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int}}, true}
    kap::SubArray{T, 0, Vector{T}, Tuple{Int}, true}

    ztsk::SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int}}, true}
    z_views::Vector{SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int}}, true}}
    s_views::Vector{SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int}}, true}}
    dual_views::Vector{SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int}}, true}}
    primal_views::Vector{SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int}}, true}}

    Point{T}() where {T <: Real} = new{T}()
end

function Point(
    model::Models.Model{T};
    ztsk_only::Bool = false,
    ) where {T <: Real}
    point = Point{T}()
    (n, p, q) = (model.n, model.p, model.q)
    tau_idx = n + p + q + 1
    dim = tau_idx + q + 1
    vec = point.vec = zeros(T, tau_idx + q + 1)

    @views point.ztsk = vec[(n + p + 1):end]
    ztsk_only && return point

    point.tau = view(vec, tau_idx)
    point.kap = view(vec, dim)
    @views begin
        point.z = vec[n + p .+ (1:q)]
        point.s = vec[tau_idx .+ (1:q)]
        point.x = vec[1:n]
        point.y = vec[n .+ (1:p)]
    end

    point.z_views = [view(point.z, idxs) for idxs in model.cone_idxs]
    point.s_views = [view(point.s, idxs) for idxs in model.cone_idxs]
    point.dual_views = [Cones.use_dual_barrier(cone_k) ? point.s_views[k] :
        point.z_views[k] for (k, cone_k) in enumerate(model.cones)]
    point.primal_views = [Cones.use_dual_barrier(cone_k) ? point.z_views[k] :
        point.s_views[k] for (k, cone_k) in enumerate(model.cones)]

    return point
end


text:
@testset "K-Means" begin

    Random.seed!(42)

    @testset "Internal functions" begin
        # data inputs
        a = Float64[1, 3, 5, 7, 11, 13]
        b = Float64[2, 4, 6, 8, 12, 14]
        A = Float64[1 3 5; 7 11 13]
        B = Float64[1 3; 5 7; 11 13]

        # euclidean
        @test ClusterAnalysis.euclidean(a, b) == 2.449489742783178

        # squared_error in vector and matrix
        @test ClusterAnalysis.squared_error(a) == 107.33333333333334
        @test ClusterAnalysis.squared_error(A[1, :]) == 8.0

        # totalwithinss
        cluster = [1, 1, 2]
        @test ClusterAnalysis.totalwithinss(B, 2, cluster) == 16.0
    end

    @testset "Kmeans Results" begin
        A = rand(100, 4)
        k = 4

        @testset "test dimensions results" begin
            res1 = kmeans(A, k, init=:kmpp)
            @test length(res1.centroids) == k
            @test length(res1.cluster) == size(A, 1)
            @test length(res1.centroids[1]) == size(A, 2)
        end

        @testset "diferent types of float" begin
            B = rand(Float32, 100, 4)
            res2 = kmeans(B, k, init=:random)
            @test eltype(B) == eltype(res2.centroids[1])
        end

        @testset "convert Int to float" begin
            C = rand(Int64, 100, 4)
            res3 = kmeans(C, k)
            @test Float64 == eltype(res3.centroids[1])
        end

        @testset "test errors" begin
            @test_throws ArgumentError kmeans(A, k, init=:wrongargument)
            @test_throws ArgumentError kmeans([1, 4, 5, 6, 7], k)
            @test_throws MethodError kmeans(A, 1.2)
        end
    end
end


text:
# documented package

err_data = @capture_err begin

  include("../packages/DocumentedCode/src/DocumentedCode.jl")

end

err_data = Boot.strip_output_color(err_data)

@test !contains(err_data, "WARNING: replacing docs for")

@test !contains(err_data, "WARNING: Method definition") &&
  !contains(err_data, "overwritten at")

@test isempty(strip(err_data))


text:
#
# Coherency
#

struct Coherency <: PairwiseStatistic; end
Base.eltype(::Coherency, X::AbstractArray{Complex{T}}) where {T<:Real} = Complex{T}

# Single input matrix
allocwork(::Coherency, X::AbstractVecOrMat{Complex{T}}) where {T<:Real} = nothing
computestat!(::Coherency, out::AbstractMatrix{Complex{T}}, work::Nothing,
             X::AbstractVecOrMat{Complex{T}}) where {T<:Real} =
    cov2coh!(out, Ac_mul_A!(out, X))

# Two input matrices
allocwork(::Coherency, X::AbstractVecOrMat{Complex{T}}, Y::AbstractVecOrMat{Complex{T}}) where {T<:Real} =
    (cov2coh_work(X), cov2coh_work(Y))
function computestat!(::Coherency, out::AbstractMatrix{Complex{T}},
             work::Tuple{Array{T}, Array{T}},
             X::AbstractVecOrMat{Complex{T}}, Y::AbstractVecOrMat{Complex{T}}) where T<:Real
    cov2coh!(out, X, Y, work[1], work[2], Ac_mul_B!(out, X, Y))
end

#
# Coherence (as the absolute value of the correlation matrix)
#

struct Coherence <: PairwiseStatistic; end
Base.eltype(::Coherence, X::AbstractArray{Complex{T}}) where {T<:Real} = T

# Single input matrix
allocwork(::Coherence, X::AbstractVecOrMat{Complex{T}}) where {T<:Real} =
    Array(Complex{T}, nchannels(X), nchannels(X))
computestat!(::Coherence, out::AbstractMatrix{T},
             work::Matrix{Complex{T}},
             X::AbstractVecOrMat{Complex{T}}) where {T<:Real} =
    cov2coh!(out, Ac_mul_A!(work, X), Base.AbsFun())

# Two input matrices
allocwork(::Coherence, X::AbstractVecOrMat{Complex{T}}, Y::AbstractVecOrMat{Complex{T}}) where {T<:Real} =
    (Array(Complex{T}, nchannels(X), nchannels(Y)), cov2coh_work(X), cov2coh_work(Y))
computestat!(::Coherence, out::AbstractMatrix{T},
             work::Tuple{Matrix{Complex{T}}, Array{T}, Array{T}},
             X::AbstractVecOrMat{Complex{T}}, Y::AbstractVecOrMat{Complex{T}}) where {T<:Real} =
    cov2coh!(out, X, Y, work[2], work[3], Ac_mul_B!(work[1], X, Y), Base.AbsFun())

#
# Jackknifing for Coherency/Coherence
#
surrogateval(::Coherence, v) = abs(v)
surrogateval(::Coherency, v) = v

# Single input matrix
allocwork(t::Union{AbstractJackknifeSurrogates{Coherency}, AbstractJackknifeSurrogates{Coherence}},
          X::AbstractVecOrMat{Complex{T}}) where {T<:Real} = (allocwork(t.transform, X), Array(T, div(size(X, 1), jnn(t)), size(X, 2)))
accumulator_array(::AbstractJackknifeSurrogates{Coherency}, work::Nothing, out::AbstractMatrix) = out
accumulator_array(::AbstractJackknifeSurrogates{Coherence}, work::AbstractMatrix, out::AbstractMatrix) = work
function computestat!(t::Union{AbstractJackknifeSurrogates{Coherency}, AbstractJackknifeSurrogates{Coherence}},
                      out::JackknifeSurrogatesOutput,
                      work::Tuple{Union{Matrix{Complex{T}}, Nothing}, Matrix{T}},
                      X::AbstractVecOrMat{Complex{T}}) where T<:Real
    stat = t.transform
    trueval = out.trueval
    surrogates = out.surrogates
    XXc = accumulator_array(t, work[1], out.trueval)
    jnspec = work[2]
    chkinput(trueval, X)
    ntrials(X) % jnn(t) == 0 || throw(DimensionMismatch("ntrials not evenly divisible by $(jnn(t))"))
    size(out.surrogates, 1) == div(ntrials(X), jnn(t)) || throw(DimensionMismatch("invalid output size"))

    Ac_mul_A!(XXc, X)

    # Compute jackknifed power spectra
    @inbounds for j = 1:size(X, 2)
        ssq = real(XXc[j, j])
        for i = 1:size(surrogates, 1)
            ssqdel = ssq
            for idel = (i-1)*jnn(t)+1:i*jnn(t)
                ssqdel -= abs2(X[idel, j])
            end
            jnspec[i, j] = 1/sqrt(ssqdel)
        end
    end

    # Surrogates
    @inbounds for k = 1:size(X, 2)
        for j = 1:k-1, i = 1:size(surrogates, 1)
            v = XXc[j, k]
            for idel = (i-1)*jnn(t)+1:i*jnn(t)
                v -= conj(X[idel, j])*X[idel, k]
            end
            surrogates[i, j, k] = surrogateval(t.transform, v)*(jnspec[i, j]*jnspec[i, k])
        end
        for i = 1:size(surrogates, 1)
            surrogates[i, k, k] = 1
        end
    end

    # Finish true value
    if isa(t.transform, Coherence)
        cov2coh!(trueval, XXc, Base.AbsFun())
    else
        cov2coh!(trueval, XXc, Base.IdFun())
    end

    out
end

# Two input matrices
allocwork(t::AbstractJackknifeSurrogates{Coherency}, X::AbstractVecOrMat{Complex{T}},
          Y::AbstractVecOrMat{Complex{T}}) where {T<:Real} =
    (nothing, cov2coh_work(X), cov2coh_work(Y), Array(T, div(size(X, 1), jnn(t)), size(X, 2)),
     Array(T, div(size(Y, 1), jnn(t)), size(Y, 2)))
allocwork(t::AbstractJackknifeSurrogates{Coherence}, X::AbstractVecOrMat{Complex{T}},
          Y::AbstractVecOrMat{Complex{T}}) where {T<:Real} =
    (Array(Complex{T}, nchannels(X), nchannels(Y)), cov2coh_work(X), cov2coh_work(Y),
     Array(T, div(size(X, 1), jnn(t)), size(X, 2)), Array(T, div(size(Y, 1), jnn(t)), size(Y, 2)))
function computestat!(t::Union{AbstractJackknifeSurrogates{Coherency}, AbstractJackknifeSurrogates{Coherence}},
                      out::JackknifeSurrogatesOutput,
                      work::Tuple{V, Array{T}, Array{T}, Matrix{T}, Matrix{T}},
                      X::AbstractVecOrMat{Complex{T}}, Y::AbstractVecOrMat{Complex{T}}) where {T<:Real,V}
    stat = t.transform
    trueval = out.trueval
    surrogates = out.surrogates
    XYc = accumulator_array(t, work[1], out.trueval)
    Xjnspec = work[4]
    Yjnspec = work[5]
    chkinput(trueval, X, Y)
    ntrials(X) % jnn(t) == 0 || throw(DimensionMismatch("ntrials not evenly divisible by $(jnn(t))"))
    size(out.surrogates, 1) == div(ntrials(X), jnn(t)) || throw(DimensionMismatch("invalid output size"))

    Ac_mul_B!(XYc, X, Y)

    # Compute jackknifed power spectra
    @inbounds for (arr, jnspec) in ((X, Xjnspec), (Y, Yjnspec))
        for j = 1:size(arr, 2)
            ssq = zero(T)
            @simd for i = 1:size(arr, 1)
                ssq += abs2(arr[i, j])
            end
            for i = 1:size(surrogates, 1)
                ssqdel = ssq
                for idel = (i-1)*jnn(t)+1:i*jnn(t)
                    ssqdel -= abs2(arr[idel, j])
                end
                jnspec[i, j] = 1/sqrt(ssqdel)
            end
        end
    end

    # Surrogates
    @inbounds for k = 1:size(Y, 2), j = 1:size(X, 2), i = 1:size(surrogates, 1)
        v = XYc[j, k]
        for idel = (i-1)*jnn(t)+1:i*jnn(t)
            v -= conj(X[idel, j])*Y[idel, k]
        end
        # XXX maybe precompute sqrt for each channel and trial?
        surrogates[i, j, k] = surrogateval(t.transform, v)*(Xjnspec[i, j]*Yjnspec[i, k])
    end

    # Finish true value
    if isa(stat, Coherence)
        cov2coh!(trueval, X, Y, work[2], work[3], XYc, Base.AbsFun())
    else
        cov2coh!(trueval, X, Y, work[2], work[3], XYc, Base.IdFun())
    end

    out
end


text:
boundingbox(d::Vector{T}) where {T <: Number} = minimum(d)..maximum(d)
boundingbox(d::Set{T}) where {T<:Number} = minimum(d)..maximum(d)

"Return the bounding box of the union of two or more bounding boxes."
unionbox(d::Domain) = d
unionbox(d1::Domain, d2::Domain) = unionbox(promote_domains(d1, d2)...)
unionbox(d1::Domain, d2::Domain, domains...) =
    unionbox(unionbox(d1,d2), domains...)

unionbox(d1::Domain{T}, d2::Domain{T}) where {T} = unionbox1(d1, d2)
unionbox1(d1, d2) = unionbox2(d1, d2)
unionbox2(d1, d2) = FullSpace{eltype(d1)}()
unionbox1(d1::EmptySpace, d2) = d2
unionbox1(d1::FullSpace, d2) = d1
unionbox2(d1, d2::EmptySpace) = d1
unionbox2(d1, d2::FullSpace) = d2

unionbox(d1::D, d2::D) where {D<:FixedInterval} = d1

function unionbox(d1::AbstractInterval{T}, d2::AbstractInterval{T}) where {T}
    a, b = endpoints(d1)
    c, d = endpoints(d2)
    A = min(a, c)
    B = max(b, d)
    isinf(A) && isinf(B) ? FullSpace{T}() : A..B
end

unionbox(d1::HyperRectangle{T}, d2::HyperRectangle{T}) where {T} =
    Rectangle{T}(map(unionbox, components(d1), components(d2)))

"Return the bounding box of the intersection of two or more bounding boxes."
intersectbox(d::Domain) = d
intersectbox(d1::Domain, d2::Domain) = intersectbox(promote_domains(d1, d2)...)
intersectbox(d1::Domain, d2::Domain, domains...) =
    intersectbox(intersectbox(d1,d2), domains...)

intersectbox(d1::Domain{T}, d2::Domain{T}) where {T} = intersectbox1(d1, d2)
intersectbox1(d1, d2) = intersectbox2(d1, d2)
intersectbox2(d1, d2) = FullSpace{eltype(d1)}()
intersectbox1(d1::EmptySpace, d2) = d1
intersectbox1(d1::FullSpace, d2) = d2
intersectbox2(d1, d2::EmptySpace) = d2
intersectbox2(d1, d2::FullSpace) = d1

intersectbox(d1::D, d2::D) where {D<:FixedInterval} = d1

intersectbox(d1::AbstractInterval{T}, d2::AbstractInterval{T}) where {T} =
    intersectdomain(d1, d2)

function intersectbox(d1::HyperRectangle{T}, d2::HyperRectangle{T}) where {T}
    d = Rectangle{T}(map(intersectbox, components(d1), components(d2)))
    isempty(d) ? EmptySpace{T}() : d
end

boundingbox(d::AbstractMappedDomain) = map_boundingbox(boundingbox(superdomain(d)), forward_map(d))

function map_boundingbox(box::Interval, fmap)
    l,r = (infimum(box),supremum(box))
    ml = fmap(l); mr = fmap(r)
    min(ml,mr)..max(ml,mr)
end

# This is a best effort implementation, it could be wrong for some maps,
# because we only map the corners. Hence we restrict to affine maps here.
map_boundingbox(box::HyperRectangle, fmap::AbstractAffineMap) =
    map_boundingbox_generic(box, fmap)

# This is a best effort implementation. It could be wrong for some maps,
# because we only map the corners.
function map_boundingbox_generic(box::HyperRectangle{T}, fmap) where {T}
    mapped_corners = map(fmap, corners(box))
    left = [minimum(x[j] for x in mapped_corners) for j in 1:dimension(box)]
    right = [maximum(x[j] for x in mapped_corners) for j in 1:dimension(box)]
    Rectangle{T}(left, right)
end


text:
import DrWatson: quickactivate
quickactivate(@__DIR__, "Chemostat_Human1")

# --------------------------------------------------------------------
using ArgParse

set = ArgParseSettings()
@add_arg_table! set begin
    "--name", "-n"
        help = "A name that uniquely identify the simulation, it is used in caching"
        default = "maxent_ep_v1"
        required = false
    "--beta0"
        help = "the lower limit of the range to search experimental beta"
        default = "0.0"
        required = false
    "--beta1"
        help = "the upper limit of the range to search experimental beta"
        default = "10000.0"
        required = false
end

if isinteractive()
    # This must uniquely identify the script
    const SIM_IDER = "maxent_ep_v1"
    const beta0 = 0.0
    const beta1 = 10000.0
else
    const parsed_args = parse_args(set)
    # This must uniquely identify the script
    const SIM_IDER = parsed_args["name"]
    const beta0 = parse(Float64, parsed_args["beta0"])
    const beta1 = parse(Float64, parsed_args["beta1"])

end
println("sim_ider: ", SIM_IDER)
println("beta0: ", beta0)
println("beta1: ", beta1)

## --------------------------------------------------------------------
import MAT
import SparseArrays

import Chemostat_Human1
const H1 = Chemostat_Human1

import Chemostat
const Ch = Chemostat
const ChU = Ch.Utils
const ChLP = Ch.LP
const ChSU = Ch.SimulationUtils
ChU.set_cache_dir(H1.CACHE_DIR)

## --------------------------------------------------------------------
# Tools
fba_objval(model) = ChLP.fba(model, H1.BIOMASS_IDER).obj_val
# data/processed/fva_pp_models/constLevel_0/HS_578T__ecModel.bson

## --------------------------------------------------------------------
# here I should load a fva_pp_model
cell_line = "HS_578T"
const_level = 0
model_type = :ecModel
model_file = H1.get_human1_model_file(cell_line, const_level, model_type)
@assert isfile(model_file)
exp_grate = H1.get_NCI60_exchange(model_type, cell_line, const_level, :exp, "biomass")
println("exp objval:   ", exp_grate)
model = ChU.read_mat(model_file);
# model = ChU.clampfields!(model, [:b, :lb, :ub]; abs_max = 100.0, zeroth = 1e-8)
println("model size:   ", size(model))
println("model objval: ", fba_objval(model))
# model = ChLP.fva_preprocess(model; check_obj = H1.BIOMASS_IDER)
# println("model objval: ", fba_objval(model))
# fvafile = joinpath(H1.FVA_PP_MODELS_DIR, "fva_pp_model.bson")
# ChU.save_data(fvafile, ChU.compressed_model(model))

## --------------------------------------------------------------------
# load cache Test
# model = ChU.load_data(model_file);
# model = ChU.uncompressed_model(model)
# model = ChU.clampfields!(model, [:b, :lb, :ub]; abs_max = 100.0, zeroth = 1e-8)
# println("model size:   ", size(model))
# println("model objval: ", fba_objval(model))

## --------------------------------------------------------------------
# epouts
const cache_id = (:EP_OUTS, SIM_IDER)
dat = ChU.load_cache(cache_id; verbose = false)
const epouts = something(dat, Dict())

## --------------------------------------------------------------------
# find exp_beta
exp_beta = ChSU.find_beta(model; beta0, beta1,
    obj_ider = H1.BIOMASS_IDER, errorth = 0.1,
    target_objval = exp_grate, epouts,
    after_maxent_ep = function(epout)
        # caching results
        ChU.save_cache(cache_id, epouts)
    end
)

## --------------------------------------------------------------------
# saving
ChU.save_data("maxent_ep_data.bson", (exp_beta, ChU.compressed_model(model), epouts))


text:
include("ispc_toolchain.jl")

# read the code from a file
code = open(f -> read(f, String), "webql.ispc")

# compile the SPMD code
lib = load_ispc(code, `--opt=fast-math --addressing=32`)

# get function pointers
const radial_spec_fptr = Libc.Libdl.dlsym(lib, "calculate_radial_spectrumF16")
const square_spec_fptr = Libc.Libdl.dlsym(lib, "calculate_square_spectrumF16")

const radial_view_fptr = Libc.Libdl.dlsym(lib, "calculate_radial_viewport_spectrumF16")
const square_view_fptr = Libc.Libdl.dlsym(lib, "calculate_square_viewport_spectrumF16")

const make_linear_video_frame_fptr = Libc.Libdl.dlsym(lib, "make_video_frameF16_linear")
const make_logistic_video_frame_fptr = Libc.Libdl.dlsym(lib, "make_video_frameF16_logistic")
const make_ratio_video_frame_fptr = Libc.Libdl.dlsym(lib, "make_video_frameF16_ratio")
const make_square_video_frame_fptr = Libc.Libdl.dlsym(lib, "make_video_frameF16_square")
const make_legacy_video_frame_fptr = Libc.Libdl.dlsym(lib, "make_video_frameF16_legacy")