text:
#auxiliary function to evaluate pl.eph. case: no jet transport
function evaleph(eph::TaylorInterpolant, t::Taylor1, ::Taylor1{T}) where {T<:Real}
    return eph(t)
end

# #auxiliary function to evaluate pl.eph. case: Taylor1 jet transport
function evaleph(eph::TaylorInterpolant, t::Taylor1, q1::Taylor1{Taylor1{T}}) where {T<:Real}
    return map(x->Taylor1( x.coeffs*one(q1[0]) ), eph(t))
end

#auxiliary function to evaluate pl.eph. case: TaylorN jet transport
function evaleph(eph::TaylorInterpolant, t::Taylor1, q1::Taylor1{TaylorN{T}}) where {T<:Real}
    return map(x->Taylor1( x.coeffs*one(q1[0]) ), eph(t))
end

#auxiliary function to evaluate pl.eph. case: Lyapunov spectrum
function evaleph(eph::TaylorInterpolant, t::Taylor1, q1::TaylorN{Taylor1{T}}) where {T<:Real}
    return one(q1)*eph(t)
end

function auxzero(a::AbstractSeries)
    return zero(a)
end

function auxzero(a::TaylorN{Taylor1{T}}) where {T<:Number}
    return TaylorN(zero.(a.coeffs))
end

# Nearth-Earth asteroid dynamical model (d=2.0)
# Bodies considered in the model are: the Sun, the eight planets, the Moon and Ceres,
# as well as the asteroid of interest as a test particle with null mass. Dynamical
# effects considered are:
# - post-Newtonian point-mass accelerations between all bodies,
# - figure-effects (oblateness) of the Earth (J2 and J3)
# - J2 effect of the Sun
# - J2 and J3 effect of the Moon
# - Kinematic model for the precession and nutation of the Earth's orientation (IAU 1976/1980 Earth orientation model)
# - Kinematic model for the Moons's orientation (Seidelmann et al., 2006)
# - Non-gravitational accelerations acting upon the asteroid
# are included (Yarkovsky effect) a_nongrav = A2*t_vec*(au/r)^d, where t_vec is the
# unit heliocentric transverse vector, au is 1 astronomical unit, r is the
# asteroid's heliocentric range, A2 is a coefficient (with units of au/day^2),
# and d = 2.0
@taylorize function RNp1BP_pN_A_J23E_J2S_ng_eph!(dq, q, params, t)
    local jd0 = params[4] # Julian date of start time
    local dsj2k = t+(jd0-JD_J2000) # days since J2000.0 = 2.451545e6
    local ss16asteph_t = evaleph(params[1], dsj2k, q[1]) # params[2](t)*one(q[1]) #ss16asteph(t)
    local acceph_t = evaleph(params[2], dsj2k, q[1]) # params[2](t)*one(q[1]) #acc_eph(t)
    local newtonianNb_Potential_t = evaleph(params[3], dsj2k, q[1]) # params[3](t)*one(q[1]) #newtonianNb_Potential(t), massive bodies
    local S = eltype(q)
    local UJ_interaction = params[5] # interaction matrix with flattened bodies
    local N = params[6] # number of bodies, including NEA
    local Nm1 = N-1
    local μ = params[7] # vector of GM's

    # parameters related to speed of light, c
    local c_p2 = 29979.063823897606 # c^2 = 29979.063823897606 au^2/d^2
    local c_m2 = 3.3356611996764786e-5 # c^-2 = 3.3356611996764786e-5 d^2/au^2

    local zero_q_1 = auxzero(q[1]) # zero(q[1])

    X = Array{S}(undef, N)
    Y = Array{S}(undef, N)
    Z = Array{S}(undef, N)

    r_p2 = Array{S}(undef, N)
    r_p3d2 = Array{S}(undef, N)
    r_p7d2 = Array{S}(undef, N)

    newtonianCoeff = Array{S}(undef, N)

    ui = Array{S}(undef, N-1)
    vi = Array{S}(undef, N-1)
    wi = Array{S}(undef, N-1)

    #post-Newtonian stuff
    U = Array{S}(undef, N)
    V = Array{S}(undef, N)
    W = Array{S}(undef, N)

    _4U_m_3X = Array{S}(undef, N)
    _4V_m_3Y = Array{S}(undef, N)
    _4W_m_3Z = Array{S}(undef, N)

    UU = Array{S}(undef, N)
    VV = Array{S}(undef, N)
    WW = Array{S}(undef, N)

    r_p1d2 = Array{S}(undef, N)

    newtonianNb_Potential = Array{S}(undef, N)
    newtonian1b_Potential = Array{S}(undef, N)
    newtonianCoeff = Array{S}(undef, N)
    newton_acc_X = Array{S}(undef, N)
    newton_acc_Y = Array{S}(undef, N)
    newton_acc_Z = Array{S}(undef, N)

    v2 = Array{S}(undef, N)
    vi_dot_vj = Array{S}(undef, N)
    pn2 = Array{S}(undef, N)
    pn3 = Array{S}(undef, N)
    _4ϕj = Array{S}(undef, N)
    ϕi_plus_4ϕj = Array{S}(undef, N)
    sj2_plus_2si2_minus_4vivj = Array{S}(undef, N)
    ϕs_and_vs = Array{S}(undef, N)
    U_t_pn2 = Array{S}(undef, N)
    V_t_pn2 = Array{S}(undef, N)
    W_t_pn2 = Array{S}(undef, N)
    pn1t1_7 = Array{S}(undef, N)

    pntempX = zero_q_1
    pntempY = zero_q_1
    pntempZ = zero_q_1
    pn1 = Array{S}(undef, N)
    X_t_pn1 = Array{S}(undef, N)
    Y_t_pn1 = Array{S}(undef, N)
    Z_t_pn1 = Array{S}(undef, N)
    pNX_t_pn3 = Array{S}(undef, N)
    pNY_t_pn3 = Array{S}(undef, N)
    pNZ_t_pn3 = Array{S}(undef, N)
    pNX_t_X = Array{S}(undef, N)
    pNY_t_Y = Array{S}(undef, N)
    pNZ_t_Z = Array{S}(undef, N)

    # J2 acceleration auxiliaries
    t31 = Array{S}(undef, N)
    t32 = Array{S}(undef, N)
    t33 = Array{S}(undef, N)
    r_sin_ϕ = Array{S}(undef, N)
    F_J2_x = Array{S}(undef, N)
    F_J2_y = Array{S}(undef, N)
    F_J2_z = Array{S}(undef, N)
    F_J2_x1 = Array{S}(undef, N)
    F_J2_y1 = Array{S}(undef, N)
    F_J2_z1 = Array{S}(undef, N)
    F_J2_x2 = Array{S}(undef, N)
    F_J2_y2 = Array{S}(undef, N)
    F_J2_z2 = Array{S}(undef, N)
    # temp_accX_j = Array{S}(undef, N)
    # temp_accY_j = Array{S}(undef, N)
    # temp_accZ_j = Array{S}(undef, N)
    temp_accX_i = Array{S}(undef, N)
    temp_accY_i = Array{S}(undef, N)
    temp_accZ_i = Array{S}(undef, N)
    sin_ϕ = Array{S}(undef, N)
    sin2_ϕ = Array{S}(undef, N)
    sin3_ϕ = Array{S}(undef, N)
    sin4_ϕ = Array{S}(undef, N)
    ϕ = Array{S}(undef, N)
    cos_ϕ = Array{S}(undef, N)
    P_2_sin_ϕ = Array{S}(undef, N)
    ∂P_2_sin_ϕ = Array{S}(undef, N)
    P_3_sin_ϕ = Array{S}(undef, N)
    ∂P_3_sin_ϕ = Array{S}(undef, N)
    m_c_ϕ_∂P_2 = Array{S}(undef, N)
    m_c_ϕ_∂P_3 = Array{S}(undef, N)
    Λ2j_div_r4 = Array{S}(undef, N)
    Λ3j_div_r5 = Array{S}(undef, N)
    F_J_ξ = Array{S}(undef, N)
    F_J_η = Array{S}(undef, N)
    F_J_ζ = Array{S}(undef, N)
    F_J2_ξ = Array{S}(undef, N)
    F_J2_η = Array{S}(undef, N)
    F_J2_ζ = Array{S}(undef, N)
    F_J3_ξ = Array{S}(undef, N)
    F_J3_η = Array{S}(undef, N)
    F_J3_ζ = Array{S}(undef, N)
    ξx = Array{S}(undef, N)
    ξy = Array{S}(undef, N)
    ξz = Array{S}(undef, N)
    ηx = Array{S}(undef, N)
    ηy = Array{S}(undef, N)
    ηz = Array{S}(undef, N)
    ηx1 = Array{S}(undef, N)
    ηy1 = Array{S}(undef, N)
    ηz1 = Array{S}(undef, N)
    ηx2 = Array{S}(undef, N)
    ηy2 = Array{S}(undef, N)
    ηz2 = Array{S}(undef, N)
    ζx = Array{S}(undef, N)
    ζy = Array{S}(undef, N)
    ζz = Array{S}(undef, N)
    ζx1 = Array{S}(undef, N)
    ζy1 = Array{S}(undef, N)
    ζz1 = Array{S}(undef, N)
    ζx2 = Array{S}(undef, N)
    ζy2 = Array{S}(undef, N)
    ζz2 = Array{S}(undef, N)

    # extended-body accelerations
    accX = zero_q_1
    accY = zero_q_1
    accZ = zero_q_1

    # rotations to and from Earth, Sun and Moon pole-oriented frames
    local M_ = Array{S}(undef, 3, 3, N)
    local M_[:,:,ea] = t2c_jpl_de430(dsj2k)

    dq[1] = q[4]
    dq[2] = q[5]
    dq[3] = q[6]

    newtonianNb_Potential[N] = zero_q_1

    #compute point-mass Newtonian accelerations, all bodies
    for i in 1:Nm1
        ui[i] = ss16asteph_t[3(N-1+i)-2]
        vi[i] = ss16asteph_t[3(N-1+i)-1]
        wi[i] = ss16asteph_t[3(N-1+i)  ]

        X[i] = ss16asteph_t[3i-2]-q[1]
        Y[i] = ss16asteph_t[3i-1]-q[2]
        Z[i] = ss16asteph_t[3i  ]-q[3]

        U[i] = ui[i]-dq[1]
        V[i] = vi[i]-dq[2]
        W[i] = wi[i]-dq[3]

        _4U_m_3X[i] = (4dq[1]) - (3ui[i])
        _4V_m_3Y[i] = (4dq[2]) - (3vi[i])
        _4W_m_3Z[i] = (4dq[3]) - (3wi[i])

        pn2x = X[i]*_4U_m_3X[i]
        pn2y = Y[i]*_4V_m_3Y[i]
        pn2z = Z[i]*_4W_m_3Z[i]

        UU[i] = ui[i]*dq[1]
        VV[i] = vi[i]*dq[2]
        WW[i] = wi[i]*dq[3]

        vi_dot_vj[i] = ( UU[i]+VV[i] ) + WW[i]

        r_p2[i] = ( (X[i]^2)+(Y[i]^2) ) + (Z[i]^2)

        r_p1d2[i] = sqrt(r_p2[i])
        r_p3d2[i] = r_p2[i]^1.5
        r_p7d2[i] = r_p2[i]^3.5

        newtonianCoeff[i] =  μ[i]/r_p3d2[i]

        pn2[i] = newtonianCoeff[i]*(( pn2x+pn2y ) + pn2z)

        newton_acc_X[i] = X[i]*newtonianCoeff[i]
        newton_acc_Y[i] = Y[i]*newtonianCoeff[i]
        newton_acc_Z[i] = Z[i]*newtonianCoeff[i]

        newtonian1b_Potential[i] = μ[i]/r_p1d2[i]
        pn3[i] = 3.5newtonian1b_Potential[i]
        U_t_pn2[i] = pn2[i]*U[i]
        V_t_pn2[i] = pn2[i]*V[i]
        W_t_pn2[i] = pn2[i]*W[i]

        #J2 accelerations, if i-th body is flattened
        if UJ_interaction[i]
            # # rotate from inertial frame to extended-body frame
            t31[i] = -X[i]*M_[1,3,i]
            t32[i] = -Y[i]*M_[2,3,i]
            t33[i] = -Z[i]*M_[3,3,i]
            r_sin_ϕ[i] = (t31[i]+t32[i])+t33[i]

            # compute cartesian coordinates of acceleration due to body figure in body frame
            sin_ϕ[i] = r_sin_ϕ[i]/r_p1d2[i] # z/r = latitude of point mass wrt body frame
            ϕ[i] = asin(sin_ϕ[i])
            cos_ϕ[i] = cos(ϕ[i])
            sin2_ϕ[i] = sin_ϕ[i]^2
            sin3_ϕ[i] = sin_ϕ[i]^3
            P_2_sin_ϕ[i] = 1.5sin2_ϕ[i] - 0.5
            ∂P_2_sin_ϕ[i] = 3sin_ϕ[i]
            P_3_sin_ϕ[i] = (-1.5sin_ϕ[i]) + (2.5sin3_ϕ[i])
            ∂P_3_sin_ϕ[i] = -1.5 + 7.5sin2_ϕ[i]
            Λ2j_div_r4[i] = (-Λ2[i])/(r_p2[i]^2)
            Λ3j_div_r5[i] = (-Λ3[i])/(r_p1d2[i]^5)
            m_c_ϕ_∂P_2[i] = (-cos_ϕ[i])*∂P_2_sin_ϕ[i]
            m_c_ϕ_∂P_3[i] = (-cos_ϕ[i])*∂P_3_sin_ϕ[i]
            F_J2_ξ[i] = ( Λ2j_div_r4[i]*(3P_2_sin_ϕ[i]) )
            #F_J2_η[i] = zero_q_1
            F_J2_ζ[i] = Λ2j_div_r4[i]*m_c_ϕ_∂P_2[i]
            F_J3_ξ[i] = ( Λ3j_div_r5[i]*(4P_3_sin_ϕ[i]) )
            #F_J3_η[i] = zero_q_1
            F_J3_ζ[i] = Λ3j_div_r5[i]*m_c_ϕ_∂P_3[i]
            F_J_ξ[i] = F_J2_ξ[i] #+ F_J3_ξ[i]
            #F_J_η[i] = zero_q_1
            F_J_ζ[i] = F_J2_ζ[i] #+ F_J3_ζ[i]
            #Compute unit vectors ξ,η,ζ
            ξx[i] = -X[i]/r_p1d2[i]
            ξy[i] = -Y[i]/r_p1d2[i]
            ξz[i] = -Z[i]/r_p1d2[i]
            #Compute η = p x ξ
            ηx1[i] = M_[2,3,i]*ξz[i]
            ηy1[i] = M_[3,3,i]*ξx[i]
            ηz1[i] = M_[1,3,i]*ξy[i]
            ηx2[i] = M_[3,3,i]*ξy[i]
            ηy2[i] = M_[1,3,i]*ξz[i]
            ηz2[i] = M_[2,3,i]*ξx[i]
            ηx[i] = ηx1[i] - ηx2[i]
            ηy[i] = ηy1[i] - ηy2[i]
            ηz[i] = ηz1[i] - ηz2[i]
            #Compute ζ = ξ x η
            ζx1[i] = ξy[i]*ηz[i]
            ζy1[i] = ξz[i]*ηx[i]
            ζz1[i] = ξx[i]*ηy[i]
            ζx2[i] = ξz[i]*ηy[i]
            ζy2[i] = ξx[i]*ηz[i]
            ζz2[i] = ξy[i]*ηx[i]
            ζx[i] = ζx1[i] - ζx2[i]
            ζy[i] = ζy1[i] - ζy2[i]
            ζz[i] = ζz1[i] - ζz2[i]
            # compute cartesian coordinates of acceleration due to body figure in inertial frame
            F_J2_x1[i] = F_J_ξ[i]*ξx[i]
            F_J2_y1[i] = F_J_ξ[i]*ξy[i]
            F_J2_z1[i] = F_J_ξ[i]*ξz[i]
            F_J2_x2[i] = F_J_ζ[i]*ζx[i]
            F_J2_y2[i] = F_J_ζ[i]*ζy[i]
            F_J2_z2[i] = F_J_ζ[i]*ζz[i]
            F_J2_x[i] = F_J2_x1[i] + F_J2_x2[i]
            F_J2_y[i] = F_J2_y1[i] + F_J2_y2[i]
            F_J2_z[i] = F_J2_z1[i] + F_J2_z2[i]
        end # if UJ_interaction[i]
        v2[i] = ( (ui[i]^2)+(vi[i]^2) ) + (wi[i]^2)
    end #for, i
    v2[N] = ( (q[4]^2)+(q[5]^2) ) + (q[6]^2)

    for i in 1:Nm1
        temp_004 = newtonian1b_Potential[i] + newtonianNb_Potential[N]
        newtonianNb_Potential[N] = temp_004
        if UJ_interaction[i]
            # # # add result to total acceleration on upon j-th body figure due to i-th point mass
            # # @show "acc",j,"+μ",i,"Λ2",j
            # temp_accX_j[i,j] = accX[j] + (μ[i]*F_J2_x[i,j])
            # accX[j] = temp_accX_j[i,j]
            # temp_accY_j[i,j] = accY[j] + (μ[i]*F_J2_y[i,j])
            # accY[j] = temp_accY_j[i,j]
            # temp_accZ_j[i,j] = accZ[j] + (μ[i]*F_J2_z[i,j])
            # accZ[j] = temp_accZ_j[i,j]

            # # reaction force on i-th body
            # @show "acc",i,"-μ",j,"Λ2",j
            temp_accX_i[i] = accX - (μ[i]*F_J2_x[i])
            accX = temp_accX_i[i]
            temp_accY_i[i] = accY - (μ[i]*F_J2_y[i])
            accY = temp_accY_i[i]
            temp_accZ_i[i] = accZ - (μ[i]*F_J2_z[i])
            accZ = temp_accZ_i[i]
        end
    end

    #post-Newtonian accelerations due to Sun, Moon and planets (Mercury through Neptune)
    #Moyer, 1971, page 7 eq. 35
    # post-Newtonian iterative procedure setup and initialization
    _4ϕj[N] = 4newtonianNb_Potential[N]
    for i in 1:10
        ϕi_plus_4ϕj[i] = newtonianNb_Potential_t[i] + _4ϕj[N]
        sj2_plus_2si2_minus_4vivj[i] = ( (2v2[i]) - (4vi_dot_vj[i]) ) + v2[N]
        ϕs_and_vs[i] = sj2_plus_2si2_minus_4vivj[i] - ϕi_plus_4ϕj[i]
        Xij_t_Ui = X[i]*ui[i]
        Yij_t_Vi = Y[i]*vi[i]
        Zij_t_Wi = Z[i]*wi[i]
        Rij_dot_Vi = ( Xij_t_Ui+Yij_t_Vi ) + Zij_t_Wi
        # the expression below inside the (...)^2 should have a minus sign in front of the numerator,
        # but upon squaring it is eliminated, so at the end of the day, it is irrelevant ;)
        pn1t7 = (Rij_dot_Vi^2)/r_p2[i]
        pn1t2_7 = ϕs_and_vs[i] - (1.5pn1t7)
        pn1t1_7[i] = c_p2 + pn1t2_7

        pNX_t_X[i] = acceph_t[3i-2]*X[i]
        pNY_t_Y[i] = acceph_t[3i-1]*Y[i]
        pNZ_t_Z[i] = acceph_t[3i  ]*Z[i]
        pn1[i] = (  pn1t1_7[i]  +  (0.5*( (pNX_t_X[i]+pNY_t_Y[i]) + pNZ_t_Z[i] ))  )

        X_t_pn1[i] = newton_acc_X[i]*pn1[i]
        Y_t_pn1[i] = newton_acc_Y[i]*pn1[i]
        Z_t_pn1[i] = newton_acc_Z[i]*pn1[i]

        pNX_t_pn3[i] = acceph_t[3i-2]*pn3[i]
        pNY_t_pn3[i] = acceph_t[3i-1]*pn3[i]
        pNZ_t_pn3[i] = acceph_t[3i  ]*pn3[i]
    end #for i
    for i in 1:10
        termpnx = ( X_t_pn1[i] + (U_t_pn2[i]+pNX_t_pn3[i]) )
        sumpnx = pntempX + termpnx
        pntempX = sumpnx
        termpny = ( Y_t_pn1[i] + (V_t_pn2[i]+pNY_t_pn3[i]) )
        sumpny = pntempY + termpny
        pntempY = sumpny
        termpnz = ( Z_t_pn1[i] + (W_t_pn2[i]+pNZ_t_pn3[i]) )
        sumpnz = pntempZ + termpnz
        pntempZ = sumpnz
    end
    # compute Newtonian accelerations due to Pluto and 16 asteroid perturbers
    for i in 11:Nm1
        X_t_pn1[i] = c_p2*newton_acc_X[i]
        Y_t_pn1[i] = c_p2*newton_acc_Y[i]
        Z_t_pn1[i] = c_p2*newton_acc_Z[i]
    end #for i
    for i in 11:Nm1
        termpnx = X_t_pn1[i]
        sumpnx = pntempX + termpnx
        pntempX = sumpnx
        termpny = Y_t_pn1[i]
        sumpny = pntempY + termpny
        pntempY = sumpny
        termpnz = Z_t_pn1[i]
        sumpnz = pntempZ + termpnz
        pntempZ = sumpnz
    end
    postNewtonX = pntempX*c_m2
    postNewtonY = pntempY*c_m2
    postNewtonZ = pntempZ*c_m2

    # compute non-gravitational acceleration
    hx = (Y[1]*W[1])-(Z[1]*V[1])
    hy = (Z[1]*U[1])-(X[1]*W[1])
    hz = (X[1]*V[1])-(Y[1]*U[1])

    #cartesian components of transversal unit vector:
    t_x = (hz*Y[1]) - (hy*Z[1]) # Note: Y[1] = y_Sun - y_Apophis, etc.
    t_y = (hx*Z[1]) - (hz*X[1])
    t_z = (hy*X[1]) - (hx*Y[1])
    t_norm = sqrt( ((t_x^2)+(t_y^2))+(t_z^2) )
    t_x_unit = t_x/t_norm
    t_y_unit = t_y/t_norm
    t_z_unit = t_z/t_norm

    # cartesian components of radial unit vector:
    r_x_unit = -(X[1]/r_p1d2[1])
    r_y_unit = -(Y[1]/r_p1d2[1])
    r_z_unit = -(Z[1]/r_p1d2[1])

    # evaluate non-grav acceleration (solar radiation pressure, Yarkovsky):
    g_r = r_p2[1]
    A2_t_g_r = q[7]/g_r
    A1_t_g_r = q[8]/g_r

    NGAx = (A2_t_g_r*t_x_unit) + (A1_t_g_r*r_x_unit)
    NGAy = (A2_t_g_r*t_y_unit) + (A1_t_g_r*r_y_unit)
    NGAz = (A2_t_g_r*t_z_unit) + (A1_t_g_r*r_z_unit)

    dq[4] = ( postNewtonX + accX ) + NGAx
    dq[5] = ( postNewtonY + accY ) + NGAy
    dq[6] = ( postNewtonZ + accZ ) + NGAz

    dq[7] = zero_q_1

    nothing
