text:
cd(@__DIR__)
using Pkg; Pkg.activate("."); Pkg.instantiate()

#This script simulates the Fisher-KPP equation and fits
#a neural PDE to the data with the growth (aka reaction) term replaced
#by a feed-forward neural network and the diffusion term with a CNN

using PyPlot, Printf
using LinearAlgebra
using Flux, DiffEqFlux, Optim, DiffEqSensitivity
using BSON: @save, @load
using Flux: @epochs
using OrdinaryDiffEq

#parameter
D = 0.01; #diffusion
r = 1.0; #reaction rate

#domain
X = 1.0; T = 5;
dx = 0.04; dt = T/10;
x = collect(0:dx:X);
t = collect(0:dt:T);
Nx = Int64(X/dx+1);
Nt = Int64(T/dt+1);

#initial conditions
Amp = 1.0;
Delta = 0.2
#IC-1
rho0 = Amp*(tanh.((x .- (0.5 - Delta/2))/(Delta/10)) - tanh.((x .- (0.5 + Delta/2))/(Delta/10)))/2
#IC-2
#rho0 = Amp*(1 .- tanh.((x .- 0.2)/(Delta/6)))/2.

save_folder = "data"

if isdir(save_folder)
    rm(save_folder, recursive=true)
end
mkdir(save_folder)

close("all")
figure()
plot(x, rho0)
title("Initial Condition")
gcf()

########################
# Generate training data
########################
reaction(u) = r * u .* (1 .- u)
lap = diagm(0 => -2.0 * ones(Nx), 1=> ones(Nx-1), -1 => ones(Nx-1)) ./ dx^2
#Periodic BC
lap[1,end] = 1.0/dx^2
lap[end,1] = 1.0/dx^2
#Neumann BC
#lap[1,2] = 2.0/dx^2
#lap[end,end-1] = 2.0/dx^2

function rc_ode(rho, p, t)
    #finite difference
    D * lap * rho + reaction.(rho)
end

prob = ODEProblem(rc_ode, rho0, (0.0, T), saveat=dt)
sol = solve(prob, Tsit5());
ode_data = Array(sol);

figure(figsize=(8,3))

subplot(121)
pcolor(x,t,ode_data')
xlabel("x"); ylabel("t");
colorbar()

subplot(122)
for i in 1:2:Nt
    plot(x, ode_data[:,i], label="t=$(sol.t[i])")
end
xlabel("x"); ylabel(L"$\rho$")
legend(frameon=false, fontsize=7, bbox_to_anchor=(1, 1), loc="upper left", ncol=1)
tight_layout()
savefig(@sprintf("%s/training_data.pdf", save_folder))
gcf()

########################
# Define the neural PDE
########################
n_weights = 3

A = [FourierBasis(n_weights)]
rx_nn = TensorLayer(A, 1)

#conv with bias with initial values as 1/dx^2
w_err = 0.0
init_w = reshape([1.1 -2.5 1.0], (3, 1, 1, 1))
diff_cnn_ = Conv(init_w, [0.], pad=(0,0,0,0))

#initialize D0 close to D/dx^2
D0 = [6.5]

p1 = 0.01randn(n_weights)
p2,re2 = Flux.destructure(diff_cnn_)
p = [p1;p2;D0]

function nn_ode(u,p,t)
    u_cnn_1   = [p[end-4] * u[end] + p[end-3] * u[1] + p[end-2] * u[2]]
    u_cnn     = [p[end-4] * u[i-1] + p[end-3] * u[i] + p[end-2] * u[i+1] for i in 2:Nx-1]
    u_cnn_end = [p[end-4] * u[end-1] + p[end-3] * u[end] + p[end-2] * u[1]]

    # Equivalent using Flux, but slower!
    #CNN term with periodic BC
    #diff_cnn_ = Conv(reshape(p[(end-4):(end-2)],(3,1,1,1)), [0.0], pad=(0,0,0,0))
    #u_cnn = reshape(diff_cnn_(reshape(u, (Nx, 1, 1, 1))), (Nx-2,))
    #u_cnn_1 = reshape(diff_cnn_(reshape(vcat(u[end:end], u[1:1], u[2:2]), (3, 1, 1, 1))), (1,))
    #u_cnn_end = reshape(diff_cnn_(reshape(vcat(u[end-1:end-1], u[end:end], u[1:1]), (3, 1, 1, 1))), (1,))

    [rx_nn([u[i]/π],p[1:length(p1)])[1] for i in 1:Nx] + p[end] * vcat(u_cnn_1, u_cnn, u_cnn_end)
end

########################
# Soving the neural PDE and setting up loss function
########################
prob_nn = ODEProblem(nn_ode, rho0, (0.0, T), p)
sol_nn = concrete_solve(prob_nn,Tsit5(), rho0, p)

function predict_rd(θ)
  # No ReverseDiff if using Flux
  Array(concrete_solve(prob_nn,Tsit5(),rho0,θ,saveat=dt,sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP())))
end

#match data and force the weights of the CNN to add up to zero
function loss_rd(p)
    pred = predict_rd(p)
    sum(abs2, ode_data .- pred) + 10^2 * abs(sum(p[end-4 : end-2])), pred
end

########################
# Training
########################

#Optimizer
opt = ADAM(0.001)

global count = 0
global save_count = 0
save_freq = 50

train_arr = Float64[]
diff_arr = Float64[]
w1_arr = Float64[]
w2_arr = Float64[]
w3_arr = Float64[]

