text:
using FixedEffectModels, CSV, DataFrames, LinearAlgebra, Test
df = CSV.read(joinpath(dirname(pathof(FixedEffectModels)), "../dataset/Cigar.csv"))
df.id1 = df.State
df.id2 = df.Year
df.pid1 = categorical(df.id1)
df.pid2 = categorical(df.id2)

df.mid1 = div.(df.id1, Ref(10))
df.y = df.Sales
df.x1 = df.Price
df.z1 = df.Pimin
df.x2 = df.NDI
df.w = df.Pop

##############################################################################
##
## coefficients
##
##############################################################################

# simple
m = @formula y ~ x1
x = reg(df, m)
@test coef(x) ≈ [139.73446,-0.22974] atol = 1e-4
m = @formula y ~ x1
x = reg(df, m, weights = :w)
@test coef(x) ≈ [137.72495428982756,-0.23738] atol = 1e-4

df.SalesInt = round.(Int64, df.Sales)
m = @formula SalesInt ~ Price
x = reg(df, m)
@test coef(x) ≈ [139.72674,-0.2296683205] atol = 1e-4

# absorb
m = @formula y ~ x1 + fe(id1)
x = reg(df, m)
@test coef(x) ≈ [-0.20984] atol = 1e-4
@test x.iterations == 1

m = @formula y ~ x1 + fe(id1) + fe(id2)
x = reg(df, m)
@test coef(x)  ≈ [-1.08471] atol = 1e-4
m = @formula y ~ x1 + fe(id1) + fe(id1)*id2
x = reg(df, m)
@test coef(x) ≈  [-0.53470, 0.0] atol = 1e-4
m = @formula y ~ x1 + fe(id1)*id2
x = reg(df, m)
@test coef(x) ≈  [-0.53470, 0.0] atol = 1e-4
#@test isempty(coef(reg(df, @formula(y ~ 0), @fe(id1*x1))))
m = @formula y ~ x1 + fe(mid1)&fe(id2)
x = reg(df, m)
@test coef(x)  ≈  [-1.44255] atol = 1e-4

m = @formula y ~ x1 + fe(id1)&id2
x = reg(df, m)
@test coef(x) ≈ [13.993028174622104,-0.5804357763515606] atol = 1e-4
m = @formula y ~ x1 + id2&fe(id1)
x = reg(df, m)
@test coef(x) ≈ [13.993028174622104,-0.5804357763515606] atol = 1e-4
m = @formula y ~ 1 + id2&fe(id1)
x = reg(df, m)
@test coef(x) ≈ [174.4084407796102] atol = 1e-4

m = @formula y ~ x1 + fe(id1)&id2 + fe(id2)&id1
x = reg(df, m)
@test coef(x) ≈ [51.2359,- 0.5797] atol = 1e-4
m = @formula y ~ x1 + x2 + fe(id1)&id2 + fe(id2)&id1
x = reg(df, m)
@test coef(x) ≈ [-46.4464,-0.2546, -0.005563] atol = 1e-4
m = @formula y ~ 0 + x1 + x2 + fe(id1)&id2 + fe(id2)&id1
x = reg(df, m)
@test coef(x) ≈  [-0.21226562244177932,-0.004775616634862829] atol = 1e-4

# recheck these two below
m = @formula y ~ z1 + x1&x2&fe(id1)
x = reg(df, m)
@test coef(x)  ≈  [122.98713, 0.30933] atol = 1e-4
m = @formula y ~ z1 + (x1&x2)*fe(id1)
x = reg(df, m)
@test coef(x) ≈   [0.421406, 0.0] atol = 1e-4

# only one intercept
m = @formula y ~ 1 + fe(id1) + fe(id2)
x = reg(df, m)

# TO DO: REPORT INTERCEPT IN CASE OF FIXED EFFFECTS, LIKE STATA
df.id3  = categorical(mod.(1:size(df, 1), Ref(3)))
df.id4  = categorical(div.(1:size(df, 1), Ref(10)))

m = @formula y ~ id3
x = reg(df, m)
@test length(coef(x)) == 3

m = @formula y ~ 0 + id3
x = reg(df, m)
@test length(coef(x)) == 3

m = @formula y ~  id3 + fe(id4)
x = reg(df, m)

m = @formula y ~ id2 + fe(id1)
x = reg(df, m)

m = @formula y ~ id3&x1
x = reg(df, m)
@test length(coef(x)) == 4
m = @formula y ~ id3&x1 + x1
x = reg(df, m)
@test length(coef(x)) == 4

m = @formula y ~ id2 + fe(id1)
x = reg(df, m)

m = @formula y ~ id2&x1 + fe(id1)
x = reg(df, m)

# absorb + weights
m = @formula y ~ x1 + fe(id1)
x = reg(df, m, weights = :w)
@test coef(x) ≈ [- 0.21741] atol = 1e-4
m = @formula y ~ x1 + fe(id1) + fe(id2)
x = reg(df, m, weights = :w)
@test coef(x) ≈ [- 0.88794] atol = 1e-3
m = @formula y ~ x1 + fe(id1) + fe(id1)&id2
x = reg(df, m, weights = :w)
@test coef(x) ≈ [- 0.461085492] atol = 1e-4

# iv
m = @formula y ~ (x1 ~ z1)
x = reg(df, m)
@test coef(x) ≈ [138.19479,- 0.20733] atol = 1e-4
m = @formula y ~ x2 + (x1 ~ z1)
x = reg(df, m)
@test coef(x) ≈  [137.45096,0.00516,- 0.76276] atol = 1e-4
m = @formula y ~ x2 + (x1 ~ z1 + w)
x = reg(df, m)
@test coef(x) ≈ [137.57335,0.00534,- 0.78365] atol = 1e-4
## multiple endogeneous variables
m = @formula y ~ (x1 + x2 ~ z1 + w)
x = reg(df, m)
@test coef(x) ≈ [139.544, .8001, -.00937] atol = 1e-4
m = @formula y ~ 1 + (x1 + x2 ~ z1 + w)
x = reg(df, m)
@test coef(x) ≈ [139.544, .8001, -.00937] atol = 1e-4
result = [196.576, 0.00490989, -2.94019, -3.00686, -2.94903, -2.80183, -2.74789, -2.66682, -2.63855, -2.52394, -2.34751, -2.19241, -2.18707, -2.09244, -1.9691, -1.80463, -1.81865, -1.70428, -1.72925, -1.68501, -1.66007, -1.56102, -1.43582, -1.36812, -1.33677, -1.30426, -1.28094, -1.25175, -1.21438, -1.16668, -1.13033, -1.03782]

m = @formula y ~ x2 + (x1&pid2 ~ z1&pid2)
x = reg(df, m)
@test coef(x) ≈ result atol = 1e-4

# iv + weight
m = @formula y ~ (x1 ~ z1)
x = reg(df, m,  weights = :w)
@test coef(x) ≈ [137.03637,- 0.22802] atol = 1e-4

# iv + weight + absorb
m = @formula y ~ (x1 ~ z1) + fe(id1)
x = reg(df, m)
@test coef(x) ≈ [-0.20284] atol = 1e-4
m = @formula y ~ (x1 ~ z1) + fe(id1)
x = reg(df, m,  weights = :w)
@test coef(x) ≈ [-0.20995] atol = 1e-4
m = @formula y ~ x2 + (x1 ~ z1)
x = reg(df, m)
@test coef(x) ≈  [137.45096580480387,0.005169677634275297,-0.7627670265757879] atol = 1e-4
m = @formula y ~ x2 + (x1 ~ z1) + fe(id1)
x = reg(df, m)
@test coef(x)  ≈  [0.0011021722526916768,-0.3216374943695231] atol = 1e-4

# non high dimensional factors
m = @formula y ~ x1 + pid2
x = reg(df, m)
m = @formula y ~ pid2 + fe(id1)
x = reg(df, m)
m = @formula y ~ x1 + pid2 + fe(id1)
x = reg(df, m)
@test coef(x)[1] ≈ -1.08471 atol = 1e-4
m = @formula y ~ x1 + pid2 + fe(id1)
x = reg(df, m,  weights = :w)
@test coef(x)[1] ≈ -0.88794 atol = 1e-4
m = @formula y ~ x2 + (x1 ~ z1) + pid2 + fe(id1)
x = reg(df, m)
@test coef(x)[1] ≈ -0.00525 atol = 1e-4

##############################################################################
##
## Programming
##
##############################################################################
using StatsModels
reg(df, Term(:Sales) ~ Term(:NDI) + fe(Term(:State)) + fe(Term(:Year)), Vcov.cluster(:State))

# old syntax
m = @model y ~ x1 fe = pid1
x = reg(df, m)

m = @model y ~ x1 vcov = cluster(State + Year)
x = reg(df, m)

##############################################################################
##
## Functions
##
##############################################################################

# function
m = @formula y ~ log(x1)
x = reg(df, m)
@test coef(x)[1] ≈184.98520688 atol = 1e-4

# function defined in user space
mylog(x) = log(x)
m = @formula y ~ mylog(x1)
x = reg(df, m)

# Function returning Inf
df.x1_zero = copy(df.x1)
df.x1_zero[1] = 0.0
m = @formula y ~ log(x1_zero)
@test_throws  "Some observations for the regressor are infinite" x = reg(df, m)
##############################################################################
##
## collinearity
## add more tests
##
##############################################################################
# ols
df.x12 = df.x1
m = @formula y ~ x1 + x12
x = reg(df, m)
@test coef(x) ≈ [139.7344639806166,-0.22974688593485126,0.0] atol = 1e-4

# iv
df.x22 = df.x2
m = @formula y ~ x22 + x2 + (x1 ~ z1)
x = reg(df, m)
@test coef(x)[2] == 0 || coef(x)[3] == 0

df.zz1 = df.z1
m = @formula y ~ zz1 + (x1 ~ x2 + z1)
x = reg(df, m)
@test coef(x)[2] != 0.0

# catch when IV underidentified
@test_throws "Model not identified. There must be at least as many ivs as endogeneneous variables" reg(df, @formula(y ~ x1 + (x2 + w ~ x2)))

# Make sure all coefficients are estimated
using Random
Random.seed!(0)
df_r = DataFrame(x1 = randn(10000) * 100)
df_r.x2 = df_r.x1.^4
result = reg(df_r, @formula(x1 ~ x2 ))
@test sum(abs.(coef(result)) .> 0)  == 2

##############################################################################
##
## std errors
##
##############################################################################

# Simple - matches stata
m = @formula y ~ x1
x = reg(df, m)
@test stderror(x) ≈ [1.521269, 0.0188963] atol = 1e-6
# Stata ivreg - ivreg2 discrepancy - matches with df_add=-2
m = @formula y ~ (x1 ~ z1)
x = reg(df, m)
@test stderror(x) ≈ [1.53661, 0.01915] atol = 1e-4
# Stata areg
m = @formula y ~ x1 + fe(id1)
x = reg(df, m)
@test stderror(x) ≈  [0.0098003] atol = 1e-7

# White
# Stata reg - matches stata
m = @formula y ~ x1
x = reg(df, m, Vcov.robust())
@test stderror(x) ≈ [1.686791, 0.0167042] atol = 1e-6
# Stata ivreg - ivreg2 discrepancy - matches with df_add=-2
m = @formula y ~ (x1 ~ z1)
x = reg(df, m, Vcov.robust())
@test stderror(x) ≈ [1.63305, 0.01674] atol = 1e-4
# Stata areg - matches areg
m = @formula y ~ x1 + fe(id1)
x = reg(df, m, Vcov.robust())
@test stderror(x) ≈ [0.0110005] atol = 1e-7

# Clustering models
# cluster - matches stata
m = @formula y ~ x1
x = reg(df, m, Vcov.cluster(:id1))
@test stderror(x)[2] ≈ 0.0379228 atol = 1e-7
# cluster with fe - matches areg & reghdfe
m = @formula y ~ x1 + fe(id1)
x = reg(df, m, Vcov.cluster(:id2))
@test stderror(x) ≈ [0.0220563] atol = 1e-7
# stata reghxe - matches reghdfe (not areg)
m = @formula y ~ x1 + fe(id1)
x = reg(df, m, Vcov.cluster(:id1))
@test stderror(x) ≈ [0.0357498] atol = 1e-7
# iv + fe + cluster - matches ivreghdfe
m = @formula y ~ x2 + (x1 ~z1) + fe(id1)
x = reg(df, m, Vcov.cluster(:id1))
@test stderror(x) ≈ [0.0019704, 0.1893396] atol = 1e-7
# iv + fe + cluster + weights - matches ivreghdfe
m = @formula y ~ x2 + (x1 ~ z1) + fe(id1)
x = reg(df, m,  Vcov.cluster(:id1), weights = :w)
@test stderror(x) ≈ [0.000759, 0.070836] atol = 1e-6
# iv + fe + cluster + weights - matches ivreghdfe
m = @formula y ~ (x1 ~ z1) + fe(id1)
x = reg(df, m,  Vcov.cluster(:id1), weights = :w)
@test stderror(x) ≈ [0.0337439] atol = 1e-7
# multiway clustering - matches reghdfe
m = @formula y ~ x1
x = reg(df, m, Vcov.cluster(:id1, :id2))
@test stderror(x) ≈ [6.196362, 0.0403469] atol = 1e-6
# multiway clustering - matches reghdfe
m = @formula y ~ x1
x = reg(df, m, Vcov.cluster(:id1, :id2))
@test stderror(x) ≈ [6.196362, 0.0403469] atol = 1e-6
# fe + multiway clustering - matches reghdfe
m = @formula y ~ x1 + fe(id1)
x = reg(df, m, Vcov.cluster(:id1, :id2))
@test stderror(x) ≈ [0.0405335] atol = 1e-7
m = @formula y ~ x1 + fe(id1)
x = reg(df, m, Vcov.cluster(:id1, :id2))
@test stderror(x) ≈ [0.0405335] atol = 1e-7
# fe + clustering on interactions - matches reghdfe
#m = @formula y ~ x1 + fe(id1) vcov = cluster(id1&id2)
#x = reg(df, m)
#@test stderror(x) ≈ [0.0110005] atol = 1e-7
# fe partially nested in interaction clusters - matches reghdfe
#m=@formula y ~ x1 + fe(id1) vcov=cluster(id1&id2)
#x = reg(df, m)
#@test stderror(x) ≈ [0.0110005] atol = 1e-7
# regressor partially nested in interaction clusters - matches reghdfe
#m=@formula y ~ x1 + pid1 vcov=cluster(id1&id2)
#x = reg(df, m)
#@test stderror(x)[1:2] ≈ [3.032187, 0.0110005] atol=1e-5