end

@taylorize function RNp1BP_pN_A_J23E_J2S_ng_eph_threads!(dq, q, params, t)
    local jd0 = params[4] # Julian date of start time
    local dsj2k = t+(jd0-JD_J2000) # days since J2000.0 = 2.451545e6
    local ss16asteph_t = evaleph(params[1], dsj2k, q[1]) # params[2](t)*one(q[1]) #ss16asteph(t)
    local acceph_t = evaleph(params[2], dsj2k, q[1]) # params[2](t)*one(q[1]) #acc_eph(t)
    local newtonianNb_Potential_t = evaleph(params[3], dsj2k, q[1]) # params[3](t)*one(q[1]) #newtonianNb_Potential(t), massive bodies
    local S = eltype(q)
    local UJ_interaction = params[5] # interaction matrix with flattened bodies
    local N = params[6] # number of bodies, including NEA
    local Nm1 = N-1
    local μ = params[7] # vector of GM's

    # parameters related to speed of light, c
    local c_p2 = 29979.063823897606 # c^2 = 29979.063823897606 au^2/d^2
    local c_m2 = 3.3356611996764786e-5 # c^-2 = 3.3356611996764786e-5 d^2/au^2

    local zero_q_1 = auxzero(q[1]) # zero(q[1])

    X = Array{S}(undef, N)
    Y = Array{S}(undef, N)
    Z = Array{S}(undef, N)

    r_p2 = Array{S}(undef, N)
    r_p3d2 = Array{S}(undef, N)
    r_p7d2 = Array{S}(undef, N)

    newtonianCoeff = Array{S}(undef, N)

    ui = Array{S}(undef, N-1)
    vi = Array{S}(undef, N-1)
    wi = Array{S}(undef, N-1)

    #post-Newtonian stuff
    U = Array{S}(undef, N)
    V = Array{S}(undef, N)
    W = Array{S}(undef, N)

    _4U_m_3X = Array{S}(undef, N)
    _4V_m_3Y = Array{S}(undef, N)
    _4W_m_3Z = Array{S}(undef, N)

    UU = Array{S}(undef, N)
    VV = Array{S}(undef, N)
    WW = Array{S}(undef, N)

    r_p1d2 = Array{S}(undef, N)

    newtonianNb_Potential = Array{S}(undef, N)
    newtonian1b_Potential = Array{S}(undef, N)
    newtonianCoeff = Array{S}(undef, N)
    newton_acc_X = Array{S}(undef, N)
    newton_acc_Y = Array{S}(undef, N)
    newton_acc_Z = Array{S}(undef, N)

    v2 = Array{S}(undef, N)
    vi_dot_vj = Array{S}(undef, N)
    pn2 = Array{S}(undef, N)
    pn3 = Array{S}(undef, N)
    _4ϕj = Array{S}(undef, N)
    ϕi_plus_4ϕj = Array{S}(undef, N)
    sj2_plus_2si2_minus_4vivj = Array{S}(undef, N)
    ϕs_and_vs = Array{S}(undef, N)
    U_t_pn2 = Array{S}(undef, N)
    V_t_pn2 = Array{S}(undef, N)
    W_t_pn2 = Array{S}(undef, N)
    pn1t1_7 = Array{S}(undef, N)

    pntempX = zero_q_1
    pntempY = zero_q_1
    pntempZ = zero_q_1
    pn1 = Array{S}(undef, N)
    X_t_pn1 = Array{S}(undef, N)
    Y_t_pn1 = Array{S}(undef, N)
    Z_t_pn1 = Array{S}(undef, N)
    pNX_t_pn3 = Array{S}(undef, N)
    pNY_t_pn3 = Array{S}(undef, N)
    pNZ_t_pn3 = Array{S}(undef, N)
    pNX_t_X = Array{S}(undef, N)
    pNY_t_Y = Array{S}(undef, N)
    pNZ_t_Z = Array{S}(undef, N)

    # J2 acceleration auxiliaries
    t31 = Array{S}(undef, N)
    t32 = Array{S}(undef, N)
    t33 = Array{S}(undef, N)
    r_sin_ϕ = Array{S}(undef, N)
    F_J2_x = Array{S}(undef, N)
    F_J2_y = Array{S}(undef, N)
    F_J2_z = Array{S}(undef, N)
    F_J2_x1 = Array{S}(undef, N)
    F_J2_y1 = Array{S}(undef, N)
    F_J2_z1 = Array{S}(undef, N)
    F_J2_x2 = Array{S}(undef, N)
    F_J2_y2 = Array{S}(undef, N)
    F_J2_z2 = Array{S}(undef, N)
    # temp_accX_j = Array{S}(undef, N)
    # temp_accY_j = Array{S}(undef, N)
    # temp_accZ_j = Array{S}(undef, N)
    temp_accX_i = Array{S}(undef, N)
    temp_accY_i = Array{S}(undef, N)
    temp_accZ_i = Array{S}(undef, N)
    sin_ϕ = Array{S}(undef, N)
    sin2_ϕ = Array{S}(undef, N)
    sin3_ϕ = Array{S}(undef, N)
    sin4_ϕ = Array{S}(undef, N)
    ϕ = Array{S}(undef, N)
    cos_ϕ = Array{S}(undef, N)
    P_2_sin_ϕ = Array{S}(undef, N)
    ∂P_2_sin_ϕ = Array{S}(undef, N)
    P_3_sin_ϕ = Array{S}(undef, N)
    ∂P_3_sin_ϕ = Array{S}(undef, N)
    m_c_ϕ_∂P_2 = Array{S}(undef, N)
    m_c_ϕ_∂P_3 = Array{S}(undef, N)
    Λ2j_div_r4 = Array{S}(undef, N)
    Λ3j_div_r5 = Array{S}(undef, N)
    F_J_ξ = Array{S}(undef, N)
    F_J_η = Array{S}(undef, N)
    F_J_ζ = Array{S}(undef, N)
    F_J2_ξ = Array{S}(undef, N)
    F_J2_η = Array{S}(undef, N)
    F_J2_ζ = Array{S}(undef, N)
    F_J3_ξ = Array{S}(undef, N)
    F_J3_η = Array{S}(undef, N)
    F_J3_ζ = Array{S}(undef, N)
    ξx = Array{S}(undef, N)
    ξy = Array{S}(undef, N)
    ξz = Array{S}(undef, N)
    ηx = Array{S}(undef, N)
    ηy = Array{S}(undef, N)
    ηz = Array{S}(undef, N)
    ηx1 = Array{S}(undef, N)
    ηy1 = Array{S}(undef, N)
    ηz1 = Array{S}(undef, N)
    ηx2 = Array{S}(undef, N)
    ηy2 = Array{S}(undef, N)
    ηz2 = Array{S}(undef, N)
    ζx = Array{S}(undef, N)
    ζy = Array{S}(undef, N)
    ζz = Array{S}(undef, N)
    ζx1 = Array{S}(undef, N)
    ζy1 = Array{S}(undef, N)
    ζz1 = Array{S}(undef, N)
    ζx2 = Array{S}(undef, N)
    ζy2 = Array{S}(undef, N)
    ζz2 = Array{S}(undef, N)

    # extended-body accelerations
    accX = zero_q_1
    accY = zero_q_1
    accZ = zero_q_1

    # rotations to and from Earth, Sun and Moon pole-oriented frames
    local M_ = Array{S}(undef, 3, 3, N)
    local M_[:,:,ea] = t2c_jpl_de430(dsj2k)

    dq[1] = q[4]
    dq[2] = q[5]
    dq[3] = q[6]

    newtonianNb_Potential[N] = zero_q_1

    #compute point-mass Newtonian accelerations, all bodies
    Threads.@threads for i in 1:Nm1
        ui[i] = ss16asteph_t[3(N-1+i)-2]
        vi[i] = ss16asteph_t[3(N-1+i)-1]
        wi[i] = ss16asteph_t[3(N-1+i)  ]

        X[i] = ss16asteph_t[3i-2]-q[1]
        Y[i] = ss16asteph_t[3i-1]-q[2]
        Z[i] = ss16asteph_t[3i  ]-q[3]

        U[i] = ui[i]-dq[1]
        V[i] = vi[i]-dq[2]
        W[i] = wi[i]-dq[3]

        _4U_m_3X[i] = (4dq[1]) - (3ui[i])
        _4V_m_3Y[i] = (4dq[2]) - (3vi[i])
        _4W_m_3Z[i] = (4dq[3]) - (3wi[i])

        pn2x = X[i]*_4U_m_3X[i]
        pn2y = Y[i]*_4V_m_3Y[i]
        pn2z = Z[i]*_4W_m_3Z[i]

        UU[i] = ui[i]*dq[1]
        VV[i] = vi[i]*dq[2]
        WW[i] = wi[i]*dq[3]

        vi_dot_vj[i] = ( UU[i]+VV[i] ) + WW[i]

        r_p2[i] = ( (X[i]^2)+(Y[i]^2) ) + (Z[i]^2)

        r_p1d2[i] = sqrt(r_p2[i])
        r_p3d2[i] = r_p2[i]^1.5
        r_p7d2[i] = r_p2[i]^3.5

        newtonianCoeff[i] =  μ[i]/r_p3d2[i]

        pn2[i] = newtonianCoeff[i]*(( pn2x+pn2y ) + pn2z)

        newton_acc_X[i] = X[i]*newtonianCoeff[i]
        newton_acc_Y[i] = Y[i]*newtonianCoeff[i]
        newton_acc_Z[i] = Z[i]*newtonianCoeff[i]

        newtonian1b_Potential[i] = μ[i]/r_p1d2[i]
        pn3[i] = 3.5newtonian1b_Potential[i]
        U_t_pn2[i] = pn2[i]*U[i]
        V_t_pn2[i] = pn2[i]*V[i]
        W_t_pn2[i] = pn2[i]*W[i]

        #J2 accelerations, if i-th body is flattened
        if UJ_interaction[i]
            # # rotate from inertial frame to extended-body frame
            t31[i] = -X[i]*M_[1,3,i]
            t32[i] = -Y[i]*M_[2,3,i]
            t33[i] = -Z[i]*M_[3,3,i]
            r_sin_ϕ[i] = (t31[i]+t32[i])+t33[i]

            # compute cartesian coordinates of acceleration due to body figure in body frame
            sin_ϕ[i] = r_sin_ϕ[i]/r_p1d2[i] # z/r = latitude of point mass wrt body frame
            ϕ[i] = asin(sin_ϕ[i])
            cos_ϕ[i] = cos(ϕ[i])
            sin2_ϕ[i] = sin_ϕ[i]^2
            sin3_ϕ[i] = sin_ϕ[i]^3
            P_2_sin_ϕ[i] = 1.5sin2_ϕ[i] - 0.5
            ∂P_2_sin_ϕ[i] = 3sin_ϕ[i]
            P_3_sin_ϕ[i] = (-1.5sin_ϕ[i]) + (2.5sin3_ϕ[i])
            ∂P_3_sin_ϕ[i] = -1.5 + 7.5sin2_ϕ[i]
            Λ2j_div_r4[i] = (-Λ2[i])/(r_p2[i]^2)
            Λ3j_div_r5[i] = (-Λ3[i])/(r_p1d2[i]^5)
            m_c_ϕ_∂P_2[i] = (-cos_ϕ[i])*∂P_2_sin_ϕ[i]
            m_c_ϕ_∂P_3[i] = (-cos_ϕ[i])*∂P_3_sin_ϕ[i]
            F_J2_ξ[i] = ( Λ2j_div_r4[i]*(3P_2_sin_ϕ[i]) )
            #F_J2_η[i] = zero_q_1
            F_J2_ζ[i] = Λ2j_div_r4[i]*m_c_ϕ_∂P_2[i]
            F_J3_ξ[i] = ( Λ3j_div_r5[i]*(4P_3_sin_ϕ[i]) )
            #F_J3_η[i] = zero_q_1
            F_J3_ζ[i] = Λ3j_div_r5[i]*m_c_ϕ_∂P_3[i]
            F_J_ξ[i] = F_J2_ξ[i] #+ F_J3_ξ[i]
            #F_J_η[i] = zero_q_1
            F_J_ζ[i] = F_J2_ζ[i] #+ F_J3_ζ[i]
            #Compute unit vectors ξ,η,ζ
            ξx[i] = -X[i]/r_p1d2[i]
            ξy[i] = -Y[i]/r_p1d2[i]
            ξz[i] = -Z[i]/r_p1d2[i]
            #Compute η = p x ξ
            ηx1[i] = M_[2,3,i]*ξz[i]
            ηy1[i] = M_[3,3,i]*ξx[i]
            ηz1[i] = M_[1,3,i]*ξy[i]
            ηx2[i] = M_[3,3,i]*ξy[i]
            ηy2[i] = M_[1,3,i]*ξz[i]
            ηz2[i] = M_[2,3,i]*ξx[i]
            ηx[i] = ηx1[i] - ηx2[i]
            ηy[i] = ηy1[i] - ηy2[i]
            ηz[i] = ηz1[i] - ηz2[i]
            #Compute ζ = ξ x η
            ζx1[i] = ξy[i]*ηz[i]
            ζy1[i] = ξz[i]*ηx[i]
            ζz1[i] = ξx[i]*ηy[i]
            ζx2[i] = ξz[i]*ηy[i]
            ζy2[i] = ξx[i]*ηz[i]
            ζz2[i] = ξy[i]*ηx[i]
            ζx[i] = ζx1[i] - ζx2[i]
            ζy[i] = ζy1[i] - ζy2[i]
            ζz[i] = ζz1[i] - ζz2[i]
            # compute cartesian coordinates of acceleration due to body figure in inertial frame
            F_J2_x1[i] = F_J_ξ[i]*ξx[i]
            F_J2_y1[i] = F_J_ξ[i]*ξy[i]
            F_J2_z1[i] = F_J_ξ[i]*ξz[i]
            F_J2_x2[i] = F_J_ζ[i]*ζx[i]
            F_J2_y2[i] = F_J_ζ[i]*ζy[i]
            F_J2_z2[i] = F_J_ζ[i]*ζz[i]
            F_J2_x[i] = F_J2_x1[i] + F_J2_x2[i]
            F_J2_y[i] = F_J2_y1[i] + F_J2_y2[i]
            F_J2_z[i] = F_J2_z1[i] + F_J2_z2[i]
        end # if UJ_interaction[i]
        v2[i] = ( (ui[i]^2)+(vi[i]^2) ) + (wi[i]^2)
    end #for, i
    v2[N] = ( (q[4]^2)+(q[5]^2) ) + (q[6]^2)

    for i in 1:Nm1
        temp_004 = newtonian1b_Potential[i] + newtonianNb_Potential[N]
        newtonianNb_Potential[N] = temp_004
        if UJ_interaction[i]
            # # # add result to total acceleration on upon j-th body figure due to i-th point mass
            # # @show "acc",j,"+μ",i,"Λ2",j
            # temp_accX_j[i,j] = accX[j] + (μ[i]*F_J2_x[i,j])
            # accX[j] = temp_accX_j[i,j]
            # temp_accY_j[i,j] = accY[j] + (μ[i]*F_J2_y[i,j])
            # accY[j] = temp_accY_j[i,j]
            # temp_accZ_j[i,j] = accZ[j] + (μ[i]*F_J2_z[i,j])
            # accZ[j] = temp_accZ_j[i,j]

            # # reaction force on i-th body
            # @show "acc",i,"-μ",j,"Λ2",j
            temp_accX_i[i] = accX - (μ[i]*F_J2_x[i])
            accX = temp_accX_i[i]
            temp_accY_i[i] = accY - (μ[i]*F_J2_y[i])
            accY = temp_accY_i[i]
            temp_accZ_i[i] = accZ - (μ[i]*F_J2_z[i])
            accZ = temp_accZ_i[i]
        end
    end

    #post-Newtonian accelerations due to Sun, Moon and planets (Mercury through Neptune)
    #Moyer, 1971, page 7 eq. 35
    # post-Newtonian iterative procedure setup and initialization
    _4ϕj[N] = 4newtonianNb_Potential[N]
    Threads.@threads for i in 1:10
        ϕi_plus_4ϕj[i] = newtonianNb_Potential_t[i] + _4ϕj[N]
        sj2_plus_2si2_minus_4vivj[i] = ( (2v2[i]) - (4vi_dot_vj[i]) ) + v2[N]
        ϕs_and_vs[i] = sj2_plus_2si2_minus_4vivj[i] - ϕi_plus_4ϕj[i]
        Xij_t_Ui = X[i]*ui[i]
        Yij_t_Vi = Y[i]*vi[i]
        Zij_t_Wi = Z[i]*wi[i]
        Rij_dot_Vi = ( Xij_t_Ui+Yij_t_Vi ) + Zij_t_Wi
        # the expression below inside the (...)^2 should have a minus sign in front of the numerator,
        # but upon squaring it is eliminated, so at the end of the day, it is irrelevant ;)
        pn1t7 = (Rij_dot_Vi^2)/r_p2[i]
        pn1t2_7 = ϕs_and_vs[i] - (1.5pn1t7)
        pn1t1_7[i] = c_p2 + pn1t2_7

        pNX_t_X[i] = acceph_t[3i-2]*X[i]
        pNY_t_Y[i] = acceph_t[3i-1]*Y[i]
        pNZ_t_Z[i] = acceph_t[3i  ]*Z[i]
        pn1[i] = (  pn1t1_7[i]  +  (0.5*( (pNX_t_X[i]+pNY_t_Y[i]) + pNZ_t_Z[i] ))  )

        X_t_pn1[i] = newton_acc_X[i]*pn1[i]
        Y_t_pn1[i] = newton_acc_Y[i]*pn1[i]
        Z_t_pn1[i] = newton_acc_Z[i]*pn1[i]

        pNX_t_pn3[i] = acceph_t[3i-2]*pn3[i]
        pNY_t_pn3[i] = acceph_t[3i-1]*pn3[i]
        pNZ_t_pn3[i] = acceph_t[3i  ]*pn3[i]
    end #for i
    for i in 1:10
        termpnx = ( X_t_pn1[i] + (U_t_pn2[i]+pNX_t_pn3[i]) )
        sumpnx = pntempX + termpnx
        pntempX = sumpnx
        termpny = ( Y_t_pn1[i] + (V_t_pn2[i]+pNY_t_pn3[i]) )
        sumpny = pntempY + termpny
        pntempY = sumpny
        termpnz = ( Z_t_pn1[i] + (W_t_pn2[i]+pNZ_t_pn3[i]) )
        sumpnz = pntempZ + termpnz
        pntempZ = sumpnz
    end
    # compute Newtonian accelerations due to Pluto and 16 asteroid perturbers
    Threads.@threads for i in 11:Nm1
        X_t_pn1[i] = c_p2*newton_acc_X[i]
        Y_t_pn1[i] = c_p2*newton_acc_Y[i]
        Z_t_pn1[i] = c_p2*newton_acc_Z[i]
    end #for i
    for i in 11:Nm1
        termpnx = X_t_pn1[i]
        sumpnx = pntempX + termpnx
        pntempX = sumpnx
        termpny = Y_t_pn1[i]
        sumpny = pntempY + termpny
        pntempY = sumpny
        termpnz = Z_t_pn1[i]
        sumpnz = pntempZ + termpnz
        pntempZ = sumpnz
    end
    postNewtonX = pntempX*c_m2
    postNewtonY = pntempY*c_m2
    postNewtonZ = pntempZ*c_m2

    # compute non-gravitational acceleration
    hx = (Y[1]*W[1])-(Z[1]*V[1])
    hy = (Z[1]*U[1])-(X[1]*W[1])
    hz = (X[1]*V[1])-(Y[1]*U[1])

    #cartesian components of transversal unit vector: t = r × h
    t_x = (hz*Y[1]) - (hy*Z[1]) # Note: Y[1] = y_Sun - y_Apophis, etc.
    t_y = (hx*Z[1]) - (hz*X[1])
    t_z = (hy*X[1]) - (hx*Y[1])
    t_norm = sqrt( ((t_x^2)+(t_y^2))+(t_z^2) )
    t_x_unit = t_x/t_norm
    t_y_unit = t_y/t_norm
    t_z_unit = t_z/t_norm

    # cartesian components of radial unit vector:
    r_x_unit = -(X[1]/r_p1d2[1])
    r_y_unit = -(Y[1]/r_p1d2[1])
    r_z_unit = -(Z[1]/r_p1d2[1])

    # evaluate non-grav acceleration (solar radiation pressure, Yarkovsky):
    g_r = r_p2[1]
    A2_t_g_r = q[7]/g_r
    A1_t_g_r = q[8]/g_r

    NGAx = (A2_t_g_r*t_x_unit) + (A1_t_g_r*r_x_unit)
    NGAy = (A2_t_g_r*t_y_unit) + (A1_t_g_r*r_y_unit)
    NGAz = (A2_t_g_r*t_z_unit) + (A1_t_g_r*r_z_unit)

    dq[4] = ( postNewtonX + accX ) + NGAx
    dq[5] = ( postNewtonY + accY ) + NGAy
    dq[6] = ( postNewtonZ + accZ ) + NGAz

    dq[7] = zero_q_1

    nothing
