text:
using SubmodularMaximization
using POMDPs
using MCTS
using PyPlot

pygui(false)

steps = 20
grid_size = 10
horizon = SubmodularMaximization.default_horizon
iterations = 100

grid = Grid(grid_size, grid_size)
sensor = RangingSensor()

robot_state = random_state(grid)
target_state = random_state(grid)

target_states = Array{State}(undef, steps)
target_states[1] = target_state

robot_states = Array{State}(undef, steps)
robot_states[1] = robot_state

histogram_filter = Filter(grid, target_state)

problem = SingleRobotTargetTrackingProblem(grid, sensor, horizon,
                                           [histogram_filter])

plot_state_space(grid)
xlim([0, grid_size+1])
ylim([0, grid_size+1])

for ii = 2:steps
  println("Step ", ii)

  # Before the target moves and the robot receives a measurement, execute robot
  # dynamics
  global robot_state = solve_single_robot(problem, robot_state,
                                          n_iterations = iterations)
  robot_states[ii] = robot_state

  # Then update the target and sample the observation
  global target_state = target_dynamics(grid, target_states[ii-1])
  target_states[ii] = target_state

  range_observation = generate_observation(grid, sensor, robot_state,
                                           target_state)

  # After updating states, update the filter (in place for now)
  process_update!(histogram_filter, transition_matrix(grid))
  measurement_update!(histogram_filter, robot_state, get_states(grid), sensor,
                      grid, range_observation)

  plots=[]
  append!(plots, plot_trajectory(robot_states[1:ii], color=:blue))
  append!(plots, plot_observation(robot_state, range_observation, color=:blue))
  append!(plots, plot_trajectory(target_states[1:ii]))
  append!(plots, visualize_filter(histogram_filter))

  savefig(string("fig/single_robot_gif_", ii, ".png"))

  foreach(x->x.remove(), plots)
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
###############################################################################
#
# Big leaf canopy model partition
#
###############################################################################
"""
    big_leaf_partition!(
                node::SPACSimple{FT},
                zenith::FT,
                r_all::FT
    ) where {FT <:AbstractFloat}

Partition the big-leaf canopy into sunlit and shaded layers, given
- `partition` Container for partition
- `zenith` Zenith angle in degree
- `r_all` Total radiation in `[W m⁻²]`
"""
function big_leaf_partition!(
            node::SPACSimple{FT},
            zenith::FT,
            r_all::FT
) where {FT <:AbstractFloat}
    # unpack values
    @unpack laba, lai = node;

    # 1. use big_leaf_partition function from CanopyLayers module
    ratio, q_slm, q_shm, e_sl, e_sh = big_leaf_partition(lai, zenith, r_all);

    # 2. update the information
    node.container2L.frac_sl = ratio;
    node.container2L.frac_sh = 1 - ratio;
    node.container2L.la_sl   = laba * ratio;
    node.container2L.la_sh   = laba * (1 - ratio);
    node.container2L.lai_sl  = lai * ratio;
    node.container2L.lai_sh  = lai * (1 - ratio);
    node.container2L.par_sl  = q_slm;
    node.container2L.par_sh  = q_shm;
    node.container2L.rad_sl  = e_sl;
    node.container2L.rad_sh  = e_sh;

    return nothing
end

ext: jl

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
# 32-bit fixed point; parameter `f` is the number of fraction bits
struct Fixed{T <: Signed,f} <: FixedPoint{T,  f}
    i::T

    # constructor for manipulating the representation;
    # selected by passing an extra dummy argument
    Fixed{T, f}(i::Integer, _) where {T,f} = new{T, f}(i % T)
end

Fixed{T, f}(x::AbstractChar) where {T,f} = throw(ArgumentError("Fixed cannot be constructed from a Char"))
Fixed{T, f}(x::Complex) where {T,f} = Fixed{T, f}(convert(real(typeof(x)), x))
Fixed{T, f}(x::Base.TwicePrecision) where {T,f} = Fixed{T, f}(convert(Float64, x))
Fixed{T,f}(x::Integer) where {T,f} = Fixed{T,f}(round(T, convert(widen1(T),x)<<f),0)
Fixed{T,f}(x::AbstractFloat) where {T,f} = Fixed{T,f}(round(T, trunc(widen1(T),x)<<f + rem(x,1)*(one(widen1(T))<<f)),0)
Fixed{T,f}(x::Rational) where {T,f} = Fixed{T,f}(x.num)/Fixed{T,f}(x.den)

reinterpret(::Type{Fixed{T,f}}, x::T) where {T <: Signed,f} = Fixed{T,f}(x, 0)

typechar(::Type{X}) where {X <: Fixed} = 'Q'
signbits(::Type{X}) where {X <: Fixed} = 1

for T in (Int8, Int16, Int32, Int64)
    for f in 0:sizeof(T)*8-1
        sym = Symbol(String(take!(showtype(_iotypealias, Fixed{T,f}))))
        @eval begin
            const $sym = Fixed{$T,$f}
            export $sym
        end
    end
end

# basic operators
-(x::Fixed{T,f}) where {T,f} = Fixed{T,f}(-x.i,0)
abs(x::Fixed{T,f}) where {T,f} = Fixed{T,f}(abs(x.i),0)

+(x::Fixed{T,f}, y::Fixed{T,f}) where {T,f} = Fixed{T,f}(x.i+y.i,0)
-(x::Fixed{T,f}, y::Fixed{T,f}) where {T,f} = Fixed{T,f}(x.i-y.i,0)

# with truncation:
#*{f}(x::Fixed32{f}, y::Fixed32{f}) = Fixed32{f}(Base.widemul(x.i,y.i)>>f,0)
# with rounding up:
*(x::Fixed{T,f}, y::Fixed{T,f}) where {T,f} = Fixed{T,f}((Base.widemul(x.i,y.i) + (one(widen(T)) << (f-1)))>>f,0)

/(x::Fixed{T,f}, y::Fixed{T,f}) where {T,f} = Fixed{T,f}(div(convert(widen(T), x.i) << f, y.i), 0)

# # conversions and promotions
function Fixed{T,f}(x::Fixed{T2,f2}) where {T <: Integer,T2 <: Integer,f,f2}
#    reinterpret(Fixed{T,f},T(reinterpret(x)<<(f-f2)))
    U = Fixed{T,f}
    y = round(((1<<f)/(1<<f2))*reinterpret(x))
    (typemin(T) <= y) & (y <= typemax(T)) || throw_converterror(U, x)
    reinterpret(U, _unsafe_trunc(T, y))
end

rem(x::Integer, ::Type{Fixed{T,f}}) where {T,f} = Fixed{T,f}(rem(x,T)<<f,0)
rem(x::Real,    ::Type{Fixed{T,f}}) where {T,f} = Fixed{T,f}(rem(Integer(trunc(x)),T)<<f + rem(Integer(round(rem(x,1)*(one(widen1(T))<<f))),T),0)

float(x::Fixed) = convert(floattype(x), x)

Base.BigFloat(x::Fixed{T,f}) where {T,f} =
    BigFloat(x.i>>f) + BigFloat(x.i&(one(widen1(T))<<f - 1))/BigFloat(one(widen1(T))<<f)
(::Type{TF})(x::Fixed{T,f}) where {TF <: AbstractFloat,T,f} =
    TF(x.i>>f) + TF(x.i&(one(widen1(T))<<f - 1))/TF(one(widen1(T))<<f)

Base.Bool(x::Fixed{T,f}) where {T,f} = x.i!=0
function Base.Integer(x::Fixed{T,f}) where {T,f}
    isinteger(x) || throw(InexactError())
    Integer(x.i>>f)
end
function (::Type{TI})(x::Fixed{T,f}) where {TI <: Integer,T,f}
    isinteger(x) || throw(InexactError())
    TI(x.i>>f)
end

