text:
# using AutomotiveDrivingModels
# using NearestNeighbors

# import AutomotiveDrivingModels: get_actions!, observe!, action_context, get_name
# import Base.rand
# import PyPlot

# export
#     HRHC,
#     curveDist,
#     wrap_to_π,
#     kdProject,
#     generateObstacleMap,
#     updateObstacleMap!,
#     generateMotionMap,
#     screenCollision,
#     tailgateAvoidance,
#     getSuccessorStates,
#     loopProjectionKD,
#     computeTrajectory,
#     screenTrajectory,
#     checkCollision,
#     calculateObjective,
#     plot_stϕ,
#     plotHRHCInfo,
#     plotObjectiveHorizon,
#     plotSplineRoadway

# utility functions
function curveDist(pt1::CurvePt, pt2::CurvePt)
    d = sqrt((pt1.pos.x - pt2.pos.x)^2 + (pt1.pos.y - pt2.pos.y)^2)
end
function wrap_to_π(θ)
    θ = θ - div(θ,2*Float64(π))*(2*Float64(π))
    θ = θ + (θ .< -π).*(2*Float64(π)) - (θ .> π).*(2*Float64(π))
end
function kdProject(x,y,θ,tree,roadway,hrhc)
    """
    project single (x,y,Θ) point to roadway spline using kdtree to find the nearest spline point
    """
    curve = roadway.segments[1].lanes[1].curve
    # Δs = roadway.segments[1].lanes[1].curve[2].s - roadway.segments[1].lanes[1].curve[1].s
    idx_list,dist = knn(tree,[x;y],1)
    idx = idx_list[1]
    idxA = idx-1
    idxB = idx+1

    if idx == length(curve)
        idxB = 1 # back to the beginning of the curve
    end
    if idx == 1
        idxA = length(curve)
    end
    dA = sqrt(sum(([curve[idxA].pos.x, curve[idxA].pos.y]-[x,y]).^2))
    dB = sqrt(sum(([curve[idxB].pos.x, curve[idxB].pos.y]-[x,y]).^2))
    if dA < dB
        idxB = idx
    else
        idxA = idx
    end

    # project
    vec1 = [curve[idxB].pos.x - curve[idxA].pos.x, curve[idxB].pos.y - curve[idxA].pos.y, 0]
    vec2 = [x - curve[idxA].pos.x, y - curve[idxA].pos.y, 0]
    idx_t = dot(vec2, vec1)/norm(vec1)^2

    pθ = curve[idxA].pos.θ + idx_t*(curve[idxB].pos.θ - curve[idxA].pos.θ)

    s = curve[idxA].s + idx_t*hrhc.Δs
    t = norm(vec2 - idx_t*vec1)*sign(sum(cross(vec1, vec2)))
    ϕ = wrap_to_π(θ - pθ)

    s,t,ϕ,idxA
end

# type Hierarchical Receding Horizion Controller
type HRHC <: DriverModel{AccelDesang}
#     action_context::IntegratedContinuous
    car_ID::Int
    v_map
    δ_map
    motion_map # state changes associated with cmd = (v_command, δ_command) - this is the motion map
    v_cmd::Int # index of v in v_cmds
    δ_cmd::Int # index of δ in δ_cmds
    #car parameters with bicycle geometry model
    car_length::Float64 # wheel base
    wheel_base::Float64
    car_width::Float64
    ellipseA::Float64
    ellipseB::Float64

    a_step::Float64 # max acceleration (m/s)
    μ::Float64 # friction coefficient
    v_range # possible velocityies
    δ_range # possible steering angles
    successor_states # array of next states
    # current v, current δ
    v::Float64
    δ::Float64
    curve_ind::Int
    Δs::Float64
    # planning horizon
    h::Int
    Δt::Float64
    # logit level
    k::Int
    # maximum deviation from center of track (if |t| > T_MAX, car is out of bounds)
    T_MAX::Float64
    # Action = Next State
    action::VehicleState
    function HRHC(car_ID::Int,roadway;
        car_length::Float64=4.8,
        wheel_base::Float64=4.4,
        car_width::Float64=2.5,
        v::Float64=0.0,
        δ::Float64=0.0,
        h::Int=10,
        Δt::Float64=1.0/24,
        ΔV₊::Float64=1.55,
        ΔV₋::Float64=3.05,
        Δδ::Float64=Float64(π)/12,
        v_min::Float64=0.0,
        v_max::Float64=100.0,
        a_step::Float64=10.0,
        a_range=[-1,0,1],
        μ::Float64=20.0,
        g::Float64=9.81,
        δ_max::Float64=Float64(π)/6,
        δ_step::Float64=Float64(π)/128,
        k::Int=1
        )

        hrhc = new()

        hrhc.T_MAX=(roadway.segments[1].lanes[1].width - car_width)/2.0

        hrhc.car_ID = car_ID
        hrhc.car_length=car_length
        hrhc.wheel_base=wheel_base
        hrhc.car_width=car_width
        hrhc.h=h
        hrhc.Δt=Δt
        hrhc.a_step=a_step
        hrhc.μ=μ

        hrhc.k=k

        hrhc.motion_map, hrhc.v_map, hrhc.δ_map, hrhc.v_range, hrhc.δ_range = generateMotionMap(v_min,
            v_max,a_step,δ_max,δ_step,wheel_base,h,μ,Δt,a_range=a_range)

        hrhc.successor_states = zeros(size(hrhc.motion_map[1],1),size(hrhc.motion_map[1],2),3)

        hrhc.v=v
        hrhc.δ=δ
        hrhc.v_cmd = 1
        hrhc.δ_cmd = Int((length(hrhc.δ_range) - 1)/2)
        hrhc.curve_ind=1
        hrhc.Δs=roadway.segments[1].lanes[1].curve[2].s-roadway.segments[1].lanes[1].curve[1].s
#         hrhc.action_context=context
        hrhc.action = VehicleState(VecSE2(0,0,0),0.0)

        # calculate semimajor axes of bounding ellipse with minimal area (for collision checking)
        W = car_width/2.0
        L = car_length/2.0

        # use quadratic formula to find B^2
        a = 2
        b = -4*L - (W^2)*2*(L^2)
        c = 2*L^4

        B = sqrt((-b + sqrt(b^2 - 4*a*c))/(2*a))
        A = sqrt(((W)^2) / (1 - (L/(B))^2))

        hrhc.ellipseA = A
        hrhc.ellipseB = B

        hrhc
    end
end
function generateObstacleMap(scene, models)
    k = maximum([driver.k for (id, driver) in models])
    n = length(scene)
    h = maximum([driver.h for (id, driver) in models]) # h should be the same for all vehicles on the track
    obstacleDict = Dict()
    for level in 0:k
        obstacleDict[level] = Dict()
        for (id, driver) in models
            obstacleDict[level][id] = zeros(h,5) # x,y,θ,v,δ
        end
    end

    return obstacleDict
end
function updateObstacleMap!(obstacleMap, level, car_ID, trajectory)
    obstacleMap[level][car_ID][1:size(trajectory,1),1:size(trajectory,2)] = trajectory
end
function generateMotionMap(v_min,v_max,a_step,δ_max,δ_step,wheel_base,h,μ,Δt;g=9.81,a_range=[-1,0,1])

    v_step = a_step*Δt
    v_range = linspace(v_min,v_max,round(v_max/v_step))
    δ_range = linspace(-δ_max,δ_max,2*round((δ_max)/δ_step)+1)
    δ_range = (1/maximum(δ_range))*δ_range.*abs.(δ_range) # concentrate toward the middle
    R_range = wheel_base./tan.(δ_range)
    # a_range = [-1,0,1]
    motion_map = Dict()
    v_map = Dict()
    δ_map = Dict()

    for v_idx in 1:length(v_range)
        v = v_range[v_idx]
        v_map[v_idx] = zeros(Int,length(a_range),length(δ_range))
        δ_map[v_idx] = zeros(Int,length(a_range),length(δ_range))
        for a_idx in 1:length(a_range)
            Δv_idx = a_range[a_idx]
            tire_force = sqrt(((v^2)./abs.(R_range)).^2 + (Δv_idx*a_step)^2)
            if v_idx + Δv_idx >= 1 && v_idx + Δv_idx <= length(v_range)
                v_map[v_idx][a_idx,:] = Int(v_idx + Δv_idx)# v_range[v_idx+Δv_idx]
            else
                v_map[v_idx][a_idx,:] = Int(v_idx) #v_range[v_idx]
            end
            if v_idx > 1
                for j in 1:length(δ_range)
                    if tire_force[j] >= μ*g
                        v_map[v_idx][a_idx,j] = Int(v_idx) - 1# v_range[v_idx-1] # decelerate
                    end
                end
            end
            #δ_map[v_idx][a_idx,:] = δ_range.*(tire_force .< μ*g) +
             #   maximum(abs.(δ_range[tire_force .< μ*g])).*sign(δ_range).*(tire_force .>= μ*g)
            δ_map[v_idx][a_idx,:] = [i for i in 1:length(δ_range)].*(tire_force .< μ*g) +
            indmax(abs.(δ_range).*(tire_force .< μ*g)).*(tire_force .>= μ*g).*-sign(δ_range) +
            (length(δ_range)+1).*(sign(δ_range) .> 0).*(tire_force .>= μ*g)
        end
    end

    # fill motion map
    for v_idx in 1:length(v_range)
        motion_map[v_idx] = zeros(length(a_range),length(δ_range),h,3)
        level_v_idx = copy(v_map[v_idx])
        level_δ_idx = copy(δ_map[v_idx])
        Δθ = zeros(length(a_range),length(δ_range))
        for i in 1:h
            radius = wheel_base./tan.(δ_range[level_δ_idx])
            Δs = ((v_range[v_idx] + v_range[level_v_idx])./2.0)*Δt
            Δθ = Δs./radius
            ΔX = abs.(radius) .* sin.(abs.(Δθ))
            ΔX[:,Int((length(δ_range)-1)/2)+1] = Δs[:,Int((length(δ_range)-1)/2)+1]
            ΔY = radius.*(1 - cos.(Δθ))
            ΔY[:,Int((length(δ_range)-1)/2)+1] = 0
            if i == 1
                motion_map[v_idx][:,:,i,1] = ΔX
                motion_map[v_idx][:,:,i,2] = ΔY
                motion_map[v_idx][:,:,i,3] = Δθ
            else
                motion_map[v_idx][:,:,i,1] = motion_map[v_idx][:,:,i-1,1] +
                    ΔX.*cos.(motion_map[v_idx][:,:,i-1,3]) - ΔY.*sin.(motion_map[v_idx][:,:,i-1,3])
                motion_map[v_idx][:,:,i,2] = motion_map[v_idx][:,:,i-1,2] +
                    ΔX.*sin.(motion_map[v_idx][:,:,i-1,3]) + ΔY.*cos.(motion_map[v_idx][:,:,i-1,3])
                motion_map[v_idx][:,:,i,3] = Δθ + motion_map[v_idx][:,:,i-1,3]  # + motion_map[v_idx][:,:,i-1,3]
            end
            for j in 1:length(level_v_idx) # update to next velocity
                level_v_idx[j] = v_map[level_v_idx[j]][j]
                level_δ_idx[j] = δ_map[level_v_idx[j]][j]
            end
        end
    end
    return motion_map, v_map, δ_map, v_range, δ_range