#callback function to observe training
cb = function (p,l,pred)
    rx_nn_p, diff_cnn_, D0 = p[1:length(p1)],re2(p[(length(p1)+1):(end-1)]),p[end]
    push!(train_arr, l)
    push!(diff_arr, p[end])

    weight = diff_cnn_.weight[:]
    push!(w1_arr, weight[1])
    push!(w2_arr, weight[2])
    push!(w3_arr, weight[3])

    println(@sprintf("Loss: %0.4f\tD0: %0.4f Weights:(%0.4f,\t %0.4f, \t%0.4f) \t Sum: %0.4f"
            ,l, D0[1], weight[1], weight[2], weight[3], sum(weight)))

    global count

    if count==0
        fig = figure(figsize=(8,2.5));
        ttl = fig.suptitle(@sprintf("Epoch = %d", count), y=1.05)
        global ttl
        subplot(131)
        pcolormesh(x,t,ode_data')
        xlabel(L"$x$"); ylabel(L"$t$"); title("Data")
        colorbar()

        subplot(132)
        img = pcolormesh(x,t,pred')
        global img
        xlabel(L"$x$"); ylabel(L"$t$"); title("Prediction")
        colorbar(); clim([0, 1]);

        ax = subplot(133); global ax
        u = collect(0:0.01:1)
        rx_line = plot(u, rx_nn.([[elem/π] for elem in u],(rx_nn_p,)), label="NN")[1];
        global rx_line
        plot(u, reaction.(u), label="True")
        title("Reaction Term")
        legend(loc="upper right", frameon=false, fontsize=8);
        ylim([0, r*0.25+0.2])

        subplots_adjust(top=0.8)
        tight_layout()
    end

    if count>0
        println("updating figure")
        img.set_array(pred[1:end-1, 1:end-1][:])
        ttl.set_text(@sprintf("Epoch = %d", count))

        u = collect(0:0.01:1)
        rx_pred = reduce(vcat,rx_nn.([[elem/π] for elem in u],(rx_nn_p,)))
        rx_line.set_ydata(rx_pred)
        u = collect(0:0.01:1)

        min_lim = min(minimum(rx_pred), minimum(reaction.(u)))-0.1
        max_lim = max(maximum(rx_pred), maximum(reaction.(u)))+0.1

        ax.set_ylim([min_lim, max_lim])
    end

    global save_count
    if count%save_freq == 0
        println("saved figure")
        savefig(@sprintf("%s/pred_%05d.png", save_folder, save_count), dpi=200, bbox_inches="tight")
        save_count += 1
    end

    display(gcf())
    count += 1

    l < 0.01 # Exit when fit to 2 decimal places
end

#train
@time begin
    res1 = DiffEqFlux.sciml_train(loss_rd, p, ADAM(0.001), cb=cb, maxiters = 100)
    res2 = DiffEqFlux.sciml_train(loss_rd, res1.minimizer, ADAM(0.001), cb=cb, maxiters = 300)
    res3 = DiffEqFlux.sciml_train(loss_rd, res2.minimizer, BFGS(initial_stepnorm=0.00001), cb=cb, maxiters = 1000, allow_f_increases=true)
end

pstar = res3.minimizer

## Save trained model
@save @sprintf("%s/model.bson", save_folder) pstar

########################
# Plot for paper
########################
@load @sprintf("%s/model.bson", save_folder) pstar
#re-defintions for newly loaded data

diff_cnn_ = Conv(reshape(pstar[(end-4):(end-2)],(3,1,1,1)), [0.0], pad=(0,0,0,0))
diff_cnn(x) = diff_cnn_(x) .- diff_cnn_.bias
D0 = res3.minimizer[end]

fig = figure(figsize=(4,4))

rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["font.size"] = 10
rcParams["text.usetex"] = true
rcParams["font.family"] = "serif"
rcParams["font.sans-serif"] = "Helvetica"
rcParams["axes.titlesize"] = 10

subplot(221)
pcolormesh(x,t,ode_data', rasterized=true)
xlabel(L"$x$"); ylabel(L"$t$"); title("Data")
yticks([0, 1, 2, 3, 4, 5])

ax = subplot(222)
cur_pred = predict_rd(pstar)[1]
img = pcolormesh(x,t,cur_pred', rasterized=true)
global img
xlabel(L"$x$"); ylabel(L"$t$"); title("Prediction")
yticks([0, 1, 2, 3, 4, 5])
cax = fig.add_axes([.48,.62,.02,.29])
colb = fig.colorbar(img, cax=cax)
colb.ax.set_title(L"$\rho$")
clim([0, 1]);
colb.set_ticks([0, 1])

subplot(223)
plot(Flux.data(w1_arr ./ w3_arr) .- 1, label=L"$w_1/w_3 - 1$")
plot(Flux.data(w1_arr .+ w2_arr .+ w3_arr), label=L"$w_1 + w_2 + w_3$")
axhline(0.0, linestyle="--", color="k")
xlabel("Epochs"); title("CNN Weights")
xticks([0, 1500, 3000]); yticks([-0.4, -0.3,-0.2, -0.1, 0.0, 0.1])
legend(loc="lower right", frameon=false, fontsize=6)

subplot(224)
u = collect(0:0.01:1)
plot(u, rx_nn.([[elem] for elem in u]), label="UPDE")[1];
plot(u, reaction.(u), linestyle="--", label="True")
xlabel(L"$\rho$")
title("Reaction Term")
legend(loc="lower center", frameon=false, fontsize=6);
ylim([0, 0.3])

tight_layout(h_pad=1)
gcf()
savefig(@sprintf("%s/fisher_kpp.pdf", save_folder))

#plot loss vs epochs and save
figure(figsize=(6,3))
plot(log.(train_arr), "k.", markersize=1)
xlabel("Epochs"); ylabel("Log(loss)")
tight_layout()
savefig(@sprintf("%s/loss_vs_epoch.pdf", save_folder))
gcf()

#=

# 3 parameters

Loss: 0.0094    D0: 5.7811 Weights:(1.1501,      -2.3004,       1.1503)         Sum: -0.0000
updating figure
237.483188 seconds (875.53 M allocations: 44.014 GiB, 3.10% gc time)

Loss: 0.0095    D0: 5.7752 Weights:(1.1496,      -2.2988,       1.1492)         Sum: 0.0000
updating figure
saved figure
240.444140 seconds (885.66 M allocations: 44.519 GiB, 3.11% gc time)

Loss: 0.0086    D0: 5.7983 Weights:(1.1519,      -2.3035,       1.1517)         Sum: -0.0000
updating figure
234.218061 seconds (917.69 M allocations: 46.116 GiB, 3.20% gc time)

Loss: 0.0090    D0: 5.1810 Weights:(1.2899,      -2.5798,       1.2899)         Sum: 0.0000
updating figure
232.159919 seconds (930.76 M allocations: 46.780 GiB, 3.25% gc time)

Loss: 0.0095    D0: 5.7827 Weights:(1.1524,      -2.3052,       1.1527)         Sum: -0.0000
updating figure
239.525352 seconds (953.27 M allocations: 47.905 GiB, 3.23% gc time)

# 5 parameters

Loss: 0.0016    D0: 5.8669 Weights:(1.0638,      -2.1274,       1.0636)         Sum: 0.0000
updating figure
243.928965 seconds (1.00 G allocations: 50.893 GiB, 3.40% gc time)

Loss: 0.0069    D0: 5.8545 Weights:(1.0472,      -2.0954,       1.0482)         Sum: -0.0000
updating figure
244.803636 seconds (1.01 G allocations: 51.578 GiB, 3.51% gc time)

Loss: 0.0069    D0: 5.8545 Weights:(1.0472,      -2.0954,       1.0482)         Sum: -0.0000
updating figure
240.307752 seconds (1.00 G allocations: 51.091 GiB, 3.43% gc time)

Loss: 0.0085    D0: 5.8546 Weights:(1.0782,      -2.1505,       1.0723)         Sum: 0.0000
updating figure
265.215108 seconds (1.14 G allocations: 57.873 GiB, 3.59% gc time)

Loss: 0.0039    D0: 5.8330 Weights:(1.0403,      -2.0784,       1.0382)         Sum: -0.0000
updating figure
246.503862 seconds (1.07 G allocations: 54.304 GiB, 3.59% gc time)

# 7 parameters

Loss: 0.0012    D0: 5.9098 Weights:(1.0402,      -2.0810,       1.0408)         Sum: -0.0000
updating figure
242.417536 seconds (1.14 G allocations: 56.570 GiB, 3.71% gc time)

Loss: 0.0055    D0: 5.9350 Weights:(1.0404,      -2.0825,       1.0421)         Sum: -0.0000
updating figure
252.128401 seconds (1.16 G allocations: 57.507 GiB, 3.88% gc time)

Loss: 0.0086    D0: 5.9322 Weights:(1.0501,      -2.1005,       1.0503)         Sum: -0.0001
updating figure
249.834412 seconds (1.14 G allocations: 56.479 GiB, 3.71% gc time)

Loss: 0.0083    D0: 5.8940 Weights:(1.0432,      -2.0869,       1.0437)         Sum: 0.0001
updating figure
255.021489 seconds (1.21 G allocations: 59.928 GiB, 3.92% gc time)

Loss: 0.0052    D0: 5.9229 Weights:(1.0481,      -2.0961,       1.0481)         Sum: 0.0000
updating figure
253.374001 seconds (1.23 G allocations: 60.636 GiB, 3.71% gc time)

# 15 parameters

Loss: 0.0032    D0: 6.1526 Weights:(1.0477,      -2.0939,       1.0462)         Sum: -0.0000
updating figure
267.148464 seconds (1.72 G allocations: 83.585 GiB, 4.00% gc time)

Loss: 0.0057    D0: 6.1925 Weights:(1.0960,      -2.1901,       1.0942)         Sum: 0.0000
updating figure
281.902850 seconds (1.91 G allocations: 92.720 GiB, 4.15% gc time)

Loss: 0.0043    D0: 6.1548 Weights:(1.0607,      -2.1218,       1.0611)         Sum: 0.0000
updating figure
262.655202 seconds (1.76 G allocations: 85.483 GiB, 4.01% gc time)

Loss: 0.0088    D0: 6.1617 Weights:(1.0751,      -2.1503,       1.0751)         Sum: -0.0001
updating figure
267.427766 seconds (1.80 G allocations: 87.616 GiB, 3.99% gc time)

Loss: 0.0095    D0: 6.1630 Weights:(1.0691,      -2.1416,       1.0724)         Sum: -0.0001
updating figure
269.782135 seconds (1.80 G allocations: 87.805 GiB, 4.00% gc time)
=#

x = [242.417536,252.128401,249.834412,255.021489,253.374001]
mean(x)
using Statistics
std(x)


text:
# following code is used to caclulate moments in SymmetricTensor form ##
"""
    blockel(X::Matrix{T}, i::Tuple, j::Tuple, b::Int)

Returns Float, the element of the block (indexed by j) of the moment's tensor
of X, at index cd i inside a block, where b is a standard blocks' size

```jldoctest
julia> M = [1. 2.  5. 6. ; 3. 4.  7. 8.];

julia> mom_el(M, (1,1), (1,1), 2)
5.0

julia> mom_el(M, (1,1), (2,2), 2)
37.0
```
"""
function blockel(data::Matrix{T}, mi::Tuple, mj::Tuple, b::Int) where T <: AbstractFloat
  ret = 0.
  t = size(data, 1)
  for l in 1:t
    temp = 1.
    for k in 1:length(mi)
      @inbounds ind = (mj[k]-1)*b+mi[k]
      @inbounds temp *= data[l,ind]
    end
    ret += temp
  end
  ret/t
end

"""
  momentblock(X::Matrix{T}, j::Tuple, dims::Tuple, b::Int)

Returns a block of a moment's tensor of X. A block is indexed by j and if size dims,
b is a standatd block size.

```jldoctest
julia> M = [1. 2.  5. 6. ; 3. 4.  7. 8.];

julia> momentblock(M, (1,1), (2,2), 2)
2×2 Array{Float64,2}:
 5.0   7.0
 7.0  10.0
```
"""
function momentblock(X::Matrix{T}, j::Tuple, dims::Tuple,
                                             b::Int) where {T <: AbstractFloat}
  ret = zeros(T, dims)
  for ind = 1:(prod(dims))
    i = Tuple(CartesianIndices(dims)[ind])
    @inbounds ret[i...] = blockel(X, i, j, b)
  end
  ret
end

"""
    usebl(bind::Tuple, n::Int, b::Int, nbar::Int)

Returns: Tuple{Int}, sizes of the last block
"""
function usebl(bind::Tuple, n::Int, b::Int, nbar::Int)
  bl = n - b*(nbar-1)
  map(i -> (i == nbar) ? (bl) : (b), bind)
end

"""
    momentn1c(X::Matrix{Float}, m::Int, b::Int)

Returns: SymmetricTensor{Float, m}, a tensor of the m'th moment of X, where b
is a block size. Uses 1 core implementation
"""
function moment1c(X::Matrix{T}, m::Int, b::Int=2) where T <: AbstractFloat
  n = size(X, 2)
  sizetest(n, b)
  nbar = mod(n,b)==0 ? n÷b : n÷b + 1
  ret = arraynarrays(T, fill(nbar, m)...,)
  for j in pyramidindices(m, nbar)
    dims = (mod(n,b) == 0 || !(nbar in j)) ? (fill(b,m)...,) : usebl(j, n, b, nbar)
    @inbounds ret[j...] = momentblock(X, j, dims, b)
  end
  SymmetricTensor(ret; testdatstruct = false)
end

"""
    momentnc(X::Matrix}, m::Int, b::Int)

Returns: SymmetricTensor{Float, m}, a tensor of the m'th moment of X, where b
is a block size. Uses multicore parallel implementation via pmap()
"""
function momentnc(x::Matrix{T}, m::Int, b::Int = 2) where T <: AbstractFloat
  t = size(x, 1)
  f(z::Matrix{T}) = moment1c(z, m, b)
  k = length(workers())
  r = mod(t,k)==0 ? t÷k : t÷k + 1
  y = [x[ind2range(i, r, t), :] for i in 1:k]
  ret = pmap(f, y)
  (r*sum(ret[1:(end-1)])+(t-(k-1)*r)*ret[end])/t
end

"""
    moment(X::Matrix}, m::Int, b::Int)

Returns: SymmetricTensor{Float, m}, a tensor of the m'th moment of X, where b
is a block size. Calls 1 core or multicore moment function.
"""
moment(X::Matrix{T}, m::Int, b::Int=2) where T <: AbstractFloat =
  (size(X,1)/10>nworkers()>1) ? momentnc(X, m, b) : moment1c(X, m, b)

# ---- following code is used to caclulate cumulants in SymmetricTensor form----
"""

Type that stores a partition of multiindex into subests, sizes of subests,
size of original multitindex and number of subsets

"""
mutable struct IndexPart
  part::Vector{Vector{Int64}}
  subsetslen::Vector{Int64}
  nind::Int
  npart::Int
  (::Type{IndexPart})(part::Vector{Vector{Int64}}, subsetslen::Vector{Int64},
nind::Int, npart::Int) = new(part, subsetslen, nind, npart)
end

"""
    indpart(nind::Int, npart::Int, e::Int = 1)

Returns vector of IndexPart type, that includes partitions of set [1, 2, ..., nind]
into npart subests of size != e, sizes of each subest, size of original set and
number of partitions

```jldoctest
julia>indpart(4,2)
3-element Array{Cumulants.IndexPart,1}:
 IndexPart(Array{Int64,1}[[1,2],[3,4]],[2,2],4,2)
 IndexPart(Array{Int64,1}[[1,3],[2,4]],[2,2],4,2)
 IndexPart(Array{Int64,1}[[1,4],[2,3]],[2,2],4,2)
```
"""
function indpart(nind::Int, npart::Int, e::Int = 1)
    part_set = IndexPart[]
    for part in partitions(1:nind, npart)
      subsetslen = map(length, part)
      if !(e in subsetslen)
        push!(part_set, IndexPart(part, subsetslen, nind, npart))
      end
    end
    part_set
end

"""
  accesscum(mulind::Tuple{Int, ...}, ::IndexPart,
    cum::SymmetricTensor{Float}...)

Returns: vector of blocks from cumulants. Each block correspond to a subests
 of partition (part) of multiindex (multiind).

 ```jldoctest
 julia> cum = SymmetricTensor([1.0 2.0 3.0; 2.0 4.0 6.0; 3.0 6.0 5.0]);

julia> accesscum((1,1,1,1), IndexPart(Array{Int64,1}[[1,2],[3,4]],[2,2],4,2), cum)
Array{Float64,N}[
[1.0 2.0; 2.0 4.0],
[1.0 2.0; 2.0 4.0]]

julia> accesscum((1,1,1,2), IndexPart(Array{Int64,1}[[1,2],[3,4]],[2,2],4,2), cum)

Array{Float64,N}[
[1.0 2.0; 2.0 4.0],
[3.0 0.0; 6.0 0.0]]

julia> accesscum((1,1,1,1), IndexPart(Array{Int64,1}[[1,4],[2,3]],[2,2],4,2), cum)

Array{Float64,N}[
[1.0 2.0; 2.0 4.0],
[1.0 2.0; 2.0 4.0]]
 ```
"""
function accesscum(mulind::Tuple, part::IndexPart,
                                  cum::SymmetricTensor{T}...) where T <: AbstractFloat
  blocks = Array{Array{T}}(undef, part.npart)
  sq = cum[1].sqr || !(cum[1].bln in mulind)
  for k in 1:part.npart
    data = getblockunsafe(cum[part.subsetslen[k]], mulind[part.part[k]])
    if sq
      @inbounds blocks[k] = data
    else
      ind = map(i -> 1:size(data,i), 1:part.subsetslen[k])
      datapadded = zeros(T, fill(cum[1].bls, part.subsetslen[k])...,)
      @inbounds datapadded[ind...] = data
      @inbounds blocks[k] = datapadded
    end
  end
  blocks
end

"""

    outprodblocks(n::Int, part::Vector{Vector{Int}}, blocks::Vector{Array{T}}

Returns: n dims Array of outer product of blocks, given partition of indices, part.

```jldoctest
julia> blocks = 2-element Array{Array{Float64,N},1}[[1.0 2.0; 2.0 4.0], [1.0 2.0; 2.0 4.0]];

julia> outprodblocks(IndexPart(Array{Int64,1}[[1,2],[3,4]],[2,2],4,2), blocks)

2×2×2×2 Array{Float64,4}:
[:, :, 1, 1] =
 1.0  2.0
 2.0  4.0

[:, :, 2, 1] =
 2.0  4.0
 4.0  8.0

[:, :, 1, 2] =
 2.0  4.0
 4.0  8.0

[:, :, 2, 2] =
 4.0   8.0
 8.0  16.0
```
"""
function outprodblocks(inp::IndexPart,
                       blocks::Vector{Array{T}}) where T <: AbstractFloat
  b = size(blocks[1], 1)
  block = zeros(T, fill(b, inp.nind)...,)
  for i = 1:(b^inp.nind)
    muli = Tuple(CartesianIndices((fill(b, inp.nind)...,))[i])
    @inbounds block[muli...] =
    mapreduce(k -> blocks[k][muli[inp.part[k]]...], *, 1:inp.npart)
  end
  block
end

"""
    outerprodcum(retd::Int, npart::Int, cum::SymmetricTensor...; exclpartlen::Int = 1)

Returns retd dims outer products of npart cumulants in SymmetricTensor form.
exclpartlen is a length of partitions to be excluded in calculations,
in this algorithm exclpartlen = 1

```jldoctest
julia> cum = SymmetricTensor([1.0 2.0 3.0; 2.0 4.0 6.0; 3.0 6.0 5.0]);

julia> outerprodcum(4,2,cum, cum)
SymmetricTensors.SymmetricTensor{Float64,4}(Union{Array{Float64,4}, Void}[[3.0 6.0; 6.0 12.0]

[6.0 12.0; 12.0 24.0]

[6.0 12.0; 12.0 24.0]

[12.0 24.0; 24.0 48.0] nothing; nothing nothing]

Union{Array{Float64,4}, Void}[nothing nothing; nothing nothing]

Union{Array{Float64,4}, Void}[[9.0 18.0; 18.0 36.0]

[18.0 36.0; 36.0 72.0] nothing; nothing nothing]

Union{Array{Float64,4}, Void}[[23.0 46.0; 46.0 92.0] [45.0; 90.0]; nothing [75.0]], 2, 2, 3, false)
```
"""
function outerprodcum(retd::Int, npart::Int,
                                 cum::SymmetricTensor{T}...;
                                 exclpartlen::Int = 1) where T <: AbstractFloat
  parts = indpart(retd, npart, exclpartlen)
  prodcum = arraynarrays(T, fill(cum[1].bln, retd)...,)
  for muli in pyramidindices(retd, cum[1].bln)
    block = zeros(T, fill(cum[1].bls, retd)...,)
    for part in parts
      blocks = accesscum(muli, part, cum...)
      @inbounds block += outprodblocks(part, blocks)
    end
    if !cum[1].sqr && cum[1].bln in muli
      ran = map(k->cum[1].bln == muli[k] ? (1:cum[1].dats% cum[1].bls) : (1:cum[1].bls), 1:retd)
      @inbounds block = block[ran...]
    end
    @inbounds prodcum[muli...] = block
  end
  SymmetricTensor(prodcum; testdatstruct = false)
end

"""
    cumulant(X::Vector{Matrix}, cum::SymmetricTensor...)

Returns: SymmetricTensor{Float, m}, a tensor of the m'th cumulant of X, given Vector
of cumulants of order 2, ..., m-2
"""
function cumulant(X::Matrix{T}, cum::SymmetricTensor{T}...) where T <: AbstractFloat
  m = length(cum) + 2
  ret =  moment(X, m, cum[1].bls)
  for sigma in 2:div(m, 2)
    ret -= outerprodcum(m, sigma, cum...)
  end
  ret
end

"""
    cumulants(X::Matrix, m::Int, b::Int)

Returns [SymmetricTensor{Float, 1}, SymmetricTensor{Float, 2}, ...,
SymmetricTensor{Float, m}], vector of cumulant tensors

```
julia> M =  [[-0.88626   0.279571];[-0.704774  0.131896]];

julia> convert(Array, cumulants(M, 3)[3])
2×2×2 Array{Float64,3}:
[:, :, 1] =
 0.0  0.0
 0.0  0.0

[:, :, 2] =
 0.0  0.0
 0.0  0.0
```
"""
function cumulants(X::Matrix{T}, m::Int = 4, b::Int = 2) where T <: AbstractFloat
  cvec = Array{SymmetricTensor{T}}(undef, m)
  cvec[1] = moment1c(X, 1, b)
  X = X .- mean(X, dims=1)
  for i = 2:m
    @inbounds cvec[i] = (i < 4) ? moment1c(X, i, b) : cumulant(X, cvec[1:(i-2)]...)
  end
  cvec
end

ext: jl


text:
"""
    DNN(d = 0,
        metric = Euclidean(),
        algorithm = :kdtree,
        leafsize = 10,
        reorder = true,
        parallel = false)

Anomaly score based on the number of neighbors in a hypersphere of radius `d`. Knorr et al. [1] directly converted the
resulting outlier scores to labels, thus this implementation does not fully reflect the approach from the paper.

Parameters
----------
    d::Real
The hypersphere radius used to calculate the global density of an instance.

$_knn_shared

Examples
--------
$(_score_unsupervised("DNN"))

References
----------
[1] Knorr, Edwin M.; Ng, Raymond T. (1998): Algorithms for Mining Distance-Based Outliers in Large Datasets.
"""
MMI.@mlj_model mutable struct DNN <: UnsupervisedDetector
    metric::DI.Metric = DI.Euclidean()
    algorithm::Symbol = :kdtree::(_ in (:kdtree, :balltree))
    leafsize::Integer = 10::(_ ≥ 0)
    reorder::Bool = true
    parallel::Bool = false
    d::Real = 0::(_ > 0) # warns if `d` is not set
end

struct DNNModel <: Model
    tree::NN.NNTree
end

function fit(detector::DNN, X::Data)::Fit
    # create the specified tree
    tree = buildTree(X, detector.metric, detector.algorithm, detector.leafsize, detector.reorder)

    # use tree to calculate distances
    scores = dnn_others(NN.inrange(tree, X, detector.d))
    Fit(DNNModel(tree), scores)
end

@score function score(detector::DNN, model::Fit, X::Data)::Result
    if detector.parallel
        # already returns scores
        return dnn_parallel(model.tree, X, detector.d)
    else
        return dnn(NN.inrange(model.tree, X, detector.d))
    end
end

@inline function dnn(idxs::AbstractVector{<:AbstractVector})::Scores
    # Helper function to reduce the instances to a global density score.
    1 ./ (length.(idxs) .+ 0.1) # min score = 0, max_score = 10
end

@inline function dnn_others(idxs::AbstractVector{<:AbstractVector})::Scores
    # Remove the (self) point previously added when fitting the tree, otherwise during `fit`, that point would always
    # be included in the density estimation
    1 ./ (length.(idxs) .- 0.9) # 1 - 0.1
end


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# knapsack.jl
#
# Solves a simple knapsack problem:
# max sum(p_j x_j)
#  st sum(w_j x_j) <= C
#     x binary
#############################################################################

using JuMP

# Maximization problem
m = Model()

@variable(m, x[1:5], Bin)

profit = [ 5, 3, 2, 7, 4 ]
weight = [ 2, 8, 4, 2, 5 ]
capacity = 10

# Objective: maximize profit
@objective(m, Max, dot(profit, x))

# Constraint: can carry all
@constraint(m, dot(weight, x) <= capacity)

# Solve problem using MIP solver
status = solve(m)

println("Objective is: ", getobjectivevalue(m))
println("Solution is:")
for i = 1:5
    print("x[$i] = ", getvalue(x[i]))
    println(", p[$i]/w[$i] = ", profit[i]/weight[i])
end


text:
using Test, Astro

println("10.1 Dynamical Time and Universal Time")
# UT is based on the rotation of the earth, but the earth is slowing down and irregular in its orbital habits..
# DynamicalTime is related to Ephemeris time, but changes year to year.
# For example, 1990 January 27, 0h UT is 57 seconds later than 1990 January 27, 0h TD.
jd = cal_to_jd(1990,1,27) + hms_to_fday(0,0,0)
@test isapprox(jd, 2.4479185e6, atol=0.1)
deltat = deltaT_seconds(jd)
@test deltat ≈  57 atol= 1# accurate to a second -— good enough?
println(" passed")

println("10.2 Dynamical Time and Universal Time")
# New Moon took place on 1977 Feburary 18 at 3h37m40s
# deltaT was equal to +48 seconds, correpsonding UT time is 3h36m52s
jd = cal_to_jd(1977,2,18) + hms_to_fday(3,37,40)
deltat = deltaT_seconds(jd)
@test isapprox(deltat, 48, atol=1)
println(" passed")


text:
struct ResourceCounter
    path::String
    id::Cint
    nresources::Cint
    undo::Bool
    tok::Cftok
    handle::Cint
    counts::Vector{Cushort}

    function ResourceCounter(path::Union{String,Tuple{String,Integer}}, nresources::Integer=1; create::Bool=true, create_exclusive::Bool=false, permission::Integer=0o660, undo::Bool=false)
        id = 0
        if isa(path, Tuple)
            path,id = path
        end
        tok = ftok(path, id)
        handle = semcreate(tok, nresources; create=create, create_exclusive=create_exclusive, permission=permission)
        new(path, id, nresources, undo, tok, handle, Vector{Cushort}(undef, nresources))
    end
end

close(sem::ResourceCounter) = nothing
delete!(sem::ResourceCounter) = semrm(sem.handle)

function count(sem::ResourceCounter)
    semget(sem.handle, sem.counts)
    sem.counts
end
count(sem::ResourceCounter, which) = semget(sem.handle, which)

reset(sem::ResourceCounter, val::T, which=0) where {T<:Integer} = semset(sem.handle, Cint(val), which)
reset(sem::ResourceCounter, vals::Vector{T}) where {T<:Integer} = semset(sem.handle, convert(Vector{Cushort},vals))

change(sem::ResourceCounter, by::T, which=0; wait::Bool=true, undo::Bool=sem.undo) where {T<:Integer} = change(sem, [SemBuf(which,by;wait=wait, undo=undo)])
change(sem::ResourceCounter, operations::Vector{SemBuf}) = semop(sem.handle, operations)


text:
struct AspectRatio
    x::Int32
    y::Int32
end


text:
################################################################################################
#
# test 1: Create a blue color image and display it
#
################################################################################################

println("\nOpenCV.jl tests: press ESC to continue")

function test1()
    println("test 1: Create a blank color image and display it")
    # blue color vector (BGR)
    color = cvScalar(255, 0, 0)

    # Create an 8-bit unsigned RGB image, rows(height) x columns(width)
    img = Mat(300, 300, CV_8UC3, color)

    # Display in window
    imdisplay(img, "Blue image")
    closeWindows(0,27,"")  # press ESC(27) to exit
end

test1()

#################################################################################################
#
# test 2: Convert an image from .png to .jpg format (with compression)
#
#################################################################################################

function test2()
    println("test 2: Convert an image from .png to .jpg format (with compression)")
    inputfile = joinpath(Pkg.dir("OpenCV"), "./test/images/lena.png")
    outfile = joinpath(homedir(), "lena_copy.jpeg")
    img = imread(inputfile)

    #  compression parameters
    compression_params = stdvec(0,cint(0))
    stdpush!(compression_params, cint(IMWRITE_JPEG_QUALITY))
    stdpush!(compression_params, cint(98))
    # jpeg-> png
    # push_back(compression_params, cint(IMWRITE_PNG_COMPRESSION))
    # push_back(compression_params, cint(9))

    !(imwrite(outfile, img, compression_params)) ? throw(ArgumentError("Can not save to disk!")) : nothing
    imdisplay(img, "Converted .jpeg image")
    closeWindows(0,27,"")  # press ESC(27) to exit
end

test2()

#################################################################################################
#
# test 3: Binary thresholding
#
#################################################################################################

function test3()
    println("test 3: Binary thresholding")
    filename = joinpath(Pkg.dir("OpenCV"), "./test/images/lena.png")
    img = imread(filename)
    dst = Mat(int(rows(img)), int(cols(img)), CV_8UC1)
    cvtColor(img, dst, COLOR_BGR2GRAY)
    gaussianBlur(dst, dst, cvSize(5,5))
    thresh = Mat()
    threshold(dst, thresh, 120, 255, THRESH_BINARY)  # thresh = 0, max = 255
    #THRESH_OTSU
    #THRESH_BINARY_INV
    #THRESH_TRUNC
    #THRESH_TOZERO
    #THRESH_TOZERO_INV
    imdisplay(img, "Original")
    imdisplay(dst, "Gaussian blur")
    imdisplay(thresh, "Thresholded")
    closeWindows(0,27,"")  # press ESC(27) to exit
end

test3()

################################################################################################
#
# test 4: Live video capture (webcam)
#
################################################################################################

# see src/OpenCV_util.jl for details

function test4()
    println("test 4: Live video capture")
    videocam()    # press ESC to stop
end

test4()

################################################################################################
#
# test 5: Basic video display: set video properties
#
################################################################################################

function test5()
   println("test 5: Basic video display: set video properties")
   fvideo = joinpath(Pkg.dir("OpenCV"), "./test/images/movie.avi")
   vid = videoCapture(fvideo)
   !isOpened(vid) ? throw(ArgumentError("Can not open video stream!")) : nothing

   setVideoId(vid, CAP_PROP_POS_MSEC, 300.0)     #300 ms into the video
   # setVideoId(vid, CAP_PROP_FRAME_WIDTH, 600.0)
   # setVideoId(vid, CAP_PROP_FRAME_HEIGHT, 400.0)
   # setVideoId(vid, CAP_PROP_FPS, 5.0)

   fps = getVideoId(vid, CAP_PROP_FPS)
   println("Frames per second: ", fps)

   namedWindow("Video player")
   frame = Mat()

   # Loop until user presses ESC or frame is empty
   while(true)
      if !(videoRead(vid, frame))
         throw(ArgumentError("Can not acquire video stream!"))
         break
      end

      imshow("Video player", frame)

      if (waitkey(30) == 27)
          destroyAllWindows()
          break
     end
    end
end

test5()


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Marcos L. Alencastro,  Estefane S. Horn (Revised Gerson B. Bicca)
#* $Id$
#*--------------------------------------------------------------------
type centrifugal_compressor
	centrifugal_compressor()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of chemical components",
				:Lower=>1
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Constant of Gases",
				:Unit=>"kJ/kmol/K",
				:Default=>8.31451,
				:Hidden=>true
			)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Molar Weight",
				:Hidden=>true
			)),(NComp)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Compressor Operation Model",
				:Valid=>["Polytropic","Isentropic"],
				:Default=>"Isentropic"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Method of Calculation",
				:Valid=>["GPSA Method","ASME Method"],
				:Default=>"GPSA Method"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Pressure Ratio",
				:Lower=>1E-6,
				:Symbol=>"P_{ratio}"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Drop",
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Increase",
				:Lower=>0,
				:DisplayUnit=>"kPa",
				:Symbol=>"P_{incr}"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Compressor efficiency - Polytropic or Isentropic (See Compressor Type)",
				:Lower=>1E-3,
				:Upper=>1
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Mechanical efficiency",
				:Lower=>1E-3,
				:Upper=>1
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Power Losses"
			)),
			energy_mass(Dict{Symbol,Any}(
				:Brief=>"Actual Head",
				:Protected=>true
			)),
			energy_mass(Dict{Symbol,Any}(
				:Brief=>"Isentropic Head",
				:Protected=>true
			)),
			energy_mass(Dict{Symbol,Any}(
				:Brief=>"Polytropic Head",
				:Protected=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Schultz Polytropic Head Correction",
				:Protected=>true
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Fluid Power",
				:Protected=>true
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Brake Power",
				:Protected=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Polytropic Coefficient",
				:Lower=>0.2,
				:Protected=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Isentropic Coefficient",
				:Lower=>0.2,
				:Protected=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Polytropic efficiency",
				:Lower=>1E-3,
				:Upper=>1,
				:Protected=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Isentropic efficiency",
				:Lower=>1E-3,
				:Upper=>1,
				:Protected=>true
			)),
			temperature(Dict{Symbol,Any}(
				:Brief=>"Isentropic Temperature",
				:Protected=>true
			)),
			enth_mol(Dict{Symbol,Any}(
				:Brief=>"Enthalpy at constant entropy",
				:Hidden=>true
			)),
			molweight(Dict{Symbol,Any}(
				:Brief=>"Mixture Molar Weight",
				:Hidden=>true
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Mass Density at inlet conditions",
				:Lower=>1E-6,
				:Protected=>true
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Mass Density at outlet conditions",
				:Lower=>1E-6,
				:Protected=>true
			)),
			dens_mass(Dict{Symbol,Any}(
				:Brief=>"Mass Density at isentropic conditions",
				:Lower=>1E-6,
				:Hidden=>true
			)),
			fraction(Dict{Symbol,Any}(
				:Brief=>"Compressibility factor at inlet",
				:Lower=>1E-3,
				:Protected=>true
			)),
			fraction(Dict{Symbol,Any}(
				:Brief=>"Compressibility factor at outlet",
				:Lower=>1E-3,
				:Protected=>true
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet stream",
				:PosX=>0.16,
				:PosY=>1,
				:Symbol=>"_{in}"
			)),
			streamPH(Dict{Symbol,Any}(
				:Brief=>"Outlet stream",
				:PosX=>0.87,
				:PosY=>0.0,
				:Symbol=>"_{out}"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Work Inlet",
				:PosX=>0,
				:PosY=>0.45,
				:Protected=>true
			)),
			[
				:(Outlet.F = Inlet.F),
				:(Outlet.z = Inlet.z),
				:(Mwm = sum(Mw*Inlet.z)),
				:(Outlet.P = Inlet.P * Pratio),
				:(Outlet.P = Inlet.P - Pdrop),
				:(Outlet.P = Inlet.P + Pincrease),
				:(rho_in = PP.VapourDensity(Inlet.T, Inlet.P, Inlet.z)),
				:(rho_out= PP.VapourDensity(Outlet.T, Outlet.P, Outlet.z)),
				:(rho_ise= PP.VapourDensity(Tisentropic, Outlet.P, Outlet.z)),
				:(hise = PP.VapourEnthalpy(Tisentropic, Outlet.P, Outlet.z)),
				:(Zfac_in = PP.VapourCompressibilityFactor(Inlet.T,Inlet.P,Inlet.z)),
				:(Zfac_out = PP.VapourCompressibilityFactor(Outlet.T,Outlet.P,Outlet.z)),
				:(IsentropicEff*(Outlet.h-Inlet.h) = (hise-Inlet.h)),
				:(Head*Mwm = (Outlet.h-Inlet.h)),
				:(PP.VapourEntropy(Tisentropic, Outlet.P, Outlet.z) = PP.VapourEntropy(Inlet.T, Inlet.P, Inlet.z)),
				:(BrakePower = -WorkIn),
				:(BrakePower*MechanicalEff = FluidPower),
				:(PowerLoss = BrakePower - FluidPower),
				:(PolytropicEff*HeadIsentropic = HeadPolytropic*IsentropicEff),
				:(FluidPower = Head*Mwm*Inlet.F),
				:(EfficiencyOperation = IsentropicEff),
				:(IseCoeff*ln(rho_ise/rho_in) = ln(Outlet.P/Inlet.P)),
				:(PolyCoeff*ln(rho_out/rho_in) = ln(Outlet.P/Inlet.P)),
				:(HeadIsentropic*rho_in = (IseCoeff/(IseCoeff-1.001))*Inlet.P*HeadCorrection*((Pratio)^((IseCoeff-1.001)/IseCoeff) - 1)),
				:(HeadPolytropic*rho_in = (PolyCoeff/(PolyCoeff-1.001))*Inlet.P*HeadCorrection*((Pratio)^((PolyCoeff-1.001)/PolyCoeff) - 1)),
				:(HeadCorrection*Mwm*(IseCoeff/(IseCoeff-1.001))*(Outlet.P/rho_ise -Inlet.P/rho_in) = (hise-Inlet.h)),
				:(EfficiencyOperation = IsentropicEff),
				:(PolytropicEff*IseCoeff*(PolyCoeff-1) = PolyCoeff*(IseCoeff-1)),
				:(HeadIsentropic*Mwm*((IseCoeff-1.001)/IseCoeff) = (0.5*Zfac_in+0.5*Zfac_out)*Rgas*Inlet.T*((Pratio)^((IseCoeff-1.001)/IseCoeff) - 1)),
				:(HeadPolytropic*Mwm*((PolyCoeff-1.001)/PolyCoeff) = (0.5*Zfac_in+0.5*Zfac_out)*Rgas*Inlet.T*((Pratio)^((PolyCoeff-1.001)/PolyCoeff) - 1)),
				:(HeadCorrection =1),
				:(HeadIsentropic = Head*IsentropicEff),
				:(EfficiencyOperation = PolytropicEff),
				:(PolytropicEff*IseCoeff*(PolyCoeff-1) = PolyCoeff*(IseCoeff-1)),
				:(HeadIsentropic*Mwm*((IseCoeff-1.001)/IseCoeff) = (0.5*Zfac_in+0.5*Zfac_out)*Rgas*Inlet.T*((Pratio)^((IseCoeff-1.001)/IseCoeff) - 1)),
				:(HeadPolytropic*Mwm*((PolyCoeff-1.001)/PolyCoeff) = (0.5*Zfac_in+0.5*Zfac_out)*Rgas*Inlet.T*((Pratio)^((PolyCoeff-1.001)/PolyCoeff) - 1)),
				:(HeadCorrection =1),
				:(HeadIsentropic = Head*IsentropicEff),
				:(EfficiencyOperation = PolytropicEff),
				:(IseCoeff*ln(rho_ise/rho_in) = ln(Outlet.P/Inlet.P)),
				:(PolyCoeff*ln(rho_out/rho_in) = ln(Outlet.P/Inlet.P)),
				:(HeadIsentropic*rho_in = (IseCoeff/(IseCoeff-1.001))*Inlet.P*HeadCorrection*((Pratio)^((IseCoeff-1.001)/IseCoeff) - 1)),
				:(HeadPolytropic*rho_in = (PolyCoeff/(PolyCoeff-1.001))*Inlet.P*HeadCorrection*((Pratio)^((PolyCoeff-1.001)/PolyCoeff) - 1)),
				:(HeadCorrection*Mwm*(IseCoeff/(IseCoeff-1.001))*(Outlet.P/rho_ise -Inlet.P/rho_in) = (hise-Inlet.h)),
			],
			[
				"Overall Molar Balance","Component Molar Balance","Average Molecular Weight","Pressure Ratio","Pressure Drop","Pressure Increase","Mass Density at inlet conditions","Mass Density at outlet conditions","Mass Density at isentropic conditions","Enthalpy at isentropic conditions","Compressibility factor at Inlet Conditions","Compressibility factor at Outlet Conditions","Isentropic Efficiency","Actual Head","Isentropic Outlet Temperature","Brake Power","Brake Power","Power Loss","Polytropic-Isentropic Relation","Fluid Power","Efficiency","Isentropic Coefficient","Polytropic Coefficient","Isentropic Head","Polytropic Head","Schultz Polytropic Head Correction","Efficiency","Polytropic Efficiency","Isentropic Coefficient","Polytropic Coefficient","Head Correction","Isentropic Head","Efficiency","Polytropic Efficiency","Isentropic Coefficient","Polytropic Coefficient","Head Correction","Isentropic Head","Efficiency","Isentropic Coefficient","Polytropic Coefficient","Isentropic Head","Polytropic Head","Schultz Polytropic Head Correction",
			],
			[:PP,:NComp,:Rgas,:Mw,:CompressorOperation,:Method,],
			[:Pratio,:Pdrop,:Pincrease,:EfficiencyOperation,:MechanicalEff,:PowerLoss,:Head,:HeadIsentropic,:HeadPolytropic,:HeadCorrection,:FluidPower,:BrakePower,:PolyCoeff,:IseCoeff,:PolytropicEff,:IsentropicEff,:Tisentropic,:hise,:Mwm,:rho_in,:rho_out,:rho_ise,:Zfac_in,:Zfac_out,:Inlet,:Outlet,:WorkIn,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Rgas::positive
	Mw::Array{molweight}
	CompressorOperation::DanaSwitcher
	Method::DanaSwitcher
	Pratio::positive
	Pdrop::press_delta
	Pincrease::press_delta
	EfficiencyOperation::positive
	MechanicalEff::positive
	PowerLoss::power
	Head::energy_mass
	HeadIsentropic::energy_mass
	HeadPolytropic::energy_mass
	HeadCorrection::positive
	FluidPower::power
	BrakePower::power
	PolyCoeff::positive
	IseCoeff::positive
	PolytropicEff::positive
	IsentropicEff::positive
	Tisentropic::temperature
	hise::enth_mol
	Mwm::molweight
	rho_in::dens_mass
	rho_out::dens_mass
	rho_ise::dens_mass
	Zfac_in::fraction
	Zfac_out::fraction
	Inlet::stream
	Outlet::streamPH
	WorkIn::power
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export centrifugal_compressor
function set(in::centrifugal_compressor)
	Mw = PP.MolecularWeight()
	 Rgas = 8.31451*"kJ/kmol/K"

end
function setEquationFlow(in::centrifugal_compressor)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	addEquation(18)
	addEquation(19)
	addEquation(20)
	let switch=CompressorOperation
		if switch=="Isentropic"
			let switch=Method
				if switch=="ASME Method"
					addEquation(21)
					addEquation(22)
					addEquation(23)
					addEquation(24)
					addEquation(25)
					addEquation(26)
				elseif switch=="GPSA Method"
					addEquation(27)
					addEquation(28)
					addEquation(29)
					addEquation(30)
					addEquation(31)
					addEquation(32)
				end
			end
		elseif switch=="Polytropic"
			let switch=Method
				if switch=="GPSA Method"
					addEquation(33)
					addEquation(34)
					addEquation(35)
					addEquation(36)
					addEquation(37)
					addEquation(38)
				elseif switch=="ASME Method"
					addEquation(39)
					addEquation(40)
					addEquation(41)
					addEquation(42)
					addEquation(43)
					addEquation(44)
				end
			end
		end
	end
end
function atributes(in::centrifugal_compressor,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/CentrifugalCompressor"
	fields[:Brief]="Model of a centrifugal compressor."
	fields[:Info]="To be documented

== References ==

[1] GPSA, 1979, Engineering Data Book, Chapter 4, 5-9 - 5-10.

[2] Bloch, Heinz P., A Practical Guide to Compressor Technology, John Wiley & Sons, Incorporate, 2006.

[3] Mark R. Sandberg, Equation Of State Influences On Compressor Performance Determination,PROCEEDINGS OF THE THIRTY-FOURTH TURBOMACHINERY SYMPOSIUM, 2005."
	drive!(fields,_)
	return fields
end
centrifugal_compressor(_::Dict{Symbol,Any})=begin
	newModel=centrifugal_compressor()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using Persistence
using Persistence.Persist
import Persistence.Lib

using BenchmarkTools

function transact(P::PersistentArray, values)
    transaction(P) do
        _copy(P, values)
    end
end

nontransact(P::PersistentArray, values) = _copy(P, values)

_copy(P::PersistentArray, values) = unsafe_copyto!(P.base, pointer(values), length(values))

function test(size, iters = 100000)
    file = "/mnt/test.pool"
    ispath(file) && rm(file)
    pool = Lib.create(file, 2^30)

    P = PersistentArray{Int64}(pool, (size,))
    values = rand(Int64, size)
    for _ in 1:iters
        _copy(P, values)
    end
end

function benchmark()
    sizes = [10 ^ i for i in 3:7]

    poolfile = "/mnt/test.pool"

    for size in sizes
        try
            println("Working with file $poolfile")
            ispath(poolfile) && rm(poolfile)
            pool = Lib.create(poolfile, 2^30)

            P = PersistentArray{Int64}(pool, (size,))
            indices = 1:size
            values = rand(Int64, size)

            results = @benchmark nontransact($P, $indices, $values)

            println("Results for size: $size")
            display(results)

            Lib.close(pool)
        finally
            rm(poolfile)
        end
    end
end


text:
# Structure:
# Generic algorithms pre_partition, pre_mc for different kinds of queries
# These are implemented by concrete algorithms which specialise on type

"Extracts properties (e.g. samples, probabilities) from random variables"
abstract type InferenceAlgorithm end

"Generate Partitions of Sets"
abstract type PartitionAlgorithm <: InferenceAlgorithm end
abstract type SamplingAlgorithm <: InferenceAlgorithm end
abstract type MCMCAlgorithm <: SamplingAlgorithm end

"A stop function used as dummy stopping criteria in while loops"
neverstop(_...) = true

"is this box small (enough to be accepted)?"
function issmall(box::Boxes, precision::Float64)::Bool
  for dim in dims(box)
    (measure(box[dim]) > precision) && return false
  end
  return true
end

include("partition.jl")
include("chain.jl")
include("mcmc.jl")

# Concrete algorithms
include("ams.jl")
include("bfs.jl")


text:
using BenchmarkTools
using Statistics
using StatsBase: sample, quantile
using Distributions
using Random
using Plots
using Revise
using EvoTrees
using EvoTrees: sigmoid, logit
# using ProfileView

# prepare a dataset
Random.seed!(12)
features = rand(10_000) .* 5
X = reshape(features, (size(features)[1], 1))
Y = sin.(features) .* 0.5 .+ 0.5
Y = logit(Y) + randn(size(Y))
Y = sigmoid(Y)
𝑖 = collect(1:size(X, 1))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace=false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1)) + 1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]

