text:
using FLOWMath
using Test
import ForwardDiff
import FiniteDiff
using LinearAlgebra: diag

@testset "FLOWMath.jl" begin

# ------ trapz --------
# tests from matlab trapz docs
x = 1:5
y = [1 4 9 16 25]
z = trapz(x, y)
@test z == 42.0

x = range(0.0, stop=pi+1e-15, step=pi/100)
y = sin.(x)
z = trapz(x, y)
@test isapprox(z, 1.999835503887444, atol=1e-15)

# -------------------------

# ------ Brent's method ------

f(x) = x^2 - 1.0

xstar, _ = brent(f, -2.0, 0)
@test xstar == -1.0
xstar, _ = brent(f, 0.0, 2)
@test xstar == 1.0

f(x) = x^3 - 1
xstar, _ = brent(f, 0, 3)
@test xstar == 1.0

f(x) = sin(x)
atol = 2e-12
xstar, _ = brent(f, 1, 4, atol=atol)
@test isapprox(xstar, pi, atol=atol)

atol = 1e-15
xstar, _ = brent(f, 1, 4, atol=atol)
@test isapprox(xstar, pi, atol=atol)

# -------------------------

# ------ abs_smooth ---------
x = -3.0
delta_x = 0.0
y = abs_smooth(x, delta_x)
@test y == 3.0

x = 3.0
delta_x = 0.1
y = abs_smooth(x, delta_x)
@test y == 3.0

x = -3.0
delta_x = 0.1
y = abs_smooth(x, delta_x)
@test y == 3.0

x = 0
delta_x = 0.1
y = abs_smooth(x, delta_x)
@test y == 0.05

x = 0.1
delta_x = 0.1
y = abs_smooth(x, delta_x)
@test y == 0.1

x = -0.1
delta_x = 0.1
y = abs_smooth(x, delta_x)
@test y == 0.1

x = 0.05
delta_x = 0.1
y = abs_smooth(x, delta_x)
@test y == 0.0625

x = -0.05
delta_x = 0.1
y = abs_smooth(x, delta_x)
@test y == 0.0625
# -------------------------

# ------ ksmax ---------
x = [0.0, 0.0]
x_max_smooth = ksmax(x)
@test isapprox(x_max_smooth, 0.013862943611198907)

# overflow
x = [1e6, 1e6]
x_max_smooth = ksmax(x)
@test isapprox(x_max_smooth, 1.0000000138629436e6)

# underflow
x = [-1e6, -1e6]
x_max_smooth = ksmax(x)
@test isapprox(x_max_smooth, -999999.9861370564)

# hardness
x = [0.0, 0.0]
hardness = 100.0
x_max_smooth = ksmax(x, hardness)
@test isapprox(x_max_smooth, 0.006931471805599453)

# ------ ksmin ---------
x = [0.0, 0.0]
x_max_smooth = ksmin(x)
@test isapprox(x_max_smooth, -0.013862943611198907)

# overflow
x = [1e6, 1e6]
x_max_smooth = ksmin(x)
@test isapprox(x_max_smooth, 999999.9861370564)

# underflow
x = [-1e6, -1e6]
x_max_smooth = ksmin(x)
@test isapprox(x_max_smooth, -1.0000000138629436e6)

# hardness
x = [0.0, 0.0]
hardness = 100.0
x_max_smooth = ksmin(x, hardness)
@test isapprox(x_max_smooth, -0.006931471805599453)

# -------------------------

# ------ ksmax_adaptive ---------
x = [0.0, 0.0]
x_max_smooth = ksmax_adaptive(x)
@test isapprox(x_max_smooth, 0.0008325546111576975)

# overflow
x = [1e6, 1e6]
x_max_smooth = ksmax_adaptive(x)
@test isapprox(x_max_smooth, 1.0000000008325547e6)

# underflow
x = [-1e6, -1e6]
x_max_smooth = ksmax_adaptive(x)
@test isapprox(x_max_smooth, -999999.9991674453)

# hardness
x = [0.0, 0.0]
hardness = 100.0
x_max_smooth = ksmax_adaptive(x, hardness)
@test isapprox(x_max_smooth, 0.0008325546111576975)

x = [-0.1, 0.0]
hardness = 100.0
x_max_smooth = ksmax_adaptive(x, hardness)
@test isapprox(x_max_smooth, 4.5398899216870535e-7)

# tolerance
x = [0.0, 0.0]
tol = 1e-3
x_max_smooth = ksmax_adaptive(x, tol=tol)
@test isapprox(x_max_smooth, 0.013862943611198907)

# smoothing_fraction
x = [-0.165, 0.0]
smoothing_fraction = 0.1
x_max_smooth = ksmax_adaptive(x, smoothing_fraction=smoothing_fraction)
@test isapprox(x_max_smooth, 5.261253727654178e-6)

x = [-0.165, 0.0]
smoothing_fraction = 0.2
x_max_smooth = ksmax_adaptive(x, smoothing_fraction=smoothing_fraction)
@test isapprox(x_max_smooth, 5.2856933329025475e-6)

# ------ ksmin_adaptive ---------
x = [0.0, 0.0]
x_max_smooth = ksmin_adaptive(x)
@test isapprox(x_max_smooth, -0.0008325546111576975)

# overflow
x = [1e6, 1e6]
x_max_smooth = ksmin_adaptive(x)
@test isapprox(x_max_smooth, 999999.9991674453)

# underflow
x = [-1e6, -1e6]
x_max_smooth = ksmin_adaptive(x)
@test isapprox(x_max_smooth, -1.0000000008325547e6)

# hardness
x = [0.0, 0.0]
hardness = 100.0
x_max_smooth = ksmin_adaptive(x, hardness)
@test isapprox(x_max_smooth, -0.0008325546111576975)

x = [0.1, 0.0]
hardness = 100.0
x_max_smooth = ksmin_adaptive(x, hardness)
@test isapprox(x_max_smooth, -4.5398899216870535e-7)

# tolerance
x = [0.0, 0.0]
tol = 1e-3
x_max_smooth = ksmin_adaptive(x, tol=tol)
@test isapprox(x_max_smooth, -0.013862943611198907)

# smoothing_fraction
x = [0.165, 0.0]
smoothing_fraction = 0.1
x_max_smooth = ksmin_adaptive(x, smoothing_fraction=smoothing_fraction)
@test isapprox(x_max_smooth, -5.261253727654178e-6)

x = [0.165, 0.0]
smoothing_fraction = 0.2
x_max_smooth = ksmin_adaptive(x, smoothing_fraction=smoothing_fraction)
@test isapprox(x_max_smooth, -5.2856933329025475e-6)

# -------------------------

# ------ sigmoid ---------
x = -0.5
y = sigmoid(x)
@test isapprox(y, 0.37754066879814546)

x = 0.0
y = sigmoid(x)
@test isapprox(y, 0.5)

x = 0.5
y = sigmoid(x)
@test isapprox(y, 0.6224593312018546)

# overflow
x = 1e6
y = sigmoid(x)
@test isapprox(y, 1.0)

# underflow
x = -1e6
y = sigmoid(x)
@test isapprox(y, 0.0)

# -------------------------

# ------ sigmoid_blend ---------
x = -0.1
f1x = x
f2x = x^2
xt = 0.0
y = sigmoid_blend(f1x, f2x, x, xt)
@test isapprox(y, -0.09926378639832867)

x = 0.0
f1x = x
f2x = x^2
xt = 0.0
y = sigmoid_blend(f1x, f2x, x, xt)
@test isapprox(y, 0.0)

x = 0.1
f1x = x
f2x = x^2
xt = 0.0
y = sigmoid_blend(f1x, f2x, x, xt)
@test isapprox(y, 0.010602356583185632)

# overflow
x = 1e6
f1x = x
f2x = x^2
xt = 0.0
y = sigmoid_blend(f1x, f2x, x, xt)
@test isapprox(y, f2x)

# underflow
x = -1e6
f1x = x
f2x = x^2
xt = 0.0
y = sigmoid_blend(f1x, f2x, x, xt)
@test isapprox(y, f1x)

# hardness
x = -0.1
f1x = x
f2x = x^2
xt = 0.0
hardness = 100
y = sigmoid_blend(f1x, f2x, x, xt, hardness)
@test isapprox(y, -0.09999500623444274)

x = 0.0
f1x = x
f2x = x^2
xt = 0.0
hardness = 100
y = sigmoid_blend(f1x, f2x, x, xt, hardness)
@test isapprox(y, 0.0)

x = 0.1
f1x = x
f2x = x^2
xt = 0.0
hardness = 100
y = sigmoid_blend(f1x, f2x, x, xt, hardness)
@test isapprox(y, 0.010004085808183225)

# vectorized
x = -0.25:0.05:0.25
f1x = x
f2x = x.^2
xt = 0.0
hardness = 100
y = sigmoid_blend.(f1x, f2x, x, xt, hardness)
ytest = [-0.24999999999566003,
         -0.19999999950532316,
         -0.14999994723186585,
         -0.09999500623444274,
         -0.04964862532647505,
          0.0,
          0.0028179104189035298,
          0.010004085808183225,
          0.022500039002533945,
          0.040000000329784596,
          0.062500000002604]
@test isapprox(y, ytest)

# -------------------------

# ------ cubic_blend ---------
x = -3.0
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = cubic_blend(f1x, f2x, x, xt, delta_x)
@test y == -3.0

x = 3.0
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = cubic_blend(f1x, f2x, x, xt, delta_x)
@test y == 9.0

x = 0.0
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = cubic_blend(f1x, f2x, x, xt, delta_x)
@test y == 0.0

x = -0.1
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = cubic_blend(f1x, f2x, x, xt, delta_x)
@test y == f1x

x = 0.1
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = cubic_blend(f1x, f2x, x, xt, delta_x)
@test y == f2x

x = -0.05
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = cubic_blend(f1x, f2x, x, xt, delta_x)
@test isapprox(y, -0.041796875000000004)

x = 0.05
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = cubic_blend(f1x, f2x, x, xt, delta_x)
@test isapprox(y, 0.009921875000000004)

# vectorized
x = -0.25:0.05:0.25
f1x = x
f2x = x.^2
xt = 0.0
delta_x = 0.1
y = cubic_blend.(f1x, f2x, x, xt, delta_x)
ytest = [-0.25,
         -0.2,
         -0.15,
         -0.1,
         -0.041796875000000004,
          0.0,
          0.009921875000000004,
          0.010000000000000002,
          0.0225,
          0.04000000000000001,
          0.0625]
@test isapprox(y, ytest)

# -------------------------
# ------ quintic_blend ---------
x = -3.0
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = quintic_blend(f1x, f2x, x, xt, delta_x)
@test y == -3.0

x = 3.0
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = quintic_blend(f1x, f2x, x, xt, delta_x)
@test y == 9.0

x = 0.0
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = quintic_blend(f1x, f2x, x, xt, delta_x)
@test y == 0.0

x = -0.1
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = quintic_blend(f1x, f2x, x, xt, delta_x)
@test y == f1x

x = 0.1
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = quintic_blend(f1x, f2x, x, xt, delta_x)
@test y == f2x

x = -0.05
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = quintic_blend(f1x, f2x, x, xt, delta_x)
@test isapprox(y, -0.044565429687500005)

x = 0.05
f1x = x
f2x = x^2
xt = 0.0
delta_x = 0.1
y = quintic_blend(f1x, f2x, x, xt, delta_x)
@test isapprox(y, 0.007416992187499999)

# vectorized
x = -0.25:0.05:0.25
f1x = x
f2x = x.^2
xt = 0.0
delta_x = 0.1
y = quintic_blend.(f1x, f2x, x, xt, delta_x)
ytest = [ -0.25,
          -0.2,
          -0.15,
          -0.1,
          -0.044565429687500005,
           0.0,
           0.007416992187499999,
           0.010000000000000002,
           0.0225,
           0.04000000000000001,
           0.0625]
@test isapprox(y, ytest)
# -------------------------

# ---------- akima interpolate -------------

x = 0:pi/4:2*pi
y = sin.(x)
xpt = 0:pi/16:2*pi

ypt = akima(x, y, xpt)

# FIXME: this is not an independent test.  right now it is just testing against refactoring.
ytest = [0.0, 0.21394356492851746, 0.40088834764831843, 0.5641656518405971, 0.7071067811865475, 0.8281808021479407, 0.920495128834866, 0.9793385864009951, 1.0, 0.9873810649285174, 0.9419417382415921, 0.8523082377305078, 0.7071067811865476, 0.5303300858899105, 0.35355339059327384, 0.17677669529663714, 1.2246467991473532e-16, -0.17677669529663698, -0.35355339059327373, -0.5303300858899104, -0.7071067811865475, -0.8523082377305078, -0.941941738241592, -0.9873810649285175, -1.0, -0.9793385864009951, -0.9204951288348657, -0.8281808021479408, -0.7071067811865477, -0.5641656518405971, -0.400888347648319, -0.21394356492851788, 0.0]

@test isapprox(ypt, ytest, atol=1e-8)

spl = Akima(x, y)

dydx = derivative(spl, pi/16)

@test isapprox(dydx, 1.0180260961104746, atol=1e-12)

dydx = gradient(spl, xpt)

dydx_test = [1.1640128599466306, 1.0180260961104746, 0.889005522603991, 0.7769511394271804, 0.6818629465800423, 0.5473879346340718, 0.3889191062220746, 0.20645646134405082, 1.8419755006770065e-16, -0.13818978335121798, -0.33430576382780836, -0.5883479414297708, -0.900316316157106, -0.9003163161571062, -0.9003163161571062, -0.9003163161571061, -0.9003163161571061, -0.9003163161571064, -0.9003163161571064, -0.9003163161571065, -0.9003163161571061, -0.5883479414297709, -0.33430576382780913, -0.13818978335121862, -1.8419755006770062e-16, 0.20645646134405057, 0.3889191062220747, 0.5473879346340714, 0.6818629465800421, 0.7769511394271802, 0.8890055226039906, 1.0180260961104746, 1.1640128599466313]

@test isapprox(dydx, dydx_test, atol=1e-12)

# differentiate with dual numbers

dydx = ForwardDiff.derivative(spl, pi/16)
@test isapprox(dydx, 1.0180260961104746, atol=1e-12)

J = ForwardDiff.jacobian(spl, xpt)
dydx = diag(J)
@test isapprox(dydx, dydx_test, atol=1e-12)

# differentiate coordinates of spline

wrapper(y) = Akima(x, y, 0.1)(xpt)
J = ForwardDiff.jacobian(wrapper, y)
J2 = FiniteDiff.finite_difference_jacobian(wrapper, y)
@test maximum(abs.(J - J2)) < 1e-6

wrapper2(y) = Akima(x, y, 0.1)(xpt)
J = ForwardDiff.jacobian(wrapper2, x)
J2 = FiniteDiff.finite_difference_jacobian(wrapper2, x)
@test maximum(abs.(J - J2)) < 1e-6

# ---------------------------

# ----- 1D linear interpolation ------
xvec = [1.0, 2.0, 4.0, 5.0]
yvec = [2.0, 3.0, 5.0, 8.0]

y = linear(xvec, yvec, 1.0)
@test y == 2.0
y = linear(xvec, yvec, 1.5)
@test y == 2.5
y = linear(xvec, yvec, 3.0)
@test y == 4.0
y = linear(xvec, yvec, 4.5)
@test y == 6.5
y = linear(xvec, yvec, 5.0)
@test y == 8.0

