text:
using Test
using LinearAlgebra
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

@testset "checkStruct" begin
	square=[([[0.575;-0.175] [0.575;0.175] [0.925;-0.175] [0.925;0.175]],
	[[0,1,2,3]])]
	@test Lar.checkStruct(square)==length(square[1][1][:,1])
	@test Lar.checkStruct(square)==2
end

@testset "2D" begin
	square = ([[0.; 0] [0; 1] [1; 0] [1; 1]], [[1, 2, 3,
			4]], [[1,2], [1,3], [2,4], [3,4]])
	@testset "apply Translation 2D" begin
		@test typeof(Lar.apply(Lar.t(-0.5,-0.5),square))==Tuple{Array{Float64,2},Array{Array{Int64,1},1},Array{Array{Int64,1},1}}
		@test Lar.apply(Lar.t(-0.5,-0.5),square)==([-0.5 -0.5 0.5 0.5; -0.5 0.5 -0.5 0.5], Array{Int64,1}[[1, 2, 3, 4]], Array{Int64,1}[[1, 2], [1, 3], [2, 4], [3, 4]])
	end
	@testset "apply Scaling 2D" begin
		@test typeof(Lar.apply(Lar.s(-0.5,-0.5),square))==Tuple{Array{Float64,2},Array{Array{Int64,1},1},Array{Array{Int64,1},1}}
		@test Lar.apply(Lar.s(-0.5,-0.5),square)==([0.0 0.0 -0.5 -0.5; 0.0 -0.5 0.0 -0.5], Array{Int64,1}[[1, 2, 3, 4]], Array{Int64,1}[[1, 2], [1, 3], [2, 4], [3, 4]])
	end
	@testset "apply Rotation 2D" begin
		@test typeof(Lar.apply(Lar.r(0),square))==Tuple{Array{Float64,2},Array{Array{Int64,1},1},Array{Array{Int64,1},1}}
		@test Lar.apply(Lar.r(0),square)==([0.0 0.0 1.0 1.0; 0.0 1.0 0.0 1.0], Array{Int64,1}[[1, 2, 3, 4]], Array{Int64,1}[[1, 2], [1, 3], [2, 4], [3, 4]])
	end
end

@testset "3D" begin
	cube=Lar.cuboid([1,1,1])
	@testset "apply Translation 3D" begin
		@test typeof(Lar.apply(Lar.t(-0.5,-0.5,-0.5),cube))==Tuple{Array{Float64,2},Array{Array{Int64,1},1}}
		@test Lar.apply(Lar.t(-0.5, -0.5, -0.5),cube) ==
		([-0.5 -0.5 -0.5 -0.5 0.5 0.5 0.5 0.5;
		-0.5 -0.5 0.5 0.5 -0.5 -0.5 0.5 0.5; -0.5 0.5 -0.5 0.5 -0.5 0.5 -0.5 0.5],
		Array{Int64,1}[[1, 2, 3, 4, 5, 6, 7, 8]])
	end
	@testset "apply Scaling 3D" begin
		@test typeof(Lar.apply(Lar.s(-0.5,-0.5,-0.5),cube))==
		Tuple{Array{Float64,2},Array{Array{Int64,1},1}}
		@test Lar.apply(Lar.s(-0.5, -0.5, -0.5),cube) ==
		([0.0 0.0 0.0 0.0 -0.5 -0.5 -0.5 -0.5;
		0.0 0.0 -0.5 -0.5 0.0 0.0 -0.5 -0.5; 0.0 -0.5 0.0 -0.5 0.0 -0.5 0.0 -0.5],
		Array{Int64,1}[[1, 2, 3, 4, 5, 6, 7, 8]])
	end
	@testset "apply Rotation 3D" begin
		@test typeof(Lar.apply(Lar.r(pi,0,0),cube))==
		Tuple{Array{Float64,2},Array{Array{Int64,1},1}}
		@test Lar.apply(Lar.r(0, 0, 0),cube) ==
		 ([0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0; 0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0;
		 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0], Array{Int64,1}[[1, 2, 3, 4, 5, 6, 7, 8]])
	end
end

@testset "cuboid Tests" begin
	square=Lar.cuboid([1,1])
	cube=Lar.cuboid([1,1,1])
	@testset "cuboid Tests 2D" begin
		@test typeof(square)==Tuple{Array{Float64,2},Array{Array{Int64,1},1}}
		@test length(square)==2
		@test size(square[1])==(2, 4)
	end
	@testset "cuboid Tests 3D" begin
		@test typeof(cube)==Tuple{Array{Float64,2},Array{Array{Int64,1},1}}
		@test length(cube)==2
		@test size(cube[1])==(3, 8)
	end
end

@testset "traversal" begin
	square=Lar.cuboid([1,1])
	dim = 2
	structure=Lar.Struct([square])
	@test typeof(structure.body) ==
	Array{Tuple{Array{Float64,2},Array{Array{Int64,1},1}},1}
	@test length(Lar.traversal(Matrix{Float64}(LinearAlgebra.I, dim+1, dim+1),[],structure,[]))==length(structure.body)
	@test typeof(Lar.traversal(Matrix{Float64}(LinearAlgebra.I, dim+1, dim+1),[],structure,[]))==Array{Any,1}
end

@testset "Struct Tests" begin
	square=Lar.cuboid([1,1])
	@test Lar.Struct([square]).body==[square]
	@test Lar.Struct([square]).dim==size(square[1],1)
	@test Lar.Struct([square]).body[1]==square
end

# to test
#@testset "embedTraversal Tests" begin
#	square = Lar.cuboid([1,1])
#	x = Lar.Struct([square])
#	@test length(Lar.embedTraversal(Lar.Struct(),x,1,"New").body[2][1][1])==
#	length(x.body[1][1][1])+1
#	#in this case n=1, but generally:
#	# length(length(embedTraversal(x,x,1,"New")=length(x.body[1][1][1])+n
#	@test length(embedTraversal(deepcopy(x),deepcopy(x),3,"New").body[2][1][1])?==
#	length(x.body[1][1][1])+3
#	@test typeof(embedTraversal(deepcopy(x),deepcopy(x),1,"New"))==Struct
#end

#@testset "embedStruct Tests" begin
#	square = Lar.cuboid([1,1])
#	x = Lar.Struct([square])
#	@test length(Lar.embedStruct(1)(x).body[1][1][1])==length(x.body[1][1][1])+1
#	#in this case n = 1, but generally:
#	#length(embedStruct(n)(x).body[1][1][1])=length(x.body[1][1][1])+n
#	@test length(Lar.embedStruct(3)(x).body[1][1][1])==length(x.body[1][1][1])+3
#	@test typeof(Lar.embedStruct(1)(x))==Struct
#end

@testset "removeDups Tests" begin
	CW1=[[0,1,2,3],[4,5,6,7],[0,1,4,5],[2,3,6,7],[0,2,4,6],[1,3,5,7],
		 [4,5,6,7],[8,9,10,11],[4,5,8,9],[6,7,10,11],[4,6,8,10],[5,7,9,11]]
	CW2=[[0,1,2,3],[4,5,6,7],[8,9,10,11],[12,13,14,15],[16,17,18,19]]
	@testset "removeDups 3D" begin
		@test length(Lar.removeDups(CW1))<= length(CW1)
		@test typeof(Lar.removeDups(CW1))==Array{Array{Int64,1},1}
	end
	@testset "removeDups 2D" begin
		@test length(Lar.removeDups(CW2))<= length(CW2)
		@test typeof(Lar.removeDups(CW2))==Array{Array{Int64,1},1}
	end
end

@testset "struct2lar" begin
   @testset "struct2lar 2D" begin
      square = Lar.cuboid([1,1])
      table = Lar.apply(Lar.t(-0.5,-0.5),square)
      structure = Lar.Struct([repeat([table,Lar.r(pi/2)],outer=2)...])
      @test typeof(Lar.struct2lar(structure))==
      		Tuple{Array{Float64,2},Array{Array{Int64,1},1}}
	  @test size(Lar.struct2lar(structure)[1])==(2, 4)
   end
   @testset "struct2lar 3D" begin
      BV = [[1,2,3,4],[5,6,7,8],[1,2,5,6],[3,4,7,8],[1,3,5,7],[2,4,6,8]]
      V = hcat([[0.0,0.0,0.0],[0.0,0.0,1.0],[0.0,1.0,0.0],[0.0,1.0,1.0],
      [1.0,0.0,0.0],[1.0,0.0,1.0],[1.0,1.0,0.0],[1.0,1.0,1.0]]...)
      block = (V,BV)
      structure = Lar.Struct(repeat([block,Lar.t(1,0,0)],outer=2));
      @test typeof(Lar.struct2lar(structure))==
      	Tuple{Array{Float64,2},Array{Array{Int64,1},1}}
      @test size(Lar.struct2lar(structure)[1],1)==3
    end
end


text:
const SpectrumOrProperties = Union{Spectrum, Dict{Symbol,Any}}

function spectrum(sop::SpectrumOrProperties)::Spectrum
    @assert sop isa Spectrum "The spectrum is not available."
    sop # Fails fast
end
properties(sop::SpectrumOrProperties)::Dict{Symbol, Any} = sop isa Spectrum ? properties(sop) : sop
hasspectrum(sop::SpectrumOrProperties) = sop isa Spectrum

"""
    ReferenceLabel

A label associated with reference spectra.  The label encapsulates the original spectrum and the range of channels
represented by this reference object.  structs that extend ReferenceLabel should have `.roi`, `.spectrum` and ".hash"
fields.
"""
abstract type ReferenceLabel <: Label
    #spectrum::SpectrumOrProperties
    #roi::UnitRange{Int}
    #xrays::Vector{CharXRay}
    #hash::UInt
end

"""
   channels(rl::ReferenceLabel)::UnitRange{Int}

The range of channels associated with the specified ReferenceLabel.
"""
channels(rl::ReferenceLabel) = rl.roi
Base.hash(rl::ReferenceLabel, h::UInt) = hash(rl.hash, h)
Base.isequal(c1::ReferenceLabel, c2::ReferenceLabel) =
    (typeof(c1)==typeof(c2)) &&
    (c1.hash == c2.hash) &&
    isequal(c1.roi, c2.roi) &&
    isequal(c1.xrays, c2.xrays) &&
    isequal(c1.spectrum, c2.spectrum)
xrays(cl::ReferenceLabel) = cl.xrays
spectrum(cl::ReferenceLabel)::Spectrum = spectrum(cl.spectrum)
properties(cl::ReferenceLabel) = properties(cl.spectrum)
hasspectrum(cl::ReferenceLabel) = hasspectrum(cl.spectrum)
composition(cl::ReferenceLabel) = get(properties(cl), :Composition, missing)
NeXLCore.element(cl::ReferenceLabel) = element(cl.xrays[1])
Base.isless(rl1::ReferenceLabel, rl2::ReferenceLabel) =
    return isequal(rl1.roi, rl2.roi) ? isless(rl1.spectrum[:Name], rl2.spectrum[:Name]) :
           (
        isequal(rl1.roi.start, rl2.roi.start) ? isless(rl1.roi.stop, rl2.roi.stop) :
        isless(rl1.roi.start, rl2.roi.start)
    )

"""
    CharXRayLabel

A ReferenceLabel that represents a reference spectrum or reference properties associated with a set of
characteristic x-rays (CharXRay) objects over a contiguous range of spectrum channels.
"""
struct CharXRayLabel <: ReferenceLabel
    spectrum::SpectrumOrProperties # The spectrum used as the reference for fitting...
    roi::UnitRange{Int}
    xrays::Vector{CharXRay}
    hash::UInt
    function CharXRayLabel(spec::SpectrumOrProperties, roi::UnitRange{Int}, xrays::Vector{CharXRay})
        @assert all(xr -> element(xr) == element(xrays[1]), xrays)
        new(spec, roi, xrays, hash(spec, hash(roi, hash(xrays))))
    end
end

function Base.show(io::IO, cl::CharXRayLabel)
    comp = composition(cl)
    compname = isnothing(comp) ? "Unspecified" : name(comp)
    print(io,"k[$(name(cl.xrays)), $compname]")
end

"""
    EscapeLabel

A ReferenceLabel<:FilteredLabel that Represents a reference spectrum associated with an escape peak from a set of
characteristic x-rays (CharXRay) objects over a contiguous range of spectrum channels.
"""
struct EscapeLabel <: ReferenceLabel
    spectrum::SpectrumOrProperties
    roi::UnitRange{Int}
    xrays::Vector{EscapeArtifact}
    hash::UInt

    EscapeLabel(spc::SpectrumOrProperties, roi::UnitRange{Int}, escs::AbstractVector{EscapeArtifact}) =
        new(spc, roi, convert(Vector{EscapeArtifact}, escs), hash(spc, hash(roi, hash(escs))))
end

Base.show(io::IO, escl::EscapeLabel) = print(io, name(escl))

NeXLCore.name(escl::EscapeLabel) = "Ecs[$(name([esc.xray for esc in escl.xrays]))]"
NeXLCore.element(escl::EscapeLabel) = element(escl.xrays[1])

"""
    UnknownLabel

A Label that represents the unknown spectrum.
"""
struct UnknownLabel <: Label
    spectrum::Union{HyperSpectrum,Spectrum}
    hash::UInt

    UnknownLabel(spec::Union{HyperSpectrum,Spectrum}) = new(spec, hash(spec))
end

Base.show(io::IO, unk::UnknownLabel) = print(io, unk.spectrum[:Name])
Base.isequal(ul1::UnknownLabel, ul2::UnknownLabel) = (ul1.hash==ul2.hash) && isequal(ul1.spectrum, ul2.spectrum)
spectrum(unkl::UnknownLabel) = unkl.spectrum
properties(unkl::UnknownLabel) = properties(unkl.spectrum)

"""
    HyperSpectrumLabel

A Label that represents a single spectrum with a HyperSpectrum.
"""
struct HyperSpectrumLabel <: Label
    hyperspectrum::HyperSpectrum
    index::CartesianIndex
    HyperSpectrumLabel(hs::HyperSpectrum, idx::Int...) = new(hs, CartesianIndex(idx...))
end

Base.show(io::IO, unk::HyperSpectrumLabel) =
    print(io, unk.hyperspectrum[:Name] * "[$(unk.index.I)]")
Base.isequal(ul1::HyperSpectrumLabel, ul2::HyperSpectrumLabel) =
    (ul1.hyperspectrum === ul2.hyperspec) && isequal(ul1.index, ul2.index)
spectrum(unkl::HyperSpectrumLabel) = unkl.hyperspectrum[unkl.index]
properties(unkl::HyperSpectrumLabel) = properties(unkl.hyperspectrum[unkl.index])

