text:
"""
`module  JAC.Pulse`
    ... a submodel of JAC that contains all structs and methods to deal with time-dependent pulses of the em field.
"""
module Pulse

    using   ..Basics, ..Defaults, ..Radial, GSL, SpecialFunctions

    export  AbstractEnvelope, AbstractBeam

    """
    `abstract type Pulse.AbstractEnvelope`
        ... defines an abstract type to comprise various envelopes of (possible) laser pulses in terms of their shape,
            pulse duration or number of cycles, etc.

        + InfiniteEnvelope         ... to represent an infinte (plane-wave) pulse.
        + RectangularEnvelope      ... to represent a finite rectangular pulse.
        + SinSquaredEnvelope       ... to represent a finite sin^2 pulse.
        + GaussianEnvelope         ... to represent a Gaussian light pulse.
    """
    abstract type  AbstractEnvelope  end

    """
    `struct  Pulse.InfiniteEnvelope  <: Pulse.AbstractEnvelope`   ... to represent an infinte (plane-wave) pulse.
    """
    struct   InfiniteEnvelope        <: Pulse.AbstractEnvelope    end

    function Base.string(env::InfiniteEnvelope)
        sa = "infinite pulse"
        return( sa )
    end

    function Base.show(io::IO, env::InfiniteEnvelope)
        sa = string(env);       print(io, sa, "\n")
    end

    """
    `struct  Pulse.RectangularEnvelope  <: Pulse.AbstractEnvelope`   ... to represent a finte rectangular pulse.

        + cycles      ::Int64     ... Number of cycles of the pulse.
    """
    struct   RectangularEnvelope       <: Pulse.AbstractEnvelope
        cycles        ::Int64
    end

    function Base.string(env::RectangularEnvelope)
        sa = "rectangular pulse of $(env.cycles) cycles"
        return( sa )
    end

    function Base.show(io::IO, env::RectangularEnvelope)
        sa = string(env);       print(io, sa, "\n")
    end

    """
    `struct  Pulse.SinSquaredEnvelope  <: Pulse.AbstractEnvelope`   ... to represent a finite sin^2 pulse.

        + cycles      ::Int64     ... Number of cycles of the pulse.
    """
    struct   SinSquaredEnvelope        <: Pulse.AbstractEnvelope
        cycles      ::Int64
    end

    function Base.string(env::SinSquaredEnvelope)
        sa = "sin^2 pulse of $(env.cycles) cycles"
        return( sa )
    end

    function Base.show(io::IO, env::SinSquaredEnvelope)
        sa = string(env);       print(io, sa, "\n")
    end

    """
    `struct  Pulse.GaussianEnvelope  <: Pulse.AbstractEnvelope`   ... to represent a (infinite) Gaussian pulse.

        + fwhm        ::Float64     ... FWHM which is often taken as pulse duration
    """
    struct   GaussianEnvelope        <: Pulse.AbstractEnvelope
        fwhm          ::Float64
    end

    function Base.string(env::GaussianEnvelope)
        sa = "Gaussian pulse with pulse duration or FWHM = $(env.fwhm)"
        return( sa )
    end

    function Base.show(io::IO, env::GaussianEnvelope)
        sa = string(env);       print(io, sa, "\n")
    end

    """
    `abstract type Pulse.AbstractBeam`
        ... defines an abstract type to comprise various basic laser pulses that are characterized in terms of their amplitude,
            frequency, carrier-envelope phase, etc. In general, the basic beam properties are independent of the (pulse) envelope
            and the polarization properties which are handled and communicated separately (to and within the program).

        + PlaneWaveBeam            ... to represent a plane-wave beam.
        + BesselBeam               ... to represent a Bessel beam (not yet).
    """
    abstract type  AbstractBeam  end

    """
    `struct  Pulse.PlaneWaveBeam  <: AbstractBeam`
        ... to represent a plane-wave beam with given (real) amplitude, frequency and carrier-envelope phase.

        + A0            ::Float64            ... (Constant) Amplitude of the light pulse.
        + omega         ::Float64            ... Central frequency.
        + cep           ::Float64            ... Carrier-envelope phase.
    """
    struct   PlaneWaveBeam  <: AbstractBeam
        A0              ::Float64
        omega           ::Float64
        cep             ::Float64
    end

    """
    `Pulse.PlaneWaveBeam()`  ... constructor for an `empty` instance of Pulse.PlaneWaveBeam().
    """
    function PlaneWaveBeam()
        PlaneWaveBeam( 0., 0., 0. )
    end

    function Base.string(beam::PlaneWaveBeam)
        sa = "Plane-wave beam/pulse with amplitude A0=$(beam.A0), frequency omega=$(beam.omega) a.u. and carrier-envelope phase cep=$(beam.cep)"
        return( sa )
    end

    function Base.show(io::IO, beam::PlaneWaveBeam)
        sa = string(beam);       print(io, sa, "\n")
    end

    #======================================================================================
    #  These data types need to be worked through
    """
    `struct  Pulse.Envelope`  ... defines a type for the envelope (function) of an em pulse.

        + timeDelay          ::Float64             ... Time delay with regard to an arbitrary time t=0, often used for the first pulse.
        + f0                 ::Float64             ... Amplitude of the envelope function f(t) = f0 * f_shape(t)
        + shapeFunction      ::Array{Float64,1}    ... (Normalized) shape function of the pulse.
    """
    struct Envelope
        timeDelay            ::Float64
        f0                   ::Float64
        shapeFunction        ::Array{Float64,1}
    end

    """
    `Pulse.Envelope()`  ... constructor for an `empty` instance of Pulse.Envelope().
    """
    function Envelope()
        Envelope( 0., 0., Float64[] )
    end

    # `Base.show(io::IO, envelope::Pulse.Envelope)`  ... prepares a proper printout of the variable envelope::Pulse.Envelope.
    function Base.show(io::IO, envelope::Pulse.Envelope)
        println(io, "timeDelay:             $(envelope.timeDelay)  ")
        println(io, "f0:                    $(envelope.f0)  ")
        println(io, "shapeFunction:         $(envelope.shapeFunction)  ")
    end

    """
    `@enum   PolarizationType`  ... defines a type of polarization of an experimentally described light pulse

    + NoType, Linear, LeftCircular, RightCircular, Elliptical      ... with obvious meaning.
    """
    @enum   PolarizationType    NoType    Linear    LeftCircular    RightCircular    Elliptical

    """
    `Pulse.PolarizationType(sa::String)`  ... constructor for a given String.
    """
    function PolarizationType(sa::String)
        if       sa == "linear"                  wa = Linear
        elseif   sa == "left-circular"           wa = LeftCircular
        elseif   sa == "right-circular"          wa = RightCircular
        elseif   sa == "elliptical"              wa = Elliptical
        else
            error("stop a")
        end
        PolarizationType(wa)
    end

    # `Base.show(io::IO, ptyp::Pulse.PolarizationType)`  ... prepares a proper printout of the variable ptyp::Pulse.PolarizationType.
    function Base.show(io::IO, ptyp::Pulse.PolarizationType)
        print(io, string(ptyp) )
    end

    """
    `Base.string(ptyp::Pulse.PolarizationType)`  ... provides a proper printout of the variable ptyp::Pulse.PolarizationType.
    """
    function Base.string(ptyp::Pulse.PolarizationType)
        if      ptyp == NoType              return( "undefined polarization" )
        elseif  ptyp == Linear              return( "linear polarization" )
        elseif  ptyp == LeftCircular        return( "left-circular polarization" )
        elseif  ptyp == RightCircular       return( "right-circular polarization" )
        elseif  ptyp == Elliptical          return( "elliptical polarization" )
        elseif  ptyp == Linear              return( "linear polarization" )
        else    error("stop a")
        end
    end

    """
    `struct  Pulse.Polarization`  ... defines a type for the polarization of an em pulse.

        + type               ::Pulse.PolarizationType    ... General type of polarization.
        + circularDegree     ::Float64                   ... Degree of circular polarization.
        + linearDegree       ::Float64                   ... Degree of linear polarization.
        + linearPolarization ::SolidAngle                ... Direction of the (linear) polarization vector as described by the
                                                             unit vector u = u(Omega') with regard to the x'-z' plane of the incident pulse.
        + has_gCoefficients  ::Bool                      ... True, of the complex gPlus and gMinus coefficients are properly defined.
        + gPlus              ::Complex{Float64}          ... g_+1 coefficient.
        + gMinus             ::Complex{Float64}          ... g_-1 coefficient.
    """
    struct Polarization
        typ                  ::Pulse.PolarizationType
        circularDegree       ::Float64
        linearDegree         ::Float64
        linearPolarization   ::SolidAngle
        has_gCoefficients    ::Bool
        gPlus                ::Complex{Float64}
        gMinus               ::Complex{Float64}
    end

    """
    `Pulse.Polarization()`  ... constructor for an `empty` instance of Pulse.Polarization().
    """
    function Polarization()
        Polarization( NoType, 0., 0., SolidAngle(0., 0.), false, Complex(0.), Complex(0.) )
    end

    # `Base.show(io::IO, polarization::Pulse.Polarization)`  ... prepares a proper printout of the variable polarization::Pulse.Polarization.
    function Base.show(io::IO, polarization::Pulse.Polarization)
        println(io, "type:                    $(polarization.typ)  ")
        println(io, "circularDegree:          $(polarization.circularDegree)  ")
        println(io, "linearDegree:            $(polarization.linearDegree)  ")
        println(io, "linearPolarization:      $(polarization.linearPolarization)  ")
        println(io, "has_gCoefficients:       $(polarization.has_gCoefficients)  ")
        println(io, "gPlus:                   $(polarization.gPlus)  ")
        println(io, "gMinus:                  $(polarization.gMinus)  ")
    end

    """
    `struct  Pulse.ExperimentalCharacterization`  ... defines a type for characterizing an experimental or physically described em pulse.

        + kind             ::Pulse.Shape            ... General shape of a given pulse.
        + propagation      ::SolidAngle             ... Propagation direction of the pulse as described by the unit vector u = u(Omega).
        + omega            ::Float64                ... Central frequency of the em pulse.
        + maxIntensity     ::Float64                ... maximum intensity in units ...
        + fwhm             ::Float64                ... approximate full-width-half maximum (fwhm) for rather long pulses; the nearest widths
                                                        with an integer No. of cycles is taken.
        + timeDelay        ::Float64                ... time delay with regard to some arbitrarely chosen time t = 0
                                                        (often taken for the first pulse).
        + NoCycles         ::Int64                  ... No. of cycles which always determines the lengths of the pulse.
        + multipoles       ::Array{EmMultipole,1}   ... Multipoles of the em field to be included in the description of the light field.
        + polarization     ::Pulse.Polarization     ... Describes the polarization properties of the em pulse.
    """
    struct ExperimentalCharacterization
        kind               ::Pulse.Shape
        propagation        ::SolidAngle
        omega              ::Float64
        maxIntensity       ::Float64
        fwhm               ::Float64
        timeDelay          ::Float64
        NoCycles           ::Int64
        multipoles         ::Array{EmMultipole,1}
        polarization       ::Pulse.Polarization
    end

    """
    `Pulse.ExperimentalCharacterization()`  ... constructor for an `empty` instance of ExperimentalCharacterization().
    """
    function ExperimentalCharacterization()
        ExperimentalCharacterization( NoShape, SolidAngle(0., 0.), 0., 0., 0., 0., 0, EmMultipole[], Polarization() )
    end

    # `Base.show(io::IO, chz::Pulse.ExperimentalCharacterization)`
    #		... prepares a proper printout of the variable  chz::Pulse.ExperimentalCharacterization.
    function Base.show(io::IO, characterization::Pulse.ExperimentalCharacterization)
        println(io, "kind:              $(chz.kind)  ")
        println(io, "propagation:       $(chz.propagation)  ")
        println(io, "omega:             $(chz.omega)  ")
        println(io, "maxIntensity:      $(chz.maxIntensity)  ")
        println(io, "fwhm:              $(chz.fwhm)  ")
        println(io, "timeDelay:         $(chz.timeDelay)  ")
        println(io, "NoCycles:          $(chz.NoCycles)  ")
        println(io, "multipoles:        $(chz.multipoles)  ")
        println(io, "polarization:      $(chz.polarization)  ")
    end

    """
    `struct  Pulse.Gaussian`  ... defines a type for a Gaussian light pulse that is used for evaluating time-dependent statistical tensors.

        + propagation      ::SolidAngle             ... Propagation direction of the pulse as described by the unit vector u = u(Omega).
        + omega            ::Float64                ... Central frequency of the em pulse.
        + multipoles       ::Array{EmMultipole,1}   ... Multipoles of the em field to be included in the description of the light field.
        + envelope         ::Pulse.Envelope         ... Envelope (function) of the light pulse.
        + polarization     ::Pulse.Polarization     ... Polarization of the light pulse with typically well-defined g_+1 and g_-1 coefficients.
    """
    struct Gaussian
        propagation        ::SolidAngle
        omega              ::Float64
        multipoles         ::Array{EmMultipole,1}
        envelope           ::Pulse.Envelope
        polarization       ::Pulse.Polarization
    end

    """
    `Pulse.Gaussian()`  ... constructor for an `empty` instance of Pulse.Gaussian().
    """
    function Gaussian()
        Gaussian( SolidAngle(0., 0.), 0., EmMultipole[], Pulse.Envelope(), Pulse.Polarization()  )
    end

    # `Base.show(io::IO, gaussian::Pulse.Gaussian)`  ... prepares a proper printout of the variable gaussian::Pulse.Gaussian.
    function Base.show(io::IO, gaussian::Pulse.Gaussian)
        println(io, "propagation:        $(gaussian.propagation)  ")
        println(io, "omega:              $(gaussian.omega)  ")
        println(io, "multipoles:         $(gaussian.multipoles)  ")
        println(io, "envelope:           $(gaussian.envelope)  ")
        println(io, "polarization:       $(gaussian.polarization)  ")
    end
    ============================================================================================#

    ###################################################################################################################
    ###################################################################################################################
    ###################################################################################################################

    """
    `Pulse.computeFieldAmplitude(intensity::Float64, omega::Float64)`
        ... compute the field amplitude from the given (maximum) intensity [in a.u.] and the central frequency [in a.u.]
            of the light field: A0 = sqrt( 8pi * alpha * intensity) / omega.
    """
    function computeFieldAmplitude(intensity::Float64, omega::Float64)
        wa = sqrt(8 * pi * Defaults.getDefaults("alpha") * intensity) / omega
        return( wa )
    end

    """
    `Pulse.envelopeVolkovIntegral(plus::Bool, envelope::Pulse.InfiniteEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                  thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)`
        ... evaluates the envelope (time) integral F^(Volkov) [+/-; omega; f^(infinite); A; angles & energies] for an infinite pulse with given
            parameters; an ntg::Complex{Float64} is retured.
    """
    function envelopeVolkovIntegral(plus::Bool, envelope::Pulse.InfiniteEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                    thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)
        wa = 0. * im;
        phiCep = beam.cep;       a = beam.A0 * sqrt(2*energyp) * sin(thetap) / sqrt(2) / beam.omega;   Up = beam.A0^2 / 4
        lambda = Basics.determinePolarizationLambda(polarization)
        #
        # Compute the summation over the Bessel functions first; start with value for s = 0
        for  s = -10:10
            if  plus   wb = Basics.diracDelta((s-1)*beam.omega + energyp - initialEn + Up, 1.0e-3)
            else       wb = Basics.diracDelta((s+1)*beam.omega + energyp - initialEn + Up, 1.0e-3)
            end
            #
            if  wb != 0.
                wa = wa + GSL.sf_bessel_Jn(s, a) * exp(im*s * (phiCep - lambda*phip)) * wb
            end
        end
        if  plus   wa = wa * 2pi * beam.A0 * exp(-im*phiCep)
        else       wa = wa * 2pi * beam.A0 * exp(im*phiCep)
        end

        return( wa )
    end

    """
    `Pulse.envelopeVolkovIntegral(plus::Bool, envelope::Pulse.RectangularEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                  thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)`
        ... evaluates the envelope (time) integral F^(Volkov) [+/-; omega; f^(rectangular); A; angles & energies] for an infinite pulse with given
            parameters; an ntg::Complex{Float64} is retured.
    """
    function envelopeVolkovIntegral(plus::Bool, envelope::Pulse.RectangularEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                    thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)
        wa = 0. * im;   Tp = 2pi * envelope.cycles / beam.omega
        phiCep = beam.cep;       a = beam.A0 * sqrt(2*energyp) * sin(thetap) / ( sqrt(2) * beam.omega );   Up = beam.A0^2 / 4
        lambda = Basics.determinePolarizationLambda(polarization)
        phaseConstant = phiCep - lambda*phip
        #
        # Compute the summation over the Bessel functions first; start with value for s = 0
        for  s = -20:20
            if  plus   wb = (s-1)*beam.omega + energyp - initialEn + Up
            else       wb = (s+1)*beam.omega + energyp - initialEn + Up
            end
            #
            if  wb != 0.
                wa = wa + GSL.sf_bessel_Jn(s, a) * exp(im*s * phaseConstant) / wb * ( exp(im*wb*Tp) - 1 )
            end
        end
        if  plus   wa = wa * (-im) * exp(-im * a * sin(phaseConstant)) * beam.A0 * exp(-im*phiCep)
        else       wa = wa * (-im) * exp(-im * a * sin(phaseConstant)) * beam.A0 * exp(im*phiCep)
        end

        return( wa )
    end

    """
    `Pulse.envelopeVolkovIntegral(plus::Bool, envelope::Pulse.SinSquaredEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                  thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)`
        ... evaluates the envelope (time) integral F^(Volkov) [+/-; omega; f^(rectangular); A; angles & energies] for a sine-squared pulse with given
            parameters; an ntg::Complex{Float64} is retured.
    """
    function envelopeVolkovIntegral(plus::Bool, envelope::Pulse.SinSquaredEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                    thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)
        wa = 0. * im;   np = envelope.cycles;   Tp = 2pi * np / beam.omega
        phiCep = beam.cep;  a = beam.A0 * sqrt(2*energyp) * sin(thetap) / ( sqrt(2) * beam.omega );   Up = beam.A0^2 / 4
        sinSqrArg = 0.5 * beam.omega / np;  volkovFac = 0.25/sinSqrArg * Up
        lambda = Basics.determinePolarizationLambda(polarization)
        phaseConstant = phiCep - lambda*phip

        kfac = [ (1-1/np), 1, (1+1/np) ]
        ack = a * [ -0.25/kfac[1], 0.5/kfac[2], -0.25/kfac[3] ]

        #Define Gauss-Legendre grid, convergence is typically good for orderGL = 100 * np (time consuming for np > 10); tested up to np = 20
        if  np <= 10     orderGL = 100*np
        else             orderGL = 1000
        end
        gaussLegendre = Radial.GridGL("Finite",0.0,Tp,orderGL)
        tgrid = gaussLegendre.t
        weights = gaussLegendre.wt

        #Sum over grid and compute Gauss-Legendre sum
        for    j = 1:orderGL
            t = tgrid[j]

            #Compute Volkov phase at gridpoint t
            SVolkov = ( energyp + 3/8 * Up ) * t - volkovFac * sin(2*sinSqrArg*t) + 0.125 * volkovFac * sin(4*sinSqrArg*t)
            for k = 1:3
                SVolkov = SVolkov + ack[k] * ( sin( beam.omega * kfac[k] * t + phaseConstant ) - sin(phaseConstant) )
            end

            #Compute integrand at gridpoint t
            if  plus    integrand = sin( sinSqrArg * t )^2 * exp( -im * ( ( initialEn + beam.omega ) * t - SVolkov ) )
            else        integrand = sin( sinSqrArg * t )^2 * exp( -im * ( ( initialEn - beam.omega ) * t - SVolkov ) )
            end

            #Gauss-Legendre sum
            wa = wa + weights[j] * integrand
        end

        #Multiply with global factor
        if  plus    wa = beam.A0 * exp(-im * phiCep) * wa
        else        wa = beam.A0 * exp(im * phiCep) * wa
        end

        return( wa )
    end

    """
    `Pulse.envelopeVolkovIntegral(plus::Bool, envelope::Pulse.GaussianEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                  thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)`
        ... evaluates the envelope (time) integral F^(Volkov) [+/-; omega; f^(rectangular); A; angles & energies] for a gaussian pulse with given
            parameters; an ntg::Complex{Float64} is retured.
    """
    function envelopeVolkovIntegral(plus::Bool, envelope::Pulse.GaussianEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                    thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64 )
        wa = 0. * im;  Tp = envelope.fwhm
        phiCep = beam.cep;   Up = beam.A0^2 / 4
        a1 = 0.25 * Up * sqrt(pi/log(4)) * Tp
        a2 = 0.125 * Tp * sqrt(pi/log(2)) * beam.A0 * sqrt(2*energyp) * sin(thetap) / sqrt(2) * exp( -beam.omega^2 * Tp^2 / log( 65536 ) )
        lambda = Basics.determinePolarizationLambda(polarization)
        phaseConstant = phiCep - lambda*phip

        # Factor from (linear) variable transformation tau -> t in order to allow Gauss-Hermite integration ( weight function e^(-t^2) )
        transformFac = Tp/(2*sqrt(log(2)))

        #Define Gauss-Hermite grid, convergence is typically good for orderGH = 10000
        orderGH = 10000
        gaussHermite = Radial.GridGH(orderGH)
        tgrid = gaussHermite.t
        weights = gaussHermite.wt

        # Sum over grid and compute Gauss-Legendre sum
        for    j = 1:orderGH
            t = transformFac * tgrid[j]

            #Compute Volkov phase at gridpoint t
            SVolkov = energyp * t + a1 * erf(2 * sqrt(log(4)) * t/Tp)
                        + a2 * ( exp( -im*phaseConstant ) * erf( (im*Tp^2*beam.omega + 8*log(2)*t)/(4*Tp*sqrt(log(2))) )
                               - exp(  im*phaseConstant ) * erf( (im*Tp^2*beam.omega - 8*log(2)*t)/(4*Tp*sqrt(log(2))) ) )

            # Compute integrand at gridpoint t
            if  plus    integrand = exp( -im * ( ( initialEn + beam.omega ) * t - SVolkov ) )
            else        integrand = exp( -im * ( ( initialEn - beam.omega ) * t - SVolkov ) )
            end

            # Gauss-Hermite sum
            wa = wa + weights[j] * integrand
        end

        # Multiply with global factor
        if  plus    wa = transformFac * beam.A0 * exp(-im * phiCep) * wa
        else        wa = transformFac * beam.A0 * exp(im * phiCep) * wa
        end

        return( wa )
    end

    """
    `Pulse.envelopeVolkovIntegral(plus::Bool, envelope::Pulse.GaussianEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                  thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64, orderGH::Int64)`
        ... evaluates the envelope (time) integral F^(Volkov) [+/-; omega; f^(rectangular); A; angles & energies] for a gaussian pulse with given
            parameters; an ntg::Complex{Float64} is retured.
    """
    function envelopeVolkovIntegral(plus::Bool, envelope::Pulse.GaussianEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                    thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64, orderGH::Int64)
        wa = 0. * im;  Tp = envelope.fwhm
        phiCep = beam.cep;   Up = beam.A0^2 / 4
        a1 = 0.25 * Up * sqrt(pi/log(4)) * Tp
        a2 = 0.125 * Tp * sqrt(pi/log(2)) * beam.A0 * sqrt(2*energyp) * sin(thetap) / sqrt(2) * exp( -beam.omega^2 * Tp^2 / log( 65536 ) )
        lambda = Basics.determinePolarizationLambda(polarization)
        phaseConstant = phiCep - lambda*phip

        # Factor from (linear) variable transformation tau -> t in order to allow Gauss-Hermite integration ( weight function e^(-t^2)
        transformFac = Tp/(2*sqrt(log(2)))

        gaussHermite = Radial.GridGH(orderGH)
        tgrid = gaussHermite.t
        weights = gaussHermite.wt

        # Sum over grid and compute Gauss-Legendre sum
        for    j = 1:orderGH
            t = transformFac * tgrid[j]

            #Compute Volkov phase at gridpoint t
            SVolkov = energyp * t + a1 * erf(2 * sqrt(log(4)) * t/Tp)
                        + a2 * ( exp( -im*phaseConstant ) * erf( (im*Tp^2*beam.omega + 8*log(2)*t)/(4*Tp*sqrt(log(2))) )
                               - exp(  im*phaseConstant ) * erf( (im*Tp^2*beam.omega - 8*log(2)*t)/(4*Tp*sqrt(log(2))) ) )

            # Compute integrand at gridpoint t
            if  plus    integrand = exp( -im * ( ( initialEn + beam.omega ) * t - SVolkov ) )
            else        integrand = exp( -im * ( ( initialEn - beam.omega ) * t - SVolkov ) )
            end

            # Gauss-Hermite sum
            wa = wa + weights[j] * integrand
        end

        #Multiply with global factor
        if  plus    wa = transformFac * beam.A0 * exp(-im * phiCep) * wa
        else        wa = transformFac * beam.A0 * exp(im * phiCep) * wa
        end

        return( wa )
    end

    """
    `Pulse.envelopeQuadVolkovIntegral(envelope::Pulse.InfiniteEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                      thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)`
        ... evaluates the envelope (time) integral F^(quad, Volkov) [f^(infinite); A; angles & energies] for an infinite pulse with given
            parameters; an ntg::Complex{Float64} is retured.
    """
    function envelopeQuadVolkovIntegral(envelope::Pulse.InfiniteEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                        thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)
        wa = 0. * im;
        phiCep = beam.cep;       a = beam.A0 * sqrt(2*energyp) * sin(thetap) / sqrt(2) / beam.omega;   Up = beam.A0^2 / 4
        lambda = Basics.determinePolarizationLambda(polarization)
        #
        # Compute the summation over the Bessel functions first; start with value for s = 0
        for  s = -10:10
            wb = Basics.diracDelta(s*beam.omega + energyp - initialEn + Up, 1.0e-3)
            #
            if  wb != 0.    wa = wa + GSL.sf_bessel_Jn(s, a) * exp(im*s * (phiCep - lambda*phip)) * wb
            end
        end
        wa = wa * 4pi * Up
        return( wa )
    end

    """
    `Pulse.envelopeQuadVolkovIntegral(envelope::Pulse.RectangularEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                      thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)`
        ... evaluates the envelope (time) integralF^(quad, Volkov) [f^(infinite); A; angles & energies]  for a rectangular pulse with given
            parameters; an ntg::Complex{Float64} is retured.
    """
    function envelopeQuadVolkovIntegral(envelope::Pulse.RectangularEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                        thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)
        wa = 0. * im;   Tp = 2pi * envelope.cycles / beam.omega
        phiCep = beam.cep;       a = beam.A0 * sqrt(2*energyp) * sin(thetap) / sqrt(2) / beam.omega;   Up = beam.A0^2 / 4
        lambda = Basics.determinePolarizationLambda(polarization)
        phaseConstant = phiCep - lambda*phip
        @show beam.omega
        # Compute the summation over the Bessel functions first; start with value for s = 0
        for  s = -20:20
            wb = s*beam.omega + energyp - initialEn + Up
            #
            if  wb != 0.    wa = wa + GSL.sf_bessel_Jn(s, a) * exp(im*s * phaseConstant) / wb * ( exp(im*wb*Tp) - 1 )
            end
        end
        wa = wa * (-im) * 2 * Up * exp(-im * a * sin(phaseConstant))

        return( wa )
    end

    """
    `Pulse.envelopeQuadVolkovIntegral(envelope::Pulse.SinSquaredEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                      thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)`
        ... evaluates the envelope (time) integral F^(quad, Volkov) [f^(infinite); A; angles & energies] for a sine-squared pulse with given
            parameters; an ntg::Complex{Float64} is retured.
    """
    function envelopeQuadVolkovIntegral(envelope::Pulse.SinSquaredEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                    thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)
        wa = 0. * im;   np = envelope.cycles;   Tp = 2pi * np / beam.omega
        phiCep = beam.cep;  a = beam.A0 * sqrt(2*energyp) * sin(thetap) / ( sqrt(2) * beam.omega );   Up = beam.A0^2 / 4
        sinSqrArg = 0.5 * beam.omega / np;  volkovFac = 0.25/sinSqrArg * Up
        lambda = Basics.determinePolarizationLambda(polarization)
        phaseConstant = phiCep - lambda*phip

        kfac = [ (1-1/np), 1, (1+1/np) ]
        ack = a * [ -0.25/kfac[1], 0.5/kfac[2], -0.25/kfac[3] ]

        # Define Gauss-Legendre grid, convergence is typically good for orderGL = 100 * np (time consuming for np > 10); tested up to np = 20
        if  np <= 10     orderGL = 100*np
        else             orderGL = 1000
        end
        gaussLegendre = Radial.GridGL("Finite",0.0,Tp,orderGL)
        tgrid = gaussLegendre.t
        weights = gaussLegendre.wt

        # Sum over grid and compute Gauss-Legendre sum
        for    j = 1:orderGL
            t = tgrid[j]

            #Compute Volkov phase at gridpoint t
            SVolkov = ( energyp + 3/8 * Up ) * t - volkovFac * sin(2*sinSqrArg*t) + 0.125 * volkovFac * sin(4*sinSqrArg*t)
            for k = 1:3
                SVolkov = SVolkov + ack[k] * ( sin( beam.omega * kfac[k] * t + phaseConstant ) - sin(phaseConstant) )
            end

            # Compute integrand at gridpoint t
            integrand = 0.5 * sin( sinSqrArg * t )^4 * exp( -im * ( initialEn * t - SVolkov ) )

            # Gauss-Legendre sum
            wa = wa + weights[j] * integrand
        end

        # Multiply with global factor
        wa = beam.A0^2 * wa

        return( wa )
    end

    """
    `Pulse.envelopeQuadVolkovIntegral(envelope::Pulse.GaussianEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                      thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)`
        ... evaluates the envelope (time) integral F^(quad, Volkov) [f^(infinite); A; angles & energies] for a gaussian pulse with given
            parameters; an ntg::Complex{Float64} is retured.
    """
    function envelopeQuadVolkovIntegral(envelope::Pulse.GaussianEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                    thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)
        wa = 0. * im;  Tp = envelope.fwhm
        phiCep = beam.cep;   Up = beam.A0^2 / 4
        a1 = 0.25 * Up * sqrt(pi/log(4)) * Tp
        a2 = 0.125 * Tp * sqrt(pi/log(2)) * beam.A0 * sqrt(2*energyp) * sin(thetap) / sqrt(2) * exp( -beam.omega^2 * Tp^2 / log( 65536 ) )
        lambda = Basics.determinePolarizationLambda(polarization)
        phaseConstant = phiCep - lambda*phip

        # Factor from (linear) variable transformation tau -> t in order to allow Gauss-Hermite integration ( weight function e^(-t^2) )
        transformFac = Tp/(2*sqrt(2*log(2)))

        # Define Gauss-Hermite grid, convergence is typically good for ......
        orderGH = 10000
        gaussHermite = Radial.GridGH(orderGH)
        tgrid = gaussHermite.t
        weights = gaussHermite.wt

        # Sum over grid and compute Gauss-Legendre sum
        for    j = 1:orderGH
            t = transformFac * tgrid[j]

            # Compute Volkov phase at gridpoint t
            SVolkov = energyp * t + a1 * erf(2 * sqrt(log(4)) * t/Tp)
                        + a2 * ( exp( -im*phaseConstant ) * erf( (im*Tp^2*beam.omega + 8*log(2)*t)/(4*Tp*sqrt(log(2))) )
                               - exp( im*phaseConstant )  * erf( (im*Tp^2*beam.omega - 8*log(2)*t)/(4*Tp*sqrt(log(2))) ) )

            # Compute integrand at gridpoint t
            integrand = exp( -im * ( initialEn * t - SVolkov ) )

            # Gauss-Hermite sum
            wa = wa + weights[j] * integrand
        end

        # Multiply with global factor
        wa = 0.5 * transformFac * beam.A0^2 * wa

        return( wa )
    end

    """
    `Pulse.envelopeQuadVolkovIntegral(envelope::Pulse.AbstractEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                      thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)`
        ... evaluates the envelope (time) integralF^(quad, Volkov) [f^(infinite); A; angles & energies]  for all pulses for which no analytical
            expression is so easily available.an ntg::Complex{Float64} is retured.
    """
    function envelopeQuadVolkovIntegral(envelope::Pulse.AbstractEnvelope, beam::AbstractBeam, polarization::Basics.AbstractPolarization,
                                        thetap::Float64, phip::Float64, energyp::Float64, initialEn::Float64)
        wa = 0. * im;
        # Collect parameters that are specific to a given pulse envelope
        if       typeof(envelope) == SinSquaredEnvelope
        elseif   typeof(envelope) == GaussianEnvelope
        end
        #
        # Determine first the integrant; the timeGrid must still be adapted to thos
        timeGrid = [0.1i for i = 1:10]
        for   t  in  timeGrid
            phase = Pulse.volkovPhase(t, envelope) ## , ...)
        end
        #
        error("Not yet implemented.")

        return( wa )
    end

