text:
@graphs begin
	minion_prio = """
	Select a minion that can be reqruited.

	Priority:
	1. Saruman
	2. Witch King
	3. Mouth of Sauron
	"""

	recruit_saruman = """
	Recruit Saruman.
	"""

	recruit_wk = """
	Recruit the Witch King, place it in a valid region with an army.

	Priority:
	1. Army is *mobile*
	2. Army's *target* is in nation at war
	3. Army becomes *mobile* if the Witch King is added
	4. Free Peoples' army at *target* or on the route to *target* does not contain Gandalf the White
	5. Free Peoples' army at *target* or on the route to *target* does not contain a hobbit
	6. Army is adjacent to a *threat*
	7. Army that is conduction a siege
	8. Army is adjacent to its *target*
	9. Highest *value* Shadow army
	10. Random
	"""

	recruit_mos = """
	Recruit Mouth of Sauron, place it in a valid region with an army or stronghold.

	Priority:
	1. Army is conducting a siege.
	2. Army is *mobile*
	3. Army becomes *mobile* if Mouth of Sauron is added
	4. Army contains Saruman
	5. Army with the highest *value*
	6. Stronghold closest to army whose *target* is in a nation at war
	7. Stronghold closest to army whose *target* is in an active nation
	8. Stronghold closest to army whose *target* is in a passive nation
	9. Random
	"""

	################################################################################
	@node muster_minion = Start() -> m_2_reserved_check
	@node m_2_reserved_check = BinaryCondition("""
											   A die has been reserved for recruiting a minion as a last action.
											   """) -> [n_true = m_2_return, n_false = m_2]
	@node m_2 = BinaryCondition("""
								Minion can be recruited.
								""") -> [n_true = m_2_1, n_false = m_2_return]
	@node m_2_return = ReturnFromGraph() -> []

	@node m_2_1 = BinaryCondition("""
								  The Free Peoples' have a Will of the West die.
								  And, Gandalf the White has not been recruited.
								  And, no minion have been recruited.
								  """) -> [n_true = m_2_1_yes, n_false = m_2_minion_selection]
	@node m_2_1_yes = UseActiveDie() -> m_2_1_reserve
	@node m_2_1_reserve = PerformAction("""
										Set aside the die (and ring if necessary). Use this die as an last action to recruit a minion. Minion selection and placement can be made from the main menu of Phase 5.
										""") -> m_2_1_return
	@node m_2_1_return = ReturnFromGraph() -> []

	@node m_2_minion_selection = MultipleChoice(minion_prio) -> [m_2_saruman_die, m_2_wk_die, m_2_mos_die]

	@node m_2_saruman_die = UseActiveDie() -> m_2_saruman_placement
	@node m_2_saruman_placement = PerformAction(recruit_saruman) -> m_2_saruman_end

	@node m_2_wk_die = UseActiveDie() -> m_2_wk_placement
	@node m_2_wk_placement = PerformAction(recruit_wk) -> m_2_wk_end

	@node m_2_mos_die = UseActiveDie() -> m_2_mos_placement
	@node m_2_mos_placement = PerformAction(recruit_mos) -> m_2_mos_end

	@node m_2_saruman_end = End() -> []
	@node m_2_wk_end = End() -> []
	@node m_2_mos_end = End() -> []

	################################################################################
	@node muster_minion_selection = Start() -> m_2_minion_selection_last
	@node m_2_minion_selection_last = MultipleChoice(minion_prio) -> [m_2_saruman_last, m_2_wk_last, m_2_mos_last]

	@node m_2_saruman_last = PerformAction(recruit_saruman) -> m_2_saruman_end_last
	@node m_2_wk_last = PerformAction(recruit_wk) -> m_2_wk_end_last
	@node m_2_mos_last = PerformAction(recruit_mos) -> m_2_mos_end_last

	@node m_2_saruman_end_last = End() -> []
	@node m_2_wk_end_last = End() -> []
	@node m_2_mos_end_last = End() -> []

	################################################################################
	@node muster_politics = Start() -> m_3
	@node m_3 = BinaryCondition("A Shadow nation is not at war.") -> [n_true = m_3_yes, n_false = m_3_return]
	@node m_3_return = ReturnFromGraph() -> []

	@node m_3_yes = UseActiveDie() -> m_3_action
	@node m_3_action = PerformAction("""
									 Move a nation down one step on the political track.

									 Priority:
									 1. Isengard
									 2. Sauron
									 3. Southrons and Easterlings
									 """) -> m_3_end
	@node m_3_end = End() -> []

	################################################################################
	@node muster_muster = Start() -> m_4
	@node m_4 = BinaryCondition("A card that musters is *playable*.") -> [n_true = m_4_die, n_false = m_5]
	@node m_4_die = UseActiveDie() -> m_4_action
	@node m_4_action = PerformAction("""
									 Play a *playable* muster card.

									 Priority:
									 1. Ascending order of initiative
									 2. Random
									 """) -> m_4_end
	@node m_4_end = End() -> []

	@node m_5 = BinaryCondition("""
								Muster is possible.
								""") -> [n_true = m_6, n_false = m_return]
	@node m_return = ReturnFromGraph() -> []

	@node m_6 = BinaryCondition("Muster can create an *exposed* region.") -> [n_true = m_6_die, n_false = m_7]
	@node m_6_die = UseActiveDie() -> m_6_action
	@node m_6_action = PerformAction("""
									 *Focus* priority:
									 1. Region which creates an *exposed* region
									 2. Random

									 Muster:
									 *Primary*: Elite
									 *Secondary*: Regular

									 If unit is unavailable, rotate as:
									 Elite -> Regular -> Nazgûl -> Elite
									 """) -> m_6_end
	@node m_6_end = End() -> []

	@node m_7 = BinaryCondition("""
								The Fellowship is adjacent to, or in, a region it is possible to muster in.
								And, the progress put the Fellowship outside Mordor.
								And, no army is adjacent to the Fellowship's current region.
								""") -> [n_true = m_7_die, n_false = m_8]
	@node m_7_die = UseActiveDie() -> m_7_action
	@node m_7_action = PerformAction("""
									 *Focus* priority:
									 1. The Fellowship's current region

									 Muster:
									 *Primary*: Regular
									 *Secondary*: Nazgûl

									 If unit is unavailable, rotate as:
									 Elite -> Regular -> Nazgûl -> Elite
									 """) -> m_7_end
	@node m_7_end = End() -> []

	@node m_8 = BinaryCondition("Muster is possible in a region containing a Shadow army.") -> [n_true = m_8_die, n_false = m_9]
	@node m_8_die = UseActiveDie() -> m_8_action
	@node m_8_action = PerformAction("""
									 *Focus* priority:
									 1. Army is conducting a siege.
									 2. Army is *mobile*
									 3. Army becomes mobile if Mouth of Sauron is added
									 4. Army contains Saruman
									 5. Army with the highest *value*
									 6. Random

									 Muster:
									 *Primary*: Regular
									 *Secondary*: Nazgûl

									 If unit is unavailable, rotate as:
									 Elite -> Regular -> Nazgûl -> Elite
									 """) -> m_8_end
	@node m_8_end = End() -> []

	@node m_9 = BinaryCondition("Less than 6 Nazgûl are in play.") -> [n_true = m_10_yes, n_false = m_10_no]
	@node m_10_yes = UseActiveDie() -> m_10_yes_action
	@node m_10_yes_action = PerformAction("""
										  *Focus* priority:
										  1. Closest to army whose *target* is in a nation at war
										  2. Closest to army whose *target* is in an active nation.
										  3. Closest to a *mobile* army
										  4. Closest to army whose *target* is in a passive nation.
										  5. Random

										  Muster:
										  *Primary*: Nazgûl
										  *Secondary*: Nazgûl

										  If unit is unavailable, rotate as:
										  Elite -> Regular -> Nazgûl -> Elite
										  """) -> m_10_end_yes
	@node m_10_end_yes = End() -> []

	@node m_10_no = UseActiveDie() -> m_10_no_action
	@node m_10_no_action = PerformAction("""
										 *Focus* priority:
										 1. Closest to army whose *target* is in a nation at war
										 2. Closest to army whose *target* is in an active nation.
										 3. Closest to a *mobile* army
										 4. Closest to army whose *target* is in a passive nation.
										 5. Random

										 Muster:
										 *Primary*: Elite
										 *Secondary*: Nazgûl

										 If unit is unavailable, rotate as:
										 Elite -> Regular -> Nazgûl -> Elite
										 """) -> m_10_end_no
	@node m_10_end_no = End() -> []

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

	@node muster_card = Start() -> m_c_6
	@node m_c_6 = BinaryCondition("Muster can create an *exposed* region.") -> [n_true = m_c_6_action, n_false = m_c_7]
	@node m_c_6_action = PerformAction("""
									   *Focus* priority:
									   1. Region which creates an *exposed* region
									   2. Random

									   Muster:
									   *Primary*: Units according to card
									   *Secondary*: Units according to card
									   """) -> m_c_6_end
	@node m_c_6_end = End() -> []

	@node m_c_7 = BinaryCondition("""
								  The Fellowship is adjacent to, or in, a region it is possible to muster in.
								  And, the progress put the Fellowship outside Mordor.
								  And, no army is adjacent to the Fellowship's current region.
								  """) -> [n_true = m_c_7_action, n_false = m_c_8]
	@node m_c_7_action = PerformAction("""
									   *Focus* priority:
									   1. The Fellowship's current region

									   Muster:
									   *Primary*: Units according to card
									   *Secondary*: Units according to card
									   """) -> m_c_7_end
	@node m_c_7_end = End() -> []

	@node m_c_8 = BinaryCondition("""
								  Muster is possible in a region containing a Shadow army.
								  """) -> [n_true = m_c_8_action, n_false = m_c_9]
	@node m_c_8_action = PerformAction("""
									   *Focus* priority:
									   1. Army is *mobile*
									   2. Army's *target* is a nation at war
									   3. Army becomes mobile if the Witch King is added
									   4. Free Peoples' army at *target* or on the route to *target* does not contain Gandalf the White
									   5. Army with the highest *value*
									   6. Random

									   Muster:
									   *Primary*: Units according to card
									   *Secondary*: Units according to card
									   """) -> m_c_8_end
	@node m_c_8_end = End() -> []

	@node m_c_9 = PerformAction("""
								*Focus* priority:
								1. Closest to army whose *target* is in a nation at war
								2. Closest to army whose *target* is in an active nation.
								3. Closest to a *mobile* army
								4. Closest to army whose *target* is in a passive nation.
								5. Random

								Muster:
								*Primary*: Units according to card
								*Secondary*: Units according to card
								""") -> m_c_9_end
	@node m_c_9_end = End() -> []

end

max_stars_repo_name: mvmorin/queller-bot

max_issues_repo_name: mvmorin/queller-bot

max_forks_repo_name: mvmorin/queller-bot


text:
using Revise
using Random
using Printf
using StatsBase
using StaticArrays
using POMDPs
using AutoViz
using AutomotiveDrivingModels
using AutomotivePOMDPs
using POMDPSimulators
using BeliefUpdaters
using POMDPPolicies
using DeepQLearning
using Flux
using RLInterface
using POMDPModelTools
# using Interact
# using Blink
using MCTS
using Reel

includet("environment.jl")
includet("generative_mdp.jl")
includet("masking.jl")
includet("cooperative_IDM.jl")
includet("belief_updater.jl")
includet("rendering.jl")
includet("overlays.jl")
includet("make_gif.jl")

function test_state(s,v=0.0, acc=0.0,)
    ego = Vehicle(vehicle_state(35.0, merge_lane(mdp.env), 5.0, mdp.env.roadway), VehicleDef(), EGO_ID)
    veh1 = Vehicle(vehicle_state(s, main_lane(mdp.env), 4.9, mdp.env.roadway), VehicleDef(), EGO_ID + 1)
    scene = Scene()
    push!(scene, ego)
    push!(scene, veh1)
    return AugScene(scene, (acc=acc,))
end

rng = MersenneTwister(1)

mdp = GenerativeMergingMDP(random_n_cars=true, dt=0.5)

mdp.driver_models[2] = CooperativeIDM(c=1.0)
set_desired_speed!(mdp.driver_models[2], 5.0)

s0 = test_state(85.0)

policy = FunctionPolicy(s->7)
for (k,m) in mdp.driver_models
    reset_hidden_state!(m)
end
hr = HistoryRecorder(rng = rng, max_steps=40)
hist = simulate(hr, mdp, policy, s0)

include("visualizer.jl");

mdp.driver_models[2].other_acc = s0.ego_info.acc
observe!(mdp.driver_models[2], s0.scene, mdp.env.roadway, 2)

make_gif(hist, mdp)

generate_s(mdp, s0, 7, rng)

ext: jl


text:
using SatisfiabilityInterface
using Symbolics

different_neighbours(E, c) = [c[i] ≠ c[j] for (i, j) in E]

"Set up a graph colouring problem for a graph
with vertex set `V`, edge set `E` (a vector of edges (i, j)),
with `k` colours"
function graph_colouring_problem(V, E, k=3)

    colours = [:red, :green, :yellow, :blue, :black][1:k]

    # c = [Num(Variable(:c, i)) for i in 1:length(V)]   # colour variables

    @variables c[1:length(V)]

    constraints =
    [
        [c[i] ∈ colours for i in 1:length(V)]

        [c[i] ≠ c[j] for (i, j) in E]  #  different_colours(E, c)
    ]

    return DiscreteCSP(constraints)
end

@testset "Graph colouring" begin

    V = [1, 2, 3]  # vertices
    E = [(1, 2), (2, 3)]  # edges

    prob = graph_colouring_problem(V, E, 2)
    status, results = solve(prob)

    @test status==:sat

    ks = sort(collect(keys(prob.varmap)))
    colours = [results[k] for k in ks]

    @test all(different_neighbours(E, colours))
end

@testset "Ring graph" begin

    function ring_graph(n=11)
        V = 1:n
        E = [(i, mod1((i+1), n)) for i in 1:n]

        return V, E
    end

    V, E = ring_graph(11)

    prob = graph_colouring_problem(V, E, 2)
    status, results = solve(prob)
    @test status==:unsat

    prob = graph_colouring_problem(V, E, 3)
    status, results = solve(prob)
    @test status==:sat

    ks = sort(collect(keys(prob.varmap)))
    colours = [results[k] for k in ks]

    @test all(different_neighbours(E, colours))

end


text:
# # Decaying barotropic QG beta-plane turbulence
#
#md # This example can be run online via [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/barotropicqg_betadecay.ipynb).
#md # Also, it can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/barotropicqg_betadecay.ipynb).
#
# An example of decaying barotropic quasi-geostrophic turbulence on a beta plane.

using FourierFlows, Plots, Printf, Random

using Statistics: mean
using FFTW: irfft

import GeophysicalFlows.BarotropicQG
import GeophysicalFlows.BarotropicQG: energy, enstrophy

# ## Choosing a device: CPU or GPU

dev = CPU()     # Device (CPU/GPU)
nothing # hide

