text:
# scratch.jl
# scratch work ray-trace
#

##
#= Functions:
	phase_refractive_index!():
		inputs:	μ, dμdψ, r, θ, χ, freq
		outputs: mutates μ, dμdψ

	ddr!():
		inputs:	μ, dμdψ, r, θ, χ, freq, dμdr
		outputs: mutates dμdr
		calls: phase_refractive_index!()

	ddθ!():
		inputs:	μ, dμdψ, r, θ, χ, freq, dμdθ
		outputs: mutates dμdθ
		calls: phase_refractive_index!()

	ddχ!():
		inputs:	μ, dμdψ, r, θ, χ, freq, dμdχ
		outputs: mutates dμdχ
		calls: phase_refractive_index!()

	ddf!():
		inputs:	μ, dμdψ, r, θ, χ, freq, dμdf
		outputs: mutates dμdf
		calls: phase_refractive_index!()

	haselgrove!():
		inputs: μ, dμdψ, r, θ, χ, freq, t, dμ, ddt
				dμdr,dμdθ,dμdχ,dμdf = dμ
				ddt = [drdt, dθdt, dχdt]
		outputs: mutates ddt

=#

#=TODO:
	1. Get functions and argument-passing working with single input (i.e. non-vectorized, non-looping)
	2. Vectorize!
	3. Run and time solver
	4. Switch to in-place functions where possible
	5. Run and time solver, optimize
=#
#const rₑ = 6378100
c = 2.99792458e8
re = 6.3712e6      	# radius of Earth in meters
B0 = 3.0696381e-5   # magnitude of B field at Earth's equator surface in Tesla
e = 1.602e-19 		# electric charge in Coulombs
me = 9.1093e-31		# electron rest mass
mp = 1.6726219e-27	# proton rest mass
eps = 8.854e-12 	# permittivity of free space in mks

using DifferentialEquations

# calculate phase refractive index μ
# DONE! todo: convert this function to mutating form
# 	QUESTION: is this possible, given that phase_refractive_index!() must be called multiple times in a single ddr!() (e.g.) call, before using the result?
# 	ANSWER: yes! see call/assignment in ddr!()

u = function phase_refractive_index(r,θ,χ,freq)
# u = [μ, dμdψ]
    # convert from radial angle to wave normal angle
    dip = atan(2.0*cot(θ))     			# dip angle: angle between horizontal and B field;
	# 11/9: fixed error tan(pi/2.0-θ) ➡ tan(θ)
	# 11/13: no! dip = 0 when θ = π/2 ➡ tan(θ - π/2) is correct!
	# 11/18: wrong again! Dip is defined as positive toward the North pole, i.e. should be positive in the Northern hemisphere!
	# 11/19: really need to check this; looks like factor of 2 was probably not correct
	#ψ = 2*π - (3π/2 - dip - χ) 		# this should be ≡ to the below, but the plotted ray paths are not all that similar!
	#ψ = π/2.0 + dip + χ
	#ψ = χ - π/2 - dip
	ψ = (3.0/2.0)*pi - dip - χ
	#ψ = χ - (3.0/2.0)*pi + dip		# wave normal angle: angle between wave vector and B; 11/9: fixed error (χ - ϕ) ➡ (ϕ - χ)
	# NOTE 11/24/2020: Rice 1997 uses ψ = χ - 3π/2 + dip, as well as a strange method of calculating χ (see FORTRAN code on p. 112)

    # convert from frequency to angular frequency
    ω = 2.0*pi*freq

    # find magnetic field at (r,θ) from dipole field model
    # Bmag = B0*((re^3)/(r^3))*sqrt(4.0-3.0*cos((pi/2.)-θ)*cos((pi/2.0)-θ))
	# 	   = B0*((re^3)/(r^3))*sqrt(4 - 3sin²θ)
	# 	   ≡ B0*((re^3)/(r^3))*sqrt(1 + 3cos²θ)
	# dipole field magnitude from Parks p. 61:
	# B(r,λ) = μ₀M/(4πr³)*(1+3sin²λ)^(1/2)
	# 		 = μ₀M/(4πr³)*(1+3sin²(π/2 - θ))^(1/2)
	# 		 = μ₀M/(4πr³)*(1+3cos²θ)^(1/2)
	# here B0*re³ = μ₀M/(4π)
	Bmag = B0*(re^3/(r^3))*sqrt(1+3*cos(θ)*cos(θ))

    # calculate electron and proton density profiles
	# isotropic ionosphere
	ne_iono = (1.8e5*exp(-4.183119*((r/re)-1.0471))) # cm^-3
	#n_p = 1.e6*(1.8e5*exp(-4.183119*((r/re)-1.0471)))

	# Carpenter and Anderson 1992 plasmasphere
	λ = π/2 - θ
	Lshell = r./(re*cos(λ)^2)

	if Lshell <= Lppi
	     log_ne = (-0.3145*Lshell + 3.9043) + (0.15*(cos((2*π*(d+9))/365) - 0.5*cos((4*π*(d+9))/365)) + 0.00127*R̄ - 0.0635)*exp((2-Lshell)/1.5)
	     ne_plasma = 10^(log_ne)
	elseif Lppi < Lshell <= Lppo
	     ne_plasma = ne_Lppi*10^((Lppi-Lshell)/0.1)
	elseif Lppo < Lshell
	     ne_plasma = (5800 + 300*t)*Lshell^(-4.5) + (1 - exp((2-Lshell)/10))
	else
	     ne_plasma = 0.0
	end

	n_e = (ne_iono + ne_plasma)*1e6
	n_p = n_e

	# electron plasma angular frequency squared
    ω_e2 = (n_e*(e^2))/(eps*me)
	# proton plasma angular frequency squared
    ω_p2 = (n_p*(e^2))/(eps*mp)

	# electron angular gyrofrequency
    Ω_e = (e*Bmag)/(me)
	# proton angular gyrofrequency
    Ω_p = (e*Bmag)/(mp)

	# define R, L, P, D, S,
	# R ≡ 1 - Σ(ωₖ²/ω²)(ω/(ω+ϵₖΩₖ))
	R = 1.0 - (ω_e2/ω^2.0)*(ω/(ω - Ω_e)) - (ω_p2/ω^2.0)*(ω/(ω + Ω_p))
	#println("R = ",R)

	# L ≡ 1 - Σ(ωₖ²/ω²)(ω/(ω-ϵₖΩₖ))
	L = 1.0 - (ω_e2/ω^2.0)*(ω/(ω + Ω_e)) - (ω_p2/ω^2.0)*(ω/(ω - Ω_p))
	#println("L = ",L)

	# P ≡ 1 - Σ(ωₖ²/ω²)
	P = 1.0 - (ω_e2/ω^2.0) - (ω_p2/ω^2.0)
	#println("P = ",P)

	# D ≡ ½(R-L); S ≡ ½(R+L)
	D = (R - L)/2.0
	S = (R + L)/2.0

	# define parts of dispersion relation
	# Aμ⁴ - Bμ² + C = 0
	# A = Ssin²ψ + Pcos²ψ
	A = S*sin(ψ)^2.0 + P*cos(ψ)^2.0	# ψ = π/2: A = S

	# B = RLsin²ψ + PS(1 + cos²ψ)
	B = R*L*sin(ψ)^2.0 + P*S*(1.0+cos(ψ)^2.0) # ψ = π/2: B = RL + PS

	# C  = PRL
	C = P*R*L

	# solve the dispersion relation for μ:
	# μ² = (B +- F)/2A
	# where F² = (RL-PS)²sin⁴ψ + 4P²D²cos²ψ
	#QUESTION: is it faster to (a) solve for F², then sqrt(); or (b) solve for F directly? (a) requires fewer sqrt() calls
	F2 = (R*L - P*S)^2.0*sin(ψ)^4.0 + 4.0*(P*D*cos(ψ))^2.0
		# ψ = π/2: F² = (RL - PS)² ➡ F = RL - PS
	F = sqrt(F2)

	# Rice 1997: typical solution to dispersion relation with quadratic formula
	μ2_minus = (B - F)/(2.0*A)
	μ2_plus = (B + F)/(2.0*A)
		# ψ = π/2: μ²₊ = (RL + PS + RL - PS)/2A
		# 			   = RL/A
		# 			   = RL/S
		# 			   = 2RL/(R + L)

	# Bortnik 2004:
	# if B > 0
	# 	μ2_minus = (B - F)/(2.0*A)
	# else
	# 	μ2_plus = (2.0*C)/(B + F)
	# end

	μ_minus = try
		sqrt(abs(μ2_minus))	# abs() is not physical! for test only

	catch err
		if isa(err,DomainError)
			println("DomainError in μ_minus!")
			println("B=", B)
			println("F=", F)
			println("A=", A)
		else
			# println(err)
		end
	end

	μ_plus = try
		sqrt(abs(μ2_plus)) # abs() is not physical! for test only

	catch err
		if isa(err,DomainError)
			# println("DomainError in μ_plus!")
			# println("B=", B)
			# println("F=", F)
			# println("A=", A)
		else
			# println(err)
		end
	end

	# Electron whistler case: ψ = 0, μ² = R; this is the μ_plus case
	#μ = μ_minus		# EMIC case
	μ = μ_plus		# electron whistler case

	# Find dA/dψ, dB/dψ, dC/dψ, dμ/dψ
	dAdψ = 2.0*(S-P)*sin(ψ)*cos(ψ) 		# ψ = π/2: dAdψ = 0
	dBdψ = 2.0*(R*L-P*S)*sin(ψ)*cos(ψ) 	# ψ = π/2: dBdψ = 0
    dCdψ = 0.0
    #dμdψ = ((μ^4.0)*dAdψ-(μ^2.0)*dBdψ+dCdψ)/(4.0*A*(μ^3.0)-2.0*B*μ)

	dFdψ = 1/(2.0*F)*((R*L-P*S)^2 * 4*sin(ψ)^3*cos(ψ) - 8*(P*D)^2*sin(ψ)*cos(ψ)) # ψ = π/2: dFdψ = 0
	#dFdψ = sqrt(abs((R*L-P*S)^2 * 4*sin(ψ)^3*cos(ψ) - 8*(P*D)^2*sin(ψ)*cos(ψ)))
	dμdψ = 1/(2.0*μ)*((dBdψ + dFdψ)/(2*A) - 2*dAdψ*(B + F)/(2*A^2))
		# ψ = π/2: dμdψ = 0
	# NOTE 11/24/2020: from Rice 1997; choose '+' for (B +- F) and (dBdψ +- dFdψ)

	#DEBUG check values
	#println("μ = ",μ)
	#println("dμdψ = ", dμdψ)

	# non-mutating output
	u = [μ,dμdψ]

end

# Derivatives of phase refractive index w.r.t. r, θ, χ,
# calculate derivative of μ w.r.t. r
dμdr = function ddr(r,θ,χ,freq)
	#μ,dμdψ = u			# unpack
	dr = 1.0e-11		# r step size

	μ_l = phase_refractive_index(r-dr/2.0,θ,χ,freq)
	#μ_l = μ
	μ_r = phase_refractive_index(r+dr/2.0,θ,χ,freq)
	#μ_r = μ

	dμdr = (μ_r[1] - μ_l[1])/dr
end

dμdθ = function ddθ(r,θ,χ,freq)
	#μ,dμdψ = u			# unpack
	dθ = 1.0e-11		# θ step size

	μ_l = phase_refractive_index(r,θ-dθ/2.0,χ,freq)
	#μ_l = μ
	#println("μ_l = ", μ)
	μ_r = phase_refractive_index(r,θ+dθ/2.0,χ,freq)
	#μ_r = μ
	#println("μ_r = ", μ)

	dμdθ = (μ_r[1] - μ_l[1])/dθ
end

dμdχ = function ddχ(r,θ,χ,freq)
	#μ,dμdψ = u			# unpack
	dχ = 1.0e-11		# χ step size

	μ_l = phase_refractive_index(r,θ,χ-dχ/2.0,freq)
	#μ_l = μ
	μ_r = phase_refractive_index(r,θ,χ+dχ/2.0,freq) # NOTE Nov 3: corrected error 'χ-dχ/2.0'➡'χ+dχ/2.0'
	#μ_r = μ

	dμdχ = (μ_r[1] - μ_l[1])/dχ
end

dμdf = function ddf(r,θ,χ,freq)
	#μ,dμdψ = u			# unpack
	df = 1.0e-5		# f step size

	μ_l = phase_refractive_index(r,θ,χ,freq-df/2.0)
	#μ_l = μ
	μ_r = phase_refractive_index(r,θ,χ,freq+df/2.0)
	#μ_r = μ

	dμdf = (μ_r[1] - μ_l[1])/df
end

# calculate derivatives w.r.t. time
function haselgrove!(du,u,p,t)
# 	du = drdt, dθdt, dχdt, dfdt
# 	u = r, θ, χ, f
#  	p = freq, μ, dμdψ, dμdr, dμdθ, dμdχ, dμdf
# 		^ this compiles and runs, but parameters
# 	t = time

	r, θ, χ, freq = u 						# unpack u
	#freq = p[1] 	# unpack p

	dμdr = ddr(r,θ,χ,freq)			# mutating form: ddr!([μ,dμdψ],r,θ,χ,freq,dμdr)
	dμdθ = ddθ(r,θ,χ,freq)			# calls phase_refractive_index(r,θ,χ,freq)
	dμdχ = ddχ(r,θ,χ,freq)
	dμdf = ddf(r,θ,χ,freq)

	v = phase_refractive_index(r,θ,χ,freq)
	μ = v[1]
	dμdψ = v[2]

	# 11/20: try dμdχ --> -1*dμdψ; these should be equal but are calculated differently
    du[1] = 1/(μ^2)*(μ*cos(χ)+dμdψ*sin(χ))
    du[2] = 1/(r*μ^2)*(μ*sin(χ)-dμdψ*cos(χ))
    du[3] = 1/(r*μ^2)*(dμdθ*cos(χ)-(r*dμdr + μ)*sin(χ))
	du[4] = 1/c*(1+(freq/μ)*dμdf)

	#du = [drdt, dθdt, dχdt, dfdt]

	#println("drdt = ", drdt)
	#println("dθdt = ", dθdt)
	#println("dχdt = ", dχdt)

	#DEBUG check values
	#println("drdt = ",drdt)
	#println("dθdt = ",dθdt)
	#println("dχdt = ",dχdt)

    # dμdψ = 1/(2*μ)*((dBdψ + dFdψ)/(2*A) - 2*dAdψ*(B+F)/(2*A^2))
    # dμdω = 1/(2*μ)*((dBdω + dFdω)/(2*A) - 2*dAdω*(B+F)/(2*A^2))
	#
    # dμdr = (μ*(r+δr/2) - μ*(r-δr/2))/δr
    # dμdθ = (μ*(θ+δθ/2) - μ*(θ-δθ/2))/δθ