y = linear(xvec, yvec, [1.0, 1.5, 3.0, 4.5, 5.0])
@test y == [2.0, 2.5, 4.0, 6.5, 8.0]

dydx = derivative(xvec, yvec, 1.0)
@test dydx == 1.0
dydx = derivative(xvec, yvec, 1.5)
@test dydx == 1.0
dydx = derivative(xvec, yvec, 3.0)
@test dydx == 1.0
dydx = derivative(xvec, yvec, 4.5)
@test dydx == 3.0
dydx = derivative(xvec, yvec, 5.0)
@test dydx == 3.0

dydx = gradient(xvec, yvec, [1.0, 1.5, 3.0, 4.5, 5.0])
@test dydx == [1.0, 1.0, 1.0, 3.0, 3.0]

# ---------------------------

# ----- 2D linear interpolation -----
xdata = [1.0, 2.0, 3.0, 4.0]
ydata = [2.0, 4.0, 6.0, 8.0]
fdata = [1.0 2.0 3.0 4.0;
        5.0 6.0 7.0 8.0;
        9.0 10.0 11.0 12.0;
        13.0 14.0 15.0 16.0]
interp1d = linear

f = interp2d(interp1d, xdata, ydata, fdata, 1.0, 2.0)
@test f[1,1] == 1.0
f = interp2d(interp1d, xdata, ydata, fdata, 1.5, 2.0)
@test f[1,1] == 3.0
f = interp2d(interp1d, xdata, ydata, fdata, 3.5, 5.0)
@test f[1,1] == 12.5

# with ForwardDiff
wrapper(x) = interp2d(interp1d, xdata, ydata, fdata, x, 5.0)
J = ForwardDiff.derivative(wrapper, 2.5)
J = ForwardDiff.jacobian(wrapper, [2.5, 3.0, 4.0])

wrapper2(y) = interp2d(interp1d, xdata, ydata, fdata, 2.5, y)
J = ForwardDiff.derivative(wrapper2, 3.0)
J = ForwardDiff.jacobian(wrapper2, [3.0, 5.0, 7.0])

# ---------------------------

# ----- 3D linear interpolation -----
xdata = [1.0, 2.0, 3.0, 4.0]
ydata = [2.0, 4.0, 6.0, 8.0]
zdata = [3.0, 5.0, 7.0, 9.0]

fdata = zeros(4, 4, 4)
fdata[:,:,1] = [1.0 2.0 3.0 4.0;
                5.0 6.0 7.0 8.0;
                9.0 10.0 11.0 12.0;
                13.0 14.0 15.0 16.0]
fdata[:,:,2] = fdata[:,:,1] .+ 1
fdata[:,:,3] = fdata[:,:,1] .+ 2
fdata[:,:,4] = fdata[:,:,1] .+ 3
interp1d = linear

f = interp3d(interp1d, xdata, ydata, zdata, fdata, 1.0, 2.0, 3.0)
@test f[1,1,1] == 1.0
f = interp3d(interp1d, xdata, ydata, zdata, fdata, 1.5, 2.0, 3.0)
@test f[1,1,1] == 3.0
f = interp3d(interp1d, xdata, ydata, zdata, fdata, 1.0, 3.0, 3.0)
@test f[1,1,1] == 1.5
f = interp3d(interp1d, xdata, ydata, zdata, fdata, 1.0, 2.0, 6.0)
@test f[1,1,1] == 2.5
f = interp3d(interp1d, xdata, ydata, zdata, fdata, 3.5, 5.0, 4.0)
@test f[1,1,1] == 13.0

# with ForwardDiff
wrapper(x) = interp3d(interp1d, xdata, ydata, zdata, fdata, x, 2.0, 3.0)
J = ForwardDiff.derivative(wrapper, 1.5)
J = ForwardDiff.jacobian(wrapper, [1.5; 2.5; 3.5])

wrapper2(y) = interp3d(interp1d, xdata, ydata, zdata, fdata, 1.0, y, 3.0)
J = ForwardDiff.derivative(wrapper2, 5.0)
J = ForwardDiff.jacobian(wrapper2, [3.0; 5.0; 7.0])

wrapper3(z) = interp3d(interp1d, xdata, ydata, zdata, fdata, 1.0, 2.0, z)
J = ForwardDiff.derivative(wrapper3, 8.0)
J = ForwardDiff.jacobian(wrapper3, [4.0; 6.0; 8.0])
# ---------------------------

# ----- 4D linear interpolation -----
xdata = [1.0, 2.0, 3.0, 4.0]
ydata = [2.0, 4.0, 6.0, 8.0]
zdata = [3.0, 5.0, 7.0, 9.0]
tdata = [0.0, 1.0, 1.5, 2.0]

fdata = zeros(4, 4, 4, 4)
fdata[:,:,1,1] = [1.0 2.0 3.0 4.0;
                  5.0 6.0 7.0 8.0;
                  9.0 10.0 11.0 12.0;
                  13.0 14.0 15.0 16.0]
for i in 1:4
    for j = 1:4
        fdata[:,:,i,j] = fdata[:,:,1,1] .+ (i-1) .+ (j-1)
    end
end
interp1d = linear

f = interp4d(interp1d, xdata, ydata, zdata, tdata, fdata, 1.0, 2.0, 3.0, 0.5)
@test f[1,1,1,1] == 1.5
f = interp4d(interp1d, xdata, ydata, zdata, tdata, fdata, 1.5, 2.0, 3.0, 1.0)
@test f[1,1,1,1] == 4.0
f = interp4d(interp1d, xdata, ydata, zdata, tdata, fdata, 1.0, 3.0, 3.0, 0.0)
@test f[1,1,1,1] == 1.5
f = interp4d(interp1d, xdata, ydata, zdata, tdata, fdata, 1.0, 2.0, 6.0, 1.5)
@test f[1,1,1,1] == 4.5
f = interp4d(interp1d, xdata, ydata, zdata, tdata, fdata, 3.5, 5.0, 4.0, 1.75)
@test f[1,1,1,1] == 15.5

# with ForwardDiff
wrapper(x) = interp4d(interp1d, xdata, ydata, zdata, tdata, fdata, x, 2.0, 3.0, 1.0)
J = ForwardDiff.derivative(wrapper, 1.5)
J = ForwardDiff.jacobian(wrapper, [1.5; 2.5; 3.5])

wrapper2(y) = interp4d(interp1d, xdata, ydata, zdata, tdata, fdata, 1.0, y, 3.0, 1.0)
J = ForwardDiff.derivative(wrapper2, 5.0)
J = ForwardDiff.jacobian(wrapper2, [3.0; 5.0; 7.0])

wrapper3(z) = interp4d(interp1d, xdata, ydata, zdata, tdata, fdata, 1.0, 2.0, z, 1.0)
J = ForwardDiff.derivative(wrapper3, 8.0)
J = ForwardDiff.jacobian(wrapper3, [4.0; 6.0; 8.0])

wrapper4(t) = interp4d(interp1d, xdata, ydata, zdata, tdata, fdata, 1.0, 2.0, 3.0, t)
J = ForwardDiff.derivative(wrapper4, 8.0)
J = ForwardDiff.jacobian(wrapper4, [3.0; 2.0; 1.0])
# ---------------------------

end


text:
struct Point1D
    x :: Float64
end

struct Point2D
    x :: Float64
    y :: Float64
end

*(a, p::Point2D) = Point2D(a*p.x,a*p.y)
+(p::Point2D, q::Point2D) = Point2D(p.x+q.x, p.y+q.y)

struct Bond1D{T<:Number}
    one :: Int
    two :: Int
    offx :: Int
    t   :: T
end

struct Bond2D{T<:Number}
    one :: Int
    two :: Int
    offx :: Int
    offy :: Int
    t   :: T
end

struct UnitCell1D{T<:Number}
    n :: Int
    sites :: Vector{Point1D}
    vs :: Vector{Point1D}
    bonds :: Vector{Bond1D{T}}
    mus :: Vector{T}
end

struct UnitCell2D{T<:Number}
    n :: Int
    sites :: Vector{Point2D}
    vs :: Vector{Point2D}
    bonds :: Vector{Bond2D{T}}
    mus :: Vector{T}
end

function enlargeunitcell(uc::UnitCell2D{T}, lx::Int, ly::Int) where{T}

    sites = Point2D[]
    bonds = Bond2D{T}[]
    mus = T[]
    for (y, x) in Iterators.product(1:ly, 1:lx)
        for idx in 1:uc.n
            push!(sites, (y-1)*uc.vs[2] + (x-1)*uc.vs[1] + uc.sites[idx])
            push!(mus, uc.mus[idx])
        end
        l = (x-1) * ly * uc.n + (y-1) * uc.n
        for b in uc.bonds
            one = l+b.one

            offx, padx = fldmod(x+b.offx-1, lx)
            offy, pady = fldmod(y+b.offy-1, ly)

            two = b.two + (padx * uc.n * ly) + (pady * uc.n)
            push!(bonds, Bond2D(one, two, offx, offy, b.t))
        end
    end
    UnitCell2D(
        uc.n*ly*lx,
        sites,
        [lx * uc.vs[1], ly * uc.vs[2]],
        bonds,
        mus,
    )
end

function chainunitcell(t::T, mu::T=zero(T)) where {T<:Number}
    uc = UnitCell1D(
        1,
        [Point1D(0,0)],
        [Point1D(+1)],
        [Bond2D(1, 1, +1, t)],
        [mu])
    uc
end

function squareunitcell(tx::T, ty::T, mu::T=zero(T)) where {T<:Number}
    uc = UnitCell2D(
        1,
        [Point2D(0,0)],
        [Point2D(0,1), Point2D(1,0)],
        [Bond2D(1, 1, +1, 0, tx), Bond2D(1, 1, 0, +1, ty)],
        [mu])
    uc
end

function triangularunitcell(t1::T, t2::T, t3::T, mu::T=zero(T)) where {T<:Number}
    uc = UnitCell2D(
        1,
        [Point2D(0,0)],
        [Point2D(1,0), Point2D(0.5,sin(pi/3))],
        [Bond2D(1, 1, +1, 0, t1), Bond2D(1, 1, 0, +1, t2), Bond2D(1, 1, -1, +1, t3)],
        [mu])
    uc
end

abstract type AbstractHamiltonian end
mutable struct BdGHamiltonian{T} <: AbstractHamiltonian
    nsites :: Int
    mat :: Matrix{T}
end
function BdGHamiltonian{T}(nsites::Int) where {T<:Number}
    BdGHamiltonian{T}(nsites, zeros(nsites, nsites))
end

function initializehamilt(::Type{T}, mode::Symbol, nsites::Int) where {T<:Number}
    if mode == :BdG
        BdGHamiltonian{T}(nsites)
    elseif mode == :SPIN
        SpinHamiltonian(nsites)
    elseif mode == :SPINU1
        U1SpinHamiltonian(nsites, div(nsites, 2))
    else
        error("unrecognized mode : $mode")
    end
end

function addterm!(H::BdGHamiltonian{T}, i::Int, mu::T) where {T<:Number}
    H.mat[i, i] += mu
    H
end

function addterm!(H::BdGHamiltonian{T}, i::Int, j::Int, t::T) where {T<:Number}
    H.mat[i, j] += t
    H.mat[j, i] += conj(t)
    H
end

mutable struct SpinHamiltonian <: AbstractHamiltonian
    nsites :: Int
    mat :: SparseMatrixCSC{Float64, Int}
end
function SpinHamiltonian(nsites::Int)
    SpinHamiltonian(nsites, spzeros(2^nsites, 2^nsites))
end

function addterm!(H::SpinHamiltonian, i::Int, mu::Float64)
    H.mat += mu * opextend(sz_half+1/2*I2, H.nsites, i)
    H
end

function addterm!(H::SpinHamiltonian, i::Int, j::Int, J::Float64)
    i == j && error("Add term for the same site!")
    idx1st = j > i ? i : j
    dis = abs(i - j)
    if dis > 1
        between = foldl(⊗, [2*sz_half for i=1:dis-1])
        H.mat += J * opextend(2*sp_half*sz_half ⊗ between ⊗ sm_half, H.nsites-dis, idx1st)
        H.mat += J * opextend(2*sz_half*sm_half ⊗ between ⊗ sp_half, H.nsites-dis, idx1st)
    else
        H.mat += J * opextend(2*sp_half*sz_half ⊗ sm_half, H.nsites-dis, idx1st)
        H.mat += J * opextend(2*sz_half*sm_half ⊗ sp_half, H.nsites-dis, idx1st)
    end
    H
end

mutable struct U1SpinHamiltonian <: AbstractHamiltonian
    nsites  :: Int
    mat     :: SparseMatrixCSC{Float64}
    I       :: Vector{Int32}
    J       :: Vector{Int32}
    V       :: Vector{Float64}
    sector  :: Int
    states  :: Vector{Int}
end

function U1SpinHamiltonian(nsites, sector)
    blocksize = binomial(nsites, sector)
    states = Vector{Int}(undef, blocksize)
    index = 0
    for state=0:2^nsites-1
        if count_ones(state) == sector
            index += 1
            states[index] = state
        end
    end
    U1SpinHamiltonian(nsites,
                      spzeros(blocksize, blocksize),
                      [], [], [],
                      sector, states)
end

function addterm!(H::U1SpinHamiltonian, i::Int, mu::Float64)
    for idx in eachindex(H.states)
        state = H.states[idx]
        if (state >> i) & 1 != 0
            push!(H.I, idx)
            push!(H.J, idx)
            push!(H.V, mu)
        end
    end
    H
end

function addterm!(H::U1SpinHamiltonian, i::Int, j::Int, J::Float64)
    i == j && error("Add term for the same site!")
    idx1st = j > i ? i : j
    dis = abs(i - j)
    for idx in eachindex(H.states)
        state = H.states[idx]
        if xor((state >> (i-1)) & 1, (state >> (j-1)) & 1) != 0
            if dis > 1
                sign = iseven(dis-1 - count_ones(state & (sum([2^n for n=0:dis-2])*2^idx1st))) ? +1 : -1
            else
                sign = +1
            end
            flippedstate = xor(state, 1 << (i-1) | 1 << (j-1))
            idx_flipped = searchsortedfirst(H.states, flippedstate)
            push!(H.I, idx)
            push!(H.J, idx_flipped)
            push!(H.V, -sign*J)
        end
    end
    H
end

function makehamiltonian(
    uc::UnitCell2D{T},
    lx::Int,
    ly::Int;
    boundary::Symbol=:OBC,
    mode::Symbol=:BdG) where {T<:Number}

    nsites = uc.n*ly*lx
    H = initializehamilt(T, mode, nsites)
    for (y,x) in Iterators.product(1:ly, 1:lx)
        l = (x-1) * ly * uc.n + (y-1) * uc.n
        for n in 1:uc.n
            addterm!(H, l+n, uc.mus[n])
        end
        for b in uc.bonds
            index1 = l+b.one
            if 0 < y+b.offy <= ly && 0 < x+b.offx <= lx
                index2 = l+b.two + (b.offx * uc.n * ly) + (b.offy * uc.n)
                addterm!(H, index1, index2, -b.t)
            elseif boundary == :PBCX && 0 < y+b.offy <= ly
                x+b.offx > lx ? lnew = l - lx * ly * uc.n : lnew = l + lx * ly * uc.n
                index2 = lnew + b.two + (b.offx * uc.n * ly) + (b.offy * uc.n)
                addterm!(H, index1, index2, -b.t)
            elseif boundary == :APBCX && 0 < y+b.offy <= ly
                x+b.offx > lx ? lnew = l - lx * ly * uc.n : lnew = l + lx * ly * uc.n
                index2 = lnew + b.two + (b.offx * uc.n * ly) + (b.offy * uc.n)
                addterm!(H, index1, index2, b.t)
            elseif !(boundary in [:OBC, :PBCX, :APBCX])
                error("unrecognized boundary condition : ", boundary)
            end
        end
    end
    H
