text:
unsafe_load(pointer.(strs) + 3)

a = "def"
b  = unsafe_string(pointer(a))

pointer(a)
pointer(b)

a = UInt8[64, 65, C_NULL, 62, 53]
pp= pointer(a)
print(pp)
@time as = String(a)
pointer(as)

unsafe_wrap(a, Ptr{String})

a = "def"
unsafe_load(pointer(a)+3)

x = "id".*string.(rand(UInt8, 1_000_000))

io = open("c:/data/io.bin", "w")
ncw = write.(Ref(io), x)
close(io)

io = open("c:/data/io.bin", "r")
buffer = Vector{UInt8}(undef, sum(ncw))
readbytes!(io, buffer, sum(ncw))
close(io)

@time aa = test(buffer)

function test(array)
    start = 1
    strings = String[]
    GC.@preserve array begin
        ptr = pointer(array) - 1
        for i in eachindex(array)
            @inbounds char = array[i]
            if char == UInt8(',')
                len = i - start
                str = unsafe_string(ptr + start, len)
                push!(strings, str)
                start = i + 1
            end
        end
    end
    strings
end

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
using MolecularDynamicsIntegrators

include("./common.jl")
include("./verlet.jl")
include("./rattle.jl")

ext: jl


text:
using Base.Test

push!(LOAD_PATH,"src")

using DalitzPlotAnalysis

fewtests = [
    let s = 1.4, mπ = 0.14, mπ2 = mπ^2, z=2rand()-1
        σ1=4mπ2+((√s-mπ)^2-4mπ2)*rand();
        τ = change_basis(σ1,1.0,0.0,z,0.0,mπ2,mπ2,mπ2,s)
        cosθ3 = cross_basis_cosθ3(σ1,z,mπ2,mπ2,mπ2,s)
        # that is the test
        τ[2] ≈ cosθ3 #
    end
    for i in 1:1000]

@test sum(fewtests) == size(fewtests,1)


text:
import ..speed_of_light_mps

const harpsn_all_orders = 1:69
const harpsn_flux_hdu = 1

const harpsn_rv_orders = vcat(collect(1:54), collect(56:60), collect(62:63), collect(66:68) )

default_blaze_degree = 4
max_blaze_degree = 6
default_pixel_buffer_lo = 64
default_pixel_buffer_hi = 8

#const harpsn_var_hdu = 5   # computed from flux and ron
#const harpsn_lambda_hdu = 8 # computed from polynomial coefficients in heard

#= Since in CSV w/ header, don't need to specify columns
const harpsn_manifest_format = [
		    ManifestFormatEntry(:UTTimestamp,String, 4:26),
                    ManifestFormatEntry(:UTDate,String,     30:37),
                    ManifestFormatEntry(:ObsNum,Int,        43:46),
                    ManifestFormatEntry(:Frame,String,      50:89),
                    ManifestFormatEntry(:ObsType,String,    97:99),
                    ManifestFormatEntry(:iTime,Float64,    100:117),
                    ManifestFormatEntry(:QProg,String,     118:145),
                    ManifestFormatEntry(:Object,String,    146:190)
                    ]
=#


text:
##### Diffusion approximation for Slab Geometry (transmittance) #####
# The transmittance for the slab is calculated from Eqn. 39 from Contini 1997 for the first 11 dipoles (-10:1:10).

function DT_transslab(t, β::Array{Float64,1}, ρ::Float64,ndet::Float64, nmed::Float64, s::Float64)
	#s is slab thickness
    n::Float64 = nmed/ndet
    μa::Float64 = β[1]
    μsp::Float64 = β[2]
    D::Float64 = 1/3μsp
	ν::Float64 = 29.9792345/nmed
	xs::UnitRange{Int64} = -10:10

	z1m = Array{Float64}(undef, length(xs))
	z2m = Array{Float64}(undef, length(xs))

    Rt1 = Array{Float64}(undef, length(t))
    Rt2 = Array{Float64}(undef, length(t), length(z1m))

	if n == 1.0
		A = 1.0
	elseif n > 1.0
		A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		 5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
	else
		A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
	end

    zs::Float64 = 1/μsp
    ze::Float64 = 2A*D

    z1m = Float64[(1 .- 2m).*s .- 4m.*ze .- zs for m in xs]
    z2m = Float64[(1 .- 2m).*s .- (4m .- 2).*ze .+ zs for m in xs]

    Rt1 = @. exp(-(ρ^2/(4D*ν*t)) - μa*ν*t)
    Rt1 = @. Rt1/(2*(4π*D*ν)^(3/2)*t^(5/2))

    Rt2 = @. z1m'*exp(-(z1m'^2 / (4D*ν*t))) - z2m'*exp(-(z2m'^2 / (4D*ν*t)))

    Rt = Rt1.*sum(Rt2, dims=2)
	replace!(Rt, NaN => 0)

    return Rt
end

##### Diffusion approximation for a turbid parallelepiped (transmittance)  #####

#=
@article{kienle2005light,
  title={Light diffusion through a turbid parallelepiped},
  author={Kienle, Alwin},
  journal={JOSA A},
  volume={22},
  number={9},
  pages={1883--1888},
  year={2005},
  publisher={Optical Society of America}
}

The transmitance is directly given by equation 8 in Kienle 2005.
=#

function DT_trans_paralpip(t, β::Array{Float64,1}, ndet::Float64, nmed::Float64, rd::Array{Float64,1}, rs::Array{Float64,1}, L::Array{Float64,1})
	# rd is location of detector [x,y]
	# rs is location of source [x,y]
	# L is dimensions of parallelpiped [Lx, Ly, Lz]

	n::Float64 = nmed/ndet
    μa::Float64 = β[1]
    μsp::Float64 = β[2]
    D::Float64 = 1/3μsp
	ν::Float64 = 29.9792345/nmed
	xs::UnitRange{Int64} = -10:10

	x::Float64 = rd[1]
	y::Float64 = rd[2]
	z::Float64 = L[3]
	xu::Float64 = rs[1]
	yu::Float64 = rs[2]
	lx::Float64 = L[1]
	ly::Float64 = L[2]
	lz::Float64 = L[3]

	x1l = Array{Float64}(undef, length(xs))
	x2l = Array{Float64}(undef, length(xs))
	y1m = Array{Float64}(undef, length(xs))
	y2m = Array{Float64}(undef, length(xs))
	z1n = Array{Float64}(undef, length(xs))
	z2n = Array{Float64}(undef, length(xs))

	Rt1 = Array{Float64}(undef, length(t))
	Rt2 = Array{Float64}(undef, length(t), length(xs))
	Rt3 = Array{Float64}(undef, length(t), length(xs))
	Rt4 = Array{Float64}(undef, length(t), length(xs))

	if n == 1.0
		A = 1.0
	elseif n > 1.0
		A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		 5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
	else
		A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
	end

	zo::Float64 = 1/μsp
	zb::Float64 = 2A*D

	x1l = Float64[2l*lx + 4l*zb + xu for l in xs]
	x2l = Float64[2l*lx + (4l-2)*zb - xu for l in xs]
	y1m = Float64[2m*ly + 4m*zb + yu for m in xs]
	y2m = Float64[2m*ly + (4m-2)*zb - yu for m in xs]
	z1n = Float64[2n*lz + 4n*zb + zo for n in xs]
	z2n = Float64[2n*lz + (4n-2)*zb - zo for n in xs]

	Rt1 = @. exp(-μa*ν*t)/(2*(4π*D*ν)^(3/2)*t^(5/2))

	Rt2 = @. exp(-(x-x1l')^2 / (4D*ν*t)) - exp(-(x-x2l')^2 / (4D*ν*t))

	Rt3 = @. exp(-(y-y1m')^2 / (4D*ν*t)) - exp(-(y-y2m')^2 / (4D*ν*t))

	Rt4 = @. (lz - z1n')*exp(-(lz-z1n')^2 /(4D*ν*t)) - (lz - z2n')*exp(-(lz-z2n')^2 /(4D*ν*t))

	Rt = Rt1.*sum(Rt2, dims=2).*sum(Rt3, dims=2).*sum(Rt4, dims=2)
	replace!(Rt, NaN =>0)
	return Rt
end

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
using .ArcadeLearningEnvironment
using .ReinforcementLearningEnvironments
using BSON
using Flux:Chain

function atari_env_factory(
    name,
    state_size,
    n_frames,
    max_episode_steps = 100_000;
    seed = nothing,
    repeat_action_probability = 0.25,
    n_replica = 1
)
    init(seed) = AtariEnv(;
        name = string(name),
        grayscale_obs = true,
        noop_max = 30,
        frame_skip = 4,
        terminal_on_life_loss = false,
        repeat_action_probability = repeat_action_probability,
        max_num_frames_per_episode = n_frames * max_episode_steps,
        color_averaging = false,
        full_action_space = false,
        seed = seed,
    ) |>
    env -> StateOverriddenEnv(
        env,
        Chain(ResizeImage(state_size...), StackFrames(state_size..., n_frames))
    ) |>
    StateCachedEnv |>
    env -> RewardOverriddenEnv(env, r -> clamp(r, -1, 1))

    if n_replica == 1
        init(seed)
    else
        envs = [init(hash(seed+i)) for i in 1:n_replica]
        states = Flux.batch(state.(envs))
        rewards = reward.(envs)
        terminals = is_terminated.(envs)
        A = Space([action_space(x) for x in envs])
        S = Space(fill(0..255, size(states)))
        MultiThreadEnv(envs, states, rewards, terminals, A, S, nothing)
    end
end

"Total reward per episode before reward reshaping"
Base.@kwdef mutable struct TotalOriginalRewardPerEpisode <: AbstractHook
    rewards::Vector{Float64} = Float64[]
    reward::Float64 = 0.0
end

function (hook::TotalOriginalRewardPerEpisode)(::PostActStage, agent, env::RewardOverriddenEnv)
    hook.reward += reward(env.env)
end

function (hook::TotalOriginalRewardPerEpisode)(::PostEpisodeStage, agent, env)
    push!(hook.rewards, hook.reward)
    hook.reward = 0
end

"Total reward of each inner env per episode before reward reshaping"
struct TotalBatchOriginalRewardPerEpisode <: AbstractHook
    rewards::Vector{Vector{Float64}}
    reward::Vector{Float64}
end

function TotalBatchOriginalRewardPerEpisode(batch_size::Int)
    TotalBatchOriginalRewardPerEpisode([Float64[] for _ in 1:batch_size], zeros(batch_size))
end

function (hook::TotalBatchOriginalRewardPerEpisode)(::PostActStage, agent, env::MultiThreadEnv{<:RewardOverriddenEnv})
    for (i, e) in enumerate(env.envs)
        hook.reward[i] += reward(e.env)
        if is_terminated(e)
            push!(hook.rewards[i], hook.reward[i])
            hook.reward[i] = 0.0
        end
    end
end

for f in readdir(@__DIR__)
    if f != splitdir(@__FILE__)[2]
        include(f)
    end
end


text:
using ContinuousPOMDPTreeSearchExperiments
using ParticleFilters
using ARDESPOT
using BasicPOMCP
using POMCPOW
using POMDPs
using DiscreteValueIteration
using QMDP
using MCTS
using SubHunt
using POMDPToolbox
using DataFrames
using CSV

file_contents = readstring(@__FILE__())

pomdp = SubHuntPOMDP()

vs = ValueIterationSolver()
if !isdefined(:vp) || vp.mdp != pomdp
    vp = solve(vs, pomdp, verbose=true)
end
qs = QMDPSolver()
qp = QMDP.create_policy(qs, pomdp)
qp.alphas[:] = vp.qmat

@show max_time = 1.0
@show max_depth = 20

solvers = Dict{String, Union{Solver,Policy}}(
    "qmdp" => qp,
    "ping_first" => PingFirst(qp),

    # parameter adjustments from 1/25 Cross Entropy
    "pomcpow" => begin
        rng = MersenneTwister(13)
        solver = POMCPOWSolver(tree_queries=10_000_000_000,
                               criterion=MaxUCB(17.0),
                               final_criterion=MaxTries(),
                               max_depth=max_depth,
                               max_time=max_time,
                               enable_action_pw=false,
                               k_observation=6.0,
                               alpha_observation=1/100.0,
                               estimate_value=FOValue(vp),
                               check_repeat_obs=false,
                               default_action=ReportWhenUsed(qp),
                               tree_in_info=false,
                               rng=rng
                              )
    end,

    "pft" => begin
        rng = MersenneTwister(13)
        m = 20
        solver = DPWSolver(n_iterations=typemax(Int),
                           exploration_constant=100.0,
                           depth=max_depth,
                           max_time=max_time,
                           k_state=2.0,
                           alpha_state=1/10,
                           check_repeat_state=false,
                           check_repeat_action=false,
                           estimate_value=RolloutEstimator(qp),
                           enable_action_pw=false,
                           # default_action=ReportWhenUsed(qp),
                           tree_in_info=false,
                           rng=rng
                          )

        # node_updater = ObsAdaptiveParticleFilter(deepcopy(pomdp),
        #                                    LowVarianceResampler(m),
        #                                    0.1, rng)
        # belief_mdp = GenerativeBeliefMDP(deepcopy(pomdp), node_updater)

        belief_mdp = MeanRewardBeliefMDP(pomdp,
                                         LowVarianceResampler(m),
                                         0.1
                                        )

        solve(solver, belief_mdp)
    end,

    "despot" => begin
        rng = MersenneTwister(13)
        b = IndependentBounds(DefaultPolicyLB(qp), 100.0, check_terminal=true)
        DESPOTSolver(lambda=0.01,
                     epsilon_0=0.0,
                     K=500,
                     D=max_depth,
                     max_trials=1_000_000,
                     T_max=max_time,
                     bounds=b,
                     default_action=ReportWhenUsed(1),
                     rng=rng)
    end,

    "d_despot" => begin
        rng = MersenneTwister(13)
        b = IndependentBounds(DefaultPolicyLB(qp), 100.0, check_terminal=true)
        dpomdp = DSubHuntPOMDP(pomdp, 1.0)
        sol = DESPOTSolver(lambda=0.01,
                     epsilon_0=0.0,
                     K=500,
                     D=max_depth,
                     max_trials=1_000_000,
                     T_max=max_time,
                     bounds=b,
                     default_action=ReportWhenUsed(1),
                     rng=rng)
        solve(sol, dpomdp)
    end,

    "pomcp" => begin
        rng = MersenneTwister(13)
        POMCPSolver(max_depth=max_depth,
                    max_time=max_time,
                    c=100.0,
                    tree_queries=typemax(Int),
                    default_action=ReportWhenUsed(qp),
                    estimate_value=FOValue(vp),
                    tree_in_info=false,
                    rng=rng
                   )
    end,

    "d_pomcp" => begin
        rng = MersenneTwister(13)
        dpomdp = DSubHuntPOMDP(pomdp, 1.0)
        sol = POMCPSolver(max_depth=max_depth,
                    max_time=max_time,
                    c=100.0,
                    tree_queries=typemax(Int),
                    default_action=ReportWhenUsed(qp),
                    estimate_value=FOValue(vp),
                    tree_in_info=false,
                    rng=rng
                   )
        solve(sol, dpomdp)
    end,

    "pomcpdpw" => begin
        rng = MersenneTwister(13)
        solver = PDPWSolver(tree_queries=10_000_000,
                            c=100.0,
                            max_depth=max_depth,
                            max_time=max_time,
                            enable_action_pw=false,
                            k_observation=1.0,
                            alpha_observation=1/10,
                            estimate_value=FOValue(vp),
                            check_repeat_obs=false,
                            # default_action=ReportWhenUsed(-1),
                            rng=rng
                           )
    end,
)