# linear
params1 = EvoTreeRegressor(T=Float32,
    loss=:linear, metric=:mse,
    nrounds=200, nbins=64,
    λ=0.5, γ=0.1, η=0.1,
    max_depth=6, min_weight=1.0,
    rowsample=0.1, colsample=1.0,
    device="gpu")

@time model = fit_evotree(params1, X_train, Y_train, print_every_n=25);
# @time model = fit_evotree(params1, X_train, Y_train, X_eval=X_eval, Y_eval=Y_eval, print_every_n=25);
# 67.159 ms (77252 allocations: 28.06 MiB)
model_cpu = convert(EvoTrees.GBTree, model);
pred_train_linear_gpu = predict(model, X_train)
pred_train_linear_cpu = predict(model_cpu, X_train)
sum(pred_train_linear_gpu .- pred_train_linear_cpu)

# @btime model = grow_gbtree($X_train, $Y_train, $params1, X_eval = $X_eval, Y_eval = $Y_eval, print_every_n = 25, metric=:mae)
@time pred_train_linear = predict(model, X_train)
mean(abs.(pred_train_linear .- Y_train))
sqrt(mean((pred_train_linear .- Y_train).^2))

# logistic / cross-entropy
params1 = EvoTreeRegressor(T=Float32,
    loss=:logistic, metric=:logloss,
    nrounds=200, nbins=64,
    λ=0.5, γ=0.1, η=0.1,
    max_depth=6, min_weight=1.0,
    rowsample=0.5, colsample=1.0,
    device="gpu")