end

##TODO: combine this stuff into a function that works for every dimension.
function makemodel(
    uc::UnitCell2D{T},
    lx::Int,
    ly::Int;
    boundary::Symbol=:OBC) where{T<:Number}

    nsites = uc.n*ly*lx
    H = zeros(T, nsites, nsites)
    for (y,x) in Iterators.product(1:ly,1:lx)
        l = (x-1) * ly * uc.n + (y-1) * uc.n
        for n in 1:uc.n
            H[l+n,l+n] = uc.mus[n]
        end
        for b in uc.bonds
            index1 = l+b.one
            if 0 < y+b.offy <= ly && 0 < x+b.offx <= lx
                index2 = l+b.two + (b.offx * uc.n * ly) + (b.offy * uc.n)
                #println("$x, $y => ($index1, $index2)")
                H[index1, index2] += -b.t
                H[index2, index1] += -conj(b.t)
            elseif boundary == :PBCX && 0 < y+b.offy <= ly
                x+b.offx > lx ? lnew = l - lx * ly * uc.n : lnew = l + lx * ly * uc.n
                index2 = lnew + b.two + (b.offx * uc.n * ly) + (b.offy * uc.n)
                #println("wrap x, $y => ($index1, $index2)")
                H[index1, index2] += -b.t
                H[index2, index1] += -conj(b.t)
            elseif !(boundary in [:OBC, :PBCX])
                error("unrecognized boundary condition : ", boundary)
            end
        end
    end
    H
end

function makemodelJW(
    uc::UnitCell2D{T},
    lx::Int,
    ly::Int;
    boundary::Symbol=:OBC) where{T<:Number}

    @assert uc.n == 1
    nsites = uc.n*ly*lx
    H = spzeros(T, 2^nsites, 2^nsites)
    for (y,x) in Iterators.product(1:ly,1:lx)
        l = (x-1) * ly * uc.n + (y-1) * uc.n
        for n in 1:uc.n
            H += uc.mus[n] .* opextend(sz_half+1/2*I2, lx*ly*uc.n, l+n)
        end
        for b in uc.bonds
            index1 = l+b.one
            if 0 < y+b.offy <= ly && 0 < x+b.offx <= lx
                index2 = l+b.two + (b.offx * uc.n * ly) + (b.offy * uc.n)
                println("$x, $y => ($index1, $index2)")
                if index2 > index1
                    idx1st = index1
                    dis = index2 - index1
                else
                    dis = index1 - index2
                    idx1st = index2
                end
                if dis > 1
                    between = foldl(⊗, [2*sz_half for i=1:dis-1])
                    H += -b.t * opextend(2*sp_half*sz_half ⊗ between ⊗ sm_half, nsites-dis, idx1st)
                    H += -conj(b.t) * opextend(2*sz_half*sm_half ⊗ between ⊗ sp_half, nsites-dis, idx1st)
                else
                    H += -b.t * opextend(2*sp_half*sz_half ⊗ sm_half, nsites-dis, idx1st)
                    H += -conj(b.t) * opextend(2*sz_half*sm_half ⊗ sp_half, nsites-dis, idx1st)
                end
            elseif boundary == :PBCX && 0 < y+b.offy <= ly
                x+b.offx > lx ? lnew = l - lx * ly * uc.n : lnew = l + lx * ly * uc.n
                index2 = lnew + b.two + (b.offx * uc.n * ly) + (b.offy * uc.n)
                println("wrap x, $y => ($index1, $index2)")
                if index2 > index1
                    idx1st = index1
                    dis = index2 - index1
                else
                    dis = index1 - index2
                    idx1st = index2
                end
                if dis > 1
                    between = foldl(⊗, [2*sz_half for i=1:dis-1])
                    H += -b.t * opextend(2*sp_half*sz_half ⊗ between ⊗ sm_half, nsites-dis, idx1st)
                    H += -conj(b.t) * opextend(2*sz_half*sm_half ⊗ between ⊗ sp_half, nsites-dis, idx1st)
                else
                    H += -b.t * opextend(2*sp_half*sz_half ⊗ sm_half, nsites-dis, idx1st)
                    H += -conj(b.t) * opextend(2*sz_half*sm_half ⊗ sp_half, nsites-dis, idx1st)
                end
            elseif !(boundary in [:OBC, :PBCX])
                error("unrecognized boundary condition : ", boundary)
            end
        end
    end
    H
end

function makemodelJW_u1sym(
    uc::UnitCell2D{T},
    lx::Int,
    ly::Int;
    boundary::Symbol=:OBC,
    zsector::Int=((lx*ly*uc.n) % 2)) where{T<:Number}

    @assert uc.n == 1
    nsites = uc.n*ly*lx
    #H = spzeros(T, 2^nsites, 2^nsites)

    M = div(nsites+zsector, 2)
    block_size = binomial(nsites, M)

    #println("$M, $block_size")
    I = Int32[]
    J = Int32[]
    V = T[]

    szblock_states = Vector{Int}(undef, block_size)
    index = 0
    for state=0:2^(nsites)-1
        if count_ones(state) == M
            #println("$state, $(string(state, base=2))")
            index += 1
            szblock_states[index] = state
        end
    end

    # for all states
    for state_index=1:block_size
        state = szblock_states[state_index]
        # for all terms in Hamiltonian
        for (y,x) in Iterators.product(1:ly,1:lx)
            l = (x-1) * ly * uc.n + (y-1) * uc.n
            for n in 1:uc.n
                if (state >> (l+n)) & 1 != 0
                    push!(I, state_index)
                    push!(J, state_index)
                    push!(V, uc.mus[n])
                end
            end
            for b in uc.bonds
                index1 = l+b.one
                if 0 < y+b.offy <= ly && 0 < x+b.offx <= lx
                    index2 = l+b.two + (b.offx * uc.n * ly) + (b.offy * uc.n)
                    #println("$x, $y => ($index1, $index2)")
                    if index2 > index1
                        idx1st = index1
                        dis = index2 - index1
                    else
                        dis = index1 - index2
                        idx1st = index2
                    end
                    if xor((state >> (index1-1)) & 1, (state >> (index2-1)) & 1) != 0
                        if dis > 1
                            #println("dis = $dis, state = $(string(state, base=2)), stuff = $(string(sum([2^i for i=1:dis-1])*2^idx1st, base=2)), and = $(state & (sum([2^i for i=1:dis-1])*2^idx1st))")
                            sign = iseven(dis-1 - count_ones(state & (sum([2^i for i=0:dis-2])*2^idx1st))) ? +1 : -1
                        else
                            sign = +1
                        end

                        flipped_state = xor(state, 1 << (index1-1) | 1 << (index2-1))
                        flipped_state_index = searchsortedfirst(szblock_states, flipped_state)

                        #println("$index1, $index2, $(string(state, base=2)), $(string(flipped_state, base=2)), indeces : $state_index, $flipped_state_index, sign = $sign")
                        if flipped_state_index > block_size
                            error("Fucked, $(string(state, base=2)), $(string(flipped_state, base=2))")
                        end
                        push!(I, state_index)
                        push!(J, flipped_state_index)
                        push!(V, sign*b.t)
                    end

                elseif boundary == :PBCX && 0 < y+b.offy <= ly
                    x+b.offx > lx ? lnew = l - lx * ly * uc.n : lnew = l + lx * ly * uc.n
                    index2 = lnew + b.two + (b.offx * uc.n * ly) + (b.offy * uc.n)
                    #println("wrap x, $y => ($index1, $index2)")
                    if index2 > index1
                        idx1st = index1
                        dis = index2 - index1
                    else
                        dis = index1 - index2
                        idx1st = index2
                    end
                    if xor((state >> (index1-1)) & 1, (state >> (index2-1)) & 1) != 0
                        if dis > 1
                            sign = iseven(dis-1-count_ones(state & (sum([2^i for i=0:dis-2])*2^idx1st))) ? +1 : -1
                        else
                            sign = +1
                        end
                        flipped_state = xor(state, 1 << (index1-1) | 1 << (index2-1))
                        flipped_state_index = searchsortedfirst(szblock_states, flipped_state)
                        push!(I, state_index)
                        push!(J, flipped_state_index)
                        push!(V, sign*b.t)
                    end
                elseif !(boundary in [:OBC, :PBCX])
                    error("unrecognized boundary condition : ", boundary)
                end
            end
        end
    end
    #println(I)
    #println(J)
    return sparse(I, J, V, block_size, block_size, +), (szblock_states .+ 1)
end

ext: jl

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
module MaxHelpingHandMoreCustomNPC

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/MoreCustomNPC" MoreCustomNPC(x::Integer, y::Integer, sprite::String="player/idle", spriteRate::Int=1, dialogId::String="", onlyOnce::Bool=true, endLevel::Bool=false,
	flipX::Bool=false, flipY::Bool=false, approachWhenTalking::Bool=false, approachDistance::Int=16, indicatorOffsetX::Int=0, indicatorOffsetY::Int=0,
	frames::String="", onlyIfFlag::String="", setFlag::String="", nodes::Array{Tuple{Integer, Integer}, 1}=Tuple{Integer, Integer}[])

@mapdef Entity "MaxHelpingHand/CustomNPCSprite" CustomNPCSprite(x::Integer, y::Integer, dialogId::String="", onlyOnce::Bool=true, endLevel::Bool=false,
	flipX::Bool=false, flipY::Bool=false, approachWhenTalking::Bool=false, approachDistance::Int=16, indicatorOffsetX::Int=0, indicatorOffsetY::Int=-16,
	spriteName::String="bird", onlyIfFlag::String="", setFlag::String="", nodes::Array{Tuple{Integer, Integer}, 1}=Tuple{Integer, Integer}[])

const npcUnion = Union{MoreCustomNPC, CustomNPCSprite}

const placements = Ahorn.PlacementDict(
    # Base placements
    "More Custom NPC\n(Everest + max480's Helping Hand)" => Ahorn.EntityPlacement(
        MoreCustomNPC
    ),
    "Custom NPC (from XML)\n(Everest + max480's Helping Hand)" => Ahorn.EntityPlacement(
        CustomNPCSprite
    ),

    # Presets for Granny and an invisible NPC that don't actually require More Custom NPC but are here for convenience/completeness
    "Custom NPC (Granny) (Everest)" => Ahorn.EntityPlacement(
        Maple.EverestCustomNPC,
        "point",
        Dict{String,Any}(
            "sprite" => "oldlady/idle",
            "spriteRate" => 7,
            "indicatorOffsetY" => -2
        )
    ),
    "Custom NPC (Invisible) (Everest)" => Ahorn.EntityPlacement(
        Maple.EverestCustomNPC,
        "point",
        Dict{String,Any}(
            "sprite" => "",
            "indicatorOffsetY" => -16
        )
    ),

    # Presets for Theo, Oshiro and Badeline Boss that do require More Custom NPC
    "More Custom NPC (Theo)\n(Everest + max480's Helping Hand)" => Ahorn.EntityPlacement(
        MoreCustomNPC,
        "point",
        Dict{String,Any}(
            "sprite" => "theo/theo",
            "frames" => "0-9",
            "spriteRate" => 10,
            "indicatorOffsetY" => -5
        )
    ),
    "More Custom NPC (Oshiro)\n(Everest + max480's Helping Hand)" => Ahorn.EntityPlacement(
        MoreCustomNPC,
        "point",
        Dict{String,Any}(
            "sprite" => "oshiro/oshiro",
            "frames" => "31-41",
            "spriteRate" => 12,
            "indicatorOffsetY" => -10
        )
    ),
    "More Custom NPC (Badeline Boss)\n(Everest + max480's Helping Hand)" => Ahorn.EntityPlacement(
        MoreCustomNPC,
        "point",
        Dict{String,Any}(
            "sprite" => "badelineBoss/boss",
            "frames" => "0-23",
            "spriteRate" => 17,
            "indicatorOffsetY" => -5
        )
    )
)

Ahorn.nodeLimits(entity::npcUnion) = 0, 2

function getSpriteName(entity::MoreCustomNPC)
    spriteName = get(entity.data, "sprite", "")

    if spriteName == "oshiro/oshiro"
        return "characters/oshiro/oshiro31"
    end

    return "characters/$(spriteName)00"
end

function getSpriteName(entity::CustomNPCSprite)
    return "ahorn/MaxHelpingHand/custom_npc_xml"
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::npcUnion)
    x, y = Ahorn.position(entity)

    spriteName = getSpriteName(entity)
    sprite = Ahorn.getSprite(spriteName, "Gameplay")

    nodes = get(entity.data, "nodes", Tuple{Int, Int}[])
    offsetCenterY = - floor(Int, sprite.height / 2)

    if length(nodes) == 2
        n1x, n1y = Int.(nodes[1])
        n2x, n2y = Int.(nodes[2])
        Ahorn.drawRectangle(ctx, n1x, n1y, n2x - n1x + 8, n2y - n1y + 8, Ahorn.colors.trigger_fc, Ahorn.colors.trigger_bc)
    end

    for node in nodes
        nx, ny = Int.(node)
        Ahorn.drawRectangle(ctx, nx, ny, 8, 8, Ahorn.colors.trigger_fc, Ahorn.colors.trigger_bc)
    end
end

function Ahorn.selection(entity::npcUnion)
    x, y = Ahorn.position(entity)

    scaleX, scaleY = get(entity.data, "flipX", false) ? -1 : 1, get(entity.data, "flipY", false) ? -1 : 1

    spriteName = getSpriteName(entity)
    sprite = Ahorn.getSprite(spriteName, "Gameplay")

    if sprite.width == 0 || sprite.height == 0
        res = Ahorn.Rectangle[Ahorn.Rectangle(x - 4, y - 4, 8, 8)]
    else
        res = Ahorn.Rectangle[Ahorn.getSpriteRectangle(spriteName, x, y, jx=0.5, jy=1.0, sx=scaleX, sy=scaleY)]
    end

    for node in get(entity.data, "nodes", ())
        nx, ny = Int.(node)
        push!(res, Ahorn.Rectangle(nx, ny, 8, 8))
    end

    return res
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::npcUnion, room::Maple.Room)
    scaleX, scaleY = get(entity.data, "flipX", false) ? -1 : 1, get(entity.data, "flipY", false) ? -1 : 1
    spriteName = getSpriteName(entity)

    Ahorn.drawSprite(ctx, spriteName, 0, 0, jx=0.5, jy=1.0, sx=scaleX, sy=scaleY)
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
module LibLMDB

import ..checked_call

using LMDB_jll
export LMDB_jll

using CEnum

const __mode_t = Cushort

const mode_t = __mode_t

mutable struct MDB_txn end

const MDB_dbi = Cuint