# ## Numerical parameters and time-stepping parameters

      n = 128            # 2D resolution = n^2
stepper = "FilteredRK4"  # timestepper
     dt = 0.05           # timestep
 nsteps = 2000           # total number of time-steps
 nsubs  = 10             # number of time-steps for intermediate logging/plotting (nsteps must be multiple of nsubs)
nothing # hide

# ## Physical parameters

L = 2π        # domain size
β = 10.0      # planetary PV gradient
μ = 0.0       # bottom drag
nothing # hide

# ## Problem setup
# We initialize a `Problem` by providing a set of keyword arguments. Not providing
# a viscosity coefficient ν leads to the module's default value: ν=0. In this
# example numerical instability due to accumulation of enstrophy in high wavenumbers
# is taken care with the `FilteredTimestepper` we picked.
prob = BarotropicQG.Problem(dev; nx=n, Lx=L, β=β, μ=μ, dt=dt, stepper=stepper)
nothing # hide

# and define some shortcuts
sol, cl, vs, pr, gr = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
x, y = gr.x, gr.y
nothing # hide

# ## Setting initial conditions

# Our initial condition consist of a flow that has power only at wavenumbers with
# $8 < \frac{L}{2\pi} \sqrt{k_x^2+k_y^2} < 10$ and initial energy $E_0$:

E0 = 0.1 # energy of initial condition

K = @. sqrt(gr.Krsq)                      # a 2D array with the total wavenumber
k = [gr.kr[i] for i=1:gr.nkr, j=1:gr.nl]  # a 2D array with the zonal wavenumber

Random.seed!(1234)
qih = randn(Complex{eltype(gr)}, size(sol))
qih[K .<  8 * 2π/L] .= 0
qih[K .> 10 * 2π/L] .= 0
qih[k .== 0] .= 0         # no power at zonal wavenumber k=0 component
Ein = energy(qih, gr)     # compute energy of qi
qih *= sqrt(E0 / Ein)     # normalize qi to have energy E0
qi  = irfft(qih, gr.nx)

BarotropicQG.set_zeta!(prob, qi)
nothing #hide

# Let's plot the initial vorticity field:

p1 = heatmap(x, y, vs.q,
         aspectratio = 1,
              c = :balance,
           clim = (-12, 12),
          xlims = (-gr.Lx/2, gr.Lx/2),
          ylims = (-gr.Ly/2, gr.Ly/2),
         xticks = -3:3,
         yticks = -3:3,
         xlabel = "x",
         ylabel = "y",
          title = "initial vorticity ζ=∂v/∂x-∂u/∂y",
     framestyle = :box)

p2 = contourf(x, y, vs.psi,
        aspectratio = 1,
             c = :viridis,
        levels = range(-0.65, stop=0.65, length=10),
          clim = (-0.65, 0.65),
         xlims = (-gr.Lx/2, gr.Lx/2),
         ylims = (-gr.Ly/2, gr.Ly/2),
        xticks = -3:3,
        yticks = -3:3,
        xlabel = "x",
        ylabel = "y",
         title = "initial streamfunction ψ",
    framestyle = :box)

l = @layout grid(1, 2)
p = plot(p1, p2, layout=l, size=(900, 400))

# ## Diagnostics

# Create Diagnostics -- `energy` and `enstrophy` functions are imported at the top.
E = Diagnostic(energy, prob; nsteps=nsteps)
Z = Diagnostic(enstrophy, prob; nsteps=nsteps)
diags = [E, Z] # A list of Diagnostics types passed to "stepforward!" will  be updated every timestep.
nothing # hide

# ## Output

# We choose folder for outputing `.jld2` files and snapshots (`.png` files).
filepath = "."
plotpath = "./plots_decayingbetaturb"
plotname = "snapshots"
filename = joinpath(filepath, "decayingbetaturb.jld2")
nothing # hide

# Do some basic file management,
if isfile(filename); rm(filename); end
if !isdir(plotpath); mkdir(plotpath); end
nothing # hide

# and then create Output.
get_sol(prob) = sol # extracts the Fourier-transformed solution
get_u(prob) = irfft(im*gr.l.*gr.invKrsq.*sol, gr.nx)
out = Output(prob, filename, (:sol, get_sol), (:u, get_u))
nothing # hide

# ## Visualizing the simulation

# We define a function that plots the vorticity and streamfunction fields and
# their corresponding zonal mean structure.

function plot_output(prob)
  ζ = prob.vars.zeta
  ψ = prob.vars.psi
  ζ̄ = mean(ζ, dims=1)'
  ū = mean(prob.vars.u, dims=1)'

  pζ = heatmap(x, y, ζ,
       aspectratio = 1,
            legend = false,
                 c = :balance,
              clim = (-12, 12),
             xlims = (-gr.Lx/2, gr.Lx/2),
             ylims = (-gr.Ly/2, gr.Ly/2),
            xticks = -3:3,
            yticks = -3:3,
            xlabel = "x",
            ylabel = "y",
             title = "vorticity ζ=∂v/∂x-∂u/∂y",
        framestyle = :box)

  pψ = contourf(x, y, ψ,
       aspectratio = 1,
            legend = false,
                 c = :viridis,
            levels = range(-0.65, stop=0.65, length=10),
              clim = (-0.65, 0.65),
             xlims = (-gr.Lx/2, gr.Lx/2),
             ylims = (-gr.Ly/2, gr.Ly/2),
            xticks = -3:3,
            yticks = -3:3,
            xlabel = "x",
            ylabel = "y",
             title = "streamfunction ψ",
        framestyle = :box)

  pζm = plot(ζ̄, y,
            legend = false,
         linewidth = 2,
             alpha = 0.7,
            yticks = -3:3,
             xlims = (-2.2, 2.2),
            xlabel = "zonal mean ζ",
            ylabel = "y")
  plot!(pζm, 0*y, y, linestyle=:dash, linecolor=:black)

  pum = plot(ū, y,
            legend = false,
         linewidth = 2,
             alpha = 0.7,
            yticks = -3:3,
             xlims = (-0.55, 0.55),
            xlabel = "zonal mean u",
            ylabel = "y")
  plot!(pum, 0*y, y, linestyle=:dash, linecolor=:black)

  l = @layout grid(2, 2)
  p = plot(pζ, pζm, pψ, pum, layout = l, size = (900, 800))

  return p
end
nothing # hide

# ## Time-stepping the `Problem` forward

# We time-step the `Problem` forward in time.

startwalltime = time()

p = plot_output(prob)

anim = @animate for j=0:Int(nsteps/nsubs)

  log = @sprintf("step: %04d, t: %d, E: %.4f, Q: %.4f, walltime: %.2f min",
    cl.step, cl.t, E.data[E.i], Z.data[Z.i], (time()-startwalltime)/60)

  if j%(1000/nsubs)==0; println(log) end

  p[1][1][:z] = Array(vs.zeta)
  p[1][:title] = "vorticity, t="*@sprintf("%.2f", cl.t)
  p[3][1][:z] = Array(vs.psi)
  p[2][1][:x] = mean(vs.zeta, dims=1)'
  p[4][1][:x] = mean(vs.u, dims=1)'

  stepforward!(prob, diags, nsubs)
  BarotropicQG.updatevars!(prob)

end

mp4(anim, "barotropicqg_betadecay.mp4", fps=14)

# ## Save

# Finally save the last snapshot.
savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), cl.step)
savefig(savename)


text:
using LinearAlgebra: isapprox
using LazyLieManifolds
using Test
using StaticArrays
using LinearAlgebra

import LazyLieManifolds: ⊕

using Manifolds

@testset "LazyLieManifolds.jl" begin

    _p0 = 1.0I(3)
    _X1 = [0.0 -pi/2 1.0;  pi/2 0.0 0.0; 0.0 0.0 1.0]
    _X2 = [0.0  pi/2 1.0; -pi/2 0.0 0.0; 0.0 0.0 1.0]

    p0 = SE([1.0 0.0; 0 1.0], [0.0, 0.0])

    M = LLM.manifold(p0)
    @test all(isapprox.(_p0, affine_matrix(M, p0[])))

    X1 = se{2}([1., 0, pi/2])
    X2 = se{2}([1., 0, -pi/2])

    @test all(isapprox.(_X1, affine_matrix(M, X1[])))
    @test all(isapprox.(_X2, affine_matrix(M, X2[])))

    p1 = p0 ⊕ X1

    p2 = p1 ⊕ X2

    p3 = p1 ∘ p2

    @time LLM.PosePose(p1, p0, p2)
    @time LLM.PosePose_x(p1, p0, p2)

    SE{3}()

    p0 = LLM._identity(SE{3})
    X1 = se{3}([1., 0, 0, pi/2, 0, 0])
    X2 = se{3}([1., 0, 0, -pi/2, 0, 0])
    p1 = p0 ⊕ X1
    p2 = p1 ⊕ X2

    # @code_warntype exp(X1)
    log(p2)
    log(p1, p2)

    ##

    p0 = SO([1.0 0.0; 0 1.0])
    X1 = so{2}([pi/2])
    X2 = se{2}([1., 0, -pi/2])

    # p1 = p0 ⊕ X1
    # p2 = p1 ⊕ X2
    # p3 = p1 ∘ p2

    SO{5}()

    p0 = LLM._identity(SO{3})
    X1 = so{3}([pi/2, 0, 0])
    X2 = so{3}([-pi/2, 0, 0])

    log(p2)
    log(p1, p2)

end

##

# Compare SArray to MArray
# tested on small function as
# M: (5 allocations: 176 bytes)
# S: (3 allocations: 112 bytes)
##################
p0 = SE{2}(ProductRepr((@MVector([0.0, 0.0]), @MMatrix([1.0 0.0; 0.0 1.0]))))

X1 = se{2}(hat(M, p0[], @MVector([1., 0.0, pi/2])))
X1 = se{2}(hat(M, p0[], @MVector([1., 0.0, -pi/2])))

p1 = p0 ⊕ X1
p2 = p1 ⊕ X2
p3 = p1 ∘ p2

@time LLM.PosePose(p1, p0, p2)
@time LLM.PosePose_x(p1, p0, p2)
@time LLM.PosePose_d(p1, p0, p2)

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

p0 = SE{2}(ProductRepr((@SVector([0.0, 0.0]), @SMatrix([1.0 0.0; 0.0 1.0]))))

X1 = se{2}(hat(M, p0[], @SVector([1., 0.0, pi/2])))
X1 = se{2}(hat(M, p0[], @SVector([1., 0.0, -pi/2])))

p1 = p0 ⊕ X1
p2 = p1 ⊕ X2
p3 = p1 ∘ p2

@time LLM.PosePose(p1, p0, p2)
@time LLM.PosePose_x(p1, p0, p2)
@time LLM.PosePose_d(p1, p0, p2)

ProductRepr((@SVector([0.0, 0.0]), @SMatrix([1.0 0.0; 0.0 1.0])))
ProductRepr((@MVector([0.0, 0.0]), @MMatrix([1.0 0.0; 0.0 1.0])))


text:
module BetheBloch

export Banana

#include("material.jl")
#include("particle.jl")
include("banana.jl")

end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# essential type definitions for interacting with C code
# (platform- or OS-dependent types are defined in c.jl)

"""
    Cuchar

Equivalent to the native `unsigned char` c-type (`UInt8`).
"""
const Cuchar = UInt8

"""
    Cshort

Equivalent to the native `signed short` c-type (`Int16`).
"""
const Cshort = Int16

"""
    Cushort

Equivalent to the native `unsigned short` c-type (`UInt16`).
"""
const Cushort = UInt16

"""
    Cint

Equivalent to the native `signed int` c-type (`Int32`).
"""
const Cint = Int32

"""
    Cuint

Equivalent to the native `unsigned int` c-type (`UInt32`).
"""
const Cuint = UInt32

"""
    Cptrdiff_t

Equivalent to the native `ptrdiff_t` c-type (`Int`).
"""
const Cptrdiff_t = Int

"""
    Csize_t

Equivalent to the native `size_t` c-type (`UInt`).
"""
const Csize_t = UInt

"""
    Cssize_t

Equivalent to the native `ssize_t` c-type.
"""
const Cssize_t = Int

"""
    Cintmax_t

Equivalent to the native `intmax_t` c-type (`Int64`).
"""
const Cintmax_t = Int64

"""
    Cuintmax_t

Equivalent to the native `uintmax_t` c-type (`UInt64`).
"""
const Cuintmax_t = UInt64

"""
    Clonglong

Equivalent to the native `signed long long` c-type (`Int64`).
"""
const Clonglong = Int64

"""
    Culonglong

Equivalent to the native `unsigned long long` c-type (`UInt64`).
"""
const Culonglong = UInt64

"""
    Cfloat

Equivalent to the native `float` c-type (`Float32`).
"""
const Cfloat = Float32

"""
    Cdouble

Equivalent to the native `double` c-type (`Float64`).
"""
const Cdouble = Float64


text:
module Classification

using LinearAlgebra
export SVC, Logistic

include("utils.jl")
include("Logistic.jl")
include("SVC.jl")
end


text:
mutable struct sem_mesh_data
    nspec::Int64
    nglob::Int64
    xyz_glob::Array{Float64,2}
    ibool::Array{Int64,4}
    idoubling::Array{Int64,1}
    ispec_is_tiso::Array{Bool,1}

    xix::Array{Float64,4}
    xiy::Array{Float64,4}
    xiz::Array{Float64,4}
    etax::Array{Float64,4}
    etay::Array{Float64,4}
    etaz::Array{Float64,4}
    gammax::Array{Float64,4}
    gammay::Array{Float64,4}
    gammaz::Array{Float64,4}

    sem_mesh_data() = new()
end

mutable struct sem_mesh_location
    stat::Int64
    eid::Int64

    uvw::Vector{Float64}
    misloc::Float64
    lagrange::Array{Float64,3}

    sem_mesh_location() = new()
end

max_stars_repo_name: ziyixi/SeisScripts

max_issues_repo_name: ziyixi/SeisScripts

max_forks_repo_name: ziyixi/SeisScripts


text:
function holtslag_local(Ri)
  if Ri>0
    y=1/(1 + 10 * Ri * (1 + 8 * Ri))
  else
    y=(1 -18 * Ri)^(0.5)
  end
  return y
end


text:
# load ttmtdb as a TaylorInterpolant saved in .jld file
const ttmtdb = JLD.load(joinpath(artifact"ttmtdb_DE430_1995_2030", "ttmtdb_DE430_1995_2030.jld"), "ttmtdb")

# read JPL ephemerides (asteroid, Solar System, TT-TDB)
function loadjpleph()
    furnsh(
        joinpath(artifact"naif0012", "naif0012.tls"),
        joinpath(artifact"TTmTDBde430", "TTmTDB.de430.19feb2015.bsp"),
        joinpath(artifact"de430", "de430_1850-2150.bsp"),
        joinpath(artifact"a99942", "a99942_s197.bsp"),
        joinpath(artifact"a99942", "a99942_s199.bsp"),
    )
end

