text:
# Collection of functions to compute transit timing variations, using the AH18 integrator.

# wrapper for testing new ics.
@inline function ttv_elements!(el::ElementsIC{T},t0::T,h::T,tmax::T,tt::Array{T,2},count::Array{Int64,1},rstar::T) where T <: Real
    return ttv_elements!(el.H,t0,h,tmax,el.elements,tt,count,0.0,0,0,rstar)
end

"""

Computes Transit Timing Variations (TTVs) as a function of orbital elements, and computes Jacobian of transit times with respect to the initial orbital elements.
"""
function ttv_elements!(H::Union{Int64,Array{Int64,1}},t0::T,h::T,tmax::T,elements::Array{T,2},tt::Array{T,2},count::Array{Int64,1},dtdq0::Array{T,4},rstar::T;fout="",iout=-1,pair=zeros(Bool,H[1],H[1])) where {T <: Real}
    #
    # Input quantities:
    # n     = number of bodies
    # t0    = initial time of integration  [days]
    # h     = time step [days]
    # tmax  = duration of integration [days]
    # elements[i,j] = 2D n x 7 array of the masses & orbital elements of the bodies (currently first body's orbital elements are ignored)
    #            elements are ordered as: mass, period, t0, e*cos(omega), e*sin(omega), inclination, longitude of ascending node (Omega)
    # tt    = array of transit times of size [n x max(ntt)] (currently only compute transits of star, so first row is zero) [days]
    # count = array of the number of transits for each body
    # dtdq0 = derivative of transit times with respect to initial x,v,m [various units: day/length (3), day^2/length (3), day/mass]
    #         4D array  [n x max(ntt) ] x [n x 7] - derivatives of transits of each planet with respect to initial positions/velocities
    #             masses of *all* bodies.  Note: mass derivatives are *after* positions/velocities, even though they are at start
    #             of the elements[i,j] array.
    #
    # Output quantity:
    # dtdelements = 4D array  [n x max(ntt) ] x [n x 7] - derivatives of transits of each planet with respect to initial orbital
    #             elements/masses of *all* bodies.  Note: mass derivatives are *after* elements, even though they are at start
    #             of the elements[i,j] array
    #
    # Example: see test_ttv_elements.jl in test/ directory
    #
    # Define initial mass, position & velocity arrays:
    n = H[1]
    m=zeros(T,n)
    x=zeros(T,NDIM,n)
    v=zeros(T,NDIM,n)
    # Fill the transit-timing & jacobian arrays with zeros:
    fill!(tt,zero(T))
    fill!(dtdq0,zero(T))
    # Create an array for the derivatives with respect to the masses/orbital elements:
    dtdelements = copy(dtdq0)
    # Counter for transits of each planet:
    fill!(count,0)
    for i=1:n
        m[i] = elements[i,1]
    end
    # Initialize the N-body problem using nested hierarchy of Keplerians:
    init = ElementsIC(t0,H,elements)
    x,v,jac_init = init_nbody(init)

    #= Why is this here??
    elements_big=big.(elements); t0big = big(t0); #jac_init_big = zeros(BigFloat,7*n,7*n)
    init_big = ElementsIC(elements_big,H,t0big)
    xbig,vbig,jac_init_big = init_nbody(init_big)
    x = convert(Array{T,2},xbig); v = convert(Array{T,2},vbig); jac_init=convert(Array{T,2},jac_init_big)
    =#

    ttv!(n,t0,h,tmax,m,x,v,tt,count,dtdq0,rstar,pair;fout=fout,iout=iout)
    # Need to apply initial jacobian TBD - convert from
    # derivatives with respect to (x,v,m) to (elements,m):
    ntt_max = size(tt)[2]
    for i=1:n, j=1:count[i]
        if j <= ntt_max
            # Now, multiply by the initial Jacobian to convert time derivatives to orbital elements:
            for k=1:n, l=1:7
                dtdelements[i,j,l,k] = zero(T)
                for p=1:n, q=1:7
                    dtdelements[i,j,l,k] += dtdq0[i,j,q,p]*jac_init[(p-1)*7+q,(k-1)*7+l]
                end
            end
        end
    end
    return dtdelements
end