end


text:
@testset "URDF parse" begin
    @testset "joint bounds" begin
        acrobot = parse_urdf(joinpath(@__DIR__, "urdf", "Acrobot.urdf"), remove_fixed_tree_joints=false)
        @test position_bounds(findjoint(acrobot, "shoulder")) == [RigidBodyDynamics.Bounds(-Inf, Inf)]
        @test velocity_bounds(findjoint(acrobot, "shoulder")) == [RigidBodyDynamics.Bounds(-Inf, Inf)]
        @test effort_bounds(findjoint(acrobot, "shoulder")) == [RigidBodyDynamics.Bounds(-Inf, Inf)]
        @test position_bounds(findjoint(acrobot, "elbow")) == [RigidBodyDynamics.Bounds(-Inf, Inf)]
        @test velocity_bounds(findjoint(acrobot, "elbow")) == [RigidBodyDynamics.Bounds(-Inf, Inf)]
        @test effort_bounds(findjoint(acrobot, "elbow")) == [RigidBodyDynamics.Bounds(-Inf, Inf)]

        acrobot_with_limits = parse_urdf(joinpath(@__DIR__, "urdf", "Acrobot_with_limits.urdf"), remove_fixed_tree_joints=false)
        @test position_bounds(findjoint(acrobot_with_limits, "shoulder")) == [RigidBodyDynamics.Bounds(-6.28, 6.28)]
        @test velocity_bounds(findjoint(acrobot_with_limits, "shoulder")) == [RigidBodyDynamics.Bounds(-10, 10)]
        @test effort_bounds(findjoint(acrobot_with_limits, "shoulder")) == [RigidBodyDynamics.Bounds(0, 0)]
        @test position_bounds(findjoint(acrobot_with_limits, "elbow")) == [RigidBodyDynamics.Bounds(-6.28, 6.28)]
        @test velocity_bounds(findjoint(acrobot_with_limits, "elbow")) == [RigidBodyDynamics.Bounds(-10, 10)]
        @test effort_bounds(findjoint(acrobot_with_limits, "elbow")) == [RigidBodyDynamics.Bounds(-5, 5)]
    end

    @testset "planar joints" begin
        robot = parse_urdf(joinpath(@__DIR__, "urdf", "planar_slider.urdf"), remove_fixed_tree_joints=false)
        # the first joint is the fixed joint between the world and the base link
        jt = joint_type(joints(robot)[1])
        @test jt isa Fixed

        # the 2nd joint is planar with <axis xyz="1 0 0"/>
        jt = joint_type(joints(robot)[2])
        @test jt isa Planar
        @test jt.x_axis ≈ SVector(0, 0, -1)
        @test jt.y_axis ≈ SVector(0, 1, 0)

        # the 3rd joint is planar with <axis xyz="0 1 0"/>
        jt = joint_type(joints(robot)[3])
        @test jt isa Planar
        @test jt.x_axis ≈ SVector(1, 0, 0)
        @test jt.y_axis ≈ SVector(0, 0, -1)

        # the 4th joint is planar with <axis xyz="0 0 1"/>
        jt = joint_type(joints(robot)[4])
        @test jt isa Planar
        @test jt.x_axis ≈ SVector(1, 0, 0)
        @test jt.y_axis ≈ SVector(0, 1, 0)
    end

    @testset "rotation handling" begin
        # https://github.com/JuliaRobotics/RigidBodyDynamics.jl/issues/429
        xml = """
            <origin rpy="0 0 0"/>
        """
        xpose = LightXML.root(LightXML.parse_string(xml))
        rot, trans = RigidBodyDynamics.parse_pose(Float64, xpose)
        @test rot ≈ SDiagonal(1, 1, 1)
        @test trans ≈ SVector(0, 0, 0)

        xml = """
            <origin rpy="0.3 0 0"/>
        """
        xpose = LightXML.root(LightXML.parse_string(xml))
        rot, trans = RigidBodyDynamics.parse_pose(Float64, xpose)
        @test rot ≈ RotMatrix(RotX(0.3))
        @test trans ≈ SVector(0, 0, 0)

        xml = """
            <origin rpy="0 0.2 0"/>
        """
        xpose = LightXML.root(LightXML.parse_string(xml))
        rot, trans = RigidBodyDynamics.parse_pose(Float64, xpose)
        @test rot ≈ RotMatrix(RotY(0.2))
        @test trans ≈ SVector(0, 0, 0)

        xml = """
            <origin rpy="0 0 0.1"/>
        """
        xpose = LightXML.root(LightXML.parse_string(xml))
        rot, trans = RigidBodyDynamics.parse_pose(Float64, xpose)
        @test rot ≈ RotMatrix(RotZ(0.1))
        @test trans ≈ SVector(0, 0, 0)

        # Comparison against ROS's tf.transformations.quaternion_from_euler
        xml = """
            <origin rpy="1 2 3"/>
        """
        xpose = LightXML.root(LightXML.parse_string(xml))
        rot, trans = RigidBodyDynamics.parse_pose(Float64, xpose)
        @test rot ≈ [0.41198225 -0.83373765 -0.36763046;
                    -0.05872664 -0.42691762  0.90238159;
                    -0.90929743 -0.35017549 -0.2248451] atol=1e-7
        @test rot ≈ RotZ(3) * RotY(2) * RotX(1)
        @test trans ≈ SVector(0, 0, 0)

        # Comparison against ROS's tf.transformations.quaternion_from_euler
        xml = """
            <origin rpy="0.5 0.1 0.2"/>
        """
        xpose = LightXML.root(LightXML.parse_string(xml))
        rot, trans = RigidBodyDynamics.parse_pose(Float64, xpose)
        @test rot ≈ [0.97517033 -0.12744012  0.18111281;
                     0.19767681  0.86959819 -0.45246312;
                    -0.09983342  0.47703041  0.8731983] atol=1e-7
        @test rot ≈ RotZ(0.2) * RotY(0.1) * RotX(0.5)
        @test trans ≈ SVector(0, 0, 0)
    end
