text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module Graphics

import Base: +, -, *, /, &, fill, norm, scale
using Colors

export
    # Part 1. 2D Geometry
    Vec2, Point, BoundingBox,
    # limits in world coordinates
    isinside, xmin, xmax, ymin, ymax, center, xrange, yrange,
    aspect_ratio, with_aspect_ratio, diagonal, shift, deform,
    # TODO: more

    # TODO: 3D geometry

    # Part 2. 2D Drawing
    # device and context
    GraphicsDevice, GraphicsContext, creategc, getgc,
    # width, height are in user (world) coordinates for geometric objects,
    # but in device coordinates for GraphicsDevice, GraphicsContext, and
    # other concrete things like windows and widgets.
    width, height,

    # drawing attribute manipulation
    save, restore, set_line_width, set_dash, set_source_rgb, set_source_rgba,
    set_source,

    # coordinate systems
    reset_transform, set_coords, rotate, scale, translate, user_to_device!,
    device_to_user!, user_to_device_distance!, device_to_user_distance!,
    user_to_device, device_to_user,

    # clipping
    clip, clip_preserve, reset_clip,

    # path primitives
    move_to, line_to, rel_line_to, rel_move_to, new_path, new_sub_path,
    close_path, arc,

    # fill and stroke
    fill, fill_preserve, paint, stroke, stroke_preserve,
    stroke_transformed, stroke_transformed_preserve,

    # derived path operations
    rectangle, circle, polygon

    # TODO: text drawing API

    # TODO: rendering pipeline API

# Part 1. geometric primitives

immutable Vec2
    x::Float64
    y::Float64
end

typealias Point Vec2

(+)(a::Vec2, b::Vec2) = Vec2(a.x + b.x, a.y + b.y)
(-)(a::Vec2, b::Vec2) = Vec2(a.x - b.x, a.y - b.y)
(*)(p::Vec2, s::Real) = Vec2(p.x*s, p.y*s)
(/)(p::Vec2, s::Real) = Vec2(p.x/s, p.y/s)
(*)(s::Real, p::Vec2) = p*s

# rotate p around o by angle
function rotate(p::Vec2, angle::Real, o::Vec2)
    c = cos(angle)
    s = sin(angle)
    d = p - o
    Vec2(o.x + c*d.x - s*d.y, o.y + s*d.x + c*d.y)
end
rotate(p::Vec2, angle::Real) = rotate(p, angle, Vec2(0.,0.))

norm(p::Vec2) = hypot(p.x, p.y)

immutable BoundingBox
    xmin::Float64
    xmax::Float64
    ymin::Float64
    ymax::Float64
end

BoundingBox() = BoundingBox(NaN, NaN, NaN, NaN)

function BoundingBox(points::Point...)
    xmin, xmax, ymin, ymax = NaN, NaN, NaN, NaN
    for p in points
        xmin = min(xmin, p.x)
        xmax = max(xmax, p.x)
        ymin = min(ymin, p.y)
        ymax = max(ymax, p.y)
    end
    return BoundingBox(xmin, xmax, ymin, ymax)
end

function BoundingBox(bboxes::BoundingBox...)
    xmin, xmax, ymin, ymax = NaN, NaN, NaN, NaN
    for bb in bboxes
        xmin = min(xmin, bb.xmin)
        xmax = max(xmax, bb.xmax)
        ymin = min(ymin, bb.ymin)
        ymax = max(ymax, bb.ymax)
    end
    return BoundingBox(xmin, xmax, ymin, ymax)
end

width(bb::BoundingBox) = bb.xmax - bb.xmin
height(bb::BoundingBox) = bb.ymax - bb.ymin
diagonal(bb) = hypot(width(bb), height(bb))
aspect_ratio(bb) = height(bb)/width(bb)

xmin(bb::BoundingBox) = bb.xmin
xmax(bb::BoundingBox) = bb.xmax
ymin(bb::BoundingBox) = bb.ymin
ymax(bb::BoundingBox) = bb.ymax

center(x) = Point((xmin(x)+xmax(x))/2, (ymin(x)+ymax(x))/2)
xrange(x) = xmin(x), xmax(x)
yrange(x) = ymin(x), ymax(x)

function (+)(bb1::BoundingBox, bb2::BoundingBox)
    BoundingBox(min(bb1.xmin, bb2.xmin),
                max(bb1.xmax, bb2.xmax),
                min(bb1.ymin, bb2.ymin),
                max(bb1.ymax, bb2.ymax))
end

function (&)(bb1::BoundingBox, bb2::BoundingBox)
    BoundingBox(max(bb1.xmin, bb2.xmin),
                min(bb1.xmax, bb2.xmax),
                max(bb1.ymin, bb2.ymin),
                min(bb1.ymax, bb2.ymax))
end

function deform(bb::BoundingBox, dl, dr, dt, db)
    BoundingBox(bb.xmin + dl, bb.xmax + dr, bb.ymin + dt, bb.ymax + db)
end

# shift center by (dx,dy), keeping width & height fixed
function shift(bb::BoundingBox, dx, dy)
    BoundingBox(bb.xmin + dx, bb.xmax + dx, bb.ymin + dy, bb.ymax + dy)
end

# scale width & height, keeping center fixed
function (*)(bb::BoundingBox, s::Real)
    dw = 0.5*(s - 1)*width(bb)
    dh = 0.5*(s - 1)*height(bb)
    deform(bb, -dw, dw, -dh, dh)
end
(*)(s::Real, bb::BoundingBox) = bb*s

function rotate(bb::BoundingBox, angle::Real, p::Point)
    a = rotate(Point(bb.xmin,bb.ymin), angle, p)
    b = rotate(Point(bb.xmax,bb.ymin), angle, p)
    c = rotate(Point(bb.xmin,bb.ymax), angle, p)
    d = rotate(Point(bb.xmax,bb.ymax), angle, p)
    BoundingBox(a, b, c, d)
end

function with_aspect_ratio(bb::BoundingBox, ratio::Real)
    if ratio < aspect_ratio(bb)
        dh = height(bb) - ratio * width(bb)
        return BoundingBox(bb.xmin, bb.xmax, bb.ymin + dh/2, bb.ymax - dh/2)
    else
        dw = width(bb) - height(bb) / ratio
        return BoundingBox(bb.xmin + dw/2, bb.xmax - dw/2, bb.ymin, bb.ymax)
    end
end

isinside(bb::BoundingBox, x, y) = (bb.xmin <= x <= bb.xmax) && (bb.ymin <= y <= bb.ymax)
isinside(bb::BoundingBox, p::Point) = isinside(bb, p.x, p.y)

# Part 2. Drawing

macro mustimplement(sig)
    fname = sig.args[1]
    arg1 = sig.args[2]
    if isa(arg1,Expr)
        arg1 = arg1.args[1]
    end
    :($(esc(sig)) = error(typeof($(esc(arg1))),
                          " must implement ", $(Expr(:quote,fname))))
end

# a graphics output device; can create GraphicsContexts
abstract GraphicsDevice

@mustimplement width(gd::GraphicsDevice)
@mustimplement height(gd::GraphicsDevice)
@mustimplement creategc(gd::GraphicsDevice)
xmin(g::GraphicsDevice) = 0
xmax(g::GraphicsDevice) = width(g)
ymin(g::GraphicsDevice) = 0
ymax(g::GraphicsDevice) = height(g)

# an object that can actually be drawn to
abstract GraphicsContext

@mustimplement width(gc::GraphicsContext)
@mustimplement height(gc::GraphicsContext)
# getgc() - get a GraphicsContext from something that might be drawable
getgc(gc::GraphicsContext) = gc

# transformations

# set coordinates in terms of left, right, top, bottom
function set_coords(c::GraphicsContext, x, y, w, h, l, r, t, b)
    reset_transform(c)
    reset_clip(c)
    rectangle(c, x, y, w, h)
    clip(c)
    if (r-l) != w || (b-t) != h || l != x || t != y
        # note: Cairo assigns integer pixel-space coordinates to the grid
        # points between sample locations, not to the centers of pixels.
        xs = w/(r-l)
        ys = h/(b-t)
        scale(c, xs, ys)

        xcent = (l+r)/2
        ycent = (t+b)/2
        translate(c, -xcent + (w/2 + x)/xs, -ycent + (h/2 + y)/ys)
    end
    c
end
set_coords(c::GraphicsContext, device::BoundingBox, user::BoundingBox) = set_coords(c, device.xmin, device.ymin, width(device), height(device), user.xmin, user.xmax, user.ymin, user.ymax)

@mustimplement save(gc::GraphicsContext)
@mustimplement restore(gc::GraphicsContext)
@mustimplement reset_transform(gc::GraphicsContext)
@mustimplement rotate(gc::GraphicsContext, ::Real)
@mustimplement scale(gc::GraphicsContext, ::Real, ::Real)
@mustimplement translate(gc::GraphicsContext, ::Real, ::Real)

user_to_device!(gc::GraphicsContext, c::Vector{Float64}) = c
device_to_user!(gc::GraphicsContext, c::Vector{Float64}) = c
user_to_device_distance!(gc::GraphicsContext, c::Vector{Float64}) = c
device_to_user_distance!(gc::GraphicsContext, c::Vector{Float64}) = c

const d2ubuf = zeros(2)
function device_to_user(gc::GraphicsContext, x::Real, y::Real)
    d2ubuf[1] = x
    d2ubuf[2] = y
    device_to_user!(gc, d2ubuf)
    d2ubuf[1], d2ubuf[2]
end
function user_to_device(gc::GraphicsContext, x::Real, y::Real)
    d2ubuf[1] = x
    d2ubuf[2] = y
    user_to_device!(gc, d2ubuf)
    d2ubuf[1], d2ubuf[2]
end

# drawing and properties

@mustimplement set_line_width(gc::GraphicsContext, ::Real)
@mustimplement set_dash(gc::GraphicsContext, ::Vector{Float64}, ::Real)
@mustimplement set_source_rgb(gc::GraphicsContext, ::Real, ::Real, ::Real)
@mustimplement set_source_rgba(gc::GraphicsContext, ::Real, ::Real, ::Real, ::Real)
@mustimplement set_source(gc::GraphicsContext, src)

# set source color as a Color
function set_source(gc::GraphicsContext, c::Color)
    rgb = convert(RGB, c)
    set_source_rgb(gc, rgb.r, rgb.g, rgb.b)
end

@mustimplement clip(gc::GraphicsContext)
@mustimplement clip_preserve(gc::GraphicsContext)
@mustimplement reset_clip(gc::GraphicsContext)

@mustimplement move_to(gc::GraphicsContext, ::Real, ::Real)
@mustimplement line_to(gc::GraphicsContext, ::Real, ::Real)
@mustimplement rel_move_to(gc::GraphicsContext, ::Real, ::Real)
@mustimplement rel_line_to(gc::GraphicsContext, ::Real, ::Real)
@mustimplement arc(gc::GraphicsContext, ::Real, ::Real, ::Real, ::Real, ::Real)
@mustimplement close_path(gc::GraphicsContext)
@mustimplement new_path(gc::GraphicsContext)
@mustimplement new_sub_path(gc::GraphicsContext)

@mustimplement fill(gc::GraphicsContext)
@mustimplement fill_preserve(gc::GraphicsContext)
@mustimplement stroke(gc::GraphicsContext)
@mustimplement stroke_preserve(gc::GraphicsContext)
@mustimplement paint(gc::GraphicsContext)
stroke_transformed(gc::GraphicsContext) = stroke(gc)
stroke_transformed_preserve(gc::GraphicsContext) = stroke_preserve(gc)

# generic path functions

function rectangle(gc::GraphicsContext, x::Real, y::Real, width::Real, height::Real)
    move_to(gc, x, y)
    rel_line_to(gc, width, 0)
    rel_line_to(gc, 0, height)
    rel_line_to(gc, -width, 0)
    close_path(gc)
end
rectangle(gc::GraphicsContext, user::BoundingBox) = rectangle(gc, user.xmin, user.ymin, width(user), height(user))

circle(ctx::GraphicsContext, x::Real, y::Real, r::Real) =
    arc(ctx, x, y, r, 0., 2pi)