end
function loopProjectionKD(hrhc,scene,roadway,tree)
    """
    projects all points in hrhc.successor_states to the kdtree representing
    the spline points along the centerline of roadway
    """
    curve = roadway.segments[1].lanes[1].curve

    s_grid = zeros(size(hrhc.successor_states,1),size(hrhc.successor_states,2))
    t_grid = zeros(size(s_grid))
    ϕ_grid = zeros(size(s_grid))
    idx_grid = zeros(Int,size(s_grid))

    pts = [reshape(hrhc.successor_states[:,:,1],size(hrhc.successor_states[:,:,1],1)*size(hrhc.successor_states[:,:,1],2),1)';
    reshape(hrhc.successor_states[:,:,2],size(hrhc.successor_states[:,:,2],1)*size(hrhc.successor_states[:,:,2],2),1)']
    idxs_list, _ = knn(tree,pts,1)
    idxs=reshape(idxs_list,size(hrhc.successor_states[:,:,2],1),size(hrhc.successor_states[:,:,2],2))

    for i in 1:size(s_grid,1)
        for j in 1:size(s_grid,2)
            idxA = idxs[i,j][1]-1
            idxB = idxs[i,j][1]+1
            if idxs[i,j][1] == length(curve)
                idxB = 1 # wrap to the beginning of the curve
            end
            if idxs[i,j][1] == 1
                idxA = length(curve) # wrap to the end of the curve
            end
            x = hrhc.successor_states[i,j,1]
            y = hrhc.successor_states[i,j,2]
            dA = sqrt(sum(([curve[idxA].pos.x, curve[idxA].pos.y]-[x,y]).^2))
            dB = sqrt(sum(([curve[idxB].pos.x, curve[idxB].pos.y]-[x,y]).^2))
            if dA < dB
                idxB = idxs[i,j][1]
            else
                idxA = idxs[i,j][1]
            end

            # project
            vec1 = [curve[idxB].pos.x - curve[idxA].pos.x, curve[idxB].pos.y - curve[idxA].pos.y, 0]
            vec2 = [x - curve[idxA].pos.x, y - curve[idxA].pos.y, 0]
            idx_t = dot(vec2, vec1)/norm(vec1)^2

            s_θ = curve[idxA].pos.θ + idx_t*(curve[idxB].pos.θ - curve[idxA].pos.θ)

            s_grid[i,j] = curve[idxA].s + idx_t*hrhc.Δs
            t_grid[i,j] = norm(vec2 - idx_t*vec1)*sign(sum(cross(vec1, vec2)))
            ϕ_grid[i,j] = wrap_to_π(hrhc.successor_states[i,j,3] - s_θ)
            idx_grid[i,j] = idxA
        end
    end
    # account for wrap-around
    s_grid[s_grid .< scene[hrhc.car_ID].state.posF.s] += curve[end].s + hrhc.Δs

    return s_grid, t_grid, ϕ_grid
end
function screenTrajectory(trajectory, obstacleMap, scene, roadway, hrhc, tree, k_level)
    out_of_bounds = false
    collision_flag = false
    # check out of bounds
    for i in 1 : size(trajectory,1)
        x = trajectory[i,1]
        y = trajectory[i,2]
        θ = trajectory[i,3]
        s,t,ϕ = kdProject(x,y,θ,tree,roadway,hrhc)

        if abs.(t) > hrhc.T_MAX
            out_of_bounds=true
            return out_of_bounds
        end
    end

    # collisionFlag = zeros(size(trajectory,1),1) # stores locations collisions
    threshold_dist = hrhc.car_length*4 # must be at least this close before we care to calculate collision cost

    # if k_level >= 1
    #     # for all other cars...
    #     for (id,traj) in obstacleMap[k_level - 1]
    #         if id != hrhc.car_ID
    #             state = scene[hrhc.car_ID].state
    #             state2 = scene[id].state
    #             diff = state.posG - state2.posG
    #             s1,_,_ = kdProject(state.posG.x,state.posG.y,state.posG.θ,tree,roadway,hrhc)
    #             s2,_,_ = kdProject(state2.posG.x,state2.posG.y,state2.posG.θ,tree,roadway,hrhc)
    #             if (norm([diff.x, diff.y]) < threshold_dist) && (s1 <= s2) # don't care if opponent is behind us
    #                 # R_idx = zeros(Int, size(trajectory,1),1) # to store the sorted indices of R
    #                 # ΔX = zeros(size(trajectory,1),1) # Δx, with opponent at origin
    #                 # ΔY = zeros(size(trajectory,1),1) # Δx, with opponent at origin
    #                 # Δθ = zeros(size(trajectory,1),1) # Δx, with opponent at origin
    #                 # for i in 1:hrhc.h
    #                     # pos = VecSE2(traj[i,1:3]) # x,y,θ of opponent at time step h
    #                 ΔX = trajectory[:,1] - traj[:,1] # Δx, with opponent at origin
    #                 ΔY = trajectory[:,2] - traj[:,2] # Δy with opponent at origin
    #                 Δθ = trajectory[:,3] - traj[:,3]  # Δθ with opponent at origin
    #                 # end
    #                 R = sqrt(ΔX.^2 + ΔY.^2)
    #                 R_idx = sortperm(R) # ordered by increasing distance
    #                 for idx in R_idx
    #                     # if collisionFlag[idx] == 1
    #                     #     continue
    #                     # end
    #                     if R[idx] > hrhc.car_length # no collision, and all other R values are greater
    #                         break
    #                     end
    #                     # R is less than hrhc.car_length
    #                     ψ = atan2(ΔY[idx],ΔX[idx]) - Δθ[idx]
    #                     r = W*L/sqrt((L*sin.(ψ))^2 + (W*cos.(ψ))^2) # radius of ellipse at given angle
    #                     if R[idx]-r < W*L/8 # collision
    #                         collision_flag = true
    #                         return collision_flag
    #                     end
    #                 end
    #             end
    #         end
    #     end
    # end

    return out_of_bounds || collision_flag
end
function getSuccessorStates(ΔXYθ, car_ID, h, scene::Scene)
    """ gets legal successor_states from motion primitives library """
    pos = scene[car_ID].state.posG # global x,y,z of car

    ΔX = ΔXYθ[:,:,h,1] * cos.(pos.θ) + ΔXYθ[:,:,h,2] * -sin.(pos.θ)
    ΔY = ΔXYθ[:,:,h,1] * sin.(pos.θ) + ΔXYθ[:,:,h,2] * cos.(pos.θ)
    Δθ = ΔXYθ[:,:,h,3]

    successor_states = zeros(size(ΔXYθ[:,:,h,:]))
    successor_states[:,:,1] = ΔX + pos.x
    successor_states[:,:,2] = ΔY + pos.y
    successor_states[:,:,3] = Δθ + pos.θ

    return successor_states
end
function computeTrajectory(ΔXYθ, car_ID, scene, v_cmd, δ_cmd, h)
    pos = scene[car_ID].state.posG

    trajectory = zeros(h,3)
    trajectory[:,1] = pos.x + ΔXYθ[v_cmd,δ_cmd,1:h,1]*cos.(pos.θ) + ΔXYθ[v_cmd,δ_cmd,1:h,2]*-sin.(pos.θ)
    trajectory[:,2] = pos.y + ΔXYθ[v_cmd,δ_cmd,1:h,1]*sin.(pos.θ) + ΔXYθ[v_cmd,δ_cmd,1:h,2]*cos.(pos.θ)
    trajectory[:,3] = pos.θ + ΔXYθ[v_cmd,δ_cmd,1:h,3]

    return trajectory
end
function checkCollisionElliptic(R,Δx,Δy,Δθ,L,W)
    b = W # buffer
    ζ = atan2(Δy,Δx)
    ψ = ζ - Δθ
    r1 = W*L/sqrt((L*sin.(ζ))^2 + (W*cos.(ζ))^2) # radius of ellipse at given angle
    r2 = W*L/sqrt((L*sin.(ψ))^2 + (W*cos.(ψ))^2)
    if R < r1 + r2 + b
        return true
    else
        return false
    end
end
function screenCollision(hrhc, obstacleMap, tree, roadway, scene, k_level)
    L = hrhc.ellipseB+2
    W = hrhc.ellipseA+2

    collisionFlag = zeros(size(hrhc.successor_states[:,:,1])) # stores locations collisions
    threshold_dist = hrhc.car_length*4 # must be at least this close before we care to calculate collision cost

    if k_level >= 1
        # for all other cars...
        for (id,trajectory) in obstacleMap[k_level - 1]
            if id != hrhc.car_ID
                state = scene[hrhc.car_ID].state
                state2 = scene[id].state
                diff = state.posG - state2.posG
                s1,_,_ = kdProject(state.posG.x,state.posG.y,state.posG.θ,tree,roadway,hrhc)
                s2,_,_ = kdProject(state2.posG.x,state2.posG.y,state2.posG.θ,tree,roadway,hrhc)
                if (norm([diff.x, diff.y]) < threshold_dist) && (s1 <= s2) # don't care if opponent is behind us
                    R_idx = zeros(Int, size(hrhc.successor_states,1)*size(hrhc.successor_states,2)) # to store the sorted indices of R
                    for i in 1:hrhc.h
                        pos = VecSE2(trajectory[i,1:3]) # x,y,θ of opponent at time step h
                        successor_states = getSuccessorStates(hrhc.motion_map[hrhc.v_cmd],hrhc.car_ID,i,scene)
                        ΔX = successor_states[:,:,1] - pos.x # Δx, with opponent at origin
                        ΔY = successor_states[:,:,2] - pos.y # Δy with opponent at origin
                        Δθ = successor_states[:,:,3] - pos.θ # Δθ with opponent at origin
                        R = sqrt(ΔX.^2 + ΔY.^2)
                        R_idx = sortperm(reshape(R,size(R,1)*size(R,2),1)[:])
                        for idx in R_idx
                            if collisionFlag[idx] == 1
                                continue
                            end
                            if R[idx] > hrhc.car_length # no collision, and all other R values are greater
                                break
                            end
                            # R is less than hrhc.car_length
                            if checkCollisionElliptic(R[idx],ΔX[idx],ΔY[idx],Δθ[idx],L,W)
                                collisionFlag[idx] = 1
                            end
                            # ψ = atan2(ΔY[idx],ΔX[idx]) - Δθ[idx]
                            # r = W*L/sqrt((L*sin.(ψ))^2 + (W*cos.(ψ))^2) # radius of ellipse at given angle
                            # if R[idx]-r < W*L/8 # collision
                            #     collisionFlag[idx] = 1
                            # end
                        end
                    end
                end
            end
        end
    end
    return collisionFlag
end
function tailgateAvoidance(hrhc, obstacleMap, tree, roadway, scene, k_level)
    L = 2*hrhc.car_length
    W = hrhc.car_width

    collisionFlag = zeros(size(hrhc.successor_states[:,:,1])) # stores locations collisions
    threshold_dist = hrhc.car_length*4 # must be at least this close before we care to calculate collision cost

    R = zeros(size(hrhc.successor_states[:,:,1]))
    ψ = zeros(size(hrhc.successor_states[:,:,1]))
    cost = zeros(size(hrhc.successor_states[:,:,1]))

    if k_level >= 1
        # for all other cars...
        for (id,trajectory) in obstacleMap[k_level - 1]
            if id != hrhc.car_ID && id == 3
                state = scene[hrhc.car_ID].state
                state2 = scene[id].state
                diff = state.posG - state2.posG
                s1,_,_ = kdProject(state.posG.x,state.posG.y,state.posG.θ,tree,roadway,hrhc)
                s2,_,_ = kdProject(state2.posG.x,state2.posG.y,state2.posG.θ,tree,roadway,hrhc)
                if (norm([diff.x, diff.y]) < threshold_dist) && (s1 <= s2) # don't care if opponent is behind us
                    pos = VecSE2(trajectory[hrhc.h,1:3]) # x,y,θ of opponent at time step h
                    successor_states = getSuccessorStates(hrhc.motion_map[hrhc.v_cmd],hrhc.car_ID,hrhc.h,scene)
                    ΔX = successor_states[:,:,1] - pos.x # Δx, with opponent at origin
                    ΔY = successor_states[:,:,2] - pos.y # Δy with opponent at origin
                    Δθ = successor_states[:,:,3] - pos.θ # Δθ with opponent at origin
                    R = sqrt(ΔX.^2 + ΔY.^2)
                    ψ = atan2(ΔY,ΔX) - Δθ

                    cost = (1./(W*R.*cos.(ψ).^2 + 2*L*R.*sin.(ψ).^2 + 1)).*(-cos.(ψ).^3 + 1.1) + 1

                end
            end
        end
    end
    return cost
end
function calculateObjective(car_ID,s₀,s,t,ϕ,T_MAX;ϕ_MAX=Float64(π),s_factor=1.0)
    s_norm = (s-s₀)/maximum(s)
    t_norm = t/T_MAX
    ϕ_norm = ϕ/ϕ_MAX

    #costs
    t_cost = abs.(t_norm).^6
    ϕ_cost = abs.(ϕ_norm).^6
    s_factor = 1
    s_cost = s_factor*(1-s_norm)
    A = [1 .5; #  [ϕ t] [a1 a2] [ϕ]
        .5 0] #         [a2 a3] [t]
    tϕ_cost = A[1,1]*(ϕ_norm).^2 + (A[1,2]+A[2,1])*(ϕ_norm).*(t_norm) + A[2,2]*(t_norm).^2

    objective = 1+s_cost+t_cost+ϕ_cost+tϕ_cost
    return objective
end
function AutomotiveDrivingModels.observe!(hrhc::HRHC, scene::Scene, roadway::Roadway, egoid::Int, tree::KDTree, obstacleMap, k_level)
    """
    Observe the current environment and select optimal action to apply at next
    time step
    """

    state = scene[hrhc.car_ID].state
    hrhc.curve_ind = state.posF.roadind.ind.i
    v = state.v # current v
    hrhc.v = v
    trajectory = zeros(hrhc.h,3)
    action_selected = false
    δ_cmd = 1
    a_cmd = 1

    i = 0
    for i in 0:(hrhc.h-1)
        if action_selected
            break # out of for loop
        end
        # calculate successor states
        hrhc.successor_states = getSuccessorStates(hrhc.motion_map[hrhc.v_cmd], hrhc.car_ID, hrhc.h-i, scene)
        # project successor states onto track
        s,t,ϕ = loopProjectionKD(hrhc, scene, roadway, tree)
        # optimization objective
        objective = calculateObjective(hrhc.car_ID, scene[hrhc.car_ID].state.posF.s,
        s, t, ϕ, hrhc.T_MAX)
        tailgateCost = tailgateAvoidance(hrhc, obstacleMap, tree, roadway, scene, k_level)
        objective = objective + tailgateCost
        collisionFlag = screenCollision(hrhc, obstacleMap, tree, roadway, scene, k_level)
        objective[collisionFlag .> 0] = Inf

        while (action_selected==false) && (minimum(objective) != Inf)
            index = indmin(objective) # find get a better method of optimizing this
            a_cmd, δ_cmd = ind2sub(s, index)
            hrhc.δ_cmd = δ_cmd
            # compute full trajectory up to horizon
            trajectory = computeTrajectory(hrhc.motion_map[hrhc.v_cmd], hrhc.car_ID, scene, a_cmd, δ_cmd, hrhc.h-i)
            # screen trajectory for collisions / validity
            out_of_bounds = screenTrajectory(trajectory, obstacleMap, scene, roadway, hrhc, tree, k_level)
            if out_of_bounds
                objective[index] = Inf
            else
                action_selected=true
                # updateObstacleMap!(obstacleMap, k_level, hrhc.car_ID, trajectory)
                if k_level != hrhc.k # only assign commands for the given k_level
                    return # do not update the action
                end
            end
        end
    end

    hrhc.δ_cmd = δ_cmd
    # hrhc.v_cmd = max(hrhc.v_cmd + a_cmd - 2,1) # assumes 3 options for acceleration
    hrhc.v_cmd = hrhc.v_map[hrhc.v_cmd][a_cmd,δ_cmd]

    hrhc.δ = hrhc.δ_range[hrhc.δ_cmd] # next δ
    hrhc.v = hrhc.v_range[hrhc.v_cmd] # next v

    next_state = VehicleState(VecSE2(trajectory[1,1:3]),roadway,hrhc.v)
    hrhc.action = next_state # action
end
AutomotiveDrivingModels.get_name(::HRHC) = "HRHC"
# AutomotiveDrivingModels.action_context(driver::HRHC) = driver.action_context # AutomotiveDrivingModels.action_context
Base.rand(hrhc::HRHC) = hrhc.action

# Plotting functions
function plotSplineRoadway(x,y,θ,lane_width)
    perp_lines1 = zeros(2,length(x))
    perp_lines2 = zeros(2,length(x))

    perp_lines1[1,:] = x + (lane_width/2.0)*sin.(θ)
    perp_lines1[2,:] = y - (lane_width/2.0)*cos.(θ)
    perp_lines2[1,:] = x - (lane_width/2.0)*sin.(θ)
    perp_lines2[2,:] = y + (lane_width/2.0)*cos.(θ)

    # PyPlot.figure()
    # PyPlot.scatter(x,y)
    PyPlot.plot(x,y)
    PyPlot.plot(perp_lines1[1,:],perp_lines1[2,:],color="green")
    PyPlot.plot(perp_lines2[1,:],perp_lines2[2,:],color="green")
    PyPlot.axis("equal")
    # PyPlot.show()
end
function plotObjectiveHorizon(hrhc,scene,roadway,tree,trajectory,obstacleMap,xR,yR,θR)
    lo=hrhc.curve_ind
    if :V_MAX in fieldnames(hrhc)
        hi = hrhc.curve_ind + Int(1+div(hrhc.V_MAX*hrhc.Δt*hrhc.h,hrhc.Δs))
    else
        hi = hrhc.curve_ind + Int(1+div(hrhc.v_range[end]*hrhc.Δt*hrhc.h,hrhc.Δs))
    end
    lane_width = roadway.segments[1].lanes[1].width

    x = zeros(hrhc.h,size(hrhc.successor_states,1),size(hrhc.successor_states,2))
    y = zeros(size(x))
    Θ = zeros(size(x))
    s = zeros(size(x))
    t = zeros(size(x))
    ϕ = zeros(size(x))
    objective = zeros(size(x))

    for i in 1:hrhc.h
        getLegalMoves!(hrhc, scene, h=i)
        getSuccessorStates!(hrhc, scene)
        x[i,:,:] = copy(hrhc.successor_states[:,:,1])
        y[i,:,:] = copy(hrhc.successor_states[:,:,2])
        Θ[i,:,:] = copy(hrhc.successor_states[:,:,3])
        s[i,:,:], t[i,:,:], ϕ[i,:,:] = loopProjectionKD(hrhc,scene,roadway,tree)
        objective[i,:,:] = calculateObjective(hrhc,scene, roadway, tree,s[i,:,:],t[i,:,:],ϕ[i,:,:],obstacleMap,hrhc.k,hrhc.h)
    end

    PyPlot.figure(figsize=[12,4])

    PyPlot.subplot(141) # ϕ
    plotSplineRoadway(xR[lo:hi],yR[lo:hi],θR[lo:hi],lane_width)
    PyPlot.scatter(x,y,c=ϕ,edgecolor="none")
    PyPlot.plot(trajectory[:,1],trajectory[:,2],color="red")
    PyPlot.axis("off")
    PyPlot.title("|phi|")

    PyPlot.subplot(142) # s
    plotSplineRoadway(xR[lo:hi],yR[lo:hi],θR[lo:hi],lane_width)
    PyPlot.scatter(x,y,c=s,edgecolor="none")
    PyPlot.plot(trajectory[:,1],trajectory[:,2],color="red")
    PyPlot.axis("off")
    PyPlot.title("s")

    PyPlot.subplot(143) # t
    plotSplineRoadway(xR[lo:hi],yR[lo:hi],θR[lo:hi],lane_width)
    PyPlot.scatter(x,y,c=t,edgecolor="none")
    PyPlot.plot(trajectory[:,1],trajectory[:,2],color="red")
    PyPlot.axis("off")
    PyPlot.title("t")

    PyPlot.subplot(144) # objective
    plotSplineRoadway(xR[lo:hi],yR[lo:hi],θR[lo:hi],lane_width)
    PyPlot.scatter(x,y,c=log(objective),edgecolor="none")
    PyPlot.plot(trajectory[:,1],trajectory[:,2],color="red")
    PyPlot.axis("off")
    PyPlot.title("log objective")
end
function plot_stϕ(hrhc,roadway,scene,x,y,θ,trajectory,s,t,ϕ,objective)
    lo=hrhc.curve_ind
    hi=hrhc.curve_ind + Int(1+2*div(hrhc.V_MAX*hrhc.Δt*hrhc.h,hrhc.Δs))
    if hi > length(roadway.segments[1].lanes[1].curve)
        lo = length(roadway.segments[1].lanes[1].curve)
        hi=hrhc.curve_ind + Int(1+2*div(hrhc.V_MAX*hrhc.Δt*hrhc.h,hrhc.Δs))
    end
    lane_width = roadway.segments[1].lanes[1].width

    PyPlot.figure(figsize=[12,4])
    PyPlot.subplot(141)
    plotSplineRoadway(x[lo:hi],y[lo:hi],θ[lo:hi],lane_width)
    # PyPlot.scatter(Pts[1,:],Pts[2,:],color="red")
    PyPlot.scatter(hrhc.successor_states[:,:,1],hrhc.successor_states[:,:,2],c=abs.(ϕ),edgecolor="none")
    PyPlot.plot(trajectory[:,1],trajectory[:,2],color="red")
    PyPlot.scatter(scene[hrhc.car_ID].state.posG.x, scene[hrhc.car_ID].state.posG.y, c="k", edgecolors="none",s=40)
    PyPlot.axis("off")
    PyPlot.title("|phi|")

    PyPlot.subplot(142)
    plotSplineRoadway(x[lo:hi],y[lo:hi],θ[lo:hi],lane_width)
    # PyPlot.scatter(Pts[1,:],Pts[2,:],color="red")
    PyPlot.scatter(hrhc.successor_states[:,:,1],hrhc.successor_states[:,:,2],c=abs.(t),edgecolor="none")
    PyPlot.plot(trajectory[:,1],trajectory[:,2],color="red")
    PyPlot.scatter(scene[hrhc.car_ID].state.posG.x, scene[hrhc.car_ID].state.posG.y, c="k", edgecolors="none",s=40)
    PyPlot.axis("off")
    PyPlot.title("|t|")

    PyPlot.subplot(143)
    plotSplineRoadway(x[lo:hi],y[lo:hi],θ[lo:hi],lane_width)
    # PyPlot.scatter(Pts[1,:],Pts[2,:],color="red")
    PyPlot.scatter(hrhc.successor_states[:,:,1],hrhc.successor_states[:,:,2],c=s,edgecolor="none")
    PyPlot.scatter(scene[hrhc.car_ID].state.posG.x, scene[hrhc.car_ID].state.posG.y, c="k", edgecolors="none",s=40)
    PyPlot.plot(trajectory[:,1],trajectory[:,2],color="red")
    PyPlot.axis("off")
    PyPlot.title("s")

    PyPlot.subplot(144)
    plotSplineRoadway(x[lo:hi],y[lo:hi],θ[lo:hi],lane_width)
    # PyPlot.scatter(Pts[1,:],Pts[2,:],color="red")
    PyPlot.scatter(hrhc.successor_states[:,:,1],hrhc.successor_states[:,:,2],c=log(objective),edgecolor="none")
    PyPlot.plot(trajectory[:,1],trajectory[:,2],color="red")
    PyPlot.scatter(scene[hrhc.car_ID].state.posG.x, scene[hrhc.car_ID].state.posG.y, c="k", edgecolors="none",s=40)
    PyPlot.axis("off")
    PyPlot.title("objective")
end
function plotHRHCInfo(hrhc,models,scene,roadway,trajectory,cmd,x,y,Θ,s,t,ϕ,objective)
    lo = hrhc.curve_ind
    hi = hrhc.curve_ind + Int(1+2*div(hrhc.V_MAX*hrhc.Δt*hrhc.h,hrhc.Δs))
    lane_width = roadway.segments[1].lanes[1].width
    if hi > length(roadway.segments[1].lanes[1].curve)
        lo = length(roadway.segments[1].lanes[1].curve)
        hi=hrhc.curve_ind + Int(1+2*div(hrhc.V_MAX*hrhc.Δt*hrhc.h,hrhc.Δs))
    end
    PyPlot.figure(figsize=[12,10])
    # Plot Raceway
    PyPlot.subplot(221)
    # plotSplineRoadway(x[lo:hi],y[lo:hi],θ[lo:hi],lane_width)
    plotSplineRoadway(x,y,Θ,lane_width)
    PyPlot.scatter(hrhc.successor_states[:,:,1],hrhc.successor_states[:,:,2],color="red")
    PyPlot.plot(trajectory[:,1],trajectory[:,2],color="red")
    PyPlot.scatter(roadway.segments[1].lanes[1].curve[hrhc.curve_ind].pos.x, roadway.segments[1].lanes[1].curve[hrhc.curve_ind].pos.y, c="white", s=40)
    for (id,car) in models
        if id == hrhc.car_ID
            PyPlot.scatter(scene[id].state.posG.x,scene[id].state.posG.y,c="red",s=20)
        else
            PyPlot.scatter(scene[id].state.posG.x,scene[id].state.posG.y,c="blue",s=20)
        end
    end
    PyPlot.axis("off")
    PyPlot.title("Raceway with Motion Primitives")

    PyPlot.subplot(222)
    plotSplineRoadway(x[lo:hi],y[lo:hi],Θ[lo:hi],lane_width)
    PyPlot.scatter(scene[hrhc.car_ID].state.posG.x, scene[hrhc.car_ID].state.posG.y, c="red", edgecolors="none",s=40)
    PyPlot.scatter(hrhc.successor_states[:,:,1], hrhc.successor_states[:,:,2],c=log(objective),edgecolors="none")
    PyPlot.scatter(hrhc.successor_states[cmd[1],cmd[2],1], hrhc.successor_states[cmd[1],cmd[2],2],c="white",s=40)
    PyPlot.plot(trajectory[:,1],trajectory[:,2],color="red")
    PyPlot.axis("off")
    PyPlot.title("Log Objective Function")
end


text:
#######################
## outer contructors ##
#######################

## concept:
## If no NA values are present, initialization could be done more
## easily with arrays.

FloatArray = Union(Array{Float64, 2}, Array{Float64, 1})
import DataFrames.convert
function convert(T::Type{DataFrame}, vals::Array{Float64, 1})
    df = convert(T, (vals')')
end

for t = (:Timedata, :Timenum, :Timematr, :Timecop)
    @eval begin

        ###########################
        ## shortcuts without NAs ##
        ###########################

        ## no names or dates (just simulated values)
        function $(t)(vals::FloatArray)
            idx = [1:size(vals, 1)]
            $(t)(convert(DataFrame, vals), idx)
        end

        ## from core elements
        function $(t)(vals::FloatArray,
                      nams::Array{Symbol, 1},
                      idx::Array)
            df = composeDataFrame(vals, nams)
            return $(t)(df, idx)
        end

        ## comprehensive constructor: very general, all elements
        ## required for Timedata type
        function $(t)(vals::Array, nams::Array, idx::Array)
            df = composeDataFrame(vals, nams)
            return $(t)(df, idx)
        end

        ## two inputs only, general form
        function $(t)(vals::Array, nams::Array)
            if isa(nams[1], Symbol)
                df = composeDataFrame(vals, nams)
                $(t)(df)
            else
                df = convert(DataFrame, vals)
                $(t)(df, nams)
            end
        end

        function $(t)(vals::DataFrame, nams::Array)
            if isa(nams[1], Symbol)
                df = composeDataFrame(vals, nams)
                $(t)(df)
            else
                df = convert(DataFrame, vals)
                $(t)(df, nams)
            end
        end

        ###################
        ## NAs in values ##
        ###################

        ## no idx
        function $(t)(vals::DataFrame)
            idx = [1:size(vals, 1)]
            tn = $(t)(vals, idx)
        end

        ########################
        ## index as DataArray ##
        ########################

        function $(t)(vals::DataFrame, idx::DataVector)
            if any(idx.na)
                error("Index values may not be missing")
            else
                return $(t)(vals, asArr(idx, eltype(idx)))
            end
        end

    end
end

ext: jl


text:
module BenchStep
using BenchmarkTools
using GrandGameOfLife
using Random

suite = BenchmarkGroup()

Random.seed!(2020)
state = create_initial_state(1000, 1000, 0.7)

suite["step"] = @benchmarkable makestep($state)

end # module

BenchDistance.suite


text:
using HCubature, StaticArrays
using Test

@testset "simple" begin
      @test hcubature(x -> cos(x[1])*cos(x[2]), [0,0], [1,1])[1] ≈ sin(1)^2 ≈
            @inferred(hcubature(x -> cos(x[1])*cos(x[2]), (0,0), (1,1)))[1]
      @test @inferred(hcubature(x -> cos(x[1]), (0,), (1,)))[1] ≈ sin(1) ≈
            @inferred(hquadrature(cos, 0, 1))[1]
      @test @inferred(hcubature(x -> cos(x[1]), (0.0f0,), (1.0f0,)))[1] ≈ sin(1.0f0)
      @test @inferred(hcubature(x -> 1.7, SVector{0,Float64}(), SVector{0,Float64}()))[1] == 1.7
end

# function wrapper for counting evaluations
const gcnt = Ref(0)
cnt(f) = x -> begin
      gcnt[] += 1
      f(x)
end

@testset "big3d" begin
      let f(x) = sin(x[1] + 3*sin(2*x[2] + 4*sin(3*x[3])))
            gcnt[] = 0
            @test hcubature(cnt(f), (0,0,0),(3,3,3), rtol=1e-6)[1] ≈ -4.78802790509727 rtol=1e-6
            @test 2_400_000 < gcnt[] < 2_500_000
            gcnt[] = 0
            hcubature(cnt(f), (0,0,0),(3,3,3), maxevals=1000)
            @test gcnt[] < 1070
      end
end

@testset "complex" begin
      @test hcubature(x -> cis(x[1]), (0,), (1,))[1] ≈ (cis(1)-1)/im
      @test hcubature(x -> cis(x[1]+x[2]), (0,-1), (1,2))[1] ≈ -(cis(1)-1)*(cis(2)-cis(-1))
end

@testset "countevals" begin
      let g = HCubature.GaussKronrod(Float64)
            @test HCubature.countevals(g) == 1 + 2length(g.w)
      end
      for n = 2:10
            let g = HCubature.GenzMalik(Val{n}(), Float64)
                  @test HCubature.countevals(g) == 1 + 4length(g.p[1]) + length(g.p[3]) + length(g.p[4])
            end
      end
      @test HCubature.countevals(HCubature.Trivial()) == 1
end

@testset "axischoosing" begin
    #Issue #4
    let f = x -> 1.0 + (x[1] * x[3] * sin(x[2]))^2
        @test hcubature(f, (0.0,0.0,-0.2), (0.2,2π,0.2), rtol=1e-6)[1] ≈ (22502//140625)*π rtol=1e-6
    end
end

@testset "initdiv" begin
      for initdiv = 1:5
            @test sin(1)^2 ≈ hcubature(x -> cos(x[1])*cos(x[2]), [0,0], [1,1], initdiv=initdiv)[1]
            @test sin(1) ≈ @inferred(hquadrature(cos, 0, 1, initdiv=initdiv))[1]
      end
end


text:
function ff_data(
    conn::Union{LibPQ.Connection, DBInterface.Connection},
    date_start::Date=Date(1926, 7, 1),
    date_end::Date=today();
    cols::Array{String}=[
        "date",
        "mktrf",
        "smb",
        "hml",
        "rf",
        "umd"
    ]
)
    col_str = join(cols, ", ")
    query = """
        SELECT $col_str FROM $(default_tables["ff_factors"])
        WHERE date BETWEEN '$date_start' AND '$date_end'
    """
    return run_sql_query(conn, query)
end


text:
abstract AbstractPlanOutMapper

type Assignment <: AbstractPlanOutMapper
    salt::AbstractString
    overrides::Dict
    data::Dict
end

Assignment(salt::AbstractString) = Assignment(salt, Dict{AbstractString, Any}())
Assignment(salt::AbstractString, overrides::Dict) = Assignment(salt, copy(overrides), copy(overrides))

function setindex!(a::Assignment, value, name)
    setindex!(a.data, getValue(a, value, name), name)
end

function setsalt!(a::Assignment, val)
    a.salt = val
end

function getindex(a::Assignment, name)
    getindex(a.data, name)
end

function delete!(a::Assignment, name)
    delete!(a.data, name)
end

evaluate(a::Assignment, value) = value
start(a::Assignment) = start(a.data)
next(a::Assignment, i) = next(a.data, i)
done(a::Assignment, i) = done(a.data, i)

length(a::Assignment) = length(a.data)


text:
using Statistics
using Random
using StatsBase
using LinearAlgebra

struct ConstellationHistogram{T1,T2,T3,T4,T5}
    nbins::Int
    dopplers::Array{Float64,1}
    doppler_rates::Array{Float64,1}
    histogram::T1
    weights::T2
    fd_rate::Vector{Float64}
    f_d::Vector{Float64}
    elevations::T3
    elev_hists::T4
    max_elevations::T5
end

function get_histogram(x; nbins)
    return normalize(fit(Histogram, x, nbins=nbins), mode=:probability)
end

function initiate_hist(dopplers, doppler_rates; nbins=1000)
    histogram = normalize(fit(Histogram, (doppler_rates, dopplers),
                              nbins=nbins), mode=:probability)
    weights = collect(Iterators.flatten(histogram.weights))
    weights = Weights(weights)
    fd_rate, f_d = meshgrid(Array(histogram.edges[1][1:end-1]),
                            Array(histogram.edges[2][1:end-1]))
    fd_rate = collect(Iterators.flatten(fd_rate))
    f_d = collect(Iterators.flatten(f_d))
    return ConstellationHistogram(nbins, dopplers, doppler_rates,
                                  histogram, weights, fd_rate, f_d,
                                  missing, missing, missing)
end

function initiate_hist(dopplers, doppler_rates, elevations, max_elevations;
                       nbins=1000, elev_bins=100)
    histogram = normalize(fit(Histogram, (doppler_rates, dopplers),
                              nbins=nbins), mode=:probability)
    histogram3d = normalize(fit(Histogram, (elevations, doppler_rates, dopplers),
                               nbins=(elev_bins, nbins, nbins)), mode=:probability)
    histogram_dopplers =  normalize(fit(Histogram, (elevations, dopplers),
                                        nbins=(elev_bins, nbins)), mode=:probability)
    histogram_doppler_rates =  normalize(fit(Histogram, (elevations, doppler_rates),
                                        nbins=(elev_bins, nbins)), mode=:probability)
    elev_hists = [histogram3d, histogram_dopplers, histogram_doppler_rates]
    weights = collect(Iterators.flatten(histogram.weights))
    weights = Weights(weights)
    fd_rate, f_d = meshgrid(Array(histogram.edges[1][1:end-1]),
                            Array(histogram.edges[2][1:end-1]))
    fd_rate = collect(Iterators.flatten(fd_rate))
    f_d = collect(Iterators.flatten(f_d))
    return ConstellationHistogram(nbins, dopplers, doppler_rates,
                                  histogram, weights, fd_rate, f_d,
                                  elevations, elev_hists, max_elevations)
end

function sample_distribution(histogram::ConstellationHistogram)
    f_d = sample(histogram.f_d, histogram.weights)
    fd_rate = sample(histogram.fd_rate, histogram.weights)
    return (f_d, fd_rate)
end

function max_doppler_curve_per_elev(histogram::ConstellationHistogram, elevs)
    elevations = histogram.elevations

    doppler_bins = histogram.histogram.edges[1]
    dopplers = histogram.dopplers
    doppler_rates = histogram.doppler_rates
    N = size(histogram.elev_hists[2].weights)[1]
    doppler_medians = zeros(N)
    doppler_rate_medians = zeros(N)
    for i in 1:N
        doppler_median_idx = argmax(histogram.elev_hists[2].weights[i,:])
        doppler_rate_median_idx = argmax(histogram.elev_hists[3].weights[i,:])
        doppler_medians[i] = histogram.elev_hists[2].edges[2][doppler_median_idx]
        doppler_rate_medians[i] = histogram.elev_hists[3].edges[2][doppler_rate_median_idx]
    end
    return (doppler_medians, doppler_rate_medians)
end


text:
using JSON3
function process_components_meta(metafile)
    metadata = JSON3.read(
        read(metafile, String)
        )
    result = []
    for (file, props) in metadata
        name = split(split(string(file), "/")[end], ".")[1]
        push!(result, make_component_meta(name, props))
    end
    return result
end

function make_component_meta(name, props)
    args = filter(filter_arg, props["props"])
    regular_args = filter(collect(keys(args))) do name
        !endswith(string(name), "-*")
    end
    wild_args =
        filter(collect(keys(args))) do name
            endswith(string(name), "-*")
        end
    return OrderedDict(
        :name => Symbol(name),
        :args => regular_args,
        :wild_args => [Symbol(replace(string(a), "-*"=>"")) for a in wild_args],
        :docstr => docstring(name, args, props["description"]),
    )
end

const _reserved_words = Set(
    ["baremodule",
    "begin",
    "break",
    "catch",
    "const",
    "continue",
    "do",
    "else",
    "elseif",
    "end",
    "export",
    "false",
    "finally",
    "for",
    "function",
    "global",
    "if",
    "import",
    "let",
    "local",
    "macro",
    "module",
    "quote",
    "return",
    "struct",
    "true",
    "try",
    "using",
    "while"]
)
function is_reserved_world(w)
    return w in _reserved_words
end

function filter_arg(argpair)
    name, props = argpair
    is_reserved_world(name) && return false
    if haskey(props, "type")
        arg_type = props["type"]["name"]
        return !in(arg_type, ["func", "symbol", "instanceOf"])
    end
    if "flowType" in props
        arg_type_name = props["flowType"]["name"]
        if arg_type_name == "signature"
            # This does the same as the PropTypes filter above, but "func"
            # is under "type" if "name" is "signature" vs just in "name"
            if !in("type", props["FlowType"]) || props["FlowType"]["type"] != "object"
                return false
            end
        end
        return true
    end
    return false
end

function docstring(name, props, description)
    article = lowercase(first(name)) in ['a', 'e', 'i', 'o', 'u'] ? "An " : "A "
    result = string(
        article, name, " component", "\n",
        description, "\n\n"
    )
    if haskey(props, :children)
        result *= arg_docstring("children", props[:children]) * "\n"
    end
    if haskey(props, :id)
         result *= arg_docstring("id", props[:id]) * "\n"
    end
    other_props = sort(
        collect(filter(v->!in(v.first, [:children, :id]), props)),
        lt = (a, b) -> a[1] < b[1]
    )
    result *= join(arg_docstring.(other_props), "\n")
    return result

end

_jl_type(::Val{:array}, type_object) = "Array"
_jl_type(::Val{:bool}, type_object) = "Bool"
_jl_type(::Val{:string}, type_object) = "String"
_jl_type(::Val{:object}, type_object) = "Dict"
_jl_type(::Val{:any}, type_object) = "Bool | Real | String | Dict | Array"
_jl_type(::Val{:element}, type_object) = "dash component"
_jl_type(::Val{:node}, type_object) = "a list of or a singular dash component, string or number"

_jl_type(::Val{:enum}, type_object) = join(
        string.(
            getindex.(type_object["value"], :value)
        ), ", "
    )

function _jl_type(::Val{:union}, type_object)
    join(
        filter(a->!isempty(a), jl_type.(type_object["value"])),
        " | "
    )
end

function _jl_type(::Val{:arrayOf}, type_object)
    result = "Array"
    if type_object["value"] != ""
        result *= string(" of ", jl_type(type_object["value"]), "s")
    end
    return result
end

_jl_type(::Val{:objectOf}, type_object) =
    string("Dict with Strings as keys and values of type ", jl_type(type_object["value"]))

function _jl_type(::Val{:shape}, type_object)
    child_names = join(string.(keys(type_object["value"])), ", ")
    result = "lists containing elements $(child_names)"
    result *= join(
        [
            arg_docstring(name, prop, prop["required"], get(prop, "description", ""), 1)
            for (name, prop) in type_object["value"]
        ]
    )
    return result
end
_jl_type(::Val{:exact}, type_object) = _jl_type(Val(:shape), type_object)
_jl_type(val, type_object) = ""

jl_type(type_object) = _jl_type(Val(Symbol(type_object["name"])), type_object)

arg_docstring(name_prop::Pair, indent_num = 0) = arg_docstring(name_prop[1], name_prop[2], indent_num)
arg_docstring(name, prop, indent_num = 0) =
    arg_docstring(
        string(name),
        haskey(prop, "type") ? prop["type"] : prop["flowType"],
        prop["required"],
        prop["description"],
        indent_num
    )

function arg_docstring(prop_name, type_object, required, description, indent_num)
    typename = jl_type(type_object)
    indent_spacing = repeat("   ", indent_num)
    if occursin("\n", typename)
        return string(
            indent_spacing, "- `", prop_name, "` ",
            "(", required ? "required" : "optional", "):",
            description, ". ",
            prop_name, " has the following type: ", typename
        )
    end

    return string(
        indent_spacing, "- `", prop_name, "` ",
            "(",
                 isempty(typename) ? "" : string(typename, "; "),
                required ? "required" : "optional",
            ")",
            isempty(description) ? "" : string(": ", description)
    )
end


text:
include("../common/ortho_gram_schmidt.jl")
include("../common/wrappers_fft.jl")
include("../common/gen_lattice_pwscf.jl")

include("gen_dr.jl")
include("init_pot_harm_3d.jl")
include("op_K.jl")
include("op_Vpot.jl")
include("op_H.jl")
include("calc_rho.jl")
include("calc_grad.jl")
include("calc_Etot.jl")
include("schsolve_Emin_sd.jl")
include("schsolve_Emin_cg.jl")
include("Kprec.jl")
include("diag_lobpcg.jl")
include("diag_davidson.jl")

function test_main( ; Ns_in=nothing, ecutwfc_Ry=nothing,
                   solution_method="diag_lobpcg" )

    LatVecs = gen_lattice_sc(6.0)

    if PWGRID_VERSION == 2
        if Ns_in != nothing
            pw = PWGrid( Ns_in, LatVecs )
        else
            @printf("PWGRID_VERSION 2 needs Ns_in\n")
            exit
        end
    elseif PWGRID_VERSION == 3
        if ecutwfc_Ry != nothing
            pw = PWGrid( 0.5*ecutwfc_Ry, LatVecs )
        else
            @printf("PWGRID_VERSION 3 needs ecutwfc_Ry\n")
            exit
        end
    else
        @printf("ERROR: Must specify Ns or ecutwfc_Ry")
        exit()
    end

    Ω  = pw.Ω
    r  = pw.r
    G  = pw.gvec.G
    G2 = pw.gvec.G2
    Ns = pw.Ns
    Npoints = prod(Ns)
    Ngwx = pw.gvecw.Ngwx

    @printf("Ns   = (%d,%d,%d)\n", Ns[1], Ns[2], Ns[3])
    @printf("Ngwx = %d\n", Ngwx)
    @printf("G2mx = %f\n", maximum(pw.gvec.G2[pw.gvecw.idx_gw2r]))

    if PWGRID_VERSION == 3
        @printf("ecutwfc_Ry = %f\n", ecutwfc_Ry)
    end

    actual = Npoints/Ngwx
    theor = 1/(4*pi*0.25^3/3)
    @printf("Compression: actual, theor: %f , %f\n", actual, theor)

    # Generate array of distances
    center = 6.0*ones(3)/2
    dr = gen_dr( r, center )

    # Setup potential
    Vpot = init_pot_harm_3d( pw, dr )
    println("sum(Vpot)*Ω/Npoints = ", sum(Vpot)*Ω/Npoints)

    Nstates = 4
    srand(2222)
    psi  = rand(ComplexF64,Ngwx,Nstates)
    psi = ortho_gram_schmidt(psi)

    if solution_method == "Emin"

        psi, Etot = schsolve_Emin_sd( pw, Vpot, psi, NiterMax=10 )
        psi, Etot = schsolve_Emin_cg( pw, Vpot, psi, NiterMax=1000 )

        Y = ortho_gram_schmidt(psi)
        mu = Y' * op_H( pw, Vpot, Y )
        evals, evecs = eigen(mu)
        Psi = Y*evecs

    else

        evals, psi = diag_lobpcg( pw, Vpot, psi, verbose=true, tol_avg=1e-10 )

    end

    # Davidson diagonalization is not working yet
    #evals, psi = diag_davidson( pw, Vpot, psi, verbose=true, tol_avg=1e-10 )

    for ist = 1:Nstates
        @printf("State # %d, Energy = %18.10f\n", ist, real(evals[ist]))
    end

end


text:
import FourierTransforms
import FourierTransforms: fft_kernel_sizes
using Test, LinearAlgebra, OffsetArrays

include("symbolic.jl")

function testdft(x::AbstractVector{T}) where T<:Union{Real,Complex}
    Tr = real(float(T))
    Tc = complex(float(T))
    y = similar(x, Tc)
    N = length(x)
    for (k1, yind) in enumerate(eachindex(y))
        fct = zero(Tr)
        yk  = zero(Tc)
        for xi in x
            yk  += xi*complex(cos(fct), -sin(fct))
            fct += one(Tr)/N*(k1 - 1)*2*π
        end
        y[yind] = yk
    end
    return y
end

@testset "Multiplication of size: $sz" for sz in [collect(fft_kernel_sizes); 511; 512; 513]
  @testset "Real plan" begin
      xr = rand(sz)
      @test FourierTransforms.fft(xr) ≈ testdft(xr)
      p = FourierTransforms.plan_fft(xr)
      @test p*xr == FourierTransforms.fft(xr)
      @test p\(p*xr) ≈ xr
  end

  @testset "Complex plan applied to $atype" for (atype, x) in (
      ("Vector", rand(Complex{Float64}, sz)),
      ("OffsetVector", OffsetArray(rand(Complex{Float64}, sz), 0:(sz - 1)))
  )
      p = FourierTransforms.plan_fft(x)
      @test mul!(similar(x), p, x) == p*x
      @test p*x ≈ testdft(x)
      @test p\(p*x) ≈ x
  end
end


text:
# Autogenerated thin wrapper around CNOVAS using Clang.jl
module LibNOVAS

using NOVAS_jll
export NOVAS_jll

using CEnum

function solarsystem(tjd, body, origin, position, velocity)
    return ccall((:solarsystem, libnovas), Cshort,
                 (Cdouble, Cshort, Cshort, Ptr{Cdouble}, Ptr{Cdouble}), tjd, body, origin,
                 position, velocity)
end

function solarsystem_hp(tjd, body, origin, position, velocity)
    return ccall((:solarsystem_hp, libnovas), Cshort,
                 (Ptr{Cdouble}, Cshort, Cshort, Ptr{Cdouble}, Ptr{Cdouble}), tjd, body,
                 origin, position, velocity)
end

function iau2000a(jd_high, jd_low, dpsi, deps)
    return ccall((:iau2000a, libnovas), Cvoid,
                 (Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}), jd_high, jd_low, dpsi,
                 deps)
end

function iau2000b(jd_high, jd_low, dpsi, deps)
    return ccall((:iau2000b, libnovas), Cvoid,
                 (Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}), jd_high, jd_low, dpsi,
                 deps)
end

function nu2000k(jd_high, jd_low, dpsi, deps)
    return ccall((:nu2000k, libnovas), Cvoid,
                 (Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}), jd_high, jd_low, dpsi,
                 deps)