# this is an auxiliary function which converts a [x,y,z,vx,vy,vz] "state" vector from km,km/sec units to au,au/day
function kmsec2auday(pv)
    pv /= au # (km, km/sec) -> (au, au/sec)
    pv[4:6] *= daysec # (au, au/sec) -> (au, au/day)
    return pv
end

# this is an auxiliary function which converts a [x,y,z,vx,vy,vz] "state" vector from au,au/day units to km,km/sec
function auday2kmsec(pv)
    pv *= au # (au, au/day) -> (km, km/day)
    pv[4:6] /= daysec # (km, km/day) -> (km, km/sec)
    return pv
end

# get [x,y,z,vx,vy,vz] geometric "state" vector at TDB instant `et` from
# SPK-formatted ephemeris file wrt J2000 frame
function getpv(target::Int, observer::Int, et)
    return spkgeo(target, et, "J2000", observer)[1] # units: km,km/sec
end

# NAIF IDs:
#0: Solar System Barycenter
#10: Sun (heliocenter)
#2099942: Apophis
#399: Earth (geocenter)
#301: Moon
#1000000001 from body 1000000000: TT-TDB
# Here, we follow the convention from the CSPICE, library, that the ephemeris
# time is referred to the J2000 frame epoch:
# https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/spk.html#Terminology
# argument `et` represents "ephemeris seconds" (TDB seconds) since J2000.0 TDB epoch
# position and velocity are assumed to be returned in km, km/sec, resp., by spkgeo
# https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/spkgeo_c.html (see: Detailed ouput section)
apophis_pv_197(et) = getpv(9904406, 0, constant_term(et)) # units: km, km/second
apophis_pv_199(et) = getpv(2099942, 0, constant_term(et)) # units: km, km/second
sun_pv(et) = getpv(10, 0, constant_term(et)) # units: km, km/second
earth_pv(et) = getpv(399, 0, constant_term(et)) # units: km, km/second
moon_pv(et) = getpv(301, 0, constant_term(et)) # units: km, km/second
tt_tdb(et) = getpv(1000000001, 1000000000, constant_term(et))[1] # units: seconds
dtt_tdb(et) = getpv(1000000001, 1000000000, constant_term(et))[4] # units: seconds/seconds

# Convert julian days to ephemeris seconds since J2000
function julian2etsecs(jd)
    return (jd-JD_J2000)*daysec
end

# Convert ephemeris seconds since J2000 to julian days
function etsecs2julian(et)
    return JD_J2000 + et/daysec
end

# Standard formula for relativistic (Shapiro) delay
function shapiro_delay(e, p, q)
    shap = 0.0 #2μ[1]/(c_au_per_day^2)
    shap_del_days = (2PlanetaryEphemeris.μ[su]/(c_au_per_day^3))*log( (e+p+q+shap)/(e+p-q+shap) ) # days
    return shap_del_days*daysec # seconds
end

function shapiro_doppler(e, de, p, dp, q, dq, F_tx)
    # shap_del_diff = 2μ[1]*( (de+dp+dq)/(e+p+q) - (de+dp-dq)/(e+p-q) )/(c_au_per_day^3) # (adim.)
    shap_del_diff = (4PlanetaryEphemeris.μ[su]/(c_au_per_day^3))*(  ( dq*(e+p) - q*(de+dp) )/( (e+p)^2 - q^2 )  ) # differential of Shapiro delay (adim.)
    shap_dop = -F_tx*shap_del_diff # ν = -F_tx*dτ/dt (units of F_tx) <-- Shapiro, Ash, Tausner (1966), footnote 10
    return shap_dop # (units of F_tx)
end

# Density of ionized electrons in interplanetary medium (ESAA 2014, p. 323, Sec. 8.7.5, Eq. 8.22)
# ESAA 2014 in turn refers to Muhleman and Anderson (1981)
# Ostro (1993) gives a reference to Anderson (1978), where this model is fitted to Mariner 9 ranging data
# Reading https://gssc.esa.int/navipedia/index.php/Ionospheric_Delay
# Helped a lot to clarify things, especially the 40.3, although they talk about Earth's ionosphere
# Another valuable source is Standish, E.M., Astron. Astrophys. 233, 252-271 (1990)
# Ne: ionized electrons density: electrons/cm^3
# p1: signal departure point (transmitter/bounce) (au)
# p2: signal arrival point (bounce/receiver) (au)
# r_s_t0: Barycentric position (au) of Sun at initial time of propagation of signal path (bounce time for down-leg; transmit time for up-leg)
# ds: current distance travelled by ray from emission point (au)
# ΔS: total distance between p1 and p2 (au)
function Ne(p1::Vector{S}, p2::Vector{S}, r_s_t0::Vector{S}, ds::U, ΔS::Real) where {S<:Number, U<:Number}
    # s: linear parametrization of ray path, such that
    # s=0 -> point on ray path is at p1; s=1 -> point on ray path is at p2
    s = ds/ΔS
    s_p2_p1 = map(x->s*x, Taylor1.(p2-p1, s.order))
    r_vec = Taylor1.(p1, s.order) + s_p2_p1 - Taylor1.(r_s_t0, s.order) # heliocentric position (au) of point on ray path at time t_tdb_jul (Julian days)
    r = sqrt( r_vec[1]^2 + r_vec[2]^2 + r_vec[3]^2 ) # heliocentric distance (au) of point on ray path at time t_tdb_jul (Julian days)
    # compute heliocentric position vector of point on ray path wrt Sun's rotation pole and equator (i.e., heliographic)
    α_p_sun_rad = deg2rad(α_p_sun)
    δ_p_sun_rad = deg2rad(δ_p_sun)
    r_vec_heliographic = inv( pole_rotation(α_p_sun_rad, δ_p_sun_rad) )*r_vec
    # compute heliographic solar latitude (Anderson, 1978) of point on ray path
    β = asin( r_vec_heliographic[3]/r ) # ecliptic solar latitude of point on ray path (rad)
    r_sr = r/R_sun
    Ne_t1 = (A_sun/r_sr^6)
    Ne_t2 = ( (a_sun*b_sun)/sqrt((a_sun*sin(β))^2 + (b_sun*cos(β))^2) )/(r_sr^2)
    Ne_val = Ne_t1 + Ne_t2
    return Ne_val
end

# Integral of Ne, evaluated with TaylorIntegration.jl
# TODO: @taylorize!
# p1: signal departure point (transmitter/bounce) (au)
# p2: signal arrival point (bounce/receiver) (au)
# r_s_t0: Barycentric position (au) of Sun at initial time of propagation of signal path (bounce time for down-leg; transmit time for up-leg)
# output is in units of electrons/cm^2
function Ne_path_integral(p1::Vector{S}, p2::Vector{S}, r_s_t0::Vector{S}) where {S<:Number}
    ΔS = (100_000au)*norm(p2-p1) # total distance between p1 and p2, in centimeters
    # kernel of path integral; distance parameter `s` and total distance `ΔS` is in cm
    function int_kernel(x, params, s)
        return Ne(p1, p2, r_s_t0, s, ΔS)
    end
    # do path integral
    i0 = zero(p1[1])
    tT = Taylor1(24)
    iT = Taylor1(i0, 24)
    TaylorIntegration.jetcoeffs!(int_kernel, tT, iT, nothing)
    return iT(ΔS)
end

# Time-delay due to thin plasma of solar corona
# Ne: ionized electrons density: electrons/cm^3
# p1: signal departure point (transmitter/bounce for up/down-link, resp.) (au)
# p2: signal arrival point (bounce/receiver for up/down-link, resp.) (au)
# t_tdb_jd1, t_tdb_jd2: Two-part Julian date (TDB) of signal path (bounce time for downlink; transmit time for uplink)
# F_tx: transmitter frequency (MHz)
# From https://gssc.esa.int/navipedia/index.php/Ionospheric_Delay it seems that
# ESAA 2014 text probably should say that in the formula for Δτ_corona,
# the expression 40.3*Ne/f^2 is adimensional, where Ne [electrons/centimeters^3] and f is in Hz
# therefore, the integral (40.3/f^2)*∫Ne*ds is in centimeters, where ds is in centimeters
# and the expression (40.3/(c*f^2))*∫Ne*ds, with c in centimeters/second, is in seconds
function corona_delay(p1::Vector{S}, p2::Vector{S}, r_s_t0::Vector{S}, F_tx::U, station_code::Int) where {S<:Number, U<:Real}
    # for the time being, we're removing the terms associated with higher-order terms in the variationals (ie, Yarkovsky)
    int_path = Ne_path_integral(p1, p2, r_s_t0) # (electrons/cm^2)
    Δτ_corona = 40.3e-6int_path/(c_cm_per_sec*(F_tx)^2) # seconds
    return Δτ_corona # seconds
end

# *VERY* elementary computation of zenith distance
# r_antenna: position of antenna at receive/transmit time in celestial frame wrt geocenter
# ρ_vec_ae: slant-range vector from antenna to asteroid
function zenith_distance(r_antenna::Vector{T}, ρ_vec_ae::Vector{S}) where {T<:Number, S<:Number}
    norm_r_antenna = sqrt(r_antenna[1]^2 + r_antenna[2]^2 + r_antenna[3]^2)
    norm_ρ_vec_ae = sqrt(ρ_vec_ae[1]^2 + ρ_vec_ae[2]^2 + ρ_vec_ae[3]^2)
    cos_antenna_slant = dot(r_antenna, ρ_vec_ae)/(norm_r_antenna*norm_ρ_vec_ae)
    return acos(cos_antenna_slant)
end

# Time delay from the Earth's troposphere for radio frequencies
# oscillates between 0.007μs and 0.225μs for z between 0 and π/2 rad
# z: zenith distance (radians)
tropo_delay(z) = (7e-9)/( cos(z) + 0.0014/(0.045+cot(z)) ) # seconds

function tropo_delay(r_antenna::Vector{T}, ρ_vec_ae::Vector{S}) where {T<:Number, S<:Number}
    zd = zenith_distance(r_antenna, ρ_vec_ae)
    return tropo_delay(zd) # seconds
end

# auxiliary function to compute (TDB-UTC); useful to convert TDB to UTC via UTC + (TDB-UTC) and viceversa
# TDB: Solar System barycentric ephemeris time
# et: TDB seconds since J2000.0 (TDB)
# TT: Terrestrial time
# TAI: International Atomic Time
# UTC: Coordinated Universal Time
# TDB-UTC = (TDB-TAI) + (TAI-UTC)
#         = (TDB-TT) + (TT-TAI) + (TAI-UTC)
#         = (TDB-TT) + (32.184 s) + ΔAT
# does not include correction due to position of measurement station v_E*(r_S.r_E)/c^2 (Folkner et al. 2014; Moyer, 2003)
function tdb_utc(et::T) where {T<:Number}
    tt_tdb_et = ttmtdb(et)
    tt_tai = 32.184
    et_00 = constant_term(constant_term(et))
    utc_secs = et_00 - deltet(et_00, "ET") # used only to determine ΔAT; no high-precision needed
    jd_utc = JD_J2000 + utc_secs/daysec
    tai_utc = get_ΔAT(jd_utc)
    return (tt_tai + tai_utc) - tt_tdb_et # TDB-UTC = (TDB-TT) + (TT-TAI) + (TAI-UTC) = (TDB-TT) + 32.184 s + ΔAT
end

# TODO: add tdb_utc(utc) method!!!
# strategy: given UTC, do UTC + (TT-TAI) + (TAI-UTC) to get TT
# then, use ttmtdb(et) function iteratively (Newton) to compute TDB

# function dtutc2et(t_utc::DateTime)
#     tt_tai = 32.184
#     jd_utc = datetime2julian(t_utc)
#     fd_utc = (jd_utc+0.5) - floor(jd_utc+0.5)
#     j, tai_utc = iauDat(year(t_utc), month(t_utc), day(t_utc), fd_utc)
#     return et
# end