@show N=1000

alldata = DataFrame()
for (k, solver) in solvers
# test = ["pft"]
# for (k, solver) in [(s, solvers[s]) for s in test]
    @show k
    if isa(solver, Solver)
        planner = solve(solver, pomdp)
    else
        planner = solver
    end
    sims = []
    for i in 1:N
        srand(planner, i+50_000)
        filter = SIRParticleFilter(deepcopy(pomdp), 100_000, rng=MersenneTwister(i+90_000))

        md = Dict(:solver=>k, :i=>i)
        sim = Sim(deepcopy(pomdp),
                  planner,
                  filter,
                  rng=MersenneTwister(i+70_000),
                  max_steps=100,
                  metadata=md
                 )

        push!(sims, sim)
    end

    data = run_parallel(sims)
    # data = run(sims)

    rs = data[:reward]
    println(@sprintf("reward: %6.3f ± %6.3f", mean(rs), std(rs)/sqrt(length(rs))))
    if isempty(alldata)
        alldata = data
    else
        alldata = vcat(alldata, data)
    end
end

datestring = Dates.format(now(), "E_d_u_HH_MM")
copyname = Pkg.dir("ContinuousPOMDPTreeSearchExperiments", "icaps_2018", "data", "subhunt_table_$(datestring).jl")
write(copyname, file_contents)
filename = Pkg.dir("ContinuousPOMDPTreeSearchExperiments", "icaps_2018", "data", "subhunt_$(datestring).csv")
println("saving to $filename...")
CSV.write(filename, alldata)
println("done.")


text:
function getLaplacianKernel{T <: FloatingPoint}(:: Type{T})
        nType = T == Float64 ? "double" : "float"

        return "
        #if defined(cl_khr_fp64)  // Khronos extension available?
        #pragma OPENCL EXTENSION cl_khr_fp64 : enable
        #elif defined(cl_amd_fp64)  // AMD extension available?
        #pragma OPENCL EXTENSION cl_amd_fp64 : enable
        #endif
        #define number $nType

        #define Ap(x,y) a[y*D2 + x]

        #define A_lap(x,y) out[y*D2 + x]

        __kernel void laplacian(
                      __global const number *a,
                      __global number *out,
                      const int D1,
                      const int D2) {

            int i = get_global_id(0);
            int j = get_global_id(1);
            int west = j-1;
            int east = j+1;
            int north = i+1;
            int south = i-1;

            if(j == 0)
                west = D2 - 1;
            if(j == D2 - 1)
                east = 0;
            if(i == D1 - 1)
                north = 0;
            if(i == 0)
                south = D1 - 1;

            // Calulate the laplacian

            A_lap(i,j) = -Ap(i,j) +(Ap(north, east) +
                                    Ap(north, j   ) +
                                    Ap(north, west) +
                                    Ap(i    , west) +
                                    Ap(south, west) +
                                    Ap(south, j   ) +
                                    Ap(south, east) +
                                    Ap(i    , east))/8;
    }
"
end

function laplacianCL!(
    a_buff,
    out_buff,
    d1 :: Int64, d2 :: Int64,
    ctx, queue, program)

    k = Kernel(program, "laplacian")

    call(queue, k, (d1,d2), nothing, a_buff, out_buff, int32(d1), int32(d2))
end


text:
tol = 1e-13

# check boundaries
N = 5
M = 3
(p,t,bel) = luteos.makecube( M, N, M )

indchecks = [3, 1, 3, 1, 2, 2]
checkval  = [0.0, 1.0, 1.0, 0.0, 0.0, 1.0]

for jj in 1:size(bel,1)
  ic   = indchecks[ bel[jj,4] ]
  cval = checkval[ bel[jj,4] ]
  for ii in 1:3
    @test p[ bel[ jj, ii ], ic ] == cval
  end
end

# check connectivity
for ii in [2,5], jj in [2,5], kk in [2,5], pp in 1:3
  mesh   = Mesh3D( "cube", pp, N = ii, M = jj, Q = kk )
  master = Master3D( pp )
  for ff in 1:size(mesh.f,3)

    tr = mesh.f[ff,4]
    tl = mesh.f[ff,5]

    if tl > 0

      itrf = find( mesh.t2f[tr,1:4] .== ff )
      itlf = find( mesh.t2f[tl,1:4] .== ff )

      nodr = master.perm[ :,itrf,abs.(mesh.t2f[tr,itrf+4]) ]
      nodl = master.perm[ :,itlf,abs.(mesh.t2f[tl,itlf+4]) ]

      for qq in 1:length(nodr), dd in 1:3
        @test abs.(mesh.nodes[ nodr[qq], dd, tr ] -
                   mesh.nodes[ nodl[qq], dd, tl ]) < tol
      end

    else
      # check that nodes are actually on boundary they should be on
      itrf =   find( mesh.t2f[tr,1:4] .== ff )
      nodr =   master.perm[ :,itrf,abs.(mesh.t2f[tr,itrf+4]) ]
      nb   = - mesh.f[ff,end]
      for qq in 1:length(nodr)
        ic   = indchecks[ nb ]
        cval = checkval[  nb ]
        @test mesh.nodes[ nodr[qq], ic, tr ] == cval
      end

    end

  end
end

# check integration through volume integrals
for ii in [2,5,7,10], jj in [2,5,7,10], kk in [2,5,7,10], pp in 1:3

  mesh = Mesh3D( "cube", pp, N = ii, M = jj, Q = kk )
  master = Master3D( pp )

  compJacob!( mesh, master )

  vol = 0

  for el in 1:size(mesh.nodes,3)
    vol += sum( mesh.jcw[:,el] )
  end

  @test abs(vol - 1) < tol

end