@time model = fit_evotree(params1, X_train, Y_train, print_every_n=25);
# @time model = fit_evotree(params1, X_train, Y_train, X_eval=X_eval, Y_eval=Y_eval, print_every_n=25);
# 218.040 ms (123372 allocations: 34.71 MiB)
# @btime model = fit_evotree($params1, $X_train, $Y_train, X_eval = $X_eval, Y_eval = $Y_eval)
@time pred_train_logistic = predict(model, X_train)
sqrt(mean((pred_train_logistic .- Y_train).^2))

x_perm = sortperm(X_train[:,1])
plot(X_train, Y_train, msize=1, mcolor="gray", mswidth=0, background_color=RGB(1, 1, 1), seriestype=:scatter, xaxis=("feature"), yaxis=("target"), legend=true, label="")
plot!(X_train[:,1][x_perm], pred_train_linear[x_perm], color="navy", linewidth=1.5, label="Linear")
plot!(X_train[:,1][x_perm], pred_train_logistic[x_perm], color="darkred", linewidth=1.5, label="Logistic")
# plot!(X_train[:,1][x_perm], pred_train_poisson[x_perm], color = "green", linewidth = 1.5, label = "Poisson")
# plot!(X_train[:,1][x_perm], pred_train_L1[x_perm], color = "pink", linewidth = 1.5, label = "L1")
savefig("figures/regression_sinus_gpu.png")

###############################
## gaussian
###############################
params1 = EvoTreeGaussian(T=Float32,
    loss=:gaussian, metric=:gaussian,
    nrounds=200, nbins=64,
    λ=1.0, γ=0.1, η=0.1,
    max_depth=6, min_weight=0.1,
    rowsample=0.5, colsample=1.0, rng=123,
    device="gpu")

@time model = fit_evotree(params1, X_train, Y_train, print_every_n=25);
# @time model = fit_evotree(params1, X_train, Y_train, X_eval=X_eval, Y_eval=Y_eval, print_every_n=25);
# @time model = fit_evotree(params1, X_train, Y_train, print_every_n = 10);
@time pred_train_gaussian = EvoTrees.predict(model, X_train)

pred_gauss = [Distributions.Normal(pred_train_gaussian[i,1], pred_train_gaussian[i,2]) for i in 1:size(pred_train_gaussian, 1)]
pred_q80 = quantile.(pred_gauss, 0.8)
pred_q20 = quantile.(pred_gauss, 0.2)

mean(Y_train .< pred_q80)
mean(Y_train .< pred_q20)