end # module


text:
function Base.cov{T <: Model}(
    ::Type{T},
    inputs::Vector,
    weights::Vector,
)
    # Construct the negative log likelihood function and its gradient.
    loss = LogLoss()
    λ = 0.0
    nll, nll_gr! = make_closures(T, loss, λ, inputs, weights)

    # Find the MLE for the model.
    res = fit(T, inputs, weights)
    Θ = Optim.minimizer(res)

    # Compute the Hessian at the MLE.
    H = ForwardDiff.hessian(nll, Θ)

    # Construct the covariance matrix.
    return inv(H)
end

ext: jl


text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# DynamicBoundspODEIneq.jl
# See https://github.com/PSORLab/DynamicBoundspODEIneq.jl
#############################################################################
# src/contractors.jl
# Defines contractors (specifically the Ig contractor of Kai Shen and Joseph
# Scott, “Exploiting Nonlinear Invariants and Path Constraints to Achieve
# Tighter Reachable Set Enclosures Using Differential Inequalities”,
# Mathematics of Control, Signals, and Systems, In Press).
#############################################################################

const POLYHEDRON_A_TOL = 1e-4
const POLYHEDRON_WIDTH_TOL = 1e-12

function mid3_I(x::Interval{Float64}, y::Interval{Float64}, z::Interval{Float64})
    a = max(min(x,y), min(max(x,y),z))
    return a