end

struct cat_entry
    starname::NTuple{51,Cchar}
    catalog::NTuple{4,Cchar}
    starnumber::Clong
    ra::Cdouble
    dec::Cdouble
    promora::Cdouble
    promodec::Cdouble
    parallax::Cdouble
    radialvelocity::Cdouble
end

struct object
    type::Cshort
    number::Cshort
    name::NTuple{51,Cchar}
    star::cat_entry
end

struct on_surface
    latitude::Cdouble
    longitude::Cdouble
    height::Cdouble
    temperature::Cdouble
    pressure::Cdouble
end

struct in_space
    sc_pos::NTuple{3,Cdouble}
    sc_vel::NTuple{3,Cdouble}
end

struct observer
    where::Cshort
    on_surf::on_surface
    near_earth::in_space
end

struct sky_pos
    r_hat::NTuple{3,Cdouble}
    ra::Cdouble
    dec::Cdouble
    dis::Cdouble
    rv::Cdouble
end

struct ra_of_cio
    jd_tdb::Cdouble
    ra_cio::Cdouble
end

function readeph(mp, name, jd, err)
    return ccall((:readeph, libnovas), Ptr{Cdouble}, (Cint, Ptr{Cchar}, Cdouble, Ptr{Cint}),
                 mp, name, jd, err)