end

#DEBUG check values
#haselgrove!(μ,dμdψ,r,θ,χ,freq,dμ,ddt,t)
#haselgrove!(1.0,1.0,10000.0,pi/4.0,0.0,5000.0,[1.0,1.0,1.0,1.0],[1.0,1.0,1.0],0.0)

# du0 = [1.0, 1.0, 1.0, 1.0]
# u0 = [8000.0, pi/4.0, 0.0, 5000.0]
# p = []
# t0 = 0.0
# haselgrove!(du0, u0, p, t0)
# println(du0)

## callbacks!
# re_cb: terminate when ray intersects Earth surface
function re_term_condition(u,t,integrator)
	u[1] - re
end

# function μ2_term_condition(u,t,integrator)
# 	body
# end

function terminate_affect!(integrator)
	terminate!(integrator)
end

re_cb = ContinuousCallback(re_term_condition,terminate_affect!)

# saveμ_cb: save μ, dμ values at each timestep
function save_func(u,t,integrator)
	r, θ, χ, f = u
	v = phase_refractive_index(r, θ, χ, f)
	dip = atan(2.0*cot(θ))
	ψ = π/2.0 + dip + χ
	vcat(v,dip,ψ)
end

saved_μ = SavedValues(Float64, Array{Float64,1})
saveμ_cb = SavingCallback(save_func, saved_μ)

cb = CallbackSet(re_cb, saveμ_cb)

## ODE solver!
#QUESTION: how does ODEProblem work?

using BenchmarkTools
using LSODA
using Sundials

u0 = [re+1.0e+6, pi/3, 0.0, 1000.0]					# r0, θ0, χ0
# r_pw = 7.90441264e6
# θ_pw = 1.848517482
# χ_pw = 0.0*π/2.0
# f_pw = 1000.0
# u0 = [r_pw, θ_pw, χ_pw, f_pw]					# r0, θ0, χ0

p = []	# f0, dμdψ, dμdr, dμdθ, dμdχ, dμdf
tspan = (0.0,5.0e+10)

hasel_prob = ODEProblem(haselgrove!,u0,tspan,p, callback=cb)
hasel_soln = solve(hasel_prob, CVODE_BDF(), reltol=1e-7) #reltol=1e-4
# LSODA_BDF()?

using Plots
plotly()
plot(hasel_soln)
plot(hasel_soln, vars=(1,2,3))

t = hasel_soln.t
r = hasel_soln[1,:]
θ = hasel_soln[2,:]
χ = hasel_soln[3,:]
f = hasel_soln[4,:]

x = r.*sin.(θ)
y = r.*cos.(θ)

plot(x,y, aspect_ratio=1, legend=:none)

#earthx = re/1.0e3.*cos([0:0.01:2*π;])
plot!(re.*cos.([0:0.01:2*π;]),re.*sin.([0:0.01:2*π;]), aspect_ratio = 1)
# plot!(re.*[0:0.01:2;],re.*[0:0.01:2;], aspect_ratio = 1)
# plot!(re.*[0:0.01:4;],re.*[0:0.005:2;], aspect_ratio = 1)

## Plot callback saved_μ

t_s = saved_μ.t
val_s = saved_μ.saveval
val_s_r = reduce(hcat, val_s)
val_s_r = val_s_r'
μ_s = val_s_r[:,1]
dμdψ_s = val_s_r[:,2]
dip_s = val_s_r[:,3]
ψ_s = val_s_r[:,4]

plot(t_s, μ_s, title="refractive index v. time")
plot(t_s, dμdψ_s, title="d/dpsi refractive index v. time")
plot(t_s, dip_s.*(180/pi), title="dip angle v. time")
plot(t_s, ψ_s.*(180/pi), title="wave normal angle v. time")

## Plot magnetic field lines
gridrange = [-pi/2.0:0.01:pi/2.0;]
#thtgrid = np.arange(0.,np.pi,0.1)
gridmag = pi/2.0.-gridrange
rmag = re.*(cos.(gridmag)).^2

#rw = 2.*ray_cmks.Re
xmag1 = 1.0.*rmag.*cos.(gridmag)
ymag1 = 1.0.*rmag.*sin.(gridmag)

xmag2 = 2.0.*rmag.*cos.(gridmag)
ymag2 = 2.0.*rmag.*sin.(gridmag)

xmag3 = 3.0.*rmag.*cos.(gridmag)
ymag3 = 3.0.*rmag.*sin.(gridmag)

xmag4 = 4.0.*rmag.*cos.(gridmag)
ymag4 = 4.0.*rmag.*sin.(gridmag)

plot!(xmag4,ymag4,color = :red)
plot!(xmag3,ymag3,color = :red)
plot!(xmag2,ymag2,color = :red)
plot!(xmag1,ymag1,color = :red)

plot!(xmag4,-ymag4,color = :red)
plot!(xmag3,-ymag3,color = :red)
plot!(xmag2,-ymag2,color = :red)
plot!(xmag1,-ymag1,color = :red)

θ_test = pi/2
plot!([-2e7:1e5:2e7;].*sin(θ_test),[-2e7:1e5:2e7;].*cos(θ_test))
x_dip = [0:1e5:1e6;].*sin(dip_test +pi/2 - θ_test)
y_dip = [0:1e5:1e6;].*cos(dip_test +pi/2 - θ_test)
x_dip4 = x_dip .+ 4.0*re*sin(θ_test)^2*sin(θ_test)
y_dip4 = y_dip .+ 4.0*re*sin(θ_test)^2*cos(θ_test)
plot!(x_dip4,y_dip4,aspect_ratio=1)

## calculate and plot refractive index surface
# initial conditions
x_test = 1*9.031E+6
y_test = 1*1.537E+6
θ_test = atan(abs(y_test/x_test))
r_test = sqrt(x_test^2 + y_test^2) #2.0*re*(sin(θ_test))^2 # on plotted B field line
#θ_test = θ_pw
#r_test = r_pw
ψ_test = [0:0.01:2π;]
dip_test = atan(2.0*cot(θ_test))
#χ_test = ψ_test .+ 3π/2 .- dip_test
χ_test = -1.0.*ψ_test .+ 3.0*π/2 .- dip_test
f_test = 1000.0

# calculate μ, dμdψ
u_test = phase_refractive_index.(r_test, θ_test, χ_test, f_test)
u_test_r = reduce(hcat,u_test)
u_test_r = u_test_r'
μ_test =  u_test_r[:,1]
dμdψ_test = u_test_r[:,2]

# xμ_test = μ_test.*sin.(ψ_test .+ π/2)
# yμ_test = μ_test.*cos.(ψ_test .+ π/2)
# # plot surface without rotation into B|| frame
xμ_test = μ_test.*sin.(χ_test)
yμ_test = μ_test.*cos.(χ_test)

xdμ_test = dμdψ_test.*sin.(ψ_test)
ydμ_test = dμdψ_test.*cos.(ψ_test)

xμ_test_fieldline = 4.0*(re*(sin(θ_test))^2)*sin(θ_test) .+ xμ_test.*1e5
yμ_test_fieldline = 4.0*(re*(sin(θ_test))^2)*cos(θ_test) .+ yμ_test.*1e5

plot(xμ_test, yμ_test, aspect_ratio = 1)
#plot(xdμ_test, ydμ_test, aspect_ratio = 1)

plot!(xμ_test_fieldline, yμ_test_fieldline, aspect_ratio = 1)|

plot!([-30:1:30;].*cos(dip_test),[-30:1:30;].*sin(dip_test))
plot!([-30:1:30;].*sin(θ_test+π/2),[-30:1:30;].*cos(θ_test+π/2))

##
function lorenz!(du, u, p, t)
    du[1] = 10.0*(u[2]-u[1])
    du[2] = u[1]*(28.0-u[3]) - u[2]
    du[3] = u[1]*u[2] - (8/3)*u[3]
	du[4] = (du[1] + du[2] + du[3])
end
# call this function in a problem
u0 = [1.0;0.0;0.0;1.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan)
sol = solve(prob)

plot(sol,vars=(1,2,3))
plot(sol,vars=4)

## B-field check

re = 6.3712e6      	# radius of Earth in meters
B0 = 3.0696381e-5	# magnetic field strentgh (where?) in Tesla
r = 5000.0
θ = pi/4.0
Bmag = B0*((re^3)/(r^3))*sqrt(4.0-3.0*cos((pi/2.)-θ)*cos((pi/2.0)-θ))

##

out = function fun(in1,in2)
	out = in1*in2
end

a = fun(1.0,2.0)

function mutefun!(out::Array{Int64},in1,in2)
	out .= in1*in2		# don't forget '.' notation for vectorized i/o!
end

function right1_fill_with_twos!(v::Vector{Int64})
	v[:] =[2 for ii in 1:length(v)]
end

a = [1:5;]

function varargs(args...)
	return args
end
varargs(1,2,3)

a = [1:5;]
sum(a)

## event/callback scratch
# see https://tutorials.sciml.ai/html/introduction/04-callbacks_and_events.html
using DifferentialEquations, ParameterizedFunctions

# define a ballistic projectile, i.e. a bouncy ball
ball! = @ode_def BallBounce begin
	dy = v	# d/dposition = velocity
	dv = -g # d/dvelocity = -1*acceleration due to gravity
end g

function condition(u,t,integrator)
	u[1]
end

# make the ball bounce!
# integrator.u : position, velocity, acceleration
# integrator.t : time
# integrator.p : parameters (g,friction coeff)
function affect!(integrator)
	integrator.u[2] = -integrator.p[2] * integrator.u[2]
	integrator.p[2] = sqrt(integrator.p[2]) # friction parameter changes with each bounce
end

# build the callback!
# when condition: y == 0, affect! occurs (velocity changes sign, minus loss from friction coeff)
bounce_cb = ContinuousCallback(condition, affect!)

# now add a discrete callback: when time = 2, some kid kicks the ball, adding 50 to its velocity
function condition_kick(u,t,integrator)
	t == 2
	#t - 2 # continuous callback version: when t - 2 --> 0, condition_kick is triggered
end

function affect_kick!(integrator)
	integrator.u[2] += 50
end

# build the discrete callback!
kick_cb = DiscreteCallback(condition_kick,affect_kick!)

# build a CallbackSet in order to combine the bounce_cb and kick_cb
cb = CallbackSet(bounce_cb,kick_cb)

# make the ODEProblem with callback
u_0 = [50.0,0.0]
t_span = (0.0,15.0)
p_ = [9.8,0.9]
prob_ = ODEProblem(ball!,u_0,t_span,p_,callback=cb)
# with bounce_cb only: prob_ = ODEProblem(ball!,u_0,t_span,p_,callback=bounce_cb)
sol_ = solve(prob_,Tsit5(),tstops=[2.0]) #note: need to specify that intergration scheme should step at exactly t == 2 in order to trigger discrete callback

plot(sol_)

# integration termination and directional handling
# example: harmonic oscillator
harmonic! = @ode_def HarmonicOscillator begin
	dv = -x
	dx = v
end

#terminate the integration when a condition is met

function terminate_condition(u,t,integrator)
	u[2]
end

function terminate_affect!(integrator)
	terminate!(integrator)
end

terminate_cb = ContinuousCallback(terminate_condition,terminate_affect!)
terminate_upcrossing_cb = ContinuousCallback(terminate_condition,terminate_affect!,nothing) # if two affect!s are given in a Callback build, the first is for "upcrossings" (i.e. when the condition trigger crosses 0 from below), and the second is for "downcrossings" (when the condition trigger crosses 0 from above)

u_0 = [1.0,0.0]
t_span = (0.0,10.0)
prob_ = ODEProblem(harmonic!,u_0,t_span,callback=terminate_upcrossing_cb)
sol_ = solve(prob_) #NOTE: can also add callbacks to solve() command! this allows us to distinguish between model feature (callbacks in problem statement) and integration commands (callbacks in solve command)!
plot(sol_)

## type/value check
r = re + 1.0e6
θ = pi/4.0
χ = 0.0
freq = 5000.0
dip = atan(2.0*cot(θ))
dip*180/pi	# dip angle: angle between horizontal and B field
ϕ = (3.0/2.0)*pi - dip              # intermediate angle -- NOT azimuth
ψ = ϕ - χ							# wave normal angle: angle between wave direction and B
ψ*180/π

# convert from frequency to angular frequency
ω = 2.0*pi*freq

# convert radius to fraction of Earth radius
(r/re) = r/re
re

# find magnetic field at (r,θ) from dipole field model
B0
Bmag = B0*((re^3)/(r^3))*sqrt(4.0-3.0*cos((pi/2.0)-θ)*cos((pi/2.0)-θ))
B0*((re^3)/(r^3))
sqrt(4.0-3.0*cos((pi/2.0)-θ)*cos((pi/2.0)-θ))
# calculate electron and proton density profiles
n_e = 1.e6*(1.8e5*exp(-4.183119*((r/re)-1.0471)))
n_p = 1.e6*(1.8e5*exp(-4.183119*((r/re)-1.0471)))

# electron plasma angular frequency squared
ω_e2 = (n_e*(e^2))/(eps*me)
# proton plasma angular frequency squared
ω_p2 = (n_p*(e^2))/(eps*mp)

# electron angular gyrofrequency
Ω_e = (e*Bmag)/(me)
# proton angular gyrofrequency
Ω_p = (e*Bmag)/(mp)