# Compute radar-astrometric round-trip time for an asteroid at
# UTC instant `t_r_utc` from tracking station with code `station_code`.
# station_code: observing station identifier (MPC nomenclature)
# t_r_utc: UTC time of echo reception (DateTime)
# t_offset: time offset wrt echo reception time, to compute Doppler shifts by range differences (seconds)
# niter: number of light-time solution iterations
# xve: Earth ephemeris wich takes TDB seconds since J2000 as input and returns Earth barycentric position in km and velocity in km/second
# xvs: Sun ephemeris wich takes TDB seconds since J2000 as input and returns Sun barycentric position in km and velocity in km/second
# xva: asteroid ephemeris wich takes TDB seconds since J2000 as input and returns asteroid barycentric position in km and velocity in km/second
function delay(station_code::Int, t_r_utc::DateTime, t_offset::Real,
        niter::Int=10; eo::Bool=true, xve=earth_pv, xvs=sun_pv,
        xva=apophis_pv_197)
    et_r_secs = str2et(string(t_r_utc)) + t_offset
    # Compute geocentric position/velocity of receiving antenna in inertial frame (au, au/day)
    R_r, V_r = observer_position(station_code, et_r_secs, eo=eo)
    # Earth's barycentric position and velocity at receive time
    rv_e_t_r = xve(et_r_secs)
    r_e_t_r = rv_e_t_r[1:3]
    v_e_t_r = rv_e_t_r[4:6]
    # Receiver barycentric position and velocity at receive time
    r_r_t_r = r_e_t_r + R_r
    # Asteroid barycentric position and velocity at receive time
    rv_a_t_r = xva(et_r_secs)
    r_a_t_r = rv_a_t_r[1:3]
    # Sun barycentric position and velocity at receive time
    rv_s_t_r = xvs(et_r_secs)
    r_s_t_r = rv_s_t_r[1:3]
    # down-leg iteration
    # τ_D first approximation: Eq. (1) Yeomans et al. (1992)
    ρ_vec_r = r_a_t_r - r_r_t_r
    ρ_r = sqrt(ρ_vec_r[1]^2 + ρ_vec_r[2]^2 + ρ_vec_r[3]^2)
    τ_D = ρ_r/clightkms # (seconds) -R_b/c, but delay is wrt asteroid Center (Brozovic et al., 2018)
    # bounce time, new estimate Eq. (2) Yeomans et al. (1992)
    et_b_secs = et_r_secs - τ_D

    Δτ_D = zero(τ_D)
    Δτ_rel_D = zero(τ_D)
    # Δτ_corona_D = zero(τ_D)
    Δτ_tropo_D = zero(τ_D)

    for i in 1:niter
        # asteroid barycentric position (in au) at bounce time (TDB)
        rv_a_t_b = xva(et_b_secs)
        r_a_t_b = rv_a_t_b[1:3]
        v_a_t_b = rv_a_t_b[4:6]
        # Eq. (3) Yeomans et al. (1992)
        ρ_vec_r = r_a_t_b - r_r_t_r
        # Eq. (4) Yeomans et al. (1992)
        ρ_r = sqrt(ρ_vec_r[1]^2 + ρ_vec_r[2]^2 + ρ_vec_r[3]^2)
        # compute down-leg Shapiro delay
        # NOTE: when using PPN, substitute 2 -> 1+γ in expressions for Shapiro delay, Δτ_rel_[D|U]
        e_D_vec  = r_r_t_r - r_s_t_r
        e_D = sqrt(e_D_vec[1]^2 + e_D_vec[2]^2 + e_D_vec[3]^2) # heliocentric distance of Earth at t_r
        rv_s_t_b = xvs(et_b_secs) # barycentric position and velocity of Sun at estimated bounce time
        r_s_t_b = rv_s_t_b[1:3]
        p_D_vec  = r_a_t_b - r_s_t_b
        p_D = sqrt(p_D_vec[1]^2 + p_D_vec[2]^2 + p_D_vec[3]^2) # heliocentric distance of asteroid at t_b
        q_D = ρ_r #signal path distance (down-leg)
        # Shapiro correction to time-delay
        Δτ_rel_D = shapiro_delay(e_D, p_D, q_D)
        # troposphere correction to time-delay
        Δτ_tropo_D = tropo_delay(R_r, ρ_vec_r) # seconds
        # Δτ_corona_D = corona_delay(constant_term.(r_a_t_b), r_r_t_r, r_s_t_r, F_tx, station_code) # seconds
        Δτ_D = Δτ_rel_D # + Δτ_tropo_D #+ Δτ_corona_D # seconds
        p_dot_23 = dot(ρ_vec_r, v_a_t_b)/ρ_r
        Δt_2 = (τ_D - ρ_r/clightkms - Δτ_rel_D)/(1.0-p_dot_23/clightkms)
        τ_D = τ_D - Δt_2
        et_b_secs = et_r_secs - τ_D
    end
    rv_a_t_b = xva(et_b_secs)
    r_a_t_b = rv_a_t_b[1:3]
    v_a_t_b = rv_a_t_b[4:6]

    # up-leg iteration
    # τ_U first estimation: Eq. (5) Yeomans et al. (1992)
    τ_U = τ_D
    # transmit time, 1st estimate Eq. (6) Yeomans et al. (1992)
    et_t_secs = et_b_secs - τ_U
    # Geocentric position and velocity of transmitting antenna in inertial frame (au, au/day)
    R_t, V_t = observer_position(station_code, et_t_secs, eo=eo)
    rv_e_t_t = xve(et_t_secs)
    r_e_t_t = rv_e_t_t[1:3]
    v_e_t_t = rv_e_t_t[4:6]
    # Transmitter barycentric position and velocity of at transmit time
    r_t_t_t = r_e_t_t + R_t
    # Eq. (7) Yeomans et al. (1992)
    ρ_vec_t = r_a_t_b - r_t_t_t
    ρ_t = sqrt(ρ_vec_t[1]^2 + ρ_vec_t[2]^2 + ρ_vec_t[3]^2)

    Δτ_U = zero(τ_U)
    Δτ_rel_U = zero(τ_U)
    # Δτ_corona_U = zero(τ_U)
    Δτ_tropo_U = zero(τ_U)

    # println("   *** DOWNLEG LOOP ***")
    for i in 1:niter
        # Geocentric position and velocity of transmitting antenna in inertial frame (au, au/day)
        # TODO: remove `constant_term` to take into account dependency of R_t, V_t wrt initial conditions variations via et_t_secs
        R_t, V_t = observer_position(station_code, et_t_secs, eo=eo)
        # Earth's barycentric position and velocity at the transmit time
        rv_e_t_t = xve(et_t_secs)
        r_e_t_t = rv_e_t_t[1:3]
        v_e_t_t = rv_e_t_t[4:6]
        # Barycentric position and velocity of the transmitter at the transmit time
        r_t_t_t = r_e_t_t + R_t
        v_t_t_t = v_e_t_t + V_t
        # Eq. (7) Yeomans et al. (1992)
        ρ_vec_t = r_a_t_b - r_t_t_t
        ρ_t = sqrt(ρ_vec_t[1]^2 + ρ_vec_t[2]^2 + ρ_vec_t[3]^2)
        # compute up-leg Shapiro delay
        # Sun barycentric position and velocity (in au, au/day) at transmit time (TDB)
        rv_s_t_t = xvs(et_t_secs)
        r_s_t_t = rv_s_t_t[1:3]
        e_U_vec = r_t_t_t - r_s_t_t
        e_U = sqrt(e_U_vec[1]^2 + e_U_vec[2]^2 + e_U_vec[3]^2) # heliocentric distance of Earth at t_t
        rv_s_t_b = xvs(et_b_secs) # barycentric position/velocity of Sun at bounce time
        r_s_t_b = rv_s_t_b[1:3]
        p_U_vec = r_a_t_b - r_s_t_b
        p_U = sqrt(p_U_vec[1]^2 + p_U_vec[2]^2 + p_U_vec[3]^2) # heliocentric distance of asteroid at t_b
        q_U_vec = r_a_t_b - r_e_t_t
        q_U = ρ_t # signal path distance (up-leg)
        Δτ_rel_U = shapiro_delay(e_U, p_U, q_U) # seconds
        Δτ_tropo_U = tropo_delay(R_t, ρ_vec_t) # seconds
        # Δτ_corona_U = corona_delay(constant_term.(r_t_t_t), constant_term.(r_a_t_b), constant_term.(r_s_t_b), F_tx, station_code) # seconds
        Δτ_U = Δτ_rel_U # + Δτ_tropo_U #+ Δτ_corona_U # seconds
        p_dot_12 = -dot(ρ_vec_t, v_t_t_t)/ρ_t
        Δt_1 = (τ_U - ρ_t/clightkms - Δτ_rel_U)/(1.0-p_dot_12/clightkms)
        τ_U = τ_U - Δt_1
        # transmit time, new estimate
        et_t_secs = et_b_secs - τ_U
    end

    # compute TDB-UTC at transmit time
    # corrections to TT-TDB from Moyer (2003) / Folkner et al. (2014) due to position of measurement station on Earth are of order 0.01μs
    # Δtt_tdb_station_t = - dot(v_e_t_t, r_t_t_t-r_e_t_t)/clightkms^2
    tdb_utc_t = tdb_utc(et_t_secs) # + Δtt_tdb_station_t
    # compute TDB-UTC at receive time
    # corrections to TT-TDB from Moyer (2003) / Folkner et al. (2014) due to position of measurement station on Earth  are of order 0.01μs
    # Δtt_tdb_station_r = - dot(v_e_t_r, r_r_t_r-r_e_t_r)/clightkms^2
    tdb_utc_r = tdb_utc(et_r_secs) # + Δtt_tdb_station_r

    # compute total time delay (UTC seconds); relativistic delay is already included in τ_D, τ_U
    # Eq. (9) Yeomans et al. (1992)
    τ = (τ_D + τ_U) + (Δτ_tropo_D + Δτ_tropo_U) + (tdb_utc_t - tdb_utc_r) # seconds

    return 1e6τ # total signal delay (μs)
end

# Compute Taylor series expansion of time-delay observable around echo reception
# time. This allows to compute dopplers via autodiff using ν = -F_tx*dτ/dt
# NOTE: Works only with TaylorInterpolant ephemeris
function delay(station_code::Int, t_r_utc::DateTime,
        niter::Int=10; eo::Bool=true, xve::TaylorInterpolant=earth_pv,
        xvs::TaylorInterpolant=sun_pv, xva::TaylorInterpolant=apophis_pv_197,
        tord::Int=xva.x[1].order)
    q1 = xva.x[1] # auxiliary to evaluate JT ephemeris
    et_r_secs_0 = str2et(string(t_r_utc))
    et_r_secs = Taylor1([et_r_secs_0,1.0].*one(q1[0]), tord)
    # Compute geocentric position/velocity of receiving antenna in inertial frame (au, au/day)
    R_r, _ = observer_position(station_code, et_r_secs, eo=eo)
    # Earth's barycentric position and velocity at receive time
    r_e_t_r = xve(et_r_secs)[1:3]
    # Receiver barycentric position and velocity at receive time
    r_r_t_r = r_e_t_r + R_r
    # Asteroid barycentric position and velocity at receive time
    r_a_t_r = xva(et_r_secs)[1:3]
    # Sun barycentric position and velocity at receive time
    r_s_t_r = xvs(et_r_secs)[1:3]
    # down-leg iteration
    # τ_D first approximation: Eq. (1) Yeomans et al. (1992)
    ρ_vec_r = r_a_t_r - r_r_t_r
    ρ_r = sqrt(ρ_vec_r[1]^2 + ρ_vec_r[2]^2 + ρ_vec_r[3]^2)
    # @show ρ_r
    τ_D = ρ_r/clightkms # (seconds) -R_b/c, but delay is wrt asteroid Center (Brozovic et al., 2018)
    # @show τ_D
    # bounce time, new estimate Eq. (2) Yeomans et al. (1992)
    # @show et_r_secs
    et_b_secs = et_r_secs - τ_D
    # @show et_b_secs et_r_secs

    Δτ_D = zero(τ_D)
    Δτ_rel_D = zero(τ_D)
    # Δτ_corona_D = zero(τ_D)
    Δτ_tropo_D = zero(τ_D)
    # @show Δτ_D

    for i in 1:niter
        # asteroid barycentric position (in au) at bounce time (TDB)
        rv_a_t_b = xva(et_b_secs)
        r_a_t_b = rv_a_t_b[1:3]
        v_a_t_b = rv_a_t_b[4:6]
        # Eq. (3) Yeomans et al. (1992)
        ρ_vec_r = r_a_t_b - r_r_t_r
        # Eq. (4) Yeomans et al. (1992)
        ρ_r = sqrt(ρ_vec_r[1]^2 + ρ_vec_r[2]^2 + ρ_vec_r[3]^2)
        # @show ρ_r
        # compute down-leg Shapiro delay
        # NOTE: when using PPN, substitute 2 -> 1+γ in expressions for Shapiro delay, Δτ_rel_[D|U]
        e_D_vec  = r_r_t_r - r_s_t_r
        # @show e_D_vec
        e_D = sqrt(e_D_vec[1]^2 + e_D_vec[2]^2 + e_D_vec[3]^2) # heliocentric distance of Earth at t_r
        # barycentric position and velocity of Sun at estimated bounce time
        r_s_t_b = xvs(et_b_secs)[1:3]
        # @show r_s_t_b
        p_D_vec  = r_a_t_b - r_s_t_b
        p_D = sqrt(p_D_vec[1]^2 + p_D_vec[2]^2 + p_D_vec[3]^2) # heliocentric distance of asteroid at t_b
        q_D = ρ_r #signal path distance (down-leg)
        # Shapiro correction to time-delay
        Δτ_rel_D = shapiro_delay(e_D, p_D, q_D)
        # @show Δτ_rel_D
        # troposphere correction to time-delay
        Δτ_tropo_D = tropo_delay(R_r, ρ_vec_r) # seconds
        # @show Δτ_tropo_D
        # Δτ_corona_D = corona_delay(constant_term.(r_a_t_b), r_r_t_r, r_s_t_r, F_tx, station_code) # seconds
        Δτ_D = Δτ_rel_D # + Δτ_tropo_D #+ Δτ_corona_D # seconds
        # @show Δτ_D
        p_dot_23 = dot(ρ_vec_r, v_a_t_b)/ρ_r
        # @show p_dot_23
        Δt_2 = (τ_D - ρ_r/clightkms - Δτ_rel_D)/(1.0-p_dot_23/clightkms)
        # @show Δt_2
        τ_D = τ_D - Δt_2
        et_b_secs = et_r_secs - τ_D
    end
    # @show τ_D
    rv_a_t_b = xva(et_b_secs)
    r_a_t_b = rv_a_t_b[1:3]
    v_a_t_b = rv_a_t_b[4:6]

    # up-leg iteration
    # τ_U first estimation: Eq. (5) Yeomans et al. (1992)
    τ_U = τ_D
    # transmit time, 1st estimate Eq. (6) Yeomans et al. (1992)
    et_t_secs = et_b_secs - τ_U
    # @show τ_U et_t_secs et_b_secs et_r_secs
    # Geocentric position and velocity of transmitting antenna in inertial frame (au, au/day)
    R_t, V_t = observer_position(station_code, et_t_secs, eo=eo)
    # @show R_t V_t
    # Earth's barycentric position and velocity at transmit time
    rv_e_t_t = xve(et_t_secs)
    r_e_t_t = rv_e_t_t[1:3]
    v_e_t_t = rv_e_t_t[4:6]
    # @show r_e_t_t norm(r_e_t_t)/au norm(v_e_t_t)
    # Transmitter barycentric position and velocity of at transmit time
    r_t_t_t = r_e_t_t + R_t
    # # Eq. (7) Yeomans et al. (1992)
    ρ_vec_t = r_a_t_b - r_t_t_t
    ρ_t = sqrt(ρ_vec_t[1]^2 + ρ_vec_t[2]^2 + ρ_vec_t[3]^2)
    # @show ρ_t

    Δτ_U = zero(τ_U)
    Δτ_rel_U = zero(τ_U)
    # Δτ_corona_U = zero(τ_U)
    Δτ_tropo_U = zero(τ_U)
    # @show Δτ_U

    # println("   *** DOWNLEG LOOP ***")
    for i in 1:niter
        # Geocentric position and velocity of transmitting antenna in inertial frame (au, au/day)
        R_t, V_t = observer_position(station_code, et_t_secs, eo=eo)
        # Earth's barycentric position and velocity at transmit time
        rv_e_t_t = xve(et_t_secs)
        r_e_t_t = rv_e_t_t[1:3]
        v_e_t_t = rv_e_t_t[4:6]
        # Barycentric position and velocity of the transmitter at the transmit time
        r_t_t_t = r_e_t_t + R_t
        v_t_t_t = v_e_t_t + V_t
        # Eq. (7) Yeomans et al. (1992)
        ρ_vec_t = r_a_t_b - r_t_t_t
        ρ_t = sqrt(ρ_vec_t[1]^2 + ρ_vec_t[2]^2 + ρ_vec_t[3]^2)
        # @show ρ_t
        # compute up-leg Shapiro delay
        # Sun barycentric position and velocity (in au, au/day) at transmit time (TDB)
        r_s_t_t = xvs(et_t_secs)[1:3]
        e_U_vec = r_t_t_t - r_s_t_t
        # @show e_U_vec
        e_U = sqrt(e_U_vec[1]^2 + e_U_vec[2]^2 + e_U_vec[3]^2) # heliocentric distance of Earth at t_t
        # @show e_U
        r_s_t_b = xvs(et_b_secs)[1:3] # barycentric position/velocity of Sun at bounce time
        p_U_vec = r_a_t_b - r_s_t_b
        p_U = sqrt(p_U_vec[1]^2 + p_U_vec[2]^2 + p_U_vec[3]^2) # heliocentric distance of asteroid at t_b
        # @show p_U
        q_U_vec = r_a_t_b - r_e_t_t
        q_U = ρ_t # signal path distance (up-leg)
        # @show q_U
        Δτ_rel_U = shapiro_delay(e_U, p_U, q_U) # seconds
        Δτ_tropo_U = tropo_delay(R_t, ρ_vec_t) # seconds
        # Δτ_corona_U = corona_delay(constant_term.(r_t_t_t), constant_term.(r_a_t_b), constant_term.(r_s_t_b), F_tx, station_code) # seconds
        Δτ_U = Δτ_rel_U # + Δτ_tropo_U #+ Δτ_corona_U # seconds
        # @show Δτ_U
        p_dot_12 = -dot(ρ_vec_t, v_t_t_t)/ρ_t
        Δt_1 = (τ_U - ρ_t/clightkms - Δτ_rel_U)/(1.0-p_dot_12/clightkms)
        # @show Δt_1
        τ_U = τ_U - Δt_1
        # transmit time, new estimate
        et_t_secs = et_b_secs - τ_U
    end

    # compute TDB-UTC at transmit time
    # corrections to TT-TDB from Moyer (2003) / Folkner et al. (2014) due to position of measurement station on Earth are of order 0.01μs
    # Δtt_tdb_station_t = - dot(v_e_t_t, r_t_t_t-r_e_t_t)/clightkms^2
    tdb_utc_t = tdb_utc(et_t_secs) # + Δtt_tdb_station_t
    # compute TDB-UTC at receive time
    # corrections to TT-TDB from Moyer (2003) / Folkner et al. (2014) due to position of measurement station on Earth  are of order 0.01μs
    # Δtt_tdb_station_r = - dot(v_e_t_r, r_r_t_r-r_e_t_r)/clightkms^2
    tdb_utc_r = tdb_utc(et_r_secs) # + Δtt_tdb_station_r
    # @show tdb_utc_t tdb_utc_r

    # # compute total time delay (UTC seconds); relativistic delay is already included in τ_D, τ_U
    # # Eq. (9) Yeomans et al. (1992)
    τ = (τ_D + τ_U) + (Δτ_tropo_D + Δτ_tropo_U) + (tdb_utc_t - tdb_utc_r) # seconds

    return 1e6τ # total signal delay (μs)
