text:
#=
# 3. Non-uniqueness and discrete circulation

There is no unique solution to the two-dimensional potential flow problem in the presence of an impenetrable surface, since one can choose any value for the uniform value $\mathfrak{s}_0$ and still enforce the no-penetration condition. Equivalently, we can specify any circulation about the body and still enforce this condition. For later use, let us write this uniform surface streamfunction as $\mathfrak{s}_0 = s_0 \mathfrak{1}$, where $s_0$ is a single scalar value and $\mathfrak{1}$ is a vector of ones on all surface points. The discrete circulation $\Gamma_b$ about the body is given by the sum of the bound vortex sheet data and can be written compactly as

$\Gamma_b = \mathfrak{1}^{T} \mathfrak{f}.$

In `GridPotentialFlow.jl`, the user cannot specify $s_0$ directly. Instead, the package solves the system

$\begin{bmatrix}
\mathsf{L} & \mathsf{R} & 0\\
\mathsf{E} &  0 & \mathfrak{1}\\
0 & \mathfrak{1}^{T} & 0\\
\end{bmatrix} \begin{pmatrix} \mathsf{s} \\ \mathfrak{f} \\ s_0 \end{pmatrix} =
\begin{pmatrix} -\mathsf{w} \\ \mathfrak{s}'_b \\ \Gamma_b \end{pmatrix}.$

and the circulation can be specified using the `Γ` keyword for each body. If it is not specified, it is set to zero. The return value of `solvesystem` has a field `ψ₀`, whose entries represent the value for $s_0$ for each body, which are in this case the Lagrange multipliers to enforce the specified circulation.
=#

# We now illustrate how to specify the circulation by simulating the flow around two cylinders that have different circulations.

#md # ```@setup 3.-Non-uniqueness-and-discrete-circulation
#md # using GridPotentialFlow
#md # using Plots
#md # ```
#!md using GridPotentialFlow
#!md using Plots
Δx = 0.01
Lx = 2.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx);

# The two cylinders are positioned using the tools from `RigidTransform` from `RigidBodyTools.jl` and the circulation of the first one is set to a non-zero value.
R = Lx/8;
b_left = PotentialFlowBody(Circle(R,Δx),Γ=1.0)
T = RigidTransform((-Lx/4,0.0),0.0)
T(b_left);
b_right = PotentialFlowBody(Circle(R,Δx))
T = RigidTransform((Lx/4,0.0),0.0)
T(b_right);

# A plot of the streamfunction then clearly shows the non-zero circulation about the left cylinder.
model = VortexModel(g,bodies=[b_left,b_right])
s = streamfunction(model)
plot(s,g)
plot!(b_left,fillcolor=:black,fillrange=0,fillalpha=0.25,linecolor=:black,linewidth=2)
plot!(b_right,fillcolor=:black,fillrange=0,fillalpha=0.25,linecolor=:black,linewidth=2)

#=
When the previously introduced saddle point system

$\begin{bmatrix}
\mathsf{L} & \mathsf{R} \\
\mathsf{E} &  0
\end{bmatrix} \begin{pmatrix} \mathsf{s} \\ \mathfrak{f} \end{pmatrix} =
\begin{pmatrix} -\mathsf{w} \\ \mathfrak{s}'_b - \mathfrak{s}_0 \end{pmatrix}.$

is solved, the resulting circulation is

$\Gamma_{b} = \mathfrak{1}^{T} \mathsf{S}^{-1}\left(\mathfrak{s}_{b}^{\prime}+\mathsf{EL}^{-1} \mathsf{w}\right)-s_{0} \mathfrak{1}^{T} \mathsf{S}^{-1} \mathfrak{1},$

where $\mathsf{S} = -\mathsf{EL}^{-1}\mathsf{R}$ is the Schur complement. The scalar factor $\mathfrak{1}^{T} \mathsf{S}^{-1} \mathfrak{1}$ in this expression is a property of the set of points and their immersion into the Cartesian grid. Part of this factor, $\mathsf{S}^{-1} \mathfrak{1}$, represents the bound vortex sheet strength associated with a uniform, unit-strength streamfunction on the surface. This sheet has a particularly important role in some of the discussion to follow, so we will denote its strength by $\mathfrak{f}_0$:

$\mathfrak{f}_0 \equiv \mathsf{S}^{-1} \mathfrak{1}.$

We demonstrate below the distribution of $\mathfrak{f}_0$ for elliptical cylinders with different aspect ratios.
=#

Δx = 0.004
Lx = 2.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx);

# Let us consider four cylinders with aspect ratios 1, 2, 3, and infinity. An ellipse with infinite aspect ratio is a flat plate, so we construct that case separately.
a = Lx/4 # semi-major axis
AR_list = [1.0,2.0,3.0]
ellipses = [Ellipse(a,a/AR,Δx) for AR in AR_list]
bodies = PotentialFlowBody[]
push!(bodies,[PotentialFlowBody(SplinedBody(hcat(e.x,e.y),3*cellsize(g))) for e in ellipses]...)
push!(bodies,PotentialFlowBody(Plate(2*a,3*Δx)));
colors = [:blue,:red,:green,:black]
plot()
for i in 1:length(bodies)
    plot!(bodies[i],fillrange=0,fillalpha=0.0,linecolor=colors[i],linewidth=2)
end
plot!(xlim=xlim,ylim=ylim,xlabel="x",ylabel="y")

# To compute $\mathfrak{f}_0$, we have to access the internal fields of the associated `system` of our vortexmodel.
f₀_list = []
for i in 1:length(bodies)
    Δs = dlength(bodies[i]);
    model = VortexModel(g,bodies=[bodies[i]]);
    ones = ScalarData(length(bodies[i]))
    ones .= 1.0
    f₀ = model.system.ibp.Sfact\ones
    push!(f₀_list,f₀)
end

# Because the flat plate is an open body, we need to mirror and append its contents to itself to compare it to the ellipses. Furthermore, its bound vortex sheet strength should be divided by two because it represents $\gamma = \gamma^{+} + \gamma^{-}$
append!(f₀_list[end],reverse!(f₀_list[end]));
Δs_plate = dlengthmid(bodies[end]);
append!(Δs_plate,reverse!(Δs_plate));
plot(ylim=(-2,0),xlabel="p/N",title="f₀ for ellipses of different aspect ratios")
for i in 1:length(bodies)-1
    plot!((1:length(f₀_list[i]))/length(f₀_list[i]),f₀_list[i]./dlength(bodies[i]),linecolor=colors[i],label="AR=$(AR_list[i])")
end
plot!((1:length(f₀_list[end]))/length(f₀_list[end]),0.5*f₀_list[end]./Δs_plate,linecolor=colors[end],label="AR=inf")


text:
using Dga
using Test

@testset "Dga.jl" begin

end

ext: jl


text:
using Random

mutable struct Chain{Int,Int} <: AbstractEnvironment
    state::Int
    reward::Float64
    done::Bool
    size::Int
    time_step::Int
	stochastic::Bool
	randchance::Float64
    function Chain(size::Int, stochastic::Bool=false, chance::Float64=0.2)
        new{Int,Int}(1, 0., false, size, 0, stochastic, chance)
    end
end

function get_gamma(env::Chain)
	return 1.0
end

function randomize_env(::Type{T}, rng::AbstractRNG, size::Int, stochastic::Bool) where {T <: Chain}
	ϵ = rand(rng, Uniform(0., 0.25))
	return Chain(size, stochastic, ϵ)
end

function get_state_desc(env::Chain)::Int
    return env.size
end

function get_action_desc(env::Chain)::Int
    return 2
end

function step!(env::Chain, action::Int, rng::AbstractRNG)::Float64
    env.time_step += 1
    if action <= 0 || action > 2
        error("Action needs to be an integer in [1, 2]")
    end

    if env.stochastic
		tmp = rand(rng)
		if tmp < env.randchance / 2.  # if random 50% chance to not move
			action = 0
		elseif tmp < env.randchance  # if random 50% chance to move in oposite direction
			action = 3 - action
		end
	end

    if action == 1
        env.state -= 1
    elseif action == 2
        env.state += 1
    end

	env.state = clamp(env.state, 1, env.size)

    env.reward = -1.0
	max_steps = 20 * env.size#^2
    env.done = (env.state == env.size) | (env.time_step ≥ max_steps)

	if env.done
		env.reward = 0.
	end

    return env.reward
end

function get_state(env::Chain)::Int
    return env.state
end

function get_reward(env::Chain)::Float64
    return env.reward
end

function is_terminal(env::Chain)::Bool
    return env.done
end

function env_reset!(env::Chain, rng)
	env.time_step = 0
    env.state = 1
    env.done = false
end

function clone(env::Chain)::Chain
    env2 = Chain(env.size, env.stochastic, env.randchance)
    env2.time_step = deepcopy(env.time_step)
    env2.state = deepcopy(env.state)
    env2.done = deepcopy(env.done)
    env2.reward = deepcopy(env.reward)
    return env2
end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
# m6.8s.jl

using Pkg, DrWatson

using MonteCarloMeasurements
using StanSample
using StatisticalRethinking

N = 100
df = DataFrame(
    :h0 => rand(Normal(10,2 ), N),
    :treatment => vcat(zeros(Int, Int(N/2)), ones(Int, Int(N/2)))
);
df[!, :fungus] = [rand(Binomial(1, 0.5 - 0.4 * df[i, :treatment]), 1)[1] for i in 1:N]
df[!, :h1] = [df[i, :h0] + rand(Normal(5 - 3 * df[i, :fungus]), 1)[1] for i in 1:N]

stan6_8 = "
data {
    int <lower=1> N;
    vector[N] h0;
    vector[N] h1;
    vector[N] treatment;
}
parameters{
    real a;
    real bt;
    real<lower=0> sigma;
}
model {
  vector[N] mu;
    vector[N] p;
    a ~ lognormal(0, 0.2);
    bt ~ normal(0, 0.5);
    sigma ~ exponential(1);
    for ( i in 1:N ) {
    p[i] = a + bt*treatment[i];
    mu[i] = h0[i] * p[i];
    }
    h1 ~ normal(mu, sigma);
}
"

m6_8_data = Dict(
    :N => nrow(df),
    :h0 => df[:, :h0],
    :h1 => df[:, :h1],
    :treatment => df[:, :treatment]
)

m6_8s = SampleModel("m6.8s", stan6_8)

rc6_8s = stan_sample(m6_8s; data=m6_8_data)

if success(rc6_8s)
    part6_8s = read_samples(m6_8s, :particles);
    part6_8s |> display
end

# End of m6.8s.jl


text:
##############################################################################
##
## Utilities for reading compound RDA items: lists, arrays etc
##
##############################################################################

function readdummy(ctx::RDAContext, fl::RDATag)
    # for reading elements without body,
    # e.g. NULL, empty environment etc
    RDummy{sxtype(fl)}()
end

function readattrs(ctx::RDAContext, fl::RDATag)
    if !hasattr(fl) return emptyhash
    else convert(Hash, readitem(ctx)) end
end

readnil = readdummy