# define R, L, P, D, S,
# R ≡ 1 - Σ(ωₖ²/ω²)(ω/(ω+ϵₖΩₖ))
R = 1.0 - (ω_e2/ω^2.0)*(ω/(ω - Ω_e)) - (ω_p2/ω^2.0)*(ω/(ω + Ω_p))
(ω_e2/ω^2.0)*(ω/(ω - Ω_e))
(ω_p2/ω^2.0)*(ω/(ω + Ω_p))
#println("R = ",R)

# L ≡ 1 - Σ(ωₖ²/ω²)(ω/(ω-ϵₖΩₖ))
L = 1.0 - (ω_e2/ω^2.0)*(ω/(ω + Ω_e)) - (ω_p2/ω^2.0)*(ω/(ω - Ω_p))
(ω_e2/ω^2.0)*(ω/(ω + Ω_e))
(ω_p2/ω^2.0)*(ω/(ω - Ω_p))
#println("L = ",L)

# P ≡ 1 - Σ(ωₖ²/ω²)
P = 1.0 - (ω_e2/ω^2.0) - (ω_p2/ω^2.0)
#println("P = ",P)

# D ≡ ½(R-L); S ≡ ½(R+L)
D = (R - L)/2.0
S = (R + L)/2.0

# define parts of dispersion relation
# Aμ⁴ - Bμ² + C = 0
# A = Ssin²ψ + Pcos²ψ
A = S*sin(ψ)^2.0 + P*cos(ψ)^2.0

# B = RLsin²ψ + PS(1 + cos²ψ)
B = R*L*sin(ψ)^2.0 + P*S*(1.0+cos(ψ)^2.0)

# C  = PRL
C = P*R*L

# solve the dispersion relation for μ:
# μ² = (B +- F)/2A
# where F² = (RL-PS)²sin⁴ψ + 4P²D²cos²ψ
#QUESTION: is it faster to (a) solve for F², then sqrt(); or (b) solve for F directly? (a) requires fewer sqrt() calls
F2 = (R*L - P*S)^2.0*sin(ψ)^4.0 + 4.0*(P*D*cos(ψ))^2.0
F = sqrt(F2)

μ2_minus = (B - F)/(2.0*A)
μ2_plus = (B + F)/(2.0*A)

##
hbi = 1.0e6
b = -π/3
re
α = acos((re*sin(b)^2.0 + cos(b)*((re + hbi)^2.0 - re^2.0 * sin(b)^2.0)^(0.5))/(re + hbi))
α*180/π
chi = b - α
chi*180/π

## plot magnetic field strength
# both plots: need to figure out how plot a log scale surface
histogram2d(randn(10000),randn(10000),nbins = 20)

plot(contourf(randn(10,20)))

# r,θ: need to figure out how to convert grid to x,y
r_gridmag = re:1e5:5e7
θ_gridmag = 0:0.01:2*π
f_bmag(r_gridmag, θ_gridmag) = begin
	B0*(re^3/(r_gridmag^3))*sqrt(1+3*cos(θ_gridmag)*cos(θ_gridmag))
end
R_gridmag = repeat(reshape(r_gridmag,1,:), length(θ_gridmag), 1)
T_gridmag = repeat(θ_gridmag, 1, length(r_gridmag))
Z_gridmag = map(f_bmag, R_gridmag, T_gridmag)
p1 = contour(r_gridmag, θ_gridmag, f_bmag, fill = true)
p2 = contour(r_gridmag, θ_gridmag, Z_gridmag)
plot(p1,p2)

# x,y: need to figure out how to exclude area inside Earth from grid
r_gridmag = re:1e5:5e7
θ_gridmag = range(0,stop=2π,length=(length(r_gridmag)))
x_gridmag = r_gridmag.*sin.(θ_gridmag) # = -5e7:1e5:5e7
sort!(x_gridmag)
y_gridmag = r_gridmag.*cos.(θ_gridmag) #x_gridmag
sort!(y_gridmag)
f_log_bmag_xy(x_gridmag, y_gridmag) = begin
	log10(B0*(re^3/(sqrt(x_gridmag^2 + y_gridmag^2)^3))*sqrt(1+3*cos(atan(x_gridmag/y_gridmag))^2))
end
X_gridmag = repeat(reshape(x_gridmag,1,:), length(y_gridmag), 1)
Y_gridmag = repeat(y_gridmag, 1, length(x_gridmag))
Z_gridmag_xy = map(f_log_bmag_xy, X_gridmag, Y_gridmag)
p1_xy = contour(x_gridmag, y_gridmag, f_log_bmag_xy, fill = (true))
p2_xy = contour(x_gridmag, y_gridmag, Z_gridmag_xy)
plot(p1,p2)

# plot various particle density distributions
r_grid = re:1e5:10*re
λ_grid = range(0,stop=2π,length=(length(r_grid)))
#x_grid = r_grid.*sin.(θ_grid) # = -5e7:1e5:5e7
#sort!(x_grid)
#y_grid = r_grid.*cos.(θ_grid) #x_gridmag
#sort!(y_grid)
f_particles(r_grid,λ_grid) = begin
	#(cos.(π/2 - θ_grid).^2)*(1.e6*(1.8e5*exp(-4.183119*((r_grid/re)-1.0471))))

end

f_Lshell(r_grid, λ_grid) = begin
	r_grid./(re*cos(λ_grid)^2)
end

R_grid = repeat(reshape(r_grid,1,:), length(λ_grid), 1)
T_grid = repeat(λ_grid, 1, length(r_grid))
Z_grid = map(f_Lshell, R_grid, T_grid)
p1 = contour(r_grid, λ_grid, f_Lshell, fill = true)
p2 = contour(r_grid, λ_grid, Z_grid)
plot(p1,p2)

#n_equator = (cos.(0).^2).*(1.e6*(1.8e5.*exp.(-4.183119*((r_grid./re).-1.0471))))
#n_45 = (cos.(π/4).^2).*(1.e6*(1.8e5.*exp.(-4.183119*((r_grid./re).-1.0471))))
#n_pole = (cos.(π/2).^2).*(1.e6*(1.8e5.*exp.(-4.183119*((r_grid./re).-1.0471))))

#plot(r_grid, n_equator)
#plot!(r_grid, n_45)
#plot!(r_grid, n_pole)

max_stars_repo_name: andersontodds/ray-trace

max_issues_repo_name: andersontodds/ray-trace

max_forks_repo_name: andersontodds/ray-trace


text:
module Options

export parse_args, get_option, do_print_usage

using Glob: glob

isoption(a) = a isa Symbol || a isa Pair

function parse_args(args)
    options = [a for a in args if isoption(a)]
    pos_args = String[]
    for (a, g) in ((a, glob(string(a))) for a in args if !isoption(a))
        if length(g) == 0
            push!(pos_args, string(a))
        else
            append!(pos_args, g)
        end
    end
    return options, pos_args
end

optionequal(option::Symbol, name) = return option === name
optionequal(option::Pair, name) = return option.first === name
anyoptionequal(option, names...) = any(optionequal(option, n) for n in names)

function get_option(options, names...; default = nothing)
    i = findnext(o -> anyoptionequal(o, names...), options, 1)
    if i === nothing
        if !(default isa Pair)
            default = names[1] => default
        end
        return default
    end
    found = options[i]
    if !(found isa Pair)
        return found => true
    else
        return found
    end
end

function do_print_usage(options, docs)
    if get_option(options, :h, :help, default = false).second
        println(docs)
        return true
    end
    return false
end

end # module

ext: jl


text:
#=
hcat.jl
provide hcat support for LinearMaps
Jeff Fessler, University of Michigan
with advice from Daniel Karrasch per
https://github.com/Jutho/LinearMaps.jl/issues/45
=#

export hcat_lm

using LinearMaps
using LinearAlgebra: I, UniformScaling
using Base: hcat
using Test