end

function test_inverse_dynamics_match(mechanism1::Mechanism, mechanism2::Mechanism; atol::Real=0, rtol::Real=atol>0 ? 0 : √eps)
    state1 = MechanismState(mechanism1)
    v̇_vec = rand(num_velocities(mechanism1))
    rand!(state1)
    v̇1 = copyto!(similar(velocity(state1)), v̇_vec)
    τ1 = inverse_dynamics(state1, v̇1)

    state2 = MechanismState(mechanism2)
    copyto!(state2, Vector(state1))
    v̇2 = copyto!(similar(velocity(state2)), v̇_vec)
    τ2 = inverse_dynamics(state2, v̇2)

    @test τ1 ≈ τ2 rtol=rtol atol=atol
end

function test_kinematic_graph_layout_match(mechanism1::Mechanism, mechanism2::Mechanism)
    for (j1, j2) in zip(joints(mechanism1), joints(mechanism2))
        @test string(j1) == string(j2)
        @test string(predecessor(j1, mechanism1)) == string(predecessor(j2, mechanism2))
        @test string(successor(j1, mechanism1)) == string(successor(j2, mechanism2))
        @test joint_type(j1) == joint_type(j2)
    end
end

function test_urdf_serialize_deserialize(mechanism1::Mechanism; remove_fixed_tree_joints::Bool)
    state1 = MechanismState(mechanism1)

    mktempdir() do dir
        urdf2 = joinpath(dir, "test.urdf")
        write_urdf(urdf2, mechanism1, robot_name="test")
        mechanism2 = parse_urdf(urdf2, remove_fixed_tree_joints=remove_fixed_tree_joints)
        test_inverse_dynamics_match(mechanism1, mechanism2, atol=1e-10)

        for joint1 in tree_joints(mechanism1)
            if remove_fixed_tree_joints && joint_type(joint1) isa Fixed
                continue
            else
                joint2 = findjoint(mechanism2, string(joint1))
                @test position_bounds(joint1) == position_bounds(joint2)
                @test velocity_bounds(joint1) == velocity_bounds(joint2)
                @test effort_bounds(joint1) == effort_bounds(joint2)
            end
        end
    end
end

@testset "URDF write" begin
    @testset "Basics" begin
        Random.seed!(124)
        urdfdir = joinpath(@__DIR__, "urdf")
        for basename in readdir(urdfdir)
            last(splitext(basename)) == ".urdf" || continue
            urdf = joinpath(urdfdir, basename)
            for floating in (true, false)
                for remove_fixed_joints_before in (true, false)
                    mechanism = parse_urdf(urdf, remove_fixed_tree_joints=remove_fixed_joints_before, floating=floating)
                    for remove_fixed_joints_after in (true, false)
                        test_urdf_serialize_deserialize(mechanism, remove_fixed_tree_joints=remove_fixed_joints_after)
                    end
                end
            end
        end
    end

    @testset "include_root" begin
        Random.seed!(1245)
        urdfdir = joinpath(@__DIR__, "urdf")
        for basename in readdir(urdfdir)
            last(splitext(basename)) == ".urdf" || continue
            urdf = joinpath(urdfdir, basename)

            unmodified_mechanism = parse_urdf(urdf, remove_fixed_tree_joints=false)
            floating_mechanism = parse_urdf(urdf, remove_fixed_tree_joints=false, floating=true)

            # Write unmodified mechanism to URDF using `include_root=false`, then parse
            # using the default `Fixed` joint type as the root joint type. Ensure we get
            # the same kinematic tree.
            mktempdir() do dir
                urdf2 = joinpath(dir, "test.urdf")
                write_urdf(urdf2, unmodified_mechanism, include_root=false)
                unmodified_mechanism_back = parse_urdf(urdf2, remove_fixed_tree_joints=false)
                test_kinematic_graph_layout_match(unmodified_mechanism, unmodified_mechanism_back)
                test_inverse_dynamics_match(unmodified_mechanism, unmodified_mechanism_back, atol=1e-10)
            end

            # Write floating-base mechanism to URDF using `include_root=true` (default), then parse
            # using the default `Fixed` joint type as the root joint type:
            mktempdir() do dir
                urdf2 = joinpath(dir, "test.urdf")
                write_urdf(urdf2, floating_mechanism)
                floating_mechanism_back = parse_urdf(urdf2)
                # Note: floating_mechanism_back has an extra fixed joint at the root, compared to floating_mechanism,
                # but is dynamically equivalent.
                test_inverse_dynamics_match(floating_mechanism, floating_mechanism_back, atol=1e-10)
            end

            # Write floating-base mechanism to URDF using `include_root=false`, then parse and compare
            # to unmodified mechanism:
            mktempdir() do dir
                urdf2 = joinpath(dir, "test.urdf")
                write_urdf(urdf2, floating_mechanism, include_root=false)
                unmodified_mechanism_back = parse_urdf(urdf2, remove_fixed_tree_joints=false)
                test_kinematic_graph_layout_match(unmodified_mechanism, unmodified_mechanism_back)
                test_inverse_dynamics_match(unmodified_mechanism, unmodified_mechanism_back, atol=1e-10)
            end
        end
    end
end

ext: jl


text:
export permutations_in_Sn, compose, cycle_structure, data, part, #Functions working with partitions and permutations
    partition, Haar, expectation, WeingartenUnitary

typealias partition Vector{Int}
#Functions working with partitions and permutations
# TODO Migrate these functions to Catalan

#Iterate over partitions of n in lexicographic order
function part(n::Integer)
    if n==0 produce([]) end
    if n<=0 return end
    for p in @task part(n-1)
        p = [p; 1]
        produce(p)
        p = p[1:end-1]
        if length(p) == 1 || (length(p)>1 && p[end]<p[end-1])
            p[end] += 1
            produce(p)
            p[end] -= 1
        end
    end
end

function permutations_in_Sn(n::Integer)
    P = permutation_calloc(n)
    while true
        produce(P)
        try permutation_next(P) catch break end
    end
end

function compose(P::Ptr{gsl_permutation}, Q::Ptr{gsl_permutation})
    #Compose the permutations
    n=convert(Int64, permutation_size(P))
    @assert n==permutation_size(Q)
    x=permutation_alloc(n)
    Pv = [x+1 for x in pointer_to_array(permutation_data(P), (n,))]
    Qv = [x+1 for x in pointer_to_array(permutation_data(Q), (n,))]
    Xv = [Qv[i] for i in Pv]
    for i=1:n
         unsafe_assign(permutation_data(x), Xv[i]-1, i)
    end
    permutation_valid(x)
    x