end

function delay_doppler(station_code::Int, t_r_utc::DateTime, F_tx::Real,
        niter::Int=10; eo::Bool=true, tc::Real=1.0, xve=earth_pv, xvs=sun_pv,
        xva=apophis_pv_197, autodiff::Bool=true, tord::Int=10)
    if autodiff
        τ = delay(station_code, t_r_utc, niter, eo=eo, xve=xve, xvs=xvs, xva=xva, tord=tord)
        return τ[0], -F_tx*τ[1]
    else
        τe = delay(station_code, t_r_utc,  tc/2, niter, eo=eo, xve=xve, xvs=xvs, xva=xva)
        τn = delay(station_code, t_r_utc,   0.0, niter, eo=eo, xve=xve, xvs=xvs, xva=xva)
        τs = delay(station_code, t_r_utc, -tc/2, niter, eo=eo, xve=xve, xvs=xvs, xva=xva)
        return τn, -F_tx*((τe-τs)/tc)
    end
end

function delay_doppler(astradarfile::String,
        niter::Int=10; eo::Bool=true, tc::Real=1.0, xve=earth_pv, xvs=sun_pv,
        xva=apophis_pv_197, autodiff::Bool=true, tord::Int=10)

    astradardata = process_radar_data_jpl(astradarfile)

    et1 = str2et(string(astradardata[1].utcepoch))
    a1_et1 = xva(et1)[1]
    S = typeof(a1_et1)

    vdelay = Array{S}(undef, length(astradardata))
    vdoppler = Array{S}(undef, length(astradardata))

    for i in eachindex(astradardata)
        vdelay[i], vdoppler[i] = delay_doppler(
            astradardata[i].rcvr,
            astradardata[i].utcepoch,
            astradardata[i].freq,
            niter,
            eo = eo,
            tc = tc,
            xve = xve,
            xvs = xvs,
            xva = xva,
            autodiff = autodiff,
            tord = tord
        )
    end

    delay_index = map(x->x.delay_units=="us", astradardata)
    doppler_index = map(x->x.doppler_units=="Hz", astradardata)

    radobs_t = DataFrame(
        (
            dt_utc_obs=utcepoch.(astradardata),
            τ_obs=delay.(astradardata),
            ν_obs=doppler.(astradardata),
            τ_comp=vdelay,
            ν_comp=vdoppler,
            σ_τ=delay_sigma.(astradardata),
            σ_ν=doppler_sigma.(astradardata),
            τ_units=delay_units.(astradardata),
            ν_units=doppler_units.(astradardata),
            freq=freq.(astradardata),
            rcvr=rcvr.(astradardata),
            xmit=xmit.(astradardata),
            bouncepoint=bouncepoint.(astradardata),
            delay_index=delay_index,
            doppler_index=doppler_index
        )
    )

    return radobs_t
end

# Compute round-trip time and Doppler shift radar astrometry for an asteroid at
# UTC instant `t_r_utc` from tracking station with code `station_code` from Earth,
# Sun and asteroid ephemerides. Dopplers are computed following Yeomans et al. (1992)
# station_code: observing station identifier (MPC nomenclature)
# et_r_secs: time of echo reception (TDB seconds since J2000.0 TDB)
# F_tx: transmitter frequency (MHz)
# niter: number of light-time solution iterations
# xve: Earth ephemeris wich takes et seconds since J2000 as input and returns Earth barycentric position in au and velocity in au/day
# xvs: Sun ephemeris wich takes et seconds since J2000 as input and returns Sun barycentric position in au and velocity in au/day
# xva: asteroid ephemeris wich takes et seconds since J2000 as input and returns asteroid barycentric position in au and velocity in au/day
function delay_doppler_yeomansetal92(station_code::Int, t_r_utc::DateTime,
        F_tx::Real, niter::Int=10; eo::Bool=true, xve=earth_pv, xvs=sun_pv,
        xva=apophis_pv_197)
    # Transform receiving time from UTC to TDB seconds since j2000
    et_r_secs = str2et(string(t_r_utc))
    # Compute geocentric position/velocity of receiving antenna in inertial frame (au, au/day)
    R_r, V_r = observer_position(station_code, et_r_secs, eo=eo)
    # Earth's barycentric position and velocity at receive time
    rv_e_t_r = xve(et_r_secs)
    r_e_t_r = rv_e_t_r[1:3]
    v_e_t_r = rv_e_t_r[4:6]
    # Receiver barycentric position and velocity at receive time
    r_r_t_r = r_e_t_r + R_r
    v_r_t_r = v_e_t_r + V_r
    # Asteroid barycentric position and velocity at receive time
    rv_a_t_r = xva(et_r_secs)
    r_a_t_r = rv_a_t_r[1:3]
    v_a_t_r = rv_a_t_r[4:6]
    # Sun barycentric position and velocity at receive time
    rv_s_t_r = xvs(et_r_secs)
    r_s_t_r = rv_s_t_r[1:3]
    v_s_t_r = rv_s_t_r[4:6]
    # down-leg iteration
    # τ_D first approximation: Eq. (1) Yeomans et al. (1992)
    ρ_vec_r = r_a_t_r - r_r_t_r
    ρ_r = sqrt(ρ_vec_r[1]^2 + ρ_vec_r[2]^2 + ρ_vec_r[3]^2)
    τ_D = ρ_r/clightkms # (seconds) -R_b/c, but delay is wrt asteroid Center (Brozovic et al., 2018)
    # bounce time, 1st estimate, Eq. (2) Yeomans et al. (1992)
    et_b_secs = et_r_secs - τ_D
    # asteroid barycentric position (in au) at bounce time (TDB)
    rv_a_t_b = xva(et_b_secs)
    r_a_t_b = rv_a_t_b[1:3]
    v_a_t_b = rv_a_t_b[4:6]

    Δτ_D = zero(τ_D)
    Δτ_rel_D = zero(τ_D)
    # Δτ_corona_D = zero(τ_D)
    Δτ_tropo_D = zero(τ_D)

    for i in 1:niter
        # Eq. (3) Yeomans et al. (1992)
        ρ_vec_r = r_a_t_b - r_r_t_r
        ρ_vec_dot_r = v_a_t_b - v_r_t_r
        # Eq. (4) Yeomans et al. (1992)
        ρ_r = sqrt(ρ_vec_r[1]^2 + ρ_vec_r[2]^2 + ρ_vec_r[3]^2)
        τ_D = ρ_r/clightkms # (seconds) -R_b/c (COM correction) + Δτ_D (relativistic, tropo, iono...)
        # bounce time, new estimate Eq. (2) Yeomans et al. (1992)
        et_b_secs = et_r_secs - τ_D
        # asteroid barycentric position (in au) at bounce time (TDB)
        rv_a_t_b = xva(et_b_secs)
        r_a_t_b = rv_a_t_b[1:3]
        v_a_t_b = rv_a_t_b[4:6]
        # compute down-leg Shapiro delay
        # NOTE: when using PPN, substitute 2 -> 1+γ in expressions for Shapiro delay, Δτ_rel_[D|U]
        e_D_vec  = r_e_t_r - r_s_t_r
        de_D_vec = v_e_t_r - v_s_t_r
        e_D = sqrt(e_D_vec[1]^2 + e_D_vec[2]^2 + e_D_vec[3]^2) # heliocentric distance of Earth at t_r
        rv_s_t_b = xvs(et_b_secs) # barycentric position and velocity of Sun at estimated bounce time
        r_s_t_b = rv_s_t_b[1:3]
        v_s_t_b = rv_s_t_b[4:6]
        p_D_vec  = constant_term.(r_a_t_b - r_s_t_b)
        dp_D_vec = constant_term.(v_a_t_b - v_s_t_b)
        p_D = sqrt(p_D_vec[1]^2 + p_D_vec[2]^2 + p_D_vec[3]^2) # heliocentric distance of asteroid at t_b
        q_D_vec  = constant_term.(r_a_t_b - r_e_t_r)
        dq_D_vec = constant_term.(v_a_t_b - v_e_t_r)
        q_D = sqrt(q_D_vec[1]^2 + q_D_vec[2]^2 + q_D_vec[3]^2) #signal path distance (down-leg)
        # Shapiro correction to time-delay
        Δτ_rel_D = shapiro_delay(e_D, p_D, q_D)
        # troposphere correction to time-delay
        Δτ_tropo_D = tropo_delay(R_r, ρ_vec_r) # seconds
        # Δτ_corona_D = corona_delay(constant_term.(r_a_t_b), r_r_t_r, r_s_t_r, F_tx, station_code) # seconds
        Δτ_D = Δτ_rel_D + Δτ_tropo_D #+ Δτ_corona_D # seconds
    end
    τ_D = τ_D + Δτ_D
    # get latest estimates of ρ_vec_r and ρ_r
    # Eq. (3) Yeomans et al. (1992)
    ρ_vec_r = r_a_t_b - r_r_t_r
    # Eq. (4) Yeomans et al. (1992)
    ρ_r = sqrt(ρ_vec_r[1]^2 + ρ_vec_r[2]^2 + ρ_vec_r[3]^2)

    # up-leg iteration
    # τ_U first estimation: Eq. (5) Yeomans et al. (1992)
    τ_U = τ_D
    # transmit time, 1st estimate Eq. (6) Yeomans et al. (1992)
    et_t_secs = et_r_secs - (τ_U+τ_D)
    # Geocentric position and velocity of transmitting antenna in inertial frame (au, au/day)
    R_t, V_t = observer_position(station_code, constant_term(et_t_secs), eo=eo)
    rv_e_t_t = xve(et_t_secs)
    r_e_t_t = rv_e_t_t[1:3]
    v_e_t_t = rv_e_t_t[4:6]
    # Transmitter barycentric position and velocity of at transmit time
    r_t_t_t = r_e_t_t + R_t
    v_t_t_t = v_e_t_t + V_t
    # Eq. (7) Yeomans et al. (1992)
    ρ_vec_t = r_a_t_b - r_t_t_t
    ρ_vec_dot_t = v_a_t_b - v_t_t_t
    ρ_t = sqrt(ρ_vec_t[1]^2 + ρ_vec_t[2]^2 + ρ_vec_t[3]^2)

    Δτ_U = zero(τ_U)
    Δτ_rel_U = zero(τ_U)
    # Δτ_corona_U = zero(τ_U)
    Δτ_tropo_U = zero(τ_U)

    for i in 1:niter
        # Eq. (8) Yeomans et al. (1992)
        τ_U = ρ_t/clightkms # (seconds) -R_b/c (COM correction) + Δτ_U (relativistic, tropo, iono...)
        # transmit time, new estimate
        et_t_secs = et_r_secs-(τ_U+τ_D)
        # Geocentric position and velocity of transmitting antenna in inertial frame (au, au/day)
        R_t, V_t = observer_position(station_code, constant_term(et_t_secs), eo=eo)
        # Earth's barycentric position and velocity at the transmit time
        rv_e_t_t = xve(et_t_secs)
        r_e_t_t = rv_e_t_t[1:3]
        v_e_t_t = rv_e_t_t[4:6]
        # Barycentric position and velocity of the transmitter at the transmit time
        r_t_t_t = r_e_t_t + R_t
        v_t_t_t = v_e_t_t + V_t
        # Eq. (7) Yeomans et al. (1992)
        ρ_vec_t = r_a_t_b - r_t_t_t
        ρ_vec_dot_t = v_a_t_b - v_t_t_t
        ρ_t = sqrt(ρ_vec_t[1]^2 + ρ_vec_t[2]^2 + ρ_vec_t[3]^2)
        # compute up-leg Shapiro delay
        # Sun barycentric position and velocity (in au, au/day) at transmit time (TDB)
        rv_s_t_t = xvs( et_t_secs )
        r_s_t_t = rv_s_t_t[1:3]
        v_s_t_t = rv_s_t_t[4:6]
        e_U_vec = constant_term.(r_e_t_t - r_s_t_t)
        de_U_vec = constant_term.(v_e_t_t - v_s_t_t)
        e_U = sqrt(e_U_vec[1]^2 + e_U_vec[2]^2 + e_U_vec[3]^2) # heliocentric distance of Earth at t_t
        rv_s_t_b = xvs(et_b_secs) # barycentric position/velocity of Sun at bounce time
        r_s_t_b = rv_s_t_b[1:3]
        v_s_t_b = rv_s_t_b[4:6]
        p_U_vec = constant_term.(r_a_t_b - r_s_t_b)
        dp_U_vec = constant_term.(v_a_t_b - v_s_t_b)
        p_U = sqrt(p_U_vec[1]^2 + p_U_vec[2]^2 + p_U_vec[3]^2) # heliocentric distance of asteroid at t_b
        q_U_vec = constant_term.(r_a_t_b - r_e_t_t)
        dq_U_vec = constant_term.(v_a_t_b - v_e_t_t)
        q_U = sqrt(q_U_vec[1]^2 + q_U_vec[2]^2 + q_U_vec[3]^2) # signal path (up-leg)
        Δτ_rel_U = shapiro_delay(e_U, p_U, q_U) # seconds
        Δτ_tropo_U = tropo_delay(R_t, ρ_vec_t) # seconds
        # Δτ_corona_U = corona_delay(constant_term.(r_t_t_t), constant_term.(r_a_t_b), constant_term.(r_s_t_b), F_tx, station_code) # seconds
        Δτ_U = Δτ_rel_U + Δτ_tropo_U #+ Δτ_corona_U # seconds
    end
    # Sun barycentric position (in au) at transmit time (TDB)
    r_s_t_t = xvs(et_t_secs)[1:3]
    τ_U = τ_U + Δτ_U

    # compute TDB-UTC at transmit time
    # corrections to TT-TDB from Moyer (2003) / Folkner et al. (2014) due to position of measurement station on Earth are of order 0.01μs
    # Δtt_tdb_station_t = - dot(v_e_t_t/daysec, r_t_t_t-r_e_t_t)/c_au_per_sec^2
    tdb_utc_t = tdb_utc(constant_term(et_t_secs)) #+ Δtt_tdb_station_t
    # compute TDB-UTC at receive time
    # corrections to TT-TDB from Moyer (2003) / Folkner et al. (2014) due to position of measurement station on Earth  are of order 0.01μs
    # Δtt_tdb_station_r = dot(v_e_t_r/daysec, r_r_t_r-r_e_t_r)/c_au_per_sec^2
    tdb_utc_r = tdb_utc(et_r_secs) #+ Δtt_tdb_station_r

    # compute total time delay (UTC seconds)
    # Eq. (9) Yeomans et al. (1992)
    τ = (τ_U + τ_D) + (tdb_utc_t - tdb_utc_r)

    # compute Doppler shift ν
    # Eq. (10) Yeomans et al. (1992)
    ρ_vec_dot_t = v_a_t_b - v_t_t_t
    ρ_vec_dot_r = v_a_t_b - v_r_t_r
    # Eq. (11) Yeomans et al. (1992)
    ρ_dot_t = dot(ρ_vec_t, ρ_vec_dot_t)/ρ_t
    ρ_dot_r = dot(ρ_vec_r, ρ_vec_dot_r)/ρ_r
    # Eq. (12) Yeomans et al. (1992)
    doppler_c = (ρ_dot_t+ρ_dot_r)/clightkms
    p_t = dot(ρ_vec_t, v_t_t_t)/ρ_t
    p_r = dot(ρ_vec_r, v_r_t_r)/ρ_r
    doppler_c2_t1 = ρ_dot_t*p_t - ρ_dot_r*p_r - ρ_dot_t*ρ_dot_r # order c^-2, 1st term
    r_ts_vec = r_t_t_t - r_s_t_t
    r_ts = sqrt(r_ts_vec[1]^2+r_ts_vec[2]^2+r_ts_vec[3]^2)
    r_rs_vec = r_r_t_r - r_s_t_r
    r_rs = sqrt(r_rs_vec[1]^2+r_rs_vec[2]^2+r_rs_vec[3]^2)
    doppler_c2_t2 = (PlanetaryEphemeris.μ[su]*((au^3)/(daysec^2)))*( (1/r_ts) - (1/r_rs) ) # order c^-2, 2nd term
    doppler_c2_t3 = (  dot(v_t_t_t, v_t_t_t) - dot(v_r_t_r, v_r_t_r)  )/2  # order c^-2, 3rd term
    doppler_c2 = (doppler_c2_t1 + doppler_c2_t2 + doppler_c2_t3)/(clightkms^2)
    ν = -F_tx*(doppler_c) # + doppler_c2)

    return 1e6τ, 1e6ν # total signal delay (μs) and Doppler shift (Hz)