function _mdb_dbi_open(txn, name, flags, dbi)
    ccall((:mdb_dbi_open, liblmdb), Cint, (Ptr{MDB_txn}, Ptr{Cchar}, Cuint, Ptr{MDB_dbi}), txn, name, flags, dbi)
end

mdb_dbi_open(txn, name, flags, dbi) = checked_call(_mdb_dbi_open, txn, name, flags, dbi)

mutable struct MDB_env end

function _mdb_dbi_close(env, dbi)
    ccall((:mdb_dbi_close, liblmdb), Cvoid, (Ptr{MDB_env}, MDB_dbi), env, dbi)
end

mdb_dbi_close(env, dbi) = checked_call(_mdb_dbi_close, env, dbi)

const mdb_mode_t = mode_t

const mdb_filehandle_t = Cint

mutable struct MDB_cursor end

struct MDB_val
    mv_size::Csize_t
    mv_data::Ptr{Cvoid}
end

# typedef int ( MDB_cmp_func ) ( const MDB_val * a , const MDB_val * b )
const MDB_cmp_func = Cvoid

# typedef void ( MDB_rel_func ) ( MDB_val * item , void * oldptr , void * newptr , void * relctx )
const MDB_rel_func = Cvoid

@cenum MDB_cursor_op::UInt32 begin
    MDB_FIRST = 0
    MDB_FIRST_DUP = 1
    MDB_GET_BOTH = 2
    MDB_GET_BOTH_RANGE = 3
    MDB_GET_CURRENT = 4
    MDB_GET_MULTIPLE = 5
    MDB_LAST = 6
    MDB_LAST_DUP = 7
    MDB_NEXT = 8
    MDB_NEXT_DUP = 9
    MDB_NEXT_MULTIPLE = 10
    MDB_NEXT_NODUP = 11
    MDB_PREV = 12
    MDB_PREV_DUP = 13
    MDB_PREV_NODUP = 14
    MDB_SET = 15
    MDB_SET_KEY = 16
    MDB_SET_RANGE = 17
    MDB_PREV_MULTIPLE = 18
end

struct MDB_stat
    ms_psize::Cuint
    ms_depth::Cuint
    ms_branch_pages::Csize_t
    ms_leaf_pages::Csize_t
    ms_overflow_pages::Csize_t
    ms_entries::Csize_t
end

struct MDB_envinfo
    me_mapaddr::Ptr{Cvoid}
    me_mapsize::Csize_t
    me_last_pgno::Csize_t
    me_last_txnid::Csize_t
    me_maxreaders::Cuint
    me_numreaders::Cuint
end

function _mdb_version(major, minor, patch)
    ccall((:mdb_version, liblmdb), Ptr{Cchar}, (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), major, minor, patch)
end

mdb_version(major, minor, patch) = checked_call(_mdb_version, major, minor, patch)

function _mdb_strerror(err)
    ccall((:mdb_strerror, liblmdb), Ptr{Cchar}, (Cint,), err)
end

mdb_strerror(err) = checked_call(_mdb_strerror, err)

function _mdb_env_create(env)
    ccall((:mdb_env_create, liblmdb), Cint, (Ptr{Ptr{MDB_env}},), env)
end

mdb_env_create(env) = checked_call(_mdb_env_create, env)

function _mdb_env_open(env, path, flags, mode)
    ccall((:mdb_env_open, liblmdb), Cint, (Ptr{MDB_env}, Ptr{Cchar}, Cuint, mdb_mode_t), env, path, flags, mode)
end

mdb_env_open(env, path, flags, mode) = checked_call(_mdb_env_open, env, path, flags, mode)

function _mdb_env_copy(env, path)
    ccall((:mdb_env_copy, liblmdb), Cint, (Ptr{MDB_env}, Ptr{Cchar}), env, path)
end

mdb_env_copy(env, path) = checked_call(_mdb_env_copy, env, path)

function _mdb_env_copyfd(env, fd)
    ccall((:mdb_env_copyfd, liblmdb), Cint, (Ptr{MDB_env}, mdb_filehandle_t), env, fd)
end

mdb_env_copyfd(env, fd) = checked_call(_mdb_env_copyfd, env, fd)

function _mdb_env_copy2(env, path, flags)
    ccall((:mdb_env_copy2, liblmdb), Cint, (Ptr{MDB_env}, Ptr{Cchar}, Cuint), env, path, flags)
end

mdb_env_copy2(env, path, flags) = checked_call(_mdb_env_copy2, env, path, flags)

function _mdb_env_copyfd2(env, fd, flags)
    ccall((:mdb_env_copyfd2, liblmdb), Cint, (Ptr{MDB_env}, mdb_filehandle_t, Cuint), env, fd, flags)
end

mdb_env_copyfd2(env, fd, flags) = checked_call(_mdb_env_copyfd2, env, fd, flags)

function _mdb_env_stat(env, stat)
    ccall((:mdb_env_stat, liblmdb), Cint, (Ptr{MDB_env}, Ptr{MDB_stat}), env, stat)
end

mdb_env_stat(env, stat) = checked_call(_mdb_env_stat, env, stat)

function _mdb_env_info(env, stat)
    ccall((:mdb_env_info, liblmdb), Cint, (Ptr{MDB_env}, Ptr{MDB_envinfo}), env, stat)
end

mdb_env_info(env, stat) = checked_call(_mdb_env_info, env, stat)

function _mdb_env_sync(env, force)
    ccall((:mdb_env_sync, liblmdb), Cint, (Ptr{MDB_env}, Cint), env, force)
end

mdb_env_sync(env, force) = checked_call(_mdb_env_sync, env, force)

function _mdb_env_close(env)
    ccall((:mdb_env_close, liblmdb), Cvoid, (Ptr{MDB_env},), env)
end

mdb_env_close(env) = checked_call(_mdb_env_close, env)

function _mdb_env_set_flags(env, flags, onoff)
    ccall((:mdb_env_set_flags, liblmdb), Cint, (Ptr{MDB_env}, Cuint, Cint), env, flags, onoff)
end

mdb_env_set_flags(env, flags, onoff) = checked_call(_mdb_env_set_flags, env, flags, onoff)

function _mdb_env_get_flags(env, flags)
    ccall((:mdb_env_get_flags, liblmdb), Cint, (Ptr{MDB_env}, Ptr{Cuint}), env, flags)
end

mdb_env_get_flags(env, flags) = checked_call(_mdb_env_get_flags, env, flags)

function _mdb_env_get_path(env, path)
    ccall((:mdb_env_get_path, liblmdb), Cint, (Ptr{MDB_env}, Ptr{Ptr{Cchar}}), env, path)
end

mdb_env_get_path(env, path) = checked_call(_mdb_env_get_path, env, path)

function _mdb_env_get_fd(env, fd)
    ccall((:mdb_env_get_fd, liblmdb), Cint, (Ptr{MDB_env}, Ptr{mdb_filehandle_t}), env, fd)
end

mdb_env_get_fd(env, fd) = checked_call(_mdb_env_get_fd, env, fd)

function _mdb_env_set_mapsize(env, size)
    ccall((:mdb_env_set_mapsize, liblmdb), Cint, (Ptr{MDB_env}, Csize_t), env, size)
end

mdb_env_set_mapsize(env, size) = checked_call(_mdb_env_set_mapsize, env, size)

function _mdb_env_set_maxreaders(env, readers)
    ccall((:mdb_env_set_maxreaders, liblmdb), Cint, (Ptr{MDB_env}, Cuint), env, readers)
end

mdb_env_set_maxreaders(env, readers) = checked_call(_mdb_env_set_maxreaders, env, readers)

function _mdb_env_get_maxreaders(env, readers)
    ccall((:mdb_env_get_maxreaders, liblmdb), Cint, (Ptr{MDB_env}, Ptr{Cuint}), env, readers)
end

mdb_env_get_maxreaders(env, readers) = checked_call(_mdb_env_get_maxreaders, env, readers)

function _mdb_env_set_maxdbs(env, dbs)
    ccall((:mdb_env_set_maxdbs, liblmdb), Cint, (Ptr{MDB_env}, MDB_dbi), env, dbs)
end

mdb_env_set_maxdbs(env, dbs) = checked_call(_mdb_env_set_maxdbs, env, dbs)

function _mdb_env_get_maxkeysize(env)
    ccall((:mdb_env_get_maxkeysize, liblmdb), Cint, (Ptr{MDB_env},), env)
end

mdb_env_get_maxkeysize(env) = checked_call(_mdb_env_get_maxkeysize, env)

function _mdb_env_set_userctx(env, ctx)
    ccall((:mdb_env_set_userctx, liblmdb), Cint, (Ptr{MDB_env}, Ptr{Cvoid}), env, ctx)
end

mdb_env_set_userctx(env, ctx) = checked_call(_mdb_env_set_userctx, env, ctx)

function _mdb_env_get_userctx(env)
    ccall((:mdb_env_get_userctx, liblmdb), Ptr{Cvoid}, (Ptr{MDB_env},), env)
end

mdb_env_get_userctx(env) = checked_call(_mdb_env_get_userctx, env)

# typedef void MDB_assert_func ( MDB_env * env , const char * msg )
const MDB_assert_func = Cvoid

function _mdb_env_set_assert(env, func)
    ccall((:mdb_env_set_assert, liblmdb), Cint, (Ptr{MDB_env}, Ptr{Cvoid}), env, func)
end

mdb_env_set_assert(env, func) = checked_call(_mdb_env_set_assert, env, func)

function _mdb_txn_begin(env, parent, flags, txn)
    ccall((:mdb_txn_begin, liblmdb), Cint, (Ptr{MDB_env}, Ptr{MDB_txn}, Cuint, Ptr{Ptr{MDB_txn}}), env, parent, flags, txn)
end

mdb_txn_begin(env, parent, flags, txn) = checked_call(_mdb_txn_begin, env, parent, flags, txn)

function _mdb_txn_env(txn)
    ccall((:mdb_txn_env, liblmdb), Ptr{MDB_env}, (Ptr{MDB_txn},), txn)
end

mdb_txn_env(txn) = checked_call(_mdb_txn_env, txn)

function _mdb_txn_id(txn)
    ccall((:mdb_txn_id, liblmdb), Csize_t, (Ptr{MDB_txn},), txn)
end

mdb_txn_id(txn) = checked_call(_mdb_txn_id, txn)

function _mdb_txn_commit(txn)
    ccall((:mdb_txn_commit, liblmdb), Cint, (Ptr{MDB_txn},), txn)
end

mdb_txn_commit(txn) = checked_call(_mdb_txn_commit, txn)

function _mdb_txn_abort(txn)
    ccall((:mdb_txn_abort, liblmdb), Cvoid, (Ptr{MDB_txn},), txn)
end

mdb_txn_abort(txn) = checked_call(_mdb_txn_abort, txn)

function _mdb_txn_reset(txn)
    ccall((:mdb_txn_reset, liblmdb), Cvoid, (Ptr{MDB_txn},), txn)
end

mdb_txn_reset(txn) = checked_call(_mdb_txn_reset, txn)

function _mdb_txn_renew(txn)
    ccall((:mdb_txn_renew, liblmdb), Cint, (Ptr{MDB_txn},), txn)
end

mdb_txn_renew(txn) = checked_call(_mdb_txn_renew, txn)

function _mdb_stat(txn, dbi, stat)
    ccall((:mdb_stat, liblmdb), Cint, (Ptr{MDB_txn}, MDB_dbi, Ptr{MDB_stat}), txn, dbi, stat)
end

mdb_stat(txn, dbi, stat) = checked_call(_mdb_stat, txn, dbi, stat)

function _mdb_dbi_flags(txn, dbi, flags)
    ccall((:mdb_dbi_flags, liblmdb), Cint, (Ptr{MDB_txn}, MDB_dbi, Ptr{Cuint}), txn, dbi, flags)
end

mdb_dbi_flags(txn, dbi, flags) = checked_call(_mdb_dbi_flags, txn, dbi, flags)

function _mdb_drop(txn, dbi, del)
    ccall((:mdb_drop, liblmdb), Cint, (Ptr{MDB_txn}, MDB_dbi, Cint), txn, dbi, del)
end

mdb_drop(txn, dbi, del) = checked_call(_mdb_drop, txn, dbi, del)

function _mdb_set_compare(txn, dbi, cmp)
    ccall((:mdb_set_compare, liblmdb), Cint, (Ptr{MDB_txn}, MDB_dbi, Ptr{Cvoid}), txn, dbi, cmp)
end

mdb_set_compare(txn, dbi, cmp) = checked_call(_mdb_set_compare, txn, dbi, cmp)

function _mdb_set_dupsort(txn, dbi, cmp)
    ccall((:mdb_set_dupsort, liblmdb), Cint, (Ptr{MDB_txn}, MDB_dbi, Ptr{Cvoid}), txn, dbi, cmp)
end

mdb_set_dupsort(txn, dbi, cmp) = checked_call(_mdb_set_dupsort, txn, dbi, cmp)

function _mdb_set_relfunc(txn, dbi, rel)
    ccall((:mdb_set_relfunc, liblmdb), Cint, (Ptr{MDB_txn}, MDB_dbi, Ptr{Cvoid}), txn, dbi, rel)
end

mdb_set_relfunc(txn, dbi, rel) = checked_call(_mdb_set_relfunc, txn, dbi, rel)

function _mdb_set_relctx(txn, dbi, ctx)
    ccall((:mdb_set_relctx, liblmdb), Cint, (Ptr{MDB_txn}, MDB_dbi, Ptr{Cvoid}), txn, dbi, ctx)
end

mdb_set_relctx(txn, dbi, ctx) = checked_call(_mdb_set_relctx, txn, dbi, ctx)

function _mdb_get(txn, dbi, key, data)
    ccall((:mdb_get, liblmdb), Cint, (Ptr{MDB_txn}, MDB_dbi, Ptr{MDB_val}, Ptr{MDB_val}), txn, dbi, key, data)
end

mdb_get(txn, dbi, key, data) = checked_call(_mdb_get, txn, dbi, key, data)

function _mdb_put(txn, dbi, key, data, flags)
    ccall((:mdb_put, liblmdb), Cint, (Ptr{MDB_txn}, MDB_dbi, Ptr{MDB_val}, Ptr{MDB_val}, Cuint), txn, dbi, key, data, flags)
end

mdb_put(txn, dbi, key, data, flags) = checked_call(_mdb_put, txn, dbi, key, data, flags)

function _mdb_del(txn, dbi, key, data)
    ccall((:mdb_del, liblmdb), Cint, (Ptr{MDB_txn}, MDB_dbi, Ptr{MDB_val}, Ptr{MDB_val}), txn, dbi, key, data)
end

mdb_del(txn, dbi, key, data) = checked_call(_mdb_del, txn, dbi, key, data)

function _mdb_cursor_open(txn, dbi, cursor)
    ccall((:mdb_cursor_open, liblmdb), Cint, (Ptr{MDB_txn}, MDB_dbi, Ptr{Ptr{MDB_cursor}}), txn, dbi, cursor)
end

mdb_cursor_open(txn, dbi, cursor) = checked_call(_mdb_cursor_open, txn, dbi, cursor)

function _mdb_cursor_close(cursor)
    ccall((:mdb_cursor_close, liblmdb), Cvoid, (Ptr{MDB_cursor},), cursor)
end

mdb_cursor_close(cursor) = checked_call(_mdb_cursor_close, cursor)