end

#Compute cycle structure, i.e. lengths of each cycle in the cycle factorization of the permutatation
function cycle_structure(P::Ptr{gsl_permutation})
    n=convert(Int64, permutation_size(P))
    Pcanon = permutation_linear_to_canonical(P)
    PCANON = data(Pcanon)
    HaveNewCycleAtPos = [PCANON[i]>PCANON[i+1] for i=1:n-1]
    cycleindices = Int[1]
    for i=1:n-1 if HaveNewCycleAtPos[i] push!(cycleindices, i+1) end end
    push!(cycleindices, n+1)
    cyclestructure = Int[cycleindices[i+1]-cycleindices[i] for i=1:length(cycleindices)-1]
    @assert sum(cyclestructure) == n
    cyclestructure
end

#Returns a vector of indices (starting from 1 in the Julia convention)
data(P::Ptr{gsl_permutation}) = [convert(Int64, x)+1 for x in
    pointer_to_array(permutation_data(P), (convert(Int64, permutation_size(P)) ,))]

type Haar <: ContinuousMatrixDistribution
    beta::Real
end

# In random matrix theory one often encounters expressions of the form
#
#X = Q * A * Q' * B
#
#where A and B are deterministic matrices with fixed numerical matrix entries
#and Q is a random matrix that does not have explicitly defined matrix
#elements. Instead, one takes an expectation over of expressions of this form
#and this "integrates out" the Qs to produce a numeric result.
#
#expectation(X) #= some number
#
#Here is a function that symbolically calculates the expectation of a product
#of matrices over the symmetric group that Q is uniform Haar over.
#It takes an expression consisting of a product of matrices and replaces it
#with an evaluated symbolic expression which is the expectation.
function expectation(X::Expr)
    if X.head != :call
        error(string("Unexpected type of expression: ", X.head))
    end

    n = length(X.args) - 1
    if n < 1 return eval(X) end #nothing to do Haar-wise

    if X.args[1] != :*
        error("Unexpected expression, only products supported")
    end

    # Parse expression involving products of matrices to extract the
    # positions of Haar matrices and their ctransposes
    Qidx=[] #Indices for Haar matrices
    Qpidx=[] #Indices for ctranspose of Haar matrices
    Others=[]
    MyQ=None
    for i=1:n
        thingy=X.args[i+1]
        if isa(thingy, Symbol)
            if isa(eval(thingy), Haar)
                if MyQ==None MyQ=thingy end
                if MyQ == thingy
                    Qidx=[Qidx; i]
                else
                    warning("only one instance of Haar supported, skipping the other guy ", thingy) end
            else
                Others = [Others; (thingy, i, i+1)]
            end
            println(i, ' ', thingy, "::", typeof(eval(thingy)))
        elseif isa(thingy, Expr)
            println(i, ' ', thingy, "::Expr")
            if thingy.head==symbol('\'') && length(thingy.args)>=1 #Maybe this is a Q'
                if isa(thingy.args[1], Symbol) && isa(eval(thingy.args[1]), Haar)
                    println("Here is a Qtranspose")
                    Qpidx=[Qpidx; i]
                end
            end
        else
            error(string("Unexpected token ", thingy ," of type ", typeof(thingy)))
        end
    end
    if length(Qidx) == length(Qpidx) == 0 return eval(X) end #nothing to do Haar-wise
    println(MyQ, " is in places ", Qidx)
    println(MyQ, "' is in places ", Qpidx)

    n = length(Qidx)
    #If there are different Qs and Q's, the answer is a big fat 0
    if n != length(Qpidx) return zeros(size(eval(X.args[2]),1)...) end

    ##################################
    # Step 2. Enumerate permutations #
    ##################################
    AllTerms = Any[]
    for sigma in @task permutations_in_Sn(n)
        for tau in @task permutations_in_Sn(n)
            sigma_inv=permutation_inverse(sigma)
            #Compose the permutations
            perm=compose(sigma_inv, tau)

            #Compute cycle structure, i.e. lengths of each cycle in the cycle
            #factorization of the permutatation
            cyclestruct = cycle_structure(perm)
            Qr = Int64[n+1 for n in Qidx]
            Qpr= Int64[Qpidx[n]+1 for n in data(tau)]
            #Consolidate deltas
            Deltas=Dict()
            for i=1:n
                Deltas[Qr[i]] = Qpidx[data(sigma)[i]]
                Deltas[Qidx[i]] =  Qpr[data(tau)[i]]
            end
            #Print deltas
            print("V(", cyclestruct, ") ")
            for k in keys(Deltas)
                print("δ(",k,",",Deltas[k],") ")
            end
            for (Symb, col_idx, row_idx) in Others
                print(Symb,"(",col_idx,",",row_idx,") ")
            end
            println()
            print("= V(", cyclestruct, ") ")
            #Evaluate deltas over the indices of the other matrices
            ReindexedSymbols = []
            for (Symb, col_idx, row_idx) in Others
                new_col, new_row = col_idx, row_idx
                if has(Deltas, col_idx) new_col = Deltas[col_idx] end
                if has(Deltas, row_idx) new_row = Deltas[row_idx] end
                print(Symb,"(",new_col,",",new_row,") ")
                ReindexedSymbols = [ReindexedSymbols; (Symb, new_col, new_row)]
            end
            println()
            #Parse coefficient
            Coefficient= WeingartenUnitary(perm)
            #Reconstruct expression
            println("START PARSING")
            println("The term is =", ReindexedSymbols[end])
            Symb, left_idx, right_idx = pop!(ReindexedSymbols)
            Expression=[Any[Symbol[Symb], left_idx, right_idx]]
            while length(ReindexedSymbols) > 0
                pop_idx = expr_idx = do_transpose = is_left = nothing
                for expr_iter in enumerate(Expression)
                    expr_idx, expr_string = expr_iter
                    _, left_idx, right_idx = expr_string
                    for iter in enumerate(ReindexedSymbols)
                        idx, data = iter
                        Symb, col_idx, row_idx = data
                        if row_idx == left_idx
                            pop_idx = idx
                            do_transpose = false
                            is_left = true
                            println("Case A")
                            break
                        elseif col_idx == right_idx
                            pop_idx = idx
                            do_transpose = false
                            is_left = false
                            println("Case B")
                            break
                        elseif row_idx == right_idx
                            pop_idx = idx
                            do_transpose = true
                            is_left = false
                            println("Case C")
                            break
                        elseif col_idx == left_idx
                            pop_idx = idx
                            do_transpose = true
                            is_left = true
                            println("Case D")
                            break
                        end
                    end
                    if pop_idx != nothing break end
                end
                println("Terms left = ", ReindexedSymbols)
                if pop_idx == nothing #Found nothing, start new expression blob
                println("New word: The term is ", ReindexedSymbols[1])
                    Symb, left_idx, right_idx = delete!(ReindexedSymbols, 1)
                    push!(Expression, [Symbol[Symb], left_idx, right_idx])
                else #Found something
                    println("The term is =", ReindexedSymbols[pop_idx])
                    Symb, col_idx, row_idx = delete!(ReindexedSymbols, pop_idx)
                    Term = do_transpose ? Expr(symbol("'"), Symb) : Symb
                    if is_left
                        insert!(Expression[expr_idx][1], 1, Term)
                        Expression[expr_idx][2] = do_transpose ? row_idx : col_idx
                    else
                        push!(Expression[expr_idx][1], Term)
                        Expression[expr_idx][3] = do_transpose ? col_idx : row_idx
                    end
                end
            end
            println("DONE PARSING TREE")

            # Evaluate closed cycles
            NewExpression=Any[]
            for ExprChunk in Expression
                ExprBlob, start_idx, end_idx = ExprChunk
                print(ExprChunk, " => ")
                if start_idx == end_idx #Have a cycle; this is a trace
                    ex= length(ExprBlob)==1 ? :(trace($(ExprBlob[1]))) :
                        Expr(:call, :trace, Expr(:call, :*, ExprBlob...))
                else #Not a cycle, regular chain of multiplications
                    ex= length(ExprBlob)==1 ? ExprBlob[1] :
                        Expr(:call, :*, ExprBlob...)
                end
                push!(NewExpression, ex)
                println(ex)
            end
            ex = Expr(:call, :*, Coefficient, NewExpression...)
            println("Final expression is: ", ex)
            push!(AllTerms, ex)
        end
    end
    X = length(AllTerms)==1 ? AllTerms[1] : Expr(:call, :+, AllTerms...)
    println("THE ANSWER IS ", X)
    eval(X)
end

#Computes the Weingarten function for permutations
function WeingartenUnitary(P::Ptr{gsl_permutation})
    C = cycle_structure(P)
    WeingartenUnitary(C)
end

#Computes the Weingarten function for partitions
function WeingartenUnitary(P::partition)
    n = sum(P)
    m = length(P)
    thesum = 0.0
    for irrep in @task part(n)
        #Character of the partition
        S = character(irrep, P)
        #Character of the identity divided by n!
        T = character_identity(irrep)
        #Denominator f_r(N) of (2.10)
        f = prod([factorial(BigInt(N + P[i] - i)) / factorial(BigInt(N - i))
                for i=1:m])
        thesum += S*T/(factorial(n)*f)
        println(irrep, " ", thesum)
    end
    thesum
end


text:
using DynamicalSystems

# Short time series of two coupled logistic maps (x drives y)
sys = logistic2_unidir(c_xy = 0.5)
npts = 40
x, y = columns(trajectory(sys, npts, Ttr = 200))

# Add some uncertainties to the time series
uvals_x = [UncertainValue(Normal, X, abs(rand(Normal(0, 0.05)))) for X in x]
uvals_y = [UncertainValue(Normal, Y, abs(rand(Normal(0, 0.05)))) for Y in y]
UVX = UncertainValueDataset(uvals_x)
UVY = UncertainValueDataset(uvals_y)

# average over a few different binning schemes, use a few different
# prediction lags and use the transfer operator grid estimator on the point
# cloud generated from the invariant measure over the triangulation.
binnings = [RectangularBinning(i) for i = 3:4]
estimator = TransferOperatorGrid()
ηs = [1]

# Perform causality test, both from x to y and y to x
test = ApproximateSimplexIntersectionTest(binning = binnings, estimator = estimator, ηs = ηs)
@test causality(uvals_x, y, test) isa Vector{T} where {T}
@test causality(uvals_x, uvals_y, test) isa Vector{T} where {T}
@test causality(x, uvals_y, test) isa Vector{T} where {T}
@test causality(UVX, y, test) isa Vector{T} where {T}
@test causality(UVX, UVY, test) isa Vector{T} where {T}
@test causality(x, UVY, test) isa Vector{T} where {T}


text:
push!(LOAD_PATH, joinpath(@__DIR__, "../src"))

using Documenter
using Actors

makedocs(
    sitename = "Actors",
    format = Documenter.HTML(),
    modules = [Actors],
    pages = [
        "Guide" => "index.md",
        "Reference" => "reference.md",
        "Contributing" => "contributing.md"
    ]
)


text:
using LinearAlgebra, DelimitedFiles, StatsBase

freqs = countmap(vec(readdlm("data_6.dat", ',', Int)))

freq_list = Int[]
for ii in 0:8
    try
        append!(freq_list, freqs[ii])
    catch
        append!(freq_list, 0)
    end
end
println(freq_list)

step = [0 1 0 0 0 0 0 0 0;
        0 0 1 0 0 0 0 0 0;
        0 0 0 1 0 0 0 0 0;
        0 0 0 0 1 0 0 0 0;
        0 0 0 0 0 1 0 0 0;
        0 0 0 0 0 0 1 0 0;
        1 0 0 0 0 0 0 1 0;
        0 0 0 0 0 0 0 0 1;
        1 0 0 0 0 0 0 0 0]

println(sum(step^80 * freq_list))
println(sum(step^256 * freq_list))


text:
using PackageTest
using Test

@testset "PackageTest.jl" begin
    @test squared(4)==16
    @test squared(0)==0
    @test squared(-2)==4
end

max_stars_repo_name: jannittayao/PackageTest

max_issues_repo_name: jannittayao/PackageTest

max_forks_repo_name: jannittayao/PackageTest


text:
function deleteResource(ctx, subscription_id::String, resource_group_name::String, diskuri::String; kwargs...)
    storage_account, key = extract_account_and_key(ctx, subscription_id, resource_group_name, diskuri)

    req = Azure.REST.ServiceRequest(storage_account, "DELETE", diskuri)
    resp = Azure.REST.execute(req, key; kwargs...)
    Azure.REST.issuccess(resp)
end

function extract_account_and_key(ctx, subscription_id::String, resource_group_name::String, diskuri::String)
    # find the storage account
    storage_account = String(split(split(diskuri, "/")[3], '.'; limit=2)[1])

    # get keys for it
    allkeys = storageAccountsListKeys(Azure.api(ctx, StorageAccountsApi), resource_group_name, storage_account, apiver(StorageAccountsApi), subscription_id)
    key = (allkeys.keys[1]).value

    storage_account, key
end


text:
using CSV
using DataFrames
using Plots

function ansatz(x,α)
	return exp(-α*x^2)
end

function probability_density(x,α)
	return (ansatz(x,α)^2)/(sqrt(pi/α))
end

function E_local(x,α)
	return α + (x^2)*(1/2 - 2*α^2)
end

function metropolis(N,α)
	L = 3/sqrt(2*α)
	x = rand(-L:0.000001:L)
	E = 0
	E2 = 0
	Eln_avg = 0
	ln_avg = 0
	rejection_ratio = 0
	#start the loop
	for i in 1:N
		x_t = x + 0.4*(2*rand()*L - L)
		if probability_density(x_t,α) >= probability_density(x,α)
			x = x_t
		else
			dummy = rand()
			if dummy < probability_density(x_t,α)/probability_density(x,α)
				x = x_t
			else
				rejection_ratio += 1/N
			end
		end

		E += E_local(x,α)/N
		E2 += E_local(x,α)^2/N
		Eln_avg += (-1*E_local(x,α)*(x^2))/N
		ln_avg += (-x^2)/N
	end
	return E,E2,Eln_avg,ln_avg,rejection_ratio