ext: jl

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
function  s1(StartPoint, Volume)
    #  S1 is a function which looks for a lower pixel location in a local
    #  neighborhood
    #
    #   EndPoint = S1(StartPoint, GradientVolume, StepSize);
    #
    #    inputs :
    #        StartPoint: 2D or 3D location in vectorfield
    #        Volume: 2D or 3D matrix with values
    #
    #   outputs :
    #        EndPoint : The new location (zero if outside image)
    #
    # Function is written by D.Kroon University of Twente (Oct 2010)
    # Ported from Matlab to Julia by Lars Hellemo (May 2014)

    StartPoint=round(StartPoint)
    S=StartPoint
    for stepsize=1:3
        sxm=max(S[1]-stepsize,1)
        sym=max(S[2]-stepsize,1)
        sxp=min(S[1]+stepsize,size(Volume,1))
        syp=min(S[2]+stepsize,size(Volume,2))

        X=sxm:sxp
        Y=sym:syp
        SubVolume=Volume[sxm:sxp ; sym:syp]

        CVolume=SubVolume<Volume[S[1] ,S[2]]
        check=any(CVolume(:))
        if check
            break
        end
    end
    if check
        (temp,ind)=findmin(SubVolume(:))
        (i,j)=ind2sub(size(CVolume),ind)
        EndPoint=[X[i];Y[j]]
    else
        EndPoint=StartPoint
    end
end

# TODO: 3D
# else
#     for stepsize=1:3
#         sxm=max(S(1)-stepsize,1);
#         sym=max(S(2)-stepsize,1);
#         szm=max(S(3)-stepsize,1);

#         sxp=min(S(1)+stepsize,size(Volume,1));
#         syp=min(S(2)+stepsize,size(Volume,2));
#         szp=min(S(3)+stepsize,size(Volume,3));

#         X=sxm:sxp;
#         Y=sym:syp;
#         Z=szm:szp;
#         SubVolume=Volume(sxm:sxp,sym:syp,szm:szp);

#         CVolume=SubVolume<Volume(S(1),S(2),S(3));
#         check=any(CVolume(:));
#         if(check), break; end
#     end
#     if(check),
#         [temp,ind]=min(SubVolume(:));
#         [i,j,k]=ind2sub(size(CVolume),ind);
#         EndPoint=[X(i);Y(j);Z(k)];
#     else
#         EndPoint=StartPoint;
#     end
# end


text:
include("ctqw/ctqw.jl")
include("szegedy/szegedy.jl")


text:
# ## GPU Computing

# https://github.com/JuliaGPU

using Plots, BenchmarkTools, FFTW, LinearAlgebra

# ### Advection equation for a rotation in two dimensional domain
#md #
#md # ```math
#md #  \\frac{d f}{dt} +  (y \\frac{d f}{dx} - x \\frac{d f}{dy}) = 0
#md # ```
#md #
#md # ``x \in [-π, π], y \in [-π, π] `` and  `` t \in [0, 200π] ``

#md # ---

struct Mesh

    nx   :: Int64
    ny   :: Int64
    x    :: Vector{Float64}
    y    :: Vector{Float64}
    kx   :: Vector{Float64}
    ky   :: Vector{Float64}

    function Mesh( xmin, xmax, nx, ymin, ymax, ny)
        ## periodic boundary condition, we remove the end point.
        x = range(xmin, stop=xmax, length=nx+1)[1:end-1]
        y = range(ymin, stop=ymax, length=ny+1)[1:end-1]
        kx  = 2π ./ (xmax-xmin) .* [0:nx÷2-1;nx÷2-nx:-1]
        ky  = 2π ./ (ymax-ymin) .* [0:ny÷2-1;ny÷2-ny:-1]
        new( nx, ny, x, y, kx, ky)
    end
end

#md # ---

function exact(time, mesh :: Mesh; shift=1.0)

    f = zeros(Float64,(mesh.nx, mesh.ny))
    for (i, x) in enumerate(mesh.x), (j, y) in enumerate(mesh.y)
        xn = cos(time)*x - sin(time)*y
        yn = sin(time)*x + cos(time)*y
        f[i,j] = exp(-(xn-shift)*(xn-shift)/0.1)*exp(-(yn-shift)*(yn-shift)/0.1)
    end

    f
end

#md # ---

function animation( tf, nt)

    mesh = Mesh( -π, π, 64, -π, π, 64)
    dt = tf / nt
    t = 0
    anim = @animate for n=1:nt

       f = exact(t, mesh)
       t += dt
       p = contour(mesh.x, mesh.y, f, axis=[], framestyle=:none)
       plot!(p[1]; clims=(0.,1.), aspect_ratio=:equal, colorbar=false, show=false)
       plot!(sqrt(2) .* cos.(-pi:0.1:pi+0.1),
             sqrt(2) .* sin.(-pi:0.1:pi+0.1), label="", show=false)
       xlims!(-π,π)
       ylims!(-π,π)

    end

    anim

end

#md # ---

anim = animation( 2π, 100)
#md gif(anim, "rotation2d.gif", fps = 20)
#md nothing # hide

#md # ![](rotation2d.gif)

#md # ---