end

function app_star(jd_tt, star, accuracy, ra, dec)
    return ccall((:app_star, libnovas), Cshort,
                 (Cdouble, Ptr{cat_entry}, Cshort, Ptr{Cdouble}, Ptr{Cdouble}), jd_tt, star,
                 accuracy, ra, dec)
end

function virtual_star(jd_tt, star, accuracy, ra, dec)
    return ccall((:virtual_star, libnovas), Cshort,
                 (Cdouble, Ptr{cat_entry}, Cshort, Ptr{Cdouble}, Ptr{Cdouble}), jd_tt, star,
                 accuracy, ra, dec)
end

function astro_star(jd_tt, star, accuracy, ra, dec)
    return ccall((:astro_star, libnovas), Cshort,
                 (Cdouble, Ptr{cat_entry}, Cshort, Ptr{Cdouble}, Ptr{Cdouble}), jd_tt, star,
                 accuracy, ra, dec)
end

function app_planet(jd_tt, ss_body, accuracy, ra, dec, dis)
    return ccall((:app_planet, libnovas), Cshort,
                 (Cdouble, Ptr{object}, Cshort, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
                 jd_tt, ss_body, accuracy, ra, dec, dis)
end

function virtual_planet(jd_tt, ss_body, accuracy, ra, dec, dis)
    return ccall((:virtual_planet, libnovas), Cshort,
                 (Cdouble, Ptr{object}, Cshort, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
                 jd_tt, ss_body, accuracy, ra, dec, dis)
end

function astro_planet(jd_tt, ss_body, accuracy, ra, dec, dis)
    return ccall((:astro_planet, libnovas), Cshort,
                 (Cdouble, Ptr{object}, Cshort, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
                 jd_tt, ss_body, accuracy, ra, dec, dis)
end

function topo_star(jd_tt, delta_t, star, position, accuracy, ra, dec)
    return ccall((:topo_star, libnovas), Cshort,
                 (Cdouble, Cdouble, Ptr{cat_entry}, Ptr{on_surface}, Cshort, Ptr{Cdouble},
                  Ptr{Cdouble}), jd_tt, delta_t, star, position, accuracy, ra, dec)
end

function local_star(jd_tt, delta_t, star, position, accuracy, ra, dec)
    return ccall((:local_star, libnovas), Cshort,
                 (Cdouble, Cdouble, Ptr{cat_entry}, Ptr{on_surface}, Cshort, Ptr{Cdouble},
                  Ptr{Cdouble}), jd_tt, delta_t, star, position, accuracy, ra, dec)
end

function topo_planet(jd_tt, ss_body, delta_t, position, accuracy, ra, dec, dis)
    return ccall((:topo_planet, libnovas), Cshort,
                 (Cdouble, Ptr{object}, Cdouble, Ptr{on_surface}, Cshort, Ptr{Cdouble},
                  Ptr{Cdouble}, Ptr{Cdouble}), jd_tt, ss_body, delta_t, position, accuracy,
                 ra, dec, dis)
end

function local_planet(jd_tt, ss_body, delta_t, position, accuracy, ra, dec, dis)
    return ccall((:local_planet, libnovas), Cshort,
                 (Cdouble, Ptr{object}, Cdouble, Ptr{on_surface}, Cshort, Ptr{Cdouble},
                  Ptr{Cdouble}, Ptr{Cdouble}), jd_tt, ss_body, delta_t, position, accuracy,
                 ra, dec, dis)
end

function mean_star(jd_tt, ra, dec, accuracy, ira, idec)
    return ccall((:mean_star, libnovas), Cshort,
                 (Cdouble, Cdouble, Cdouble, Cshort, Ptr{Cdouble}, Ptr{Cdouble}), jd_tt, ra,
                 dec, accuracy, ira, idec)
end

function place(jd_tt, cel_object, location, delta_t, coord_sys, accuracy, output)
    return ccall((:place, libnovas), Cshort,
                 (Cdouble, Ptr{object}, Ptr{observer}, Cdouble, Cshort, Cshort,
                  Ptr{sky_pos}), jd_tt, cel_object, location, delta_t, coord_sys, accuracy,
                 output)
end

function equ2gal(rai, deci, glon, glat)
    return ccall((:equ2gal, libnovas), Cvoid,
                 (Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}), rai, deci, glon, glat)
end

function equ2ecl(jd_tt, coord_sys, accuracy, ra, dec, elon, elat)
    return ccall((:equ2ecl, libnovas), Cshort,
                 (Cdouble, Cshort, Cshort, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}),
                 jd_tt, coord_sys, accuracy, ra, dec, elon, elat)
end

function equ2ecl_vec(jd_tt, coord_sys, accuracy, pos1, pos2)
    return ccall((:equ2ecl_vec, libnovas), Cshort,
                 (Cdouble, Cshort, Cshort, Ptr{Cdouble}, Ptr{Cdouble}), jd_tt, coord_sys,
                 accuracy, pos1, pos2)
end

function ecl2equ_vec(jd_tt, coord_sys, accuracy, pos1, pos2)
    return ccall((:ecl2equ_vec, libnovas), Cshort,
                 (Cdouble, Cshort, Cshort, Ptr{Cdouble}, Ptr{Cdouble}), jd_tt, coord_sys,
                 accuracy, pos1, pos2)
end

function equ2hor(jd_ut1, delta_t, accuracy, xp, yp, location, ra, dec, ref_option, zd, az,
                 rar, decr)
    return ccall((:equ2hor, libnovas), Cvoid,
                 (Cdouble, Cdouble, Cshort, Cdouble, Cdouble, Ptr{on_surface}, Cdouble,
                  Cdouble, Cshort, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
                 jd_ut1, delta_t, accuracy, xp, yp, location, ra, dec, ref_option, zd, az,
                 rar, decr)
end

function gcrs2equ(jd_tt, coord_sys, accuracy, rag, decg, ra, dec)
    return ccall((:gcrs2equ, libnovas), Cshort,
                 (Cdouble, Cshort, Cshort, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}),
                 jd_tt, coord_sys, accuracy, rag, decg, ra, dec)
end

function sidereal_time(jd_high, jd_low, delta_t, gst_type, method, accuracy, gst)
    return ccall((:sidereal_time, libnovas), Cshort,
                 (Cdouble, Cdouble, Cdouble, Cshort, Cshort, Cshort, Ptr{Cdouble}), jd_high,
                 jd_low, delta_t, gst_type, method, accuracy, gst)
end

function era(jd_high, jd_low)
    return ccall((:era, libnovas), Cdouble, (Cdouble, Cdouble), jd_high, jd_low)
end

function ter2cel(jd_ut_high, jd_ut_low, delta_t, method, accuracy, option, xp, yp, vec1,
                 vec2)
    return ccall((:ter2cel, libnovas), Cshort,
                 (Cdouble, Cdouble, Cdouble, Cshort, Cshort, Cshort, Cdouble, Cdouble,
                  Ptr{Cdouble}, Ptr{Cdouble}), jd_ut_high, jd_ut_low, delta_t, method,
                 accuracy, option, xp, yp, vec1, vec2)
end

function cel2ter(jd_ut_high, jd_ut_low, delta_t, method, accuracy, option, xp, yp, vec1,
                 vec2)
    return ccall((:cel2ter, libnovas), Cshort,
                 (Cdouble, Cdouble, Cdouble, Cshort, Cshort, Cshort, Cdouble, Cdouble,
                  Ptr{Cdouble}, Ptr{Cdouble}), jd_ut_high, jd_ut_low, delta_t, method,
                 accuracy, option, xp, yp, vec1, vec2)
end

function spin(angle, pos1, pos2)
    return ccall((:spin, libnovas), Cvoid, (Cdouble, Ptr{Cdouble}, Ptr{Cdouble}), angle,
                 pos1, pos2)
end

function wobble(tjd, direction, xp, yp, pos1, pos2)
    return ccall((:wobble, libnovas), Cvoid,
                 (Cdouble, Cshort, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}), tjd,
                 direction, xp, yp, pos1, pos2)