end

function iterate_alpha()
	α = 0.1
	alpha_iter = 20
	N = 500
	random_walkers = 200
	γ = 0.9
	energy = Vector{Float64}(undef,alpha_iter)
	alpha = Vector{Float64}(undef,alpha_iter)
	variance = Vector{Float64}(undef,alpha_iter)
	En_analytical = Vector{Float64}(undef,alpha_iter)
	Var_analytical = Vector{Float64}(undef,alpha_iter)
	for i in 1:alpha_iter
		E = 0
		E2 = 0
		dE_dalpha = 0
		Eln = 0
		ln = 0
		rejection_ratio = 0
		for j in 1:random_walkers
			E_met, E2_met, Eln_met, ln_met, rejections_met = metropolis(N,α)
			E += E_met/random_walkers
			E2 += E2_met/random_walkers
			Eln += Eln_met/random_walkers
			ln += ln_met/random_walkers
			rejection_ratio += rejections_met/random_walkers
		end
		E_analytical = α/2 + 1/8/α
		var_analytical = (1-4*α^2)^2/32/α^2

		# Define the next α
		# dE_dalpha = 2*(Eln-E*ln)
		α = α + 0.05

		# Update the arrays for plotting
		energy[i] = E
		alpha[i] = α
		variance[i] = E2 - E^2
		En_analytical[i] = E_analytical
		Var_analytical[i] = var_analytical
	end
	return energy,alpha,variance,En_analytical,Var_analytical
end

energy,alpha,variance,E_analytical,var_analytical = iterate_alpha()
CSV.write(datadir("sims/harmonic_oscillator.csv"),(Energy = energy, Variance = variance, Alpha = alpha))

energy_plot = plot()
plot!(energy_plot,energy,
	label = "calculated expectation of energy",
	title = "Energy Plots",
	xlabel = "timestep",
	ylabel = "<E>",
	linecolor = :blue,
)
plot!(energy_plot,E_analytical,label = "Analytical Energy")

safesave(plotsdir()*"/energy_ho.png",energy_plot)

variance_plot = plot()
plot(variance_plot,variance,
	title = "Variance in energy",
	label = false,
	ylabel = "Var(E)",
	xlabel = "timestep",
	linecolor = :blue
     )

safesave(plotsdir()*"/variance_ho.png",variance_plot)

energy_vs_alpha = plot()
scatter(energy_vs_alpha,alpha,energy,
	label = "Calculated energy",
	title = "Energy vs α",
	markercolor = :blue,
    xlabel = "α",
	ylabel = "<E>"
)

scatter!(energy_vs_alpha,alpha,
	E_analytical,
	label = "Calculated energy",
	title = "Energy vs α",
	markercolor = :red
)
safesave(plotsdir()*"/energy_vs_alpha_ho.png",energy_vs_alpha)

max_stars_repo_name: Ved-Mahajan/Variational-Monto-Carlo

max_issues_repo_name: Ved-Mahajan/Variational-Monto-Carlo

max_forks_repo_name: Ved-Mahajan/Variational-Monto-Carlo


text:
# TODO: finish implementing AbstractArray interface.
# TODO: better bounds checking for getting and setting

"""Sparse array with a band of nonzeroes."""
mutable struct BandedArray{T} <: AbstractArray{T,2}
    data::Array{T,2}
    nrows::Int
    ncols::Int
    bandwidth::Int
    # offsets of band
    h_offset::Int  # how many more columns than rows
    v_offset::Int  # rows than columns
    # limits for determining if [i, j] is in band
    lower::Int
    upper::Int
    # used for initial creation and resizing
    row_padding::Int  # set by user
    col_padding::Int  # set by user
    default::T
    initialize::Bool

    function BandedArray{T}(data::Array{T,2}, shape::Tuple{Int,Int},
                            bandwidth::Int,
                            row_padding::Int, col_padding::Int,
                            default::T, initialize::Bool) where T
        if bandwidth < 1
            error("bandwidth must be positive")
        end
        nrows, ncols = shape
        drows = ndatarows(nrows, ncols, bandwidth)
        if size(data) != (drows + row_padding, ncols + col_padding)
            error("data is wrong shape")
        end
        h_offset = max(ncols - nrows, 0)
        v_offset = max(nrows - ncols, 0)
        lower, upper = bandlimits(nrows, ncols, bandwidth)
        new(data, nrows, ncols, bandwidth,
            h_offset, v_offset,
            lower, upper, row_padding, col_padding,
            default, initialize)
    end
end

function bandlimits(nrows, ncols, bandwidth)
    if ncols > nrows
        lower = nrows - ncols - bandwidth
        upper = bandwidth
    else
        lower = -bandwidth
        upper = nrows - ncols + bandwidth
    end
    return lower, upper
end

function BandedArray(T::Type, shape::Tuple{Int,Int}, bandwidth::Int;
                     row_padding::Int=0, col_padding::Int=0,
                     default=zero(T), initialize::Bool=true)
    nrows, ncols = shape
    data = allocate_data(T, nrows, ncols, bandwidth, row_padding, col_padding, initialize)
    return BandedArray{T}(data, shape, bandwidth, row_padding, col_padding, default, initialize)
end

function allocate_data(T::Type, nrows::Int, ncols::Int, bandwidth::Int,
                       row_padding::Int, col_padding::Int, initialize::Bool)
    ndrows = ndatarows(nrows, ncols, bandwidth)
    if initialize
        data = zeros(T, ndrows + row_padding, ncols + col_padding)
    else
        data = Array{T}(ndrows + row_padding, ncols + col_padding)
    end
    return data
end

function reallocate!(A::BandedArray{T}) where {T}
    # TODO: copy if requested
    A.data = allocate_data(T, A.nrows, A.ncols, A.bandwidth,
                           A.row_padding, A.col_padding, A.initialize)
end

function Base.resize!(A::BandedArray{T}, shape::Tuple{Int,Int}) where {T}
    nrows, ncols = shape
    A.nrows = nrows
    A.ncols = ncols
    A.h_offset = max(ncols - nrows, 0)
    A.v_offset = max(nrows - ncols, 0)
    A.lower, A.upper = bandlimits(nrows, ncols, A.bandwidth)

    drows, dcols = size(A.data)
    new_ndrows = ndatarows(nrows, ncols, A.bandwidth)
    if new_ndrows > drows || ncols > dcols
        reallocate!(A)
    end
end

function newbandwidth!(A::BandedArray{T}, bandwidth::Int) where {T}
    A.bandwidth = bandwidth
    reallocate!(A)
end

"""Number of used data rows"""
function ndatarows(nrows::Int, ncols::Int, bandwidth::Int)
    length_diff = abs(nrows - ncols)
    return 2 * bandwidth + length_diff + 1
end

Base.size(A::BandedArray) = (A.nrows, A.ncols)

"""The row in `data` which contains element [i, j]"""
function data_row(A::BandedArray, i::Int, j::Int)
    if !inband(A, i, j)
        error("[$i, $j] is not in band")
    end
    return (i - j) + A.h_offset + A.bandwidth + 1
end

function Base.getindex(A::BandedArray{T}, i::Int, j::Int) where {T}
    if inband(A, i, j)
        return A.data[data_row(A, i, j), j]
    else
        return A.default
    end
end

function Base.setindex!(A::BandedArray{T}, v, i::Int, j::Int) where {T}
    if inband(A, i, j)
        A.data[data_row(A, i, j), j] = v
    else
        error("Cannot set out-of-band element [$i, $j].")
    end
end

"""The rows in A's column `j` that are dense"""
function row_range(A::BandedArray, j::Int)
    start = max(1, j - A.h_offset - A.bandwidth)
    stop = min(j + A.v_offset + A.bandwidth, A.nrows)
    return start, stop
end

"""The rows in data's column `j` that can be used"""
function data_row_range(A::BandedArray, j::Int)
    a, b = row_range(A, j)
    return data_row(A, a, j), data_row(A, b, j)
end

"""Get a view of the in-band elements of column `j`."""
function sparsecol(A::BandedArray, j::Int)
    start, stop = data_row_range(A, j)
    return view(A.data, start:stop, j)
end

"""Is [i, j] in the banded region?"""
function inband(A::BandedArray, i::Int, j::Int)
    if i < 1 || j < 1 || i > A.nrows || j > A.ncols
        return false
    end
    return A.lower <= i - j <= A.upper
end

"""Return a dense representation of A"""
function Base.full(A::BandedArray{T}) where {T}
    result = zeros(T, size(A))
    for j = 1:A.ncols
        start, stop = row_range(A, j)
        dstart, dstop = data_row_range(A, j)
        result[start:stop, j] = A.data[dstart:dstop, j]
    end
    return result
end

"""Reverse the rows and the columns of `A`.

If `A` has shape `m` by `n`, then element [i, j] goes to
position [m - i, n - j].

"""
function flip!(A::BandedArray{T}) where {T}
    nrows = ndatarows(A.nrows, A.ncols, A.bandwidth)
    a, b = divrem(A.ncols, 2)
    for j in 1:a
        for i in 1:nrows
            first = A.data[i, j]
            second = A.data[nrows - i + 1, A.ncols - j + 1]
            A.data[i, j] = second
            A.data[nrows - i + 1, A.ncols - j + 1] = first
        end
    end
    # handle the middle column
    if b == 1
        c = div(nrows, 2)
        j = a + 1
        for i in 1:c
            first = A.data[i, j]
            second = A.data[nrows - i + 1, A.ncols - j + 1]
            A.data[i, j] = second
            A.data[nrows - i + 1, A.ncols - j + 1] = first
        end
    end
end

"""
    equal_ranges(a_range, b_range)

Takes the true ranges of two sub-columns of a banded array and returns
the start:stop range for their overlapping elements.

# Examples
```julia-repl
julia> equal_ranges((3, 5), (4, 6))
((2, 3), (1, 2))

julia> equal_ranges((1, 5), (1, 2))
((1, 2), (1, 2))

julia> equal_ranges((1, 5), (4, 5))
((4, 5), (1, 2))

```

"""
function equal_ranges(a_range::Tuple{Int,Int},
                      b_range::Tuple{Int,Int})
    a_start, a_stop = a_range
    b_start, b_stop = b_range
    alen = a_stop - a_start + 1
    blen = b_stop - b_start + 1
    amin = max(b_start - a_start + 1, 1)
    amax = alen - max(a_stop - b_stop, 0)
    bmin = max(a_start - b_start + 1, 1)
    bmax = blen - max(b_stop - a_stop, 0)
    return (amin, amax), (bmin, bmax)
end


text:
# Test Tools
# ==========

using Test: Test
using Random: seed!, randstring

TEST_RANDOM_SEED = 12345

function test_roundtrip_read(encoder, decoder)
    seed!(TEST_RANDOM_SEED)
    for n in vcat(0:30, sort!(rand(500:100_000, 30))), alpha in (0x00:0xff, 0x00:0x0f)
        data = rand(alpha, n)
        file = IOBuffer(data)
        stream = decoder(encoder(file))
        Test.@test hash(read(stream)) == hash(data)
        close(stream)
    end
end

function test_roundtrip_write(encoder, decoder)
    seed!(TEST_RANDOM_SEED)
    for n in vcat(0:30, sort!(rand(500:100_000, 30))), alpha in (0x00:0xff, 0x00:0x0f)
        data = rand(alpha, n)
        file = IOBuffer()
        stream = encoder(decoder(file))
        write(stream, data, TOKEN_END); flush(stream)
        Test.@test hash(take!(file)) == hash(data)
        close(stream)
    end
end

function test_roundtrip_transcode(encode, decode)
    seed!(TEST_RANDOM_SEED)
    encoder = encode()
    initialize(encoder)
    decoder = decode()
    initialize(decoder)
    for n in vcat(0:30, sort!(rand(500:100_000, 30))), alpha in (0x00:0xff, 0x00:0x0f)
        data = rand(alpha, n)
        Test.@test hash(transcode(decode, transcode(encode, data))) == hash(data)
        Test.@test hash(transcode(decoder, transcode(encoder, data))) == hash(data)
    end
    finalize(encoder)
    finalize(decoder)
end

function test_roundtrip_lines(encoder, decoder)
    seed!(TEST_RANDOM_SEED)
    lines = String[]
    buf = IOBuffer()
    stream = encoder(buf)
    for i in 1:100_000
        line = randstring(rand(0:1000))
        println(stream, line)
        push!(lines, line)
    end
    write(stream, TOKEN_END)
    flush(stream)
    seekstart(buf)
    Test.@test hash(lines) == hash(readlines(decoder(buf)))
end

function test_roundtrip_fileio(Encoder, Decoder)
    data = b"""
    Lorem ipsum dolor sit amet, consectetur adipiscing elit. Nulla sit amet tempus felis. Etiam molestie urna placerat iaculis pellentesque. Maecenas porttitor et dolor vitae posuere. Lorem ipsum dolor sit amet, consectetur adipiscing elit. Nunc eget nibh quam. Nullam aliquet interdum fringilla. Duis facilisis, lectus in consectetur varius, lorem sem tempor diam, nec auctor tellus nibh sit amet sapien. In ex nunc, elementum eget facilisis ut, luctus eu orci. Sed sapien urna, accumsan et elit non, auctor pretium massa. Phasellus consectetur nisi suscipit blandit aliquam. Nulla facilisi. Mauris pellentesque sem sit amet mi vestibulum eleifend. Nulla faucibus orci ac lorem efficitur, et blandit orci interdum. Aenean posuere ultrices ex sed rhoncus. Donec malesuada mollis sem, sed varius nunc sodales sed. Curabitur lobortis non justo non tristique.
    """
    mktemp() do filename, file
        stream = TranscodingStream(Encoder(), file)
        write(stream, data)
        close(stream)
        stream = TranscodingStream(Decoder(), open(filename))
        Test.@test hash(read(stream)) == hash(data)
        close(stream)
    end
end