x_perm = sortperm(X_train[:,1])
plot(X_train[:, 1], Y_train, ms=1, mcolor="gray", mswidth=0, background_color=RGB(1, 1, 1), seriestype=:scatter, xaxis=("feature"), yaxis=("target"), legend=true, label="")
plot!(X_train[:,1][x_perm], pred_train_gaussian[x_perm, 1], color="navy", linewidth=1.5, label="mu")
plot!(X_train[:,1][x_perm], pred_train_gaussian[x_perm, 2], color="darkred", linewidth=1.5, label="sigma")
plot!(X_train[:,1][x_perm], pred_q20[x_perm, 1], color="green", linewidth=1.5, label="q20")
plot!(X_train[:,1][x_perm], pred_q80[x_perm, 1], color="green", linewidth=1.5, label="q80")
savefig("figures/gaussian_sinus_gpu.png")


text:
using SpatialDatasets
using Documenter

DocMeta.setdocmeta!(SpatialDatasets, :DocTestSetup, :(using SpatialDatasets); recursive=true)

makedocs(;
    modules=[SpatialDatasets],
    authors="Javier Barbero and contributors",
    repo="https://github.com/javierbarbero/SpatialDatasets.jl/blob/{commit}{path}#{line}",
    sitename="SpatialDatasets.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://javierbarbero.github.io/SpatialDatasets.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "List of Spatial Datasets" => "sdatasets.md",
        "API" => "api.md",
    ],
)

deploydocs(;
    repo="github.com/javierbarbero/SpatialDatasets.jl",
    devbranch = "main"
)


text:
eltype_nullable(x::Nullable) = eltype(x)
eltype_nullable(x) = typeof(x)
eltype_nullable{T<:Nullable}(::Type{T}) = eltype(T)
eltype_nullable{T}(::Type{T}) = T

eltypes(x) = Tuple{eltype_nullable(x)}
eltypes(x, xs...) = Tuple{eltype_nullable(x), eltypes(xs...).parameters...}

"""
    lift(f, xs...)

Lift function `f`, passing it arguments `xs...`, using standard lifting semantics:
for a function call `f(xs...)`, return null if any `x` in `xs` is null; otherwise,
return `f` applied to values of `xs`.
"""
@inline @generated function lift{F}(f::F, xs::Tuple)
    N = nfields(xs)
    args = (:(unsafe_get(xs[$i])) for i in 1:N)
    checknull = (:(!isnull(xs[$i])) for i in 1:N)
    if isdefined(f, :instance) && null_safe_op(f.instance, map(eltype_nullable, xs.parameters)...)
        return quote
            val = f($(args...))
            nonull = (&)($(checknull...))
            @compat Nullable(val, nonull)
        end
    elseif VERSION >= v"0.6.0-dev"
      return quote
          U = Core.Inference.return_type(f, eltypes(xs...))
          if (&)($(checknull...))
              return Nullable(f($(args...)))
          else
              return isleaftype(U) ? Nullable{U}() : Nullable()
          end
      end
    else # Inference fails with the previous branch on Julia 0.5
      if isdefined(f, :instance)
          U = Core.Inference.return_type(f.instance, map(eltype_nullable, xs.parameters))
          isleaftype(U) || (U = Union{})
      elseif F === DataType # Function is a constructor
          U = f.parameters[1]
      else
          U = Union{}
      end
      return quote
          if (&)($(checknull...))
              return Nullable(f($(args...)))
          else
              return Nullable{$U}()
          end
      end
    end
end

lift(f) = Nullable(f())


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

"""
    abstract type DifferentiationAlgorithm

*Experimental feature, not part of stable public API.*

Abstract type for integrated autocorrelation length estimation algorithms.
"""
abstract type DifferentiationAlgorithm end
export DifferentiationAlgorithm

"""
    bat_valgrad(f::Function, [algorithm::DifferentiationAlgorithm])

*Experimental feature, not part of stable public API.*

Generated a function that calculates both value and gradient of `f` at given
points.

The function `f` must support `ValueShapes.varshape(f)` and
`ValueShapes.unshaped(f)`.

Returns a NamedTuple of the shape

```julia
(result = valgrad_f, ...)
```

with

```julia
f_at_x, grad_of_f_at_x = valgrad_f(x)

grad_of_f_at_x = zero(x)
f_at_x = valgrad_f(!, grad_of_f_at_x, x)
```

Result properties not listed here are algorithm-specific and are not part
of the stable public API.

!!! note

    Do not add add methods to `bat_valgrad`, add methods to
    `bat_valgrad_impl` instead.
"""
function bat_valgrad end
export bat_valgrad

function bat_valgrad_impl end

function bat_valgrad(
    f::Function,
    algorithm = bat_default_withdebug(bat_valgrad, Val(:algorithm), f)
)
    r = bat_valgrad_impl(f, algorithm)
    result_with_args(r, (algorithm = algorithm,))
end

function argchoice_msg(::typeof(bat_valgrad), ::Val{:algorithm}, x::DifferentiationAlgorithm)
    "Using automiatic differentiation algorithm $x"
end

# ToDo, move to BAT.jl:
#=
gradient_shape(vs::AbstractValueShape) = replace_const_shapes(ValueShapes.const_zero_shape, vs)

abstract type GradientFunction end
=#


text:
peek(v::Vector{T}) where {T} = v[end]

function shove!(a::Vector{T}, i::Integer) where {T}
  ln = length(a)
  i = min(max(1, ln - i), ln + 1)
  pop!(insert!(a, i, a[end]))
  return a
end

function yank!(a::Vector{T}, i::Integer) where {T}
  ln = length(a)
  i = min(max(1, ln - i), ln)
  push!(a, splice!(a, i))
end

function yankdup!(a::Vector{T}, i::Integer) where {T}
  ln = length(a)
  i = min(max(1, ln - i), ln)
  push!(a, a[i])
end


text:
##
## Angle Conversions
##

const MINUTES_TO_WHOLE = inv(60)
const SECONDS_TO_WHOLE = inv(3600)
const HOURS_PER_DEGREE = 24 / 360
const HOURS_PER_RADIAN = 12 / π
const DEGREES_PER_HOUR = 360 / 24
const RADIANS_PER_HOUR = π / 12

### rad2xxx

"""
    rad2ha(angle)

Convert radians to hour angles
"""
rad2ha(angle) = angle * HOURS_PER_RADIAN

"""
    rad2dms(angle)

Convert radians to (degrees, arcminutes, arcseconds) tuple
"""
rad2dms(angle) = rad2deg(angle) |> deg2dms

"""
    rad2hms(angle)

Convert radians to (hours, minutes, seconds) tuple
"""
rad2hms(angle) = rad2ha(angle) |> ha2hms

### deg2xxx

"""
    deg2ha(angle)

Convert degrees to hour angles
"""
deg2ha(angle) = angle * HOURS_PER_DEGREE

"""
    deg2dms(angle)

Convert degrees to (degrees, arcminutes, arcseconds) tuple
"""
function deg2dms(angle)
    remain_degrees, degrees = modf(angle)
    fraction_arcmin = abs(remain_degrees) * 60
    remain_arcmin, arcmin = modf(fraction_arcmin)
    arcsec = remain_arcmin * 60
    return degrees, arcmin, arcsec
end

"""
    deg2hms(angle)

Convert degrees to (hours, minutes, seconds) tuple
"""
deg2hms(angle) = deg2ha(angle) |> ha2hms

### ha2xxx

"""
    ha2rad(angle)

Convert hour angles to radians
"""
ha2rad(angle) = angle * RADIANS_PER_HOUR

"""
    ha2deg(angle)

Convert hour angles to degrees
"""
ha2deg(angle) = angle * DEGREES_PER_HOUR

"""
    ha2hms(angle)

Convert hour angles to (hours, minutes, seconds) tuple
"""
function ha2hms(angle)
    remain_hours, hours = modf(angle)
    fraction_min = abs(remain_hours) * 60
    remain_min, minutes = modf(fraction_min)
    seconds = remain_min * 60
    return hours, minutes, seconds
end

"""
    ha2dms(angle)

Convert hour angles to (degrees, arcminutes, arcseconds) tuple
"""
ha2dms(angle) = ha2deg(angle) |> deg2dms

### dms2xxx

"""
    dms2deg(degrees, arcmin, arcsec)
    dms2deg(parts)
    dms2deg(input::AbstractString)

Convert (degrees, arcminutes, arcseconds) tuple to degrees. If a string is given, will parse with [`parse_dms`](@ref) first. If an angle is input will treat as a no-op.
"""
function dms2deg(degrees, arcminutes, arcseconds)
    frac = arcminutes * MINUTES_TO_WHOLE + arcseconds * SECONDS_TO_WHOLE
    return signbit(degrees) ? degrees - frac : degrees + frac
end

"""
    dms2rad(degrees, arcmin, arcsec)
    dms2rad(parts)
    dms2rad(input::AbstractString)

Convert (degrees, arcminutes, arcseconds) tuple to radians. If a string is given, will parse with [`parse_dms`](@ref) first. If an angle is input will treat as a no-op.
"""
dms2rad(degrees, arcminutes, arcseconds) = dms2deg(degrees, arcminutes, arcseconds) |> deg2rad

"""
    dms2ha(degrees, arcmin, arcsec)
    dms2ha(parts)
    dms2ha(input::AbstractString)

Convert (degrees, arcminutes, arcseconds) tuple to hour angles. If a string is given, will parse with [`parse_dms`](@ref) first. If an angle is input will treat as a no-op.
"""
dms2ha(degrees, arcminutes, arcseconds) = dms2deg(degrees, arcminutes, arcseconds) |> deg2ha

# code-gen for string inputs and no-ops
for func in (:dms2deg, :dms2rad, :dms2ha)
    @eval $func(input::AbstractString) = parse_dms(input) |> $func
    @eval $func(input::Number) = input
end

### hms2xxx

"""
    hms2ha(hours, mins, secs)
    hms2ha(parts)
    hms2ha(input::AbstractString)

Convert (hours, minutes, seconds) tuple to hour angles. If a string is given, will parse with [`parse_hms`](@ref) first. If an angle is input will treat as a no-op.
"""
function hms2ha(hours, minutes, seconds)
    frac = minutes * MINUTES_TO_WHOLE + seconds * SECONDS_TO_WHOLE
    return signbit(hours) ? hours - frac : hours + frac
end

"""
    hms2deg(hours, mins, secs)
    hms2deg(parts)
    hms2deg(input::AbstractString)

Convert (hours, minutes, seconds) tuple to degrees. If a string is given, will parse with [`parse_hms`](@ref) first. If an angle is input will treat as a no-op.
"""
hms2deg(hours, minutes, seconds) = hms2ha(hours, minutes, seconds) |> ha2deg

"""
    hms2rad(hours, mins, secs)
    hms2rad(parts)
    hms2rad(input::AbstractString)

Convert (hours, minutes, seconds) tuple to radians. If a string is given, will parse with [`parse_hms`](@ref) first. If an angle is input will treat as a no-op.
"""
hms2rad(hours, minutes, seconds) = hms2ha(hours, minutes, seconds) |> ha2rad

# code-gen for string inputs and no-ops
for func in (:hms2deg, :hms2rad, :hms2ha)
    @eval $func(input::AbstractString) = parse_hms(input) |> $func
    @eval $func(input::Number) = input
end

# code-gen for accepting inputs separate or in a collection
for func in (:dms2rad, :dms2deg, :dms2ha, :hms2rad, :hms2deg, :hms2ha)
    @eval $func(parts) = $func(parts...)
end


text:
using LinearAlgebra
using CompScienceMeshes
using SauterSchwab3D

const pI = point(3,4,2)
const pII = point(3,7,4)
const pIII = point(8,3,5)

const P = simplex(pI, pII, pIII)
const Q = simplex(pI, pII, pIII)

const sing = SauterSchwab3D.singularity_detection(P,Q)

Accuracy2 = 24
cf_ref = CommonFace4D(sing,SauterSchwab3D._legendre(Accuracy2,0.0,1.0))