end

function terra(location, st, pos, vel)
    return ccall((:terra, libnovas), Cvoid,
                 (Ptr{on_surface}, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}), location, st, pos,
                 vel)
end

function e_tilt(jd_tdb, accuracy, mobl, tobl, ee, dpsi, deps)
    return ccall((:e_tilt, libnovas), Cvoid,
                 (Cdouble, Cshort, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
                  Ptr{Cdouble}), jd_tdb, accuracy, mobl, tobl, ee, dpsi, deps)
end

function cel_pole(tjd, type, dpole1, dpole2)
    return ccall((:cel_pole, libnovas), Cshort, (Cdouble, Cshort, Cdouble, Cdouble), tjd,
                 type, dpole1, dpole2)
end

function ee_ct(jd_high, jd_low, accuracy)
    return ccall((:ee_ct, libnovas), Cdouble, (Cdouble, Cdouble, Cshort), jd_high, jd_low,
                 accuracy)
end

function frame_tie(pos1, direction, pos2)
    return ccall((:frame_tie, libnovas), Cvoid, (Ptr{Cdouble}, Cshort, Ptr{Cdouble}), pos1,
                 direction, pos2)
end

function proper_motion(jd_tdb1, pos, vel, jd_tdb2, pos2)
    return ccall((:proper_motion, libnovas), Cvoid,
                 (Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Ptr{Cdouble}), jd_tdb1, pos,
                 vel, jd_tdb2, pos2)