function test_chunked_read(Encoder, Decoder)
    seed!(TEST_RANDOM_SEED)
    alpha = b"色即是空"
    encoder = Encoder()
    initialize(encoder)
    for _ in 1:500
        chunks = [rand(alpha, rand(0:100)) for _ in 1:rand(1:100)]
        data = mapfoldl(x->transcode(encoder, x), vcat, chunks, init=UInt8[])
        buffer = NoopStream(IOBuffer(data))
        ok = true
        for chunk in chunks
            stream = TranscodingStream(Decoder(), buffer, stop_on_end=true)
            ok &= hash(read(stream)) == hash(chunk)
            ok &= eof(stream)
            close(stream)
        end
        Test.@test ok
    end
    finalize(encoder)
end

function test_chunked_write(Encoder, Decoder)
    seed!(TEST_RANDOM_SEED)
    alpha = b"空即是色"
    encoder = Encoder()
    initialize(encoder)
    for _ in 1:500
        chunks = [rand(alpha, rand(0:100)) for _ in 1:2]
        data = map(x->transcode(encoder, x), chunks)
        buffer = IOBuffer()
        stream = TranscodingStream(Decoder(), buffer, stop_on_end=true)
        write(stream, vcat(data...))
        flush(stream)
        ok = true
        ok &= hash(take!(buffer)) == hash(chunks[1])
        ok &= buffersize(stream.state.buffer1) == length(data[2])
        Test.@test ok
        close(stream)
    end
    finalize(encoder)
end


text:
__precompile__()

module BilevelTests

export test_bilevel_learn

using Printf
using SparseArrays
using FileIO
using ColorTypes: Gray
import ColorVectorSpace
using ImageContrastAdjustment

using AlgTools.Util
using AlgTools.LinOps
using AlgTools.LinkedLists
using ImageTools.Visualise

using VariationalImaging

using VariationalImaging.Util
using VariationalImaging.GradientOps
using VariationalImaging.OpDenoise

using VariationalImaging.TestDatasets
using VariationalImaging.Bilevel

# Parameters
const default_save_prefix="bilevel_result_"

const default_params = (
    verbose_iter = 1,
    maxiter = 20,
    save_results = true,
    dataset_name = "cameraman128_5",
    save_iterations = false
)

const bilevel_params = (
    η₁ = 0.25,
    η₂ = 0.75,
    β₁ = 0.25,
    β₂ = 2.0,
    Δ₀ = 1.0,
    x₀ = 0.001
)

function save_results(params, b, b_data, x, opt_img,st)
    if params.save_results
        perffile = params.save_prefix * ".txt"
        println("Saving " * perffile)
        write_log(perffile, st.log, "# params = $(params), x = $x\n")
        fn = (t, ext) -> "$(params.save_prefix)_$(t).$(ext)"
        save(File(format"PNG", fn("true", "png")), grayimg(b))
        save(File(format"PNG", fn("data", "png")), grayimg(b_data))
        save(File(format"PNG", fn("reco", "png")), grayimg(opt_img))
    end
end

###############
# Learning function
###############
function learning_function(x,data)
    op = FwdGradientOp()
    u = denoise(data[2],x,op)
    cost = 0.5*norm₂²(u-data[1])
    grad = gradient(x,u,data[1],op)
    return u,cost,grad
end

function denoise(data,x,op::LinOp)
    denoise_params = (
        ρ = 0,
        α = x,
        op=op,
        # PDPS
        τ₀ = 5,
        σ₀ = 0.99/5,
        accel = true,
        save_results = false,
        maxiter = 1000,
        verbose_iter = 1001,
        save_iterations = false
    )
    st_opt, iterate_opt = initialise_visualisation(false)
    opt_img = op_denoise_pdps(data; iterate=iterate_opt, params=denoise_params)
    finalise_visualisation(st_opt)
    return opt_img
end

function gradient(α,u::AbstractArray{T,3},ū::AbstractArray{T,3},op::LinOp) where T

	M,N,O = size(u)
	grad = 0
	for i = 1:O
		u1 = @view u[:,:,i]
		u2 = @view ū[:,:,i]
		g = gradient(α,u1,u2,op)
		grad += g
	end
	return grad
end

function gradient(α,u::AbstractArray{T,2},ū::AbstractArray{T,2},op::LinOp) where T
	u = Float64.(Gray{Float64}.(u))
	ū = Float64.(Gray{Float64}.(ū))
	# Obtain Active and inactive sets
	n = size(u,1)

	# Generate centered gradient matrix
    G = matrix(op,n)
	#G = createDivMatrix(n)
	Gu = G*u[:]
	nGu = xi(Gu)
	act = nGu .< 1e-12
	inact = 1 .- act
	Act = spdiagm(0=>act)
	Inact = spdiagm(0=>inact)

	# Vector with grad norm in inactive components and one in the active
	den = Inact*nGu+act
	Den = spdiagm(0=>1 ./den)

	# prod KuKuᵗ/norm³
	prodKuKu = prodesc(Gu ./den.^3,Gu)

	Adj = [spdiagm(0=>ones(n^2)) α*G';
			Act*G+Inact*(prodKuKu-Den)*G Inact+eps()*Act]

	Track=[(u[:]-ū[:]);zeros(2*n^2)]
	mult = Adj\Track
	p = @view mult[1:n^2]
	return -p'*(G'*Inact*Den*Gu)
end

###############
# Bilevel learn test
###############

function test_bilevel_learn(;
    visualise=true,
    save_prefix=default_save_prefix,
    kwargs...)

    # Parameters for this experiment
    params = default_params ⬿ bilevel_params ⬿ kwargs
    params = params ⬿ (save_prefix = save_prefix * "denoise_" * params.dataset_name,)

    # Load dataset
    b,b_noisy = TestDatasets.testdataset(params.dataset_name)
    b = Float64.(Gray{Float64}.(b))
    b_noisy = Float64.(Gray{Float64}.(b_noisy))

    # Launch (background) visualiser
    st, iterate = initialise_visualisation(visualise)

    # Run algorithm
    x, u, st = bilevel_learn((b,b_noisy),learning_function; xinit=params.x₀,iterate=iterate, params=params)

    adjust_histogram!(u,LinearStretching())

    save_results(params, b, b_noisy, x, u, st)

    # Exit background visualiser
    finalise_visualisation(st)
end

end

max_stars_repo_name: dvillacis/BPLDenoising

max_issues_repo_name: dvillacis/BPLDenoising

max_forks_repo_name: dvillacis/BPLDenoising


text:
__precompile__(true)

"""
Main module for `Documenter.jl` -- a documentation generation package for Julia.

Two functions are exported from this module for public use:

- [`makedocs`](@ref). Generates documentation from docstrings and templated markdown files.
- [`deploydocs`](@ref). Deploys generated documentation from *Travis-CI* to *GitHub Pages*.

Additionally it provides the unexported [`Documenter.generate`](@ref), which can be used to
generate documentation stubs for new packages.

$(EXPORTS)

"""
module Documenter

using Compat, DocStringExtensions

#
# Submodules
#
# All submodules of this package are declared in the following loop.
#
# They can either be a single file in the "modules" folder or a subfolder of
# "modules" containing a file with the same name. Large submodules should be
# split into several files in a subfolder.
#
for mod in [
    "Utilities",
    "DocSystem",
    "Selectors",
    "Formats",
    "Anchors",
    "Documents",
    "Builder",
    "Expanders",
    "Walkers",
    "CrossReferences",
    "DocChecks",
    "Writers",
    "Deps",
    "Generator",
]
    dir = dirname(@__FILE__)
    file = joinpath(dir, mod * ".jl")
    isfile(file) ? include(file) : include(joinpath(dir, mod, mod * ".jl"))
end

# User Interface.
# ---------------

export Deps, makedocs, deploydocs

"""
    makedocs(
        root    = "<current-directory>",
        source  = "src",
        build   = "build",
        clean   = true,
        doctest = true,
        modules = Module[],
        repo    = "",
    )

Combines markdown files and inline docstrings into an interlinked document.
In most cases [`makedocs`](@ref) should be run from a `make.jl` file:

```julia
using Documenter
makedocs(
    # keywords...
)
```

which is then run from the command line with:

```sh
\$ julia make.jl
```

The folder structure that [`makedocs`](@ref) expects looks like:

    docs/
        build/
        src/
        make.jl

# Keywords

**`root`** is the directory from which `makedocs` should run. When run from a `make.jl` file
this keyword does not need to be set. It is, for the most part, needed when repeatedly
running `makedocs` from the Julia REPL like so:

    julia> makedocs(root = Pkg.dir("MyPackage", "docs"))

**`source`** is the directory, relative to `root`, where the markdown source files are read
from. By convention this folder is called `src`. Note that any non-markdown files stored
in `source` are copied over to the build directory when [`makedocs`](@ref) is run.

**`build`** is the directory, relative to `root`, into which generated files and folders are
written when [`makedocs`](@ref) is run. The name of the build directory is, by convention,
called `build`, though, like with `source`, users are free to change this to anything else
to better suit their project needs.

**`clean`** tells [`makedocs`](@ref) whether to remove all the content from the `build`
folder prior to generating new content from `source`. By default this is set to `true`.

**`doctest`** instructs [`makedocs`](@ref) on whether to try to test Julia code blocks
that are encountered in the generated document. By default this keyword is set to `true`.
Doctesting should only ever be disabled when initially setting up a newly developed package
where the developer is just trying to get their package and documentation structure correct.
After that, it's encouraged to always make sure that documentation examples are runnable and
produce the expected results. See the [Doctests](@ref) manual section for details about
running doctests.

**`modules`** specifies a vector of modules that should be documented in `source`. If any
inline docstrings from those modules are seen to be missing from the generated content then
a warning will be printed during execution of [`makedocs`](@ref). By default no modules are
passed to `modules` and so no warnings will appear. This setting can be used as an indicator
of the "coverage" of the generated documentation.
For example Documenter's `make.jl` file contains:

```julia
makedocs(
    modules = [Documenter],
    # ...
)
```

and so any docstring from the module `Documenter` that is not spliced into the generated
documentation in `build` will raise a warning.

**`repo`** specifies a template for the "link to source" feature. If you are
using GitHub, this is automatically generated from the remote. If you are using
a different host, you can use this option to tell Documenter how URLs should be
generated. The following placeholders will be replaced with the respective
value of the generated link:

  - `{commit}` Git commit id
  - `{path}` Path to the file in the repository
  - `{line}` Line (or range of lines) in the source file

For example if you are using GitLab.com, you could use

```julia
makedocs(repo = "https://gitlab.com/user/project/blob/{commit}{path}#L{line}")
```

# See Also

A guide detailing how to document a package using Documenter's [`makedocs`](@ref) is provided
in the [Usage](@ref) section of the manual.
"""
function makedocs(; debug = false, args...)
    document = Documents.Document(; args...)
    cd(document.user.root) do
        Selectors.dispatch(Builder.DocumentPipeline, document)
    end
    debug ? document : nothing
end