function hcat_lm(As...) # helper routine for hcat of LinearMaps
	nrow = 0
	for i=1:length(As) # find first non-I to specify number of rows
		if As[i] != I
			nrow = size(As[i],1)
			break
		end
	end
	nrow == 0 && throw(DimensionMismatch("need one non-I")) # [I I] illegal

	for i=1:length(As) # check row consistency
		if As[i] != I
			ndims(As[i]) != 2 && throw(DimensionMismatch("need ndims=2"))
			size(As[i],1) != nrow &&
				throw(DimensionMismatch("need same # of rows, " *
					"got (nrow = $nrow, $i => $(size(As[i],1)))"))
		end
	end
	ncol = map(A -> A == I ? nrow : size(A, 2), As)
	ncol = [0; ncol...] # trick for ranger()

	ranger = i -> (1+sum(ncol[1:i])) : sum(ncol[1:(i+1)]) # ith input range
	f = x -> sum([As[i] * (@view x[ranger(i)]) for i=1:length(As)])
	fc = y -> vcat(map(A -> A'*y, As)...)

	lacks_fc = A -> isa(A, LinearMap) && A.fc == nothing
	if any(map(A -> lacks_fc(A), As)) # any LinearMap missing adjoint?
		return LinearMap(f, nrow, sum(ncol)) # no fc
	end
	fc = y -> vcat(map(A -> A'*y, As)...)
	return LinearMap(f, fc, nrow, sum(ncol))
end

# overload hcat for LinearMap
Base.hcat(A::LinearMap, B::LinearMap) = hcat_lm(A, B)

# [A I] [A I I] [A I B C] etc.
Base.hcat(A::LinearMap, Bs...) = hcat_lm(A, Bs...)

Base.hcat(I::UniformScaling{<:Number}, A::LinearMap, Bs...) =
	hcat_lm(I, A, Bs...)

# this will support a single leading matrix, but not more than one
Base.hcat(A::AbstractMatrix{<:Number}, B::LinearMap, Cs...) =
	hcat_lm(A, B, Cs...)

function hcat_lm_test()
	a = ones(3,2)
	b = zeros(3,4)

	A = LinearMap(x -> a*x, y -> a'*y, 3, 2)
	B = LinearMap(x -> b*x, y -> b'*y, 3, 4)

	C = [A B I] # now returns a new LinearMap :)
	@test size(C) == (3,2+4+3)
	@test C * ones(size(C,2)) == 3. * ones(3)
	@test C' * ones(size(C,1)) == [3.0, 3.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0]
	@test Matrix(C)' == Matrix(C')

	D = [I A] # LinearMap
	@test Matrix(D)' == Matrix(D')

	E = [a B] # LinearMap
	@test Matrix(E)' == Matrix(E')

	F = LinearMap(x -> b*x, size(b)...)
	G = [A F] # LinearMap
	@test Matrix(G) == [a b]
	true
end

"""
`hcat_lm(:test)`
self test
"""
function hcat_lm(test::Symbol)
	test != :test && throw("hcat_lm test")
	@test hcat_lm_test()
	true
end


text:
import StatsBase: RealVector, RealMatrix
using KernelEstimator, Optim

"""
    bwlscv(xdata::RealVector, kernel::Function)

Optimal bandwidth estimation using least square cross validation.
Adapted from https://github.com/panlanfeng/KernelEstimator.jl/blob/master/src/bandwidth.jl

# Examples
```julia-repl
julia> bar([1, 2], [1, 2])
1
```
"""
function bwlscv(xdata::RealVector, kernel::Function)
    n=length(xdata)
    w  = zeros(n)
    h0=bwnormal(xdata)
    #when n is large, leaveoneout is more expensive than numeric integration
    if kernel == gaussiankernel && n<200
        return Optim.minimizer(Optim.optimize(h -> KernelEstimator.Jh(xdata, h, w, n), 0.01*h0, 10*h0, iterations=200, abs_tol=h0/n))
    end

    xlb, xub = extrema(xdata)
    h0=KernelEstimator.midrange(xdata)
    if h0==0 #Algorithm returns 0 when h0=0. This happens when there are many ties, mixdrange can return 0
        h0 = median(xdata)
    end
    hlb = h0/n
    hub = h0
    if kernel == betakernel
        xlb = 0.0
        xub = 1.0
        hub = 0.25
        logxdata = log.(xdata)
        log1_xdata = log.(1.0 .- xdata)
        return Optim.minimizer(Optim.optimize(h -> KernelEstimator.Jh(xdata, logxdata, log1_xdata, kernel, h, w, n, xlb,xub), hlb, hub, iterations=200,abs_tol=h0/n^2))
    elseif kernel == gammakernel
        xlb = 0.0
        logxdata = log.(xdata)
        return Optim.minimizer(Optim.optimize(h -> KernelEstimator.Jh(xdata, logxdata, kernel, h, w, n, xlb,xub), hlb, hub, iterations=200,abs_tol=h0/n^2))
    end
    return Optim.minimizer(Optim.optimize(h -> KernelEstimator.Jh(xdata, kernel, h, w, n, xlb,xub), hlb, hub, iterations=200,abs_tol=h0/n^2))
end

"""
    create_gamma_kde_mixture_parameters(smp::Array{Float64,1})

This is using Chen, Song Xi. "Probability density function estimation using gamma kernels." Annals of the Institute of Statistical Mathematics 52, no. 3 (2000): 471-480.
Practically, given the bandwidth bw, the kernels are Gamma(α = x/bw+1, β = 1/bw) in the shape, rate parameterisation.
Returns the list of α for all kernels and the common β.

# Examples
```julia-repl
julia> bar([1, 2], [1, 2])
1
```
"""
function create_gamma_kde_mixture_parameters_one_β(smp::Array{T,1}) where T <: Real
    bw = bwlscv(smp, gammakernel)
    # if bw == 0
    #     bw = bwlcv(smp, gammakernel)
    # end
    # bw = bwlcv(smp, gammakernel)
    if bw==0
        error("bandwidth estimation by least square cross validation failed")
    else
        return smp ./ bw .+ 1, 1/bw
    end
end

function create_gamma_kde_mixture_parameters(smp::Array{T,1}) where T <: Real
    α_list, β = create_gamma_kde_mixture_parameters_one_β(smp)

    return α_list, repeat([β], length(α_list))
end

function create_gamma_mixture_density_αβ(α_list::AbstractArray{T,1}, β_list::AbstractArray{U,1}, wms::AbstractArray{V,1}) where {T <: Real, U <: Real, V <: Real}
    #use 1/θ because of the way the Gamma distribution is parameterised in Julia Distributions.jl
    function res(x::Real)
        sum(wms[i] * pdf(Gamma(α_list[i], 1/β_list[i]), x) for i in eachindex(wms))
    end
    return res
end

"""
    create_gamma_kde_mixture_parameters(xdata::RealMatrix)

This is using Aitchison, J., & Lauder, I. J. (1985). Kernel density estimation for compositional data. Applied Statistics, 129–137.
Optimal bandwidth selection is performed using a likelihood criterion, as in the article, which may not be the best choice.
The type RealMatrix is inherited from package StatsFuns

# Examples
```julia-repl
julia> using Distributions, Random
julia> Random.seed!(0);
julia> xdata = rand(Dirichlet([0.3,5.,2.3]), 5)
3×5 Array{Float64,2}:
 0.103154  0.000797253  0.0274887  0.0004612  0.069083
 0.525477  0.48484      0.620501   0.220089   0.500526
 0.371369  0.514362     0.352011   0.77945    0.430391
 julia> create_Dirichlet_kde_mixture_parameters(xdata)
 3×5 Array{Float64,2}:
 1.02705  1.00021  1.00721  1.00012  1.01812
 1.13781  1.12716  1.16273  1.05772  1.13127
 1.0974   1.1349   1.09232  1.20442  1.11288
```
"""
function create_Dirichlet_kde_mixture_parameters(xdata::RealMatrix)
    λ = DualOptimalFiltering.bwlcv(xdata, dirichletkernel)
    # if bw == 0
    #     bw = bwlcv(smp, gammakernel)
    # end
    # bw = bwlcv(smp, gammakernel)
    if λ==0
        error("bandwidth estimation by least square cross validation failed")
    else
        return 1 .+ xdata ./ λ
    end
end

function create_Dirichlet_kde_mixture_parameters(xdata::Array{Array{Float64,1},1})
    λ = DualOptimalFiltering.bwlcv(xdata, dirichletkernel)
    # if bw == 0
    #     bw = bwlcv(smp, gammakernel)
    # end
    # bw = bwlcv(smp, gammakernel)
    if λ==0
        error("bandwidth estimation by least square cross validation failed")
    else
        return Array{Float64,1}[1 .+ x ./ λ for x in xdata]
    end
end

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
""""

 NOTE: Here we don't implement the interface for cost transformation (e.g. ∂l/∂ξ).
 Instead, the problem write should provide the transfomed cost directly.

 Feedback linearizable systems shall implement this interface and join the trait.
"""

# Holy traits for linearization style dispatch
abstract type LinearizationStyle end
# systems without a special structure (e.g. use jacobian linearization)
struct JacobianLinearization <: LinearizationStyle end
# systems that are feedback linearizable (i.e. differentially flat systems)
struct FeedbackLinearization <: LinearizationStyle end
# systems that are trivially linearized (e.g. already linear)
struct TrivialLinearization <: LinearizationStyle end

LinearizationStyle(cs) = JacobianLinearization()

"--------------------- The Feedback Linearization Interface ---------------------"

"""
    $(FUNCTIONNAME)(cs::ControlSystem)

Returns the number of ξ states (linearized states).
"""
function n_linstates end

"""
$(FUNCTIONNAME)(cs::ControlSystem)

Returns the indices of the cartesian x and y coordinates in the linearized state
vector ξ.
"""
# TODO, technically this could be computed by propagating xyindex(cs) through
# ξ_from.
function ξxyindex end

"""
    $(FUNCTIONNAME)(cs::ControlSystem)

Returns a discretized LinearSystem that descibes the dynamics of the linearized
state coordinate.

NOTE: this must not be used with the non-linear state coordinates.
"""
#TODO: maybe wrap state type for safety.
function feedbacklin end

## State transformation λ

"""
    $(FUNCTIONNAME)(cs::ControlSystem, ξ)

The state conversion map `λ(ξ)`. Transforms linear states `ξ` to nonlinear states
`x`.
"""

function x_from end

"""
    $(FUNCTIONNAME)(cs::ControlSystem, x)

Inverse state conversion map `λ⁻¹(x)`. Transforms nonlinear states `x` to linear
states `ξ`.
"""
function ξ_from end

"""
    $(FUNCTIONNAME)(cs::ControlSystem, ξ)

Returns true if the state converstion map is singular at the given *linear* system
states `ξ`.
"""
function λ_issingular end

# Cost transformation

"""
    $(FUNCTIONNAME)(cs::ControlSystem, c::PlayerCost, np::Int)

Computes the transformed cost for the linear ξ coordinates.

NOTE: This may either be an exact transformation (by applying the chain-rule to the
x-domain cost) or an approximation based on a new cost directly formulated in ξ (as
suggested in https://arxiv.org/abs/1910.00681)
"""
function transformed_cost end

"-------------------------- Optional Interface Extension --------------------------"

## Affine input transformation: u = Minv(x) * (z - m(x))

"""
    $(FUNCTIONNAME)(cs::ControlSystem, x)

Returns the inverse decoupling matrix M⁻¹(x) for input transformation.
"""
function inverse_decoupling_matrix end

"""
    $(FUNCTIONNAME)(cs::ControlSystem, x)

Returns the drift term m(x) for input transformation.
"""
function decoupling_drift_term end


text:
function updateJacobian!(JF, rShamᵢ, rSham₀, ArmijoFail, ∇ₓF, τᵢ, xᵢ; update_factors=true)
    if update_factors && rShamᵢ ≤ rSham₀ && ~ArmijoFail
        JF.age += 1
        JF.facage +=1
    else
        JF.fac = factorize(JF.fac, ∇ₓF(xᵢ) - I / τᵢ, update_factors=update_factors)
        JF.age = 0
        update_factors ? JF.facage = 0 : JF.facage += 1
    end
    return JF
end

function NewtonChordShamanskii2(F, ∇ₓF, nrm, xinit, τstop; preprint="", maxItNewton=50)
    if preprint ≠ ""
        println(preprint * "(No initial Jacobian factors fed to Newton solver)")
    end
    # Initial Jacobian
    Nx = nrm(xinit)
    NF = nrm(F(xinit))
    τ = Nx / NF
    J = ∇ₓF(xinit) - I / τ
    Jfinit = factorize(J) # construct JF
    return NewtonChordShamanskii(F, ∇ₓF, nrm, xinit, τstop, Jfinit; preprint=preprint, maxItNewton=maxItNewton)
end

function NewtonChordShamanskii2(F, ∇ₓF, nrm, xinit, τstop, Jfinit; preprint="", maxItNewton=50, update_factors=true, firstJ=false)
    if preprint ≠ ""
        println(preprint * "Solving F(x) = 0 (using Shamanskii Method)")
        preprint_end = preprint * "└─> "
        preprint = preprint * "│   "
    end

    maxItArmijo = 20 # Max iterations in Armijo line search
     # Max quasi Newton iterations
    rSham₀ = 0.5     # Shamanskii minimum reduction to keep old J

    # Initial values for while loop
    Fᵢ = F(xinit)      ;  Fᵢ₋₁ = copy(Fᵢ)
    solType = eltype(Fᵢ) # The type of the solution will be that of J (or f)
    xᵢ₋₁ = convert(Vector{solType}, xinit)
    xᵢ = copy(xᵢ₋₁)
    Nxᵢ = nrm(xᵢ)   ;  Nxᵢ₋₁ = Nxᵢ
    NFᵢ = nrm(Fᵢ)   ;  NFᵢ₋₁ = NFᵢ
    δxᵢ = copy(xᵢ₋₁)

    τᵢ = Nxᵢ / NFᵢ
    rShamᵢ = 0

    # Initialize the first Jacobian
    if firstJ
        JF = AgedJacobianFactors(Jfinit, 10, 10)
    else
        JF = initialize_AgedJacobianFactor(∇ₓF(xinit) - I / τᵢ, Jfinit, update_factors)
    end

    i = 0 # counter of quasi-Newton δxⱼs
    if preprint ≠ ""
        print(preprint)
        @printf "%3s   %8s   %8s   %7s   %7s" "iteration" "|F(x)|" "|δx|/|x|" "Jac age" "fac age"
        if solType == Complex{Float64}
            @printf "%12s" "|δix|/|ix|"
        elseif solType == Dual{Float64}
            @printf "%12s" "|δεx|/|εx|"
        elseif solType == Hyper{Float64}
            @printf "%12s" "|δhx|/|hx|"
        end
        println("")
        print(preprint)
        @printf "%5d       %8.1e   %8s   %7s   %7s\n" i NFᵢ "" "" ""
    end

    # Condition on relative tolerance for imaginary part of xᵢ
    nonrealtolx = true
    ArmijoFail = false

    # main iteration loop
    while (τᵢ < τstop || nonrealtolx) && (i < maxItNewton)
        i += 1
        if preprint ≠ ""
            print(preprint)
            @printf "%5d       " i
        end
        # fill containers with current iteration values
        xᵢ₋₁ .= xᵢ
        Fᵢ₋₁ .= Fᵢ
        Nxᵢ₋₁ = Nxᵢ
        NFᵢ₋₁ = NFᵢ

        # Update Jacobian (or rather its factors, JF)
        # Shamanski: only update J if |F(xᵢ)| / |F(xᵢ₋₁)| > rSham₀
        JF = updateJacobian!(JF, rShamᵢ, rSham₀, ArmijoFail, ∇ₓF, τᵢ, xᵢ, update_factors=update_factors)

        # Newton Step
        δxᵢ .= JF.fac \ -Fᵢ

        # Update δxᵢ, xᵢ, and Fᵢ with an Armijo line search
        δxᵢ, xᵢ, Fᵢ, nArmijo, ArmijoFail = searchLineArmijo!(δxᵢ, xᵢ, Fᵢ, F, maxItArmijo, nrm, preprint)

        Nxᵢ = nrm(xᵢ)
        NFᵢ = nrm(Fᵢ)
        τᵢ = Nxᵢ / NFᵢ
        rShamᵢ = NFᵢ / NFᵢ₋₁

        if ArmijoFail
            if JF.age > 0 # If Armijo fails and J is old, come back and update J
                xᵢ .= xᵢ₋₁
                Fᵢ .= Fᵢ₋₁
            else # else if J is fresh it's a complete failure
                preprint ≠ "" ? println("Complete Failure") : nothing
                return xᵢ
            end
        end

        Nδxᵢ = nrm(δxᵢ)
        RNδxᵢ = Nδxᵢ / Nxᵢ₋₁

        # Check that the non-real part is converging
        # Here I chose to only look at the relative step size, |δx|/|x|
        if solType == Float64
            nonrealtolx = false
        elseif solType == Complex{Float64}
            RNimδxᵢ = nrm(imag.(δxᵢ)) / nrm(imag.(xᵢ₋₁))
            nonrealtolx = (RNimδxᵢ > 1e-7)
        elseif solType == Dual{Float64}
            RNεδxᵢ = nrm(dualpart.(δxᵢ)) / nrm(dualpart.(xᵢ₋₁))
            nonrealtolx = (RNεδxᵢ > 1e-7)
        elseif solType == Hyper{Float64}
            RNhδxᵢ = nrm(ε₁part.(δxᵢ)) / nrm(ε₁part.(xᵢ₋₁)) +
                     nrm(ε₂part.(δxᵢ)) / nrm(ε₂part.(xᵢ₋₁)) +
                     nrm(ε₁ε₂part.(δxᵢ)) / nrm(ε₁ε₂part.(xᵢ₋₁))
            nonrealtolx = (RNhδxᵢ > 1e-7)
        end

        if preprint ≠ ""
            if ArmijoFail
                if JF.age > 0
                    print(preprint * "Armijo Failure, but old J")
                else
                    print(preprint * "Complete Armijo failure (new J): ")
                    @printf("|x|/|F(x)| = %.2g years\n", τᵢ/(365*24*60*60))
                end
            end
            nArmijo == 0 ? nothing : print(preprint * "    └─────> ") # alignment thing
            @printf "%8.1e   %8.1e   " NFᵢ RNδxᵢ
            print_marker(JF.age)
            print_marker(JF.facage)
            if solType == Complex{Float64}
                @printf "%8.1e" RNimδxᵢ
            elseif solType == Dual{Float64}
                @printf "%8.1e" RNεδxᵢ
            elseif solType == Hyper{Float64}
                @printf "%8.1e" RNhδxᵢ
            end
            println("")
        end
    end

    if preprint ≠ ""
        (τᵢ < τstop) && print(preprint_end * "Newton has reached max iterations, ")
        (τᵢ > τstop) && print(preprint_end * "Newton has converged, ")
        @printf("|x|/|F(x)| = %.2g years\n", τᵢ/(365*24*60*60))
    end

    return xᵢ
end


text:
# automatically generated by the FlatBuffers compiler, do not modify

NamespaceC.eval(quote

import ..NamespaceA

FlatBuffers.@with_kw mutable struct TableInC
    refer_to_a1::Union{NamespaceA.TableInFirstNS, Nothing} = nothing
    refer_to_a2::Union{NamespaceA.SecondTableInA, Nothing} = nothing
end
FlatBuffers.@ALIGN(TableInC, 1)
FlatBuffers.slot_offsets(::Type{T}) where {T<:TableInC} = [
    0x00000004, 0x00000006
]

TableInC(buf::AbstractVector{UInt8}) = FlatBuffers.read(TableInC, buf)
TableInC(io::IO) = FlatBuffers.deserialize(io, TableInC)

end)

max_stars_repo_name: rjkat/flatbuffers-julia

max_issues_repo_name: rjkat/flatbuffers-julia

max_forks_repo_name: rjkat/flatbuffers-julia


text:
using SimpleDifferentialOperators
using Test, LinearAlgebra, DualNumbers

@testset "Truncation rules when indices are explicitly given" begin
    x̄ = 1:5
    jumps = [0; -1; 0]
    @test JumpProcess(x̄, jumps).jumps == jumps
    @test JumpProcess(x̄, jumps, (:interior, :boundary)).jumps == jumps
    @test JumpProcess(x̄, jumps, (:boundary, :interior)).jumps == jumps
    @test JumpProcess(x̄, jumps, (:boundary, :boundary)).jumps == jumps

    jumps = [0; 1; 0]
    @test JumpProcess(x̄, jumps).jumps == jumps
    @test JumpProcess(x̄, jumps, (:interior, :boundary)).jumps == jumps
    @test JumpProcess(x̄, jumps, (:boundary, :interior)).jumps == jumps
    @test JumpProcess(x̄, jumps, (:boundary, :boundary)).jumps == jumps

    jumps = [0; -2; 0]
    @test JumpProcess(x̄, jumps).jumps == [0; -1; 0]
    @test JumpProcess(x̄, jumps, (:interior, :boundary)).jumps == [0; -1; 0]
    @test JumpProcess(x̄, jumps, (:boundary, :interior)).jumps == jumps
    @test JumpProcess(x̄, jumps, (:boundary, :boundary)).jumps == jumps

    jumps = [0; 2; 0]
    @test JumpProcess(x̄, jumps).jumps == [0; 1; 0]
    @test JumpProcess(x̄, jumps, (:interior, :boundary)).jumps == jumps
    @test JumpProcess(x̄, jumps, (:boundary, :interior)).jumps == [0; 1; 0]
    @test JumpProcess(x̄, jumps, (:boundary, :boundary)).jumps == jumps

    jumps = [-1; -2; 0]
    @test JumpProcess(x̄, jumps).jumps == [0; -1; 0]
    @test JumpProcess(x̄, jumps, (:interior, :boundary)).jumps == [0; -1; 0]
    @test JumpProcess(x̄, jumps, (:boundary, :interior)).jumps == jumps
    @test JumpProcess(x̄, jumps, (:boundary, :boundary)).jumps == jumps

    jumps = [0; 2; 1]
    @test JumpProcess(x̄, jumps).jumps == [0; 1; 0]
    @test JumpProcess(x̄, jumps, (:interior, :boundary)).jumps == jumps
    @test JumpProcess(x̄, jumps, (:boundary, :interior)).jumps == [0; 1; 0]
    @test JumpProcess(x̄, jumps, (:boundary, :boundary)).jumps == jumps

    jumps = [-1; -2; 1]
    @test JumpProcess(x̄, jumps).jumps == [0; -1; 0]
    @test JumpProcess(x̄, jumps, (:interior, :boundary)).jumps == [0; -1; 1]
    @test JumpProcess(x̄, jumps, (:boundary, :interior)).jumps == [-1; -2; 0]
    @test JumpProcess(x̄, jumps, (:boundary, :boundary)).jumps == jumps

    jumps = [-1; 2; 1]
    @test JumpProcess(x̄, jumps).jumps == [0; 1; 0]
    @test JumpProcess(x̄, jumps, (:interior, :boundary)).jumps == [0; 2; 1]
    @test JumpProcess(x̄, jumps, (:boundary, :interior)).jumps == [-1; 1; 0]
    @test JumpProcess(x̄, jumps, (:boundary, :boundary)).jumps == jumps
end

@testset "Truncation rules when jump functions are provided" begin
    x̄ = 1:5
    jumps = [0; -1; 0]
    jumpf = x -> (x == 3) ? -1 : 0
    @test JumpProcess(x̄, jumpf).jumps == jumps
    @test JumpProcess(x̄, jumpf, (:interior, :boundary)).jumps == jumps
    @test JumpProcess(x̄, jumpf, (:boundary, :interior)).jumps == jumps
    @test JumpProcess(x̄, jumpf, (:boundary, :boundary)).jumps == jumps

    jumps = [0; 1; 0]
    jumpf = x -> (x == 3) ? 1 : 0
    @test JumpProcess(x̄, jumpf).jumps == jumps
    @test JumpProcess(x̄, jumpf, (:interior, :boundary)).jumps == jumps
    @test JumpProcess(x̄, jumpf, (:boundary, :interior)).jumps == jumps
    @test JumpProcess(x̄, jumpf, (:boundary, :boundary)).jumps == jumps

    jumps = [0; -2; 0]
    jumpf = x -> (x == 3) ? -2 : 0
    @test JumpProcess(x̄, jumpf).jumps == [0; -1; 0]
    @test JumpProcess(x̄, jumpf, (:interior, :boundary)).jumps == [0; -1; 0]
    @test JumpProcess(x̄, jumpf, (:boundary, :interior)).jumps == jumps
    @test JumpProcess(x̄, jumpf, (:boundary, :boundary)).jumps == jumps

    jumps = [0; 2; 0]
    jumpf = x -> (x == 3) ? 2 : 0
    @test JumpProcess(x̄, jumpf).jumps == [0; 1; 0]
    @test JumpProcess(x̄, jumpf, (:interior, :boundary)).jumps == jumps
    @test JumpProcess(x̄, jumpf, (:boundary, :interior)).jumps == [0; 1; 0]
    @test JumpProcess(x̄, jumpf, (:boundary, :boundary)).jumps == jumps

    jumps = [-1; -2; 0]
    jumpf = x -> (x <= 3) ? -(x-1) : 0
    @test JumpProcess(x̄, jumpf).jumps == [0; -1; 0]
    @test JumpProcess(x̄, jumpf, (:interior, :boundary)).jumps == [0; -1; 0]
    @test JumpProcess(x̄, jumpf, (:boundary, :interior)).jumps == jumps
    @test JumpProcess(x̄, jumpf, (:boundary, :boundary)).jumps == jumps

    jumps = [0; 2; 1]
    jumpf = x -> (x >= 3) ? (5-x) : 0
    @test JumpProcess(x̄, jumpf).jumps == [0; 1; 0]
    @test JumpProcess(x̄, jumpf, (:interior, :boundary)).jumps == jumps
    @test JumpProcess(x̄, jumpf, (:boundary, :interior)).jumps == [0; 1; 0]
    @test JumpProcess(x̄, jumpf, (:boundary, :boundary)).jumps == jumps

    jumps = [-1; -2; 1]
    jumpf = x -> (x <= 3) ? -(x-1) : 1
    @test JumpProcess(x̄, jumpf).jumps == [0; -1; 0]
    @test JumpProcess(x̄, jumpf, (:interior, :boundary)).jumps == [0; -1; 1]
    @test JumpProcess(x̄, jumpf, (:boundary, :interior)).jumps == [-1; -2; 0]
    @test JumpProcess(x̄, jumpf, (:boundary, :boundary)).jumps == jumps

    jumps = [-1; 2; 1]
    jumpf = x -> (x >= 3) ? (5-x) : -1
    @test JumpProcess(x̄, jumpf).jumps == [0; 1; 0]
    @test JumpProcess(x̄, jumpf, (:interior, :boundary)).jumps == [0; 2; 1]
    @test JumpProcess(x̄, jumpf, (:boundary, :interior)).jumps == [-1; 1; 0]
    @test JumpProcess(x̄, jumpf, (:boundary, :boundary)).jumps == jumps

    # on irregular grids
    x̄ = [1; 2.5; 3; 4.9; 5]
    jumps = [-1; 2; 1]
    jumpf = x -> (x >= 3) ? (5-x) : -1
    @test JumpProcess(x̄, jumpf).jumps == [0; 1; 0]
    @test JumpProcess(x̄, jumpf, (:interior, :boundary)).jumps == [0; 2; 1]
    @test JumpProcess(x̄, jumpf, (:boundary, :interior)).jumps == [-1; 1; 0]
    @test JumpProcess(x̄, jumpf, (:boundary, :boundary)).jumps == jumps
end

@testset "Accuracy test for jump diffusion extension operators" begin
    x̄ = [1; 2.3; 3.4; 4.5; 5.0]
    jumps = [0; -1; 0]
    method = JumpProcess(x̄, jumps)
    L̄ = ExtensionDifferentialOperator(x̄, method)
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 1. -1. 0. 0.; 0. 0. 0. 0. 0.]

    jumps = [0; 1; 0]
    method = JumpProcess(x̄, jumps)
    L̄ = ExtensionDifferentialOperator(x̄, method)
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 0. -1. 1. 0.; 0. 0. 0. 0. 0.]

    jumps = [0; -2; 0]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 1. -1. 0. 0.; 0. 0. 0. 0. 0.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:interior, :boundary)))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 1. -1. 0. 0.; 0. 0. 0. 0. 0.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:boundary, :interior)))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 1. 0. -1. 0. 0.; 0. 0. 0. 0. 0.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:boundary, :boundary)))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 1. 0. -1. 0. 0.; 0. 0. 0. 0. 0.]

    jumps = [0; 2; 0]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 0. -1. 1. 0.; 0. 0. 0. 0. 0.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:interior, :boundary)))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 0. -1. 0. 1.; 0. 0. 0. 0. 0.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:boundary, :interior)))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 0. -1. 1. 0.; 0. 0. 0. 0. 0.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:boundary, :boundary)))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 0. -1. 0. 1.; 0. 0. 0. 0. 0.]

    jumps = [-1; -2; 0]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 1. -1. 0. 0.; 0. 0. 0. 0. 0.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:interior, :boundary)))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 1. -1. 0. 0.; 0. 0. 0. 0. 0.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:boundary, :interior)))
    @test Array(L̄) == [1. -1. 0. 0. 0.; 1. 0. -1. 0. 0.; 0. 0. 0. 0. 0.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:boundary, :boundary)))
    @test Array(L̄) == [1. -1. 0. 0. 0.; 1. 0. -1. 0. 0.; 0. 0. 0. 0. 0.]

    jumps = [0; 2; 1]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 0. -1. 1. 0.; 0. 0. 0. 0. 0.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:interior, :boundary)))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 0. -1. 0. 1.; 0. 0. 0. -1. 1.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:boundary, :interior)))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 0. -1. 1. 0.; 0. 0. 0. 0. 0.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:boundary, :boundary)))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 0. -1. 0. 1.; 0. 0. 0. -1. 1.]

    jumps = [-1; -2; 1]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 1. -1. 0. 0.; 0. 0. 0. 0. 0.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:interior, :boundary)))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 1. -1. 0. 0.; 0. 0. 0. -1. 1.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:boundary, :interior)))
    @test Array(L̄) == [1. -1. 0. 0. 0.; 1. 0. -1. 0. 0.; 0. 0. 0. 0. 0.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:boundary, :boundary)))
    @test Array(L̄) == [1. -1. 0. 0. 0.; 1. 0. -1. 0. 0.; 0. 0. 0. -1. 1.]

    jumps = [-1; 2; 1]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 0. -1. 1. 0.; 0. 0. 0. 0. 0.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:interior, :boundary)))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 0. -1. 0. 1.; 0. 0. 0. -1. 1.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:boundary, :interior)))
    @test Array(L̄) == [1. -1. 0. 0. 0.; 0. 0. -1. 1. 0.; 0. 0. 0. 0. 0.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:boundary, :boundary)))
    @test Array(L̄) == [1. -1. 0. 0. 0.; 0. 0. -1. 0. 1.; 0. 0. 0. -1. 1.]

    jumps = [-1; -1; -1]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 1. -1. 0. 0.; 0. 0. 1. -1. 0.]
    @test L̄ == ExtensionDifferentialOperator(x̄, JumpProcess(x̄, -1))
    @test L̄ == ExtensionDifferentialOperator(x̄, JumpProcess(x̄, -1.0))
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:interior, :boundary)))
    @test Array(L̄) == [0. 0. 0. 0. 0.; 0. 1. -1. 0. 0.; 0. 0. 1. -1. 0.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:boundary, :interior)))
    @test Array(L̄) == [1. -1. 0. 0. 0.; 0. 1. -1. 0. 0.; 0. 0. 1. -1. 0.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:boundary, :boundary)))
    @test Array(L̄) == [1. -1. 0. 0. 0.; 0. 1. -1. 0. 0.; 0. 0. 1. -1. 0.]

    jumps = [1; 1; 1]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps))
    @test L̄ == ExtensionDifferentialOperator(x̄, JumpProcess(x̄, 1))
    @test L̄ == ExtensionDifferentialOperator(x̄, JumpProcess(x̄, 1.0))
    @test Array(L̄) == [0. -1. 1. 0. 0.; 0. 0. -1. 1. 0.; 0. 0. 0. 0. 0.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:interior, :boundary)))
    @test Array(L̄) == [0. -1. 1. 0. 0.; 0. 0. -1. 1. 0.; 0. 0. 0. -1. 1.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:boundary, :interior)))
    @test Array(L̄) == [0. -1. 1. 0. 0.; 0. 0. -1. 1. 0.; 0. 0. 0. 0. 0.]
    L̄ = ExtensionDifferentialOperator(x̄, JumpProcess(x̄, jumps, (:boundary, :boundary)))
    @test Array(L̄) == [0. -1. 1. 0. 0.; 0. 0. -1. 1. 0.; 0. 0. 0. -1. 1.]