end

function bary2obs(pos, pos_obs, pos2, lighttime)
    return ccall((:bary2obs, libnovas), Cvoid,
                 (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), pos, pos_obs,
                 pos2, lighttime)
end

function geo_posvel(jd_tt, delta_t, accuracy, obs, pos, vel)
    return ccall((:geo_posvel, libnovas), Cshort,
                 (Cdouble, Cdouble, Cshort, Ptr{observer}, Ptr{Cdouble}, Ptr{Cdouble}),
                 jd_tt, delta_t, accuracy, obs, pos, vel)
end

function light_time(jd_tdb, ss_object, pos_obs, tlight0, accuracy, pos, tlight)
    return ccall((:light_time, libnovas), Cshort,
                 (Cdouble, Ptr{object}, Ptr{Cdouble}, Cdouble, Cshort, Ptr{Cdouble},
                  Ptr{Cdouble}), jd_tdb, ss_object, pos_obs, tlight0, accuracy, pos, tlight)
end

function d_light(pos1, pos_obs)
    return ccall((:d_light, libnovas), Cdouble, (Ptr{Cdouble}, Ptr{Cdouble}), pos1, pos_obs)
end

function grav_def(jd_tdb, loc_code, accuracy, pos1, pos_obs, pos2)
    return ccall((:grav_def, libnovas), Cshort,
                 (Cdouble, Cshort, Cshort, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
                 jd_tdb, loc_code, accuracy, pos1, pos_obs, pos2)
end

function grav_vec(pos1, pos_obs, pos_body, rmass, pos2)
    return ccall((:grav_vec, libnovas), Cvoid,
                 (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Ptr{Cdouble}), pos1,
                 pos_obs, pos_body, rmass, pos2)
end

function aberration(pos, ve, lighttime, pos2)
    return ccall((:aberration, libnovas), Cvoid,
                 (Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Ptr{Cdouble}), pos, ve, lighttime,
                 pos2)
end

function rad_vel(cel_object, pos, vel, vel_obs, d_obs_geo, d_obs_sun, d_obj_sun, rv)
    return ccall((:rad_vel, libnovas), Cvoid,
                 (Ptr{object}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cdouble, Cdouble,
                  Cdouble, Ptr{Cdouble}), cel_object, pos, vel, vel_obs, d_obs_geo,
                 d_obs_sun, d_obj_sun, rv)
end

function precession(jd_tdb1, pos1, jd_tdb2, pos2)
    return ccall((:precession, libnovas), Cshort,
                 (Cdouble, Ptr{Cdouble}, Cdouble, Ptr{Cdouble}), jd_tdb1, pos1, jd_tdb2,
                 pos2)
end

function nutation(jd_tdb, direction, accuracy, pos, pos2)
    return ccall((:nutation, libnovas), Cvoid,
                 (Cdouble, Cshort, Cshort, Ptr{Cdouble}, Ptr{Cdouble}), jd_tdb, direction,
                 accuracy, pos, pos2)
end

function nutation_angles(t, accuracy, dpsi, deps)
    return ccall((:nutation_angles, libnovas), Cvoid,
                 (Cdouble, Cshort, Ptr{Cdouble}, Ptr{Cdouble}), t, accuracy, dpsi, deps)
end

function fund_args(t, a)
    return ccall((:fund_args, libnovas), Cvoid, (Cdouble, Ptr{Cdouble}), t, a)
end

function mean_obliq(jd_tdb)
    return ccall((:mean_obliq, libnovas), Cdouble, (Cdouble,), jd_tdb)
end

function vector2radec(pos, ra, dec)
    return ccall((:vector2radec, libnovas), Cshort,
                 (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), pos, ra, dec)
end

function radec2vector(ra, dec, dist, vector)
    return ccall((:radec2vector, libnovas), Cvoid,
                 (Cdouble, Cdouble, Cdouble, Ptr{Cdouble}), ra, dec, dist, vector)
end

function starvectors(star, pos, vel)
    return ccall((:starvectors, libnovas), Cvoid,
                 (Ptr{cat_entry}, Ptr{Cdouble}, Ptr{Cdouble}), star, pos, vel)
end

function tdb2tt(tdb_jd, tt_jd, secdiff)
    return ccall((:tdb2tt, libnovas), Cvoid, (Cdouble, Ptr{Cdouble}, Ptr{Cdouble}), tdb_jd,
                 tt_jd, secdiff)
end

function cio_ra(jd_tt, accuracy, ra_cio)
    return ccall((:cio_ra, libnovas), Cshort, (Cdouble, Cshort, Ptr{Cdouble}), jd_tt,
                 accuracy, ra_cio)
end

function cio_location(jd_tdb, accuracy, ra_cio, ref_sys)
    return ccall((:cio_location, libnovas), Cshort,
                 (Cdouble, Cshort, Ptr{Cdouble}, Ptr{Cshort}), jd_tdb, accuracy, ra_cio,
                 ref_sys)
end

function cio_basis(jd_tdb, ra_cio, ref_sys, accuracy, x, y, z)
    return ccall((:cio_basis, libnovas), Cshort,
                 (Cdouble, Cdouble, Cshort, Cshort, Ptr{Cdouble}, Ptr{Cdouble},
                  Ptr{Cdouble}), jd_tdb, ra_cio, ref_sys, accuracy, x, y, z)
end

function cio_array(jd_tdb, n_pts, cio)
    return ccall((:cio_array, libnovas), Cshort, (Cdouble, Clong, Ptr{ra_of_cio}), jd_tdb,
                 n_pts, cio)
end

function ira_equinox(jd_tdb, equinox, accuracy)
    return ccall((:ira_equinox, libnovas), Cdouble, (Cdouble, Cshort, Cshort), jd_tdb,
                 equinox, accuracy)
end

function ephemeris(jd, cel_obj, origin, accuracy, pos, vel)
    return ccall((:ephemeris, libnovas), Cshort,
                 (Ptr{Cdouble}, Ptr{object}, Cshort, Cshort, Ptr{Cdouble}, Ptr{Cdouble}),
                 jd, cel_obj, origin, accuracy, pos, vel)
end

function transform_hip(hipparcos, hip_2000)
    return ccall((:transform_hip, libnovas), Cvoid, (Ptr{cat_entry}, Ptr{cat_entry}),
                 hipparcos, hip_2000)
end

function transform_cat(option, date_incat, incat, date_newcat, newcat_id, newcat)
    return ccall((:transform_cat, libnovas), Cshort,
                 (Cshort, Cdouble, Ptr{cat_entry}, Cdouble, Ptr{Cchar}, Ptr{cat_entry}),
                 option, date_incat, incat, date_newcat, newcat_id, newcat)
end

function limb_angle(pos_obj, pos_obs, limb_ang, nadir_ang)
    return ccall((:limb_angle, libnovas), Cvoid,
                 (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), pos_obj, pos_obs,
                 limb_ang, nadir_ang)
end

function refract(location, ref_option, zd_obs)
    return ccall((:refract, libnovas), Cdouble, (Ptr{on_surface}, Cshort, Cdouble),
                 location, ref_option, zd_obs)
end

function julian_date(year, month, day, hour)
    return ccall((:julian_date, libnovas), Cdouble, (Cshort, Cshort, Cshort, Cdouble), year,
                 month, day, hour)
end

function cal_date(tjd, year, month, day, hour)
    return ccall((:cal_date, libnovas), Cvoid,
                 (Cdouble, Ptr{Cshort}, Ptr{Cshort}, Ptr{Cshort}, Ptr{Cdouble}), tjd, year,
                 month, day, hour)
end

function norm_ang(angle)
    return ccall((:norm_ang, libnovas), Cdouble, (Cdouble,), angle)
end

function make_cat_entry(star_name, catalog, star_num, ra, dec, pm_ra, pm_dec, parallax,
                        rad_vel_, star)
    return ccall((:make_cat_entry, libnovas), Cshort,
                 (Ptr{Cchar}, Ptr{Cchar}, Clong, Cdouble, Cdouble, Cdouble, Cdouble,
                  Cdouble, Cdouble, Ptr{cat_entry}), star_name, catalog, star_num, ra, dec,
                 pm_ra, pm_dec, parallax, rad_vel_, star)
end

function make_object(type, number, name, star_data, cel_obj)
    return ccall((:make_object, libnovas), Cshort,
                 (Cshort, Cshort, Ptr{Cchar}, Ptr{cat_entry}, Ptr{object}), type, number,
                 name, star_data, cel_obj)
end

function make_observer(where, obs_surface, obs_space, obs)
    return ccall((:make_observer, libnovas), Cshort,
                 (Cshort, Ptr{on_surface}, Ptr{in_space}, Ptr{observer}), where,
                 obs_surface, obs_space, obs)
end

function make_observer_at_geocenter(obs_at_geocenter)
    return ccall((:make_observer_at_geocenter, libnovas), Cvoid, (Ptr{observer},),
                 obs_at_geocenter)
end

function make_observer_on_surface(latitude, longitude, height, temperature, pressure,
                                  obs_on_surface)
    return ccall((:make_observer_on_surface, libnovas), Cvoid,
                 (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{observer}), latitude,
                 longitude, height, temperature, pressure, obs_on_surface)
end

function make_observer_in_space(sc_pos, sc_vel, obs_in_space)
    return ccall((:make_observer_in_space, libnovas), Cvoid,
                 (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{observer}), sc_pos, sc_vel, obs_in_space)
end

function make_on_surface(latitude, longitude, height, temperature, pressure, obs_surface)
    return ccall((:make_on_surface, libnovas), Cvoid,
                 (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{on_surface}), latitude,
                 longitude, height, temperature, pressure, obs_surface)
end

function make_in_space(sc_pos, sc_vel, obs_space)
    return ccall((:make_in_space, libnovas), Cvoid,
                 (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{in_space}), sc_pos, sc_vel, obs_space)
end

function ephem_open(ephem_name, jd_begin, jd_end, de_number)
    return ccall((:ephem_open, libnovas), Cshort,
                 (Ptr{Cchar}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cshort}), ephem_name,
                 jd_begin, jd_end, de_number)