end


text:
immutable BicubicSpline <: AbstractInterpolation
    BicubicSpline() = BicubicSpline
end

immutable BicubicSplineInterpolation{T}
    x::Vector{T}
    y::Vector{T}
    z::Array{T,2}
    A::Array{T,4}
end

function default_settings{D,T}(::Type{BicubicSpline},axes::NTuple{D,T})
    nx = length(axes[1])
    ny = length(axes[2])
    return Irregular{2,Val{(nx,ny)}}, Nearest
end

function deriv(x,y,i)
    nx = length(x)
    if i == 1
        dy = y[i+1] - y[i]
        return dy
    end
    if i == nx
        dy = y[i] - y[i-1]
        return dy
    end
    dy1 = y[i+1] - y[i]
    dy2 = y[i] - y[i-1]
    return 0.5*(dy1 + dy2)
end

function BicubicSpline(x::AbstractVector,y::AbstractVector, z::AbstractArray)
    issorted(x) || throw(ArgumentError("x points must be in ascending order"))
    issorted(y) || throw(ArgumentError("y points must be in ascending order"))
    x = convert(Array{Float64,1},x)
    y = convert(Array{Float64,1},y)
    z = convert(Array{Float64,2},z)
    nx = length(x)
    ny = length(y)
    nx == size(z,1) || throw(ArgumentError("length(x) != size(z,1)"))
    ny == size(z,2) || throw(ArgumentError("length(y) != size(z,2)"))

    px = similar(z)
    py = similar(z)
    pxy = similar(z)

    for i=1:nx
        for j = 1:ny
            px[i,j] = deriv(x,z[:,j],i)
            py[i,j] = deriv(y,z[i,:],j)
        end
        for j=1:ny
            pxy[i,j] = deriv(y,px[i,:],j)
        end
    end
    C = [1 0 0 0; 0 0 1 0; -3 3 -2 -1; 2 -2 1 1]
    A = zeros(4,4,nx,ny)
    for i=1:nx-1
        for j=1:ny-1
            F = [z[i,j]    z[i,j+1]    py[i,j]    py[i,j+1];
                 z[i+1,j]  z[i+1,j+1]  py[i+1,j]  py[i+1,j+1];
                 px[i,j]   px[i,j+1]   pxy[i,j]   pxy[i,j+1];
                 px[i+1,j] px[i+1,j+1] pxy[i+1,j] pxy[i+1,j+1]]
            A[:,:,i,j] = C*(F*C')
        end
    end

    x,y,z, A = promote(x,y,z,A)
    return BicubicSplineInterpolation{eltype(x)}(x,y,z,A)
end

@inline function preprocess(::Type{BicubicSpline}, axes, f)
    return BicubicSpline(axes...,f)
end

function apply_bc(::Type{BicubicSpline},BC::Type{Nearest}, I::Symbol, x::Symbol)
    xbc = apply_nearest(I, 1, x, 1)
    ybc = apply_nearest(I, 2, x, 2)
    quote
        $(xbc)
        $(ybc)
    end
end

function apply_bc(::Type{BicubicSpline},BC::Type{Reflect}, I::Symbol, x::Symbol)
    xbc = apply_reflect(I,1,x,1)
    ybc = apply_reflect(I,2,x,2)
    quote
        $(xbc)
        $(ybc)
    end
end

function apply_bc{V}(::Type{BicubicSpline},BC::Type{Value{V}}, I::Symbol, x::Symbol)
    xbc = apply_value(BC,I,1,x,1)
    ybc = apply_value(BC,I,2,x,2)
    quote
        $(xbc)
        $(ybc)
    end
end

function apply_bc(::Type{BicubicSpline},BC::Type{None}, I::Symbol, x::Symbol)
    apply_none()
end

function apply_bc(::Type{BicubicSpline},BC::Type{Error}, I::Symbol, x::Symbol)
    xbc = apply_error(I,1,x,1)
    ybc = apply_error(I,2,x,2)
    quote
        $(xbc)
        $(ybc)
    end
end

@inline function get_indices{S}(::Type{BicubicSpline}, G::Type{Uniform{2,S}}, I, x)
    dx = I.x[2]-I.x[1]
    ix = clamp(Int(floor((x[1] - I.x[1])/dx)) + 1, 1, S.parameters[1][1]-1)
    dy = I.y[2]-I.y[1]
    iy = clamp(Int(floor((x[2] - I.y[1])/dy)) + 1, 1, S.parameters[1][2]-1)
    return (ix, iy)
end

@inline function get_indices{S}(::Type{BicubicSpline}, G::Type{Irregular{2,S}}, I, x)
    xr = searchsortedlast(I.x,x[1],Base.Order.Forward)
    ix = clamp(xr, 1, S.parameters[1][1]-1)
    yr = searchsortedlast(I.y,x[2],Base.Order.Forward)
    iy = clamp(yr, 1, S.parameters[1][2]-1)
    return (ix, iy)
end

@inline function get_value(::Type{BicubicSpline}, I, x, inds)
    ix, iy = inds
    A = I.A[:,:,ix,iy]
    dx = I.x[ix+1] - I.x[ix]
    xx = (x[1] - I.x[ix])/dx
    dy = I.y[iy+1] - I.y[iy]
    yy = (x[2] - I.y[iy])/dy
    v = (1.0) *(A[1,1] + A[1,2]*yy + A[1,3]*(yy^2) + A[1,4]*(yy^3)) +
        (xx)  *(A[2,1] + A[2,2]*yy + A[2,3]*(yy^2) + A[2,4]*(yy^3)) +
        (xx^2)*(A[3,1] + A[3,2]*yy + A[3,3]*(yy^2) + A[3,4]*(yy^3)) +
        (xx^3)*(A[4,1] + A[4,2]*yy + A[4,3]*(yy^2) + A[4,4]*(yy^3))
    return v
end


text:
"""
    struct DemoPage <: Any
    DemoPage(root::String)

Constructs a demo page object.

# Fields

Besides the root path to the demo page folder `root`, this struct has some other fields:

* `title`: page title
* `template`: template content of the demo page.
* `sections`: demo sections found in `root`

# Configuration

You can manage an extra `config.json` file to customize rendering of a demo page.
Supported items are:

* `order`: specify the sections order. By default, it's case-insensitive alphabetic order.
* `template`: path to template filename. The content of the template file should has one and only one `{{{democards}}}`.
* `title`: specify the title of this demo page. By default, it's the folder name of `root`. Will be override by `template`.

The following is an example of `config.json`:

```json
{
    "template": "template.md",
    "order": [
        "basic",
        "advanced"
    ]
}
```

# Examples

The following is the simplest folder structure of a `DemoPage`:

```text
demos
└── basic
    ├── demo_1.md
    ├── demo_2.md
    ├── demo_3.md
    ├── demo_4.md
    ├── demo_5.md
    ├── demo_6.md
    ├── demo_7.md
    └── demo_8.md
```

!!! note

    A `DemoPage` doesn't manage demo files directly, so here you'll need
    a DemoSection `basic` to manage them.

The following is a typical folder structure of a `DemoPage`:

```text
demos
├── advanced
│   ├── advanced_demo_1.md
│   └── advanced_demo_2.md
├── basic
│   ├── part1
│   │   ├── basic_demo_1_1.md
│   │   └── basic_demo_1_2.md
│   └── part2
│       ├── config.json
│       ├── basic_demo_2_1.md
│       └── basic_demo_2_2.md
├── config.json
└── template.md
```

!!! warning

    A section should only hold either subsections or demo files. A folder that has both subfolders and demo files (e.g., `*.md`) is invalid.

See also: [`MarkdownDemoCard`](@ref DemoCards.MarkdownDemoCard), [`DemoSection`](@ref DemoCards.DemoSection)
"""
struct DemoPage
    root::String
    sections::Vector{DemoSection}
    template::String
    title::String
end

basename(page::DemoPage) = basename(page.root)

function DemoPage(root::String)::DemoPage
    root = replace(root, r"[/\\]" => Base.Filesystem.path_separator) # windows compatibility
    isdir(root) || throw(ArgumentError("page root should be a valid dir, instead it's $(root)"))
    root = rstrip(root, '/') # otherwise basename(root) will returns `""`

    path = joinpath.(root, filter(x->!startswith(x, "."), readdir(root))) # filter out hidden files
    section_paths = filter(x->isdir(x)&&!(basename(x) in ignored_dirnames),
                           path)
    sections = map(DemoSection, section_paths)

    # first consturct an incomplete page
    # then load the config and reconstruct a new one
    page = DemoPage(root, sections, "", "")

    section_paths = joinpath.(root, load_config(page, "order"))
    ordered_sections = map(DemoSection, section_paths) # TODO: technically, we don't need to regenerate sections here

    title = load_config(page, "title")
    page = DemoPage(root, ordered_sections, "", title)

    # default template requires a title
    template = load_config(page, "template")
    DemoPage(root, ordered_sections, template, title)
end

function load_config(page::DemoPage, key)
    path = joinpath(page.root, config_filename)
    config = isfile(path) ? JSON.parsefile(path) : Dict()

    if key == "order"
        haskey(config, key) || return get_default_order(page)

        order = config[key]
        validate_order(order, page)
        return order
    elseif key == "title"
        return load_title(page, config)
    elseif key == "template"
        haskey(config, key) || return get_default_template(page)

        template_file = replace(config[key],
                                r"[/\\]" => Base.Filesystem.path_separator) # windows compatibility
        template_path = joinpath(page.root, template_file)

        check_ext(template_path, :markdown)
        if 1 != sum(occursin.("{{{democards}}}", readlines(template_path)))
            throw(ArgumentError("invalid template file $(template_path): it should has one and only one {{{democards}}}"))
        end

        return read(template_path, String)
    else
        throw(ArgumentError("Unrecognized key $(key) for DemoPage"))
    end
end

"""
    load_title(page::DemoPage, config)

load title from config using the following priority rules:

1. parse the title out from the template file
2. use config["title"]
3. fallback: use the folder name
"""
function load_title(page::DemoPage, config)
    key = "title"
    title = load_template_config(page, key)

    if haskey(config, key)
        if isnothing(title)
            return config[key]
        else
            path = joinpath(page.root, config_filename)
            @warn("config item $(key) in $(path) is suppressed by \"template\"")
            return title
        end
    else
        if isnothing(title)
            return basename(page)
        else
            return title
        end
    end
end

get_default_order(page::DemoPage) =
    sort(basename.(page.sections), by = x->lowercase(x))

function get_default_template(page::DemoPage)
    header = "# $(page.title)\n\n"
    # TODO: by doing this we loss the control on section-level template
    content = "{{{democards}}}" # render by Mustache
    footer = ""
    return header * content * footer
end

"""
    load_template_config(page::DemoPage, key)

parse out and return value of `key` from the template file of `page`.

If `page` doesn't have a template file, it returns `nothing`.
"""
function load_template_config(page::DemoPage, key)
    path = joinpath(page.root, config_filename)
    config = isfile(path) ? JSON.parsefile(path) : Dict()
    template_path = joinpath(page.root, get(config, "template", ""))
    isfile(template_path) || return nothing

    config = parse_markdown(template_path)
    get(config, key, nothing)
end


text:
function fpt(mu::Float64, bound::Float64, Δt::Float64, tmax::Float64)
  # mu::Float64=1.0
  # bound::Float64=1.0
  # dt::Float64=0.01
  # tmax::Float64=4

  nsteps = Int(tmax/Δt)
  g1 = Array{Float64}(undef, nsteps)
  g2 = Array{Float64}(undef, nsteps)

  err = ccall((:fpt, "./ddm_fpt_lib"),
                Cint, (Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}),
                mu, bound, Δt, tmax, g1, g2)

  return g1, g2
end