end

@testset "Accuracy test for jump diffusion operators on the interior" begin
    x̄ = [1; 2.3; 3.4; 4.5; 5.0]
    jumps = [0; -1; 0]
    bc = (Absorbing(), Absorbing())
    method = JumpProcess(x̄, jumps)
    L = Lₙbc(x̄,bc,method)
    @test Array(L) == [0. 0. 0.; 1. -1. 0.; 0. 0. 0.]

    jumps = [0; 1; 0]
    method = JumpProcess(x̄, jumps)
    L = Lₙbc(x̄,bc,method)
    @test Array(L) == [0. 0. 0.; 0. -1. 1.; 0. 0. 0.]

    jumps = [0; -1; 0]
    L = Lₙbc(x̄,bc,JumpProcess(x̄, jumps))
    @test Array(L) == [0. 0. 0.; 1. -1. 0.; 0. 0. 0. ]
    L = Lₙbc(x̄,bc,JumpProcess(x̄, jumps, (:interior, :boundary)))
    @test Array(L) == [0. 0. 0.; 1. -1. 0.; 0. 0. 0. ]
    L = Lₙbc(x̄,bc,JumpProcess(x̄, jumps, (:boundary, :interior)))
    @test Array(L) == [0. 0. 0.; 1. -1. 0.; 0. 0. 0. ]
    L = Lₙbc(x̄,bc,JumpProcess(x̄, jumps, (:boundary, :boundary)))
    @test Array(L) == [0. 0. 0.; 1. -1. 0.; 0. 0. 0. ]

    jumps = [0; -2; 0]
    L = Lₙbc(x̄,bc,JumpProcess(x̄, jumps))
    @test Array(L) == [0. 0. 0.; 1. -1. 0.; 0. 0. 0. ]
    L = Lₙbc(x̄,bc,JumpProcess(x̄, jumps, (:interior, :boundary)))
    @test Array(L) == [0. 0. 0.; 1. -1. 0.; 0. 0. 0. ]
    L = Lₙbc(x̄,bc,JumpProcess(x̄, jumps, (:boundary, :interior)))
    @test Array(L) == [0. 0. 0.; 0. -1. 0.; 0. 0. 0. ]
    L = Lₙbc(x̄,bc,JumpProcess(x̄, jumps, (:boundary, :boundary)))
    @test Array(L) == [0. 0. 0.; 0. -1. 0.; 0. 0. 0. ]

    bc = (Reflecting(), Reflecting())
    jumps = [0; -2; 0]
    L = Lₙbc(x̄,bc,JumpProcess(x̄, jumps))
    @test Array(L) == [0. 0. 0.; 1. -1. 0.; 0. 0. 0. ]
    L = Lₙbc(x̄,bc,JumpProcess(x̄, jumps, (:interior, :boundary)))
    @test Array(L) == [0. 0. 0.; 1. -1. 0.; 0. 0. 0. ]
    @test_throws ErrorException Lₙbc(x̄,bc,JumpProcess(x̄, jumps, (:boundary, :interior)))
    @test_throws ErrorException Lₙbc(x̄,bc,JumpProcess(x̄, jumps, (:boundary, :boundary)))