function _mdb_cursor_renew(txn, cursor)
    ccall((:mdb_cursor_renew, liblmdb), Cint, (Ptr{MDB_txn}, Ptr{MDB_cursor}), txn, cursor)
end

mdb_cursor_renew(txn, cursor) = checked_call(_mdb_cursor_renew, txn, cursor)

function _mdb_cursor_txn(cursor)
    ccall((:mdb_cursor_txn, liblmdb), Ptr{MDB_txn}, (Ptr{MDB_cursor},), cursor)
end

mdb_cursor_txn(cursor) = checked_call(_mdb_cursor_txn, cursor)

function _mdb_cursor_dbi(cursor)
    ccall((:mdb_cursor_dbi, liblmdb), MDB_dbi, (Ptr{MDB_cursor},), cursor)
end

mdb_cursor_dbi(cursor) = checked_call(_mdb_cursor_dbi, cursor)

function _mdb_cursor_get(cursor, key, data, op)
    ccall((:mdb_cursor_get, liblmdb), Cint, (Ptr{MDB_cursor}, Ptr{MDB_val}, Ptr{MDB_val}, MDB_cursor_op), cursor, key, data, op)
end

mdb_cursor_get(cursor, key, data, op) = checked_call(_mdb_cursor_get, cursor, key, data, op)

function _mdb_cursor_put(cursor, key, data, flags)
    ccall((:mdb_cursor_put, liblmdb), Cint, (Ptr{MDB_cursor}, Ptr{MDB_val}, Ptr{MDB_val}, Cuint), cursor, key, data, flags)
end

mdb_cursor_put(cursor, key, data, flags) = checked_call(_mdb_cursor_put, cursor, key, data, flags)

function _mdb_cursor_del(cursor, flags)
    ccall((:mdb_cursor_del, liblmdb), Cint, (Ptr{MDB_cursor}, Cuint), cursor, flags)
end

mdb_cursor_del(cursor, flags) = checked_call(_mdb_cursor_del, cursor, flags)

function _mdb_cursor_count(cursor, countp)
    ccall((:mdb_cursor_count, liblmdb), Cint, (Ptr{MDB_cursor}, Ptr{Csize_t}), cursor, countp)
end

mdb_cursor_count(cursor, countp) = checked_call(_mdb_cursor_count, cursor, countp)

function _mdb_cmp(txn, dbi, a, b)
    ccall((:mdb_cmp, liblmdb), Cint, (Ptr{MDB_txn}, MDB_dbi, Ptr{MDB_val}, Ptr{MDB_val}), txn, dbi, a, b)
end

mdb_cmp(txn, dbi, a, b) = checked_call(_mdb_cmp, txn, dbi, a, b)

function _mdb_dcmp(txn, dbi, a, b)
    ccall((:mdb_dcmp, liblmdb), Cint, (Ptr{MDB_txn}, MDB_dbi, Ptr{MDB_val}, Ptr{MDB_val}), txn, dbi, a, b)
end

mdb_dcmp(txn, dbi, a, b) = checked_call(_mdb_dcmp, txn, dbi, a, b)

# typedef int ( MDB_msg_func ) ( const char * msg , void * ctx )
const MDB_msg_func = Cvoid

function _mdb_reader_list(env, func, ctx)
    ccall((:mdb_reader_list, liblmdb), Cint, (Ptr{MDB_env}, Ptr{Cvoid}, Ptr{Cvoid}), env, func, ctx)
end

mdb_reader_list(env, func, ctx) = checked_call(_mdb_reader_list, env, func, ctx)

function _mdb_reader_check(env, dead)
    ccall((:mdb_reader_check, liblmdb), Cint, (Ptr{MDB_env}, Ptr{Cint}), env, dead)
end

mdb_reader_check(env, dead) = checked_call(_mdb_reader_check, env, dead)

const MDB_VERSION_MAJOR = 0

const MDB_VERSION_MINOR = 9

const MDB_VERSION_PATCH = 27

#const MDB_VERSION_FULL = MDB_VERINT(MDB_VERSION_MAJOR, MDB_VERSION_MINOR, MDB_VERSION_PATCH)

const MDB_VERSION_DATE = "October 26, 2020"

#const MDB_VERSION_STRING = MDB_VERFOO(MDB_VERSION_MAJOR, MDB_VERSION_MINOR, MDB_VERSION_PATCH, MDB_VERSION_DATE)

const MDB_FIXEDMAP = 0x01

const MDB_NOSUBDIR = 0x4000

const MDB_NOSYNC = 0x00010000

const MDB_RDONLY = 0x00020000

const MDB_NOMETASYNC = 0x00040000

const MDB_WRITEMAP = 0x00080000

const MDB_MAPASYNC = 0x00100000

const MDB_NOTLS = 0x00200000

const MDB_NOLOCK = 0x00400000

const MDB_NORDAHEAD = 0x00800000

const MDB_NOMEMINIT = 0x01000000

const MDB_REVERSEKEY = 0x02

const MDB_DUPSORT = 0x04

const MDB_INTEGERKEY = 0x08

const MDB_DUPFIXED = 0x10

const MDB_INTEGERDUP = 0x20

const MDB_REVERSEDUP = 0x40

const MDB_CREATE = 0x00040000

const MDB_NOOVERWRITE = 0x10

const MDB_NODUPDATA = 0x20

const MDB_CURRENT = 0x40

const MDB_RESERVE = 0x00010000

const MDB_APPEND = 0x00020000

const MDB_APPENDDUP = 0x00040000

const MDB_MULTIPLE = 0x00080000

const MDB_CP_COMPACT = 0x01

const MDB_SUCCESS = 0

const MDB_KEYEXIST = -30799

const MDB_NOTFOUND = -30798

const MDB_PAGE_NOTFOUND = -30797

const MDB_CORRUPTED = -30796

const MDB_PANIC = -30795

const MDB_VERSION_MISMATCH = -30794

const MDB_INVALID = -30793

const MDB_MAP_FULL = -30792

const MDB_DBS_FULL = -30791

const MDB_READERS_FULL = -30790

const MDB_TLS_FULL = -30789

const MDB_TXN_FULL = -30788

const MDB_CURSOR_FULL = -30787

const MDB_PAGE_FULL = -30786

const MDB_MAP_RESIZED = -30785

const MDB_INCOMPATIBLE = -30784

const MDB_BAD_RSLOT = -30783

const MDB_BAD_TXN = -30782

const MDB_BAD_VALSIZE = -30781

const MDB_BAD_DBI = -30780

const MDB_LAST_ERRCODE = MDB_BAD_DBI

# exports
const PREFIXES = ["mdb_", "MDB_","_mdb_"]
for name in names(@__MODULE__; all=true), prefix in PREFIXES
    if startswith(string(name), prefix)
        @eval export $name
    end
end

end # module


text:
module ReadHTML

using DataFrames, HTTP, Gumbo, AbstractTrees
export gettables, gettablerow, extractdata, makerowdata, read_html, read_table
const StringMissing = Union{String, Missing}