#check palue printed is correct
m = @formula y ~ x1 + fe(id1)
x = reg(df, m, Vcov.cluster(:id1))
@test coeftable(x).mat[1, 4] ≈ 4.872723900371927e-7 atol = 1e-7

##############################################################################
##
## subset
##
##############################################################################
m = @formula y ~ x1 + pid1
x0 = reg(df[df.State .<= 30, :], m)

m = @formula y ~ x1 + pid1
x1 = reg(df, m, subset = df.State .<= 30)
@test length(x1.esample) == size(df, 1)
@test coef(x0) ≈ coef(x1) atol = 1e-4
@test vcov(x0) ≈ vcov(x1) atol = 1e-4

df.id1_missing = ifelse.(df.id1 .<= 30, df.id1, missing)
df.pid1_missing = categorical(df.id1_missing)
m = @formula y ~ x1 + pid1_missing
x2 = reg(df, m)
@test length(x2.esample) == size(df, 1)
@test coef(x0) ≈ coef(x2) atol = 1e-4
@test vcov(x0) ≈ vcov(x2) atol = 1e-2

m = @formula y ~ x1 + fe(id1)
x3 = reg(df, m, subset = df.State .<= 30)
@test length(x3.esample) == size(df, 1)
@test coef(x0)[2] ≈ coef(x3)[1] atol = 1e-4

m = @formula y ~ x1 + fe(id1_missing)
x4 = reg(df, m)
@test coef(x0)[2] ≈ coef(x4)[1] atol = 1e-4

# categorical variable as fixed effects
m = @formula y ~ x1 + fe(id1)
x5 = reg(df, m, subset = df.State .>= 30)

#Error reported by Erik
m = @formula y ~ z1 + CPI
x = reg(df, m, Vcov.cluster(:id1), subset = df.State .>= 30)
@test diag(x.vcov) ≈ [130.7464887, 0.0257875, 0.0383939] atol = 1e-4

##############################################################################
##
## R2
##
##############################################################################
m = @formula y ~ x1
x = reg(df, m)
@test r2(x) ≈ 0.0969 atol = 1e-4
@test adjr2(x) ≈ 0.09622618 atol = 1e-4

##############################################################################
##
## F Stat
##
##############################################################################
m = @formula y ~ x1
x = reg(df, m)
@test x.F  ≈ 147.82425 atol = 1e-4
m = @formula y ~ x1 + fe(id1)
x = reg(df, m)
@test x.F  ≈ 458.45825 atol = 1e-4
m = @formula y ~ (x1 ~ z1)
x = reg(df, m)
@test x.F  ≈ 117.17329 atol = 1e-4
m = @formula y ~ x1 + w
x = reg(df, m)
@test x.F ≈ 79.091576 atol = 1e-4
m = @formula y ~ x1
x = reg(df, m, Vcov.cluster(:id1))
@test x.F    ≈ 36.70275 atol = 1e-4
m = @formula y ~ (x1 ~ z1)
x = reg(df, m, Vcov.cluster(:id1))
@test x.F  ≈ 39.67227 atol = 1e-4
#  xtivreg2
m = @formula y ~ (x1 ~ z1) + fe(id1)
x = reg(df, m)
@test x.F  ≈ 422.46444 atol = 1e-4

# p value
m = @formula w ~ z1
x = reg(df, m)
@test x.p ≈ 0.003998718283554043 atol = 1e-4
m = @formula w ~ z1 + x1
x = reg(df, m)
@test x.p ≈ 3.369423076033613e-14 atol = 1e-4
##############################################################################
##
## F_kp r_kp statistics for IV. Difference degrees of freedom.
##
##############################################################################
m = @formula y ~ (x1 ~ z1)
x = reg(df, m)
@test x.F_kp ≈ 52248.79247 atol = 1e-4
m = @formula y ~ x2 + (x1 ~ z1)
x = reg(df, m)
@test x.F_kp ≈ 5159.812208193612 atol = 1e-4
m = @formula y ~ (x1 ~ z1 + CPI)
x = reg(df, m)
@test x.F_kp   ≈ 27011.44080 atol = 1e-4
m = @formula y ~ (x1 ~ z1) + fe(id1)
x = reg(df, m)
@test x.F_kp ≈ 100927.75710 atol = 1e-4

# exactly same with ivreg2
m = @formula y ~ (x1 ~ z1)
x = reg(df, m, Vcov.robust())
@test x.F_kp ≈ 23160.06350 atol = 1e-4
m = @formula y ~ (x1 ~ z1) + fe(id1)
x = reg(df, m, Vcov.robust())
@test x.F_kp  ≈ 37662.82808 atol = 1e-4
m = @formula y ~ CPI + (x1 ~ z1)
x = reg(df, m, Vcov.robust())
@test x.F_kp ≈ 2093.46609 atol = 1e-4
m = @formula y ~ (x1 ~ z1 + CPI)
x = reg(df, m, Vcov.robust())
@test x.F_kp  ≈ 16418.21196 atol = 1e-4

# like in ivreg2 but += 5 difference. combination iv difference, degrees of freedom difference?
# iv + cluster - F_kp varies from ivreg2 about 7e-4 (SEs off by more than 2 df)
m = @formula y ~ (x1 ~ z1)
x = reg(df, m, Vcov.cluster(:id1))
@test x.F_kp     ≈ 7249.88606 atol = 1e-4
# iv + cluster - F_kp varies from ivreg2 about 5e-6 (SEs off by more than 2 df)
m = @formula y ~ CPI + (x1 ~ z1)
x = reg(df, m, Vcov.cluster(:id1))
@test x.F_kp   ≈ 538.40393 atol = 1e-4
# Modified test values below after multiway clustering update
# iv + 2way clustering - F_kp matches ivreg2 (SEs match with df_add=-2)
m = @formula y ~ CPI + (x1 ~ z1)
x = reg(df, m, Vcov.cluster(:id1, :id2))
@test x.F_kp   ≈ 421.9651 atol = 1e-4
# multivariate iv + clustering - F_kp varies from ivreg2 about 3 (SEs off by more than 2 df)
m = @formula y ~ (x1 ~ z1 + CPI)
x = reg(df, m, Vcov.cluster(:id1))
@test x.F_kp      ≈ 4080.66081 atol = 1e-4
# multivariate iv + multiway clustering - F_kp varies from ivreg2 about 2 (SEs off by more than 2 df)
m = @formula y ~ (x1 ~ z1 + CPI)
x = reg(df, m, Vcov.cluster(:id1, :id2))
@test x.F_kp  ≈ 2873.1405 atol = 1e-4

##############################################################################
##
## Test singleton
##
##
##############################################################################
df.n = max.(1:size(df, 1), 60)
df.pn = categorical(df.n)
m = @formula y ~ x1 + fe(pn)
x = reg(df, m, Vcov.cluster(:id1))
@test x.nobs == 60

m = @formula y ~ x1 + fe(pn)
x = reg(df, m, Vcov.cluster(:id1), drop_singletons = false)
@test x.nobs == 1380

##############################################################################
##
## Test unbalanced panel
##
## corresponds to abdata in Stata, for instance reghxe wage emp [w=indoutpt], a(id year)
##
##############################################################################
df = CSV.read(joinpath(dirname(pathof(FixedEffectModels)), "../dataset/EmplUK.csv"))
df.id1 = df.Firm
df.id2 = df.Year
df.y = df.Wage
df.x1 = df.Emp
df.w = df.Output

methods_vec = [:lsmr, :lsmr_threads, :lsmr_cores]
if isdefined(FixedEffectModels.FixedEffects, :FixedEffectSolverLSMRGPU)
	push!(methods_vec, :lsmr_gpu)
end

for method in methods_vec
	@show method
	m = @formula y ~ x1 + fe(id1)
	x = reg(df, m, method = method)
	@test coef(x) ≈ [- 0.11981270017206136] atol = 1e-4
	m = @formula y ~ x1 + fe(id1)&id2
	x = reg(df, m, method = method)
	@test coef(x)  ≈ [-315.0000747500431,- 0.07633636891202833] atol = 1e-4
	m = @formula y ~ x1 + id2&fe(id1)
	x = reg(df, m, method = method)
	@test coef(x) ≈  [-315.0000747500431,- 0.07633636891202833] atol = 1e-4
	m = @formula y ~ 1 + id2&fe(id1)
	x = reg(df, m, method = method)
	@test coef(x) ≈  [- 356.40430526316396] atol = 1e-4
	m = @formula y ~ x1 + fe(id1)
	x = reg(df, m, weights = :w,  method = method)
	@test coef(x) ≈ [- 0.11514363590574725] atol = 1e-4

	# absorb + weights
	m = @formula y ~ x1 + fe(id1) + fe(id2)
	x = reg(df, m, method = method)
	@test coef(x)  ≈  [- 0.04683333721137311] atol = 1e-4
	m = @formula y ~ x1 + fe(id1) + fe(id2)
	x = reg(df, m, weights = :w, method = method)
	@test coef(x) ≈  [- 0.043475472188120416] atol = 1e-3

	## the last two ones test an ill conditioned model matrix
	m = @formula y ~ x1 + fe(id1) + fe(id1)&id2
	x = reg(df, m, method = method)
	@test coef(x)  ≈   [- 0.122354] atol = 1e-4
	@test x.iterations <= 30

	m = @formula y ~ x1 + fe(id1) + fe(id1)&id2
	x = reg(df, m, weights = :w, method = method)
	@test coef(x) ≈ [- 0.11752306001586807] atol = 1e-4
	@test x.iterations <= 50

	# same thing with float32 precision
	m = @formula y ~ x1 + fe(id1)
	x = reg(df, m, method = method, double_precision = false)
	@test coef(x) ≈ [- 0.11981270017206136] rtol = 1e-4
	m = @formula y ~ x1 + fe(id1)&id2
	x = reg(df, m, method = method, double_precision = false)
	@test coef(x)  ≈ [-315.0000747500431,- 0.07633636891202833] rtol = 1e-4
	m = @formula y ~ x1 + id2&fe(id1)
	x = reg(df, m, method = method, double_precision = false)
	@test coef(x) ≈  [-315.0000747500431,- 0.07633636891202833] rtol = 1e-4
	m = @formula y ~ 1 + id2&fe(id1)
	x = reg(df, m, method = method, double_precision = false)
	@test coef(x) ≈  [- 356.40430526316396] rtol = 1e-4
	m = @formula y ~ x1 + fe(id1)
	x = reg(df, m, weights = :w, method = method, double_precision = false)
	@test coef(x) ≈ [- 0.11514363590574725] rtol = 1e-4
	m = @formula y ~ x1 + fe(id1) + fe(id2)
	x = reg(df, m, method = method, double_precision = false)
	@test coef(x)  ≈  [- 0.04683333721137311] rtol = 1e-4
	m = @formula y ~ x1 + fe(id1) + fe(id2)
	x = reg(df, m, weights = :w, method = method, double_precision = false)
	@test coef(x) ≈  [- 0.043475472188120416] atol = 1e-3

end

# add tests with missing fixed effects
df.id1_missing = ifelse.(df.id1 .<= 30, missing, df.id1)

## test with missing fixed effects
m = @formula y ~ x1 + fe(id1_missing)
x = reg(df, m)
@test coef(x) ≈ [-.1093657] atol = 1e-4
@test stderror(x) ≈  [.032949 ] atol = 1e-4
@test r2(x) ≈ 0.8703 atol = 1e-2
@test adjr2(x) ≈ 0.8502 atol = 1e-2
@test x.nobs == 821

## test with missing interaction
df.id3 = df.id2 .>= 1980
m = @formula y ~ x1 + fe(id1_missing) & fe(id3)
x = reg(df, m)
@test coef(x) ≈ [-0.100863] atol = 1e-4
@test stderror(x) ≈  [0.04149] atol = 1e-4
@test x.nobs == 821

##############################################################################
##
## Missing dependent / independent variables
##
##############################################################################