# check integration through surface integrals
for ii in [5,10], jj in [2,10], kk in [2,10], pp in 1:3

  mesh = Mesh3D( "cube", pp, N = ii, M = jj, Q = kk )
  master = Master3D( pp )
  compJacob!( mesh, master )

  vol = 0

  for el in 1:size(mesh.nodes,3)

    # using divergence theorem: ∭ᵥ (∇⋅F) dV = ∯ₛ(F⋅n)dS
    # In this case F = 1/3 [x,y,z]

    voltemp = 0

    for qq in 1:4

      indF = mesh.t2f[el,qq]

      (ϕ2d,p2d,nods,normal,jcw2D) = luteos.compJacobFace( mesh, master, el, qq )

      volttemp = 0.0
      for ll in 1:size(p2d,1)
        volttemp += 1/3 * jcw2D[ll] * (normal[ll,:]' * p2d[ll,:])
      end
      voltemp += volttemp
    end

    vol += voltemp

  end

  @test abs(vol - 1) < tol

end

### Check skewed mesh
# check integration through volume integrals
for ii in [2,5,7,10], jj in [2,5,7,10], kk in [2,5,7,10], pp in 1:3

  mesh = Mesh3D( "skewcube", pp, N = ii, M = jj, Q = kk )
  master = Master3D( pp )

  compJacob!( mesh, master )

  vol = 0

  for el in 1:size(mesh.nodes,3)
    vol += sum( mesh.jcw[:,el] )
  end

  @test abs(vol - 1) < tol

end

# check integration through surface integrals
for ii in [5,10], jj in [2,10], kk in [2,10], pp in 1:3

  mesh = Mesh3D( "skewcube", pp, N = ii, M = jj, Q = kk )
  master = Master3D( pp )
  compJacob!( mesh, master )

  vol = 0

  for el in 1:size(mesh.nodes,3)

    # using divergence theorem: ∭ᵥ (∇⋅F) dV = ∯ₛ(F⋅n)dS
    # In this case F = 1/3 [x,y,z]

    voltemp = 0

    for qq in 1:4

      indF = mesh.t2f[el,qq]

      (ϕ2d,p2d,nods,normal,jcw2D) = luteos.compJacobFace( mesh, master, el, qq )

      volttemp = 0.0
      for ll in 1:size(p2d,1)
        volttemp += 1/3 * jcw2D[ll] * (normal[ll,:]' * p2d[ll,:])
      end
      voltemp += volttemp
    end

    vol += voltemp

  end

  @test abs(vol - 1) < tol

end

# check that basis functions and gauss points are correct
# for ii in [2,5,7,10], jj in [2,5,7,10], kk in [2,5,7,10], pp in 1:3
#
#   mesh = Mesh2D( "square", pp, N = ii )
#   master = Master2D( pp )
#
#   compJacob!( mesh, master )
#
#   barycent = fill( 0.0, size(master.gpts,1), 3 )
#   barycent[:,1]   = sum( master.gpts, 2 )
#   barycent[:,2:3] = master.gpts
#
#   for el in 1:size(mesh.nodes,3), ig in 1:size(master.gpts,1), dd in 1:2
#     cgpt = barycent[ig,1] * mesh.p[mesh.t[el,1],dd] +
#            barycent[ig,2] * mesh.p[mesh.t[el,2],dd] +
#            barycent[ig,3] * mesh.p[mesh.t[el,3],dd]
#     cgpt_ϕ = master.ϕ[:,ig]' * mesh.nodes[:,dd,el]
#     @test abs(cgpt_ϕ - cgpt) < tol
#   end
#
# end

# check integration through volume integrals using basis functions
# ∭ᵥ x dV
for ii in [2,5,7,10], jj in [2,5,7,10], kk in [2,5,7,10], pp in 1:3

  mesh = Mesh3D( "skewcube", pp, N = ii, M = jj, Q = kk )
  master = Master3D( pp )

  compJacob!( mesh, master )

  vol1 = 0
  vol2 = 0
  vol3 = 0

  for el in 1:size(mesh.nodes,3)
    vol1 += sum( diagm(mesh.jcw[:,el]) * (master.ϕ' * mesh.nodes[:,1,el] ) )
    vol2 += sum( diagm(mesh.jcw[:,el]) * (master.ϕ' * mesh.nodes[:,2,el] ) )
    vol3 += sum( diagm(mesh.jcw[:,el]) * (master.ϕ' * mesh.nodes[:,3,el] ) )
  end

  @test abs(vol1 - 1)   < tol
  @test abs(vol2 - 1)   < tol
  @test abs(vol3 - 0.5) < tol

end

# check that derivatives of basis functions are correct
for ii in [2,5,7,10], pp in 1:3
  mesh = Mesh2D( "square", pp, N = ii )
  master = Master2D( pp )

  compJacob!( mesh, master )

  for el in 1:size(mesh.nodes,3)
    ∇ϕc   = luteos.getderbfel( master, mesh.∂ξ∂x[:,el,:] )
    func  = mesh.nodes[:,1,el]
    res1  = ∇ϕc[:,:,1]' * func
    res2  = master.ϕ' * mesh.nodes[:,1,el].^0
    for cc in 1:length(res2)
      @test abs(res1[cc] - res2[cc]) < tol
    end

    func  = mesh.nodes[:,2,el]#0.5 * mesh.nodes[:,2,el].^2
    res1  = ∇ϕc[:,:,2]' * func
    res2  = master.ϕ' * mesh.nodes[:,2,el].^0
    for cc in 1:length(res2)
      @test abs(res1[cc] - res2[cc]) < tol
    end
  end

end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

module AstroLib
using Compat
using Base.Dates
import Compat.String

# Note on function definitions in this package.  Most functions are defined as
# follows:
#
#     _funcname{T<:AbstractFloat}(x::T, y::T) = ...
#     funcname(x::Real, y::Real) = _funcname(promote(float(x), float(y))...)
#
# "_funcname" is the function with the implementation of the full algorithm.
# "_funcname" doesn't have optional argument nor keywords.  Instead, "funcname"
# is the function exposed to the user which ensures all input arguments are
# converted to the same floating point type, in order to improve both
# type-stability and precision of calculations.  "funcname" may have optional
# arguments or keywords.  In addition, the use of AbstractFloat type in
# "_funcname" enables compatibility with Measurements.jl package, which defines
# a new type Measurement{T<:AbstractFloat} <: AbstractFloat.

# Core functions
include("misc.jl")
include("utils.jl")
# New types
include("types.jl")
# Common constants
include("common.jl")

end # module


text:
struct RandomProjectionTreeNode{T <: Number,
                                V <: AbstractVector{T}}
    indices
    isleaf::Bool
    hyperplane::Union{V, Nothing}
    offset::Union{T, Nothing}
    leftchild::Union{RandomProjectionTreeNode{T, V}, Nothing}
    rightchild::Union{RandomProjectionTreeNode{T, V}, Nothing}
end

struct RandomProjectionTree{T <: Number, V <: AbstractVector{T}}
    root::RandomProjectionTreeNode{T, V}
    leafsize
end

"""
    RandomProjectionTree(data; leafsize = 30)

A RandomProjectionTree is a binary tree whose non-leaf nodes correspond to random
hyperplanes in ℜᵈ and whose leaf nodes represent subregions that partition this
space.
"""
function RandomProjectionTree(data; leafsize = 30)
    root = build_rptree(data, 1:length(data), leafsize)
    return RandomProjectionTree(root, leafsize)
end

"""
    build_rptree(data::AbstractVector, indices, leafsize) -> RandomProjectionTreeNode

Recursively construct the RP tree by randomly splitting the data into nodes.
"""
function build_rptree(data::U, indices, leafsize) where {T <: Number,
                                                         V <: AbstractVector{T},
                                                         U <: AbstractVector{V}}
    if length(indices) <= leafsize
        return RandomProjectionTreeNode{T, V}(indices,
                                              true,
                                              nothing,
                                              nothing,
                                              nothing,
                                              nothing)
    end

    leftindices, rightindices, hyperplane, offset = rp_split(data, indices)
    leftchild = build_rptree(data, leftindices, leafsize)
    rightchild = build_rptree(data, rightindices, leafsize)

    return RandomProjectionTreeNode(nothing,
                                    false,
                                    hyperplane,
                                    offset,
                                    leftchild,
                                    rightchild)
end

"""
    search_tree(tree, point) -> node::RandomProjectionTreeNode

Search a random projection tree for the node to which a point belongs.
"""
function search_rptree(tree::RandomProjectionTree{T, V},
                       point::V
                       ) where {T, V}

    node = tree.root
    while !node.isleaf
        node = select_side(node.hyperplane, node.offset, point) ? node.leftchild : node.rightchild
    end
    return node

end

#############################
# Random Projection Forests
#############################

struct RandomProjectionForest{T, V}
    trees::Vector{RandomProjectionTree{T, V}}
end

"""
    RandomProjectionForest(data, args...; n_trees, kwargs...) -> RandomProjectionForest

Create a collection of `n_trees` random projection trees built using `data`.
`args...` and `kwargs...` are passed to the RandomProjectionTree constructor.
"""
function RandomProjectionForest(data, args...;
                                n_trees, kwargs...)
    # TODO: threading
    trees = [RandomProjectionTree(data, args...; kwargs...) for _ in 1:n_trees]
    return RandomProjectionForest(trees)
end

"""
    approx_knn(tree, data, point, n_neighbors) -> indices, distances

Find approximate nearest neighbors to `point` in `data` using a random
projection tree built on this data. Returns the indices and distances to the
approximate knns as arrays, sorted by distance.
"""
function approx_knn(tree::RandomProjectionTree,
                    data,
                    point,
                    n_neighbors)
    # TODO: handle case when n_neighbors > tree.leafsize
    candidates = search_rptree(tree, point).indices

    return _approx_knn(data, point, candidates, n_neighbors)
end

"""

"""
function approx_knn(forest::RandomProjectionForest,
                    data,
                    point,
                    n_neighbors)
    # TODO: threading
    candidates = []
    for tree in forest.trees
        union!(candidates, search_rptree(tree, point).indices)
    end
    return _approx_knn(data, point, candidates, n_neighbors)
end

"""

"""
function _approx_knn(data, point, candidates, n_neighbors)
    length(candidates) >= n_neighbors || @warn "Fewer candidates than n_neighbors!"
    # TODO: threading
    distances = [norm(data[i] - point) for i in candidates]
    perm = sortperm(distances)
    indices = candidates[perm]
    distances = distances[perm]

    if n_neighbors < length(indices)
        indices = indices[1:n_neighbors]
        distances = distances[1:n_neighbors]
    end
    return indices, distances
end


text:
#===================================================================================================
  Kernel Functions Module
===================================================================================================#

module MLKernels

import Base: call, convert, eltype, print, show, string, ==, *, /, +, -, ^, besselk, exp, gamma,
             tanh

export
    # Functions
    ∘,
    kernel,
    kernelmatrix,
    centerkernelmatrix!,
    centerkernelmatrix,

    ismercer,
    isnegdef,
    ispositive,
    isnonnegative,

    # Hyper Parameters
    Bound,
    Interval,
    leftbounded,
    rightbounded,
    unbounded,
    checkbounds,
    Variable,
    fixed,
    Argument,
    HyperParameter,

    # Kernel Type
    Kernel,
        StandardKernel,
            PairwiseKernel,
                AdditiveKernel,
                    SquaredDistanceKernel,
                    SineSquaredKernel,
                    ChiSquaredKernel,
                    ScalarProductKernel,
                ARD,
            KernelComposition,
        KernelOperation,
            KernelAffinity,
            KernelProduct,
            KernelSum,

    # Composition Classes
    CompositionClass,
        PositiveMercerClass,
            ExponentialClass,
            GammaExponentialClass,
            RationalClass,
            GammaRationalClass,
            MaternClass,
            ExponentiatedClass,
        PolynomialClass,
        NonNegativeNegativeDefiniteClass,
            PowerClass,
            LogClass,
            GammaLogClass,
        SigmoidClass,

    # Kernel Constructors
        GaussianKernel,
            SquaredExponentialKernel,
            RadialBasisKernel,
        LaplacianKernel,
        PeriodicKernel,
        RationalQuadraticKernel,
        MaternKernel,
            MatérnKernel,
        PolynomialKernel,
        LinearKernel,
        SigmoidKernel

include("hyperparameter.jl")
using MLKernels.HyperParameters: Bound, Interval, leftbounded, rightbounded, unbounded, checkbounds,
    Variable, fixed, Argument, HyperParameter

include("meta.jl")
include("kernel.jl")
include("pairwisefunction.jl")
include("kernelfunction.jl")

end # MLKernels


text:
# TODO: Make this more comprehensive. For now look at eval.jl or the add_command! docs.

module Commands

using Discord

function main()
    c = Client(ENV["DISCORD_TOKEN"])
    add_command!(c, :echo; pattern=r"^echo (.+)") do c, msg, noprefix
        reply(c, msg, noprefix)
    end
    open(c)
    return c
end

end

if abspath(PROGRAM_FILE) == @__FILE__
    c = Commands.main()
    wait(c)
end


text:
@testset "optimization" begin
    function f_lm(x)
      [x[1], 2.0 - x[2]]
    end
    function g_lm(x)
      [1.0 0.0; 0.0 -1.0]
    end

    initial_x = [100.0, 100.0]

    R1 = OnceDifferentiable(f_lm, g_lm, zeros(2), zeros(2); inplace = false)
    results = LsqFit.levenberg_marquardt(R1, initial_x)
    @assert norm(OptimBase.minimizer(results) - [0.0, 2.0]) < 0.01

    function rosenbrock_res(r, x)
        r[1] = 10.0 * (x[2] - x[1]^2 )
        r[2] =  1.0 - x[1]
        return r
    end

    function rosenbrock_jac(j, x)
        j[1, 1] = -20.0 * x[1]
        j[1, 2] =  10.0
        j[2, 1] =  -1.0
        j[2, 2] =   0.0
        return j
    end

    initial_xrb = [-1.2, 1.0]

    R2 = OnceDifferentiable(rosenbrock_res, rosenbrock_jac, zeros(2), zeros(2); inplace = true)

    results = LsqFit.levenberg_marquardt(R2, initial_xrb)

    @assert norm(OptimBase.minimizer(results) - [1.0, 1.0]) < 0.01

    # check estimate is within the bound PR #278
    result = LsqFit.levenberg_marquardt(R2, [150.0, 150.0]; lower = [10.0, 10.0], upper = [200.0, 200.0])
    @test OptimBase.minimizer(result)[1] >= 10.0
    @test OptimBase.minimizer(result)[2] >= 10.0

    # tests for #178, taken from LsqFit.jl, but stripped
    let
        Random.seed!(12345)

        # TODO: Change to `.-x` when 0.5 support is dropped
        model(x, p) = p[1] .* exp.(-x .* p[2])
        _nobs = 20
        xdata = range(0, stop = 10, length = _nobs)
        ydata = model(xdata, [1.0 2.0]) + 0.01*randn(_nobs)

        f_lsq = p -> model(xdata, p) - ydata
        g_lsq = p -> NLSolversBase.FiniteDiff.finite_difference_jacobian(f_lsq, p)

        R3 = OnceDifferentiable(f_lsq, g_lsq, zeros(2), zeros(_nobs); inplace = false)
        results = LsqFit.levenberg_marquardt(R3, [0.5, 0.5])

        @assert norm(OptimBase.minimizer(results) - [1.0, 2.0]) < 0.05
    end

    let
        Random.seed!(12345)

        # TODO: Change to `.-x` when 0.5 support is dropped
        model(x, p) = p[1] .* exp.(-x ./ p[2]) .+ p[3]

        xdata = 1:100
        p_generate = [10.0, 10.0, 10.0]
        ydata = model(xdata, p_generate) + 0.1*randn(length(xdata))

        # TODO: Change to `model.(xdata, p) .- ydata` when 0.4 support is dropped
        f_lsq = p -> model(xdata, p) - ydata
        g_lsq = p -> NLSolversBase.FiniteDiff.finite_difference_jacobian(f_lsq, p)

        R4 = OnceDifferentiable(f_lsq, g_lsq, similar(p_generate), zeros(length(xdata)); inplace = false)

        # tests for box constraints, PR #196
        @test_throws ArgumentError LsqFit.levenberg_marquardt(R4, [15.0, 15.0, 15.0], lower=[5.0, 11.0])
        @test_throws ArgumentError LsqFit.levenberg_marquardt(R4, [5.0, 5.0, 5.0], upper=[15.0, 9.0])
        @test_throws ArgumentError LsqFit.levenberg_marquardt(R4, [15.0, 10.0, 15.0], lower=[5.0, 11.0, 5.0])
        @test_throws ArgumentError LsqFit.levenberg_marquardt(R4, [5.0, 10.0, 5.0], upper=[15.0, 9.0, 15.0])

        lower=[5.0, 11.0, 5.0]
        results = LsqFit.levenberg_marquardt(R4, [15.0, 15.0, 15.0], lower=lower)
        OptimBase.minimizer(results)
        @test OptimBase.converged(results)
        @test all(OptimBase.minimizer(results) .>= lower)

        upper=[15.0, 9.0, 15.0]
        results = LsqFit.levenberg_marquardt(R4, [5.0, 5.0, 5.0], upper=upper)
        OptimBase.minimizer(results)
        @test OptimBase.converged(results)
        @test all(OptimBase.minimizer(results) .<= upper)

        # tests for PR #267
        LsqFit.levenberg_marquardt(R4, [15.0, 15.0, 15.0], show_trace=true)
    end
end


text:
"""

Carry out AHL21 mapping and calculates both the Jacobian and derivative with respect to the time step.
"""
function ahl21!(s::State{T},d::Derivatives{T},h::T) where T<:AbstractFloat
    zilch = zero(T); half::T = 0.5; two::T = 2.0;
    h2::T = half*h; n = s.n; sevn = 7*n; h6::T = h/6.0
    zero_out!(d)
    fill!(s.dqdt,zilch)
    kickfast!(s,d,h6)
    d.dqdt_kick ./= 6.0 # Since step is h/6
    # Since I removed identity from kickfast, need to add in dqdt:
    mul!(d.tmp7n, d.jac_kick, s.dqdt)
    s.dqdt .+= d.dqdt_kick .+ d.tmp7n #*(d.jac_kick,s.dqdt)
    # Multiply Jacobian from kick step:
    mul!(d.jac_copy, d.jac_kick, s.jac_step)
    drift_grad!(s,h2)
    # Compute time derivative of drift step:
    @inbounds for i=1:n, k=1:3
        s.dqdt[(i-1)*7+k] = half*s.v[k,i] + h2*s.dqdt[(i-1)*7+3+k]
    end
    # Add back in the identity portion of the Jacobian with compensated summation:
    comp_sum_matrix!(s.jac_step,s.jac_error,d.jac_copy)
    indi = 0; indj = 0
    @inbounds for i=1:s.n-1
        indi = (i-1)*7
        @inbounds for j=i+1:s.n
            indj = (j-1)*7
            if ~s.pair[i,j]  # Check to see if kicks have not been applied
                kepler_driftij_gamma!(s,d,i,j,h2,true)
                copy_submatrix!(s,d,indi,indj,sevn)
                # Carry out multiplication on the i/j components of matrix:
                mul!(d.jac_tmp2, d.jac_ij, d.jac_tmp1)
                # Add back in the Jacobian with compensated summation:
                comp_sum_matrix!(d.jac_tmp1,d.jac_err1,d.jac_tmp2)
                # Add in partial derivatives with respect to time:
                # Need to multiply by 1/2 since we're taking 1/2 time step:
                d.dqdt_ij .*= half
                mul!(d.tmp14, d.jac_ij, d.dqdt_tmp1)
                d.dqdt_ij .+= d.dqdt_tmp1 .+ d.tmp14#*(d.jac_ij,d.dqdt_tmp1)
                # Copy back time derivatives:
                ypoc_submatrix!(s,d,indi,indj,sevn)
            end
        end
    end
    phic!(s,d,h)
    phisalpha!(s,d,h,two)
    mul!(d.jac_copy, d.jac_phi, s.jac_step)
    # Add in time derivative with respect to prior parameters:
    # Copy result to dqdt:
    mul!(d.tmp7n, d.jac_phi, s.dqdt)
    s.dqdt .+= d.dqdt_phi .+ d.tmp7n #*(d.jac_phi,s.dqdt)
    # Add back in the identity portion of the Jacobian with compensated summation:
    comp_sum_matrix!(s.jac_step,s.jac_error,d.jac_copy)
    indi=0; indj=0
    @inbounds for i=s.n-1:-1:1
        indi=(i-1)*7
        @inbounds for j=s.n:-1:i+1
            indj=(j-1)*7
            if ~s.pair[i,j]  # Check to see if kicks have not been applied
                kepler_driftij_gamma!(s,d,i,j,h2,false)
                # Pick out indices for bodies i & j:
                # Carry out multiplication on the i/j components of matrix:
                copy_submatrix!(s,d,indi,indj,sevn)
                # Carry out multiplication on the i/j components of matrix:
                mul!(d.jac_tmp2,d.jac_ij,d.jac_tmp1)
                # Add back in the Jacobian with compensated summation:
                comp_sum_matrix!(d.jac_tmp1,d.jac_err1,d.jac_tmp2)
                # Add in partial derivatives with respect to time:
                # Need to multiply by 1/2 since we're taking 1/2 time step:
                d.dqdt_ij .*= half
                mul!(d.tmp14, d.jac_ij, d.dqdt_tmp1)
                d.dqdt_ij .+= d.dqdt_tmp1 .+ d.tmp14#*(d.jac_ij,d.dqdt_tmp1)
                # Copy back derivatives:
                ypoc_submatrix!(s,d,indi,indj,sevn)
            end
        end
    end
    drift_grad!(s,h2)
    # Compute time derivative of drift step:
    @inbounds for i=1:n, k=1:3
        s.dqdt[(i-1)*7+k] += half*s.v[k,i] + h2*s.dqdt[(i-1)*7+3+k]
    end
    fill!(d.dqdt_kick,zilch)
    kickfast!(s,d,h6)
    d.dqdt_kick ./= 6.0 # Since step is h/6
    # Copy result to dqdt:
    mul!(d.tmp7n, d.jac_kick, s.dqdt)
    s.dqdt .+= d.dqdt_kick .+ d.tmp7n#*(d.jac_kick,s.dqdt)
    # Multiply Jacobian from kick step:
    mul!(d.jac_copy,d.jac_kick,s.jac_step)
    # Add back in the identity portion of the Jacobian with compensated summation:
    comp_sum_matrix!(s.jac_step,s.jac_error,d.jac_copy)
    # Edit this routine to do compensated summation for Jacobian [x]
end

function ahl21!(s::State{T},d::Jacobian{T},h::T) where T<:AbstractFloat
    zilch = zero(T); uno = one(T); half = convert(T,0.5); two = convert(T,2.0); h2 = half*h; sevn = 7*s.n
    zero_out!(d)

    #drift!(s.x,s.v,s.xerror,s.verror,h2,s.n,s.jac_step,s.jac_error)
    #kickfast!(s.x,s.v,s.xerror,s.verror,h/6,s.m,s.n,d.jac_kick,d.dqdt_kick,pair)
    kickfast!(s,d,h/6)
    # Multiply Jacobian from kick step:
    if T == BigFloat
        d.jac_copy .= *(d.jac_kick,s.jac_step)
    else
        start = time()
        #BLAS.gemm!('N','N',uno,d.jac_kick,s.jac_step,zilch,d.jac_copy)
        mul!(d.jac_copy,d.jac_kick,s.jac_step)
        d.ctime[1] += time()-start
    end
    # Add back in the identity portion of the Jacobian with compensated summation:
    comp_sum_matrix!(s.jac_step,s.jac_error,d.jac_copy)
    drift_grad!(s,h2)
    indi = 0; indj = 0
    @inbounds for i=1:s.n-1
        indi = (i-1)*7
        @inbounds for j=i+1:s.n
            indj = (j-1)*7
            if ~s.pair[i,j]  # Check to see if kicks have not been applied
                kepler_driftij_gamma!(s,d,i,j,h2,true)
                #kepler_driftij_gamma!(s.m,s.x,s.v,s.xerror,s.verror,i,j,h2,d.jac_ij,d.dqdt_ij,true)
                # Pick out indices for bodies i & j:
                @inbounds for k2=1:sevn, k1=1:7
                    d.jac_tmp1[k1,k2] = s.jac_step[ indi+k1,k2]
                    d.jac_err1[k1,k2] = s.jac_error[indi+k1,k2]
                end
                @inbounds for k2=1:sevn, k1=1:7
                    d.jac_tmp1[7+k1,k2] = s.jac_step[ indj+k1,k2]
                    d.jac_err1[7+k1,k2] = s.jac_error[indj+k1,k2]
                end
                # Carry out multiplication on the i/j components of matrix:
                if T == BigFloat
                    d.jac_tmp2 .= *(d.jac_ij,d.jac_tmp1)
                else
                    start = time()
                    #BLAS.gemm!('N','N',uno,d.jac_ij,d.jac_tmp1,zilch,d.jac_tmp2)
                    mul!(d.jac_tmp2,d.jac_ij,d.jac_tmp1)
                    d.ctime[1] += time()-start
                end
                # Add back in the Jacobian with compensated summation:
                comp_sum_matrix!(d.jac_tmp1,d.jac_err1,d.jac_tmp2)
                # Copy back to the Jacobian:
                @inbounds for k2=1:sevn, k1=1:7
                    s.jac_step[ indi+k1,k2]=d.jac_tmp1[k1,k2]
                    s.jac_error[indi+k1,k2]=d.jac_err1[k1,k2]
                end
                @inbounds for k2=1:sevn, k1=1:7
                    s.jac_step[ indj+k1,k2]=d.jac_tmp1[7+k1,k2]
                    s.jac_error[indj+k1,k2]=d.jac_err1[7+k1,k2]
                end
            end
        end
    end
    #phic!(s.x,s.v,s.xerror,s.verror,h,s.m,s.n,d.jac_phi,d.dqdt_phi,pair)
    #phisalpha!(s.x,s.v,s.xerror,s.verror,h,s.m,two,s.n,d.jac_phi,d.dqdt_phi,pair) # 10%
    phic!(s,d,h)
    phisalpha!(s,d,h,two)
    if T == BigFloat
        d.jac_copy .= *(d.jac_phi,s.jac_step)
    else
        start = time()
        #BLAS.gemm!('N','N',uno,d.jac_phi,s.jac_step,zilch,d.jac_copy)
        mul!(d.jac_copy,d.jac_phi,s.jac_step)
        d.ctime[1] += time()-start
    end
    # Add back in the identity portion of the Jacobian with compensated summation:
    comp_sum_matrix!(s.jac_step,s.jac_error,d.jac_copy)
    indi=0; indj=0
    @inbounds for i=s.n-1:-1:1
        indi=(i-1)*7
        @inbounds for j=s.n:-1:i+1
            indj=(j-1)*7
            if ~s.pair[i,j]  # Check to see if kicks have not been applied
                #kepler_driftij_gamma!(s.m,s.x,s.v,s.xerror,s.verror,i,j,h2,d.jac_ij,d.dqdt_ij,false)
                kepler_driftij_gamma!(s,d,i,j,h2,false)
                # Pick out indices for bodies i & j:
                # Carry out multiplication on the i/j components of matrix:
                @inbounds for k2=1:sevn, k1=1:7
                    d.jac_tmp1[k1,k2] = s.jac_step[ indi+k1,k2]
                    d.jac_err1[k1,k2] = s.jac_error[indi+k1,k2]
                end
                @inbounds for k2=1:sevn, k1=1:7
                    d.jac_tmp1[7+k1,k2] = s.jac_step[ indj+k1,k2]
                    d.jac_err1[7+k1,k2] = s.jac_error[indj+k1,k2]
                end
                # Carry out multiplication on the i/j components of matrix:
                if T == BigFloat
                    d.jac_tmp2 .= *(d.jac_ij,d.jac_tmp1)
                else
                    start = time()
                    #BLAS.gemm!('N','N',uno,d.jac_ij,d.jac_tmp1,zilch,d.jac_tmp2)
                    mul!(d.jac_tmp2,d.jac_ij,d.jac_tmp1)
                    d.ctime[1] += time()-start
                end
                # Add back in the Jacobian with compensated summation:
                comp_sum_matrix!(d.jac_tmp1,d.jac_err1,d.jac_tmp2)
                # Copy back to the Jacobian:
                @inbounds for k2=1:sevn, k1=1:7
                    s.jac_step[ indi+k1,k2]=d.jac_tmp1[k1,k2]
                    s.jac_error[indi+k1,k2]=d.jac_err1[k1,k2]
                end
                @inbounds for k2=1:sevn, k1=1:7
                    s.jac_step[ indj+k1,k2]=d.jac_tmp1[7+k1,k2]
                    s.jac_error[indj+k1,k2]=d.jac_err1[7+k1,k2]
                end
            end
        end
    end
    drift_grad!(s,h2)
    #kickfast!(x,v,h2,m,n,jac_kick,dqdt_kick,pair)
    #kickfast!(s.x,s.v,s.xerror,s.verror,h/6,s.m,s.n,d.jac_kick,d.dqdt_kick,pair)
    kickfast!(s,d,h/6)
    # Multiply Jacobian from kick step:
    if T == BigFloat
        d.jac_copy .= *(d.jac_kick,s.jac_step)
    else
        start = time()
        #BLAS.gemm!('N','N',uno,d.jac_kick,s.jac_step,zilch,d.jac_copy)
        mul!(d.jac_copy,d.jac_kick,s.jac_step)
        d.ctime[1] += time()-start
    end
    # Add back in the identity portion of the Jacobian with compensated summation:
    comp_sum_matrix!(s.jac_step,s.jac_error,d.jac_copy)
    # Edit this routine to do compensated summation for Jacobian [x]
    #drift!(s.x,s.v,s.xerror,s.verror,h2,s.n,s.jac_step,s.jac_error)
    return
end

function ahl21!(s::State{T},d::dTime{T},h::T) where T<:AbstractFloat
    # [Currently this routine is not giving the correct dqdt values. -EA 8/12/2019]
    n = s.n
    zilch = zero(T); uno = one(T); half = convert(T,0.5); two = convert(T,2.0); h2 = half*h; sevn = 7*n
    zero_out!(d)
    fill!(s.dqdt,zilch)
    kickfast!(s,d,h/6)
    d.dqdt_kick ./= 6 # Since step is h/6
    # Since I removed identity from kickfast, need to add in dqdt:
    s.dqdt .+= d.dqdt_kick + *(d.jac_kick,s.dqdt)
    drift!(s,h2)
    # Compute time derivative of drift step:
    @inbounds for i=1:n, k=1:3
        s.dqdt[(i-1)*7+k] = half*s.v[k,i] + h2*s.dqdt[(i-1)*7+3+k]
    end
    @inbounds for i=1:n-1
        indi = (i-1)*7
        @inbounds for j=i+1:n
            indj = (j-1)*7
            if ~s.pair[i,j]  # Check to see if kicks have not been applied
                kepler_driftij_gamma!(s,d,i,j,h2,true)
                # Copy current time derivatives for multiplication purposes:
                @inbounds for k1=1:7
                    d.dqdt_tmp1[  k1] = s.dqdt[indi+k1]
                    d.dqdt_tmp1[7+k1] = s.dqdt[indj+k1]
                end
                # Add in partial derivatives with respect to time:
                # Need to multiply by 1/2 since we're taking 1/2 time step:
                d.dqdt_ij .*= half
                d.dqdt_ij .+= *(d.jac_ij,d.dqdt_tmp1) + d.dqdt_tmp1
                # Copy back time derivatives:
                @inbounds for k1=1:7
                    s.dqdt[indi+k1] = d.dqdt_ij[  k1]
                    s.dqdt[indj+k1] = d.dqdt_ij[7+k1]
                end
            end
        end
    end
    # Looks like we are missing phic! here: [ ]
    # Since I haven't added dqdt to phic yet, for now, set jac_phi equal to identity matrix
    # (since this is commented out within phisalpha):
    phic!(s,d,h)
    phisalpha!(s,d,h,two)
    # Add in time derivative with respect to prior parameters:
    # Copy result to dqdt:
    s.dqdt .+= d.dqdt_phi + *(d.jac_phi,s.dqdt)
    indi=0; indj=0
    @inbounds for i=n-1:-1:1
        indi=(i-1)*7
        @inbounds for j=n:-1:i+1
            if ~s.pair[i,j]  # Check to see if kicks have not been applied
                indj=(j-1)*7
                kepler_driftij_gamma!(s,d,i,j,h2,false)
                # Copy current time derivatives for multiplication purposes:
                @inbounds for k1=1:7
                    d.dqdt_tmp1[  k1] = s.dqdt[indi+k1]
                    d.dqdt_tmp1[7+k1] = s.dqdt[indj+k1]
                end
                # Add in partial derivatives with respect to time:
                # Need to multiply by 1/2 since we're taking 1/2 time step:
                d.dqdt_ij .*= half
                d.dqdt_ij .+= *(d.jac_ij,d.dqdt_tmp1) + d.dqdt_tmp1
                # Copy back time derivatives:
                @inbounds for k1=1:7
                    s.dqdt[indi+k1] = d.dqdt_ij[  k1]
                    s.dqdt[indj+k1] = d.dqdt_ij[7+k1]
                end
            end
        end
    end
    drift_grad!(s,h2)
    # Compute time derivative of drift step:
    @inbounds for i=1:n, k=1:3
        s.dqdt[(i-1)*7+k] += half*s.v[k,i] + h2*s.dqdt[(i-1)*7+3+k]
    end
    fill!(d.dqdt_kick,zilch)
    kickfast!(s,d,h/6)
    d.dqdt_kick ./= 6 # Since step is h/6
    # Copy result to dqdt:
    s.dqdt .+= d.dqdt_kick + *(d.jac_kick,s.dqdt)
    return
end

"""

AHL21 drift step. Drifts all particles with compensated summation. (with Jacobian)
"""
function drift_grad!(s::State{T},h::T) where {T <: Real}
    indi::Int64 = 0
    @inbounds for i=1:s.n
        indi = (i-1)*7
        @inbounds for j=1:NDIM
            s.x[j,i],s.xerror[j,i] = comp_sum(s.x[j,i],s.xerror[j,i],h*s.v[j,i])
        end
        # Now for Jacobian:
        @inbounds for k=1:7*s.n, j=1:NDIM
            s.jac_step[indi+j,k],s.jac_error[indi+j,k] = comp_sum(s.jac_step[indi+j,k],s.jac_error[indi+j,k],h*s.jac_step[indi+3+j,k])
        end
    end
    return
end

"""

AHL21 kick step. Computes "fast" kicks for pairs of bodies (in lieu of -drift+Kepler). Include Jacobian and compensated summation.
"""
function kickfast!(s::State{T},d::AbstractDerivatives{T},h::T) where {T <: Real}
    n::Int64 = s.n
    s.rij .= 0.0
    # Getting rid of identity since we will add that back in in calling routines:
    d.jac_kick .= 0.0
    @inbounds for i=1:n-1
        indi = (i-1)*7
        for j=i+1:n
            indj = (j-1)*7
            if s.pair[i,j]
                for k=1:3
                    s.rij[k] = s.x[k,i] - s.x[k,j]
                end
                r2inv::T = 1.0/dot_fast(s.rij)
                r3inv::T = r2inv*sqrt(r2inv)
                fac2::T  = h*GNEWT*r3inv
                for k=1:3
                    fac::T = fac2*s.rij[k]
                    # Apply impulses:
                    s.v[k,i],s.verror[k,i] = comp_sum(s.v[k,i],s.verror[k,i],-s.m[j]*fac)
                    s.v[k,j],s.verror[k,j] = comp_sum(s.v[k,j],s.verror[k,j], s.m[i]*fac)
                    # Compute time derivative:
                    d.dqdt_kick[indi+3+k] -= s.m[j]*fac/h
                    d.dqdt_kick[indj+3+k] += s.m[i]*fac/h
                    # Computing the derivative
                    # Mass derivative of acceleration vector (10/6/17 notes):
                    # Impulse of ith particle depends on mass of jth particle:
                    d.jac_kick[indi+3+k,indj+7] -= fac
                    # Impulse of jth particle depends on mass of ith particle:
                    d.jac_kick[indj+3+k,indi+7] += fac
                    # x derivative of acceleration vector:
                    fac *= 3.0*r2inv
                    # Dot product x_ij.\delta x_ij means we need to sum over components:
                    for p=1:3
                        d.jac_kick[indi+3+k,indi+p] += fac*s.m[j]*s.rij[p]
                        d.jac_kick[indi+3+k,indj+p] -= fac*s.m[j]*s.rij[p]
                        d.jac_kick[indj+3+k,indj+p] += fac*s.m[i]*s.rij[p]
                        d.jac_kick[indj+3+k,indi+p] -= fac*s.m[i]*s.rij[p]
                    end
                    # Final term has no dot product, so just diagonal:
                    d.jac_kick[indi+3+k,indi+k] -= fac2*s.m[j]
                    d.jac_kick[indi+3+k,indj+k] += fac2*s.m[j]
                    d.jac_kick[indj+3+k,indj+k] -= fac2*s.m[i]
                    d.jac_kick[indj+3+k,indi+k] += fac2*s.m[i]
                end
            end
        end
    end
    return
end

"""

Computes correction for pairs which are kicked, with Jacobian, dq/dt, and compensated summation.
"""
function phic!(s::State{T},d::AbstractDerivatives{T},h::T) where {T <: Real}
    s.a .= 0.0
    s.rij .= 0.0
    s.aij .= 0.0
    d.dadq .= 0.0  # There is no velocity dependence
    d.dotdadq .= 0.0 # There is no velocity dependence
    # Set jac_step to zeros:
    fill!(d.jac_phi,zero(T))
    fac::T = 0.0; fac1::T = 0.0; fac2::T = 0.0; fac3::T = 0.0; r1::T = 0.0; r2::T = 0.0; r3::T = 0.0
    coeff::T = h^3/36*GNEWT
    n::Int64 = s.n
    @inbounds for i=1:n-1
        indi = (i-1)*7
        for j=i+1:n
            if s.pair[i,j]
                indj = (j-1)*7
                for k=1:3
                    s.rij[k] = s.x[k,i] - s.x[k,j]
                end
                r2inv::T = inv(dot_fast(s.rij))
                r3inv::T = r2inv*sqrt(r2inv)
                for k=1:3
                    # Apply impulses:
                    fac = GNEWT*s.rij[k]*r3inv
                    facv = fac*2*h/3
                    s.v[k,i],s.verror[k,i] = comp_sum(s.v[k,i],s.verror[k,i],-s.m[j]*facv)
                    s.v[k,j],s.verror[k,j] = comp_sum(s.v[k,j],s.verror[k,j],s.m[i]*facv)
                    # Compute time derivative:
                    d.dqdt_phi[indi+3+k] -= 1/h*s.m[j]*facv
                    d.dqdt_phi[indj+3+k] += 1/h*s.m[i]*facv
                    s.a[k,i] -= s.m[j]*fac
                    s.a[k,j] += s.m[i]*fac
                    # Impulse of ith particle depends on mass of jth particle:
                    d.jac_phi[indi+3+k,indj+7] -= facv
                    # Impulse of jth particle depends on mass of ith particle:
                    d.jac_phi[indj+3+k,indi+7] += facv
                    # x derivative of acceleration vector:
                    facv *= 3.0*r2inv
                    # Dot product x_ij.\delta x_ij means we need to sum over components:
                    for p=1:3
                        d.jac_phi[indi+3+k,indi+p] += facv*s.m[j]*s.rij[p]
                        d.jac_phi[indi+3+k,indj+p] -= facv*s.m[j]*s.rij[p]
                        d.jac_phi[indj+3+k,indj+p] += facv*s.m[i]*s.rij[p]
                        d.jac_phi[indj+3+k,indi+p] -= facv*s.m[i]*s.rij[p]
                    end
                    # Final term has no dot product, so just diagonal:
                    facv = 2h/3*GNEWT*r3inv
                    d.jac_phi[indi+3+k,indi+k] -= facv*s.m[j]
                    d.jac_phi[indi+3+k,indj+k] += facv*s.m[j]
                    d.jac_phi[indj+3+k,indj+k] -= facv*s.m[i]
                    d.jac_phi[indj+3+k,indi+k] += facv*s.m[i]
                    # Mass derivative of acceleration vector (10/6/17 notes):
                    # Since there is no velocity dependence, this is fourth parameter.
                    # Acceleration of ith particle depends on mass of jth particle:
                    d.dadq[k,i,4,j] -= fac
                    d.dadq[k,j,4,i] += fac
                    # x derivative of acceleration vector:
                    fac *= 3.0*r2inv
                    # Dot product x_ij.\delta x_ij means we need to sum over components:
                    for p=1:3
                        d.dadq[k,i,p,i] += fac*s.m[j]*s.rij[p]
                        d.dadq[k,i,p,j] -= fac*s.m[j]*s.rij[p]
                        d.dadq[k,j,p,j] += fac*s.m[i]*s.rij[p]
                        d.dadq[k,j,p,i] -= fac*s.m[i]*s.rij[p]
                    end
                    # Final term has no dot product, so just diagonal:
                    fac = GNEWT*r3inv
                    d.dadq[k,i,k,i] -= fac*s.m[j]
                    d.dadq[k,i,k,j] += fac*s.m[j]
                    d.dadq[k,j,k,j] -= fac*s.m[i]
                    d.dadq[k,j,k,i] += fac*s.m[i]
                end
            end
        end
    end
    # Next, compute g_i acceleration vector.
    # Note that jac_step[(i-1)*7+k,(j-1)*7+p] is the derivative of the kth coordinate
    # of planet i with respect to the pth coordinate of planet j.
    indi::Int64 = 0; indj::Int64 = 0; indd::Int64 = 0
    @inbounds for i=1:n-1
        indi = (i-1)*7
        for j=i+1:n
            if s.pair[i,j] # correction for Kepler pairs
                indj = (j-1)*7
                for k=1:3
                    s.aij[k] = s.a[k,i] - s.a[k,j]
                    s.rij[k] = s.x[k,i] - s.x[k,j]
                end
                # Compute dot product of r_ij with \delta a_ij:
                fill!(d.dotdadq,0.0)
                @inbounds for di=1:n, p=1:4, k=1:3
                    d.dotdadq[p,di] += s.rij[k]*(d.dadq[k,i,p,di]-d.dadq[k,j,p,di])
                end
                r2 = dot_fast(s.rij)
                r1 = sqrt(r2)
                ardot = dot_fast(s.aij,s.rij)
                fac1 = coeff/(r2*r2*r1)
                fac2 = 3*ardot
                for k=1:3
                    fac = fac1*(s.rij[k]*fac2- r2*s.aij[k])
                    s.v[k,i],s.verror[k,i] = comp_sum(s.v[k,i],s.verror[k,i],s.m[j]*fac)
                    s.v[k,j],s.verror[k,j] = comp_sum(s.v[k,j],s.verror[k,j],-s.m[i]*fac)
                    # Compute time derivative of 4th order correction
                    d.dqdt_phi[indi+3+k] += 3/h*s.m[j]*fac
                    d.dqdt_phi[indj+3+k] -= 3/h*s.m[i]*fac
                    # Mass derivative (first part is easy):
                    d.jac_phi[indi+3+k,indj+7] += fac
                    d.jac_phi[indj+3+k,indi+7] -= fac
                    # Position derivatives:
                    fac *= 5.0/r2
                    for p=1:3
                        d.jac_phi[indi+3+k,indi+p] -= fac*s.m[j]*s.rij[p]
                        d.jac_phi[indi+3+k,indj+p] += fac*s.m[j]*s.rij[p]
                        d.jac_phi[indj+3+k,indj+p] -= fac*s.m[i]*s.rij[p]
                        d.jac_phi[indj+3+k,indi+p] += fac*s.m[i]*s.rij[p]
                    end
                    # Diagonal position terms:
                    fac = fac1*fac2
                    d.jac_phi[indi+3+k,indi+k] += fac*s.m[j]
                    d.jac_phi[indi+3+k,indj+k] -= fac*s.m[j]
                    d.jac_phi[indj+3+k,indj+k] += fac*s.m[i]
                    d.jac_phi[indj+3+k,indi+k] -= fac*s.m[i]
                    # Dot product \delta rij terms:
                    fac = -2*fac1*s.aij[k]
                    for p=1:3
                        fac3 = fac*s.rij[p] + fac1*3.0*s.rij[k]*s.aij[p]
                        d.jac_phi[indi+3+k,indi+p] += s.m[j]*fac3
                        d.jac_phi[indi+3+k,indj+p] -= s.m[j]*fac3
                        d.jac_phi[indj+3+k,indj+p] += s.m[i]*fac3
                        d.jac_phi[indj+3+k,indi+p] -= s.m[i]*fac3
                    end
                    # Diagonal acceleration terms:
                    fac = -fac1*r2
                    # Duoh.  For dadq, have to loop over all other parameters!
                    @inbounds for di=1:n
                        indd = (di-1)*7
                        for p=1:3
                            d.jac_phi[indi+3+k,indd+p] += fac*s.m[j]*(d.dadq[k,i,p,di]-d.dadq[k,j,p,di])
                            d.jac_phi[indj+3+k,indd+p] -= fac*s.m[i]*(d.dadq[k,i,p,di]-d.dadq[k,j,p,di])
                        end
                        # Don't forget mass-dependent term:
                        d.jac_phi[indi+3+k,indd+7] += fac*s.m[j]*(d.dadq[k,i,4,di]-d.dadq[k,j,4,di])
                        d.jac_phi[indj+3+k,indd+7] -= fac*s.m[i]*(d.dadq[k,i,4,di]-d.dadq[k,j,4,di])
                    end
                    # Now, for the final term:  (\delta a_ij . r_ij ) r_ij
                    fac = 3.0*fac1*s.rij[k]
                    @inbounds for di=1:n
                        indd = (di-1)*7
                        for p=1:3
                            d.jac_phi[indi+3+k,indd+p] += fac*s.m[j]*d.dotdadq[p,di]
                            d.jac_phi[indj+3+k,indd+p] -= fac*s.m[i]*d.dotdadq[p,di]
                        end
                        d.jac_phi[indi+3+k,indd+7] += fac*s.m[j]*d.dotdadq[4,di]
                        d.jac_phi[indj+3+k,indd+7] -= fac*s.m[i]*d.dotdadq[4,di]
                    end
                end
            end
        end
    end
    return
end

"""

Computes the 4th-order correction, with Jacobian, dq/dt, and compensated summation.
"""
function phisalpha!(s::State{T},d::AbstractDerivatives{T},h::T,alpha::T) where {T <: Real}
    s.a .= 0.0
    d.dadq .= 0.0  # There is no velocity dependence
    d.dotdadq .= 0.0  # There is no velocity dependence
    s.rij .= 0.0
    s.aij .= 0.0
    coeff::T = alpha*h^3/96*2*GNEWT
    fac::T = 0.0; fac1::T = 0.0; fac2::T = 0.0; fac3::T = 0.0; r1::T = 0.0; r2::T = 0.0; r3::T = 0.0
    n::Int64 = s.n
    @inbounds for i=1:n-1
        for j=i+1:n
            if ~s.pair[i,j] # correction for Kepler pairs
                for k=1:3
                    s.rij[k] = s.x[k,i] - s.x[k,j]
                end
                r2 = dot_fast(s.rij)
                r3 = r2*sqrt(r2)
                fac2 = GNEWT/r3
                for k=1:3
                    fac = fac2*s.rij[k]
                    s.a[k,i] -= s.m[j]*fac
                    s.a[k,j] += s.m[i]*fac
                    # Mass derivative of acceleration vector (10/6/17 notes):
                    # Since there is no velocity dependence, this is fourth parameter.
                    # Acceleration of ith particle depends on mass of jth particle:
                    d.dadq[k,i,4,j] -= fac
                    d.dadq[k,j,4,i] += fac
                    # x derivative of acceleration vector:
                    fac *= 3.0/r2
                    # Dot product x_ij.\delta x_ij means we need to sum over components:
                    for p=1:3
                        d.dadq[k,i,p,i] += fac*s.m[j]*s.rij[p]
                        d.dadq[k,i,p,j] -= fac*s.m[j]*s.rij[p]
                        d.dadq[k,j,p,j] += fac*s.m[i]*s.rij[p]
                        d.dadq[k,j,p,i] -= fac*s.m[i]*s.rij[p]
                    end
                    # Final term has no dot product, so just diagonal:
                    d.dadq[k,i,k,i] -= fac2*s.m[j]
                    d.dadq[k,i,k,j] += fac2*s.m[j]
                    d.dadq[k,j,k,j] -= fac2*s.m[i]
                    d.dadq[k,j,k,i] += fac2*s.m[i]
                end
            end
        end
    end
    # Next, compute \tilde g_i acceleration vector (this is rewritten
    # slightly to avoid reference to \tilde a_i):
    # Note that jac_step[(i-1)*7+k,(j-1)*7+p] is the derivative of the kth coordinate
    # of planet i with respect to the pth coordinate of planet j.
    indi = 0; indj=0; indd = 0
    @inbounds for i=1:n-1
        indi = (i-1)*7
        for j=i+1:n
            if ~s.pair[i,j] # correction for Kepler pairs
                indj = (j-1)*7
                for k=1:3
                    s.aij[k] = s.a[k,i] - s.a[k,j]
                    s.rij[k] = s.x[k,i] - s.x[k,j]
                end
                # Compute dot product of r_ij with \delta a_ij:
                fill!(d.dotdadq,0.0)
                @inbounds for di=1:n, p=1:4, k=1:3
                    d.dotdadq[p,di] += s.rij[k]*(d.dadq[k,i,p,di]-d.dadq[k,j,p,di])
                end
                r2 = dot_fast(s.rij)
                r1 = sqrt(r2)
                ardot = dot_fast(s.aij, s.rij)
                fac1 = coeff/(r2*r2*r1)
                fac2 = (2*GNEWT*(s.m[i]+s.m[j])/r1 + 3*ardot)
                for k=1:3
                    fac = fac1*(s.rij[k]*fac2-r2*s.aij[k])
                    s.v[k,i],s.verror[k,i] = comp_sum(s.v[k,i],s.verror[k,i], s.m[j]*fac)
                    s.v[k,j],s.verror[k,j] = comp_sum(s.v[k,j],s.verror[k,j],-s.m[i]*fac)
                    # Compute time derivative:
                    d.dqdt_phi[indi+3+k] += 3/h*s.m[j]*fac
                    d.dqdt_phi[indj+3+k] -= 3/h*s.m[i]*fac
                    # Mass derivative (first part is easy):
                    d.jac_phi[indi+3+k,indj+7] += fac
                    d.jac_phi[indj+3+k,indi+7] -= fac
                    # Position derivatives:
                    fac *= 5.0/r2
                    for p=1:3
                        d.jac_phi[indi+3+k,indi+p] -= fac*s.m[j]*s.rij[p]
                        d.jac_phi[indi+3+k,indj+p] += fac*s.m[j]*s.rij[p]
                        d.jac_phi[indj+3+k,indj+p] -= fac*s.m[i]*s.rij[p]
                        d.jac_phi[indj+3+k,indi+p] += fac*s.m[i]*s.rij[p]
                    end
                    # Second mass derivative:
                    fac = 2*GNEWT*fac1*s.rij[k]/r1
                    d.jac_phi[indi+3+k,indi+7] += fac*s.m[j]
                    d.jac_phi[indi+3+k,indj+7] += fac*s.m[j]
                    d.jac_phi[indj+3+k,indj+7] -= fac*s.m[i]
                    d.jac_phi[indj+3+k,indi+7] -= fac*s.m[i]
                    #  (There's also a mass term in dadq [x]. See below.)
                    # Diagonal position terms:
                    fac = fac1*fac2
                    d.jac_phi[indi+3+k,indi+k] += fac*s.m[j]
                    d.jac_phi[indi+3+k,indj+k] -= fac*s.m[j]
                    d.jac_phi[indj+3+k,indj+k] += fac*s.m[i]
                    d.jac_phi[indj+3+k,indi+k] -= fac*s.m[i]
                    # Dot product \delta rij terms:
                    fac = -2*fac1*(s.rij[k]*GNEWT*(s.m[i]+s.m[j])/(r2*r1)+s.aij[k])
                    for p=1:3
                        fac3 = fac*s.rij[p] + fac1*3.0*s.rij[k]*s.aij[p]
                        d.jac_phi[indi+3+k,indi+p] += s.m[j]*fac3
                        d.jac_phi[indi+3+k,indj+p] -= s.m[j]*fac3
                        d.jac_phi[indj+3+k,indj+p] += s.m[i]*fac3
                        d.jac_phi[indj+3+k,indi+p] -= s.m[i]*fac3
                    end
                    # Diagonal acceleration terms:
                    fac = -fac1*r2
                    # Duoh.  For dadq, have to loop over all other parameters!
                    @inbounds for di=1:n
                        indd = (di-1)*7
                        for p=1:3
                            d.jac_phi[indi+3+k,indd+p] += fac*s.m[j]*(d.dadq[k,i,p,di]-d.dadq[k,j,p,di])
                        end
                        for p=1:3
                            d.jac_phi[indj+3+k,indd+p] -= fac*s.m[i]*(d.dadq[k,i,p,di]-d.dadq[k,j,p,di])
                        end
                        # Don't forget mass-dependent term:
                        d.jac_phi[indi+3+k,indd+7] += fac*s.m[j]*(d.dadq[k,i,4,di]-d.dadq[k,j,4,di])
                        d.jac_phi[indj+3+k,indd+7] -= fac*s.m[i]*(d.dadq[k,i,4,di]-d.dadq[k,j,4,di])
                    end
                    # Now, for the final term:  (\delta a_ij . r_ij ) r_ij
                    fac = 3.0*fac1*s.rij[k]
                    @inbounds for di=1:n
                        indd = (di-1)*7
                        for p=1:3
                            d.jac_phi[indi+3+k,indd+p] += fac*s.m[j]*d.dotdadq[p,di]
                        end
                        for p=1:3
                            d.jac_phi[indj+3+k,indd+p] -= fac*s.m[i]*d.dotdadq[p,di]
                        end
                        d.jac_phi[indi+3+k,indd+7] += fac*s.m[j]*d.dotdadq[4,di]
                        d.jac_phi[indj+3+k,indd+7] -= fac*s.m[i]*d.dotdadq[4,di]
                    end
                end
            end
        end
    end
    return
end

"""

Carries out a Kepler step and reverse drift for bodies i & j, and computes Jacobian. Uses new version of the code with gamma in favor of s.
"""
function kepler_driftij_gamma!(s::State{T},d::AbstractDerivatives{T},i::Int64,j::Int64,h::T,drift_first::Bool) where {T <: Real}
    # Initial state:
    @inbounds for k=1:NDIM
        s.x0[k] = s.x[k,i] - s.x[k,j] # x0 = positions of body i relative to j
        s.v0[k] = s.v[k,i] - s.v[k,j] # v0 = velocities of body i relative to j
    end
    gm = GNEWT*(s.m[i]+s.m[j])
    if gm == 0; return; end
    # jac_ij should be the Jacobian for going from (x_{0,i},v_{0,i},m_i) &  (x_{0,j},v_{0,j},m_j)
    # to  (x_i,v_i,m_i) &  (x_j,v_j,m_j), a 14x14 matrix for the 3-dimensional case.:
    fill!(d.jac_ij,zero(T))
    s.delxv .= 0.0
    d.jac_kepler .= 0.0
    d.jac_mass .= 0.0
    params::NTuple{22,T} = jac_delxv_gamma!(s,gm,h,drift_first)
    compute_jacobian_gamma!(params...,s.x0,s.v0,d.jac_kepler,d.jac_mass,drift_first,false)
    mijinv::T = one(T)/(s.m[i] + s.m[j])
    mi::T = s.m[i]*mijinv # Normalize the masses
    mj::T = s.m[j]*mijinv
    @inbounds for k=1:3
        # Add kepler-drift differences, weighted by masses, to start of step:
        s.x[k,i],s.xerror[k,i] = comp_sum(s.x[k,i],s.xerror[k,i], mj*s.delxv[k])
        s.x[k,j],s.xerror[k,j] = comp_sum(s.x[k,j],s.xerror[k,j],-mi*s.delxv[k])
    end
    @inbounds for k=1:3
        s.v[k,i],s.verror[k,i] = comp_sum(s.v[k,i],s.verror[k,i], mj*s.delxv[3+k])
        s.v[k,j],s.verror[k,j] = comp_sum(s.v[k,j],s.verror[k,j],-mi*s.delxv[3+k])
    end
    # Compute Jacobian:
    @inbounds for l=1:6, k=1:6
        # Compute derivatives of x_i,v_i with respect to initial conditions:
        d.jac_ij[  k,  l] += mj*d.jac_kepler[k,l]
        d.jac_ij[  k,7+l] -= mj*d.jac_kepler[k,l]
        # Compute derivatives of x_j,v_j with respect to initial conditions:
        d.jac_ij[7+k,  l] -= mi*d.jac_kepler[k,l]
        d.jac_ij[7+k,7+l] += mi*d.jac_kepler[k,l]
    end
    @inbounds for k=1:6
        # Compute derivatives of x_i,v_i with respect to the masses:
        d.jac_ij[   k, 7] = d.jac_mass[k]*s.m[j]
        d.jac_ij[   k,14] =  mi*s.delxv[k]*mijinv + GNEWT*mj*d.jac_kepler[  k,7]
        # Compute derivatives of x_j,v_j with respect to the masses:
        d.jac_ij[ 7+k, 7] = -mj*s.delxv[k]*mijinv - GNEWT*mi*d.jac_kepler[  k,7]
        d.jac_ij[ 7+k,14] = -d.jac_mass[k]*s.m[i]
    end
    # The following lines are meant to compute dq/dt for kepler_driftij,
    # but they currently contain an error (test doesn't pass in test_ahl21.jl). [ ]
    @inbounds for k=1:6
        # Position/velocity derivative, body i:
        d.dqdt_ij[  k] =  mj*d.jac_kepler[k,8]
        # Position/velocity derivative, body j:
        d.dqdt_ij[7+k] = -mi*d.jac_kepler[k,8]
    end
    return
end

"""

Computes Jacobian of delx and delv with respect to x0, v0, k, and h.
"""
function jac_delxv_gamma!(s::State{T},k::T,h::T,drift_first::Bool;debug::Bool=false,grad::Bool=true) where {T <: Real}
    # Compute r0:
    r0 = zero(T)
    s.rtmp[1] = s.x0[1]-h*s.v0[1]
    s.rtmp[2] = s.x0[2]-h*s.v0[2]
    s.rtmp[3] = s.x0[3]-h*s.v0[3]
    drift_first ?  r0 = norm(s.rtmp) : r0 = norm(s.x0)
    # And its inverse:
    r0inv::T = inv(r0)
    # Compute beta_0:
    beta0::T = 2k*r0inv-dot_fast(s.v0,s.v0)
    beta0inv::T = inv(beta0)
    signb::T = sign(beta0)
    sqb::T = sqrt(signb*beta0)
    zeta::T = k-r0*beta0
    gamma_guess = zero(T)
    # Compute \eta_0 = x_0 . v_0:
    eta = zero(T)
    #drift_first ?  eta = dot(s.x0 .- h .* s.v0, s.v0) : eta = dot(s.x0,s.v0)
    if drift_first
        eta = dot_fast(s.rtmp,s.v0)
    else
        eta = dot_fast(s.x0,s.v0)
    end
    if zeta != zero(T)
        # Make sure we have a cubic in gamma (and don't divide by zero):
        zinv = 6/zeta
        gamma_guess = cubic1(0.5*eta*sqb*zinv,r0*signb*beta0*zinv,-h*signb*beta0*sqb*zinv)
    else
        # Check that we have a quadratic in gamma (and don't divide by zero):
        if eta != zero(T)
            reta = r0/eta
            disc = reta^2+2h/eta
            disc > zero(T) ?  gamma_guess = sqb*(-reta+sqrt(disc)) : gamma_guess = h*r0inv*sqb
        else
            gamma_guess = h*r0inv*sqb
        end
    end
    gamma  = copy(gamma_guess)
    # Make sure prior two steps differ:
    gamma1::T = 2*copy(gamma)
    gamma2::T = 3*copy(gamma)
    iter = 0
    ITMAX = 20
    # Compute coefficients: (8/28/19 notes)
    c1 = k; c2 = -2zeta; c3 = 2*eta*signb*sqb; c4 = -sqb*h*beta0; c5 = 2eta*signb*sqb
    # Solve for gamma:
    while iter < ITMAX
        gamma2 = gamma1
        gamma1 = gamma
        xx = 0.5*gamma
        if beta0 > 0
            sx,cx = sincos(xx);
        else
            sx = sinh(xx); cx = exp(-xx)+sx
        end
        gamma -= (k*gamma+c2*sx*cx+c3*sx^2+c4)/(2signb*zeta*sx^2+c5*sx*cx+r0*beta0)
        iter +=1
        if gamma == gamma2 || gamma == gamma1
            break
        end
    end
    # Set up a single output array for delx and delv:
    if debug
        s.delxv = zeros(T,12)
    else
        s.delxv .= zero(T)
    end
    # Since we updated gamma, need to recompute:
    xx = 0.5*gamma
    if beta0 > 0
        sx,cx = sincos(xx)
    else
        sx = sinh(xx); cx = exp(-xx)+sx
    end
    # Now, compute final values.  Compute Wisdom/Hernandez G_i^\beta(s) functions:
    g1bs = 2sx*cx/sqb
    g2bs = 2signb*sx^2*beta0inv
    g0bs = one(T)-beta0*g2bs
    g3bs = G3(gamma,beta0,sqb)
    h1 = zero(T); h2 = zero(T)
    # Compute r from equation (35):
    r = r0*g0bs+eta*g1bs+k*g2bs
    rinv = inv(r)
    dfdt = -k*g1bs*rinv*r0inv # [x]
    if drift_first
        # Drift backwards before Kepler step: (1/22/2018)
        fm1 = -k*r0inv*g2bs # [x]
        # This is given in 2/7/2018 notes: g-h*f
        gmh = k*r0inv*(h*g2bs-r0*g3bs)
    else
        # Drift backwards after Kepler step: (1/24/2018)
        # The following line is f-1-h fdot:
        h1= H1(gamma,beta0); h2= H2(gamma,beta0,sqb)
        fm1 =  k*rinv*(g2bs-k*r0inv*h1)
        # This is g-h*dgdt
        gmh = k*rinv*(r0*h2+eta*h1)
    end
    # Compute velocity component functions:
    if drift_first
        # This is gdot - h fdot - 1:
        dgdtm1 = k*r0inv*rinv*(h*g1bs-r0*g2bs)
    else
        # This is gdot - 1:
        dgdtm1 = -k*rinv*g2bs # [x]
    end
    @inbounds for j=1:3
        # Compute difference vectors (finish - start) of step:
        s.delxv[  j] = fm1*s.x0[j]+gmh*s.v0[j]        # position x_ij(t+h)-x_ij(t) - h*v_ij(t) or -h*v_ij(t+h)
    end
    @inbounds for j=1:3
        s.delxv[3+j] = dfdt*s.x0[j]+dgdtm1*s.v0[j]    # velocity v_ij(t+h)-v_ij(t)
    end
    if debug
        s.delxv[7] = gamma
        s.delxv[8] = r
        s.delxv[9] = fm1
        s.delxv[10] = dfdt
        s.delxv[11] = gmh
        s.delxv[12] = dgdtm1
    end
    if grad
        return gamma,g0bs,g1bs,g2bs,g3bs,h1,h2,dfdt,fm1,gmh,dgdtm1,r0,r,r0inv,rinv,k,h,beta0,beta0inv,eta,sqb,zeta
    end
end

"""

Computes the gradient analytically.
"""
function compute_jacobian_gamma!(gamma::T,g0::T,g1::T,g2::T,g3::T,h1::T,h2::T,dfdt::T,fm1::T,gmh::T,dgdtm1::T,
    r0::T,r::T,r0inv::T,rinv::T,k::T,h::T,beta::T,betainv::T,eta::T,sqb::T,zeta::T,x0::Array{T,1},v0::Array{T,1},
    delxv_jac::Array{T,2},jac_mass::Array{T,1},drift_first::Bool,debug::Bool) where {T <: Real}
    # Computes Jacobian:
    r0inv2 = r0inv^2
    r0inv3 = r0inv2*r0inv
    rinv2 = rinv^2
    rinv3 = rinv2*rinv
    hsq = h^2
    ksq = k^2
    if drift_first
        # First, the diagonal terms:
        # Now the off-diagonal terms:
        d   = (h + eta*g2 + 2*k*g3)*betainv
        c1 = d-r0*g3
        c2 = eta*g0+g1*zeta
        c3  = d*k+g1*r0^2
        c4 = eta*g1+2*g0*r0
        c13 = g1*h-g2*r0
        c9 = 2*g2*h-3*g3*r0
        c10 = k*r0inv2^2*(-g2*r0*h+k*c9*betainv-c3*c13*rinv)
        c24 = r0inv3*(r0*(2*k*r0inv-beta)*betainv-g1*c3*rinv/g2)
        h6 = H6(gamma,beta)
        # Derivatives of \delta x with respect to x0, v0, k & h:
        dfm1dxx = fm1*c24
        dfm1dxv = -fm1*(g1*rinv+h*c24)
        dfm1dvx = dfm1dxv
        dfm1dvv = fm1*rinv*(-r0*g2 + k*h6*betainv/g2 + h*(2*g1+h*r*c24))
        dfm1dh  = fm1*(g1*rinv*(1/g2+2*k*r0inv-beta)-eta*c24)
        dfm1dk  = fm1*(1/k+g1*c1*rinv*r0inv/g2-2*betainv*r0inv)
        h4 = -H1(gamma,beta)*beta
        h5 = H5(gamma,beta,sqb)
        dfm1dk2  = (r0*h4+k*h6)
        dgmhdxx = c10
        dgmhdxv =  -g2*k*c13*rinv*r0inv-h*c10
        dgmhdvx =  dgmhdxv
        h3 = H3(gamma,beta,sqb)
        h8 = -2h3+3h5
        dgmhdvv =  2*g2*h*k*c13*rinv*r0inv+hsq*c10+
        k*betainv*rinv*r0inv*(r0^2*h8-beta*h*r0*g2^2 + (h*k+eta*r0)*h6)
        dgmhdh  =  g2*k*r0inv+k*c13*rinv*r0inv+g2*k*(2*k*r0inv-beta)*c13*rinv*r0inv-eta*c10
        dgmhdk  =  r0inv*(k*c1*c13*rinv*r0inv+g2*h-g3*r0-k*c9*betainv*r0inv)
        dgmhdk2 = (h6*g3*ksq+eta*r0*(h6+g2*h4)+r0^2*g0*h5+k*eta*g2*h6+(g1*h6+g3*h4)*k*r0)
        @inbounds for j=1:3
            # First, compute the \delta x-derivatives:
            delxv_jac[  j,  j] = fm1
            delxv_jac[  j,3+j] = gmh
            @inbounds for i=1:3
                delxv_jac[j,  i] += (dfm1dxx*x0[i]+dfm1dxv*v0[i])*x0[j] + (dgmhdxx*x0[i]+dgmhdxv*v0[i])*v0[j]
                delxv_jac[j,3+i] += (dfm1dvx*x0[i]+dfm1dvv*v0[i])*x0[j] + (dgmhdvx*x0[i]+dgmhdvv*v0[i])*v0[j]
            end
            delxv_jac[  j,  7] = dfm1dk*x0[j] + dgmhdk*v0[j]
            delxv_jac[  j,  8] = dfm1dh*x0[j] + dgmhdh*v0[j]
            # Compute the mass jacobian separately since otherwise cancellations happen in kepler_driftij_gamma:
            jac_mass[  j] = (GNEWT*r0inv)^2*betainv*rinv*(dfm1dk2*x0[j]-dgmhdk2*v0[j])
        end
        # Derivatives of \delta v with respect to x0, v0, k & h:
        c5 = (r0-k*g2)*rinv/g1
        c6 = (r0*g0-k*g2)*betainv
        c7 = g2*(1/g1+c2*rinv)
        c8 = (k*c6+r*r0+c3*c5)*r0inv3
        c12 = g0*h-g1*r0
        c17 = r0-r-g2*k
        c18 = eta*g1+2*g2*k
        c20 = k*(g2*k+r)-g0*r0*zeta
        c21 = (g2*k-r0)*(beta*c3-k*g1*r)*betainv*rinv2*r0inv3/g1+eta*g1*rinv*r0inv2-2r0inv2
        c22 = rinv*(-g1-g0*g2/g1+g2*c2*rinv)
        c25 = k*rinv*r0inv2*(-g2+k*(c13-g2*r0)*betainv*r0inv2-c13*r0inv-c12*c3*rinv*r0inv2+
                              c13*c2*c3*rinv2*r0inv2-c13*(k*(g2*k+r)-g0*r0*zeta)*betainv*rinv*r0inv2)
        c26 = k*rinv2*r0inv*(-g2*c12-g1*c13+g2*c13*c2*rinv)
        ddfdtdxx = dfdt*c21
        ddfdtdxv = dfdt*(c22-h*c21)
        ddfdtdvx = ddfdtdxv
        c34 = (-beta*eta^2*g2^2-eta*k*h8-h6*ksq-2beta*eta*r0*g1*g2+(g2^2-3*g1*g3)*beta*k*r0
               - beta*g1^2*r0^2)*betainv*rinv2+(eta*g2^2)*rinv/g1 + (k*h8)*betainv*rinv/g1
        ddfdtdvv = dfdt*(c34 - 2*h*c22 +hsq*c21)
        ddfdtdk  = dfdt*(1/k-betainv*r0inv-c17*betainv*rinv*r0inv-c1*(g1*c2-g0*r)*rinv2*r0inv/g1)
        ddfdtdk2 = -(g2*k-r0)*(beta*r0*(g3-g1*g2)-beta*eta*g2^2+k*h3)*betainv*rinv2*r0inv
        ddfdtdh  = dfdt*(g0*rinv/g1-c2*rinv2-(2*k*r0inv-beta)*c22-eta*c21)
        dgdtmhdfdtm1dxx = c25
        dgdtmhdfdtm1dxv = c26-h*c25
        dgdtmhdfdtm1dvx = c26-h*c25
        h2 = H2(gamma,beta,sqb)
        c33 = d*k*rinv3*r0inv*k*(h*g2- r0*g3)+k*(-eta*k*g1*g2^2-g1*g2*g3*ksq-r0*eta*beta*g1*g2^2-r0*k*g1*h2 - beta*g2^2*g0*r0^2)*betainv*rinv2*r0inv
        dgdtmhdfdtm1dvv = c33-2*h*c26+hsq*c25
        dgdtmhdfdtm1dk = rinv*r0inv*(-k*(c13-g2*r0)*betainv*r0inv+c13-k*c13*c17*betainv*rinv*r0inv+k*c1*c12*rinv*r0inv-k*c1*c2*c13*rinv2*r0inv)
        dgdtmhdfdtm1dk2 = k*betainv*rinv2*r0inv*(-beta*eta^2*g2^4+eta*g2*(g1*g2^2+g1^2*g3-5*g2*g3)*k+g2*g3*h3*ksq+
                                                  2eta*r0*beta*g2^2*(g3-g1*g2)+(4g3-g0*g3-g1*g2)*(g3-g1*g2)*r0*k+beta*(2g1*g3*g2-g1^2*g2^2-g3^2)*r0^2)
        dgdtmhdfdtm1dh = g1*k*rinv*r0inv+k*c12*rinv2*r0inv-k*c2*c13*rinv3*r0inv-(2*k*r0inv-beta)*c26-eta*c25
        @inbounds for j=1:3
            # Next, compute the \delta v-derivatives:
            delxv_jac[3+j,  j] = dfdt
            delxv_jac[3+j,3+j] = dgdtm1
            @inbounds for i=1:3
                delxv_jac[3+j,  i] += (ddfdtdxx*x0[i]+ddfdtdxv*v0[i])*x0[j] + (dgdtmhdfdtm1dxx*x0[i]+dgdtmhdfdtm1dxv*v0[i])*v0[j]
                delxv_jac[3+j,3+i] += (ddfdtdvx*x0[i]+ddfdtdvv*v0[i])*x0[j] + (dgdtmhdfdtm1dvx*x0[i]+dgdtmhdfdtm1dvv*v0[i])*v0[j]
            end
            delxv_jac[3+j,  7] = ddfdtdk*x0[j] + dgdtmhdfdtm1dk*v0[j]
            delxv_jac[3+j,  8] = ddfdtdh*x0[j] + dgdtmhdfdtm1dh*v0[j]
            # Compute the mass jacobian separately since otherwise cancellations happen in kepler_driftij_gamma:
            jac_mass[3+j] = GNEWT^2*r0inv*rinv*(ddfdtdk2*x0[j]+dgdtmhdfdtm1dk2*v0[j])
        end
        if debug
            # Now include derivatives of gamma, r, fm1, dfdt, gmh, and dgdtmhdfdtm1:
            @inbounds for i=1:3
                delxv_jac[ 7,i] = -sqb*rinv*((g2-h*c3*r0inv3)*v0[i]+c3*x0[i]*r0inv3); delxv_jac[7,3+i] = sqb*rinv*((-d+2*g2*h-hsq*c3*r0inv3)*v0[i]+(-g2+h*c3*r0inv3)*x0[i])
                delxv_jac[ 8,i] = (c20*betainv-c2*c3*rinv)*r0inv3*x0[i]+((eta*g2+g1*r0)*rinv+h*r0inv3*(c2*c3*rinv-c20*betainv))*v0[i]
                drdv0x0 = (beta*g1*g2+((eta*g2+k*g3)*eta*g0*c3)*rinv*r0inv3 + (g1*g0*(2k*eta^2*g2+3eta*ksq*g3))*betainv*rinv*r0inv2-
                           k*betainv*r0inv3*(eta*g1*(eta*g2+k*g3)+g3*g0*r0^2*beta+2h*g2*k)+(g1*zeta)*rinv*((h*c3)*r0inv3 - g2) -
                           (eta*(beta*g2*g0*r0+k*g1^2)*(eta*g1+k*g2))*betainv*rinv*r0inv2)
                delxv_jac[8,3+i] = drdv0x0*x0[i] - (k*betainv*rinv*(eta*(beta*g2*g3-h8) - h6*k + (g2^2 - 2*g1*g3)*beta*r0) + h*drdv0x0)*v0[i]
                delxv_jac[ 9,i] = dfm1dxx*x0[i]+dfm1dxv*v0[i]; delxv_jac[ 9,3+i]=dfm1dvx*x0[i]+dfm1dvv*v0[i]
                delxv_jac[10,i] = ddfdtdxx*x0[i]+ddfdtdxv*v0[i]; delxv_jac[10,3+i]=ddfdtdvx*x0[i]+ddfdtdvv*v0[i]
                delxv_jac[11,i] = dgmhdxx*x0[i]+dgmhdxv*v0[i]; delxv_jac[11,3+i]=dgmhdvx*x0[i]+dgmhdvv*v0[i]
                delxv_jac[12,i] = dgdtmhdfdtm1dxx*x0[i]+dgdtmhdfdtm1dxv*v0[i]; delxv_jac[12,3+i]=dgdtmhdfdtm1dvx*x0[i]+dgdtmhdfdtm1dvv*v0[i]
            end
            delxv_jac[ 7,7] = sqb*c1*r0inv*rinv; delxv_jac[7,8] = sqb*rinv*(1+eta*c3*r0inv3+g2*(2k*r0inv-beta))
            delxv_jac[ 8,7] = betainv*r0inv*rinv*(-g2*r0^2*beta-eta*g1*g2*(k+beta*r0)+eta*g0*g3*(2*k+zeta)-
                                                  g2^2*(beta*eta^2+2*k*zeta)+g1*g3*zeta*(3*k-beta*r0)); delxv_jac[8,8] = c2*rinv
            delxv_jac[ 8,8] = ((r0*g1+eta*g2)*rinv)*(beta-2*k*r0inv)+c2*rinv+eta*r0inv3*(c2*c3*rinv-c20*betainv)
            delxv_jac[ 9,7] = dfm1dk; delxv_jac[ 9,8] = dfm1dh
            delxv_jac[10,7] = ddfdtdk; delxv_jac[10,8] = ddfdtdh
            delxv_jac[11,7] = dgmhdk; delxv_jac[11,8] = dgmhdh
            delxv_jac[12,7] = dgdtmhdfdtm1dk; delxv_jac[12,8] = dgdtmhdfdtm1dh
        end
    else
        # Now compute the Kepler-Drift Jacobian terms:
        # First, the diagonal terms:
        # Now the off-diagonal terms:
        d   = (h + eta*g2 + 2*k*g3)*betainv
        c1 = d-r0*g3
        c2 = eta*g0+g1*zeta
        c3  = d*k+g1*r0^2
        c4 = eta*g1+2*g0*r0
        c6 = (r0*g0-k*g2)*betainv
        c9  = g2*r-h1*k
        c14 = r0*g2-k*h1
        c15 = eta*h1+h2*r0
        c16 = eta*h2+g1*gamma*r0/sqb
        c17 = r0-r-g2*k
        c18 = eta*g1+2*g2*k
        c19 = 4*eta*h1+3*h2*r0
        c23 = h2*k-r0*g1
        h6 = H6(gamma,beta)
        h3 = H3(gamma,beta,sqb)
        h5 = H5(gamma,beta,sqb)
        h8 = -2h3+3h5
        # Derivatives of \delta x with respect to x0, v0, k & h:
        dfm1dxx = k*rinv3*betainv*r0inv^4*(k*h1*r^2*r0*(beta-2*k*r0inv)+beta*c3*(r*c23+c14*c2)+c14*r*(k*(r-g2*k)+g0*r0*zeta))
        dfm1dxv = k*rinv2*r0inv*(k*(g2*h2+g1*h1)-2g1*g2*r0+g2*c14*c2*rinv)
        dfm1dvx = dfm1dxv
        dfm1dvv = k*r0inv*rinv2*betainv*(2eta*k*(g2*g3-g1*h1)+(3g3*h2-4h1*g2)*ksq +
                                          beta*g2*r0*(3h1*k-g2*r0)+c14*rinv*(-beta*g2^2*eta^2+eta*k*(2g0*g3-h2)-
                                                                             h6*ksq+(-2eta*g1*g2+k*(h1-2g1*g3))*beta*r0-beta*g1^2*r0^2))
        dfm1dh  = (g1*k-h2*ksq*r0inv-k*c14*c2*rinv*r0inv)*rinv2
        dfm1dk  = rinv*r0inv*(4*h1*ksq*betainv*r0inv-k*h1-2*g2*k*betainv+c14-k*c14*c17*betainv*rinv*r0inv+
                              k*(g1*r0-k*h2)*c1*rinv*r0inv-k*c14*c1*c2*rinv2*r0inv)
        # New expression for d(f-1-h \dot f)/dk with cancellations of higher order terms in gamma is:
        dfm1dk2  = betainv*r0inv*rinv2*(r*(2eta*k*(g1*h1-g3*g2)+(4g2*h1-3g3*h2)*ksq-eta*r0*beta*g1*h1 + (g3*h2-4g2*h1)*beta*k*r0 + g2*h1*beta^2*r0^2) -
                                         # In the following line I need to replace 3g0*g3-g1*g2 by -H8:
                                         c14*(-eta^2*beta*g2^2 - k*eta*h8 - ksq*h6 - eta*r0*beta*(g1*g2 + g0*g3) + 2*(h1 - g1*g3)*beta*k*r0 - (g2 - beta*g1*g3)*beta*r0^2))
        dgmhdxx = k*rinv*r0inv*(h2+k*c19*betainv*r0inv2-c16*c3*rinv*r0inv2+c2*c3*c15*(rinv*r0inv)^2-c15*(k*(g2*k+r)-g0*r0*zeta)*betainv*rinv*r0inv2)
        dgmhdxv = k*rinv2*(h1*r-g2*c16-g1*c15+g2*c2*c15*rinv)
        dgmhdvx = dgmhdxv
        dgmhdvv = k*betainv*rinv2*(2*eta^2*(g1*h1-g2*g3)+eta*k*(4g2*h1-3h2*g3)+r0*eta*(4g0*h1-2g1*g3)+
                                    # In the following lines I need to replace g1*g2-3g0*g3-g1*g2 by H8:
                                    3r0*k*((g1+beta*g3)*h1-g3*g2)+(g0*h8-beta*g1*(g2^2+g1*g3))*r0^2 -
                                    c15*rinv*(beta*g2^2*eta^2+eta*k*h8+h6*ksq+(2eta*g1*g2-k*(g2^2-3g1*g3))*beta*r0+beta*g1^2*r0^2))
        dgmhdk  = rinv*(k*c1*c16*rinv*r0inv+c15-k*c15*c17*betainv*rinv*r0inv-k*c19*betainv*r0inv-k*c1*c2*c15*rinv2*r0inv)
        h7 = beta*g1*g2^2-g0*h8
        dgmhdk2 =  betainv*rinv2*(r*(2eta^2*(g3*g2-g1*h1) + eta*k*(3g3*h2 - 4g2*h1) +
                                      r0*eta*(beta*g3*(g1*g2 + g0*g3) - 2g0*h6) + (-h6*(g1 + beta*g3) + g2*(2g3 - h2))*r0*k +
                                      (h7 - beta^2*g1*g3^2)*r0^2)- c15*(-beta*eta^2*g2^2 + eta*k*(-h2 + 2g0*g3) - h6*ksq -
                                                                        r0*eta*beta*(h2 + 2g0*g3) + 2beta*(2*h1 - g2^2)*r0*k + beta*(beta*g1*g3 - g2)*r0^2))
        dgmhdh  = k*rinv3*(r*c16-c2*c15)
        @inbounds for j=1:3
            # First, compute the \delta x-derivatives:
            delxv_jac[  j,  j] = fm1
            delxv_jac[  j,3+j] = gmh
            @inbounds for i=1:3
                delxv_jac[j,  i] += (dfm1dxx*x0[i]+dfm1dxv*v0[i])*x0[j] + (dgmhdxx*x0[i]+dgmhdxv*v0[i])*v0[j]
                delxv_jac[j,3+i] += (dfm1dvx*x0[i]+dfm1dvv*v0[i])*x0[j] + (dgmhdvx*x0[i]+dgmhdvv*v0[i])*v0[j]
            end
            delxv_jac[  j,  7] = dfm1dk*x0[j] + dgmhdk*v0[j]
            delxv_jac[  j,  8] = dfm1dh*x0[j] + dgmhdh*v0[j]
            jac_mass[  j] = GNEWT^2*rinv*r0inv*(dfm1dk2*x0[j]+dgmhdk2*v0[j])
        end
        # Derivatives of \delta v with respect to x0, v0, k & h:
        c5 = (r0-k*g2)*rinv/g1
        c7 = g2*(1/g1+c2*rinv)
        c8 = (k*c6+r*r0+c3*c5)*r0inv3
        c12 = g0*h-g1*r0
        c20 = k*(g2*k+r)-g0*r0*zeta
        ddfdtdxx = dfdt*(eta*g1*rinv-2-g0*c3*rinv*r0inv/g1+c2*c3*r0inv*rinv2-k*(k*g2-r0)*betainv*rinv*r0inv)*r0inv2
        ddfdtdxv = -dfdt*(g0*g2/g1+(r0*g1+eta*g2)*rinv)*rinv
        ddfdtdvx = ddfdtdxv
        ddfdtdvv = -k*rinv3*r0inv*betainv*((beta*eta*g2^2+k*h8)*(r0*g0+k*g2)+
                                            g1*(- h6*ksq + (-2eta*g1*g2+(h1-2g1*g3)*k)*beta*r0 - beta*g1^2*r0^2))
        ddfdtdk  = dfdt*(1/k+c1*(r0-g2*k)*r0inv*rinv2/g1-betainv*r0inv*(1+c17*rinv))
        ddfdtdk2  = (r0-g2*k)*betainv*r0inv*rinv2*(-eta*beta*g2^2+h3*k+(g3-g1*g2)*beta*r0)
        ddfdtdh  = dfdt*(r0-g2*k)*rinv2/g1
        dgdotm1dxx = rinv2*r0inv3*((eta*g2+g1*r0)*k*c3*rinv+g2*k*(k*(g2*k-r)-g0*r0*zeta)*betainv)
        dgdotm1dxv = k*g2*rinv3*(r*g1+r0*g1+eta*g2)
        dgdotm1dvx = dgdotm1dxv
        dgdotm1dvv = k*betainv*rinv3*(eta^2*beta*g2^3-eta*k*g2*h3+3r0*eta*beta*g1*g2^2 +
                                       r0*k*(-g0*h6+3beta*g1*g2*g3)+beta*g2*(g0*g2+g1^2)*r0^2)
        dgdotm1dk = rinv*r0inv*(-r0*g2+g2*k*(r+r0-g2*k)*betainv*rinv-k*g1*c1*rinv+k*g2*c1*c2*rinv2)
        dgdotm1dk2 = betainv*rinv2*(-beta*eta^2*g2^3+eta*k*g2*h3+eta*r0*beta*g2*(g3-2g1*g2)+
                                     (h6-beta*g2^3)*r0*k + beta*g1*(g3 - g1*g2)*r0^2)
        dgdotm1dh = k*rinv3*(g2*c2-r*g1)
        @inbounds for j=1:3
            # Next, compute the \delta v-derivatives:
            delxv_jac[3+j,  j] = dfdt
            delxv_jac[3+j,3+j] = dgdtm1
            @inbounds for i=1:3
                delxv_jac[3+j,  i] += (ddfdtdxx*x0[i]+ddfdtdxv*v0[i])*x0[j] + (dgdotm1dxx*x0[i]+dgdotm1dxv*v0[i])*v0[j]
                delxv_jac[3+j,3+i] += (ddfdtdvx*x0[i]+ddfdtdvv*v0[i])*x0[j] + (dgdotm1dvx*x0[i]+dgdotm1dvv*v0[i])*v0[j]
            end
            delxv_jac[3+j,  7] = ddfdtdk*x0[j] + dgdotm1dk*v0[j]
            delxv_jac[3+j,  8] = ddfdtdh*x0[j] + dgdotm1dh*v0[j]
            jac_mass[3+j] = GNEWT^2*rinv*r0inv*(ddfdtdk2*x0[j]+dgdotm1dk2*v0[j])
        end
        if debug
            # Now include derivatives of gamma, r, fm1, gmh, dfdt, and dgdtmhdfdtm1:
            @inbounds for i=1:3
                delxv_jac[ 7,i] = -sqb*rinv*(g2*v0[i]+c3*x0[i]*r0inv3); delxv_jac[7,3+i] = -sqb*rinv*(d*v0[i]+g2*x0[i])
                delxv_jac[ 8,i] = (c20*betainv-c2*c3*rinv)*r0inv3*x0[i]+((r0*g1+eta*g2)*rinv)*v0[i]
                delxv_jac[8,3+i] = (c18*betainv-d*c2*rinv)*v0[i]+((r0*g1+eta*g2)*rinv)*x0[i]
                delxv_jac[ 9,i] = dfm1dxx*x0[i]+dfm1dxv*v0[i]; delxv_jac[ 9,3+i]=dfm1dvx*x0[i]+dfm1dvv*v0[i]
                delxv_jac[10,i] = ddfdtdxx*x0[i]+ddfdtdxv*v0[i]; delxv_jac[10,3+i]=ddfdtdvx*x0[i]+ddfdtdvv*v0[i]
                delxv_jac[11,i] = dgmhdxx*x0[i]+dgmhdxv*v0[i]; delxv_jac[11,3+i]=dgmhdvx*x0[i]+dgmhdvv*v0[i]
                delxv_jac[12,i] = dgdotm1dxx*x0[i]+dgdotm1dxv*v0[i]; delxv_jac[12,3+i]=dgdotm1dvx*x0[i]+dgdotm1dvv*v0[i]
            end
            delxv_jac[ 7,7] = sqb*c1*r0inv*rinv; delxv_jac[7,8] = sqb*rinv
            delxv_jac[ 8,7] = betainv*r0inv*rinv*(-g2*r0^2*beta-eta*g1*g2*(k+beta*r0)+eta*g0*g3*(2*k+zeta)-
                                                  g2^2*(beta*eta^2+2*k*zeta)+g1*g3*zeta*(3*k-beta*r0)); delxv_jac[8,8] = c2*rinv
            delxv_jac[ 9,7] = dfm1dk; delxv_jac[ 9,8] = dfm1dh
            delxv_jac[10,7] = ddfdtdk; delxv_jac[10,8] = ddfdtdh
            delxv_jac[11,7] = dgmhdk; delxv_jac[11,8] = dgmhdh
            delxv_jac[12,7] = dgdotm1dk; delxv_jac[12,8] = dgdotm1dh
        end
    end
    return
end


text:
mutable struct FSpaceXD
	ndim       :: Integer                                  # number of dimensions
	coeff      :: Array{Float64}		                   # coefficient vector/
	basis      :: Dict{Integer,BSpline}	                   # dict with spline objects
	curr_basis :: Dict{Int,SparseMatrixCSC{Float64,Int64}} # dict with current spline evals

	function FSpaceXD(ndim::Integer,coeff::Array{Float64},bs::Dict{Integer,BSpline})
		if ndim != length(bs)
			throw(ArgumentError("number of basis functions must be equal to ndim"))
		end
		cols = 1
		for (k,v) in bs
			cols *= getNumCoefs(bs[k])
		end
		if cols != size(coeff,1)
			throw(ArgumentError("Dimension mismatch: size(coeff,1) must be equal to the product of all spline coefficients"))
		end
		cb = Dict{Int,SparseMatrixCSC{Float64,Int64}}()
		for i in ndim
			cb[i] = spzeros(getNumCoefs(bs[i]),1)
		end

		return new(ndim,coeff,bs,cb)
	end
end

function getValue(fx::FSpaceXD,x::Array{Float64,1})

	# 0) check length of x
	if length(x) != fx.ndim
		throw(ArgumentError("length of x must be equal to fx.ndim"))
	end

	# 1) evaluate each spline[j] at x[j]
	for i in 1:fx.ndim
		fx.curr_basis[i] = getBasis(x[i],fx.basis[i])
	end

	# 2) use evalTensorXX to evaluate at correct coeff vector (in case there are many)
	evalTensor(fx.curr_basis,fx.coeff)
end

function show(io::IO,f::FSpaceXD)
	print(io,"FSpaceXD with $(f.ndim) dimensions\n")
end


text:
using Test
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

@testset "GLColorBuffer.jl" begin

   # function_name()
   @testset "function_name()" begin
      @test
      @test
      @test
      @test
   end

end


text:
using LickometerAnalysis
using Test

@testset "LickometerAnalysis.jl" begin
    # Write your tests here.
end