"""
    gettables(url::String)

Extract all
~
 contents.
"""
function gettables(url::String)
    data = String(HTTP.get(url).body)
    data = replace(data, "\n"=>"")
    table_list = String[]
    offset = 1
    while true
        m = match(r"
(.*?)
"s, data, offset)
        if isnothing(m)
            break
        else
            push!(table_list, m.match)
            offset = m.offset + 1
        end
    end
    return table_list
end

"""
    gettablerow(table)

Extract all row data of table.
"""
function gettablerow(table, header)
    l = StringMissing[]
    offset = 1
    while true
        m = match(r"
(.*?)
"s, table, offset)
        if isnothing(m)
            break
        else
            offset = m.offset + 1
            if header
                push!(l, m.captures[1])
            elseif !occursin(r"", m.captures[1])
                push!(l, m.captures[1])
            end
        end
    end
    return l
end

function extractdata(data)
    data = replace(data, "\n"=>"")
    t = parsehtml(data)
    t = t.root[2]

    str = String[]
    for elem in PreOrderDFS(t)
        if isa(elem, HTMLText)
            push!(str, elem.text)
        end
    end
    return join(str, ' ')
end

"""
    makerowdata(tr)

Extract row elements.
"""
function makerowdata(tr)
    l = StringMissing[]
    offset = 1
    while true
        m = match(r"<t[dh].*?>(.*?)</t[dh]>"s, tr, offset)
        if isnothing(m)
            break
        else
            offset = m.offset + 1
            push!(l, extractdata(m.match))
        end
    end
    return l
end

"""
    checknumrows(trs)

Check how many rows there are.
"""
function checknumrows(trs)
    tr = trs[1]
    return length(makerowdata(tr))
end

"""
    read_table(table; header=true)

Read HTML table into DataFrame.
"""
function read_table(table; header=true)
    containheader = occursin("	", table)
    trs = gettablerow(table, header)
    isempty(trs) && return nothing
    ncol = checknumrows(trs)
    nrow = size(trs, 1)
    arr = Matrix{StringMissing}(undef, nrow, ncol)
    for (rowidx, tr) in enumerate(trs)
        rowdata = makerowdata(tr)
        for (colidx, colitem) in enumerate(rowdata)
            arr[rowidx, colidx] = colitem
        end
    end
    try
        if containheader && header
            DataFrame(arr[2:end, :], Symbol.(arr[1,:]))
        else
            DataFrame(arr)
        end
    catch
        return nothing
    end
end

"""
    read_html(url; header=true)

Obtain HTML tables from `url` and read them into DataFrame.
"""
function read_html(url; header=true)
    tables = gettables(url)
    dfs = DataFrame[]
    for table in tables
        tb = read_table(table, header=header)
        !isnothing(tb) && push!(dfs, tb)
    end
    return dfs
end

end # module


text:
using Documenter, MonteCarloMeasurements, Unitful

using Plots
plotly()

makedocs(
      sitename = "MonteCarloMeasurements Documentation",
      doctest = false,
      modules = [MonteCarloMeasurements],
      pages = [
            "Home" => "index.md",
            "Supporting new functions" => "overloading.md",
            "Examples" => "examples.md",
            "Linear vs. Monte-Carlo uncertainty propagation" => "comparison.md",
            "Performance tips" => "performance.md",
            "Advanced usage" => "advanced_usage.md",
            "API" => "api.md",
      ],
      format = Documenter.HTML(prettyurls = haskey(ENV, "CI")),
) # Due to lots of plots, this will just have to be run on my local machine

deploydocs(
      deps = Deps.pip("pygments", "mkdocs", "python-markdown-math", "mkdocs-cinder"),
      repo = "github.com/baggepinnen/MonteCarloMeasurements.jl.git",
)


text:
function resultConverter_v3(x, V, anchoLado, matConexionVertices, vecVertices, vecAlturas)

    template = Int(x[end])
    alt = min(x[1], maximum(vecAlturas))
    theta = x[2]

    # Genera Poligono que corta VolTeor a altura alt
    idAlt = findfirst(y -> y >= alt, vecAlturas) - 1
    vert_poly_0 = Int.(vecVertices[idAlt])
    vert_poly_1 = Int.(matConexionVertices[vert_poly_0,2])
    poly_0 = copy(V[vert_poly_0, 1:2])[:,1,:]
    poly_1 = copy(V[vert_poly_1, 1:2])[:,1,:]
    alfa_ = (vecAlturas[idAlt + 1] - alt) / (vecAlturas[idAlt + 1] - vecAlturas[idAlt])
    polyCorte_alt = alfa_ .* poly_0 .+ (1 - alfa_) .* poly_1

    if template == 1

        alfa = x[3]
        pos_x = x[4]
        pos_y = x[5]

        largo1 = anchoLado #
        largo2 = anchoLado #

        R1 = poly2D.rotationMatrix(theta);
        cr1 = [pos_x; pos_y]
        p1_1 = [pos_x; pos_y];
        p2_1 = R1 * ([pos_x + largo1; pos_y] - cr1) + cr1;
        p3_1 = R1 * ([pos_x + largo1; pos_y + anchoLado] - cr1) + cr1;
        p4_1 = R1 * ([pos_x ; pos_y + anchoLado] - cr1) + cr1;
        R2 = poly2D.rotationMatrix(alfa + theta);
        cr2 = [pos_x; pos_y]
        p1_2 = [pos_x; pos_y];
        p2_2 = R2 * ([pos_x + anchoLado; pos_y] - cr2) + cr2;
        p3_2 = R2 * ([pos_x + anchoLado; pos_y + largo2] - cr2) + cr2;
        p4_2 = R2 * ([pos_x; pos_y + largo2] - cr2) + cr2;
        V1 = [p1_1';p2_1';p3_1';p4_1']
        V2 = [p1_2';p2_2';p3_2';p4_2'];

        ps1 = PolyShape([V1], 1)
        ps2 = PolyShape([V2], 1)
        psCorte = PolyShape([polyCorte_alt], 1)

        flagV1 = polyShape.isPolyInPoly(ps1, psCorte)
        flagV2 = polyShape.isPolyInPoly(ps2, psCorte)

        if flagV1 && flagV2
            V1_ = poly2D.expandPolygonSide(V1, 500, 2)
            pV1 = V1[[2, 3],:]
            VX1 = poly2D.intersectPoly2d(polyCorte_alt, V1_)
            DM1 = poly2D.distanceMat(pV1, VX1[:,[3,4]])
            minD1 = minimum(DM1)
            V1__ = poly2D.expandPolygonSide(V1, minD1, 2)
            ps1max = PolyShape([V1__], 1)

            V2_ = poly2D.expandPolygonSide(V2, 500, 3)
            pV2 = V2[[3, 4],:]
            VX2 = poly2D.intersectPoly2d(polyCorte_alt, V2_)
            DM2 = poly2D.distanceMat(pV2, VX2[:,[3,4]])
            minD2 = minimum(DM2)
            V2__ = poly2D.expandPolygonSide(V2, minD2, 3)
            ps2max = PolyShape([V2__], 1)

            ps_base = polyShape.polyUnion_v2(ps1max, ps2max)
            areaBasal = polyShape.polyArea_v2(ps_base)
            ps_baseSeparada = PolyShape([V1__, V2__], 2)

        else
            areaBasal = -1000
            ps_base = []
            ps_baseSeparada = []
        end

    elseif template == 2
        phi1 = x[3]
        phi2 = x[4]
        pos_x0 = x[5]
        pos_y0 = x[6]
        largo0 = x[7]

        largo1 = anchoLado #
        largo2 = anchoLado #

        R_theta = poly2D.rotationMatrix(theta);
        cr_theta  = [pos_x0; pos_y0];

        p1_0  = (R_theta * ([pos_x0; pos_y0] - cr_theta) + cr_theta)'
        p2_0  = (R_theta * ([pos_x0 + largo0; pos_y0] - cr_theta) + cr_theta)'
        p3_0  = (R_theta * ([pos_x0 + largo0; pos_y0 + anchoLado] - cr_theta) + cr_theta)'
        p4_0  = (R_theta * ([pos_x0; pos_y0 + anchoLado] - cr_theta) + cr_theta)'

        R_phi1 = poly2D.rotationMatrix(phi1);
        cr_phi1  = [pos_x0; pos_y0];
        p1_1  = (R_theta * ((R_phi1 * ([pos_x0; pos_y0] - cr_phi1) + cr_phi1) - cr_theta) + cr_theta)';
        p2_1  = (R_theta * ((R_phi1 * ([pos_x0 + anchoLado; pos_y0] - cr_phi1) + cr_phi1) - cr_theta) + cr_theta)';
        p3_1  = (R_theta * ((R_phi1 * ([pos_x0 + anchoLado; pos_y0 + largo1] - cr_phi1) + cr_phi1) - cr_theta) + cr_theta)';
        p4_1  = (R_theta * ((R_phi1 * ([pos_x0; pos_y0 + largo1] - cr_phi1) + cr_phi1) - cr_theta) + cr_theta)';

        R_phi2 = poly2D.rotationMatrix(-phi2);
        cr_phi2  = [pos_x0 + largo0; pos_y0];
        p1_2  = (R_theta * ((R_phi2 * ([pos_x0 + largo0 - anchoLado; pos_y0] - cr_phi2) + cr_phi2) - cr_theta) + cr_theta)';
        p2_2  = (R_theta * ((R_phi2 * ([pos_x0 + largo0; pos_y0] - cr_phi2) + cr_phi2) - cr_theta) + cr_theta)';
        p3_2  = (R_theta * ((R_phi2 * ([pos_x0 + largo0; pos_y0 + largo2] - cr_phi2) + cr_phi2) - cr_theta) + cr_theta)';
        p4_2  = (R_theta * ((R_phi2 * ([pos_x0 + largo0 - anchoLado; pos_y0 + largo2] - cr_phi2) + cr_phi2) - cr_theta) + cr_theta)';

        V0 = [p1_0;p2_0;p3_0;p4_0]
        V1 = [p1_1;p2_1;p3_1;p4_1]
        V2 = [p1_2;p2_2;p3_2;p4_2]

        ps0 = PolyShape([V0], 1)
        ps1 = PolyShape([V1], 1)
        ps2 = PolyShape([V2], 1)
        psCorte = PolyShape([polyCorte_alt], 1)

        flagV0 = polyShape.isPolyInPoly(ps0, psCorte)
        flagV1 = polyShape.isPolyInPoly(ps1, psCorte)
        flagV2 = polyShape.isPolyInPoly(ps2, psCorte)

        if flagV0 && flagV1 && flagV2
            V1_ = poly2D.expandPolygonSide(V1, 500, 3)
            pV1 = V1[[3, 4],:]
            VX1 = poly2D.intersectPoly2d(polyCorte_alt, V1_)
            DM1 = poly2D.distanceMat(pV1, VX1[:,[3,4]])
            minD1 = minimum(DM1)
            V1__ = poly2D.expandPolygonSide(V1, minD1, 3)
            ps1max = PolyShape([V1__], 1)

            V2_ = poly2D.expandPolygonSide(V2, 500, 3)
            pV2 = V2[[3, 4],:]
            VX2 = poly2D.intersectPoly2d(polyCorte_alt, V2_)
            DM2 = poly2D.distanceMat(pV2, VX2[:,[3,4]])
            minD2 = minimum(DM2)
            V2__ = poly2D.expandPolygonSide(V2, minD2, 3)
            ps2max = PolyShape([V2__], 1)

            ps_base = polyShape.polyUnion_v2(ps0, ps1max)
            ps_base = polyShape.polyUnion_v2(ps_base, ps2max)
            areaBasal = polyShape.polyArea_v2(ps_base)
            ps_baseSeparada = PolyShape([V0, V1__, V2__], 3)

        else
            areaBasal = -1000
            ps_base = []
            ps_baseSeparada = []
        end

    elseif template == 3
        pos_x0 = x[3]
        pos_y0 = x[4]
        largo0 = x[5]

        R_theta = poly2D.rotationMatrix(theta);
        cr_theta  = [pos_x0; pos_y0];

        p1_0  = (R_theta * ([pos_x0; pos_y0] - cr_theta) + cr_theta)'
        p2_0  = (R_theta * ([pos_x0 + anchoLado; pos_y0] - cr_theta) + cr_theta)'
        p3_0  = (R_theta * ([pos_x0 + anchoLado; pos_y0 + largo0] - cr_theta) + cr_theta)'
        p4_0  = (R_theta * ([pos_x0; pos_y0 + largo0] - cr_theta) + cr_theta)'

        p1_1  = (R_theta * ([pos_x0; pos_y0 - anchoLado] - cr_theta) + cr_theta)'
        p2_1  = (R_theta * ([pos_x0 + anchoLado; pos_y0 - anchoLado] - cr_theta) + cr_theta)'
        p3_1  = (R_theta * ([pos_x0 + anchoLado; pos_y0] - cr_theta) + cr_theta)'
        p4_1  = (R_theta * ([pos_x0; pos_y0] - cr_theta) + cr_theta)'

        p1_2  = (R_theta * ([pos_x0; pos_y0 + largo0] - cr_theta) + cr_theta)'
        p2_2  = (R_theta * ([pos_x0 + anchoLado; pos_y0 + largo0] - cr_theta) + cr_theta)'
        p3_2  = (R_theta * ([pos_x0 + anchoLado; pos_y0 + largo0 + anchoLado] - cr_theta) + cr_theta)'
        p4_2  = (R_theta * ([pos_x0; pos_y0 + largo0 + anchoLado] - cr_theta) + cr_theta)'

        V0 = [p1_0;p2_0;p3_0;p4_0]
        V1 = [p1_1;p2_1;p3_1;p4_1]
        V2 = [p1_2;p2_2;p3_2;p4_2]

        ps0 = PolyShape([V0], 1)
        ps1 = PolyShape([V1], 1)
        ps2 = PolyShape([V2], 1)
        psCorte = PolyShape([polyCorte_alt], 1)

        flagV0 = polyShape.isPolyInPoly(ps0, psCorte)
        flagV1 = polyShape.isPolyInPoly(ps1, psCorte)
        flagV2 = polyShape.isPolyInPoly(ps2, psCorte)

        if flagV0 && flagV1 && flagV2
            V1_ = poly2D.expandPolygonSide(V1, 500, 2)
            VX1 = poly2D.intersectPoly2d(polyCorte_alt, V1_)
            pV1 = V1[[2, 3],:]
            DM1 = poly2D.distanceMat(pV1, VX1[:,[3,4]])
            minD1 = minimum(DM1)
            V1 = poly2D.expandPolygonSide(V1, minD1, 2)
            V1_ = poly2D.expandPolygonSide(V1, 500, 4)
            VX1 = poly2D.intersectPoly2d(polyCorte_alt, V1_)
            pV1 = V1[[4, 1],:]
            DM1 = poly2D.distanceMat(pV1, VX1[:,[3,4]])
            minD1 = minimum(DM1)
            V1 = poly2D.expandPolygonSide(V1, minD1, 4)
            ps1max = PolyShape([V1], 1)

            V2_ = poly2D.expandPolygonSide(V2, 500, 2)
            VX2 = poly2D.intersectPoly2d(polyCorte_alt, V2_)
            pV2 = V2[[2, 3],:]
            DM2 = poly2D.distanceMat(pV2, VX2[:,[3,4]])
            minD2 = minimum(DM2)
            V2 = poly2D.expandPolygonSide(V2, minD2, 2)
            V2_ = poly2D.expandPolygonSide(V2, 500, 4)
            VX2 = poly2D.intersectPoly2d(polyCorte_alt, V2_)
            pV2 = V2[[4, 1],:]
            DM2 = poly2D.distanceMat(pV2, VX2[:,[3,4]])
            minD2 = minimum(DM2)
            V2 = poly2D.expandPolygonSide(V2, minD2, 4)
            ps2max = PolyShape([V2], 1)

            ps_base = polyShape.polyUnion_v2(ps0, ps1max)
            ps_base = polyShape.polyUnion_v2(ps_base, ps2max)
            areaBasal = polyShape.polyArea_v2(ps_base)
            ps_baseSeparada = PolyShape([V0, V1, V2], 3)

        else
            areaBasal = -1000
            ps_base = []
            ps_baseSeparada = []
        end

    elseif template == 4

        alfa = x[3]
        pos_x = x[4]
        pos_y = x[5]

        largo1 = anchoLado #
        largo2 = anchoLado #

        R1 = poly2D.rotationMatrix(theta);
        cr1 = [pos_x; pos_y]
        p1_1 = [pos_x; pos_y];
        p2_1 = R1 * ([pos_x + largo1; pos_y] - cr1) + cr1;
        p3_1 = R1 * ([pos_x + largo1; pos_y + anchoLado] - cr1) + cr1;
        p4_1 = R1 * ([pos_x ; pos_y + anchoLado] - cr1) + cr1;
        R2 = poly2D.rotationMatrix(theta - alfa);
        cr2 = [pos_x; pos_y]
        p1_2 = R2 * ([pos_x - anchoLado; pos_y] - cr2) + cr2;
        p2_2 = R2 * ([pos_x; pos_y] - cr2) + cr2;
        p3_2 = R2 * ([pos_x; pos_y + largo2] - cr2) + cr2;
        p4_2 = R2 * ([pos_x - anchoLado; pos_y + largo2] - cr2) + cr2;
        V1 = [p1_1';p2_1';p3_1';p4_1']
        V2 = [p1_2';p2_2';p3_2';p4_2'];

        ps1 = PolyShape([V1], 1)
        ps2 = PolyShape([V2], 1)
        psCorte = PolyShape([polyCorte_alt], 1)

        flagV1 = polyShape.isPolyInPoly(ps1, psCorte)
        flagV2 = polyShape.isPolyInPoly(ps2, psCorte)

        V1_ = poly2D.expandPolygonSide(V1, 500, 2)
        VX1 = poly2D.intersectPoly2d(polyCorte_alt, V1_)

        V2_ = poly2D.expandPolygonSide(V2, 500, 3)
        VX2 = poly2D.intersectPoly2d(polyCorte_alt, V2_)

        if flagV1 && flagV2
            pV1 = V1[[2, 3],:]
            DM1 = poly2D.distanceMat(pV1, VX1[:,[3,4]])
            minD1 = minimum(DM1)
            V1__ = poly2D.expandPolygonSide(V1, minD1, 2)
            ps1max = PolyShape([V1__], 1)

            pV2 = V2[[3, 4],:]
            DM2 = poly2D.distanceMat(pV2, VX2[:,[3,4]])
            minD2 = minimum(DM2)
            V2__ = poly2D.expandPolygonSide(V2, minD2, 3)
            ps2max = PolyShape([V2__], 1)

            ps_base = polyShape.polyUnion_v2(ps1max, ps2max)
            areaBasal = polyShape.polyArea_v2(ps_base)
            ps_baseSeparada = PolyShape([V1__, V2__], 2)

        else
            areaBasal = -1000
            ps_base = []
            ps_baseSeparada = []
        end

    end

    return alt, areaBasal, ps_base, ps_baseSeparada, polyCorte_alt

end


text:
using Revise

import KRLS
using Test

# filter_param = KRLS.FilterParam(2, 0.01, 10)

a = 1


text:
using StatsBase
using Base.Test

srand(10)
v1 = rand(3)
v2 = ["Good", "Great", "Bad"]
v3 = rand(Int8, 3)
v4 = [StatsBase.PValue(rand()./10000) for i in 1:3]
m = rand(3,4)
@test sprint(show, CoefTable(Any[v1, v2, v3, v4],
    ["Estimate", "Comments", "df", "p"],
    ["x1", "x2", "x3"])) == """\
     Estimate Comments  df     p
x1   0.112582     Good  88 <1e-4
x2   0.368314    Great -90 <1e-4
x3   0.344454      Bad -80 <1e-4
"""

@test sprint(show, CoefTable(m, ["Estimate", "Stderr", "df", "p"],
    ["x1", "x2", "x3"], 4)) == """\
     Estimate   Stderr       df      p
x1   0.819778 0.844007 0.923676 0.1717
x2   0.669931  0.67919 0.066098 0.4204
x3   0.453058  0.72525 0.999172 0.5567
"""

@test sprint(show, StatsBase.PValue(1.0)) == "1.0000"
@test sprint(show, StatsBase.PValue(1e-1)) == "0.1000"
@test sprint(show, StatsBase.PValue(1e-5)) == "<1e-4"
@test sprint(show, StatsBase.PValue(NaN)) == "NaN"
@test_throws ErrorException StatsBase.PValue(-0.1)
@test_throws ErrorException StatsBase.PValue(1.1)


text:
function scaleto8(G)
	# Read a grid with GMT and return an image object of it scaled to [1 255]
	add_off = 1;
	range = 254. / (G.hdr[6] - G.hdr[5])
	img = zeros(UInt8,  G.n_rows, G.n_columns)
	GMT_grid_flip_vertical(convert(Ptr{Void}, pointer(G.z)), G.n_columns, G.n_rows, 0, sizeof(G.z[1]))
	for (k = 1:G.n_columns * G.n_rows)
		img[k] = round(UInt8, (((G.z[k] - G.hdr[5]) * range) + add_off))
	end

	P = gmt("makecpt -T0/255/1")
	nColors = size(P.colormap,1)
	cm = zeros(UInt32, nColors, 3)
	for (k = 1:nColors*3)
		cm[k] = round(UInt32, P.colormap[k]*255)
	end
	pal = zeros(UInt32, nColors)
	for (k = 1:nColors)
		pal[k] =  (UInt32(0) << 24) | (UInt32(cm[k]) << 16) | (UInt32(cm[k+nColors]) << 8) | UInt32(cm[k+2*nColors])
	end
	pPalette = pointer(pal)

	n = G.n_columns * G.n_rows
	data = pointer([pointer(img)])
	im = imImage(G.n_columns, G.n_rows, IM_MAP, 0, 0, 1, G.n_columns, n,
	             n, n, data, pPalette, nColors, C_NULL)
	image = pointer([im])

	return image
end


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

struct DensityProduct{
    N,
    D<:NTuple{N,AbstractDensity},
    B<:AbstractParamBounds
} <: AbstractDensity
    densities::D
    bounds::B
end

# ToDo: Check for equal number of parameters in DensityProduct ctor.

DensityProduct(densities::D, bounds::B) where {N,D<:NTuple{N,AbstractDensity},B<:AbstractParamBounds} =
    BAT.DensityProduct{N,D,B}(densities, bounds)

Base.parent(density::DensityProduct) = density.densities

param_bounds(density::DensityProduct) = density.bounds

nparams(density::DensityProduct) = nparams(density.bounds)

import Base.*
function *(a::AbstractDensity, b::AbstractDensity)
    nparams(a) != nparams(b) && throw(ArgumentError("Can't multiply densities with different number of arguments"))
    new_bounds = param_bounds(a) ∩ param_bounds(b)
    _unsafe_prod(a, b, new_bounds)
end

_unsafe_prod(a::AbstractDensity, b::AbstractDensity, new_bounds::AbstractParamBounds) =
    DensityProduct((a,b), new_bounds)

_unsafe_prod(a::AbstractDensity, b::DensityProduct, new_bounds::AbstractParamBounds) =
    DensityProduct((a, b.densities...), new_bounds)

_unsafe_prod(a::DensityProduct, b::DensityProduct, new_bounds::AbstractParamBounds) =
    DensityProduct((a.densities..., b.densities...), new_bounds)

_unsafe_prod(a::DensityProduct, b::AbstractDensity, new_bounds::AbstractParamBounds) =
    DensityProduct((a.densities...,b), new_bounds)

function unsafe_density_logval(density::DensityProduct, params::AbstractVector{<:Real},
    exec_context::ExecContext
)
    ds = density.densities
    isempty(ds) && throw(ArgumentError("Can't evaluate density_logval on empty DensityProduct"))
    sum(map(d -> unsafe_density_logval(d, params, exec_context), ds))
end

exec_capabilities(::typeof(unsafe_density_logval), density::DensityProduct, params::AbstractVector{<:Real}) =
    ∩(map(d -> exec_capabilities(density_logval, d, params), density.densities)...)

function unsafe_density_logval!(r::AbstractArray{<:Real}, density::DensityProduct, params::AbstractMatrix{<:Real},
    exec_context::ExecContext
)
    ds = density.densities
    isempty(ds) && throw(ArgumentError("Can't evaluate density_logval! on empty DensityProduct"))
    fill!(r, 0)
    tmp = similar(r)  # ToDo: Avoid memory allocation
    for d in ds
        unsafe_density_logval!(tmp, d, params, exec_context)
        r .+= tmp
    end
    r
end

exec_capabilities(::typeof(unsafe_density_logval!), r::AbstractArray{<:Real}, density::DensityProduct, params::AbstractMatrix{<:Real}) =
    ∩(map(d -> exec_capabilities(unsafe_density_logval!, r, d, params), density.densities)...)


text:
Base.@kwdef struct RoomParams
    max_allowed_objects = 8
    bottom_left_lim::Vector{Float64} = [-5, -5, 0]
    top_right_lim::Vector{Float64} = [5, 5, 3]
end

Base.@kwdef struct CameraParams
    image_dim_x = 256
    image_dim_y = 256
    horizontal_fov::Float64 = 55
    vertical_fov::Float64 = 55
end

export RoomParams
export CameraParams
export LineSegment

max_stars_repo_name: zhangir-azerbayev/MetaGen-New

max_issues_repo_name: zhangir-azerbayev/MetaGen-New

max_forks_repo_name: zhangir-azerbayev/MetaGen-New


text:
"""
    map_column_pairs(f, X::Matrix{Float64}[, Y::Matrix{Float64} = X])

Create a matrix by applying function `f` to each pair of columns of input matrices
`X` and `Y`.
"""
function map_column_pairs(f, X::AbstractMatrix, Y::AbstractMatrix)
    nobs1 = size(X,2)
    nobs2 = size(Y,2)
    D = Array{eltype(Y)}(undef, nobs1, nobs2)
    map_column_pairs!(D, f, X, Y)
    return D
end

function map_column_pairs(f, X::AbstractMatrix)
    dim, nobs = size(X)
    D = Array{eltype(X)}(undef, nobs, nobs)
    for i in 1:nobs
        for j in 1:i
            @inbounds D[i,j] = f(X[:,i], X[:,j])
            if i != j; @inbounds D[j,i] = D[i,j]; end;
        end
    end
    return D
end

"""
    map_column_pairs!(D::Matrix{Float64}, f, X::Matrix{Float64}[, Y::Matrix{Float64} = X])

Like [`map_column_pairs`](@ref), but stores the result in `D` rather than a new matrix.
"""
function map_column_pairs!(D::AbstractMatrix, f, X::AbstractMatrix, Y::AbstractMatrix)
    dim, nobs1 = size(X)
    nobs2 = size(Y,2)
    dim == size(Y,1) || throw(ArgumentError("Input observation matrices must have consistent dimensions"))
    size(D,1) == nobs1 || throw(ArgumentError(@sprintf("D has %d rows, while X has %d columns (should be same)",
                                                       size(D,1), nobs1)))
    size(D,2) == nobs2 || throw(ArgumentError(@sprintf("D has %d columns, while Y has %d columns (should be same)",
                                                       size(D,2), nobs2)))
    for i in 1:nobs1, j in 1:nobs2
        @inbounds D[i,j] = f(view(X,:,i), view(Y,:,j))
    end
    return D
end

function map_column_pairs!(D::AbstractMatrix, f, X::AbstractMatrix)
    dim, nobs = size(X)
    size(D,1) == nobs || throw(ArgumentError(@sprintf("D has %d rows, while X has %d columns (should be same)",
                                                       size(D,1), nobs)))
    size(D,2) == nobs || throw(ArgumentError(@sprintf("D has %d columns, while X has %d columns (should be same)",
                                                       size(D,2), nobs)))
    @inbounds for i in 1:nobs
        for j in 1:i
            D[i,j] = f(view(X,:,i), view(X,:,j))
            if i != j; D[j,i] = D[i,j]; end;
        end
    end
    return D
end


text:
module LatLongDF

using DataFrames

export distanceKm, LatLongS

const pid180 = pi / 180

struct LatLongS
	lat
	long
	LatLongS(lat, long) = new(lat * pid180, long * pid180)
end

# see https://spatialthoughts.com/2013/07/06/calculate-distance-spreadsheet/

function distanceKm(ll1::LatLongS, ll2::LatLongS)
	df = 2 * 6371 * asin(sqrt((sin((ll2.lat-ll1.lat)/2))^2+cos(ll2.lat)*cos(ll1.lat)*sin(((ll2.long-ll1.long)/2))^2))
end

function distanceKm(startdf::DataFrame, enddf::DataFrame)
	startp = eachrow(startdf[!, [:Latitude, :Longitude]])[1]
	endp = eachrow(enddf[!, [:Latitude, :Longitude]])[1]
	distanceKm(LatLongS(startp[1], startp[2]), LatLongS(endp[1], endp[2]))
end

function distanceKm(lat, long, enddf::DataFrame)
	endp = eachrow(enddf[!, [:Latitude, :Longitude]])[1]
	distanceKm(LatLongS(lat, long), LatLongS(endp[1], endp[2]))
end

function distanceKm(lat1, long1, lat2, long2)
	distanceKm(LatLongS(lat1, long1), LatLongS(lat2, long2))
end

###
end


text:
module types

FL = Float64

CustomSquareSparseMatrixRawPairs = Array{Pair{String, FL}}

mutable struct SquareSparseMatrix

    dict::Dict

    """
    ### Initializes new `CustomSquareSparseMatrix` object from an array of pairs.

    """
    function SquareSparseMatrix(init::CustomSquareSparseMatrixRawPairs)
        return new(Dict(init))
    end

    function SquareSparseMatrix(init::Dict)
        return new(init)
    end

end

# define some utility methods for this new type
function Base.getindex(obj::SquareSparseMatrix, i, j)
    return get(obj.dict, string(i, " ", j), 0)
end

function Base.setindex!(obj::SquareSparseMatrix, v, i, j)
    key = string(i, " ", j)
    # if v == 0
    #     delete!(obj.dict, key)
    # else
        obj.dict[key] = v
    # end
end

function Base.copy(obj::SquareSparseMatrix)
    dict = copy(obj.dict)
    return SquareSparseMatrix(dict)
end

export FL, SquareSparseMatrix, CustomSquareSparseMatrixRawPairs

end

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
export iauTttcg
"""
Time scale transformation:  Terrestrial Time, TT, to Geocentric
Coordinate Time, TCG.

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

Status:  canonical.

Given:
   tt1,tt2    double    TT as a 2-part Julian Date

Returned:
   tcg1,tcg2  double    TCG as a 2-part Julian Date

Returned (function value):
               int       status:  0 = OK

Note:

   tt1+tt2 is Julian Date, apportioned in any convenient way between
   the two arguments, for example where tt1 is the Julian Day Number
   and tt2 is the fraction of a day.  The returned tcg1,tcg2 follow
   suit.

References:

   McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
   IERS Technical Note No. 32, BKG (2004)

   IAU 2000 Resolution B1.9

This revision:  2013 June 18

SOFA release 2018-01-30

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

# int iauTttcg(double tt1, double tt2, double *tcg1, double *tcg2)
function iauTttcg(tt1::Real, tt2::Real)
   ref_tcg1 = Ref{Float64}(0.0)
   ref_tcg2 = Ref{Float64}(0.0)

   status = ccall((:iauTttcg, libsofa_c), Cint,
                   (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
                   convert(Float64, tt1), convert(Float64, tt2), ref_tcg1, ref_tcg2)

   return status, ref_tcg1[], ref_tcg2[]
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 2cf61e2a-ff3f-11ea-2d5b-77a2b52fd6e7
using CSV, DataFrames

# ╔═╡ 78701408-ff37-11ea-04a8-07805ef3cc1c
md"""
## Outline

In this project, we will implement a sowing plan for a mixed culture of three types of crops. At our disposal, we have fifteen plots of land, to be planted using a limited supply of seeds. The different lots have varying soil compositions, influencing plant growth. Furthermore, some combinations of the crops grow well together; other combinations show competition. These aspects make for a challenging optimization problem!
"""

# ╔═╡ bea5f54e-ff37-11ea-2c37-8707c75a9b5c
n = 15

# ╔═╡ 9b5ec9c8-ff37-11ea-1816-0365aa76c976
md"""
## A model for plant growth

We have $n=15$ plots of land. The three types of plants are denoted with x, y, and z, respectively. Indices are used to indicate the amount of seed distributed to a particular plot, e.g., $x_i$ is the amount of seed assigned to the $i$-th plot. So this optimization problem can be solved in a $15\times 3=45$-dimensional space.

Small-capped variables indicate the amount of seed assigned to a field. Large-cap variables represent the corresponding yield of that field, i.e., $X_i$ is the yield (in kg) for field $i$. The yield for each type of seed can be computed using the following equations:

$$X_i = \frac{A_i^x x_i^2}{120 + x_i^2 + 2y_i - 0.8x_iz_i + z_i^2}\,,$$

$$Y_i = \frac{A_i^y y_i^2}{30 + 4x_i + y_i^2 + 7z_i}\,,$$

$$Z_i = \frac{A_i^z z_i^2}{80 + 0.4x_i^2 + 0.2x_iz_i +0.6y_i + z_i^2}\,.$$
"""

# ╔═╡ d327e0c2-ff37-11ea-3006-53be838b7dc8
md"""
Here, $A_i^x$, $A_i^y$, $A_i^z$ represent the maximum in yields in field $i$ for the different seed types. It depends on the amount of nitrogen $u_i$ and the water status $v_i$ of the soil. The following equation can be used to compute these coefficients:

$$\begin{bmatrix}
A_i^x \\
A_i^y\\
A_i^z
\end{bmatrix}
=
\begin{bmatrix}
4 & 2  \\
1 & 0.3 \\
-0.5 & 4
\end{bmatrix}
\begin{bmatrix}
u_i \\
v_j\\
\end{bmatrix}
+
\begin{bmatrix}
100 \\
300 \\
210
\end{bmatrix}\,.$$
"""

# ╔═╡ c9897d6a-ff47-11ea-1212-9d5eb14cd0d7
S, b = [4 2; 1 0.3; -0.5 4], [100, 300, 210]  # soil effects

# ╔═╡ e4a29d80-ff37-11ea-0d0b-59545bce9bd2
md"""
The growth model has several interesting facets:
- sowing more seed results in a larger yield, however, the effect saturates, increasing quantities have diminishing returns;
- the plants show competition;
- z positively influences x, while x negatively influences z (z produces nitrogen while x requires much nitrogen);
- y is a good producer but very sensitive to adverse conditions.
"""

# ╔═╡ ea977544-ff37-11ea-310b-79307ef120d3
md"Below are the concentrations of nitrogen ($u$) and water status ($v$) for every field."

# ╔═╡ eba6a9dc-ff37-11ea-233d-1b9aca686ea4
u = [36.6776, 36.9967, 83.033, 50.3725, 43.4616, 55.5842, 44.8919, 99.6519, 20.158, 102.325, 96.8896, 33.7957, 26.6129, 38.7194, 60.1461]

# ╔═╡ 0e9a6d52-ff38-11ea-2472-810a521cb8f8
v = [34.5773,  24.3003,  24.3952,  28.462,  37.2912,  38.196,  36.4821,  30.1988,  20.9124,  35.207,  38.0924,  24.438,  28.3169,  20.3022,  24.8884]

# ╔═╡ 100a63f4-ff38-11ea-1fbf-05c136011e28
md"The total amount of seed for every type is fixed by $c_x$, $c_y$, $c_z$:"

# ╔═╡ 17088b36-ff38-11ea-3519-85055dd00929
cx, cy, cz = 250, 175, 325

# ╔═╡ 25b72854-ff38-11ea-0817-0f44e6420882
md"Finally, the yield for each crop can be sold at different prices (EUR/kg) $w_x$, $w_y$, $w_z$:"

# ╔═╡ 223cec90-ff38-11ea-040f-47c9ee1c9036
wx, wy, wz = 0.7, 0.85, 0.6

# ╔═╡ 3282d3c6-ff38-11ea-0458-af74eb6b42dc
md"""
So the objective is

$$f(\mathbf{x}, \mathbf{y}, \mathbf{z}) = \sum_{i=1}^{15}w_xX_i(x_i,y_i,z_i) + w_yY_i(x_i,y_i,z_i)+w_zZ_i(x_i,y_i,z_i)\,,$$

but the constraints are

$$\sum_{i=1}^{15} x_i \le c_x \quad\text{and}\quad x_i\geq 0\text{ for }i=1,\ldots,15\,,$$

$$\sum_{i=1}^{15} y_i \le c_y \quad\text{and}\quad y_i\geq 0\text{ for }i=1,\ldots,15\,,$$

$$\sum_{i=1}^{15} z_i \le c_z \quad\text{and}\quad z_i\geq 0\text{ for }i=1,\ldots,15\,.$$

Note: the sum constraints are inequalities because you don't need to use all the seed. However, the optimal solution will likely use all the available seed.
"""

# ╔═╡ 625ffadc-ff48-11ea-1849-39df2c3ddb85
md"""
## An example

Let us compute the yield for plot 1 when we sow $x_1=10,y_1=6, z_1=8$. First, we compute the maximum yields, then the resulting concrete yields and, finally the value.
"""

# ╔═╡ 57dc25e6-ff38-11ea-342c-f19e90a2c732
x1, y1, z1 = 10, 6, 8  # amount of seed per plot

# ╔═╡ 6137fea8-ff38-11ea-34b7-631de98c3ad4
(Ax1, Ay1, Az1) = S * [u[1], v[1]] .+ b

# ╔═╡ 7e33d806-ff38-11ea-2297-512188693a6a
X1 = (Ax1 * x1^2) / (120 + x1^2 + 2y1 - 0.8x1 * z1 + z1^2)

# ╔═╡ 84e7171e-ff38-11ea-07fd-5fb12ff4f039
Y1 = (Ay1 * y1^2) / (30 + 4x1 + y1^2 + 7z1)

# ╔═╡ 85d3cd6e-ff38-11ea-2d7d-539601eec8aa
Z1 = (Az1 * z1^2) / (80 + 0.4x1 + 0.2x1*z1 +0.6y1 + z1^2)

# ╔═╡ 8a574622-ff38-11ea-15dc-1b4bb63497df
revenue_plot1 = wx * X1 + wy * Y1 + wz * Z1  # in EUR

# ╔═╡ 92a481f4-ff39-11ea-1cc0-855d93b3967b
md"""
# Assignments

1. Give the formal optimization problem. Also provide the Lagrangian formulation.
2. Is the optimization problem of the optimal sowing plan concave? (since it is a maximization problem, I mean is minimizing the negative total revenue of a plan convex). You don't have to prove this formally, but you can make visual arguments. For a given field, make a the contour plots for (x,y), (x,z) and (y,z), always setting the third variable to 0.
3. Give a good/optimal sowing plan. You may solve this either using custom code or using a Julia package [Optim.jl](https://github.com/JuliaNLSolvers/Optim.jl), [Convex.jl](https://github.com/jump-dev/Convex.jl), or [JuMP.jl](https://github.com/jump-dev/JuMP.jl). Make plots of your solution and discuss it.
4. Given an optimal sowing plan found in the previous assignment, what is the price per kg for every seed you would be willing to buy, as to be break-even? How would you use this additional marginal quantity of seed. (HINT: you can obtain this from the Lagrangian).
5. Provide three alternative solutions, each in which you only use one of the types of seed, e.g. $\max_\mathbf{x}f(\mathbf{x}, 0, 0)$. Show and discuss the difference.
"""

# ╔═╡ ea2c0d3e-ff43-11ea-14f4-932602b2fecd
md"""
## Submission

Hand in the solved Jupyter notebook by **14 November 2020**. Send your notebook to [me](michiel.stock@ugent.be) both as Pluto notebook file (.jl) and as **PDF or HTML file**. Hand in a CSV file with your best **valid** solution for question 3. In this file, every row is a plot of land and the three columns represent x, y, z, respectively.
"""

# ╔═╡ 64049c58-ff37-11ea-1ae8-f907d5e33509
student_name = missing  # fill in your name(s)

# ╔═╡ 2568b67a-ff37-11ea-37fd-87a3d569fc31
md"""
# Project: Sowing plan

**STMO**

2020-2021

Project by: $student_name
"""

# ╔═╡ ddab3252-07bc-11eb-390a-750512ee8c90
if ismissing(student_name)
	md"fill in your name(s) below"
end

# ╔═╡ 562bc43a-ff43-11ea-0451-d9e7005c56a1
md"""
### 1. Formal description of the problem

COMPLETE
"""

# ╔═╡ 5232b83c-ff43-11ea-36e4-7d08693203dc
md"""
### 2. Concavity

COMPLETE
"""

# ╔═╡ 8cb9cf68-ff43-11ea-399d-df193400c077
md"""
### 3. Making a sowing plan

Here, you have the space to solve the problem. Save your final solution in the variable `solution` and save it in a csv file using `save_solution`. This is done automatically when your name is filled in.
"""

# ╔═╡ 9d2eabae-ff43-11ea-3103-0f22b5c59e6b

# ╔═╡ 9cd521d8-ff43-11ea-174e-71556bcdebac

# ╔═╡ 67882c04-ff3f-11ea-19b5-890bd3bd15a9
solution = fill(1.0, 15, 3)  # valid, but can be improved...

# ╔═╡ 0a220a2a-ff45-11ea-2b35-c70cdb18c0a8
md"Check if the solution is valid."

# ╔═╡ 1411424e-ff45-11ea-2b3e-c75c7d39917a
md"Save the solution."

# ╔═╡ 07575a02-ff45-11ea-2794-63681f41a7bf
md"""
### 4. Price of seed

Compute how much you would pay for each additional unit of seed in your solution.
"""

# ╔═╡ 79df7272-ff44-11ea-01bb-716d0bf52464

# ╔═╡ 79a661b0-ff44-11ea-10ef-fbc61315162b
md"""
### 5. One type of seed

Make three solutions in which you only use one type of seed. Analyse your solutions!
"""

# ╔═╡ 79897d9a-ff44-11ea-269c-bba522b74d66

# ╔═╡ 796e75c2-ff44-11ea-2a85-0b3fbdc675b5

# ╔═╡ 0ba3f5f6-ff46-11ea-0449-29b505430cf9

# ╔═╡ 7ade6b7e-ff44-11ea-130a-bdd36b48b8af
md"## Functions"

# ╔═╡ b96c0828-ff3b-11ea-3ee3-b979608fe705
"""Check if a solution is valid."""
function isvalidsolution(solution::Matrix)
	return size(solution) == (15, 3) && all(solution .≥ 0.0) && all(sum(solution, dims=1) .≤ [cx cy cz])
end

# ╔═╡ 6f98672e-ff3f-11ea-0198-5fbef05027da
isvalidsolution(solution)

# ╔═╡ a9ca6492-ff44-11ea-3b51-3905bf9280cf
"""Save the solution to a file for submission."""
function save_solution(fname, solution::Matrix)
	@assert isvalidsolution(solution) "Oh no, your solution is invalid!"
	CSV.write(fname, DataFrame(solution))
end

# ╔═╡ 928fd032-ff3f-11ea-13a0-07140325ab9f
!ismissing(student_name) && save_solution("solution_$(student_name).csv", solution);

# ╔═╡ Cell order:
# ╟─2568b67a-ff37-11ea-37fd-87a3d569fc31
# ╟─ddab3252-07bc-11eb-390a-750512ee8c90
# ╟─78701408-ff37-11ea-04a8-07805ef3cc1c
# ╠═bea5f54e-ff37-11ea-2c37-8707c75a9b5c
# ╟─9b5ec9c8-ff37-11ea-1816-0365aa76c976
# ╟─d327e0c2-ff37-11ea-3006-53be838b7dc8
# ╠═c9897d6a-ff47-11ea-1212-9d5eb14cd0d7
# ╟─e4a29d80-ff37-11ea-0d0b-59545bce9bd2
# ╟─ea977544-ff37-11ea-310b-79307ef120d3
# ╠═eba6a9dc-ff37-11ea-233d-1b9aca686ea4
# ╠═0e9a6d52-ff38-11ea-2472-810a521cb8f8
# ╟─100a63f4-ff38-11ea-1fbf-05c136011e28
# ╠═17088b36-ff38-11ea-3519-85055dd00929
# ╟─25b72854-ff38-11ea-0817-0f44e6420882
# ╠═223cec90-ff38-11ea-040f-47c9ee1c9036
# ╟─3282d3c6-ff38-11ea-0458-af74eb6b42dc
# ╟─625ffadc-ff48-11ea-1849-39df2c3ddb85
# ╠═57dc25e6-ff38-11ea-342c-f19e90a2c732
# ╠═6137fea8-ff38-11ea-34b7-631de98c3ad4
# ╠═7e33d806-ff38-11ea-2297-512188693a6a
# ╠═84e7171e-ff38-11ea-07fd-5fb12ff4f039
# ╠═85d3cd6e-ff38-11ea-2d7d-539601eec8aa
# ╠═8a574622-ff38-11ea-15dc-1b4bb63497df
# ╟─92a481f4-ff39-11ea-1cc0-855d93b3967b
# ╟─ea2c0d3e-ff43-11ea-14f4-932602b2fecd
# ╠═64049c58-ff37-11ea-1ae8-f907d5e33509
# ╠═562bc43a-ff43-11ea-0451-d9e7005c56a1
# ╠═5232b83c-ff43-11ea-36e4-7d08693203dc
# ╟─8cb9cf68-ff43-11ea-399d-df193400c077
# ╠═9d2eabae-ff43-11ea-3103-0f22b5c59e6b
# ╠═9cd521d8-ff43-11ea-174e-71556bcdebac
# ╠═67882c04-ff3f-11ea-19b5-890bd3bd15a9
# ╟─0a220a2a-ff45-11ea-2b35-c70cdb18c0a8
# ╠═6f98672e-ff3f-11ea-0198-5fbef05027da
# ╟─1411424e-ff45-11ea-2b3e-c75c7d39917a
# ╠═928fd032-ff3f-11ea-13a0-07140325ab9f
# ╟─07575a02-ff45-11ea-2794-63681f41a7bf
# ╠═79df7272-ff44-11ea-01bb-716d0bf52464
# ╟─79a661b0-ff44-11ea-10ef-fbc61315162b
# ╠═79897d9a-ff44-11ea-269c-bba522b74d66
# ╠═796e75c2-ff44-11ea-2a85-0b3fbdc675b5
# ╠═0ba3f5f6-ff46-11ea-0449-29b505430cf9
# ╟─7ade6b7e-ff44-11ea-130a-bdd36b48b8af
# ╠═2cf61e2a-ff3f-11ea-2d5b-77a2b52fd6e7
# ╠═b96c0828-ff3b-11ea-3ee3-b979608fe705
# ╠═a9ca6492-ff44-11ea-3b51-3905bf9280cf

max_stars_repo_name: tjvneste/STMO

max_issues_repo_name: tjvneste/STMO

max_forks_repo_name: tjvneste/STMO


text:
module QuadraticAssignmentProblem

include("solvers/exact.jl")
include("solvers/lp.jl")
include("solvers/faq.jl")
include("solvers/mms.jl")
include("solvers/ns.jl")
include("solvers/random.jl")

include("util/generators.jl")
include("util/graphmatching.jl")
include("util/metrics.jl")
include("util/util.jl")

export qap_exact
export qap_linearization, qap_lprounding
export faq
export qap_mms
export qap_ns, qap_starpacking, qap_densemapping
export qap_random

export uniform_matrix
export metric_matrix
export zeroone_matrix
export generate_qap

export qap_objective

export qap_set_optimizer!

end


text:
import Test
import JSON

include("../tools/travis/travis-rhf.jl")
#include("s22_gamess_values.jl")

#== select input files ==#
directory = joinpath(@__DIR__, "../example_inputs/S22/")
inputs = readdir(directory)
inputs .= directory .* inputs

display(inputs)

#== initialize JuliaChem ==#
JuliaChem.initialize()

#== read in GAMESS values for comparison ==#
S22_GAMESS_file = open(joinpath(@__DIR__, "s22_gamess_values.json"))
  S22_GAMESS_string = read(S22_GAMESS_file, String)
close(S22_GAMESS_file)

S22_GAMESS = JSON.parse(S22_GAMESS_string)

#== run S22 calculations ==#
molecules = collect(2:3)

s22_test_results = Dict([])
for imol in molecules
  s22_test_results[imol] = travis_rhf(inputs[imol])
end

#== check energies ==#
Test.@testset "S22 Energy" begin
  for imol in molecules
    Test.@test s22_test_results[imol][:Energy]["Energy"] ≈ S22_GAMESS["$imol"]["Energy"]
  end
end

#== check dipole moments ==#
Test.@testset "S22 Dipoles" begin
  for imol in molecules
    if S22_GAMESS["$imol"]["Dipole"] == 1.0E-6
      Test.@test abs(s22_test_results[imol][:Properties]["Dipole"][:moment]) <= S22_GAMESS["$imol"]["Dipole"] #check if approximately zero
    else
      Test.@test s22_test_results[imol][:Properties]["Dipole"][:moment] ≈ S22_GAMESS["$imol"]["Dipole"] atol=5.0E-5
    end
  end
end

#== check HOMO-LUMO gaps ==#
Test.@testset "S22 HOMO-LUMO Gaps" begin
  for imol in molecules
    Test.@test s22_test_results[imol][:Properties]["MO Energies"][:homo_lumo] ≈ S22_GAMESS["$imol"]["HOMO-LUMO Gap"] atol=5.0E-4
  end
end

#== check Mulliken charges ==#
Test.@testset "S22 Mulliken Charges" begin
  for imol in molecules
    Test.@test s22_test_results[imol][:Properties]["Mulliken Population"] ≈
      S22_GAMESS["$imol"]["Mulliken Population"] atol=5.0E-6
  end
end

#== finalize JuliaChem ==#
JuliaChem.finalize()


text:
using LinearAlgebra
using QuasiTriangular
using Random
using Test

Aorig = [1 3; 0 2]
A1 = [Aorig Aorig; zeros(2,2) Aorig]
A = QuasiUpperTriangular(A1)
B = randn(4,4)
C = zeros(4,4)
#A_mul_B!(C, 1, B, 1, 4, 4, A, 1, 4)

#@test C ≈ B*A1

Aorig = [1 3; 0.5 2]
A = QuasiUpperTriangular(Aorig)
B = Matrix{Float64}(I(2))
X = zeros(2,2)
#A_ldiv_B!(A,B)

#@test B ≈ inv(Aorig)

B = Matrix{Float64}(I(2))
A_rdiv_Bt!(B,A)
@test B ≈ inv(Aorig)'

B = Matrix{Float64}(I(2))
A_rdiv_B!(B,A)

@test B ≈ inv(Aorig)

Random.seed!(123)
n = 7
a = randn(n,n)
S = schur(a)
t = S.T
b = randn(n,n)
c = similar(b)

@test t*b ≈ A_mul_B!(QuasiUpperTriangular(t),b)
@test t'*b ≈ At_mul_B!(QuasiUpperTriangular(t),b)
@test b*t ≈ A_mul_B!(b,QuasiUpperTriangular(t))
@test b*t' ≈ A_mul_Bt!(b,QuasiUpperTriangular(t))

@test t*b ≈ A_mul_B!(c,QuasiUpperTriangular(t),b)
@test t'*b ≈ At_mul_B!(c,QuasiUpperTriangular(t),b)
@test b*t ≈ A_mul_B!(c,b,QuasiUpperTriangular(t))
@test b*t' ≈ A_mul_Bt!(c,b,QuasiUpperTriangular(t))

b1 = copy(b)
x = zeros(n,n)
A_ldiv_B!(QuasiUpperTriangular(t),b1)
@test t\b ≈ b1
b1 = copy(b)
A_rdiv_B!(b1,QuasiUpperTriangular(t))
@test b/t ≈ b1
b1 = copy(b)
A_rdiv_Bt!(b1,QuasiUpperTriangular(t))
@test b/t' ≈ b1

b = rand(n)
b1 = copy(b)
r = rand()
I_plus_rA_ldiv_B!(r,QuasiUpperTriangular(t),b1)
@test b1 ≈ (Matrix{Float64}(I(n)) + r*t)\b
b1 = copy(b)
s = rand()
I_plus_rA_plus_sB_ldiv_C!(r,s,QuasiUpperTriangular(t),QuasiUpperTriangular(t*t),b1)
@test b1 ≈ (Matrix{Float64}(I(n)) + r*t + s*t*t)\b

b = rand(n,n)
b1 = copy(b)
r = rand()
I_plus_rA_ldiv_B!(r,QuasiUpperTriangular(t),b1)
@test b1 ≈ (Matrix{Float64}(I(n)) + r*t)\b
b1 = copy(b)
s = rand()
I_plus_rA_plus_sB_ldiv_C!(r,s,QuasiUpperTriangular(t),QuasiUpperTriangular(t*t),b1)
@test b1 ≈ (Matrix{Float64}(I(n)) + r*t + s*t*t)\b


text:
# load DQMC framework
include(joinpath(ENV["JULIA_DQMC"], "src/dqmc_framework.jl"))

using HDF5
c = h5read("conf_near_qcp_L_10_B_40.h5", "conf")

p = Params()
p.output_file = "live.out.h5.running"
xml2parameters!(p, "dqmc.in.xml")
mc = DQMC(p)

h5open("scipost_sdw.h5", isfile("scipost_sdw.h5") ? "r+" : "w") do file
    for β in (40)#range(5,40,step=5)
        mc.p.beta = β
        mc.p.slices = Int(mc.p.beta/mc.p.delta_tau)
        init!(mc)
        p.hsfield .= c
        for l in 1:mc.p.slices
            file["B_$β/$l"] = slice_matrix(mc,l)
        end
    end
end