# Missing
df1 = DataFrame(a=[1.0, 2.0, 3.0, 4.0], b=[5.0, 7.0, 11.0, 13.0])
df2 = DataFrame(a=[1.0, missing, 3.0, 4.0], b=[5.0, 7.0, 11.0, 13.0])
x = reg(df1,@formula(a ~ b))
@test coef(x) ≈ [ -0.6500000000000004, 0.35000000000000003] atol = 1e-4
x = reg(df1,@formula(b ~ a))
@test coef(x) ≈ [2.0, 2.8] atol = 1e-4
x = reg(df2,@formula(a ~ b))
@test coef(x) ≈ [ -0.8653846153846163, 0.3653846153846155] atol = 1e-4
x = reg(df2,@formula(b ~ a))
@test coef(x) ≈ [ 2.4285714285714253, 2.7142857142857157] atol = 1e-4

# Works with Integers
df1 = DataFrame(a=[1, 2, 3, 4], b=[5, 7, 11, 13], c = categorical([1, 1, 2, 2]))
x = reg(df1,@formula(a ~ b))
@test coef(x) ≈ [-0.65, 0.35] atol = 1e-4
x = reg(df1,@formula(a ~ b + fe(c)))
@test coef(x) ≈ [0.5] atol = 1e-4


text:
### A Pluto.jl notebook ###
# v0.17.0

using Markdown
using InteractiveUtils

# ╔═╡ 2d6d0588-f4e4-4cad-8353-9c0ba020fec6
begin
	using FileIO
	using ImageShow
	using ImageIO
	using PlutoUI
	using Cbc
	using DataFrames
	using Distances
	using JuMP
	using LinearAlgebra
end

# ╔═╡ d44f3f70-705b-43a2-85c0-a441e08fbbe6
load("Pydata_Jeddah_22_new-06.png")

# ╔═╡ a12d1135-ce89-414c-866e-14b4415e1990
md"""
# Geographical Clustering of Saudi Cities With Additional Constraint
**Originally Contributed by**: Matthew Helm ([with help from Mathieu Tanneau on Julia Discourse](https://discourse.julialang.org/t/which-jump-jl-solver-for-this-problem/43350/17?u=mthelm85))
**Customized by**: Hassan Alsawadi
"""

# ╔═╡ d35c109a-e5b6-4c01-a9c2-333b883783f2
md"""
#### Intializing packages

_When running this notebook for the first time, this could take up to 5 minutes. Hang in there!_
"""

# ╔═╡ 2877f488-1ab3-46da-8a49-7db65c3d8a45
md"""
## Objective
The goal of this exercise is to cluster 𝑛 cities into 𝑘 groups, minimizing the total pairwise distance between cities and ensuring that the variance in the total populations of each group is relatively small.

For this example, we'll use the 20 most populous cities in Saudi Arabia.
"""

# ╔═╡ cf6d3f77-781d-4cd0-ad45-9c5c58380354
cities = DataFrame(
    city=["Riyadh", "Jeddah", "Mecca", "Medina", "Dammam", "Abha", "Ha'il", "Hofuf", "Al-Mubarraz", "Sakakah", "Taif", "Khamis Mushait", "Jubail", "Buraydah", "Tabuk", "Najran", "Qatif", "Khobar", "Al Bahah", "Hafr Al-Batin"],
    population=[7676654, 4697000, 2042000, 1488782, 1252523, 1093705, 936465, 858395, 837000, 779480, 695000, 667000, 684531, 669000, 667000, 505652, 474573, 455541, 366000, 359000],
    lat=[24.633333, 21.543333, 21.4225, 24.466667, 26.433333, 18.216944, 27.516667, 25.383333, 25.383333, 29.969722, 21.275094, 18.3, 27, 26.333333, 28.397222, 17.491667, 26.556, 26.283333, 20.0125, 28.434167],
    lon=[46.716667, 39.172778, 39.823333, 39.6, 50.1, 42.505278, 41.683333, 49.583333, 49.6, 40.2, 40.406156, 42.733333, 49.666667, 43.966667, 36.578889, 44.132222, 49.996, 50.2, 41.465278, 45.963611])

# ╔═╡ 11c8323a-9cd5-4580-93d9-ebcd085b4551
md"""
## Model Specifics
We will cluster these 20 cities into 3 different groups and we will assume that the ideal or target population  𝑃  for a group is simply the total population of the 20 cities divided by 3:
"""

# ╔═╡ 0a9212cd-4c79-4752-930b-45dce6a1f3e8
begin
	n = size(cities,1)
	k = 3
	P = sum(cities.population) / k
end

# ╔═╡ 7f7439d3-3d3b-4098-8be5-ca9aa59373ce
md"""### Obtaining the distances between each city
Let's leverage the *Distances.jl* package to compute the pairwise Haversine distance between each of the cities in our data
set and store the result in a variable we'll call `dm`:
"""