"""
    deploydocs(
        root   = "<current-directory>",
        target = "site",
        repo   = "<required>",
        branch = "gh-pages",
        latest = "master",
        osname = "linux",
        julia  = "nightly",
        deps   = <Function>,
        make   = <Function>,
    )

Converts markdown files generated by [`makedocs`](@ref) to HTML and pushes them to `repo`.
This function should be called from within a package's `docs/make.jl` file after the call to
[`makedocs`](@ref), like so

```julia
using Documenter, PACKAGE_NAME
makedocs(
    # options...
)
deploydocs(
    repo = "github.com/..."
)
```

# Keywords

**`root`** has the same purpose as the `root` keyword for [`makedocs`](@ref).

**`target`** is the directory, relative to `root`, where generated HTML content should be
written to. This directory **must** be added to the repository's `.gitignore` file. The
default value is `"site"`.

**`repo`** is the remote repository where generated HTML content should be pushed to. This
keyword *must* be set and will throw an error when left undefined. For example this package
uses the following `repo` value:

```julia
repo = "github.com/JuliaDocs/Documenter.jl.git"
```

**`branch`** is the branch where the generated documentation is pushed. By default this
value is set to `"gh-pages"`.

**`latest`** is the branch that "tracks" the latest generated documentation. By default this
value is set to `"master"`.

**`osname`** is the operating system which will be used to deploy generated documentation.
This defaults to `"linux"`. This value must be one of those specified in the `os:` section
of the `.travis.yml` configuration file.

**`julia`** is the version of Julia that will be used to deploy generated documentation.
This defaults to `"nightly"`. This value must be one of those specified in the `julia:`
section of the `.travis.yml` configuration file.

**`deps`** is the function used to install any dependencies needed to build the
documentation. By default this function installs `pygments` and `mkdocs` using the
[`Deps.pip`](@ref) function:

```julia
deps = Deps.pip("pygments", "mkdocs")
```

**`make`** is the function used to convert the markdown files to HTML. By default this just
runs `mkdocs build` which populates the `target` directory.

# See Also

The [Hosting Documentation](@ref) section of the manual provides a step-by-step guide to
using the [`deploydocs`](@ref) function to automatically generate docs and push then to
GitHub.
"""
function deploydocs(;
        root   = Utilities.currentdir(),
        target = "site",

        repo   = error("no 'repo' keyword provided."),
        branch = "gh-pages",
        latest = "master",

        osname = "linux",
        julia  = "nightly",

        deps   = Deps.pip("pygments", "mkdocs"),
        make   = () -> run(`mkdocs build`),
    )
    # Get environment variables.
    github_api_key      = get(ENV, "GITHUB_API_KEY",       "")
    travis_branch       = get(ENV, "TRAVIS_BRANCH",        "")
    travis_pull_request = get(ENV, "TRAVIS_PULL_REQUEST",  "")
    travis_repo_slug    = get(ENV, "TRAVIS_REPO_SLUG",     "")
    travis_tag          = get(ENV, "TRAVIS_TAG",           "")
    travis_osname       = get(ENV, "TRAVIS_OS_NAME",       "")
    travis_julia        = get(ENV, "TRAVIS_JULIA_VERSION", "")

    # Other variables.
    sha          = readchomp(`git rev-parse --short HEAD`)
    ssh_key_file = abspath(joinpath(root, ".documenter.enc"))
    has_ssh_key  = isfile(ssh_key_file)

    # When should a deploy be attempted?
    should_deploy =
        contains(repo, travis_repo_slug) &&
        travis_pull_request == "false"   &&
        (
            # Support token and ssh key deployments.
            github_api_key != "" ||
            has_ssh_key
        ) &&
        travis_osname == osname &&
        travis_julia  == julia  &&
        (
            travis_branch == latest ||
            travis_tag    != ""
        )

    if get(ENV, "DOCUMENTER_DEBUG", "") == "true"
        Utilities.debug("TRAVIS_REPO_SLUG       = \"$travis_repo_slug\"")
        Utilities.debug("TRAVIS_PULL_REQUEST    = \"$travis_pull_request\"")
        Utilities.debug("TRAVIS_OS_NAME         = \"$travis_osname\"")
        Utilities.debug("TRAVIS_JULIA_VERSION   = \"$travis_julia\"")
        Utilities.debug("TRAVIS_BRANCH          = \"$travis_branch\"")
        Utilities.debug("TRAVIS_TAG             = \"$travis_tag\"")
        Utilities.debug("git commit SHA         = $sha")
        Utilities.debug("GITHUB_API_KEY exists  = $(github_api_key != "")")
        Utilities.debug(".documenter.enc path   = $(ssh_key_file)")
        Utilities.debug(".documenter.enc exists = $(has_ssh_key)")
        Utilities.debug("should_deploy          = $should_deploy")
    end

    if should_deploy
        # Add local bin path if needed.
        Deps.updatepath!()
        # Install dependencies when applicable.
        if deps !== nothing
            Utilities.log("installing dependencies.")
            deps()
        end
        # Change to the root directory and try to deploy the docs.
        cd(root) do
            Utilities.log("setting up target directory.")
            isdir(target) || mkpath(target)
            # Run extra build steps defined in `make` if required.
            if make !== nothing
                Utilities.log("running extra build steps.")
                make()
            end
            Utilities.log("pushing new documentation to remote: $repo:$branch.")
            mktempdir() do temp
                # Versioned docs directories.
                latest_dir = joinpath(temp, "latest")
                stable_dir = joinpath(temp, "stable")
                tagged_dir = joinpath(temp, travis_tag)

                keyfile, _ = splitext(ssh_key_file)
                target_dir = abspath(target)

                # The upstream URL to which we push new content and the ssh decryption commands.
                upstream, ssh_script =
                    if has_ssh_key
                        key = getenv(r"encrypted_(.+)_key")
                        iv  = getenv(r"encrypted_(.+)_iv")
                        "git@$(replace(repo, "github.com/", "github.com:"))",
                        """
                        openssl aes-256-cbc -K $key -iv $iv -in $keyfile.enc -out $keyfile -d
                        chmod 600 $keyfile
                        eval `ssh-agent -s`
                        ssh-add $keyfile
                        """
                    else
                        "https://$github_api_key@$repo", ""
                    end

                # On non-tagged builds we just build `latest`,
                # otherwise build the `stable` and `version` builds.
                copy_script =
                    if travis_tag == ""
                        """
                        rm -rf $latest_dir
                        cp -r  $target_dir $latest_dir
                        """
                    else
                        """
                        rm -rf $stable_dir
                        cp -r  $target_dir $stable_dir
                        rm -rf $tagged_dir
                        cp -r  $target_dir $tagged_dir
                        """
                    end

                # Auto authorise SSH authentication requests for github.com.
                open(joinpath(homedir(), ".ssh", "config"), "a") do io
                    println(io,
                        """
                        Host github.com
                            StrictHostKeyChecking no
                        """
                    )
                end

                # Write, run, and delete the deploy script.
                mktemp() do path, io
                    script = buildscript(
                        temp,
                        upstream,
                        branch,
                        ssh_script,
                        copy_script,
                        sha,
                    )
                    println(io, script); flush(io) # `flush`, otherwise `path` is empty.
                    run(`sh $path`)
                    # Remove the unencrypted private key.
                    isfile(keyfile) && rm(keyfile)
                end
            end
        end
    else
        Utilities.log("skipping docs deployment.")
    end
end

buildscript(dir, upstream, branch, ssh_script, copy_script, sha) =
    """
    $ssh_script

    cd $dir

    git init

    git config user.name  "autodocs"
    git config user.email "autodocs"

    git remote add upstream "$upstream"

    git fetch upstream

    git checkout -b $branch upstream/$branch

    $copy_script

    git add -A .
    git commit -m "build based on $sha"

    git push -q upstream HEAD:$branch
    """

function getenv(k::Regex)
    found = collect(filter(s -> ismatch(k, s), keys(ENV)))
    length(found) === 1 ? ENV[found[1]] : error("no keys found in ENV 'key/iv' pair.")
end

export Travis

"""
Package functions for interacting with Travis.

$(EXPORTS)
"""
module Travis

using Compat, DocStringExtensions

export genkeys

"""
$(SIGNATURES)

Generate ssh keys for package `package` to automatically deploy docs from Travis to GitHub
pages. Uses the `remote` information to get the user and repository values.
Requires the following command lines programs to be installed:

- `which`
- `git`
- `travis`
- `ssh-keygen`

# Examples

```jlcon
julia> using Documenter

julia> Travis.genkeys("MyPackageName")
[ ... output ... ]

julia> Travis.genkeys("MyPackageName", remote="organization")
[ ... output ... ]
```
"""
function genkeys(package; remote="origin")
    # Error checking. Do the required programs exist?
    isdir(Pkg.dir(package))     || error("'$package' could not be found in '$(Pkg.dir())'.")
    success(`which which`)      || error("'which' not found.")
    success(`which git`)        || error("'git' not found.")
    success(`which travis`)     || error("'travis' not found.")
    success(`which ssh-keygen`) || error("'ssh-keygen' not found.")

    directory = "docs"
    filename  = ".documenter"

    cd(Pkg.dir(package, directory)) do
        # Get remote details.
        user, repo =
            let r = readchomp(`git config --get remote.$remote.url`)
                m = match(isdefined(Base, :LibGit2) ?
                    Base.LibGit2.GITHUB_REGEX :
                    Pkg.Git.GITHUB_REGEX, r)
                m === nothing && error("no remote repo named '$remote' found.")
                m[2], m[3]
            end

        run(`ssh-keygen -N "" -f $filename`)
        run(`travis login --auto`)
        run(`travis encrypt-file $filename -r $user/$repo`)

        warn("removing private key.")
        rm(filename)

        println(
            """

            Add the following public deploy key to '$user/$repo' with write access

            $(strip(readstring("$filename.pub")))

            on the following page:

                https://github.com/$user/$repo/settings/keys

            Then commit the '$filename.enc' file. Do not edit '.travis.yml'.
            """
        )
        warn("removing public key.")
        rm("$filename.pub")
    end
end

end

"""
$(SIGNATURES)

Creates a documentation stub for a package called `pkgname`. The location of
the documentation is assumed to be `<package directory>/docs`, but this can
be overriden with the keyword argument `dir`.

It creates the following files

```
docs/
    .gitignore
    src/index.md
    make.jl
    mkdocs.yml
```

# Arguments

**`pkgname`** is the name of the package (without `.jl`). It is used to
determine the location of the documentation if `dir` is not provided.

# Keywords

**`dir`** defines the directory where the documentation will be generated.
It defaults to `<package directory>/docs`. The directory must not exist.

# Examples

```jlcon
julia> using Documenter

julia> Documenter.generate("MyPackageName")
[ ... output ... ]
```
"""
function generate(pkgname::AbstractString; dir=nothing)
    # TODO:
    #   - set up deployment to `gh-pages`
    #   - fetch url and username automatically (e.g from git remote.origin.url)

    # Check the validity of the package name
    if length(pkgname) == 0
        error("Package name can not be an empty string.")
    end
    # Determine the root directory where we wish to generate the docs and
    # check that it is a valid directory.
    docroot = if dir === nothing
        pkgdir = Pkg.dir(pkgname)
        if !isdir(pkgdir)
            error("Unable to find package $(pkgname).jl at $(pkgdir).")
        end
        joinpath(pkgdir, "docs")
    else
        dir
    end

    if ispath(docroot)
        error("Directory $(docroot) already exists.")
    end

    # deploy the stub
    try
        info("Deploying documentation to $(docroot)")
        mkdir(docroot)

        # create the root doc files
        Generator.savefile(docroot, ".gitignore") do io
            write(io, Generator.gitignore())
        end
        Generator.savefile(docroot, "make.jl") do io
            write(io, Generator.make(pkgname))
        end
        Generator.savefile(docroot, "mkdocs.yml") do io
            write(io, Generator.mkdocs(pkgname))
        end

        # Create the default documentation source files
        Generator.savefile(docroot, "src/index.md") do io
            write(io, Generator.index(pkgname))
        end
    catch
        rm(docroot, recursive=true)
        rethrow()
    end
    nothing
end

end


text:
export @linq, linq

##############################################################################
##
## @linq - general macro that creates a mini DSL for chaining and macro and
##         function calls
##
##############################################################################
"""
    @linq df ...

!!! note
    `@linq` is deprecated. Use `@chain` instead. See `? @chain` for details.

General macro that creates a mini DSL for chaining and macro calls.

### Details

The following embedded function calls are equivalent to their macro version:

- `with`
- `where`
- `select`
- `transform`
- `by`
- `groupby`
- `orderby`
- `combine`

### Examples

```jldoctest
julia> using DataFramesMeta, Statistics

julia> df = DataFrame(
            a = repeat(1:4, outer = 2),
            b = repeat(2:-1:1, outer = 4),
            x = 1:8);

julia> x1 = @linq transform(where(df, :a .> 2, :b .!= "c"), :y = 10 .* :x);

julia> x1 = @linq by(x1, :b, :meanX = mean(:x), :meanY = mean(:y));

julia> @linq select(orderby(x1, :b, -:meanX), :var = :b, :meanX, :meanY)
2×3 DataFrame
│ Row │ var   │ meanX   │ meanY   │
│     │ Int64 │ Float64 │ Float64 │
├─────┼───────┼─────────┼─────────┤
│ 1   │ 1     │ 6.0     │ 60.0    │
│ 2   │ 2     │ 5.0     │ 50.0    │

julia> @linq df |>
           transform(y = 10 .* :x) |>
           where(:a .> 2) |>
           by(:b, :meanX = mean(:x), :meanY = mean(:y)) |>
           orderby(:meanX) |>
           select(:meanX, :meanY, var = :b)
2×3 DataFrame
│ Row │ meanX   │ meanY   │ var   │
│     │ Float64 │ Float64 │ Int64 │
├─────┼─────────┼─────────┼───────┤
│ 1   │ 5.0     │ 50.0    │ 2     │
│ 2   │ 6.0     │ 60.0    │ 1     │

```
"""
macro linq(arg)
    msg = "@linq is deprecated. Use @chain instead. See ? @chain for details"
    Base.depwarn(msg, Symbol("@linq"))
    esc(replacefuns(replacechains(arg)))
end

# Snippet from Calculus.jl
struct SymbolParameter{T} end
SymbolParameter(s::Symbol) = SymbolParameter{s}()

replacefuns(x) = x  # default for non-expression stuff
function replacefuns(e::Expr)
    for i in 1:length(e.args)
        e.args[i] = replacefuns(e.args[i])
    end
    if e.head == :call && isa(e.args[1], Symbol)
        return linq(SymbolParameter(e.args[1]), e.args[2:end]...)
    else
        return e
    end
end

replacechains(x) = x
function replacechains(e::Expr)
    for i in 1:length(e.args)
        e.args[i] = replacechains(e.args[i])
    end
    if e.head == :call && e.args[1] == :|> && isa(e.args[3], Expr)
        newe = e.args[3]
        insert!(newe.args, 2, e.args[2])
        return newe
    else
        return e
    end
end

##############################################################################
##
## Various linq helper definitions
##
##############################################################################

## Default, no-op:
linq(::SymbolParameter{s}, args...) where {s} = Expr(:call, s, args...)

function linq(::SymbolParameter{:with}, d, body)
    with_helper(d, body)
end

function linq(::SymbolParameter{:where}, d, args...)
    where_helper(d, args...)
end

function linq(::SymbolParameter{:orderby}, d, args...)
    orderby_helper(d, args...)
end

function linq(::SymbolParameter{:transform}, x, args...)
    transform_helper(x, args...)
end

function linq(::SymbolParameter{:combine}, x, args...)
    combine_helper(x, args...)
end

function linq(::SymbolParameter{:by}, x, what, args...)
    by_helper(x, what, args...)
end

function linq(::SymbolParameter{:select}, x, args...)
    select_helper(x, args...)
end


text:
using Revise
using POMDPs
using POMDPModelTools
using Compose
import Cairo

using DroneSurveillance

pomdp = DroneSurveillancePOMDP(size=(5,5))

s0 = DSState([2,2],[3,3])
c = render(pomdp, (s=s0,))

c |> PDF("dronesurveillance.pdf")

using RockSample
using Random

pomdp = RockSamplePOMDP{3}(rocks_positions=[[4,2], [2,3], [4, 4]])

s0 = initialstate(pomdp, MersenneTwister(1))

c = render(pomdp, (s=s0,a=1))

c |> PDF("rocksample.pdf")

using POMDPModels
using POMDPModelChecking
using Colors
include("blind_gridworld.jl")

mdp = SimpleGridWorld(size=(5,5), rewards=Dict(GWPos(9,3)=>10.0), terminate_from=Set([]))

const LABELLED_STATES = Dict(GWPos(2,3) => :C,  GWPos(4,4)=>:A, GWPos(5,2)=>:B)

function POMDPModelChecking.labels(mdp::SimpleGridWorld, s, a)
    if haskey(LABELLED_STATES, s)
        return tuple(LABELLED_STATES[s])
    else
        return ()
    end
end

pomdp = BlindGridWorld(exit=GWPos(4,3), simple_gw = mdp)

c = render(mdp, (s=GWPos(1,1),), valuecolor=nothing, landmark = s -> s in keys(LABELLED_STATES));

c |> PDF("gridworld.pdf")


text:
# Autogenerated wrapper script for LoopTools_jll for i686-w64-mingw32-libgfortran3
export libooptools

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("LoopTools")
JLLWrappers.@declare_library_product(libooptools, "libooptools.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libooptools,
        "bin\\libooptools.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()