end

"""
$(FUNCTIONNAME)

Contracts the state variable bounds based on the `PolyhedralConstraint`.
"""
function polyhedral_contact!(d::PolyhedralConstraint, Xin::Vector{Interval{Float64}},
                             Xp::Vector{Interval{Float64}}, nx::Int, nm::Int)
    X = Xin
    A = d.A
    b = d.b
    Aik = 0.0
    Aij = 0.0
    zL = Interval{Float64}(0.0)
    zU = Interval{Float64}(0.0)
    alphaL = Interval{Float64}(0.0)
    alphaU = Interval{Float64}(0.0)
    lambda = Interval{Float64}(0.0)
    gamma = Interval{Float64}(0.0)
    Xtemp = Interval{Float64}(0.0)
    for i = 1:nm
        alphaL = @interval(b[i])
        alphaU = @interval(b[i])
        for k = 1:nx
            Aik = A[i, k]
            zL = @interval(X[k].lo)
            zU = @interval(X[k].hi)
            alphaL -= max(Aik*zL, Aik*zU)
            alphaU -= min(Aik*zL, Aik*zU)
        end
        for j = 1:nx
            Aij = A[i,j]
            if abs(Aij) > POLYHEDRON_A_TOL
                zL = @interval(X[j].lo)
                zU = @interval(X[j].hi)
                alphaL += max(Aij*zL, Aij*zU)
                alphaU += min(Aij*zL, Aij*zU)
                lambda = min(alphaL/Aij, alphaU/Aij)
                gamma = max(alphaL/Aij, alphaU/Aij)
                zL = mid3_I(zL, zU, lambda)
                zU = mid3_I(zL, zU, gamma)
                Xtemp = @interval(zL, zU)
                if diam(Xtemp) > POLYHEDRON_WIDTH_TOL
                    X[j] = Xtemp
                end
                alphaL -= max(Aij*zL, Aij*zU)
                alphaU -= min(Aij*zL, Aij*zU)
            end
        end
    end
    Xin[:] = X[:]
    return