# ╔═╡ 3e708ffa-05d9-4cf7-ba4c-11ca5a3e43e7
dm = Distances.pairwise(Haversine(6372.8), Matrix(cities[:, [3,4]])', dims=2)

# ╔═╡ 54082ac4-c31d-4d27-9cc3-5248018c6838
md"""Our distance matrix is symmetric so we'll convert it to a `LowerTriangular` matrix so that we can better interpret the
objective value of our model (if we don't do this the total distance will be doubled):"""

# ╔═╡ 6bc42cc7-e3d4-4470-bece-6af950280372
dm2 = LowerTriangular(dm)

# ╔═╡ a5f3abe7-a843-4ea5-8e77-dfd075b243a5
md"""
## Build the model
Now that we have the basics taken  care of, we can set up our model, create decision variables, add constraints, and then
solve.

First, we'll set up a model that leverages the [Cbc](https://github.com/coin-or/Cbc) solver. Next, we'll set up a binary
variable $x_{i,k}$ that takes the value $1$ if city $i$ is in group $k$ and $0$ otherwise. Each city must be in a group, so
we'll add the constraint $\sum_kx_{i,k} = 1$ for every $i$.
"""

# ╔═╡ ee41a6b3-fb13-4269-ab22-ad2a85bf26e5
begin
	model = Model(Cbc.Optimizer)

	@variable(model, x[1:n, 1:k], Bin)

	for i in 1:n
	    @constraint(model, sum(x[i,:]) == 1)
	end
end

# ╔═╡ dd262e8b-4697-4aa7-9bc7-a32ac25c96dc
md"""
The total population of a group $k$ is $Q_k = \sum_ix_{i,k}q_i$ where $q_i$ is simply the $i$th value from the `population`
column in our `cities` DataFrame. Let's add constraints so that $\alpha \leq (Q_k - P) \leq \beta$. We'll set $\alpha$
equal to -2,500,000 and $\beta$ equal to 2,500,000. By adjusting these thresholds you'll find that there is a tradeoff
between having relatively even populations between groups and having geographically close cities within each group. In
other words, the larger the absolute values of $\alpha$ and $\beta$, the closer together the cities in a group will be but
the variance between the group populations will be higher.
"""

# ╔═╡ d5fc7ca7-33cd-4f8a-aecf-8897bcd40e15
begin
	α = -2_500_000
	β = 2_500_000
	for i in 1:k
	    @constraint(model, (x' * cities.population)[i] - P <= β)
	    @constraint(model, (x' * cities.population)[i] - P >= α)
	end
end

# ╔═╡ 4ab84877-f5d8-4c1b-b0d2-346f0048bb20
md"""
Now we need to add one last binary variable $z_{i,j}$ to our model that we'll use to compute the total distance between the
cities in our groups, defined as  $\sum_{i,j}d_{i,j}z_{i,j}$. Variable $z_{i,j}$ will equal $1$ if cities $i$ and $j$ are
in the same group, and $0$ if they are not in the same group.

To ensure that $z_{i,j} = 1$ if and only if cities $i$ and $j$ are in the same group, we add the constraints $z_{i,j} \geq
x_{i,k} + x_{j,k} - 1$ for every pair $i,j$ and every $k$:
"""

# ╔═╡ 1054b6a4-1c28-46b6-88a4-6f7cae84edee
begin
	@variable(model, z[1:n,1:n], Bin)

	for k in 1:k, i in 1:n, j in 1:n
	    @constraint(model, z[i,j] >= x[i,k] + x[j,k] - 1)
	end
end

# ╔═╡ 07f58656-4826-4ce5-8c97-bf5570f0fc60
md"""
We can now add an objective to our model which will simply be to minimize the dot product of $z$ and our distance matrix,
`dm`. We can then call `optimize!` and review the results.
"""

# ╔═╡ 17088a0b-a568-431c-8d01-13f957e2a4c0
begin
	@objective(model, Min, dot(z,dm));

	optimize!(model)
end

# ╔═╡ 4ded521e-c41f-4927-854b-86f8912269db
md"""## Reviewing the Results
Now that we have results, we can add a column to our `cities` DataFrame for the group and then loop through our $x$
variable to assign each city to its group. Once we have that, we can look at the total population for each group and also
plot the cities and their groups to verify visually that they are grouped by geographic proximity.
"""

# ╔═╡ dab4e559-9dc6-4a72-8d38-9a0487b6251a
PlutoUI.with_terminal() do
	cities.group = zeros(n)

	for i in 1:n, j in 1:k
	    if round(value.(x)[i,j]) == 1.0
	        cities.group[i] = j
	    end
	end

	for group in groupby(cities, :group)
	    @show group
	    println("")
	    @show sum(group.population)
	    println("")
	end
end

# ╔═╡ aa340dc7-8833-4e40-b63e-c602096e40e3
md"""
The populations of each group are fairly even and we can see from the plot below that the groupings look good in terms of
geographic proximity:
"""

# ╔═╡ a2fad4bc-66fa-4466-ae29-b495590f035f
load("output.png")

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19"
ImageShow = "4e3cecfd-b093-5904-9786-8bbb286a6a31"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
Cbc = "~0.9.1"
DataFrames = "~1.3.2"
Distances = "~0.10.7"
FileIO = "~1.13.0"
ImageIO = "~0.6.1"
ImageShow = "~0.3.3"
JuMP = "~0.22.3"
PlutoUI = "~0.7.34"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ASL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "6252039f98492252f9e47c312c8ffda0e3b9e78d"
uuid = "ae81ac8f-d209-56e5-92de-9978fef736f9"
version = "0.1.3+0"

[[AbstractFFTs]]
deps = ["ChainRulesCore", "LinearAlgebra"]
git-tree-sha1 = "6f1d9bc1c08f9f4a8fa92e3ea3cb50153a1b40d4"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.1.0"

[[AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.3"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[BenchmarkTools]]
deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
git-tree-sha1 = "4c10eee4af024676200bc7752e536f858c6b8f93"
uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
version = "1.3.1"

[[BinaryProvider]]
deps = ["Libdl", "Logging", "SHA"]
git-tree-sha1 = "ecdec412a9abc8db54c0efc5548c64dfce072058"
uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
version = "0.5.10"

[[Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.8+0"

[[CEnum]]
git-tree-sha1 = "215a9aa4a1f23fbd05b92769fdd62559488d70e9"
uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
version = "0.4.1"

[[Calculus]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad"
uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
version = "0.5.1"

[[Cbc]]
deps = ["BinaryProvider", "CEnum", "Cbc_jll", "Libdl", "MathOptInterface", "SparseArrays"]
git-tree-sha1 = "6656166f484075dd146c9f452b1428116eaf76d4"
uuid = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
version = "0.9.1"

[[Cbc_jll]]
deps = ["ASL_jll", "Artifacts", "Cgl_jll", "Clp_jll", "CoinUtils_jll", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "OpenBLAS32_jll", "Osi_jll", "Pkg"]
git-tree-sha1 = "a3c5986d7713bce4260d9826deead060a17c8e2d"
uuid = "38041ee0-ae04-5750-a4d2-bb4d0d83d27d"
version = "200.1000.501+0"

[[Cgl_jll]]
deps = ["Artifacts", "Clp_jll", "CoinUtils_jll", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Osi_jll", "Pkg"]
git-tree-sha1 = "11eb7b7688925e9751b5d7a187aaa4291eae2664"
uuid = "3830e938-1dd0-5f3e-8b8e-b3ee43226782"
version = "0.6000.300+0"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "7dd38532a1115a215de51775f9891f0f3e1bac6a"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.12.1"

[[ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[Clp_jll]]
deps = ["Artifacts", "CoinUtils_jll", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "METIS_jll", "MUMPS_seq_jll", "OpenBLAS32_jll", "Osi_jll", "Pkg"]
git-tree-sha1 = "b1031dcfbb44553194c9e650feb5ab65e372504f"
uuid = "06985876-5285-5a41-9fcb-8948a742cc53"
version = "100.1700.601+0"

[[CodecBzip2]]
deps = ["Bzip2_jll", "Libdl", "TranscodingStreams"]
git-tree-sha1 = "2e62a725210ce3c3c2e1a3080190e7ca491f18d7"
uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd"
version = "0.7.2"

[[CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[CoinUtils_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "OpenBLAS32_jll", "Pkg"]
git-tree-sha1 = "44173e61256f32918c6c132fc41f772bab1fb6d1"
uuid = "be027038-0da8-5614-b30d-e42594cb92df"
version = "200.1100.400+0"

[[ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[ColorVectorSpace]]
deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "TensorCore"]
git-tree-sha1 = "3f1f500312161f1ae067abe07d13b40f78f32e07"
uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4"
version = "0.9.8"

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[CommonSubexpressions]]
deps = ["MacroTools", "Test"]
git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
version = "0.3.0"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "44c37b4636bc54afac5c574d2d02b625349d6582"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.41.0"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

[[DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "ae02104e835f219b8930c7664b8012c93475c340"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.3.2"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[DiffResults]]
deps = ["StaticArrays"]
git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805"
uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
version = "1.0.3"

[[DiffRules]]
deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
git-tree-sha1 = "dd933c4ef7b4c270aacd4eb88fa64c147492acf0"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.10.0"

[[Distances]]
deps = ["LinearAlgebra", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "3258d0659f812acde79e8a74b11f17ac06d0ca04"
uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
version = "0.10.7"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "80ced645013a5dbdc52cf70329399c35ce007fae"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.13.0"

[[FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[ForwardDiff]]
deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
git-tree-sha1 = "1bd6fc0c344fc0cbee1f42f8d2e7ec8253dda2d2"
uuid = "f6369f11-7733-5829-9624-2563aa707210"
version = "0.10.25"

[[Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[Graphics]]
deps = ["Colors", "LinearAlgebra", "NaNMath"]
git-tree-sha1 = "1c5a84319923bea76fa145d49e93aa4394c73fc2"
uuid = "a2bd30eb-e257-5431-a919-1863eab51364"
version = "1.1.1"

[[Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[ImageBase]]
deps = ["ImageCore", "Reexport"]
git-tree-sha1 = "b51bb8cae22c66d0f6357e3bcb6363145ef20835"
uuid = "c817782e-172a-44cc-b673-b171935fbb9e"
version = "0.1.5"

[[ImageCore]]
deps = ["AbstractFFTs", "ColorVectorSpace", "Colors", "FixedPointNumbers", "Graphics", "MappedArrays", "MosaicViews", "OffsetArrays", "PaddedViews", "Reexport"]
git-tree-sha1 = "9a5c62f231e5bba35695a20988fc7cd6de7eeb5a"
uuid = "a09fc81d-aa75-5fe9-8630-4744c3626534"
version = "0.9.3"

[[ImageIO]]
deps = ["FileIO", "JpegTurbo", "Netpbm", "OpenEXR", "PNGFiles", "QOI", "Sixel", "TiffImages", "UUIDs"]
git-tree-sha1 = "464bdef044df52e6436f8c018bea2d48c40bb27b"
uuid = "82e4d734-157c-48bb-816b-45c225c6df19"
version = "0.6.1"

[[ImageShow]]
deps = ["Base64", "FileIO", "ImageBase", "ImageCore", "OffsetArrays", "StackViews"]
git-tree-sha1 = "d0ac64c9bee0aed6fdbb2bc0e5dfa9a3a78e3acc"
uuid = "4e3cecfd-b093-5904-9786-8bbb286a6a31"
version = "0.3.3"

[[Imath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "87f7662e03a649cffa2e05bf19c303e168732d3e"
uuid = "905a6f67-0a94-5f89-b386-d35d92009cd1"
version = "3.1.2+0"

[[IndirectArrays]]
git-tree-sha1 = "012e604e1c7458645cb8b436f8fba789a51b257f"
uuid = "9b13fd28-a010-5f03-acff-a1bbcff69959"
version = "1.0.0"

[[Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "a7254c0acd8e62f1ac75ad24d5db43f5f19f3c65"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.2"

[[InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[JpegTurbo]]
deps = ["CEnum", "FileIO", "ImageCore", "JpegTurbo_jll", "TOML"]
git-tree-sha1 = "a77b273f1ddec645d1b7c4fd5fb98c8f90ad10a5"
uuid = "b835a17e-a41a-41e7-81f0-2f016b05efe0"
version = "0.1.1"

[[JpegTurbo_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b53380851c6e6664204efb2e62cd24fa5c47e4ba"
uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
version = "2.1.2+0"

[[JuMP]]
deps = ["Calculus", "DataStructures", "ForwardDiff", "JSON", "LinearAlgebra", "MathOptInterface", "MutableArithmetics", "NaNMath", "OrderedCollections", "Printf", "Random", "SparseArrays", "SpecialFunctions", "Statistics"]
git-tree-sha1 = "fe0f87cc077fc6a23c21e469318993caf2947d10"
uuid = "4076af6c-e467-56ae-b986-b466b2749572"
version = "0.22.3"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "e5718a00af0ab9756305a0392832c8952c7426c1"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.6"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[METIS_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "1d31872bb9c5e7ec1f618e8c4a56c8b0d9bddc7e"
uuid = "d00139f3-1899-568f-a2f0-47f597d42d70"
version = "5.1.1+0"

[[MUMPS_seq_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "METIS_jll", "OpenBLAS32_jll", "Pkg"]
git-tree-sha1 = "29de2841fa5aefe615dea179fcde48bb87b58f57"
uuid = "d7ed1dd3-d0ae-5e8e-bfb4-87a502085b8d"
version = "5.4.1+0"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[MappedArrays]]
git-tree-sha1 = "e8b359ef06ec72e8c030463fe02efe5527ee5142"
uuid = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900"
version = "0.4.1"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MathOptInterface]]
deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "JSON", "LinearAlgebra", "MutableArithmetics", "OrderedCollections", "Printf", "SparseArrays", "Test", "Unicode"]
git-tree-sha1 = "e8c9653877adcf8f3e7382985e535bb37b083598"
uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
version = "0.10.9"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MosaicViews]]
deps = ["MappedArrays", "OffsetArrays", "PaddedViews", "StackViews"]
git-tree-sha1 = "b34e3bc3ca7c94914418637cb10cc4d1d80d877d"
uuid = "e94cdb99-869f-56ef-bcf0-1ae2bcbe0389"
version = "0.3.3"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[MutableArithmetics]]
deps = ["LinearAlgebra", "SparseArrays", "Test"]
git-tree-sha1 = "842b5ccd156e432f369b204bb704fd4020e383ac"
uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
version = "0.3.3"

[[NaNMath]]
git-tree-sha1 = "b086b7ea07f8e38cf122f5016af580881ac914fe"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.7"

[[Netpbm]]
deps = ["FileIO", "ImageCore"]
git-tree-sha1 = "18efc06f6ec36a8b801b23f076e3c6ac7c3bf153"
uuid = "f09324ee-3d7c-5217-9330-fc30815ba969"
version = "1.0.2"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OffsetArrays]]
deps = ["Adapt"]
git-tree-sha1 = "043017e0bdeff61cfbb7afeb558ab29536bbb5ed"
uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
version = "1.10.8"

[[OpenBLAS32_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "ba4a8f683303c9082e84afba96f25af3c7fb2436"
uuid = "656ef2d0-ae68-5445-9ca0-591084a874a2"
version = "0.3.12+1"

[[OpenEXR]]
deps = ["Colors", "FileIO", "OpenEXR_jll"]
git-tree-sha1 = "327f53360fdb54df7ecd01e96ef1983536d1e633"
uuid = "52e1d378-f018-4a11-a4be-720524705ac7"
version = "0.3.2"

[[OpenEXR_jll]]
deps = ["Artifacts", "Imath_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "923319661e9a22712f24596ce81c54fc0366f304"
uuid = "18a262bb-aa17-5467-a713-aee519bc75cb"
version = "3.1.1+0"

[[OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Osi_jll]]
deps = ["Artifacts", "CoinUtils_jll", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "OpenBLAS32_jll", "Pkg"]
git-tree-sha1 = "28e0ddebd069f605ab1988ab396f239a3ac9b561"
uuid = "7da25872-d9ce-5375-a4d3-7a845f58efdd"
version = "0.10800.600+0"

[[PNGFiles]]
deps = ["Base64", "CEnum", "ImageCore", "IndirectArrays", "OffsetArrays", "libpng_jll"]
git-tree-sha1 = "eb4dbb8139f6125471aa3da98fb70f02dc58e49c"
uuid = "f57f5aa1-a3ce-4bc8-8ab9-96f992907883"
version = "0.3.14"

[[PaddedViews]]
deps = ["OffsetArrays"]
git-tree-sha1 = "03a7a85b76381a3d04c7a1656039197e70eda03d"
uuid = "5432bcbf-9aad-5242-b902-cca2824c8663"
version = "0.5.11"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "13468f237353112a01b2d6b32f3d0f80219944aa"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.2"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PkgVersion]]
deps = ["Pkg"]
git-tree-sha1 = "a7a7e1a88853564e551e4eba8650f8c38df79b37"
uuid = "eebad327-c553-4316-9ea0-9fa01ccd7688"
version = "0.1.1"

[[PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "8979e9802b4ac3d58c503a20f2824ad67f9074dd"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.34"

[[PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "db3a23166af8aebf4db5ef87ac5b00d36eb771e2"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.4.0"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "2cf929d64681236a2e074ffafb8d568733d2e6af"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.3"

[[PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "dfb54c4e414caa595a1f2ed759b160f5a3ddcba5"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.3.1"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Profile]]
deps = ["Printf"]
uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"

[[ProgressMeter]]
deps = ["Distributed", "Printf"]
git-tree-sha1 = "afadeba63d90ff223a6a48d2009434ecee2ec9e8"
uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
version = "1.7.1"

[[QOI]]
deps = ["ColorTypes", "FileIO", "FixedPointNumbers"]
git-tree-sha1 = "18e8f4d1426e965c7b532ddd260599e1510d26ce"
uuid = "4b34888f-f399-49d4-9bb3-47ed5cae4e65"
version = "1.0.0"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sixel]]
deps = ["Dates", "FileIO", "ImageCore", "IndirectArrays", "OffsetArrays", "REPL", "libsixel_jll"]
git-tree-sha1 = "8fb59825be681d451c246a795117f317ecbcaa28"
uuid = "45858cf5-a6b0-47a3-bbea-62219f50df47"
version = "0.1.2"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "8d0c8e3d0ff211d9ff4a0c2307d876c99d10bdf1"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.1.2"

[[StackViews]]
deps = ["OffsetArrays"]
git-tree-sha1 = "46e589465204cd0c08b4bd97385e4fa79a0c770c"
uuid = "cae243ae-269e-4f55-b966-ac2d0dc13c15"
version = "0.1.1"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "95c6a5d0e8c69555842fc4a927fc485040ccc31c"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.3.5"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StatsAPI]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "c3d8ba7f3fa0625b062b82853a7d5229cb728b6b"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.2.1"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "bb1064c9a84c52e277f1096cf41434b675cd368b"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.1"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[TensorCore]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "1feb45f88d133a655e001435632f019a9a1bcdb6"
uuid = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50"
version = "0.1.1"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TiffImages]]
deps = ["ColorTypes", "DataStructures", "DocStringExtensions", "FileIO", "FixedPointNumbers", "IndirectArrays", "Inflate", "OffsetArrays", "PkgVersion", "ProgressMeter", "UUIDs"]
git-tree-sha1 = "991d34bbff0d9125d93ba15887d6594e8e84b305"
uuid = "731e570b-9d59-4bfa-96dc-6df516fadf69"
version = "0.5.3"

[[TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[libpng_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
version = "1.6.38+0"

[[libsixel_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "78736dab31ae7a53540a6b752efc61f77b304c5b"
uuid = "075b6546-f08a-558a-be8f-8157d0f608a5"
version = "1.8.6+1"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─d44f3f70-705b-43a2-85c0-a441e08fbbe6
# ╟─a12d1135-ce89-414c-866e-14b4415e1990
# ╟─d35c109a-e5b6-4c01-a9c2-333b883783f2
# ╠═2d6d0588-f4e4-4cad-8353-9c0ba020fec6
# ╟─2877f488-1ab3-46da-8a49-7db65c3d8a45
# ╠═cf6d3f77-781d-4cd0-ad45-9c5c58380354
# ╟─11c8323a-9cd5-4580-93d9-ebcd085b4551
# ╠═0a9212cd-4c79-4752-930b-45dce6a1f3e8
# ╟─7f7439d3-3d3b-4098-8be5-ca9aa59373ce
# ╠═3e708ffa-05d9-4cf7-ba4c-11ca5a3e43e7
# ╟─54082ac4-c31d-4d27-9cc3-5248018c6838
# ╠═6bc42cc7-e3d4-4470-bece-6af950280372
# ╟─a5f3abe7-a843-4ea5-8e77-dfd075b243a5
# ╠═ee41a6b3-fb13-4269-ab22-ad2a85bf26e5
# ╟─dd262e8b-4697-4aa7-9bc7-a32ac25c96dc
# ╠═d5fc7ca7-33cd-4f8a-aecf-8897bcd40e15
# ╟─4ab84877-f5d8-4c1b-b0d2-346f0048bb20
# ╠═1054b6a4-1c28-46b6-88a4-6f7cae84edee
# ╟─07f58656-4826-4ce5-8c97-bf5570f0fc60
# ╠═17088a0b-a568-431c-8d01-13f957e2a4c0
# ╟─4ded521e-c41f-4927-854b-86f8912269db
# ╠═dab4e559-9dc6-4a72-8d38-9a0487b6251a
# ╟─aa340dc7-8833-4e40-b63e-c602096e40e3
# ╟─a2fad4bc-66fa-4466-ae29-b495590f035f
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl


text:
export
    Field, f_inner_type, f_from,
    f_zero, f_one, f_add, f_neg, f_sub, f_mul, f_inv, f_div, f_pow,
    PrimeField, RealField, RationalField

# Dynamic PLs can model this more concisely. An example in Python:
# <https://jeremykun.com/2014/03/13/programming-with-finite-fields/>.

#===== Field =====#
abstract type Field{T} end

function f_inner_type(f::Field{T})::Type where T
    T
end

@inline function f_from(f::Field{T}, x::T) where T
    x
end

@inline function f_zero(f::Field{T})::T where T
    zero(T)
end

@inline function f_one(f::Field{T})::T where T
    one(T)
end

@inline function f_add(_f::Field{T}, a::T, b::T) where T
    a + b
end

@inline function f_neg(_f::Field{T}, a::T) where T
    -a
end

@inline function f_sub(f::Field{T}, a::T, b::T) where T
    f_add(f, a, f_neg(f, b))
end

@inline function f_mul(_f::Field{T}, a::T, b::T)::T where T
    a * b
end

@inline function f_inv(_f::Field{T}, a::T) where T
    inv(a)
end

@inline function f_div(f::Field{T}, a::T, b::T) where T
    f_mul(f, a, f_inv(f, b))
end

function f_pow(f::Field{T}, a::T, p::Int) where T
    acc = f_one(f)
    while p != 0
        if p & 1 != 0
            acc = f_mul(f, acc, a)
        end
        a = f_mul(f, a, a)
        p >>= 1
    end
    acc
end

#===== Prime Field =====#
"GF(p)"
struct PrimeField{T<:Integer} <: Field{T}
    p::T
end

function f_from(f::PrimeField{T}, x::T) where T
    mod(x, f.p)
end

function f_add(f::PrimeField{T}, a::T, b::T)::T where T
    mod(a + b, f.p)
end

function f_neg(f::PrimeField{T}, a::T)::T where T
    mod(-a, f.p)
end

function f_mul(f::PrimeField{T}, a::T, b::T)::T where T
    # Uses Egyptian multiplication to prevent overflow.
    res = zero(T)
    m = f.p
    while a != 0
        if a & 1 != 0
            if b >= m - res
                res -= m
            end
            res += b
        end
        a >>= 1

        temp_b = b
        if b >= m - b
            temp_b -= m
        end
        b += temp_b
    end
    res
end

function f_inv(f::PrimeField{T}, a::T)::T where T
    m = f.p
    mod(gcdx(a, m)[2], m)
end

function f_pow(f::PrimeField{T}, a::T, p::Int) where T
    powermod(a, p, f.p)
end

#===== Real Field =====#
struct RealField{T<:AbstractFloat} <: Field{T}
end

#===== Rational Field =====#
struct RationalField{T<:Rational} <: Field{T}
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
function StaticArrays._one(s::Size, ::Type{SM}) where {Sel, N, SM <: HybridArrays.SSubArray{Tuple{N,N}, Sel}}
    return StaticArrays._one(s, SMatrix{N,N,Sel})
end

Base.one(A::HA) where {HA<:HybridMatrix} = HA(one(A.data))

@inline function Base.zero(a::HA) where {S, Sel, HA <: SSubArray{S, Sel}}
    return StaticArrays.zeros(SArray{S, Sel})
end

Base.fill!(A::HybridArray, x) = fill!(A.data, x)

function Base.zero(a::HybridArray{S}) where {S}
    return HybridArray{S}(zero(a.data))
end


text:
# compare exceptions to returning a type

function run_exception(n)
    count = 0
    for i in 1:n
        try
            random_exception()
        catch
            count += 1
        end
    end
    println("exceptions: $count")
end

function random_exception()
    if rand(1:2) == 1
        error()
    end
end

function run_type(n)
    count = 0
    for i in 1:n
        if isa(random_type(), A)
            count += 1
        end
    end
    println("types: $count")
end

type A end
type B end

function random_type()
    if rand(1:2) == 1
        A()
    else
        B()
    end
end

run_exception(10)
run_type(10)
@time run_exception(1000)
@time run_type(1000)
@time run_exception(1000)
@time run_type(1000)

# exceptions are much slower

# exceptions: 521
#  33.651 milliseconds (703 allocations: 20269 bytes)
# types: 519
# 225.958 microseconds (31 allocations: 1248 bytes)
# exceptions: 518
#  33.506 milliseconds (549 allocations: 9552 bytes)
# types: 498
# 205.418 microseconds (30 allocations: 1232 bytes)


text:
module ITF1788tests

include("add.jl")

export regenerate_test

end # module


text:
using Test, ImageFeatures, Images

@testset "GLCM" begin
    img = [ 0 0 1 1
            0 0 1 1
            0 2 2 2
            2 2 3 3 ]
    img_gray = Gray.(img / 3)

    glcm_mat = glcm(img, 1, 0, 4)
    expected_1 = [ 2 2 1 0
                   0 2 0 0
                   0 0 3 1
                   0 0 0 1 ]
    @test all(expected_1 .== glcm_mat)
    glcm_mat = glcm(img_gray, 1, 0, 4)
    @test all(expected_1 .== glcm_mat)
    glcm_mat = glcm_norm(img, 1, 0, 4)
    @test all(expected_1 / sum(expected_1) .== glcm_mat)
    glcm_mat = glcm_symmetric(img, 1, 0, 4)
    @test all(expected_1 + expected_1' .== glcm_mat)

    glcm_mat = glcm(img, 1, pi/2, 4)
    expected_2 = [ 3 0 2 0
                   0 2 2 0
                   0 0 1 2
                   0 0 0 0 ]
    @test all(expected_2 .== glcm_mat)

    glcm_mats = glcm(img, [1, 1], [0, pi/2], 4)
    @test all(glcm_mats[1, 1] .== expected_1)
    @test all(glcm_mats[1, 2] .== expected_2)
    @test all(glcm_mats[2, 1] .== expected_1)
    @test all(glcm_mats[2, 2] .== expected_2)

    glcm_mats = glcm(img_gray, [1, 1], [0, pi/2], 4)
    @test all(glcm_mats[1, 1] .== expected_1)
    @test all(glcm_mats[1, 2] .== expected_2)
    @test all(glcm_mats[2, 1] .== expected_1)
    @test all(glcm_mats[2, 2] .== expected_2)

    glcm_mats = glcm_symmetric(img, [1, 1], [0, pi/2], 4)
    @test all(glcm_mats[1, 1] .== expected_1 + expected_1')
    @test all(glcm_mats[1, 2] .== expected_2 + expected_2')
    @test all(glcm_mats[2, 1] .== expected_1 + expected_1')
    @test all(glcm_mats[2, 2] .== expected_2 + expected_2')

    glcm_mats = glcm_norm(img, [1, 1], [0, pi/2], 4)
    @test all(glcm_mats[1, 1] .== expected_1 / sum(expected_1))
    @test all(glcm_mats[1, 2] .== expected_2 / sum(expected_2))
    @test all(glcm_mats[2, 1] .== expected_1 / sum(expected_1))
    @test all(glcm_mats[2, 2] .== expected_2 / sum(expected_2))

    glcm_mat = glcm(img, 2, 0, 4)
    expected_3 = [ 0 4 1 0
                   0 0 0 0
                   0 0 1 2
                   0 0 0 0 ]
    @test all(expected_3 .== glcm_mat)

    glcm_mat = glcm(img, 1, 0, 2)
    expected_4 = [ 6 1
                   0 5 ]
    @test all(expected_4 .== glcm_mat)

    expected_5 = [ 1 0 3 0
                   0 0 2 2
                   0 0 0 0
                   0 0 0 0 ]
    glcm_mats = glcm(img, 2, [0, pi/2], 4)
    @test all(glcm_mats[1] .== expected_3)
    @test all(glcm_mats[2] .== expected_5)

    glcm_mats = glcm_symmetric(img, 2, [0, pi/2], 4)
    @test all(glcm_mats[1] .== expected_3 + expected_3')
    @test all(glcm_mats[2] .== expected_5 + expected_5')

    glcm_mats = glcm_norm(img, 2, [0, pi/2], 4)
    @test all(glcm_mats[1] .== expected_3 / sum(expected_3))
    @test all(glcm_mats[2] .== expected_5 / sum(expected_5))

    glcm_mats = glcm(img, [1, 2], 0, 4)
    @test all(glcm_mats[1] .== expected_1)
    @test all(glcm_mats[2] .== expected_3)

    glcm_mats = glcm_symmetric(img, [1, 2], 0, 4)
    @test all(glcm_mats[1] .== expected_1 + expected_1')
    @test all(glcm_mats[2] .== expected_3 + expected_3')

    glcm_mats = glcm_norm(img, [1, 2], 0, 4)
    @test all(glcm_mats[1] .== expected_1 / sum(expected_1))
    @test all(glcm_mats[2] .== expected_3 / sum(expected_3))
end

@testset "Properties" begin
    img = convert(Array{Int}, reshape(1:1:30, 5, 6))
    @test glcm_prop(img, max_prob) == maxfinite(img)
    @test glcm_prop(img, contrast) == 2780
    @test glcm_prop(img, dissimilarity) == 930
    @test glcm_prop(img, ASM) == 9455
    @test isapprox(glcm_prop(img, energy), 97.2368, rtol = 0.001)
    @test isapprox(glcm_prop(img, glcm_entropy), -1357.0889, rtol = 0.001)
    @test glcm_prop(img, glcm_mean_ref) == 1455
    @test isapprox(glcm_prop(img, glcm_var_ref), 31307.95505, rtol = 0.001)
    @test glcm_prop(img, glcm_mean_neighbour) == 2065
    @test isapprox(glcm_prop(img, glcm_var_neighbour), 44433.61666, rtol = 0.001)
    @test isapprox(glcm_prop(img, correlation), 0.99999, rtol = 0.001)
    @test isapprox(glcm_prop(img, IDM), 165.5176, rtol = 0.001)

    glcm_props = glcm_prop(img, 3, max_prob)
    expected = [  7.0  12.0  17.0  22.0  27.0  27.0
                  8.0  13.0  18.0  23.0  28.0  28.0
                  9.0  14.0  19.0  24.0  29.0  29.0
                 10.0  15.0  20.0  25.0  30.0  30.0
                 10.0  15.0  20.0  25.0  30.0  30.0 ]
    @test all(glcm_props .== expected)
    glcm_props = glcm_prop(img, 4, 3, max_prob)
    expected = [  8.0  13.0  18.0  23.0  28.0  28.0
                  9.0  14.0  19.0  24.0  29.0  29.0
                 10.0  15.0  20.0  25.0  30.0  30.0
                 10.0  15.0  20.0  25.0  30.0  30.0
                 10.0  15.0  20.0  25.0  30.0  30.0 ]
    @test all(glcm_props .== expected)

    # Test Correlation when Variance is Zero

    glcm_mat = [ 1 0 0 0
                 0 0 0 0
                 0 0 0 0
                 0 0 0 0 ]
    @test glcm_prop(glcm_mat, correlation) == 1
end


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

x="println(\"x=\$(repr(x))\\n\$x\")"
println("x=$(repr(x))\n$x")

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
using FMIndexes
using Combinatorics
using Base.Test

srand(12345)

# simple DNA sequence type
@enum Nuc A C G T
type DNASeq
    data::Vector{Nuc}
end
Base.getindex(seq::DNASeq, i::Integer) = seq.data[i]
Base.length(seq::DNASeq) = length(seq.data)
Base.endof(seq::DNASeq)  = length(seq.data)

@testset "construct" begin
    @testset "one" begin
        σ = 2
        seq = [0x00]
        index = FMIndex(seq, σ)
        @test typeof(index) == FMIndex{1,UInt8}
    end

    @testset "short" begin
        σ = 4
        seq = [0x00, 0x01, 0x02, 0x03]
        index = FMIndex(seq, σ)
        @test typeof(index) == FMIndex{2,UInt8}
    end

    @testset "long" begin
        σ = 2
        seq = rand(0x00:0x01, 2^24)
        index = FMIndex(seq, σ)
        @test typeof(index) == FMIndex{1,UInt32}
    end

    @testset "dna" begin
        σ = 4
        seq = DNASeq([A, C, G, T, A, T])
        index = FMIndex(seq, σ)
        @test typeof(index) == FMIndex{2,UInt8}
    end

    @testset "boundaries" begin
        σ = 4
        # 8 bits
        n = 2^8
        seq = rand(0x00:0x03, n)
        index = FMIndex(seq, σ)
        @test typeof(index) == FMIndex{2,UInt8}
        # 16 bits
        n = 2^16
        seq = rand(0x00:0x03, n)
        index = FMIndex(seq, σ)
        @test typeof(index) == FMIndex{2,UInt16}
        # 32 bits
        # too long to test!
    end

    @testset "mmap" begin
        σ = 4
        seq = rand(0x00:0x03, 2^10)
        index = FMIndex(seq, σ, mmap=true)
        @test typeof(index) == FMIndex{2,UInt16}
    end

    @testset "use pSAscan" begin
        σ = 4
        seq = rand(0x00:0x03, 2^24)
        if haskey(ENV, "PSASCAN")
            index = FMIndex(seq, σ, program=:psascan, psascan=ENV["PSASCAN"], mmap=true)
            @test typeof(index) == FMIndex{2,UInt32}
        else
            info("Skipped a test")
            #@pending typeof(index) --> FMIndex{2,UInt32}
        end
    end
end

@testset "restore" begin
    @testset "examples" begin
        σ = 2
        seq = [0x01]
        index = FMIndex(seq, σ)
        @test restore(index) == seq

        σ = 2
        seq = [0x00, 0x00, 0x01, 0x01]
        index = FMIndex(seq, σ)
        @test restore(index) == seq

        seq = [0x01, 0x01, 0x00, 0x00]
        index = FMIndex(seq, σ)
        @test restore(index) == seq

        σ = 128
        seq = Vector{UInt8}("abracadabra")
        index = FMIndex(seq, σ)
        @test restore(index) == seq
    end

    @testset "dna" begin
        σ = 4
        seq = DNASeq([A, C, G, T, A, T])
        index = FMIndex(seq, σ)
        @test restore(index) == [0x00, 0x01, 0x02, 0x03, 0x0, 0x03]
    end

    @testset "random" begin
        σ = 128
        for n in [2, 15, 100], _ in 1:100
            seq = randstring(n)
            index = FMIndex(Vector{UInt8}(seq), σ)
            @test restore(index) == Vector{UInt8}(seq)
        end
    end
end

@testset "count" begin
    @testset "[0x00]" begin
        σ = 2
        index = FMIndex([0x00], σ)
        @test count([0x00], index) == 1
        @test count([0x01], index) == 0
    end

    @testset "\"abracadabra\"" begin
        σ = 128
        seq = Vector{UInt8}("abracadabra")
        index = FMIndex(seq, σ)

        @test count("a", index) == 5
        @test count("ab", index) == 2
        @test count("abr", index) == 2
        @test count("abra", index) == 2
        @test count("abrac", index) == 1

        @test count("d", index) == 1
        @test count("da", index) == 1
        @test count("dab", index) == 1

        @test count("r", index) == 2
        @test count("ra", index) == 2
        @test count("rac", index) == 1

        # not existing
        @test count("x", index) == 0
        @test count("ax", index) == 0
        @test count("aa", index) == 0
        @test count("aaa", index) == 0
        @test count("braa", index) == 0
    end

    @testset "[0x00, 0x00, 0x01, 0x01]" begin
        σ = 2
        seq = [0x00, 0x00, 0x01, 0x01]
        index = FMIndex(seq, σ)

        @test count([0x00], index) == 2
        @test count([0x00, 0x00], index) == 1
        @test count([0x00, 0x00, 0x01], index) == 1
        @test count([0x00, 0x01, 0x01], index) == 1
        @test count([0x00, 0x00, 0x01, 0x01], index) == 1

        @test count([0x01], index) == 2
        @test count([0x01, 0x01], index) == 1

        @test count([0x01, 0x00], index) == 0
        @test count([0x01, 0x01, 0x00], index) == 0
        @test count([0x00, 0x01, 0x01, 0x00], index) == 0
        @test count([0x01, 0x01, 0x00, 0x00], index) == 0
    end
end

function linear_search(query, seq)
    locs = Int[]
    for i in 1:endof(seq)-length(query)+1
        j = 1
        while j ≤ length(query) && query[j] == seq[i+j-1]
            j += 1
        end
        if j > length(query)
            push!(locs, i)
        end
    end
    return locs
end

@testset "locate/locateall" begin
    @testset "[0x00]" begin
        σ = 2
        index = FMIndex([0x00], σ)
        @test locateall([0x00], index) == [1]
        @test isempty(locateall([0x01], index))
    end

    @testset "\"abracadabra\"" begin
        σ = 128
        seq = Vector{UInt8}("abracadabra")
        index = FMIndex(seq, σ)

        @test locate("a", index) |> collect |> sort == [1, 4, 6, 8, 11]

        @test locateall("a", index) |> sort == [1, 4, 6, 8, 11]
        @test locateall("ab", index) |> sort == [1, 8]
        @test locateall("abr", index) |> sort == [1, 8]
        @test locateall("abra", index) |> sort == [1, 8]
        @test locateall("abrac", index) |> sort == [1]

        @test locateall("d", index) |> sort == [7]
        @test locateall("da", index) |> sort == [7]
        @test locateall("dab", index) |> sort == [7]

        @test locateall("r", index) |> sort == [3, 10]
        @test locateall("ra", index) |> sort == [3, 10]
        @test locateall("rac", index) |> sort == [3]

        @test isempty(locate("x", index))

        @test isempty(locateall("x", index))
        @test isempty(locateall("ax", index))
        @test isempty(locateall("aa", index))
        @test isempty(locateall("aaa", index))
        @test isempty(locateall("braa", index))
    end

    @testset "random" begin
        σ = 4
        seq = rand(0x00:0x03, 10_000)
        index = FMIndex(seq, σ)
        for m in [1, 2, 3, 5, 10]
            query = rand(0x00:0x03, m)
            @test locateall(query, index) |> sort == linear_search(query, seq)
        end
    end

    @testset "boundaries" begin
        σ = 4
        for n in [2^8, 2^16]
            seq = rand(0x00:0x03, n)
            index = FMIndex(seq, σ)
            for set in ([0x00],
                        [0x03],
                        [0x00, 0x00],
                        [0x00, 0x01],
                        [0x02, 0x03],
                        [0x00, 0x01, 0x02],
                        [0x01, 0x02, 0x03],
                        [0x00, 0x01, 0x02, 0x03])
                for query in permutations(set)
                    @test locateall(query, index) |> sort == linear_search(query, seq)
                end
            end
        end
    end
end

@testset "full-text search" begin
    function linear_search(query)
        locs = Int[]
        loc = 0
        while (loc = searchindex(text, Vector{UInt8}(query), loc + 1)) > 0
            push!(locs, loc)
        end
        return locs
    end

    text = open(read, joinpath(dirname(@__FILE__), "lorem_ipsum.txt"))
    index = FMIndex(text, r=2)

    @test count("Lorem", index) == 1
    @test locateall("Lorem", index) == [1]
    @test count("hoge", index) == 0
    @test isempty(locateall("hoge", index))
    for query in ["a", ".",
                  "ex", "In",
                  "non", "Sed",
                  "odio", "Cras",
                  "sollicitudin"]
        locs = linear_search(query)
        @test count(query, index) == length(locs)
        @test locateall(query, index) |> sort == locs
    end

    @test String(restore(index)) == String(text)
end


text:
#---------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Model of a tank basic
#*----------------------------------------------------------------------
#*
#*   Description:
#*       Generic model for a dynamic tank.
#*
#*   Assumptions:
#*		* single- and two-phases involved
#*       * dynamic
#*
#*----------------------------------------------------------------------
#* Author: Rodolfo Rodrigues
#* $Id$
#*--------------------------------------------------------------------
type tank_basic
	tank_basic()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External physical properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of components",
				:Default=>1
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet stream",
				:PosX=>0,
				:PosY=>0,
				:Symbol=>"_{in}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Intermediary outlet stream",
				:Symbol=>"_{outm}"
			)),
			vol_tank(Dict{Symbol,Any}(
				:Brief=>"Routine to volume tank calculation",
				:Symbol=>"_{tank}"
			)),
			fill(mol(Dict{Symbol,Any}(
				:Brief=>"Component molar holdup"
			)),(NComp)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Total component molar holdup"
			)),
			energy(Dict{Symbol,Any}(
				:Brief=>"Internal energy"
			)),
			heat_rate(Dict{Symbol,Any}(
				:Brief=>"Reactor duty",
				:Default=>0
			)),
			[
				:(diff(M) = Inlet.F*Inlet.z - Outletm.F*Outletm.z),
				:(M = Mt*Outletm.z),
				:(sum(Outletm.z) = 1),
				:(diff(E) = Inlet.F*Inlet.h - Outletm.F*Outletm.h + Q),
			],
			[
				"Component molar balance","Component molar","Mole fraction normalisation","Energy balance",
			],
			[:PP,:NComp,],
			[:Inlet,:Outletm,:Tank,:M,:Mt,:E,:Q,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Inlet::stream
	Outletm::stream
	Tank::vol_tank
	M::Array{mol}
	Mt::mol
	E::energy
	Q::heat_rate
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export tank_basic
function setEquationFlow(in::tank_basic)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
end
function atributes(in::tank_basic,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Brief]="Basic model for a dynamic tank"
	drive!(fields,_)
	return fields
end
tank_basic(_::Dict{Symbol,Any})=begin
	newModel=tank_basic()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
include("Util_single.jl")

@everywhere module ReversibleJumpSingleParticle
using StatsBase
using KernelDensity
using Distributions
using SingleParticle

export bounds_check, judge, MCstep, MCstep_parallel, update_stepsize
export ReversibleParticle, init_particle

mutable struct ReversibleParticle
    x_present::Array{Float64, 1}
    obj_func::Function # return relative weight
    obj_func_val::Float64
    weight::Float64
    cov::Array{Float64, 2}
    model_id::Int
end

function init_particle(lbs, ubs, obj_func, weight, model_id)
    """
    obj_func: must return -1 as a flag if has an error.
    """
    @assert length(lbs) == length(ubs)
    num_param = length(lbs)
    obj_func_val = -1
    kernel = Uniform(0, 10)
    x_present = zeros(length(lbs))   # generate parameter
    while obj_func_val < 0
        for i in 1:num_param
            x_present[i] = rand(kernel)
        end
        obj_func_val = obj_func(x_present)
    end
    return ReversibleParticle(x_present, obj_func, obj_func_val, weight, eye(num_param), model_id)
end

end


text:
cd(@__DIR__); include("setups/distROT_distTSD_ratio.jl")
gr(dpi = 200)

## (a) P64
include("setups/path64.jl")
ρ1 = generate_ROT_TSD_ratio(nsim, 𝚽, ∇𝚽, 𝛌, Q)
# plt = ROT_TSD_ratio_histogram(ρ1)
# savefig(plt, "../figs/Path64_ROT_TSD.png")

## (b) P₇ x P₃
include("setups/grid7x3.jl")
ρ2 = generate_ROT_TSD_ratio(nsim, 𝚽, ∇𝚽, 𝛌, Q)
# plt = ROT_TSD_ratio_histogram(ρ2)
# savefig(plt, "../figs/Grid7x3_ROT_TSD.png")

## (c) Erdos Rényi
include("setups/er.jl")
ρ3 = generate_ROT_TSD_ratio(nsim, 𝚽, ∇𝚽, 𝛌, Q)
# plt = ROT_TSD_ratio_histogram(ρ3)
# savefig(plt, "../figs/Erdos_Renyi_ROT_TSD.png")

## (d) weighted RGC100
include("setups/rgc100.jl")
ρ4 = generate_ROT_TSD_ratio(nsim, 𝚽, ∇𝚽, 𝛌, Q; edge_length = edge_length)
# plt = ROT_TSD_ratio_histogram(ρ4)
# savefig(plt, "../figs/wRGC100_ROT_TSD.png")

## boxplot
plt = boxplot(["(a) Path" "(b) Grid" "(c) ER" "(d) RGC100"], [ρ1, ρ2, ρ3, ρ4];
    legend = false, frame = :box, ylim = [0.9, 1.9], tickfontsize = 11,
    outliers = true, grid = false, range = 3, lw = 1, size = (800, 600),
    ylab = "ρ", yguidefontsize = 14, xtickfontsize = 12, color = :white)
savefig(plt, "../figs/ROT_TSD_boxplots.png")

## Table 4.1
display_basic_stats([ρ1, ρ2, ρ3, ρ4])


text:
using NeuralNetworks

net = Serial(Float64, (28 * 28, 1), (10, 1))
add_layer!(net, Linear(Float64, 1, 28 * 28, 300))
add_layer!(net, Sigmoid(Float64, (300, 1)))
add_layer!(net, Linear(Float64, 1, 300, 10))
add_layer!(net, Sigmoid(Float64, (10, 1)))
validate(net)

initWeights!(net)

const num_iters = 1000
const num_images = 60000

using MNIST

target = zeros(Float64, 10, 1)
criterion = MSECriterion(Float64, 1, 10)

sgd_state = ObjectIdDict()
function sgd(x, dfdx)
    local lr = 0.1
    local mom = 0.9
    local cur_dfdx
    if haskey(sgd_state, x)
        cur_dfdx = sgd_state[x]
    else
        cur_dfdx = zeros(x)
    end
    cur_dfdx[:] = mom * cur_dfdx + (1.0 - mom) * dfdx
    x[:] = x - lr * cur_dfdx
end

for iter = 1:num_iters
    local loss
    loss = 0.0
    for img = 1:num_images
        local image_data = reshape(MNIST.trainfeatures(img), 28 * 28, 1)
        local label = Int(MNIST.trainlabel(img)) + 1
        target[label,1] = 1.0

        local out = forward!(net, image_data)
        local l = forward!(criterion, out, target)[1]
        loss += l
        local gi = backward!(criterion, out, target)
        zeroGradParameters!(net)
        backward!(net, image_data, gi)

        target[label,1] = 0.0

        for (x, dfdx) in getParameters(net)
            sgd(x, dfdx)
        end
        println("  ", l)
    end
    loss /= num_images
    println(iter, " ", loss)
end


text:
__precompile__()
module MD5

export md5

mutable struct MD5_CTX
    state::Vector{UInt32}
    bytecount::UInt64
    buffer::Vector{UInt8}
    M::Vector{UInt32}
end

function md5(data::Array{UInt8,1})
    ctx = MD5_CTX()
    update!(ctx, data)
    return digest!(ctx)
end

md5(str::AbstractString) = md5(convert(Array{UInt8,1}, str))
md5(io::IO) = md5(read(io))

digestlen(::Type{MD5_CTX}) = 16
state_type(::Type{MD5_CTX}) = UInt32
blocklen(::Type{MD5_CTX}) = UInt64(64)

const MD5_initial_hash_value = UInt32[
    0x67452301,
    0xefcdab89,
    0x98badcfe,
    0x10325476,]

const S_MD5 = UInt32[
    7, 12, 17, 22,  7, 12, 17, 22,  7, 12, 17, 22,  7, 12, 17, 22,
    5,  9, 14, 20,  5,  9, 14, 20,  5,  9, 14, 20,  5,  9, 14, 20,
    4, 11, 16, 23,  4, 11, 16, 23,  4, 11, 16, 23,  4, 11, 16, 23,
    6, 10, 15, 21,  6, 10, 15, 21,  6, 10, 15, 21,  6, 10, 15, 21,]

@assert length(S_MD5) == 64
@assert S_MD5[33] == 4
@assert last(S_MD5) == 21

MD5_CTX() = MD5_CTX(copy(MD5_initial_hash_value), 0, zeros(UInt8, blocklen(MD5_CTX)),
                    zeros(UInt32,16))

import Base.copy
copy(ctx::MD5_CTX) = MD5_CTX(copy(ctx.state), ctx.bytecount, copy(ctx.buffer))

include("sha.jl")

const K_MD5 = UInt32[floor(UInt32, 2^32 * abs(sin(i))) for i in 1:64]
@assert last(K_MD5) == 0xeb86d391
@assert K_MD5[29] == 0xa9e3e905

# transform!(ctx::MD5_CTX) = transform_baseline!(ctx)
transform!(ctx::MD5_CTX) = transform_unrolled!(ctx)

@generated function transform_unrolled!(context::MD5_CTX)
    ret = quote
        pbuf = Ptr{UInt32}(pointer(context.buffer))
        M = context.M
    end
    for i in 1:16
        ex = :(M[$i] = unsafe_load(pbuf,$i))
        push!(ret.args, ex)
    end
    ex  = quote
        A = context.state[1]
        B = context.state[2]
        C = context.state[3]
        D = context.state[4]
    end
    push!(ret.args, ex)
    for i in 0:63
        if 0 ≤ i ≤ 15
            ex = :(F = (B & C) | ((~B) & D))
            g = i
        elseif 16 ≤ i ≤ 31
            ex = :(F = (D & B) | ((~D) & C))
            g = 5i + 1
        elseif 32 ≤ i ≤ 47
            ex = :(F = B ⊻ C ⊻ D)
            g = 3i + 5
        elseif 48 ≤ i ≤ 63
            ex = :(F = C ⊻ (B | (~D)))
            g = 7i
        end
        push!(ret.args, ex)
        g = (g % 16) + 1
        ex = quote
            temp = D
            D = C
            C = B
            inner = A + F + $(K_MD5[i+1]) + M[$g]
            rot_inner = lrot($(S_MD5[i+1]), inner, 32)
            B = B + rot_inner
            A = temp
        end
        push!(ret.args, ex)
    end

    ex = quote
        context.state[1] += A
        context.state[2] += B
        context.state[3] += C
        context.state[4] += D
    end
    push!(ret.args, ex)
    quote
        @inbounds $ret
    end
end

function transform_baseline!(context::MD5_CTX)
    pbuf = Ptr{UInt32}(pointer(context.buffer))
    for i in 1:16
        context.M[i] = unsafe_load(pbuf,i)
    end
    A = context.state[1]
    B = context.state[2]
    C = context.state[3]
    D = context.state[4]
    for i in 0:63
        if 0 ≤ i ≤ 15
            F = (B & C) | ((~B) & D)
            g = i
        elseif 16 ≤ i ≤ 31
            F = (D & B) | ((~D) & C)
            g = (5i + 1)
        elseif 32 ≤ i ≤ 47
            F = xor(B,C,D)
            g = (3i + 5)
        elseif 48 ≤ i ≤ 63
            F = C ⊻ (B | (~D))
            g = (7i)
        end
        g = g % 16

        temp = D
        D = C
        C = B
        inner = A + F + K_MD5[i+1] + context.M[g+1]
        rot_inner = lrot(S_MD5[i+1], inner, 32)
        B = B + rot_inner
        A = temp
    end

    context.state[1] += A
    context.state[2] += B
    context.state[3] += C
    context.state[4] += D
end

function digest!(context::MD5_CTX)
    T = typeof(context)
    usedspace = context.bytecount % blocklen(T)
    usedspace += 1
    context.buffer[usedspace] = 0x80
    if usedspace <= 56
        while usedspace < 56
            usedspace += 1
            context.buffer[usedspace] = 0x00
        end
    else
        context.buffer[usedspace+1:end] = 0x00
        transform!(context)
        fill!(context.buffer, 0x00)
        usedspace = 56
    end
    @assert usedspace == 56
    B = typeof(context.bytecount)
    bitcount = context.bytecount * B(8)
    pbuf = Ptr{B}(pointer(context.buffer))
    index = 8
    unsafe_store!(pbuf, bitcount, index)
    transform!(context)
    reinterpret(UInt8, context.state)
end

end # module


text:
"""
    ecdf(y::Vector{<:Real})::Tuple{Vector{<:Real}, Vector{<:Real}}

Compute the empirical cumulative distribution function using the Gumbel formula.

The empirical quantiles are computed using the Gumbel plotting positions as
as recommended by [Makkonen (2006)](https://journals.ametsoc.org/jamc/article/45/2/334/12668/Plotting-Positions-in-Extreme-Value-Analysis).

# Example
```julia-repl
julia> (x, F̂) = Extremes.ecdf(y)
```
# Reference
Makkonen, L. (2006). Plotting positions in extreme value analysis. Journal of
Applied Meteorology and Climatology, 45(2), 334-340.
"""
function ecdf(y::Vector{<:Real})::Tuple{Vector{<:Real}, Vector{<:Real}}
    ys = sort(y)
    n = length(ys)
    p = collect(1:n)/(n+1)

    return ys, p
end

"""
    pp_analysis(sepp::SEPP)

Compute the rescaled time transformed of the realizations of the point process.

It should behave like order statistics from a standard uniform distribution.
"""
function pp_analysis(sepp::SEPP)
    ts=sepp.data
    isnothing(ts) && error("No data in model")

    mts = MarkedTimeSeries(ts)
    marks = mts.marks
    times = first(mts.times) isa TimeType ? Dates.value.(Date.(mts.times)) : mts.times

    starttime = first(times)

    μ = sepp.μ
    ϕ = sepp.ϕ
    γ = sepp.γ
    δ = sepp isa SEMPPExpKern ? sepp.δ : 0
    impact_func = sepp isa SEMPPExpKern ? sepp.impact_function : (_ -> 1)

    aux(i) = sum(impact_func.(δ .* marks[1:i]) .* (1 .- exp.(-γ .*(times[i] .- times[1:i]))))

    τs = μ .* (times .- starttime) .+ ϕ/γ .* aux.(collect(1:length(times)))

    scaled = τs ./ (last(τs))

    n = length(scaled)

    return scaled, (1/n):(1/n):1
end

"""
    transformed_marks_ecdf(sempp::SEMPPExpKern)

Compute the transformed marks ecdf for ploting.

See Li2020 4.2..
"""
function transformed_marks_ecdf(sempp::SEMPPExpKern)
    mts = sempp.data
    isnothing(mts) && error("No data in model")

    times = mts.times
    marks = mts.marks

    γ = sempp.γ
    δ = sempp.δ
    impact_func = sempp.impact_function
    markdens = sempp.markdens
    ξ = sempp.ξ
    α = sempp.α
    β = sempp.β
    κ = sempp.κ

    vol = volfunc(times, mts, γ, δ, impact_func)
    σ = β .+ α .* vol

    sigmarks = hcat(σ, marks)

    if markdens == Distributions.GeneralizedPareto
        u = (sigmark -> cdf(Distributions.GeneralizedPareto(0, sigmark[1], ξ), sigmark[2])).(eachrow(sigmarks))
    else    # EGPD case
        u = (sigmark -> cdf(EGPD.EGPpower(sigmark[1], ξ, κ), sigmark[2])).(eachrow(sigmarks))
    end

    return ecdf(u)
end

"""
    marks_unit_exponential_qq(sempp::SEMPPExpKern)
Compute the quantile plot data based on the unit exponential distribution.

See 4.2 in Li2020.
"""
function marks_unit_exponential_qq(sempp::SEMPPExpKern)

    u, p = transformed_marks_ecdf(sempp)

    model_qs = -log.(1 .- p)

    emp_qs = -log.(1 .- u)

    return emp_qs, model_qs
end


text:
let
    global reader, active_reader_id, start_reading, stop_reading, read_wav
    _readers::Dict{String, Union{Base.Process,PyObject,IOBuffer}}               = Dict{String, Union{Base.Process,PyObject,IOBuffer}}()
    _active_reader_id::String                                                   = ""
    reader(id::String=DEFAULT_READER_ID)::Union{Base.Process,PyObject,IOBuffer} = _readers[id]
    active_reader_id()::String                                                  = (if (_active_reader_id=="") @APIUsageError("no reader is active.") end; return _active_reader_id)

    function start_reading(filename::String; id::String=DEFAULT_READER_ID) # NOTE: here could be started multiple readers.
        _readers[id] = Wave.open(filename, "rb")
        check_wav(_readers[id])
        _active_reader_id = id
        return _readers[id]
    end

    function start_reading(audio::AbstractVector{UInt8}; id::String=DEFAULT_READER_ID) # NOTE: here could be started multiple readers.
        _readers[id] = IOBuffer(audio)
        _active_reader_id = id
        return _readers[id]
    end

    function start_reading(audio_input_cmd::Cmd; id::String=DEFAULT_READER_ID) # NOTE: here could be started multiple readers.
        _readers[id] = open(audio_input_cmd)
        _active_reader_id = id
        return _readers[id]
    end

    function stop_reading(; id::String=DEFAULT_READER_ID)
        close(_readers[id])
    end

    function check_wav(reader::PyObject)
        if reader.getnchannels() != 1 || reader.getsampwidth() != 2 || reader.getcomptype() != "NONE" || reader.getframerate() != SAMPLERATE
            @FileError("the audio file must be WAV format mono PCM @ $SAMPLERATE Hz.")
        end
    end

    function read_wav(filepath::String)
        wav_reader = start_reading(filepath)
        audio = zeros(UInt8, wav_reader.getnframes()*sizeof(AUDIO_ELTYPE))
        bytes_read = readbytes!(wav_reader, audio)
        close(wav_reader)
        if (bytes_read == 0) @FileError("file $filepath could not be read.") end
        return audio
    end
end


text:
SQ._collect{R<:Relation}(_src::Tuple{Grouped{R}}, q::SQ.Node{SQ.Summarize}) =
    Relations._collect(first(_src), first(_src).src, q)

function _collect{R<:Relation,S<:Attributed,T}(gr::Grouped{R}, src::Relation{S,T}, q::SQ.Node{SQ.Summarize})
    res, fs, gs, _res = pre_apply(gr, q)
    grp_levels, grp_indices = gr.group_levels, gr.group_indices
    for level in grp_levels
        empty!(_res)
        indices = grp_indices[level]
        grow_nonnull!(_res, fs, src, indices)
        grow_summaries!(res, gs, _res, level)
    end
    return res
end

function pre_apply(gr::Grouped, q::SQ.Node{SQ.Summarize})
    _cols = Vector{Any}()
    _attrs = [ d.res_field for d in SQ.dos(q) ]
    fs = tuple([ d.f for d in SQ.dos(q) ]...)
    nrow = length(gr.src)
    T = typeof(gr.src).parameters[2]
    for f in fs
        U = Core.Inference.return_type(f, Tuple{T})
        if isleaftype(U)
            U <: Nullable ? push!(_cols, Vector{eltype(U)}()) :
                            push!(_cols, Vector{U}())
        # else TODO: make this a real code path for when inferred type is non-concrete
        end
    end
    _res = Relation(tuple(_attrs...), _cols...)
    # attributes, columns of res include groupbys and summaries
    cols = Any[ similar(getfield(gr.src.src, groupby), 0) for groupby in gr.groupbys ]
    attrs = tuple(vcat(gr.groupbys, _attrs)...)
    gs = tuple([ d.g for d in SQ.dos(q) ]...)
    for (_attr,g) in zip(_attrs, gs)
        U = Core.Inference.return_type(g, Tuple{typeof(getfield(_res.src, _attr))})
        if isleaftype(U)
            push!(cols, Vector{U}())
        else # TODO: make this a real code path
            error()
        end
    end
    res = Relation(attrs, cols...)
    return res, fs, gs, _res
end

# gs are aggregating functions
@generated function grow_summaries!(res, gs, _res, level)
    m = length(level.parameters) # number of groupbys
    n = length(gs.parameters) # number of summarizations
    v_exs, push!_exs = Expr(:block), Expr(:block)
    for j in 1:m
        push!(v_exs.args,
              Expr(:(=), Symbol("v_$j"), Expr(:call, :getfield, :level, j)))
        push!(push!_exs.args,
              Expr(:call, :push!,
                   Expr(:call, :getfield, :(res.src), j), Symbol("v_$j")))
    end
    for j in 1:n
        push!(v_exs.args,
              Expr(:(=), Symbol("v_$(j+m)"),
              Expr(:call, Expr(:call, :getfield, :gs, j),
                   Expr(:call, :getfield, :(_res.src), j))))
        push!(push!_exs.args,
              Expr(:call, :push!,
                   Expr(:call, :getfield, :(res.src), j+m), Symbol("v_$(j+m)")))
    end
    return quote
        $v_exs; $push!_exs
    end
end

@generated function grow_nonnull!(_res, fs, src, indices)
    n = length(fs.parameters)
    v_exs, push!_exs = Expr(:block), Expr(:block)
    for j in 1:n
        push!(v_exs.args,
              Expr(:(=), Symbol("v_$j"),
              Expr(:call, Expr(:call, :getfield, :fs, j), :row)))
        # NOTE: This may result in a Relation with unequal-length columns!!
        #       Should really be a "Pseudo-Relation"
        push!(push!_exs.args,
              Expr(:if, Expr(:call, :!, Expr(:call, :isnull, Symbol("v_$j"))),
                   Expr(:call, :push!,
                        Expr(:call, :getfield, :(_res.src), j),
                             Expr(:call, :unsafe_get, Symbol("v_$j")))))
    end
    return quote
        for row in src[indices]
            $v_exs; $push!_exs
        end
        return
    end
end


text:
# These do most of the work in the package.
# They are all @generated or recusive functions for performance.

const UnionAllTupleOrVector = Union{Vector{UnionAll},Tuple{UnionAll,Vararg}}

"""
Sort dimensions into the order they take in the array.

Missing dimensions are replaced with `nothing`
"""
@inline Base.permutedims(tosort::AbDimTuple, perm::Union{Vector{<:Integer},Tuple{<:Integer,Vararg}}) =
    map(p -> tosort[p], Tuple(perm))
@inline Base.permutedims(tosort::Union{Vector{<:Integer},Tuple{<:Integer,Vararg}}, perm::AbDimTuple) =
    tosort

@inline Base.permutedims(tosort::AbDimTuple, order::UnionAllTupleOrVector) =
    permutedims(tosort, Tuple(map(u -> u(), order)))
@inline Base.permutedims(tosort::UnionAllTupleOrVector, order::AbDimTuple) =
    permutedims(Tuple(map(u -> u(), tosort)), order)
@inline Base.permutedims(tosort::AbDimTuple, order::AbDimVector) =
    permutedims(tosort, Tuple(order))
@inline Base.permutedims(tosort::AbDimVector, order::AbDimTuple) =
    permutedims(Tuple(tosort), order)
@inline Base.permutedims(tosort::AbDimTuple, order::AbDimTuple) =
    Base.permutedims(tosort, order, map(d -> dims(grid(d)), order))
@generated Base.permutedims(tosort::AbDimTuple, order::AbDimTuple, griddims) =
    permutedims_inner(tosort, order, griddims)

permutedims_inner(tosort::Type, order::Type, griddims::Type) = begin
    indexexps = []
    for (i, dim) in enumerate(order.parameters)
        dimindex = findfirst(d -> basetypeof(d) <: basetypeof(dim), tosort.parameters)
        if dimindex == nothing
            # The grid may allow dimensions not in dims
            # that will be transformed to the actual dimensions (ie for rotations).
            if griddims <: Nothing
                push!(indexexps, :(nothing))
            else
                gridindex = findfirst(d -> basetypeof(d) <: basetypeof(griddims.parameters[i]), tosort.parameters)
                if gridindex == nothing
                    push!(indexexps, :(nothing))
                else
                    push!(indexexps, :(tosort[$gridindex]))
                end
            end
        else
            push!(indexexps, :(tosort[$dimindex]))
        end
    end
    Expr(:tuple, indexexps...)
end

"""
Convert a tuple of AbstractDimension to indices, ranges or Colon.
"""
@inline dims2indices(A, lookup, emptyval=Colon()) =
    dims2indices(dims(A), lookup, emptyval)
@inline dims2indices(dims::Tuple, lookup, emptyval=Colon()) =
    dims2indices(dims, (lookup,), emptyval)
@inline dims2indices(dims::Tuple, lookup::Tuple, emptyval=Colon()) =
    dims2indices(map(grid, dims), dims, permutedims(lookup, dims), emptyval)

# Deal with irregular grid types that need multiple dimensions indexed together
@inline dims2indices(grids::Tuple{DependentGrid,Vararg}, dims::Tuple, lookup::Tuple, emptyval) = begin
    (irregdims, irreglookup), (regdims, reglookup) = splitgridtypes(grids, dims, lookup)
    (irreg2indices(map(grid, irregdims), irregdims, irreglookup, emptyval)...,
     dims2indices(map(grid, regdims), regdims, reglookup, emptyval)...)
end

@inline dims2indices(grids::Tuple, dims::Tuple, lookup::Tuple, emptyval) = begin
    (dims2indices(grids[1], dims[1], lookup[1], emptyval),
     dims2indices(tail(grids), tail(dims), tail(lookup), emptyval)...)
end
@inline dims2indices(grids::Tuple{}, dims::Tuple{}, lookup::Tuple{}, emptyval) = ()

@inline dims2indices(grid, dim::AbDim, lookup::Type{<:AbDim}, emptyval) = Colon()
@inline dims2indices(grid, dim::AbDim, lookup::Nothing, emptyval) = emptyval
@inline dims2indices(grid, dim::AbDim, lookup::AbDim, emptyval) = val(lookup)
@inline dims2indices(grid, dim::AbDim, lookup::AbDim{<:Selector}, emptyval) =
    sel2indices(grid, dim, val(lookup))

# Selectors select on grid dimensions
@inline irreg2indices(grids::Tuple, dims::Tuple, lookup::Tuple{AbDim{<:Selector},Vararg}, emptyval) =
    sel2indices(grids, dims, map(val, lookup))
# Other dims select on regular dimensions
@inline irreg2indices(grids::Tuple, dims::Tuple, lookup::Tuple, emptyval) = begin
    (dims2indices(grids[1], dims[1], lookup[1], emptyval),
     dims2indices(tail(grids), tail(dims), tail(lookup), emptyval)...)
end

@inline splitgridtypes(grids::Tuple{DependentGrid,Vararg}, dims, lookup) = begin
    (irregdims, irreglookup), reg = splitgridtypes(tail(grids), tail(dims), tail(lookup))
    irreg = (dims[1], irregdims...), (lookup[1], irreglookup...)
    irreg, reg
end
@inline splitgridtypes(grids::Tuple{IndependentGrid,Vararg}, dims, lookup) = begin
    irreg, (regdims, reglookup) = splitgridtypes(tail(grids), tail(dims), tail(lookup))
    reg = (dims[1], regdims...), (lookup[1], reglookup...)
    irreg, reg
end
@inline splitgridtypes(grids::Tuple{}, dims, lookup) = ((), ()), ((), ())

"""
Slice the dimensions to match the axis values of the new array

All methods returns a tuple conatining two tuples: the new dimensions,
and the reference dimensions. The ref dimensions are no longer used in
the new struct but are useful to give context to plots.

Called at the array level the returned tuple will also include the
previous reference dims attached to the array.
"""
@inline slicedims(A, dims::AbDimTuple) = slicedims(A, dims2indices(A, dims))
@inline slicedims(dims2slice::AbDimTuple, dims::AbDimTuple) =
    slicedims(dims2slice, dims2indices(dims2slice, dims))
@inline slicedims(A, I::Tuple) = slicedims(dims(A), refdims(A), I)
@inline slicedims(dims::Tuple, refdims::Tuple, I::Tuple) = begin
    newdims, newrefdims = slicedims(dims, I)
    newdims, (refdims..., newrefdims...)
end
@inline slicedims(dims::Tuple{}, I::Tuple) = ((), ())
@inline slicedims(dims::AbDimTuple, I::Tuple) = begin
    d = slicedims(dims[1], I[1])
    ds = slicedims(tail(dims), tail(I))
    (d[1]..., ds[1]...), (d[2]..., ds[2]...)
end
@inline slicedims(dims::Tuple{}, I::Tuple{}) = ((), ())
@inline slicedims(d::AbDim, i::Colon) = ((rebuild(d, val(d)),), ())
@inline slicedims(d::AbDim, i::Number) = ((), (rebuild(d, val(d)[i]),))
@inline slicedims(d::AbDim{<:AbstractArray}, i::AbstractArray) =
    ((rebuild(d, val(d)[i]),), ())
@inline slicedims(d::AbDim{<:LinRange}, i::UnitRange) = begin
    range = val(d)
    start, stop, len = range[first(i)], range[last(i)], length(i) ÷ step(i)
    ((rebuild(d, LinRange(start, stop, len)),), ())
end

"""
Get the number of an AbstractDimension as ordered in the array
"""
@inline dimnum(A, lookup) = dimnum(typeof(dims(A)), lookup)
@inline dimnum(dimtypes::Type, lookup::AbstractArray) = dimnum(dimtypes, (lookup...,))
@inline dimnum(dimtypes::Type, lookup::Number) = lookup
@inline dimnum(dimtypes::Type, lookup::Tuple) =
    (dimnum(dimtypes, lookup[1]), dimnum(dimtypes, tail(lookup))...,)
@inline dimnum(dimtypes::Type, lookup::Tuple{}) = ()
@inline dimnum(dimtypes::Type, lookup::AbDim) = dimnum(dimtypes, typeof(lookup))
@generated dimnum(dimtypes::Type{DTS}, lookup::Type{D}) where {DTS,D} = begin
    index = findfirst(dt -> D <: basetypeof(dt), DTS.parameters)
    if index == nothing
        :(throw(ArgumentError("No $lookup in $dimtypes")))
    else
        :($index)
    end
end

"""
Format the dimension to match internal standards.

Mostily this means converting tuples and UnitRanges to LinRange,
which is easier to handle. Errors are thrown if dims don't match the array dims or size.
"""
@inline formatdims(A::AbstractArray{T,N}, dims::Tuple) where {T,N} = begin
    dimlen = length(dims)
    dimlen == N || throw(ArgumentError("dims ($dimlen) don't match array dimensions $(N)"))
    formatdims(size(A), dims)
end
@inline formatdims(size::Tuple, dims::Tuple) = map(formatdims, size, dims)
@inline formatdims(len::Integer, dim::AbDim{<:AbstractArray}) =
    if length(val(dim)) == len
        dim
    else
        throw(ArgumentError("length of $dim $(length(val(dim))) does not match
                             size of array dimension $len"))
    end
@inline formatdims(len::Integer, dim::AbDim{<:Union{UnitRange,NTuple{2}}}) = linrange(dim, len)
@inline formatdims(len::Integer, dim::AbDim) = dim

linrange(dim, len) = begin
    range = val(dim)
    start, stop = first(range), last(range)
    rebuild(dim, LinRange(start, stop, len))
end

"""
Replace the specified dimensions with an index of 1 to match
a new array size where the dimension has been reduced to a length
of 1, but the number of dimensions has not changed.

Used in mean, reduce, etc.
"""
@inline reducedims(A, dimstoreduce) = reducedims(A, (dimstoreduce,))
@inline reducedims(A, dimstoreduce::Tuple) = reducedims(dims(A), dimstoreduce)
@inline reducedims(dims::AbDimTuple, dimstoreduce::Tuple) =
    map(reducedims, dims, permutedims(dimstoreduce, dims))
@inline reducedims(dims::AbDimTuple, dimstoreduce::Tuple{Vararg{Int}}) =
    map(reducedims, dims, permutedims(map(i -> dims[i], dimstoreduce), dims))

@inline reducedims(dim::AbDim, dimtoreduce::AbDim) = basetypeof(dim)(first(val(dim)))
@inline reducedims(dim::AbDim, dimtoreduce::Nothing) = dim

"""
Get the dimension(s) matching the type(s) of the lookup dimension.
"""
@inline dims(A, lookup) = dims(dims(A), lookup)
@inline dims(ds::AbDimTuple, lookup::Integer) = ds[lookup]
@inline dims(ds::AbDimTuple, lookup::Tuple) =
    (dims(ds, lookup[1]), dims(ds, tail(lookup))...)
@inline dims(ds::AbDimTuple, lookup::Tuple{}) = ()
@inline dims(ds::AbDimTuple, lookup) = dims(ds, basetypeof(lookup))
@generated dims(ds::DT, lookup::Type{L}) where {DT<:AbDimTuple,L} = begin
    index = findfirst(dt -> dt <: L, DT.parameters)
    if index == nothing
        :(throw(ArgumentError("No $lookup in $dims")))
    else
        :(ds[$index])
    end
end