end

function ephem_close()
    return ccall((:ephem_close, libnovas), Cshort, ())
end

function planet_ephemeris(tjd, target, center, position, velocity)
    return ccall((:planet_ephemeris, libnovas), Cshort,
                 (Ptr{Cdouble}, Cshort, Cshort, Ptr{Cdouble}, Ptr{Cdouble}), tjd, target,
                 center, position, velocity)
end

function state(jed, target, target_pos, target_vel)
    return ccall((:state, libnovas), Cshort,
                 (Ptr{Cdouble}, Cshort, Ptr{Cdouble}, Ptr{Cdouble}), jed, target,
                 target_pos, target_vel)
end

function interpolate(buf, t, ncm, na, position, velocity)
    return ccall((:interpolate, libnovas), Cvoid,
                 (Ptr{Cdouble}, Ptr{Cdouble}, Clong, Clong, Ptr{Cdouble}, Ptr{Cdouble}),
                 buf, t, ncm, na, position, velocity)
end

function split(tt, fr)
    return ccall((:split, libnovas), Cvoid, (Cdouble, Ptr{Cdouble}), tt, fr)
end

const SIZE_OF_OBJ_NAME = 51

const SIZE_OF_CAT_NAME = 4

const BARYC = 0

const HELIOC = 1

end


text:
# Carries along the `u` which is an allocation to save when no callbacks
function handle_callbacks!(integrator,eval_sol_fcn)
  discrete_callbacks = integrator.opts.callback.discrete_callbacks
  continuous_callbacks = integrator.opts.callback.continuous_callbacks
  atleast_one_callback = false

  continuous_modified = false
  discrete_modified = false
  saved_in_cb = false
  if !(typeof(continuous_callbacks)<:Tuple{})
    time,upcrossing,event_occured,idx,counter =
              DiffEqBase.find_first_continuous_callback(integrator,continuous_callbacks...)
    if event_occured
      integrator.event_last_time = idx
      continuous_modified,saved_in_cb = DiffEqBase.apply_callback!(integrator,continuous_callbacks[idx],time,upcrossing)
    else
      integrator.event_last_time = 0
    end
  end
  if !(typeof(discrete_callbacks)<:Tuple{})
    discrete_modified,saved_in_cb = DiffEqBase.apply_discrete_callback!(integrator,discrete_callbacks...)
  end
  if !saved_in_cb
    savevalues!(integrator)
  end

  integrator.u_modified = continuous_modified || discrete_modified
end

function DiffEqBase.savevalues!(integrator::ODEInterfaceIntegrator,force_save=false)::Tuple{Bool,Bool}
  saved, savedexactly = false, false
  !integrator.opts.save_on && return saved, savedexactly
  uType = eltype(integrator.sol.u)

  if integrator.opts.save_everystep || force_save
    saved = true
    push!(integrator.sol.t,integrator.t)
    save_value!(integrator.sol.u,copy(integrator.u),uType,integrator.sizeu)
  end

  while !isempty(integrator.opts.saveat) &&
    integrator.tdir*top(integrator.opts.saveat) < integrator.tdir*integrator.t
    saved = true
    curt = pop!(integrator.opts.saveat)
    tmp = integrator(curt)::Vector{Float64}
    push!(integrator.sol.t,curt)
    save_value!(integrator.sol.u,tmp,uType,integrator.sizeu)
  end
  savedexactly = last(integrator.sol.t) == integrator.t
  return saved, savedexactly
end

function DiffEqBase.change_t_via_interpolation!(integrator::ODEInterfaceIntegrator,t)
  integrator.t = t
  tmp = integrator(integrator.t)::Vector{Float64}
  if eltype(integrator.sol.u) <: Vector
      integrator.u .= tmp
  else
      integrator.u .= reshape(tmp,integrator.sizeu)
  end
  nothing
end
DiffEqBase.get_tmp_cache(i::ODEInterfaceIntegrator,args...) = nothing

@inline function Base.getproperty(integrator::ODEInterfaceIntegrator, sym::Symbol)
  if sym == :dt
    return integrator.t-integrator.tprev
  else
    return getfield(integrator, sym)
  end
end


text:
module TestUtils

using Test
using RLEVectors

@testset begin

    # rep
    @test RLEVectors.rep([4, 5, 6], each = 2) == [4, 4, 5, 5, 6, 6]
    @test RLEVectors.rep([4, 5, 6], times = 2) == [4, 5, 6, 4, 5, 6]
    @test RLEVectors.rep([4, 5, 6], each = 3, times = 2) ==
          [4, 4, 4, 5, 5, 5, 6, 6, 6, 4, 4, 4, 5, 5, 5, 6, 6, 6]
    @test RLEVectors.rep(["p", "d", "q"], each = 3, times = 2) ==
          ["p", "p", "p", "d", "d", "d", "q", "q", "q", "p", "p", "p", "d", "d", "d", "q", "q", "q"]
    @test_throws ArgumentError RLEVectors.rep([1, 2, 3], each = [4, 5])
    @test rep(3, each = 2) == [3, 3]

    # searchsortedfirst
    @test RLEVectors.searchsortedfirst([3, 6, 9, 12, 15], 6, 1, 5) == 2
    @test RLEVectors.searchsortedfirst([3, 6, 9, 12, 15], 7, 1, 5) == 3
    @test RLEVectors.searchsortedfirst([3, 6, 9, 12, 15], 7, 0, 22) == 3
    @test RLEVectors.searchsortedfirst([3, 6, 9, 12, 15], 40, 1, 5) == 6
    @test RLEVectors.searchsortedfirst([3, 6, 9, 12, 15], -2, 1, 5) == 1
    @test RLEVectors.searchsortedfirst([3, 6, 9, 12, 15], 14, 4, 5) == 5
    @test RLEVectors.searchsortedfirst([3, 6, 9, 12, 15], 3, 4, 5) == 4

    @test RLEVectors.searchsortedfirst([0, 5, 10, 15], [-3, 2, 3, 7, 22]) == [1, 2, 2, 3, 5]
    @test RLEVectors.searchsortedfirst([0, 5, 10, 15], [-3, 2, -3, 7, 22]) == [1, 2, 1, 3, 5]

    # growat!
    x = collect(1:10)
    RLEVectors.growat!(x, 3, 2)
    @test x == [1, 2, 3, 4, 3, 4, 5, 6, 7, 8, 9, 10]
    x = RLEVector([1, 1, 2, 2, 3, 3])
    RLEVectors.growat!(x, 2, 1)
    @test x.runvalues[3:4] == [2, 3]
    @test x.runends[3:4] == [4, 6]

end # testset

end # module


text:
export spgline

"""
Non-monotone linesearch
"""
function spgline(A::TA,
                 f::Tf,
                 d::AbstractVector{<:Number},
                 gtd_in::Number,
                 x::AbstractArray{ETx},
                 fMax::Tf,
                 funForward::Function,
                 funPenalty::Function,
                 params::Dict{String,Any},
                 b::AbstractVector{ETb},
                 feasSrchIt::Int64,
                 linear::Bool,
                 options::spgOptions,
                 timeProject::Float64) where
                   {TA<:Union{joAbstractLinearOperator, AbstractArray},
                    Tf<:Number, ETx<:Number, ETb<:Number}

    (options.verbosity > 1) && println("Script entered spgline for A::Function")

    localProdA = 0
    EXIT_CONVERGED = 0
    EXIT_LINEITERATIONS = 1
    maxIts = feasSrchIt
    step = one(real(ETx))
    iter = zero(Int64)
    gamma = 1e-4
    gtd = -abs(gtd_in)

    if linear
        Ad = funForward(A, d, Array{ETx,1}(), params)[1]
        localProdA += 1
        r = copy(x)
    end

    # Define outputs outside of while loop scope
    rNew::Array{ETb,1} = Array{ETb,1}()
    fNew::Tf = zero(Tf)
    err = -1

    while true

        if linear
            rNew = r .- step*Ad
        else
            rNew = b - funForward(A, x + step*d, Array{ETx,1}(), params)[1]
            localProdA += 1
        end

        fNew = funPenalty(rNew, params)[1]

        # Check exit conditions
        (options.verbosity > 1) && println(""" fNew: $fNew
                    fMax: $(fMax)
                    add $(gamma*step*gtd)""")
        if fNew <= fMax + gamma*step*gtd
            err = EXIT_CONVERGED
            break
        elseif iter >= maxIts
            err = EXIT_LINEITERATIONS
            break
        end

        iter += 1
        (options.verbosity > 1) && println("Line Search Iter: $iter")
        if step <= 0.1
            step /= 2
        else
            tmp = (-gtd*step^2)/(2*(fNew-f-step*gtd))
            if (tmp < 0.1) | (tmp > 0.9*step) | (isnan(tmp))
                tmp = step/2
            end
            step = copy(tmp)
        end

    end # while

    return fNew, step, rNew, iter, err, localProdA

end # spgline