function polygon(gc::GraphicsContext, verts::Matrix, idx::Vector)
    move_to(gc, verts[1,idx[1]], verts[2,idx[1]])
    for i=2:length(idx)
        n = idx[i]
        line_to(gc, verts[1,n], verts[2,n])
    end
    close_path(gc)
end

function polygon(self::GraphicsContext, points::AbstractVector)
    move_to(self, points[1].x, points[1].y)
    for i in 2:length(points)
        line_to(self, points[i].x, points[i].y)
    end
    close_path(self)
end

end # module


text:
#
module Spectral

using LinearAlgebra
using FastGaussQuadrature

abstract type AbstractSpectralSpace end

end # module

ext: jl


text:
# Continuation of ex 7.13, Least Squares WRM

using NumericalMethodsforEngineers
using Test

@sym begin
  ClearAll(xi, yi, N1, Y, Ydotdot, C11, ytilde1, ytilde2)
  xi = [0, 1//3, 2//3, 1]
  yi = [0, a, b, 1]
  Y(x_) := LagrangePolynomial(xi, yi)
  #
  # Can be re-formulated as Y(x) = F(x) + C1(a) * Ψ1(x) + C2(a) * Ψ2(x)
  #
  # F(x_) := (1/2)*(9x^3 - 9x^2 + 2x)
  # C1 = (27/2)*(a-b)
  Ψ1(x_) := x^3 - x^2
  # C2 = -(9/2)*(2a-b)
  Ψ2(x_) := x^2 - x
  #
  # Y(x_) := F(x) + C1(a)*Ψ1(x) + C2(a)*Ψ2(x)
  #
  Ydotdot(x_) = D(Y(x), x, 2)
  R(x_) := Ydotdot(x) - 3*x - 4*Y(x)
  Rdota(x_) = D(R(x), a)
  Rdotb(x_) = D(R(x), b)
  sol = Solve([Integrate(R(x)*Rdota(x), [x, 0, 1]), Integrate(R(x)*Rdotb(x), [x, 0, 1])], [a, b])
  SetJ(r, ToString(Simplify(R(x))))
  SetJ(a, ToString(sol[1][1][2]))
  SetJ(b, ToString(sol[1][2][2]))
  ytilde1 = Simplify(Y(x) ./ (a => sol[1][1][2]))
  ytilde2 = Simplify(ytilde1 ./ (b => sol[1][2][2]))
  SetJ(y, ToString(Simplify(Expand(ytilde2))));
end

println("\n\nExample 7.17b: y'' = 3x + 4y, y(0)=0, y(1)=1")
println("by 2-point least squares Weighted Residual Method")
@sym Println("\nY(x) = $(Simplify(Expand(Y(x))))\n")
@sym Println("R(x) = $(R(x))\n")
@eval a = $(Meta.parse(a))
@eval b = $(Meta.parse(b))
println("(a,b) = $((a, b))\n")
println("C1 = $((27/2*(a-b)))\n")
println("C2 = $(-(9/2*(2a-b)))\n")
@eval rf_2pt_leastsquares(x, a, b) = $(Meta.parse(r))
@eval ytilde_2pt_leastsquares(x) = $(Meta.parse(y))
println()

rf_2pt_leastsquares_1(x) = rf_2pt_leastsquares(x, a, b)
println()

@test r == "-9 - 45a + 36b + 20x - 63b*x + 45a*x + 18*x^2 - 72b*x^2 + 90a*x^2 - 18x^3 - 54a*x^3 + 54b*x^3"
@test y == "(1/9693)*x*(2420 - 3311x + 10584*x^2)"
#@test (quadgk(rf_2pt_leastsquares_1, 0, 1))[1] < 5*eps()


text:
module ReinterpretTests

using Test
using Gridap.TensorValues

v = VectorValue(1.3,2.1)
V = fill(v,5)
A = reinterpret(V)
R = [1.3 1.3 1.3 1.3 1.3; 2.1 2.1 2.1 2.1 2.1]
@test A == R

v = TensorValue(2.0,3.5,1.0,4.0)
V = fill(v,5)
A = reinterpret(V)
R = [2.0 2.0 2.0 2.0 2.0; 3.5 3.5 3.5 3.5 3.5; 1.0 1.0 1.0 1.0 1.0; 4.0 4.0 4.0 4.0 4.0]
@test A == R

v = SymFourthOrderTensorValue(1,2,3, 4,5,6, 7,8,9)
V = fill(v,5)
A = reinterpret(V)
R = [1 1 1 1 1; 2 2 2 2 2; 3 3 3 3 3; 4 4 4 4 4; 5 5 5 5 5; 6 6 6 6 6; 7 7 7 7 7; 8 8 8 8 8; 9 9 9 9 9]
@test A == R

v = SymTensorValue(1,2,3)
V = fill(v,5)
A = reinterpret(V)
R = [1 1 1 1 1; 2 2 2 2 2; 3 3 3 3 3]
@test A == R

end # module ReinterpretTests


text:
using Test
using StatsBase, Statistics, Random

using CommonLH.StatsLH

## Weighted density
## should be tested more substantively
function density_weighted_test()
	n = 1000;
	xV = collect(range(5.0, 2.0, length = n));
	wtV = collect(range(0.5, 0.8, length = n));
	densityV, midV = StatsLH.density_weighted(xV, wtV, true);
	@testset "Density weighted" begin
		@test length(xV) == length(densityV) + 1;
	end
end

## Weighted std deviation
function std_w_test()
	n = Int64(1e4);
	rng = MersenneTwister(1234);
	z = randn(rng, Float64, 1, n);
	wtV = ones(1, n);
	zStd, zMean = std_w(z, wtV)
	println("zMean = $zMean,  zStd = $zStd")
	@testset "std_w" begin
		@test zMean ≈ 0.0 atol=0.01
		@test zStd ≈ 1.0 atol = 1e-3

		zStd2 = std(z)
		@test zStd2 ≈ zStd atol = 1e-4
	end
end

@testset "Statistics" begin
	include("statistics/std_std_test.jl")
	density_weighted_test();
	std_w_test();
	discretize_test();
	bin_edges_test();
	bin_edges_weighted_test();
	discretize_given_percentiles_test();
	discretize_given_pct_weighted_test();
end

#@testset "StatsLH" begin
#    std_w_test()
#    std_std_test()
#    discretize_test()
#    discretize_given_pct_weighted_test();
#    bin_edges_test();
#    bin_edges_weighted_test();
#    discretize_given_percentiles_test()
#    density_weighted_test()
#    # return true
#end

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


text:
export Ban

"""
A [`User`](@ref) ban.
More details [here](https://discordapp.com/developers/docs/resources/guild#ban-object).
"""
struct Ban
    reason::Nullable{String}
    user::User
end
@boilerplate Ban :constructors :docs :lower :merge :mock


text:
"""
# module DumpInput

# Examples

```jldoctest
julia>
```
"""
module DumpInput

using StaticArrays: SMatrix
using JuliennedArrays
using FilePaths

export dumps,
    dump

dumps(v::NamelistVariable) = string(v.name, " = ", fstring(v.value))

dumps(n::Namelist) = "\&$(uppercase(n.name))\n" * join([dumps(v) * ",\n" for v in n.values]) * "/\n"

dumps(c::Card) = "$(uppercase(c.name))\n" * dumps(c.values)

dumps(c::OptionedCard) = "$(uppercase(c.name)) \{ $(c.option) \}\n" * dumps(c.values)

dumps(input::NamelistsInput) = join([dumps(namelist) for namelist in input.namelists])

dumps(input::NamelistsCardsInput) = (join ∘ append!)([dumps(namelist) for namelist in input.namelists], [dumps(card) for card in input.cards])

dump(io::IO, input::AbstractInput) = write(io, dumps(input))

dump(path::AbstractPath, input::AbstractInput) = open(path, "w") do f
    dump(f, input)
end

dumps(a::AtomicSpecies) = " $(a.symbol)  $(a.mass)  $(a.pseudopotential)"

function dumps(c::CellParameters, delim::String=" ")
    m = SMatrix{3, 4, String}(hcat((x -> string(x, delim)).(c.coordinates), ["\n", "\n", "\n"]))
    join(map(Reduce(*), julienne(m, (*, :))))
end

end


text:
function nonrisk_clustering(dist::Sampleable{Multivariate, Continuous}, Ω::AbstractRiskRegion,
                            num_risk::Int64, num_non_risk::Int64)
    max_non_risk = (num_risk + num_non_risk) * 10
    dim = length(dist)
    scenarios = Array{Float64}(undef, dim, num_risk + num_non_risk)
    non_risk_scen = Array{Float64}(undef, dim, max_non_risk)
    r = 0                  # Counter for risk scenarios
    nr = 0   # Counter for non-risk scenarios

    while r < num_risk
        scen = rand(dist)
        if scen ∈ Ω
            scenarios[:,r+1] = scen
            r += 1
        else
            non_risk_scen[:, nr+1] = scen
            nr += 1
        end
    end

    if nr <= num_non_risk
        scenarios[:,r+1:r+nr] = non_risk_scen[:,1:nr]
        probs = fill(1.0/(r+nr), r+nr)
        return scenarios, probs
    else
        num_clusters = num_non_risk
        cluster_results = kmeans(non_risk_scen[:,1:nr], num_clusters)
        scenarios[:, (num_risk+1):(num_risk+num_non_risk)] = cluster_results.centers
        probs = Array{Float64}(undef, num_risk + num_clusters)
        probs[1:r] .= 1.0/(r + nr)
        probs[r+1:r + num_clusters] .= cluster_results.cweights/(r + nr)
    end
    return scenarios, probs
end

function nonrisk_clustering(scenarios::Matrix{Float64}, Ω::AbstractRiskRegion, num_non_risk::Int64)
    dim, num_scen = size(scenarios)
    new_scenarios = Array{Float64}(undef, dim, num_scen)
    non_risk_scenarios = Array{Float64}(undef, dim, num_scen)
    r = 0   # Risk scenario counter
    nr = 0  # Non-risk scenario counter
    for s in 1:num_scen
        scen = scenarios[:,s]
        if scen ∈ Ω
            new_scenarios[:,r+1] = scen
            r += 1
        else
            non_risk_scenarios[:, nr+1] = scen
            nr += 1
        end
    end
    non_risk_scenarios = non_risk_scenarios[:, 1:nr]
    if nr > num_non_risk
        cluster_results = kmeans(non_risk_scenarios, num_non_risk)
        new_scenarios[:,(r+1):(r+num_non_risk)] = cluster_results.centers
        # The following line causes a garbage-collection related bug, so we use something less efficient...
        # new_scenarios = pointer_to_array(pointer(new_scenarios), (dim, r + num_non_risk))
        new_scenarios = new_scenarios[:, 1:r + num_non_risk]
        new_probs = Array{Float64}(undef, r + num_non_risk)
        new_probs[1:r] = fill(1.0/(r + nr), r)
        new_probs[r+1:r + num_non_risk] = cluster_results.cweights/(r + nr)
    else
        new_scenarios[:,(r+1):(r+nr)] = non_risk_scenarios
        # new_scenarios = pointer_to_array(pointer(new_scenarios), (dim, r + nr)) # Bug
        new_scenarios = new_scenarios[:, 1:(r+nr)]
        new_probs = fill(1.0/(r+nr))
    end
    return new_scenarios, new_probs
end


text:
baremodule PublicAPI

export @public

macro public end
macro strict end

function of end

module Internal

using ..PublicAPI: PublicAPI
import ..PublicAPI: @public, @strict

include("registry.jl")
include("query.jl")
include("imports.jl")

# Import README as a docstring for testing it via doctest.
# Not calling `@eval` here to avoid "eval from module PublicAPI to Internal".
function define_docstring()
    path = joinpath(@__DIR__, "..", "README.md")
    include_dependency(path)
    doc = read(path, String)
    doc = replace(doc, r"^```julia"m => "```jldoctest README")
    @eval PublicAPI $Base.@doc $doc PublicAPI
end

end  # module Internal

Internal.define_docstring()

@public @public
@public @strict
@public of

end  # baremodule PublicAPI


text:
function gausshermite( n::Int )
# GAUSSHERMITE(n) COMPUTE THE GAUSS-HERMITE NODES AND WEIGHTS IN O(n) time.

if n < 0
    x = (Float64[],Float64[])
    return x
elseif n == 0
    x = (Float64[],Float64[])
    return x
elseif n == 1
    x = ([0.0],[sqrt(pi)])
    return x
elseif n <= 20
   # GW algorithm
   x = hermpts_gw( n )
elseif n <= 200
   # REC algorithm
   x = hermpts_rec( n )
else
   # ASY algorithm
   x = hermpts_asy( n )
end

if mod(n,2) == 1                              # fold out
    w = [flipdim(x[2][:],1); x[2][2:end]]
    w = (sqrt(pi)/sum(w))*w
    x = ([-flipdim(x[1],1) ; x[1][2:end]], w)
else
    w = [flipdim(x[2][:],1); x[2][:]]
    w = (sqrt(pi)/sum(w))*w
    x = ([-flipdim(x[1],1) ; x[1]], w)
end

end

function hermpts_asy( n::Int )
# Compute Hermite nodes and weights using asymptotic formula

x0 = HermiteInitialGuesses( n ) # get initial guesses
t0 = x0./sqrt(2n+1)
theta0 = acos(t0)               # convert to theta-variable
val = x0;
for k = 1:20
    val = hermpoly_asy_airy(n, theta0);
    dt = -val[1]./(sqrt(2)*sqrt(2n+1)*val[2].*sin(theta0))
    theta0 = theta0 - dt;                        # Newton update
    if norm(dt,Inf) < sqrt(eps(Float64))/10
       break
    end
end
t0 = cos(theta0)
x = sqrt(2n+1)*t0                          #back to x-variable
ders = x.*val[1] + sqrt(2)*val[2]
w = (exp(-x.^2)./ders.^2)';            # quadrature weights

x = (x, w)
end

function hermpts_rec( n::Int )
# Compute Hermite nodes and weights using recurrence relation.

x0 = HermiteInitialGuesses( n )
x0 = x0.*sqrt(2)
val = x0
for kk = 1:10
    val = hermpoly_rec(n, x0)
    dx = val[1]./val[2]
    dx[ isnan( dx ) ] = 0
    x0 = x0 - dx
    if norm(dx, Inf)<sqrt(eps(Float64))
        break
    end
end
x = x0/sqrt(2)
w = exp(-x.^2)./val[2].^2           # quadrature weights

x = (x, w)
end

function hermpoly_rec( n::Int, x0)
# HERMPOLY_rec evaluation of scaled Hermite poly using recurrence

# evaluate:
Hold = exp(-x0.^2/4)
H = x0.*exp(-x0.^2/4)
for k = 1:n-1
    Hold, H = H, (x0.*H./sqrt(k+1) - Hold./sqrt(1+1/k))
end
# return (value, derivative):
val = (H, (-x0.*H + sqrt(n)*Hold))
end

function hermpoly_asy_airy(n::Int, theta)
# HERMPOLY_ASY evaluation hermite poly using Airy asymptotic formula in
# theta-space.

musq = 2n+1;
cosT = cos(theta)
sinT = sin(theta)
sin2T = 2*cosT.*sinT
eta = .5*theta - .25*sin2T
chi = -(3*eta/2).^(2/3)
phi = (-chi./sinT.^2).^(1/4)
C = 2*sqrt(pi)*musq^(1/6)*phi
Airy0 = real(airy(musq.^(2/3)*chi))
Airy1 = real(airy(1,musq.^(2/3)*chi))

# Terms in (12.10.43):
a0 = 1; b0 = 1
a1 = 15/144; b1 = -7/5*a1
a2 = 5*7*9*11/2/144^2; b2 = -13/11*a2
a3 = 7*9*11*13*15*17/6/144^3
b3 = -19/17*a3

# u polynomials in (12.10.9)
u0 = 1; u1 = (cosT.^3-6*cosT)/24
u2 = (-9*cosT.^4 + 249*cosT.^2 + 145)/1152
u3 = (-4042*cosT.^9+18189*cosT.^7-28287*cosT.^5-151995*cosT.^3-259290*cosT)/414720

#first term
A0 = 1
val = A0*Airy0

#second term
B0 = -(a0*phi.^6.*u1+a1*u0)./chi.^2
val = val + B0.*Airy1./musq.^(4/3)

# third term
A1 = (b0*phi.^12.*u2 + b1*phi.^6.*u1 + b2*u0)./chi.^3
val = val + A1.*Airy0/musq.^2

# fourth term
B1 = -(phi.^18.*u3 + a1*phi.^12.*u2 + a2*phi.^6.*u1 + a3*u0)./chi.^5
val = val + B1.*Airy1./musq.^(4/3+2)

val = C.*val

## Derivative

eta = .5*theta - .25*sin2T
chi = -(3*eta/2).^(2/3)
phi = (-chi./sinT.^2).^(1/4)
C = sqrt(2*pi)*musq^(1/3)./phi

# v polynomials in (12.10.10)
v0 = 1; v1 = (cosT.^3+6*cosT)/24
v2 = (15*cosT.^4-327*cosT.^2-143)/1152
v3 = (259290*cosT + 238425*cosT.^3 - 36387*cosT.^5 + 18189*cosT.^7 - 4042*cosT.^9)/414720

# first term
C0 = -(b0*phi.^6.*v1 + b1.*v0)./chi
dval = C0.*Airy0/musq.^(2/3)

# second term
D0 =  a0*v0
dval = dval + D0*Airy1

# third term
C1 = -(phi.^18.*v3 + b1*phi.^12.*v2 + b2*phi.^6.*v1 + b3*v0)./chi.^4
dval = dval + C1.*Airy0/musq.^(2/3+2)

#fourth term
D1 = (a0*phi.^12.*v2 + a1*phi.^6.*v1 + a2*v0)./chi.^3
dval = dval + D1.*Airy1/musq.^2

dval = C.*dval

val = (val, dval)
end

function HermiteInitialGuesses( n::Int )
#HERMITEINTITIALGUESSES(N), Initial guesses for Hermite zeros.
#
# [1] L. Gatteschi, Asymptotics and bounds for the zeros of Laguerre
# polynomials: a survey, J. Comput. Appl. Math., 144 (2002), pp. 7-27.
#
# [2] F. G. Tricomi, Sugli zeri delle funzioni di cui si conosce una
# rappresentazione asintotica, Ann. Mat. Pura Appl. 26 (1947), pp. 283-300.

# Error if n < 20 because initial guesses are based on asymptotic expansions:
@assert n>=20

# Gatteschi formula involving airy roots [1].
# These initial guess are good near x = sqrt(n+1/2);
if mod(n,2) == 1
    m = (n-1)>>1
    bess = collect(1:m)*pi
    a = .5
else
    m = n>>1
    bess = (collect(0:m-1)+.5)*pi
    a = -.5
end
nu = 4*m + 2*a + 2
T(t) = t.^(2/3).*(1+5/48*t.^(-2)-5/36*t.^(-4)+(77125/82944)*t.^(-6) -108056875/6967296*t.^(-8)+162375596875/334430208*t.^(-10));
airyrts = -T(3/8*pi*(4*collect(1:m)-1))

airyrts_exact = [-2.338107410459762           # Exact Airy roots.
    -4.087949444130970
    -5.520559828095555
    -6.786708090071765
    -7.944133587120863
    -9.022650853340979
    -10.040174341558084
    -11.008524303733260
    -11.936015563236262
    -12.828776752865757]
airyrts[1:10] = airyrts_exact  # correct first 10.

x_init = sqrt(abs(nu + 2^(2/3)*airyrts*nu^(1/3) + 1/5*2^(4/3)*airyrts.^2*nu^(-1/3) +
    (11/35-a^2-12/175*airyrts.^3)/nu + (16/1575*airyrts+92/7875*airyrts.^4)*2^(2/3)*nu^(-5/3) -
    (15152/3031875*airyrts.^5+1088/121275*airyrts.^2)*2^(1/3)*nu^(-7/3)))
x_init_airy = real( flipdim(x_init,1) )

# Tricomi initial guesses. Equation (2.1) in [1]. Originally in [2].
# These initial guesses are good near x = 0 . Note: zeros of besselj(+/-.5,x)
# are integer and half-integer multiples of pi.
# x_init_bess =  bess/sqrt(nu).*sqrt((1+ (bess.^2+2*(a^2-1))/3/nu^2) );
Tnk0 = pi/2*ones(m)
nu = (4*m+2*a+2)
rhs = (4*m-4*collect(1:m)+3)./nu*pi

for k = 1:7
    val = Tnk0 - sin(Tnk0) - rhs
    dval = 1 - cos(Tnk0)
    dTnk0 = val./dval
    Tnk0 = Tnk0 - dTnk0
end

tnk = cos(Tnk0/2).^2
x_init_sin = sqrt(nu*tnk - (5./(4*(1-tnk).^2) - 1./(1-tnk)-1+3*a^2)/3/nu)

# Patch together
p = 0.4985+eps(Float64)
x_init = [x_init_sin[1:convert(Int,floor(p*n))] ;
x_init_airy[convert(Int,ceil(p*n)):end]]

if mod(n, 2) == 1
    x_init = [0 ; x_init]
    x_init = x_init[1:m+1]
else
    x_init = x_init[1:m]
end

return x_init
end

function hermpts_gw( n::Int )
# Golub--Welsch algorithm. Used here for n<=20.

    beta = sqrt(.5*(1:n-1))              # 3-term recurrence coeffs
    T = diagm(beta, 1) + diagm(beta, -1)   # Jacobi matrix
    (D, V) = eig(T)                      # Eigenvalue decomposition
    indx = sortperm(D)                  # Hermite points
    x = D[indx]
    w = sqrt(pi)*V[1,indx].^2            # weights

    # Enforce symmetry:
    ii = floor(Int, n/2)+1:n
    x = x[ii]
    w = w[ii]
    return (x,w)
end


text:
"""
    expval1(basis::SingleBlockBasis, op::AbstractMatrix, i::Int, wf::AbstractVector{Float64})

Calculate the expectation value of the one-body operator `op` at site `i` in
the state `wf`.
"""
function expval1(basis::SingleBlockBasis{N}, op::AbstractMatrix, i::Int, wf::AbstractVector{Float64}) where {N}
    w = Array{Int}(undef, 2N)

    result = 0.0

    for col in 1:basis.block.size
        v = basis.block.vectors[col]

        li = v[2i-1]
        mi = v[2i]
        coli = li*(li+1) + mi + 1

        rowi = 0
        for lip in 0:basis.l_max
            for mip in (-lip):lip
                rowi += 1
                op[rowi, coli] == 0.0 && continue

                w .= v
                w[2i-1] = lip
                w[2i] = mip

                haskey(basis.block.lookup, w) || continue

                result += wf[basis.block.lookup[w]] * wf[col] * op[rowi, coli]
            end
        end
    end

    result
end

"""
    expval2(basis::SingleBlockBasis, op_i::AbstractMatrix, i::Int, op_j::AbstractMatrix, j::Int, wf::AbstractVector{Float64})

Calculate the expectation value of the two-body operator composed of `op_i` at
site `i` and `op_j` at site `j` in the state `wf`.
"""
function expval2(basis::SingleBlockBasis{N}, op_i::AbstractMatrix, i::Int, op_j::AbstractMatrix, j::Int, wf::AbstractVector{Float64}) where {N}
    w = Array{Int}(undef, 2N)

    result = 0.0

    for col in 1:basis.block.size
        v = basis.block.vectors[col]

        li = v[2i-1]
        mi = v[2i]
        coli = li*(li+1) + mi + 1
        lj = v[2j-1]
        mj = v[2j]
        colj = lj*(lj+1) + mj + 1

        rowi = 0
        for lip in 0:basis.l_max
            for mip in (-lip):lip
                rowi += 1
                op_i[rowi, coli] == 0.0 && continue

                rowj = 0
                for ljp in 0:basis.l_max
                    for mjp in (-ljp):ljp
                        rowj += 1
                        op_j[rowj, colj] == 0.0 && continue

                        w .= v
                        w[2i-1] = lip
                        w[2i] = mip
                        w[2j-1] = ljp
                        w[2j] = mjp

                        haskey(basis.block.lookup, w) || continue

                        result += wf[basis.block.lookup[w]] * wf[col] * op_i[rowi, coli] * op_j[rowj, colj]
                    end
                end
            end
        end
    end

    result
end


text:
# Conversion & Promotion
# ======================
#
# Conversion methods for biological sequences.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

# Promotion
# ---------
for alph in (DNAAlphabet, RNAAlphabet)
    @eval function Base.promote_rule(::Type{BioSequence{A}}, ::Type{BioSequence{B}}) where {A<:$alph,B<:$alph}
        return BioSequence{promote_rule(A, B)}
    end
end

# Conversion
# ----------

# Create a 4 bit DNA/RNA sequence from a 2 bit DNA/RNA sequence, and vice-versa.
for (alpha, alphb) in [(DNAAlphabet{4}, DNAAlphabet{2}), # DNA to DNA
                       (DNAAlphabet{2}, DNAAlphabet{4}),
                       (RNAAlphabet{4}, RNAAlphabet{2}), # RNA to RNA
                       (RNAAlphabet{2}, RNAAlphabet{4}),
                       (DNAAlphabet{2}, RNAAlphabet{4}), # DNA to RNA
                       (DNAAlphabet{4}, RNAAlphabet{2}),
                       (DNAAlphabet{2}, RNAAlphabet{2}),
                       (DNAAlphabet{4}, RNAAlphabet{4}),
                       (RNAAlphabet{4}, DNAAlphabet{2}), # RNA to DNA
                       (RNAAlphabet{2}, DNAAlphabet{4}),
                       (RNAAlphabet{2}, DNAAlphabet{2}),
                       (RNAAlphabet{4}, DNAAlphabet{4})]

    @eval function Base.convert(::Type{BioSequence{$alpha}}, seq::BioSequence{$alphb})
        return BioSequence{$alpha}(seq)
    end
end

# Convert from a BioSequence to to a DNA or RNA vector
function Base.convert(::Type{BioSequence{A}}, seq::Vector) where A<:Alphabet
    return BioSequence{A}(seq)
end
Base.convert(::Type{Vector}, seq::BioSequence) = collect(seq)
function Base.convert(::Type{Vector{DNA}}, seq::BioSequence{<:DNAAlphabet})
    return collect(seq)
end
function Base.convert(::Type{Vector{RNA}}, seq::BioSequence{<:RNAAlphabet})
    return collect(seq)
end
Base.convert(::Type{Vector{AminoAcid}}, seq::AminoAcidSequence) = collect(seq)

# Covert from a string to a BioSequence and _vice versa_.
function Base.convert(::Type{S}, seq::BioSequence) where {S<:AbstractString}
    return S([Char(x) for x in seq])
end
Base.String(seq::BioSequence) = convert(String, seq)
Base.convert(::Type{BioSequence{A}}, seq::AbstractString) where A = BioSequence{A}(seq)


text:
export get_hash_types

# This is a mirror of the nettle-meta.h:nettle_hash struct
struct NettleHash
    name::Ptr{UInt8}
    context_size::Cuint
    digest_size::Cuint
    block_size::Cuint
    init::Ptr{Cvoid}     # nettle_hash_init_func
    update::Ptr{Cvoid}   # nettle_hash_update_func
    digest::Ptr{Cvoid}   # nettle_hash_digest_func
end

# We convert from the above NettleHash to a HashType for two reasons:
#   First, we need a way to keep the actual pointer address around
#   Secondly, it's nice to convert the name from a pointer to an actual string
# This cuts down on the amount of work we have to do for some operations
# (especially HMAC operations) at the cost of a bit of memory.  I'll take it.
struct HashType
    name::String
    context_size::Cuint
    digest_size::Cuint
    block_size::Cuint
    init::Ptr{Cvoid}     # nettle_hash_init_func
    update::Ptr{Cvoid}   # nettle_hash_update_func
    digest::Ptr{Cvoid}   # nettle_hash_digest_func

    # This pointer member not actually in the original nettle struct
    ptr::Ptr{Cvoid}
end

# The function that maps from a NettleHash to a HashType
function HashType(nh::NettleHash, nhptr::Ptr{Cvoid})
    HashType( uppercase(unsafe_string(nh.name)),
                    nh.context_size, nh.digest_size, nh.block_size,
                    nh.init, nh.update, nh.digest, nhptr)
end

# The global dictionary of hash types we know how to construct
const _hash_types = Dict{AbstractString,HashType}()

# We're going to load in each NettleHash struct individually, deriving
# HashAlgorithm types off of the names we find, and storing the output
# and context size from the data members in the C structures
function get_hash_types()
    # If we have already gotten the hash types from libnettle, don't query again
    if isempty(_hash_types)
        hash_idx = 1
        # nettle_hashes is an array of pointers ended by a NULL pointer, continue reading hash types until we hit it
        while( true )
            nhptr = unsafe_load(cglobal(("nettle_hashes",libnettle),Ptr{Ptr{Cvoid}}),hash_idx)
            if nhptr == C_NULL
                break
            end
            hash_idx += 1

            nh = unsafe_load(convert(Ptr{NettleHash}, nhptr))
            hash_type = HashType(nh, convert(Ptr{Cvoid},nhptr))
            _hash_types[hash_type.name] = hash_type
        end
    end
    return _hash_types
end


text:
using SuperTB, PyPlot, DelimitedFiles
close("all")

# set up a model (amorphous system with AII-class hopping)
hop = (ri, rj, per) -> hopAgarwala(ri, rj, per)
tb = TBModel(hop)
tb.pos         = rand(24^2,2) #rand(576,2)
tb.periodicity = trues(2) # periodic boundary conditions
tb.basisdim    = 2

# solve the system
ε, eigvec = solve(tb; eigvecflag=true)

# plot the eigenenergies
figure()
plot(ε, ".")
# compare with manually read-off data from Agarwala PRL paper (Fig. 2a)
fig2adata = readdlm(join(pkgdir(SuperTB)*"/data/fig2a_agarwala.txt"), ',', Float64)
plot(fig2adata[:,1], fig2adata[:,2], ".r")
ylim(-4.5, 4.5)

# plot the amorphous lattice (and its neighbors, if for-loops are uncommted)
plot(tb)

# calculate the density of states (dos) and plot it
dos(eeig, espan) = map( (x)-> sum(imag(inv.((x - 0.01im) .- eeig))), espan)
εˢ = range(minimum(ε), maximum(ε), length=250)
dosˢ = dos(ε, εˢ)

figure()
plot(εˢ, dosˢ, "-k")

# compute Bott indices
bottatindices = size(tb.pos,1) .+ (-150:150) # focus on the middle part of the spectrum...
β = bott(eigvec, tb, bottatindices)
figure()
# TODO: Does this actually work? Seems like we broke it? 8 June, 2020
plot(ε[bottatindices], β, ".-")

max_issues_repo_name: thchr/super-tb

max_forks_repo_name: thchr/super-tb


text:
"""

```
mc(sys::System{T},opt=Options=Options()) where T
```

Peforms a Monte-Carlos simulation.

"""
function mc(sys::System{T},opt::Options=Options()) where T

  println("""
  -----------------------------------------------------------------------------------------
  Monte-Carlo simulation
  -----------------------------------------------------------------------------------------
  """)

  # aliases to simplify code
  @unpack x0 = sys
  @unpack nsteps, kT, alpha = opt
  u(x) = potential(x,sys,opt)

  # Copy the initial point, to preserve it
  x = copy(x0)
  xtrial = copy(x0)

  # Open trajectory file for writting
  trajectory_file = open(opt.trajectory_file,"w")

  ucurrent = u(x)
  println(" Potential energy at initial point: " , ucurrent)

  # log vector will contain potential energy
  out = zeros(1)

  # Write coordinates to trajectory file, and update log vectors
  printxyz(0,x,sys,trajectory_file)
  out[1] = ucurrent

  # Running simulation

  last_print = 0
  last_printxyz = 0
  naccepted = 0
  for istep in 1:nsteps

    # Creating new positions
    @. xtrial = x + alpha*randn(T)

    # Computing the potential energy of this new configuration
    utrial = u(xtrial)

    # Check if this conformation is accepted
    if (utrial < ucurrent) ||
       (exp( -(utrial-ucurrent) / kT) > rand())

      naccepted += 1
      ucurrent = utrial
      @. x = xtrial
      push!(out,ucurrent)

      if istep - last_print >= opt.iprint
        @printf(" STEP= %10i U = %12.5f AR = %12.5f\n",
                istep, ucurrent, naccepted/istep)
        last_print = istep
      end
      if istep - last_printxyz >= opt.iprintxyz
        printxyz(naccepted,x,sys,trajectory_file)
        last_printxyz = istep
      end

    end

  end
  println(" Final acceptance rate = ", naccepted/nsteps)

  close(trajectory_file)
  println(" Wrote trajectory file: ", opt.trajectory_file)
  return out
end


text:
good_initialstate =
Complex{Float64}[0.12985363466416291 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.09182038560852626 - 0.11779111254135884im 0.20809568265079803 + 0.06492681733208146im 0.0 + 0.0im -0.20776581546266068 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.09182038560852626 + 0.11779111254135884im 0.0 + 0.0im 0.0 + 0.0im 0.17177593422562445 + 0.0im 0.08825031202852393 + 0.23467519763362607im 0.0 + 0.0im -0.14691261697364202 - 0.18846578006617412im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.20809568265079803 - 0.06492681733208146im 0.0 + 0.0im 0.0 + 0.0im 0.08825031202852393 - 0.23467519763362607im 0.3659451263699554 + 0.0im 0.0 + 0.0im -0.3329530922412769 + 0.10388290773133034im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 -0.20776581546266068 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.14691261697364202 + 0.18846578006617412im -0.3329530922412769 - 0.10388290773133034im 0.0 + 0.0im 0.3324253047402571 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im]
good_hamiltonian =
Complex{Float64}[-2.1766638633710218e-24 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 2.1766638633710218e-24 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 6.294766119337809e-27 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 6.529991590113067e-24 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 8.902143617934187e-27 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.0883319316855107e-23 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.0902854740456259e-26 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.523664704359715e-23 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.2589532238675618e-26 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.9589974770339197e-23 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.4075524945301896e-26 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 2.3943302497081238e-23 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.5418965042537032e-26 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 2.8296630223823285e-23 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.6654385713082973e-26 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 3.264995795056533e-23 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.7804287235868373e-26 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 3.700328567730737e-23 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.8884298358013428e-26 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 4.135661340404942e-23 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.9905798275166757e-26 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 4.570994113079146e-23 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 2.0877377360885034e-26 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 6.294766119337809e-27 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 2.1766638633710218e-24 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 8.902143617934187e-27 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 6.529991590113065e-24 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.0902854740456259e-26 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.088331931685511e-23 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.2589532238675618e-26 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.523664704359715e-23 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.4075524945301896e-26 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.9589974770339197e-23 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.5418965042537032e-26 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 2.3943302497081244e-23 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.6654385713082973e-26 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 2.829663022382328e-23 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.7804287235868373e-26 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 3.264995795056533e-23 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.8884298358013428e-26 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 3.700328567730738e-23 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 1.9905798275166757e-26 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 4.135661340404942e-23 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 2.0877377360885034e-26 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 4.5709941130791465e-23 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 5.0063268857533505e-23 + 0.0im]
good_time_vec = 0.0:1.6666666666666665e-9:4.983333333333333e-7
good_time_evo_array_10 =
Complex{Float64}[-0.15643446504026298 + 0.9876883405951233im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im -0.09780950046421157 - 0.6175448817059057im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.7708220139742106 + 0.1220862132731213im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.13614606094595783 + 0.26720168949539713im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.8499975783949111 - 0.4330953982011464im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.01413968726204487 - 0.014139687262319484im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.7069653946581456 + 0.7069653946581942im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.19436705129183957 - 0.09903495917643126im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.44305693845700483 - 0.869548201952292im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.41288354107963293 + 0.06539432857306797im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.14211025171633407 + 0.8972488170253525im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.5757802750662213 + 0.09119463661639748im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.12710342286844417 - 0.8024994286553415im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.637980531437569 - 0.32506731704889646im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.31692051565025847 + 0.6219915334261098im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.5799230301987155 + 0.5799230301990679im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.40458531738701853 - 0.40458531738677045im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.407695987898866 - 0.8001484289040002im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.39198867058774545 + 0.19972820357760443im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.14896657393809504 + 0.9405379318421944im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.30152389111727645 - 0.047756692738848514im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.15412512269687714 - 0.9731077268758683im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.1690846239077706 - 0.026780373525144868im 0.0 + 0.0im
                 0.0 + 0.0im -0.7708220139742108 + 0.12208621327312136im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.09780950046423909 - 0.6175448817054787im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.8499975783949112 - 0.43309539820114606im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.1361460609457362 + 0.26720168949515516im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.7069653946581458 + 0.7069653946581937im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.014139687262247902 - 0.014139687261971068im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.4430569384570048 - 0.8695482019522922im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.19436705129203227 - 0.0990349591765151im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.14211025171633396 + 0.8972488170253523im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.41288354107958547 + 0.06539432857307387im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.12710342286844423 - 0.8024994286553415im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.5757802750662158 + 0.09119463661629737im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.31692051565025825 + 0.6219915334261097im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.6379805314374246 - 0.32506731704903513im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.40458531738701825 - 0.4045853173867705im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.5799230301987033 + 0.5799230301990606im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.3919886705877454 + 0.19972820357760443im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.4076959878988859 - 0.8001484289039934im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.30152389111727657 - 0.047756692738848514im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.14896657393807192 + 0.9405379318422152im 0.0 + 0.0im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.16908462390777063 - 0.02678037352514489im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.154125122696876 - 0.973107726875868im 0.0 + 0.0im
                 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -0.45399049974022765 - 0.8910065241878162im]
good_excited_prob_15 = 0.14437365772081995 - 1.734723475976807e-18im
good_photons =
[0.12985363466406216,8.956826296345213e-13,7.247624648163295e-13,0.1717759342254107,0.3659451263677979,2.3265359513421833e-12,0.33242530473921916,9.93579359938434e-13,1.1193469441562706e-12,1.3635633644591424e-12,5.1266325158194194e-14,8.03699601310089e-13]
good_densitymatrix =
Complex{Float64}[0.1298536346651647 - 1.063402748901653e-16im -1.7153270010721216e-14 - 4.7080086215367077e-14im -2.5511702622259665e-13 + 7.973876287596511e-14im 0.09182038560763076 - 0.11779111253945646im 0.20809568264892228 + 0.06492681733149624im -1.158362235351889e-16 - 5.498501201275925e-17im -0.20776581546116304 + 3.361036859530484e-17im 4.953566800215108e-17 + 6.265009337479234e-17im 3.5455427460540837e-17 + 2.5810088937138977e-17im 1.4762382136831412e-16 + 1.6050771250793775e-16im 1.1261023746769581e-18 + 2.9579136115209247e-18im 7.982687469890087e-17 - 1.2790945920844224e-16im
                 -1.7118263337229884e-14 + 4.705674843303964e-14im -2.0352354545838354e-12 + 1.648559917164815e-16im -4.9370724338264355e-14 - 1.3196058194158109e-13im 4.0899909115615276e-14 - 1.3075291889174604e-14im 1.2218157887270876e-13 + 4.3976505370705745e-13im 1.0583794865968257e-12 + 3.2957789631292955e-13im 7.31710584060233e-16 + 1.3155556139910176e-16im 7.244914381687481e-14 + 5.723046759709995e-17im -2.389171122024771e-17 - 6.283865273786088e-18im 6.009266392897577e-17 + 1.6348262995047207e-16im -8.481472206351771e-17 - 2.0047546904841818e-17im -3.4864993962514174e-17 - 2.4688176266123025e-17im
                 -2.5514920669267246e-13 - 7.970153510347655e-14im -4.9368658802276435e-14 + 1.3211442495155999e-13im 3.0612368128103658e-12 - 1.4735383355504342e-16im 3.434719886292726e-13 + 9.129111675454094e-13im 4.58898415235114e-13 - 1.4322907136522458e-13im 4.0589002529434126e-13 + 1.9881327240580967e-13im 5.920334699672492e-13 + 1.8426645940518877e-13im -3.6217396381620984e-16 + 6.003194483150313e-17im 8.256016414892149e-13 + 2.0604496425047978e-17im -4.1155207727329517e-17 - 1.2466132731311331e-17im -9.243042501685913e-17 - 8.450204737344317e-17im -1.8985803388237902e-16 - 2.479165965109378e-16im
                 0.09182038560763067 + 0.11779111253945646im 4.103629386651778e-14 + 1.3028637055914218e-14im 3.436258312840232e-13 - 9.127528977639034e-13im 0.17177593422589837 - 5.649274873295089e-18im 0.08825031202766084 + 0.23467519763133007im 9.26768015962802e-13 - 2.8934872567172806e-13im -0.14691261697227043 - 0.18846578006476883im 5.838248543713399e-13 + 1.8164562089848732e-13im -5.629904401075741e-16 - 1.1826901620755228e-16im -6.693418122664687e-13 - 8.180647511318491e-19im -6.186831321090148e-17 + 2.537646542941238e-17im -1.883186327933527e-17 - 4.2603628942146454e-16im
                 0.20809568264892228 - 0.0649268173314963im 1.2216034774529736e-13 - 4.3977108600891917e-13im 4.58800786616887e-13 + 1.4328643031850113e-13im 0.08825031202766041 - 0.23467519763133016im 0.3659451263674329 + 1.1091706209568025e-17im 5.042200048614163e-13 + 1.3409704009201915e-12im -0.3329530922382139 + 0.10388290773037483im -5.586014781054598e-13 - 7.054303621848114e-13im 3.4365511475385343e-13 + 1.0692923154846483e-13im 3.156774514293581e-16 - 1.5137710439407656e-16im -7.087357621669534e-13 - 4.3924856111280685e-17im -2.3978419241998925e-16 - 6.530950000269174e-17im
                 -2.443230873360892e-16 + 5.2619774273194625e-17im 1.0583229786972691e-12 - 3.29567428007537e-13im 4.060365831825773e-13 - 1.9878018690022158e-13im 9.267299771968195e-13 + 2.8947633464788145e-13im 5.051238568217575e-13 - 1.341249973016051e-12im 1.3964824662438406e-12 - 7.269634753277736e-17im -5.546934912105006e-13 - 1.4759935463199806e-12im -2.6146522652801576e-13 + 8.143105546033151e-14im 2.5133988622170374e-13 + 1.1221177216147583e-12im -1.1175097776044959e-13 - 3.4546250749197675e-14im -5.722833332529649e-16 - 1.221902249601922e-16im 4.134868356213999e-13 - 1.50488548529027e-16im
                 -0.20776581546116268 + 2.231738183402708e-17im 8.637671347709314e-16 - 1.0916593202017007e-16im 5.921774016806236e-13 - 1.8429609240500084e-13im -0.14691261697227043 + 0.18846578006476883im -0.33295309223821407 - 0.1038829077303748im -5.557041827830988e-13 + 1.4763033612257727e-12im 0.33242530473666504 + 1.1534861443919397e-16im 2.1568646451770002e-13 + 5.727926283688743e-13im 5.774304579016301e-13 - 1.799191690481103e-13im -6.461307907547694e-13 - 8.61905138001205e-13im -1.1040179370807719e-12 - 3.435390031559379e-13im -6.338819574842895e-16 - 1.4335769632626904e-16im
                 4.1737699961512194e-17 + 4.161593859558891e-17im 7.25343722055598e-14 + 4.791293838540227e-18im -2.9924200186438947e-16 + 1.3468008688960724e-17im 5.839947110225024e-13 - 1.8166398377969424e-13im -5.585480818370964e-13 + 7.052941868019906e-13im -2.612995166529622e-13 - 8.153711976288852e-14im 2.158971806052743e-13 - 5.728705589444684e-13im 1.2511725298017272e-12 + 4.8327125807884367e-17im -5.80397596075875e-13 - 1.5416792060658365e-12im 3.4379442388972335e-12 - 1.0718480340214739e-12im -1.297848061482234e-13 - 1.4908952608450638e-12im -2.4414211604204425e-13 - 7.548340767884563e-14im
                 2.0827480792435328e-17 + 2.030934958759866e-17im -2.050047953882304e-17 - 1.6256031442946778e-17im 8.255761170691514e-13 - 8.572096773785579e-17im -5.937236734758678e-16 + 1.1404448691793852e-16im 3.4355324570743183e-13 - 1.069460604807609e-13im 2.5145686215362744e-13 - 1.1221249201336459e-12im 5.773345963537267e-13 + 1.7994760273432847e-13im -5.798769379066649e-13 + 1.541781149656775e-12im -8.890167439919724e-13 + 3.951276407509798e-17im 2.930231683153928e-12 + 7.781836377983164e-12im -4.4240392386742905e-13 + 1.380970979107165e-13im 7.787267825068578e-13 + 1.395840966221983e-12im
                 -1.258884020671135e-17 - 2.234519134958331e-16im 3.266096576634333e-17 - 1.1556079303386561e-16im -3.1311535369870906e-17 - 3.768968923856418e-17im -6.693914056672307e-13 - 4.9415774169004345e-17im 3.786368218245795e-16 + 1.2477987167800442e-16im -1.1174283028992002e-13 + 3.4542479190522555e-14im -6.460062108332094e-13 + 8.619036854195803e-13im 3.437688682244186e-12 + 1.0719071185335953e-12im 2.9306133779801788e-12 - 7.782869387688646e-12im -7.768176195559152e-13 + 3.698125280805289e-17im -1.9542452812739042e-12 - 5.189271507488401e-12im -2.0686531118990907e-12 + 6.451322263756805e-13im
                 -5.058718123938594e-17 - 4.335028580376187e-18im -1.1647123398088106e-17 + 3.9659466887952465e-18im -5.3974129171671386e-17 + 6.457194035725409e-17im -2.954381322232078e-17 + 3.247830119932039e-17im -7.085904624360996e-13 + 7.681780670429755e-17im -5.751798283480544e-16 + 6.328595611970431e-17im -1.103898698738213e-12 + 3.435741526345953e-13im -1.2991216120287792e-13 + 1.4908051086184604e-12im -4.4244106321630977e-13 - 1.3809524094327247e-13im -1.9546846477882953e-12 + 5.190176215497891e-12im 1.796982682355512e-12 - 2.8527195347994325e-18im -1.9645512586678855e-12 - 5.221179904755555e-12im
                 5.96412491095817e-17 + 1.4078380886051428e-16im -1.0560874054472164e-16 + 1.6611720717388e-17im -8.915669771545307e-17 + 1.9510939584517166e-16im -9.073770267743712e-18 + 4.0977528209479005e-16im -2.262373695328381e-16 + 3.929262120992453e-17im 4.135600121853594e-13 + 1.4251641231779014e-16im -6.99035587005359e-16 + 9.361023885046045e-17im -2.441371098470797e-13 + 7.550343245870475e-14im 7.786477517647994e-13 - 1.3958304287897092e-12im -2.0684691422657994e-12 - 6.452110705042345e-13im -1.9652449607624674e-12 + 5.221907879517976e-12im 2.652791702468874e-12 - 7.845214038328464e-17im]


text:
module AWSCorePrototype

using XMLDict

include("S3/S3.jl")

export Endpoint

bucket_result = S3.list_buckets()["Buckets"].x
last_bucket = parse_xml(string(bucket_result.lastelement))
bucket_name = Dict(last_bucket)["Name"]

objects_result = S3.list_objects(bucket_name)

for object in Dict(objects_result)["Contents"]
    println(object["Key"])
end

end  # End Module

max_stars_repo_name: mattBrzezinski/AWSCore-Prototype

max_issues_repo_name: mattBrzezinski/AWSCore-Prototype

max_forks_repo_name: mattBrzezinski/AWSCore-Prototype


text:
"""
```
FEMParameters(;
    type::AbstractMesh,
    order::AbstractElementOrder,
    model::AbstractModel,
    dx, dy, dz, mx, my, mz
)
```
Creates [`FEMParameters`](@ref).
"""
function FEMParameters(;
    type::AbstractMesh,
    order::AbstractElementOrder,
    model::AbstractModel,
    dx,
    dy,
    dz,
    mx,
    my,
    mz,
)
    return FEMParameters(type, order, model, dx, dy, dz, mx, my, mz)
end

"""
```
ForceDisplacement(; uTilde, uHat, u, fTilde, fHat, f)
```
Creates [`ForceDisplacement`](@ref).
"""
function ForceDisplacement(; uTilde, uHat, u, fTilde, fHat, f)
    return ForceDisplacement(uTilde, uHat, u, fTilde, fHat, f)
end

"""
```
buildMesh(
    matParams::MaterialParameters,
    femParams::FEMParameters{F1,F2,F3,F4,F5} where {F1<:DispatchRegularCuboidMesh,F2,F3,F4,F5}
)
```
Creates a [`RegularCuboidMesh`](@ref).
"""
function buildMesh(
    matParams::MaterialParameters,
    femParams::FEMParameters{
        F1,
        F2,
        F3,
        F4,
        F5,
    } where {F1 <: DispatchRegularCuboidMesh, F2, F3, F4, F5},
)
    return RegularCuboidMesh(matParams, femParams)
end

"""
```
RegularCuboidMesh(
    matParams::MaterialParameters,
    femParams::FEMParameters{F1,F2,F3,F4,F5} where {F1<:DispatchRegularCuboidMesh,F2<:LinearElement,F3,F4,F5}
)
```

FE domain with linear hexahedral elements with `2 × 2 × 2` Gauss nodes. Uses the canonical node, face and edge ordering of [https://www.mcs.anl.gov/uploads/cels/papers/P1573A.pdf](https://www.mcs.anl.gov/uploads/cels/papers/P1573A.pdf).

```
Corners:

      7----------8
     /|         /|
    / |        / |
   /  |       /  |
  /   4------/---3
 /   /      /   /
5----------6   /
|  /       |  /
| /        | /
|/         |/
1----------2

Edges:

      .----11-----.
     /|          /|
    / 8         / 7
  12  |       10  |
  /   .-----3-/---.
 /   /       /   /
.-----9-----.   /
|  4        |  2
5 /         6 /
|/          |/
.-----1-----.

Faces:

      .----------.
     /|         /|
    / |    3   / |
   /  | 6     /  |
  /   .------/---.
 / 4 /      / 2 /
.----------.   /
|  /     5 |  /
| /  1     | /
|/         |/
.----------.

Z
^   y
|  ^
| /
.-----> x
```

"""
function RegularCuboidMesh(
    matParams::MaterialParameters,
    femParams::FEMParameters{
        F1,
        F2,
        F3,
        F4,
        F5,
    } where {F1 <: DispatchRegularCuboidMesh, F2 <: LinearElement, F3, F4, F5},
)
    μ = matParams.μ
    ν = matParams.ν
    νomνInv = matParams.νopνInv # ν/(1+ν)
    order = femParams.order
    dx = femParams.dx
    dy = femParams.dy
    dz = femParams.dz
    mx = femParams.mx
    my = femParams.my
    mz = femParams.mz

    dxType = typeof(dx)
    mxType = typeof(mx)

    # Element width, height, depth
    w = dx / mx
    h = dy / my
    d = dz / mz

    numElem = mx * my * mz
    mx1 = mx + 1
    my1 = my + 1
    mz1 = mz + 1
    mxm1 = mx - 1
    mym1 = my - 1
    mzm1 = mz - 1
    mx1mz1 = mx1 * mz1
    mymx1mz1 = my * mx1mz1
    mzmx1 = mz * mx1
    numNode = mx1 * my1 * mz1

    # Length scale.
    scale = SVector{3, dxType}(dx, dy, dz) * cbrt(dx * dy * dz)

    # For a regular cuboid mesh this is predefined. Making a volumetric polytope allows us to check if points lie inside the volume simply by doing [x, y, z] ∈ vertices, or checking if they do not belong by doing [x, y, z] ∉ vertices. The symbols are typed as \in and \notin + Tab.
    vtx = SMatrix{3, 8, dxType}(
        0,
        0,
        0,
        dx,
        0,
        0,
        dx,
        dy,
        0,
        0,
        dy,
        0,
        0,
        0,
        dz,
        dx,
        0,
        dz,
        dx,
        dy,
        dz,
        0,
        dy,
        dz,
    )

    vertices = VPolytope(vtx)

    # Faces as defined by the vertices.
    faces = SMatrix{4, 6, mxType}(
        1,
        2,
        6,
        5, # xz plane @ min y
        2,
        3,
        7,
        6, # yz plane @ max x
        3,
        4,
        8,
        7, # xz plane @ max y
        4,
        1,
        5,
        8, # yz plane @ min x
        4,
        3,
        2,
        1, # xy plane @ min z
        5,
        6,
        7,
        8, # xy plane @ max z
    )

    faceMidPt = mean(vtx[:, faces], dims = 2)[:, 1, :]

    # Face normal of the corresponding face.
    faceNorm = SMatrix{3, 6, dxType}(
        0,
        -1,
        0,   # xz plane @ min y
        1,
        0,
        0,    # yz plane @ max x
        0,
        1,
        0,    # xz plane @ max y
        -1,
        0,
        0,   # yz plane @ min x
        0,
        0,
        -1,   # xy plane @ min z
        0,
        0,
        1,    # xy plane @ max z
    )

    coord = zeros(dxType, 3, numNode)           # Node coordinates.
    connectivity = zeros(mxType, 8, numElem)    # Element connectivity.

    cornerNode = SVector{8, Symbol}(
        :x0y0z0,
        :x1y0z0,
        :x1y1z0,
        :x0y1z0,
        :x0y0z1,
        :x1y0z1,
        :x1y1z1,
        :x0y1z1,
    )
    edgeNode = SVector{12, Symbol}(
        :x_y0z0,
        :y_x1z0,
        :x_y1z0,
        :y_x0z0,
        :x_y0z1,
        :y_x1z1,
        :x_y1z1,
        :y_x0z1,
        :z_x0y0,
        :z_x1y0,
        :z_x1y1,
        :z_x0y1,
    )
    faceNode = SVector{6, Symbol}(:xz_y0, :yz_x1, :xz_y1, :yz_x0, :xy_z0, :xy_z1)
    surfNode = (
        # Corners
        x0y0z0 = 1,
        x1y0z0 = mx1,
        x1y1z0 = mx1 * my1,
        x0y1z0 = mx1 * my + 1,
        x0y0z1 = 1 + mx1 * my1 * mz,
        x1y0z1 = mx1 + mx1 * my1 * mz,
        x1y1z1 = mx1 * my1 * mz1,
        x0y1z1 = mx1 * my + 1 + mx1 * my1 * mz,
        # Edges
        x_y0z0 = zeros(mxType, mxm1),
        y_x1z0 = zeros(mxType, mym1),
        x_y1z0 = zeros(mxType, mxm1),
        y_x0z0 = zeros(mxType, mym1),
        x_y0z1 = zeros(mxType, mxm1),
        y_x1z1 = zeros(mxType, mym1),
        x_y1z1 = zeros(mxType, mxm1),
        y_x0z1 = zeros(mxType, mym1),
        z_x0y0 = zeros(mxType, mzm1),
        z_x1y0 = zeros(mxType, mzm1),
        z_x1y1 = zeros(mxType, mzm1),
        z_x0y1 = zeros(mxType, mzm1),
        # Faces
        xz_y0 = zeros(mxType, mxm1 * mzm1),
        yz_x1 = zeros(mxType, mym1 * mzm1),
        xz_y1 = zeros(mxType, mxm1 * mzm1),
        yz_x0 = zeros(mxType, mym1 * mzm1),
        xy_z0 = zeros(mxType, mxm1 * mym1),
        xy_z1 = zeros(mxType, mxm1 * mym1),
    )
    surfNodeArea = (
        # Corners
        x0y0z0 = (h * d + w * d + w * h) / 4,
        x1y0z0 = (h * d + w * d + w * h) / 4,
        x1y1z0 = (h * d + w * d + w * h) / 4,
        x0y1z0 = (h * d + w * d + w * h) / 4,
        x0y0z1 = (h * d + w * d + w * h) / 4,
        x1y0z1 = (h * d + w * d + w * h) / 4,
        x1y1z1 = (h * d + w * d + w * h) / 4,
        x0y1z1 = (h * d + w * d + w * h) / 4,
        # Edges
        x_y0z0 = (h + d) * w / 2,
        y_x1z0 = (w + d) * h / 2,
        x_y1z0 = (h + d) * w / 2,
        y_x0z0 = (w + d) * h / 2,
        x_y0z1 = (h + d) * w / 2,
        y_x1z1 = (w + d) * h / 2,
        x_y1z1 = (h + d) * w / 2,
        y_x0z1 = (w + d) * h / 2,
        z_x0y0 = (h + w) * d / 2,
        z_x1y0 = (h + w) * d / 2,
        z_x1y1 = (h + w) * d / 2,
        z_x0y1 = (h + w) * d / 2,
        # Surfaces
        xz_y0 = w * d,
        yz_x1 = h * d,
        xz_y1 = w * d,
        yz_x0 = h * d,
        xy_z0 = w * h,
        xy_z1 = w * h,
    )

    surfNodeNorm = (
        # Corners
        x0y0z0 = normalize(faceNorm[:, 4] + faceNorm[:, 1] + faceNorm[:, 5]),
        x1y0z0 = normalize(faceNorm[:, 2] + faceNorm[:, 1] + faceNorm[:, 5]),
        x1y1z0 = normalize(faceNorm[:, 2] + faceNorm[:, 3] + faceNorm[:, 5]),
        x0y1z0 = normalize(faceNorm[:, 4] + faceNorm[:, 3] + faceNorm[:, 5]),
        x0y0z1 = normalize(faceNorm[:, 4] + faceNorm[:, 1] + faceNorm[:, 6]),
        x1y0z1 = normalize(faceNorm[:, 2] + faceNorm[:, 1] + faceNorm[:, 6]),
        x1y1z1 = normalize(faceNorm[:, 2] + faceNorm[:, 3] + faceNorm[:, 6]),
        x0y1z1 = normalize(faceNorm[:, 4] + faceNorm[:, 3] + faceNorm[:, 6]),
        # Edges
        x_y0z0 = normalize(faceNorm[:, 1] + faceNorm[:, 5]),
        y_x1z0 = normalize(faceNorm[:, 2] + faceNorm[:, 5]),
        x_y1z0 = normalize(faceNorm[:, 3] + faceNorm[:, 5]),
        y_x0z0 = normalize(faceNorm[:, 4] + faceNorm[:, 5]),
        x_y0z1 = normalize(faceNorm[:, 1] + faceNorm[:, 6]),
        y_x1z1 = normalize(faceNorm[:, 2] + faceNorm[:, 6]),
        x_y1z1 = normalize(faceNorm[:, 3] + faceNorm[:, 6]),
        y_x0z1 = normalize(faceNorm[:, 4] + faceNorm[:, 6]),
        z_x0y0 = normalize(faceNorm[:, 4] + faceNorm[:, 1]),
        z_x1y0 = normalize(faceNorm[:, 2] + faceNorm[:, 1]),
        z_x1y1 = normalize(faceNorm[:, 2] + faceNorm[:, 3]),
        z_x0y1 = normalize(faceNorm[:, 4] + faceNorm[:, 3]),
        # Faces
        xz_y0 = faceNorm[:, 1],
        yz_x1 = faceNorm[:, 2],
        xz_y1 = faceNorm[:, 3],
        yz_x0 = faceNorm[:, 4],
        xy_z0 = faceNorm[:, 5],
        xy_z1 = faceNorm[:, 6],
    )

    nodeEl = 1:8 # Local node numbers.
    dofLocal = Tuple(
        Iterators.flatten((
            3 * (nodeEl .- 1) .+ 1,
            3 * (nodeEl .- 1) .+ 2,
            3 * (nodeEl .- 1) .+ 3,
        )),
    )

    E = 2 * μ * (1 + ν)                 # Young's modulus
    Lame = E * νomνInv / (1 - 2 * ν)    # Lamé's relation
    CDiag = Lame + 2 * μ

    # Stiffness tensor.
    C = SMatrix{6, 6, dxType}(
        CDiag,
        Lame,
        Lame,
        0,
        0,
        0,
        Lame,
        CDiag,
        Lame,
        0,
        0,
        0,
        Lame,
        Lame,
        CDiag,
        0,
        0,
        0,
        0,
        0,
        0,
        μ,
        0,
        0,
        0,
        0,
        0,
        0,
        μ,
        0,
        0,
        0,
        0,
        0,
        0,
        μ,
    )

    # Local nodes Gauss Quadrature. Using 8 nodes per element.
    p = 1 / sqrt(3)

    gaussNodes = SMatrix{3, 8, dxType}(
        -p,
        -p,
        -p, # x0y0z0
        p,
        -p,
        -p,  # x1y0z0
        p,
        p,
        -p,   # x1y1z0
        -p,
        p,
        -p,  # x0y1z0
        -p,
        -p,
        p,  # x0y0z1
        p,
        -p,
        p,   # x1y0z1
        p,
        p,
        p,    # x1y1z1
        -p,
        p,
        p,   # x0y1z1
    )

    # Shape functions and their derivatives.
    N = shapeFunction(
        LinearQuadrangle3D(),
        gaussNodes[1, :],
        gaussNodes[2, :],
        gaussNodes[3, :],
    )
    dNdS = shapeFunctionDeriv(
        LinearQuadrangle3D(),
        gaussNodes[1, :],
        gaussNodes[2, :],
        gaussNodes[3, :],
    )

    # Fill node coordinates.
    @inbounds @simd for k in 1:mz1
        km1 = k - 1
        for j in 1:my1
            jm1 = j - 1
            for i in 1:mx1
                globalNode = i + jm1 * mx1 + km1 * mx1 * my1
                coord[1, globalNode] = (i - 1) * w
                coord[2, globalNode] = jm1 * h
                coord[3, globalNode] = km1 * d
            end
        end
    end

    # Fill element connectivity.
    @inbounds @simd for k in 1:mz
        km1 = k - 1
        for j in 1:my
            jm1 = j - 1
            for i in 1:mx
                globalElem = i + jm1 * mx + km1 * mx * my
                connectivity[1, globalElem] = i + jm1 * mx1 + km1 * mx1 * my1
                connectivity[2, globalElem] = connectivity[1, globalElem] + 1
                connectivity[4, globalElem] = i + j * mx1 + km1 * mx1 * my1
                connectivity[3, globalElem] = connectivity[4, globalElem] + 1
                connectivity[5, globalElem] = i + jm1 * mx1 + k * mx1 * my1
                connectivity[6, globalElem] = connectivity[5, globalElem] + 1
                connectivity[8, globalElem] = i + j * mx1 + k * mx1 * my1
                connectivity[7, globalElem] = connectivity[8, globalElem] + 1
            end
        end
    end

    # Edge node along x
    @inbounds @simd for i in 1:mxm1
        surfNode[:x_y0z0][i] = i + 1                              # x_y0z0
        surfNode[:x_y1z0][i] = i + 1 + mx1 * my                   # x_y1z0
        surfNode[:x_y0z1][i] = i + 1 + mx1 * my1 * mz             # x_y0z1
        surfNode[:x_y1z1][i] = i + 1 + mx1 * my1 * mz + mx1 * my  # x_y1z1
    end
    # Edge node along y
    @inbounds @simd for i in 1:mym1
        surfNode[:y_x0z0][i] = i * mx1 + 1                      # y_x0z0
        surfNode[:y_x1z0][i] = i * mx1 + mx1                    # y_x1z0
        surfNode[:y_x0z1][i] = i * mx1 + 1 + mx1 * my1 * mz     # y_x0z1
        surfNode[:y_x1z1][i] = i * mx1 + mx1 + mx1 * my1 * mz   # y_x1z1
    end
    # Edge node along z
    @inbounds @simd for i in 1:mzm1
        surfNode[:z_x0y0][i] = i * mx1 * my1 + 1                # z_x0y0
        surfNode[:z_x1y0][i] = i * mx1 * my1 + mx1              # z_x1y0
        surfNode[:z_x0y1][i] = i * mx1 * my1 + 1 + mx1 * my     # z_x0y1
        surfNode[:z_x1y1][i] = i * mx1 * my1 + mx1 + mx1 * my   # z_x1y1
    end

    # Face node xy
    @inbounds @simd for j in 1:mym1
        jm1 = j - 1
        for i in 1:mxm1
            surfNode[:xy_z0][i + mxm1 * jm1] = j * mx1 + 1 + i                  # xy_z0
            surfNode[:xy_z1][i + mxm1 * jm1] = j * mx1 + 1 + i + mx1 * my1 * mz # xy_z1
        end
    end

    # Face node xz
    @inbounds @simd for j in 1:mzm1
        jm1 = j - 1
        for i in 1:mxm1
            surfNode[:xz_y0][i + mxm1 * jm1] = j * mx1 * my1 + 1 + i            # xz_y0
            surfNode[:xz_y1][i + mxm1 * jm1] = j * mx1 * my1 + 1 + i + mx1 * my # xz_y1
        end
    end

    # Face node yz
    @inbounds @simd for j in 1:mzm1
        jm1 = j - 1
        for i in 1:mym1
            surfNode[:yz_x0][i + mym1 * jm1] = j * mx1 * my1 + 1 + i * mx1  # yz_x0
            surfNode[:yz_x1][i + mym1 * jm1] = j * mx1 * my1 + mx1 + i * mx1  # yz_x0
        end
    end

    localCoord = SMatrix{3, 8, dxType}(
        coord[1, connectivity[1, 1]],
        coord[2, connectivity[1, 1]],
        coord[3, connectivity[1, 1]],
        coord[1, connectivity[2, 1]],
        coord[2, connectivity[2, 1]],
        coord[3, connectivity[2, 1]],
        coord[1, connectivity[3, 1]],
        coord[2, connectivity[3, 1]],
        coord[3, connectivity[3, 1]],
        coord[1, connectivity[4, 1]],
        coord[2, connectivity[4, 1]],
        coord[3, connectivity[4, 1]],
        coord[1, connectivity[5, 1]],
        coord[2, connectivity[5, 1]],
        coord[3, connectivity[5, 1]],
        coord[1, connectivity[6, 1]],
        coord[2, connectivity[6, 1]],
        coord[3, connectivity[6, 1]],
        coord[1, connectivity[7, 1]],
        coord[2, connectivity[7, 1]],
        coord[3, connectivity[7, 1]],
        coord[1, connectivity[8, 1]],
        coord[2, connectivity[8, 1]],
        coord[3, connectivity[8, 1]],
    )

    localdNdS = SMatrix{8, 3, dxType}(
        dNdS[1][1, 1],
        dNdS[1][1, 2],
        dNdS[1][1, 3],
        dNdS[1][1, 4],
        dNdS[1][1, 5],
        dNdS[1][1, 6],
        dNdS[1][1, 7],
        dNdS[1][1, 8],
        dNdS[1][2, 1],
        dNdS[1][2, 2],
        dNdS[1][2, 3],
        dNdS[1][2, 4],
        dNdS[1][2, 5],
        dNdS[1][2, 6],
        dNdS[1][2, 7],
        dNdS[1][2, 8],
        dNdS[1][3, 1],
        dNdS[1][3, 2],
        dNdS[1][3, 3],
        dNdS[1][3, 4],
        dNdS[1][3, 5],
        dNdS[1][3, 6],
        dNdS[1][3, 7],
        dNdS[1][3, 8],
    )

    J = localCoord * localdNdS
    detJ = det(J)
    invJ = inv(J)'

    B = zeros(dxType, 6, 24, 8)     # Jacobian matrix.
    # Gauss quadrature nodes.
    @inbounds @simd for q in nodeEl
        nx = invJ * dNdS[q]
        # Nodes in element.
        for a in nodeEl
            idx = (a - 1) * 3
            idx1 = idx + 1
            idx2 = idx + 2
            idx3 = idx + 3
            B[1, idx1, q] = nx[1, a]
            # B[1, idx2, q] = 0
            # B[1, idx3, q] = 0

            # B[2, idx1, q] = 0
            B[2, idx2, q] = nx[2, a]
            # B[2, idx3, q] = 0

            # B[3, idx1, q] = 0
            # B[3, idx2, q] = 0
            B[3, idx3, q] = nx[3, a]

            B[4, idx1, q] = nx[2, a]
            B[4, idx2, q] = nx[1, a]
            # B[4, idx3, q] = 0

            B[5, idx1, q] = nx[3, a]
            # B[5, idx2, q] = 0
            B[5, idx3, q] = nx[1, a]

            # B[6, idx1, q] = 0
            B[6, idx2, q] = nx[3, a]
            B[6, idx3, q] = nx[2, a]
        end
    end

    localK = zeros(dxType, 24, 24)  # K for an element.
    @inbounds @simd for q in nodeEl
        localK .+= B[:, :, q]' * C * B[:, :, q]
    end

    localK = localK * detJ
    cntr = 0
    V1 = zeros(dxType, numElem * 24^2)
    V2 = zeros(dxType, numElem * 24^2)
    V3 = zeros(dxType, numElem * 24^2)
    @inbounds @simd for p in 1:numElem
        globalNode = connectivity[nodeEl, p] # Global node numbers
        dofGlobal = Tuple(
            Iterators.flatten((
                3 * (globalNode .- 1) .+ 1,
                3 * (globalNode .- 1) .+ 2,
                3 * (globalNode .- 1) .+ 3,
            )),
        ) # Global degree of freedom
        for i in 1:24
            for j in 1:24
                cntr += 1
                V1[cntr] = dofGlobal[i]
                V2[cntr] = dofGlobal[j]
                V3[cntr] = localK[dofLocal[i], dofLocal[j]]
            end
        end
    end

    numNode3 = 3 * numNode
    globalK = droptol!(sparse(V1, V2, V3, numNode3, numNode3), eps(dxType))
    droptol!(globalK, eps(dxType))
    globalK = sparse(Hermitian(globalK))

    # Surface elements
    mxmy = mx * my
    mxmz = mx * mz
    mymz = my * mz
    surfElemNode = zeros(mxType, 4, 2 * (mxmy + mxmz + mymz))
    cntr = 0

    @inbounds @simd for i in 1:size(faces, 2)
        label = vec(connectivity[faces[:, i], :]')

        if i == 1 || i == 3
            xyz = 2
            dimCoord = @view coord[xyz, label]
            i == 1 ? lim = minimum(dimCoord) : lim = maximum(dimCoord)
        elseif i == 2 || i == 4
            xyz = 1
            dimCoord = @view coord[xyz, label]
            i == 4 ? lim = minimum(dimCoord) : lim = maximum(dimCoord)
        else
            xyz = 3
            dimCoord = @view coord[xyz, label]
            i == 5 ? lim = minimum(dimCoord) : lim = maximum(dimCoord)
        end

        idx = findall(x -> x ≈ lim, dimCoord)
        n = div(length(idx), 4)

        label = label[idx]

        surfElemNode[:, (1:n) .+ cntr] .= reshape(label, :, 4)'
        cntr += n
    end

    return RegularCuboidMesh(
        order,
        dx,
        dy,
        dz,
        mx,
        my,
        mz,
        w,
        h,
        d,
        scale,
        numElem,
        numNode,
        C,
        vertices,
        faces,
        faceNorm,
        faceMidPt,
        cornerNode,
        edgeNode,
        faceNode,
        surfNode,
        surfNodeArea,
        surfNodeNorm,
        surfElemNode,
        coord,
        connectivity,
        globalK,
    )
end

"""
```
BoundaryNode(; index, node)
```
Create [`BoundaryNode`](@ref).
"""
function BoundaryNode(; index, node)
    return BoundaryNode(index, node)
end

"""
```
Boundaries(
    ::FEMParameters{T1,T2,T3,T4,T5} where {T1,T2,T3<:CantileverLoad,T4,T5},
    femMesh::RegularCuboidMesh;
    kw...
)
```
Creates [`Boundaries`](@ref) for loading a hexahedral cantilever.
"""
function Boundaries(
    ::FEMParameters{T1, T2, T3, T4, T5} where {T1, T2, T3 <: CantileverLoad, T4, T5},
    femMesh::RegularCuboidMesh;
    kw...,
)
    numNode = femMesh.numNode
    numNode3 = numNode * 3
    faceNorm = femMesh.faceNorm
    surfNode = femMesh.surfNode
    K = femMesh.K

    if !haskey(kw, :uGamma)
        uIndex = SVector{9, Symbol}(
            [
                :x0y0z0
                :x0y1z0
                :x0y0z1
                :x0y1z1
                :y_x0z0
                :y_x0z1
                :z_x0y0
                :z_x0y1
                :yz_x0
            ],
        )
        uNodes = collect(Iterators.flatten([surfNode[uIndex[i]] for i in 1:length(uIndex)]))
        uGamma = BoundaryNode(; index = uIndex, node = uNodes)
    else
        uGamma = kw[:uGamma]
    end

    if !haskey(kw, :mGamma)
        mIndex = SVector{3, Symbol}([:x1y0z1; :x1y1z1; :y_x1z1])
        mNodes = collect(Iterators.flatten([surfNode[mIndex[i]] for i in 1:length(mIndex)]))
        mGamma = BoundaryNode(; index = mIndex, node = mNodes)
    else
        mGamma = kw[:mGamma]
    end

    if !haskey(kw, :tGamma)
        tIndex = SVector{14, Symbol}(
            Symbol.(setdiff(string.(keys(surfNode)), string.([uIndex; mIndex]))),
        )
        tNodes = collect(Iterators.flatten([surfNode[tIndex[i]] for i in 1:length(tIndex)]))
        tGamma = BoundaryNode(; index = tIndex, node = tNodes)
    else
        tGamma = kw[:tGamma]
    end

    uGammaNode = uGamma.node
    mGammaNode = mGamma.node
    tGammaNode = tGamma.node

    !haskey(kw, "uDofs") ?
    uDofs = [3 * uGammaNode .- 2; 3 * uGammaNode .- 1; 3 * uGammaNode; 3 * mGammaNode] :
    uDofs = kw["uDofs"]

    !haskey(kw, "tDofs") ? tDofs = setdiff(1:numNode3, uDofs) : tDofs = kw["tDofs"]

    !haskey(kw, "noExit") ? noExit = 4 : noExit = kw["noExit"]

    C = try
        factorize(K[tDofs, tDofs])
    catch
        K[tDofs, tDofs]
    end

    boundaries = Boundaries(;
        model = CantileverLoad(),
        noExit = noExit,
        uGamma = uGamma,
        tGamma = tGamma,
        mGamma = mGamma,
        uDofs = uDofs,
        tDofs = tDofs,
        mDofs = nothing,
        tK = C,
    )

    forceDisplacement = ForceDisplacement(;
        uTilde = spzeros(numNode3),
        uHat = spzeros(numNode3),
        u = spzeros(numNode3),
        fTilde = spzeros(numNode3),
        fHat = spzeros(numNode3),
        f = spzeros(numNode3),
    )
    return boundaries, forceDisplacement
end

function Boundaries(
    ::FEMParameters{T1, T2, T3, T4, T5} where {T1, T2, T3 <: PillarLoad, T4, T5},
    femMesh::RegularCuboidMesh;
    kw...,
)
    numNode = femMesh.numNode
    numNode3 = numNode * 3
    faceNorm = femMesh.faceNorm
    surfNode = femMesh.surfNode
    K = femMesh.K

    if !haskey(kw, :uGamma)
        uIndex = SVector{9, Symbol}(
            [
                :x0y0z0
                :x0y1z0
                :x0y0z1
                :x0y1z1
                :y_x0z0
                :y_x0z1
                :z_x0y0
                :z_x0y1
                :yz_x0
            ],
        )
        uNodes = collect(Iterators.flatten([surfNode[uIndex[i]] for i in 1:length(uIndex)]))
        uGamma = BoundaryNode(; index = uIndex, node = uNodes)
    else
        uGamma = kw[:uGamma]
    end

    if !haskey(kw, :mGamma)
        mIndex = SVector{9, Symbol}(
            [
                :x1y0z0
                :x1y1z0
                :x1y0z1
                :x1y1z1
                :y_x1z0
                :y_x1z1
                :z_x1y0
                :z_x1y1
                :yz_x1
            ],
        )
        mNodes = collect(Iterators.flatten([surfNode[mIndex[i]] for i in 1:length(mIndex)]))
        mGamma = BoundaryNode(; index = mIndex, node = mNodes)
    else
        mGamma = kw[:mGamma]
    end

    if !haskey(kw, :tGamma)
        tIndex = SVector{8, Symbol}(
            Symbol.(setdiff(string.(keys(surfNode)), string.([uIndex; mIndex]))),
        )
        tNodes = collect(Iterators.flatten([surfNode[tIndex[i]] for i in 1:length(tIndex)]))
        tGamma = BoundaryNode(; index = tIndex, node = tNodes)
    else
        tGamma = kw[:tGamma]
    end

    uGammaNode = uGamma.node
    mGammaNode = mGamma.node
    tGammaNode = tGamma.node

    !haskey(kw, "uDofs") ?
    uDofs =
        uDofs = [
            3 * uGammaNode .- 2 # x-coodinate for the yz plane at x = 0
            collect(
                Iterators.flatten([3 * surfNode[:x0y1z0] .- 1; 3 * surfNode[:x0y1z1] .- 1]),
            ) # y-coordinate
            collect(Iterators.flatten(3 * surfNode[:z_x0y1] .- 1)) # y-coordinate
            collect(Iterators.flatten([3 * surfNode[:x0y0z0]; 3 * surfNode[:x0y1z0]])) # z-coordinate
            collect(Iterators.flatten(3 * surfNode[:y_x0z0])) # z-coordinate
            3 * mGammaNode .- 2
        ] : uDofs = kw["uDofs"]

    !haskey(kw, "tDofs") ? tDofs = setdiff(1:numNode3, uDofs) : tDofs = kw["tDofs"]

    !haskey(kw, "noExit") ? noExit = SVector{2, Int}(2, 4) : noExit = kw["noExit"]

    C = try
        factorize(K[tDofs, tDofs])
    catch
        K[tDofs, tDofs]
    end

    boundaries = Boundaries(;
        model = PillarLoad(),
        noExit = noExit,
        uGamma = uGamma,
        tGamma = tGamma,
        mGamma = mGamma,
        uDofs = uDofs,
        tDofs = tDofs,
        mDofs = nothing,
        tK = C,
    )

    forceDisplacement = ForceDisplacement(;
        uTilde = spzeros(numNode3),
        uHat = spzeros(numNode3),
        u = spzeros(numNode3),
        fTilde = spzeros(numNode3),
        fHat = spzeros(numNode3),
        f = spzeros(numNode3),
    )
    return boundaries, forceDisplacement
end

"""
```
Boundaries(; model, noExit, uGamma, tGamma, mGamma, uDofs, tDofs, mDofs, tK)
```
Creates [`Boundaries`](@ref).
"""
function Boundaries(; model, noExit, uGamma, tGamma, mGamma, uDofs, tDofs, mDofs, tK)
    uGammaDln = [uGamma.node; mGamma.node]
    tGammaDln = [tGamma.node; mGamma.node]
    uDofsDln = [3 * uGammaDln .- 2; 3 * uGammaDln .- 1; 3 * uGammaDln]
    tDofsDln = [3 * tGammaDln .- 2; 3 * tGammaDln .- 1; 3 * tGammaDln]
    return Boundaries(
        model,
        noExit,
        uGammaDln,
        tGammaDln,
        uDofsDln,
        tDofsDln,
        uGamma,
        tGamma,
        mGamma,
        uDofs,
        tDofs,
        mDofs,
        tK,
    )
end

function getTGammaDlnNormsArea(boundaries::Boundaries, mesh::AbstractMesh)
    surfNode = mesh.surfNode
    surfNorm = mesh.surfNodeNorm
    surfArea = mesh.surfNodeArea
    tGamma = boundaries.tGamma
    mGamma = boundaries.mGamma
    index = [tGamma.index; mGamma.index]

    lenIdx = length(index)

    N = reshape(
        collect(
            Iterators.flatten(
                Iterators.flatten([
                    Iterators.repeated(surfNorm[index[i]], length(surfNode[index[i]]))
                    for i in 1:lenIdx
                ]),
            ),
        ),
        3,
        :,
    )
    A = collect(
        Iterators.flatten([
            Iterators.repeated(surfArea[index[i]], length(surfNode[index[i]])) for
            i in 1:lenIdx
        ]),
    )
    return N, A
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
function read_file(filename)
    data = read(open(filename, "r"), String)

    parsed_data = Vector{Vector{Int}}()
    for line in split(data, "\n")
        parsed_line = []
        for number in split(line, " ")
            if number != ""
                number = parse(Int, number)
                push!(parsed_line, number)
            end
        end
        push!(parsed_data, parsed_line)
    end
    parsed_data
end

function parse_file(d::Vector{Vector{Int}}, ds, id)
    i = 1
    (num_agents, num_jobs) = d[i]
    # return num_agents, num_jobs
    agent_job_cost_mat = hcat(d[i+1:i+num_agents]...)
    i += num_agents
    # return agent_job_cost_mat
    agent_job_resource_mat = hcat(d[i+1:i+num_agents]...)
    i += num_agents
    agents_resource_cap = d[i+1]
    GAPProb(
        GAPID(ds, num_agents, num_jobs),
        agent_job_cost_mat,
        agent_job_resource_mat,
        agents_resource_cap)
end

function flat_parse_file(d::Vector{Int}, ds, id)
    i = 0
    function ni(l)
        i += 1
        l[i]
    end

    num_agents = ni(d)
    num_jobs = ni(d)

    agent_job_costs = []
    for agent in 1:num_agents
        row = []
        for job in 1:num_jobs
            push!(row, ni(d))
        end
        push!(agent_job_costs, row)
    end

    agent_job_resources = []
    for agent in 1:num_agents
        row = []
        for job in 1:num_jobs
            push!(row, ni(d))
        end
        push!(agent_job_resources, row)
    end

    agents_resource_cap = []
    for agent in 1:num_agents
        push!(agents_resource_cap, ni(d))
    end

    prep(l) = hcat(l...)
    GAPProb(
        GAPID(id, ds, num_agents, num_jobs),
        prep(agent_job_costs),
        prep(agent_job_resources),
        agents_resource_cap)
end

function load_folder(folder="GAP/benchmark_problems/")
    gaps = Vector{GAPProb}()
    id = 0
    for file in readdir(folder)
        id += 1
        dataset = file[1]
        data = read_file(joinpath(folder, file))
        gap = flat_parse_file(vcat(data...), dataset, id)
        push!(gaps, gap)
    end
    sort(gaps, by=x->x.id.name)
end

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
using RecipesBase

@recipe function f(z::Array{Polar{T}}) where T
    projection --> :polar
    angle.(z),abs.(z)
end

@recipe function f(z::Array{Spherical{T}};sphere=true) where T
    delete!(plotattributes,:sphere)
    markersize --> 1
    aspect_ratio --> 1
    xlims --> (-1,1)
    ylims --> (-1,1)
    zlims --> (-1,1)

    @series begin
        x = [ cos(z.lat)*cos(z.lon) for z in z ]
        y = [ cos(z.lat)*sin(z.lon) for z in z ]
        z = [ sin(z.lat) for z in z ]
        x,y,z
    end

    function latcurves(n)
        lats = π*(1-n:2:n-1)/(2n+2)
        ϕ = π*(-200:200)/200
        [(cos(θ)*cos.(ϕ),cos(θ)*sin.(ϕ),fill(sin(θ),length(ϕ))) for θ in lats]
    end

    function loncurves(n)
        θ = π*(-100:100)/200
        longs = 2π*(0:n-1)/n
        [(cos.(θ)*cos(ϕ),cos.(θ)*sin(ϕ),sin.(θ)) for ϕ in longs]
    end

    if isa(sphere,Tuple)
        nlat,nlon = sphere
        sphere = true
    elseif isa(sphere,Int)
        nlat = nlon = sphere
        sphere = true
    else
        nlon = 12
        nlat = 7
    end
    if sphere
        for c in latcurves(nlat)
            @series begin
                group := 2
                seriestype := :path3d
                color := :gray
                linestyle := :dot
                linealpha := 0.5
                markershape := :none
                linewidth := 0.5
                label := ""
                c[1],c[2],c[3]
            end
        end

        for c in loncurves(nlon)
            @series begin
                group := 2
                seriestype := :path3d
                color := :gray
                linealpha := 0.5
                linestyle := :dot
                markershape := :none
                linewidth := 0.5
                label := ""
                c[1],c[2],c[3]
            end
        end
    end

end