function integrand(x,y)
      return ((x-pI)'*(y-pII))*exp(-im*1*norm(x-y))/(4pi*norm(x-y))
end

function INTEGRAND(u,v)
   n1 = neighborhood(P,u)
   n2 = neighborhood(Q,v)
   x = cartesian(n1)
   y = cartesian(n2)
   output = integrand(x,y)*jacobian(n1)*jacobian(n2)
   return(output)
end

print("Ref: ")
ref = sauterschwab_parameterized(INTEGRAND, cf_ref)
println(ref)
println()

res_tp =[]
res_sp =[]
res_gm =[]

n1 = []
n2 = []
n3 = []
for i in 2:1:14
   Accuracy = i
   cf = CommonFace4D(sing,SauterSchwab3D._legendre(Accuracy,0.0,1.0))

   int_tp = sauterschwab_parameterized(INTEGRAND, cf)

   num_pts = 6*length(cf.qps)^4
   push!(n1,num_pts)

   push!(res_tp,int_tp)
end

for i in 2:1:7
   Accuracy = i
   cf_s = CommonFace4D_S(sing,(SauterSchwab3D._legendre(Accuracy,0.0,1.0),
                          SauterSchwab3D._shunnham3D(Accuracy)))

   int_sp = sauterschwab_parameterized(INTEGRAND, cf_s)

   num_pts = 6*length(cf_s.qps[1])*length(cf_s.qps[2])
   push!(n2,num_pts)
   push!(res_sp,int_sp)

end

for i in 2:1:12
   Accuracy = i

   cf_gm = CommonFace4D_S(sing,(SauterSchwab3D._legendre(Accuracy,0.0,1.0),
                          SauterSchwab3D._grundmannMoeller3D(2*(Accuracy-1)+1)))

   int_gm = sauterschwab_parameterized(INTEGRAND, cf_gm)

   num_pts = 6*length(cf_gm.qps[1])*length(cf_gm.qps[2])
   push!(n3,num_pts)
   push!(res_gm,int_gm)
end

err_tp = norm.(res_tp.-ref)/norm(ref)
err_sp = norm.(res_sp.-ref)/norm(ref)
err_gm = norm.(res_gm.-ref)/norm(ref)

println(err_tp)
println(err_sp)
println(err_gm)

using Plots
plot( yaxis=:log, xaxis=:log, fontfamily="Times")
plot!(n1,err_tp, label="Gauss Tensor-Product",markershape=:circle)
plot!(n2,err_sp, label="Simplex Tensor-Product",markershape=:rect)
#plot!(n3,err_gm, label="Simplex-Product GM",markershape=:x)
plot!(xlims=(1e1,1e4),ylims=(1e-6,1))
plot!(xlabel="#Quad. pts/Func. evals", ylabel="Rel. Error.", title="Common Face 4D",legend=:bottomleft)

using BenchmarkTools

ref = 9.199951912100538 - 5.619050336211777im

cf = CommonFace4D(sing,SauterSchwab3D._legendre(5,0.0,1.0))

int_tp = sauterschwab_parameterized(INTEGRAND, cf)

num_pts = 6*length(cf.qps)^4
println("#Pts: ",num_pts)
err_tp = norm.(int_tp.-ref)/norm(ref)
println("Rel. Err: ",err_tp)
@time for i in 1:100 sauterschwab_parameterized(INTEGRAND, cf) end

cf_s = CommonFace4D_S(sing,(SauterSchwab3D._legendre(6,0.0,1.0),
SauterSchwab3D._shunnham3D(6)))

int_sp = sauterschwab_parameterized(INTEGRAND, cf_s)

num_pts = 6*length(cf_s.qps[1])*length(cf_s.qps[2])

println("#Pts: ",num_pts)
err_sp = norm.(int_sp.-ref)/norm(ref)
println("Rel Err: ",err_sp)
@time for i in 1:100 sauterschwab_parameterized(INTEGRAND, cf_s) end

#=
cf_gm = CommonFace4D_S((SauterSchwab3D._legendre(6,0.0,1.0),
SauterSchwab3D._grundmannMoeller3D(2*6-1)))

int_gm = sauterschwab_parameterized(INTEGRAND, cf_gm)

num_pts = 6*length(cf_gm.qps[1])*length(cf_gm.qps[2])
println("#Pts: ",num_pts)
err_gm = norm.(int_gm.-ref)/norm(ref)
println("Rel Err: ",err_gm)
=#


text:
module FluxMPI

include("mpi_extensions.jl")

import ChainRulesCore: @non_differentiable
import ComponentArrays: ComponentArray, getdata, getaxes
import CUDA
import Dates: now
import Functors: fmap
import LearnBase
import MLDataUtils
import MLUtils: getobs, numobs
import MPI
import MPI: Barrier, Comm_rank, Comm_size, COMM_WORLD
import .MPIExtensions: allreduce!, bcast!, cpu, gpu, reduce!, Iallreduce!, Ibcast!
import Optimisers: Leaf, init, apply!
import Setfield: @set!

# General Utilities -- Init, Clean Printing
include("common.jl")

# synchronize! Model Parameters -- Works with both Zygote.Params which Flux uses and NamedTuples for ExplicitFluxLayers
include("synchronize.jl")

# Implementation of Distributed Optimiser
include("optimiser.jl")

# Extends LearnBase & MLDataUtils API for Distributed Datasets -- compatible with DataLoaders.jl
include("data.jl")

export local_rank, total_workers, DistributedOptimiser, clean_print, clean_println, DistributedDataContainer, allreduce_gradients

end


text:
baremodule ZeroBasedIndexing
using Main: Base
export zerobased

for n = 1:5
    Base.eval(ZeroBasedIndexing, :(getindex(A, $([Symbol("i$i") for i = 1:n]...)) =
                                  Base.getindex(A, $([:($(Symbol("i$i"))+1) for i = 1:n]...))))
end
getindex(A, inds...) = Base.getindex(A, map(x->x+1, inds)...)

for n = 1:5
    Base.eval(ZeroBasedIndexing, :(setindex!(A, X, $([Symbol("i$i") for i = 1:n]...)) =
                                  Base.setindex!(A, X, $([:($(Symbol("i$i"))+1) for i = 1:n]...))))
end
setindex!(A, X, inds...) = Base.setindex!(A, X, map(x->x+1, inds)...)

function zerobased(mod=Base.current_module())
    Base.eval(mod, quote
        const getindex = $getindex
        const setindex! = $setindex!
    end)
    nothing
end
end # module


text:
module EMST

export compute_emst , verify_emst

include("emst_dual_boruvka.jl")
include("emst_validation.jl")

end # module


text:
import AbstractAlgebra.Generic: Perm, SymmetricGroup

# disambiguation
GroupsCore.order(::Type{I}, G::SymmetricGroup) where {I<:Integer} =
    convert(I, factorial(G.n))

# disambiguation
GroupsCore.order(::Type{I}, g::Perm) where {I<:Integer} =
    convert(I, foldl(lcm, length(c) for c in AbstractAlgebra.cycles(g)))

# correct the AA length:
Base.length(G::SymmetricGroup) = order(Int, G)

# genuinely new methods:
Base.IteratorSize(::Type{<:AbstractAlgebra.AbstractPermutationGroup}) = Base.HasLength()

function GroupsCore.gens(G::SymmetricGroup{I}) where {I}
    G.n == 1 && return eltype(G)[]
    if G.n == 2
        a = one(G)
        a.d[1], a.d[2] = 2, 1
        return [a]
    end
    a, b = one(G), one(G)
    circshift!(a.d, b.d, -1)
    b.d[1], b.d[2] = 2, 1
    return [a, b]
end

Base.deepcopy_internal(g::Perm, stackdict::IdDict) =
    Perm(Base.deepcopy_internal(g.d, stackdict), false)


text:
#!/usr/bin/env julia

function process(case, dgs)
    stack = Int[0]
    ret = Union{Int,String}[]
    for digit in dgs
        diff = digit - stack[end]

        if diff != 0
            pop!(stack)
            push!(stack, digit)
        end

        diff > 0 && push!(ret, '('^diff)
        diff < 0 && push!(ret, ')'^-diff)
        push!(ret, digit)
    end

    push!(ret, ')'^stack[end])
    nested = join(ret, "")
    println("Case #$(case): $(nested)")
end

function main()
    T = parse(Int, readline(stdin))
    for case in 1:T
        dgs = parse.(Int, split(readline(stdin), ""))
        process(case, dgs)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
module Painter
using Printf
using Markdown
using Cairo
using ComplexVisual
@ComplexVisual.import_huge

import Main.DocGenerator: DocSource, DocCreationEnvironment, DocContext,
        Document, substitute_marker_in_markdown, create_doc_icon, append_md

"""
# [![./Painter_docicon.png]({image_from_canvas: get_doc_icon()}) Painter](./Painter.md)

Painters have the ability to "draw"/"paint" something inside objects
with math coordinate systems (e.g. `CV_Math2DCanvas`).

Typically a `CV_ContextStyle` is "attached" to a painter (via
`↦:Tuple{CV_ContextStyle, CV_Painter}` or `CV_StyledPainter`) to govern the
appearance (color, font, line width, etc.) of the painting operations.

## Quick links

|  "area" painters     |    curve  painters    |      other painters   |
|:---------------------|:----------------------|:----------------------|
| ![./Painter_fillpainter_icon.png]({image_from_canvas: create_icon(example_fill_painter())}) `CV_FillPainter`             | ![./Painter_linepainter_icon.png]({image_from_canvas: create_icon(example_line_painter())}) `CV_LinePainter`    | ![./Painter_markpainter_icon.png]({image_from_canvas: create_icon(example_mark_painter())}) `CV_ValueMarkPainter` |
| ![./Painter_portraitpainter_icon.png]({image_from_canvas: create_icon(example_portrait_painter())}) `CV_PortraitPainter` | ![./Painter_dirpainter_icon.png]({image_from_canvas: create_icon(example_dir_painter())}) `CV_DirectionPainter` | ![./Painter_gridpainter_icon.png]({image_from_canvas: create_icon(example_grid_painter())}) `CV_GridPainter`      |
| ![./Painter_m2dcanvaspainter_icon.png]({image_from_canvas: create_icon(example_m2d_canvas_painter())}) `CV_Math2DCanvasPainter` |                                                                                | `CV_CombiPainter`   `→:Tuple{T, S} where {T<:CV_Painter, S<:CV_Painter}`                                                                                |
| ![./Painter_windingpainter_icon.png]({image_from_canvas: create_icon(example_winding_painter())}) `CV_WindingPainter`           |                       | `CV_StyledPainter`, `↦:Tuple{CV_ContextStyle, CV_Painter}`    |

construction of line segments: `cv_parallel_lines`  `cv_arc_lines`  `cv_star_lines`

## How painters work

For the  painting operation of a subtype of `CV_Painter` the function
`cv_paint` is called:

```
cv_paint(context::C, painter::P)
    context     C  <: CV_Context    (often C <: CV_2DCanvasContext)
    painter     P  <: CV_Painter    (often P <: CV_2DCanvasPainter)
```

## `doc: CV_Painter`

## `doc: CV_CanvasPainter`

## `doc: CV_2DCanvasPainter`

"""
painter_intro() = nothing

"""
## The `trafo` (or `dst_trafo`) argument

Many painters support a `trafo` argument (which is sometimes called `dst_trafo`
when other transformations are involved) in order to "transform" the objects
before painting them. Such a transformation has always the form

```julia
    trafo(z::ComplexF64) :: ComplexF64
```

If this `trafo` depends on other parameters then this can be used to
change the output of the painters by chaning the parameters (of the `trafo`).

Here is an example.

![./Painter_trafoarg.png]({image_from_canvas: example_trafo()})

```julia
{func: example_trafo}
```
"""
painter_trafo() = nothing

function example_trafo()
    math_canvas = CV_Math2DCanvas(-1.0 + 1.0im, 3.0 - 1.0im, 110)

    bg_fill = cv_white ↦ CV_FillPainter()  # for background
    grid_style = cv_color(0.8, 0.8, 0.8) → cv_linewidth(1)
    grid = grid_style ↦ CV_GridPainter(-1.0:0.2:3.0, -1.0:0.2:1.0)

    n = 0
    trafo = z -> z^n + n - 1

    lp = CV_LinePainter(trafo, [cv_arc_lines(0, 2π, (0.8,))[1] .+ 0.1 .+ 0.1im])
    style1 = cv_color(1,0,0) → cv_linewidth(2)
    style2 = cv_color(0,1,0) → cv_linewidth(2)
    style3 = cv_color(0,0,1) → cv_linewidth(2)

    cv_create_context(math_canvas) do canvas_context
        cv_paint(canvas_context, bg_fill)
        cv_paint(canvas_context, grid)

        n = 1; cv_paint(canvas_context, style1 ↦ lp)
        n = 2; cv_paint(canvas_context, style2 ↦ lp)
        n = 3; cv_paint(canvas_context, style3 ↦ lp)
    end

    return math_canvas
end

"""
## `doc: CV_FillPainter`

## Example for `CV_FillPainter`

![./Painter_fillpainter.png]({image_from_canvas: example_fill_painter()})

```julia
{func: example_fill_painter}
```
"""
help_fill_painter() = nothing

function example_fill_painter()
    math_canvas = CV_Math2DCanvas(0.0 + 1.0im, 1.0 + 0.0im, 220)

    fill_painter = CV_FillPainter()
    styled_painter = cv_color(0.7, 0.4, 0.4) ↦ fill_painter

    cv_create_context(math_canvas) do canvas_context
        cv_paint(canvas_context, styled_painter)
    end

    return math_canvas
end

"""
## `doc: CV_ValueMarkPainter`

## Example for `CV_ValueMarkPainter`

![./Painter_markpainter.png]({image_from_canvas: example_mark_painter()})

```julia
{func: example_mark_painter}
```
"""
help_mark_painter() = nothing

function example_mark_painter()
    math_canvas = CV_Math2DCanvas(0.0 + 1.0im, 1.0 + 0.0im, 220)

    x_pos, y_pos = CV_TranslateByOffset(Float64), CV_TranslateByOffset(Float64)
    x_pos.value, y_pos.value = 0.7, 0.3

    bg_fill = cv_white ↦ CV_FillPainter()  # for background
    grid_style = cv_color(0.8, 0.8, 0.8) → cv_linewidth(1)
    grid = grid_style ↦ CV_GridPainter(0.0:0.1:1.0, 0.0:0.1:1.0)
    horiz_mark = CV_ValueMarkPainter(x_pos, 0.5, 0.1, true)
    vert_mark = CV_ValueMarkPainter(y_pos, 0.2, 0.1, false)

    h_painter = (cv_color(1,0,0) → cv_linewidth(2)) ↦ horiz_mark
    v_painter = (cv_color(0,1,0) → cv_linewidth(2)) ↦ vert_mark

    cv_create_context(math_canvas) do canvas_context
        cv_paint(canvas_context, bg_fill)
        cv_paint(canvas_context, grid)
        cv_paint(canvas_context, h_painter)
        cv_paint(canvas_context, v_painter)
    end

    return math_canvas
end

"""
## `doc: CV_GridPainter`

## Example for `CV_GridPainter`

![./Painter_gridpainter.png]({image_from_canvas: example_grid_painter()})

```julia
{func: example_grid_painter}
```
"""
help_grid_painter() = nothing

function example_grid_painter()
    math_canvas = CV_Math2DCanvas(0.0 + 1.0im, 1.0 + 0.0im, 220)

    bg_fill = cv_white ↦ CV_FillPainter()  # for background

    style1 = cv_color(0.7, 0.3, 0.3) → cv_linewidth(2)
    grid1 = style1 ↦ CV_GridPainter(0.0:0.2:1.0, 0.0:0.2:1.0)

    style2 = cv_color(0.3, 0.7, 0.3) → cv_linewidth(1)
    grid2 = style2 ↦ CV_GridPainter(0.1:0.2:0.9, 0.1:0.2:0.9)

    cv_create_context(math_canvas) do canvas_context
        cv_paint(canvas_context, bg_fill)
        cv_paint(canvas_context, grid2)
        cv_paint(canvas_context, grid1)
    end

    return math_canvas
end

"""
## `doc: CV_LinePainter`

## Example for `CV_LinePainter`

![./Painter_linepainter.png]({image_from_canvas: example_line_painter()})

```julia
{func: example_line_painter}
```
"""
help_line_painter() = nothing

function example_line_painter()
    math_canvas = CV_Math2DCanvas(-1.0 + 1.0im, 1.0 + -1.0im, 110)
    bg_fill = cv_white ↦ CV_FillPainter()  # for background

    grid_style = cv_color(0.7, 0.7, 0.7) → cv_linewidth(1)
    grid = grid_style ↦ CV_GridPainter(-1.0:0.2:1.0, -1.0:0.2:1.0)

    segment = [0.2im + exp(ϕ*2im)*ϕ/7 for ϕ in LinRange(0, 2*π, 200)]
    style = cv_color(0, 0, 1) → cv_linewidth(2)

    seg_painter = style ↦ CV_LinePainter([segment])

    cv_create_context(math_canvas) do canvas_context
        cv_paint(canvas_context, bg_fill)
        cv_paint(canvas_context, grid)
        cv_paint(canvas_context, seg_painter)
    end

    return math_canvas
end

"""
## `doc: CV_DirectionPainter`

## Example for `CV_DirectionPainter`

![./Painter_dirpainter.png]({image_from_canvas: example_dir_painter()})

```julia
{func: example_dir_painter}
```
"""
help_dir_painter() = nothing

function example_dir_painter()
    math_canvas = CV_Math2DCanvas(-1.0 + 1.0im, 1.0 + -1.0im, 110)
    bg_fill = cv_white ↦ CV_FillPainter()  # for background

    grid_style = cv_color(0.7, 0.7, 0.7) → cv_linewidth(1)
    grid = grid_style ↦ CV_GridPainter(-1.0:0.2:1.0, -1.0:0.2:1.0)

    segment = [0.2im + exp(ϕ*2im)*ϕ/7 for ϕ in LinRange(0, 2*π, 200)]
    style = cv_color(0.8, 0.8, 1) → cv_linewidth(1)  # light blue for curve
    seg_painter = style ↦ CV_LinePainter([segment])

    dir_style = cv_color(0.9, 0, 0)
    dir_painter = dir_style ↦ CV_DirectionPainter(identity,
        [segment]; every_len=0.5, arrow=0.1*exp(1im*π*8/9))

    cv_create_context(math_canvas) do canvas_context
        cv_paint(canvas_context, bg_fill)
        cv_paint(canvas_context, grid)
        cv_paint(canvas_context, seg_painter)
        cv_paint(canvas_context, dir_painter)
    end

    return math_canvas
end

"""
## `doc: CV_PortraitPainter`

## Example for `CV_PortraitPainter`

![./Painter_portraitpainter.png]({image_from_canvas: example_portrait_painter()})

```julia
{func: example_portrait_painter}
```
"""
help_portrait_painter() = nothing

function example_portrait_painter()
    math_canvas = CV_Math2DCanvas(0.0 + 1.0im, 1.0 + 0.0im, 220)

    trafo = z -> (z - 0.6 - 0.2im)^2 + 0.15*exp(z)
    painter = CV_PortraitPainter(trafo)

    cv_create_context(math_canvas) do canvas_context
        cv_paint(canvas_context, painter)
    end

    return math_canvas
end

"""
## `doc: CV_CombiPainter`

## `doc: →:Tuple{T, S} where {T<:CV_Painter, S<:CV_Painter}`
"""
help_combi_painter() = nothing

"""
## `doc: CV_StyledPainter`

## `doc: ↦:Tuple{CV_ContextStyle, CV_Painter}`
"""
help_styled_painter() = nothing

"""
## `doc: CV_Math2DCanvasPainter`

## Example for `CV_Math2DCanvasPainter`

![./Painter_m2dcanvaspainter.png]({image_from_canvas: example_m2d_canvas_painter()})

```julia
{func: example_m2d_canvas_painter}
```
"""
help_m2d_canvas_painter() = nothing

function example_m2d_canvas_painter()
    math_canvas = CV_Math2DCanvas(0.0 + 1.0im, 1.0 + 0.0im, 220)

    bg_fill = cv_white ↦ CV_FillPainter()
    grid_style = cv_color(0.8, 0.8, 0.8) → cv_linewidth(1)
    grid = grid_style ↦ CV_GridPainter(0.0:0.1:1.0, 0.0:0.1:1.0)

    can_letter = cv_example_image_letter(; letter="E",
        canvas=CV_Math2DCanvas(-0.5 + 0.5im, 0.5 -0.5im, 400))

    trafo = z -> (0.8*z - 0.6 - 0.2im)^2 + 0.10*exp(0.8*z) + 0.2im
    painter = CV_Math2DCanvasPainter(trafo, can_letter)

    cv_create_context(math_canvas) do canvas_context
        cv_paint(canvas_context, bg_fill)
        cv_paint(canvas_context, grid)
        cv_paint(canvas_context, painter)
    end

    return math_canvas
end

"""
## `doc: CV_WindingPainter`

## Example for `CV_WindingPainter`

![./Painter_windingpainter.png]({image_from_canvas: example_winding_painter()})

```julia
{func: example_winding_painter}
```
"""
help_winding_painter() = nothing

function example_winding_painter()
    math_canvas = CV_Math2DCanvas(0.0 + 1.0im, 1.0 + 0.0im, 220)

    trafo = z -> (0.7*z)^2 + 0.5 + 0.5im
    curve = cv_arc_lines(0, 2π, (0.8,))[1][1:end-1] .+ 0.1 .+ 0.1im

    bg_fill = cv_white ↦ CV_FillPainter()
    grid_style = cv_color(0.8, 0.8, 0.8) → cv_linewidth(1)
    grid = grid_style ↦ CV_GridPainter(0.0:0.1:1.0, 0.0:0.1:1.0)

    painter = CV_WindingPainter(trafo, [curve])

    cv_create_context(math_canvas) do canvas_context
        cv_paint(canvas_context, bg_fill)
        cv_paint(canvas_context, grid)
        cv_paint(canvas_context, painter)
    end

    return math_canvas
end

"""
## `doc: cv_parallel_lines`

## Example for `cv_parallel_lines`

![./Painter_parallel_lines.png]({image_from_canvas: example_parallel_lines()})

```julia
{func: example_parallel_lines}
```
"""
help_parallel_lines() = nothing

function example_parallel_lines()
    math_canvas = CV_Math2DCanvas(0.0 + 1.0im, 1.0 + 0.0im, 220)

    bg_fill = cv_white ↦ CV_FillPainter()
    grid_style = cv_color(0.8, 0.8, 0.8) → cv_linewidth(1)
    grid = grid_style ↦ CV_GridPainter(0.0:0.1:1.0, 0.0:0.1:1.0)

    trafo_red = z -> z + 0.4 + 0.4im
    red_lines = cv_parallel_lines(0.2 + 0.2im; width=0.5)
    red_style = cv_color(1, 0, 0) → cv_linewidth(2) → cv_antialias_best
    p_red = red_style ↦ CV_LinePainter(trafo_red, red_lines)

    trafo_blue = z -> z + 0.8 + 0.8im
    blue_lines = cv_parallel_lines(-0.05 + 0.2im; width=0.2, lines=10)
    blue_style = cv_color(0, 0, 1) → cv_linewidth(1) → cv_antialias_best
    p_blue = blue_style ↦ CV_LinePainter(trafo_blue, blue_lines)

    cv_create_context(math_canvas) do canvas_context
        cv_paint(canvas_context, bg_fill)
        cv_paint(canvas_context, grid)
        cv_paint(canvas_context, p_red)
        cv_paint(canvas_context, p_blue)
    end

    return math_canvas
end

"""
## `doc: cv_arc_lines`

## Example for `cv_arc_lines`

![./Painter_arc_lines.png]({image_from_canvas: example_arc_lines()})

```julia
{func: example_arc_lines}
```
"""
help_arc_lines() = nothing

function example_arc_lines()
    math_canvas = CV_Math2DCanvas(0.0 + 1.0im, 1.0 + 0.0im, 220)

    bg_fill = cv_white ↦ CV_FillPainter()
    grid_style = cv_color(0.8, 0.8, 0.8) → cv_linewidth(1)
    grid = grid_style ↦ CV_GridPainter(0.0:0.1:1.0, 0.0:0.1:1.0)

    trafo_red = z -> z + 0.4 + 0.4im
    red_lines = cv_arc_lines(-π/4, 3π/2, (0.1, 0.2, 0.3, 0.4))
    red_style = cv_color(1, 0, 0) → cv_linewidth(2) → cv_antialias_best
    p_red = red_style ↦ CV_LinePainter(trafo_red, red_lines)

    trafo_blue = z -> z + 0.8 + 0.8im
    blue_lines = cv_arc_lines(π, 2π, tuple(0.03:0.03:0.15...))
    blue_style = cv_color(0, 0, 1) → cv_linewidth(1) → cv_antialias_best
    p_blue = blue_style ↦ CV_LinePainter(trafo_blue, blue_lines)

    cv_create_context(math_canvas) do canvas_context
        cv_paint(canvas_context, bg_fill)
        cv_paint(canvas_context, grid)
        cv_paint(canvas_context, p_red)
        cv_paint(canvas_context, p_blue)
    end

    return math_canvas
end

"""
## `doc: cv_star_lines`

## Example for `cv_star_lines`

![./Painter_star_lines.png]({image_from_canvas: example_star_lines()})

```julia
{func: example_star_lines}
```
"""
help_star_lines() = nothing

function example_star_lines()
    math_canvas = CV_Math2DCanvas(0.0 + 1.0im, 1.0 + 0.0im, 220)

    bg_fill = cv_white ↦ CV_FillPainter()
    grid_style = cv_color(0.8, 0.8, 0.8) → cv_linewidth(1)
    grid = grid_style ↦ CV_GridPainter(0.0:0.1:1.0, 0.0:0.1:1.0)

    trafo_red = z -> z + 0.4 + 0.4im
    red_lines = cv_star_lines(0.05, 0.4, tuple(0:π/4:2π...))
    red_style = cv_color(1, 0, 0) → cv_linewidth(2) → cv_antialias_best
    p_red = red_style ↦ CV_LinePainter(trafo_red, red_lines)

    trafo_blue = z -> z + 0.8 + 0.8im
    blue_lines = cv_star_lines(0, 0.15, tuple(π:π/8:2π...))
    blue_style = cv_color(0, 0, 1) → cv_linewidth(1) → cv_antialias_best
    p_blue = blue_style ↦ CV_LinePainter(trafo_blue, blue_lines)

    cv_create_context(math_canvas) do canvas_context
        cv_paint(canvas_context, bg_fill)
        cv_paint(canvas_context, grid)
        cv_paint(canvas_context, p_red)
        cv_paint(canvas_context, p_blue)
    end

    return math_canvas
end

"""
## `doc: cv_star_arc_lines`
"""
help_star_arc_lines() = nothing

create_icon(can) = create_doc_icon(can, cv_rect_blwh(Int32, 50, 70, 70, 70))

function get_doc_icon()
    math_canvas = CV_Math2DCanvas(0.0 + 1.0im, 1.0 + 0.0im, 220)

    trafo = z -> (z - 0.6 - 0.2im)^2 + 0.15*exp(z)
    ppainter = CV_PortraitPainter(trafo)

    arc_base = cv_arc_lines(-π, π/2, (0.4,))
    arc = [arc_base[1] .+ 0.4 .+ 0.6im]

    arc_style = cv_white → cv_linewidth(4)
    arc_painter = arc_style ↦ CV_LinePainter(arc)
    dir_painter = arc_style ↦ CV_DirectionPainter(identity, arc;
        every_len=0.5, arrow=0.15*exp(1im*π*8/9))

    cv_create_context(math_canvas) do canvas_context
        cv_paint(canvas_context, ppainter)
        cv_paint(canvas_context, arc_painter)
        cv_paint(canvas_context, dir_painter)
    end

    icon = create_doc_icon(math_canvas, cv_rect_blwh(Int32, 10, 10, 200, 200))
    return icon
end

function create_document(doc_env::DocCreationEnvironment)
    doc_source = DocSource("Painter", @__MODULE__)
    context = DocContext(doc_env, doc_source)

    md = Markdown.MD()
    for part in (painter_intro, painter_trafo,
            help_fill_painter, help_mark_painter,
            help_grid_painter, help_line_painter, help_dir_painter,
            help_portrait_painter, help_combi_painter, help_styled_painter,
            help_m2d_canvas_painter, help_winding_painter,
            help_parallel_lines, help_arc_lines, help_star_lines,
            help_star_arc_lines)
        part_md = Base.Docs.doc(part)
        substitute_marker_in_markdown(context, part_md)

        append_md(md, part_md)
    end

    doc = Document(doc_source, md)
    return doc
end

end

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
using Documenter, Comonicon

makedocs(;
    modules = [Comonicon],
    format = Documenter.HTML(prettyurls = !("local" in ARGS)),
    pages = [
        "Home" => "index.md",
        "Command Types" => "types.md",
        "Command Parsing" => "parse.md",
        "Code Generation" => "codegen.md",
        "Build and Install CLI" => "build.md",
    ],
    repo = "https://github.com/Roger-luo/Comonicon.jl",
    sitename = "Comonicon.jl",
)

deploydocs(; repo = "github.com/Roger-luo/Comonicon.jl")


text:
using ADVUMPS
using ADVUMPS:qrpos,lqpos,num_grad,Ising
using Test
using Zygote
using KrylovKit
using LinearAlgebra
using OMEinsum
using ChainRulesCore
using ChainRulesTestUtils
using Random

@testset "Zygote" begin
    a = randn(10,10)
    @test Zygote.gradient(norm, a)[1] ≈ num_grad(norm, a)

    foo1 = x -> sum(Float64[x 2x; 3x 4x])
    @test Zygote.gradient(foo1, 1)[1] ≈ num_grad(foo1, 1)

    function foo3(x)
        x = 2*x + x
        return x
    end
    @test Zygote.gradient(foo3, 1)[1] ≈ num_grad(foo3, 1)
end

@testset "Zygote.@ignore" begin
    function foo2(x)
        return x^2
    end
    function foo3(x)
        return x^2 + Zygote.@ignore x^3
    end
    @test foo2(1) != foo3(1)
    @test Zygote.gradient(foo2,1)[1] ≈ Zygote.gradient(foo3,1)[1]
end

@testset "QR factorization" begin
    M = rand(10,10)
    function foo5(x)
        A = M.*x
        Q, R = qrpos(A)
        return norm(Q) + norm(R)
    end
    @test isapprox(Zygote.gradient(foo5, 1)[1], num_grad(foo5, 1), atol = 1e-5)
end

@testset "LQ factorization" begin
    M = rand(10,10)
    function foo6(x)
        A = M .*x
        L, Q = lqpos(A)
        return norm(L) + norm(Q)
    end
    @test isapprox(Zygote.gradient(foo6, 1)[1], num_grad(foo6, 1), atol = 1e-5)
end

@testset "eigsolve and linsolve" begin
    Random.seed!(100)
    N = 1000
    A = rand(N,N)
    λLs, Ls, info = eigsolve(L -> ein"a,ab -> b"(L,A), rand(N), 2, :LM)
    λL = real(λLs[1])
    L = real(Ls[1])
    # @show λLs[1] λLs[2]
    λRs, Rs, info = eigsolve(R -> ein"ab,b -> a"(A,R), rand(N), 2, :LM)
    λR = real(λRs[1])
    R = real(Rs[1])
    @test λL ≈ λR
    @test ein"a,a -> "(L,L)[] ≈ 1
    @test ein"a,a -> "(R,R)[] ≈ 1

    dL = rand(N)
    dL -= ein"a,a -> "(L,dL)[] .* L
    @test isapprox(ein"a,a -> "(L,dL)[], 0, atol = 1e-9)
    ξL, info = linsolve(R -> ein"ab,b -> a"(A,R), dL, -λL, 1)
    @test isapprox(ein"a,a -> "(ξL,L)[], 0, atol = 1e-9)

    dR = rand(N)
    dR -= ein"a,a -> "(R,dR)[] .* R
    @test isapprox(ein"a,a -> "(R,dR)[], 0, atol = 1e-9)
    ξR, info = linsolve(L -> ein"a, ab -> b"(L,A), dR, -λR, 1)
    @test isapprox(ein"a,a -> "(ξR,R)[], 0, atol = 1e-9)
end


text:
#Look up Table
include("lut/mc.jl")

"""
`marching_cubes(sdf::SignedDistanceField, [iso = 0.0,] [MT = HomogenousMesh{Point{3,Float64},Face{3,Int}}])`

Construct a `HomogenousMesh` from a `SignedDistanceField` using the
marching cubes algorithm. This method is faster than Marching Tetrahedra
and generates fewer vertices and faces (about 1/4 as many).
However it may generate non-manifold meshes, while Marching
Tetrahedra guarentees a manifold mesh.
"""
function marching_cubes(sdf::SignedDistanceField{3,ST,FT},
                               iso=0.0,
                               MT::Type{M}=SimpleMesh{Point{3,Float64},Face{3,Int}},
                               eps=0.00001, reduceverts=Val(true)) where {ST,FT,M<:AbstractMesh}
    nx, ny, nz = size(sdf)
    h = HyperRectangle(sdf)
    w = widths(h)
    orig = origin(HyperRectangle(sdf))

    # we subtract one from the length along each axis because
    # an NxNxN SDF has N-1 cells on each axis
    s = Point{3,Float64}(w[1]/(nx-1), w[2]/(ny-1), w[3]/(nz-1))

    # arrays for vertices and faces
    vts = Point{3,Float64}[]
    fcs = Face{3,Int}[]
    mt = max(nx,ny,nz)
    typeof(reduceverts) == Val{true} && sizehint!(vts, mt*mt*5)
    typeof(reduceverts) == Val{false} && sizehint!(vts, mt*mt*6)
    sizehint!(fcs, mt*mt*2)
    vertlist = Vector{Point{3,Float64}}(undef, 12)
    @inbounds for xi = 1:nx-1, yi = 1:ny-1, zi = 1:nz-1

        iso_vals = (sdf[xi,yi,zi],
                    sdf[xi+1,yi,zi],
                    sdf[xi+1,yi+1,zi],
                    sdf[xi,yi+1,zi],
                    sdf[xi,yi,zi+1],
                    sdf[xi+1,yi,zi+1],
                    sdf[xi+1,yi+1,zi+1],
                    sdf[xi,yi+1,zi+1])

        #Determine the index into the edge table which
        #tells us which vertices are inside of the surface
        cubeindex = _get_cubeindex(iso_vals, iso)

        # Cube is entirely in/out of the surface
        (cubeindex == 0x00 || cubeindex == 0xff) && continue

        points = (Point{3,Float64}(xi-1,yi-1,zi-1) .* s .+ orig,
                  Point{3,Float64}(xi,yi-1,zi-1) .* s .+ orig,
                  Point{3,Float64}(xi,yi,zi-1) .* s .+ orig,
                  Point{3,Float64}(xi-1,yi,zi-1) .* s .+ orig,
                  Point{3,Float64}(xi-1,yi-1,zi) .* s .+ orig,
                  Point{3,Float64}(xi,yi-1,zi) .* s .+ orig,
                  Point{3,Float64}(xi,yi,zi) .* s .+ orig,
                  Point{3,Float64}(xi-1,yi,zi) .* s .+ orig)

        # Find the vertices where the surface intersects the cube
        find_vertices_interp!(vertlist, points, iso_vals, cubeindex, iso, eps)

        # Create the triangle
        typeof(reduceverts) == Val{true} && _mc_unique_triangles!(vts, fcs, vertlist, cubeindex)
        typeof(reduceverts) == Val{false} && _mc_create_triangles!(vts, fcs, vertlist, cubeindex)
    end
    MT(vts,fcs)
end

function marching_cubes(f::Function,
                        bounds::HyperRectangle,
                        samples::NTuple{3,Int}=(256,256,256),
                        iso=0.0,
                        MT::Type{M}=SimpleMesh{Point{3,Float64},Face{3,Int}},
                        eps=0.00001, reduceverts=Val(true)) where {ST,FT,M<:AbstractMesh}
    nx, ny, nz = samples[1], samples[2], samples[3]
    w = widths(bounds)
    orig = origin(bounds)

    # we subtract one from the length along each axis because
    # an NxNxN SDF has N-1 cells on each axis
    s = Point{3,Float64}(w[1]/(nx-1), w[2]/(ny-1), w[3]/(nz-1))

    # arrays for vertices and faces
    vts = Point{3,Float64}[]
    fcs = Face{3,Int}[]
    mt = max(nx,ny,nz)
    typeof(reduceverts) == Val{true} && sizehint!(vts, mt*mt*5)
    typeof(reduceverts) == Val{false} && sizehint!(vts, mt*mt*6)
    sizehint!(fcs, mt*mt*2)
    vertlist = Vector{Point{3,Float64}}(undef, 12)
    iso_vals = Vector{Float64}(undef,8)
    @inbounds for xi = 1:nx-1, yi = 1:ny-1, zi = 1:nz-1

        points = (Point{3,Float64}(xi-1,yi-1,zi-1) .* s .+ orig,
                  Point{3,Float64}(xi,yi-1,zi-1) .* s .+ orig,
                  Point{3,Float64}(xi,yi,zi-1) .* s .+ orig,
                  Point{3,Float64}(xi-1,yi,zi-1) .* s .+ orig,
                  Point{3,Float64}(xi-1,yi-1,zi) .* s .+ orig,
                  Point{3,Float64}(xi,yi-1,zi) .* s .+ orig,
                  Point{3,Float64}(xi,yi,zi) .* s .+ orig,
                  Point{3,Float64}(xi-1,yi,zi) .* s .+ orig)

        if zi == 1
            for i = 1:8
                iso_vals[i] = f(points[i])
            end
        else
            iso_vals[1] = iso_vals[5]
            iso_vals[2] = iso_vals[6]
            iso_vals[3] = iso_vals[7]
            iso_vals[4] = iso_vals[8]
            iso_vals[5] = f(points[5])
            iso_vals[6] = f(points[6])
            iso_vals[7] = f(points[7])
            iso_vals[8] = f(points[8])
        end

        #Determine the index into the edge table which
        #tells us which vertices are inside of the surface
        cubeindex = _get_cubeindex(iso_vals, iso)

        # Cube is entirely in/out of the surface
        (cubeindex == 0x00 || cubeindex == 0xff) && continue

        # Find the vertices where the surface intersects the cube
        # TODO this can use the underlying function to find the zero.
        # The underlying space is non-linear so there will be error otherwise
        find_vertices_interp!(vertlist, points, iso_vals, cubeindex, iso, eps)

        # Create the triangle
        typeof(reduceverts) == Val{true} && _mc_unique_triangles!(vts, fcs, vertlist, cubeindex)
        typeof(reduceverts) == Val{false} && _mc_create_triangles!(vts, fcs, vertlist, cubeindex)

    end
    MT(vts,fcs)
end

@inline function _mc_create_triangles!(vts, fcs, vertlist, cubeindex)
    fct = length(vts) + 3

    push!(vts, vertlist[tri_table[cubeindex][1]],
               vertlist[tri_table[cubeindex][2]],
               vertlist[tri_table[cubeindex][3]])
    push!(fcs, Face{3,Int}(fct, fct-1, fct-2))

    iszero(tri_table[cubeindex][4]) && return
    fct += 3
    push!(vts, vertlist[tri_table[cubeindex][4]],
               vertlist[tri_table[cubeindex][5]],
               vertlist[tri_table[cubeindex][6]])
    push!(fcs, Face{3,Int}(fct, fct-1, fct-2))

    iszero(tri_table[cubeindex][7]) && return
    fct += 3
    push!(vts, vertlist[tri_table[cubeindex][7]],
               vertlist[tri_table[cubeindex][8]],
               vertlist[tri_table[cubeindex][9]])
    push!(fcs, Face{3,Int}(fct, fct-1, fct-2))

    iszero(tri_table[cubeindex][10]) && return
    fct += 3
    push!(vts, vertlist[tri_table[cubeindex][10]],
               vertlist[tri_table[cubeindex][11]],
               vertlist[tri_table[cubeindex][12]])
    push!(fcs, Face{3,Int}(fct, fct-1, fct-2))

    iszero(tri_table[cubeindex][13]) && return
    fct += 3
    push!(vts, vertlist[tri_table[cubeindex][13]],
               vertlist[tri_table[cubeindex][14]],
               vertlist[tri_table[cubeindex][15]])
    push!(fcs, Face{3,Int}(fct, fct-1, fct-2))
end

"""
Create triangles by only adding unique vertices within the voxel.
Each face may share a reference to a vertex with another face.
"""
@inline function _mc_unique_triangles!(vts, fcs, vertlist, cubeindex)
    fct = length(vts)

    vert_to_add = _mc_verts[cubeindex]
    # Each vertex list will have atleast 3 elements so we can
    # add them to the list immediately
    push!(vts, vertlist[vert_to_add[1]])
    push!(vts, vertlist[vert_to_add[2]])
    push!(vts, vertlist[vert_to_add[3]])

    for i = 4:12
        elt = vert_to_add[i]
        iszero(elt) && break
        push!(vts, vertlist[elt])
    end
    offsets = _mc_connectivity[cubeindex]

    # There is atleast one face so we can push it immediately
    push!(fcs, Face{3,Int}(fct+offsets[3], fct+offsets[2], fct+offsets[1]))

    iszero(offsets[4]) && return
    push!(fcs, Face{3,Int}(fct+offsets[6], fct+offsets[5], fct+offsets[4]))

    iszero(offsets[7]) && return
    push!(fcs, Face{3,Int}(fct+offsets[9], fct+offsets[8], fct+offsets[7]))

    iszero(offsets[10]) && return
    push!(fcs, Face{3,Int}(fct+offsets[12], fct+offsets[11], fct+offsets[10]))

    iszero(offsets[13]) && return
    push!(fcs, Face{3,Int}(fct+offsets[15], fct+offsets[14], fct+offsets[13]))

end

@inline function find_vertices_interp!(vertlist, points, iso_vals, cubeindex, iso, eps)
    if !iszero(edge_table[cubeindex] & 0x001)
        vertlist[1] =
            vertex_interp(iso,points[1],points[2],iso_vals[1],iso_vals[2], eps)
    end
    if !iszero(edge_table[cubeindex] & 0x002)
        vertlist[2] =
            vertex_interp(iso,points[2],points[3],iso_vals[2],iso_vals[3], eps)
    end
    if !iszero(edge_table[cubeindex] & 0x004)
        vertlist[3] =
            vertex_interp(iso,points[3],points[4],iso_vals[3],iso_vals[4], eps)
    end
    if !iszero(edge_table[cubeindex] & 0x008)
        vertlist[4] =
            vertex_interp(iso,points[4],points[1],iso_vals[4],iso_vals[1], eps)
    end
    if !iszero(edge_table[cubeindex] & 0x010)
        vertlist[5] =
            vertex_interp(iso,points[5],points[6],iso_vals[5],iso_vals[6], eps)
    end
    if !iszero(edge_table[cubeindex] & 0x020)
        vertlist[6] =
            vertex_interp(iso,points[6],points[7],iso_vals[6],iso_vals[7], eps)
    end
    if !iszero(edge_table[cubeindex] & 0x040)
        vertlist[7] =
            vertex_interp(iso,points[7],points[8],iso_vals[7],iso_vals[8], eps)
    end
    if !iszero(edge_table[cubeindex] & 0x080)
        vertlist[8] =
            vertex_interp(iso,points[8],points[5],iso_vals[8],iso_vals[5], eps)
    end
    if !iszero(edge_table[cubeindex] & 0x100)
        vertlist[9] =
            vertex_interp(iso,points[1],points[5],iso_vals[1],iso_vals[5], eps)
    end
    if !iszero(edge_table[cubeindex] & 0x200)
        vertlist[10] =
            vertex_interp(iso,points[2],points[6],iso_vals[2],iso_vals[6], eps)
    end
    if !iszero(edge_table[cubeindex] & 0x400)
        vertlist[11] =
            vertex_interp(iso,points[3],points[7],iso_vals[3],iso_vals[7], eps)
    end
    if !iszero(edge_table[cubeindex] & 0x800)
        vertlist[12] =
            vertex_interp(iso,points[4],points[8],iso_vals[4],iso_vals[8], eps)
    end
end

# Linearly interpolate the position where an isosurface cuts
# an edge between two vertices, each with their own scalar value
function vertex_interp(iso, p1, p2, valp1, valp2, eps = 0.00001)

    abs(iso - valp1) < eps && return p1
    abs(iso - valp2) < eps && return p2
    abs(valp1-valp2) < eps && return p1
    mu = (iso - valp1) / (valp2 - valp1)
    p = p1 + mu * (p2 - p1)

    return p
end

struct MarchingCubes{T} <: AbstractMeshingAlgorithm
    iso::T
    eps::T
    reduceverts::Bool
end

MarchingCubes(;iso::T1=0.0, eps::T2=1e-3, reduceverts::Bool=true) where {T1, T2} = MarchingCubes{promote_type(T1, T2)}(iso, eps, reduceverts)
MarchingCubes(iso) = MarchingCubes(iso=iso)
MarchingCubes(iso,eps) = MarchingCubes(iso=iso,eps=eps)

function (::Type{MT})(df::SignedDistanceField, method::MarchingCubes)::MT where {MT <: AbstractMesh}
    marching_cubes(df, method.iso, MT, method.eps, Val(method.reduceverts))
end

function (::Type{MT})(f::Function, h::HyperRectangle, size::NTuple{3,Int}, method::MarchingCubes)::MT where {MT <: AbstractMesh}
    marching_cubes(f, h, size, method.iso, MT, method.eps, Val(method.reduceverts))
end

function (::Type{MT})(f::Function, h::HyperRectangle, method::MarchingCubes; size::NTuple{3,Int}=(128,128,128))::MT where {MT <: AbstractMesh}
    marching_cubes(f, h, size, method.iso, MT, method.eps, Val(method.reduceverts))
end