"""
Non-monotone linesearch
"""
function spgline(A::TA,
                 f::Tf,
                 d::AbstractVector{<:Number},
                 gtd_in::Number,
                 x::AbstractArray{ETx},
                 fMax::AbstractFloat,
                 funForward::Function,
                 funPenalty::Function,
                 params::Dict{String,Any},
                 b::AbstractVector{ETb},
                 feasSrchIt::Int64,
                 linear::Bool,
                 options::spgOptions,
                 timeProject::Float64) where
                   {TA<:Function, Tf<:Number, ETx<:Number, ETb<:Number}

    (options.verbosity > 1) && println("Script entered spgline for A::Function")

    localProdA = 0
    EXIT_CONVERGED = 0
    EXIT_LINEITERATIONS = 1
    maxIts = feasSrchIt
    step::real(ETx) = one(real(ETx))
    iter = zero(Int64)
    gamma = 1e-4
    gtd = -abs(gtd_in)

    if linear
        Ad = funForward(A, d, Arrat{ETx,1}(), params)[1]
        localProdA += 1
        r = copy(x)
    end

    # Define outputs outside of while loop scope
    rNew = Array{ETb,1}()
    fNew = zero(Tf)
    err = -1

    while true

        if linear
            rNew = r - step*Ad
        else
            rNew = b - funForward(A, x + step*d, Array{ETx,1}(), params)[1]
            localProdA += 1
        end

        fNew = funPenalty(rNew, params)[1]

        # Check exit conditions
        (options.verbosity > 1) && println(""" fNew: $fNew
                    fMax: $(fMax)
                    add $(gamma*step*gtd)""")
        if fNew <= fMax + gamma*step*gtd
            err = EXIT_CONVERGED
            break
        elseif iter >= maxIts
            err = EXIT_LINEITERATIONS
            break
        end

        iter += 1
        (options.verbosity > 1) && println("Line Search Iter: $iter")
        if step <= 0.1
            step /= 2
        else
            tmp = (-gtd*step^2)/(2*(fNew-f-step*gtd))
            if (tmp < 0.1) | (tmp > 0.9*step) | (isnan(tmp))
                tmp = step/2
            end
            step = copy(tmp)
        end

    end # while

    return fNew, step, rNew, iter, err, localProdA

end # spgline


text:
export propensitygrad_sparsity_pattern

function propensitygrad_sparsity_pattern(speciescount, parametercount, propensities, parameters)
    propensitycount = length(propensities)
    ModelingToolkit.@variables t, x[1:speciescount](t), alphas[1:propensitycount](t)
    ModelingToolkit.@parameters p[1:parametercount]
    eqs = Vector{ModelingToolkit.Equation}([])
    for r in 1:propensitycount
        if istimevarying(propensities[r])
            push!(eqs, alphas[r]~propensities[r](t, x, p))
        else
            push!(eqs, alphas[r]~propensities[r](x, p))
        end
    end
    ModelingToolkit.@named sys = ModelingToolkit.ODESystem(eqs,t,[x;alphas],p)
    eqdep = ModelingToolkit.equation_dependencies(sys; variables=[p...])
    depgraph = ModelingToolkit.asgraph(eqdep, Dict([s=>i for (i,s) in enumerate(p)]))

    I = Vector{UInt32}()
    J = Vector{UInt32}()
    vals = Vector{Bool}()
    for r in 1:propensitycount
        for ip in depgraph.fadjlist[r]
            push!(I, r)
            push!(J, ip)
            push!(vals, true)
        end
    end

    return sparse(I,J,vals)
end


text:
const MKL_DSS_DEFAULTS = 0
const MKL_DSS_NON_SYMMETRIC = 536871104
const MKL_DSS_SUCCESS = 0
const MKL_DSS_AUTO_ORDER = 268435520
const MKL_DSS_POSITIVE_DEFINITE = 134217792

function dss_analyze(A::SparseMatrixCSC)
  handle = Int[0]
  opt = MKL_DSS_DEFAULTS
  error = ccall((:dss_create, Sparso.LIB_PATH), Cint,
        (Ptr{Void}, Ptr{Cint}),
        handle, &opt)
  if error != MKL_DSS_SUCCESS
    println("dss_create returned error code $error")
  end

  opt = MKL_DSS_NON_SYMMETRIC
  m = size(A, 1)
  nnz = A.colptr[m + 1] - 1
  error = ccall((:dss_define_structure, Sparso.LIB_PATH), Cint,
        (Ptr{Void}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
        handle, &opt, A.colptr, &m, &m, A.rowval, &nnz)
  if error != MKL_DSS_SUCCESS
    println("dss_define_structure returned error code $error")
  end

  opt = MKL_DSS_AUTO_ORDER
  error = ccall((:dss_reorder, Sparso.LIB_PATH), Cint,
                (Ptr{Void}, Ptr{Cint}, Ptr{Cint}),
                handle, &opt, C_NULL)
  if error != MKL_DSS_SUCCESS
    println("dss_reorder returned error code $error")
  end

  handle
end

function dss_factor(handle, A::SparseMatrixCSC)
  opt = MKL_DSS_DEFAULTS
  error = ccall((:dss_factor_real, Sparso.LIB_PATH), Cint,
                (Ptr{Void}, Ptr{Cint}, Ptr{Cdouble}),
                handle, &opt, A.nzval)
  if error != MKL_DSS_SUCCESS
    println("dss_factor_real returned error code $error")
  end
end

# solve A*sol = rhs
function dss_solve!(handle, rhs::Vector, sol::Vector)
  opt = MKL_DSS_DEFAULTS
  nrhs = 1
  error = ccall((:dss_solve_real, Sparso.LIB_PATH), Cint,
                (Ptr{Void}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cdouble}),
                handle, &opt, rhs, &nrhs, sol)
  ccall((:dss_delete, Sparso.LIB_PATH), Cint,
        (Ptr{Void}, Ptr{Cint}), handle, &opt)
end

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
using Turing: gradient_logp_forward, gradient_logp_reverse
using Test

function test_ad(f, at = 0.5; rtol = 1e-8, atol = 1e-8)
    isarr = isa(at, AbstractArray)
    reverse = Tracker.data(Tracker.gradient(f, at)[1])
    if isarr
        forward = ForwardDiff.gradient(f, at)
        @test isapprox(reverse, forward, rtol=rtol, atol=atol)
    else
        forward = ForwardDiff.derivative(f, at)
        finite_diff = central_fdm(5,1)(f, at)
        @test isapprox(reverse, forward, rtol=rtol, atol=atol)
        @test isapprox(reverse, finite_diff, rtol=rtol, atol=atol)
    end
end

"""
    test_reverse_mode_ad(forward, f, ȳ, x...; rtol=1e-8, atol=1e-8)

Check that the reverse-mode sensitivities produced by an AD library are correct for `f`
at `x...`, given sensitivity `ȳ` w.r.t. `y = f(x...)` up to `rtol` and `atol`.
`forward` should be either `Tracker.forward` or `Zygote.forward`.
"""
function test_reverse_mode_ad(forward, f, ȳ, x...; rtol=1e-8, atol=1e-8)

    # Perform a regular forwards-pass.
    y = f(x...)

    # Use tracker to compute reverse-mode sensitivities.
    y_tracker, back = forward(f, x...)
    x̄s_tracker = back(ȳ)

    # Use finite differencing to compute reverse-mode sensitivities.
    x̄s_fdm = FDM.j′vp(central_fdm(5, 1), f, ȳ, x...)
    if length(x) == 1
        x̄s_fdm = (x̄s_fdm,)
    end

    # Check that forwards-pass produces the correct answer.
    @test y ≈ y_tracker

    # Check that reverse-mode sensitivities are correct.
    @test all([x̄_tracker ≈ x̄_fdm for (x̄_tracker, x̄_fdm) in zip(x̄s_tracker, x̄s_fdm)])
end

# See `test_reverse_mode_ad` for details.
function test_tracker_ad(f, ȳ, x...; rtol=1e-8, atol=1e-8)
    return test_reverse_mode_ad(Tracker.forward, f, ȳ, x...; rtol=rtol, atol=atol)
end

function test_model_ad(model, f, syms::Vector{Symbol})
    # Set up VI.
    vi = Turing.VarInfo(model)

    # Collect symbols.
    vnms = Vector(undef, length(syms))
    vnvals = Vector{Float64}()
    for i in 1:length(syms)
        s = syms[i]
        vnms[i] = getfield(vi.metadata, s).vns[1]

        vals = getval(vi, vnms[i])
        for i in eachindex(vals)
            push!(vnvals, vals[i])
        end
    end

    spl = SampleFromPrior()
    _, ∇E = gradient_logp_forward(vi[spl], vi, model)
    grad_Turing = sort(∇E)

    # Call ForwardDiff's AD
    grad_FWAD = sort(ForwardDiff.gradient(f, vec(vnvals)))

    # Compare result
    @test grad_Turing ≈ grad_FWAD atol=1e-9
end


text:
using Random

import Random: rand, SamplerTrivial, CloseOpen12_64, CloseOpen01_64, BitInteger, UInt52Raw, CloseOpen12, SamplerUnion, SamplerType

"""
    MutableRandomNumbers([rng_base=GLOBAL_RNG], size, mode:=static)

Create a `MutableRandomNumbers` RNG object with a vector of size `size`
containing floating-point random_numbers initially generated with `rng_base`. Random
numbers are then generated sequentially from this vector. The RNG object can be
initialized in two modes:
* :static - then there is an exception thrown once all random numbers are used
* :dynamic - then there are new random numbers generated on the flow from the (copied) rng

# Examples
```jldoctest
julia> rng = MutableRandomNumbers(MersenneTwister(1234),100);
julia> x1 = rand(rng, 2)
2-element Array{Float64,1}:
 0.5908446386657102
 0.7667970365022592
julia> rng = MersenneTwister(1234);
julia> x2 = rand(rng, 2)
2-element Array{Float64,1}:
 0.5908446386657102
 0.7667970365022592
julia> x1 == x2
true
```

Importantly, the random numbers can be accessed and manipulated as in ordinary
array objects.

```jldoctest
julia> rng[1]
0.5908446386657102
julia> rng[3] = 0.2
julia> rand(rng)
0.2
```

Use `reset` in order to rerun the (manipulated) `random` number sequence.

```jldoctest
julia> reset(rng)

```

"""
mutable struct MutableRandomNumbers <: AbstractRNG
    random_numbers::Vector{Float64}
    index_current::Int
    mode::Symbol
    rng_base::AbstractRNG

    function MutableRandomNumbers(rng_base::AbstractRNG, size::Int;
                                  mode::Symbol=:static)
        new(rand(rng_base, size), 0, mode, copy(rng_base))
    end
end

MutableRandomNumbers(rng_base::AbstractRNG; mode=:static) = MutableRandomNumbers(rng_base, 1, mode=mode)
MutableRandomNumbers(size::Int; mode=:static) = MutableRandomNumbers(Random.GLOBAL_RNG, size, mode=mode)
MutableRandomNumbers(; mode=:static) = MutableRandomNumbers(Random.GLOBAL_RNG, 1, mode=mode)

"""
    reset(rng::MutableRandomNumbers, [index::Int=0])

Reset the state of a MutableRandomNumbers object `rng` to `index`. Default
resets the RNG object to the pre-initial index (0)
"""
function reset(r::MutableRandomNumbers, index::Int)
    r.index_current = index
end

reset(r::MutableRandomNumbers) = reset(r,0)

# easy access overloads
Base.length(r::MutableRandomNumbers) = Base.length(r.random_numbers)
Base.getindex(r::MutableRandomNumbers, i::Int) = Base.getindex(r.random_numbers, i)

function Base.setindex!(r::MutableRandomNumbers, value::Real, i::Int)
    if 0 <= value < 1
        Base.setindex!(r.random_numbers, value, i)
    else
        throw(DomainError(value, "random numbers have to be in the domain [0,1)"))
    end
end

# take the next random number from the generator
function gen_rand01(r::MutableRandomNumbers)
    r.index_current += 1
    if checkbounds(Bool, r.random_numbers, r.index_current)
        return @inbounds r.random_numbers[r.index_current]
    else
        if r.mode == :dynamic
            push!(r.random_numbers, rand(r.rng_base))
            return @inbounds r.random_numbers[r.index_current]
        else
            throw(DomainError(r.index_current, "MutableRandomNumbers object has no further random numbers left and is configured in mode `:static`"))
        end
    end
end

#### helper functions

# random numbers work in the range [1,2) because of floating point
# precision, so we have to shift our stored random numbers again
rand_inbounds(r::MutableRandomNumbers, ::CloseOpen12_64) = gen_rand01(r) + 1
rand_inbounds(r::MutableRandomNumbers, ::CloseOpen01_64=CloseOpen01()) =
    rand_inbounds(r, CloseOpen12()) - 1.0

rand_inbounds(r::MutableRandomNumbers, ::UInt52Raw{T}) where {T<:BitInteger} =
    reinterpret(UInt64, rand_inbounds(r, CloseOpen12())) % T

#### generation

function rand(r::MutableRandomNumbers, x::SamplerTrivial{UInt52Raw{UInt64}})
    rand_inbounds(r, x[])
end

rand(r::MutableRandomNumbers, sp::SamplerTrivial{CloseOpen12_64}) = rand_inbounds(r, sp[])

# currently no support for integer cast of random numbers ...
#rand(r::MersenneTwister, T::SamplerUnion(Bool, Int8, UInt8, Int16, Int64, Int128, UInt16, Int32, UInt64, UInt32, UInt128)) =
#    rand(r, UInt52Raw()) % T[]