end

function polyhedral_contact!(d::Nothing, Xin::Vector{Interval{Float64}}, Xp::Vector{Interval{Float64}}, nx::Int, nm::Int)
    return
end


text:
module JuliaToJupyter

# Write your package code here.
using JSON
f = open("test\\Test_1.jl","r")
file = open("writing files\\Test_1.ipynb","w")
start = "{ \n \"cells\": [\n"
finish = "],\n"
write(file,start)
while !eof(f)
    s = readline(f)
    if (s[1]=='#')
        println("Its a comment")
        line = Dict(
            "cell_type"=>"markdown",
            "metadata"=>"{}",
            "source"=>"[$s]"
            )
        json_object = json(line,4)
        write(file,json_object)
        write(file,",\n")
    else
        println("Its a code line")
        line = Dict(
            "cell_type"=>"code",
            "execution_count"=>"null",
            "metadata"=>"{}",
            "outputs"=>"[]",
            "source"=>"[$s]"
        )
        json_object = json(line,2)
        write(file,json_object)
        write(file,",\n")
    end
end
write(file,finish)
metadata_template = open("templates\\metadata.json")
while !eof(metadata_template)
    s = readline(metadata_template)
    write(file,s)
    write(file,"\n")
end
write(file,"}")
close(f)
close(file)
end


text:
using Mocking

userEndpoint = "api/$curApiVersion/users/{1}"
configEndpoint = "api/$curApiVersion/users/{1}/config"

"""
$(SIGNATURES)
Add a user to Synchrony.
Return: Returns the created user.
"""
function addUser(user::UserRequest)::UserResponse
    config = getGraffConfig()
    if config == nothing
        error("Graff config is not set, please call setGraffConfig with a valid configuration.")
    end
    url = "$(config.apiEndpoint)/$(format(userEndpoint, user.id))"
    response = @mock _sendRestRequest(config, HTTP.post, url, data=JSON.json(user), debug=true)
    if(response.status != 200)
        error("Error creating user, received $(response.status) with body '$(String(response.body))'.")
    end
    return _unmarshallWithLinks(String(response.body), UserResponse)
end

"""
$(SIGNATURES)
Gets a user given the user ID.
Return: The user for the given user ID.
"""
function getUser(userId::String)::UserResponse
    config = getGraffConfig()
    if config == nothing
        error("Graff config is not set, please call setGraffConfig with a valid configuration.")
    end
    url = "$(config.apiEndpoint)/$(format(userEndpoint, userId))"
    response = @mock _sendRestRequest(config, HTTP.get, url)
    body = String(response.body)
    if(response.status != 200)
        error("Error getting user, received $(response.status) with body '$(body)'.")
    end
    return JSON2.read(body, UserResponse)
end

"""
$(SIGNATURES)
Update a user.
Return: The updated user from the service.
"""
function updateUser(user::UserRequest)::UserResponse
    config = getGraffConfig()
    if config == nothing
        error("Graff config is not set, please call setGraffConfig with a valid configuration.")
    end
    url = "$(config.apiEndpoint)/$(format(userEndpoint, user.id))"
    response = @mock _sendRestRequest(config, HTTP.put, url, data=JSON.json(user))
    if(response.status != 200)
        error("Error updating user, received $(response.status) with body '$(String(response.body))'.")
    end
    return _unmarshallWithLinks(String(response.body), UserResponse)
end

"""
$(SIGNATURES)
Delete a user given a user ID.
NOTE: All robots must be deleted first, the call will fail if robots are still associated to the user.
Return: The deleted user.
"""
function deleteUser(userId::String)::UserResponse
    config = getGraffConfig()
    if config == nothing
        error("Graff config is not set, please call setGraffConfig with a valid configuration.")
    end
    url = "$(config.apiEndpoint)/$(format(userEndpoint, userId))"
    response = @mock _sendRestRequest(config, HTTP.delete, url)
    if(response.status != 200)
        error("Error deleting user, received $(response.status) with body '$(String(response.body))'.")
    end
    return _unmarshallWithLinks(String(response.body), UserResponse)
end


text:
######################################################
##The ABC-PMC algorithms and functions they all share.
##See abc_pmc_interace.jl for interface to most of these
##(The exception is Algorithm 4 which doesn't require an interace).
######################################################

"
Perform a version of ABC-PMC proposed in the paper, Algorithm 4.
This chooses the distance function for iteration t at the end of iteration t-1.

Its arguments are as follows:

* `abcinput`
An `ABCInput` variable.

* `N`
Number of accepted particles in each iteration.

* `α`
A tuning parameter between 0 and 1 determining how fast the acceptance threshold is reduced. (In more detai, the acceptance threshold in iteration t is the α quantile of distances from particles which meet the acceptance criteria of the previous iteration.)

* `maxsims`
The algorithm terminates once this many simulations have been performed.

* `nsims_for_init`
How many simulations are stored to initialise the distance function (by default 10,000).