function rand_asym(mu::Float64, bound_lo::Float64, bound_hi::Float64, Δt_rand::Float64, n::Int, seed::Int)
  # mu::Float64=1.0
  # bound_lo::Float64=-1.5
  # bound_hi::Float64=1.0
  # dt_rand::Float64=0.001
  # n::Int32=10000
  # seed::Int32=0

  t = Array{Float64}(undef,n)
  bound_cond = Array{Int}(undef,n)

  err=ccall((:rand_asym,"./ddm_fpt_lib"),Cint,(Cdouble,Cdouble,Cdouble,Cdouble,Cint,Cint,Ptr{Cdouble},Ptr{Cint}),
            mu,bound_lo,bound_hi, Δt_rand, n, seed, t, bound_cond)

  return t, bound_cond
end

# mu::Float64=1.0
# bound::Float64=1.0
# dt::Float64=0.01
# tmax::Float64=4

# bound_lo::Float64=-1.5
# bound_hi::Float64=1.0
# Δt_rand::Float64=0.001
# n::Int=10000
# seed::Int=0

g1, g2 = fpt(1.0, 1.0, 0.01, 10.0)
rt, choice = rand_asym(1.0, 1.5, 1.0, 0.01, 1, 123)

max_stars_repo_name: awellis/dm

max_issues_repo_name: awellis/dm

max_forks_repo_name: awellis/dm


text:
using HDFS

host = isempty(ARGS) ? "localhost" : ARGS[1]
port = isempty(ARGS) ? 9000 : int(ARGS[2])

# connect
fs = hdfs_connect(host, port)
println("connected to hdfs:://$(host):$(port)")

# print capacity, used
capacity = hdfs_get_capacity(fs)
used = hdfs_get_used(fs)
println("$(used) used from total $(capacity) bytes")
@assert 0 < used < capacity

blk_sz = hdfs_get_default_block_size(fs)
println("default block sz: $(blk_sz)")
@assert blk_sz > 0

# get current working directory
pwd = hdfs_pwd(fs)
println("current directory: $(pwd)")
@assert length(pwd) > 0

# make directory
const TEST_DIR = "test_dir"
const TEST_FILE = "test_file"
hdfs_exists(fs, TEST_DIR) && hdfs_delete(fs, TEST_DIR)
@assert 0 == hdfs_mkdir(fs, TEST_DIR)
println("created new directory $(TEST_DIR)")
@assert true == hdfs_is_directory(fs, TEST_DIR)
@assert 0 == hdfs_cd(fs, TEST_DIR)
pwd1 = hdfs_pwd(fs)
println("current directory: $(pwd1)")
@assert pwd1 == string(pwd, "/", TEST_DIR)

@assert false == hdfs_exists(fs, TEST_FILE)                         # try exists on a file
fi = hdfs_open(fs, TEST_FILE, "w")                                  # create the file
println("created new file $(TEST_FILE)")
@assert 600 == hdfs_write(fi, repeat("hello ", 100))                # write to the file
@assert 0 == hdfs_flush(fi)                                         # flush writes
@assert 0 == hdfs_close(fi)                                         # close the file
println("wrote 600 bytes to file")
@assert false == hdfs_is_directory(fs, TEST_FILE)

@assert true == hdfs_exists(fs, TEST_FILE)                          # try exists on a file
fi = hdfs_open(fs, TEST_FILE, "r")                                  # read from the file
avlb = hdfs_available(fi)
println("bytes available: $(avlb)")
finfo = hdfs_get_path_info(fs, TEST_FILE)
println("file size bytes: $(finfo.size)")
@assert 600 == finfo.size
@assert 0 <= avlb <= finfo.size

println("reading...")
(buff, bytes) = hdfs_read(fi, 5)
readstr = bytestring(buff[1:bytes])
@assert "hello" == readstr

@assert 5 == hdfs_tell(fi)
@assert 0 == hdfs_seek(fi, 2)

(buff, bytes) = hdfs_read(fi, 5)
readstr = bytestring(buff[1:bytes])
@assert "llo h" == readstr

# test read at end of file
@assert 0 == hdfs_seek(fi, finfo.size-10)
(buff, bytes) = hdfs_read(fi, 50)
@assert 50 == length(buff)
@assert 10 == bytes

println("closing...")
@assert 0 == hdfs_close(fi)

blocks = hdfs_blocks(fs, TEST_FILE, 0, finfo.size)
@assert length(blocks) == 1

const NEW_TEST_FILE = string(TEST_FILE, ".new")
const NEW_TEST_FILE_2 = string(TEST_FILE, ".new2")
@assert 0 == hdfs_copy(fs, TEST_FILE, fs, NEW_TEST_FILE)
@assert true == hdfs_exists(fs, NEW_TEST_FILE)
@assert 0 == hdfs_move(fs, NEW_TEST_FILE, fs, NEW_TEST_FILE_2)
@assert true == hdfs_exists(fs, NEW_TEST_FILE_2)
@assert false == hdfs_exists(fs, NEW_TEST_FILE)

println("listing the directory...")
entries = hdfs_list_directory(fs, ".")
@assert 2 == length(entries)
for finfo in entries
    println("\t$(finfo.name)")
end

println("cleaning up...")
@assert 0 == hdfs_cd(fs, "..")
@assert 0 == hdfs_delete(fs, TEST_DIR)

println("passed")


text:
struct PSHestonAssetPricingModelParameters <: PSAbstractAssetReturnModel

    # data -
    μ::Float64
    𝝹::Float64
    ϴ::Float64
    σ::Float64
    𝜌::Float64

    function PSHestonAssetPricingModelParameters(μ::Float64, 𝝹::Float64, ϴ::Float64, σ::Float64, 𝜌::Float64)
        this = new(μ, 𝝹, ϴ, σ, 𝜌)
    end
end

struct PSGeometricBrownianMotionModelParameters <: PSAbstractAssetReturnModel

    # model parameters -
    μ::Float64
    σ::Float64

    function PSGeometricBrownianMotionModelParameters(μ=0.0, σ=0.0)
        this = new(μ, σ)
    end
end

struct PSSingleIndexModelParameters <: PSAbstractAssetReturnModel

    # model parameters -
    ⍺::Float64
    β::Float64
    riskFreeRate::Float64
    𝝐::ContinuousUnivariateDistribution

    function PSSingleIndexModelParameters(⍺, β, riskFreeRate, 𝝐)
        this = new(⍺, β, riskFreeRate, 𝝐)
    end
end

struct PSRandomWalkModelParameters <: PSAbstractAssetReturnModel

    # model parameters -
    𝝙::UnivariateKDE

    function PSRandomWalkModelParameters(𝝙::UnivariateKDE)
        this = new(𝝙)
    end
end


text:
mutable struct Vertical <: Geometry
    label::String
    shape::Vector{Box}
    nsample::Int
end

function Vertical( ; label, shape, nsample)
    return Vertical(label, shape, nsample)
end


text:
immutable HomotheticHyperbolic{T} <: Model
    a::T
    α::T
end

function Base.show(io::IO, ::Type{HomotheticHyperbolic})
    return print(io, "Homothetic Hyperbolic")
end

function Base.show(io::IO, m::HomotheticHyperbolic)
    @printf(
        io,
        "Homothetic Hyperbolic (a = %.3f, α = %.3f)",
        m.a,
        m.α,
    )
end

function initial_parameters(::Type{HomotheticHyperbolic})
    return Float64[log(1.0), log(1.0)]
end

canonical_parameters(m::HomotheticHyperbolic) = (m.a, m.α)

function HomotheticHyperbolic(parameters::Vector)
    a = exp(parameters[1])
    α = exp(parameters[2])
    return HomotheticHyperbolic(a, α)
end

@inline function predict(
    m::HomotheticHyperbolic,
    x2::Real,
    t2::Real,
    x1::Real,
    t1::Real,
)
    a, α = m.a, m.α
    u1 = log(x1 * (1 / (1 + α * t1)))
    u2 = log(x2 * (1 / (1 + α * t2)))
    z = a * (u2 - u1)
    p = invlogit(z)
    return p
end

@inline function gradient_component(
    m::HomotheticHyperbolic,
    x2::Real,
    t2::Real,
    x1::Real,
    t1::Real,
)
    a, α = m.a, m.α
    d_t2 = 1 / (1 + α * t2)
    d_t1 = 1 / (1 + α * t1)
    gr = (
        (log(x2 * d_t2) - log(x1 * d_t1)) * a,
        a * (t1 * d_t1 - t2 * d_t2) * α,
    )
    return gr
end


text:
using GViz

gp = GenomicPos(;height=1cm, margin_top = 8mm)
dt = DataTrack(x = [0,1,2,3,4,5],
               y = [7, 10, 20, 30, 40, 42],
               fill = ["red", "green", "blue","red", "green", "blue"],
               size = [3,1,2,3,4,5], alpha = 0.3,
               layers = [Dict(:geom => :line, :size => 1, :color => "black", :alpha => 1.0),
                         Dict(:geom => :point)])

#doc = RefFrame(;width=20cm, height = 3cm, defaultparams...)
#p0 = RefFrame(transform = Translation(0cm, 0cm), children=[gp])
#p1 = RefFrame(transform=Translation(0cm, 2cm), children=[dt], title="Track1")
#p2 = RefFrame(transform=Translation(0cm, 5.5cm), children=[dt], title="Track2")
#
#push!(doc.children, p0)
#push!(doc.children, p1)
#push!(doc.children, p2)
#
#pdf(doc, "/tmp/dt.pdf", 0, 5; absolute=true)

plotTracks([gp, dt, dt, gp], "/tmp/dt.pdf", -0.1, 5.2; defaultparams..., ymin=0, ymax=50)

dt = DataTrack()

p5 = RefFrame(children=[dt], title = "Track 5",
              shape="diamond", size=1, symbol_line_width = 0,
              fill = "red",
              x = range(0; stop = 5, length = 500), y = rand(500))

plotTracks([gp, dt, p5], "/tmp/dt.pdf", 0, 5; track_spacing = 2mm,defaultparams..., width=20cm, title="default", x = [0,1,2,3,4,5],y = [7, 10, 20, 30, 40, 42])

# chr18:80,139,089-80,265,529
plotTracks([gp], "/tmp/dt.pdf", 80139089, 80265529; defaultparams...)

# ---- GenomicFeatureTrack

chrom_ = "chr20"
bpStart_ = 31400000
bpEnd_ = 31550229

gp = GenomicPos(;height=1cm, margin_top = 8mm)

gfs = genes(GTFFile(tabix(GViz.pkgpath("data", "gencode.v35.annotation.gtf.bgz"),
                          chrom_, bpStart_, bpEnd_)))
gt = GenomicFeatureTrack(;features = gfs, title = "Genes", feature_height=5mm, feature_fill = "darkgreen",
                         feature_label_spacing=1mm, feature_label_color = "black")

gt = GeneTrack(chrom_, bpStart_, bpEnd_;)
               title = "Genes",
               feature_height=5mm,
               feature_fill = "darkgreen",
               feature_label_spacing=1mm,
               feature_label_color = "black")

dt = DataTrack(x = range(bpStart_, bpEnd_, length=1000), y = rand(1000),
               ymin = -0.05, ymax = 1.05, ylab = "Intensity", title = "Random data", height=5cm )

plotTracks([gp, dt, gt, gp], "/tmp/dt.pdf", bpStart_, bpEnd_; defaultparams..., width=20cm)

# ---- gorfile

using GorJulia

meth = GorFile(GViz.pkgpath("data", "methylationFraction.gor.tsv");
               first = (chrom_, bpStart_), last = (chrom_, bpEnd_)) |> Tables.columns

mt = DataTrack(x = meth[:Pos], y = meth[:methFrac], ymin=0, ymax = 1, title="Methylation",
               margin_top=0mm, color="red", fill="red", size=1)
plotTracks([gp, mt, gt], "/tmp/dt.pdf", bpStart_, bpEnd_; defaultparams..., width=20cm)

# ---- multiple layer dicts in DataTrack

hmeth = GorFile(GViz.pkgpath("data", "PofO_ASM.gor");
                first = (chrom_, bpStart_), last = (chrom_, bpEnd_)) |> Tables.columns

hmt = DataTrack(x = hmeth[:Pos],
                layers = [Dict(:y => hmeth[:mfrac], :fill => "blue"),
                          Dict(:y => hmeth[:pfrac], :fill => "red")],
                ymin=0, ymax = 1, title="ASM",
                size=1, margin_top=0mm, color="white")

zt = DataTrack(x = hmeth[:Pos], y = -log10.(hmeth[:pval]), ymin=0, ymax = 20, title="ASM -log10(p)",
               size=1, margin_top=0mm, geom = :col, color = "darkblue")

plotTracks([gp, zt,mt, hmt, gt], "/tmp/dt.pdf", bpStart_, bpEnd_; defaultparams..., width=25cm)


text:
module GridapP4estTests

using GridapP4est
using MPI
using ArgParse
using Test

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--image-file", "-i"
        help = "Path to the image file that one can use in order to accelerate MPI tests"
        arg_type = String
        default="GridapDistributed.so"
    end
    return parse_args(s)
end

parsed_args = parse_commandline()
image_file_path=parsed_args["image-file"]
image_file_exists=isfile(image_file_path)

nprocs_str = get(ENV, "JULIA_GRIDAP_P4EST_TEST_NPROCS","")
nprocs = nprocs_str == "" ? clamp(Sys.CPU_THREADS, 2, 4) : parse(Int, nprocs_str)
#mpiexec_args = Base.shell_split("--allow-run-as-root --tag-output") #Base.shell_split(get(ENV, "JULIA_MPIEXEC_TEST_ARGS", ""))
testdir = @__DIR__
istest(f) = endswith(f, ".jl") && !(f=="runtests.jl")
testfiles = sort(filter(istest, readdir(testdir)))
@time @testset "$f" for f in testfiles
  MPI.mpiexec() do cmd
     if f in ["UniformlyRefinedForestOfOctreesDiscreteModelsTests.jl"]
       np = 4
       extra_args = "-s 2 2 -r 2"
     else
       np = nprocs
       extra_args = ""
     end
     if ! image_file_exists
       cmd = `$cmd -n $(np) --allow-run-as-root --oversubscribe $(Base.julia_cmd()) --project=. $(joinpath(testdir, f)) $(split(extra_args))`
     else
       cmd = `$cmd -n $(np) --allow-run-as-root --oversubscribe $(Base.julia_cmd()) -J$(image_file_path) --project=. $(joinpath(testdir, f)) $(split(extra_args))`
     end
     @show cmd
     run(cmd)
     @test true
  end
end

end # module


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

include("abstract_variate_transform.jl")
include("distribution_transform.jl")


text:
export Evaluation,ivp,bvp

## Evaluation constructors

abstract Evaluation{T}<:Operator{T}

@functional Evaluation

# M = Bool if endpoint
immutable ConcreteEvaluation{S,M,OT,T} <: Evaluation{T}
    space::S
    x::M
    order::OT
end
Evaluation{T}(::Type{T},sp::UnivariateSpace,x::Bool,order::Integer) =
    ConcreteEvaluation{typeof(sp),typeof(x),typeof(order),T}(sp,x,order)