function readnumeric(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == REALSXP
    RNumericVector(readfloatorNA(ctx.io, readlength(ctx.io)),
                   readattrs(ctx, fl))
end

function readinteger(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == INTSXP
    RIntegerVector(readintorNA(ctx.io, readlength(ctx.io)),
                   readattrs(ctx, fl))
end

function readlogical(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == LGLSXP # excluding this check, the method is the same as readinteger()
    RLogicalVector(readintorNA(ctx.io, readlength(ctx.io)),
                   readattrs(ctx, fl))
end

function readcomplex(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == CPLXSXP
    n = readlength(ctx.io)
    v = Vector{ComplexF64}(undef, n)
    readfloatorNA!(ctx.io, reinterpret(Float64, v))
    RComplexVector(v, readattrs(ctx, fl))
end

function readstring(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == STRSXP
    RStringVector(readcharacter(ctx.io, readlength(ctx.io))...,
                  readattrs(ctx, fl))
end

function readlist(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == VECSXP
    n = readlength(ctx.io)
    RList(RSEXPREC[readitem(ctx) for i in 1:n],
          readattrs(ctx, fl))
end

function readref(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == REFSXP
    ix = fl >> 8
    ix = ix != 0 ? ix : readint32(ctx.io)
    (ix <= length(ctx.ref_tab)) || throw(BoundsError("undefined reference index=$ix"))
    return ctx.ref_tab[ix]
end

function readsymbol(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == SYMSXP
    registerref!(ctx, RSymbol(readcharacter(ctx.io)))
end

function readS4(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == S4SXP
    # S4 objects doesn't contain any data other than its attributes
    # FIXME S4 is read just as a named hash of its attributes
    RS4Obj(readattrs(ctx, fl))
end

function readenv(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == ENVSXP
    is_locked = readint32(ctx.io)
    res = registerref!(ctx, REnvironment()) # registering before reading the contents
    res.enclosed = readitem(ctx)
    res.frame = readitem(ctx)
    res.hashtab = readitem(ctx)
    attr = readitem(ctx)
    if isa(attr, Hash)
        res.attr = attr
    end
    return res
end

function readname(ctx)
    (readint32(ctx.io) == zero(Int32)) ||
        error("Names in persistent strings not supported")
    n = readint32(ctx.io)
    return RString[readcharacter(ctx.io) for i in 1:n]
end

function readnamespace(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == NAMESPACESXP
    registerref!(ctx, RNamespace(readname(ctx)))
end

function readpackage(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == PACKAGESXP
    registerref!(ctx, RPackage(readname(ctx)))
end

# reads single-linked lists R objects
function readpairedobjects(ctx::RDAContext, fl::RDATag)
    res = RPairList(readattrs(ctx, fl))
    ifl = fl # RDATag for the list item
    while true
        if sxtype(ifl) == sxtype(fl)
            if hastag(ifl)
                tag = readitem(ctx)
                if isa(tag, RSymbol)
                    nm = tag.displayname
                else
                    nm = emptyhashkey
                end
            else
                nm = emptyhashkey
            end
            item = readitem(ctx)
            push!(res, item, nm)
            ifl = readuint32(ctx.io)
            # read the item attributes
            # FIXME what to do with the attributes?
            (sxtype(ifl) == sxtype(fl)) && readattrs(ctx, ifl)
        elseif sxtype(ifl) == NILVALUE_SXP # end of list
            break
        else # end of list (not a single-linked list item)
            # it's not clear whether it's an error of handling AltReps
            # or a feature of AltReps (it only occurs within AltReps)
            # normally pairlists should be terminated by NILVALUE_SXP
            @warn "0x$(string(ifl, base=16)) element in a 0x$(string(fl, base=16)) list, assuming it's the last element"
            item = readitem(ctx, ifl)
            push!(res, item, emptyhashkey)
            break
        end
    end

    return res
end

function readpairlist(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == LISTSXP
    readpairedobjects(ctx, fl)
end

function readlang(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == LANGSXP
    readpairedobjects(ctx, fl)
end

function readclosure(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == CLOSXP
    res = RClosure(readattrs(ctx, fl))
    hastag(fl) && (res.env = readitem(ctx))
    res.formals = readitem(ctx)
    res.body = readitem(ctx)
    return res
end

function readpromise(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == PROMSXP
    res = RPromise(readattrs(ctx, fl))
    hastag(fl) && (res.env = readitem(ctx))
    res.value = readitem(ctx)
    res.expr = readitem(ctx)
    return res
end

function readraw(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == RAWSXP
    n = readlength(ctx.io)
    RRaw(readuint8(ctx.io, n), readattrs(ctx, fl))
end

# reads built-in R objects
function readbuiltin(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == BUILTINSXP || sxtype(fl) == SPECIALSXP
    RBuiltin(readnchars(ctx.io, readint32(ctx.io)))
end

function readextptr(ctx::RDAContext, fl::RDATag)
    @assert sxtype(fl) == EXTPTRSXP
    res = registerref!(ctx, RExtPtr()) # registering before reading the contents
    res.protected = readitem(ctx)
    res.tag = readitem(ctx)
    res.attr = readattrs(ctx, fl)
    return res
end

"""
Context for reading R bytecode.
"""
struct BytecodeContext
    ctx::RDAContext                            # parent RDA context
    ref_tab::Vector{Union{RBytecode, Missing}} # table of bytecode references

    BytecodeContext(ctx::RDAContext, nrefs::Int32) =
        new(ctx, Vector{Union{RBytecode, Missing}}(missing, nrefs))
end

const BYTECODELANG_Types = Set([BCREPREF, BCREPDEF, LANGSXP, LISTSXP, ATTRLANGSXP, ATTRLISTSXP])

function readbytecodelang(bctx::BytecodeContext, bctype::Int32)
    if bctype == BCREPREF # refer to an already defined bytecode
        res = bctx.ref_tab[readint32(bctx.ctx.io)+1]
        @assert !ismissing(res)
        return res
    elseif bctype ∈ BYTECODELANG_Types
        pos = 0
        hasattr = false
        if bctype == BCREPDEF # define a reference
            pos = readint32(bctx.ctx.io)+1
            bctype = readint32(bctx.ctx.io)
        end # now bctype might be updated
        if bctype == ATTRLANGSXP
            bctype = LANGSXP # it's an R expression
            hasattr = true
        elseif bctype == ATTRLISTSXP
            bctype = LISTSXP # it's pairlist
            hasattr = true
        end
        res = RBytecode() # FIXME create a RPairlist if LANG/LISTSXP?
        hasattr && (res.attr = readitem(bctx.ctx))
        res.tag = readitem(bctx.ctx)
        res.car = readbytecodelang(bctx, readint32(bctx.ctx.io))
        res.cdr = readbytecodelang(bctx, readint32(bctx.ctx.io))
        (pos > 0) && setindex!(bctx.ref_tab, res, pos)
        return res
    else
        return readitem(bctx.ctx)
    end
end

function readbytecodeconsts(bctx::BytecodeContext)
    nconsts = readint32(bctx.ctx.io)
    v = fill!(Vector{RSEXPREC}(undef, nconsts), RDummy{NILVALUE_SXP}())
    @inbounds for i = 1:nconsts
        bctype = readint32(bctx.ctx.io)
        v[i] = if bctype == BCODESXP
            readbytecodecontents(bctx)
        elseif bctype ∈ BYTECODELANG_Types
            readbytecodelang(bctx, bctype)
        else
            readitem(bctx.ctx)
        end
    end
    return RList(v)
end

function readbytecodecontents(bctx::BytecodeContext)
    RBytecode(readitem(bctx.ctx),
              readbytecodeconsts(bctx))
end

function readbytecode(ctx::RDAContext, fl::RDATag)
    @assert fl == BCODESXP
    res = readbytecodecontents(BytecodeContext(ctx, readint32(ctx.io)))
    res.attr = readattrs(ctx, fl)
    return res
end

function readaltrep(ctx::RDAContext, fl::RDATag)
    info = readitem(ctx)
    state = readitem(ctx)
    attr = readitem(ctx)
    return RAltRep(info, state, isa(attr, RDummy) ? emptyhash : attr)
end

function readunsupported(ctx::RDAContext, fl::RDATag)
    throw(UnsupportedROBJ(sxtype(fl), "Reading SEXPREC of type $(sxtype(fl)) ($(SXTypes[sxtype(fl)].name)) is not supported"))
end

"""
Definition of R type.
"""
struct SXTypeInfo
    name::String         # R type name
    reader::Function     # function to deserialize R type from RDA stream
end

"""
Maps R type id (`SXType`) to its `SXTypeInfo`.
"""
const SXTypes = Dict{SXType, SXTypeInfo}(
    NILSXP     => SXTypeInfo("NULL",readdummy),
    SYMSXP     => SXTypeInfo("Symbol",readsymbol),
    LISTSXP    => SXTypeInfo("Pairlist",readpairlist),
    CLOSXP     => SXTypeInfo("Closure",readclosure),
    ENVSXP     => SXTypeInfo("Environment",readenv),
    PROMSXP    => SXTypeInfo("Promise",readpromise),
    LANGSXP    => SXTypeInfo("Lang",readlang),
    SPECIALSXP => SXTypeInfo("Special",readbuiltin),
    BUILTINSXP => SXTypeInfo("Builtin",readbuiltin),
    CHARSXP    => SXTypeInfo("Char",readunsupported),
    LGLSXP     => SXTypeInfo("Logical",readlogical),
    INTSXP     => SXTypeInfo("Integer",readinteger),
    REALSXP    => SXTypeInfo("Real",readnumeric),
    CPLXSXP    => SXTypeInfo("Complex",readcomplex),
    STRSXP     => SXTypeInfo("String",readstring),
    DOTSXP     => SXTypeInfo("Dot",readunsupported),
    ANYSXP     => SXTypeInfo("Any",readunsupported),
    VECSXP     => SXTypeInfo("List",readlist),
    EXPRSXP    => SXTypeInfo("Expr",readunsupported),
    BCODESXP   => SXTypeInfo("ByteCode",readbytecode),
    EXTPTRSXP  => SXTypeInfo("XPtr",readextptr),
    WEAKREFSXP => SXTypeInfo("WeakRef",readunsupported),
    RAWSXP     => SXTypeInfo("Raw",readraw),
    S4SXP      => SXTypeInfo("S4",readS4),
    NEWSXP     => SXTypeInfo("New",readunsupported),
    FREESXP    => SXTypeInfo("Free",readunsupported),
    FUNSXP     => SXTypeInfo("Function",readunsupported),
    BASEENV_SXP       => SXTypeInfo("BaseEnv",readdummy),
    EMPTYENV_SXP      => SXTypeInfo("EmptyEnv",readdummy),
    BCREPREF          => SXTypeInfo("BCREPREF",readunsupported), # handled within readbytecode()
    BCREPDEF          => SXTypeInfo("BCREPDEF",readunsupported), # handled within readbytecode()
    GENERICREFSXP     => SXTypeInfo("GenericRef",readunsupported),
    CLASSREFSXP       => SXTypeInfo("ClassRef",readunsupported),
    PERSISTSXP        => SXTypeInfo("Persist",readunsupported),
    PACKAGESXP        => SXTypeInfo("Package",readpackage),
    NAMESPACESXP      => SXTypeInfo("Namespace",readnamespace),
    BASENAMESPACE_SXP => SXTypeInfo("BaseNamespace",readdummy),
    MISSINGARG_SXP    => SXTypeInfo("MissingArg",readdummy),
    UNBOUNDVALUE_SXP  => SXTypeInfo("UnboundValue",readdummy),
    GLOBALENV_SXP     => SXTypeInfo("GlobalEnv",readdummy),
    NILVALUE_SXP      => SXTypeInfo("NilValue",readnil),
    REFSXP            => SXTypeInfo("Ref",readref),
    ALTREP_SXP        => SXTypeInfo("AltRep",readaltrep)
)

function readitem(ctx::RDAContext, fl::RDATag)
    sxt = sxtype(fl)
    haskey(SXTypes, sxt) || throw(UnsupportedROBJ(sxt, "encountered unknown SEXPREC type $sxt"))
    sxtinfo = SXTypes[sxt]
    return sxtinfo.reader(ctx, fl)
### Should not occur at the top level
###    if sxt == NILVALUE_SXP return nothing end      # terminates dotted pair lists
###    if sxt == CHARSXP return readcharacter(ctx.io, ff) end
end

readitem(ctx::RDAContext) = readitem(ctx, readuint32(ctx.io))


text:
"""
    @alias new = original

Define `const new = original` and attach `original`'s docstring to `new`.
"""
macro alias(expr)
    expr.head == :(=) || error("must be an assignment expression")
    new, original = expr.args
    return quote
        @doc (@doc $original)
        const $(esc(new)) = $(esc(original))
    end
end

"""
    @withfb "Reticulating splines" slow_function()

Give user feedback: what is happening when the program hangs; and when is it done.
"""
macro withfb(description, expr)
    return quote
        print($(esc(description)), " … ")
        flush(stdout)
        t0 = time()
        $(esc(expr))
        dt = time() - t0
        print("done")
        if dt > 0.1  # seconds
            @printf " (%.1f s)" dt
        end
        println()
    end
end


text:
using Revise
using AdFem
using PyCall
using LinearAlgebra
using PyPlot
using SparseArrays
using MAT
using Statistics
np = pyimport("numpy")

model_id = 2
if length(ARGS)==1
  global model_id = parse(Int64, ARGS[1])
end

n = 15
m = 2n
h = 0.01
NT = 20
Δt = 2/NT

bdedge = bcedge("right", m, n, h)
bdnode = bcnode("lower", m, n, h)

λ = constant(2.0)
μ = constant(0.2)

M = compute_fem_mass_matrix1(m, n, h)
Zero = spzeros((m+1)*(n+1), (m+1)*(n+1))
M = SparseTensor([M Zero;Zero M])

## alpha-scheme
β = 1/4; γ = 1/2

function eta_fun(σ)
  if model_id==1
    return constant(10*ones(4*m*n)) + 5.0/(1+1000*sum(σ[:,1:2]^2, dims=2))
  else
    return constant(10*ones(4*m*n)) + relu(50.0/(1+1000*sum(σ[:,1:2]^2, dims=2)) - 10.0)
  end
end

# invη is a 4*m*n array
function make_matrix(invη)
  a = b = 0.1
  fn_G = invη->begin
    G = tensor([1/Δt+2/3*μ*invη -μ/3*invη 0.0
      -μ/3*invη 1/Δt+2/3*μ*invη 0.0
      0.0 0.0 1/Δt+μ*invη])
    invG = inv(G)
  end
  invG = map(fn_G, invη)
  S = tensor([2μ/Δt+λ/Δt λ/Δt 0.0
      λ/Δt 2μ/Δt+λ/Δt 0.0
      0.0 0.0 μ/Δt])

  H = invG*S
  K = compute_fem_stiffness_matrix(H, m, n, h)
  C = a*M + b*K # damping matrix
  L = M + γ*Δt*C + β*Δt^2*K
  # L = pl[1]*L
  L, Lbd = fem_impose_Dirichlet_boundary_condition_experimental(L, bdnode, m, n, h)

  return C, K, L, S, invG
end

function simulate(FORCE_SCALE)

  a = TensorArray(NT+1); a = write(a, 1, zeros(2(m+1)*(n+1))|>constant)
  v = TensorArray(NT+1); v = write(v, 1, zeros(2(m+1)*(n+1))|>constant)
  d = TensorArray(NT+1); d = write(d, 1, zeros(2(m+1)*(n+1))|>constant)
  U = TensorArray(NT+1); U = write(U, 1, zeros(2(m+1)*(n+1))|>constant)
  Sigma = TensorArray(NT+1); Sigma = write(Sigma, 1, zeros(4*m*n, 3)|>constant)
  Varepsilon = TensorArray(NT+1); Varepsilon = write(Varepsilon, 1,zeros(4*m*n, 3)|>constant)

  Forces = zeros(NT, 2(m+1)*(n+1))
  for i = 1:NT
    T = eval_f_on_boundary_edge((x,y)->0.1*FORCE_SCALE, bdedge, m, n, h)

    T = [-T T]
    rhs = compute_fem_traction_term(T, bdedge, m, n, h)

    Forces[i, :] = rhs
  end
  Forces = constant(Forces)

  function condition(i, tas...)
    i <= NT
  end

  function body(i, tas...)
    a_, v_, d_, U_, Sigma_, Varepsilon_ = tas
    a = read(a_, i)
    v = read(v_, i)
    d = read(d_, i)
    U = read(U_, i)
    Sigma = read(Sigma_, i)
    Varepsilon = read(Varepsilon_, i)

    invη = eta_fun(Sigma)
    C, K, L, S, invG = make_matrix(invη)

    res = batch_matmul(invG/Δt, Sigma)
    F = compute_strain_energy_term(res, m, n, h) - K * U
    rhs = Forces[i] - F

    td = d + Δt*v + Δt^2/2*(1-2β)*a
    tv = v + (1-γ)*Δt*a
    rhs = rhs - C*tv - K*td
    rhs = scatter_update(rhs, constant([bdnode; bdnode.+(m+1)*(n+1)]), constant(zeros(2*length(bdnode))))

    ## alpha-scheme
    a = L\rhs # bottleneck
    d = td + β*Δt^2*a
    v = tv + γ*Δt*a
    U_new = d

    Varepsilon_new = eval_strain_on_gauss_pts(U_new, m, n, h)

    res2 = batch_matmul(invG * S, Varepsilon_new-Varepsilon)
    Sigma_new = res +  res2

    i+1, write(a_, i+1, a), write(v_, i+1, v), write(d_, i+1, d), write(U_, i+1, U_new),
          write(Sigma_, i+1, Sigma_new), write(Varepsilon_, i+1, Varepsilon_new)
  end

  i = constant(1, dtype=Int32)
  _, _, _, _, u, sigma, varepsilon = while_loop(condition, body,
                    [i, a, v, d, U, Sigma, Varepsilon])

  U = stack(u)
  Sigma = stack(sigma)
  Varepsilon = stack(varepsilon)

  U = set_shape(U, (NT+1, size(U,2)))
  Sigma = set_shape(Sigma, (NT+1, 4*m*n, 3))

  return U, Sigma
end

function visualize(U0, U_, Sigma0, Sigma_)
  close("all")
  figure(figsize=(13,4))
  subplot(121)
  plot(LinRange(0, 20, NT+1), U0[:,1], "r--")  # reference
  plot(LinRange(0, 20, NT+1), U_[:,1], "r")
  plot(LinRange(0, 20, NT+1), U0[:,1+(n+1)*(m+1)], "g--")
  plot(LinRange(0, 20, NT+1), U_[:,1+(n+1)*(m+1)], "g")
  xlabel("Time")
  ylabel("Displacement")
  subplot(122)
  plot(LinRange(0, 20, NT+1), mean(Sigma0[:,1:4,1], dims=2)[:],"r--", label="\$\\sigma_{xx}\$")
  plot(LinRange(0, 20, NT+1), mean(Sigma0[:,1:4,2], dims=2)[:],"b--", label="\$\\sigma_{yy}\$")
  plot(LinRange(0, 20, NT+1), mean(Sigma0[:,1:4,3], dims=2)[:],"g--", label="\$\\sigma_{xy}\$")
  plot(LinRange(0, 20, NT+1), mean(Sigma_[:,1:4,1], dims=2)[:],"r-")
  plot(LinRange(0, 20, NT+1), mean(Sigma_[:,1:4,2], dims=2)[:],"b-")
  plot(LinRange(0, 20, NT+1), mean(Sigma_[:,1:4,3], dims=2)[:],"g-")
  legend()
  legend()
  xlabel("Time")
  ylabel("Stress")
  savefig("disp$i.png")
  matwrite("nn$i.mat", Dict("U"=>U_, "S"=>Sigma_))
end

U = Array{PyObject}(undef, 4)
Sigma = Array{PyObject}(undef, 4)

for (k,FORCE_SCALE) in enumerate([0.5, 0.8, 1.5, 1.0])
  U[k], Sigma[k] = simulate(FORCE_SCALE)
end
sess = Session(); init(sess)

Uval = run(sess, U)
Sigmaval = run(sess, Sigma)
matwrite("data$model_id.mat", Dict("U"=>Uval, "S"=>Sigmaval))

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
function en(an,anm1,Ln)
	return sqrt((1-abs2(an))/((one(an)-abs2(anm1))*(one(an)-abs2(Ln))))
end

function ln(z::Vector{Complex128},an::Complex128,anm1::Complex128,phinm1::Vector{Complex128},phinm1star::Vector{Complex128})
    f1 = [((z[i]-anm1)/(one(Complex128)-conj(an)*z[i]))*conj(phinm1[i]) for i=1:length(z)]
    f2 = [((one(Complex128)-conj(anm1)*z[i])/(one(Complex128)-conj(an)*z[i]))*conj(phinm1star[i]) for i=1:length(z)]
    num = sum(f1)
    den = sum(f2)
    return -num/den
end

function ln(z::Vector{Complex128},a::Vector{Complex128},phinm1::Vector{Complex128})
    if length(a)==1
        an = copy(a[end])
        anm1 = zero(Complex128)
        phinm1star = phinm1
        f1 = [((z[i]-anm1)/(one(Complex128)-conj(an)*z[i]))*conj(phinm1[i]) for i=1:length(z)]
        f2 = [((one(Complex128)-conj(anm1)*z[i])/(one(Complex128)-conj(an)*z[i]))*phinm1star[i] for i=1:length(z)]
        num = sum(f1)
        den = sum(f2)
        r = -num/den
    else
        an = copy(a[end])
        anm1 = copy(a[end-1])
        bProd = length(a)>=3 ? conjBlaschkeProd(z,[a[1:end-2]]) : ones(Complex128,length(z))
        f1 = [((z[i]-anm1)/(one(Complex128)-conj(an)*z[i]))*conj(phinm1[i]) for i=1:length(z)]
        f2 = [((z[i]-anm1)/(one(Complex128)-conj(an)*z[i]))*phinm1[i]*bProd[i] for i=1:length(z)]
        num = sum(f1)
        den = sum(f2)
        r = -num/den
    end
    return r
end

function phiNext(z::Vector{Complex128},a::Vector{Complex128},phinm1::Vector{Complex128})
    phin = zeros(Complex128,length(phinm1))
    if length(a)==1
        an = copy(a[end])
        anm1 = zero(Complex128)
        L = ln(z,a,phinm1)
        c = en(an,anm1,L)
        for k=1:length(z)
            q = (z[k]-anm1)
            p = (one(Complex128)-conj(anm1)*z[k])
            C = c/(z[k]-an)
            PHI = p*phinm1[k]
            PHISTAR = q*phinm1[k]
            phin[k] = C*(PHI+conj(L)*PHISTAR)
        end
    else
        an = a[end]
        anm1 = a[end-1]
        L = ln(z,a,phinm1)
        c = en(an,anm1,L)
        bProd = length(a)>=3 ? blaschkeProd(z,[a[1:end-2]]) : ones(Complex128,length(z))
        for k=1:length(z)
            q = (z[k]-an)
            p = (one(Complex128)-conj(anm1)*z[k])
            C = c*p/q
            phin[k] = C*(phinm1[k]+conj(L)*bProd[k]*conj(phinm1[k]))
        end
    end
    return (phin, L)
end

function phiNext(z::Vector{Complex128},an::Complex128,anm1::Complex128,phinm1::Vector{Complex128},phinm1star::Vector{Complex128})
    phin = zeros(Complex128,length(phinm1))
    phinstar = zeros(Complex128,length(phinm1))
    L = ln(z,an,anm1,phinm1,phinm1star)
    c = en(an,anm1,L)
    for k=1:length(z)
        q = (z[k]-anm1)
        p = (one(Complex128)-conj(anm1)*z[k])
        C = c/(z[k]-an)
        PHI = p*phinm1[k]
        PHISTAR = q*phinm1star[k]
        phin[k] = C*(PHI+conj(L)*PHISTAR)
        phinstar[k] = C*(L*PHI+PHISTAR)
    end
    return (phin, phinstar, L)
end

function phiNext(z::Vector{Complex128},an::Complex128,anm1::Complex128,phinm1::Vector{Complex128},phinm1star::Vector{Complex128},L::Complex128)
    phin = zeros(Complex128,length(phinm1))
    phinstar = zeros(Complex128,length(phinm1))
    c = en(an,anm1,L)
    for k=1:length(z)
        q = (z[k]-anm1)
        p = (one(Complex128)-conj(anm1)*z[k])
        C = c/(z[k]-an)
        PHI = p*phinm1[k]
        PHISTAR = q*phinm1star[k]
        phin[k] = C*(PHI+conj(L)*PHISTAR)
        phinstar[k] = C*(L*PHI+PHISTAR)
    end
    return (phin, phinstar, L)
end

function blaschkeProd(z::Vector{Complex128},aj::Vector{Complex128})
    r = ones(Complex128,length(z))
    for i=1:length(aj)
        t = [(one(Complex128)-conj(aj[i])*z[k])./(z[k]-aj[i]) for k=1:length(z)]
        r = r.*t
    end
    return r
end

function conjBlaschkeProd(z::Vector{Complex128},aj::Vector{Complex128})
    r = ones(Complex128,length(z))
    for i=1:length(aj)
        t = [(z[k]-aj[i])/(one(Complex128)-conj(aj[i])*z[k]) for k=1:length(z)]
        r = r.*t
    end
    return r
end

function qrUpdate(z::Vector{Complex128},p::Vector{Complex128},B::Array{Complex128,2},Bstar::Array{Complex128,2})
    (L,C) = size(B)
    (Q,R) = qr(B)
    Lk = zeros(Complex128,C-1)
    Bnew = copy(Q)
    BstarNew = zeros(Complex128,(L,C))
    for k=2:C
        BstarNew[:,k] = conj(Bnew[:,k]).*blaschkeProd(z,[p[1:k]])
        Lk[k-1] = ln(z,p[k],p[k-1],Bnew[:,k-1],BstarNew[:,k-1])
    end
    return (Bnew[:,2:end],BstarNew[:,2:end],Lk)
end

function genZ(N::Int)
	return exp((im*2*pi/N)*[0:N-1])
end

function genZArc(N::Int,Div::Int)
    return exp((im*2*pi/(Div*N))*[0:N-1])
end

function genBasis(N::Int,z::Vector,a::Vector{Complex128})
    M = length(a)+1
    upfreq = 500;
    a = cat(1,zero(Complex128),a)
    B = ones(Complex128,(N,M))/sqrt(N)
    Bstar = ones(Complex128,(N,M))/sqrt(N)
    #BstarOld = ones(Complex128,N)/sqrt(N)
    Lk = zeros(Complex128,M-1)
    for k=2:M
        (B[:,k],Bstar[:,k],Lk[k-1]) = phiNext(z,a[k],a[k-1],B[:,k-1],Bstar[:,k-1])
        #(B[:,k],BstarNew,Lk[k-1]) = phiNext(z,a[k],a[k-1],B[:,k-1],BstarOld)
        #BstarOld = copy(BstarNew)
        if mod(k+1,upfreq)==0
            println("QR updating")
            (B[:,(k-upfreq+1):k],Bstar[:,(k-upfreq+1):k],Lk[(k-upfreq+1):k]) = qrUpdate(z,a,B[:,(k-upfreq):k],Bstar[:,(k-upfreq):k])
        end
    end
    return (B[:,2:end],Lk,Bstar[:,2:end])
end

function genBasis(N::Int,z::Vector,a::Vector{Complex128},Lk::Vector{Complex128},InitVal)
    M = length(a)+1
    a = cat(1,zero(Complex128),a)
    B = ones(Complex128,(N,M))*InitVal
    BstarOld = ones(Complex128,N)*InitVal
    for k=2:M
        (B[:,k],BstarNew) = phiNext(z,a[k],a[k-1],B[:,k-1],BstarOld,Lk[k-1])
        BstarOld = copy(BstarNew)
    end
    return B[:,2:end]
end

function genWBasis(N::Int,a::Vector{Complex128})
    z = genZ(N)
    B = genBasis(N,z,a)
    return fftshift(B,1)
end

function genTBasis(N::Int,a::Vector{Complex128})
    z = genZ(N)
	(B,L) = genBasis(N,z,a)
	T = fft(B,1)/sqrt(N)
    T = flipud(T)
    return T
end

function unitDiskGrid(n::Int,delta::Float64)
    N = nextpow2(n)
    M = N>>1
    v = linspace(delta,one(Float64)-delta,M)
    Zjk = zeros(Complex128,N,N)
    for j=1:(M)
            #bottom, running counter-clockwise
            Zjk[(M+j),((M-(j-1)):(M+1+(j-1)))] = v[j]*exp(im*linspace(5*pi/4,7*pi/4,2*j))
            #top, running clockwise
            Zjk[(M-j+1),((M-(j-1)):(M+1+(j-1)))] = v[j]*exp(im*linspace(3*pi/4,1*pi/4,2*j))
    end
    for j=1:(M)
            #right, running clockwise
            Zjk[((M-(j-1)):(M+1+(j-1))),M+j] = v[j]*exp(im*linspace(pi/4,-pi/4,2*j)) #right
            #left, running counter-clockwise
            Zjk[((M-(j-1)):(M+1+(j-1))),M-j+1] = v[j]*exp(im*linspace(3*pi/4,5*pi/4,2*j)) #left
    end
    return (reshape(Zjk,N^2))
end

function ompTMFreqDomain(D::Array{Complex128,2},Z::Vector{Complex128},Npoles::Int,Npad::Int)
    D = cat(1,zeros(Npad,size(D,2)),D,zeros(Npad,size(D,2)))
    N = size(D,1)
    Q = size(D,2)
    M = length(Z)
    z = genZ(N)
    B = zeros(Complex128,(N,M))
    Bstar = zeros(Complex128,(N,M))
    Bnm1 = ones(Complex128,(N,M))/sqrt(N)
    Bstarnm1 = ones(Complex128,(N,M))/sqrt(N)
    poles = zeros(Complex128,Npoles+1)
    amps = zeros(Complex128,((Npoles+1),Q))
    jopt = one(Int);
    R = copy(D)
    for k=2:(Npoles+1)
        for j=1:M
            (B[:,j],Bstar[:,j]) = phiNext(z,Z[j],poles[k-1],Bnm1[:,jopt],Bstarnm1[:,jopt])
        end
        Bnm1 = copy(B)
        Bstarnm1 = copy(Bstar)
        B = fftshift(B,1)
        P = B'*R #projection on to column space of data
        O = sum(abs2(P),2) #this is the vector we wish to maximize
        jopt = indmax(O)
        poles[k] = Z[jopt]
        amps[k,:] = P[jopt,:]
        R -= B[:,jopt]*amps[k,:]
    end
    return (poles[2:end], amps[2:end,:])
end

function ompTMTimeDomain(D::Array{Complex128,2},Z::Vector{Complex128},Npoles::Int,Npad::Int)
    D = cat(1,zeros(Npad,size(D,2)),D,zeros(Npad,size(D,2)))
    N = size(D,1)
    Q = size(D,2)
    M = length(Z)
    z = genZ(N)
    B = zeros(Complex128,(N,M))
    Bstar = zeros(Complex128,(Npoles,M))
    Ltemp = zeros(Complex128,M)
    Lopt = zeros(Complex128,Npoles)
    Bnm1 = ones(Complex128,(N,M))/sqrt(N)
    Bstarnm1 = ones(Complex128,(N,M))/sqrt(N)
    poles = zeros(Complex128,Npoles+1)
    amps = zeros(Complex128,((Npoles+1),Q))
    jopt = ones(Int,Npoles+1);
    R = copy(D)
    for k=2:(Npoles+1)
        for j=1:M
            (B[:,j],Bstar[:,j],Ltemp[j]) = phiNext(z,Z[j],poles[k-1],Bnm1[:,jopt[k-1]],Bstarnm1[:,jopt[k-1]])
        end
        Bnm1 = copy(B)
        Bstarnm1 = copy(Bstar)
        B = flipud(fft(B,1))/sqrt(N)
        P = B'*R #projection on to column space of data
        O = sum(abs2(P),2) #this is the vector we wish to maximize
        O[jopt[1:(k-1)]] = 0 #prevent repeat pole selection
        jopt[k] = indmax(O)
        poles[k] = Z[jopt[k]]
        amps[k,:] = P[jopt[k],:]
        R -= B[:,jopt[k]]*amps[k,:]
        println(Ltemp[jopt[k]])
        Lopt[k-1] = Ltemp[jopt[k]]
    end
    return (poles[2:end], amps[2:end,:], Lopt)
end

function ompSepTMDecomp(Npoles::Int,Nsing::Int,Ngrid::Int,Ntarget::Int,D::Array{Complex128,2})
    (u,s,v) = svd(D)
    vp = v[:,1:Nsing]
    sp = s[1:Nsing]
    up  = u[:,1:Nsing]
    Z = unitDiskGrid(int(floor(sqrt(Ngrid))),0.001)
    (polesx,ampsx,Lk) = ompTMTimeDomain(vp,Z,Npoles,0)
    Bx = genTBasis(Ntarget,polesx)
    R = up*diagm(sp)*(Bx*ampsx)'
    return (polesx,ampsx,R,Lk)
end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

using ACE, JuLIP, BenchmarkTools

#---

basis = ACE.Utils.rpi_basis(; species=:Si, N = 5, maxdeg = 14)
@show length(basis)
V = ACE.Random.randcombine(basis)

#---

at = bulk(:Si, cubic=true) * 5
nlist = neighbourlist(at, cutoff(V))
tmp = JuLIP.alloc_temp(V, maxneigs(nlist))
tmp_d = JuLIP.alloc_temp_d(V, maxneigs(nlist))
F = zeros(JVecF, length(at))
@code_warntype JuLIP.energy!(tmp, V, at)
@code_warntype JuLIP.forces!(F, tmp_d, V, at)

@btime JuLIP.forces!($F, $tmp_d, $V, $at)
@btime JuLIP.forces($V, $at)


text:
using BenchmarkTools
using Distributions

d = Normal()
X = rand(d, 5000)

@btime pdf.(d, X)
@btime logpdf.(d, X)

using Statistikk

@btime Statistikk.pdf(d, X)
@btime Statistikk.logpdf(d, X)

max_stars_repo_name: maxmouchet/Statistikk

max_issues_repo_name: maxmouchet/Statistikk

max_forks_repo_name: maxmouchet/Statistikk


text:
struct BernoulliArray <: Distribution{Array{Bool}} end
const bernoulli_array = BernoulliArray()

function Gen.logpdf(::BernoulliArray, x::Array{Bool}, weights::Array{Float64})
    heads = x
    tails =.!x
    sum(log.(weights[heads])) + sum(log.(1.0 .- weights[tails]))
end

function Gen.random(::BernoulliArray, weights::Array{Float64})
    [random(bernoulli,p) for p in weights]
end

export bernoulli_array


text:
using IntervalContractors
using IntervalArithmetic

@testset "forward_backward_contractor" begin

    vars = @variables x, y
    @variables a

    ex = x^2 + y^2

    C = forward_backward_contractor(ex, vars)
    @test C(IntervalBox(-10..10, -10..10), 0..1) == ( (-1..1, -1..1), 0..200 )

    ex = x^2 + a * y^2

    C = forward_backward_contractor(ex, vars, [a])
    @test C(IntervalBox(-10..10, -10..10), 0..1, 1..1) == ( (-1..1, -1..1), 0..200 )

end


text:
using Molecules.PointGroup

mols = ["water", "ammonia", "benzene", "ethanol", "formaldehyde", "glycine", "methane", "phosphaethene", "acetic_acid", "nitrogen"]
Ds = Any[Molecules.PointGroup.DistanceMatrix(Any[Atom{Int16, Float64}(8, 15.999, [1.2091536548, 1.7664118189, -0.0171613972]), Atom{Int16, Float64}(1, 1.008, [2.1984800075, 1.7977100627, 0.0121161719]), Atom{Int16, Float64}(1, 1.008, [0.9197881882, 2.458018557, 0.6297938832])], [0.0 0.9902542038608343 0.9902542038706021; 0.9902542038608343 0.0 1.5660733799622173; 0.9902542038706021 1.5660733799622173 0.0]), Molecules.PointGroup.DistanceMatrix(Any[Atom{Int16, Float64}(7, 14.007, [0.0, 0.0, 0.117347]), Atom{Int16, Float64}(1, 1.008, [0.0, 0.932649, -0.273809]), Atom{Int16, Float64}(1, 1.008, [0.807698, -0.466325, -0.273809]), Atom{Int16, Float64}(1, 1.008, [-0.807698, -0.466325, -0.273809])], [0.0 1.0113541286498018 1.011354577368887 1.011354577368887; 1.0113541286498018 0.0 1.6153960232339313 1.6153960232339313; 1.011354577368887 1.6153960232339313 0.0 1.615396; 1.011354577368887 1.6153960232339313 1.615396 0.0]), Molecules.PointGroup.DistanceMatrix(Any[Atom{Int16, Float64}(6, 12.011, [0.0, 1.391673, 0.0]), Atom{Int16, Float64}(6, 12.011, [1.205224, 0.695836, 0.0]), Atom{Int16, Float64}(6, 12.011, [1.205224, -0.695836, 0.0]), Atom{Int16, Float64}(6, 12.011, [0.0, -1.391673, 0.0]), Atom{Int16, Float64}(6, 12.011, [-1.205224, -0.695836, 0.0]), Atom{Int16, Float64}(6, 12.011, [-1.205224, 0.695836, 0.0]), Atom{Int16, Float64}(1, 1.008, [0.0, 2.469588, 0.0]), Atom{Int16, Float64}(1, 1.008, [2.138726, 1.234794, 0.0]), Atom{Int16, Float64}(1, 1.008, [2.138726, -1.234794, 0.0]), Atom{Int16, Float64}(1, 1.008, [0.0, -2.469588, 0.0]), Atom{Int16, Float64}(1, 1.008, [-2.138726, -1.234794, 0.0]), Atom{Int16, Float64}(1, 1.008, [-2.138726, 1.234794, 0.0])], [0.0 1.391673101250793 2.4104478246286516 2.783346 2.4104478246286516 1.391673101250793 1.077915 2.1444719451923357 3.3871046345167724 3.861261 3.3871046345167724 2.1444719451923357; 1.391673101250793 0.0 1.391672 2.4104478246286516 2.7833452025014793 2.410448 2.144472207252871 1.0779154464836285 2.1444715388421454 3.387104071319923 3.861260648984992 3.3871045641172643; 2.4104478246286516 1.391672 0.0 1.391673101250793 2.410448 2.7833452025014793 3.387104071319923 2.1444715388421454 1.0779154464836285 2.144472207252871 3.3871045641172643 3.861260648984992; 2.783346 2.4104478246286516 1.391673101250793 0.0 1.391673101250793 2.4104478246286516 3.861261 3.3871046345167724 2.1444719451923357 1.077915 2.1444719451923357 3.3871046345167724; 2.4104478246286516 2.7833452025014793 2.410448 1.391673101250793 0.0 1.391672 3.387104071319923 3.861260648984992 3.3871045641172643 2.144472207252871 1.0779154464836285 2.1444715388421454; 1.391673101250793 2.410448 2.7833452025014793 2.4104478246286516 1.391672 0.0 2.144472207252871 3.3871045641172643 3.861260648984992 3.387104071319923 2.1444715388421454 1.0779154464836285; 1.077915 2.144472207252871 3.387104071319923 3.861261 3.387104071319923 2.144472207252871 0.0 2.4695880477342778 4.277451917321806 4.939176 4.277451917321806 2.4695880477342778; 2.1444719451923357 1.0779154464836285 2.1444715388421454 3.3871046345167724 3.861260648984992 3.3871045641172643 2.4695880477342778 0.0 2.469588 4.277451917321806 4.9391760954685555 4.277452; 3.3871046345167724 2.1444715388421454 1.0779154464836285 2.1444719451923357 3.3871045641172643 3.861260648984992 4.277451917321806 2.469588 0.0 2.4695880477342778 4.277452 4.9391760954685555; 3.861261 3.387104071319923 2.144472207252871 1.077915 2.144472207252871 3.387104071319923 4.939176 4.277451917321806 2.4695880477342778 0.0 2.4695880477342778 4.277451917321806; 3.3871046345167724 3.861260648984992 3.3871045641172643 2.1444719451923357 1.0779154464836285 2.1444715388421454 4.277451917321806 4.9391760954685555 4.277452 2.4695880477342778 0.0 2.469588; 2.1444719451923357 3.3871045641172643 3.861260648984992 3.3871046345167724 2.1444715388421454 1.0779154464836285 2.4695880477342778 4.277452 4.9391760954685555 4.277451917321806 2.469588 0.0]), Molecules.PointGroup.DistanceMatrix(Any[Atom{Int16, Float64}(6, 12.011, [1.161583, -0.406755, 0.0]), Atom{Int16, Float64}(6, 12.011, [0.0, 0.552718, 0.0]), Atom{Int16, Float64}(8, 15.999, [-1.187114, -0.21286, 0.0]), Atom{Int16, Float64}(1, 1.008, [-1.932434, 0.383817, 0.0]), Atom{Int16, Float64}(1, 1.008, [2.10286, 0.13584, 0.0]), Atom{Int16, Float64}(1, 1.008, [1.122347, -1.039829, 0.881134]), Atom{Int16, Float64}(1, 1.008, [1.122347, -1.039829, -0.881134]), Atom{Int16, Float64}(1, 1.008, [0.056147, 1.193553, 0.880896]), Atom{Int16, Float64}(1, 1.008, [0.056147, 1.193553, -0.880896])], [0.0 1.5066066187356275 2.3566868414861575 3.1934221899825586 1.0864675442708818 1.08568838951515 1.08568838951515 2.135170299478709 2.135170299478709; 1.5066066187356275 0.0 1.4125683427997386 1.9398012042879549 2.143783442067785 2.1382857296381137 2.1382857296381137 1.0907794211709352 1.0907794211709352; 2.3566868414861575 1.4125683427997386 0.0 0.9547383645423493 3.3084015189628966 2.606508193625717 2.606508193625717 2.073565336203805 2.073565336203805; 3.1934221899825586 1.9398012042879549 0.9547383645423493 0.0 4.042906164996289 3.4835114492754293 3.4835114492754293 2.3207982562198293 2.3207982562198293; 1.0864675442708818 2.143783442067785 3.3084015189628966 4.042906164996289 0.0 1.7663523053700245 1.7663523053700245 2.466529679033682 2.466529679033682; 1.08568838951515 2.1382857296381137 2.606508193625717 3.4835114492754293 1.7663523053700245 0.0 1.762268 2.4748288131844594 3.0380137127445623; 1.08568838951515 2.1382857296381137 2.606508193625717 3.4835114492754293 1.7663523053700245 1.762268 0.0 3.0380137127445623 2.4748288131844594; 2.135170299478709 1.0907794211709352 2.073565336203805 2.3207982562198293 2.466529679033682 2.4748288131844594 3.0380137127445623 0.0 1.761792; 2.135170299478709 1.0907794211709352 2.073565336203805 2.3207982562198293 2.466529679033682 3.0380137127445623 2.4748288131844594 1.761792 0.0]), Molecules.PointGroup.DistanceMatrix(Any[Atom{Int16, Float64}(8, 15.999, [0.0, 0.0, 0.674493]), Atom{Int16, Float64}(6, 12.011, [0.0, 0.0, -0.529724]), Atom{Int16, Float64}(1, 1.008, [0.0, 0.934728, -1.108799]), Atom{Int16, Float64}(1, 1.008, [0.0, -0.934728, -1.108799])], [0.0 1.2042169999999999 2.0134166958799167 2.0134166958799167; 1.2042169999999999 0.0 1.0995655003723062 1.0995655003723062; 2.0134166958799167 1.0995655003723062 0.0 1.869456; 2.0134166958799167 1.0995655003723062 1.869456 0.0]), Molecules.PointGroup.DistanceMatrix(Any[Atom{Int16, Float64}(6, 12.011, [0.0, 0.550615, 0.0]), Atom{Int16, Float64}(8, 15.999, [1.177607, 0.811873, 0.0]), Atom{Int16, Float64}(8, 15.999, [-0.969424, 1.492729, 0.0]), Atom{Int16, Float64}(6, 12.011, [-0.584762, -0.856381, 0.0]), Atom{Int16, Float64}(7, 14.007, [0.400658, -1.92777, 0.0]), Atom{Int16, Float64}(1, 1.008, [-0.5071, 2.345833, 0.0]), Atom{Int16, Float64}(1, 1.008, [-1.245659, -0.945604, 0.881315]), Atom{Int16, Float64}(1, 1.008, [-1.245659, -0.945604, -0.881315]), Atom{Int16, Float64}(1, 1.008, [1.018457, -1.781229, 0.803234]), Atom{Int16, Float64}(1, 1.008, [1.018457, -1.781229, -0.803234])], [0.0 1.2062396067999923 1.351799423276989 1.5236746177120626 2.5105614991847935 1.8654645741809195 2.1370666272877408 2.1370666272877408 2.668320812784887 2.668320812784887; 1.2062396067999923 0.0 2.2524002769705476 2.4267294659019987 2.847682130443986 2.2784360771039855 3.1205223347237876 3.1205223347237876 2.7193182166969723 2.7193182166969723; 1.351799423276989 2.2524002769705476 0.0 2.380395481079562 3.6846896878468614 0.9703246445350132 2.607391363669635 2.607391363669635 3.9135222900452478 3.9135222900452478; 1.5236746177120626 2.4267294659019987 2.380395481079562 0.0 1.4556534497334865 3.2031556140843356 1.1051975920906632 1.1051975920906632 2.017632237009758 2.017632237009758; 2.5105614991847935 2.847682130443986 3.6846896878468614 1.4556534497334865 0.0 4.368948064256772 2.1099113358788326 2.1099113358788326 1.02388218454957 1.02388218454957; 1.8654645741809195 2.2784360771039855 0.9703246445350132 3.2031556140843356 4.368948064256772 0.0 3.4865087194319475 3.4865087194319475 4.472711679825673 4.472711679825673; 2.1370666272877408 3.1205223347237876 2.607391363669635 1.1051975920906632 2.1099113358788326 3.4865087194319475 0.0 1.76263 2.4146608549943407 2.9431608409127086; 2.1370666272877408 3.1205223347237876 2.607391363669635 1.1051975920906632 2.1099113358788326 3.4865087194319475 1.76263 0.0 2.9431608409127086 2.4146608549943407; 2.668320812784887 2.7193182166969723 3.9135222900452478 2.017632237009758 1.02388218454957 4.472711679825673 2.4146608549943407 2.9431608409127086 0.0 1.606468; 2.668320812784887 2.7193182166969723 3.9135222900452478 2.017632237009758 1.02388218454957 4.472711679825673 2.9431608409127086 2.4146608549943407 1.606468 0.0]), Molecules.PointGroup.DistanceMatrix(Any[Atom{Int16, Float64}(6, 12.011, [0.0, 0.0, 0.0]), Atom{Int16, Float64}(1, 1.008, [0.626891, 0.626891, 0.626891]), Atom{Int16, Float64}(1, 1.008, [-0.626891, -0.626891, 0.626891]), Atom{Int16, Float64}(1, 1.008, [-0.626891, 0.626891, -0.626891]), Atom{Int16, Float64}(1, 1.008, [0.626891, -0.626891, -0.626891])], [0.0 1.085807062807661 1.085807062807661 1.085807062807661 1.085807062807661; 1.085807062807661 0.0 1.7731155086592638 1.7731155086592638 1.7731155086592638; 1.085807062807661 1.7731155086592638 0.0 1.7731155086592638 1.7731155086592638; 1.085807062807661 1.7731155086592638 1.7731155086592638 0.0 1.7731155086592638; 1.085807062807661 1.7731155086592638 1.7731155086592638 1.7731155086592638 0.0]), Molecules.PointGroup.DistanceMatrix(Any[Atom{Int16, Float64}(6, 12.011, [0.056665, 1.080332, 0.0]), Atom{Int16, Float64}(15, 30.9737619985, [0.056665, -0.600673, 0.0]), Atom{Int16, Float64}(1, 1.008, [-0.834835, 1.695163, 0.0]), Atom{Int16, Float64}(1, 1.008, [1.00253, 1.611132, 0.0]), Atom{Int16, Float64}(1, 1.008, [-1.357659, -0.778191, 0.0])], [0.0 1.681005 1.0829540195968617 1.0846240077672076 2.335469999487255; 1.681005 0.0 2.462851028563441 2.4055647894517413 1.4254209965129598; 1.0829540195968617 2.462851028563441 0.0 1.839285554824481 2.528008098146048; 1.0846240077672076 2.4055647894517413 1.839285554824481 0.0 3.35847532580634; 2.335469999487255 1.4254209965129598 2.528008098146048 3.35847532580634 0.0]), Molecules.PointGroup.DistanceMatrix(Any[Atom{Int16, Float64}(6, 12.011, [0.984363, -0.979655, 0.0]), Atom{Int16, Float64}(6, 12.011, [0.0, 0.141435, 0.0]), Atom{Int16, Float64}(8, 15.999, [0.295443, 1.309606, 0.0]), Atom{Int16, Float64}(1, 1.008, [1.969363, -0.539961, 0.0]), Atom{Int16, Float64}(1, 1.008, [0.849207, -1.603789, 0.878633]), Atom{Int16, Float64}(1, 1.008, [0.849207, -1.603789, -0.878633]), Atom{Int16, Float64}(8, 15.999, [-1.26951, -0.279192, 0.0]), Atom{Int16, Float64}(1, 1.008, [-1.781416, 0.533546, 0.0])], [0.0 1.4919159841857719 2.390674944972863 1.078682443370615 1.0861889075943465 1.0861889075943465 2.360210142020833 3.1526672427076727; 1.4919159841857719 0.0 1.20495230257882 2.0839124584744435 2.1304824072763426 2.1304824072763426 1.3373790462053006 1.8240597581704938; 2.390674944972863 1.20495230257882 0.0 2.4945753654457903 3.0929792987360907 3.0929792987360907 2.2301024588599065 2.2171180459057656; 1.078682443370615 2.0839124584744435 2.4945753654457903 0.0 1.7772100119594758 1.7772100119594758 3.249353594407663 3.901379292748912; 1.0861889075943465 2.1304824072763426 3.0929792987360907 1.7772100119594758 0.0 1.757266 2.6486817262908358 3.5014817176508286; 1.0861889075943465 2.1304824072763426 3.0929792987360907 1.7772100119594758 1.757266 0.0 2.6486817262908358 3.5014817176508286; 2.360210142020833 1.3373790462053006 2.2301024588599065 3.249353594407663 2.6486817262908358 2.6486817262908358 0.0 0.9605159079786239; 3.1526672427076727 1.8240597581704938 2.2171180459057656 3.901379292748912 3.5014817176508286 3.5014817176508286 0.9605159079786239 0.0]), Molecules.PointGroup.DistanceMatrix(Any[Atom{Int16, Float64}(7, 14.007, [0.0, 0.0, 0.550296]), Atom{Int16, Float64}(7, 14.007, [0.0, 0.0, -0.550296])], [0.0 1.100592; 1.100592 0.0])]
SEAs = Any[Molecules.PointGroup.SEA[Molecules.PointGroup.SEA{Int64, Float64}("None", [1], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [2, 3], [0.0, 0.0, 0.0])], Molecules.PointGroup.SEA[Molecules.PointGroup.SEA{Int64, Float64}("None", [1], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [2, 3, 4], [0.0, 0.0, 0.0])], Molecules.PointGroup.SEA[Molecules.PointGroup.SEA{Int64, Float64}("None", [1, 2, 3, 4, 5, 6], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [7, 8, 9, 10, 11, 12], [0.0, 0.0, 0.0])], Molecules.PointGroup.SEA[Molecules.PointGroup.SEA{Int64, Float64}("None", [1], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [2], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [3], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [4], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [5], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [6, 7], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [8, 9], [0.0, 0.0, 0.0])], Molecules.PointGroup.SEA[Molecules.PointGroup.SEA{Int64, Float64}("None", [1], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [2], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [3, 4], [0.0, 0.0, 0.0])], Molecules.PointGroup.SEA[Molecules.PointGroup.SEA{Int64, Float64}("None", [1], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [2], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [3], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [4], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [5], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [6], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [7, 8], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [9, 10], [0.0, 0.0, 0.0])], Molecules.PointGroup.SEA[Molecules.PointGroup.SEA{Int64, Float64}("None", [1], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [2, 3, 4, 5], [0.0, 0.0, 0.0])], Molecules.PointGroup.SEA[Molecules.PointGroup.SEA{Int64, Float64}("None", [1], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [2], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [3], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [4], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [5], [0.0, 0.0, 0.0])], Molecules.PointGroup.SEA[Molecules.PointGroup.SEA{Int64, Float64}("None", [1], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [2], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [3], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [4], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [5, 6], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [7], [0.0, 0.0, 0.0]), Molecules.PointGroup.SEA{Int64, Float64}("None", [8], [0.0, 0.0, 0.0])], Molecules.PointGroup.SEA[Molecules.PointGroup.SEA{Int64, Float64}("None", [1, 2], [0.0, 0.0, 0.0])]]
COMs = Any[[1.2483188267782848, 1.8068607904101412, 0.02067611111507077], [0.0, -5.9186189890467085e-8, 0.04789390012330456], [0.0, -5.685142354124265e-18, 0.0], [-0.0541194764592242, -0.017770381536391034, 0.0], [0.0, 0.0, 0.07304868310797315], [0.012654881758962023, 0.040814231526502945, 0.0], [0.0, 0.0, 0.0], [0.02686983463465735, -0.06696346354690948, 0.0], [-0.030964198561246933, 0.05292127193099316, 0.0], [0.0, 0.0, 0.0]]
COMshifts = Any[Any[Atom{Int16, Float64}(8, 15.999, [-0.03916517197828484, -0.04044897151014126, -0.03783750831507077]), Atom{Int16, Float64}(1, 1.008, [0.9501611807217154, -0.009150727710141249, -0.00855993921507077]), Atom{Int16, Float64}(1, 1.008, [-0.3285306385782848, 0.6511577665898587, 0.6091177720849292])], Any[Atom{Int16, Float64}(7, 14.007, [0.0, 5.9186189890467085e-8, 0.06945309987669546]), Atom{Int16, Float64}(1, 1.008, [0.0, 0.9326490591861898, -0.3217029001233046]), Atom{Int16, Float64}(1, 1.008, [0.807698, -0.4663249408138101, -0.3217029001233046]), Atom{Int16, Float64}(1, 1.008, [-0.807698, -0.4663249408138101, -0.3217029001233046])], Any[Atom{Int16, Float64}(6, 12.011, [0.0, 1.391673, 0.0]), Atom{Int16, Float64}(6, 12.011, [1.205224, 0.695836, 0.0]), Atom{Int16, Float64}(6, 12.011, [1.205224, -0.695836, 0.0]), Atom{Int16, Float64}(6, 12.011, [0.0, -1.391673, 0.0]), Atom{Int16, Float64}(6, 12.011, [-1.205224, -0.695836, 0.0]), Atom{Int16, Float64}(6, 12.011, [-1.205224, 0.695836, 0.0]), Atom{Int16, Float64}(1, 1.008, [0.0, 2.469588, 0.0]), Atom{Int16, Float64}(1, 1.008, [2.138726, 1.234794, 0.0]), Atom{Int16, Float64}(1, 1.008, [2.138726, -1.234794, 0.0]), Atom{Int16, Float64}(1, 1.008, [0.0, -2.469588, 0.0]), Atom{Int16, Float64}(1, 1.008, [-2.138726, -1.234794, 0.0]), Atom{Int16, Float64}(1, 1.008, [-2.138726, 1.234794, 0.0])], Any[Atom{Int16, Float64}(6, 12.011, [1.2157024764592241, -0.38898461846360893, 0.0]), Atom{Int16, Float64}(6, 12.011, [0.0541194764592242, 0.570488381536391, 0.0]), Atom{Int16, Float64}(8, 15.999, [-1.132994523540776, -0.19508961846360895, 0.0]), Atom{Int16, Float64}(1, 1.008, [-1.8783145235407759, 0.40158738153639106, 0.0]), Atom{Int16, Float64}(1, 1.008, [2.1569794764592243, 0.15361038153639103, 0.0]), Atom{Int16, Float64}(1, 1.008, [1.176466476459224, -1.022058618463609, 0.881134]), Atom{Int16, Float64}(1, 1.008, [1.176466476459224, -1.022058618463609, -0.881134]), Atom{Int16, Float64}(1, 1.008, [0.1102664764592242, 1.211323381536391, 0.880896]), Atom{Int16, Float64}(1, 1.008, [0.1102664764592242, 1.211323381536391, -0.880896])], Any[Atom{Int16, Float64}(8, 15.999, [0.0, 0.0, 0.6014443168920268]), Atom{Int16, Float64}(6, 12.011, [0.0, 0.0, -0.6027726831079732]), Atom{Int16, Float64}(1, 1.008, [0.0, 0.934728, -1.1818476831079732]), Atom{Int16, Float64}(1, 1.008, [0.0, -0.934728, -1.1818476831079732])], Any[Atom{Int16, Float64}(6, 12.011, [-0.012654881758962023, 0.509800768473497, 0.0]), Atom{Int16, Float64}(8, 15.999, [1.164952118241038, 0.771058768473497, 0.0]), Atom{Int16, Float64}(8, 15.999, [-0.982078881758962, 1.451914768473497, 0.0]), Atom{Int16, Float64}(6, 12.011, [-0.597416881758962, -0.8971952315265029, 0.0]), Atom{Int16, Float64}(7, 14.007, [0.388003118241038, -1.968584231526503, 0.0]), Atom{Int16, Float64}(1, 1.008, [-0.519754881758962, 2.305018768473497, 0.0]), Atom{Int16, Float64}(1, 1.008, [-1.258313881758962, -0.986418231526503, 0.881315]), Atom{Int16, Float64}(1, 1.008, [-1.258313881758962, -0.986418231526503, -0.881315]), Atom{Int16, Float64}(1, 1.008, [1.005802118241038, -1.822043231526503, 0.803234]), Atom{Int16, Float64}(1, 1.008, [1.005802118241038, -1.822043231526503, -0.803234])], Any[Atom{Int16, Float64}(6, 12.011, [0.0, 0.0, 0.0]), Atom{Int16, Float64}(1, 1.008, [0.626891, 0.626891, 0.626891]), Atom{Int16, Float64}(1, 1.008, [-0.626891, -0.626891, 0.626891]), Atom{Int16, Float64}(1, 1.008, [-0.626891, 0.626891, -0.626891]), Atom{Int16, Float64}(1, 1.008, [0.626891, -0.626891, -0.626891])], Any[Atom{Int16, Float64}(6, 12.011, [0.02979516536534265, 1.1472954635469095, 0.0]), Atom{Int16, Float64}(15, 30.9737619985, [0.02979516536534265, -0.5337095364530905, 0.0]), Atom{Int16, Float64}(1, 1.008, [-0.8617048346346573, 1.7621264635469094, 0.0]), Atom{Int16, Float64}(1, 1.008, [0.9756601653653426, 1.6780954635469094, 0.0]), Atom{Int16, Float64}(1, 1.008, [-1.3845288346346574, -0.7112275364530904, 0.0])], Any[Atom{Int16, Float64}(6, 12.011, [1.0153271985612469, -1.0325762719309932, 0.0]), Atom{Int16, Float64}(6, 12.011, [0.030964198561246933, 0.08851372806900684, 0.0]), Atom{Int16, Float64}(8, 15.999, [0.32640719856124695, 1.2566847280690068, 0.0]), Atom{Int16, Float64}(1, 1.008, [2.000327198561247, -0.5928822719309932, 0.0]), Atom{Int16, Float64}(1, 1.008, [0.8801711985612469, -1.6567102719309932, 0.878633]), Atom{Int16, Float64}(1, 1.008, [0.8801711985612469, -1.6567102719309932, -0.878633]), Atom{Int16, Float64}(8, 15.999, [-1.238545801438753, -0.33211327193099316, 0.0]), Atom{Int16, Float64}(1, 1.008, [-1.7504518014387531, 0.4806247280690068, 0.0])], Any[Atom{Int16, Float64}(7, 14.007, [0.0, 0.0, 0.550296]), Atom{Int16, Float64}(7, 14.007, [0.0, 0.0, -0.550296])]]
MOITs = Any[[0.8506310823932961 0.1990554320457495 0.18620403155090115; 0.1990554320457495 1.4403375118431023 -0.42437006076117806; 0.18620403155090115 -0.42437006076117806 1.4970247630712796], [1.6957183090708572 0.0 0.0; 0.0 1.6957182586275732 -3.2427652335575985e-7; 0.0 -3.2427652335575985e-7 2.630380321153812], [88.2299814706577 0.0 0.0; 0.0 88.22999596081817 0.0; 0.0 0.0 176.45997743147586], [14.715300174100918 4.353853093951238 0.0; 4.353853093951238 52.51463093256641 0.0; 0.0 0.0 60.97075555513861], [14.728706546062075 0.0 0.0; 0.0 12.967294215150332 0.0; 0.0 0.0 1.761412330911744], [127.18693509903127 15.180098551948333 0.0; 15.180098551948333 51.91095466126222 0.0; 0.0 0.0 173.36478497675412], [3.1690901159043836 0.0 0.0; 0.0 3.1690901159043836 0.0; 0.0 0.0 3.1690901159043836], [31.11102406202894 -1.0303999752233461 0.0; -1.0303999752233461 3.678418822081161 0.0; 0.0 0.0 34.7894428841101], [47.608404223837276 4.399008427594539 0.0; 4.399008427594539 48.88054006235208 0.0; 0.0 0.0 93.3762566210753], [8.483358812874624 0.0 0.0; 0.0 8.483358812874624 0.0; 0.0 0.0 0.0]]
EMOITs = Any[([0.6578934196149924, 1.236103259038848, 1.8939966786538391], [-0.81649546927384 0.5773518413024179 3.045108609711633e-11; 0.4216331768098601 0.5962769215475041 -0.6831392954903539; 0.3944117301355493 0.5577801396379573 0.7302880958614506]), ([1.695718258627461, 1.6957183090708572, 2.6303803211539245], [0.0 1.0 0.0; 0.9999999999999398 0.0 -3.4694520763925603e-7; 3.4694520763925603e-7 0.0 0.9999999999999398]), ([88.2299814706577, 88.22999596081817, 176.45997743147586], [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]), ([14.220291330478627, 53.009639776188706, 60.97075555513861], [0.9935987805836436 0.1129666465054909 0.0; -0.11296664650549092 0.9935987805836436 0.0; 0.0 0.0 1.0]), ([1.761412330911744, 12.967294215150332, 14.728706546062075], [0.0 0.0 1.0; 0.0 1.0 0.0; 1.0 0.0 0.0]), ([48.96503551000708, 130.1328542502864, 173.36478497675412], [-0.19051030587628406 0.9816851956482407 0.0; 0.9816851956482407 0.19051030587628406 0.0; 0.0 0.0 1.0]), ([3.1690901159043836, 3.1690901159043836, 3.1690901159043836], [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]), ([3.639770272689056, 31.149672611421042, 34.7894428841101], [0.03748194117262715 0.9992973051529418 0.0; 0.9992973051529418 -0.03748194117262715 0.0; 0.0 0.0 1.0]), ([43.79971595811038, 52.68922832807898, 93.3762566210753], [0.756010993372136 0.654558918586002 0.0; -0.654558918586002 0.756010993372136 0.0; 0.0 0.0 1.0]), ([0.0, 8.483358812874624, 8.483358812874624], [0.0 0.0 1.0; 0.0 1.0 0.0; 1.0 0.0 0.0])]
tol = 1E-5

@testset "SymTools" begin
    for i = eachindex(mols)
        path = joinpath(@__DIR__, "xyz/"*mols[i]*".xyz")
        mol = Molecules.parse_file(path)

        D = buildD(mol)
        @test isapprox(D.distances, Ds[i].distances, rtol=tol)

        sea = findSEA(D, 5)
        @test isapprox([sea[j].set for j = 1:size(sea)[1]], [SEAs[i][j].set for j = 1:size(SEAs[i])[1]], rtol=tol)
        moit = calcmoit(mol)
        @test isapprox(moit, MOITs[i], rtol=tol)

        emoit = eigenmoit(moit)
        @test isapprox(emoit[1], EMOITs[i][1], rtol=tol)
        @test isapprox(emoit[2], EMOITs[i][2], rtol=tol)

    end
end


text:
module YAJL

export @yajl

using Libdl

# Thrown when an invalid callback is registered.
const invalid_sig = ArgumentError("Invalid callback type signature")

"""
Base type for YAJL contexts.
To implement a custom `Context`'s behaviour, see [`@yajl`](@ref) and [`collect`](@ref).
For a full example, see `minifier.jl`.
"""
abstract type Context end

# YAJL callback types.
const CALLBACKS = (:null, :boolean, :integer, :double, :number, :string, :map_start,
                   :map_key, :map_end, :array_start, :array_end)

# Default callbacks.
for s in [:null, :boolean, :integer, :double, :number, :string, :map_start, :map_key,
          :map_end, :array_start, :array_end]
    f = Symbol(:cb_, s)
    @eval $f(::Context) = C_NULL
end

"""
    collect(ctx::Context)

Override this function for your custom [`Context`](@ref) to specify what is returned from [`run`](@ref).
By default, `ctx` itself is returned.
"""
collect(ctx::Context) = ctx

struct Callbacks
    null::Ptr{Cvoid}
    boolean::Ptr{Cvoid}
    integer::Ptr{Cvoid}
    double::Ptr{Cvoid}
    number::Ptr{Cvoid}
    string::Ptr{Cvoid}
    map_start::Ptr{Cvoid}
    map_key::Ptr{Cvoid}
    map_end::Ptr{Cvoid}
    array_start::Ptr{Cvoid}
    array_end::Ptr{Cvoid}

    Callbacks(ctx::Context) = new(
        cb_null(ctx),
        cb_boolean(ctx),
        cb_integer(ctx),
        cb_double(ctx),
        cb_number(ctx),
        cb_string(ctx),
        cb_map_start(ctx),
        cb_map_key(ctx),
        cb_map_end(ctx),
        cb_array_start(ctx),
        cb_array_end(ctx),
    )
end

# Parse options.
const ALLOW_COMMENTS = 0x01
const DONT_VALIDATE_STRINGS = 0x02
const ALLOW_TRAILING_GARBAGE = 0x04
const ALLOW_MULTIPLE_VALUES = 0x08
const ALLOW_PARTIAL_VALUES = 0x10
const OPTIONS = [
    ALLOW_COMMENTS, DONT_VALIDATE_STRINGS, ALLOW_TRAILING_GARBAGE, ALLOW_MULTIPLE_VALUES,
    ALLOW_PARTIAL_VALUES,
]

# Check if a context type has defined a function.
# We do this specific check because cb_* is implemented for all Context types.
hasmeth(T::Type{<:Context}, f::Function) = any(methods(f)) do m
    arg1 = m.sig isa UnionAll ? m.sig.body.types[2] : m.sig.types[2]
    arg1 !== Context && T <: arg1
end

# Check that a callback's type signature is valid.
function checktypes(f::Symbol, Ts::Type...)
    # Drop the leading on_.
    f = Symbol(string(f)[4:end])
    err = invalid_sig
    if f in (:boolean, :integer)
        length(Ts) == 1 && Ts[1] <: Integer || throw(err)
    elseif f === :double
        length(Ts) == 1 && Ts[1] <: AbstractFloat || throw(err)
    elseif f in (:number, :string, :map_key)
        length(Ts) == 2 && Ts[1] in (Ptr{UInt8}, Cstring) && Ts[2] <: Integer || throw(err)
    else
        isempty(Ts) || throw(err)
    end
end

"""
Register a callback for a specific data type.
If the callback function returns a value via the `return` keyword,
then it must be an `Integer` value.
If the value is `false` or `0`, any further processing is cancelled.

The callbacks to be overridden are as follows, where `T` is a `Context` subtype:

- `null(ctx::T)`: Called on `null` values.
- `boolean(ctx::T, v::Bool)`: Called on boolean values.
- `integer(ctx::T, v::Int)`: Called on integer values (see note below).
- `double(ctx::T, v::Float64)`: Called on float values (see note below).
- `number(ctx::T, v::Ptr{UInt8}, len::Int)`: Called on numeric values (see note below).
- `string(ctx::T, v::Ptr{UInt8}, len::Int)`: Called on string values.
- `map_start(ctx::T)`: Called when an object begins (`{`).
- `map_key(ctx::T, v::Ptr{UInt8}, len::Int)`: Called on object keys.
- `map_end(ctx::T)`: Called when an object ends (`}`).
- `array_start(ctx::T)`: Called when an array begins (`[`).
- `array_end(ctx::T)`: Called when an array ends (`]`).

For string arguments which appear as `Ptr{UInt8}`, `Cstring` can also be used.
However, `Ptr{UInt8}` is usually better if you want to use `unsafe_string(v, len)`.

!!! note
    To handle numbers, implement either `number` or both `integer` and `double`.
    Usually, `number` is a better choice because `integer` and `double have limited precision.
    See [here](https://lloyd.github.io/yajl/yajl-2.1.0/structyajl__callbacks.html) for more details.

For a full example, see `minifier.jl`.
"""
macro yajl(ex::Expr)
    # We start with an expression that should be a function definition. One of:
    # - f(...) = ...
    # - function f(...) ... end
    ex.head in (:(=), :function) ||
        throw(ArgumentError("Expression must be a function definition"))

    # If the context type for which the callback is being implemented is parametric,
    # then the function definition might have a "where" clause that changes the structure.
    # The value is either:
    # - A list of symbols or expressions: f(...) where {T, U <: IO, ...} = ...
    # - nothing: f(...) = ...
    where = ex.args[1].head === :where ? ex.args[1].args[2:end] : nothing

    # Next, we extract the function signature, whose position can change depending on
    # whether a where clause is present, and whether a type conversion is present.
    # It looks like:
    # - f(...)
    sig = where === nothing ? ex.args[1] : ex.args[1].args[1]
    sig.head === :(::) && (sig = sig.args[1])

    # Now we make sure that the function returns Cint.
    # TODO: There must be a more elegant way...
    push!(ex.args[2].args, true)
    if where === nothing && ex.args[1].head === :(::)
        ex.args[1].args[end] = :Cint
    elseif where === nothing
        ex.args[1] = Expr(:(::), ex.args[1], :Cint)
    elseif ex.args[1].args[1].head === :(::)
        ex.args[1].args[1].args[end] = :Cint
    else
        ex.args[1].args[1] = Expr(:(::), ex.args[1].args[1], :Cint)
    end

    # Unmodified function name. This should be :null, :number, :string, etc.
    f = sig.args[1]
    f in CALLBACKS || throw(ArgumentError("Invalid callback name"))

    # Name of the cb_* function, which returns a function pointer.
    cb = Expr(:., :YAJL, QuoteNode(Symbol(:cb_, f)))

    # Rename the function to on_* to avoid any Base conflicts.
    f = sig.args[1] = Symbol(:on_, f)

    # Argument expressions.
    # There must be at least one, and each one should be typed with an optional name.
    # For example:
    # - foo::T
    # - ::T
    args = sig.args[2:end]
    isempty(args) && throw(invalid_sig)
    any(ex -> ex isa Symbol, args) &&
        throw(ArgumentError("Callback arguments must be typed"))

    # Argument types, which will be used for @cfunction.
    # Using the last expression argument gets the type for both named/unnamed arguments.
    # If all of the types were simple symbols (e.g. :Int), then Ts has eltype Symbol.
    # However, we're going to insert an Expr so we need to do a conversion.
    Ts = convert(Vector{Any}, map(ex -> ex.args[end], args))

    # Change the first argument to be a reference, since that's what @cfunction needs.
    # In this case, we need no type parameters at all because @cfunction doesn't like them.
    T = Ts[1]
    Tnp = T isa Expr && T.head === :curly ? T.args[1] : T
    Ts[1] = :(Ref{$Tnp})

    # Sometimes we also need the whole type with where.
    Twh = where === nothing ? T : Expr(:where, T, where...)

    # We'll check the context type manually, but we still need the later ones to
    # make sure that the signature is correct.
    tocheck = map(esc, Ts[2:end])

    # Create the callback function, which will create a C function pointer.
    cbfun = Expr(:call, cb, :(::$T))
    where === nothing || (cbfun = Expr(:where, cbfun, where...))
    cbfun = Expr(:(=), cbfun, :(@cfunction $f Cint ($(Ts...),)))

    # Warn if a useless or destructive callback is being added.
    # For more info, see the note in the docstring.
    warn1 = if f === :on_number
        quote
            (YAJL.hasmeth($Twh, YAJL.cb_integer) || YAJL.hasmeth($Twh, YAJL.cb_double)) &&
                @warn "Implementing number callback for $($T) disables both integer and double callbacks"
        end
    end
    warn2 = if f in (:on_integer, :on_double)
        quote
            YAJL.hasmeth($Twh, YAJL.cb_number) &&
                @warn "Implementing integer or double callback for $($T) has no effect because number callback is already implemented"
        end
    end

    quote
        $(esc(Twh)) <: Context || throw(invalid_sig)
        checktypes($(QuoteNode(f)), $(tocheck...))
        $(esc(warn1))
        $(esc(warn2))
        $(esc(ex))
        $(esc(cbfun))
    end
end

const ST_OK = 0
const ST_CLIENT_CANCELLED = 1
const ST_ERROR = 2

# A YAJL parser error.
struct ParseError <: Exception
    reason::String
end

# Check the parser status and throw an exception if there's an error.
function checkstatus(handle::Ptr{Cvoid}, status::Cint, text::Vector{UInt8}, len::Int)
    return if status == ST_OK
        true
    elseif status == ST_CLIENT_CANCELLED
        false
    elseif status == ST_ERROR
        err = ccall(yajl[:get_error], Cstring, (Ptr{Cvoid}, Cint, Ptr{Cuchar}, Csize_t),
                    handle, 1, text, len)
        reason = unsafe_string(err)
        ccall(yajl[:free_error], Cvoid, (Ptr{Cvoid}, Cstring), handle, err)
        throw(ParseError(reason))
    else
        @warn "yajl_parse returned unknown status: $status"
        true
    end
end

"""
    run(io::IO, ctx::Context; chunk::Integer=2^16, options::Integer=0x0)

Parse the JSON data from `io` and process it with `ctx`'s callbacks.
The return value is determined by the implementation of [`collect`](@ref) for `ctx`.

## Keywords
- `chunk::Integer=2^16`: Number of bytes to read from `io` at a time.
- `options::Integer=0x0`: YAJL parser options, ORed together.
"""
function run(io::IO, ctx::T; chunk::Integer=2^16, options::Integer=0x0) where T <: Context
    cbs = Callbacks(ctx)
    handle = ccall(yajl[:alloc], Ptr{Cvoid}, (Ref{Callbacks}, Ptr{Cvoid}, Ref{T}),
                   cbs, C_NULL, ctx)

    for o in OPTIONS
        if options & o == o
            ccall(yajl[:config], Cint, (Ptr{Cvoid}, Cuint), handle, o)
        end
    end

    cancelled = false
    text = Vector{UInt8}(undef, chunk)
    while !eof(io)
        n = readbytes!(io, text)
        status = ccall(yajl[:parse], Cint, (Ptr{Cvoid}, Ptr{Cuchar}, Csize_t),
                       handle, text, n)
        if !checkstatus(handle, status, text, n)
            cancelled = true
            break
        end
    end

    if !cancelled
        status = ccall(yajl[:complete_parse], Cint, (Ptr{Cvoid},), handle)
        checkstatus(handle, status, UInt8[], 0)
    end

    ccall(yajl[:free], Cvoid, (Ptr{Cvoid},), handle)

    return collect(ctx)
end

# Container for function pointers.
const yajl = Dict{Symbol, Ptr{Cvoid}}()

const depsfile = joinpath(dirname(@__DIR__), "deps", "deps.jl")
isfile(depsfile) ? include(depsfile) : error("""Run Pkg.build("YAJL")""")

# Load functions at runtime.
function __init__()
    check_deps()
    lib = Libdl.dlopen(libyajl)
    for f in [:alloc, :complete_parse, :config, :free, :free_error, :get_error, :parse]
        yajl[f] = Libdl.dlsym(lib, Symbol(:yajl_, f))
    end
end

include("minifier.jl")

end


text:
module Rocket

using Unrolled

include("utils.jl")
include("helpers/list.jl")
include("helpers/mstorage.jl")
include("helpers/updates.jl")

include("teardown.jl")
include("teardown/void.jl")

include("scheduler.jl")
include("actor.jl")
include("subscribable.jl")
include("operator.jl")
include("subject.jl")

include("schedulers/asap.jl")
include("schedulers/async.jl")
include("schedulers/threads.jl")
include("schedulers/pending.jl")
include("schedulers/postponed.jl")

include("actor/function.jl")
include("actor/lambda.jl")
include("actor/logger.jl")
include("actor/void.jl")
include("actor/sync.jl")
include("actor/keep.jl")
include("actor/buffer.jl")
include("actor/server.jl")
include("actor/storage.jl")
include("actor/test.jl")

include("subjects/subject.jl")
include("subjects/behavior.jl")
include("subjects/replay.jl")
include("subjects/pending.jl")
include("subjects/recent.jl")

@generate_subscribe! Subject AbstractSubject
@generate_subscribe! BehaviorSubjectInstance AbstractSubject
@generate_subscribe! ReplaySubjectInstance AbstractSubject
@generate_subscribe! PendingSubjectInstance AbstractSubject
@generate_subscribe! RecentSubjectInstance AbstractSubject

include("observable/generate.jl")
include("observable/single.jl")
include("observable/array.jl")
include("observable/iterable.jl")
include("observable/faulted.jl")
include("observable/never.jl")
include("observable/completed.jl")
include("observable/proxy.jl")
include("observable/timer.jl")
include("observable/interval.jl")
include("observable/function.jl")
include("observable/file.jl")
include("observable/network.jl")
include("observable/combined.jl")
include("observable/combined_updates.jl")
include("observable/collected.jl")
include("observable/race.jl")
include("observable/merged.jl")
include("observable/concat.jl")
include("observable/lazy.jl")
include("observable/connectable.jl")
include("observable/scheduled.jl")
include("observable/defer.jl")
include("observable/zipped.jl")

include("operators/map.jl")
include("operators/map_to.jl")
include("operators/reduce.jl")
include("operators/scan.jl")
include("operators/filter.jl")
include("operators/filter_type.jl")
include("operators/find.jl")
include("operators/find_index.jl")
include("operators/some.jl")
include("operators/count.jl")
include("operators/enumerate.jl")
include("operators/take.jl")
include("operators/take_until.jl")
include("operators/first.jl")
include("operators/last.jl")
include("operators/tap.jl")
include("operators/tap_on_subscribe.jl")
include("operators/tap_on_unsubscribe.jl")
include("operators/tap_on_complete.jl")
include("operators/sum.jl")
include("operators/max.jl")
include("operators/min.jl")
include("operators/delay.jl")
include("operators/uppercase.jl")
include("operators/lowercase.jl")
include("operators/to_array.jl")
include("operators/tuple_with.jl")
include("operators/switch_map.jl")
include("operators/switch_map_to.jl")
include("operators/merge_map.jl")
include("operators/concat_map.jl")
include("operators/concat_map_to.jl")
include("operators/exhaust_map.jl")
include("operators/multicast.jl")
include("operators/ref_count.jl")
include("operators/publish.jl")
include("operators/share.jl")
include("operators/catch_error.jl")
include("operators/rerun.jl")
include("operators/safe.jl")
include("operators/noop.jl")
include("operators/default_if_empty.jl")
include("operators/error_if_empty.jl")
include("operators/error_if.jl")
include("operators/error_if_not.jl")
include("operators/debounce_time.jl")
include("operators/skip_next.jl")
include("operators/skip_error.jl")
include("operators/skip_complete.jl")
include("operators/schedule_on.jl")
include("operators/async.jl")
include("operators/parallel.jl")
include("operators/discontinue.jl")
include("operators/with_latest.jl")
include("operators/ignore.jl")
include("operators/start_with.jl")
include("operators/accumulated.jl")
include("operators/pairwise.jl")
include("operators/limit_subscribers.jl")
include("operators/substitute.jl")
include("operators/override.jl")

include("extensions/observable/single.jl")

end # module