function rotation_on_cpu( mesh :: Mesh, nt :: Int64, tf :: Float64)

    dt = tf / nt

    f   = zeros(ComplexF64,(mesh.nx,mesh.ny))
    f  .= exact( 0.0, mesh )

    exky = exp.( 1im*tan(dt/2) .* mesh.x  .* mesh.ky')
    ekxy = exp.(-1im*sin(dt)   .* mesh.y' .* mesh.kx )

    for n = 1:nt

        fft!(f, 2)
        f .= exky .* f
        ifft!(f,2)

        fft!(f, 1)
        f .= ekxy .* f
        ifft!(f, 1)

        fft!(f, 2)
        f .= exky .* f
        ifft!(f, 2)

    end

    real(f)

end

#md # ---

mesh = Mesh( -π, π, 1024, -π, π, 1024)
nt, tf = 100, 20.
rotation_on_cpu(mesh, 1, 0.1) # trigger building
etime = @time norm( rotation_on_cpu(mesh, nt, tf) .- exact( tf, mesh))
println(etime)

@test true #src

#md # ---

using Pkg

GPU_ENABLED = haskey(Pkg.installed(), "CUDAdrv")

if GPU_ENABLED

    using CUDAdrv, CuArrays, CuArrays.CUFFT

    println(CUDAdrv.name(CuDevice(0)))

end

#md # ---

if GPU_ENABLED

    function rotation_on_gpu( mesh :: Mesh, nt :: Int64, tf :: Float64)

        dt  = tf / nt
        f   = zeros(ComplexF64,(mesh.nx, mesh.ny))
        f  .= exact( 0.0, mesh)

        d_f    = CuArray(f) # allocate f on GPU

        p_x    = plan_fft!(d_f,  [1]) # Create fft plans on GPU
        pinv_x = plan_ifft!(d_f, [1])
        p_y    = plan_fft!(d_f,  [2])
        pinv_y = plan_ifft!(d_f, [2])

        d_exky = CuArray(exp.( 1im*tan(dt/2) .* mesh.x  .* mesh.ky'))
        d_ekxy = CuArray(exp.(-1im*sin(dt)   .* mesh.y' .* mesh.kx ))

        for n = 1:nt

            p_y * d_f
            d_f .*= d_exky
            pinv_y * d_f

            p_x * d_f
            d_f .*= d_ekxy
            pinv_x * d_f

            p_y * d_f
            d_f .*= d_exky
            pinv_y * d_f

        end

        f .= collect(d_f) # Transfer f from GPU to CPU
        real(f)

    end

end

#md # ---

if GPU_ENABLED

    nt, tf = 100, 20.
    rotation_on_gpu(mesh, 1, 0.1)
    etime = @time norm( rotation_on_gpu(mesh, nt, tf) .- exact( tf, mesh))
    println(etime)
    @test true #src

end

#md # Next: [Metaprogrammimg and Particle Group](/Numkin2019/04/build/index.html)


text:
using ViewerGL, LinearAlgebraicRepresentation
GL = ViewerGL
Lar = LinearAlgebraicRepresentation

#Tetrahemihexahedron
V, (VV,EV,FV,CV) = Lar.simplex(3, true)
tetra = V, EV,FV,CV
model = Lar.Struct([ tetra,
					Lar.r(π,0,0), tetra,
					Lar.r(0,π,0), tetra,
					Lar.r(0,π,0), Lar.r(0,0,π), tetra ]);

V, EV,FV,CV = Lar.struct2lar(model)
GL.VIEW([ GL.GLGrid(V,FV), GL.GLFrame ]);

EV = collect(Set(EV)) # removes duplicates in EV
cop_EV = Lar.coboundary_0(EV::Lar.Cells);
cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
W = convert(Lar.Points, V');

# arrangement don't work: it loops in minimal_3cycles! TODO: fix
V, copEV, copFE, copCF = Lar.Arrangement.spatial_arrangement( W::Lar.Points, cop_EV::Lar.ChainOp, cop_FE::Lar.ChainOp)
χ = Lar.euler_characteristic(W, copEV, copFE)

FVs = [ [FV[k]] for k=1:length(FV)]
GL.VIEW(GL.GLExplode(V,FVs,1.5,1.5,1.5,99, .5));

EVs = [ [EV[k]] for k=1:length(EV)]
GL.VIEW(GL.GLExplode(V,EVs,1.2,1.2,1.2,99,1));


text:
Dict{Int32, String}(
    004 => "AFG",
    248 => "ALA",
    008 => "ALB",
    012 => "DZA",
    016 => "ASM",
    020 => "AND",
    024 => "AGO",
    660 => "AIA",
    010 => "ATA",
    028 => "ATG",
    032 => "ARG",
    051 => "ARM",
    533 => "ABW",
    036 => "AUS",
    040 => "AUT",
    031 => "AZE",
    044 => "BHS",
    048 => "BHR",
    050 => "BGD",
    052 => "BRB",
    112 => "BLR",
    056 => "BEL",
    084 => "BLZ",
    204 => "BEN",
    060 => "BMU",
    064 => "BTN",
    068 => "BOL",
    535 => "BES",
    070 => "BIH",
    072 => "BWA",
    074 => "BVT",
    076 => "BRA",
    086 => "IOT",
    096 => "BRN",
    100 => "BGR",
    854 => "BFA",
    108 => "BDI",
    132 => "CPV",
    116 => "KHM",
    120 => "CMR",
    124 => "CAN",
    136 => "CYM",
    140 => "CAF",
    148 => "TCD",
    152 => "CHL",
    156 => "CHN",
    162 => "CXR",
    166 => "CCK",
    170 => "COL",
    174 => "COM",
    178 => "COG",
    180 => "COD",
    184 => "COK",
    188 => "CRI",
    384 => "CIV",
    191 => "HRV",
    192 => "CUB",
    531 => "CUW",
    196 => "CYP",
    203 => "CZE",
    208 => "DNK",
    262 => "DJI",
    212 => "DMA",
    214 => "DOM",
    218 => "ECU",
    818 => "EGY",
    222 => "SLV",
    226 => "GNQ",
    232 => "ERI",
    233 => "EST",
    748 => "SWZ",
    231 => "ETH",
    238 => "FLK",
    234 => "FRO",
    242 => "FJI",
    246 => "FIN",
    250 => "FRA",
    254 => "GUF",
    258 => "PYF",
    260 => "ATF",
    266 => "GAB",
    270 => "GMB",
    268 => "GEO",
    276 => "DEU",
    288 => "GHA",
    292 => "GIB",
    300 => "GRC",
    304 => "GRL",
    308 => "GRD",
    312 => "GLP",
    316 => "GUM",
    320 => "GTM",
    831 => "GGY",
    324 => "GIN",
    624 => "GNB",
    328 => "GUY",
    332 => "HTI",
    334 => "HMD",
    336 => "VAT",
    340 => "HND",
    344 => "HKG",
    348 => "HUN",
    352 => "ISL",
    356 => "IND",
    360 => "IDN",
    364 => "IRN",
    368 => "IRQ",
    372 => "IRL",
    833 => "IMN",
    376 => "ISR",
    380 => "ITA",
    388 => "JAM",
    392 => "JPN",
    832 => "JEY",
    400 => "JOR",
    398 => "KAZ",
    404 => "KEN",
    296 => "KIR",
    408 => "PRK",
    410 => "KOR",
    414 => "KWT",
    417 => "KGZ",
    418 => "LAO",
    428 => "LVA",
    422 => "LBN",
    426 => "LSO",
    430 => "LBR",
    434 => "LBY",
    438 => "LIE",
    440 => "LTU",
    442 => "LUX",
    446 => "MAC",
    807 => "MKD",
    450 => "MDG",
    454 => "MWI",
    458 => "MYS",
    462 => "MDV",
    466 => "MLI",
    470 => "MLT",
    584 => "MHL",
    474 => "MTQ",
    478 => "MRT",
    480 => "MUS",
    175 => "MYT",
    484 => "MEX",
    583 => "FSM",
    498 => "MDA",
    492 => "MCO",
    496 => "MNG",
    499 => "MNE",
    500 => "MSR",
    504 => "MAR",
    508 => "MOZ",
    104 => "MMR",
    516 => "NAM",
    520 => "NRU",
    524 => "NPL",
    528 => "NLD",
    540 => "NCL",
    554 => "NZL",
    558 => "NIC",
    562 => "NER",
    566 => "NGA",
    570 => "NIU",
    574 => "NFK",
    580 => "MNP",
    578 => "NOR",
    512 => "OMN",
    586 => "PAK",
    585 => "PLW",
    275 => "PSE",
    591 => "PAN",
    598 => "PNG",
    600 => "PRY",
    604 => "PER",
    608 => "PHL",
    612 => "PCN",
    616 => "POL",
    620 => "PRT",
    630 => "PRI",
    634 => "QAT",
    638 => "REU",
    642 => "ROU",
    643 => "RUS",
    646 => "RWA",
    652 => "BLM",
    654 => "SHN",
    659 => "KNA",
    662 => "LCA",
    663 => "MAF",
    666 => "SPM",
    670 => "VCT",
    882 => "WSM",
    674 => "SMR",
    678 => "STP",
    682 => "SAU",
    686 => "SEN",
    688 => "SRB",
    690 => "SYC",
    694 => "SLE",
    702 => "SGP",
    534 => "SXM",
    703 => "SVK",
    705 => "SVN",
    090 => "SLB",
    706 => "SOM",
    710 => "ZAF",
    239 => "SGS",
    728 => "SSD",
    724 => "ESP",
    144 => "LKA",
    729 => "SDN",
    740 => "SUR",
    744 => "SJM",
    752 => "SWE",
    756 => "CHE",
    760 => "SYR",
    158 => "TWN",
    762 => "TJK",
    834 => "TZA",
    764 => "THA",
    626 => "TLS",
    768 => "TGO",
    772 => "TKL",
    776 => "TON",
    780 => "TTO",
    788 => "TUN",
    792 => "TUR",
    795 => "TKM",
    796 => "TCA",
    798 => "TUV",
    800 => "UGA",
    804 => "UKR",
    784 => "ARE",
    826 => "GBR",
    840 => "USA",
    581 => "UMI",
    858 => "URY",
    860 => "UZB",
    548 => "VUT",
    862 => "VEN",
    704 => "VNM",
    092 => "VGB",
    850 => "VIR",
    876 => "WLF",
    732 => "ESH",
    887 => "YEM",
    894 => "ZMB",
    716 => "ZWE",
    typemax(Int32) => "MISSNG"
)


text:
using Mods
using Primes

struct PrimeGroup <: CyclicGroup
    G::Mod
    q::Integer
    t::Int # security of the group
    function PrimeGroup(G,q,t)
        @assert isprime(G.mod)
        # I could also test by calculating Euler totient
        new(G,q,t)
    end
end

PrimeGroup(g,p,q,t) = PrimeGroup(Mod(g,p),q,t)
PrimeGroup(val::Integer,G::PrimeGroup) = PrimeGroup(val,G.G.mod,G.q,G.t)

value(G::PrimeGroup) = G.G.val
order(G::PrimeGroup) = G.q
security(G::PrimeGroup) = G.t

import Base.*
function *(X::PrimeGroup,Y::PrimeGroup)
    # Parametrizing group with p and q is not an option as that would recompile the code.
    qx,qy = order(X), order(Y)
    px,py = X.G.mod, Y.G.mod

    if px!=py || qx!=qy
        error("Groups are not equal")
    else
        G = X.G * Y.G
        return PrimeGroup(G,X.q,X.t)
    end
end

import Base.==
function ==(X::PrimeGroup,Y::PrimeGroup)
    qx,qy = order(X), order(Y)
    px,py = X.G.mod, Y.G.mod

    if px!=py || qx!=qy
        error("Groups are not equal")
    else
        return X.G==Y.G
    end
end

### Some prime group generation algorithms. References:
# + https://crypto.stackexchange.com/questions/820/how-does-one-calculate-a-primitive-root-for-diffie-hellman
# + https://math.stackexchange.com/questions/124408/finding-a-primitive-root-of-a-prime-number

### I could manage to get 512 security
function SophieGermainGroup(rng::AbstractRNG,g::Integer,t::Int)
    @assert g!=1
    while true
        q = rngprime(rng,2*t)
        p = 2*q + 1
        if g!=p-1 && isprime(p)
            return PrimeGroup(g,p,2*q,t)
        end
    end
end

function generateqp(rng::AbstractRNG,t)
    while true
        q = rngprime(rng,t)
        r = rngint(rng,t)

        p = q*r + 1
        if isprime(p)
            return q,p
        end
    end
end

function DSAStandartGroup(rng::AbstractRNG,tp::Int,tg::Int)
    q,p = generateqp(rng,2*tp)
    while true
        u = Mod(rngint(rng,tg),p)
        g = u^div((p-1),q)
        if g.val!=1 || g.val!=0
            return PrimeGroup(g,q,tp)
        end
    end
end


text:
struct xReLU{N,Aθ,Aγ,AΔ,Aξ} <: AbstractLayer{N}
    θ::Aθ
    γ::Aγ
    Δ::AΔ
    ξ::Aξ
    function xReLU(θ::AbstractArray, γ::AbstractArray, Δ::AbstractArray, ξ::AbstractArray)
        @assert size(θ) == size(γ) == size(Δ) == size(ξ)
        return new{ndims(θ), typeof(θ), typeof(γ), typeof(Δ), typeof(ξ)}(θ, γ, Δ, ξ)
    end
end

function xReLU(::Type{T}, n::Int...) where {T}
    θ = zeros(T, n)
    γ = ones(T, n)
    Δ = zeros(T, n)
    ξ = zeros(T, n)
    return xReLU(θ, γ, Δ, ξ)
end

xReLU(n::Int...) = xReLU(Float64, n...)

Base.repeat(l::xReLU, n::Int...) = xReLU(
    repeat(l.θ, n...), repeat(l.γ, n...), repeat(l.Δ, n...), repeat(l.ξ, n...)
)

energies(layer::xReLU, x::AbstractArray) = energies(dReLU(layer), x)
free_energies(layer::xReLU, inputs::Union{Real,AbstractArray} = 0) = free_energies(dReLU(layer), inputs)
sample_from_inputs(layer::xReLU, inputs::Union{Real,AbstractArray} = 0) = sample_from_inputs(dReLU(layer), inputs)
mode_from_inputs(layer::xReLU, inputs::Union{Real,AbstractArray} = 0) = mode_from_inputs(dReLU(layer), inputs)
mean_from_inputs(layer::xReLU, inputs::Union{Real,AbstractArray} = 0) = mean_from_inputs(dReLU(layer), inputs)
var_from_inputs(layer::xReLU, inputs::Union{Real,AbstractArray} = 0) = var_from_inputs(dReLU(layer), inputs)
meanvar_from_inputs(layer::xReLU, inputs::Union{Real,AbstractArray} = 0) = meanvar_from_inputs(dReLU(layer), inputs)
std_from_inputs(layer::xReLU, inputs::Union{Real,AbstractArray} = 0) = sqrt.(var_from_inputs(layer, inputs))
mean_abs_from_inputs(layer::xReLU, inputs::Union{Real,AbstractArray} = 0) = mean_abs_from_inputs(dReLU(layer), inputs)

function ∂free_energies(layer::xReLU, inputs::Union{Real,AbstractArray} = 0)
    drelu = dReLU(layer)

    lp = ReLU( drelu.θp, drelu.γp)
    ln = ReLU(-drelu.θn, drelu.γn)
    Fp = free_energies(lp,  inputs)
    Fn = free_energies(ln, -inputs)
    F = -logaddexp.(-Fp, -Fn)
    pp = exp.(F - Fp)
    pn = exp.(F - Fn)
    μp, νp = meanvar_from_inputs(lp,  inputs)
    μn, νn = meanvar_from_inputs(ln, -inputs)
    μ2p = @. νp + μp^2
    μ2n = @. νn + μn^2

    η = @. layer.ξ / (1 + abs(layer.ξ))

    ∂θ = @. -(pp * μp - pn * μn)
    ∂γ = @. (pp * μ2p / (1 + η) + pn * μ2n / (1 - η)) / 2
    ∂γ .*= sign.(layer.γ)
    ∂Δ = @. -(pp * μp / (1 + η) + pn * μn / (1 - η))
    abs_γ = abs.(layer.γ)
    ∂ξ = @. (
        pp * (-abs_γ/2 * μ2p + layer.Δ * μp) / (1 + layer.ξ + abs(layer.ξ))^2 +
        pn * ( abs_γ/2 * μ2n - layer.Δ * μn) / (1 - layer.ξ + abs(layer.ξ))^2
    )
    return (θ = ∂θ, γ = ∂γ, Δ = ∂Δ, ξ = ∂ξ)
end

const xReLUStats = pReLUStats

function ∂energy(layer::xReLU, stats::xReLUStats)
    @assert size(layer) == size(stats)
    η = @. layer.ξ / (1 + abs(layer.ξ))
    ∂γ = @. (stats.xp2 / (1 + η) + stats.xn2 / (1 - η)) / 2
    ∂γ .*= sign.(layer.γ)
    ∂Δ = @. -stats.xp1 / (1 + η) + stats.xn1 / (1 - η)
    abs_γ = abs.(layer.γ)
    ∂ξ = @. (
        (-abs_γ/2 * stats.xp2 + layer.Δ * stats.xp1) / (1 + layer.ξ + abs(layer.ξ))^2 +
        ( abs_γ/2 * stats.xn2 + layer.Δ * stats.xn1) / (1 - layer.ξ + abs(layer.ξ))^2
    )
    return (θ = -stats.x, γ = ∂γ, Δ = ∂Δ, ξ = ∂ξ)
end

suffstats(layer::xReLU, data::AbstractArray; wts = nothing) = xReLUStats(layer, data; wts)


text:
using MAT
using MatrixMarket

function countz(i)
    local count = 0
    for e in i
        if e != 0
            count += 1
        end
    end
    count
end

function write_graph(A, filename)
    open(filename, "w") do out
        x,y,z=findnz(A)

        A = sparse(x,y,map(a->round(Int, abs(a)), z))

        count = 0

        rows = rowvals(A)
        vals = nonzeros(A)
        x,y = size(A)
        d = countz(diag(A))
        write(out, "$x $(round(Int, (nnz(A)-d)/2))\n")
        for i in 1:y
            for j in nzrange(A, i)
                if rows[j] != i
                    count += 1
                    write(out, " $(rows[j])")
                end
            end
            write(out, "\n")
        end
    end
end

file = ARGS[1]
ext = splitext(file)[2]
A = if ext == ".mtx"
        mmread(file)
    else
        P = matread(file)["Problem"]
        A = P["A"]
    end
if all(x -> x == 0, sum(A, 1))
    A = spdiagm(diag(A)) - A
end
write_graph(A, ARGS[2])

max_stars_repo_name: ligmg/ligmg

max_issues_repo_name: ligmg/ligmg

max_forks_repo_name: ligmg/ligmg


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using NNFEM
Random.seed!(233)

function rivlin_saunders(strain,c1,c2)
    rivlin_saunders_ = load_op_and_grad("./build/libRivlinSaunders","rivlin_saunders")
    strain,c1,c2 = convert_to_tensor([strain,c1,c2], [Float64,Float64,Float64])
    rivlin_saunders_(strain,c1,c2)
end

c1 = 2.0
c2 = 3.0
N = 100
ε = rand(N,3)
σ = zeros(N,3)
for i = 1:N
    σ[i,:] = NNFEM.PlaneStressIncompressibleRivlinSaundersStress(ε[i,1],ε[i,2],ε[i,3],c1,c2)
end

# TODO: specify your input parameters
u = rivlin_saunders(ε,c1,c2)
sess = Session(); init(sess)
@show run(sess, u)-σ

# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(m)
    # return sum(rivlin_saunders(ε,c1,c2)^2)
    # return sum(rivlin_saunders(m,c1,c2)^2)
    return sum(rivlin_saunders(ε,c1,m)^2)
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(rand())
v_ = rand()
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
include(joinpath("symplectic_euler", "integrator.jl"))
include(joinpath("symplectic_euler", "body.jl"))
# include(joinpath("symplectic_euler", "joint.jl"))
include(joinpath("symplectic_euler", "bound.jl"))
include(joinpath("symplectic_euler", "test.jl"))


text:
#Pkg.add("ParallelDataTransfer")

addprocs(3) #number of additional processors

@everywhere using ParallelDataTransfer
@everywhere using ExoplanetsSysSim

import DataArrays.skipmissing

@everywhere include("clusters.jl")
@everywhere include("planetary_catalog.jl")
@everywhere include("optimization.jl")

@everywhere sim_param = setup_sim_param_model()

##### To start saving the model iterations in the optimization into a file:

@everywhere add_param_fixed(sim_param,"num_targets_sim_pass_one",78005)

model_name = "Clustered_P_R_broken_R_optimization"
optimization_number = "_random"*ARGS[1] #if want to run on the cluster with random initial active parameters: "_random"*ARGS[1]
use_KS_or_AD = "KS" #'KS' or 'AD'
AD_mod = true
Kep_or_Sim = "Sim" #'Kep' or 'Sim'
max_evals = 3000
num_evals_weights = 20
Pop_per_param = 4

file_name = model_name*optimization_number*"_targs"*string(get_int(sim_param,"num_targets_sim_pass_one"))*"_evals"*string(max_evals)

sendto(workers(), max_evals=max_evals, file_name=file_name)

@everywhere f = open(file_name*"_worker"*string(myid())*".txt", "w")
println(f, "# All initial parameters:")
write_model_params(f, sim_param)

##### To run the same model multiple times to see how it compares to a simulated catalog with the same parameters:

srand(1234) #to have the same reference catalog and simulated catalogs for calculating the weights

#To generate a simulated catalog to fit to:
cat_phys = generate_kepler_physical_catalog(sim_param)
cat_phys_cut = ExoplanetsSysSim.generate_obs_targets(cat_phys,sim_param)
cat_obs = observe_kepler_targets_single_obs(cat_phys_cut,sim_param)
summary_stat_ref = calc_summary_stats_model(cat_obs,sim_param)

@passobj 1 workers() summary_stat_ref #to send the 'summary_stat_ref' object to all workers

#To simulate more observed planets for the subsequent model generations:
@everywhere add_param_fixed(sim_param,"max_incl_sys",80.0) #degrees; 0 (deg) for isotropic system inclinations; set closer to 90 (deg) for more transiting systems

active_param_true, weights, target_fitness, target_fitness_std = compute_weights_target_fitness_std_perfect_model(num_evals_weights, use_KS_or_AD ; AD_mod=AD_mod, weight=true, save_dist=true)

##### To draw the initial values of the active parameters randomly within a search range:

active_param_keys = ["break_radius", "log_rate_clusters", "log_rate_planets_per_cluster", "power_law_P", "power_law_r1", "power_law_r2"]
    #["break_radius", "log_rate_clusters", "log_rate_planets_per_cluster", "mr_power_index", "num_mutual_hill_radii", "power_law_P", "power_law_r1", "power_law_r2", "sigma_hk", "sigma_incl", "sigma_incl_near_mmr", "sigma_log_radius_in_cluster", "sigma_logperiod_per_pl_in_cluster"]
active_params_box = [(0.5*ExoplanetsSysSim.earth_radius, 10.*ExoplanetsSysSim.earth_radius), (log(1.), log(5.)), (log(1.), log(5.)), (-0.5, 1.5), (-6., 0.), (-6., 0.)] #search ranges for all of the active parameters
    #[(0.5*ExoplanetsSysSim.earth_radius, 10.*ExoplanetsSysSim.earth_radius), (log(1.), log(5.)), (log(1.), log(5.)), (1., 4.), (3., 20.), (-0.5, 1.5), (-6., 0.), (-6., 0.), (0., 0.1), (0., 5.), (0., 5.), (0.1, 1.0), (0., 0.3)] #search ranges for all of the active parameters

#To randomly draw (uniformly) a value for each active model parameter within its search range:

srand() #to have a random set of initial parameters and optimization run

for (i,param_key) in enumerate(active_param_keys)
    active_param_draw = active_params_box[i][1] + (active_params_box[i][2] - active_params_box[i][1])*rand(1)
    add_param_active(sim_param,param_key,active_param_draw[1])
end
active_param_start = make_vector_of_sim_param(sim_param)

PopSize = length(active_param_true)*Pop_per_param

println("# Active parameters: ", make_vector_of_active_param_keys(sim_param))
println(f, "# Active parameters: ", make_vector_of_active_param_keys(sim_param))
println(f, "# Starting active parameter values: ", active_param_start)
println(f, "# Optimization active parameters search bounds: ", active_params_box)
println(f, "# Method: dxnes")
println(f, "# PopulationSize: ", PopSize)
println(f, "# Format: Active_params: [active parameter values]")
println(f, "# Format: Dist: [distances][total distance]")
println(f, "# Format: Dist_weighted: [weighted distances][total weighted distance]")
println(f, "# Distances used: ", use_KS_or_AD)
println(f, "# AD_mod: ", AD_mod)

target_function(active_param_start, use_KS_or_AD, Kep_or_Sim ; AD_mod=AD_mod, weights=weights, all_dist=false, save_dist=true) #to simulate the model once with the drawn parameters before starting the optimization

##### To use automated optimization routines to optimize the active model parameters:

# Pkg.add("BlackBoxOptim")       # only need to do these once
# Pkg.checkout("BlackBoxOptim")  # needed to get the lastest version
using BlackBoxOptim              # see https://github.com/robertfeldt/BlackBoxOptim.jl for documentation

tic()
opt_result = bboptimize(active_params -> target_function(active_params, use_KS_or_AD, Kep_or_Sim ; AD_mod=AD_mod, weights=weights, all_dist=false, save_dist=true); SearchRange = active_params_box, NumDimensions = length(active_param_true), Method = :dxnes, PopulationSize = PopSize, MaxFuncEvals = max_evals, TargetFitness = target_fitness, FitnessTolerance = target_fitness_std, TraceMode = :verbose, Workers = workers())
t_elapsed = toc()

println(f, "# best_candidate: ", best_candidate(opt_result))
println(f, "# best_fitness: ", best_fitness(opt_result))
println(f, "# elapsed time: ", t_elapsed, " seconds")
@everywhere close(f)


text:
function get_nous1()
m = Model()

# ----- Variables ----- #
@variable(m, objvar)
x_Idx = Any[2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49]
@variable(m, x[x_Idx])
b_Idx = Any[50, 51]
@variable(m, b[b_Idx], Bin)
JuMP.set_lower_bound(x[16], 0.0)
JuMP.set_lower_bound(x[14], 0.0)
JuMP.set_lower_bound(x[38], 0.0)
JuMP.set_lower_bound(x[42], 0.0)
JuMP.set_lower_bound(x[22], 0.0)
JuMP.set_lower_bound(x[2], 0.0)
JuMP.set_lower_bound(x[9], 0.0)
JuMP.set_lower_bound(x[8], 0.0)
JuMP.set_lower_bound(x[43], 0.0)
JuMP.set_lower_bound(x[36], 0.0)
JuMP.set_lower_bound(x[4], 0.0)
JuMP.set_lower_bound(x[32], 0.0)
JuMP.set_lower_bound(x[27], 0.0)
JuMP.set_lower_bound(x[3], 0.0)
JuMP.set_lower_bound(x[25], 0.0)
JuMP.set_lower_bound(x[30], 0.0)
JuMP.set_lower_bound(x[11], 0.0)
JuMP.set_lower_bound(x[29], 0.0)
JuMP.set_lower_bound(x[5], 0.0)
JuMP.set_lower_bound(x[37], 0.0)
JuMP.set_lower_bound(x[24], 0.0)
JuMP.set_lower_bound(x[41], 0.0)
JuMP.set_lower_bound(x[18], 0.0)
JuMP.set_lower_bound(x[7], 0.0)
JuMP.set_lower_bound(x[13], 0.0)
JuMP.set_lower_bound(x[21], 0.0)
JuMP.set_lower_bound(x[10], 0.0)
JuMP.set_lower_bound(x[26], 0.0)
JuMP.set_lower_bound(x[45], 0.0)
JuMP.set_lower_bound(x[12], 0.0)
JuMP.set_lower_bound(x[40], 0.0)
JuMP.set_lower_bound(x[44], 0.0)
JuMP.set_lower_bound(x[31], 0.0)
JuMP.set_lower_bound(x[33], 0.0)
JuMP.set_lower_bound(x[28], 0.0)
JuMP.set_lower_bound(x[35], 0.0)
JuMP.set_lower_bound(x[6], 0.0)
JuMP.set_lower_bound(x[17], 0.0)
JuMP.set_lower_bound(x[23], 0.0)
JuMP.set_lower_bound(x[34], 0.0)
JuMP.set_lower_bound(x[19], 0.0)
JuMP.set_lower_bound(x[20], 0.0)
JuMP.set_lower_bound(x[39], 0.0)
JuMP.set_lower_bound(x[15], 0.0)
JuMP.set_upper_bound(x[2], 300.0)
JuMP.set_upper_bound(x[3], 300.0)
JuMP.set_upper_bound(x[4], 300.0)
JuMP.set_upper_bound(x[5], 300.0)
JuMP.set_upper_bound(x[6], 300.0)
JuMP.set_upper_bound(x[7], 300.0)
JuMP.set_upper_bound(x[8], 300.0)
JuMP.set_upper_bound(x[9], 300.0)
JuMP.set_upper_bound(x[10], 300.0)
JuMP.set_upper_bound(x[11], 300.0)
JuMP.set_upper_bound(x[12], 300.0)
JuMP.set_upper_bound(x[13], 300.0)
JuMP.set_upper_bound(x[14], 300.0)
JuMP.set_upper_bound(x[15], 300.0)
JuMP.set_upper_bound(x[16], 300.0)
JuMP.set_upper_bound(x[17], 300.0)
JuMP.set_upper_bound(x[18], 300.0)
JuMP.set_upper_bound(x[19], 300.0)
JuMP.set_upper_bound(x[22], 100.0)
JuMP.set_upper_bound(x[23], 100.0)
JuMP.set_upper_bound(x[24], 100.0)
JuMP.set_upper_bound(x[25], 100.0)
JuMP.set_upper_bound(x[26], 100.0)
JuMP.set_upper_bound(x[27], 100.0)
JuMP.set_upper_bound(x[28], 1.0)
JuMP.set_upper_bound(x[29], 1.0)
JuMP.set_upper_bound(x[30], 1.0)
JuMP.set_upper_bound(x[31], 1.0)
JuMP.set_upper_bound(x[32], 1.0)
JuMP.set_upper_bound(x[33], 1.0)
JuMP.set_upper_bound(x[34], 1.0)
JuMP.set_upper_bound(x[35], 1.0)
JuMP.set_upper_bound(x[36], 1.0)
JuMP.set_upper_bound(x[37], 1.0)
JuMP.set_upper_bound(x[38], 1.0)
JuMP.set_upper_bound(x[39], 1.0)
JuMP.set_upper_bound(x[40], 1.0)
JuMP.set_upper_bound(x[41], 1.0)
JuMP.set_upper_bound(x[42], 1.0)
JuMP.set_upper_bound(x[43], 1.0)
JuMP.set_upper_bound(x[44], 1.0)
JuMP.set_upper_bound(x[45], 1.0)
JuMP.set_lower_bound(x[46], 0.85)
JuMP.set_upper_bound(x[46], 1.0)
JuMP.set_lower_bound(x[47], 0.85)
JuMP.set_upper_bound(x[47], 1.0)
JuMP.set_lower_bound(x[48], 0.85)
JuMP.set_upper_bound(x[48], 1.0)
JuMP.set_lower_bound(x[49], 0.85)
JuMP.set_upper_bound(x[49], 1.0)

# ----- Constraints ----- #
@NLconstraint(m, e1, -((0.0042656*x[29]-0.0005719*x[28]+0.0093514*x[46]+0.0077308*x[47]-0.0139904)*x[4]+(0.0016371*x[31]+0.0288996*x[32]+0.0338147*x[48]+0.0373349*x[49]-0.0661588)*x[5])+objvar-0.23947*b[50]-0.75835*b[51] == 0.0)
@constraint(m, e2, x[2]+x[3]+x[20]+x[21] == 300.0)
@constraint(m, e3, x[6]-x[12]-x[13] == 0.0)
@constraint(m, e4, x[7]-x[11]-x[14]-x[15] == 0.0)
@constraint(m, e5, x[8]-x[10]-x[16]-x[17] == 0.0)
@constraint(m, e6, x[9]-x[18]-x[19] == 0.0)
@NLconstraint(m, e7, -x[10]*x[40]-0.333333333333333*x[2]+x[22] == 0.0)
@NLconstraint(m, e8, -x[10]*x[41]-0.333333333333333*x[2]+x[23] == 0.0)
@NLconstraint(m, e9, -x[10]*x[42]-0.333333333333333*x[2]+x[24] == 0.0)
@NLconstraint(m, e10, -x[11]*x[37]-0.333333333333333*x[3]+x[25] == 0.0)
@NLconstraint(m, e11, -x[11]*x[38]-0.333333333333333*x[3]+x[26] == 0.0)
@NLconstraint(m, e12, -x[11]*x[39]-0.333333333333333*x[3]+x[27] == 0.0)
@NLconstraint(m, e13, -(x[6]*x[34]+x[7]*x[37])+x[22] == 0.0)
@NLconstraint(m, e14, -(x[6]*x[35]+x[7]*x[38])+x[23] == 0.0)
@NLconstraint(m, e15, -(x[6]*x[36]+x[7]*x[39])+x[24] == 0.0)
@NLconstraint(m, e16, -(x[8]*x[40]+x[9]*x[43])+x[25] == 0.0)
@NLconstraint(m, e17, -(x[8]*x[41]+x[9]*x[44])+x[26] == 0.0)
@NLconstraint(m, e18, -(x[8]*x[42]+x[9]*x[45])+x[27] == 0.0)
@NLconstraint(m, e19, x[22]*x[46]-x[6]*x[34] == 0.0)
@NLconstraint(m, e20, x[23]*x[47]-x[7]*x[38] == 0.0)
@NLconstraint(m, e21, x[26]*x[48]-x[8]*x[41] == 0.0)
@NLconstraint(m, e22, x[27]*x[49]-x[9]*x[45] == 0.0)
@NLconstraint(m, e23, x[12]*x[34]+x[14]*x[37]+x[16]*x[40]+x[18]*x[43]+0.333333333333333*x[20] == 30.0)
@NLconstraint(m, e24, x[12]*x[35]+x[14]*x[38]+x[16]*x[41]+x[18]*x[44]+0.333333333333333*x[20] == 50.0)
@NLconstraint(m, e25, x[12]*x[36]+x[14]*x[39]+x[16]*x[42]+x[18]*x[45]+0.333333333333333*x[20] == 30.0)
@NLconstraint(m, e26, x[13]*x[34]+x[15]*x[37]+x[17]*x[40]+x[19]*x[43]+0.333333333333333*x[21] == 70.0)
@NLconstraint(m, e27, x[13]*x[35]+x[15]*x[38]+x[17]*x[41]+x[19]*x[44]+0.333333333333333*x[21] == 50.0)
@NLconstraint(m, e28, x[13]*x[36]+x[15]*x[39]+x[17]*x[42]+x[19]*x[45]+0.333333333333333*x[21] == 70.0)
@NLconstraint(m, e29, x[4]*x[28]-x[22] == 0.0)
@NLconstraint(m, e30, x[4]*x[29]-x[23] == 0.0)
@NLconstraint(m, e31, x[4]*x[30]-x[24] == 0.0)
@NLconstraint(m, e32, x[5]*x[31]-x[25] == 0.0)
@NLconstraint(m, e33, x[5]*x[32]-x[26] == 0.0)
@NLconstraint(m, e34, x[5]*x[33]-x[27] == 0.0)
@constraint(m, e35, x[34]+x[35]+x[36] == 1.0)
@constraint(m, e36, x[37]+x[38]+x[39] == 1.0)
@constraint(m, e37, x[40]+x[41]+x[42] == 1.0)
@constraint(m, e38, x[43]+x[44]+x[45] == 1.0)
@constraint(m, e39, x[28]+x[29]+x[30] == 1.0)
@constraint(m, e40, x[31]+x[32]+x[33] == 1.0)
@constraint(m, e41, x[36] == 0.0)
@constraint(m, e42, x[43] == 0.0)
@constraint(m, e43, x[4]-300*b[50] <= 0.0)
@constraint(m, e44, x[5]-300*b[51] <= 0.0)

# ----- Objective ----- #
@objective(m, Min, objvar)
return m
end


text:
using ProgressLogging
using NBodySimulator, OrdinaryDiffEq, StaticArrays
using Plots, DataFrames, StatsPlots

function setup(t)
    T = 120.0 # K
    kb = 1.38e-23 # J/K
    ϵ = T * kb # J
    σ = 3.4e-10 # m
    ρ = 1374 # kg/m^3
    m = 39.95 * 1.6747 * 1e-27 # kg
    N = 350
    L = (m*N/ρ)^(1/3)
    R = 3.5σ
    v_dev = sqrt(kb * T / m) # m/s

    _L = L / σ
    _σ = 1.0
    _ϵ = 1.0
    _m = 1.0
    _v = v_dev / sqrt(ϵ / m)
    _R = R / σ

    bodies = generate_bodies_in_cell_nodes(N, _m, _v, _L)
    lj_parameters = LennardJonesParameters(_ϵ, _σ, _R)
    pbc = CubicPeriodicBoundaryConditions(_L)
    lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => lj_parameters));
    simulation = NBodySimulation(lj_system, (0.0, t), pbc, _ϵ/T)

    return simulation
end

function benchmark(energyerr, rts, bytes, allocs, nt, nf, t, configs)
    simulation = setup(t)
    prob = SecondOrderODEProblem(simulation)
    for config in configs
        alg = config.alg
        sol, rt, b, gc, memalloc = @timed solve(prob, alg();
            save_everystep=false, progress=true, progress_name="$alg", config...)
        result = NBodySimulator.SimulationResult(sol, simulation)
        ΔE = total_energy(result, t) - total_energy(result, 0)
        energyerr[alg] = ΔE
        rts[alg] = rt
        bytes[alg] = b
        allocs[alg] = memalloc
        nt[alg] = sol.destats.naccept
        nf[alg] = sol.destats.nf + sol.destats.nf2
    end
end

function run_benchmark!(results, t, integrators, tol...; c=ones(length(integrators)))
    @progress "Benchmark at t=$t" for τ in zip(tol...)
        runtime = Dict()
        ΔE = Dict()
        nt = Dict()
        nf = Dict()
        b = Dict()
        allocs = Dict()
        cfg = config(integrators, c, τ...)

        GC.gc()
        benchmark(ΔE, runtime, b, allocs, nt, nf, t, cfg)
        get_tol(idx) = haskey(cfg[idx], :dt) ? cfg[idx].dt : (cfg[idx].abstol, cfg[idx].rtol)

        for (idx,i) in enumerate(integrators)
            push!(results, [string(i), runtime[i], get_tol(idx)..., abs(ΔE[i]), nt[i], nf[i], c[idx]])
        end
    end
    return results
end

symplectic_integrators = [
    VelocityVerlet,
    VerletLeapfrog,
    PseudoVerletLeapfrog,
    McAte2,
    CalvoSanz4,
    McAte5,
    Yoshida6,
    KahanLi8,
    SofSpa10
];

config(integrators, c, τ) = [ (alg=a, dt=τ*cₐ) for (a,cₐ) in zip(integrators, c)]

t = 35.0
τs = 1e-3

# warmup
c_symplectic = ones(length(symplectic_integrators))
benchmark(Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), 10.,
    config(symplectic_integrators, c_symplectic, τs))

# results = DataFrame(:integrator=>String[], :runtime=>Float64[], :τ=>Float64[],
#     :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
# run_benchmark!(results, t, symplectic_integrators, τs)

# c_symplectic .= results[!, :runtime] ./ results[!, :timesteps]
# c_Verlet = c_symplectic[1]
# c_symplectic /= c_Verlet

c_symplectic = [
    1.00,   # VelocityVerlet
    1.05,   # VerletLeapfrog
    0.98,   # PseudoVerletLeapfrog
    1.02,   # McAte2
    2.38,   # CalvoSanz4
    2.92,   # McAte5
    3.74,   # Yoshida6
    8.44,   # KahanLi8
    15.76   # SofSpa10
]

t = 40.0
τs = 10 .^range(-4, -3, length=10)

results = DataFrame(:integrator=>String[], :runtime=>Float64[], :τ=>Float64[],
    :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
run_benchmark!(results, t, symplectic_integrators, τs, c=c_symplectic)

@df results plot(:EnergyError, :runtime, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")

@df results plot(:timesteps, :runtime, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Number of timesteps", ylabel="Runtime (s)")

function benchmark(energyerr, rts, ts, t, configs)
    simulation = setup(t)
    prob = SecondOrderODEProblem(simulation)
    for config in configs
        alg = config.alg
        sol, rt = @timed solve(prob, alg(); progress=true, progress_name="$alg", config...)
        result = NBodySimulator.SimulationResult(sol, simulation)
        ΔE(t) = total_energy(result, t) - total_energy(result, 0)
        energyerr[alg] = [ΔE(t) for t in sol.t[2:10^2:end]]
        rts[alg] = rt
        ts[alg] = sol.t[2:10^2:end]
    end
end

ΔE = Dict()
rt = Dict()
ts = Dict()
configs = config(symplectic_integrators, c_symplectic, 2.3e-4)
benchmark(ΔE, rt, ts, 40., configs)

plt = plot(xlabel="Rescaled Time", ylabel="Energy error", legend=:bottomleft);
for c in configs
    plot!(plt, ts[c.alg], abs.(ΔE[c.alg]), label="$(c.alg), $(rt[c.alg])s")
end
plt

adaptive_integrators=[
    # Non-stiff ODE methods
    Tsit5,
    Vern7,
    Vern9,
    # DPRKN
    DPRKN6,
    DPRKN8,
    DPRKN12,
];

config(integrators, c, at, rt) = [ (alg=a, abstol=at*2^cₐ, rtol=rt*2^cₐ) for (a,cₐ) in zip(integrators, c)]

t = 35.0
ats = 10 .^range(-14, -4, length=10)
rts = 10 .^range(-14, -4, length=10)

# warmup
c_adaptive = ones(length(adaptive_integrators))
benchmark(Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), 10.,
    config(adaptive_integrators, 1, ats[1], rts[1]))

# results = DataFrame(:integrator=>String[], :runtime=>Float64[], :abstol=>Float64[],
#    :reltol=>Float64[], :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
# run_benchmark!(results, t, adaptive_integrators, ats[1], rts[1])

# c_adaptive .= results[!, :runtime] ./ results[!, :timesteps]
# c_adaptive /= c_Verlet

c_adaptive = [
    3.55,   # Tsit5,
    7.84,   # Vern7,
    11.38,  # Vern9
    3.56,   # DPRKN6,
    5.10,   # DPRKN8,
    8.85    # DPRKN12,
]

t = 40.0

results = DataFrame(:integrator=>String[], :runtime=>Float64[], :abstol=>Float64[],
    :reltol=>Float64[], :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
run_benchmark!(results, t, adaptive_integrators, ats, rts, c=c_adaptive)

@df results plot(:EnergyError, :runtime, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")

@df results plot(:EnergyError, :f_evals, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Number of f evals")

t = 40.0

symplectic_integrators = [
    VelocityVerlet,
    VerletLeapfrog,
    PseudoVerletLeapfrog,
    McAte2,
    CalvoSanz4
]

c_symplectic = [
    1.00,   # VelocityVerlet
    1.05,   # VerletLeapfrog
    0.98,   # PseudoVerletLeapfrog
    1.02,   # McAte2
    2.38,   # CalvoSanz4
]

results1 = DataFrame(:integrator=>String[], :runtime=>Float64[], :τ=>Float64[],
    :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
run_benchmark!(results1, t, symplectic_integrators, τs, c=c_symplectic)

adaptive_integrators=[
    DPRKN6,
    DPRKN8,
    DPRKN12,
]

c_adaptive = [
    3.56,   # DPRKN6,
    5.10,   # DPRKN8,
    8.85    # DPRKN12,
]

results2 = DataFrame(:integrator=>String[], :runtime=>Float64[], :abstol=>Float64[],
    :reltol=>Float64[], :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
run_benchmark!(results2, t, adaptive_integrators, ats, rts, c=c_adaptive)

append!(results1, results2, cols=:union)
results1

@df results1 plot(:EnergyError, :runtime, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")


text:
module DataHelpers

using JSON
using GraffSDK
using DocStringExtensions
using Base64

"""
$(SIGNATURES)
Encode data and return the data request.
"""
function encodeJsonData(id::String, description::String, data::Any)::BigDataElementRequest
    return BigDataElementRequest(id, "Mongo", description, JSON.json(data), "application/json");
end

function encodeBinaryData(id::String, description::String, data::Vector{UInt8}; mimeType="application/octet-stream")
    return BigDataElementRequest(id, "Mongo", description, base64encode(data), mimeType);
end

"""
$(SIGNATURES)
Read an image from a file and encode it as a Graff data element.
"""
function readFileIntoDataRequest(file::String, id::String, description::String, mimeType)::BigDataElementRequest
    try
        fid = open(file,"r");
        imgBytes = read(fid);
        close(fid);
        return BigDataElementRequest(id, "Mongo", description, base64encode(imgBytes), mimeType);
    catch ex
        showerror(stderr, ex)
        error("Unable to read the image from $file - $ex")
    end
end

"""
$(SIGNATURES)
Check if a datastructure is safe for JSON serialization.
"""
function isSafeToJsonSerialize(data::Any)
    try
        JSON.json(data)
        return true
    catch ex
        showerror(stderr, ex)
        return false
    end
end

end


text:
module GridapODEs

include("ODETools/ODETools.jl")

include("TransientFETools/TransientFETools.jl")

include("DiffEqsWrappers/DiffEqsWrappers.jl")

end #module


text:
using BinDeps

@BinDeps.setup

libSDL = library_dependency("libSDL", aliases = ["libSDL", "SDL"])
SDLgfx = library_dependency("libSDL_gfx", aliases = ["libSDL_gfx"], depends = [libSDL], os = :Unix)
SDLimage = library_dependency("libSDL_image", aliases = ["libSDL_image"], depends = [libSDL], os = :Unix)
SDLmixer = library_dependency("libSDL_mixer", aliases = ["libSDL_mixer"], depends = [libSDL], os = :Unix)
SDLttf = library_dependency("libSDL_ttf", aliases = ["libSDL_ttf"], depends = [libSDL], os = :Unix)

@windows_only begin
	using WinRPM
	provides(WinRPM.RPM, "libSDL", libSDL, os = :Windows)
end

@osx_only begin
	if Pkg.installed("Homebrew") === nothing
		error("Hombrew package not installed, please run Pkg.add(\"Homebrew\")")
	end
	using Homebrew
	provides(Homebrew.HB, "sdl", libSDL, os = :Darwin)
	provides(Homebrew.HB, "sdl_gfx", SDLgfx, os = :Darwin)
	provides(Homebrew.HB, "sdl_image", SDLimage, os = :Darwin)
	provides(Homebrew.HB, "sdl_mixer", SDLmixer, os = :Darwin)
	provides(Homebrew.HB, "sdl_ttf", SDLttf, os = :Darwin)
end

@linux_only begin
    provides(AptGet,
	    	{"libsdl1.2-dev" => libSDL,
	    	 "libsdl-gfx1.2-dev" => SDLgfx,
	    	 "libsdl-image1.2-dev" => SDLimage,
	    	 "libsdl-mixer1.2-dev" => SDLmixer,
	    	 "libsdl-ttf2.0-dev" => SDLttf})

    provides(Yum,
    		{"SDL-devel" => libSDL,
    		 "SDL_gfx-devel" => SDLgfx,
    		 "SDL_image-devel" => SDLimage,
    		 "SDL_mixer-devel" => SDLmixer,
    		 "SDL_ttf-devel" => SDLttf})
end

@BinDeps.install [:libSDL => :libSDL]


text:
@diffusion_process Lorenz96 begin
    :dimensions
    process --> DiffusionDefinition._TEMP_DIMENSION_PROCESS
    wiener --> DiffusionDefinition._TEMP_DIMENSION_WIENER

    :parameters
    (θ, σ) --> Float64

    :additional
    constdiff --> true
    diagonaldiff --> true
end

#-------------------------------------------------------------------------------
#                Some helper functions for static indexing

function si1(::DiffusionDefinition.DiffusionProcess{T,DP}) where {T,DP}
    SVector{DP,Int64}(mod1.(2:DP+1, DP))
end

function si2(::DiffusionDefinition.DiffusionProcess{T,DP}) where {T,DP}
    SVector{DP,Int64}(mod1.(-1:DP-2, DP))
end

function si3(::DiffusionDefinition.DiffusionProcess{T,DP}) where {T,DP}
    SVector{DP,Int64}(mod1.(0:DP-1, DP))
end

function si4(::DiffusionDefinition.DiffusionProcess{T,DP}) where {T,DP}
    SVector{DP,Int64}(1:DP)
end

function sid_mat(::DiffusionDefinition.DiffusionProcess{T,DP,DW}) where {T,DP,DW}
    SDiagonal{DW,T}(I)
end

#-------------------------------------------------------------------------------
#                Some helper functions for dynamic indexing

function i1(::DiffusionDefinition.DiffusionProcess{T,DP}) where {T,DP}
    mod1.(2:DP+1, DP)
end

function i2(::DiffusionDefinition.DiffusionProcess{T,DP}) where {T,DP}
    mod1.(-1:DP-2, DP)
end

function i3(::DiffusionDefinition.DiffusionProcess{T,DP}) where {T,DP}
    mod1.(0:DP-1, DP)
end

function i4(::DiffusionDefinition.DiffusionProcess{T,DP}) where {T,DP}
    1:DP
end

function id_mat(::DiffusionDefinition.DiffusionProcess{T,DP,DW}) where {T,DP,DW}
    Diagonal{T}(I, DW)
end

#-------------------------------------------------------------------------------
#                 Drift and diffusion coefficients

function b(t, x, P::Lorenz96)
    (
        (view(x, si1(P)) - view(x, si2(P))) .* view(x, si3(P))
        - view(x, si4(P)) .+ P.θ
    )
end

σ(t, x, P::Lorenz96) = P.σ * sid_mat(P)

function b!(buffer, t, x, P::Lorenz96)
    buffer.b .= (
        (view(x, i1(P)) .- view(x, i2(P))) .* view(x, i3(P))
        .- view(x, i4(P))
        .+ P.θ
    )
end

σ!(buffer, t, x, P::Lorenz96) = (buffer.σ.diag .= P.σ)


text:
module RefFEs

using Gridap
using Gridap.Helpers
using Gridap.DOFBases: LagrangianDOFBasis

export RefFE
export LagrangianRefFE
export shfbasis
export polytope
export nfacedofs
export dofbasis
export field_type

import Gridap: evaluate, evaluate!

# Abstract types and interfaces

"""
Abstract Reference Finite Element
"""
abstract type RefFE{D,T} end

dofbasis(this::RefFE{D,T} where {D,T})::DOFBasis{D,T} = @abstractmethod

polytope(this::RefFE{D,T} where {D,T})::Polytope{D} = @abstractmethod

shfbasis(this::RefFE{D,T} where {D,T})::Basis{D,T} = @abstractmethod

nfacedofs(this::RefFE{D,T} where {D,T})::Vector{Vector{Int}} = @abstractmethod

function field_type(this::RefFE{D,T} ) where {D,T}
  return T
end

# @fverdugo this abstract interface needs a tester

"""
Extract the lists of dofs, but only for the nfaces of dimension dim
"""
function nfacedofs(reffe::RefFE{D},dim::Integer) where D
  @assert 0 <= dim
  @assert dim <= D
  nface_to_dofs = nfacedofs(reffe)
  p = polytope(reffe)
  nface_to_dofs[nf_dim(p)[end][dim+1]]
end

"""
Reference Finite Element a la Ciarlet, i.e., it relies on a local function
(polynomial) space, an array of nodes (DOFs), and a polytope (cell topology).
The rank of the approximating field can be arbitrary. The current implementation
relies on the prebasis (e.g., monomial basis of polynomials) and a
change-of-basis (using the node array) to generate the canonical basis, i.e.,
the shape functions.
"""
struct LagrangianRefFE{D,T} <: RefFE{D,T}
  polytope::Polytope{D}
  dofbasis::LagrangianDOFBasis{D,T}
  shfbasis::Basis{D,T}
  nfacedofs::Vector{Vector{Int}}
  nfacenodes::Vector{Vector{Int}}
  # this type is unstable
end

"""
Constructor that given a vector of orders, generates high-order Lagrangian
RefFEs on n-cubes or n-tets.
"""
function LagrangianRefFE(::Type{T}, p::Polytope{D}, orders::Vector{Int}) where {D,T}
  @assert length(orders) == D
  @assert D > 0
  if !(all(extrusion(p).array.== HEX_AXIS) || all(orders.==orders[1]))
    error("One can consider anisotropic orders on n-cubes only")
  end
  if (all(orders.==1))
    nodes, nfacenodes = _linear_lagrangian_nodes_polytope(p)
  else
    nodes, nfacenodes = _high_order_lagrangian_nodes_polytope(p,orders)
  end
  dofsb = Gridap.RefFEs.LagrangianDOFBasis(T,nodes)
  # @santiagobadia : This part is missing for anisotropic orders !
  prebasis = Gridap.RefFEs._monomial_basis(p,T,orders[1])
  aux = zeros(Float64,numlocaldofs(dofsb),numlocaldofs(dofsb))
  @assert numlocaldofs(dofsb) == length(prebasis)
  evaluate!(dofsb,prebasis,aux)
  changeofbasis=inv(aux)
  basis = change_basis(prebasis, changeofbasis)
  nfacedofs = _genereate_nface_to_dofs(nfacenodes,dofsb.node_and_comp_to_dof)
  return LagrangianRefFE{D,T}(p, dofsb, basis, nfacedofs, nfacenodes)
end

"""
Version of the constructor for a scalar order
"""
function LagrangianRefFE(::Type{T}, p::Polytope{D}, order::Int) where {D,T}
  _order = order*ones(Int,D)
  return LagrangianRefFE(T,p,_order)
end

dofbasis(this::LagrangianRefFE{D,T} where {D,T}) = this.dofbasis

polytope(this::LagrangianRefFE{D,T} where {D,T}) = this.polytope

shfbasis(this::LagrangianRefFE{D,T} where {D,T}) = this.shfbasis

nfacedofs(this::LagrangianRefFE{D,T} where {D,T}) = this.nfacedofs

function closurenfacedofs(this::LagrangianRefFE{D,T} where {D,T})
  cv1 = CellValueFromArray(this.polytope.nf_nfs) # cell to index
  cv2 = CellValueFromArray(this.nfacedofs)
  return collect(CellVectorByComposition(cv1,cv2))
end

# Generate the linear nodes by computing the polytope vertices. Only the
# vertices have nodes.
function _linear_lagrangian_nodes_polytope(p::Polytope)
  function _linear_nfacedofs(p)
    nfacedofs = Vector{Int}[]
    for (inf,nf) in enumerate(nfaces(p)[nf_dim(p)[end][1]])
      push!(nfacedofs, [inf])
    end
    for id in 2:length(nf_dim(p)[end])
      for (inf,nf) in enumerate(nfaces(p)[nf_dim(p)[end][id]])
        push!(nfacedofs, Int[])
      end
    end
    return nfacedofs
  end
  nodes = Gridap.Polytopes.vertices_coordinates(p)
  nfacedofs = _linear_nfacedofs(p)
  return nodes, nfacedofs
end

# Determine the high order nodes of the polytope by traversing the n-faces,
# and generate their open nodes using a reference polytope of the n-face and
# a linear reference FE on top of it
function _high_order_lagrangian_nodes_polytope(p::Polytope, order)
  vs_p = Gridap.Polytopes.vertices_coordinates(p)
  ns_float_p = [i for i in vs_p]
  ref_ps = Gridap.Polytopes.nface_ref_polytopes(p)
  # rfe_p = Gridap.RefFEs._high_order_lagrangian_reffe(p,Float64,1)
  rfe_p = LagrangianRefFE(Float64,p,1)
  nfacedofs = copy(rfe_p.nfacedofs)
  nfs = nfaces(p)
  k = length(vs_p)
  for inf_dim = 1:length(nf_dim(p)[end])-1
    nfs_dim = nf_dim(p)[end][inf_dim+1]
    nfs_vs = Gridap.Polytopes._dimfrom_fs_dimto_fs(p,inf_dim,0)
    for (i_nf_dim,i_nf) in enumerate(nf_dim(p)[end][inf_dim+1])
      ref_p = ref_ps[i_nf]
      # _order = Tuple(order*ones(Int,length(ref_p.extrusion)))
      _order = _extract_nonzeros(nfs[i_nf].extrusion,order)
      ns = Gridap.Polytopes._interior_nodes_int_coords(ref_p, _order)
      ns_float = Gridap.Polytopes._interior_nodes_int_to_real_coords(ns,_order)
      rfe = LagrangianRefFE(Float64,ref_p,1)
      nf_vs = nfs_vs[i_nf_dim]
      vs = vs_p[nf_vs]
      if ( length(ns_float) > 0 )
        a = evaluate(rfe.shfbasis,ns_float)
        nfacedofs[i_nf] = [k+i for i in 1:length(ns)]
        k += length(ns)
        hons = Gridap.RefFEs._map_high_order_lagrangian_nodes(a, vs)
        for i in hons
          push!(ns_float_p, i)
        end
      else
        nfacedofs[i_nf] = Int[]
      end
    end
  end
  return ns_float_p, nfacedofs
end

# Given the evaluation of linear shape functions of the reference
# polytope of a given n-face evaluated on the open high order nodes
# in the reference space, map these high order nodes in the n-face
# by combining the shape functions with the n-face vertices
function _map_high_order_lagrangian_nodes(shfs_ns, vs)
  a = shfs_ns
  T = typeof(vs[1].array)
  ndofs, npoints = size(a)
  v = Point{length(vs[1]),Float64}[]
  for j in 1:npoints
    aux = zero(T)
    for i in 1:ndofs
      aux += outer(a[i,j],vs[i]).array
    end
    cs = convert(Vector{Float64}, [aux...])
    push!(v, aux)
  end
  return v
end

# @santiagobadia : The basis for prisms, pyramids, etc, not implemented
# It would require to enrich the polynomial machinery, too compicated for a
# filter, I think it would not be efficient
function _monomial_basis(p::Polytope{D}, T, order) where D
  if (_is_hex(p))
    orders = order*ones(Int,D)
    prebasis = MonomialBasis(T,orders)
  elseif (_is_tet(p))
    filter(e,order) = sum(e) <= order
    prebasis = MonomialBasis(D, T, filter, order)
  else
    @notimplemented
  end
end

function _is_hex(p)
  for i in 2:length(extrusion(p))
    if extrusion(p)[i] != 1 return false && break end
  end
  return true
end

function _is_tet(p)
  for i in 2:length(extrusion(p))
    if extrusion(p)[i] != 2 return false && break end
  end
  return true
end

function _genereate_nface_to_dofs(nface_to_nodes,node_and_comp_to_dof)
  nnode, ncomp = size(node_and_comp_to_dof)
  nface_to_dofs = [zeros(Int,length(nodes)*ncomp) for nodes in nface_to_nodes]
  nnface = length(nface_to_nodes)
  nface_to_i = ones(Int,nnface)
  for (nface,nodes) in enumerate(nface_to_nodes)
    for comp in 1:ncomp
      for node in nodes
        dof = node_and_comp_to_dof[node,comp]
        i = nface_to_i[nface]
        nface_to_dofs[nface][i] = dof
        nface_to_i[nface] += 1
      end
    end
  end
  nface_to_dofs
end

function _extract_nonzeros(mask,values)
  b = Int[]
  for (m,n) in zip(mask,values)
    if (m != 0)
      push!(b, n)
    end
  end
  return Tuple(b)
end

end # module RefFEs


text:
using Compose: LinePrimitive, CirclePrimitive, Form, SimplePolygonPrimitive,
    RectanglePrimitive, TextPrimitive, ArcPrimitive, CurvePrimitive
const PointT = Tuple{Float64,Float64}
const TextT = Tuple{Float64,Float64,String}
const EdgeT = Tuple{PointT, PointT}
const NODE_CACHE = MyOrderedDict{Context, Vector{PointT}}()
const EDGE_CACHE = MyOrderedDict{Context, Vector{EdgeT}}()
const TEXT_CACHE = MyOrderedDict{Context, Vector{TextT}}()
const TextBrush = Pair{Context,String}

export flush!, canvas

function put_node!(brush::Context, x::PointT)
    if haskey(NODE_CACHE, brush)
        push!(NODE_CACHE[brush], x)
    else
        NODE_CACHE[brush] = [x]
    end
    return x
end

function put_edge!(brush::Context, x::PointT, y::PointT)
    if haskey(EDGE_CACHE, brush)
        push!(EDGE_CACHE[brush], (x, y))
    else
        EDGE_CACHE[brush] = [(x, y)]
    end
    return (x, y)
end

function put_text!(brush::Context, x::PointT, y::String)
    if haskey(TEXT_CACHE, brush)
        push!(TEXT_CACHE[brush], (x..., y))
    else
        TEXT_CACHE[brush] = [(x..., y)]
    end
    return (x..., y)
end

function empty_cache!()
    empty!(NODE_CACHE)
    empty!(EDGE_CACHE)
    empty!(TEXT_CACHE)
end

nnode() = isempty(NODE_CACHE) ? 0 : sum(length, values(NODE_CACHE))
nedge() = isempty(EDGE_CACHE) ? 0 : sum(length, values(EDGE_CACHE))
ntext() = isempty(TEXT_CACHE) ? 0 : sum(length, values(TEXT_CACHE))

function inner_most_containers(f, c::Context)
    f(c)
    for child in c.container_children
        inner_most_containers(f, child)
    end
end

function Base.:>>(brush::Context, x::Tuple{<:NTuple{2,T}, <:NTuple{2,T}}) where T<:Real
    put_edge!(brush, Float64.(x[1]), Float64.(x[2]))
end

function Base.:>>(brush::Context, x::Tuple{<:NTuple{2,T}, TS}) where {T<:Real,TS<:AbstractString}
    put_text!(brush, Float64.(x[1]), string(x[2]))
end

function Base.:>>(brush::Context, position::NTuple{2,Real})
    put_node!(brush, Float64.(position))
end

function update_locs!(c, locs)
    isempty(c) && return
    c.head = similar(c.head, locs)
    update_locs!(c.tail, locs)
    return
end

function flush!(d::MyOrderedDict)
    lst = Context[]
    for (brush, lines) in d
        b = deepcopy(brush)
        inner_most_containers(b) do c
            update_locs!(c.form_children, lines)
        end
        push!(lst, b)
    end
    empty!(d)
    return lst
end

function flush!()
    compose(context(),
        flush!(TEXT_CACHE)...,
        flush!(NODE_CACHE)...,
        flush!(EDGE_CACHE)...,
        )
end

"""
    canvas(f)

paint on global canvas.
"""
function canvas(f)
    empty_cache!()
    f()
    flush!()
end

similar(e::Form{<:LinePrimitive}, point_arrays) = line(point_arrays)

function similar(e::Form{<:LinePrimitive}, point_arrays::AbstractVector{Tuple{Float64,Float64}})
    c = first(e.primitives)
    newpoints = map(x->[x .+ _value.(p) for p in c.points], point_arrays)
    line(newpoints)
end

function similar(e::Form{<:ArcPrimitive}, point_arrays)
    c = first(e.primitives)
    as = getindex.(point_arrays,1)
    bs = getindex.(point_arrays,2)
    arc(as, bs, [c.radius], [c.angle1], [c.angle2], [c.sector])
end

function similar(e::Form{<:ArcPrimitive}, point_arrays::AbstractVector{Tuple{Float64,Float64}})
    c = first(e.primitives)
    center = map(x->x .+ _value.(c.center), point_arrays)
    arc(getindex.(center, 1), getindex.(center, 2), [c.radius], [c.angle1], [c.angle2], [c.sector])
end

function similar(e::Form{<:CurvePrimitive}, point_arrays)
    c = first(e.primitives)
    as = getindex.(point_arrays,1)
    bs = getindex.(point_arrays,2)
    ctrl0s = map(x->x[1] .- _value.(c.anchor0) .+ _value.(c.ctrl0), point_arrays)
    ctrl1s = map(x->x[2] .- _value.(c.anchor1) .+ _value.(c.ctrl1), point_arrays)
    curve(as, ctrl0s, ctrl1s, bs)
end

function similar(e::Form{<:CurvePrimitive}, point_arrays::AbstractVector{Tuple{Float64,Float64}})
    c = first(e.primitives)
    anchor0s = map(x->x .+ _value.(c.anchor0), point_arrays)
    anchor1s = map(x->x .+ _value.(c.anchor1), point_arrays)
    ctrl0s = map(x->x .+ _value.(c.ctrl0), point_arrays)
    ctrl1s = map(x->x .+ _value.(c.ctrl1), point_arrays)
    curve(anchor0s, ctrl0s, ctrl1s, anchor1s)
end

function similar(n::Form{<:CirclePrimitive}, xys)
    c = first(n.primitives)
    xs = map(p->_value(c.center[1]) + p[1], xys)
    ys = map(p->_value(c.center[2]) + p[2], xys)
    circle(xs, ys, [c.radius])
end

function similar(p::Form{<:SimplePolygonPrimitive}, xys)
    c = first(p.primitives)
    polygon(map(x-> map(pi->(x[1] + _value(pi[1]), x[2] + _value(pi[2])), c.points), xys))
end

function similar(n::Form{<:RectanglePrimitive}, xys)
    c = first(n.primitives)
    xs = map(p->_value(c.corner[1]) + p[1], xys)
    ys = map(p->_value(c.corner[2]) + p[2], xys)
    rectangle(xs, ys, [c.width], [c.height])
end

function similar(n::Form{<:TextPrimitive}, xys)
    c = first(n.primitives)
    xs = map(p->_value(c.position[1]) + p[1], xys)
    ys = map(p->_value(c.position[2]) + p[2], xys)
    ts = getindex.(xys, 3)
    text(xs, ys, ts, [c.halign], [c.valign], [c.rot], [c.offset])
end

function similar(n::Form{<:MathJaxPrimitive}, xys)
    c = first(n.primitives)
    xs = map(p->_value(c.position[1]) + p[1], xys)
    ys = map(p->_value(c.position[2]) + p[2], xys)
    ts = getindex.(xys, 3)
    mathjax(xs, ys, [c.size[1]], [c.size[2]], ts, [c.rot], [c.offset])
end


text:
using URIParser

struct Address
    scheme::String
    host::String
    port::Integer
end

function from_uri(uri::String, default_port::Integer=0)
    parsed = URI(uri)
    return Address(parsed.scheme, parsed.host, if parsed.port == 0 default_port else parsed.port end)
end

Base.:(==)(x::Address, y::Address) = x.scheme == y.scheme && x.host == y.host && x.port == y.port

max_stars_repo_name: virtualgraham/JuliaBolt

max_issues_repo_name: virtualgraham/JuliaBolt

max_forks_repo_name: virtualgraham/JuliaBolt


text:
abstract type AbstractList{T} end

### Interface

# Insert an item at the current position.
insert!(l::AbstractList, item) = nothing

# Append an item at the end.
append!(l::AbstractList, item) = nothing

# Remove and return the current item.
remove!(l::AbstractList) = nothing

# Set current position to the start/end/position.
movetostart!(l::AbstractList) = nothing
movetoend!(l::AbstractList) = nothing
movetoposition!(l::AbstractList, pos) = nothing

# Move current position one step left/right.
prev!(l::AbstractList) = nothing
next!(l::AbstractList) = nothing

# Return the position of the current element.
currentposition(l::AbstractList) = nothing

# Return the current element.
getvalue(l::AbstractList) = nothing

# Return the number of elements in the list.
length(l::AbstractList) = nothing

# Check if list is empty.
isempty(l::AbstractList) = nothing

### Generic

eltype(l::AbstractList{T}) where T = T


text:
"""
    MeshHexahedronModule

Module for generation of  hexahedral meshes.
"""
module MeshHexahedronModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FESetModule: FESet, FESetQ4, FESetH8, FESetH20, FESetH27, subset, bfun, connasarray, setlabel!, updateconn!
import FinEtools.FENodeSetModule: FENodeSet, count, xyz3
import FinEtools.MeshQuadrilateralModule: Q4elliphole, Q4circlen
import FinEtools.MeshUtilModule: makecontainer, addhyperface!, findhyperface!, linearspace
import FinEtools.MeshModificationModule: meshboundary, mergemeshes
import FinEtools.MeshSelectionModule: selectelem, connectednodes
import LinearAlgebra: norm
import Statistics: mean

"""
    H8block(Length::FFlt, Width::FFlt, Height::FFlt, nL::FInt, nW::FInt, nH::FInt)

Make  a mesh  of a 3D block consisting of  eight node hexahedra.

Length, Width, Height= dimensions of the mesh in Cartesian coordinate axes,
smallest coordinate in all three directions is  0 (origin)
nL, nW, nH=number of elements in the three directions
"""
function H8block(Length::FFlt, Width::FFlt, Height::FFlt, nL::FInt, nW::FInt, nH::FInt)
    return H8blockx(collect(linearspace(0.0, Length, nL+1)), collect(linearspace(0.0, Width, nW+1)), collect(linearspace(0.0, Height, nH+1)));
end

"""
    H8blockx(xs::FFltVec, ys::FFltVec, zs::FFltVec)

Graded mesh of a 3-D block of H8 finite elements.
"""
function H8blockx(xs::FFltVec, ys::FFltVec, zs::FFltVec)
    nL =length(xs)-1;
    nW =length(ys)-1;
    nH =length(zs)-1;

    nnodes=(nL+1)*(nW+1)*(nH+1);
    xyz =zeros(FFlt, nnodes, 3);
    ncells=(nL)*(nW)*(nH);
    conns =zeros(FInt, ncells, 8);

    # first go along Length,  then Width,  and finally Height
    # function  break(i, j, k, nL, nW, nH)
    #     f=(k-1)*((nL+1)*(nW+1))+(j-1)*(nL+1)+i;
    #     nn=[f (f+1)  f+(nL+1)+1 f+(nL+1)];
    #     nn=[nn broadcast(+, nn, (nL+1)*(nW+1))];
    # end

    f=1;
    for k=1:(nH+1)
        for j=1:(nW+1)
            for i=1:(nL+1)
                xyz[f, 1] = xs[i]
                xyz[f, 2] = ys[j]
                xyz[f, 3] = zs[k];
                f=f+1;
            end
        end
    end

    gc=1;
    for i=1:nL
        for j=1:nW
            for k=1:nH
                f=(k-1)*((nL+1)*(nW+1))+(j-1)*(nL+1)+i;
                conns[gc, 1] = f
                conns[gc, 2] = (f+1)
                conns[gc, 3] = f+(nL+1)+1
                conns[gc, 4] = f+(nL+1)
                f = f + (nL+1)*(nW+1)
                conns[gc, 5] = f
                conns[gc, 6] = (f+1)
                conns[gc, 7] = f+(nL+1)+1
                conns[gc, 8] = f+(nL+1)
                gc=gc+1;
            end
        end
    end
    # create the nodes
    fens = FENodeSet(xyz);
    # Create the finite elements
    fes = FESetH8(conns);

    return fens, fes;
end

"""
    H8sphere(radius::FFlt, nrefine::FInt)

    Create a mesh of 1/8 of the sphere of "radius". The  mesh will consist of
    four hexahedral elements if "nrefine==0",  or more if "nrefine>0".
    "nrefine" is the number of bisections applied  to refine the mesh.
"""
function H8sphere(radius::FFlt, nrefine::FInt)
    a=sqrt(2.0)/2.0;
    b=1.0/sqrt(3.0);
    c=0.6*a;
    d=0.6*b;
    z= 0.0;
    h= 0.5;
    o= 1.0;
    xyz= [z  z  z;
        h  z  z;
        c  c  z;
        z  h  z;
        z  z  h;
        c  z  c;
        d  d  d;
        z  c  c;
        z  z  o;
        a  z  a;
        o  z  z;
        a  a  z;
        z  o  z;
        z  a  a;
        b  b  b]*radius;
    conns=[1  2  3  4  5  6  7  8;
        2  11  12  3  6  10  15  7;
        4  3  12  13  8  7  15  14;
        5  6  7  8  9  10  15  14];

    fens = FENodeSet(xyz);
    fes = FESetH8(conns);
    for i = 1:nrefine
        fens, fes = H8refine(fens, fes);
        bg = meshboundary(fes);
        l = selectelem(fens, bg, facing=true,
            direction=[1, 1, 1], dotmin= 0.01);
        cn = connectednodes(subset(bg, l))   ;
        for j=1:length(cn)
            fens.xyz[cn[j], :]=fens.xyz[cn[j], :]*radius/norm(fens.xyz[cn[j], :]);
        end
    end
    return fens,  fes
end

"""
    H8refine(fens::FENodeSet,  fes::FESetH8)

Refine a mesh of H8 hexahedrals by octasection.
"""
function H8refine(fens::FENodeSet,  fes::FESetH8)
    fens, fes = H8toH27(fens, fes);
    conn=connasarray(fes);
    nconn=zeros(FInt, 8*size(conn, 1), 8);
    nc=1;
    for i= 1:size(conn, 1)
        conn27=conn[i, :];
        nconn[nc, :] =conn27[[1, 9, 21, 12, 17, 22, 27, 25]];        nc= nc+ 1;
        nconn[nc, :] =conn27[[9, 2, 10, 21, 22, 18, 23, 27]];        nc= nc+ 1;
        nconn[nc, :] =conn27[[21, 10, 3, 11, 27, 23, 19, 24]];        nc= nc+ 1;
        nconn[nc, :] =conn27[[12, 21, 11, 4, 25, 27, 24, 20]];        nc= nc+ 1;
        nconn[nc, :] =conn27[[17, 22, 27, 25, 5, 13, 26, 16]];        nc= nc+ 1;
        nconn[nc, :] =conn27[[22, 18, 23, 27, 13, 6, 14, 26]];        nc= nc+ 1;
        nconn[nc, :] =conn27[[27, 23, 19, 24, 26, 14, 7, 15]];        nc= nc+ 1;
        nconn[nc, :] =conn27[[25, 27, 24, 20, 16, 26, 15, 8]];        nc= nc+ 1;
    end

    fes=FESetH8(nconn);
    return fens,  fes
end

"""
    H8toH27(fens::FENodeSet,  fes::FESetH8)

    Convert a mesh of hexahedra H8 to hexahedra H27.
"""
function   H8toH27(fens::FENodeSet,  fes::FESetH8)
    nedges=12;
    nfaces=6;
    ec = [1   2; 2   3; 3   4; 4   1; 5   6; 6   7; 7   8; 8   5; 1   5; 2   6; 3   7; 4   8;];
    fc = [1     4     3     2;
        1     2     6     5;
        2     3     7     6;
        3     4     8     7;
        4     1     5     8;
        6     7     8     5];
    conns = connasarray(fes);
    labels = deepcopy(fes.label)
    # Additional node numbers are numbered from here
    newn=count(fens)+1;
    # make a search structure for edges
    edges=makecontainer();
    for i= 1:size(conns, 1)
        conn = conns[i, :];
        for J = 1:nedges
            ev=conn[ec[J, :]];
            newn = addhyperface!(edges,  ev,  newn);
        end
    end
    # make a search structure for faces
    faces=makecontainer();
    for i= 1:size(conns, 1)
        conn = conns[i, :];
        for J = 1:nfaces
            fv=conn[fc[J, :]];
            newn = addhyperface!(faces,  fv,  newn);
        end
    end
    # make a search structure for volumes
    volumes=makecontainer();
    for i= 1:size(conns, 1)
      conn = conns[i, :];
      newn = addhyperface!(volumes,  conn,  newn);
    end
    xyz1 =fens.xyz;             # Pre-existing nodes
   # Allocate for vertex nodes plus edge nodes plus face nodes
    xyz =zeros(FFlt, newn-1, 3);
    xyz[1:size(xyz1, 1), :] = xyz1; # existing nodes are copied over
    # calculate the locations of the new nodes
    # and construct the new nodes
    for i in keys(edges)
        C=edges[i];
        for J = 1:length(C)
          ix = vec([item for item in C[J].o])
          push!(ix,  i)
          xyz[C[J].n, :]=mean(xyz[ix, :], dims = 1);
        end
    end
   # calculate the locations of the new nodes
    # and construct the new nodes
    for i in keys(faces)
        C=faces[i];
        for J = 1:length(C)
          ix = vec([item for item in C[J].o])
          push!(ix,  i)
          xyz[C[J].n, :] = mean(xyz[ix, :], dims = 1);
        end
    end
    # calculate the locations of the new nodes
    # and construct the new nodes
    for i in keys(volumes)
       C=volumes[i];
        for J = 1:length(C)
          ix = vec([item for item in C[J].o])
          push!(ix,  i)
          xyz[C[J].n, :] = mean(xyz[ix, :], dims = 1);
        end
    end
     # construct new geometry cells
    nconns =zeros(FInt, size(conns, 1), 27);
    nc=1;
    for i= 1:size(conns, 1)
        conn = conns[i, :];
        econn=zeros(FInt, 1, nedges);
        for J = 1:nedges
            ev=conn[ec[J, :]];
            h, n=findhyperface!(edges,  ev);
            econn[J]=n;
        end
        fconn=zeros(FInt, 1, nfaces);
        for J = 1:nfaces
            fv=conn[fc[J, :]];
            h, n=findhyperface!(faces,  fv);
            fconn[J]=n;
        end
        h, n=findhyperface!(volumes,  conn);
        vconn=n;
        nconns[nc, :] =vcat(vec(conn),  vec(econn),  vec(fconn),  vec([vconn]))
        nc= nc+ 1;
    end
    fens =FENodeSet(xyz);
    fes = FESetH27(nconns) ;
    setlabel!(fes, labels);
    return fens, fes;
end

"""
    H8hexahedron(xyz::FFltMat, nL::FInt, nW::FInt, nH::FInt; blockfun=nothing)

Mesh of a general hexahedron given by the location of the vertices.

`xyz` = One vertex location per row; Either two rows (for a rectangular
     block given by the its corners),  or eight rows (general hexahedron).
`nL`,  `nW`,  `nH` = number of elements in each direction
`blockfun` = Optional argument: function of the block-generating mesh function
     (having the signature of the function `H8block()`).
"""
function H8hexahedron(xyz::FFltMat, nL::FInt, nW::FInt, nH::FInt; blockfun=nothing)
    npts=size(xyz, 1);
    @assert (npts == 2) || (npts == 8) "Need 2 or 8 points"
    if npts == 2
        lo=minimum(xyz, dims = 1);
        hi=maximum(xyz, dims = 1);
        xyz=[lo[1]  lo[2]  lo[3];
            hi[1]  lo[2]  lo[3];
            hi[1]  hi[2]  lo[3];
            lo[1]  hi[2]  lo[3];
            lo[1]  lo[2]  hi[3];
            hi[1]  lo[2]  hi[3];
            hi[1]  hi[2]  hi[3];
            lo[1]  hi[2]  hi[3]];
    end

    if blockfun == nothing
        blockfun = H8block; # default block type
    end

    fens, fes= blockfun(2.0, 2.0, 2.0, nL, nW, nH);

    dummy = FESetH8(reshape(collect(1:8), 1, 8))
    pxyz=fens.xyz;
    for i=1:count(fens)
        N = bfun(dummy, broadcast(-, pxyz[i, :], 1.0));# shift coordinates by -1
        pxyz[i, :] =N'*xyz;
    end
    fens.xyz = pxyz;

    return fens, fes;
end

"""
    H27block(Length::FFlt, Width::FFlt, Height::FFlt, nL::FInt, nW::FInt, nH::FInt)

Create mesh of a 3-D block of H27 finite elements.
"""
function H27block(Length::FFlt, Width::FFlt, Height::FFlt, nL::FInt, nW::FInt, nH::FInt)
    fens, fes = H8block(Length, Width, Height, nL, nW, nH);
    fens, fes = H8toH27(fens, fes);
    return fens, fes
end

"""
    H27blockx(xs::FFltVec, ys::FFltVec, zs::FFltVec)

Graded mesh of a 3-D block of H27 finite elements.
"""
function H27blockx(xs::FFltVec, ys::FFltVec, zs::FFltVec)
    fens, fes = H8blockx(xs, ys, zs);
    fens, fes = H8toH27(fens, fes);
end

function doextrude(fens, fes::FESetQ4, nLayers, extrusionh)
    nn1=count(fens);
    nnt=nn1*nLayers;
    ngc=count(fes)*nLayers;
    hconn=zeros(FInt, ngc, 8);
    conn = connasarray(fes)
    xyz =zeros(FFlt, nn1*(nLayers+1), 3);
    for j=1:nn1
        xyz[j, :] =extrusionh(fens.xyz[j, :], 0);
    end
    for k=1:nLayers
        for j=1:nn1
            f=j+k*nn1;
            xyz[f, :] =extrusionh(fens.xyz[j, :], k);
        end
    end

    gc=1;
    for k=1:nLayers
        for i=1:count(fes)
            hconn[gc, :]=[broadcast(+, conn[i, :], (k-1)*nn1) broadcast(+, conn[i, :], k*nn1)];
            gc=gc+1;
        end
    end
    efes = FESetH8(hconn);
    efens = FENodeSet(xyz);
    return efens, efes
end

"""
    H8extrudeQ4(fens::FENodeSet,  fes::FESetQ4, nLayers::FInt, extrusionh::Function)

Extrude a mesh of quadrilaterals into a mesh of hexahedra (H8).
"""
function H8extrudeQ4(fens::FENodeSet,  fes::FESetQ4, nLayers::FInt, extrusionh::F) where {F<:Function}
    id = vec([i for i in 1:count(fens)])
    cn=connectednodes(fes);
    id[cn[:]]=vec([i for i in 1:length(cn)]);
    q4fes= deepcopy(fes);
    updateconn!(q4fes, id);
    q4fens = FENodeSet(fens.xyz[cn[:], :]);
    h8fens, h8fes= doextrude(q4fens, q4fes, nLayers, extrusionh);
    return h8fens, h8fes
end

"""
    H8spheren(radius::FFlt, nperradius::FInt)

Create a solid mesh of 1/8 of sphere.

Create a solid mesh of 1/8 of the sphere of "radius",  with nperradius
elements per radius.

"""
function H8spheren(radius::FFlt, nperradius::FInt)
    if (mod(nperradius, 2) != 0)
        nperradiu = nperradius+1;
    end
    nL = ceil(FInt, nperradius/2); nW=nL; nH=nL;

    a=sqrt(2.0)/2.0;
    b=1.0/sqrt(3.0);
    c=0.6*a;
    d=0.6*b;
    z= 0.0;
    h= 0.5;
    o= 1.0;
    xyz =  [z  z  z;
            h  z  z;
            c  c  z;
            z  h  z;
            z  z  h;
            c  z  c;
            d  d  d;
            z  c  c;
            z  z  o;
            a  z  a;
            o  z  z;
            a  a  z;
            z  o  z;
            z  a  a;
            b  b  b]*radius;
    conns = [1  2  3  4  5  6  7  8;
            2  11  12  3  6  10  15  7;
            4  3  12  13  8  7  15  14;
            5  6  7  8  9  10  15  14];
    tolerance=radius*1.0e-6;

    fens, fes = H8hexahedron(xyz[conns[1, :][:], :], nL, nW, nH);
    fens1, fes1 = H8hexahedron(xyz[conns[2, :][:], :], nL, nW, nH);
    fens, fes1, fes2 = mergemeshes(fens1,  fes1,  fens,  fes,  tolerance);
    fes=cat(fes1, fes2);
    fens1, fes1 = H8hexahedron(xyz[conns[3, :][:], :], nL, nW, nH);
    fens, fes1, fes2 = mergemeshes(fens1,  fes1,  fens,  fes,  tolerance);
    fes=cat(fes1, fes2);
    fens1, fes1 = H8hexahedron(xyz[conns[4, :][:], :], nL, nW, nH);
    fens, fes1, fes2 = mergemeshes(fens1,  fes1,  fens,  fes,  tolerance);
    fes=cat(fes1, fes2);

    xyz = deepcopy(fens.xyz);
    layer = broadcast(+, zeros(FFlt, size(xyz,  1), 1), oftype(1.0, Inf));
    conn = deepcopy(connasarray(fes));
    bg = meshboundary(fes);
    l = selectelem(fens, bg; facing=true,  direction=[1. 1. 1.], dotmin = 0.01);
    cn = connectednodes(subset(bg, l))   ;
    layer[cn] .= 1;
    for j = 1:nperradius-1
        for k = 1:size(conn, 1)
            ll = layer[conn[k, :]];
            ml = minimum(ll);
            if (ml==j)
                ix = isinf.(ll);
                ll[ix] .= j+1;
                layer[conn[k, :]] = ll;
            end
        end
    end
    nxyz = deepcopy(xyz);
    for j = 1:size(xyz, 1)
        if (!isinf.(layer[j]))
            nxyz[j, :] = nxyz[j, :]*(nperradius-layer[j]+1)/nperradius*radius/norm(nxyz[j, :]);
        end
    end
    s =  collect(linearspace(0.,  1., nperradius));
    for j = 1:size(xyz, 1)
        ell = layer[j]
        # show(ell)
        if (!isinf.(ell))
            ell = Int(ell)
            nxyz[j, :] = s[ell]*xyz[j, :] + (1-s[ell])*nxyz[j, :];
        end
    end
    fens.xyz = deepcopy(nxyz);
    return fens, fes
end

"""
    H20block(Length::FFlt, Width::FFlt, Height::FFlt, nL::FInt, nW::FInt, nH::FInt)

Create mesh of a 3-D block of H20 finite elements.
"""
function H20block(Length::FFlt, Width::FFlt, Height::FFlt, nL::FInt, nW::FInt, nH::FInt)
    fens, fes = H8block(Length, Width, Height, nL, nW, nH);
    fens, fes = H8toH20(fens, fes);
end

"""
    H20blockx(xs::FFltVec, ys::FFltVec, zs::FFltVec)

Graded mesh of a 3-D block of H20 finite elements.
"""
function H20blockx(xs::FFltVec, ys::FFltVec, zs::FFltVec)
    fens, fes = H8blockx(xs, ys, zs);
    fens, fes = H8toH20(fens, fes);
end

"""
    H8toH20(fens::FENodeSet,  fes::FESetH8)

Convert a mesh of hexahedra H8 to hexahedra H20.
"""
function   H8toH20(fens::FENodeSet,  fes::FESetH8)
    nedges=12;
    ec = [1   2; 2   3; 3   4; 4   1; 5   6; 6   7; 7   8; 8   5; 1   5; 2   6; 3   7; 4   8;];
    conns = connasarray(fes);
    labels = deepcopy(fes.label)
    # Additional node numbers are numbered from here
    newn=count(fens)+1;
    # make a search structure for edges
    edges=makecontainer();
    for i= 1:size(conns, 1)
        conn = conns[i, :];
        for J = 1:nedges
            ev=conn[ec[J, :]];
            newn = addhyperface!(edges,  ev,  newn);
        end
    end
    xyz1 =fens.xyz;             # Pre-existing nodes
    # Allocate for vertex nodes plus edge nodes plus face nodes
    xyz =zeros(FFlt, newn-1, 3);
    xyz[1:size(xyz1, 1), :] = xyz1; # existing nodes are copied over
    # calculate the locations of the new nodes
    # and construct the new nodes
    for i in keys(edges)
        C=edges[i];
        for J = 1:length(C)
            ix = vec([item for item in C[J].o])
            push!(ix,  i) # Add the anchor point as well
            xyz[C[J].n, :] = mean(xyz[ix, :], dims = 1);
        end
    end
    # construct new geometry cells
    nconns =zeros(FInt, size(conns, 1), 20);
    nc=1;
    for i= 1:size(conns, 1)
        conn = conns[i, :];
        econn=zeros(FInt, 1, nedges);
        for J = 1:nedges
            ev=conn[ec[J, :]];
            h, n=findhyperface!(edges,  ev);
            econn[J]=n;
        end
        nconns[nc, :] =vcat(vec(conn),  vec(econn))
        nc= nc+ 1;
    end
    fens =FENodeSet(xyz);
    fes = FESetH20(nconns) ;
    setlabel!(fes, labels);
    return fens, fes;
end

# Construct arrays to describe a hexahedron mesh created from voxel image.
#
# img = 3-D image (array),  the voxel values  are arbitrary
# voxval =range of voxel values to be included in the mesh,
# voxval =  [minimum value,  maximum value].  Minimum value == maximum value is
# allowed.
# Output:
# t = array of hexahedron connectivities,  one hexahedron per row
# v =Array of vertex locations,  one vertex per row
function H8voximggen(img::Array{DataT, 3},
  voxval::Array{DataT, 1}) where {DataT<:Number}
   M=size(img,  1); N=size(img,  2); P=size(img,  3);

     function find_nonempty(minvoxval, maxvoxval)
        Nvoxval=0
        for I= 1:M
            for J= 1:N
                for K= 1:P
                    if (img[I, J, K]>=minvoxval) && (img[I, J, K]<=maxvoxval)
                        Nvoxval=Nvoxval+1
                    end
                end
            end
        end
        return Nvoxval
    end
    minvoxval= minimum(voxval)  # include voxels at or above this number
    maxvoxval= maximum(voxval)  # include voxels at or below this number
    Nvoxval =find_nonempty(minvoxval, maxvoxval) # how many "full" voxels are there?

    # Allocate output arrays
    h =zeros(FInt, Nvoxval, 8);
    v =zeros(FInt, (M+1)*(N+1)*(P+1), 3);
    hmid =zeros(FInt, Nvoxval);

    Slice =zeros(FInt, 2, N+1, P+1); # auxiliary buffer
    function find_vertex(I, IJK)
        vidx = zeros(FInt, 1, size(IJK, 1));
        for r= 1:size(IJK, 1)
            if (Slice[IJK[r, 1], IJK[r, 2], IJK[r, 3]]==0)
                nv=nv+1;
                v[nv, :] =IJK[r, :]; v[nv, 1] += I-1;
                Slice[IJK[r, 1], IJK[r, 2], IJK[r, 3]] =nv;
            end
            vidx[r] =Slice[IJK[r, 1], IJK[r, 2], IJK[r, 3]];
        end
        return vidx
    end
    function store_hex(I, J, K)
        locs =[1 J K;1+1 J K;1+1 J+1 K;1 J+1 K;1 J K+1;1+1 J K+1;1+1 J+1 K+1;1 J+1 K+1];
        vidx = find_vertex(I, locs);
        nh =nh +1;
        h[nh, :] = vidx;
        hmid[nh] = convert(FInt, img[I, J, K]);
    end

    nv =0;                      # number of vertices
    nh =0;                      # number of elements
    for I= 1:M
        for J= 1:N
            for K= 1:P
                if  (img[I, J, K]>=minvoxval) && (img[I, J, K]<=maxvoxval)
                    store_hex(I, J, K);
                end
            end
        end
        Slice[1, :, :] = Slice[2, :, :] ;
        Slice[2, :, :] .= 0;
    end
    # Trim output arrays
    v=v[1:nv, :];
    h=h[1:nh, :] ;
    hmid=hmid[1:nh] ;

    return h, v, hmid
end

"""
    H8voximg(img::Array{DataT, 3}, voxdims::FFltVec,
        voxval::Array{DataT, 1}) where {DataT<:Number}

Generate a hexahedral mesh  from three-dimensional image.
"""
function H8voximg(img::Array{DataT, 3}, voxdims::FFltVec,
    voxval::Array{DataT, 1}) where {DataT<:Number}
    h, v, hmid = H8voximggen(img, voxval)
    xyz = zeros(FFlt, size(v, 1), 3)
    for k=1:3
        for j=1:size(v, 1)
            xyz[j, k] = v[j, k]*voxdims[k]
        end
    end
    fens  = FENodeSet(xyz);
    fes = FESetH8(h);
    setlabel!(fes, hmid)
    return fens, fes;
end

"""
    H8layeredplatex(xs::FFltVec, ys::FFltVec, ts::FFltVec, nts::FIntVec)

H8 mesh for a layered block (composite plate) with specified in plane coordinates.

xs,ys =Locations of the individual planes of nodes.
ts= Array of layer thicknesses,
nts= array of numbers of elements per layer

The finite elements of each layer are labeled with the layer number, starting
from 1.
"""
function H8layeredplatex(xs::FFltVec, ys::FFltVec, ts::FFltVec, nts::FIntVec)
    tolerance = minimum(abs.(ts))/maximum(nts)/10.;
    @assert length(ts) >= 1
    layer = 1
    zs = collect(linearspace(0.0,ts[layer],nts[layer]+1))
    fens, fes = H8blockx(xs, ys, zs);
    setlabel!(fes, layer);
    for layer = 2:length(ts)
        zs = collect(linearspace(0.0,ts[layer],nts[layer]+1))
        fens1, fes1 = H8blockx(xs, ys, zs);
        setlabel!(fes1, layer);
        fens1.xyz[:, 3] = fens1.xyz[:, 3] .+ sum(ts[1:layer-1]);
        fens, fes1, fes2 = mergemeshes(fens1, fes1, fens, fes, tolerance);
        fes = cat(fes1,fes2);
    end
    return fens,fes
end

"""
    H8elliphole(xradius::FFlt, yradius::FFlt, L::FFlt, H::FFlt, T::FFlt,
        nL::FInt, nH::FInt, nW::FInt, nT::FInt)

Mesh of one quarter of a rectangular plate with an elliptical hole.

`xradius`,`yradius` = radii of the ellipse,
`L`,`H` = dimensions of the plate,
`T` = thickness of the plate
`nL`,`nH`= numbers of edges along the side of the plate; this is also
  the number of edges along the circumference of the elliptical hole
`nW` = number of edges along the remaining straight edge (from the hole
  in the radial direction),
"""
function H8elliphole(xradius::FFlt, yradius::FFlt, L::FFlt, H::FFlt, T::FFlt,
    nL::FInt, nH::FInt, nW::FInt, nT::FInt)
    fens,fes  =  Q4elliphole(xradius, yradius, L, H, nL, nH, nW)
    fens.xyz = xyz3(fens)
    ex(xyz, layer) = xyz + reshape(layer/nT*T*[0.0 0.0 1.0], size(xyz, 1), size(xyz, 2))
    fens,fes  =  H8extrudeQ4(fens, fes, nT, ex);
    return fens,fes
end

"""
    H8cylindern(Radius::FFlt, Length::FFlt, nperradius, nL)

H8 mesh of a solid  cylinder with given number of edges per radius (`nperradius`)
and per length (`nL`).
"""
function H8cylindern(Radius::FFlt, Length::FFlt, nperradius, nL)
    tol = min((Length/nL), (Radius/2/nperradius))/100;
    fens,fes = Q4circlen(Radius, nperradius);
    fens2 = deepcopy(fens)
    for i = 1:count(fens2)
    	fens2.xyz[i, 1], fens2.xyz[i, 2] = -fens2.xyz[i, 2], fens2.xyz[i, 1]
    end
    fens,fes,fes2 = mergemeshes(fens, fes, fens2, fes, tol);
    fes = cat(fes,fes2);
    fens2 = deepcopy(fens)
    for i = 1:count(fens2)
    	fens2.xyz[i, 1], fens2.xyz[i, 2] = -fens2.xyz[i, 1], -fens2.xyz[i, 2]
    end
    fens,fes,fes2 = mergemeshes(fens, fes, fens2, fes, tol);
    fes = cat(fes,fes2);
    fens,fes = H8extrudeQ4(fens,fes,nL,(x,k) -> [x[1],x[2],k*Length/nL]);
    return fens,fes
end

end


text:
# numerical stochastic climate-carbon cycle model
# from Arnscheidt and Rothman - Asymmetry of extreme Cenozoic climate-carbon cycle events (2021)

using DifferentialEquations
using Plots
using BenchmarkTools
using DelimitedFiles

# unit convention: SI, except for mass (Pg) and time (years)

# PARAMETERS

# 1. carbonate chemistry approximation parameters
const γ = 6.5
const Iᵦ= 58000.0

# 2. isotopic parameters

const vpdb = 0.01123720 	# ratio of 13C/12C in VPDB standard
const δᵪ = 1.0 			# isotopic composition of carbonate carbon
const Rχ = (δᵪ/1000 + 1)*vpdb 	# ratio of 13C/12C in carbonate carbon
const δₒ = -25
const Rₒ = (δₒ/1000 + 1)*vpdb

# 4. other parameters
const Myr = 1.0*10^6
const kyr = 1000.0
const sinyr = 365.25*24*3600.0	# seconds in a year
const P₀ = 400.0*10^(-6)  	# pCO₂ in reference state
const I₀ = 38000.0		# reference surficial carbon stock
const τ = 100kyr		# timescale of weathering feedback
const λ = 5.0 			# Earth System Sensitivity (K)
const a = 2.2*sinyr		# slope of OLR with T (J/yr/m²/K)
const C = 2.0*10^(8)		# heat capacity of the Earth (J/m²/K)

end

# FUNCTION SPECIFICATION

function model!(du,u,par,t) # the model equations

	# define explicitly for ease of understanding
	ΔT = u[1]
	I = u[2]
	I₁₃ = u[3]

	if I<0
		I = 0.0
	end

	P = ((I^γ)/(I^γ+Iᵦ^γ))*10^(-6)*7000.0

	du[1] = (1/C)*(-a*(ΔT-par[5](t))+a*λ*log(P/P₀)/log(2))
	du[2] = - (I-I₀)/τ
	du[3] = - (I₁₃-I₀*Rχ)/τ
	du[4] = 0

end

function noise!(du,u,par,t)

	du[1,1] = par[1]*(u[1]+par[3])
	du[2,1] = par[2]*(u[1]+par[3])
	du[3,1] = par[2]*(u[1]+par[3])*Rₒ
	du[4,1] = par[2]*(u[1]+par[3])

	du[1,2] = par[4]
	du[2,2] = 0
	du[3,2] = 0
	du[4,2] = 0
end

# ENSEMBLE RUN - periodic forcing
u0 = [0.0,I₀,Rχ*I₀,0.0]
tspan = (0.0,10Myr)
tstep = 10000
t = range(tspan[1],tspan[2],step=tstep)

F = t->0.6(sin(2*π*t/(400000)))

νₜ= 0.2
νᵪ= 1.0
c = 1.0
μ = 0.4

prob_st = SDEProblem(model!,noise!,u0,tspan,[νₜ,νᵪ,c,μ,t->F(t)],noise_rate_prototype=zeros(4,2))
ensemble_prob = EnsembleProblem(prob_st)
sol = solve(ensemble_prob,dt=1,saveat=tstep,EM(),trajectories=10)

# convert sol[3,:,:] to be time series of d13c
sol[3,:,:] = ((sol[3,:,:]./sol[2,:,:])/vpdb.-1)*1000

# write outputs to file
#writedlm("ens_400_T.csv",sol[1,:,:],',')
#writedlm("ens_400_dc.csv",sol[3,:,:],',')