function Evaluation{T}(::Type{T},sp::UnivariateSpace,x::Number,order::Integer)
    d=domain(sp)
    if isa(d,IntervalDomain) && isapprox(first(d),x)
        Evaluation(T,sp,false,order)
    elseif isa(d,IntervalDomain) && isapprox(last(d),x)
        Evaluation(T,sp,true,order)
    else
        ConcreteEvaluation{typeof(sp),typeof(x),typeof(order),T}(sp,x,order)
    end
end

Evaluation(sp::UnsetSpace,x::Bool,k::Integer) =
    ConcreteEvaluation{UnsetSpace,Bool,typeof(k),UnsetNumber}(sp,x,k)
Evaluation(sp::Space{ComplexBasis},x,order::Integer) =
    Evaluation(Complex{real(eltype(domain(sp)))},sp,x,order)
Evaluation(sp::Space,x,order::Integer) = Evaluation(eltype(domain(sp)),sp,x,order)

#Evaluation(sp::UnsetSpace,x::Bool)=Evaluation(sp,x,0)
Evaluation(d::Space,x::Union{Number,Bool}) = Evaluation(d,x,0)

Evaluation(d::Domain,x::Union{Number,Bool},n...) = Evaluation(Space(d),x,n...)
Evaluation(x::Union{Number,Bool}) = Evaluation(UnsetSpace(),x,0)
Evaluation(x::Union{Number,Bool},k::Integer) = Evaluation(UnsetSpace(),x,k)
Evaluation{T<:Number}(d::Vector{T},x::Union{Number,Bool},o::Integer) = Evaluation(Interval(d),x,o)

rangespace{S<:AmbiguousSpace}(E::ConcreteEvaluation{S,Bool}) = ConstantSpace()
rangespace{S}(E::ConcreteEvaluation{S,Bool}) = ConstantSpace(Point(E.x?last(domain(E)):first(domain(E))))
rangespace(E::ConcreteEvaluation) = ConstantSpace(Point(E.x))

function Base.convert{T}(::Type{Operator{T}},E::ConcreteEvaluation)
    if T == eltype(E)
        E
    else
        ConcreteEvaluation{typeof(E.space),typeof(E.x),typeof(E.order),T}(E.space,E.x,E.order)
    end
end

## default getindex
getindex(D::ConcreteEvaluation,k::Integer) =
    differentiate(Fun([zeros(eltype(D),k-1);one(eltype(D))],D.space),D.order)(D.x)

function getindex{S}(D::ConcreteEvaluation{S,Bool},k::Integer)
    T=eltype(D)
    if !D.x
        first(differentiate(Fun([zeros(T,k-1);one(T)],D.space),D.order))
    else
        last(differentiate(Fun([zeros(T,k-1);one(T)],D.space),D.order))
    end
end

## EvaluationWrapper

immutable EvaluationWrapper{S<:Space,M,FS<:Operator,OT,T<:Number} <: Evaluation{T}
    space::S
    x::M
    order::OT
    functional::FS
end

EvaluationWrapper(sp::Space,x,order,func::Operator) =
    EvaluationWrapper{typeof(sp),typeof(x),typeof(func),typeof(order),eltype(sp)}(sp,x,order,func)
getindex(E::EvaluationWrapper,k) = E.functional[k]

domainspace(E::Evaluation) = E.space
domain(E::Evaluation) = domain(E.space)
promotedomainspace{T}(E::Evaluation{T},sp::Space) =
    Evaluation(promote_type(T,eltype(sp)),sp,E.x,E.order)
Base.stride(E::EvaluationWrapper)=stride(E.functional)

## Convenience routines

evaluate(d::Domain,x) = Evaluation(d,x)
ldiffbc(d,k) = Evaluation(d,false,k)
rdiffbc(d,k) = Evaluation(d,true,k)
diffbcs(d,k) = [ldiffbc(d,k),rdiffbc(d,k)]

ldirichlet(d) = ldiffbc(d,0)
rdirichlet(d) = rdiffbc(d,0)
lneumann(d) = ldiffbc(d,1)
rneumann(d) = rdiffbc(d,1)

dirichlet(d) = diffbcs(d,0)
neumann(d) = diffbcs(d,1)

ivp(d,k) = Operator{eltype(d)}[ldiffbc(d,i) for i=0:k-1]
bvp(d,k) = vcat(Operator{eltype(d)}[ldiffbc(d,i) for i=0:div(k,2)-1],
                Operator{eltype(d)}[rdiffbc(d,i) for i=0:div(k,2)-1])

periodic(d,k) = Operator{eltype(d)}[Evaluation(d,false,i)-Evaluation(d,true,i) for i=0:k]

for op in (:rdirichlet,:ldirichlet,:dirichlet,:lneumann,:rneumann,:neumann,:ivp,:bvp)
    @eval begin
        $op()=$op(UnsetSpace())
        $op(::PeriodicDomain)=[]
    end
end

for op in (:ldiffbc,:rdiffbc,:diffbcs,:ivp,:bvp,:periodic)
    @eval $op(k::Integer)=$op(UnsetSpace(),k)
end

immutable Dirichlet{S,T} <: Operator{T}
    space::S
    order::Int
end
Dirichlet(sp::Space)=Dirichlet{typeof(sp),BandedMatrix{eltype(sp)}}(sp,0)
Dirichlet(d::Domain)=Dirichlet(Space(d))
Neumann(sp::Space)=Dirichlet{typeof(sp),BandedMatrix{eltype(sp)}}(sp,1)
Neumann(d::Domain)=Dirichlet(Space(d))

domainspace(S::Dirichlet)=S.space
rangespace(B::Dirichlet)=Space(∂(domain(B)))


text:
import ArchGDAL, GDAL;
const AG = ArchGDAL
using Test

@testset "test_gdalutilities.jl" begin
    AG.read("data/utmsmall.tif") do ds_small
        @testset "GDAL Error" begin
            @test_throws GDAL.GDALError AG.gdalinfo(
                ds_small,
                ["-novalidoption"],
            )
            @test_throws GDAL.GDALError AG.unsafe_gdaltranslate(
                ds_small,
                ["-novalidoption"],
            )
            @test_throws GDAL.GDALError AG.unsafe_gdalbuildvrt(
                [ds_small],
                ["-novalidoption"],
            )
            @test_throws GDAL.GDALError AG.unsafe_gdaldem(
                ds_small,
                "hillshade",
                ["-novalidoption"],
            )
            @test_throws GDAL.GDALError AG.unsafe_gdalnearblack(
                ds_small,
                ["-novalidoption"],
            )
            @test_throws GDAL.GDALError AG.unsafe_gdalwarp(
                [ds_small],
                ["-novalidoption"],
            )
        end

        @testset "GDAL Info" begin
            infostr = AG.gdalinfo(ds_small, ["-checksum"])
            @test occursin("Checksum=50054", infostr)
            info_default = AG.gdalinfo(ds_small)
            @test occursin("Driver: GTiff/GeoTIFF", info_default)
        end

        AG.gdaltranslate(
            ds_small, # resample to a 5×5 ascii grid
            ["-of", "AAIGrid", "-r", "cubic", "-tr", "1200", "1200"],
        ) do ds_tiny
            @testset "GDAL Translate" begin
                @test AG.read(ds_tiny, 1) == [
                    128 171 127 93 83
                    126 164 148 114 101
                    161 175 177 164 140
                    185 206 205 172 128
                    193 205 209 181 122
                ]
            end

            @testset "GDAL Build VRT" begin
                AG.gdalbuildvrt([ds_tiny]) do ds_vrt
                    @test AG.read(ds_vrt, 1) == [
                        128 171 127 93 83
                        126 164 148 114 101
                        161 175 177 164 140
                        185 206 205 172 128
                        193 205 209 181 122
                    ]
                end
            end

            @testset "GDAL DEM Processing" begin
                AG.gdaldem(ds_tiny, "hillshade", ["-of", "AAIGrid"]) do ds_dempr
                    @test AG.read(ds_dempr, 1) == [
                        0 0 0 0 0
                        0 183 180 181 0
                        0 184 182 181 0
                        0 183 181 177 0
                        0 0 0 0 0
                    ]
                end
                AG.gdaldem(
                    ds_tiny,
                    "color-relief",
                    colorfile = "data/color_relief.txt",
                ) do ds_dempr
                    @test AG.read(ds_dempr, 1) == [
                        0x80 0x87 0x80 0x7b 0x7a
                        0x80 0x86 0x83 0x7e 0x7d
                        0x85 0x87 0x88 0x86 0x82
                        0x89 0x8c 0x8c 0x87 0x80
                        0x8a 0x8c 0x8c 0x88 0x80
                    ]
                end
            end

            @testset "GDAL Near Black" begin
                AG.gdalnearblack(
                    ds_tiny,
                    ["-of", "GTiff", "-color", "0"],
                ) do ds_nearblack
                    @test AG.read(ds_nearblack, 1) == [
                        0 0 0 0 0
                        0 0 0 0 0
                        0 0 177 0 0
                        0 0 0 0 0
                        0 0 0 0 0
                    ]
                end
            end
        end

        @testset "GDAL Warp" begin
            AG.gdalwarp(
                [ds_small],
                ["-of", "MEM", "-t_srs", "EPSG:4326"],
            ) do ds_warped
                @test AG.width(ds_small) == 100
                @test AG.height(ds_small) == 100
                @test AG.width(ds_warped) == 109
                @test AG.height(ds_warped) == 91
            end
        end
        @testset "GDAL Warp" begin
            AG.gdalwarp([ds_small], ["-of", "MEM"]) do ds_warped
                @test AG.width(ds_small) == 100
                @test AG.height(ds_small) == 100
                @test AG.shortname(AG.getdriver(ds_small)) == "GTiff"
                @test AG.width(ds_warped) == 100
                @test AG.height(ds_warped) == 100
                @test AG.shortname(AG.getdriver(ds_warped)) == "MEM"
            end
        end
    end

    AG.read("data/point.geojson") do ds_point
        @testset "GDAL Grid" begin
            AG.gdalgrid(
                ds_point,
                [
                    "-of",
                    "MEM",
                    "-outsize",
                    "3",
                    "10",
                    "-txe",
                    "100",
                    "100.3",
                    "-tye",
                    "0",
                    "0.1",
                ],
            ) do ds_grid
                @test AG.getgeotransform(ds_grid) ≈
                      [100.0, 0.1, 0.0, 0.0, 0.0, 0.01]
            end
        end

        @testset "GDAL Rasterize" begin
            AG.gdalrasterize(
                ds_point,
                ["-of", "MEM", "-tr", "0.05", "0.05"],
            ) do ds_rasterize
                @test AG.getgeotransform(ds_rasterize) ≈
                      [99.975, 0.05, 0.0, 0.1143, 0.0, -0.05]
            end
        end

        @testset "GDAL Vector Translate" begin
            AG.destroy(
                AG.unsafe_gdalvectortranslate(
                    [ds_point],
                    ["-f", "CSV", "-lco", "GEOMETRY=AS_XY"],
                    dest = "data/point.csv",
                ),
            )
            @test replace(read("data/point.csv", String), "\r" => "") == """
            X,Y,FID,pointname
            100,0,2,point-a
            100.2785,0.0893,3,point-b
            100,0,0,a
            100.2785,0.0893,3,b
            """
            rm("data/point.csv")
        end
    end
end

@testset "Interactive data/utmsmall.tif" begin
    ds_small = AG.read("data/utmsmall.tif")
    @testset "GDAL Error" begin
        @test_throws GDAL.GDALError AG.gdalinfo(ds_small, ["-novalidoption"])
        @test_throws GDAL.GDALError AG.unsafe_gdaltranslate(
            ds_small,
            ["-novalidoption"],
        )
        @test_throws GDAL.GDALError AG.unsafe_gdalbuildvrt(
            [ds_small],
            ["-novalidoption"],
        )
        @test_throws GDAL.GDALError AG.unsafe_gdaldem(
            ds_small,
            "hillshade",
            ["-novalidoption"],
        )
        @test_throws GDAL.GDALError AG.unsafe_gdalnearblack(
            ds_small,
            ["-novalidoption"],
        )
        @test_throws GDAL.GDALError AG.unsafe_gdalwarp(
            [ds_small],
            ["-novalidoption"],
        )
    end

    @testset "GDAL Info" begin
        infostr = AG.gdalinfo(ds_small, ["-checksum"])
        @test occursin("Checksum=50054", infostr)
        info_default = AG.gdalinfo(ds_small)
        @test occursin("Driver: GTiff/GeoTIFF", info_default)
    end

    ds_tiny = AG.unsafe_gdaltranslate(
        ds_small, # resample to a 5×5 ascii grid
        ["-of", "AAIGrid", "-r", "cubic", "-tr", "1200", "1200"],
    )
    @test typeof(ds_tiny) == AG.Dataset
    @testset "GDAL Translate" begin
        @test AG.read(ds_tiny, 1) == [
            128 171 127 93 83
            126 164 148 114 101
            161 175 177 164 140
            185 206 205 172 128
            193 205 209 181 122
        ]
    end

    @testset "GDAL Build VRT" begin
        ds_vrt = AG.unsafe_gdalbuildvrt([ds_tiny])
        @test AG.read(ds_vrt, 1) == [
            128 171 127 93 83
            126 164 148 114 101
            161 175 177 164 140
            185 206 205 172 128
            193 205 209 181 122
        ]
    end

    @testset "GDAL DEM Processing" begin
        ds_dempr = AG.unsafe_gdaldem(ds_tiny, "hillshade", ["-of", "AAIGrid"])
        @test AG.read(ds_dempr, 1) == [
            0 0 0 0 0
            0 183 180 181 0
            0 184 182 181 0
            0 183 181 177 0
            0 0 0 0 0
        ]
    end

    @testset "GDAL Near Black" begin
        ds_nearblack =
            AG.unsafe_gdalnearblack(ds_tiny, ["-of", "GTiff", "-color", "0"])
        @test AG.read(ds_nearblack, 1) == [
            0 0 0 0 0
            0 0 0 0 0
            0 0 177 0 0
            0 0 0 0 0
            0 0 0 0 0
        ]
    end

    @testset "GDAL Warp" begin
        AG.gdalwarp(
            [ds_small],
            ["-of", "MEM", "-t_srs", "EPSG:4326"],
        ) do ds_warped
            @test AG.width(ds_small) == 100
            @test AG.height(ds_small) == 100
            @test AG.width(ds_warped) == 109
            @test AG.height(ds_warped) == 91
        end
    end

    @testset "GDAL Warp #2" begin
        ds_warped = AG.unsafe_gdalwarp([ds_small], ["-of", "MEM"])
        @test AG.width(ds_small) == 100
        @test AG.height(ds_small) == 100
        @test AG.shortname(AG.getdriver(ds_small)) == "GTiff"
        @test AG.width(ds_warped) == 100
        @test AG.height(ds_warped) == 100
        @test AG.shortname(AG.getdriver(ds_warped)) == "MEM"
    end
end