end


text:
module Interp1d

"""
    interp(xi::Vector{T}, yi::VecOrMat{T}, ::Type{Nearest}; copy=true, sorted=false) where T
    interp(xi::Vector{T}, yi::VecOrMat{T}, ::Type{Previous}; copy=true, sorted=false) where T

Returns a nearest interpolation function `f(x)` based on interpolated points `xi` and `yi`.

If `xi` has a duplicated value, `ArgumentError` is thrown.

If `copy` is true, this function makes internal copies of `xi` and `yi`, otherwise it uses references.

If `sorted` is false, this function sorts `xi` first, otherwise it assumes `xi` is monotonically increasing.

# Examples
This is a simple interpolation example:
```jldoctest
julia> using Interp1d

julia> xi = [-1.0, 0.0, 3.0/2.0];

julia> yi = [2.0, 1.0, 3.0];

julia> f = interp(xi, yi, Nearest);

julia> f(1.0)
3.0
```

`yi` can be matrix for multi univerate interpolations.
The number of rows of `yi` must be equal to the length of `xi`.

```jldoctest
julia> using Interp1d

julia> xi = [-1.0, 0.0, 3.0/2.0];

julia> yi = [ 2.0  1.0 3.0;
             -1.0 -4.4 2.0];

julia> f = interp(xi, yi, Nearest);

julia> f(0.0)
2-element Vector{Float64}:
  1.0
 -4.4
```

"""
function interp end

# Interp mode
abstract type InterpMode end

function _prepare_input(xi, yi, copy, sorted)

    # check duplicate in xi
    length(xi) != length(unique(xi)) && throw(ArgumentError("xi has duplicated value: $xi"))

    _xi, _yi = copy ? (deepcopy(xi), deepcopy(yi)) : (xi, yi)

    if !sorted
        idxs = sortperm(_xi)
        _xi = view(_xi, idxs)
        _yi = ndims(_yi) == 1 ? view(_yi, idxs) : view(_yi, :, idxs)
    end

    return _xi, _yi
end

include("Nearest.jl")
include("Previous.jl")

# Supported Interpolation mode list
const INTERP_MODE_LIST = [Nearest, Previous]

# export all public APIs
const _EXCLUDE_SYMBOLS = [Symbol(@__MODULE__), :eval, :include]
for sym in names(@__MODULE__, all = true)
    sym_string = string(sym)
    if sym in _EXCLUDE_SYMBOLS ||
       !Base.isidentifier(sym) ||
       startswith(sym_string, "_") ||
       startswith(sym_string, "@_")
        continue
    end
    @eval export $sym
end

end


text:
using Omega
using OmegaTestModels
import ForwardDiff
using Flux

const ΩTs = [LinearΩ{Vector{Int}, UnitRange{Int64}, Vector{Float64}},
             Omega.SimpleΩ{Vector{Int}, Float64},
             LinearΩ{Vector{Int}, UnitRange{Int64}, TrackedArray{Float64, 1, Array{Float64,1}}}]

const gradalgs = [Omega.TrackerGrad,
                  Omega.ForwardDiffGrad,
                  # Omega.ZygoteGrad
                  ]

function testgrad(OT = Omega.SimpleΩ{Vector{Int}, Float64})
  μ = uniform(0.0, 1.0)
  x = normal(μ, 1.0) + normal(μ, 1.0)
  y = logerr(x ==ₛ 1.0)
  ω = OT()
  y(ω)
  Omega.lineargradient(y, ω, Omega.ForwardDiffGrad)
end

testgrad(LinearΩ{Vector{Int}, UnitRange{Int64}, Vector{Float64}})

# Test them all
istracked(::Type{LinearΩ{I, SEG, T}}) where {I, SEG, T <: TrackedArray} = true
istracked(ΩT) = false

function testgrad2(; modelcond, ΩT, gradalg)
  @show ΩT
  println()
  Omega.lineargradient(logerr(modelcond), ΩT(), gradalg)
end

# Filter to differentiable models
models = filter(hascond ∧ isdiff, allmodels)

# Dont use TrackerGrad on wrong Omegas, etc
f1(m, ΩT, gradalg) = gradalg == Omega.TrackerGrad ?  istracked(ΩT) : true
f2(m, ΩT, gradalg) = gradalg == Omega.ForwardDiffGrad ? !istracked(ΩT) : true

function runtests()
  for model in models, gradalg in gradalgs, ΩT in ΩTs
    if f1(model, ΩT, gradalg) && f2(model, ΩT, gradalg)
      @show model.name
      @show ΩT
      @show gradalg
      testgrad2(; modelcond = model.cond, ΩT = ΩT, gradalg = gradalg)
    end
  end
end

runtests()

function zygotetest()
  x = normal(0, 1, (10,))
  y = sum(x)
  lineargradient(y, rand(defΩ()), ZygoteGrad)
end


text:
module TSP

import Random
using ProgressMeter
import Base
using ..Storage

# Util function to generate random instances
function generateRandomTSPData(numCities, dims=2)
    Random.rand(numCities, dims)
end

function readTourFile(path)
    file = open(path)
    lines = readlines(file)
    lines = filter(x -> isdigit(x[1]), lines)
    data = [parse(Float32, split(x)[j]) for x=lines, j=2:3]
    return data
end

# The state of a solution to the problem
struct TSPSolution{N}
    order::SizedArray{Int32, N} # I doubt people will have more than 2**31 cities
    v::Int
end

# The struct containing the problem data
# This is constant and should NEVER be changed during optimization
# TODO We should probably enforce that
struct TSPData{ND}
    coords::SizedArray{Float32, ND}
end

function generateProblemForData(problemData::AbstractArray{T, 2}) where {T}
    numCities, dims = size(problemData)

    # Make julia create a compiled version of the code for these specific problem sizes
    # Because they are passed as types and not values it tells julia to compile
    # a different version for every pair of these values
    # It makes the loss calculation stupidly fast for example
    function superFastGen(::Val{dims}, ::Val{numCities}) where {dims, numCities}

        randomInit(state, _) = state.order .= Random.randperm(numCities)

        @Base.pure function rawDistance(a, b, coords)
            result = 0.0
            for j in 1:dims
                @inbounds result += (coords[a, j] - coords[b, j]) ^2
            end
            result
        end

        function nnInit(state, problemData)
            coords = problemData.coords
            order = state.order
            doneCities = zeros(Bool, numCities)
            startIx = Random.rand(1:numCities)
            order[1] = startIx
            doneCities[startIx] = true
            prev = startIx

            @showprogress "Init state" for i in 2:numCities
                bestOther = 0
                bestDistance = Inf32
                for j in 1:numCities
                    d = rawDistance(prev, j, coords)
                    if !doneCities[j] && d < bestDistance
                        bestOther = j
                        bestDistance = d
                    end
                end
                doneCities[i] = true
                order[i] = bestOther
                prev = bestOther
            end
        end

        # We assume proper indices !!!!
        @Base.pure @inline function distance(a, b, order, coords)
            result = 0.0
            for j in 1:dims
                @inbounds result += (coords[order[a], j] - coords[order[b], j]) ^2
            end
            result
        end

        @Base.pure @inline checkBounds(n) = mod1(n, numCities)

        @inline @Base.pure function loss(state, problemData)
            order = state.order
            coords = problemData.coords
            result = distance(1, numCities, order, coords) # The first is the link from the end to start
            for i in 2:numCities
                result += distance(i , i - 1, order, coords)
            end
            result
        end

        ### Swap2
        @inline function swap2Propose(currentState, problemData, dir)
                a, b = dir[1], dir[2]
                if b < a
                    a, b = b, a
                end
                @inline @Base.pure d(a, b) = distance(checkBounds(a), checkBounds(b),
                                                      currentState.order, problemData.coords)
                if (a == b) return 0.0 end #swapping the same cities
                return (
                        -(d(a - 1, a) + d(a, a + 1) + d(b - 1, b) + d(b, b + 1)) # old distances
                        +(d(a - 1, b) + d(b, a + 1) + d(b - 1, a) + d(a, b + 1)) # new distances
                       )
        end
        @inline function swap2Apply(currentState, problemData, dir)
            @inbounds a, b = dir[1], dir[2]
            if b < a; a, b = b, a end
            if (a == b) return end
            order = currentState.order
            @inbounds order[a], order[b] = order[b], order[a]
        end
        swap2 = (propose = swap2Propose, apply = swap2Apply)

        @inline function opt2Propose(currentState, problemData, dir)
                @inbounds a, b = dir[1], dir[2]
                if b < a; a, b = b, a end
                @inline @Base.pure d(a, b) = distance(checkBounds(a), checkBounds(b),
                                                      currentState.order, problemData.coords)
                if (a == b) || (a == 1 && b == numCities) return 0.0 end
                (d(a - 1, b) + d(a, b + 1) - d(a - 1, a) - d(b, b + 1))
        end
        @inline function opt2Apply(currentState, problemData, dir)
            @inbounds a, b = dir[1], dir[2]
            if b < a; a, b = b, a end
            if (a == b) || (a == 1 && b == numCities) return end
            order = currentState.order
            for i in 0:div((b - a + 1), 2) - 1
                @inbounds order[a + i], order[b - i] = order[b - i], order[a + i]
            end
        end

        opt2 = (propose = opt2Propose, apply = opt2Apply)

        (
         # Problem Data
         dataType=TSPData{(numCities, dims)}, # How to represent the data
         data=(problemData, ), # What is the actual data we are solving for
         # Solution
         solType=TSPSolution{numCities}, # How to represent a solution
         init=randomInit, # How to init the solution
         # Loss
         loss=loss,
         # Exploration
         neighborSpace=CartesianIndices((numCities, numCities)), # What space to explore
         opt2...
        )
    end
    superFastGen(Val(dims), Val(numCities))
end

export generateRandomTSPData, generateProblemForData

end


text:
"""
`CropRatio <: ImageTransformation`

Description
============

Crops out the biggest area around the center of the given image
such that said sub-image satisfies the specified aspect ratio
(i.e. width divided by height).

Usage
======

    CropRatio()

    CropRatio(ratio = 1)

Arguments
==========

- **`ratio`** : The ratio of image height to image width that
the output image should satisfy

Methods
========

- **`transform`** : Applies the transformation to the given Image
or set of images

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

Examples
========

    using Augmentor
    using TestImages

    # load an example image
    img = testimage("lena")

    # Crop to image to a 2:1 aspect ratio
    img_new = transform(CropRatio(2), img)

see also
=========

`ImageTransformation`, `transform`
"""
immutable CropRatio <: ImageTransformation
    ratio::Float64

    function CropRatio(ratio::Real)
        ratio > 0 || throw(ArgumentError("ratio has to be greater than 0"))
        new(Float64(ratio))
    end
end

CropRatio(; ratio = 1.) = CropRatio(ratio)

Base.showcompact(io::IO, op::CropRatio) = print(io, "Crop to $(op.ratio) aspect ratio.")
multiplier(::CropRatio) = 1

function transform{T<:AbstractImage}(op::CropRatio, img::T)
    w = width(img)
    h = height(img)

    nw = floor(Int, h * op.ratio)
    nh = floor(Int, w / op.ratio)
    nw = nw > 1 ? nw : 1
    nh = nh > 1 ? nh : 1

    result = if nw == w || nh == h
        img
    elseif nw < w
        i = 1 + floor(Int, (w - nw) / 2)
        @assert i > 0
        crop(img, i:(i+nw-1), 1:h)
    elseif nh < h
        i = 1 + floor(Int, (h - nh) / 2)
        @assert i > 0
        crop(img, 1:w, i:(i+nh-1))
    end::T
    _log_operation!(result, op)
end


text:
using Pumas, MCMCChains, StatsPlots, Random
plotly()

# Define Population
population = Population(
    map(i -> Subject(
        id = i),
        1:8))