* `store_init` (optional)
Whether to store the parameters and simulations used to initialise the distance function in each iteration (useful to produce many of the paper's figures).
These are stored for every iteration even if `adaptive` is false.
By default false.

* `diag_perturb` (optional)
Whether perturbation should be based on a diagonalised estimate of the covariance.
By default this is false, giving the perturbation described in Section 2.2 of the paper.

* `silent` (optional)
When true no progress bar or messages are shown during execution.
By default false.

* `h1` (optional)
The acceptance threshold for the first iteration, by default `Inf` (accept everything).
If this is not `Inf`, then `abcinput.abcdist` will be used in the first iteration, and so must not require initialisation or be pre-initialised.

The output is a `ABCPMCOutput` object.
"
function abcPMC4(abcinput::ABCInput, N::Integer, α::Float64, maxsims::Integer, nsims_for_init=10000; store_init=false, diag_perturb=false, silent=false, h1=Inf)
    if !silent
        prog = Progress(maxsims, 1) ##Progress meter
    end
    k::Int = ceil(N*α)
    nparameters = length(abcinput.prior)
    itsdone = 0
    simsdone = 0
    firstit = true
    ##We record a sequence of distances and thresholds
    ##(all distances the same but we record a sequence for consistency with other algorithm)
    dists = ABCDistance[abcinput.abcdist]
    thresholds = Float64[h1]
    rejOutputs = ABCRejOutput[]
    cusims = Int[]
    ##Main loop
    while (simsdone < maxsims)
        samplefromprior = (firstit || thresholds[itsdone]==Inf)
        if !samplefromprior
            perturbdist = getperturbdist(curroutput, diag_perturb)
        end
        ##Initialise new reference table
        newparameters = Array{Float64}(nparameters, N)
        newsumstats = Array{Float64}(abcinput.nsumstats, N)
        newpriorweights = Array{Float64}(N)
        successes_thisit = 0
        nextparticle = 1
        ##Initialise storage of simulated parameter/summary pairs
        sumstats_forinit = Array{Float64}(abcinput.nsumstats, nsims_for_init)
        pars_forinit = Array{Float64}(nparameters, nsims_for_init)
        ##Loop to fill up new reference table
        while (nextparticle <= N && simsdone<maxsims)
            ##Sample parameters from importance density
            if samplefromprior
                proppars = rand(abcinput.prior)
            else
                proppars = rimportance(curroutput, perturbdist)
            end
            ##Calculate prior weight and reject if zero
            priorweight = pdf(abcinput.prior, proppars)
            if (priorweight == 0.0)
                continue
            end
            ##Draw summaries
            (success, propstats) = abcinput.sample_sumstats(proppars)
            simsdone += 1
            if !silent
                next!(prog)
            end
            if (!success)
                ##If rejection occurred during simulation
                continue
            end
            if (successes_thisit < nsims_for_init)
                successes_thisit += 1
                sumstats_forinit[:,successes_thisit] = propstats
                pars_forinit[:,successes_thisit] = proppars
            end
            ##Accept if all previous distances less than corresponding thresholds
            accept = propgood(propstats, dists, thresholds)
            if (accept)
                newparameters[:,nextparticle] = copy(proppars)
                newsumstats[:,nextparticle] = copy(propstats)
                newpriorweights[nextparticle] = priorweight
                nextparticle += 1
            end
        end
        ##Stop if not all sims required to continue have been done (because simsdone==maxsims)
        if nextparticle<=N
            continue
        end
        ##Update counters
        itsdone += 1
        push!(cusims, simsdone)
        ##Trim pars_forinit and sumstats_forinit to correct size
        if (successes_thisit < nsims_for_init)
            sumstats_forinit = sumstats_forinit[:,1:successes_thisit]
            pars_forinit = pars_forinit[:,1:successes_thisit]
        end

        currdist = dists[itsdone]
        if firstit && h1==Inf
            currdistances = zeros(N)
        else
            currdistances = Float64[ evaldist(currdist, newsumstats[:,i]) for i in 1:N ]
        end
        if samplefromprior
            newweights = ones(N)
        else
            oldoutput = copy(curroutput)
            newweights = getweights(newparameters, newpriorweights, oldoutput, perturbdist)
        end
        curroutput = ABCRejOutput(nparameters, abcinput.nsumstats, N, N, newparameters, newsumstats, currdistances, newweights, currdist, sumstats_forinit, pars_forinit) ##n.b. abcinput.abcdist may be uninitialised, but this is not problematic

        ##Create and store new distance
        newdist = init(abcinput.abcdist, sumstats_forinit, pars_forinit)
        push!(dists, newdist)

        ##Calculate and store threshold for next iteration
        newdistances = Float64[ evaldist(newdist, newsumstats[:,i]) for i in 1:N ]
        newthreshold = select!(newdistances, k)
        push!(thresholds, newthreshold)

        ##Record output
        push!(rejOutputs, curroutput)
        ##Report status
        if !silent
            print("\n Iteration $itsdone, $simsdone sims done\n")
            if firstit
                accrate = k/simsdone
            else
            accrate = k/(simsdone-cusims[itsdone-1])
            end
            @printf("Acceptance rate %.1e percent\n", 100*accrate)
            print("Output of most recent stage:\n")
            print(curroutput)
            print("Next threshold: $(convert(Float32, newthreshold))\n") ##Float64 shows too many significant figures
            ##TO DO: make some plots as well?
        end
        ##TO DO: consider alternative stopping conditions? (e.g. zero threshold reached)
        firstit = false
    end

    ##Put results into ABCPMCOutput object
    parameters = Array{Float64}(nparameters, N, itsdone)
    sumstats = Array{Float64}(abcinput.nsumstats, N, itsdone)
    distances = Array{Float64}(N, itsdone)
    weights = Array{Float64}(N, itsdone)
    for i in 1:itsdone
        parameters[:,:,i] = rejOutputs[i].parameters
        sumstats[:,:,i] = rejOutputs[i].sumstats
        distances[:,i] = rejOutputs[i].distances
        weights[:,i] = rejOutputs[i].weights
    end
    if (store_init)
        init_sims = Array{Array{Float64, 2}}(itsdone)
        init_pars = Array{Array{Float64, 2}}(itsdone)
        for i in 1:itsdone
            init_sims[i] = rejOutputs[i].init_sims
            init_pars[i] = rejOutputs[i].init_pars
        end
    else
        init_sims = Array{Array{Float64, 2}}(0)
        init_pars = Array{Array{Float64, 2}}(0)
    end
    output = ABCPMCOutput(nparameters, abcinput.nsumstats, itsdone, simsdone, cusims, parameters, sumstats, distances, weights, dists[1:itsdone], thresholds[1:itsdone], init_sims, init_pars)
end

"
Perform a version of ABC-PMC: either Algorithm 5 or a variant of Algorithm 3.
See `abcPMC3V` and `abcPMC5` documentation for details.
The `adaptive` argument determines which algorithm is used (`true` gives Algorithm 5).
"
function abcPMC_3V5(abcinput::ABCInput, N::Integer, α::Float64, maxsims::Integer, adaptive::Bool, nsims_for_init=10000; store_init=false, diag_perturb=false, silent=false)
    if !silent
        prog = Progress(maxsims, 1) ##Progress meter
    end
    M::Int = ceil(N/α)
    nparameters = length(abcinput.prior)
    itsdone = 0
    simsdone = 0
    firstit = true
    ##We record a sequence of distances and thresholds
    ##(for non-adaptive case all distances the same, and we only use most recent threshold)
    dists = ABCDistance[]
    thresholds = Float64[]
    rejOutputs = ABCRejOutput[]
    cusims = Int[]
    ##Main loop
    while (simsdone < maxsims)
        samplefromprior = firstit
        if !samplefromprior
            perturbdist = getperturbdist(curroutput, diag_perturb)
        end
        ##Initialise new reference table
        newparameters = Array{Float64}(nparameters, M)
        newsumstats = Array{Float64}(abcinput.nsumstats, M)
        newpriorweights = Array{Float64}(M)
        successes_thisit = 0
        if (firstit || adaptive || store_init)
            ##Initialise storage of simulated parameter/summary pairs
            sumstats_forinit = Array{Float64}(abcinput.nsumstats, nsims_for_init)
            pars_forinit = Array{Float64}(nparameters, nsims_for_init)
        end
        nextparticle = 1
        ##Loop to fill up new reference table
        while (nextparticle <= M && simsdone<maxsims)
            ##Sample parameters from importance density
            if (firstit)
                proppars = rand(abcinput.prior)
            else
                proppars = rimportance(curroutput, perturbdist)
            end
            ##Calculate prior weight and reject if zero
            priorweight = pdf(abcinput.prior, proppars)
            if (priorweight == 0.0)
                continue
            end
            ##Draw summaries
            (success, propstats) = abcinput.sample_sumstats(proppars)
            simsdone += 1
            if !silent
                next!(prog)
            end
            if (!success)
                ##If rejection occurred during simulation
                continue
            end
            if ((firstit || adaptive || store_init) && successes_thisit < nsims_for_init)
                successes_thisit += 1
                sumstats_forinit[:,successes_thisit] = propstats
                pars_forinit[:,successes_thisit] = proppars
            end
            if (firstit)
                ##No rejection at this stage in first iteration
                accept = true
            elseif (adaptive)
                ##Accept if all prev distances less than corresponding thresholds
                accept = propgood(propstats, dists, thresholds)
            else
                ##Accept if distance less than current threshold
                accept = propgood(propstats, dists[itsdone], thresholds[itsdone])
            end
            if (accept)
                newparameters[:,nextparticle] = copy(proppars)
                newsumstats[:,nextparticle] = copy(propstats)
                newpriorweights[nextparticle] = priorweight
                nextparticle += 1
            end
        end
        ##Stop if not all sims required to continue have been done (because simsdone==maxsims)
        if nextparticle<=M
            continue
        end
        ##Update counters
        itsdone += 1
        push!(cusims, simsdone)
        ##Trim pars_forinit and sumstats_forinit to correct size
        if (firstit || adaptive || store_init)
            if (successes_thisit < nsims_for_init)
                sumstats_forinit = sumstats_forinit[:,1:successes_thisit]
                pars_forinit = pars_forinit[:,1:successes_thisit]
            end
        else
            ##Create some empty arrays to use as arguments
            sumstats_forinit = Array{Float64}(0,0)
            pars_forinit = Array{Float64}(0,0)
        end
        ##Create new distance if needed
        if (firstit || adaptive)
            newdist = init(abcinput.abcdist, sumstats_forinit, pars_forinit)
        else
            newdist = dists[1]
        end
        push!(dists, newdist)

        ##Calculate distances
        distances = Float64[ evaldist(newdist, newsumstats[:,i]) for i=1:M ]
        if !samplefromprior
            oldoutput = copy(curroutput)
        end
        curroutput = ABCRejOutput(nparameters, abcinput.nsumstats, M, N, newparameters, newsumstats, distances, newpriorweights, newdist, sumstats_forinit, pars_forinit) ##Temporarily use prior weights
        sortABCOutput!(curroutput)
        ##Calculate, store and use new threshold
        newthreshold = curroutput.distances[N]
        push!(thresholds, newthreshold)
        curroutput.parameters = curroutput.parameters[:,1:N]
        curroutput.sumstats = curroutput.sumstats[:,1:N]
        curroutput.distances = curroutput.distances[1:N]
        if samplefromprior
            curroutput.weights = ones(N)
        else
            curroutput.weights = getweights(curroutput.parameters, curroutput.weights, oldoutput, perturbdist)
        end

        ##Record output
        push!(rejOutputs, curroutput)
        ##Report status
        if !silent
            print("\n Iteration $itsdone, $simsdone sims done\n")
            if firstit
                accrate = N/simsdone
            else
            accrate = N/(simsdone-cusims[itsdone-1])
            end
            @printf("Acceptance rate %.1e percent\n", 100*accrate)
            print("Output of most recent stage:\n")
            print(curroutput)
            ##TO DO: make some plots as well?
        end
        ##TO DO: consider alternative stopping conditions? (e.g. zero threshold reached)
        firstit = false
    end

    ##Put results into ABCPMCOutput object
    parameters = Array{Float64}(nparameters, N, itsdone)
    sumstats = Array{Float64}(abcinput.nsumstats, N, itsdone)
    distances = Array{Float64}(N, itsdone)
    weights = Array{Float64}(N, itsdone)
    for i in 1:itsdone
        parameters[:,:,i] = rejOutputs[i].parameters
        sumstats[:,:,i] = rejOutputs[i].sumstats
        distances[:,i] = rejOutputs[i].distances
        weights[:,i] = rejOutputs[i].weights
    end
    if (store_init)
        init_sims = Array{Array{Float64, 2}}(itsdone)
        init_pars = Array{Array{Float64, 2}}(itsdone)
        for i in 1:itsdone
            init_sims[i] = rejOutputs[i].init_sims
            init_pars[i] = rejOutputs[i].init_pars
        end
    else
        init_sims = Array{Array{Float64, 2}}(0)
        init_pars = Array{Array{Float64, 2}}(0)
    end
    output = ABCPMCOutput(nparameters, abcinput.nsumstats, itsdone, simsdone, cusims, parameters, sumstats, distances, weights, dists, thresholds, init_sims, init_pars)
end

"
Perform a version of ABC-PMC: either Algorithm 2 or 3.
See `abcPMC2` and `abcPMC3` documentation for details.
The `initialise_dist` argument specifies whether the distance will be initialised at the end of the first iteration.
If `true` the code implements Algorithm 3.
If `false` a distance must be specified which does not need initialisation and Algorithm 2 is run.
"
function abcPMC_23(abcinput::ABCInput, N::Integer, α::Float64, maxsims::Integer, initialise_dist::Bool, nsims_for_init=10000; store_init=false, diag_perturb=false, silent=false, h1=Inf)
    if initialise_dist && h1<Inf
        error("To initialise distance during the algorithm the first threshold must be Inf")
    end
    if !silent
        prog = Progress(maxsims, 1) ##Progress meter
    end
    k::Int = ceil(N*α)
    nparameters = length(abcinput.prior)
    itsdone = 0
    simsdone = 0
    firstit = true
    ##We record a sequence of distances and thresholds
    ##(all distances the same but we record a sequence for consistency with other algorithm)
    dists = ABCDistance[abcinput.abcdist]
    thresholds = Float64[h1]
    rejOutputs = ABCRejOutput[]
    cusims = Int[]
    ##Main loop
    while (simsdone < maxsims)
        samplefromprior = (firstit || thresholds[itsdone]==Inf)
        if !samplefromprior
            perturbdist = getperturbdist(curroutput, diag_perturb)
        end
        ##Initialise new reference table
        newparameters = Array{Float64}(nparameters, N)
        newsumstats = Array{Float64}(abcinput.nsumstats, N)
        newpriorweights = Array{Float64}(N)
        successes_thisit = 0
        if (firstit || store_init)
            ##Initialise storage of simulated parameter/summary pairs
            sumstats_forinit = Array{Float64}(abcinput.nsumstats, nsims_for_init)
            pars_forinit = Array{Float64}(nparameters, nsims_for_init)
        end
        nextparticle = 1
        ##Loop to fill up new reference table
        while (nextparticle <= N && simsdone<maxsims)
            ##Sample parameters from importance density
            if samplefromprior
                proppars = rand(abcinput.prior)
            else
                proppars = rimportance(curroutput, perturbdist)
            end
            ##Calculate prior weight and reject if zero
            priorweight = pdf(abcinput.prior, proppars)
            if (priorweight == 0.0)
                continue
            end
            ##Draw summaries
            (success, propstats) = abcinput.sample_sumstats(proppars)
            simsdone += 1
            if !silent
                next!(prog)
            end
            if (!success)
                ##If rejection occurred during simulation
                continue
            end
            if (((firstit && initialise_dist) || store_init) && successes_thisit < nsims_for_init)
                successes_thisit += 1
                sumstats_forinit[:,successes_thisit] = propstats
                pars_forinit[:,successes_thisit] = proppars
            end
            if (firstit && initialise_dist)
                ##No rejection at this stage in first iteration if we want to initialise distance
                accept = true
            else
                ##Accept if distance less than current threshold
                accept = propgood(propstats, dists[1], thresholds[itsdone+1])
            end
            if (accept)
                newparameters[:,nextparticle] = copy(proppars)
                newsumstats[:,nextparticle] = copy(propstats)
                newpriorweights[nextparticle] = priorweight
                nextparticle += 1
            end
        end
        ##Stop if not all sims required to continue have been done (because simsdone==maxsims)
        if nextparticle<=N
            continue
        end
        ##Update counters
        itsdone += 1
        push!(cusims, simsdone)
        ##Trim pars_forinit and sumstats_forinit to correct size
        if ((firstit && initialise_dist) || store_init)
            if (successes_thisit < nsims_for_init)
                sumstats_forinit = sumstats_forinit[:,1:successes_thisit]
                pars_forinit = pars_forinit[:,1:successes_thisit]
            end
        end
        ##Create new distance if needed
        if (firstit && initialise_dist)
            newdist = init(dists[1], sumstats_forinit, pars_forinit)
        else
            newdist = dists[1]
        end
        ##Store new distance
        if firstit
            dists[1] = newdist
        else
            push!(dists, newdist)
        end

        ##Calculate distances
        distances = Float64[ evaldist(newdist, newsumstats[:,i]) for i=1:N ]
        if !samplefromprior
            oldoutput = copy(curroutput)
        end
        curroutput = ABCRejOutput(nparameters, abcinput.nsumstats, N, N, newparameters, newsumstats, distances, newpriorweights, newdist, sumstats_forinit, pars_forinit) ##Temporarily use prior weights
        sortABCOutput!(curroutput)
        ##Calculate and store threshold for next iteration
        newthreshold = curroutput.distances[k]
        push!(thresholds, newthreshold)
        if samplefromprior
            curroutput.weights = ones(N)
        else
            curroutput.weights = getweights(curroutput.parameters, curroutput.weights, oldoutput, perturbdist)
        end

        ##Record output
        push!(rejOutputs, curroutput)
        ##Report status
        if !silent
            print("\n Iteration $itsdone, $simsdone sims done\n")
            if firstit
                accrate = k/simsdone
            else
            accrate = k/(simsdone-cusims[itsdone-1])
            end
            @printf("Acceptance rate %.1e percent\n", 100*accrate)
            print("Output of most recent stage:\n")
            print(curroutput)
            print("Next threshold: $(convert(Float32, newthreshold))\n") ##Float64 shows too many significant figures
            ##TO DO: make some plots as well?
        end
        ##TO DO: consider alternative stopping conditions? (e.g. zero threshold reached)
        firstit = false
    end

    ##Put results into ABCPMCOutput object
    parameters = Array{Float64}(nparameters, N, itsdone)
    sumstats = Array{Float64}(abcinput.nsumstats, N, itsdone)
    distances = Array{Float64}(N, itsdone)
    weights = Array{Float64}(N, itsdone)
    for i in 1:itsdone
        parameters[:,:,i] = rejOutputs[i].parameters
        sumstats[:,:,i] = rejOutputs[i].sumstats
        distances[:,i] = rejOutputs[i].distances
        weights[:,i] = rejOutputs[i].weights
    end
    if (store_init)
        init_sims = Array{Array{Float64, 2}}(itsdone)
        init_pars = Array{Array{Float64, 2}}(itsdone)
        for i in 1:itsdone
            init_sims[i] = rejOutputs[i].init_sims
            init_pars[i] = rejOutputs[i].init_pars
        end
    else
        init_sims = Array{Array{Float64, 2}}(0)
        init_pars = Array{Array{Float64, 2}}(0)
    end
    output = ABCPMCOutput(nparameters, abcinput.nsumstats, itsdone, simsdone, cusims, parameters, sumstats, distances, weights, dists, thresholds, init_sims, init_pars)
end

"
Return the perturbation distribution to be used given output `x` from the previous iteration.
Argument `diag` specifies whether the perturbation should have a diagonal covariance.
"
function getperturbdist(x::ABCRejOutput, diag::Bool)
    wv = Weights(x.weights)
    if (diag)
        ##Calculate diagonalised variance of current weighted particle approximation
        diagvar = Float64[var(x.parameters[i,:], wv, corrected=false) for i in 1:x.nparameters]
        perturbdist = MvNormal(2.0 .* diagvar)
    else
        ##Calculate variance of current weighted particle approximation
        currvar = cov(x.parameters, wv, 2, corrected=false)
        perturbdist = MvNormal(2.0 .* currvar)
    end
    return perturbdist
end

##Check if summary statistics meet acceptance requirement
function propgood(s::Array{Float64, 1}, dist::ABCDistance, threshold::Float64)
    ##n.b. threshold==Inf gives acceptance even if dist is not initialised
    threshold == Inf || evaldist(dist, s)<=threshold
end

##Check if summary statistics meet all of previous acceptance requirements
function propgood(s::Array{Float64, 1}, dists::Array{ABCDistance, 1}, thresholds::Array{Float64, 1})
    for i in length(dists):-1:1 ##Check the most stringent case first
        if !propgood(s, dists[i], thresholds[i])
            return false
        end
    end
    return true
end

##Samples from importance density defined by prev output
function rimportance(out::ABCRejOutput, dist::MvNormal)
    i = sample(Weights(out.weights))
    out.parameters[:,i] + rand(dist)
end

"
Calculate importance weight of parameters `x` given `priorweight`, their prior density, `old`, the output of the previous iteration and `perturbdist`, the perturbation distribution used in the importance density.
"
function get1weight(x::Array{Float64,1}, priorweight::Float64, old::ABCRejOutput, perturbdist::MvNormal)
    nparticles = size(old.parameters)[2]
    temp = Float64[pdf(perturbdist, x-old.parameters[:,i]) for i in 1:nparticles]
    priorweight / sum(old.weights .* temp)
end

"
Calculates importance weights for parameters `pars` (columns are parameter vectors) given `priorweights`, their priors densities, `old`, the output of the previous iteration and `perturbdist`, the perturbation distribution used in the importance density.
"
function getweights(pars::Array{Float64, 2}, priorweights::Array{Float64,1}, old::ABCRejOutput, perturbdist::MvNormal)
    nparticles = size(pars)[2]
    weights = Float64[get1weight(pars[:,i], priorweights[i], old, perturbdist) for i in 1:nparticles]
    weights ./ sum(weights)
end


text:
#################################################################
# This module provides implementations of MDP and POMDP problems
# in the POMDPs.jl format. These implementations serve as
# examples for new users and as benchmark problems for testing
# and valdiation.
#################################################################

module POMDPModels

using POMDPs
using POMDPModelTools
using BeliefUpdaters
using Distributions
using StaticArrays
using StatsBase
using Random
using Printf
using Parameters

# for visualization
using Compose
using ColorSchemes
import POMDPModelTools.render

using POMDPs

import Base: ==, hash
import Random: rand, rand!
import Distributions: pdf

import POMDPs: gen, support, discount, isterminal
import POMDPs: actions, actionindex, action, dimensions
import POMDPs: states, stateindex, transition
import POMDPs: observations, observation, obsindex
import POMDPs: initialstate, initialstate_distribution, initialobs
import POMDPs: updater, update
import POMDPs: reward
import POMDPs: convert_s, convert_a, convert_o

include("TigerPOMDPs.jl")
export
    TigerPOMDP,
    TigerDistribution,
    TigerStateSpace,
    TigerActionSpace,
    TigerObservationSpace,
    TigerBeliefUpdater,
    TIGER_LEFT,
    TIGER_RIGHT,
    TIGER_LISTEN,
    TIGER_OPEN_LEFT,
    TIGER_OPEN_RIGHT

include("gridworld.jl")
include("gridworld_visualization.jl")
export
    GWPos,
    SimpleGridWorld

include("CryingBabies.jl")
export
    BabyPOMDP,
    BabyBeliefUpdater,
    Starve,
    AlwaysFeed,
    FeedWhenCrying

include("MountainCar.jl")
export
    MountainCar,
    Energize

include("InvertedPendulum.jl")
export
    InvertedPendulum

include("Tabular.jl")
export
    TabularMDP,
    TabularPOMDP

include("Random.jl")
export
    RandomMDP,
    RandomPOMDP

include("TMazes.jl")
export
    TMaze,
    TMazeState,
    TMazeSpace,
    TMazeStateSpace,
    MazeOptimal

include("LightDark.jl")
export
    LightDark1D,
    LightDark1DState,
    LightDark1DActionSpace,
    DummyHeuristic1DPolicy,
    SmartHeuristic1DPolicy

# Legacy

include("legacy/GridWorlds.jl")
export
    LegacyGridWorld,
    GridWorldState,
    GridWorldAction,
    GridWorldActionSpace,
    GridWorldStateSpace,
    GridWorldDistribution,
    static_reward

end # module


text:
# compatibility of Julia, CUDA and LLVM

const lowest = v"0"
const highest = v"999"

verlist(vers) = join(map(ver->"$(ver.major).$(ver.minor)", sort(collect(vers))), ", ", " and ")

## version range

struct VersionRange
    lower::VersionNumber
    upper::VersionNumber
end

Base.in(v::VersionNumber, r::VersionRange) = (v >= r.lower && v <= r.upper)

import Base.(:)
(:)(a::VersionNumber, b::VersionNumber) = VersionRange(a, b)

Base.intersect(v::VersionNumber, r::VersionRange) =
    v < r.lower ? (r.lower:v) :
    v > r.upper ? (v:r.upper) : (v:v)

## devices supported by the CUDA toolkit

# Source:
# - https://en.wikipedia.org/wiki/CUDA#GPUs_supported
# - ptxas |& grep -A 10 '\--gpu-name'
const cuda_cap_db = Dict(
    v"1.0" => lowest:v"6.5",
    v"1.1" => lowest:v"6.5",
    v"1.2" => lowest:v"6.5",
    v"1.3" => lowest:v"6.5",
    v"2.0" => lowest:v"8.0",
    v"2.1" => lowest:v"8.0",
    v"3.0" => v"4.2":v"10.2",
    v"3.2" => v"6.0":v"10.2",
    v"3.5" => v"5.0":highest,   # deprecated with v11
    v"3.7" => v"6.5":highest,   # deprecated with v11
    v"5.0" => v"6.0":highest,   # deprecated with v11
    v"5.2" => v"7.0":highest,
    v"5.3" => v"7.5":highest,
    v"6.0" => v"8.0":highest,
    v"6.1" => v"8.0":highest,
    v"6.2" => v"8.0":highest,
    v"7.0" => v"9.0":highest,
    v"7.2" => v"9.2":highest,
    v"7.5" => v"10.0":highest,
    v"8.0" => v"11.0":highest,
    v"8.6" => v"11.1":highest,
)

function cuda_cap_support(ver::VersionNumber)
    caps = Set{VersionNumber}()
    for (cap,r) in cuda_cap_db
        if ver in r
            push!(caps, cap)
        end
    end
    return caps
end

## PTX ISAs supported by the CUDA toolkit

# Source:
# - PTX ISA document, Release History table
const cuda_ptx_db = Dict(
    v"1.0" => v"1.0":highest,
    v"1.1" => v"1.1":highest,
    v"1.2" => v"2.0":highest,
    v"1.3" => v"2.1":highest,
    v"1.4" => v"2.2":highest,
    v"1.5" => v"2.2":highest,
    v"2.0" => v"3.0":highest,
    v"2.1" => v"3.1":highest,
    v"2.2" => v"3.2":highest,
    v"2.3" => v"4.2":highest,
    v"3.0" => v"4.1":highest,
    v"3.1" => v"5.0":highest,
    v"3.2" => v"5.5":highest,
    v"4.0" => v"6.0":highest,
    v"4.1" => v"6.5":highest,
    v"4.2" => v"7.0":highest,
    v"4.3" => v"7.5":highest,
    v"5.0" => v"8.0":highest,
    v"6.0" => v"9.0":highest,
    v"6.1" => v"9.1":highest,
    v"6.2" => v"9.2":highest,
    v"6.3" => v"10.0":highest,
    v"6.4" => v"10.1":highest,
    v"6.5" => v"10.2":highest,
    v"7.0" => v"11.0":highest,
    v"7.1" => v"11.1":highest,
    v"7.2" => v"11.2":highest,
)

function cuda_ptx_support(ver::VersionNumber)
    caps = Set{VersionNumber}()
    for (cap,r) in cuda_ptx_db
        if ver in r
            push!(caps, cap)
        end
    end
    return caps
end

## devices supported by the LLVM NVPTX back-end

# Source: LLVM/lib/Target/NVPTX/NVPTX.td
const llvm_cap_db = Dict(
    v"2.0" => v"3.2":highest,
    v"2.1" => v"3.2":highest,
    v"3.0" => v"3.2":highest,
    v"3.2" => v"3.7":highest,
    v"3.5" => v"3.2":highest,
    v"3.7" => v"3.7":highest,
    v"5.0" => v"3.5":highest,
    v"5.2" => v"3.7":highest,
    v"5.3" => v"3.7":highest,
    v"6.0" => v"3.9":highest,
    v"6.1" => v"3.9":highest,
    v"6.2" => v"3.9":highest,
    v"7.0" => v"6.0":highest,
    v"7.2" => v"7.0":highest,
    v"7.5" => v"8.0":highest,
    v"8.0" => v"11.0":highest,
)

function llvm_cap_support(ver::VersionNumber)
    caps = Set{VersionNumber}()
    for (cap,r) in llvm_cap_db
        if ver in r
            push!(caps, cap)
        end
    end
    return caps
end

## PTX ISAs supported by the LVM NVPTX back-end

# Source: LLVM/lib/Target/NVPTX/NVPTX.td
const llvm_ptx_db = Dict(
    v"3.0" => v"3.2":v"3.5",
    v"3.1" => v"3.2":v"3.5",
    v"3.2" => v"3.5":highest,
    v"4.0" => v"3.5":highest,
    v"4.1" => v"3.7":highest,
    v"4.2" => v"3.7":highest,
    v"4.3" => v"3.9":highest,
    v"5.0" => v"3.9":highest,
    v"6.0" => v"6.0":highest,
    v"6.1" => v"7.0":highest,
    v"6.3" => v"8.0":highest,
    v"6.4" => v"9.0":highest,
    v"6.5" => v"11.0":highest,
    v"7.0" => v"11.0":highest,
)

function llvm_ptx_support(ver::VersionNumber)
    caps = Set{VersionNumber}()
    for (cap,r) in llvm_ptx_db
        if ver in r
            push!(caps, cap)
        end
    end
    return caps
end

## high-level functions that return target and isa support

function llvm_compat(version=LLVM.version())
    InitializeNVPTXTarget()

    cap_support = sort(collect(llvm_cap_support(version)))
    ptx_support = sort(collect(llvm_ptx_support(version)))

    return (cap=cap_support, ptx=ptx_support)
end

function cuda_compat(driver_release=release(), toolkit_release=toolkit_release())
    driver_cap_support = cuda_cap_support(driver_release)
    toolkit_cap_support = cuda_cap_support(toolkit_release)
    cap_support = sort(collect(driver_cap_support ∩ toolkit_cap_support))

    driver_ptx_support = cuda_ptx_support(driver_release)
    toolkit_ptx_support = cuda_ptx_support(toolkit_release)
    ptx_support = sort(collect(driver_ptx_support ∩ toolkit_ptx_support))

    return (cap=cap_support, ptx=ptx_support)
end

# select the highest capability that is supported by both the toolchain and device
@memoize function supported_capability(dev::CuDevice)
    dev_cap = capability(dev)
    compat_caps = filter(cap -> cap <= dev_cap, target_support())
    isempty(compat_caps) &&
        error("Device capability v$dev_cap not supported by available toolchain")

    return maximum(compat_caps)
end

## initialization

const __target_support = Ref{Vector{VersionNumber}}()
const __ptx_support = Ref{Vector{VersionNumber}}()

target_support() = @after_init(__target_support[])
ptx_support() = @after_init(__ptx_support[])

function __init_compatibility__()
    llvm_support = llvm_compat()
    cuda_support = cuda_compat()

    __target_support[] = sort(collect(llvm_support.cap ∩ cuda_support.cap))
    isempty(__target_support[]) && error("Your toolchain does not support any device capability")

    __ptx_support[] = sort(collect(llvm_support.ptx ∩ cuda_support.ptx))
    isempty(__ptx_support[]) && error("Your toolchain does not support any PTX ISA")

    @debug("Toolchain with LLVM $(LLVM.version()), CUDA driver $(version()) and toolkit $(toolkit_version()) supports devices $(verlist(__target_support[])); PTX $(verlist(__ptx_support[]))")
end


text:
"""
    cv_iht(d, l, x, z, y, J, path, q)
For each model specified in `path`, performs `q`-fold cross validation and
returns the (averaged) deviance residuals.
The purpose of this function is to find the best sparsity level `k`, judiciously obtained
from selecting the model with the minimum out-of-sample error. Automatically finds the
correct version of `L0_reg` to use depending on the type of `x`. By default, each CPU runs
a different model for a given fold. To use this function, start julia using 4 processors
(the more the better) by:
    julia> using Distributed
    julia> addprocs(4)
# Warning
Do not remove files with random file names when you run this function. There are
memory mapped files that will be deleted automatically once they are no longer needed.
# Arguments
- `d`: A distribution (e.g. Normal, Bernoulli)
- `l`: A link function (e.g. Loglink, ProbitLink)
- `x`: A SnpArray, which can be memory mapped to a file. Does not engage in any linear algebra
- `z`: Matrix of non-genetic covariates. The first column usually denotes the intercept.
- `y`: Response vector
- `J`: The number of maximum groups
- `path`: Vector storing different model sizes
- `q`: Number of cross validation folds. For large data do not set this to be greater than 5
# Optional Arguments:
- `group` vector storing group membership
- `weight` vector storing vector of weights containing prior knowledge on each SNP
- `folds`: Vector that separates the sample into q disjoint subsets
- `destin`: Directory where intermediate files will be generated. Directory name must end with `/`.
- `init`: Boolean indicating whether we should initialize IHT algorithm at a good starting guess
- `use_maf`: Boolean indicating we should scale the projection step by a weight vector
- `debias`: Boolean indicating whether we should debias at each IHT step
- `showinfo`: Whether we want IHT to print meaningful intermediate steps
- `parallel`: Whether we want to run cv_iht using multiple CPUs (highly recommended)
"""
function cv_iht(
    d        :: UnivariateDistribution,
    l        :: Link,
    x        :: Union{AbstractMatrix{T}, SnpArray},
    z        :: AbstractMatrix{T},
    y        :: AbstractVector{T},
    J        :: Int64,
    path     :: DenseVector{Int},
    q        :: Int64;
    group    :: AbstractVector{Int} = Int[],
    weight   :: AbstractVector{T} = T[],
    folds    :: DenseVector{Int} = rand(1:q, size(x, 1)),
    destin   :: String = "./",
    init     :: Bool = false,
    use_maf  :: Bool = false,
    debias   :: Bool = false,
    showinfo :: Bool = true,
    parallel :: Bool = false
) where {T <: Float}

    # preallocate mean squared error matrix
    nmodels = length(path)
    mses = zeros(nmodels, q)

    for fold in 1:q
        # find entries that are for test sets and train sets
        test_idx  = folds .== fold
        train_idx = .!test_idx

        # validate trained models on test data by computing deviance residuals
        mses[:, fold] .= train_and_validate(train_idx, test_idx, d, l, x, z, y, J, path, fold, group=group, weight=weight, destin=destin, init=init, use_maf=use_maf, debias=debias, showinfo=false, parallel=parallel)
    end

    #weight mses for each fold by their size before averaging
    mse = meanloss(mses, q, folds)

    # find best model size and print cross validation result
    k = path[argmin(mse)] :: Int
    showinfo && print_cv_results(mse, path, k)

    return mse
end

"""
Performs q-fold cross validation for Iterative hard thresholding to
determine the best model size `k`. The function is the same as `cv_iht`
except here each `fold` is distributed to a different CPU as opposed
to each `path` to a different CPU.
"""
function cv_iht_distribute_fold(
    d        :: UnivariateDistribution,
    l        :: Link,
    x        :: Union{AbstractMatrix{T}, SnpArray},
    z        :: AbstractMatrix{T},
    y        :: AbstractVector{T},
    J        :: Int64,
    path     :: DenseVector{Int},
    q        :: Int64;
    group    :: AbstractVector{Int} = Int[],
    weight   :: AbstractVector{T} = T[],
    folds    :: DenseVector{Int} = rand(1:q, size(x, 1)),
    destin   :: String = "./",
    init     :: Bool = false,
    use_maf  :: Bool = false,
    debias   :: Bool = false,
    showinfo :: Bool = true,
    parallel :: Bool = false
) where {T <: Float}

    # for each fold, allocate train/test set, train the model, and test the model
    mses = (parallel ? pmap : map)(1:q) do fold
        test_idx  = folds .== fold
        train_idx = .!test_idx
        betas, cs = pfold_train(train_idx, x, z, y, J, d, l, path, fold, group=group, weight=weight, destin=destin, init=init, use_maf=use_maf, debias=debias, showinfo=false)
        return pfold_validate(test_idx, betas, cs, x, z, y, J, d, l, path, fold, group=group, weight=weight, destin=destin, init=init, use_maf=use_maf, debias=debias, showinfo=false)
    end

    #weight mses for each fold by their size before averaging
    mse = meanloss(mses, q, folds)

    # find best model size and print cross validation result
    k = path[argmin(mse)] :: Int
    showinfo && print_cv_results(mse, path, k)

    return mse
end

"""
Runs IHT across many different model sizes specifed in `path`.
This is basically the same as `cv_iht` except we **DO NOT** validate each model
in a holdout set, meaning that this will definitely induce overfitting as we increase
model size. Use this to perform a quick estimate a range of feasible model sizes before
engaging in full cross validation.
"""
function iht_run_many_models(
    d        :: UnivariateDistribution,
    l        :: Link,
    x        :: Union{AbstractMatrix{T}, SnpArray},
    z        :: AbstractMatrix{T},
    y        :: AbstractVector{T},
    J        :: Int64,
    path     :: DenseVector{Int};
    group    :: AbstractVector{Int} = Int[],
    weight   :: AbstractVector{T} = T[],
    init     :: Bool = false,
    use_maf  :: Bool = false,
    debias   :: Bool = false,
    showinfo :: Bool = true,
    parallel :: Bool = false
) where {T <: Float}

    # for each k, run L0_reg and store the loglikelihoods
    typeof(x) == SnpArray && (xbm = SnpBitMatrix{T}(x, model=ADDITIVE_MODEL, center=true, scale=true);)
    results = (parallel ? pmap : map)(path) do k
        if typeof(x) == SnpArray
            return L0_reg(x, xbm, z, y, 1, k, d, l, group=group, weight=weight, init=init, use_maf=use_maf, debias=debias, show_info=false)
        else
            return L0_reg(x, z, y, 1, k, d, l, group=group, weight=weight, init=init, use_maf=use_maf, debias=debias, show_info=false)
        end
    end

    loglikelihoods = zeros(size(path, 1))
    for i in 1:length(results)
        loglikelihoods[i] = results[i].logl
    end

    #display result and then return
    showinfo && print_a_bunch_of_path_results(loglikelihoods, path)
    return loglikelihoods
end

"""
This function trains a bunch of models, where each model has a different sparsity
parameter, k, which is specified in the variable `path`. Then each trained model is used to
compute the deviance residuals (i.e. mean squared error for normal response) on the test set.
on the test set. This deviance residuals vector is returned
"""
function train_and_validate(train_idx::BitArray, test_idx::BitArray, d::UnivariateDistribution,
                    l::Link, x::SnpArray, z::AbstractMatrix{T}, y::AbstractVector{T}, J::Int64,
                    path::DenseVector{Int}, fold::Int; group::AbstractVector{Int}=Int[],
                    weight::AbstractVector{T}=T[], init::Bool=false, destin::String = "./",
                    use_maf::Bool=false, debias::Bool=false, showinfo::Bool=true,
                    parallel::Bool=false) where {T <: Float}

    # create directory for memory mapping
    train_file = destin * randstring(100) * ".bed"
    test_file  = destin * randstring(100) * ".bed"

    # first allocate arrays needed for computing deviance residuals
    p, q = size(x, 2), size(z, 2)
    nmodels = length(path)
    test_size = sum(test_idx)
    xb = zeros(T, test_size)
    zc = zeros(T, test_size)
    μ  = zeros(T, test_size)

    # allocate train model
    x_train = SnpArray(train_file, sum(train_idx), p)
    copyto!(x_train, @view(x[train_idx, :]))
    x_trainbm = SnpBitMatrix{T}(x_train, model=ADDITIVE_MODEL, center=true, scale=true);
    z_train = z[train_idx, :]
    y_train = y[train_idx]

    # allocate test model
    x_test = SnpArray(test_file, test_size, p)
    copyto!(x_test, @view(x[test_idx, :]))
    x_testbm = SnpBitMatrix{T}(x_test, model=ADDITIVE_MODEL, center=true, scale=true);
    z_test = z[test_idx, :]
    y_test = y[test_idx]

    # allocate group and weight vectors if supplied
    group_train = (group == Int[] ? Int[] : group[train_idx])
    weight_train = (weight == T[] ? T[] : weight[train_idx])

    # for each k in path, run L0_reg and compute mse
    try
        mses = (parallel ? pmap : map)(path) do k

            #run IHT on training model with given k
            result = L0_reg(x_train, x_trainbm, z_train, y_train, 1, k, d, l, group=group_train, weight=weight_train, init=init, use_maf=use_maf, debias=debias, show_info=showinfo)

            # compute estimated response Xb: [xb zc] = [x_test z_test] * [b; c] and update mean μ = g^{-1}(xb)
            A_mul_B!(xb, zc, x_testbm, z_test, result.beta, result.c)
            update_μ!(μ, xb .+ zc, l)

            # compute sum of squared deviance residuals. For normal, this is equivalent to out-of-sample error
            return deviance(d, y_test, μ)
        end
    finally
        #clean up
        rm(train_file, force=true)
        rm(test_file, force=true)
    end

    return mses
end

function train_and_validate(train_idx::BitArray, test_idx::BitArray, d::UnivariateDistribution,
                    l::Link, x::AbstractMatrix{T}, z::AbstractMatrix{T}, y::AbstractVector{T},
                    J::Int64, path::DenseVector{Int}, fold::Int; group::AbstractVector{Int}=Int[],
                    weight::AbstractVector{T}=T[], destin::String = "./", init::Bool=false,
                    use_maf::Bool=false, debias::Bool=false, showinfo::Bool=true,
                    parallel::Bool=false) where {T <: Float}

    # first allocate arrays needed for computing deviance residuals
    p, q = size(x, 2), size(z, 2)
    nmodels = length(path)
    test_size = sum(test_idx)
    xb = zeros(T, test_size)
    zc = zeros(T, test_size)
    μ  = zeros(T, test_size)

    # allocate train model and test models
    x_train = x[train_idx, :]
    x_test = x[test_idx, :]

    # allocate group and weight vectors if supplied
    group_train = (group == Int[] ? Int[] : group[train_idx])
    weight_train = (weight == T[] ? T[] : weight[train_idx])

    # for each k in path, run L0_reg and compute mse
    mses = (parallel ? pmap : map)(path) do k

        #run IHT on training model with given k
        result = L0_reg(x_train, z[train_idx, :], y[train_idx], 1, k, d, l, group=group_train, weight=weight_train, init=init, use_maf=use_maf, debias=debias, show_info=showinfo)

        # compute estimated response Xb: [xb zc] = [x_test z_test] * [b; c] and update mean μ = g^{-1}(xb)
        A_mul_B!(xb, zc, x_test, z[test_idx, :], result.beta, result.c)
        update_μ!(μ, xb .+ zc, l)

        # compute sum of squared deviance residuals. For normal, this is equivalent to out-of-sample error
        return deviance(d, y[test_idx], μ)
    end

    return mses
end

function pfold_train(train_idx::BitArray, x::SnpArray, z::AbstractMatrix{T},
    y::AbstractVector{T}, J::Int64, d::UnivariateDistribution, l::Link,
    path::DenseVector{Int}, fold::Int64; group::AbstractVector{Int}=Int[],
    weight::AbstractVector{T}=T[], destin::String = "./", init::Bool=false,
    use_maf::Bool =false, max_iter::Int = 100, max_step::Int = 3,
    debias::Bool = false, showinfo::Bool = false) where {T <: Float}

    # create directory for memory mapping
    train_file = destin * randstring(100) * ".bed"

    #preallocate arrays
    p, q = size(x, 2), size(z, 2)
    nmodels = length(path)
    betas = zeros(T, p, nmodels)
    cs = zeros(T, q, nmodels)

    # allocate training data
    x_train = SnpArray(train_file, sum(train_idx), p)
    copyto!(x_train, view(x, train_idx, :))
    x_trainbm = SnpBitMatrix{T}(x_train, model=ADDITIVE_MODEL, center=true, scale=true);

    try
        for i in 1:length(path)
            k = path[i]
            result = L0_reg(x_train, x_trainbm, z[train_idx, :], y[train_idx], 1, k, d, l, group=group, weight=weight, init=init, use_maf=use_maf, debias=debias, show_info=false)
            betas[:, i] .= result.beta
            cs[:, i] .= result.c
        end
    finally
        rm(train_file, force=true) #clean up
    end

    return betas, cs
end

"""
This function takes a trained model, and returns the mean squared error (mse) of that model
on the test set. A vector of mse is returned, where each entry corresponds to the training
set on each fold with different sparsity parameter.
"""
function pfold_validate(test_idx::BitArray, betas::AbstractMatrix{T}, cs::AbstractMatrix{T},
    x::SnpArray, z::AbstractMatrix{T}, y::AbstractVector{T}, J::Int64, d::UnivariateDistribution,
    l::Link, path::DenseVector{Int}, fold::Int64; group::AbstractVector{Int}=Int[],
    weight::AbstractVector{T}=T[], destin::String = "./", init::Bool=false, use_maf::Bool = false,
    max_iter::Int = 100, max_step::Int = 3, debias::Bool = false, showinfo::Bool = false
    ) where {T <: Float}

    # create directory for memory mapping
    test_file = destin * randstring(100) * ".bed"

    # preallocate arrays
    p, q = size(x, 2), size(z, 2)
    test_size = sum(test_idx)
    mse = zeros(T, length(path))
    xb = zeros(T, test_size)
    zc = zeros(T, test_size)
    μ  = zeros(T, test_size)

    # allocate test model
    x_test = SnpArray(test_file, sum(test_idx), p)
    y_test = y[test_idx]
    copyto!(x_test, @view(x[test_idx, :]))
    x_testbm = SnpBitMatrix{T}(x_test, model=ADDITIVE_MODEL, center=true, scale=true);

    # for each computed model stored in betas, compute the deviance residuals (i.e. generalized mean squared error) on test set
    try
        for i = 1:size(betas, 2)
            # compute estimated response Xb: [xb zc] = [x_test z_test] * [b; c] and update mean μ = g^{-1}(xb)
            A_mul_B!(xb, zc, x_testbm, z[test_idx, :], @view(betas[:, i]), @view(cs[:, i]))
            update_μ!(μ, xb .+ zc, l)

            # compute sum of squared deviance residuals. For normal, this is equivalent to out-of-sample error
            mse[i] = deviance(d, y_test, μ)
        end
    finally
        rm(test_file, force=true) #clean up
    end

    return mse
end

"""
This function scale mean squared errors (deviance residuals) for each fold by its own fold size.
"""
function meanloss(fitloss::Matrix{T}, q::Int64,
                  folds::DenseVector{Int}) where {T <: Float}
    ninfold = zeros(Int, q)
    for fold in folds
        ninfold[fold] += 1
    end

    loss = zeros(T, size(fitloss, 1))
    for j = 1:size(fitloss, 2)
        wfold = convert(T, ninfold[j]/length(folds))
        for i = 1:size(fitloss, 1)
            loss[i] += fitloss[i, j]*wfold
        end
    end

    return loss
end

function meanloss(mses::Vector{Vector{T}}, num_fold::Int64,
                  folds::DenseVector{Int}) where {T <: Float}

    fitloss = hcat(mses...) :: Matrix{T}
    ninfold = zeros(Int, num_fold)
    for fold in folds
        ninfold[fold] += 1
    end

    loss = zeros(size(fitloss, 1))
    for j = 1:size(fitloss, 2)
        wfold = convert(T, ninfold[j]/length(folds))
        for i = 1:size(fitloss, 1)
            loss[i] += fitloss[i, j]*wfold
        end
    end

    return loss
end


text:
### Manual MuseProblem

struct MuseProblem{X,S,Gθ,Pθ,GZ} <: AbstractMuseProblem
    x :: X
    sample_x_z :: S
    ∇θ_logLike :: Gθ
    logLike_and_∇z_logLike :: GZ
    logPriorθ :: Pθ
end

function MuseProblem(x, sample_x_z, logLike, logPriorθ=(θ->0), ad::AD.AbstractBackend=AD.ForwardDiffBackend())
    MuseProblem(
        x,
        sample_x_z,
        (x,z,θ) -> first(AD.gradient(ad, θ -> logLike(x,z,θ), θ)),
        (x,z,θ) -> first.(AD.value_and_gradient(ad, z -> logLike(x,z,θ), z)),
        logPriorθ
    )
end

∇θ_logLike(prob::MuseProblem, x, z, θ) = prob.∇θ_logLike(x, z, θ)
logPriorθ(prob::MuseProblem, θ) = prob.logPriorθ(θ)
logLike_and_∇z_logLike(prob::MuseProblem, x, z, θ) = prob.logLike_and_∇z_logLike(x, z, θ)
sample_x_z(prob::MuseProblem, rng::AbstractRNG, θ) = prob.sample_x_z(rng, θ)