"""

Computes TTVs as a function of initial x,v,m, and derivatives (dtdq0).
"""
function ttv!(n::Int64,t0::T,h::T,tmax::T,m::Array{T,1},x::Array{T,2},v::Array{T,2},tt::Array{T,2},count::Array{Int64,1},dtdq0::Array{T,4},rstar::T,pair::Array{Bool,2};fout="",iout=-1) where {T <: Real}
    xprior = copy(x)
    vprior = copy(v)
    xtransit = copy(x)
    vtransit = copy(v)
    xerror = zeros(T,size(x)); verror=zeros(T,size(v))
    xerr_trans = zeros(T,size(x)); verr_trans =zeros(T,size(v))
    xerr_prior = zeros(T,size(x)); verr_prior =zeros(T,size(v))
    # Set the time to the initial time:
    t = t0
    # Define error estimate based on Kahan (1965):
    s2 = zero(T)
    # Set step counter to zero:
    istep = 0
    # Jacobian for each step (7- 6 elements+mass, n_planets, 7 - 6 elements+mass, n planets):
    jac_prior = zeros(T,7*n,7*n)
    jac_error_prior = zeros(T,7*n,7*n)
    jac_transit = zeros(T,7*n,7*n)
    jac_trans_err= zeros(T,7*n,7*n)
    # Initialize matrix for derivatives of transit times with respect to the initial x,v,m:
    dtdq = zeros(T,1,7,n)
    # Initialize the Jacobian to the identity matrix:
    #jac_step = eye(T,7*n)
    jac_step = Matrix{T}(I,7*n,7*n)
    # Initialize Jacobian error array:
    jac_error = zeros(T,7*n,7*n)
    if fout != ""
        # Open file for output:
        file_handle =open(fout,"a")
    end
    # Output initial conditions:
    if mod(istep,iout) == 0 && iout > 0
        # Write to file:
        writedlm(file_handle,[convert(Float64,t);convert(Array{Float64,1},reshape(x,3n));convert(Array{Float64,1},reshape(v,3n))]') # Transpose to write each line
    end

    # Save the g function, which computes the relative sky velocity dotted with relative position
    # between the planets and star:
    gsave = zeros(T,n)
    for i=2:n
        # Compute the relative sky velocity dotted with position:
        gsave[i]= g!(i,1,x,v)
    end
    # Loop over time steps:
    dt = zero(T)
    gi = zero(T)
    ntt_max = size(tt)[2]
    param_real = all(isfinite.(x)) && all(isfinite.(v)) && all(isfinite.(m)) && all(isfinite.(jac_step))
    while t < (t0+tmax) && param_real
        # Carry out a ah18 mapping step:
        ah18!(x,v,xerror,verror,h,m,n,jac_step,jac_error,pair)
        param_real = all(isfinite.(x)) && all(isfinite.(v)) && all(isfinite.(m)) && all(isfinite.(jac_step))
        # Check to see if a transit may have occured.  Sky is x-y plane; line of sight is z.
        # Star is body 1; planets are 2-nbody (note that this could be modified to see if
        # any body transits another body):
        for i=2:n
            # Compute the relative sky velocity dotted with position:
            gi = g!(i,1,x,v)
            ri = sqrt(x[1,i]^2+x[2,i]^2+x[3,i]^2)  # orbital distance
            # See if sign of g switches, and if planet is in front of star (by a good amount):
            # (I'm wondering if the direction condition means that z-coordinate is reversed? EA 12/11/2017)
            if gi > 0 && gsave[i] < 0 && x[3,i] > 0.25*ri && ri < rstar
                # A transit has occurred between the time steps - integrate ah18! between timesteps
                count[i] += 1
                if count[i] <= ntt_max
                    dt0 = -gsave[i]*h/(gi-gsave[i])  # Starting estimate
                    xtransit .= xprior; vtransit .= vprior; jac_transit .= jac_prior; jac_trans_err .= jac_error_prior
                    xerr_trans .= xerr_prior; verr_trans .= verr_prior
                    dt = findtransit3!(1,i,n,h,dt0,m,xtransit,vtransit,xerr_trans,verr_trans,jac_transit,jac_trans_err,dtdq,pair) # 20%
                    tt[i,count[i]],stmp = comp_sum(t,s2,dt)
                    # Save for posterity:
                    for k=1:7, p=1:n
                        dtdq0[i,count[i],k,p] = dtdq[1,k,p]
                    end
                end
            end
            gsave[i] = gi
        end
        # Save the current state as prior state:
        xprior .= x
        vprior .= v
        jac_prior .= jac_step
        jac_error_prior .= jac_error
        xerr_prior .= xerror
        verr_prior .= verror
        # Increment time by the time step using compensated summation:
        t,s2 = comp_sum(t,s2,h)
        if mod(istep,iout) == 0 && iout > 0
            # Write to file:
            writedlm(file_handle,[convert(Float64,t);convert(Array{Float64,1},reshape(x,3n));convert(Array{Float64,1},reshape(v,3n));convert(Array{Float64,1},reshape(jac_step,49n^2))]') # Transpose to write each line
        end
        # t += h <- this leads to loss of precision
        # Increment counter by one:
        istep +=1
    end
    if fout != ""
        # Close output file:
        close(file_handle)
    end
    return
end

"""

Finds the transit by taking a partial ah18 step from prior time step, computes timing Jacobian, dtbvdq, with respect to initial cartesian coordinates and masses.
"""
function findtransit3!(i::Int64,j::Int64,n::Int64,h::T,tt::T,m::Array{T,1},x1::Array{T,2},v1::Array{T,2},xerror::Array{T,2},verror::Array{T,2},jac_step::Array{T,2},jac_error::Array{T,2},dtbvdq::Array{T,3},pair::Array{Bool,2}) where {T <: Real}
    # Computes the transit time, approximating the motion as a fraction of a AH17 step forward in time.
    # Also computes the Jacobian of the transit time with respect to the initial parameters, dtbvdq[1-3,7,n].
    # This version is same as findtransit2, but uses the derivative of AH17 step with respect to time
    # rather than the instantaneous acceleration for finding transit time derivative (gdot).
    # Initial guess using linear interpolation:
    dt = one(T)
    iter = 0
    r3 = zero(T)
    gdot = zero(T)
    gsky = zero(T)
    x = copy(x1)
    v = copy(v1)
    xerr_trans = copy(xerror); verr_trans = copy(verror)
    dqdt = zeros(T,7*n)
    TRANSIT_TOL = 10*eps(dt)
    tt1 = tt + 1
    tt2 = tt + 2
    ITMAX = 20
    #while abs(dt) > TRANSIT_TOL && iter < 20
    while true # while true is kind of a no-no...
        tt2 = tt1
        tt1 = tt
        x .= x1; v .= v1; xerr_trans .= xerror; verr_trans .= verror
        # Advance planet state at start of step to estimated transit time:
        ah18!(x,v,xerr_trans,verr_trans,tt,m,n,dqdt,pair)
        # Compute time offset:
        gsky = g!(i,j,x,v)
        #  # Compute derivative of g with respect to time:
        gdot = ((x[1,j]-x[1,i])*(dqdt[(j-1)*7+4]-dqdt[(i-1)*7+4])+(x[2,j]-x[2,i])*(dqdt[(j-1)*7+5]-dqdt[(i-1)*7+5])
                +  (v[1,j]-v[1,i])*(dqdt[(j-1)*7+1]-dqdt[(i-1)*7+1])+(v[2,j]-v[2,i])*(dqdt[(j-1)*7+2]-dqdt[(i-1)*7+2]))
        # Refine estimate of transit time with Newton's method:
        dt = -gsky/gdot
        # Add refinement to estimated time:
        tt += dt
        iter +=1
        # Break out if we have reached maximum iterations, or if
        # current transit time estimate equals one of the prior two steps:
        if (iter >= ITMAX) || (tt == tt1) || (tt == tt2)
            break
        end
    end
    if iter >= 20
        #  println("Exceeded iterations: planet ",j," iter ",iter," dt ",dt," gsky ",gsky," gdot ",gdot)
    end
    # Compute time derivatives:
    x .= x1; v .= v1; xerr_trans .= xerror; verr_trans .= verror
    # Compute dgdt with the updated time step.
    ah18!(x,v,xerr_trans,verr_trans,tt,m,n,jac_step,jac_error,pair)
    # Need to reset to compute dqdt:
    x .= x1; v .= v1; xerr_trans .= xerror; verr_trans .= verror
    ah18!(x,v,xerr_trans,verr_trans,tt,m,n,dqdt,pair)
    # Compute time offset:
    gsky = g!(i,j,x,v)
    # Compute derivative of g with respect to time:
    gdot  = ((x[1,j]-x[1,i])*(dqdt[(j-1)*7+4]-dqdt[(i-1)*7+4])+(x[2,j]-x[2,i])*(dqdt[(j-1)*7+5]-dqdt[(i-1)*7+5])
             +  (v[1,j]-v[1,i])*(dqdt[(j-1)*7+1]-dqdt[(i-1)*7+1])+(v[2,j]-v[2,i])*(dqdt[(j-1)*7+2]-dqdt[(i-1)*7+2]))
    # Set dtbvdq to zero:
    fill!(dtbvdq,zero(T))
    indj = (j-1)*7+1
    indi = (i-1)*7+1
    @inbounds for p=1:n
        indp = (p-1)*7
        @inbounds for k=1:7
            # Compute derivatives:
            dtbvdq[1,k,p] = -((jac_step[indj  ,indp+k]-jac_step[indi  ,indp+k])*(v[1,j]-v[1,i])+(jac_step[indj+1,indp+k]-jac_step[indi+1,indp+k])*(v[2,j]-v[2,i])+
                              (jac_step[indj+3,indp+k]-jac_step[indi+3,indp+k])*(x[1,j]-x[1,i])+(jac_step[indj+4,indp+k]-jac_step[indi+4,indp+k])*(x[2,j]-x[2,i]))/gdot
        end
    end
    # Note: this is the time elapsed *after* the beginning of the timestep:
    ntbv = size(dtbvdq)[1]
    if ntbv == 3
        # Compute the impact parameter and sky velocity:
        vsky = sqrt((v[1,j]-v[1,i])^2 + (v[2,j]-v[2,i])^2)
        bsky2 = (x[1,j]-x[1,i])^2 + (x[2,j]-x[2,i])^2
        # partial derivative v_{sky} with respect to time:
        dvdt = ((v[1,j]-v[1,i])*(dqdt[(j-1)*7+4]-dqdt[(i-1)*7+4])+(v[2,j]-v[2,i])*(dqdt[(j-1)*7+5]-dqdt[(i-1)*7+5]))/vsky
        # (note that \partial b/\partial t = 0 at mid-transit since g_{sky} = 0 mid-transit).
        @inbounds for p=1:n
            indp = (p-1)*7
            @inbounds for k=1:7
                # Compute derivatives:
                #v_{sky} = sqrt((v[1,j]-v[1,i])^2+(v[2,j]-v[2,i])^2)
                dtbvdq[2,k,p] = ((jac_step[indj+3,indp+k]-jac_step[indi+3,indp+k])*(v[1,j]-v[1,i])+(jac_step[indj+4,indp+k]-jac_step[indi+4,indp+k])*(v[2,j]-v[2,i]))/vsky + dvdt*dtbvdq[1,k,p]
                #b_{sky}^2 = (x[1,j]-x[1,i])^2+(x[2,j]-x[2,i])^2
                dtbvdq[3,k,p] = 2*((jac_step[indj  ,indp+k]-jac_step[indi  ,indp+k])*(x[1,j]-x[1,i])+(jac_step[indj+1,indp+k]-jac_step[indi+1,indp+k])*(x[2,j]-x[2,i]))
            end
        end
        # return the transit time, impact parameter, and sky velocity:
        return tt::T,vsky::T,bsky2::T
    else
        return tt::T
    end
end

"""

Finds the transit by taking a partial ah18 step from prior times step, computes timing Jacobian, dtdq, wrt initial cartesian coordinates, masses:
"""
function findtransit3!(i::Int64,j::Int64,n::Int64,h::T,tt::T,m::Array{T,1},x1::Array{T,2},v1::Array{T,2},xerror::Array{T,2},verror::Array{T,2},pair::Array{Bool,2};calcbvsky=false) where {T <: Real}
    # Computes the transit time, approximating the motion as a fraction of a DH17 step forward in time.
    # Also computes the Jacobian of the transit time with respect to the initial parameters, dtdq[7,n].
    # This version is same as findtransit2, but uses the derivative of dh17 step with respect to time
    # rather than the instantaneous acceleration for finding transit time derivative (gdot).
    # Initial guess using linear interpolation:
    dt = one(T)
    iter = 0
    r3 = zero(T)
    gdot = zero(T)
    gsky = gdot
    x = copy(x1); v = copy(v1); xerr_trans = copy(xerror); verr_trans = copy(verror)
    dqdt = zeros(T,7*n)
    #TRANSIT_TOL = 10*sqrt(eps(dt)
    TRANSIT_TOL = 10*eps(dt)
    ITMAX = 20
    tt1 = tt + 1
    tt2 = tt + 2
    #while abs(dt) > TRANSIT_TOL && iter < 20
    while true
        tt2 = tt1
        tt1 = tt
        x .= x1; v .= v1; xerr_trans .= xerror; verr_trans .= verror
        # Advance planet state at start of step to estimated transit time:
        #  dh17!(x,v,tt,m,n,pair)
        #dh17!(x,v,xerror,verror,tt,m,n,dqdt,pair)
        ah18!(x,v,xerr_trans,verr_trans,tt,m,n,dqdt,pair)
        # Compute time offset:
        gsky = g!(i,j,x,v)
        #  # Compute derivative of g with respect to time:
        gdot = ((x[1,j]-x[1,i])*(dqdt[(j-1)*7+4]-dqdt[(i-1)*7+4])+(x[2,j]-x[2,i])*(dqdt[(j-1)*7+5]-dqdt[(i-1)*7+5])
                +  (v[1,j]-v[1,i])*(dqdt[(j-1)*7+1]-dqdt[(i-1)*7+1])+(v[2,j]-v[2,i])*(dqdt[(j-1)*7+2]-dqdt[(i-1)*7+2]))
        # Refine estimate of transit time with Newton's method:
        dt = -gsky/gdot
        # Add refinement to estimated time:
        tt += dt
        iter +=1
        # Break out if we have reached maximum iterations, or if
        # current transit time estimate equals one of the prior two steps:
        if (iter >= ITMAX) || (tt == tt1) || (tt == tt2)
            break
        end
    end
    if iter >= 20
        #  println("Exceeded iterations: planet ",j," iter ",iter," dt ",dt," gsky ",gsky," gdot ",gdot)
    end
    # Note: this is the time elapsed *after* the beginning of the timestep:
    if calcbvsky
        # Compute the sky velocity and impact parameter:
        vsky = sqrt((v[1,j]-v[1,i])^2 + (v[2,j]-v[2,i])^2)
        bsky2 = (x[1,j]-x[1,i])^2 + (x[2,j]-x[2,i])^2
        # return the transit time, impact parameter, and sky velocity:
        return tt::T,vsky::T,bsky2::T
    else
        return tt::T
    end
end

"""Used in computing transit time inside `findtransit3`."""
function g!(i::Int64,j::Int64,x::Array{T,2},v::Array{T,2}) where {T <: Real}
    # See equation 8-10 Fabrycky (2008) in Seager Exoplanets book
    g = (x[1,j]-x[1,i])*(v[1,j]-v[1,i])+(x[2,j]-x[2,i])*(v[2,j]-v[2,i])
    return g
end

"""

Computes TTVs as a function of orbital elements, allowing for a single log perturbation of dlnq for body jq and element iq. (NO GRADIENT)
"""
function ttv_elements!(H::Union{Int64,Array{Int64,1}},t0::T,h::T,tmax::T,elements::Array{T,2},tt::Array{T,2},count::Array{Int64,1},dlnq::T,iq::Int64,jq::Int64,rstar::T;fout="",iout=-1,pair = zeros(Bool,H[1],H[1])) where {T <: Real}
    #
    # Input quantities:
    # n     = number of bodies
    # t0    = initial time of integration  [days]
    # h     = time step [days]
    # tmax  = duration of integration [days]
    # elements[i,j] = 2D n x 7 array of the masses & orbital elements of the bodies (currently first body's orbital elements are ignored)
    #            elements are ordered as: mass, period, t0, e*cos(omega), e*sin(omega), inclination, longitude of ascending node (Omega)
    # tt    = pre-allocated array to hold transit times of size [n x max(ntt)] (currently only compute transits of star, so first row is zero) [days]
    #         upon output, set to transit times of planets.
    # count = pre-allocated array of the number of transits for each body upon output
    #
    # dlnq  = fractional variation in initial parameter jq of body iq for finite-difference calculation of
    #         derivatives [this is only needed for testing derivative code, below].
    #
    # Example: see test_ttv_elements.jl in test/ directory
    #
    #fcons = open("fcons.txt","w");
    # Set up mass, position & velocity arrays.  NDIM =3
    n = H[1]
    m=zeros(T,n)
    x=zeros(T,NDIM,n)
    v=zeros(T,NDIM,n)
    # Fill the transit-timing array with zeros:
    fill!(tt,0.0)
    # Counter for transits of each planet:
    fill!(count,0)
    # Insert masses from the elements array:
    for i=1:n
        m[i] = elements[i,1]
    end
    # Allow for perturbations to initial conditions: jq labels body; iq labels phase-space element (or mass)
    # iq labels phase-space element (1-3: x; 4-6: v; 7: m)
    dq = zero(T)
    if iq == 7 && dlnq != 0.0
        dq = m[jq]*dlnq
        m[jq] += dq
    end
    # Initialize the N-body problem using nested hierarchy of Keplerians:
    init = ElementsIC(t0,H,elements)
    x,v,_ = init_nbody(init)
    #elements_big=big.(elements); t0big = big(t0)
    #init_big = ElementsIC(elements_big,H,t0big)
    #xbig,vbig,_ = init_nbody(init_big)
    #if T != BigFloat
    #  println("Difference in x,v init: ",x-convert(Array{T,2},xbig)," ",v-convert(Array{T,2},vbig)," (dlnq version)")
    #end
    #x = convert(Array{T,2},xbig); v = convert(Array{T,2},vbig)
    # Perturb the initial condition by an amount dlnq (if it is non-zero):
    if dlnq != 0.0 && iq > 0 && iq < 7
        if iq < 4
            if x[iq,jq] != 0
                dq = x[iq,jq]*dlnq
            else
                dq = dlnq
            end
            x[iq,jq] += dq
        else
            # Same for v
            if v[iq-3,jq] != 0
                dq = v[iq-3,jq]*dlnq
            else
                dq = dlnq
            end
            v[iq-3,jq] += dq
        end
    end
    ttv!(n,t0,h,tmax,m,x,v,tt,count,fout,iout,rstar,pair)
    return dq
end

"""

Computes TTVs as a function of initial x,v,m. (NO GRADIENT).
"""
function ttv!(n::Int64,t0::T,h::T,tmax::T,m::Array{T,1},x::Array{T,2},v::Array{T,2},tt::Array{T,2},count::Array{Int64,1},fout::String,iout::Int64,rstar::T,pair::Array{Bool,2}) where {T <: Real}
    # Make some copies to allocate space for saving prior step and computing coordinates at the times of transit.
    xprior = copy(x)
    vprior = copy(v)
    xtransit = copy(x)
    vtransit = copy(v)
    #xerror = zeros(x); verror = zeros(v)
    xerror = copy(x); verror = copy(v)
    fill!(xerror,zero(T)); fill!(verror,zero(T))
    #xerr_prior = zeros(x); verr_prior = zeros(v)
    xerr_prior = copy(xerror); verr_prior = copy(verror)
    # Set the time to the initial time:
    t = t0
    # Define error estimate based on Kahan (1965):
    s2 = zero(T)
    # Set step counter to zero:
    istep = 0
    # Jacobian for each step (7 elements+mass, n_planets, 7 elements+mass, n planets):
    # Save the g function, which computes the relative sky velocity dotted with relative position
    # between the planets and star:
    gsave = zeros(T,n)
    gi  = 0.0
    dt::T = 0.0
    # Loop over time steps:
    ntt_max = size(tt)[2]
    param_real = all(isfinite.(x)) && all(isfinite.(v)) && all(isfinite.(m))
    if fout != ""
        # Open file for output:
        file_handle =open(fout,"a")
    end
    # Output initial conditions:
    if mod(istep,iout) == 0 && iout > 0
        # Write to file:
        writedlm(file_handle,[convert(Float64,t);convert(Array{Float64,1},reshape(x,3n));convert(Array{Float64,1},reshape(v,3n))]') # Transpose to write each line
    end
    while t < t0+tmax && param_real
        # Carry out a phi^2 mapping step:
        #  phi2!(x,v,h,m,n)
        ah18!(x,v,xerror,verror,h,m,n,pair)
        #  xbig = big.(x); vbig = big.(v); hbig = big(h); mbig = big.(m)
        #dh17!(x,v,xerror,verror,h,m,n,pair)
        #  dh17!(xbig,vbig,hbig,mbig,n,pair)
        #  x .= convert(Array{Float64,2},xbig); v .= convert(Array{Float64,2},vbig)
        param_real = all(isfinite.(x)) && all(isfinite.(v)) && all(isfinite.(m))
        # Check to see if a transit may have occured.  Sky is x-y plane; line of sight is z.
        # Star is body 1; planets are 2-nbody:
        for i=2:n
            # Compute the relative sky velocity dotted with position:
            gi = g!(i,1,x,v)
            ri = sqrt(x[1,i]^2+x[2,i]^2+x[3,i]^2)
            # See if sign switches, and if planet is in front of star (by a good amount):
            if gi > 0 && gsave[i] < 0 && x[3,i] > 0.25*ri && ri < rstar
                # A transit has occurred between the time steps.
                # Approximate the planet-star motion as a Keplerian, weighting over timestep:
                count[i] += 1
                #      tt[i,count[i]]=t+findtransit!(i,h,gi,gsave[i],m,xprior,vprior,x,v,pair)
                if count[i] <= ntt_max
                    dt0 = -gsave[i]*h/(gi-gsave[i])
                    xtransit .= xprior
                    vtransit .= vprior
                    #        dt = findtransit2!(1,i,n,h,dt0,m,xtransit,vtransit,pair)
                    #hbig = big(h); dt0big=big(dt0); mbig=big.(m); xtbig = big.(xtransit); vtbig = big.(vtransit)
                    #dtbig = findtransit2!(1,i,n,hbig,dt0big,mbig,xtbig,vtbig,pair)
                    #dt = convert(Float64,dtbig)
                    dt = findtransit3!(1,i,n,h,dt0,m,xtransit,vtransit,xerr_prior,verr_prior,pair)
                    #tt[i,count[i]]=t+dt
                    tt[i,count[i]],stmp = comp_sum(t,s2,dt)
                end
                #      tt[i,count[i]]=t+findtransit2!(1,i,n,h,gi,gsave[i],m,xprior,vprior,pair)
            end
            gsave[i] = gi
        end
        # Save the current state as prior state:
        xprior .= x
        vprior .= v
        xerr_prior .= xerror
        verr_prior .= verror
        # Increment time by the time step using compensated summation:
        #s2 += h; tmp = t + s2; s2 = (t - tmp) + s2
        #t = tmp
        t,s2 = comp_sum(t,s2,h)
        if mod(istep,iout) == 0 && iout > 0
            # Write to file:
            writedlm(file_handle,[convert(Float64,t);convert(Array{Float64,1},reshape(x,3n));convert(Array{Float64,1},reshape(v,3n))]') # Transpose to write each line
        end
        # t += h  <- this leads to loss of precision
        # Increment counter by one:
        istep +=1
        #  println("xerror: ",xerror)
        #  println("verror: ",verror)
    end
    #println("xerror: ",xerror)
    #println("verror: ",verror)
    if fout != ""
        # Close output file:
        close(file_handle)
    end
    return
end


text:
@testset "Cipher" begin
  @testset "Cipher: atbash" begin
        input = "test"
        @test encode(input) == "gvhg"
        input = "abcdefghijklmnopqrstuvwxyz"
        @test encode(input) == "zyxwvutsrqponmlkjihgfedcba"
        input = "hello"
        @test encode(input) == "svool"
    end
end

ext: jl

max_stars_repo_name: standardgalactic/Julia-algorithms


text:
#=
Tokenize openfoam polymesh dictionary
    - Currently only implemented for ascii files

TODO: Implement function converting binary and gz files to Vector{Token}
=#

const TokenTypeID_Dict = Dict(
    :eof => 0,
    :emptyspace => 1,
    :newline => 2,
    :string => 3,
    :word => 4,
    :int_number => 5,
    :float_number => 6,
    :comment => 7,
    :bar => Int64('|'),
    :bracket_close => Int64(')'),
    :bracket_open => Int64('('),
    :colon => Int64(':'),
    :comma => Int64(','),
    :curly_brace_close => Int64('}'),
    :curly_brace_open => Int64('{'),
    :dot => Int64('.'),
    :doublequote => Int64('"'),
    :equal => Int64('='),
    :greater => Int64('>'),
    :minus => Int64('-'),
    :plus => Int64('+'),
    :semicolon => Int64(';'),
    :slash_bw => Int64('\\'),
    :slash_fw => Int64('/'),
    :smaller => Int64('<'),
    :star => Int64('*'),
   )

const TokenSymbol_Dict = Dict(value => key for (key, value) in TokenTypeID_Dict)

struct Token
    type_id::Int64 # see TokenTypeID_Dict
    string::String
end

mutable struct Lexer
    buffer::String
    pos::Int64
    Lexer(buffer::String) = new(buffer, 1)
end

#= ---------------------- =#
#= Lexer iterator methods =#
#= ---------------------- =#

"""
    start!(lex::Lexer)

reset postion of lexer to 1
"""
function start!(lex::Lexer)
    lex.pos = 1
end

done(lex::Lexer) = length(lex.buffer) < lex.pos
iteratorsize(lex::Lexer) =  Base.SizeUnknown()

"""
    next_ascii!(lex::Lexer)

divide lexer into tokens -> return next token and set lex.pos after token
"""
function next_ascii!(lex::Lexer)

    if done(lex)
        return (Token(TokenTypeID_Dict[:eof], ""), lex.pos)
    end

    ch = lex.buffer[lex.pos]

    # while empty space skip forward to next possible token
    while ch in " \t\n"
        lex.pos += 1
        if done(lex)
            return Token(TokenTypeID_Dict[:eof], "")
        end
        ch = lex.buffer[lex.pos]
    end

    if in(ch, "/")
        # Looking for comments
        if lex.buffer[lex.pos+1] == '/' # line comment
            spos = lex.pos + 2
            npos = findnext("//", lex.buffer,spos)
            lex.pos =  last(npos) + 1
            return Token(TokenTypeID_Dict[:comment],
                                lex.buffer[spos:(first(npos)-1)])

        elseif lex.buffer[lex.pos+1] == '*' # multi line comment
            spos = lex.pos + 2
            npos = findnext("*/", lex.buffer, spos)
            lex.pos =  last(npos) + 1
            return Token(TokenTypeID_Dict[:comment],
                                lex.buffer[spos:(first(npos)-1)])

        end

        lex.pos += 1
        return Token(Int(ch), string(ch))
    end

    if ch == '"'
        # looking for quoted strings
        spos = lex.pos + 1
        npos = findnext('"', lex.buffer, spos) # TODO: Does not handle escaped quotes
        lex.pos =  npos + 1
        return Token(TokenTypeID_Dict[:string], lex.buffer[spos : (first(npos)-1)])
    elseif  in(ch, "+-/*{}();\\:=,|.<>")
          # Looking for general syntax tokens
          lex.pos += 1
          return Token(Int(ch), string(ch))

    elseif isdigit(ch) || ch == '-'
        # numbers musst start with digit or - no support for .1 etc ...
        #looking for numbers
        # range = findnext(r"\d+", lex.buffer, lex.pos)
        # lex.pos =  last(range) + 1
        # return Token(TokenTypeID_Dict[:number], lex.buffer[range])
        return tokenize_ascii_number!(lex)
    elseif isletter(ch)
        # looking for word
        range = findnext(r"\w+", lex.buffer, lex.pos)
        lex.pos =  last(range) + 1
        return Token(TokenTypeID_Dict[:word], lex.buffer[range])

    else
        error("Unknown character '$ch' at position: $(lex.pos)")
    end
end

"""
    tokenize_ascii_number!(lex::Lexer)

find strings that represent numbers and return token, also set lex.pos after
number string
"""
function tokenize_ascii_number!(lex::Lexer)
    regexstr = r"[+\-]?(?:0|[1-9]\d*)(?:\.\d*)?(?:[eE][+\-]?\d+)?"

    range = findnext(regexstr, lex.buffer, lex.pos)

    if (!isempty(range) && first(range) == lex.pos)

        lex.pos = last(range) + 1

        if (contains(lex.buffer[range], '.') ||
            contains(lex.buffer[range], 'E') ||
            contains(lex.buffer[range], 'e') )

            return Token(TokenTypeID_Dict[:float_number], lex.buffer[range])
        else
            return Token(TokenTypeID_Dict[:int_number], lex.buffer[range])
        end

    else
        lex.pos += 1
        return Token(Int(lex.buffer[lex.pos]), lex.buffer[range])
    end
end

#= -------------------------- =#
#= End lexer iterator methods =#
#= -------------------------- =#

"""
    tokenize_ascii(lex::Lexer)

tokenize of lexer (buffer)
"""
function tokenize_ascii(lex::Lexer)
    tokens = Vector{Token}(undef,0)
    start!(lex) # just to be shure lex.pos = 1

    while true
        token = next_ascii!(lex)

        push!(tokens, token)

        if token.type_id == TokenTypeID_Dict[:eof]
            break
        end
    end

    return tokens
end

"""
    tokenize_ofmeshfile(filename::String)

tokenize openfoam (ascii) dictionary (filename) from lex buffer
"""
function tokenize_ascii_ofmeshfile(filename::String)
    file = open(filename,"r")
    lexer = Lexer(read(file, String))
    close(file)

    tokens = tokenize_ascii(lexer)
    remove_comments!(tokens)
end

"""
    remove_comments!(tokens::Vector{Token})

remove comments from tokens
"""
function remove_comments!(tokens::Vector{Token})
    filter!(t->t.type_id != TokenTypeID_Dict[:comment], tokens)
end

"""
    printtokens(tokens::Vector{Token})

simple printing of tokens
"""
function printtokens(tokens::Vector{Token})
    for (i,t) in enumerate(tokens)
        println(i," ",TokenSymbol_Dict[t.type_id], " \n\t\t: ", t.string, " \n")
    end
end

max_stars_repo_name: c-schubert/linearMeshIO

max_issues_repo_name: c-schubert/linearMeshIO

max_forks_repo_name: c-schubert/linearMeshIO


text:
import Term: Panel, TextBox, cleantext

@testset "\e[34mPanel - panel creation" begin
    p1 = Panel("this panel has fixed width"; width = 44, justify = :right)
    @test string(p1) ==
        "╭──────────────────────────────────────────╮\n│               this panel has fixed width │\n╰──────────────────────────────────────────╯"

    p2 = Panel(
        "this one too, but the text is at the center!"; width = 66, justify = :center
    )
    @test string(p2) ==
        "╭────────────────────────────────────────────────────────────────╮\n│          this one too, but the text is at the center!          │\n╰────────────────────────────────────────────────────────────────╯"

    p3 = Panel("this one fits its content"; width = :fit)
    @test string(p3) ==
        "╭───────────────────────────╮\n│ this one fits its content │\n╰───────────────────────────╯"

    p4 = Panel(
        "[red]This is the panel's first line.[/red]",
        "[bold green]and this is another, panel just stacks all inputs into one piece of content[/bold green]";
        width = :fit,
    )
    @test string(p4) ==
        "╭─────────────────────────────────────────────────────────────────────────────╮\n│ \e[31mThis is the panel's first line.\e[39m                                             │\n│ \e[1m\e[32mand this is another, panel just stacks all inputs into one piece of content\e[22m\e[39m │\n╰─────────────────────────────────────────────────────────────────────────────╯"

    p5 = Panel(
        "content "^10;
        subtitle = "another panel",
        subtitle_style = "dim underline",
        subtitle_justify = :right,
        width = 44,
    )
    @test string(p5) ==
        "╭──────────────────────────────────────────╮\n│ content content content content content  │\n│ content content content content content  │\n╰─────────────────────────\e[0m \e[2m\e[4manother panel\e[22m\e[24m ──╯"

    p6 = Panel("content "^10; box = :DOUBLE, style = "blue", width = 44)
    @test string(p6) ==
        "\e[34m╔══════════════════════════════════════════╗\e[39m\n\e[34m║\e[39m content content content content content  \e[34m║\e[39m\n\e[34m║\e[39m content content content content content  \e[34m║\e[39m\n\e[34m╚══════════════════════════════════════════╝\e[39m"

    # test with wider chars
    p7 = Panel("나랏말싸미 듕귁에 달아")
    p8 = Panel("こんにちは(わ)")

    @test p7.measure.w == 26
    @test p8.measure.w == 18
    @test p8.measure.h == 3

    @test string(p7) ==
        "╭────────────────────────╮\n│ 나랏말싸미 듕귁에 달아 │\n╰────────────────────────╯"
    @test string(p8) == "╭────────────────╮\n│ こんにちは(わ) │\n╰────────────────╯"

    # create empty panels with varying size
    p9 = Panel(; height = 20, width = 4)
    p10 = Panel(""; height = 5)
    p11 = Panel(; width = 12, height = 6)

    # check all  lines in each panel have exactly the same size
    for panel in (p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11)
        _p = string(panel)
        widths = textwidth.(cleantext.(split(_p, "\n")))
        @test length(unique(widths)) == 1
    end
end

@testset "\e[34mPanel - panel layout  " begin
    pright = Panel(
        "content [red]with style[/red] "^26;
        title = "My Panel",
        title_style = "bold red",
        width = 44,
    )
    pleft = Panel("content "^30; box = :DOUBLE, style = "blue", width = 66)
    pp = Panel(
        pleft / pright;
        style = "green dim",
        title_style = "green",
        title = "vertically stacked!",
        title_justify = :center,
        subtitle = "styled by Term.jl",
        subtitle_justify = :right,
        subtitle_style = "dim",
        justify = :center,
    )

    @test pright.measure.w == 44
    @test pright.measure.h == 15

    @test pleft.measure.w == 66
    @test pleft.measure.h == 6

    @test pp.measure.w == 70
    @test pp.measure.h == 23

    pv = pright / pleft
    @test pv.measure.w == 66
    @test pv.measure.h == 21

    ph = pright * pleft
    @test ph.measure.w == 110
    @test ph.measure.h == 15
end

@testset "\e[34mPanel - tbox creation " begin
    my_long_text = "Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat. Duis aute irure dolor in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur. Excepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mollit anim id est laborum."
    tb = TextBox(my_long_text; width = 44)

    @test tb.measure.w == 44
    @test tb.measure.h == 14
    @test string(tb) ==
        "\e[8m╭──────────────────────────────────────────╮\e[28m\n\e[8m│\e[28m Lorem ipsum dolor sit amet, consectetur  \e[8m│\e[28m\n\e[8m│\e[28m adipiscing elit, sed do eiusmod tempor i \e[8m│\e[28m\n\e[8m│\e[28m ncididunt ut labore et dolore magna aliq \e[8m│\e[28m\n\e[8m│\e[28m ua. Ut enim ad minim veniam, quis nostru \e[8m│\e[28m\n\e[8m│\e[28m d exercitation ullamco laboris nisi ut a \e[8m│\e[28m\n\e[8m│\e[28m liquip ex ea commodo consequat. Duis aut \e[8m│\e[28m\n\e[8m│\e[28m e irure dolor in reprehenderit in volupt \e[8m│\e[28m\n\e[8m│\e[28m ate velit esse cillum dolore eu fugiat n \e[8m│\e[28m\n\e[8m│\e[28m ulla pariatur. Excepteur sint occaecat c \e[8m│\e[28m\n\e[8m│\e[28m upidatat non proident, sunt in culpa qui \e[8m│\e[28m\n\e[8m│\e[28m officia deserunt mollit anim id est lab  \e[8m│\e[28m\n\e[8m│\e[28m orum.                                    \e[8m│\e[28m\n\e[8m╰──────────────────────────────────────────╯\e[28m"

    tb2 = TextBox(
        my_long_text;
        title = "This is my long text",
        title_style = "bold red",
        title_justify = :center,
        subtitle = "styled by Term.jl",
        subtitle_justify = :right,
        subtitle_style = "dim",
        width = 88,
    )
    @test string(tb2) ==
        "\e[8m╭────────────────────────────────\e[28m\e[0m \e[1m\e[31mThis is my long text\e[22m\e[39m \e[8m────────────────────────────────╮\e[28m\n\e[8m│\e[28m Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incid \e[8m│\e[28m\n\e[8m│\e[28m idunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud exerci \e[8m│\e[28m\n\e[8m│\e[28m tation ullamco laboris nisi ut aliquip ex ea commodo consequat. Duis aute irure dolo \e[8m│\e[28m\n\e[8m│\e[28m r in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur. E \e[8m│\e[28m\n\e[8m│\e[28m xcepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mo \e[8m│\e[28m\n\e[8m│\e[28m llit anim id est laborum.                                                            \e[8m│\e[28m\n\e[8m╰─────────────────────────────────────────────────────────────────\e[28m\e[8m\e[0m \e[2mstyled by Term.jl\e[22m \e[8m──╯\e[28m\e[28m"
end

@testset "\e[34mPannel - tbox markup  " begin
    another_long_one =
        "This is a [red bold]very[/red bold] piece of [green italic]content[/green italic]. But TextBox can handle [underline]anything[/underline]!! "^10
    tb = TextBox(another_long_one; width = 44)

    @test tb.measure.w == 44
    @test tb.measure.h == 19

    @test string(tb) ==
        "\e[8m╭──────────────────────────────────────────╮\e[28m\n\e[8m│\e[28m This is a \e[1m\e[31mvery\e[22m\e[39m piece of \e[3m\e[32mcontent\e[23m\e[39m. But Tex \e[8m│\e[28m\n\e[8m│\e[28m tBox can handle \e[4manything\e[24m!! This is a \e[1m\e[31mver\e[22m\e[39m \e[8m│\e[28m\n\e[8m│\e[28m \e[1m\e[31my\e[22m\e[39m piece of \e[3m\e[32mcontent\e[23m\e[39m. But TextBox can hand \e[8m│\e[28m\n\e[8m│\e[28m le \e[4manything\e[24m!! This is a \e[1m\e[31mvery\e[22m\e[39m piece of \e[3m\e[32mco\e[23m\e[39m \e[8m│\e[28m\n\e[8m│\e[28m \e[3m\e[32mntent\e[23m\e[39m. But TextBox can handle \e[4manything\e[24m!! \e[8m│\e[28m\n\e[8m│\e[28m This is a \e[1m\e[31mvery\e[22m\e[39m piece of \e[3m\e[32mcontent\e[23m\e[39m. But Te  \e[8m│\e[28m\n\e[8m│\e[28m xtBox can handle \e[4manything\e[24m!! This is a \e[1m\e[31mve\e[22m\e[39m \e[8m│\e[28m\n\e[8m│\e[28m \e[1m\e[31mry\e[22m\e[39m piece of \e[3m\e[32mcontent\e[23m\e[39m. But TextBox can han \e[8m│\e[28m\n\e[8m│\e[28m dle \e[4manything\e[24m!! This is a \e[1m\e[31mvery\e[22m\e[39m piece of \e[3m\e[32mc\e[23m\e[39m \e[8m│\e[28m\n\e[8m│\e[28m \e[3m\e[32montent\e[23m\e[39m. But TextBox can handle \e[4manything\e[24m! \e[8m│\e[28m\n\e[8m│\e[28m ! This is a \e[1m\e[31mvery\e[22m\e[39m piece of \e[3m\e[32mcontent\e[23m\e[39m. But T \e[8m│\e[28m\n\e[8m│\e[28m extBox can handle \e[4manything\e[24m!! This is a \e[1m\e[31mv\e[22m\e[39m \e[8m│\e[28m\n\e[8m│\e[28m \e[1m\e[31mery\e[22m\e[39m piece of \e[3m\e[32mcontent\e[23m\e[39m. But TextBox can ha \e[8m│\e[28m\n\e[8m│\e[28m ndle \e[4manything\e[24m!! This is a \e[1m\e[31mvery\e[22m\e[39m piece of  \e[8m│\e[28m\n\e[8m│\e[28m \e[3m\e[32mcontent\e[23m\e[39m. But TextBox can handle \e[4manything\e[24m \e[8m│\e[28m\n\e[8m│\e[28m \e[4m\e[24m!! This is a \e[1m\e[31mvery\e[22m\e[39m piece of \e[3m\e[32mcontent\e[23m\e[39m. But  \e[8m│\e[28m\n\e[8m│\e[28m TextBox can handle \e[4manything\e[24m!!            \e[8m│\e[28m\n\e[8m╰──────────────────────────────────────────╯\e[28m"
end


text:
using Random, Statistics, Plots, Printf
Random.seed!(430)
N = 10000;

# generate 10000 realizations of a Gaussian random variable Xn with mean 1 and variance 4
X_n = 2 * randn(N) .+ 1;
@printf("%d Gaussian random variable Xn with mean: %.5f and variance: %.5f\n", N, mean(X_n), var(X_n))

# Calculate S_n
S_n = 1/length(X_n[1:N]) * sum(X_n[1:N])
@printf("Obtained Sn: %.5f from Xn samples\n", S_n)

# Define the running average
function running_avg(n)
    values = []
    for i in 1:n
        S_n = 1/length(X_n[1:i]) * sum(X_n[1:i])
        push!(values, S_n)
    end
    return values
end

# Samples for each n values from 1 to 10,000
X_nSamples = running_avg(N);

# Plot |Sn − 1| as a function of n
S_nMinus1 = abs.(X_nSamples .- 1);
display(plot(1:N, S_nMinus1, labels="|Sn − 1|", xlabel="N", ylabel="|Sn − 1|"))

# Define the time correlation
function time_correlation(X_n, n)
    Kn_list = []
    for n in 1:N-1
        Kn=0
        for k in 1:n
            Kn += X_n[k]*X_n[k+1]
        end
        val = 1/n * Kn;
        push!(Kn_list, val)
    end
    return Kn_list
end

# Plot Kn
T = time_correlation(X_n, N);
display(plot(1:length(T), T, labels="Kn", xlabel="N", ylabel="Kn"))


text:
module MarkdownFormatTests

using Test
using Random

using Documenter

# Documenter package docs
@info("Building Documenter's docs with Markdown.")
const Documenter_root = normpath(joinpath(@__DIR__, "..", "..", "docs"))
build_dir_relpath = relpath(joinpath(@__DIR__, "builds/markdown"), Documenter_root)
doc = makedocs(
    format = :markdown,
    debug   = true,
    root    = Documenter_root,
    modules = Documenter,
    build   = build_dir_relpath,
)

@testset "Markdown" begin
    @test isa(doc, Documenter.Documents.Document)

    let build_dir  = joinpath(Documenter_root, build_dir_relpath),
        source_dir = joinpath(Documenter_root, "src")
        @test isdir(build_dir)
        @test isdir(joinpath(build_dir, "assets"))
        @test isdir(joinpath(build_dir, "lib"))
        @test isdir(joinpath(build_dir, "man"))

        @test isfile(joinpath(build_dir, "index.md"))
        @test isfile(joinpath(build_dir, "assets", "mathjaxhelper.js"))
        @test isfile(joinpath(build_dir, "assets", "Documenter.css"))
    end

    @test doc.user.root   == Documenter_root
    @test doc.user.source == "src"
    @test doc.user.build  == build_dir_relpath
    @test doc.user.clean  == true
end

end


text:
# Xpress model

#################################################
#
#  model type & constructors
#
#################################################

"""
    Model

Type to hold an Xpress model
"""
mutable struct Model
    ptr_model::Ptr{Nothing}
    callback::Array{Any}
    finalize_env::Bool

    time::Float64

    function Model(p::Ptr{Nothing}; finalize_env::Bool=true)
        model = new(p, Any[], finalize_env, 0.0)
        if VERSION >= v"0.7-"
            if finalize_env
                finalizer(model) do m
                    free_model(m)
                end
            end
        else
            if finalize_env
                finalizer(model, m -> free_model(m))
            end
        end
        model
    end
end

"""
    Model(; finalize_env::Bool=false)

Xpress model constructor (autimatically sets OUTPUTLOG = 1)
"""
function Model(; finalize_env::Bool=true)

    a = Array{Ptr{Nothing}}(undef, 1)
    ret = @xprs_ccall(createprob, Cint, ( Ptr{Ptr{Nothing}},), a )
    if ret != 0
        error("It was not possible to create a model, try running Env() and then create the model again.")
    end

    m = Model(a[1]; finalize_env = finalize_env)

    # turn off default printing on unix
    setparam!(m, XPRS_OUTPUTLOG, 0)
    setparam!(m, XPRS_CALLBACKFROMMASTERTHREAD, 1) #cannot be changed in julia

    load_empty(m)

    return m
end
Model(env; finalize_env::Bool=true) = Model(finalize_env=true)

function Model(name::String, sense::Symbol = :minimize)
    model = Model()
    if sense != :minimize
        set_sense!(model, sense)
    end
    model
end

#################################################
#
#  model error handling
#
#################################################
function get_error_msg(m::Model)
    #@assert env.ptr_env == 1
    out = Array{Cchar}(undef,  512)
    ret = @xprs_ccall(getlasterror, Cint, (Ptr{Nothing},Ptr{Cchar}),
        m.ptr_model, out)

    error( unsafe_string(pointer(out))  )
end
function get_error_msg(m::Model, ret::Int)
    #@assert env.ptr_env == 1
    out = Array{Cint}(undef, 1)
    out2 = @xprs_ccall(getintattrib, Cint,(Ptr{Nothing}, Cint, Ptr{Cint}),
        m.ptr_model, ret , out)
end

##############################################
# Error handling
##############################################
mutable struct XpressError
    code::Int
    msg::String
end
function XpressError(m::Model)#, code::Integer)
    XpressError( 0, get_error_msg(m) )#convert(Int, code), get_error_msg(env))
end
function XpressError(m::Model, ret::Int)
    XpressError( get_error_msg(m, ret), "get message from optimizer manual" )
end

#################################################
#
#  model manipulation
#
#################################################

Base.unsafe_convert(ty::Type{Ptr{Nothing}}, model::Model) = model.ptr_model::Ptr{Nothing}

"""
    free_model(model::Model)

Free all memory allocated in C related to Model
"""
function free_model(model::Model)
    if model.ptr_model != C_NULL
        ret = @xprs_ccall(destroyprob, Cint, (Ptr{Nothing},), model.ptr_model)
        if ret != 0
            throw(XpressError(model))
        end
        model.ptr_model = C_NULL
    end
    return nothing
end

"""
    copy(model_src::Model)

Creates copy of the Xpress Model

Callbacks and setting are not copied
"""
function copy(model_src::Model)
    # only copies problem not callbacks and controls
    if model_src.ptr_model != C_NULL
        model_dest = Model()#env)
        ret = @xprs_ccall(copyprob, Cint, (Ptr{Nothing},Ptr{Nothing},Ptr{UInt8}),
            model_dest.ptr_model, model_src.ptr_model, "")
        if ret != 0
            throw(XpressError(model_src))
        end
    end
    model_dest
end
function copycontrols!(mdest::Model, msrc::Model)
    ret = @xprs_ccall(copycontrols, Cint, (Ptr{Nothing}, Ptr{Nothing}),
    mdest.ptr_model, msrc.ptr_model)
    if ret != 0
        throw(XpressError(mdest))
    end
    nothing
end
function dumpcontrols(model::Model)
    ret = @xprs_ccall(dumpcontrols, Cint, (Ptr{Nothing},), model.ptr_model)
    if ret != 0
        throw(XpressError(mdest))
    end
    nothing
end
function setdefaults(model::Model)
    ret = @xprs_ccall(setdefaults, Cint, (Ptr{Nothing},), model.ptr_model)
    if ret != 0
        throw(XpressError(mdest))
    end
    nothing
end

addcolnames(m::Model, names::Vector) = addnames(m, names, Int32(2))
addrownames(m::Model, names::Vector) = addnames(m, names, Int32(1))
function addnames(m::Model, names::Vector, nametype::Int32)
    # XPRSaddnames(prob, int type, char name[], int first, int last)

    NAMELENGTH = 64

    #nametype = 2
    first = 0
    last = length(names)-1

    cnames = ""
    for str in names
        cnames = string(cnames, join(Base.Iterators.take(str,NAMELENGTH)), "\0")
    end
    ret = @xprs_ccall(addnames, Cint, (Ptr{Nothing}, Cint,Ptr{Cchar}, Cint, Cint),
        m.ptr_model, nametype, cnames, first, last)

    if ret != 0
        throw(XpressError(m))
    end

    nothing
end

# read / write file
#=
XPRSprob a-> Ptr{Nothing}
const char *a -> Ptr{UInt8}  -> "" : String
int -> Cint
const char a[] -> Ptr{Cchar} -> Cchar[]
=#
function load_empty(model::Model)
#=int XPRS_CC XPRSloadlp(XPRSprob prob, const char *probname, int ncol, int
      nrow, const char qrtype[], const double rhs[], const double range[],
      const double obj[], const int mstart[], const int mnel[], const int
      mrwind[], const double dmatval[], const double dlb[], const double
      dub[]);
=#
    ret = @xprs_ccall(loadlp, Cint,
        (Ptr{Nothing},
            Ptr{UInt8},
            Cint,
            Cint,
            Ptr{Cchar},#type
            Ptr{Float64},#rhs
            Ptr{Float64},#range
            Ptr{Float64},#obj
            Ptr{Cint},#mastart
            Ptr{Cint},#mnel
            Ptr{Cint},#mrwind
            Ptr{Float64},#dmat
            Ptr{Float64},#lb
            Ptr{Float64}#ub
            ),model.ptr_model,"",0,0,Cchar[],Float64[],Float64[],Float64[],
                Cint[],Cint[],Cint[],Float64[],Float64[],Float64[]
    )
    if ret != 0
        throw(XpressError(model))
    end
    nothing
end

"""
    read_model(model::Model, filename::String)

Read file and add its informationt o the model
"""
function read_model(model::Model, filename::String)
    #@assert is_valid(model.env)
    flags = ""
    ret = @xprs_ccall(readprob, Cint,
        (Ptr{Nothing}, Ptr{UInt8}, Ptr{UInt8}),
        model.ptr_model, filename, flags)
    if ret != 0
        throw(XpressError(model))
    end
    nothing
end

"""
    write_model(model::Model, filename::String, flags::String="")

Writes a model into file.
For flags setting see the manual (writeprob)
"""
function write_model(model::Model, filename::String, flags::String="")
    ret = @xprs_ccall(writeprob, Cint, (Ptr{Nothing}, Ptr{UInt8}, Ptr{UInt8}),
        model.ptr_model, filename, flags)
    if ret != 0
        throw(XpressError(model))
    end
    nothing
end

"""
    writesol(model::Model, filename::String, flags::String="")

Writes solution into file.
For flags setting see the manual (writesol)
"""
function writesol(model::Model, filename::String, flags::String="")
    # int XPRS_CC XPRSwritesol(XPRSprob prob, const char *filename, const char *flags)
    ret = @xprs_ccall(writesol, Cint, (Ptr{Nothing}, Ptr{UInt8}, Ptr{UInt8}),
        model.ptr_model, filename, flags)
    if ret != 0
        throw(XpressError(model))
    end
    nothing
end

"""
    writeptrsol(model::Model, filename::String, flags::String="")

Writes solution into a ptr file.
For flags setting see the manual (writeptrsol)
"""
function writeptrsol(model::Model, filename::String, flags::String="")
    # int XPRS_CC XPRSwriteptrsol(XPRSprob prob, const char *filename, const char *flags)
    ret = @xprs_ccall(writeptrsol, Cint, (Ptr{Nothing}, Ptr{UInt8}, Ptr{UInt8}),
        model.ptr_model, filename, flags)
    if ret != 0
        throw(XpressError(model))
    end
    nothing
end

"""
    getprobname(model::Model)::String

Query internal model name
"""
function getprobname(model::Model)
    name = " "^100
    ret = @xprs_ccall(getprobname, Cint,
        (Ptr{Nothing},Ptr{Cchar}), model.ptr_model, name)
    if ret != 0
        throw(XpressError(model))
    end
    out = strip(name)
end

"""
    setprobname(model::Model, name::String)::Nothing

Set internal model name
"""
function setprobname(model::Model, name::String)
    ret = @xprs_ccall(setprobname, Cint,
        (Ptr{Nothing},Ptr{Cchar}), model.ptr_model, name)
    if ret != 0
        throw(XpressError(model))
    end
    nothing
end

"""
    writesvf(model::Model)::Nothing

Save serialized internal model to svf format (Xpress internal format).
File is saved to `pwd()`
"""
function writesvf(model::Model)
    ret = @xprs_ccall(save, Cint,
        (Ptr{Nothing},), model.ptr_model)
    if ret != 0
        throw(XpressError(model))
    end
    nothing
end
"""
    readsvf(model::Model, name::String, force::String = "")::Nothing

Read serialized internal model to svf format (Xpress internal format)
"""
function readsvf(model::Model, name::String, force::String = "")
    ret = @xprs_ccall(restore, Cint,
        (Ptr{Nothing}, Ptr{UInt8}, Ptr{UInt8}), model.ptr_model, name, force)
    if ret != 0
        throw(XpressErrocontr(model))
    end
    nothing
end

"""
    fixglobals(model::Model, round::Bool)

Fix globals entities after a first Global search.
"""
function fixglobals(model::Model, round::Bool)
    flag = 0
    if round
        flag = 1
    end
    ret = @xprs_ccall(fixglobals, Cint, (Ptr{Nothing}, Cint),
        model.ptr_model, flag)
    if ret != 0
        throw(XpressError(model))
    end
    nothing
end

"""
    setlogfile(model::Model, filename::String)

Attach a log file to the model
"""
function setlogfile(model::Model, filename::String)
    #int XPRS_CC XPRSsetlogfile(XPRSprob prob, const char *filename);

    flags = ""
    ret = @xprs_ccall(setlogfile, Cint, (Ptr{Nothing}, Ptr{UInt8}),
        model.ptr_model, filename)
    if ret != 0
        throw(XpressError(model))
    end
    nothing
end


text:
z(x) = 3x

h(x) = 5x


text:
"""
    Discontinuity(t, order)

Object of discontinuity of order `order` at time `t`, i.e. discontinuity of `order`th
derivative at time `t`.
"""
struct Discontinuity{tType,O}
    t::tType
    order::O
end

# ordering of discontinuities
Base.:<(a::Discontinuity, b::Discontinuity) = a.t < b.t || (a.t == b.t && a.order < b.order)
Base.isless(a::Discontinuity, b::Discontinuity) = isless(a.t, b.t) || (isequal(a.t, b.t) && isless(a.order, b.order))
Base.isequal(a::Discontinuity, b::Discontinuity) = isequal(a.t, b.t) && isequal(a.order, b.order)
Base.:(==)(a::Discontinuity, b::Discontinuity) = a.t == b.t && a.order == b.order

# ordering with numbers
Base.:<(a::Discontinuity, b::Number) = a.t < b
Base.:<(a::Number, b::Discontinuity) = a < b.t
Base.isless(a::Discontinuity, b::Number) = isless(a.t, b)
Base.isless(a::Number, b::Discontinuity) = isless(a, b.t)
Base.:(==)(a::Discontinuity, b::Number) = a.t == b
Base.:(==)(a::Number, b::Discontinuity) = a == b.t
Base.isequal(a::Discontinuity, b::Number) = isequal(a.t, b)
Base.isequal(a::Number, b::Discontinuity) = isequal(a, b.t)

# multiplication with numbers
Base.:*(a::Discontinuity, b::Number) = Discontinuity(a.t * b, a.order)
Base.:*(a::Number, b::Discontinuity) = Discontinuity(a * b.t, b.order)

# conversion to numbers
Base.convert(::Type{T}, d::Discontinuity) where {T<:Number} = convert(T, d.t)
Base.convert(::Type{T}, d::Discontinuity{T}) where {T<:Number} = d.t

# conversion to discontinuity
function Base.convert(::Type{Discontinuity{T}}, x::Number) where {T<:Number}
    return convert(Discontinuity{T,Int}, x)
end
function Base.convert(::Type{Discontinuity{T,O}}, x::Number) where {T<:Number,O}
    return Discontinuity{T,O}(convert(T, x), zero(O))
end
function Base.convert(::Type{Discontinuity{T,O}}, x::T) where {T<:Number,O}
    return Discontinuity{T,O}(x, zero(O))
end

# display
Base.show(io::IO, d::Discontinuity) = print(io, d.t, " (order ", d.order, ")")
Base.show(io::IO, ::MIME"text/plain", d::Discontinuity) = print(io, typeof(d), ":\n   ", d)

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
export EpsilonGreedyExplorer, GreedyExplorer

using Random
using Distributions: DiscreteNonParametric

"""
    EpsilonGreedyExplorer{T}(;kwargs...)
    EpsilonGreedyExplorer(ϵ) -> EpsilonGreedyExplorer{:linear}(; ϵ_stable = ϵ)

> Epsilon-greedy strategy: The best lever is selected for a proportion `1 - epsilon` of the trials, and a lever is selected at random (with uniform probability) for a proportion epsilon . [Multi-armed_bandit](https://en.wikipedia.org/wiki/Multi-armed_bandit)

Two kinds of epsilon-decreasing strategy are implmented here (`linear` and `exp`).

> Epsilon-decreasing strategy: Similar to the epsilon-greedy strategy, except that the value of epsilon decreases as the experiment progresses, resulting in highly explorative behaviour at the start and highly exploitative behaviour at the finish.  - [Multi-armed_bandit](https://en.wikipedia.org/wiki/Multi-armed_bandit)

# Keywords

- `T::Symbol`: defines how to calculate the epsilon in the warmup steps. Supported values are `linear` and `exp`.
- `step::Int = 1`: record the current step.
- `ϵ_init::Float64 = 1.0`: initial epsilon.
- `warmup_steps::Int=0`: the number of steps to use `ϵ_init`.
- `decay_steps::Int=0`: the number of steps for epsilon to decay from `ϵ_init` to `ϵ_stable`.
- `ϵ_stable::Float64`: the epsilon after `warmup_steps + decay_steps`.

# Example

```julia
s = EpsilonGreedyExplorer{:linear}(ϵ_init=0.9, ϵ_stable=0.1, warmup_steps=100, decay_steps=100)
plot([RL.get_ϵ(s, i) for i in 1:500], label="linear epsilon")
```
![](https://github.com/JuliaReinforcementLearning/ReinforcementLearning.jl/raw/master/docs/src/assets/img/linear_epsilon_greedy_selector.png)

```julia
s = EpsilonGreedyExplorer{:exp}(ϵ_init=0.9, ϵ_stable=0.1, warmup_steps=100, decay_steps=100)
plot([RL.get_ϵ(s, i) for i in 1:500], label="exp epsilon")
```
![](https://github.com/JuliaReinforcementLearning/ReinforcementLearning.jl/raw/master/docs/src/assets/img/exp_epsilon_greedy_selector.png)
"""
mutable struct EpsilonGreedyExplorer{Kind,IsBreakTie,R} <: AbstractExplorer
    ϵ_stable::Float64
    ϵ_init::Float64
    warmup_steps::Int
    decay_steps::Int
    step::Int
    rng::R
end

function Base.copy(p::EpsilonGreedyExplorer{Kind,IsBreakTie,R}) where {Kind,IsBreakTie,R}
    EpsilonGreedyExplorer{Kind,IsBreakTie,R}(
        p.ϵ_stable,
        p.ϵ_init,
        p.warmup_steps,
        p.decay_steps,
        p.step,
        copy(p.rng),
    )
end

function EpsilonGreedyExplorer(;
    ϵ_stable,
    kind = :linear,
    ϵ_init = 1.0,
    warmup_steps = 0,
    decay_steps = 0,
    step = 1,
    is_break_tie = false,
    seed = nothing,
)
    rng = MersenneTwister(seed)
    EpsilonGreedyExplorer{kind,is_break_tie,typeof(rng)}(
        ϵ_stable,
        ϵ_init,
        warmup_steps,
        decay_steps,
        step,
        rng,
    )
end

EpsilonGreedyExplorer(ϵ; kwargs...) = EpsilonGreedyExplorer(; ϵ_stable = ϵ, kwargs...)

function get_ϵ(s::EpsilonGreedyExplorer{:linear}, step)
    if step <= s.warmup_steps
        s.ϵ_init
    elseif step >= (s.warmup_steps + s.decay_steps)
        s.ϵ_stable
    else
        steps_left = s.warmup_steps + s.decay_steps - step
        s.ϵ_stable + steps_left / s.decay_steps * (s.ϵ_init - s.ϵ_stable)
    end
end

function get_ϵ(s::EpsilonGreedyExplorer{:exp}, step)
    if step <= s.warmup_steps
        s.ϵ_init
    else
        n = step - s.warmup_steps
        scale = s.ϵ_init - s.ϵ_stable
        s.ϵ_stable + scale * exp(-1.0 * n / s.decay_steps)
    end
end

get_ϵ(s::EpsilonGreedyExplorer) = get_ϵ(s, s.step)

"""
    (s::EpsilonGreedyExplorer)(values; step) where T

!!! note
    If multiple values with the same maximum value are found.
    Then a random one will be returned!

    `NaN` will be filtered unless all the values are `NaN`.
    In that case, a random one will be returned.
"""
function (s::EpsilonGreedyExplorer{<:Any,true})(values)
    ϵ = get_ϵ(s)
    s.step += 1
    rand(s.rng) >= ϵ ? rand(s.rng, find_all_max(values)[2]) : rand(s.rng, 1:length(values))
end

function (s::EpsilonGreedyExplorer{<:Any,false})(values)
    ϵ = get_ϵ(s)
    s.step += 1
    rand(s.rng) >= ϵ ? find_max(values)[2] : rand(s.rng, 1:length(values))
end

function (s::EpsilonGreedyExplorer{<:Any,true})(values, mask)
    ϵ = get_ϵ(s)
    s.step += 1
    rand(s.rng) >= ϵ ? rand(s.rng, find_all_max(values, mask)[2]) :
    rand(s.rng, findall(mask))
end

function (s::EpsilonGreedyExplorer{<:Any,false})(values, mask)
    ϵ = get_ϵ(s)
    s.step += 1
    rand(s.rng) >= ϵ ? find_max(values, mask)[2] : rand(s.rng, findall(mask))
end

Random.seed!(s::EpsilonGreedyExplorer, seed) = Random.seed!(s.rng, seed)

"""
    get_prob(s::EpsilonGreedyExplorer, values) -> DiscreteNonParametric
    get_prob(s::EpsilonGreedyExplorer, values, mask) -> DiscreteNonParametric

Return the probability of selecting each action given the estimated `values` of each action.
"""
function RLBase.get_prob(s::EpsilonGreedyExplorer{<:Any,true}, values)
    ϵ, n = get_ϵ(s), length(values)
    probs = fill(ϵ / n, n)
    max_val_inds = find_all_max(values)[2]
    for ind in max_val_inds
        probs[ind] += (1 - ϵ) / length(max_val_inds)
    end
    DiscreteNonParametric(1:length(probs), probs)
end

function RLBase.get_prob(s::EpsilonGreedyExplorer{<:Any,false}, values)
    ϵ, n = get_ϵ(s), length(values)
    probs = fill(ϵ / n, n)
    probs[find_max(values)[2]] += 1 - ϵ
    DiscreteNonParametric(1:length(probs), probs)
end

function RLBase.get_prob(s::EpsilonGreedyExplorer{<:Any,true}, values, mask)
    ϵ, n = get_ϵ(s), length(values)
    probs = zeros(n)
    probs[mask] .= ϵ / sum(mask)
    max_val_inds = find_all_max(values, mask)[2]
    for ind in max_val_inds
        probs[ind] += (1 - ϵ) / length(max_val_inds)
    end
    DiscreteNonParametric(1:length(probs), probs)
end

function RLBase.get_prob(s::EpsilonGreedyExplorer{<:Any,false}, values, mask)
    ϵ, n = get_ϵ(s), length(values)
    probs = zeros(n)
    probs[mask] .= ϵ / sum(mask)
    probs[find_max(values, mask)[2]] += 1 - ϵ
    DiscreteNonParametric(1:length(probs), probs)
end

RLBase.reset!(s::EpsilonGreedyExplorer) = s.step = 1

struct GreedyExplorer <: AbstractExplorer end

(s::GreedyExplorer)(values) = find_max(values)[2]
(s::GreedyExplorer)(values, mask) = find_max(values, mask)[2]

function RLBase.get_prob(s::GreedyExplorer, values)
    prob = zeros(length(values))
    prob[find_max(values)[2]] = 1.0
    DiscreteNonParametric(1:length(prob), prob)
end

function RLBase.get_prob(s::GreedyExplorer, values, mask)
    prob = zeros(length(values))
    prob[find_max(values, mask)[2]] = 1.0
    DiscreteNonParametric(1:length(prob), prob)
end


text:
using JSON
using Logging
using ArgParse

include("src/config.jl")
include("src/network.jl")
include("src/datasets.jl")
include("src/classify.jl")

argtable = ArgParseSettings()
@add_arg_table(
    argtable,
    "--seed", arg_type = Int, default = 0,
    "--log", arg_type = String, default = "iris.log",
    "--data", arg_type = String, default = "iris",
    "--supervised", arg_type = Bool, default = true,
    "--n_epochs", arg_type = Int64, default = 1,
    "--connectivity", arg_type = Float64, default = 0.3,
    "--tinput", arg_type = Int64, default = 4,
    "--toutput", arg_type = Int64, default = 4,
    "--trest", arg_type = Int64, default = 7,
    "--ma_rate", arg_type = Float64, default = 0.8,
    "--fr", arg_type = Float64, default = 0.5,
    "--n", arg_type = Int64, default = 12,
    "--nratio", arg_type = Float64, default = 2.8,
)

args, cfg = get_config(argtable)
srand(args["seed"])

Logging.configure(filename=args["log"], level=INFO)
Logging.info(@sprintf("Arguments: %s", JSON.json(args)))
Logging.info(@sprintf("Config: %s", JSON.json(cfg)))

# nin, nout, nn
X, Y = get_data(args["data"])
nin = size(X, 1) * args["n"]
nout = length(unique(Y)) * args["n"]
nn = Int64(round((nin + nout) * args["nratio"]))

n = Network(nn, nin, nout, args["connectivity"], cfg)

pcfg = Dict()
pcfg["n_classes"] = length(unique(Y))
for k in ["tinput", "toutput", "trest"]
    pcfg[k] = args[k] * 10
end
for k in ["ma_rate", "fr", "n_epochs", "seed", "supervised", "data"]
    pcfg[k] = args[k]
end

fit = run_classify(n, X, Y, pcfg)
Logging.info(@sprintf("E%0.6f", -fit))


text:
# Test the default and optional arguments
@test ismeanfeasible(m)
@test ismeanfeasible(m, toler = 1e-8)
@test ismeanfeasible(m, solver = GLPKSolverLP())

# Test an infeasible case
setmean(m, ones(flex_data.numRVs) * -50)
@test !ismeanfeasible(m)

# Reset the mean back to the original value
setmean(m, means)
@test ismeanfeasible(m)


text:
# Visualization of the interference of two point sources in 2D

using PyPlot

const L = 6
const sep = 2
const N = 500
const k = 5

V = linspace(-L, L, N)
X = V
Y = X'
#X, Y = meshgrid(V, V)
R1 = sqrt( (X.-sep).^2 .+ Y.^2)
R2 = sqrt( (X.+sep).^2 .+ Y.^2)

# Sum of Green's functions for two point sources
Z = exp(1im*k*R1)./R1 + exp(1im*k*R2)./R2

M = 10
T = linspace(0.0, 2*pi, M)
T = T[1:(M-1)]
cut = 0.8
scale = 55/(2*cut)

for p=1:1
  for iter=1:length(T)

    figure(1)
    clf()
    hold()

    W = real(Z*exp(-1im * T[iter]))
    W = max(W, -cut)
    W = min(W, cut)

    imshow(scale*(W+cut), vmin=nothing, vmax=nothing, aspect="equal")

    filename = "Frame$(1000+iter)"
    print("Saving to $(filename)")
    savefig(filename)
  end
 end

 # convert to gif
 run(`convert -density 100 -loop 1000 -delay 10 Frame100* interference_two_sources.gif`)


text:
function insert_or_inc(d::Dictionary{Int,Int}, ind::Int)
    set!(d, ind, isassigned(d, ind) ? d[ind] + 1 : 1)
end


text:
######################################################################
# cache2.jl
#
# like infotheoretic.jl, but allows caching
######################################################################

export make_cache
export make_cache2

function make_cache(df::DF)
    n  = 2*df.n - 1
    cache = zeros(n, n, length(df.obs_list))

    for dxi = 1:n, dyi = 1:n
        dx = dxi - df.n
        dy = dyi - df.n

        # generate a theta, xp that fit this
        # just put vehicle at zero?
        xp = (0.0,0.0,0.0)
        tx = dx * df.cell_size
        ty = dy * df.cell_size
        theta = (tx, ty, 0.0, 0.0)

        for (oi,o) in enumerate(df.obs_list)
            cache[dxi, dyi, oi] = O(df.sensor, theta, xp, o)
        end
    end

    return cache
end

function mutual_information(df::DF, cache::Array{Float64,3})
    mut_info = zeros(df.n, df.n)
    for xi = 1:df.n, yi = 1:df.n
        # TODO: really, should have some other heading
        xp = (xi, yi, 0)
        mut_info[xi,yi] = mutual_information(df, xp, cache)
    end
    return mut_info
end

function mutual_information(df::DF, xp::NTuple{3,Int}, cache::Array{Float64,3})
	H_o = 0.0
	H_o_t = 0.0

    for (oi,o) in enumerate(df.obs_list)
        po = 0.0
		for txi = 1:df.n, tyi = 1:df.n
            if df.b[txi, tyi] > 0.0
                pot = cache[txi-xp[1]+df.n, tyi-xp[2]+df.n, oi]
                po += df.b[txi, tyi] * pot
                if pot > 0.0
                    H_o_t -= pot * df.b[txi, tyi] * log(pot)
                end
            end
		end
		if po > 0.0
			H_o -= po * log(po)
		end
	end
	return H_o - H_o_t
end

# for the R2 case
function p_obs(df::DF, xp::NTuple{3,Int}, oi::Int, cache::Array{Float64,3})
    prob = 0.0
    for txi = 1:df.n, tyi = 1:df.n
        if df.b[txi,tyi] > 0.0
            prob += df.b[txi, tyi]*cache[txi-xp[1]+df.n, tyi-xp[2]+df.n, oi]
        end
    end
    return prob
end

######################################################################
# SE(2)
######################################################################

function make_cache2(df::DF)
    n  = 2*df.n - 1
    cache = zeros(n, n, 36, length(df.obs_list))

    for dxi = 1:n, dyi = 1:n
        dx = dxi - df.n
        dy = dyi - df.n

        for hi = 1:36

            h = (hi - 0.5) * 10.0
            # generate a theta, xp that fit this
            # just put vehicle at zero?
            xp = (0.0,0.0,h)
            tx = dx * df.cell_size
            ty = dy * df.cell_size
            theta = (tx, ty, 0.0, 0.0)

            for (oi,o) in enumerate(df.obs_list)
                cache[dxi, dyi, hi, oi] = O(df.sensor, theta, xp, o)
            end
        end
    end

    return cache
end

function make_cache2(df::DF, n_thetas)
    n  = 2*df.n - 1
    cache = zeros(n, n, n_thetas, length(df.obs_list))

    d_theta = 360.0 / n_thetas

    for dxi = 1:n, dyi = 1:n
        dx = dxi - df.n
        dy = dyi - df.n

        for hi = 1:36

            h = (hi - 0.5) * n_thetas
            # generate a theta, xp that fit this
            # just put vehicle at zero?
            xp = (0.0,0.0,h)
            tx = dx * df.cell_size
            ty = dy * df.cell_size
            theta = (tx, ty, 0.0, 0.0)

            for (oi,o) in enumerate(df.obs_list)
                cache[dxi, dyi, hi, oi] = O(df.sensor, theta, xp, o)
            end
        end
    end

    return cache
end

function mutual_information(df::DF, cache::Array{Float64,4})
    mut_info = zeros(df.n, df.n, 36)
    for xi = 1:df.n, yi = 1:df.n, hi = 1:36
    #for hi = 1:36, yi = 1:df.n, xi = 1:df.n
        xp = (xi, yi, hi)
        mut_info[xi,yi, hi] = mutual_information(df, xp, cache)
    end
    return mut_info
end

function mutual_information(df::DF, xp::NTuple{3,Int}, cache::Array{Float64,4})
	H_o = 0.0
	H_o_t = 0.0

    for (oi,o) in enumerate(df.obs_list)
        po = 0.0

		# sum over possible jammer locations
		for txi = 1:df.n, tyi = 1:df.n
            if df.b[txi, tyi] > 0.0
                pot = cache[txi-xp[1]+df.n, tyi-xp[2]+df.n, xp[3], oi]
                po += df.b[txi,tyi] * pot
                if pot > 0.0
                    H_o_t -= pot * df.b[txi, tyi] * log(pot)
                end
            end
		end
		if po > 0.0
			H_o -= po * log(po)
		end
	end
	return H_o - H_o_t
end

function p_obs(df::DF, xp::NTuple{3,Int}, oi::Int, cache::Array{Float64,4})
    prob = 0.0
    for txi = 1:df.n, tyi = 1:df.n
        if df.b[txi,tyi] > 0.0
            dxi = txi - xp[1] + df.n    # cache x index
            dyi = tyi - xp[2] + df.n    # cache y index
            prob += df.b[txi, tyi]*cache[dxi, dyi, xp[3], oi]
        end
    end
    return prob
end


text:
export AtomicMeasure, WeightedDiracMeasure

"""
    struct WeightedDiracMeasure{T}
        center::Vector{T}
        weight::T
    end

Represents a weighted dirac measure with centered at `center` and with weight `weight`.
"""
struct WeightedDiracMeasure{T, AT<:AbstractVector{T}} # Cannot subtype AbstractMeasureLike, since it does not contain the variables corresponding to the coefficients of the center
    center::AT
    weight::T
end
function WeightedDiracMeasure(center::AbstractVector{S}, weight::T) where {S, T}
    U = promote_type(S, T)
    return WeightedDiracMeasure(convert(AbstractVector{U}, center),
                                convert(U, weight))
end

"""
    struct AtomicMeasure{T, AT, V} <: AbstractMeasureLike{T}
        variables::V                           # Vector/Tuple of variables
        atoms::Vector{WeightedDiracMeasure{T, AT}} # Atoms of the measure
    end

An measure is said to be *atomic* if it is a sum of weighted dirac measures.
For instance, ``\\eta = 2 \\delta_{(1, 0)} + 3 \\delta_{(1/2, 1/2)}`` is an atomic measure since it is a sum of the diracs centered at ``(1, 0)`` and ``(1/2, 1/2)`` and weighted respectively by 2 and 3.
That is, ``\\mathbb{E}_{\\eta}[p] = 2p(1, 0) + 3p(1/2, 1/2)``.

The `AtomicMeasure` struct stores an atomic measure where `variables` contains the variables and `atoms` contains atoms of the measure.
"""
struct AtomicMeasure{T, AT, V} <: AbstractMeasureLike{T}
    variables::V                           # Vector/Tuple of variables
    atoms::Vector{WeightedDiracMeasure{T, AT}} # Atoms of the measure
end

function Base.show(io::IO, η::AtomicMeasure)
    print(io, "Atomic measure on the variables $(join(η.variables, ", ")) with $(length(η.atoms)) atoms:")
    for atom in η.atoms
        println(io)
        print(io, " at $(atom.center) with weight $(atom.weight)")
    end
end

measure(η::AtomicMeasure, x) = Measure(η, x)
function Measure(η::AtomicMeasure{T}, x::AbstractVector{TT}) where {T, TT}
    Measure{T, MB.MonomialBasis{monomialtype(TT), monovectype(x)}}(η, x)
end
function Measure{T, BT}(η::AtomicMeasure{T}, x) where {T, BT}
    X = monovec(x)
    sum(atom.weight * dirac(X, η.variables => atom.center) for atom in η.atoms)
end

expectation(η::AtomicMeasure, p::APL) = sum(δ -> δ.weight * p(η.variables => δ.center), η.atoms)
expectation(p::APL, η::AtomicMeasure) = expectation(η, p)

function permcomp(f, m)
    picked = BitSet()
    for i in 1:m
        k = 0
        for j in 1:m
            if !(j in picked) && f(i, j)
                k = j
                break
            end
        end
        if k == 0
            return false
        end
        push!(picked, k)
    end
    true
end

Base.isapprox(η1::WeightedDiracMeasure, η2::WeightedDiracMeasure; kws...) = isapprox(η1.weight, η2.weight; kws...) && isapprox(η1.center, η2.center; kws...)
function Base.isapprox(η1::AtomicMeasure, η2::AtomicMeasure; kws...)
    m = length(η1.atoms)
    length(η2.atoms) == m && permcomp((i, j) -> isapprox(η1.atoms[i], η2.atoms[j]; kws...), m)
end


text:
using Plots, DataFrames, CSV, LinearAlgebra
include("utility.jl")

const G = 100.0
const S = 0.05

mutable struct Particle
    m::Float64
    r::Vector{Float64}
    v::Vector{Float64}
end

mutable struct System
    particles::Vector{Particle}
    #cutoff_distance::Float64
    #a::Float64
    #b::Float64
    #c::Float64
    #alpha::Float64
    #beta::Float64
    #gamma::Float64
end

# returns acceleration of p1 due to the force between p1 and p2
function a_g(p1::Particle, p2::Particle)
    r = p2.r - p1.r
    return G * p2.m * r / (norm(r)^2 + S^2)^1.5
end

# returns the net acceleration of a particle due to the forces from all other
# particles within the cutoff distance
function net_a_g(p1::Particle, s::System)
    a = [0.0, 0.0, 0.0]

    for p2 in s.particles
        if p1 != p2
            a += a_g(p1, p2)
        end
    end

    return a
end

# move system s forward by dt using velocity-verlet method
function step(s::System, dt::Float64)
    for p in s.particles
        a1 = net_a_g(p, s)
        p.r += p.v * dt + 0.5 * a1 * dt^2
        a2 = net_a_g(p, s)
        p.v += 0.5 * (a1 + a2) * dt
    end
end

# numerically integrates system particle positions from 0 to t
function integrate(s::System, t::Float64, dt::Float64, save_interval::Float64)
    rows = Int(t / save_interval)
    columns = 3 * length(s.particles) + 1
    data = Array{Float64}(undef, rows, columns)
    data[1, :] = vcat([0.0], [s.particles[j].r[k] for
        j = 1:length(s.particles) for k = 1:3])
    save_counter = 0.0
    i = 2

    for t_i in 0:dt:t
        step(s, dt)

        if save_counter >= save_interval
            data[i, :] = vcat([t_i], [s.particles[j].r[k] for
                j = 1:length(s.particles) for k = 1:3])
            i += 1
            save_counter = 0.0
        else
            save_counter += dt
        end
    end

    data = data[1:i - 1, :]
    return data
end

function main()
    system, t, dt, save_interval = load("input.txt")
    data = integrate(system, t, dt, save_interval)
    save(data, "output")
    create_plot(data, "plot")
end

main()


text:
export setϵᵣ! ,  setμᵣ! , setmaterial! , setbackground!
export Material
export ϵ₁ , ϵ₂ , n , κ

function  setϵᵣ! end
function  setμᵣ! end
function  setmaterial! end
function  setbackground! end

struct Material{T1,T2,T3,T4,T5,T6}
    ϵr_xx::T1
    ϵr_yy::T2
    ϵr_zz::T3
    μr_xx::T4
    μr_yy::T5
    μr_zz::T6
end

function Material( ; ϵᵣ::Number)
    Material(ϵᵣ,ϵᵣ,ϵᵣ,1,1,1)
end

function (sim::Simulation)(::typeof(setmaterial!), mat::Material, region)
    sim(setϵᵣ!,x̂,x -> mat.ϵr_xx,region)
    sim(setϵᵣ!,ŷ,x -> mat.ϵr_yy,region)
    sim(setϵᵣ!,ẑ,x -> mat.ϵr_zz,region)

    sim(setμᵣ!,x̂,x -> mat.μr_xx,region)
    sim(setμᵣ!,ŷ,x -> mat.μr_yy,region)
    sim(setμᵣ!,ẑ,x -> mat.μr_zz,region)
end

function (sim::Simulation)(::typeof(setbackground!), mat::Material)
    sim(setϵᵣ!,x̂,x -> mat.ϵr_xx,EVERYWHERE)
    sim(setϵᵣ!,ŷ,x -> mat.ϵr_yy,EVERYWHERE)
    sim(setϵᵣ!,ẑ,x -> mat.ϵr_zz,EVERYWHERE)

    sim(setμᵣ!,x̂,x -> mat.μr_xx,EVERYWHERE)
    sim(setμᵣ!,ŷ,x -> mat.μr_yy,EVERYWHERE)
    sim(setμᵣ!,ẑ,x -> mat.μr_zz,EVERYWHERE)
end

function get_ϵ_comp(sim::Simulation{Dim},  dir::Direction{D}) where {D,Dim}
	if D == 1 return sim.ϵᵣ_xx end
	if D == 2 return sim.ϵᵣ_yy end
	if D == 3 return sim.ϵᵣ_zz end
end

function get_μ_comp(sim::Simulation{Dim},  dir::Direction{D}) where {D,Dim}
	if D == 1 return sim.μᵣ_xx end
	if D == 2 return sim.μᵣ_yy end
	if D == 3 return sim.μᵣ_zz end
end

function (sim::Simulation)(::typeof(setϵᵣ!),dir::Direction,val::Function, reg = EVERYWHERE , gridtype::GridType= p̂)
sim(set!, get_ϵ_comp(sim,  dir),val,reg,gridtype)
end

function (sim::Simulation)(::typeof(setμᵣ!),dir::Direction,val::Function, reg = EVERYWHERE , gridtype::GridType = d̂)
sim(set!, get_μ_comp(sim,  dir),val,reg,gridtype)
end

ϵ₁(n::Real,κ::Real) = n^2 - κ^2
ϵ₂(n::Real,κ::Real) = 2*n*κ

n(ϵ₁::Real,ϵ₂::Real) = 1/sqrt(2)  * sqrt(ϵ₁  + sqrt(ϵ₁^2 + ϵ₂^2))
κ(ϵ₁::Real,ϵ₂::Real) = 1/sqrt(2)  * sqrt(-ϵ₁ + sqrt(ϵ₁^2 + ϵ₂^2))


text:
# auxiliary model: HAR-RV
# Corsi, Fulvio. "A simple approximate long-memory model
# of realized volatility." Journal of Financial Econometrics 7,
# no. 2 (2009): 174-196.
function HAR(y)
    ylags = lags(y,10)
    X = [ones(size(y,1)) ylags[:,1]  mean(ylags[:,1:4],dims=2) mean(ylags[:,1:10],dims=2)]
    # drop missings
    y = y[11:end]
    X = X[11:end,:]
    βhat = X\y
    σhat = std(y-X*βhat)
    vcat(βhat,σhat)
end

max_stars_repo_name: mcreel/SV

max_issues_repo_name: mcreel/SV

max_forks_repo_name: mcreel/SV


text:
# test GenericWrapParam
using Base: Test
using Distributions
using NLsolve
using KernelDensityEstimate
using IncrementalInference

import IncrementalInference: getSample

mutable struct FunctorWorks
  a::Array{Float64,2}
end

@testset "FunctorWorks" begin

function (fw::FunctorWorks)(x)
  fw.a[1,1] = -1.0
  nothing
end

A = rand(2,3)
At = deepcopy(A)
At[1,1] = -1.0

fvar = FunctorWorks(A)
fvar(0.0)
@test At == A

end

mutable struct FunctorArray{T}
  # This was a tuple, but array will like work better in the long term
  fnc!::Function
  a::Array{T, 1}
end

@testset "FunctorArray" begin

function testarray!(a1::Array{Float64, 2}, a2::Array{Float64,2})
  a1[1,1] = -1.0
  @show a1
  nothing
end

function (fi::FunctorArray)(x)
  fi.fnc!(fi.a...)
end

A = rand(2,3)
B = rand(2,3)
t = Array{Array{Float64,2},1}()
push!(t,A)
push!(t,B)
At = deepcopy(A)
At[1,1] = -1.0

fvar = FunctorArray(testarray!, t)
fvar([0.0])
@test At == A

end

struct Test2 <: FunctorInferenceType
end

@testset "CommonConbWrapper test" begin

function (tt::Test2)(res::Vector{Float64}, userdata::FactorMetadata, idx::Int, meas::Tuple{Array{Float64,2}}, tp1::Array{Float64,2}, tp2::Array{Float64,2})
  tp1[1,1]=-2.0;
  res[:] = 1.0
  nothing;
end

N = 3
tst2 = Test2()
A = rand(2,N)
B = rand(2,N)
At = deepcopy(A)
t = Array{Array{Float64,2},1}()
push!(t,A)
push!(t,B)
t[1][1,1] = -10.0
@test A[1,1] == -10
# @show typeof(t)
# generalwrapper = GenericWrapParam{Test2}(tst2, t, 1, 1)
ccw = CommonConvWrapper(tst2, t[1], 2, t) # generalwrapper.measurement = rand(1,1)

# x, res = zeros(2), zeros(2)
# @show t[1][:,1]
# # @time generalwrapper(res, x)
#
# At[1,1] = -2.0
# At[2,1] = 0.0
# @test A == At

ccw.cpt

ccw.cpt[Threads.threadid()].factormetadata
ccw.cpt[Threads.threadid()].particleidx
ccw.measurement = (randn(1,N),)
ccw.params[ccw.cpt[Threads.threadid()].activehypo]

x, res = zeros(2), zeros(2)
@time ccw(res, x)

At[1,1] = -2.0
At[2,1] = 0.0
@test A == At

# resource requirements

# fm = FactorMetadata()
# meas = (randn(2,3), )
#
# @time tst2(res, fm, 1, meas, A, B)

end

# abstract Nonparametric <: Function
# This is what the internmediate user would be contributing
mutable struct Pose1Pose1Test{T} <: FunctorPairwise
  Dx::T
  Pose1Pose1Test{T}() where T = new()
  Pose1Pose1Test{T}(a::T) where T = new(a)
  # Pose1Pose1Test(a::T) where T = new(a)
end
Pose1Pose1Test(a::T) where T = Pose1Pose1Test{T}(a)

getSample{T}(pp1t::Pose1Pose1Test{T}, N::Int=1) = (reshape(rand(pp1t.Dx,N),1,N),)

#proposed standardized parameter list, does not have to be functor
function (Dp::Pose1Pose1Test)(res::Array{Float64},
            userdata::FactorMetadata,
            idx::Int,
            meas::Tuple{Array{Float64,2}},
            p1::Array{Float64,2},
            p2::Array{Float64,2} )
  #
  res[1] = meas[1][1,idx] - (p2[1,idx] - p1[1,idx])
  nothing
end

@testset "Test in factor graph setting..." begin

N = 101
p1 = rand(1,N)
p2 = rand(1,N)
t = Array{Array{Float64,2},1}()
push!(t,p1)
push!(t,p2)

odo = Pose1Pose1Test{Normal}(Normal(100.0,1.0))
# generalwrapper = GenericWrapParam{Pose1Pose1Test}(odo, t, 1, 1, getSample(odo, N), getSample)
# generalwrapper.measurement = getSample(odo, N)

# x, res = zeros(1), zeros(1)
# @time for generalwrapper.particleidx in 1:N
#     # nlsolve( generalwrapper, x0? .. )
#   generalwrapper(res, x)
#   # each point should be near 100.0
#   @test res[1] > 50.0
# end

# println("Test with NLsolve for root finding using generalwrapper functor.")
#
# generalwrapper.varidx = 2
# @time for i in 1:N
#   generalwrapper.particleidx = i
#     # generalwrapper(x, res)
#   r = nlsolve( generalwrapper, generalwrapper.params[generalwrapper.varidx][:,generalwrapper.particleidx] )
#   generalwrapper.params[generalwrapper.varidx][1,generalwrapper.particleidx] = r.zero[1]
# end
#
# @test abs(Base.mean(p1)-0.0) < 4.0
# @test abs(Base.mean(p2)-100.0) < 4.0

## Repeat for new CommonConvWrapper

N = 101
p1 = rand(1,N)
p2 = rand(1,N)
t = Array{Array{Float64,2},1}()
push!(t,p1)
push!(t,p2)

odo = Pose1Pose1Test{Normal}(Normal(100.0,1.0))
ccw = CommonConvWrapper(odo, t[1], 1, t, measurement=getSample(odo, N))

freshSamples!(ccw, N)
x, res = zeros(1), zeros(1)

@time for n in 1:N
  ccw.cpt[Threads.threadid()].particleidx = n
    # nlsolve( generalwrapper, x0? .. )
  ccw(res, x)
  # each point should be near 100.0
  @test res[1] > 50.0
end

ccw.varidx = 2
@time for i in 1:N
  ccw.cpt[Threads.threadid()].particleidx = i
    # ccw(x, res)
  r = nlsolve( ccw, ccw.params[ccw.varidx][:,ccw.cpt[Threads.threadid()].particleidx] )
  ccw.params[ccw.varidx][1,ccw.cpt[Threads.threadid()].particleidx] = r.zero[1]
end

@test abs(Base.mean(p1)-0.0) < 4.0
@test abs(Base.mean(p2)-100.0) < 4.0

end

@testset "Test with CommonConvWrapper for un-permuted root finding..." begin

N = 110
p1 = rand(1,N)
p2 = rand(1,N)
t = Array{Array{Float64,2},1}()
push!(t,p1)
push!(t,p2)

odo = Pose1Pose1Test(Normal(100.0,1.0))
# varidx=2 means we are solving for p2 relative to p1

measurement = getSample(odo, N)
@show zDim = size(measurement[1],1)

solvefor = 2

ccw = CommonConvWrapper(odo, t[solvefor], zDim, t, measurement=measurement)
@show ccw.varidx = solvefor
# gwp = GenericWrapParam{Pose1Pose1Test}(odo, t, 2, 1, (zeros(0,1),) , getSample) #getSample(odo, N)

freshSamples!(ccw, N)
# and return complete fr/gwp
@time for n in 1:N
  # gwp(x, res)
  ccw.cpt[Threads.threadid()].particleidx = n
  numericRootGenericRandomizedFnc!( ccw )
end

# @show gwp.params

@test 90.0 < Base.mean(ccw.params[ccw.varidx]) < 110.0
@test -10.0 < Base.mean(ccw.params[1]) < 10.0

println("and in the reverse direction, achieved by simply changing CommonConvWrapper.varidx to 1...")

solvefor = 1
@show ccw.varidx = solvefor
ccw.params[solvefor][:,:] = -100.0*ones(size(ccw.params[solvefor]))
ccw.cpt[Threads.threadid()].X = ccw.params[solvefor]

# @show gwp.params

# fr = FastRootGenericWrapParam{Pose1Pose1Test}(gwp.params[gwp.varidx], zDim, gwp)

freshSamples!(ccw, N)
@time for n in 1:N
  # gwp(x, res)
  ccw.cpt[Threads.threadid()].particleidx = n
  numericRootGenericRandomizedFnc!( ccw )
end

# @show gwp.params

@test -10.0 < Base.mean(ccw.params[1]) < 10.0
@test 90.0 < Base.mean(ccw.params[2]) < 110.0

end

@testset "Test with FastRootGenericWrapParam for permuted root finding..." begin

warn("test not implemented yet")

end

# use the range only example, should give a circle with nothing in the middle

@testset "Generic convolution testing in factor graph context..." begin

N=101
p1 = randn(1,N)
d1 = kde!(p1)
p2 = randn(1,N)
t = Array{Array{Float64,2},1}()
push!(t,p1)
push!(t,p2)

fg = emptyFactorGraph()
# fg.registeredModuleFunctions[:Main] = getSample

v1=addNode!(fg, :x1, ContinuousScalar, N=N)
v2=addNode!(fg, :x2, ContinuousScalar, N=N)
bws = getBW(d1)[:,1]
f1 = addFactor!(fg, [v1], Obsv2(p1, reshape(bws, 1, length(bws)), [1.0]))

odo = Pose1Pose1Test(Normal(100.0,1.0))
f2 = addFactor!(fg, [v1;v2], odo)

tree = wipeBuildNewTree!(fg)

pts = getVal(getVert(fg,:x1))
@test abs(Base.mean(pts)-0.0) < 10.0
pts = getVal(getVert(fg,:x2))
@test abs(Base.mean(pts)-0.0) < 10.0

inferOverTreeR!(fg, tree, N=N)
inferOverTreeR!(fg, tree)
# @time [inferOverTreeR!(fg, tree, N=N) for i in 1:3];

# using Gadfly
# plot(y=rand(10))
# plotKDE(getVertKDE(fg,:x2))

pts = getVal(getVert(fg,:x1))
@test abs(Base.mean(pts)-0.0) < 10.0

pts = getVal(getVert(fg,:x2))
@test abs(Base.mean(pts)-100.0) < 10.0

end


text:
using LightGraphs, SimpleWeightedGraphs, GraphPlot

g = SimpleWeightedGraph(3)

add_edge!(g, 1, 2, 0.5)
add_edge!(g, 1, 3, 2.5)
add_edge!(g, 2, 3, 0.8)

gplot(g, nodelabel=vertices(g), edgelabel=weight.(collect(edges(g))))


text:
module Frameworks

using Printf: @printf
using TimerOutputs: TimerOutputs,TimerOutput,@timeit
using ..Terms: Parameters
using ...Prerequisites: decimaltostr
using ...Prerequisites.CompositeStructures: NamedContainer
using ...Prerequisites.TypeTraits: efficientoperations

import ...Interfaces: id,update!,prepare!,register!,run!,add!

export App,Engine
export Assignment,Algorithm
export dependences,rundependences!

"""
    App

Abstract type for all apps.
"""
abstract type App end

"""
    ==(app1::App,app2::App) -> Bool
    isequal(app1::App,app2::App) -> Bool

Judge whether two apps are equivalent to each other.
"""
Base.:(==)(app1::App,app2::App)= ==(efficientoperations,app1,app2)
Base.isequal(app1::App,app2::App)=isequal(efficientoperations,app1,app2)

"""
    update!(app::App;kwargs...) -> App

Update the status of an app.
"""
update!(app::App;kwargs...)=app

"""
    Engine

Abstract type for all engines.
"""
abstract type Engine end

"""
    ==(engine1::Engine,engine2::Engine) -> Bool
    isequal(engine1::Engine,engine2::Engine) -> Bool

Judge whether two engines are equivalent to each other.
"""
Base.:(==)(engine1::Engine,engine2::Engine)= ==(efficientoperations,engine1,engine2)
Base.isequal(engine1::Engine,engine2::Engine)=isequal(efficientoperations,engine1,engine2)

"""
    update!(engine::Engine;kwargs...) -> Engine

Update the status of an engine.
"""
update!(engine::Engine;kwargs...)=engine

"""
    Assignment( id::Symbol,app::App,parameters::Parameters;
                map::Function=identity,
                dependences::Tuple{Vararg{Symbol}}=(),
                data::Any=nothing,
                savedata::Bool=true,
                virgin::Bool=true,
                kwargs...
                )

An assignment associated with an app.
"""
mutable struct Assignment{A<:App,P<:Parameters,M<:Function,D<:Tuple{Vararg{Symbol}},R<:Any,I}
    app::A
    parameters::P
    map::M
    dependences::D
    data::R
    savedata::Bool
    virgin::Bool
end
function Assignment(id::Symbol,app::App,parameters::Parameters;
                    map::Function=identity,
                    dependences::Tuple{Vararg{Symbol}}=(),
                    data::Any=nothing,
                    savedata::Bool=true,
                    virgin::Bool=true,
                    kwargs...
                    )
    A,P,M,D=app|>typeof,parameters|>typeof,map|>typeof,dependences|>typeof
    R=data===nothing ? Any : data|>typeof
    return Assignment{A,P,M,D,R,id}(app,parameters,map,dependences,data,savedata,virgin)
end

"""
    ==(assign1::Assignment,assign2::Assignment) -> Bool
    isequal(assign1::Assignment,assign2::Assignment) -> Bool

Judge whether two assignments are equivalent to each other.
"""
Base.:(==)(assign1::Assignment,assign2::Assignment)= ==(efficientoperations,assign1,assign2)
Base.isequal(assign1::Assignment,assign2::Assignment)=isequal(efficientoperations,assign1,assign2)

"""
    valtype(assign::Assignment)
    valtype(::Type{<:Assignment{<:App,<:Parameters,<:Function,<:Tuple{Vararg{Symbol}},R}}) where R

The type of the data(result) of an assignment.
"""
Base.valtype(assign::Assignment)=assign|>typeof|>valtype
Base.valtype(::Type{<:Assignment{<:App,<:Parameters,<:Function,<:Tuple{Vararg{Symbol}},R}}) where R=R

"""
    id(assign::Assignment) -> Symbol
    id(::Type{<:Assignment{<:App,<:Parameters,<:Function,<:Tuple{Vararg{Symbol}},<:Any,I}}) where I -> Symbol

The id of an assignment.
"""
id(assign::Assignment)=assign|>typeof|>id
id(::Type{<:Assignment{<:App,<:Parameters,<:Function,<:Tuple{Vararg{Symbol}},<:Any,I}}) where I=I

"""
    update!(assign::Assignment;kwargs...) -> Assignment

Update the parameters of an assignment and the status of its associated app.
"""
@generated function update!(assign::Assignment;kwargs...)
    exprs=[]
    names=fieldnames(fieldtype(assign,:parameters))
    for (i,name) in enumerate(names)
        name=QuoteNode(name)
        push!(exprs,:(get(kwargs,$name,getfield(assign.parameters,$i))))
    end
    return quote
        assign.parameters=Parameters{$names}($(exprs...))
        update!(assign.app;assign.map(assign.parameters)...)
        return assign
    end
end

"""
    Algorithm(  name::String,engine::Engine;
                din::String=".",
                dout::String=".",
                parameters::Union{Parameters,Nothing}=nothing,
                map::Function=identity,
                assignments::Tuple{Vararg{Assignment}}=(),
                kwargs...
                )

An algorithm associated with an engine.
"""
mutable struct Algorithm{E<:Engine,P<:Parameters,M<:Function,S<:NamedContainer{Assignment}}
    name::String
    engine::E
    din::String
    dout::String
    parameters::P
    map::M
    sassignments::S
    dassignments::Dict{Symbol,Assignment}
    timer::TimerOutput
end
function Algorithm( name::String,engine::Engine;
                    din::String=".",
                    dout::String=".",
                    parameters::Union{Parameters,Nothing}=nothing,
                    map::Function=identity,
                    assignments::Tuple{Vararg{Assignment}}=(),
                    kwargs...
                    )
    parameters===nothing && (parameters=Parameters(engine))
    assignments=namedassignments(assignments)
    return Algorithm(name,engine,din,dout,parameters,map,assignments,Dict{Symbol,Assignment}(),TimerOutput())
end
@generated function namedassignments(assignments::Tuple{Vararg{Assignment}})
    names,values=[],[]
    for i=1:fieldcount(assignments)
        push!(names,fieldtype(assignments,1)|>id)
        push!(values,:(assignments[$i]))
    end
    names=NTuple{fieldcount(assignments),Symbol}(names)
    values=Expr(:tuple,values...)
    return :(NamedContainer{$names}($values))
end

"""
    update!(alg::Algorithm;kwargs...) -> Algorithm

Update the parameters of an algorithm and its associated engine.
"""
@generated function update!(alg::Algorithm;kwargs...)
    exprs=[]
    names=fieldnames(fieldtype(alg,:parameters))
    for (i,name) in enumerate(names)
        name=QuoteNode(name)
        push!(exprs,:(get(kwargs,$name,getfield(alg.parameters,$i))))
    end
    return quote
        alg.parameters=Parameters{$names}($(exprs...))
        update!(alg.engine;alg.map(alg.parameters)...)
        return alg
    end
end

"""
    repr(alg::Algorithm,mask::Tuple{Vararg{Symbol}}=();ndecimal::Int=10) -> String

Get the repr representation of an algorithm.

Optionally, some parameters of the algorithm can be masked. Besides, the maximum number of decimals of the parameters can also be specified.
"""
function Base.repr(alg::Algorithm,mask::Tuple{Vararg{Symbol}}=();ndecimal::Int=10)
    result=[string(alg.name),repr(alg.engine)]
    for (name,value) in pairs(alg.parameters)
        name∉mask && push!(result,decimaltostr(value,ndecimal))
    end
    return join(result,"_")
end

"""
    show(io::IO,alg::Algorithm)

Show an algorithm.
"""
function Base.show(io::IO,alg::Algorithm)
    @printf io "%s_%s" alg.name alg.engine
    for (name,value) in pairs(alg.parameters)
        @printf io "_%s" decimaltostr(value,10)
    end
end

"""
    get(alg::Algorithm,id::Symbol) -> Assignment
    get(alg::Algorithm,::Val{id}) where id -> Assignment

Find the assignment registered on a algorithm by its id.
"""
Base.get(alg::Algorithm,id::Symbol)=id∈keys(alg.sassignments) ? getfield(alg.sassignments,id) : alg.dassignments[id]
@generated Base.get(alg::Algorithm,::Val{id}) where id=id∈fieldnames(fieldtype(alg,:sassignments)) ? :(getfield(alg.sassignments,id)) : :(alg.dassignments[id])

"""
    summary(alg::Algorithm)

Provide a summary of an algorithm.
"""
function Base.summary(alg::Algorithm)
    @info "Summary of $(alg.name)($(nameof(typeof(alg.engine)))):"
    @info string(alg.timer)
end

"""
    prepare!(alg::Algorithm,assign::Assignment) -> Nothing

Prepare an assignment registered on a algorithm.
"""
prepare!(alg::Algorithm,assign::Assignment)=nothing

"""
    run!(alg::Algorithm,assign::Assignment) -> Nothing

Run an assignment registered on a algorithm.
"""
run!(alg::Algorithm,assign::Assignment)=nothing

"""
    register!(alg::Algorithm,id::Symbol,app::App;kwargs...) -> Algorithm

Add an assignment on a algorithm by providing the contents of the assignment, and run this assignment.
"""
function register!(alg::Algorithm,id::Symbol,app::App;kwargs...)
    add!(alg,id,app;kwargs...)
    run!(alg,id,true)
end

"""
    add!(alg::Algorithm,id::Symbol,app::App;kwargs...) -> Algorithm

Add an assignment on a algorithm by providing the contents of the assignment.

The difference between `add!` and `register!` is that the `add!` function does not run the newly added assignment but the `register!` function does.
"""
function add!(alg::Algorithm,id::Symbol,app::App;kwargs...)
    @assert id ∉ keys(alg.sassignments) "add! error: id($id) conflict."
    alg.dassignments[id]=Assignment(id,app,merge(alg.parameters,get(kwargs,:parameters,Parameters{()}()));kwargs...)
    return alg
end

"""
    run!(alg::Algorithm,id::Symbol,timing::Bool=true) -> Algorithm
    run!(alg::Algorithm,::Val{id},timing::Bool=true) where id -> Algorithm

Run an assignment with the given id registered on an algorithm. Optionally, the run process can be timed by setting the `timing` argument to be `true`.
"""
run!(alg::Algorithm,id::Symbol,timing::Bool=true)=run!(alg,Val(id),timing)
function run!(alg::Algorithm,::Val{id},timing::Bool=true) where id
    assign=get(alg,Val(id))
    if timing
        @timeit alg.timer string(id) algrunassignment!(alg,assign)
        @info "App $id($(nameof(assign.app|>typeof))): time consumed $(TimerOutputs.time(alg.timer[string(id)])/10^9)s."
    else
        algrunassignment!(alg,assign)
    end
    return alg
end
function algrunassignment!(alg::Algorithm,assign::Assignment)
    ismatched=match(assign.parameters,alg.parameters)
    !assign.virgin && ismatched && return
    !ismatched && update!(alg;assign.parameters...)
    prepare!(alg,assign)
    run!(alg,assign)
    assign.virgin=false
end

"""
    dependences(alg::Algorithm,assign::Assignment,::Tuple{}=()) -> Tuple{Vararg{Symbol}}
    dependences(alg::Algorithm,assign::Assignment,mask::Tuple{Vararg{Symbol}}) -> Tuple{Vararg{Symbol}}

Get the dependences of an assignment and return their ids.
"""
dependences(alg::Algorithm,assign::Assignment,::Tuple{}=())=assign.dependences
dependences(alg::Algorithm,assign::Assignment,mask::Tuple{Vararg{Symbol}})=Tuple(filter(x->x∉mask,collect(assign.dependences)))

"""
    rundependences!(alg::Algorithm,assign::Assignment,mask::Tuple{Vararg{Symbol}}=()) -> Algorithm

Run the dependences of an assignment. Optionally, some dependences can be jumped by specifying the `mask` argument.
"""
function rundependences!(alg::Algorithm,assign::Assignment,mask::Tuple{Vararg{Symbol}}=())
    for id in dependences(alg,assign,mask)
        assign=get(alg,id)
        algrundependence!(alg,assign)
    end
    return alg
end
function algrundependence!(alg::Algorithm,assign::Assignment)
    ismatched=match(alg.parameters,assign.parameters)
    !assign.virgin && ismatched && return
    !ismatched && update!(assign;alg.parameters...)
    prepare!(alg,assign)
    run!(alg,assign)
    assign.virgin=false
end

end  #module