# Build model
tumour_growth_model = @model begin
    @param   begin
        # mean_initial_volume ~ Constrained(
        #     Normal(0.2, 1), lower=0.001, upper=Inf, init=0.2)
        # mean_critical_volume ~ Constrained(
        #     Normal(0.2, 1), lower=0.001, upper=Inf, init=0.2)
        # mean_growth_rate ~ Constrained(
        #     Normal(0.2, 1), lower=0.001, upper=Inf, init=0.2)
        # std ~ Constrained(
        #     MvNormal([1 , 1, 1]),
        #     lower=[0, 0, 0],
        #     upper=[Inf, Inf, Inf],
        #     init=[0.09,0.09, 0.09])
        # σ_base ~ Constrained(
        #     Normal(0.2, 1), lower=0.001, upper=Inf, init=0.2)
        # σ_rel ~ Constrained(
        #     Normal(0.2, 1), lower=0.001, upper=Inf, init=0.2)
        mean_initial_volume ∈ RealDomain(lower=0, init = 0.2)
        mean_critical_volume ∈ RealDomain(lower=0, init = 1)
        mean_growth_rate ∈ RealDomain(lower = 0, init= 0.1)
        std ∈ PDiagDomain(init=[0.09,0.09, 0.09])
        σ_base ∈ RealDomain(lower=0, init=0.04)
        σ_rel ∈ RealDomain(lower=0, init=0.04)
    end

    @random begin
        η ~ MvNormal(std)
    end

    @pre begin
        initial_volume = mean_initial_volume * exp(η[1])
        critical_volume = mean_critical_volume * exp(η[2])
        growth_rate = mean_growth_rate * exp(η[3])
    end

    @init begin
        volume = initial_volume
    end

    @dynamics begin
       volume' = (growth_rate * critical_volume * volume) / (volume + critical_volume)
    end

  @derived begin
      σ_tot := σ_base .+ volume * σ_rel
      Measurement ~ @. Normal(volume, σ_tot)
    end
end

# Simulate Population
param = init_param(tumour_growth_model)
obs = simobs(tumour_growth_model, population, param, obstimes=0:1:30)
scatter(obs)

# Convert to dataset
simdf = DataFrame(obs)
data = read_pumas(
    simdf, observations=[:Measurement], event_data=false)

# Infer parameters
param = init_param(tumour_growth_model)
result = fit(
    tumour_growth_model, data, param, Pumas.BayesMCMC();
    nsamples=2000, nadapts=1000)

# Visualise population parameters
chains = Pumas.Chains(result)
plot(chains[1000:end])

tumour_growth_model.prob

max_stars_repo_name: DavAug/erlotinib

max_issues_repo_name: DavAug/erlotinib

max_forks_repo_name: DavAug/erlotinib


text:
using MyExample
using Test

@testset "MyExample.jl" begin
    @test my_f(2, 1) == 7
    @test my_f(2, 3) == 13
    @test new_f(1, 2) == 0.5
end

max_stars_repo_name: Mikejmnez/MyExample

max_issues_repo_name: Mikejmnez/MyExample

max_forks_repo_name: Mikejmnez/MyExample


text:
type SnapShot
     isz:: Int64
     isx:: Int64
     nz :: Int64
     nx :: Int64
     ext:: Int64
     iflag :: Int64
     dt :: Float64
     it :: Int64
     vz :: Array{Float64, 1}
     vx :: Array{Float64, 1}
     pz :: Array{Float64, 1}
     px :: Array{Float64, 1}
end

function InitSnapShot(isz::Int64, isx::Int64, nz::Int64, nx::Int64, ext::Int64, iflag::Int64, dt::Float64, it::Int64)
    if iflag == 1
       Nz = nz + 2*ext
       Nx = nx + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
       Nx = nx + 2*ext
    end
    spt = SnapShot(isz, isx, nz, nx, ext, iflag, dt, it, zeros(Nz*Nx), zeros(Nz*Nx), zeros(Nz*Nx), zeros(Nz*Nx))
    return spt
end

import Base.copy
function copy(snapShot1::SnapShot, snapShot2::SnapShot)
    if snapShot2.nz != snapShot1.nz || snapShot2.nx != snapShot1.nx || snapShot2.ext != snapShot1.ext || snapShot2.dt != snapShot1.dt || snapShot2.iflag != snapShot1.iflag
       error("the two snapShot are different")
    end
    snapShot1.isz = snapShot2.isz
    snapShot1.isx = snapShot2.isx
    snapShot1.it  = snapShot2.it
    snapShot1.vz[:] = snapShot2.vz[:]
    snapShot1.vx[:] = snapShot2.vx[:]
    snapShot1.pz[:] = snapShot2.pz[:]
    snapShot1.px[:] = snapShot2.px[:]
    return nothing
end

import Base.norm
function norm(snapShot::SnapShot)
    L = norm(snapShot.vz)
    L = L + norm(snapShot.vx)
    L = L + norm(snapShot.vz)
    L = L + norm(snapShot.pz)
    L = L + norm(snapShot.px)
    return L
end

import Base.-
function -(snapShot1::SnapShot, snapShot2::SnapShot)
    isz= snapShot1.isz
    isx= snapShot1.isx
    nz = snapShot1.nz
    nx = snapShot1.nx
    ext = snapShot1.ext
    iflag = snapShot1.iflag
    dt = snapShot1.dt
    it = snapShot1.it
    if nz != snapShot2.nz || nx != snapShot2.nx || ext != snapShot2.ext || iflag != snapShot2.iflag
       error("size dismatch")
    end
    snapShot3 = InitSnapShot(isz, isx, nz, nx, ext, iflag, dt, it)
    snapShot3.vz = snapShot1.vz - snapShot2.vz
    snapShot3.vx = snapShot1.vx - snapShot2.vx
    snapShot3.pz = snapShot1.pz - snapShot2.pz
    snapShot3.px = snapShot1.px - snapShot2.px
    return snapShot3
end

import Base.+
function +(snapShot1::SnapShot, snapShot2::SnapShot)
    isz= snapShot1.isz
    isx= snapShot1.isx
    nz = snapShot1.nz
    nx = snapShot1.nx
    ext = snapShot1.ext
    iflag = snapShot1.iflag
    dt = snapShot1.dt
    it = snapShot1.it
    if nz != snapShot2.nz || nx != snapShot2.nx || ext != snapShot2.ext || iflag != snapShot2.iflag
       error("size dismatch")
    end
    snapShot3 = InitSnapShot(isz, isx, nz, nx, ext, iflag, dt, it)
    snapShot3.vz = snapShot1.vz + snapShot2.vz
    snapShot3.vx = snapShot1.vx + snapShot2.vx
    snapShot3.pz = snapShot1.pz + snapShot2.pz
    snapShot3.px = snapShot1.px + snapShot2.px
    return snapShot3
end

import Base.write
function write(path::ASCIIString, spt::SnapShot, ot::Float64, dt::Float64)
    fid = open(path, "w")
    write(fid, spt.isz); write(fid, spt.isx  );
    write(fid, spt.nz ); write(fid, spt.nx   );
    write(fid, spt.ext); write(fid, spt.iflag);
    write(fid, ot     ); write(fid, dt       );
    write(fid, spt.pz+spt.px)
    flush(fid)
    return fid
end

function write(fid::IOStream, spt::SnapShot)
    write(fid, vec(spt.pz+spt.px))
    return nothing
end

import Base.read
function read(path::ASCIIString, it::Int64)
    fid = open(path, "r") ; seek(fid, sizeof(Int64)*2);
    nz  = read(fid, Int64); nx = read(fid, Int64);
    ext = read(fid, Int64); iflag = read(fid, Int64);
    if iflag == 1
       Nz = nz + 2*ext
       Nx = nx + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
       Nx = nx + 2*ext
    end
    position = sizeof(Int64)*8 + (it-1)*sizeof(Float64)*Nz*Nx
    seek(fid, position)
    d = read(fid, Float64, Nz*Nx)
    d = reshape(d, Nz, Nx)
    close(fid)
    return d
end

function readSnapShot(path::ASCIIString, it::Int64)
    fid = open(path, "r") ;
    isz = read(fid, Int64); isx   = read(fid, Int64)
    nz  = read(fid, Int64); nx    = read(fid, Int64);
    ext = read(fid, Int64); iflag = read(fid, Int64);
    ot= read(fid, Float64); dt  = read(fid, Float64);
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    position = sizeof(Int64)*8 + (it-1)*sizeof(Float64)*Nz*Nx*4
    seek(fid, position)
    vz = read(fid, Float64, Nz*Nx)
    vx = read(fid, Float64, Nz*Nx)
    pz = read(fid, Float64, Nz*Nx)
    px = read(fid, Float64, Nz*Nx)
    spt = SnapShot(isz, isx, nz, nx, ext, iflag, dt, it, vz, vx, pz, px)
    close(fid)
    return spt
end

function InfoWfd(path::ASCIIString; print_flag = false)
    fid = open(path, "r")
    isz = read(fid, Int64); isx = read(fid, Int64);
    nz  = read(fid, Int64); nx  = read(fid, Int64);
    ext = read(fid, Int64); iflag=read(fid, Int64);
    ot  = read(fid, Float64); dt =read(fid, Float64);
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    nt = round(Int64, (filesize(fid) - sizeof(Int64)*8) / (Nz*Nx*sizeof(Float64)))
    close(fid)
    if print_flag
       println("shot z: $isz, shot x: $isx")
       println("zlength: $nz, xlength: $nx")
       println("padding: $ext, surface: $iflag")
       println("ot: $ot, dt: $dt, nt: $nt")
    end
    return isz, isx, nz, nx, ext, iflag, ot, dt, nt
end

function InfoFullWfd(path::ASCIIString; print_flag = false)
    fid = open(path, "r")
    isz = read(fid, Int64); isx = read(fid, Int64);
    nz  = read(fid, Int64); nx  = read(fid, Int64);
    ext = read(fid, Int64); iflag=read(fid, Int64);
    ot  = read(fid, Float64); dt =read(fid, Float64);
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    nt = round(Int64, (filesize(fid) - sizeof(Int64)*8) / (Nz*Nx*4*sizeof(Float64)))
    close(fid)
    if print_flag
       println("shot z: $isz, shot x: $isx")
       println("zlength: $nz, xlength: $nx")
       println("padding: $ext, surface: $iflag")
       println("ot: $ot, dt: $dt, nt: $nt")
    end
    return isz, isx, nz, nx, ext, iflag, ot, dt, nt
end

function numberOfSnapShots(path::ASCIIString)
    fid = open(path, "r")
    isz = read(fid, Int64); isx = read(fid, Int64);
    nz  = read(fid, Int64); nx  = read(fid, Int64);
    ext = read(fid, Int64); iflag=read(fid, Int64);
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    nt = round(Int64, (filesize(fid) - sizeof(Int64)*8) / (Nz*Nx*4*sizeof(Float64)))
    close(fid)
    return nt
end

function writeSnapShot(path::ASCIIString, spt::SnapShot)
    fid = open(path, "w")
    write(fid, spt.isz); write(fid, spt.isx  );
    write(fid, spt.nz ); write(fid, spt.nx   );
    write(fid, spt.ext); write(fid, spt.iflag);
    write(fid, 0.0    ); write(fid, spt.dt   );
    write(fid, spt.vz ); write(fid, spt.vx   );
    write(fid, spt.pz ); write(fid, spt.px   );
    flush(fid)
    return fid
end

function writeSnapShot(fid::IOStream, spt::SnapShot)
    write(fid, spt.vz ); write(fid, spt.vx   );
    write(fid, spt.pz ); write(fid, spt.px   );
    flush(fid)
    return nothing
end

function reverseSnapShotsOrder(path::ASCIIString, path_tmp::ASCIIString)
    (isx, isy, nz, nx, ext, iflag, ot, dt, nt) = InfoFullWfd(path_tmp)
    if iflag == 1
       Nz = nz + 2*ext
    elseif iflag == 2
       Nx = nx +   ext
    end
    Nx = nx + 2*ext
    fid = open(path, "w")
    write(fid, isx, isy, nz, nx, ext, iflag, ot, dt)
    fid_tmp = open(path_tmp, "r")
    for it = nt:-1:1
        position = sizeof(Float64)*8 + sizeof(Float64)*Nz*Nx*4*(it-1)
        seek(fid_tmp, position)
        d = read(fid_tmp, Float64, Nz*Nx*4);
        write(fid, d)
    end
    close(fid_tmp); close(fid)
    rm(path_tmp)
    return nothing