(::Type{TR})(x::Fixed{T,f}) where {TR <: Rational,T,f} =
    TR(x.i>>f + (x.i&(1<<f-1))//(one(widen1(T))<<f))

promote_rule(ft::Type{Fixed{T,f}}, ::Type{TI}) where {T,f,TI <: Integer} = Fixed{T,f}
promote_rule(::Type{Fixed{T,f}}, ::Type{TF}) where {T,f,TF <: AbstractFloat} = TF
promote_rule(::Type{Fixed{T,f}}, ::Type{Rational{TR}}) where {T,f,TR} = Rational{TR}

@generated function promote_rule(::Type{Fixed{T1,f1}}, ::Type{Fixed{T2,f2}}) where {T1,T2,f1,f2}
    f = max(f1, f2)  # ensure we have enough precision
    T = promote_type(T1, T2)
    # make sure we have enough integer bits
    i1, i2 = 8*sizeof(T1)-f1, 8*sizeof(T2)-f2  # number of integer bits for each
    i = 8*sizeof(T)-f
    while i < max(i1, i2)
        Tw = widen1(T)
        T == Tw && break
        T = Tw
        i = 8*sizeof(T)-f
    end
    :(Fixed{$T,$f})
end

# TODO: Document and check that it still does the right thing.
decompose(x::Fixed{T,f}) where {T,f} = x.i, -f, 1


text:
# Illustrate functionality of the code with a test suite, heuristics and Gaussian quadrature
# About
#   Author       - Peter Opsomer (Peter.Opsomer@cs.kuleuven.be)
#   History      - Created October 2014, last edit October 2015

## Test suite: test all functionality of the code.
using PyPlot

include("validation.jl")
include("precomputations.jl")
include("asy_expansions.jl")
include("UExplicit.jl")

alpha = 2*rand()-0.9; beta = 3*rand()-0.7; cc = rand(); dd = rand();
function h(x) return exp(cc*x) + sin(dd*x).^2+1; end
function dh(x) return cc*exp(cc*x) +2.*sin(dd*x).*dd*cos(dd*x); end

maxOrder = rand(2:10,1)[1];
# maxP2 lower-> maybe no convergence, higher-> roundoff errors in exactPolys
maxP2 = 7;
shift = (rand(0:3,1))[1];
print("alpha=", alpha, ", beta=", beta, ", cc=", cc, ", dd=", dd, ", maxOrder=", maxOrder, ", shift=", shift, '\n')
(P,gammaP,alphaP,betaP) = exactPolys(alpha,beta,h,2^maxP2+shift);

rs = rand(1:27,1,3); # The indexes of the regions to be tested
zr = (0+0im)*zeros(27,1); # At most 3 of these will be used
th = (0+0im)*zeros(27,1);
#  Taking rand iso randn for abs(zr) since it make points not lie inside the Bernstein ellipse
for ri = rs
	if ri <= 3
	elseif ri <= 7 zr[ri] = 0.5*rand()*exp(2im*rand()*pi); # lens
	elseif ri <= 11 zr[ri] = 1.1*rand()*exp(2im*rand()*pi); # outer
	elseif ri <= 15 zr[ri] = 1 + 0.4*rand()*exp(2im*rand()*pi); # right disk
	elseif ri <= 19 zr[ri] = -1 + 0.4*rand()*exp(2im*rand()*pi); # left disk
	elseif ri <= 23 th[ri] = eps()^0.4*rand()*exp(2im*rand()*pi); zr[ri] = cos(th[ri]); # right disk Q
	else th[ri] = eps()^0.4*rand()*exp(2im*rand()*pi); zr[ri] = -cos(th[ri]); # left disk Q
	end
end

legs = ["Recurrence coefficients alpha_n","Recurrence coefficients beta_n","Leading order coefficients gamma_n of the orthonormal polynomials",# = 3
	"Lens monic","Derivative lens monic","Lens orthonormal","Derivative lens orthonormal", # = 7
	"Outer monic","Derivative outer monic","Outer orthonormal","Derivative outer orthonormal", # = 11
	"Right disk monic","Derivative right disk monic","Right disk orthonormal","Derivative right disk orthonormal", # = 15
	"Left disk monic","Derivative Left disk monic","Left disk orthonormal","Derivative left disk orthonormal", # = 19
	"Right monic with Q-s","Derivative right monic with Q-s","orthonormal Right with Q-s", "Derivative orthonormal Right with Q-s", # = 23
	"Left monic with Q-s","Derivative left monic with Q-s","orthonormal left with Q-s","derivative orthonormal left with Q-s"]; # = 27

tic()
nrQ = ceil(Integer, 2+2*maxOrder*ones(1,maxOrder-1));
rho = 3; # Trapezium rules may become inaccurate if rho is (much) higher
if (maximum(abs(imag(zr))) > (rho/2 - 1/rho/2)*0.9 ) || (maximum(real(zr)) > (rho/2 + 1/rho/2)*0.9 ) ||	(minimum(real(zr)) < -(rho/2 + 1/rho/2)*0.9 )
	warn("Possible errors when evaluating far from the interval")
end

# All zr should lie within the contour. For Q-s, size(c&d) = ceil(nrT/2) so get higher number of Q-s by
(c, d, Dinf, psi, dpsi, contpsi) = contour_integrals(alpha,beta,h,4*maximum(nrQ+1),rho,2000);
(Uright,Uleft,Qright,Qleft) = UQ(alpha,beta,Dinf,c,d,maxOrder,"UQW",nrQ);
timePrecompute = toc()

(Ure,Ule) = UExplicit(alpha,beta,Dinf,c,d);
rk = min(size(Uright,3), size(Ure,3), 6); rm = min(size(Uright,4), size(Ure,4), 3);
lk = min(size(Uleft,3), size(Ule,3), 6); lm = min(size(Uleft,4), size(Ule,4), 3);
errorUright = sum(sum(sum(sum((Ure[:,:,1:rk,1:rm]-Uright[:,:,1:rk,1:rm]).^2))))
errorUleft = sum(sum(sum(sum((Ule[:,:,1:lk,1:lm]-Uleft[:,:,1:lk,1:lm]).^2))))
print("Error in Uright = ", errorUright, " and in Uleft = ", errorUleft, '\n')

exact = (0+0im)*zeros(maxP2,1,27);
asy = (0+0im)*zeros(maxP2,maxOrder,27);
gammas = (0+0im)*zeros(maxP2,1);
hh = sqrt(eps());
# Using complex finite differences (imag(P(zr[ri]+1i*hh,nt))/hh;) as explained in "William Squire and George Trapp. Using complex variables to estimate derivatives fo real functions. SIAM Review, 40(1) pp. 110-112, 1998." could make that hh may be smaller and thus better estimate of derivative, but not if zr is not real
for tn = 1:maxP2
	nt = convert(Int64,2^tn+shift);
	gammas[tn] = gammaP[nt+1];
	    for ri = rs
		if ri == 1 exact[tn,1,ri] = alphaP[nt+1];
			for i = 1:maxOrder   asy[tn,i,ri] = alphan(nt,i,Uright,Uleft); end
		elseif ri == 2 exact[tn,1,ri] = betaP[nt+1];
			for i = 1:maxOrder	asy[tn,i,ri] = betan(nt,i,Dinf,Uright,Uleft); end
		elseif ri == 3 exact[tn,1,ri] = gammaP[nt+1];
			for i = 1:maxOrder   asy[tn,i,ri] = gamman(nt,i,Dinf,Uright,Uleft); end
		elseif ri == 4	exact[tn,1,ri] = P(zr[ri],nt)/gammas[tn];

			for i = 1:maxOrder   asy[tn,i,ri] = asy_lens(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"m"); end
		elseif ri == 5	exact[tn,1,ri] = (P(zr[ri],nt)-P(zr[ri]-hh,nt))/hh/gammas[tn];
			for i = 1:maxOrder
				(_, t) = asy_lens(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"m",dh,dpsi); asy[tn,i,ri] = t;
			end
		elseif ri == 6	exact[tn,1,ri] = P(zr[ri],nt);
			for i = 1:maxOrder   asy[tn,i,ri] = asy_lens(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"o"); end
		elseif ri == 7	exact[tn,1,ri] = (P(zr[ri],nt)-P(zr[ri]-hh,nt))/hh;
			for i = 1:maxOrder
				(_, t) = asy_lens(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"o",dh,dpsi); asy[tn,i,ri] = t;
			end
		elseif ri == 8	exact[tn,1,ri] = P(zr[ri],nt)/gammas[tn];
			for i = 1:maxOrder   asy[tn,i,ri] = asy_outer(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"m"); end
		elseif ri == 9	exact[tn,1,ri] = (P(zr[ri],nt)-P(zr[ri]-hh,nt))/hh/gammas[tn];
			for i = 1:maxOrder
				(_, t) = asy_outer(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"m",dh,dpsi);	asy[tn,i,ri] = t;
			end
		elseif ri == 10	exact[tn,1,ri] = P(zr[ri],nt);
			for i = 1:maxOrder   asy[tn,i,ri] = asy_outer(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"o"); end
		elseif ri == 11	exact[tn,1,ri] = (P(zr[ri],nt)-P(zr[ri]-hh,nt))/hh;
			for i = 1:maxOrder
				(_, t) = asy_outer(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"o",dh,dpsi);	asy[tn,i,ri] = t;
			end
		elseif ri == 12	exact[tn,1,ri] = P(zr[ri],nt)/gammas[tn];
			for i = 1:maxOrder    asy[tn,i,ri] = asy_right(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"m"); end
		elseif ri == 13	exact[tn,1,ri] = (P(zr[ri],nt)-P(zr[ri]-hh,nt))/hh/gammas[tn];
			for i = 1:maxOrder
				(_, t) = asy_right(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"m",dh,dpsi);
				asy[tn,i,ri] = t;
			end
		elseif ri == 14	exact[tn,1,ri] = P(zr[ri],nt);
			for i = 1:maxOrder    asy[tn,i,ri] = asy_right(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"o"); end
		elseif ri == 15	exact[tn,1,ri] = (P(zr[ri],nt)-P(zr[ri]-hh,nt))/hh;
			for i = 1:maxOrder
				(_, t) = asy_right(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"o",dh,dpsi);
				asy[tn,i,ri] = t;
			end
		elseif ri == 16	exact[tn,1,ri] = P(zr[ri],nt)/gammas[tn];
			for i = 1:maxOrder    asy[tn,i,ri] = asy_left(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"m"); end
		elseif ri == 17
			exact[tn,1,ri] = (P(zr[ri],nt)-P(zr[ri]-hh,nt))/hh/gammas[tn];
			for i = 1:maxOrder
				(_ , t) = asy_left(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"m",dh,dpsi);
				asy[tn,i,ri] = t;
			end
		elseif ri == 18	exact[tn,1,ri] = P(zr[ri],nt);
			for i = 1:maxOrder    asy[tn,i,ri] = asy_left(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"o"); end
		elseif ri == 19	exact[tn,1,ri] = (P(zr[ri],nt)-P(zr[ri]-hh,nt))/hh;
			for i = 1:maxOrder
				(_, t) = asy_left(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"o",dh,dpsi);
				asy[tn,i,ri] = t;
			end
		elseif ri == 20	exact[tn,1,ri] = P(zr[ri],nt)/gammas[tn];
			for i = 1:maxOrder
				asy[tn,i,ri] = asy_right(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"mQ",false,false,Qright,th[ri],contpsi);
			end
		elseif ri == 21	exact[tn,1,ri] = (P(zr[ri],nt)-P(zr[ri]-hh,nt))/hh/gammas[tn];
			for i = 1:maxOrder
				(_, t) = asy_right(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"mQ",dh,dpsi,Qright,th[ri],contpsi);
				asy[tn,i,ri] = t;
			end
		elseif ri == 22	exact[tn,1,ri] = P(zr[ri],nt);
			for i = 1:maxOrder
				asy[tn,i,ri] = asy_right(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"oQ",false,false,Qright,th[ri],contpsi); end
		elseif ri == 23	exact[tn,1,ri] = (P(zr[ri],nt)-P(zr[ri]-hh,nt))/hh;
			for i = 1:maxOrder
				(_, t) = asy_right(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"oQ",dh,dpsi,Qright,th[ri],contpsi);
				asy[tn,i,ri] = t;
			end

		elseif ri == 24	exact[tn,1,ri] = P(zr[ri],nt)/gammas[tn];
			for i = 1:maxOrder
				asy[tn,i,ri] = asy_left(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"mQ",false,false,Qleft,th[ri],contpsi); end
		elseif ri == 25	exact[tn,1,ri] = (P(zr[ri],nt)-P(zr[ri]-hh,nt))/hh/gammas[tn];
			for i = 1:maxOrder
				(_, t) = asy_left(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"mQ",dh,dpsi,Qleft,th[ri],contpsi);
				asy[tn,i,ri] = t;
			end
		elseif ri == 26	exact[tn,1,ri] = P(zr[ri],nt);
			for i = 1:maxOrder
				asy[tn,i,ri] = asy_left(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"oQ",false,false,Qleft,th[ri],contpsi); end
		elseif ri == 27	exact[tn,1,ri] = (P(zr[ri],nt)-P(zr[ri]-hh,nt))/hh;
			for i = 1:maxOrder
				(_, t) = asy_left(nt,zr[ri],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"oQ",dh,dpsi,Qleft,th[ri],contpsi);
				asy[tn,i,ri] = t;
			end
		end
	end
end
# Testing order of convergence by plotting the errors in the (3 randomly) chosen regions, take into account the remarks in the article when interpreting these
for ri = rs
	te = (0+0im)*zeros(maxP2,1);
	ta = (0+0im)*zeros(maxP2,maxOrder);
	for tn = 1:maxP2
		te[tn,1] = exact[tn,1,ri];
		for i = 1:maxOrder
			ta[tn,i] = asy[tn,i,ri];
		end
	end
	if ri <= 3 plotConv(te,ta, legs[ri]);
	elseif ri <= 19 plotConv(te,ta, string(legs[ri], " at z= ", zr[ri]) );
	else plotConv(te, ta, string(legs[ri], " at theta= ", th[ri]));
	end
end

## Heuristics: see which expansion is more accurate where.
workspace();
using PyPlot
include("validation.jl"); include("precomputations.jl"); include("asy_expansions.jl");

alpha = rand(); beta = -1+3*rand(); function h(x) 1; end
print("Heuristics: alpha=", alpha, ", beta=", beta, '\n')
n = 100;
maxOrder = 3;

(P,gammaP,alphaP,betaP) = exactPolys(alpha,beta,h,n);
# Computing relative errors on nrP points on this line
zs = linspace(-1,1,25) +0.02im;

(c, d, Dinf, psi, _, contpsi) = contour_integrals(alpha,beta,h,2*(maxOrder+3) );
(Uright,Uleft,Qright,Qleft) = UQ(alpha,beta,Dinf,c,d,maxOrder,"UQV");

relerrs = zeros(length(zs),maxOrder,6);
# Will give many warnings when evaluating the asymptotic expansions in the disks with Q-s far from endpoint. "Currently the only fix is to eliminate the warning", so do $ julia illustrateFunctionality.jl 2> junk
for zi = 1:length(zs)
	exact = P(zs[zi],n);
	for i = 1:maxOrder
		relerrs[zi,i,1] =  abs(asy_lens( n,zs[zi],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"o")-exact)/abs(exact);
		relerrs[zi,i,2] =  abs(asy_outer(n,zs[zi],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"o")-exact)/abs(exact);
		relerrs[zi,i,3] =  abs(asy_right(n,zs[zi],alpha,beta,h,psi,i,Dinf,Uright,Uleft,"o")-exact)/abs(exact);
		relerrs[zi,i,4] =  abs(asy_left( n,zs[zi],alpha,beta,h,psi,i,Dinf, Uright,Uleft,"o")-exact)/abs(exact);
		relerrs[zi,i,5] = abs(asy_right(n,zs[zi],alpha,beta,h,psi,i,Dinf, Uright,Uleft,"oQ",false,false,Qright,acos(zs[zi]),contpsi)-exact)/abs(exact);
		relerrs[zi,i,6] = abs(asy_left(n,zs[zi],alpha,beta,h,psi,i,Dinf, Uright,Uleft,"oQ",false,false,Qleft,acos(-zs[zi]),contpsi)-exact)/abs(exact);
	end
end

colors = "bgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmyk";
semilogy(real(zs),relerrs[:,1,4],string("o",colors[1]),label=string("1 term") )
for i = 2:maxOrder
	semilogy(real(zs),relerrs[:,i,4],string("o",colors[i]),label=string(i," terms") )
end
legend(loc=0);
xlabel(string("Real part of z, imaginary part is ", mean(imag(zs))) );
ylabel("Relative error");
title("Relative errors of expansion in the left disk for varying z and number of terms");
show()

T = 3;
legs =  ["lens","outer","right","left","right Q","left Q"];
for i = 1:6
	semilogy(real(zs),relerrs[:,T,i],string("o",colors[i]),label=legs[i])
end
legend(loc=10);
xlabel(string("Real part of z, imaginary part is ", mean(imag(zs))));
ylabel("Relative error");
title(string("Relative errors of different expansion for varying z and T=", T))
show()

## Gaussian quadrature: Newton procedure to compute nodes and weights
workspace(); include("precomputations.jl"); include("asy_expansions.jl");

cheb = rand(Bool);
if cheb
	alpha = -1/2; beta = -1/2;
	function h(x) 1; end;
	function dh(x) 0; end
else
	alpha = 1/sqrt(3); beta = -1/pi;
	function h(x) exp(x); end;
	function dh(x) exp(x); end
end
n = 200;
nrT = 3;

(c, d, Dinf, psi, dpsi, _) = contour_integrals(alpha,beta,h,nrT);
(Uright,Uleft) = UQ(alpha,beta,Dinf,c,d,nrT,"UW");

zw = zeros(n,3);
for ni = 1:n
	if ni <=3 # (Approximation of) zero of bessel function to approximate first zeros
		zw[ni,1] = -1 + (2*sqrt(beta+1) + (ni-1)*pi)^2/2/n^2;
		if zw[ni,1] +1 < eps()^(1/3)
			warn("Use Q-s");
		end
	else
		zw[ni,1] = 3*zw[ni-1] -3*zw[ni-2] +zw[ni-3]; # Quadratic extrapolation
	end
	dpoly = "saveFromLoop"; polyM1 = "sfl"; polyP1 = "sfl";
	# Could improve the following by stopping Newton iterates when already tried that zero and then going over all (max. 5) closest floats.
	for attempt = 1:10
		if zw[ni,1]+1 < 0.2 # This bound is not general for all weights!
			(poly, dpoly) = asy_left(n,zw[ni,1],alpha,beta,h,psi,nrT,Dinf,Uright,Uleft,"o",dh,dpsi);
			polyM1 = asy_left(n-1,zw[ni,1],alpha,beta,h,psi,nrT,Dinf,Uright,Uleft,"o");
			polyP1 = asy_left(n+1,zw[ni,1],alpha,beta,h,psi,nrT,Dinf,Uright,Uleft,"o");
		elseif 1-zw[ni,1] < 0.2
			(poly, dpoly) = asy_right(n,zw[ni,1],alpha,beta,h,psi,nrT,Dinf,Uright,Uleft,"o",dh,dpsi);
			polyM1 = asy_right(n-1,zw[ni,1],alpha,beta,h,psi,nrT,Dinf,Uright,Uleft,"o");
			polyP1 = asy_right(n+1,zw[ni,1],alpha,beta,h,psi,nrT,Dinf,Uright,Uleft,"o");
		else
			(poly, dpoly) = asy_lens(n,zw[ni,1],alpha,beta,h,psi,nrT,Dinf, Uright,Uleft,"o",dh,dpsi);
			polyM1 = asy_lens(n-1,zw[ni,1],alpha,beta,h,psi,nrT,Dinf, Uright,Uleft,"o");
			polyP1 = asy_lens(n+1,zw[ni,1],alpha,beta,h,psi,nrT,Dinf, Uright,Uleft,"o");
		end
		zw[ni,1] = zw[ni,1] - real(poly)/real(dpoly); # We know the polynomial is real on the interval
	end
	if (zw[ni,1] < -1) || (zw[ni,1] > 1)
		error("Outside of the interval");
	elseif (ni > 1) && (zw[ni,1] < zw[ni-1,1] )
		error("Not going right");
	end
	zw[ni,2] = gamman(n,nrT,Dinf,Uright,Uleft)/real(dpoly)/gamman(n-1,nrT,Dinf,Uright,Uleft)/real(polyM1);
	zw[ni,3] = -gamman(n+1,nrT,Dinf,Uright,Uleft)/real(dpoly)/gamman(n,nrT,Dinf,Uright,Uleft)/real(polyP1);
end
errzw = norm(zw[:,2]-zw[:,3])/norm(zw[:,2]) # Formulas should give the same result

# Test by integrating (x^2-1) and (x^3) over the interval
ae1 = "saveFromIf"; ae2 = "sfi";
if cheb # % Weight 1/sqrt(1-x^2)
	chebyroots = cos((2*(1:n)-1)/2/n*pi);
	errorOnRoots = norm(zw[:,1]+chebyroots)
	print("Error in roots = ", errorOnRoots, '\n')
	ae1 = sum(zw[:,2].*(zw[:,1].^2-1) ) - (-pi/2)
	ae2 = sum(zw[:,2].*(zw[:,1].^3) )
else # w(x) = (1-x)^( 1/sqrt(3))*(1+x)^(-1/pi)*exp(x)
	ae1 = sum(zw[:,2].*(zw[:,1].^2-1) ) -(-1.293225037)
	ae2 = sum(zw[:,2].*(zw[:,1].^3) ) - (-.1659585353)
end
# Integrate a random polynomial with a degree low enough to be integrated exactly
r = [randn(2,1); rand(0:(2*n-3),1); rand(0:5,1)];
function fct(x) return 7*r[1]*x.^r[3] -23*r[2]*x.^r[4]; end
function w(x) return (1-x).^alpha.*(1+x).^beta.*h(x); end
# Better use the Cubature package for approximating the integral
xs = linspace(1/20000-1,1-1/20000,10000)
ae3 = sum(zw[:,2].*fct(zw[:,1]) ) -sum(w(xs).*fct(xs))/5000
print("Error in weights = ", errzw, ", errors in first, second and third integral are ", ae1, ", ", ae2, " and ", ae3, '\n')


text:
## SMA performance: MarketTechnicals.jl TimeArray version
using TimeSeries, MarketTechnicals
using Reactive

nma = 50
ohlc_BA = TimeSeries.readtimearray("../../data/OHLC_BA.csv")
sma_ta = sma(ohlc_BA["Close"], nma)

s_ohlc = Reactive.Input(ohlc_BA[1:nma])
s_close = lift(s -> values(s["Close"])[end], s_ohlc, typ=Float64)
s_sma = lift(s -> values(sma(s["Close"], nma))[end], s_ohlc, typ=Float64)

vsig = [s_sma.value]
function vsigupd!(vsig::Vector{Float64}, nma::Int64, nmax::Int64)
    for i in (nma + 1):nmax
        push!(s_ohlc, ohlc_BA[i-nma:i])
        push!(vsig, s_sma.value)
    end
end

println("SMA MarketTechnicals TimeArray version:")
println("reactive signal updates timing")
expr_run = parse("@time vsigupd!(vsig, nma, length(ohlc_BA))")
for i = 1:5
    eval(expr_run)
end

### results
# elapsed time: 1.771798644 seconds (577773184 bytes allocated, 29.26% gc time)
### way too long, partially due to inefficient TimeArray-slicing updates (see ../Reactive_push_lift_ohlc/)

Profile.clear()  # in case we have any previous profiling data
@profile vsigupd!(vsig, nma, length(ohlc_BA))
Profile.print()


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: sso
using AWS.Compat
using AWS.UUIDs

"""
    get_role_credentials(account_id, role_name, x-amz-sso_bearer_token)
    get_role_credentials(account_id, role_name, x-amz-sso_bearer_token, params::Dict{String,<:Any})

Returns the STS short-term credentials for a given role name that is assigned to the user.

# Arguments
- `account_id`: The identifier for the AWS account that is assigned to the user.
- `role_name`: The friendly name of the role that is assigned to the user.
- `x-amz-sso_bearer_token`: The token issued by the CreateToken API call. For more
  information, see CreateToken in the AWS SSO OIDC API Reference Guide.

"""
function get_role_credentials(
    account_id,
    role_name,
    x_amz_sso_bearer_token;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return sso(
        "GET",
        "/federation/credentials",
        Dict{String,Any}(
            "account_id" => account_id,
            "role_name" => role_name,
            "headers" =>
                Dict{String,Any}("x-amz-sso_bearer_token" => x_amz_sso_bearer_token),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_role_credentials(
    account_id,
    role_name,
    x_amz_sso_bearer_token,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return sso(
        "GET",
        "/federation/credentials",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "account_id" => account_id,
                    "role_name" => role_name,
                    "headers" => Dict{String,Any}(
                        "x-amz-sso_bearer_token" => x_amz_sso_bearer_token
                    ),
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_account_roles(account_id, x-amz-sso_bearer_token)
    list_account_roles(account_id, x-amz-sso_bearer_token, params::Dict{String,<:Any})

Lists all roles that are assigned to the user for a given AWS account.

# Arguments
- `account_id`: The identifier for the AWS account that is assigned to the user.
- `x-amz-sso_bearer_token`: The token issued by the CreateToken API call. For more
  information, see CreateToken in the AWS SSO OIDC API Reference Guide.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"max_result"`: The number of items that clients can request per page.
- `"next_token"`: The page token from the previous response output when you request
  subsequent pages.
"""
function list_account_roles(
    account_id, x_amz_sso_bearer_token; aws_config::AbstractAWSConfig=global_aws_config()
)
    return sso(
        "GET",
        "/assignment/roles",
        Dict{String,Any}(
            "account_id" => account_id,
            "headers" =>
                Dict{String,Any}("x-amz-sso_bearer_token" => x_amz_sso_bearer_token),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_account_roles(
    account_id,
    x_amz_sso_bearer_token,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return sso(
        "GET",
        "/assignment/roles",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "account_id" => account_id,
                    "headers" => Dict{String,Any}(
                        "x-amz-sso_bearer_token" => x_amz_sso_bearer_token
                    ),
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_accounts(x-amz-sso_bearer_token)
    list_accounts(x-amz-sso_bearer_token, params::Dict{String,<:Any})

Lists all AWS accounts assigned to the user. These AWS accounts are assigned by the
administrator of the account. For more information, see Assign User Access in the AWS SSO
User Guide. This operation returns a paginated response.

# Arguments
- `x-amz-sso_bearer_token`: The token issued by the CreateToken API call. For more
  information, see CreateToken in the AWS SSO OIDC API Reference Guide.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"max_result"`: This is the number of items clients can request per page.
- `"next_token"`: (Optional) When requesting subsequent pages, this is the page token from
  the previous response output.
"""
function list_accounts(
    x_amz_sso_bearer_token; aws_config::AbstractAWSConfig=global_aws_config()
)
    return sso(
        "GET",
        "/assignment/accounts",
        Dict{String,Any}(
            "headers" =>
                Dict{String,Any}("x-amz-sso_bearer_token" => x_amz_sso_bearer_token),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_accounts(
    x_amz_sso_bearer_token,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return sso(
        "GET",
        "/assignment/accounts",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "headers" => Dict{String,Any}(
                        "x-amz-sso_bearer_token" => x_amz_sso_bearer_token
                    ),
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    logout(x-amz-sso_bearer_token)
    logout(x-amz-sso_bearer_token, params::Dict{String,<:Any})

Removes the client- and server-side session that is associated with the user.

# Arguments
- `x-amz-sso_bearer_token`: The token issued by the CreateToken API call. For more
  information, see CreateToken in the AWS SSO OIDC API Reference Guide.

"""
function logout(x_amz_sso_bearer_token; aws_config::AbstractAWSConfig=global_aws_config())
    return sso(
        "POST",
        "/logout",
        Dict{String,Any}(
            "headers" =>
                Dict{String,Any}("x-amz-sso_bearer_token" => x_amz_sso_bearer_token),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function logout(
    x_amz_sso_bearer_token,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return sso(
        "POST",
        "/logout",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "headers" => Dict{String,Any}(
                        "x-amz-sso_bearer_token" => x_amz_sso_bearer_token
                    ),
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
module HyperHessians

using SIMD: SIMD, Vec, vstore
using CommonSubexpressions: cse

include("rules.jl")

# This is currently "artifically" restricted to "square" HyperDuals where
# ϵ1 and ϵ2 have the same length.
struct HyperDual{N, T} <: Real
    v::T
    ϵ1::Vec{N,T} # Vec{M,T}
    ϵ2::Vec{N,T}
    ϵ12::NTuple{N, Vec{N,T}} # NTuple{M, Vec{N,T}}
end
HyperDual(v::T, ϵ1::Vec{N,T}, ϵ2::Vec{N,T}) where {N, T} = HyperDual(v, ϵ1, ϵ2, ntuple(i -> zero(Vec{N,T}), Val(N)))
HyperDual{N}(v::T) where {N, T} = HyperDual(v, zero(Vec{N,T}),zero(Vec{N,T}))
HyperDual{N,T}(v) where {N, T} = HyperDual{N}(T(v))

function HyperDual(v::T1, ϵ1::Vec{N,T2}, ϵ2::Vec{N,T2}, ϵ12::NTuple{N, Vec{N,T2}}) where {N, T1, T2}
    T = promote_type(T1, T2)
    HyperDual(T(v), convert(Vec{N,T}, ϵ1), convert(Vec{N,T}, ϵ2), convert.(Vec{N,T}, ϵ12))
end

Base.promote_rule(::Type{HyperDual{N,T1}}, ::Type{HyperDual{N,T2}}) where {N,T1,T2} = HyperDual{N, promote_type(T1, T2)}
Base.convert(::Type{HyperDual{N,T1}}, h::HyperDual{N,T2}) where {N,T1,T2} =
    HyperDual{N,T1}(T1(h.v), convert(Vec{N,T1}, h.ϵ1), convert(Vec{N,T1}, h.ϵ2), convert.(Vec{N,T1}, h.ϵ12))

# Make this look nicer
function Base.show(io::IO, h::HyperDual)
    print(io, h.v, " + ", Tuple(h.ϵ1), "ϵ1", " + ", Tuple(h.ϵ2), "ϵ2", " + ",  map(Tuple, h.ϵ12), "ϵ12")
end

Base.one(::Type{HyperDual{N, T}}) where {N, T} = HyperDual{N}(one(T))
Base.zero(::Type{HyperDual{N, T}}) where {N, T} = HyperDual{N}(zero(T))
Base.one(::HyperDual{N, T}) where {N, T} = one(HyperDual{N,T})
Base.zero(::HyperDual{N, T}) where {N, T} = zero(HyperDual{N,T})
Base.float(h::HyperDual{N,T}) where {N,T} = convert(HyperDual{N,float(T)}, h)

# Unary
@inline Base.:(-)(h::HyperDual{N}) where {N} = HyperDual(-h.v, -h.ϵ1, -h.ϵ2, ntuple(i -> -h.ϵ12[i], Val(N)))
@inline Base.:(+)(h::HyperDual) = h

# Binary
for f in (:+, :-)
    @eval begin
        @inline Base.$f(h1::HyperDual{N,T}, h2::HyperDual{N,T}) where {N,T} =
            HyperDual($f(h1.v, h2.v), $f(h1.ϵ1, h2.ϵ1), $f(h1.ϵ2, h2.ϵ2), ntuple(i -> $f(h1.ϵ12[i], h2.ϵ12[i]), Val(N)))
        @inline Base.$f(h1::HyperDual{N,T1}, h2::HyperDual{N,T2}) where {N,T1,T2} = $f(promote(h1, h2)...)
        @inline Base.$f(h::HyperDual{N}, r::Real) where {N} =
            HyperDual($f(h.v, r), h.ϵ1, h.ϵ2, h.ϵ12)
        @inline Base.$f(r::Real, h::HyperDual{N}) where {N} =
            HyperDual($f(r, h.v), h.ϵ1, h.ϵ2, h.ϵ12)
    end
end

for f in (:*, :/)
    @eval @inline Base.$f(h::HyperDual{N}, r::Real) where {N} =
        HyperDual($f(h.v, r), $f(h.ϵ1, r), $f(h.ϵ2, r), ntuple(i -> $f(h.ϵ12[i], r), Val(N)))
    @eval @inline Base.$f(r::Real, h::HyperDual{N}) where {N} =
        HyperDual($f(r, h.v), $f(r, h.ϵ1), $f(r, h.ϵ2), ntuple(i -> $f(r, h.ϵ12[i]), Val(N)))
end

@inline Base.:(/)(h1::HyperDual{N,T}, h2::HyperDual{N,T}) where {N,T} = h1 * inv(h2)

# Get's rid of a bunch of boundserror throwing in the LLVM IR, not sure it matters for runtime performance...
unsafe_getindex(v::SIMD.Vec, i::SIMD.IntegerTypes) = SIMD.Intrinsics.extractelement(v.data, i-1)
@inline ⊗(v1::Vec{N}, v2::Vec{N}) where {N} = ntuple(i -> unsafe_getindex(v1, i) * v2, Val(N))

@inline Base.:(*)(h1::HyperDual{N,T1}, h2::HyperDual{N,T2}) where {N,T1,T2} = *(promote(h1, h2)...)
@inline function Base.:(*)(h1::HyperDual{N,T}, h2::HyperDual{N,T}) where {N,T}
    r = h1.v * h2.v
    ϵ1 = muladd(h1.v, h2.ϵ1, h1.ϵ1 * h2.v)
    ϵ2 = muladd(h1.v, h2.ϵ2, h1.ϵ2 * h2.v)
    ϵ12_1 = h1.ϵ1 ⊗ h2.ϵ2
    ϵ12_2 = h2.ϵ1 ⊗ h1.ϵ2
    ϵ12 = ntuple(i -> muladd(h1.v, h2.ϵ12[i], muladd(h1.ϵ12[i], h2.v, ϵ12_1[i] + ϵ12_2[i])), Val(N))
    return HyperDual(r, ϵ1, ϵ2, ϵ12)
end
@inline Base.literal_pow(::typeof(^), x::HyperDual, ::Val{0}) = one(typeof(x))
@inline Base.literal_pow(::typeof(^), x::HyperDual, ::Val{1}) = x
@inline Base.literal_pow(::typeof(^), x::HyperDual, ::Val{2}) = x*x
@inline Base.literal_pow(::typeof(^), x::HyperDual, ::Val{3}) = x*x*x

for (f, f′, f′′) in DIFF_RULES
    expr = quote
        # Verify that the cse still works properly when changing this.
        v = $f(x)
        f′ = $f′
        f′′ = $f′′
        x23 = (f′′ * h.ϵ1) ⊗ h.ϵ2
        return HyperDual(v, h.ϵ1 * f′, h.ϵ2 * f′, ntuple(i -> h.ϵ12[i]*f′ + x23[i], Val(N)))
    end
    cse_expr = cse(expr; warn=false)
    @eval @inline function Base.$f(h::HyperDual{N,T}) where {N,T}
        x = h.v
        $cse_expr
    end
end

@inline function Base.sin(h::HyperDual{N}) where {N}
    s, c = sincos(h.v)
    f′, f′′ = c, -s
    x23 = (f′′ * h.ϵ1) ⊗ h.ϵ2
    return HyperDual(s, h.ϵ1 * f′, h.ϵ2 * f′, ntuple(i -> h.ϵ12[i]*f′ + x23[i], Val(N)))
end

@inline function Base.cos(h::HyperDual{N}) where {N}
    s, c = sincos(h.v)
    f′,  f′′ = -s, -c
    x23 = (f′′ * h.ϵ1) ⊗ h.ϵ2
    return HyperDual(c, h.ϵ1 * f′,  h.ϵ2 * f′, ntuple(i -> h.ϵ12[i]*f′ + x23[i], Val(N)))
end

##########
# Chunks #
##########

const DEFAULT_CHUNK_THRESHOLD = 8

struct Chunk{N} end

const CHUNKS = [Chunk{i}() for i in 1:DEFAULT_CHUNK_THRESHOLD]

function Chunk(input_length::Integer)
    N = pickchunksize(input_length)
    0 < N <= DEFAULT_CHUNK_THRESHOLD && return CHUNKS[N] # TODO: Check if this matters
    return Chunk{N}()
end

Chunk(x::AbstractArray) = Chunk(length(x))

# TODO: MOre testing needed here.
function pickchunksize(input_length)
    input_length == 1 && return 1
    input_length == 2 && return 2
    3 <= input_length <= 4 && return 4 # 3 seems to be very bad for some reason
    return 8
end

(chunksize(::Chunk{N})::Int) where {N} = N

##################
# Hessian Config #
##################

abstract type AbstractHessianConfig end

struct HessianConfig{D <: AbstractVector{<:HyperDual}, S} <: AbstractHessianConfig
    duals::D
    seeds::S
end
(chunksize(cfg::HessianConfig)::Int) = length(cfg.seeds)

function HessianConfig(x::AbstractVector{T}, chunk=Chunk(x)::Chunk) where {T}
    N = chunksize(chunk)
    @assert 0 < N
    duals = similar(x, HyperDual{N,T}) # not Vector
    seeds = collect(construct_seeds(NTuple{N,T}))
    return HessianConfig(duals, seeds)
end

struct HessianConfigThreaded{D <: AbstractVector{<:HyperDual}, S} <: AbstractHessianConfig
    duals::Vector{D}
    seeds::S
end

function HessianConfigThreaded(x::AbstractVector{T}, chunk=Chunk(x)::Chunk) where {T}
    N = chunksize(chunk)
    @assert 0 < N
    duals = [similar(x, HyperDual{N,T}) for i in 1:Threads.nthreads()] # not Vector
    seeds = collect(construct_seeds(NTuple{N,T}))
    return HessianConfigThreaded(duals, seeds)
end

(chunksize(cfg::HessianConfigThreaded)::Int) = length(cfg.seeds)
HessianConfig(cfg::HessianConfigThreaded) = HessianConfig(first(cfg.duals), cfg.seeds)

###########
# Seeding #
###########

@generated function single_seed(::Type{NTuple{N,T}}, ::Val{i}) where {N,T,i}
    ex = Expr(:tuple, [ifelse(i === j, :(one(T)), :(zero(T))) for j in 1:N]...)
    return :(Vec($(ex)))
end

@generated construct_seeds(::Type{NTuple{N,T}}) where {N,T} =
    Expr(:tuple, [:(single_seed(NTuple{N,T}, Val{$i}())) for i in 1:N]...)

seed_epsilon_1(d::HyperDual{N,T}, ϵ1) where {N,T} = HyperDual{N,T}(d.v, ϵ1, d.ϵ2, d.ϵ12)
seed_epsilon_2(d::HyperDual{N,T}, ϵ2) where {N,T} = HyperDual{N,T}(d.v, d.ϵ1, ϵ2, d.ϵ12)

function seed!(d::AbstractVector{<:HyperDual{N}}, x, seeds, block_i, block_j) where {N}
    d .= HyperDual{N}.(x)
    index_i = (block_i-1) * N + 1
    index_j = (block_j-1) * N + 1
    range_i = index_i : min(length(x), (index_i + N -1))
    range_j = index_j : min(length(x), (index_j + N -1))
    chunks_i = length(range_i)
    chunks_j = length(range_j)

    d[range_i] .= seed_epsilon_1.(view(d, range_i), view(seeds, 1:chunks_i))
    d[range_j] .= seed_epsilon_2.(view(d, range_j), view(seeds, 1:chunks_j))
    return d
end

# Hessian
@noinline check_scalar(x) =
    x isa Number || throw(error("expected a scalar to be returned from function passed to `hessian`"))

function extract_hessian!(H::AbstractMatrix, v::HyperDual)
    Base.require_one_based_indexing(H)
    @inbounds for i in 1:size(H, 2)
        H[:, i] .= Tuple(v.ϵ12[i])
    end
    return H
end

function symmetrize!(H::AbstractMatrix)
    Base.require_one_based_indexing(H)
    for i in 1:size(H, 1)
        for j in i:size(H,2)
            H[j,i] = H[i,j]
        end
    end
    return H
end

function extract_hessian!(H::AbstractMatrix, v::HyperDual{N}, block_i::Int, block_j::Int) where {N}
    index_i = (block_i-1) * N + 1
    index_j = (block_j-1) * N + 1
    range_i = index_i : (index_i + N -1)
    range_j = index_j : (index_j + N -1)

    for (I, i) in enumerate(range_i)
        for (J, j) in enumerate(range_j)
            if checkbounds(Bool, H, i, j)
                H[i, j] = v.ϵ12[I][J]
            end
        end
    end
    return H
end

function extract_gradient!(G::AbstractVector, v::HyperDual, block_i::Int)
    index_i = (block_i-1) * N + 1
    range_i = index_i : (index_i + N -1)
    for (I, i) in enumerate(range_i)
        G[i] .= v.ϵ1[I]
    end
    return G
end

###############
# Hessian API #
###############

# Scalar
function hessian(f, x::Real)
    dual = HyperDual(x, Vec(one(x)), Vec(one(x)))
    v = f(dual)
    check_scalar(v)
    return @inbounds v.ϵ12[1][1]
end

hessian(f, x::AbstractVector)                             = hessian!(similar(x, axes(x,1), axes(x,1)), f, x, HessianConfig(x))
hessian(f, x::AbstractVector, cfg::AbstractHessianConfig) = hessian!(similar(x, axes(x,1), axes(x,1)), f, x, cfg)

function hessian!(H::AbstractMatrix, f, x::AbstractVector{T}, cfg::AbstractHessianConfig) where {T,N}
    if chunksize(cfg) == length(x)
        if cfg isa HessianConfigThreaded
            cfg = HessianConfig(cfg)
        end
        return hessian_vector!(H, f, x, cfg)
    else
        return hessian_chunk!(H, f, x, cfg)
    end
end

function hessian_vector!(H::AbstractMatrix, f, x::AbstractVector, cfg::HessianConfig)
    @assert size(H,1) == size(H,2) == length(x)
    cfg.duals .= HyperDual.(x, cfg.seeds, cfg.seeds)
    v = f(cfg.duals)
    check_scalar(v)
    return extract_hessian!(H, v)
end

function hessian_chunk!(H::AbstractMatrix, f, x::AbstractVector{T}, cfg::HessianConfig) where {T,N}
    @assert size(H,1) == size(H,2) == length(x)
    n_chunks = ceil(Int, length(x) / chunksize(cfg))
    for i in 1:n_chunks
        for j in i:n_chunks
            seed!(cfg.duals, x, cfg.seeds, i, j)
            v = f(cfg.duals)
            check_scalar(v)
            extract_hessian!(H, v, i, j)
        end
    end
    symmetrize!(H)
    return H
end

# Threaded chunk mode
function linear_to_cartesian(n::Int, k::Int)
    k -= 1
    n += 1
    i = n - 2 - floor(sqrt(-8*k + 4*n*(n-1)-7)/2 - 0.5)
    j = k + i + 1 - n*(n-1)/2 + (n-i)*((n-i)-1)/2
    return round(Int, i+1), round(Int, j)
end

function hessian_chunk!(H::AbstractMatrix, f, x::AbstractVector{T}, cfg::HessianConfigThreaded) where {T,N}
    @assert size(H,1) == size(H,2) == length(x)
    n_chunks = ceil(Int, length(x) / chunksize(cfg))
    Threads.@threads for k in 1:((n_chunks*(n_chunks+1)) ÷ 2)
        i, j = linear_to_cartesian(n_chunks, k)
        duals = cfg.duals[Threads.threadid()]
        seed!(duals, x, cfg.seeds, i, j)
        v = f(duals)
        check_scalar(v)
        extract_hessian!(H, v, i, j)
    end
    symmetrize!(H)
    return H
end

end # module


text:
module MixedModelsSim

using LinearAlgebra
using MixedModels
using PooledArrays
using PrettyTables
using Random
using Statistics
using Tables

using MixedModels: replicate

export
    create_re,
    create_theta,
    createθ,
    cyclicshift,
    factorproduct,
    flatlowertri,
    nlevels,
    nlevstbl,
    #withinitem,
    pooled!,
    power_table,
    simdat_crossed,
    update!

export pretty_table, @pt # re-exports

include("utilities.jl")
include("power.jl")
include("simdat.jl")

end # module


text:
# Submodule to compute the coefficients of shearlet in serial in 2D and its reconstruction

# Type for the shearlet serial prepared filters

struct Serialpreparedfilters
		Xfreq
		Xrec
		preparedFilters
		dualFrameWeightsCurr
		shearletIdxs
end

##############################################################################
## Function that prepare the serial shearlet transform

"""
...
prepareserial2D(X,nScales,shearLevels = ceil.((1:nScales)/2), full = 0,
directionalFilter = filt_gen("directional_shearlet") , quadratureMirrorFilter = filt_gen("scaling_shearlet"))
...
"""
function prepareserial2D(X,nScales,shearLevels = ceil.((1:nScales)/2), full = 0,
directionalFilter = filt_gen("directional_shearlet") , quadratureMirrorFilter = filt_gen("scaling_shearlet"))
		shearLeveles = map(Int,shearLevels)
		# Compute the fft transform
		Xfreq = fftshift(fft(ifftshift(X)));
		# Compute prepared filters
		preparedFilters = preparefilters2D(size(X,1),size(X,2),nScales,shearLevels,directionalFilter,quadratureMirrorFilter)
		# Compute dual frames
		dualFrameWeightsCurr = zeros(size(X));
		Xrec = zeros(size(X));
		# Compute the indices
		shearletIdxs = getshearletidxs2D(shearLevels,full);

		# Assign values to the type
		return Serialpreparedfilters(Xfreq, Xrec, preparedFilters, dualFrameWeightsCurr, shearletIdxs)
end #prepareserial2D

# Type for the 2D serial sheralet decompostion
struct Sheardecserial2D
		coeffs
		shearlet
		dualFrameWeightsNew
		RMS
end

##############################################################################
## Function that decompose a signal in frequency domain in serial
"""
...
sheardecserial2D(Xfreq, ShearletIdx, preparedFilters, dualFrameWeightsCurr) decompose a signal in frequency domain in serial
...
"""
function sheardecserial2D(Xfreq, shearletIdx, preparedFilters, dualFrameWeightsCurr)
		# Compute the serial decomposition
		Shearlets2D = Shearlab.getshearlets2D(preparedFilters,shearletIdx);
		shearlet = Shearlets2D.shearlets;
		RMS = Shearlets2D.RMS;
		dualFrameWeights = Shearlets2D.dualFrameWeights;
		dualFrameWeightsNew = dualFrameWeightsCurr+dualFrameWeights;
		coeffs = fftshift(ifft(ifftshift(conj(shearlet).*Xfreq)));
		return Sheardecserial2D(coeffs, shearlet, dualFrameWeightsNew, RMS)
end # sheardecserial2D

##############################################################################
## Function that makes the serial recovery
"""
...
shearrecserial2D(coeffs, shearlet, Xrec) serial recovery of coefficients with certain shearlet
...
"""
function shearrecserial2D(coeffs, shearlet, Xrec)
		return Xrec+fftshift(fft(ifftshift(coeffs))).*shearlet;
end # shearrecserial2D

##############################################################################
## Function that finishes the serial recovery
"""
...
finishserial2D(Xrec,dualFrameWeightsCurr) serial recovery of coefficients with certain shearlet
...
"""
function finishserial2D(Xcurr, dualFrameWeights)
		return fftshift(ifft(ifftshift((1 ./dualFrameWeights).*Xcurr)))
end # shearrecserial2D


text:
using Random
import Random.rand

function rand(rng::AbstractRNG, ::Random.SamplerTrivial{Random.CloseOpen01{DoubleFloat{T}}}) where {T<:IEEEFloat}
    hi, lo  = rand(rng, T, 2)
    while hi == zero(T) # fairly common with Float16 or Float32
        hi = rand(rng, T)
    end
    frlo, xplo  = frexp(lo)
    xplo = Base.exponent(hi) - min(1, fld(xplo,4)) - abs(Base.exponent(eps(hi)))
    lo = ldexp(frlo, xplo)
    lo = rand(Bool) ? lo : -lo

    DoubleFloat(hi, lo)
end


text:
using Cxx
include(Pkg.dir("Cxx","test","llvmincludes.jl"))
function collectSymbolsForExport(RD::pcpp"clang::CXXRecordDecl")
    C = Cxx.instance(Cxx.__default_compiler__)
    f = Cxx.CreateFunctionWithPersonality(C, Cxx.julia_to_llvm(Void),[Cxx.julia_to_llvm(Ptr{Ptr{Void}})])
    state = Cxx.setup_cpp_env(C,f)
    builder = Cxx.irbuilder(C)
    names = Symbol[]
    nummethods = icxx"""
    auto *CGM = $(C.CGM);
    auto *RD = clang::cast<clang::CXXRecordDecl>($RD->getDefinition());
    assert(RD);
    RD->dump();
    CGM->EmitTopLevelDecl(RD);
    llvm::SmallVector<llvm::Constant *, 0> addresses;
    unsigned i = 0;
    clang::CodeGen::CGBuilderTy *builder = $builder;
    for (auto method : RD->methods()) {
        clang::GlobalDecl GD;
        if (clang::isa<clang::CXXConstructorDecl>(method))
            GD = clang::GlobalDecl(clang::cast<clang::CXXConstructorDecl>(method),clang::Ctor_Complete);
        else if (clang::isa<clang::CXXDestructorDecl>(method))
            GD = clang::GlobalDecl(clang::cast<clang::CXXDestructorDecl>(method),clang::Dtor_Complete);
        else
            GD = clang::GlobalDecl(method);
        auto name = CGM->getMangledName(GD);
        $:(push!(names,symbol(bytestring(icxx"return name.data();"))));
        auto llvmf = CGM->GetAddrOfGlobal(GD,false);
        if (clang::cast<llvm::Function>(llvmf)->isDeclaration())
            continue;
        llvm::Value *Val =
          builder->CreateBitCast(llvmf,$(Cxx.julia_to_llvm(Ptr{Void})));
        llvm::Value *Addr = builder->CreateConstGEP1_32($(Cxx.julia_to_llvm(Ptr{Void})),
          &$f->getArgumentList().front(), i++);
        builder->CreateStore(Val,
          clang::CodeGen::Address(Addr,clang::CharUnits(8)));
    }
    builder->CreateRetVoid();
    i;
    """
    Cxx.cleanup_cpp_env(C,state)
    addresses = Array(Ptr{Void},nummethods)
    eval(:(Core.Intrinsics.llvmcall($(f.ptr),Void,Tuple{Ptr{Ptr{Void}}},$(pointer(addresses)))))
    symbols = Dict{Symbol,Ptr{Void}}()
    for (k,v) in zip(names,addresses)
        symbols[k] = v
    end
    symbols
end

export_class(RD::pcpp"clang::CXXRecordDecl") = export_symbols(collectSymbolsForExport(RD))


text:
using ProgressMeter

function couplingTime(model::SMCModel, lM::F, ccsmcio::CCSMCIO{Particle},
  algorithm::Symbol = :BS, independentInitialization::Bool = false,
  rngCouple::Bool = false, maxit::Int64 = typemax(Int64)) where
  {F<:Function, Particle}

  initializeCCSMC(model, lM, ccsmcio, algorithm, independentInitialization)

  ref1 = ccsmcio.ref1
  ref2 = ccsmcio.ref2

  for i in 1:maxit
    ccXpf!(model, lM, ccsmcio, algorithm, rngCouple)
    checkEqual(ref1, ref2) && return i
  end
  return 0
end

function couplingTime(model::SMCModel, ccsmcio::CCSMCIO{Particle},
  independentInitialization::Bool=false, rngCouple::Bool = false,
  maxit::Int64 = typemax(Int64)) where Particle
  return couplingTime(model, error, ccsmcio, :AT, independentInitialization,
    rngCouple, maxit)
end

function couplingTimes(model::SMCModel, lM::F, N::Int64, n::Int64,
  m::Int64, algorithm::Symbol = :BS, independentInitialization::Bool=false,
  rngCouple::Bool = false, maxit::Int64 = typemax(Int64)) where
  {F<:Function, Particle}
  nt = Threads.nthreads()
  ccsmcios::Vector{CCSMCIO{model.particle, model.pScratch}} =
    Vector{CCSMCIO{model.particle, model.pScratch}}(undef, nt)
  Threads.@threads for i in 1:nt
    ccsmcios[i] = CCSMCIO{model.particle, model.pScratch}(N, n)
  end
  vs::Vector{Int64} = Vector{Int64}(undef, m)
  p = Progress(div(m, nt), 10)
  abort = false
  Threads.@threads for i in 1:m
    tid = Threads.threadid()
    abort && break
    vs[i] = couplingTime(model, lM, ccsmcios[tid], algorithm,
      independentInitialization, rngCouple, maxit)
    if vs[i] == 0
      abort = true
    end
    tid == 1 && update!(p, i)
  end
  if abort
    return [missing]
  else
    return vs
  end
end

function couplingTimes(model::SMCModel, N::Int64, n::Int64, m::Int64,
  independentInitialization::Bool = false, rngCouple::Bool = false,
  maxit::Int64 = typemax(Int64)) where Particle
  return couplingTimes(model, error, N, n, m, :AT, independentInitialization,
    rngCouple, maxit)
end


text:
module Obfs

using Dates
using Sockets

import Sockets.accept
import Base.read
import Base.write
import Base.readbytes!

using Crypto.Chacha20: Chacha20Encrypt, Chacha20Decrypt
using Crypto.HKDF: hkdf
using ..Common: ioCopy

import ..Common.read
import ..Common.write

const INFO = Array{UInt8}("obfs-tls")
const TIMEOUT = Millisecond(1000)
const ZERO = Millisecond(0)
const MaxSize = 0x5000

mutable struct TLSConfig
    host::Union{String, Sockets.IPAddr}
    port::UInt16
    key::Array{UInt8}
end

mutable struct TLSConnection
    conn::TCPSocket
    cache::IOStream
end

@inline function eof(io::TLSConnection)
    eof(io.conn)
end

@inline function isopen(io::TLSConnection)
    isopen(io.conn)
end

@inline function readbytes!(io::TLSConnection, buff::Array{UInt8}, bytes::Integer)
    readbytes!(io.conn, buff, nbytes)
end

@inline function write(io::TLSConnection, buff::Array{UInt8})
    write(io.conn, buff)
end

@inline function read(io::TLSConnection, buff::Array{UInt8})
    try
        eof(io)
    catch err
        return nothing, err
    end

    isopen(io) ? try
        readbytes(io, buff, 0x05)
    catch err
        return nothing, err
    end : return nothing, Error("Connection Closed")

    nbytes = UInt16(buff[4]) * 0x0100 + buff[5]

    isopen(io) ? try
        readbytes!(io, buff, nbytes)
    catch err
        return nothing, err
    end : return nothing, Error("Connection Closed")

    return nbytes, nothing
end

@inline function read(io::TLSConnection, buff::Array{UInt8}, nbytes::Integer)
    while true
        if nbytes <= bytesavailable(io.cache)
            readbytes!(io.cache, buff, nbytes)
            break
        end

        nread, err = read(io, buff)
        if err != nothing
            return nothing, err
        end

        write(io.cache, unsafe_wrap(Array{UInt8}, pointer(buff), nread))
    end

    return nbytes, nothing
end

@inline function write(io::TLSConnection, buff::Array{UInt8}, nbytes::Integer)
    isopen(io) ? try
        write(io, [[0x17, 0x03, 0x03, ]; UInt8[nbytes >> 8, nbytes & 0xff, ]])
        write(io, unsafe_wrap(Array{UInt8}, pointer(buff), nbytes))
    catch err
        return err
    end : return Error("Connection Closed")

    return nothing
end

mutable struct RandomSet
    lock::ReentrantLock
    set::Set{Array{UInt8}}
end

struct Error <: Exception
    msg::AbstractString
end

function add(random::Array{UInt8}, set::RandomSet)
    lock(random.lock)
    union!(set, [random, ])
    unlock(random.lock)
end

function clear(set::RandomSet)
    lock(random.lock)
    set.set = Set{Array{UInt8}}()
    unlock(random.lock)
end

function hasRandom(random::Array{UInt8}, set::RandomSet)::Bool
    lock(random.lock)
    x = random in set.set
    unlock(random.lock)

    x
end

function genRandom(config::TLSConfig)::Array{UInt8}
    random = Array{UInt8}(undef, 32)
    random[24:32] = rand(UInt8, 9)
    kv, err = hkdf("SHA2-256", random[24:32], config.key, INFO, 44)
    Chacha20Encrypt(random, kv[1:32], kv[33:44], Array{UInt8}(string(Dates.now())))

    random
end

function checkRandom(random::Array{UInt8}, config::TLSConfig)::Bool
    x = Array{UInt8}(undef, 23)
    kv, err = hkdf("SHA2-256", random[24:32], config.key, INFO, 44)
    Chacha20Decrypt(x, kv[1:32], kv[33:44], random[1:23])

    now = try
        Dates.now() - parse(Dates.DateTime, String(x))
    catch err
        return false
    end

    if ZERO < now < TIMEOUT
        true
    else
        false
    end
end

function parseClientHello(buff::Array{UInt8})::Tuple{Union{Array{UInt8}, Nothing}, Union{Array{UInt8}, Nothing}, Union{Exception, Nothing}}
    while true
        pointer = 6;
        if buff[pointer] != 0x01
            break
        end

        pointer += 6;
        random = buff[pointer:pointer+32]

        pointer += 33;
        sessionId = buff[pointer:pointer+32]

        return random, sessionId, nothing
    end

    return nothing, nothing, Error("Not a Valid Client Hello Message")
end

function composeExtension()
    [0x00, ]
end

function composeClientHello(config::TLSConfig)::Array{UInt8}
    [[0x16, ]; [0x03, 0x03, ]; [0x00, 0x00];
        [
            [0x01, ];
            [0x00, 0x01, 0xfc];
            [0x03, 0x03];
            genRandom(config);
            [0x20, ];
            rand(UInt8, 0x20);
            [0x00, 0x1c];
            [0x2a, 0x2a, 0xc0, 0x2b, 0xc0, 0x2f, 0xc0, 0x2c, 0xc0, 0x30, 0xcc, 0xa9, 0xcc, 0xa8, 0xc0, 0x13, 0xc0, 0x14, 0x00, 0x9c, 0x00, 0x9d, 0x00, 0x2f, 0x00, 0x35, 0x00, 0x0a, ];
            [0x01, ];
            [0x00, ];
            [0x01, 0x97, ];
            composeExtension()
        ]
    ]
end

function composeServerHello(sessionId::Array{UInt8})::Array{UInt8}
    [[0x16, ]; [0x03, 0x03, ]; [];
        [
            [0x02, ];
            [0x00, 0x00, 0x4d, ];
            [0x03, 0x03];
            rand(UInt8, 0x20);
            [0x20, ];
            sessionId;
            [0xc0, 0x30, ];
            [0x00, ];
            [0x00, 0x05, ];
            [0xff, 0x01, 0x00, 0x01, 0x00, ]
        ]
    ]
end

function composeChangeCipher()
    [[0x14, ]; [0x03, 0x03, ]; [0x00, 0x01, ]; [0x01, ]]
end

function composeFinished()
    [[0x16, ]; [0x03, 0x03, ]; [0x00, 0x40, ]; rand(UInt8, 0x28); Array{UInt8}(undef, 0x18)]
end

function accept(server::Sockets.TCPServer, config::TLSConfig)::TCPSocket
    if @isdefined set

    else
        global set = RandomSet(ReentrantLock(), Set{Array{UInt8}}())
        @async while true
            sleep(3600)
            clear(set)
        end
    end

    buff = Array{UInt8}(undef, 1024)

    while true
        conn = try
            accept(server)
        catch err
            throw(err)
        end

        nbytes, err = read(conn, buff)
        if err != nothing
            close(conn)
            continue
        end

        random, sessionId, err = parseClientHello(buff)
        if err != nothing
            close(conn)
            continue
        end

        if checkRandom(random, config) && !hasRandom(random, set)
            add(random, set)
            read(conn, buff, 75)
            return conn
        else
            fake = try
                connect(config.host, config.port)
            catch err
            end

            err = write(fake, buff, nbytes)
            if err != nothing
            end

            @async ioCopy(conn, fake)
            @async ioCopy(fake, conn)
        end
    end
end

end # module

if false

using Shadowsocks
using Sockets

config = Shadowsocks.Obfs.TLSConfig(IPv4(13,35,15,164), 0x01bb, Array{UInt8}("julia"))
server = listen(2000)

Shadowsocks.Obfs.accept(server, config)

end


text:
# coupeler_driver

# import packages
using Pkg
import SciMLBase: step!
using OrdinaryDiffEq
using OrdinaryDiffEq: ODEProblem, solve, SSPRK33, savevalues!, Euler
using LinearAlgebra
import Test: @test

# import coupler utils
include("coupler_utils/flux_calculator.jl")
include("coupler_utils/conservation_checker.jl")
include("coupler_utils/regridder.jl")

# initiate spatial and temporal info
t_end = 2592000 # 100e2
tspan = (0, t_end)
Δt_cpl = 2e2
saveat = Δt_cpl * 1000

# init MPI
include("mpi/mpi_init.jl")

# init model components
include("atmos/atmos_init.jl")
atmos_sim = atmos_init(FT, Y, spaces, integrator, params = params);

include("slab/slab_init.jl")
boundary_space = ClimaCore.Fields.level(atmos_sim.domain.center_space, 1)
slab_sim = slab_init(FT, tspan, dt = Δt_cpl, space = boundary_space, saveat = saveat);

# init boundary fields for regridding (TODO: technically this can be bypassed by directly rigridding on model grids)
T_S = ClimaCore.Fields.zeros(boundary_space)
F_S = ClimaCore.Fields.zeros(boundary_space)

# init conservation info collector
CS = ConservationCheck([], [])

# coupling loop
@show "Starting coupling loop"
walltime = @elapsed for t in (tspan[1]:Δt_cpl:tspan[end])
    @show t

    ## Atmos
    # calculate surface fluxes for Atmos BCs
    T_S .= FT(0)
    dummmy_remap!(T_S, slab_sim.integrator.u.T_sfc)

    #atmos_sim.integrator.p.dif_flux_energy .= ClimaCore.Geometry.WVector.(ClimaCore.Fields.zeros(axes(atmos_sim.integrator.p.dif_flux_energy)))
    parent(atmos_sim.integrator.p.dif_flux_energy) .= FT(0)
    calculate_surface_fluxes_atmos_grid!(atmos_sim.integrator, T_S)

    # run
    step!(atmos_sim.integrator, t - atmos_sim.integrator.t, true) # NOTE: use (t - integ_atm.t) here instead of Δt_cpl to avoid accumulating roundoff error in our timestepping.

    ## Slab
    # pre: get accumulated flux from atmos
    F_S .= ClimaCore.Fields.zeros(boundary_space)
    dummmy_remap!(F_S, atmos_sim.integrator.p.dif_flux_energy)

    # save the accumulated flux
    slab_F_sfc = slab_sim.integrator.p.F_sfc
    @. slab_F_sfc = -F_S

    # run
    step!(slab_sim.integrator, t - slab_sim.integrator.t, true)

    # conservation info "callback"
    if !is_distributed
        check_conservation_callback(CS, atmos_sim, slab_sim)
    end

end
@show walltime

@show "Postprocessing"
# collect solutions
sol_atm = atmos_sim.integrator.sol
sol_slab = slab_sim.integrator.sol

include("mpi/mpi_postprocess.jl")

# conservation  check
if !is_distributed || (is_distributed && ClimaComms.iamroot(comms_ctx))
    conservation_plot(atmos_sim, slab_sim, solu_atm, solu_slab, "conservation.png")
end

max_stars_repo_name: CliMA/ClimaCoupler


text:
# simple script to print the status (i.e. version) for
# all packages that we use

used_packages =  ["Distributions";
                  "DataFrames";
                  "StatsBase";
                  "Optim";
                  "Lasso";
                  "StatsFuns";
                  "JLD";
                  "HDF5";
	"Yeppp"]

for pack in used_packages
  Pkg.status(pack)
end

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text:
"""
Full type hierarchy intended to be mostly a 1:1 relationship with the functions listed in
the AlphaVantage documentation: https://www.alphavantage.co/documentation/.

Forex and Crypto have been grouped together under AVCurrency since they are very similar.
"""

# Abstracts
abstract type AVFunction end
abstract type AVTimeSeries <: AVFunction end
abstract type AVFundamental <: AVFunction end
abstract type AVCurrency <: AVFunction end
abstract type AVEconomicIndicator <: AVFunction end
abstract type AVTechnicalIndicator <: AVFunction end

# Time series
struct TIME_SERIES_INTRADAY <: AVTimeSeries end
struct TIME_SERIES_INTRADAY_EXTENDED <: AVTimeSeries end
struct TIME_SERIES_DAILY <: AVTimeSeries end
struct TIME_SERIES_DAILY_ADJUSTED <: AVTimeSeries end
struct TIME_SERIES_WEEKLY <: AVTimeSeries end
struct TIME_SERIES_WEEKLY_ADJUSTED <: AVTimeSeries end
struct TIME_SERIES_MONTHLY <: AVTimeSeries end
struct TIME_SERIES_MONTHLY_ADJUSTED <: AVTimeSeries end
struct GLOBAL_QUOTE <: AVTimeSeries end
struct SYMBOL_SEARCH <: AVTimeSeries end

# Fundamentals
struct INCOME_STATEMENT <: AVFundamental end
struct BALANCE_SHEET <: AVFundamental end
struct CASH_FLOW <: AVFundamental end
struct EARNINGS <: AVFundamental end
struct OVERVIEW <: AVFundamental end
struct LISTING_STATUS <: AVFundamental end
struct EARNINGS_CALENDAR <: AVFundamental end
struct IPO_CALENDAR <: AVFundamental end

# Forex / Crypto
struct CURRENCY_EXCHANGE_RATE <: AVCurrency end
struct FX_INTRADAY <: AVCurrency end
struct FX_DAILY <: AVCurrency end
struct FX_WEEKLY <: AVCurrency end
struct FX_MONTHLY <: AVCurrency end
struct AVC_PRICE <: AVCurrency end
struct CRYPTO_INTRADAY <: AVCurrency end
struct DIGITAL_CURRENCY_DAILY <: AVCurrency end
struct DIGITAL_CURRENCY_WEEKLY <: AVCurrency end
struct DIGITAL_CURRENCY_MONTHLY <: AVCurrency end

# Economic Indicators
struct REAL_GDP <: AVEconomicIndicator end
struct REAL_GDP_PER_CAPITA <: AVEconomicIndicator end
struct TREASURY_YIELD <: AVEconomicIndicator end
struct FEDERAL_FUNDS_RATE <: AVEconomicIndicator end
struct INFLATION <: AVEconomicIndicator end
struct INFLATION_EXPECTATION <: AVEconomicIndicator end
struct CONSUMER_SENTIMENT <: AVEconomicIndicator end
struct RETAIL_SALES <: AVEconomicIndicator end
struct DURABLES <: AVEconomicIndicator end
struct UNEMPLOYMENT <: AVEconomicIndicator end
struct NONFARM_PAYROLL <: AVEconomicIndicator end

# Technical Indicators
struct SMA <: AVTechnicalIndicator end
struct EMA <: AVTechnicalIndicator end
struct WMA <: AVTechnicalIndicator end
struct DEMA <: AVTechnicalIndicator end
struct TEMA <: AVTechnicalIndicator end
struct TRIMA <: AVTechnicalIndicator end
struct KAMA <: AVTechnicalIndicator end
struct MAMA <: AVTechnicalIndicator end
struct VWAP <: AVTechnicalIndicator end
struct T3 <: AVTechnicalIndicator end
struct MACD <: AVTechnicalIndicator end
struct MACDEXT <: AVTechnicalIndicator end
struct STOCH <: AVTechnicalIndicator end
struct STOCHF <: AVTechnicalIndicator end
struct RSI <: AVTechnicalIndicator end
struct STOCHRSI <: AVTechnicalIndicator end
struct WILLR <: AVTechnicalIndicator end
struct ADX <: AVTechnicalIndicator end
struct ADXR <: AVTechnicalIndicator end
struct APO <: AVTechnicalIndicator end
struct PPO <: AVTechnicalIndicator end
struct MOM <: AVTechnicalIndicator end
struct BOP <: AVTechnicalIndicator end
struct CCI <: AVTechnicalIndicator end
struct CMO <: AVTechnicalIndicator end
struct ROC <: AVTechnicalIndicator end
struct ROCR <: AVTechnicalIndicator end
struct AROON <: AVTechnicalIndicator end
struct AROONOSC <: AVTechnicalIndicator end
struct MFI <: AVTechnicalIndicator end
struct TRIX <: AVTechnicalIndicator end
struct ULTOSC <: AVTechnicalIndicator end
struct DX <: AVTechnicalIndicator end
struct MINUS_DI <: AVTechnicalIndicator end
struct PLUS_DI <: AVTechnicalIndicator end
struct MINUS_DM <: AVTechnicalIndicator end
struct PLUS_DM <: AVTechnicalIndicator end
struct BBANDS <: AVTechnicalIndicator end
struct MIDPOINT <: AVTechnicalIndicator end
struct MIDPRICE <: AVTechnicalIndicator end
struct SAR <: AVTechnicalIndicator end
struct TRANGE <: AVTechnicalIndicator end
struct ATR <: AVTechnicalIndicator end
struct NATR <: AVTechnicalIndicator end
struct AD <: AVTechnicalIndicator end
struct ADOSC <: AVTechnicalIndicator end
struct OBV <: AVTechnicalIndicator end
struct HT_TRENDLINE <: AVTechnicalIndicator end
struct HT_SINE <: AVTechnicalIndicator end
struct HT_TRENDMODE <: AVTechnicalIndicator end
struct HT_DCPERIOD <: AVTechnicalIndicator end
struct HT_DCPHASE <: AVTechnicalIndicator end
struct HT_PHASOR <: AVTechnicalIndicator end

Base.string(x::A where A <: AVFunction) = string(typeof(x))


text:
using CEDICT

using Test

@testset "unit tests" begin
    @testset "pinyin fuzzy matching" begin
        re = CEDICT._convert_pinyin_regex("jue2 dai4 shuang1 jiao1")
        @test match(re, "jue2 dai4 shuang1 jiao1") !== nothing
        @test match(re, "jue2 dai4 shuang1 jiao2") === nothing
        @test match(re, "wu2 jue2 dai4 shuang1 jiao1") === nothing
        @test match(re, "jue2 shuang1 jiao1") === nothing
        @test match(re, "jue2 dai4 shuang1 jiao1 ji4") === nothing
    end

    @testset "dictionary" begin
        dict = ChineseDictionary("res/small_dict.txt")

        @testset "loading" begin
            @test length(dict) == 744
        end

        @testset "contains" for key in ["代數拓撲", "做事", "優惠券", "优惠券"]
            @test haskey(dict, key)
        end

    end

    @testset "dictionary headword search" begin

    end

    @testset "dictionary sense/meaning search" begin

    end
end


text:
# This file was generated, do not modify it. # hide
LogisticClassifier = @load LogisticClassifier pkg=MLJLinearModels
X2 = select(X, Not([:Year, :Today]))
classif = machine(LogisticClassifier(), X2, y)

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
# Plot the Objective over the course of the solve

using TrajOptSOCPs, Plots

function plotObjective(obj::objectiveFunc, trajList)
    fList = [TrajOptSOCPs.fObjQP(obj, traj)[1] for traj in trajList]

    pltObj = plot(fList, markershape = :square)
    title!("Cost Function")
    ylabel!("Cost Function")
    xlabel!("Newton Step")

    return pltObj
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
module FloatTiming

using StaticArrays, DoubleFloats, Test, BenchmarkTools

widemult(a::Float64, b::Real) = Double64(a)*Double64(b)
widemult(a::Union{Float32,Int}, b::Float64) = Double64(a)*Double64(b)
widemult(a, b) = widemul(a,b)

""" det2(p,q, a,b)
2x2 same precision
"""
det2(p,q, a,b) = (p[a]*q[b]) - (p[b]*q[a])
@inline det2d(p,q, a,b) = (p[a]*(q[b]-p[b])) - (p[b]*(q[a]-p[a])) # more accurate for nearby soints.
det2f(p,q, a,b) = muladd(p[a],q[b], -(p[b]*q[a]))
det2w(p,q, a,b) = widemult(p[a],q[b]) - widemult(p[b],q[a]) # 2 mult, add, neg
det2wm(p,q, a,b) = muladd(p[a],q[b], -widemult(p[b],q[a])) # 2 mult, add, neg
@inline det2i(p,q, a,b) =  @inbounds (p[a]*q[b] - p[b]*q[a])
@inline det2fi(p,q, a,b) = @inbounds muladd(p[a],q[b], -(p[b]*q[a]))
@inline det2wi(p,q, a,b) = @inbounds widemult(p[a],q[b]) - widemult(p[b],q[a]) # 2 mult, add, neg
@inline det2mi(p,q, a,b) = @inbounds muladd(p[a],q[b], -widemult(p[b],q[a])) # 2 mult, add, neg

function test(alg,pts, y)
	for p in pts
		for q in pts
			@inbounds	y[1] += alg(p,q,1,2)
			@inbounds	y[2] += alg(p,q,1,3)
			@inbounds	y[3] += alg(p,q,2,3)
		end
	end
	y
end

function check(pts)
	for i in 1:length(pts)
		for j in i+1:length(pts)
			d1 = det2d(pts[i],pts[j], 1,2)
			d2 = det2w(pts[i],pts[j], 1,2)
			delta = sign(d1) - sign(d2)
			if delta ≠ zero(delta)
			 	println((i,j, delta, d1, d2, d1-d2) )
			end
		end
	end
end

function baddets(pts)
	ind = Int64[]
	for i in 1:length(pts)
		for j in i+1:length(pts)
			d1 = det2d(pts[i],pts[j], 1,2)
			d2 = det2w(pts[i],pts[j], 1,2)
			delta = sign(d1) - sign(d2)
			if delta ≠ zero(delta)
				# println((i,j, delta, d1, d2, d1-d2) )
				push!(ind, i)
				push!(ind, j)
			end
		end
	end
	unique!(ind)
end

function test(pts)
	println(size(pts))
	for alg in det2algs
		y = @MVector [alg(pts[1],pts[2],1,2), alg(pts[1],pts[2],1,3), alg(pts[1],pts[2],2,3)]
		test(alg, pts[1:2], y) # make sure it compiles first
		@time test(alg, pts, y)
	end
end

det2algs = [det2, det2f, det2w, det2wm, det2i,det2fi,det2wi,det2mi]
println("Float32")
pts = [rand(Float32,3).+1.0f0 for i in 1:100_000]
#test(pts)
pt = pts[baddets(pts)]
check(pt)

#using InteractiveUtils
#@code_warntype test(det2, pts[1:2], y)
println("Double32")
ptd = map(p->Double32.(p), pts[1:1000])
#test(ptd)
#check(ptd)

println("Int32")
pti = [rand(Int32,3) for i in 1:10_000]
#test(pti)
#check(pti)

det3(p,q,r, a,b,c) = muladd(p[a], det2(q,r,b,c), # 9 mult, 5 add, 3 neg,
					 muladd(p[b], det2(r,q,a,c),
						    p[c]*det2(q,r,a,b)))
det4ck(p,q,r,s, a,b,c,d) = muladd(p[a], det3(q,r,s, b,c,d), # 40 mult, 23 add, 12 neg
						 muladd(p[b], det3(q,r,s, a,d,c),
						 muladd(p[c], det3(q,r,s, a,b,d),
								p[d]*det3(q,r,s, a,c,b))))

# 2x2 minors is better
det4(p,q,r,s, a,b,c,d) = Double64(det2(p,q, a,b))*Double64(det2(r,s, c,d)) - # 30 mult, 17 add, 12 neg
						 Double64(det2(p,q, a,c))*Double64(det2(r,s, b,d)) +
						 Double64(det2(p,q, a,d))*Double64(det2(r,s, b,c)) +
						 Double64(det2(p,q, b,c))*Double64(det2(r,s, a,d)) -
						 Double64(det2(p,q, b,d))*Double64(det2(r,s, a,c)) +
						 Double64(det2(p,q, c,d))*Double64(det2(r,s, a,b))

function Plucker3L(p,q)
	diff = q.-p
	(q.-p, det2(p,d, 1,2), det2(p,d, 1,3), det2(p,d, 1,4), det2(p,d, 2,3), det2(p,d, 2,4), det2(p,d, 3,4))

det4dbl(p,q,r,s, a,b,c,d) = muladd(det2(p,q, a,b), det2(r,s, c,d), # 30 mult, 17 add, 12 neg
						 muladd(det2(p,q, c,a), det2(r,s, b,d),
						 muladd(det2(p,q, a,d), det2(r,s, b,c),
						 muladd(det2(p,q, b,c), det2(r,s, a,d),
						 muladd(det2(p,q, d,b), det2(r,s, a,c), det2(p,q, c,d))*det2(r,s, a,b)) ) ) )

end

max_stars_repo_name: snoeyink/CGinJulia

max_issues_repo_name: snoeyink/CGinJulia

max_forks_repo_name: snoeyink/CGinJulia


text:
# Load modules
using Plots, LaTeXStrings
pyplot()
using ThinFilmsTools

# Wavelength range [nm]
λ = 200:1100
# Angle of incidence [degrees]
θ = [0.]
# Polarisation (1.0 = p, 0.0 = s, between 0.0 and 1.0 = average)
pol = 0.5
beam = PlaneWave(λ, θ; p=pol)

# Refractive indices of incident (0) and substrate (2)
incident = RIdb.air(beam.λ)
emergent = RIdb.fused_silica_uv(beam.λ)

# Define the RI model to use
layers = [
    LayerTMMO(incident),
    ModelFit(:tauclorentz),
    LayerTMMO(emergent),
]

# Measured absolute transmittance
Texp = SpectraDB.hafnia_spectrum(beam.λ)
plot(beam.λ, Texp)
gui()

## Optimization using BBO

# One oscillator
seed = [vcat(175.0, # thickness
             [3.0, 5.5,
              545, 6.2, 13.2], # osc 1
        ),
]
lb = 0.01.*seed
ub = 2.0.*seed
lb[1][1] = 150
ub[1][1] = 200
lb[1][3] = 4.0
ub[1][3] = 6.0

solOptim = fit_tmm_optics(
    Transmittance(Texp), seed, beam, layers;
    alg=:BBO,
    SearchRange=Utils.unfoldbnd(lb,ub),
    Method=:de_rand_1_bin,
)

plot(FitSpectrum(),
    solOptim.beam.λ, solOptim.spectrumExp, solOptim.spectrumFit,
    xaxis=("Wavelength [nm]"), yaxis=("Transmittance"),
)
gui()


text:
using NLopt,ForwardDiff,DiffResults

minigerg = GERG2008([19,1]) #H2S + Ch4
aire = GERG2008([2,16])
xaire = [0.79,0.21]
P0 = 4.53e6
T0 = 190
x0 = [0.05,0.95]
x01 = [1.731e-2,0.98269]
x02 = [0.06618,0.93382]

function fmin_2phases(x::Vector,model::AbstractHelmholtzModel,x0::Vector,T0)
    nt = length(x0)
    return _helmholtzn(model,x[nt+1],T0,x[1:(nt)])+_helmholtzn(model,x[end],T0,x[(nt+2):(end-1)]) -
    _helmholtzn(model,x0[end],T0,x0[1:end-1])
end

function nlopt_form(f,xx,gg)
    if length(gg) > 0
        df = DiffResults.GradientResult(xx)
        df = ForwardDiff.gradient!(df,f,xx)
        gg .= DiffResults.gradient(df)
        return DiffResults.value(df)
    else
        return f(xx)
    end
end

#uses x = [x V], #x0 is feed composition (one phase)
function sum_restriction(result::Vector,x::Vector,grad::Matrix,x0::Vector)
    nt = length(x0)
    Π = Int32(length(x)/nt)

    xzero = zero(eltype(x))
    xone = one(eltype(x))
    if length(grad) > 0
        grad .= ones(Π*nt,nt)
        for i = 1:nt
            sum_n = xzero
            for j = 1:Π
                nn = nt*(j-1)+i
                sum_n += x[nn]
            end
            result[i]=sum_n-x0[i]
        end
    else
        for i = 1:nt
            sum_n = xzero
            for j = 1:Π
                sum_n += x[nt*(j-1)+i]
            end
            result[i]=sum_n-x0[i]
        end
    end
    return result
end

function sum_restriction2(result,x,grad,x0)

    nt = length(x0)
    Π = Int32(length(x)/nt)

    xzero = zero(eltype(x))
    for i = 1:nt
        sum_n = xzero
        for j = 1:Π
            sum_n += x[nt*(j-1)+i]
        end
        result[i]=sum_n-x0[i]
    end

    if length(grad)>0
        grad=ones(size(grad))
    end
    return result
end

function _solver2(model,v0,T0,x0)
    r = rand()

    x00 = vcat(vcat(r*x0,r*v0),vcat((1-r)*x0,(1-r)*v0))
    nn = length(x00)
    lower = zeros(nn)
    upper = ones(nn)
    tol = upper*1e-06

    opt = NLopt.Opt(:LN_COBYLA, nn)
    opt.xtol_rel = 1e-6
    opt.lower_bounds = lower
    opt.upper_bounds = upper
    equality_constraint!(opt, (r,x,g) -> sum_restriction2(r,x,g,x0),tol)
    fmin1 = (x,g) -> nlopt_form(vec -> fmin_2phases(vec,model,x0,T0),x,g)
    opt.min_objective = fmin1

    (min_f,min_x,status) = NLopt.optimize(opt, x00)

    return (min_f,min_x,status)
    return nn
end


text:
"""
Episode

    A container to keep track of environment feedback during an episode.
"""
mutable struct Episode
    total_reward
    niter
end

Episode() = Episode(0.0, 0)

"""
Environment

    Abstract stateful environment.

    Keeps track of the environment state and implements update eqn. s', r = ρ(s, a).
"""
abstract type AbstractEnvironment end
function reset!(env::AbstractEnvironment) end
function step!(env::AbstractEnvironment, action) end
function available_actions(env::AbstractEnvironment) end

"""
RandomWalk

    Random walk with random rewards and meaningless actions.
"""
mutable struct RandomWalk <: AbstractEnvironment
    state
    action_space
end

RandomWalk() = RandomWalk([0.0], [1, 2])

function reset!(env::RandomWalk)
    env.state = [0.0]
    env.action_space = [1, 2]
    return env.state, false  # state, done
end

function step!(env::RandomWalk, action)
    @assert(action in env.action_space)
    env.state[1] += randn()
    r = randn()
    done = abs(env.state[1]) > 5.0
    info = nothing
    if !done
        env.action_space = available_actions(env)
    end
    return env.state, r, done, info
end

function available_actions(env::RandomWalk)
    return [1, 2]
end


text:
using BinaryProvider

products = Product[
  LibraryProduct(
    joinpath(expanduser("~"), "usr", "local", "lib"), "libpoloapi", :libpolo
  ),
]

if any(!satisfied(p) for p in products)
  error("POLO C-API is not installed properly on your system.")
end

write_deps_file(joinpath(@__DIR__, "deps.jl"), products)


text:
using ProjectiveGeometricAlgebra3d
using ProjectiveGeometricAlgebra3d: sandwich, MultiVectorEven, MultiVector2, geomul, sandwich, MultiVector3, MultiVectorAll
const PGA = ProjectiveGeometricAlgebra3d
using Test

e = MultiVectorAll(e=1)
e0 = MultiVectorAll(e0=1)
e1 = MultiVectorAll(e1=1)
e2 = MultiVectorAll(e2=1)
e3 = MultiVectorAll(e3=1)

e01 = MultiVectorAll(e01=1)
e02 = MultiVectorAll(e02=1)
e03 = MultiVectorAll(e03=1)

e12 = MultiVectorAll(e12=1)
e13 = MultiVectorAll(e13=1)
e23 = MultiVectorAll(e23=1)

e012 = MultiVectorAll(e012=1)
e013 = MultiVectorAll(e013=1)
e023 = MultiVectorAll(e023=1)
e123 = MultiVectorAll(e123=1)

e0123=MultiVectorAll(e0123=1)
z=MultiVectorAll()

@test e0*e0 === z
@test e0*e1 === e01
@test e0*e2 === e02
@test e2*e0 === -e02
@test e0*e3 === e03
@test e0*e1*e2*e2 === e01

@test e0∧e1 === e01
@test e1∧e0 === -e01
@test e1∧e2∧e3 === e123

@test rev(e) === e
@test rev(e3) === e3
@test rev(e13) === -e13
@test rev(e123) === -e123
@test rev(e0123) === e0123

@test e*dual(e) === e0123
@test e0*dual(e0) === e0123
@test e1*dual(e1) === e0123
@test e2*dual(e2) === e0123

@test e ≈ e
@test !(e≈ e1)
@test e ≈ e1 atol=2

@testset "laws" begin
    for _ in 1:1000
        a = PGA.random(PGA.MultiVector{Float64})
        b = PGA.random(PGA.MultiVector{Float64})
        c = PGA.random(PGA.MultiVector{Float64})

        ### associativity
        @test (a∧b)∧c ≈ a∧(b∧c)
        @test (a*b)*c ≈ a*(b*c)
        @test (a∨b)∨c ≈ a∨(b∨c)

        ### involutions
        @test dual(dual(a)) === a
        @test dual(a∧b) ≈ dual(a) ∨ dual(b)
        @test rev(rev(a)) === a

        λ = randn()
        μ = randn()
        @test (λ*a + μ*b)*c ≈ (λ*a*c) + (μ*b*c)
        @test (λ*a + μ*b)∧c ≈ (λ*a∧c) + (μ*b∧c)
        @test (λ*a + μ*b)∨c ≈ (λ*a∨c) + (μ*b∨c)
    end
end

include("test_euclid.jl")


text:
module Aircraft_module

push!(LOAD_PATH, "./")

using CAS_module, Constants_module

export Dynamics, Observation, Sensor, Aircraft
export Nearest_Intruder_Tracker, Multi_Threat_Tracker

"""
	Define the Dynamics object
"""
mutable struct Dynamics
	x::Float64
	y::Float64
	dt::Float64
	v::Float64
	heading::Float64
	turn_rate::Float64
end

"""
	Define the Observation object
"""
mutable struct Observation
	rho::Float64 # Distance from ownship to intruder
	theta::Float64 # Angle to intruder relative to ownship heading direction
	psi::Float64 # Heading angle of intruder relative to ownship heading direction
	sos::Float64 # Speed of ownship
	soi::Float64 # Speed of intruder
	tau::Float64 # Time until loss of vertical separation
	pRA::Int64 # Previous advisory
	intruder_pos_angle::Float64
end

"""
	Define the Sensor object
"""
abstract type Sensor end

"""
	Define the tracking mode: nearest intruder
"""
mutable struct Nearest_Intruder_Tracker <: Sensor
	sensing_range::Float64
	obs::Array{Observation, 1}
	# default constructor
	function Nearest_Intruder_Tracker(; r::Float64=sensing_range)
		this = new()
		this.sensing_range = r
		this.obs = []
		return this
	end
end

"""
	Define the tracking mode: multi-threat
"""
mutable struct Multi_Threat_Tracker <: Sensor
	sensing_range::Float64
	obs::Array{Observation, 1}
	# default constructor
	function Multi_Threat_Tracker(; r::Float64=sensing_range)
		this = new()
		this.sensing_range = r
		this.obs = [] # empty array
		return this
	end
end

"""
	Define the Aircraft object
"""
mutable struct Aircraft
	id::Base.Random.UUID
	dyn::Dynamics

	x_init::Float64
	y_init::Float64

	x_dest::Float64 # within airspace
	y_dest::Float64 # within airspace

	colli::Bool
	desti::Bool

	advisory::Symbol # Current resolution advisory
	cas::CAS
	sensor::Sensor

	created_time::Int64
	nominal_route_len::Float64
	actual_route_len::Float64

	dist_to_dest::Float64

	# statistics for each individual agent (not global stats)
	t_to_nominal_t::Float64
	d_to_nominal_d::Float64
	cumulative_num_ac_encountered::Int64
	num_ac_encountered_arr::Array{Int64,1}
end

end # module

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
using JuMP, AmplNLWriter, Bonmin_jll

J = 11;
B = 10.0;
θ = 0.0024;
M = zeros(J).+0.1;
Gt = 1.0;
g = 9.81;
L = 1000.0;
n0 = 0.033;
nt = 0.0133;
Q = 20;
ηu = 0.0;
ηd = 0.0;
f = 1;

I = tan(θ);
Δx = L/(J-1);
x = LinRange(0,L,J);
E = (L.-x)*I;

z = zeros(J)

result = zeros(J)

function explicity_enumeration(J::Int,result::Array)
    for k in 2:J-1
        model = Model();
        set_optimizer(model,Ipopt.Optimizer)

        @variable(model, 0.0 <= H[1:J] <= L, start = 2.0);
        @variable(model, 0.0 <= A[1:J] <= L, start = 20.0);
        @variable(model, 0.0 <= SM[1:J] <= L, start = 0.0);
        @variable(model, 0.0 <= V[1:J] <= L, start = 1.0);
        @variable(model, 0.0 <= W[1:J] <= L, start = 0.0);

        @constraint(model, [j in 1:J-1], V[j]/(2*g) + H[j] + E[j] == V[j+1]/(2*g) + H[j+1] + E[j+1] + SM[j]*Δx + W[j]*z[j]);
        @constraint(model,[j in 1:J], A[j] == B*H[j]);
        @constraint(model, H[J] == 2);
        @NLexpression(model, R[j in 1:J], A[j]/(B+2*H[j]));
        @constraint(model, [j in 1:J], V[j] <= g*H[j]);
        @NLexpression(model, S[j in 1:J], V[j]*(n0^2+nt^2*z[j]+2*n0*nt*z[j])/(R[j]^(4/3)));
        @NLconstraint(model,[j in 1:J-1], SM[j] == (S[j] + S[j+1])/2);
        @constraint(model, [j in 1:J], H[j] >= Gt*z[j]);
        @constraint(model, [j in 1:J-1], -M[j]*H[j] <= H[j+1] - H[j]);
        @constraint(model, [j in 1:J-1], H[j+1] - H[j] <= M[j]*H[j]);

        @NLobjective(model, Max, -H[1]*ηu + H[J]*ηd + sum(W[j]*z[j] for j in 1:J) - sum((V[j]*A[j]^2 - Q^2)^2 for j in 1:J)/J)

        status = optimize!(model)
        result[k] = objective_value(model)
    end
    return result
end

max_stars_repo_name: Thiago-NovaesB/alpheus

max_issues_repo_name: Thiago-NovaesB/alpheus

max_forks_repo_name: Thiago-NovaesB/alpheus


text:
using DataFrames, CSV

@testset "Order-preserving hash for samples" begin
    increasing_array_1d = range(1,100,length=100)
    increasing_array_1d_smaller = range(-50,49,length=100)
    increasing_array_2d = reshape(range(1,100,length=100), (10,10))
    increasing_array_2d_smaller = reshape(range(-50,49,length=100), (10,10))

    @test orderinghash(10) < orderinghash(10000)
    @test orderinghash(-10) < orderinghash(0)
    @test orderinghash(increasing_array_1d_smaller) < orderinghash(increasing_array_1d)
    @test orderinghash(increasing_array_2d_smaller) < orderinghash(increasing_array_2d)
    @test orderinghash("abcd") < orderinghash("abce")
    @test orderinghash("abcd") < orderinghash("abcdf")
    @test orderinghash("Abcdf") < orderinghash("abcd")
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
include("../base.jl")
include("fitnessFunctions.jl")
import Random
println("hola")
println("hola")
println("hola")
"""
rng = Random.MersenneTwister(6)

setIndividualType(CGPGenotype)
setRNG(rng)
setGPExperimentInfo()
setGenerator(rampedHalfHalfGenerator, popSize = 10, generateOneByOne=false)
genPopulation!(GenJ)
gen1, gen2 = getGenotype(GenJ._population[1]), getGenotype(GenJ._population[6])
phen1 = getPhenotype(gen1)
phen2 = getPhenotype(gen2)
_,_, pointDepth = selectCrossPoint(gen1, rng)
selectCrossPoint(gen2, rng, minDepth=pointDepth, maxDepth=3)
#children = subtreeCross(gen1, gen2, GenJ._experimentInfo._GPExperimentInfo, rng)
#mutatedParent = subtreeMutation(parent1, GenJ._experimentInfo._GPExperimentInfo._functionSet, GenJ._experimentInfo._GPExperimentInfo._terminalSet, 3, GenJ._experimentInfo._rng, canonicalFullGenerator)
#fenotipo = getPhenotype(parent1, GenJ._experimentInfo._GPExperimentInfo)

"""

f(x, y) = (y/5.5)+x^2

println("hola")
println("hola")
x = collect(-5:0.5:5)
y = collect(0:0.5:10)
nvalues = length(x)
objs = Array{Number}(undef, nvalues)
for i=1:nvalues
    objs[i] = f(x[i], y[i])
end

setIndividualType(CGPGenotype)
setRandomSeed(2)
setGPExperimentInfo()
setAlgorithm(SPEA, 10)
setStopCondition(maxIterations=50)
setEvaluator([FitnessFunction(compareFunctions, objs, x, y, weight=-1), FitnessFunction(nNodes, weight=-0.2)],
            globalFitnessFunction=FitnessFunction(pareto), compareFunctionArgs="global")
setGenerator(rampedHalfHalfGenerator, popSize = 50, generateOneByOne = false)
setSelector(tournamentSelector, 2)
setCrossoverOperator(subtreeCross, probability=0.9)
setMutationOperator(subtreeMutation, probability=0.2)
setReplacementOperator(replaceAllPopulation, eliteSize=5)
setExperimentSummary(batchSize=20, printDuringExperiment=true)

@time runGenJ(verbose=false)
"""

f(x, y) = (y/5.5)+x^2

println("hola")
println("hola")
x = collect(-5:0.5:5)
y = collect(0:0.5:10)
nvalues = length(x)
objs = Array{Number}(undef, nvalues)
for i=1:nvalues
    objs[i] = f(x[i], y[i])
end

setIndividualType(CGPGenotype)
setRandomSeed(54)
setGPExperimentInfo()
setStopCondition(maxIterations=50)
setEvaluator([FitnessFunction(compareFunctions, objs, x, y, weight=-1)])
setGenerator(rampedHalfHalfGenerator, popSize = 50, generateOneByOne = false)
setSelector(tournamentSelector, 4)
setCrossoverOperator(subtreeCross, probability=0.9)
setMutationOperator(pointMutation, 0.2, probability=0.1)
setReplacementOperator(replaceAllPopulation, eliteSize=5)
setExperimentSummary(batchSize=5, printDuringExperiment=true)

@time runGenJ(verbose=false)
"""

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
using .FluxTraining

import .FluxTraining: Loggables, _combinename, cpu

struct WandbBackend <: FluxTraining.LoggerBackend
    logger::WandbLogger

    WandbBackend(; kwargs...) = new(WandbLogger(; kwargs...))
end

function Base.show(io::IO, backend::WandbBackend)
    print(io, "WandbBackend(")
    print(io, backend.logger)
    return print(io, ")")
end

for func in (:log, :close)
    @eval begin
        Base.$(func)(wa::WandbBackend, args...; kwargs...) = $(func)(wa.logger, args...; kwargs...)
    end
end

for func in (:update_config!, :finish, :save, :get_config)
    @eval begin
        $(func)(wa::WandbBackend, args...; kwargs...) = $(func)(wa.logger, args...; kwargs...)
    end
end

function FluxTraining.log_to(backend::WandbBackend, value::Loggables.Value, name, i; group=())
    name = _combinename(name, group)
    if length(group) == 0
        log(backend.logger, Dict(name => cpu(value.data)))
    else
        log(backend.logger, Dict(name => cpu(value.data), group[1] => i))
    end
end

function FluxTraining.log_to(backend::WandbBackend, image::Loggables.Image, name, i; group=())
    name = _combinename(name, group)
    imgs = Image.(collect(image.data))
    if length(group) == 0
        log(backend.logger, Dict([name * "_$i" => img for (i, img) in enumerate(imgs)]...))
    else
        log(backend.logger, Dict([name * "_$i" => img for (i, img) in enumerate(imgs)]..., group[1] => i))
    end
end

function FluxTraining.log_to(backend::WandbBackend, text::Loggables.Text, name, i; group=())
    name = _combinename(name, group)
    if length(group) == 0
        log(backend.logger, Dict(name => text.data))
    else
        log(backend.logger, Dict(name => text.data, group[i] => i))
    end
end

function FluxTraining.log_to(backend::WandbBackend, hist::Loggables.Histogram, name, i; group=())
    name = _combinename(name, group)
    if length(group) == 0
        log(backend.logger, Dict(name => Histogram(cpu(hist.data))))
    else
        log(backend.logger, Dict(name => Histogram(cpu(hist.data)), group[1] => i))
    end
end

export WandbBackend


text:
using BenchmarkTools, LazyArtifacts

t = @elapsed using Vlasiator
println("Julia version is $VERSION")
println(string("Vlasiator.jl loading time: \e[33;1;1m$t\e[m seconds"))
println()
println("Benchmarking Vlasiator.jl...")
println()

directory = artifact"testdata"
files = ("bulk.1d.vlsv", "bulk.2d.vlsv", "bulk.amr.vlsv")

const SUITE = BenchmarkGroup()

SUITE["read"] = BenchmarkGroup(["variable"])
file = joinpath(directory, files[2])
meta = load(file)
SUITE["read"]["meta"] = @benchmarkable load($file)
SUITE["read"]["DCCRG"] = @benchmarkable readvariable($meta, "proton/vg_rho")
ids = 100:110
SUITE["read"]["DCCRG_select_small"] =
   @benchmarkable readvariable($meta, "proton/vg_rho", $ids)
ids = 3000:6300
SUITE["read"]["DCCRG_select_large"] =
   @benchmarkable readvariable($meta, "proton/vg_rho", $ids)
SUITE["read"]["FG"] = @benchmarkable readvariable($meta, "fg_b")
file = joinpath(directory, files[1])
meta = load(file)
SUITE["read"]["VDF"] = @benchmarkable readvcells($meta, 5; species="proton")

SUITE["VTK"] = BenchmarkGroup(["conversion"])
file = joinpath(directory, files[3])
meta = load(file)
SUITE["VTK"]["AMR"] = @benchmarkable Vlasiator.fillmesh($meta,
   $["proton/vg_rho", "proton/vg_v", "fg_b", "vg_boundarytype"])


text:
export iauNut00b
"""
Nutation, IAU 2000B model.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  canonical model.

Given:
   date1,date2   double    TT as a 2-part Julian Date (Note 1)

Returned:
   dpsi,deps     double    nutation, luni-solar + planetary (Note 2)

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The nutation components in longitude and obliquity are in radians
   and with respect to the equinox and ecliptic of date.  The
   obliquity at J2000.0 is assumed to be the Lieske et al. (1977)
   value of 84381.448 arcsec.  (The errors that result from using
   this function with the IAU 2006 value of 84381.406 arcsec can be
   neglected.)

   The nutation model consists only of luni-solar terms, but
   includes also a fixed offset which compensates for certain long-
   period planetary terms (Note 7).

   3. This function is an implementation of the IAU 2000B abridged
   nutation model formally adopted by the IAU General Assembly in
   2000.  The function computes the MHB_2000_SHORT luni-solar
   nutation series (Luzum 2001), but without the associated
   corrections for the precession rate adjustments and the offset
   between the GCRS and J2000.0 mean poles.

   4. The full IAU 2000A (MHB2000) nutation model contains nearly 1400
   terms.  The IAU 2000B model (McCarthy & Luzum 2003) contains only
   77 terms, plus additional simplifications, yet still delivers
   results of 1 mas accuracy at present epochs.  This combination of
   accuracy and size makes the IAU 2000B abridged nutation model
   suitable for most practical applications.

   The function delivers a pole accurate to 1 mas from 1900 to 2100
   (usually better than 1 mas, very occasionally just outside
   1 mas).  The full IAU 2000A model, which is implemented in the
   function iauNut00a (q.v.), delivers considerably greater accuracy
   at current dates;  however, to realize this improved accuracy,
   corrections for the essentially unpredictable free-core-nutation
   (FCN) must also be included.

   5. The present function provides classical nutation.  The
   MHB_2000_SHORT algorithm, from which it is adapted, deals also
   with (i) the offsets between the GCRS and mean poles and (ii) the
   adjustments in longitude and obliquity due to the changed
   precession rates.  These additional functions, namely frame bias
   and precession adjustments, are supported by the SOFA functions
   iauBi00  and iauPr00.

   6. The MHB_2000_SHORT algorithm also provides "total" nutations,
   comprising the arithmetic sum of the frame bias, precession
   adjustments, and nutation (luni-solar + planetary).  These total
   nutations can be used in combination with an existing IAU 1976
   precession implementation, such as iauPmat76,  to deliver GCRS-
   to-true predictions of mas accuracy at current epochs.  However,
   for symmetry with the iauNut00a  function (q.v. for the reasons),
   the SOFA functions do not generate the "total nutations"
   directly.  Should they be required, they could of course easily
   be generated by calling iauBi00, iauPr00 and the present function
   and adding the results.

   7. The IAU 2000B model includes "planetary bias" terms that are
   fixed in size but compensate for long-period nutations.  The
   amplitudes quoted in McCarthy & Luzum (2003), namely
   Dpsi = -1.5835 mas and Depsilon = +1.6339 mas, are optimized for
   the "total nutations" method described in Note 6.  The Luzum
   (2001) values used in this SOFA implementation, namely -0.135 mas
   and +0.388 mas, are optimized for the "rigorous" method, where
   frame bias, precession and nutation are applied separately and in
   that order.  During the interval 1995-2050, the SOFA
   implementation delivers a maximum error of 1.001 mas (not
   including FCN).

References:

   Lieske, J.H., Lederle, T., Fricke, W., Morando, B., "Expressions
   for the precession quantities based upon the IAU /1976/ system of
   astronomical constants", Astron.Astrophys. 58, 1-2, 1-16. (1977)

   Luzum, B., private communication, 2001 (Fortran code
   MHB_2000_SHORT)

   McCarthy, D.D. & Luzum, B.J., "An abridged model of the
   precession-nutation of the celestial pole", Cel.Mech.Dyn.Astron.
   85, 37-49 (2003)

   Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
   Francou, G., Laskar, J., Astron.Astrophys. 282, 663-683 (1994)

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauNut00b(double date1, double date2, double *dpsi, double *deps)

function iauNut00b(date1::Real, date2::Real)
   # Allocate return values
   ref_dpsi = Ref{Float64}(0.0)
   ref_deps = Ref{Float64}(0.0)

   ccall((:iauNut00b, libsofa_c), Cvoid,
        (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
        convert(Float64, date1), convert(Float64, date2),
        ref_dpsi, ref_deps)

   return ref_dpsi[], ref_deps[]
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule DirectQhullHelper_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("DirectQhullHelper")
JLLWrappers.@generate_main_file("DirectQhullHelper", UUID("af82f447-cd1e-5518-b874-3ab508da728e"))
end  # module DirectQhullHelper_jll