end

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
function plot(X::Embedding{T}, kwargs...) where T
    d, n = size(X)

    if d == 1
        plot(X', kwargs...)
    elseif
        d == 2
        plot(X[:,1], X[:,2])
    end
end

function plot!(X::Embedding{T}, kwargs...) where T
    d, n = size(X)

    if d == 1
        plot!(X, kwargs...)
    elseif
        d == 2
        plot!(X[:,1], X[:,2], kwargs...)
    end
end


text:
export DescentMethod, CoreMethod, Wrapper

abstract type DescentMethod end

abstract type CoreMethod <: DescentMethod end

abstract type Wrapper <: DescentMethod end


text:
function HandleDefinition(ps, cg)
    r = ParseDefinition(ps)
    ir = codegen(cg, r)
    push!(cg.jit, cg.mod)
    initialize_module_and_pass_manager!(cg)
    return ir
end

function HandleExtern(ps, cg)
    protoAST = ParseExtern(ps)
    ir = codegen(cg, protoAST)
    cg.function_protos[protoAST.name] = protoAST
    return ir
end

function HandleTopLevelExpression(ps, cg)
    r = ParseTopLevelExpr(ps)
    ir = codegen(cg, r)
    oldmod = cg.mod
    push!(cg.jit, cg.mod) do
        initialize_module_and_pass_manager!(cg)
        if !haskey(LLVM.functions(cg.jit), "__anon_expr")
            error("Did not find function in execution engine")
        end
        f = LLVM.functions(cg.jit)["__anon_expr"]
        res = LLVM.run(cg.jit, f)
        println("Got ", convert(Float64, res, LLVM.DoubleType()))
        LLVM.dispose(res)
    end
    LLVM.dispose(oldmod)
    return ir
end

function repl()
    cg = CodeGen()
    initialize_module_and_pass_manager!(cg)
    while true
        print("ready> ")
        str = readline()
        _generate_ir(str, cg)
    end
    # TODO: dispose
end

function generate_ir(str::String)
    cg = CodeGen()
    initialize_module_and_pass_manager!(cg)
    _generate_ir(str, cg)
end

function _generate_ir(str, cg::CodeGen)
    initialize_module_and_pass_manager!(cg)
    ps = Parser(str)
    tok = next_token!(ps)
    while true
        tok = current_token(ps)
        if tok.kind == Kinds.EOF
            break
        elseif tok.kind == Kinds.SEMICOLON
            next_token!(ps)
        elseif tok.kind == Kinds.DEF
            HandleDefinition(ps, cg)
        elseif tok.kind == Kinds.EXTERN
            HandleExtern(ps, cg)
        else
            HandleTopLevelExpression(ps, cg)
        end
    end
    return cg.mod
end

function emit_objectfile(str::String, path::String)
    cg = CodeGen()
    initialize_module_and_pass_manager!(cg)
    host_triple = LLVM.triple()
    host_t = LLVM.Target(host_triple)
    _generate_ir(str, cg)
    LLVM.TargetMachine(host_t, host_triple) do tm
        LLVM.emit(tm, cg.mod, LLVM.API.LLVMObjectFile, path)
    end
end

str = """
def fib(x) {
    if (x < 3) then
        1
    else
        fib(x-1)+fib(x-2)
}

fib(10)
"""

Kaleidoscope.generate_ir(str)

str = """
def fib(x) {
    var a = 1, b = 1
    for i = 3, i < x, 1.0 in {
        var c = a + b
        a = b
        b = c
    }
    b
}

fib(10)
"""

Kaleidoscope.generate_ir(str)


text:
T = Float32
_pkey = :z
tbl = table((x=rand(T, 10), y=rand(T, 10), z=collect(Int, 1:10)), pkey=_pkey)
nds = ndsparse(deepcopy(tbl))

@testset "DB: id_key checks, push!/pop! primitives" begin
    # db_create_schema
    @test sort(Garamond.db_create_schema(tbl), by=x->x.column)  ==
          sort(Garamond.db_create_schema(nds), by=x->x.column)  ==
          [(column=:x, coltype=T, pkey=false),
           (column=:y, coltype=T, pkey=false),
           (column=_pkey, coltype=Int, pkey=true)]

    # db_get_primary_keys
    @test Garamond.db_get_primary_keys(tbl) == Garamond.db_get_primary_keys(nds) == (_pkey,)

    # db_check_id_key
    @test Garamond.db_check_id_key(tbl, _pkey) == Garamond.db_check_id_key(nds, _pkey) == nothing
    @test_throws ErrorException Garamond.db_check_id_key(tbl, :some_random_key)
    @test_throws ErrorException Garamond.db_check_id_key(nds, :some_random_key)
    @test_throws ErrorException Garamond.db_check_id_key(tbl, :x)
    @test_throws ErrorException Garamond.db_check_id_key(nds, :x)

    # db_select_entry
    @test Garamond.db_select_entry(tbl, 1; id_key=_pkey) == first(rows(tbl))
    @test Garamond.db_select_entry(nds, 1; id_key=_pkey) == first(rows(nds))
    _tbl = Garamond.db_select_entry(tbl, 1_000; id_key=_pkey)
    @test _tbl isa IndexedTable && isempty(_tbl)
    _nds = Garamond.db_select_entry(nds, 1_000; id_key=_pkey)
    @test _nds isa NDSparse && isempty(_nds)

    # db_push!
    _len = length(tbl)
    data = (x=T(-100.0), y=T(-11.1), z=_len+1)
    @test Garamond.db_push!(tbl, data; id_key=_pkey) == nothing && length(tbl) == _len + 1
    @test Garamond.db_push!(nds, data; id_key=_pkey) == nothing && length(nds) == _len + 1
    #@test_throws ErrorException Garamond.db_push!(tbl, data)
    #@test_throws ErrorException Garamond.db_push!(nds, data)
    wrong_data = [(x=T(0.0), y=-T(11.0), z=100),         # wrong pkey
                  (x="a", y=T(1.0), z=_len+1),        # wrong eltype
                  (x=T(1.0), y=T(2.0), z=_len+1, γ="a"), # additional column
                  (x=T(1.0), y=T(2.0))]                  # missing column
    for wrong in wrong_data
        @test_throws ErrorException Garamond.db_push!(tbl, wrong; id_key=_pkey)
        @test_throws ErrorException Garamond.db_push!(nds, wrong; id_key=_pkey)
    end

    # db_pushfirst!
    _len = length(tbl)
    data = (x=T(-100.0), y=T(-11.1), z=1)
    @test Garamond.db_pushfirst!(tbl, data; id_key=_pkey) == nothing && length(tbl) == _len + 1
    @test getproperty(columns(tbl), _pkey) == collect(1:_len+1)
    @test Garamond.db_pushfirst!(nds, data; id_key=_pkey) == nothing && length(nds) == _len + 1
    @test getproperty(columns(nds), _pkey) == collect(1:_len+1)
    #@test_throws ErrorException Garamond.db_pushfirst!(tbl, data)
    #@test_throws ErrorException Garamond.db_pushfirst!(nds, data)
    wrong_data = [(x=T(0.0), y=-T(11.0), z=100),         # wrong pkey
                  (x="a", y=T(1.0), z=_len+1),        # wrong eltype
                  (x=T(1.0), y=T(2.0), z=_len+1, γ="a"), # additional column
                  (x=T(1.0), y=T(2.0))]                  # missing column
    for wrong in wrong_data
        @test_throws ErrorException Garamond.db_pushfirst!(tbl, wrong; id_key=_pkey)
        @test_throws ErrorException Garamond.db_pushfirst!(nds, wrong; id_key=_pkey)
    end

    # db_pop!
    _len = length(tbl)
    last_row_tbl = last(rows(tbl))
    last_row_nds = last(rows(nds))
    @test Garamond.db_pop!(tbl) == last_row_tbl && length(tbl) == _len - 1
    @test Garamond.db_pop!(nds) == last_row_nds && length(nds) == _len - 1

    # db_popfirst!
    _len = length(tbl)
    first_row_tbl = first(rows(tbl))
    first_row_nds = first(rows(nds))
    @test Garamond.db_popfirst!(tbl; id_key=_pkey) == first_row_tbl && length(tbl) == _len - 1
    @test getproperty(columns(tbl), _pkey) == collect(1:_len-1)
    @test Garamond.db_popfirst!(nds; id_key=_pkey) == first_row_nds && length(nds) == _len - 1
    @test getproperty(columns(nds), _pkey) == collect(1:_len-1)

    # db_deleteat!
    _len = length(tbl)
    to_delete = sort(unique(rand(1:_len, 3)))
    @test Garamond.db_deleteat!(tbl, to_delete; id_key=_pkey) == nothing &&
        length(tbl) == _len - length(to_delete)
    @test getproperty(columns(tbl), _pkey) == collect(1:(_len-length(to_delete)))
    @test Garamond.db_deleteat!(nds, to_delete; id_key=_pkey) == nothing &&
        length(nds) == _len - length(to_delete)
    @test getproperty(columns(nds), _pkey) == collect(1:(_len-length(to_delete)))

    # db_drop_columns
    @test Garamond.db_drop_columns(tbl, [:does_not_exist]) == tbl
    @test Garamond.db_drop_columns(nds, [:does_not_exist]) == nds
    @test !in(_pkey, colnames(Garamond.db_drop_columns(tbl, [_pkey])))
    @test Garamond.db_drop_columns(nds, [_pkey]) == nds  # cannot drop index
    dropped_x = Garamond.db_drop_columns(tbl, [:x])
    @test !(:x in colnames(dropped_x))
    dropped_x = Garamond.db_drop_columns(tbl, [:x])
    @test !(:x in colnames(dropped_x))
end


text:
@everywhere Random.seed!(42)

weights = rand(400)
weights = weights ./ sum(weights)

test_sys_resample    = SMC.resample(weights, method = :systematic)
test_multi_resample  = SMC.resample(weights, method = :multinomial)
test_poly_resample   = SMC.resample(weights, method = :polyalgo)

saved_sys_resample   = load("reference/resample.jld2", "sys")
saved_multi_resample = load("reference/resample.jld2", "multi")
saved_poly_resample  = load("reference/resample.jld2", "poly")

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

@testset "Resampling methods" begin
    @test test_sys_resample   == saved_sys_resample
    @test test_multi_resample == saved_multi_resample
    @test test_poly_resample  == saved_poly_resample
end


text:
using Oceananigans.Grids: cpu_face_constructor_x, cpu_face_constructor_y, cpu_face_constructor_z, default_indices
using Oceananigans.BoundaryConditions: CBC, PBC

struct YPartition{N} <: AbstractPartition
    div :: N
    function YPartition(sizes)
        if length(sizes) > 1 && all(y -> y == sizes[1], sizes)
            sizes = length(sizes)
        end
        return new{typeof(sizes)}(sizes)
    end
end

const EqualYPartition = YPartition{<:Number}

length(p::EqualYPartition) = p.div
length(p::YPartition)      = length(p.div)

Base.summary(p::EqualYPartition) = "Equal partitioning in Y ($(p.div) regions)"
Base.summary(p::YPartition)      = "partitioning in Y [$(["$(p.div[i]) " for i in 1:length(p)]...)]"

function partition_size(p::EqualYPartition, grid)
    Nx, Ny, Nz = size(grid)
    @assert mod(Ny, p.div) == 0
    return Tuple((Nx, Ny ÷ p.div, Nz) for i in 1:length(p))
end

function partition_size(p::YPartition, grid)
    Nx, Ny, Nz = size(grid)
    @assert sum(p.div) != Ny
    return Tuple((Nx, p.div[i], Nz) for i in 1:length(p))
end

function partition_extent(p::YPartition, grid)
    x = cpu_face_constructor_x(grid)
    y = cpu_face_constructor_y(grid)
    z = cpu_face_constructor_z(grid)

    y = divide_direction(y, p)
    return Tuple((x, y = y[i], z = z) for i in 1:length(p))
end

function partition_topology(p::YPartition, grid)
    TX, TY, TZ = topology(grid)

    return Tuple((TX, (TY == Periodic ?
                       FullyConnected :
                       i == 1 ?
                       RightConnected :
                       i == length(p) ?
                       LeftConnected :
                       FullyConnected), TZ) for i in 1:length(p))
end

divide_direction(x::Tuple, p::EqualYPartition) =
    Tuple((x[1]+(i-1)*(x[2] - x[1])/length(p), x[1]+i*(x[2] - x[1])/length(p)) for i in 1:length(p))

function divide_direction(x::AbstractArray, p::EqualYPartition)
    nelem = (length(x)-1)÷length(p)
    return Tuple(x[1+(i-1)*nelem:1+i*nelem] for i in 1:length(p))
end

function partition_global_array(a::AbstractArray, ::EqualYPartition, local_size, region, arch)
    idxs = default_indices(length(size(a)))
    return arch_array(arch, a[idxs[1], local_size[2]*(region-1)+1:local_size[2]*region, idxs[3:end]...])
end

function partition_global_array(a::OffsetArray, ::EqualYPartition, local_size, region, arch)
    idxs    = default_indices(length(size(a)))
    offsets = (0, a.offsets[2], Tuple(0 for i in 1:length(idxs)-2)...)
    return arch_array(arch, OffsetArray(a[idxs[1], local_size[2]*(region-1)+1+offsets[2]:local_size[2]*region-offsets[2], idxs[3:end]...], offsets...))
end

####
#### Global reconstruction utils
####

function reconstruct_size(mrg, p::YPartition)
    Nx = mrg.region_grids[1].Nx
    Ny = sum([grid.Ny for grid in mrg.region_grids.regions])
    Nz = mrg.region_grids[1].Nz
    return (Nx, Ny, Nz)
end

function reconstruct_extent(mrg, p::YPartition)
    switch_device!(mrg.devices[1])
    x = cpu_face_constructor_x(mrg.region_grids.regions[1])
    z = cpu_face_constructor_z(mrg.region_grids.regions[1])

    if cpu_face_constructor_y(mrg.region_grids.regions[1]) isa Tuple
        y = (cpu_face_constructor_y(mrg.region_grids.regions[1])[1],
             cpu_face_constructor_y(mrg.region_grids.regions[length(p)])[end])
    else
        y = [cpu_face_constructor_y(mrg.region_grids.regions[1])...]
        for (idx, grid) in enumerate(mrg.region_grids.regions[2:end])
            switch_device!(mrg.devices[idx])
            y = [y..., cpu_face_constructor_y(grid)[2:end]...]
        end
    end
    return (; x = x, y = y, z = z)
end

function reconstruct_global_array(ma::ArrayMRO{T, N}, p::EqualYPartition, arch) where {T, N}
    local_size = size(first(ma.regions))
    global_Ny  = local_size[2] * length(p)
    idxs = default_indices(length(local_size))
    arr_out = zeros(eltype(first(ma.regions)), local_size[1], global_Ny, local_size[3:end]...)
    n = local_size[2]
    for r = 1:length(p)
        init = Int(n * (r - 1) + 1)
        fin  = Int(n * r)
        arr_out[idxs[1], init:fin, idxs[3:end]...] .= arch_array(CPU(), ma[r])[idxs[1], 1:fin-init+1, idxs[3:end]...]
    end

    return arch_array(arch, arr_out)
end

function compact_data!(global_field, global_grid, data::MultiRegionObject, p::EqualYPartition)
    Nx, Ny, Nz = size(global_grid)
    n = Ny / length(p)
    for r = 1:length(p)
        init = Int(n * (r - 1) + 1)
        fin  = Int(n * r)
        interior(global_field)[:, init:fin, :] .= data[r][:, 1:fin-init+1, :]
    end
end

#####
##### Boundary specific Utils
#####

inject_west_boundary(region, p::YPartition, bc) = bc
inject_east_boundary(region, p::YPartition, bc) = bc

function inject_south_boundary(region, p::YPartition, global_bc)
    if region == 1
        typeof(global_bc) <: Union{CBC, PBC} ?
                bc = CommunicationBoundaryCondition((rank = region, from_rank = length(p))) :
                bc = global_bc
    else
        bc = CommunicationBoundaryCondition((rank = region, from_rank = region - 1))
    end
    return bc
end

function inject_north_boundary(region, p::YPartition, global_bc)
    if region == length(p)
        typeof(global_bc) <: Union{CBC, PBC} ?
                bc = CommunicationBoundaryCondition((rank = region, from_rank = 1)) :
                bc = global_bc
    else
        bc = CommunicationBoundaryCondition((rank = region, from_rank = region + 1))
    end
    return bc
end

####
#### Global index flattening
####

@inline function displaced_xy_index(i, j, grid, region, p::YPartition)
    j′ = j + grid.Ny * (region - 1)
    t  = i + (j′ - 1) * grid.Nx
    return t
end