text:
# Autogenerated wrapper script for LevelDB_jll for x86_64-apple-darwin
export leveldb

JLLWrappers.@generate_wrapper_header("LevelDB")
JLLWrappers.@declare_library_product(leveldb, "@rpath/libleveldb.1.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        leveldb,
        "lib/libleveldb.1.23.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@hl struct Person
    name::String
end

@hl struct Developer <: Person
    salary::Int
end

@hl struct SpecializedDeveloper <: Developer
    language::String
end

function sumsalaries(first::Developer, second::Developer)
    return first.salary + second.salary
end

@hl struct Job
    nb_hours::Int
    assigned_dev::Developer
end

@concretify @hl struct ConcreteJob
    nb_hours::Int
    assigned_dev::Developer
end

ext: jl


text:
# UT

function datetime2gpsws(dt::DateTime)
    d = diff2UTC(dt, tstGPS)
    jd = datetime2julian(dt - d)
    jd -= GPS_EPOCH_JD
    w = floor(Int64, jd) ÷ 7
    s = round(Int64, (jd - w * 7) * 86400)
    return w, s
end

function gpsws2datetime(w::Int64, s::Union{Int64,Float64})
    dow =  s ÷ 86400
    jd = w * 7 + dow + GPS_EPOCH_JD +
            (s - dow * 86400) / 86400.0
    dt = julian2datetime(jd)
    dt + diff2UTC(dt, tstGPS)
end


text:
export sources

"""
`sources(f::DiscreteFunction)` returns a list of all inputs to `f`
that are not the output of `f`. That is, this is a list of the vertices
with out-degree 0 in the graph of `f`.
"""
function sources(f::DiscreteFunction)
    n = length(f)
    S = setdiff( Set(1:n), image(f))
    return sort(collect(S))
end


text:
export Entity, update!
export ECS, XYZ, accumulate_XYZ
export getEntityRow, getEntityById, getEntityRowById, getDfRowProp, setDfRowProp!
export ECSIteratorState, Level
export instantiate!, destroy!
export Scene, scene_view

abstract type Entity end

update!(e, Δ) = nothing

# magic symbols that getproperty and setproperty! (and end users) care about
const ENT = :ent
const TYPE = :type
const ID = :id
const PARENT = :parent
const CHILDREN = :children
const POSITION = :pos
const ROTATION = :rot
const ABSOLUTE_POSITION = :abs_pos
const ABSOLUTE_ROTATION = :abs_rot
const HIDDEN = :hidden
const PROPS = :props

# TODO: replace XYZ with something more elegant

# position, rotation, velocity, acceleration, whatever
mutable struct XYZ
  x::Number
  y::Number
  z::Number
  XYZ(x=0, y=0, z=0) = new(x, y, z)
end

import Base.+, Base.-, Base.*, Base.÷, Base./, Base.%
+(a::XYZ, b::XYZ) = XYZ(a.x+b.x,a.y+b.y,a.z+b.z)
-(a::XYZ, b::XYZ) = XYZ(a.x-b.x,a.y-b.y,a.z-b.z)
*(a::XYZ, b::Number) = XYZ(a.x*b,a.y*b,a.z*b)
÷(a::XYZ, b::Number) = XYZ(a.x÷b,a.y÷b,a.z÷b)
/(a::XYZ, b::Number) = XYZ(a.x/b,a.y/b,a.z/b)
%(a::XYZ, b::Number) = XYZ(a.x%b,a.y%b,a.z%b)

# columns of the in-memory database
const components = Dict(
  ENT=>Entity,
  TYPE=>DataType,
  ID=>Number,
  PARENT=>Number,
  CHILDREN=>Set{Number},
  POSITION=>XYZ,
  ROTATION=>XYZ,
  HIDDEN=>Bool,
  PROPS=>Dict{Symbol, Any}
)

mutable struct ECS
  df::DataFrame
  awoken::Bool
  lock::ReentrantLock
  next_id::Number
  function ECS()
    df = DataFrame(
      NamedTuple{Tuple(keys(components))}(
        t[] for t in values(components)
      ))
    return new(df, false, ReentrantLock(), 0)
  end
end

# internally using Base.getproperty directly so
# as to not break if the symbol values change
getEntityRow(ent::Entity) = @view ecs().df[getproperty(ecs().df, ENT) .== [ent], :]
function getEntityById(id::Number)
  arr = ecs().df[getproperty(ecs().df, ID) .== [id], ENT]
  if length(arr) > 0
    return arr[1]
  end
  return nothing
end
getEntityRowById(id::Number) = @view ecs().df[getproperty(ecs().df, ID) .== [id], :]
getDfRowProp(r, s) = r[!, s][1]
setDfRowProp!(r, s, x) = r[!, s][1] = x

function Base.propertynames(ent::Entity)
  return (
    keys(components)...,
    ABSOLUTE_POSITION,
    ABSOLUTE_ROTATION,
    [n for n in keys(getproperty(ent, PROPS))]...
  )
end

function Base.hasproperty(ent::Entity, s::Symbol)
  return s in Base.propertynames(ent)
end

function accumulate_XYZ(r, s)
  acc = XYZ()
  while true
    inc = getDfRowProp(r, s)
    acc += inc
    r = getEntityRowById(getDfRowProp(r, PARENT))
    getDfRowProp(r, PARENT) != 0 || return acc
  end
end

function Base.getproperty(ent::Entity, s::Symbol)
  e = getEntityRow(ent)
  if s == ABSOLUTE_POSITION return accumulate_XYZ(e, POSITION)
  elseif s == ABSOLUTE_ROTATION return accumulate_XYZ(e, ROTATION)
  elseif s in keys(components) return getDfRowProp(e, s)
  elseif s in keys(getDfRowProp(e, PROPS)) return getDfRowProp(e, PROPS)[s]
  else return getfield(ent, s)
  end
end

function Base.setproperty!(ent::Entity, s::Symbol, x)
  e = getEntityRow(ent)
  if s in [
    ENT, # immutable
    TYPE, # automatically set
    ID, # automatically set
    ABSOLUTE_POSITION, # computed
    ABSOLUTE_ROTATION # computed
  ]
    error("cannot set property $(s) on Entity")
  end

  lock(ecs().lock)
  if s == PARENT
    par = getEntityById(getDfRowProp(e, PARENT))
    push!(getproperty(par, CHILDREN), getDfRowProp(e, ID))
  elseif s in keys(components) && s != PROPS
    setDfRowProp!(e, s, x)
  else
    getDfRowProp(e, PROPS)[s] = x
  end
  unlock(ecs().lock)
end

Base.length(e::ECS) = size(e.df)[1]

# TODO: the type/struct approach feels clunky, fix with Vals?
abstract type ECSIterator end

mutable struct ECSIteratorState
  root::Number
  q::Queue{Number}
  root_visited::Bool
  index::Number
  ECSIteratorState(; root=0, q=Queue{Number}(),
  root_visited=false, index=1) = new(root, q, root_visited, index)
end

# refers to tree level, i.e. breadth-first,
# nothing special to see here
struct Level end
Base.length(l::Level) = length(ecs())

function Base.iterate(l::Level, state::ECSIteratorState=ECSIteratorState())
  if isempty(state.q)
    if !state.root_visited # just started
      enqueue!(state.q, state.root)
      state.root_visited = true
    else # just finished
      return nothing
    end
  end

  ent = getEntityById(dequeue!(state.q))

  for c in getproperty(ent, CHILDREN)
    enqueue!(state.q, c)
  end

  return (ent, state)
end

function handleMessage!(e::ECS, m::TICK)
  Log.@debug "ECS tick"
  try
    map((ent) -> update!(ent, m.Δ), Level()) # TODO investigate parallelization
  catch
    handleException()
  end
end

function awake!(e::ECS)
  Log.@debug "ECS awake!"
  e.awoken = true
  map(awake!, Level())
  listenFor(e, TICK)
end

function shutdown!(e::ECS)
  Log.@debug "ECS shutdown!"
  unlistenFrom(e, TICK)
  map(shutdown!, Level())
  e.awoken = false
end

function instantiate!(e::Entity; kw...)
  lock(ecs().lock)

  id = ecs().next_id
  ecs().next_id += 1

  # update ecs
  # allows invalid parents and children for now
  push!(ecs().df, Dict(
    ENT=>e,
    TYPE=>typeof(e),
    ID=>id,
    CHILDREN=>get(kw, CHILDREN, Set{Number}()),
    PARENT=>get(kw, PARENT, 0),
    POSITION=>get(kw, POSITION, XYZ()),
    ROTATION=>get(kw, ROTATION, XYZ()),
    HIDDEN=>get(kw, HIDDEN, false),
    PROPS=>merge(get(kw, PROPS, Dict{Symbol, Any}()),
      Dict(k=>v for (k,v) in kw if k ∉
      [CHILDREN, PARENT, POSITION, ROTATION, HIDDEN, PROPS]))
  ))

  if id != 0 # root has no parent but itself
    par = getEntityById(get(kw, PARENT, 0))
    push!(getproperty(par, CHILDREN), id)
  end

  unlock(ecs().lock)

  if ecs().awoken awake!(e) end

  return e
end

function destroy!(e::Entity)
  shutdown!(e)

  lock(ecs().lock)

  p = getEntityById(getproperty(e, PARENT))

  # if not root
  if getproperty(p, ID) != getproperty(e, ID)
    # update parent
    delete!(getproperty(p, CHILDREN), getproperty(e, ID))
    # update dataframe
    deleteat!(ecs().df, getproperty(ecs().df, ENT) .== [e])
  end

  unlock(ecs().lock)
end

function destroy!(es...)
  map(destroy!, es) # TODO investigate parallelization
end

# this is the scene graph, ladies and gentlemen,
# which we can traverse and mutate however we want
# during iteration, and initialize however we want
# and destroy however we want...sorry, it took me
# a long time to come up with this design, and i'm
# a little bit psyched about it. :)
scene_view() = ecs().df[(ecs().df.type .<: [Renderable]) .& (ecs().df.hidden .== [false]), :]

struct Scene end

Base.length(s::Scene) = size(scene_view())[1]

function Base.iterate(s::Scene, state::ECSIteratorState=ECSIteratorState())
  # does reverse-z order for now, only
  # suitable for simple 2d drawing
  if state.index > length(s) return nothing end
  ent = scene_view()[!, ENT][state.index]
  state.index += 1
  return (ent, state)
end

function awake!(s::Scene)
  Log.@debug "Scene awake!"
  listenFor(s, TICK)
end

function shutdown!(s::Scene)
  Log.@debug "Scene shutdown!"
  unlistenFrom(s, TICK)
end

function handleMessage!(s::Scene, m::TICK)
  # sort just once per tick rather than every time we iterate
  Log.@debug "Scene tick"
  try
    sort!(ecs().df, [order(POSITION, rev=true, by=(pos)->pos.z)])
  catch
    handleException()
  end
end


text:
using FieldTracer, PyPlot
using Documenter

makedocs(;
    modules=[FieldTracer],
    authors="Hongyang Zhou <hyzhou@umich.edu>",
    repo="https://github.com/henry2004y/FieldTracer.jl/blob/{commit}{path}#L{line}",
    sitename="FieldTracer.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://henry2004y.github.io/FieldTracer.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Example" => "example.md",
        "Log" => "log.md",
        "Internal" => "internal.md",
    ],
)

deploydocs(;
    repo="github.com/henry2004y/FieldTracer.jl",
)


text:
using WAV
using SpeechAugment

dB = (5.0, 15.0)
period = 10
# read a wav file
wav = wavread("xxx.wav")[1]
# init a addnoise function with dir only having noise audios
addnoise = initAddNoise("XXPathFullOfNoiseWAVs", period, dB)
# every period times the addnoise function would change a background noise audio
noisy = addnoise(wav);

max_stars_repo_name: sonosole/SpeechAugment

max_issues_repo_name: sonosole/SpeechAugment

max_forks_repo_name: sonosole/SpeechAugment


text:
##
# Generic base overloads
Base.extrema(primitive::GeometryPrimitive) = (minimum(primitive), maximum(primitive))
function widths(x::AbstractRange)
    mini, maxi = Float32.(extrema(x))
    return maxi - mini
end

##
# conversion & decompose
convert_simplex(::Type{T}, x::T) where {T} = (x,)

function convert_simplex(NFT::Type{NgonFace{N,T1}},
                         f::Union{NgonFace{N,T2}}) where {T1,T2,N}
    return (convert(NFT, f),)
end

function convert_simplex(NFT::Type{NgonFace{3,T}}, f::NgonFace{3,T2}) where {T,T2}
    return (convert(NFT, f),)
end
function convert_simplex(NFT::Type{NgonFace{2,T}}, f::NgonFace{2,T2}) where {T,T2}
    return (convert(NFT, f),)
end

"""
    convert_simplex(::Type{Face{3}}, f::Face{N})

Triangulate an N-Face into a tuple of triangular faces.
"""
@generated function convert_simplex(::Type{TriangleFace{T}},
                                    f::Union{SimplexFace{N},NgonFace{N}}) where {T,N}
    3 <= N || error("decompose not implemented for N <= 3 yet. N: $N")# other wise degenerate
    v = Expr(:tuple)
    for i in 3:N
        push!(v.args, :(TriangleFace{T}(f[1], f[$(i - 1)], f[$i])))
    end
    return v
end

"""
    convert_simplex(::Type{Face{2}}, f::Face{N})

Extract all line segments in a Face.
"""
@generated function convert_simplex(::Type{LineFace{T}},
                                    f::Union{SimplexFace{N},NgonFace{N}}) where {T,N}
    2 <= N || error("decompose not implented for N <= 2 yet. N: $N")# other wise degenerate

    v = Expr(:tuple)
    for i in 1:(N - 1)
        push!(v.args, :(LineFace{$T}(f[$i], f[$(i + 1)])))
    end
    # connect vertices N and 1
    push!(v.args, :(LineFace{$T}(f[$N], f[1])))
    return v
end

to_pointn(::Type{T}, x) where {T<:Point} = convert_simplex(T, x)[1]

# disambiguation method overlords
convert_simplex(::Type{Point}, x::Point) = (x,)
convert_simplex(::Type{Point{N,T}}, p::Point{N,T}) where {N,T} = (p,)
function convert_simplex(::Type{Point{N,T}}, x) where {N,T}
    N2 = length(x)
    return (Point{N,T}(ntuple(i -> i <= N2 ? T(x[i]) : T(0), N)),)
end

function convert_simplex(::Type{Vec{N,T}}, x) where {N,T}
    N2 = length(x)
    return (Vec{N,T}(ntuple(i -> i <= N2 ? T(x[i]) : T(0), N)),)
end

collect_with_eltype(::Type{T}, vec::Vector{T}) where {T} = vec
collect_with_eltype(::Type{T}, vec::AbstractVector{T}) where {T} = collect(vec)

function collect_with_eltype(::Type{T}, iter) where {T}
    # TODO we could be super smart about allocating the right length
    # but its kinda annoying, since e.g. T == Triangle and first(iter) isa Quad
    # will need double the length etc - but could all be figured out ;)
    result = T[]
    for element in iter
        # convert_simplex always returns a tuple,
        # so that e.g. convert(Triangle, quad) can return 2 elements
        for telement in convert_simplex(T, element)
            push!(result, telement)
        end
    end
    return result
end

"""
The unnormalized normal of three vertices.
"""
function orthogonal_vector(v1, v2, v3)
    a = v2 - v1
    b = v3 - v1
    return cross(a, b)
end

"""
```
normals{VT,FD,FT,FO}(vertices::Vector{Point{3, VT}},
                    faces::Vector{Face{FD,FT,FO}},
                    NT = Normal{3, VT})
```
Compute all vertex normals.
"""
function normals(vertices::AbstractVector{<:AbstractPoint{3,T}}, faces::AbstractVector{F};
                 normaltype=Vec{3,T}) where {T,F<:NgonFace}
    return normals(vertices, faces, normaltype)
end

function normals(vertices::AbstractVector{<:AbstractPoint{3,T}}, faces::AbstractVector{F},
                 ::Type{N}) where {T,F<:NgonFace,N}
    normals_result = zeros(N, length(vertices))
    for face in faces
        v = metafree.(vertices[face])
        # we can get away with two edges since faces are planar.
        n = orthogonal_vector(v[1], v[2], v[3])
        for i in 1:length(F)
            fi = face[i]
            normals_result[fi] = normals_result[fi] + n
        end
    end
    normals_result .= normalize.(normals_result)
    return normals_result
end


text:
#=

Hard to test these functions other than visually.  This script simply
runs then all to make sure that they at least run

Set the variable 'disp' to true to display the processed images
=#

disp = false

using Test, ImagePhaseCongruency, TestImages

if disp
    using PyPlot
end

println("Testing phase congruency functions...")

img = testimage("lena_gray")
img = Float64.(img)
disp ? imshow(img, ColorMap("gray")) : nothing

println("Phase preserving dynamic range compresion")
dimg = ppdrc(img, 50; clip=0.01, n=2)
disp ? imshow(dimg, ColorMap("gray")) : nothing

dimg = ppdrc(img, geoseries((20,60),3))
disp ? imshow(dimg[1], ColorMap("gray")) : nothing

println("phasecongmono")
 (PC, or, ft, T) =
         phasecongmono(img; nscale=4, minwavelength=3, mult=2,
                        sigmaonf=0.55, k=3, cutoff=0.5, g=10,
                        deviationgain=1.5, noisemethod=-1)

PC = phasecongmono(img, nscale=3)[1]
disp ? imshow(PC, ColorMap("gray")) : nothing

println("phasesymmono")
 (phaseSym, symmetryEnergy, T) =
           phasesymmono(img; nscale=3, minwavelength=3, mult=2,
                        sigmaonf=0.55, k=2, polarity=1, noisemethod=-1)

disp ? imshow(phaseSym, ColorMap("gray")) : nothing

phaseSym = phasesymmono(img)[1]
disp ? imshow(phaseSym, ColorMap("gray")) : nothing

println("phasecong3")
(M, m, or, ft, EO, T) = phasecong3(img; nscale=3, norient=6, minwavelength=3,
                          mult=2, sigmaonf=0.55, k=2, cutoff=0.5,
                          g = 10, noisemethod=-1)

disp ? imshow(M, ColorMap("gray")) : nothing
disp ? imshow(m, ColorMap("gray")) : nothing

println("phasesym")
(phaseSym, orient, totalEnergy, T) =
            phasesym(img; nscale = 5, norient = 6, minwavelength = 3, mult= 2,
                     sigmaonf = 0.55, k = 2, polarity = 0, noisemethod = -1)
disp ? imshow(phaseSym, ColorMap("gray")) : nothing

println("ppdenoise")
cleanimage = ppdenoise(img,  nscale = 5, norient = 6,
                              mult = 2.5, minwavelength = 2, sigmaonf = 0.55,
                              dthetaonsigma = 1.0, k = 3, softness = 1.0)
disp ? imshow(cleanimage, ColorMap("gray")) : nothing

println("monofilt")
nscale = 4
norient = 6
minWaveLength = 3
mult = 2
sigmaOnf = 0.55
dThetaOnSigma = 1.3
Lnorm = 0
orientWrap = false

(f, h1f, h2f, A, theta, psi) =
            monofilt(img, nscale, minWaveLength, mult, sigmaOnf, orientWrap)

println("gaborconvolve")
(EO, BP) = gaborconvolve(img,  nscale, norient, minWaveLength, mult,
                                 sigmaOnf, dThetaOnSigma, Lnorm)

disp ? imshow(real.(EO[nscale,1]), ColorMap("gray")) : nothing
disp ? imshow(BP[nscale], ColorMap("gray")) : nothing

println("highpassmonogenic and bandpassmonogenic")
minwavelength = 4
maxwavelength = 20
n = 4
(ph, orient, E) = highpassmonogenic(img, maxwavelength, n)
disp ? imshow(ph, ColorMap("gray")) : nothing
disp ? imshow(orient, ColorMap("gray")) : nothing
disp ? imshow(E, ColorMap("gray")) : nothing
(ph, orient, E) = bandpassmonogenic(img, minwavelength, maxwavelength, n)
disp ? imshow(ph, ColorMap("gray")) : nothing
disp ? imshow(orient, ColorMap("gray")) : nothing
disp ? imshow(E, ColorMap("gray")) : nothing

nothing


text:
"""
   extract(modelfile)

Return a [`ONNX.Types.Model`](@ref) and a Dict mapping input variables to size tuples (in Flux order).

Beware that missing/variable size data for a dimension results in a random size for that dimension. Therefore sizes should mostly be used to determine the number of dimensions.
"""
extract(modelfile::AbstractString) = open(io -> extract(io), modelfile)
extract(io::IO) = ONNX.readproto(io, ONNX.ModelProto())

NaiveNASlib.name(n::OnnxNode) = name(n.proto)
NaiveNASlib.name(n::ONNX.NodeProto) = n.name
NaiveNASlib.name(vip::ONNX.ValueInfoProto) = vip.name
NaiveNASlib.name(tp::ONNX.TensorProto) = tp.name

load(filename::String, insizes...; kwargs...) = open(io -> load(io, insizes...; kwargs...), filename)
load(io::IO, insizes...; kwargs...) = load(extract(io), insizes...; kwargs...)
load(m::ONNX.ModelProto, insizes...; kwargs...) = load(m.graph, insizes...; kwargs...)
load(g::ONNX.GraphProto, insizes...; kwargs...) = CompGraph(g, insizes...; kwargs...)
NaiveNASlib.CompGraph(g::ONNX.GraphProto, insizes...; kwargs...) = CompGraph(CompGraphBuilder(g, insizes...); kwargs...)
function NaiveNASlib.CompGraph(gb::CompGraphBuilder; vfun = create_vertex_default, infer_shapes=true)
    outputs::Vector{AbstractVertex} = vertex.(gb, node.(name.(gb.g.output), gb), vfun)
    graph = CompGraph(gb.inputs, outputs)
    if infer_shapes
        try_infer_sizes!(graph, (get(gb.sizes, n, (missing,)) for n in name.(inputs(graph)))...)
    end
    return graph
end

"""
   vertex(gb::CompGraphBuilder, n::ONNX.Types.Node, vfun = create_vertex_default)

Return an `AbstractVertex` created from `n`.

Inputs to the returned vertex are created recursively based on state in `gb`.
"""
function vertex(gb::CompGraphBuilder, n::OnnxNode, vfun = create_vertex_default)
      return get!(gb.created, n) do
         n_create, ins = check_combine(gb, n)
         invertices = map(ni -> vertex(gb, ni, vfun), ins)
         v = vfun(n_create, invertices)
         if is_input(v)
            push!(gb.inputs, v)
         end
         return v
      end
end

is_input(v::AbstractVertex) = is_input(base(v))
is_input(v::InputVertex) = true
is_input(v::CompVertex) = false
is_input(v::SourceVertex) = false

create_vertex_default(n::OnnxNode, inputs::Array; kwargs...) = verts[optype(n)](name(n), inputs, n.attribute, params(n)...; kwargs...)


text:
# This file is available under the terms of the MIT License

using LinearAlgebra
using Random
using Printf
using Optim
using Munkres

@doc """
Utility functions for evaluating exponential basis,
generating synthetic examples, evaluating derivatives,
etc. Many depend on BLAS for speed
"""

###########################################################
# return the column view (complex or real) of a matrix
function col_view(X)
    ndims(X) == 1 && (return X);
    m,n = size(X);
    T   = eltype(X);
    s   = sizeof(T);
    p   = pointer(X);
    x   = Array{Vector{T}, 1}(UndefInitializer(), n);
    for i = 1:n
        x[i] = unsafe_wrap(Array, p, m);
        p += m*s;
    end
    return x
end

function col_view_real(X)
    # assert X is complex
    T = typeof(real(X[1]));
    p = convert(Ptr{T}, pointer(X));
    s = sizeof(T);
    if ndims(X) == 1
        m = length(X);
        x = unsafe_wrap(Array, p, 2*m);
        return x
    else
        m,n = size(X);
        x   = Array{Vector{T},1}(UndefInitializer(), n);
        for i = 1:n
            x[i] = unsafe_wrap(Array, p, 2*m);
            p += 2*m*s;
        end
        return x
    end
end

###########################################################
# populate exponential matrix
function update_P!(params)
    T  = eltype(params.X);
    BLAS.gemm!('N', 'T', T(1.0), params.t, params.a, T(0.0), params.P);
    map!(exp, params.P, params.P);
end

function update_P_general!(P, t, a, f)
    T  = eltype(a);
    c0 = zero(T);
    c1 = one(T);
    BLAS.gemm!('N', 'T', c1, t, a, c0, P);
    map!(f, P, P);
end

# ############################################################
# # update residuals
# function update_R!(vars, params)
#     T  = eltype(params.X);
#     c0 = zero(T);
#     c1 = one(T);

#     BLAS.gemm!('N', 'N', c1, vars.P, vars.B, c0, vars.R);
#     broadcast!(-, vars.R, vars.R, params.X);
# end

# function update_r!(vars, params, id)
#     T  = eltype(params.X);
#     c0 = zero(T);
#     c1 = one(T);

#     BLAS.gemm!('N', 'N', c1, vars.P, vars.b[id], c0, vars.r[id]);
#     broadcast!(-, vars.r[id], vars.r[id], params.x[id]);
# end

# ############################################################
# # update the QR factorization of P
# function update_PQR!(vars, params, svars)
#     P = vars.P;
#     T = eltype(P);
#     PQ = svars.PQ;
#     PR = svars.PR;
#     tP = svars.tP;

#     c1 = one(T);
#     c0 = zero(T);

#     # calculate QR decomposition
#     copy!(PQ, P);
#     LAPACK.geqrf!(PQ, tP);
#     LAPACK.orgqr!(PQ, tP);
#     BLAS.gemm!('C', 'N', c1, PQ, P, c0, PR);
# end

# # solve upper triangular linear system
# function upper_solve!(PR, b)
#     k = length(b);
#     # backsubtitution
#     b[k] = b[k]/PR[k,k];
#     for i = k-1:-1:1
#         # calculate the rhs
#         for j = i+1:k
#             b[i] -= PR[i,j]*b[j];
#         end
#         b[i] = b[i]/PR[i,i];
#     end
# end

# ############################################################
# # update B, b
# function update_B!(vars, params, svars)
#     T  = eltype(params.X);
#     c0 = zero(T);
#     c1 = one(T);
#     BLAS.gemm!('C', 'N', c1, svars.PQ, params.X, c0, vars.B);
#     for i = 1:params.n
#         upper_solve!(svars.PR, vars.b[i]);
#     end
# end

# function update_b!(vars, params, svars, id)
#     T  = eltype(params.X);
#     c0 = zero(T);
#     c1 = one(T);
#     BLAS.gemv!('C', c1, svars.PQ, params.x[id], c0, vars.b[id]);
#     upper_solve!(svars.PR, vars.b[id]);
# end

# function grad_ar!(gr, vars, params, svars)
#     #
#     # Helper routine: following the inner solve,
#     # this routine computes the gradient w.r.t alpha
#     #
#     # This routine assumes that alphar, B, and R
#     # are up-to-date
#     #
#     # NOTE: this function overwrites R to save space
#     #

#     t = params.t;
#     P = vars.P;
#     B = vars.B;
#     R = vars.R;
#     tM = svars.tM;

#     T = eltype(params.X);

#     # wrap complex array around galphar
#     pr = pointer(gr);
#     pc = convert(Ptr{T}, pr);
#     gc = unsafe_wrap(Array, pc, params.k);
#     # compute complex gradient
#     c0 = zero(T);
#     c1 = T(0.5);
#     conj!(R);
#     broadcast!(*, R, R, t);
#     BLAS.gemm!('T', 'N', c1, P, R, c0, tM);
#     broadcast!(*, tM, tM, B);
#     BLAS.sum!(gc, tM);

#     gc2gr!(gr, params.k);
# end

# function grad_ar!(gr, vars, params, svars, id)
#     #
#     # Helper routine: following the inner solve,
#     # this routine computes the gradient w.r.t alpha
#     #
#     # This routine assumes that alphar, B, and R
#     # are up-to-date
#     #
#     # NOTE: this function overwrites R to save space
#     #

#     t = params.t;
#     P = vars.P;
#     b = vars.b;
#     r = vars.r;

#     T = eltype(params.X);

#     # wrap complex array around galphar
#     pr = pointer(gr);
#     pc = convert(Ptr{T}, pr);
#     gc = unsafe_wrap(Array, pc, params.k);
#     # compute complex gradient
#     c0 = zero(T);
#     c1 = T(0.5);
#     conj!(r[id]);
#     broadcast!(*, r[id], r[id], t);
#     BLAS.gemv!('T', c1, P, r[id], c0, gc);
#     broadcast!(*, gc, gc, b[id]);

#     gc2gr!(gr, params.k);
# end

# function gc2gr!(gr, k)
#     T = eltype(gr);
#     scale!(gr, T(-2.0));
#     BLAS.scal!(k, T(-1.0), gr, 2);
# end

function vc2vr(vc)
    T = eltype(real(vc[1]));
    k = length(vc);
    #
    pc = pointer(vc);
    pr = convert(Ptr{T}, pc);
    vr = unsafe_wrap(Array, pr, k<<1)
    #
    return vr
end

function vr2vc(vr)
    T = Complex{eltype(vr)};
    k = length(vr);
    #
    pr = pointer(vr);
    pc = convert(Ptr{T}, pr);
    vc = unsafe_wrap(Array, pc, k>>1)
    return vc
end

# simulate a simple synthetic example
#------------------------------------------------------------------------------
function simDMD(m, n, k, T; seed=123)
    Random.seed!(seed);
    #
    # time and space vector
    t = complex(collect(range(T(0.0), stop=T(1.0), length=m)));
    s = complex(collect(range(T(-pi), stop=T(pi), length=n)));
    #
    # data matrix
    c0 = zero(Complex{T});
    c1 = one(Complex{T});
    #
    at = im*randn(T,k);       # temporal modes
    bt = complex((rand(T,k))*n/10.0);   # spatial modes
    # phit = exp(t⋅alphatᵀ), psit = sin(s⋅betatᵀ)
    Pt = zeros(Complex{T},m,k); update_P_general!(Pt, t, at, exp);
    Bt = zeros(Complex{T},k,n); update_P_general!(Bt, bt, s, sin);
    # X = Pt⋅Bt
    X = zeros(Complex{T}, m, n);
    BLAS.gemm!('N', 'N', c1, Pt, Bt, c0, X);
    #
    return X, t, at, Bt
end

# ###########################################################
# # closed form solution of B for least squares problem

function dmdl2B!(B, alpha, m, n, k, X, t)
    T  = eltype(X);
    epsmin = 1e2*eps(typeof(real(X[1])));
    c0 = zero(T);
    c1 = one(T);
    phi = zeros(T,m,k);
    update_P_general!(phi, t, alpha, exp);

    # stabilized least squares solution

    F = svd(phi,full=false)

    s1 = maximum(F.S)
    k2 = sum(F.S .> s1*epsmin)

    Y = zeros(T,k2,n)
    U = view(F.U,:,1:k2)
    Vt = view(F.Vt,1:k2,:)
    BLAS.gemm!('C','N',c1,U,X,c0,Y)
    Y ./= F.S[1:k2]
    BLAS.gemm!('C','N',c1,Vt,Y,c0,B)

end

# ###########################################################
# # exact and trapezoidal dmd --- for generating initial guess
function dmdexactestimate(m,n,k,X,t;dmdtype="exact",
                          uselowrank=false,niter=0)
    # use the trapezoidal rule and exact DMD
    # to estimate eigenvalues
    # Assumes that the times are in order,
    # i.e. that t[i] < t[i+1]
    T  = eltype(X);
    Tr = typeof(real(X[1]));

    if (dmdtype == "exact")
        dt = t[2]-t[1]

        if uselowrank
            x1op = transpose(LinearOperator(X[1:end-1,:]))
            F = psvdfact(x1op,rank=k,sketch = :randn,
                         sketch_randn_niter = niter);
        else
            x1 = transpose(X[1:end-1,:])
            F = svd(x1,full=false);
        end
        u = F.U[:,1:k]; s = Diagonal(Array{T}(F.S[1:k]));
        v = F.Vt[1:k,:]';

        atilde = u'*(transpose(X[2:end,:])*v)/s
        a = eigvals(atilde)
        a = log.(a)/dt

        B = zeros(T,k,n)
        dmdl2B!(B,a,m,n,k,X,t)

    elseif (dmdtype == "trap")

        if uselowrank
            x1op = LinearOperator(X[1:end-1,:])
            x2op = LinearOperator(X[2:end,:])
            midop = transpose(x1op)+transpose(x2op)
            F = psvdfact(midop,rank=k,sketch=:randn,
                         sketch_randn_niter=niter)
        else
            mid = transpose(X[1:end-1,:] + X[2:end,:])
            F = svd(mid, full=false)
        end

        u = F.U[:,1:k]
        s = Diagonal(Array{T}(F.S[1:k]))
        dt = Diagonal(Array{T}(2.0./(t[2:end]-t[1:end-1])))
        v = dt*(F.Vt[1:k,:]')

        atilde = u'*(transpose(X[2:end,:])*v - transpose(X[1:end-1,:])*v)/s
        a = eigvals(atilde)
        B = zeros(T,k,n)

        dmdl2B!(B,a,m,n,k,X,t)

    else

        error("unknown value for dmdtype")

    end

    return a, B
end

###########################################################
# error measure
function besterrperm(v1,v2)
    n = length(v1)
    A = Array{typeof(abs(v1[1]))}(undef,n,n)
    for j = 1:n
        for i = 1:n
            A[i,j] = abs(v1[i]-v2[j])
        end
    end
    if any(isnan,A)
        println("besterrperm: NaN in A, abort")
        return Inf
    end

    p = munkres(A)
    err = 0.0
    for i = 1:n
        err = err + abs(v1[i]-v2[p[i]])
    end
    return err
end

function besterrperm_wi(v1,v2)
    n = length(v1)
    A = Array{typeof(abs(v1[1]))}(undef,n,n)
    for j = 1:n
        for i = 1:n
            A[i,j] = abs(v2[i]-v1[j])
        end
    end
    p = munkres(A)
    err = 0.0
    for i = 1:n
        err = err + abs(v1[p[i]]-v2[i])
    end
    return err, p
end


text:
immutable MExpr{Head}
  args::Vector
end
MExpr(head, args...) = MExpr{head}([args...])

import Base.show

function show{T}(io::IO, e::MExpr{T})
  print(io, T)
  if length(e.args) >= 1
    print(io, "["); show(io, e.args[1])
    for x in e.args[2:end] print(io, ", "); show(io, x) end
    print(io, "]")
  else
    print(io, "[]")
  end
end

show(io::IO, r::MExpr{:Rule}) =
  (show(io, r.args[1]); print(io, "→"); show(io, r.args[2]))

# Conversion

const aliases =
  Dict(:*   => :Times,
       :/   => :Divide,
       :^   => :Power,
       :+   => :Plus,
       :-   => :Subtract,
       :%   => :Mod,
       :log => :Log,
       :exp => :Exp,
       :sin => :Sin,
       :cos => :Cos,
       :tan => :Tan)

from_mma(x) = x
const symbols = Dict(:True => true, :False => false, :Null => nothing)
from_mma(s::Symbol) = haskey(symbols, s) ? symbols[s] : s

to_mma{T<:Union{Int64,Int32,Float64,Float32,Symbol,AbstractString}}(x::T) = x

function to_mma(x::Expr)
  if x.head == :call
    head = haskey(aliases, x.args[1]) ? aliases[x.args[1]] : x.args[1]
    MExpr{head}(map(to_mma, x.args[2:end]))
  elseif x.head == :block
    MExpr{:CompoundExpression}(map(to_mma, x.args))
  elseif x.head == :ref
    MExpr{:Part}(map(to_mma, x.args))
  elseif x.head == :cell1d
    MExpr{:List}(x.args)
  else
    error("Unsupported $(x.head) expression.")
  end
end

function to_mma(x::QuoteNode)
  typeof(x.value) == Symbol && return x.value
  error("Cannot call to_mma on QuoteNode($(x.value))")
end

# Other types

to_mma(x::Bool) = x ? :True : :False

to_mma{T}(x::MExpr{T}) = MExpr{T}(map(to_mma, x.args))
from_mma{T}(x::MExpr{T}) = MExpr{T}(map(from_mma, x.args))

to_mma(x::Rational) = MExpr(:Rational, x.num, x.den)
from_mma(f::MExpr{:Rational}) = f.args[1]//f.args[2]

to_mma(x::Complex) = MExpr(:Complex, to_mma(real(x)), to_mma(imag(x)))
from_mma(x::MExpr{:Complex}) = Complex(x.args[1], x.args[2])

to_mma(xs::Vector) = MExpr{:List}(map(to_mma, xs))
from_mma(l::MExpr{:List}) = map(from_mma, l.args)

# Julia Expression Conversion

to_expr(x) = x
to_expr{T}(x::MExpr{T}) = Expr(:call, T, map(to_expr, x.args)...)

for (j, m) in aliases
  j = Expr(:quote, j)
  m = Expr(:quote, m)
  @eval to_expr(x::MExpr{$m}) = Expr(:call, $j, map(to_expr,x.args)...)
end

import Base.convert
convert(::Type{Expr}, x::MExpr) = to_expr(from_mma(x))
convert(::Type{MExpr}, x::Expr) = to_mma(x)


text:
import .Frames # The dot is due to the fact that Frames is a submodule
using RecipesBase
# default(markersize=1)

@userplot TrajPlot
"""`trajplot(T,args...)` Plots a trajectory of T-matrices in a single plot"""
trajplot
@recipe function trajplot(t::TrajPlot)
    title --> "Trajectory"
    layout --> 1
    T = t.args[1]
    label --> ["x" "y" "z"]
    T[1:3,4,:]'
end

@userplot TrajPlot3
"""`plot_traj3(T, ls=".", plotFrame = 0)` Plots a trajectory of T-matrices in a 3D plot, with otional frames drawn of length `plotFrame`"""
trajplot3
@recipe function trajplot3(t::TrajPlot3, plotFrame = 0)
    T = t.args[1]
    xguide --> "x"
    yguide --> "y"
    zguide --> "z"
    if ndims(T) == 3
        @series begin
            seriestype := :scatter3d
            (T[1,4,:],T[2,4,:],T[3,4,:])
        end
        if plotFrame > 0
            for i = 1:size(T,3)
                @series begin
                    plot!(T[:,:,i],plotFrame)
                end
            end
        end
    else
        @series begin
            seriestype := :scatter3d
            TT = reduce(hcat, Robotlib.T2t.(T))
            TT[1,:],TT[2,:],TT[3,:]
        end
        if plotFrame > 0
            for i = 1:size(T,3)
                @series begin
                    plot!(T[i],plotFrame)
                end
            end
        end
    end
end

max_forks_repo_name: baggepinnen/Robotlib


text:
function emldescriptors(x, q, t, a, b, c, pbc, eta, kappa, emdescriptors, mldescriptors=nothing)

    if length(emdescriptors) > 0
        rcutmax = 0.0
        for i = 1:length(emdescriptors)
            rcutmax = max(rcutmax, emdescriptors[i].rcut)
        end
        ed, edd, edv = empiricaldescriptors(x, q, t, a, b, c, pbc, rcutmax, eta, kappa, emdescriptors)
        sz = size(edd)
        if length(sz)==3
            edd = reshape(edd, (sz[1]*sz[2], sz[3]))
        else
            edd = reshape(edd, (sz[1]*sz[2], 1))
        end
    end

    # bispectrum, bispectrum derivatives, bispectrum virial
    if mldescriptors !== nothing
        if mldescriptors.name == "SNAP"
            rcutmax = 0.0
            for ielem = 1:mldescriptors.ntypes
                rcutmax = max(2.0*mldescriptors.radelem[1+ielem]*mldescriptors.rcutfac,rcutmax)
            end
            mldescriptors.rcutmax = rcutmax
            mld, mldd, mldv = snapdescriptors(x, t, a, b, c, pbc, mldescriptors)
            sz = size(mldd)
            mldd = reshape(mldd, (sz[1]*sz[2], sz[3]))
        end
    end

    if (length(emdescriptors) > 0) & (mldescriptors !== nothing)
        d = [ed[:]; mld[:]]
        dd = cat(edd, mldd, dims=2);
        dv = cat(edv, mldv, dims=2);
        return d, dd, dv
    elseif length(emdescriptors) > 0
        return ed, edd, edv
    elseif (mldescriptors !== nothing)
        return mld, mldd, mldv
    else
        error("empirical descriptors and ML descriptors are empty")
    end
end

function emldescriptors(config, indices, eta, kappa, normalizeenergy, emdescriptors,
                mldescriptors=nothing, potential=nothing)

normalizestress = 0
# cumalative sum of numbers of atoms
nconfigs = length(config.natom)
natom = [0; cumsum(config.natom[:])];

nf = zeros(Int64, length(indices))
be = []
bf = []
bs  = []
Ae = []
Af = []
As = []
for i = 1:length(indices)
    ci = indices[i]
    if ci > nconfigs
        error("Indices are greater than nconfigs")
    end

    x = config.x[:,(natom[ci]+1):natom[ci+1]] # atom positions of configuration ci
    t = config.t[:,(natom[ci]+1):natom[ci+1]] # atom types of configuration ci
    if length(config.q)>0
        q = config.q[:,(natom[ci]+1):natom[ci+1]]
    else
        q = []
    end
    if size(config.a,2) > 1
        a = config.a[:,ci]
        b = config.b[:,ci]
        c = config.c[:,ci]
    else
        a = config.a[:,1]
        b = config.b[:,1]
        c = config.c[:,1]
    end
    if size(config.pbc,2) > 1
        pbc = config.pbc[:,ci]
    else
        pbc = config.pbc[:,1]
    end

    d, dd, dv = emldescriptors(x, q, t, a[:], b[:], c[:], pbc[:], eta, kappa, emdescriptors, mldescriptors)

    if normalizestress==1
        vol = latticevolume(a[:], b[:], c[:])
        dv = (1.6021765e6/vol) * dv
    end

    if potential !== nothing
        rcutmax = 0.0
        for j = 1:length(potential)
            rcutmax = max(rcutmax, potential[j].rcut)
        end
        eref, e, fref, v, sref = empiricalpotential(x, q, t, a[:], b[:], c[:], pbc[:], rcutmax, eta, kappa, potential)
        if normalizestress==1
            vol = latticevolume(a[:], b[:], c[:])
            sref = (1.6021765e6/vol) * sref
        end
    else
        eref = 0.0;
        fref = [0.0];
        sref = [0.0];
    end

    if normalizeenergy==1
        d = reshape(d, (1,length(d)))/config.natom[ci]
    else
        d = reshape(d, (1,length(d)))
    end

    nf[i] = size(dd,1)
    if i == 1
        Ae = d
        Af = dd
        As = dv
    else
        Ae = cat(Ae, d, dims=1)
        Af = cat(Af, dd, dims=1)
        As = cat(As, dv, dims=1)
    end

    if length(config.e)>0
        e = config.e[indices[i]] - eref
        if normalizeenergy==1
            e = e/config.natom[indices[i]]
        end
        be = [be; e]
    end
    if length(config.f)>0
        f = config.f[:,(natom[ci]+1):natom[ci+1]]
        bf = [bf; f[:] .- fref[:]]
    end
    if length(config.stress)>0
        if size(config.stress,1)==9
            ind = [1, 5, 9, 6, 3, 2]
            s = config.stress[ind,ci]
        else
            s = config.stress[:,ci]
        end
        bs = [bs; s[:] .- sref[:]]
    end
end

return Ae, Af, As, be, bf, bs, nf

end

function emldescriptors2(config, indices, eta, kappa, normalizeenergy, normalizestress,
        descriptors, potential=nothing)

# getdescriptors
emdescriptors, mldescriptors = getdescriptors(descriptors)

# cumalative sum of numbers of atoms
nconfigs = length(config.natom)
natom = [0; cumsum(config.natom[:])];

nf = zeros(Int64, length(indices))
be = []
bf = []
bs  = []
Ae = []
Af = []
As = []
for i = 1:length(indices)
    ci = indices[i]
    if ci > nconfigs
        error("Indices are greater than nconfigs")
    end

    x = config.x[:,(natom[ci]+1):natom[ci+1]] # atom positions of configuration ci
    t = config.t[:,(natom[ci]+1):natom[ci+1]] # atom types of configuration ci
    if length(config.q)>0
        q = config.q[:,(natom[ci]+1):natom[ci+1]]
    else
        q = []
    end
    if size(config.a,2) > 1
        a = config.a[:,ci]
        b = config.b[:,ci]
        c = config.c[:,ci]
    else
        a = config.a[:,1]
        b = config.b[:,1]
        c = config.c[:,1]
    end
    if size(config.pbc,2) > 1
        pbc = config.pbc[:,ci]
    else
        pbc = config.pbc[:,1]
    end

    d, dd, dv = emldescriptors(x, q, t, a[:], b[:], c[:], pbc[:], eta, kappa, emdescriptors, mldescriptors)

    if normalizestress==1
        vol = latticevolume(a[:], b[:], c[:])
        dv = (1.6021765e6/vol) * dv
    end

    if (potential !== nothing)
        if (length(potential) > 0)
            rcutmax = 0.0
            for j = 1:length(potential)
                rcutmax = max(rcutmax, potential[j].rcut)
            end
            eref, e, fref, v, sref = empiricalpotential(x, q, t, a[:], b[:], c[:], pbc[:], rcutmax, eta, kappa, potential)
            if normalizestress==1
                vol = latticevolume(a[:], b[:], c[:])
                sref = (1.6021765e6/vol) * sref
            end
        else
            eref = 0.0; fref = [0.0]; sref = [0.0];
        end
    else
        eref = 0.0; fref = [0.0]; sref = [0.0];
    end

    if normalizeenergy==1
        d = reshape(d, (1,length(d)))/config.natom[ci]
    else
        d = reshape(d, (1,length(d)))
    end

    nf[i] = size(dd,1)
    if i == 1
        Ae = d
        Af = dd
        As = dv
    else
        Ae = cat(Ae, d, dims=1)
        Af = cat(Af, dd, dims=1)
        As = cat(As, dv, dims=1)
    end

    if length(config.e)>0
        e = config.e[indices[i]] - eref
        if normalizeenergy==1
            e = e/config.natom[indices[i]]
        end
        be = [be; e]
    end
    if length(config.f)>0
        f = config.f[:,(natom[ci]+1):natom[ci+1]]
        bf = [bf; f[:] .- fref[:]]
    end
    if length(config.stress)>0
        if size(config.stress,1)==9
            ind = [1, 5, 9, 6, 3, 2]
            s = config.stress[ind,ci]
        else
            s = config.stress[:,ci]
        end
        bs = [bs; s[:] .- sref[:]]
    end

end

return Ae, Af, As, be, bf, bs, nf

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
# The filter version below would be great if the dict keys were already
# Symbols, but they're not 'cause python =\
function sub_dict(dict, keys)
    Dict(
        Symbol(index)=>value for (index, value) in pairs(dict) if index ∈ keys
    )
end

# d_gpts = filter(
#   p -> p.first ∈ ["p", "t0", "P", "rho", "inc", "b", "aRs", "q1"],
#   d_pkl_gpts
# )


text:
using Documenter, BioMedQuery

makedocs()

deploydocs(
    deps   = Deps.pip("mkdocs", "python-markdown-math", "mkdocs-material"),
    repo = "github.com/bcbi/BioMedQuery.jl.git",
    julia  = "0.4",
    osname = "osx"
)


text:
pkgdir = Pkg.dir("SLEEF")
libext = is_linux() ? ".so" :
		 is_apple()  ? ".dynlib" :
		 is_windows() ? ".dll" : error("OS couldn't be determined")
lib_dir = joinpath(pkgdir, "deps","usr","lib")

cp(joinpath(pkgdir, "deps"), joinpath(pkgdir, "_deps"); remove_destination=true)

try
	for arch in ["sse2", "avx", "avx2", "fma4", "purec"]
		println("-"^38)
		println("Building for $arch")
		println()

		rm(joinpath(pkgdir, "deps", "src"); force=true, recursive=true)    #Git rid of any compiled files (except the lib)
		rm(joinpath(lib_dir, "sleef"*libext); force=true) #get rid of old one so it rebuilds

		ENV["JLSLEEF_ARCH"] = arch
		buildscriptpath = joinpath(pkgdir, "deps", "build.jl")
		julia = joinpath(JULIA_HOME, "julia") #Deal with julia under aliases
		run(`$julia $buildscriptpath`)

		lib_path = joinpath(splitdir(@__FILE__)[1], "usr", "lib", arch*"_sleef"*libext)
		mv(joinpath(lib_dir,"sleef"*libext), lib_path)

		sleef_inner_path = joinpath(pkgdir,"src", "sleef_inner.jl")

		module_name = "SLEEF_"*arch
		#module_path = joinpath(pkgdir,"src", module_name*".jl")
		module_path = joinpath(splitdir(@__FILE__)[1], "arch_modules", module_name*".jl")

		open(module_path, "w") do fp
			println(fp, """
module $(module_name)
# This is a AutoGenerated file. Do not edit directly. Generated by $(@__FILE__)

# Macro to load a library
macro checked_lib(libname, path)
	((VERSION >= v"0.4.0-dev+3844" ? Base.Libdl.dlopen_e : Base.dlopen_e)(path) == C_NULL) && error("Unable to load \n\n\$libname (\$path)\n\nPlease re-run Pkg.build(package), and restart Julia.")
	quote const \$(esc(libname)) = \$path end
end

# Load dependencies
@checked_lib jl_sleef "$(lib_path)"

# Load-hooks

const SIMD_ARCH = "$arch"

include(\"$(sleef_inner_path)\")

end #module
""")
		end

	end
	println("Done")

finally
	rm(joinpath(pkgdir, "deps"); force=true, recursive=true)
	cp(joinpath(pkgdir, "_deps"), joinpath(pkgdir, "deps"); remove_destination=true)

	println("Deps restored")
end


text:
χNmapAB = Dict(Set((:A,:B))=>20.0)
χNmapABS = Dict(Set((:A,:B))=>20.0, Set((:A, :S))=>80.0, Set((:B, :S))=>120.0)
χNmapABS_bad = Dict(Set((:A,:B))=>20.0, Set((:A, :S))=>80.0)
χNmapABS_singular = Dict(Set((:A,:B))=>20.0, Set((:A, :S))=>0.0, Set((:B, :S))=>20.0)

@testset "chiN.jl: _species" begin
    @test Polymer._species(χNmapAB) == [:A, :B]
    @test Polymer._species(χNmapABS) == [:A, :B, :S]
end

@testset "chiN.jl: _check_consistency" begin
    @test Polymer._check_consistency(χNmapAB)
    @test Polymer._check_consistency(χNmapABS)
    @test !Polymer._check_consistency(χNmapABS_bad)
end

@testset "chiN.jl: _χNmap_to_matrix" begin
    mat = Polymer._χNmap_to_matrix(χNmapAB)
    @test mat == [0 20.0; 20.0 0]
    mat = Polymer._χNmap_to_matrix(χNmapABS)
    @test mat == [0 20.0 80.0; 20.0 0 120.0; 80.0 120.0 0]

    @test_throws ErrorException Polymer._χNmap_to_matrix(χNmapABS_singular)
end

@testset "chiN.jl: χNMatrix" begin
    m = χNMatrix(χNmapAB)
    @test m == [0 20.0; 20.0 0]
    m = χNMatrix(χNmapABS)
    @test m == [0 20.0 80.0; 20.0 0 120.0; 80.0 120.0 0]
    @test size(m) == (3,3)

    @test_throws ErrorException χNMatrix(χNmapABS_bad)
    @test_throws ErrorException χNMatrix(χNmapABS_singular)
end

@testset "chiN.jl: getindex" begin
    m = χNMatrix(χNmapABS)
    @test m[1] == 0
    @test m[7] == 80.0
    @test m[1, 2] == 20.0
    @test m[2, 1] == 20.0
    @test m[1, 3] == 80.0
    @test m[2, 3] == 120.0
    @test m[:A] == 0
    @test m[:B, :B] == 0
    @test m[:B, :A] == 20.0
    @test m[:A, :S] == 80.0
    @test m[:S, :B] == 120.0
end

@testset "chiN.jl: setindex" begin
    m = χNMatrix(χNmapABS)
    m[1] = 10.0
    @test m[1] == 0
    @test m[:A] == 0
    @test m[:A, :A] == 0

    m[7] = 88.0
    @test m[7] == 88.0
    @test m[1, 3] == 88.0
    @test m[3, 1] == 88.0

    m[:A, :A] = 99.0
    @test m[1] == 0
    @test m[1, 1] == 0
    @test m[:A] == 0
    @test m[:A, :A] == 0

    m[:S, :B] = 102.0
    @test m[:B, :S] == 102.0
    @test m[2, 3] == 102.0
    @test m[3, 2] == 102.0
end


text:
struct SAFTVRQMieParam <: EoSParam
    segment::SingleParam{Float64}
    sigma::PairParam{Float64}
    lambda_a::PairParam{Float64}
    lambda_r::PairParam{Float64}
    epsilon::PairParam{Float64}
    Mw::PairParam{Float64}
end

abstract type SAFTVRQMieModel <: SAFTVRMieModel end
@newmodel SAFTVRQMie SAFTVRQMieModel SAFTVRQMieParam

export SAFTVRQMie
function SAFTVRQMie(components; idealmodel=BasicIdeal, userlocations=String[], ideal_userlocations=String[], verbose=false)
    params = getparams(components, ["SAFT/SAFTVRQMie"]; userlocations=userlocations, verbose=verbose)

    params["Mw"].values .*= 1E-3
    Mw = sigma_LorentzBerthelot(params["Mw"])
    segment = params["m"]
    params["sigma"].values .*= 1E-10
    sigma = sigma_LorentzBerthelot(params["sigma"])
    epsilon = epsilon_HudsenMcCoubrey(params["epsilon"], sigma)
    lambda_a = lambda_LorentzBerthelot(params["lambda_a"])
    lambda_r = lambda_LorentzBerthelot(params["lambda_r"])

    packagedparams = SAFTVRQMieParam(segment, sigma, lambda_a, lambda_r, epsilon, Mw)
    references = ["todo"]

    model = SAFTVRQMie(packagedparams, idealmodel; ideal_userlocations=ideal_userlocations, references=references, verbose=verbose)
    return model
end

function a_mono(model::SAFTVRQMieModel, V, T, z,_data = @f(data))
    _,_,vrdata = _data
    return @f(a_hs,vrdata) + @f(a_disp,_data)
end

function a_res(model::SAFTVRQMieModel, V, T, z)
    @f(a_mono)
end

function data(model::SAFTVRQMieModel, V, T, z)
    _σeff = @f(σeff)
    _ϵff = @f(ϵeff)
    _d = @f(d,_σeff)
    ζi = @f(ζ0123,_d)
    _ζst = @f(ζst,_σeff)
    _ζ_X,_ = @f(ζ_X_σ3,_d)
    _ρ_S = @f(ρ_S)
    σ3x = _ζst/(_ρ_S*π/6)
    vrdata = (_d,_ρ_S,ζi,_ζ_X,_ζst,σ3x)
    return (_σeff,_ϵff,vrdata)
end

function Q1(model::SAFTVRQMieModel, V, T, z, λ)
    return λ*(λ-1)
end

function Q2(model::SAFTVRQMieModel, V, T, z, λ)
    return 1/2*(λ+2)*(λ+1)*λ*(λ-1)
end

function d(model::SAFTVRQMieModel, V, T, z, i,_data)
    ϵ,σ,λr,λa,Mw,σeff = _data

    Di = ħ^2/(12*k_B*T*Mw/N_A*σ^2)
    βi = ϵ/T

    function udud2u(x)
        _C = @f(Cλ,λa,λr)
        Q1λr = @f(Q1,λr)
        Q1λa = @f(Q1,λa)
        Q2λr = @f(Q2,λr)
        Q2λa = @f(Q2,λa)
        _u= _C*((x^-λr-x^-λa)+
            Di*(Q1λr*x^-(λr+2)-Q1λa*x^-(λa+2))+
            Di^2*(Q2λr*x^-(λr+4)-Q2λa*x^-(λa+4)))
        _du = -_C*((λr*x^-(λr+1)-λa*x^-(λa+1))+
            Di*(Q1λr*(λr+2)*x^-(λr+3)-Q1λa*(λa+2)*x^-(λa+3))+
            Di^2*(Q2λr*(λr+4)*x^-(λr+5)-Q2λa*(λa+4)*x^-(λa+5)))
        _d2u = _C*((λr*(λr+1)*x^-(λr+2)-λa*(λa+1)*x^-(λa+2))+
            Di*(Q1λr*(λr+2)*(λr+3)*x^-(λr+4)-Q1λa*(λa+2)*(λa+3)*x^-(λa+4))+
            Di^2*(Q2λr*(λr+4)*(λr+5)*x^-(λr+6)-Q2λa*(λa+4)*(λa+5)*x^-(λa+6)))

        return _u,_du,_d2u
    end

    function f(x)
        _C = @f(Cλ,λa,λr)
        Q1λr = @f(Q1,λr)
        Q1λa = @f(Q1,λa)
        Q2λr = @f(Q2,λr)
        Q2λa = @f(Q2,λa)
        _u= _C*((x^-λr-x^-λa)+
            Di*(Q1λr*x^-(λr+2)-Q1λa*x^-(λa+2))+
            Di^2*(Q2λr*x^-(λr+4)-Q2λa*x^-(λa+4)))
        return exp(-βi*_u)
    end

    function fgh(x)
        _u,_du,_d2u = udud2u(x)
        _f = exp(-βi*_u)
        _g = -βi*_f*_du
        _h = βi*_f*(βi*_du^2-_d2u)
        return _f,_g,_h
    end
    x_min = Solvers.halley(fgh,one(T))
    σ_effi = σeff/σ
    x = SAFTVRQMieconsts.x
    w = SAFTVRQMieconsts.w
    return σ*(σ_effi-(σ_effi-x_min)/2*sum(w[i]*f((σ_effi-x_min)/2*x[i]+(σ_effi+x_min)/2) for i in 1:21))
end

function d(model::SAFTVRQMieModel, V, T, z,_σeff = @f(σeff))
    _ϵ = model.params.epsilon.values
    _σ = model.params.sigma.values
    _λr = model.params.lambda_r.values
    _λa = model.params.lambda_a.values
    _Mwij = model.params.Mw.values

    _d = zeros(typeof(T),length(z))
    for i ∈ 1:length(_d)
        _data = (_ϵ[i,i],_σ[i,i],_λr[i,i],_λa[i,i],_Mwij[i,i],_σeff[i,i])
        _d = @f(d,i,_data)
    end
    _d
end

function σeff(model::SAFTVRQMieModel, V, T, z)
    _σ = model.params.sigma.values
    _λr = model.params.lambda_r.values
    _λa = model.params.lambda_a.values
    Mwij = model.params.Mw.values
    #Dij = ħ^2/(12*k_B*T*Mwij/N_A*σ[i,j]^2)

    #f(x) = Dij^2*@f(Q2,λr[i,j])+Dij*@f(Q1,λr[i,j])*x^2+x^4-x^(λr[i,j]-λa[i,j])*(Dij^2*@f(Q2,λa[i,j])+Dij*@f(Q1,λa[i,j])*x^2+x^4)
    #g(x) = 2*Dij*@f(Q1,λr[i,j])*x+4x^3-x^(λr[i,j]-λa[i,j]-1)*((λr[i,j]-λa[i,j])*Dij^2*@f(Q2,λa[i,j])+(λr[i,j]-λa[i,j]+2)*Dij*@f(Q1,λa[i,j])*x^2+(λr[i,j]-λa[i,j]+4)*x^4)
    #h(x) = 2*Dij*@f(Q1,λr[i,j])+12x^2-x^(λr[i,j]-λa[i,j]-2)*((λr[i,j]-λa[i,j]-1)*(λr[i,j]-λa[i,j])*Dij^2*@f(Q2,λa[i,j])+((λr[i,j]-λa[i,j]-1)*(λr[i,j]-λa[i,j]+2)+2*(λr[i,j]-λa[i,j]+2))*Dij*@f(Q1,λa[i,j])*x^2+((λr[i,j]-λa[i,j]-1)*(λr[i,j]-λa[i,j]+4)+4*(λr[i,j]-λa[i,j]+4))*x^4)

    function fgh(x,λa,λr,σ,T,Mw)
        Q1λr = @f(Q1,λr)
        Q1λa = @f(Q1,λa)
        Q2λr = @f(Q2,λr)
        Q2λa = @f(Q2,λa)
        Dij = ħ^2/(12*k_B*T*Mw/N_A*σ^2)
        xΔλ =x^(λr-λa)
        x2 = x*x
        x4 = x2*x2
        f = Dij^2*Q2λr+Dij*Q1λr*x2+x4-xΔλ*(Dij^2*Q2λa+Dij*Q1λa*x2+x4)
        g = 2*Dij*Q1λr*x+4x^3-(xΔλ/x)*((λr-λa)*Dij^2*Q2λa+(λr-λa+2)*Dij*Q1λa*x2+(λr-λa+4)*x4)
        h = 2*Dij*Q1λr+12x^2-x^(λr-λa-2)*((λr-λa-1)*(λr-λa)*Dij^2*Q2λa+((λr-λa-1)*(λr-λa+2)+2*(λr-λa+2))*Dij*Q1λa*x2+((λr-λa-1)*(λr-λa+4)+4*(λr-λa+4))*x4)
        return f,g,h
    end

    _σeff = zeros(typeof(T),size(_σ))
    n1,n2 = size(_σ)
    for i in 1:n1
        f0 = x -> fgh(x,_λa[i,i],_λr[i,i],_σ[i,i],T,Mwij[i,i])
        _σeff[i,i] =  _σ[i,i]*Solvers.halley(f0,one(T))
        for j in 1:i-1
        f0 = x -> fgh(x,_λa[i,j],_λr[i,j],_σ[i,j],T,Mwij[i,j])
        _σeff[i,j] =  _σ[i,j]*Solvers.halley(f0,one(T))
        _σeff[j,i] = _σeff[i,j]
        end
    end
    return _σeff
end

function ϵeff(model::SAFTVRQMieModel, V, T, z)
    ϵ = model.params.epsilon.values
    _σ = model.params.sigma.values
    _λr = model.params.lambda_r.values
    _λa = model.params.lambda_a.values
    Mwij = model.params.Mw.values#[i] + model.params.Mw.values[j])/2 # check
    #Dij = ħ^2/(12*k_B*T*Mwij/N_A*σ[i,j]^2)

    function fgh(x,λa,λr,σ,T,Mw)
        Q1λr = @f(Q1,λr)
        Q1λa = @f(Q1,λa)
        Q2λr = @f(Q2,λr)
        Q2λa = @f(Q2,λa)
        Dij = ħ^2/(12*k_B*T*Mw/N_A*σ^2)

        xΔλ =x^(λr-λa)
        x2 = x*x
        x4 = x2*x2
        f = (λr+4)*Dij^2*Q2λr+(λr+2)*Dij*Q1λr*x2+λr*x4-x^(λr-λa)*((λa+4)*Dij^2*Q2λa+(λa+2)*Dij*Q1λa*x2+λa*x4)
        g = 2*(λr+2)*Dij*Q1λr*x+4λr*x^3-x^(λr-λa-1)*((λa+4)*(λr-λa)*Dij^2*Q2λa+(λa+2)*(λr-λa+2)*Dij*Q1λa*x2+λa*(λr-λa+4)*x4)
        h = 2*(λr+2)*Dij*Q1λr+12λr*x2-x^(λr-λa-2)*((λa+4)*(λr-λa-1)*(λr-λa)*Dij^2*Q2λa+(λa+2)*((λr-λa-1)*(λr-λa+2)+2*(λr-λa+2))*Dij*Q1λa*x2+λa*((λr-λa-1)*(λr-λa+4)+4*(λr-λa+4))*x4)
    return f,g,h
    end

    function u(x,λa,λr,σ,T,Mw)
        Q1λr = @f(Q1,λr)
        Q1λa = @f(Q1,λa)
        Q2λr = @f(Q2,λr)
        Q2λa = @f(Q2,λa)
        _C = @f(Cλ,λa,λr)
        Dij = ħ^2/(12*k_B*T*Mw/N_A*σ^2)
        return _C*((x^-λr-x^-λa)+
        Dij*(Q1λr*x^-(λr+2)-Q1λa*x^-(λa+2))+
        Dij^2*(Q2λr*x^-(λr+4)-Q2λa*x^-(λa+4)))
    end
    #f(x) = (λr[i,j]+4)*Dij^2*@f(Q2,λr[i,j])+(λr[i,j]+2)*Dij*@f(Q1,λr[i,j])*x^2+λr[i,j]*x^4-x^(λr[i,j]-λa[i,j])*((λa[i,j]+4)*Dij^2*@f(Q2,λa[i,j])+(λa[i,j]+2)*Dij*@f(Q1,λa[i,j])*x^2+λa[i,j]*x^4)
    #g(x) = 2*(λr[i,j]+2)*Dij*@f(Q1,λr[i,j])*x+4λr[i,j]*x^3-x^(λr[i,j]-λa[i,j]-1)*((λa[i,j]+4)*(λr[i,j]-λa[i,j])*Dij^2*@f(Q2,λa[i,j])+(λa[i,j]+2)*(λr[i,j]-λa[i,j]+2)*Dij*@f(Q1,λa[i,j])*x^2+λa[i,j]*(λr[i,j]-λa[i,j]+4)*x^4)
    #h(x) = 2*(λr[i,j]+2)*Dij*@f(Q1,λr[i,j])+12λr[i,j]*x^2-x^(λr[i,j]-λa[i,j]-2)*((λa[i,j]+4)*(λr[i,j]-λa[i,j]-1)*(λr[i,j]-λa[i,j])*Dij^2*@f(Q2,λa[i,j])+(λa[i,j]+2)*((λr[i,j]-λa[i,j]-1)*(λr[i,j]-λa[i,j]+2)+2*(λr[i,j]-λa[i,j]+2))*Dij*@f(Q1,λa[i,j])*x^2+λa[i,j]*((λr[i,j]-λa[i,j]-1)*(λr[i,j]-λa[i,j]+4)+4*(λr[i,j]-λa[i,j]+4))*x^4)

   # σ_min = @f(Halley,f,g,h,(λr[i,j]/λa[i,j])^(1/(λr[i,j]-λa[i,j])))
    #return -ϵ[i,j]*u(σ_min)

    _ϵeff = zeros(typeof(T),size(_σ))
    n1,n2 = size(_σ)
    for i in 1:n1
        f0 = x -> fgh(x,_λa[i,i],_λr[i,i],_σ[i,i],T,Mwij[i,i])
        x0 = (_λr[i,i]/_λa[i,i])^(1/(_λr[i,i]-_λa[i,i]))
        _σmin = Solvers.halley(f0,one(T)*x0)
        uij =u(_σmin,_λa[i,i],_λr[i,i],_σ[i,i],T,Mwij[i,i])
        _ϵeff[i,i] =  -ϵ[i,i]*uij
        for j in 1:i-1
            f0 = x -> fgh(x,_λa[i,j],_λr[i,j],_σ[i,j],T,Mwij[i,j])
            x0 = (_λr[i,j]/_λa[i,j])^(1/(_λr[i,j]-_λa[i,j]))
            _σmin = Solvers.halley(f0,one(T)*x0)
            uij =u(_σmin,_λa[i,j],_λr[i,j],_σ[i,j],T,_Mw[i,j])
            _ϵeff[i,j] =  -ϵ[i,j]*uij
            _ϵeff[j,i] = _ϵeff[i,j]
        end
    end
    return _ϵeff
end

#=
function a_1(model::SAFTVRQMieModel, V, T, z, i, j,_data = @f(data))
    _σeff,_ϵff,_d,_ρ_S,ζi,_ζ_X,_ζst = _data
    ϵ = model.params.epsilon.values
    σ = model.params.sigma.values
    λr = model.params.lambda_r.values
    λa = model.params.lambda_a.values
    Mwij = (model.params.Mw.values[i] + model.params.Mw.values[j])/2 # check
    Dij = ħ^2/(12*k_B*T*Mwij/N_A*σ[i,j]^2)
    x_0ij = @f(x_0,i,j)
    dij = 0.5*(_d[i] + _d[j])
    x_0effij = _σeff[i,j]/dij
    return 2*π*ϵ[i,j]*dij^3*@f(Cλ,λa[i,j],λr[i,j])*@f(ρ_S)*
    ( (x_0ij^λa[i,j]*(@f(aS_1,λa[i,j])+@f(B,λa[i,j],x_0effij))-
        x_0ij^λr[i,j]*(@f(aS_1,λr[i,j])+@f(B,λr[i,j],x_0effij)))+
        (x_0ij^(λa[i,j]+2)*@f(Q1,λa[i,j])*(@f(aS_1,λa[i,j]+2)+@f(B,λa[i,j]+2,x_0effij))-
         x_0ij^(λr[i,j]+2)*@f(Q1,λr[i,j])*(@f(aS_1,λr[i,j]+2)+@f(B,λr[i,j]+2,x_0effij)))*Dij+
        (x_0ij^(λa[i,j]+4)*@f(Q2,λa[i,j])*(@f(aS_1,λa[i,j]+4)+@f(B,λa[i,j]+4,x_0effij))-
         x_0ij^(λr[i,j]+4)*@f(Q2,λr[i,j])*(@f(aS_1,λr[i,j]+4)+@f(B,λr[i,j]+4,x_0effij)))*Dij^2 )
end
=#

#=
function x_0eff(model::SAFTVRQMieModel, V, T, z, i, j)
    return @f(σeff,i,j)/@f(d,i,j)
end

function a_2(model::SAFTVRQMieModel, V, T, z, i, j,_data = @f(data))
    ϵ = model.params.epsilon.values
    σ = model.params.sigma.values
    λr = model.params.lambda_r.values[i,j]
    λa = model.params.lambda_a.values[i,j]
    Mwij = (model.params.Mw.values[i] + model.params.Mw.values[j])/2 # check
    Dij = ħ^2/(12*k_B*T*Mwij/N_A*σ[i,j]^2)
    x_0ij = @f(x_0,i,j)
    x_0effij = @f(x_0eff,i,j)
    return π*@f(KHS)*(1+@f(χ,i,j))*@f(ρ_S)*ϵ[i,j]^2*@f(d,i,j)^3*@f(C,i,j)^2*(x_0ij^(2*λa)*(@f(aS_1,2*λa)+@f(B,2*λa,x_0effij))-
    x_0ij^(λa+λr)*2*(@f(aS_1,λa+λr)+@f(B,λa+λr,x_0effij))+
    x_0ij^(2*λr)*(@f(aS_1,2*λr)+@f(B,2*λr,x_0effij))+
    x_0ij^(2*λa+2)*2*@f(Q1,λa)*(@f(aS_1,2*λa+2)+@f(B,2*λa+2,x_0effij))*Dij+
    x_0ij^(2*λr+2)*2*@f(Q1,λr)*(@f(aS_1,2*λr+2)+@f(B,2*λr+2,x_0effij))*Dij-
    x_0ij^(λa+λr+2)*2*(@f(Q1,λa)+@f(Q1,λr))*(@f(aS_1,λa+λr+2)+@f(B,λa+λr+2,x_0effij))*Dij+
    x_0ij^(2*λa+4)*@f(Q1,λa)^2*(@f(aS_1,2*λa+4)+@f(B,2*λa+4,x_0effij))*Dij^2+
    x_0ij^(2*λr+4)*@f(Q1,λr)^2*(@f(aS_1,2*λr+4)+@f(B,2*λr+4,x_0effij))*Dij^2-
    x_0ij^(λa+λr+4)*(2*@f(Q1,λa)*@f(Q1,λr))*(@f(aS_1,λa+λr+4)+@f(B,λa+λr+4,x_0effij))*Dij^2+
    x_0ij^(2*λa+4)*2*@f(Q2,λa)*(@f(aS_1,2*λa+4)+@f(B,2*λa+4,x_0effij))*Dij^2+
    x_0ij^(2*λr+4)*2*@f(Q2,λr)*(@f(aS_1,2*λr+4)+@f(B,2*λr+4,x_0effij))*Dij^2-
    x_0ij^(λa+λr+4)*2*(@f(Q2,λa)+@f(Q2,λr))*(@f(aS_1,λa+λr+4)+@f(B,λa+λr+4,x_0effij))*Dij^2+
    x_0ij^(2*λa+6)*2*(@f(Q1,λa)*@f(Q2,λa))*(@f(aS_1,2*λa+6)+@f(B,2*λa+6,x_0effij))*Dij^3+
    x_0ij^(2*λr+6)*2*(@f(Q1,λr)*@f(Q2,λr))*(@f(aS_1,2*λr+6)+@f(B,2*λr+6,x_0effij))*Dij^3-
    x_0ij^(λa+λr+6)*2*(@f(Q1,λr)*@f(Q2,λa)+@f(Q1,λa)*@f(Q2,λr))*(@f(aS_1,λa+λr+6)+@f(B,λa+λr+6,x_0effij))*Dij^3+
    x_0ij^(2*λa+8)*@f(Q2,λa)^2*(@f(aS_1,2*λa+8)+@f(B,2*λa+8,x_0effij))*Dij^4+
    x_0ij^(2*λr+8)*@f(Q2,λr)^2*(@f(aS_1,2*λr+8)+@f(B,2*λr+8,x_0effij))*Dij^4-
    x_0ij^(λa+λr+8)*(2*@f(Q2,λa)*@f(Q2,λr))*(@f(aS_1,λa+λr+8)+@f(B,λa+λr+8,x_0effij))*Dij^4)
end

function χ(model::SAFTVRQMieModel, V, T, z,i,j)
    λr = model.params.lambda_r.values
    λa = model.params.lambda_a.values
    σ = model.params.sigma.values
    ϵ = model.params.epsilon.values
    σeffij = @f(σeff,i,j)
    ϵeffij = @f(ϵeff,i,j)
    Mwij = (model.params.Mw.values[i] + model.params.Mw.values[j])/2 # check
    Dij = ħ^2/(12*k_B*T*Mwij/N_A*σ[i,j]^2)

    ζst_ = ζst(model, V, T, z)
    α = @f(C,i,j)*ϵ[i,j]/ϵeffij*
    (((σ[i,j]/σeffij)^λa[i,j]/(λa[i,j]-3)-(σ[i,j]/σeffij)^λr[i,j]/(λr[i,j]-3))+
        Dij*((σ[i,j]/σeffij)^(2+λa[i,j])*@f(Q1,λa[i,j])/(λa[i,j]-1)-
             (σ[i,j]/σeffij)^(2+λr[i,j])*@f(Q1,λr[i,j])/(λr[i,j]-1))+
        Dij^2*((σ[i,j]/σeffij)^(4+λa[i,j])*@f(Q2,λa[i,j])/(λa[i,j]+1)-
               (σ[i,j]/σeffij)^(4+λr[i,j])*@f(Q2,λr[i,j])/(λr[i,j]+1)))
    return @f(f,α,1)*ζst_+@f(f,α,2)*ζst_^5+@f(f,α,3)*ζst_^8
end

function a_3(model::SAFTVRQMieModel, V, T, z, i, j)
    λr = model.params.lambda_r.values
    λa = model.params.lambda_a.values
    σ = model.params.sigma.values
    ϵ = model.params.epsilon.values
    σeffij = @f(σeff,i,j)
    ϵeffij = @f(ϵeff,i,j)
    Mwij = model.params.Mw.values[i,j]
    Dij = ħ^2/(12*k_B*T*Mwij/N_A*σ[i,j]^2)

    ζst_ = ζst(model, V, T, z)
    α = @f(C,i,j)*ϵ[i,j]/ϵeffij*
    (((σ[i,j]/σeffij)^λa[i,j]/(λa[i,j]-3)-(σ[i,j]/σeffij)^λr[i,j]/(λr[i,j]-3))+
        Dij*((σ[i,j]/σeffij)^(2+λa[i,j])*@f(Q1,λa[i,j])/(λa[i,j]-1)-
             (σ[i,j]/σeffij)^(2+λr[i,j])*@f(Q1,λr[i,j])/(λr[i,j]-1))+
        Dij^2*((σ[i,j]/σeffij)^(4+λa[i,j])*@f(Q2,λa[i,j])/(λa[i,j]+1)-
               (σ[i,j]/σeffij)^(4+λr[i,j])*@f(Q2,λr[i,j])/(λr[i,j]+1)))
    return -ϵeffij^3*@f(f,α,4)*ζst_*exp(@f(f,α,5)*ζst_+@f(f,α,6)*ζst_^2)
end
=#
function a_disp(model::SAFTVRQMieModel, V, T, z,_data = @f(data))
    _σeff,_ϵff,vrdata= _data
    _d,_ρ_S,ζi,_ζ_X,_ζst  = vrdata
    comps = @comps
    l = length(comps)
    ∑z = ∑(z)
    m = model.params.segment.values
    _ϵ = model.params.epsilon.values
    _λr = model.params.lambda_r.values
    _λa = model.params.lambda_a.values
    _σ = model.params.sigma.values
    Mw = model.params.Mw.values
    m̄ = dot(z, m)
    m̄inv = 1/m̄
    a₁ = zero(V+T+first(z))
    a₂ = a₁
    a₃ = a₁
    _ζst5 = _ζst^5
    _ζst8 = _ζst^8
    _KHS = @f(KHS,_ζ_X,_ρ_S)
    for i ∈ comps
        j = i
        x_Si = z[i]*m[i]*m̄inv
        x_Sj = x_Si
        ϵ,λa,λr,σ,di= _ϵ[i,j],_λa[i,i],_λr[i,i],_σ[i,i],_d[i]
        σeff,ϵff = _σeff[i,j],_ϵff[i,j]
        _C = @f(Cλ,λa,λr)
        dij3 = di^3
        x_0ij = σ/di
        x_0effij = σeff/di
        Mwij = Mw[i,j]
        Dij = ħ^2/(12*k_B*T*Mwij/N_A*σ^2)
        Q1λr = @f(Q1,λr)
        Q1λa = @f(Q1,λa)
        Q2λr = @f(Q2,λr)
        Q2λa = @f(Q2,λa)

        #calculations for a1 - diagonal
        a1_ij = 2*π*ϵ*dij3*_C*_ρ_S*
            ( (x_0ij^λa*(@f(aS_1,λa,_ζ_X)+@f(B,λa,x_0effij,_ζ_X))-
            x_0ij^λr*(@f(aS_1,λr,_ζ_X)+@f(B,λr,x_0effij,_ζ_X)))+
            (x_0ij^(λa+2)*Q1λa*(@f(aS_1,λa+2,_ζ_X)+@f(B,λa+2,x_0effij,_ζ_X))-
            x_0ij^(λr+2)*Q1λr*(@f(aS_1,λr+2,_ζ_X)+@f(B,λr+2,x_0effij,_ζ_X)))*Dij+
            (x_0ij^(λa+4)*Q2λa*(@f(aS_1,λa+4,_ζ_X)+@f(B,λa+4,x_0effij,_ζ_X))-
            x_0ij^(λr+4)*Q2λr*(@f(aS_1,λr+4,_ζ_X)+@f(B,λr+4,x_0effij,_ζ_X)))*Dij^2 )
        #calculations for a2 - diagonal
        σcoeff =σ/σeff
        α = _C*ϵ/ϵff*
            ((σcoeff^λa/(λa-3)-σcoeff^λr/(λr-3))+
            Dij*(σcoeff^(2+λa)*Q1λa/(λa-1)-
            σcoeff^(2+λr)*Q1λr/(λr-1))+
            Dij^2*(σcoeff^(4+λa)*Q2λa/(λa+1)-
            σcoeff^(4+λr)*Q2λr/(λr+1)))
        f1,f2,f3,f4,f5,f6 = @f(f123456,α)
         _χ = f1*_ζst+f2*_ζst5+f3*_ζst8
        a2_ij = π*_KHS*(1+_χ)*_ρ_S*ϵ^2*dij3*_C^2*(x_0ij^(2*λa)*(@f(aS_1,2*λa)+@f(B,2*λa,x_0effij,_ζ_X))-
            x_0ij^(λa+λr)*2*(@f(aS_1,λa+λr)+@f(B,λa+λr,x_0effij,_ζ_X))+
            x_0ij^(2*λr)*(@f(aS_1,2*λr)+@f(B,2*λr,x_0effij,_ζ_X))+
            x_0ij^(2*λa+2)*2*Q1λa*(@f(aS_1,2*λa+2)+@f(B,2*λa+2,x_0effij,_ζ_X))*Dij+
            x_0ij^(2*λr+2)*2*Q1λr*(@f(aS_1,2*λr+2)+@f(B,2*λr+2,x_0effij,_ζ_X))*Dij-
            x_0ij^(λa+λr+2)*2*(Q1λa+Q1λr)*(@f(aS_1,λa+λr+2)+@f(B,λa+λr+2,x_0effij,_ζ_X))*Dij+
            x_0ij^(2*λa+4)*Q1λa^2*(@f(aS_1,2*λa+4,_ζ_X)+@f(B,2*λa+4,x_0effij,_ζ_X))*Dij^2+
            x_0ij^(2*λr+4)*Q1λr^2*(@f(aS_1,2*λr+4,_ζ_X)+@f(B,2*λr+4,x_0effij,_ζ_X))*Dij^2-
            x_0ij^(λa+λr+4)*(2*Q1λa*Q1λr)*(@f(aS_1,λa+λr+4,_ζ_X)+@f(B,λa+λr+4,x_0effij,_ζ_X))*Dij^2+
            x_0ij^(2*λa+4)*2*Q2λa*(@f(aS_1,2*λa+4,_ζ_X)+@f(B,2*λa+4,x_0effij,_ζ_X))*Dij^2+
            x_0ij^(2*λr+4)*2*Q2λr*(@f(aS_1,2*λr+4,_ζ_X)+@f(B,2*λr+4,x_0effij,_ζ_X))*Dij^2-
            x_0ij^(λa+λr+4)*2*(Q2λa+Q2λr)*(@f(aS_1,λa+λr+4,_ζ_X)+@f(B,λa+λr+4,x_0effij,_ζ_X))*Dij^2+
            x_0ij^(2*λa+6)*2*(Q1λa*Q2λa)*(@f(aS_1,2*λa+6)+@f(B,2*λa+6,x_0effij,_ζ_X))*Dij^3+
            x_0ij^(2*λr+6)*2*(Q1λr*Q2λr)*(@f(aS_1,2*λr+6)+@f(B,2*λr+6,x_0effij,_ζ_X))*Dij^3-
            x_0ij^(λa+λr+6)*2*(Q1λr*Q2λa+Q1λa*Q2λr)*(@f(aS_1,λa+λr+6)+@f(B,λa+λr+6,x_0effij,_ζ_X))*Dij^3+
            x_0ij^(2*λa+8)*Q2λa^2*(@f(aS_1,2*λa+8,_ζ_X)+@f(B,2*λa+8,x_0effij,_ζ_X))*Dij^4+
            x_0ij^(2*λr+8)*Q2λr^2*(@f(aS_1,2*λr+8,_ζ_X)+@f(B,2*λr+8,x_0effij,_ζ_X))*Dij^4-
            x_0ij^(λa+λr+8)*(2*Q2λa*Q2λr)*(@f(aS_1,λa+λr+8,_ζ_X)+@f(B,λa+λr+8,x_0effij,_ζ_X))*Dij^4)

        #calculations for a3 - diagonal
        a3_ij = -ϵff^3*f4*_ζst * exp(f5*_ζst+f6*_ζst^2)
        #adding - diagonal
        a₁ += a1_ij*x_Si*x_Si
        a₂ += a2_ij*x_Si*x_Si
        a₃ += a3_ij*x_Si*x_Si
        for j ∈ (i+1):l
            x_Sj = z[j]*m[j]*m̄inv
            ϵ,λa,λr,σ,dj= _ϵ[i,j],_λa[i,j],_λr[i,j],_σ[i,j],_d[j]
            σeff,ϵff = _σeff[i,j],_ϵff[i,j]
            _C = @f(Cλ,λa,λr)
            dij = 0.5*(di+dj)
            dij3 = dij^3
            x_0ij = σ/dij
            x_0effij = σeff/dij
            Q1λr = @f(Q1,λr)
            Q1λa = @f(Q1,λa)
            Q2λr = @f(Q2,λr)
            Q2λa = @f(Q2,λa)
            Mwij = Mw[i,j]
            Dij = ħ^2/(12*k_B*T*Mwij/N_A*σ^2)

            #calculations for a1
            a1_ij = 2*π*ϵ*dij3*_C*ρ_S*
            ( (x_0ij^λa*(@f(aS_1,λa,_ζ_X)+@f(B,λa,x_0effij,_ζ_X))-
                x_0ij^λr*(@f(aS_1,λr,_ζ_X)+@f(B,λr,x_0effij,_ζ_X)))+
                (x_0ij^(λa+2)*Q1λa*(@f(aS_1,λa+2,_ζ_X)+@f(B,λa+2,x_0effij,_ζ_X))-
                 x_0ij^(λr+2)*Q1λr*(@f(aS_1,λr+2,_ζ_X)+@f(B,λr+2,x_0effij,_ζ_X)))*Dij+
                (x_0ij^(λa+4)*Q2λa*(@f(aS_1,λa+4,_ζ_X)+@f(B,λa+4,x_0effij,_ζ_X))-
                 x_0ij^(λr+4)*Q2λr*(@f(aS_1,λr+4,_ζ_X)+@f(B,λr+4,x_0effij,_ζ_X)))*Dij^2 )

            #calculations for a2
            σcoeff =σ/σeff
            α = _C*ϵ/ϵff*
                ((σcoeff^λa/(λa-3)-σcoeff^λr/(λr-3))+
                Dij*(σcoeff^(2+λa)*Q1λa/(λa-1)-
                σcoeff^(2+λr)*Q1λr/(λr-1))+
                Dij^2*(σcoeff^(4+λa)*Q2λa/(λa+1)-
                σcoeff^(4+λr)*Q2λr/(λr+1)))
            f1,f2,f3,f4,f5,f6 = @f(f123456,α)
             _χ = f1*_ζst+f2*_ζst5+f3*_ζst8
             a2_ij = π*_KHS*(1+_χ)*_ρ_S*ϵ^2*dij3*_C^2*(x_0ij^(2*λa)*(@f(aS_1,2*λa)+@f(B,2*λa,x_0effij,_ζ_X))-
             x_0ij^(λa+λr)*2*(@f(aS_1,λa+λr)+@f(B,λa+λr,x_0effij,_ζ_X))+
             x_0ij^(2*λr)*(@f(aS_1,2*λr)+@f(B,2*λr,x_0effij,_ζ_X))+
             x_0ij^(2*λa+2)*2*Q1λa*(@f(aS_1,2*λa+2)+@f(B,2*λa+2,x_0effij,_ζ_X))*Dij+
             x_0ij^(2*λr+2)*2*Q1λr*(@f(aS_1,2*λr+2)+@f(B,2*λr+2,x_0effij,_ζ_X))*Dij-
             x_0ij^(λa+λr+2)*2*(Q1λa+Q1λr)*(@f(aS_1,λa+λr+2)+@f(B,λa+λr+2,x_0effij,_ζ_X))*Dij+
             x_0ij^(2*λa+4)*Q1λa^2*(@f(aS_1,2*λa+4,_ζ_X)+@f(B,2*λa+4,x_0effij,_ζ_X))*Dij^2+
             x_0ij^(2*λr+4)*Q1λr^2*(@f(aS_1,2*λr+4,_ζ_X)+@f(B,2*λr+4,x_0effij,_ζ_X))*Dij^2-
             x_0ij^(λa+λr+4)*(2*Q1λa*Q1λr)*(@f(aS_1,λa+λr+4,_ζ_X)+@f(B,λa+λr+4,x_0effij,_ζ_X))*Dij^2+
             x_0ij^(2*λa+4)*2*Q2λa*(@f(aS_1,2*λa+4,_ζ_X)+@f(B,2*λa+4,x_0effij,_ζ_X))*Dij^2+
             x_0ij^(2*λr+4)*2*Q2λr*(@f(aS_1,2*λr+4,_ζ_X)+@f(B,2*λr+4,x_0effij,_ζ_X))*Dij^2-
             x_0ij^(λa+λr+4)*2*(Q2λa+Q2λr)*(@f(aS_1,λa+λr+4,_ζ_X)+@f(B,λa+λr+4,x_0effij,_ζ_X))*Dij^2+
             x_0ij^(2*λa+6)*2*(Q1λa*Q2λa)*(@f(aS_1,2*λa+6)+@f(B,2*λa+6,x_0effij,_ζ_X))*Dij^3+
             x_0ij^(2*λr+6)*2*(Q1λr*Q2λr)*(@f(aS_1,2*λr+6)+@f(B,2*λr+6,x_0effij,_ζ_X))*Dij^3-
             x_0ij^(λa+λr+6)*2*(Q1λr*Q2λa+Q1λa*Q2λr)*(@f(aS_1,λa+λr+6)+@f(B,λa+λr+6,x_0effij,_ζ_X))*Dij^3+
             x_0ij^(2*λa+8)*Q2λa^2*(@f(aS_1,2*λa+8,_ζ_X)+@f(B,2*λa+8,x_0effij,_ζ_X))*Dij^4+
             x_0ij^(2*λr+8)*Q2λr^2*(@f(aS_1,2*λr+8,_ζ_X)+@f(B,2*λr+8,x_0effij,_ζ_X))*Dij^4-
             x_0ij^(λa+λr+8)*(2*Q2λa*Q2λr)*(@f(aS_1,λa+λr+8,_ζ_X)+@f(B,λa+λr+8,x_0effij,_ζ_X))*Dij^4)

            #calculations for a3
            a3_ij = -ϵff^3*f4*_ζst * exp(f5*_ζst+f6*_ζst^2)
            #adding
            a₁ += 2*a1_ij*x_Si*x_Sj
            a₂ += 2*a2_ij*x_Si*x_Sj
            a₃ += 2*a3_ij*x_Si*x_Sj
        end
    end
    a₁ = a₁*m̄/T/∑z
    a₂ = a₂*m̄/(T*T)/∑z
    a₃ = a₃*m̄/(T*T*T)/∑z
    adisp =  a₁ + a₂ + a₃
    return adisp
end

const SAFTVRQMieconsts = (
    x = [-0.9956571630258080807355,-0.973906528517171720078,-0.9301574913557082260012,-0.8650633666889845107321,-0.7808177265864168970637,-0.6794095682990244062343,-0.562757134668604683339,
        -0.4333953941292471907993,-0.294392862701460198131,-0.1488743389816312108848,0,0.1488743389816312108848,0.2943928627014601981311,0.4333953941292471907993,
        0.562757134668604683339,0.6794095682990244062343,0.7808177265864168970637,0.865063366688984510732,0.9301574913557082260012,0.973906528517171720078,0.9956571630258080807355],
    w = [0.0116946388673718742781,0.0325581623079647274788,0.0547558965743519960314,0.075039674810919952767,0.093125454583697605535,0.1093871588022976418992,0.123491976262065851078,
        0.134709217311473325928,0.142775938577060080797,0.1477391049013384913748,0.149445554002916905665,0.1477391049013384913748,0.1427759385770600807971,0.134709217311473325928,
        0.123491976262065851078,0.109387158802297641899,0.093125454583697605535,0.075039674810919952767,0.05475589657435199603138,0.032558162307964727479,0.0116946388673718742781],
)


text:
# Should LoopVectorization.lower and ReverseDiffExpressions.lower be distinct?
#import LoopVectorization: lower, lower!
# Because they mean the same thing -- turn obj into Expr -- I'm importing LoopVectorization.lower for now.

# Algorithm for lowering:
# 1. Pick an unpicked func.
# 2. Check if all vars are initialized. If so, go to step 4.
# 3. For all that are not, check if they have a parentfunc. Error if not, if so go to step 1 (picking that parentfunc); when done, resume here.
# 4. Lower func, initializing any vars it defines.
# 5. Return up callstack to 3, or iterate to next unpicked function.
#
# Eventually, this should be improved by picking the order more intelligently.

lower(m::Model) = lower!(Expr(:block), m)
function lower!(q::Expr, m::Model)
    @unpack vars, targets = m
    reset_funclowered!(m)
    reset_varinitialized!(m)
    foreach(t -> lower!(q, m, vars[t]), Iterators.reverse(m.targets))
    push!(q.args, Expr(:call, (:vsum), Symbol("##TARGET##")))
    q
end

function asexpr(instr::LoopVectorization.Instruction)
    if instr.mod === Symbol("")
        instr.instr
    else
        Expr(:(.), instr.mod, QuoteNode(instr.instr))
    end
end

function getprop_call!(q::Expr, m::Model, func::Func, retvar::Variable)
    @unpack vars = m
    vparents = parents(func)
    call = Expr(:(.), name(vars[vparents[1]]), QuoteNode(vars[vparents[2]].ref))
    lowernum = num_lowered(retvar)
    retname = return_name(retvar, lowernum)
    if iszero(lowernum)
        push!(q.args, Expr(:(=), retname, call))
    else
        retnametemp = gensym(retname)
        push!(q.args, Expr(:(=), retnametemp, call))
        push!(q.args, Expr(:(=), retname, Expr(:call, :(ReverseDiffExpressions.vadd!), return_name(retvar, lowernum - 1), retnametemp)))
    end
    func.lowered[] = true
    retvar.lowered_count += 1
    nothing
end
function indfunc_call!(q::Expr, m::Model, func::Func, retvar::Variable, check_paireddeps::Bool = true)
    @unpack vars, funcs = m
    vparents = parents(func)
    call = Expr(:call, asexpr(func.instr))
    for vpid ∈ parents(func)
        p = vars[vpid]
        isinitialized(p) || lower!(q, m, p, check_paireddeps)
        push!(call.args, p)
    end
    # indfuncs may get called repeatedly as a result of trying to lower others,
    # due to the immediate lowering performed by constrain_pullbacks
    # so once parents are lowered, we may find that this func has
    # been lowered in that process.
    # Therefore, check and return if so.
    if islowered(func)
        # however, this means some paired deps were skipped
        lower_paired_deps!(q, m, retvar)#vars[first(parents(func))])
        return
    end
    lowernum = num_lowered(retvar)
    retname = return_name(retvar, lowernum)
    if iszero(lowernum)
        push!(q.args, Expr(:(=), retname, call))
    else
        retnametemp = gensym(retname)
        push!(q.args, Expr(:(=), retnametemp, call))
        push!(q.args, Expr(:(=), retname, Expr(:call, :(ReverseDiffExpressions.vadd!), return_name(retvar, lowernum - 1), retnametemp)))
    end
    func.lowered[] = true
    retvar.lowered_count += 1
    nothing
end
function lower_paired_deps!(q::Expr, m::Model, v::Variable)
    for pv ∈ v.paireddeps
        lower!(q, m, m.vars[pv])
    end
end
function lower!(q::Expr, m::Model, v::Variable, check_paireddeps::Bool = true)
    isinitialized(v) && return
    @unpack vars, funcs = m
    # @assert name(v) !== Symbol("L##BAR##")
    vparents = parents(v)
    # check if any constrained
    if iszero(num_lowered(v))
        constrainid = findfirst(pf -> isconstrainttransform(funcs[pf]), vparents)
        skipid = if isnothing(constrainid)
            notupdateingid = findfirst(vparents) do pf
                func = nonindex_parent(m, funcs[pf])
                !func.probdistapi && notupdating(func)
            end
            something(notupdateingid, 1)
        else
            constrainid
        end
        skippf = vparents[skipid]
        lower!(q, m, funcs[skippf], v)
    end
    length(vparents) > 1 && foreach(pf -> lower!(q, m, funcs[pf], v), vparents)
    if check_paireddeps
        lower_paired_deps!(q, m, v)
    end
    # for pf ∈ v.useids
    #     func = funcs[pf]
    #     isindexfunc(func) && lower!(q, m, vars[func.output[]])
    # end
    nothing
end
function return_name(retvar::Variable, lowernum::Int = num_lowered(retvar))
    lowernum + 1 == length(parents(retvar)) ? name(retvar) : Symbol(name(retvar), '#', lowernum, '#')
end
function probdist_call(m::Model, func::Func, retvar::Variable)
    funcupdating = m.gradmodel
    f = funcupdating ? :∂logdensity! : :logdensity
    @unpack vars, funcs = m
    disttup = Expr(:tuple); append!(disttup.args, map(id -> istracked(vars[id]), parents(func)))
    dist = Expr(:call, Expr(:curly, func.instr.instr, disttup))
    call = Expr(:call, Expr(:(.), :ReverseDiffExpressions, QuoteNode(:stack_pointer_call)), f, STACK_POINTER_NAME)
    push!(call.args, dist)
    call
end

function add_probdist_partials!(q::Expr, call::Expr, m::Model, func::Func, retvar::Variable)
    @unpack vars, funcs = m
    ∂tup = Expr(:tuple)
    insert!(call.args, 4, ∂tup)
    for id ∈ parents(func)
        pvar = vars[id]
        if !istracked(pvar)
            push!(∂tup.args, nothing)
            continue
        end
        ∂pvar = getvar!(m, diffsym(pvar))
        lowernum = num_lowered(∂pvar)
        # It isn't fully initialized, but perhaps it isn't partially initialized either.
        if iszero(lowernum)
            # If it hasn't been initialized, we look if one of them is simply missing
            # a needed getindex call to actually extract it.
            found = false
            for pfuncid ∈ parents(∂pvar)
                immediate_parent_func = funcs[pfuncid]
                pfunc = nonindex_parent(m, immediate_parent_func)
                if islowered(pfunc) # the actual parent has been lowered: (∂pvar,) = pfunc(...)
                    found = true
                    # so we extract it: ∂pvar = (∂pvar,)[1]
                    lower!(q, m, immediate_parent_func, ∂pvar)
                    break
                end
            end
            if found
                push!(∂tup.args, return_name(∂pvar, 0))
            else
                push!(∂tup.args, nothing)
            end
        else
            push!(∂tup.args, return_name(∂pvar, lowernum - 1))
        end
    end
end
function lower!(q::Expr, m::Model, func::Func, retvar::Variable, check_paireddeps::Bool = true)
    func.lowered[] && return
    iszero(func.loopsetid) || return lower_loopset!(q, func, m)
    @unpack vars, funcs = m
    if func.probdistapi
        funcupdating = m.gradmodel
        call = probdist_call(m, func, retvar)
    elseif func.instr === Instruction(:Base, :getproperty)
        return getprop_call!(q, m, func, retvar)
    elseif isindexfunc(func)
        return indfunc_call!(q, m, func, retvar, check_paireddeps)
    else
        lowernum = num_lowered(retvar)
        funcupdating = !(iszero(lowernum) || notupdating(func))
        instr = funcupdating ? ReverseDiffExpressionsBase.MAKEUPDATING[func] : func.instr
        if instr === Instruction(Symbol(""),:constrain_reverse!)
            lower_paired_deps!(q, m, vars[first(parents(func))])
        end
        call = Expr(:call, Expr(:(.), :ReverseDiffExpressions, QuoteNode(:stack_pointer_call)), asexpr(instr), STACK_POINTER_NAME)
        funcupdating && push!(call.args, return_name(retvar, lowernum - 1))
    end
    for vpid ∈ parents(func)
        p = vars[vpid]
        lower!(q, m, p)
        push!(call.args, p)
    end
    if func.probdistapi && funcupdating
        add_probdist_partials!(q, call, m, func, retvar)
    end

    retvarid = func.output[]
    ret = Expr(:tuple, STACK_POINTER_NAME)
    call = Expr(:(=), ret, call)
    push!(q.args, call)
    if retvarid ≥ 0
        lowernum = num_lowered(retvar)
        retname = return_name(retvar, lowernum)
        if iszero(lowernum) | funcupdating
            push!(ret.args, retname)
        else
            retnametemp = gensym(retname)
            push!(ret.args, retnametemp)
            push!(q.args, Expr(:(=), retname, Expr(:call, :(ReverseDiffExpressions.vadd!), return_name(retvar, lowernum - 1), retnametemp)))
        end
    end
    func.lowered[] = true
    retvar.lowered_count += 1
    if func.instr.instr === :constrain_pullback!
        let fid1 = retvar.useids[2]
            let fid2 = (vars[funcs[fid1].output[]]).useids[2]
                func2 = funcs[fid2]
                varout = vars[func2.output[]]
                lower!(q, m, func2, varout, false)
                lower_paired_deps!(q, m, varout)
            end
        end
    elseif func.probdistapi && funcupdating

        # Need to force lowering of
    end
    nothing
end

function lower_loopset!(q::Expr, func::Func, m::Model)
    push!(q.args, setup_call(m.loops[func.loopsetid]))
end


text:
export tryuparse, dict_to_splatting, params_from_dict

function tryuparse(val::Any)
  try
    return uparse(val)
  catch e
    return val
  end
end

function tryuparse(val::Vector{T}) where {T}
  try
    return tryuparse.(val)
  catch e
    return val
  end
end

function dict_to_splatting(dict::Dict)
  splattingDict = Dict{Symbol, Any}()
  for (key, value) in dict
    if value isa Dict # Do we need recursion here?
      specializedType = nothing
      doSpecialize = true
      for (subkey, subvalue) in value
        value[subkey] = tryuparse.(subvalue) # Convert with Unitful if applicable

        if isnothing(specializedType)
          specializedType = typeof(value[subkey])
        elseif typeof(value[subkey]) != specializedType
          doSpecialize = false
        end
      end

      # Check if the types are equal and we can convert the dict to a more specialized version
      if doSpecialize
        value = convert(Dict{String, specializedType}, value)
      end
    else
      value = tryuparse.(value) # Convert with Unitful if applicable
    end
    splattingDict[Symbol(key)] = value
  end

  return splattingDict
end

function params_from_dict(type::DataType, dict::Dict)
  splattingDict = dict_to_splatting(dict)

  try
    return type(;splattingDict...)
  catch e
    if e isa UndefKeywordError
      throw(ScannerConfigurationError("The required field `$(e.var)` is missing in your configuration for a device with the params type `$type`."))
    else
      rethrow()
    end
  end
end

function stringToEnum(value::AbstractString, enumType::Type{T}) where {T <: Enum}
  stringInstances = string.(instances(enumType))
  # If lowercase is not sufficient one could try Unicode.normalize with casefolding
  index = findfirst(isequal(lowercase(value)), lowercase.(stringInstances))
  if isnothing(index)
    throw(ArgumentError("$value cannot be resolved to an instance of $(typeof(enumType)). Possible instances are: " * join(stringInstances, ", ", " and ")))
  end
  return instances(enumType)[index]
end


text:
mutable struct XYPlot <: ESDLPlot
  xaxis
  group
end
plotAxVars(p::XYPlot)=[FixedAx(p.xaxis,"X Axis",true,true,1),FixedAx(p.group,"Group",false,false,2)]

function plotCall(::XYPlot,d,ixaxis,igroup,otherinds...)

  axlist = caxes(d)
  inds1 = ntuple(  i->in(i,(ixaxis,igroup)) ? (:) : axVal2Index(axlist[i],otherinds[i],fuzzy=true), length(otherinds))
  x1 = d[inds1...]

  all(ismissing(x1)) && (x1=fill(NaN,length(x1)))
  replace!(x1,missing=>NaN)
  x1 = Float64.(x1)

  if igroup > 0
    igroup < ixaxis && (x1=transpose(x1))
    if isa(axlist[ixaxis],CategoricalAxis)
      plotf = groupedbar
    else
      plotf = plot
    end
    labs = reshape(string.(axlist[igroup].values),(1,length(axlist[igroup])))
    xlabel = axname(axlist[ixaxis])
    p=plotf(axlist[ixaxis].values,x1,
    lab=labs,
    xlabel=xlabel)
  else
    plotf = isa(axlist[ixaxis],CategoricalAxis) ? bar : plot
    p=plotf(axlist[ixaxis].values,x1,xlabel=axname(axlist[ixaxis]),fmt=:png)
  end
  p
end
nplotCubes(::XYPlot)=1

"""
    plotXY(cube; group=0, xaxis=-1, kwargs...)

Generic plotting tool for cube objects, can be called on any type of cube data.

### Keyword arguments

* `xaxis` which axis is to be used as x axis. Can be either an axis Datatype or a string. Short versions of axes names are possible as long as the axis can be uniquely determined.
* `group` it is possible to group the plot by a categorical axis. Can be either an axis data type or a string.
* `dim=value` can set other dimensions to certain values, for example `lon=51.5` will fix the longitude for the resulting plot

If a dimension is not the x axis or group variable and is not fixed through an additional keyword, a slider or dropdown menu will appear to select the axis value.
"""
function plotXY(cube;group=nothing,xaxis=nothing,kwargs...)

  return plotGeneric(XYPlot(xaxis,group),cube;kwargs...)

end

mutable struct ScatterPlot <: ESDLPlot
  vsaxis
  alongaxis
  group
  c_1
  c_2
end
plotAxVars(p::ScatterPlot)=[
  FixedAx(p.vsaxis,"VS Axis",true,true,-1),
  FixedAx(p.alongaxis,"Along",true,false,1),
  FixedAx(p.group,"Group",false,false,2),
  FixedVar(p.vsaxis,p.c_1,"X Axis",true),
  FixedVar(p.vsaxis,p.c_2,"Y Axis",true)
  ]

function plotCall(::ScatterPlot,d,ivsaxis,ialongaxis,igroup,c1,c2,otherinds...)

  axlist = caxes(d)
  inds1 = ntuple(  i->(i==ivsaxis)    ? axVal2Index(axlist[i],c1,fuzzy=true)  : (i==ialongaxis) ? (:) : (i==igroup)     ? (:) : axVal2Index(axlist[i],otherinds[i],fuzzy=true), length(otherinds))
  inds2 = ntuple(  i->(i==ivsaxis)    ? axVal2Index(axlist[i],c2,fuzzy=true)  : (i==ialongaxis) ? (:) : (i==igroup)     ? (:) : axVal2Index(axlist[i],otherinds[i],fuzzy=true), length(otherinds))
  x1 = d[inds1...]
  x2 = d[inds2...]

  all(ismissing(x1)) && (x1=fill(NaN,length(x1)))
  all(ismissing(x2)) && (x2=fill(NaN,length(x2)))
  replace!(x1,missing=>NaN)
  replace!(x2,missing=>NaN)

  goodinds = map((m1,m2)->!ismissing(m1) && !ismissing(m2),x1,x2)
  a_1 = x1[goodinds]
  a_2 = x2[goodinds]
  if isempty(a_1)
    push!(a_1,1);push!(a_2,1)
  end
  nPoints = length(a_1)
  pointSize = min(5000/nPoints,3)
  msw=pointSize > 2 ? 1 : 0
  fmt=nPoints>20000 ? :png : :svg
  if igroup > 0 && (igroup != ialongaxis)
    plotf = isa(axlist[ivsaxis],CategoricalAxis) ? groupedbar : plot
    igroup < ialongaxis && (a_1=transpose(a_1);a_2=transpose(a_2))
    p=plotf(a_1,a_2,
      xlabel=string(axlist[ivsaxis].values[inds1[ivsaxis]]),
      ylabel=string(axlist[ivsaxis].values[inds2[ivsaxis]]),
      lab=reshape(string.(axlist[igroup].values),(1,length(axlist[igroup]))),
      fmt=fmt,
      ms=pointSize,
      markerstrokewidth=msw,
      )
  else
    p = scatter(a_1,a_2,
      xlabel=string(axlist[ivsaxis].values[inds1[ivsaxis]]),
      ylabel=string(axlist[ivsaxis].values[inds2[ivsaxis]]),
      lab="",
      fmt=fmt,
      ms=pointSize,
      markerstrokewidth=msw,
      )
  end
  p
end
"""
    plotScatter(cube; vsaxis=VariableAxis, alongaxis=0, group=0, xaxis=0, yaxis=0, kwargs...)

Generic plotting tool for cube objects to generate scatter plots, like variable `A` against variable `B`. Can be called on any type of cube data.

### Keyword arguments

* `vsaxis` determines the axis from which the x and y variables are drawn.
* `alongaxis` determines the axis along which the variables are plotted. E.g. if you choose `TimeAxis`, a dot will be plotted for each time step.
* `xaxis` index or value of the variable to plot on the x axis
* `yaxis` index or value of the variable to plot on the y axis
* `group` it is possible to group the plot by an axis. Can be either an axis data type or a string. *Caution: This will increase the number of plotted data points*
* `dim=value` can set other dimensions to certain values, for example `lon=51.5` will fix the longitude for the resulting plot

If a dimension is not the `vsaxis` or `alongaxis` or `group` and is not fixed through an additional keyword, a slider or dropdown menu will appear to select the axis value.
"""
function plotScatter(cube;group=nothing,vsaxis="Variable",xaxis=nothing,yaxis=nothing,alongaxis=nothing,kwargs...)

  return plotGeneric(ScatterPlot(vsaxis,alongaxis,group,xaxis,yaxis),cube;kwargs...)

end

export plotHist
plotHist(c;kwargs...)=plotScatter(c;vsaxis="Hist",xaxis="MidPoints",yaxis="Frequency",alongaxis="Bin",kwargs...)


text:
function momentforbox(mass::Real, width::Real, height::Real)
    ccall(dlsym(libchipmunk, :cpMomentForBox), Cdouble, (Cdouble, Cdouble, Cdouble,), mass, width, height)
end

export momentforbox


text:
using CLIMAParameters.Planet: Omega
export Source, Gravity, RayleighSponge, Subsidence, GeostrophicForcing, Coriolis

# kept for compatibility
# can be removed if no functions are using this
function atmos_source!(
    f::Function,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    f(atmos, source, state, diffusive, aux, t, direction)
end
function atmos_source!(
    ::Nothing,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
) end
# sources are applied additively
@generated function atmos_source!(
    stuple::Tuple,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    N = fieldcount(stuple)
    return quote
        Base.Cartesian.@nexprs $N i -> atmos_source!(
            stuple[i],
            atmos,
            source,
            state,
            diffusive,
            aux,
            t,
            direction,
        )
        return nothing
    end
end

abstract type Source end

struct Gravity <: Source end
function atmos_source!(
    ::Gravity,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    if atmos.ref_state isa HydrostaticState
        source.ρu -= (state.ρ - aux.ref_state.ρ) * aux.orientation.∇Φ
    else
        source.ρu -= state.ρ * aux.orientation.∇Φ
    end
end

struct Coriolis <: Source end
function atmos_source!(
    ::Coriolis,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    FT = eltype(state)
    _Omega::FT = Omega(atmos.param_set)
    # note: this assumes a SphericalOrientation
    source.ρu -= SVector(0, 0, 2 * _Omega) × state.ρu
end

struct Subsidence{FT} <: Source
    D::FT
end

function atmos_source!(
    subsidence::Subsidence,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    ρ = state.ρ
    z = altitude(atmos, aux)
    w_sub = subsidence_velocity(subsidence, z)
    k̂ = vertical_unit_vector(atmos, aux)

    source.ρe -= ρ * w_sub * dot(k̂, diffusive.∇h_tot)
    source.moisture.ρq_tot -= ρ * w_sub * dot(k̂, diffusive.moisture.∇q_tot)
end

subsidence_velocity(subsidence::Subsidence{FT}, z::FT) where {FT} =
    -subsidence.D * z

struct GeostrophicForcing{FT} <: Source
    f_coriolis::FT
    u_geostrophic::FT
    v_geostrophic::FT
end
function atmos_source!(
    s::GeostrophicForcing,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    u_geo = SVector(s.u_geostrophic, s.v_geostrophic, 0)
    ẑ = vertical_unit_vector(atmos, aux)
    fkvector = s.f_coriolis * ẑ
    source.ρu -= fkvector × (state.ρu .- state.ρ * u_geo)
end

"""
    RayleighSponge{FT} <: Source

Rayleigh Damping (Linear Relaxation) for top wall momentum components
Assumes laterally periodic boundary conditions for LES flows. Momentum components
are relaxed to reference values (zero velocities) at the top boundary.
"""
struct RayleighSponge{FT} <: Source
    "Maximum domain altitude (m)"
    z_max::FT
    "Altitude at with sponge starts (m)"
    z_sponge::FT
    "Sponge Strength 0 ⩽ α_max ⩽ 1"
    α_max::FT
    "Relaxation velocity components"
    u_relaxation::SVector{3, FT}
    "Sponge exponent"
    γ::FT
end
function atmos_source!(
    s::RayleighSponge,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    z = altitude(atmos, aux)
    if z >= s.z_sponge
        r = (z - s.z_sponge) / (s.z_max - s.z_sponge)
        β_sponge = s.α_max * sinpi(r / 2)^s.γ
        source.ρu -= β_sponge * (state.ρu .- state.ρ * s.u_relaxation)
    end
end


text:
using Query
using DataFrames
using Statistics

df = DataFrame(name=repeat(["John", "Sally", "Kirk"],inner=[1],outer=[2]),
     age=vcat([10., 20., 30.],[10., 20., 30.].+3),
     children=repeat([3,2,2],inner=[1],outer=[2]),state=[:a,:a,:a,:b,:b,:b])

x = @from i in df begin
    @group i by i.state into g
    @select {group=key(g),mage=mean(g.age), oldest=maximum(g.age), youngest=minimum(g.age)}
    @collect DataFrame
end

println(x)


text:
module ExprTools

export args_tuple_expr, combinedef, parameters, signature, splitdef

include("function.jl")
include("method.jl")
include("type_utils.jl")
include("def_tools.jl")

end  # module


text:
# # Multinomial logistic regression

# Import the necessary packages:

using DataFrames
using Distributions
using MLJBase
using NNlib
using RDatasets
using Soss
using SossMLJ
using Statistics

# In this example, we fit a Bayesian multinomial logistic regression model
# with the canonical link function.

# Suppose that we are given a matrix of features `X` and a column vector of
# labels `y`. `X` has `n` rows and `p` columns. `y` has `n` elements. We assume
# that our observation vector `y` is a realization of a random variable `Y`.
# We define `μ` (mu) as the expected value of `Y`, i.e. `μ := E[Y]`. Our model
# comprises three components:
#
# 1. The probability distribution of `Y`. We assume that each `Yᵢ` follows a multinomial distribution with `k` categories, mean `μᵢ`, and one trial.
# 2. The systematic component, which consists of linear predictor `η` (eta), which we define as `η := Xβ`, where `β` is the column vector of `p` coefficients.
# 3. The link function `g`, which provides the following relationship: `g(E[Y]) = g(μ) = η = Xβ`. It follows that `μ = g⁻¹(η)`, where `g⁻¹` denotes the inverse of `g`. In multinomial logistic regression, the canonical link function is the generalized logit function. The inverse of the generalized logit function is the softmax function. Therefore, when using the canonical link function, `μ = g⁻¹(η) = softmax(η)`.

# A multinomial distribution with one trial is equivalent to the categorical distribution.
# Therefore, the following two statements are equivalent:
# - `Yᵢ` follows a multinomial distribution with `k` categories, mean `μᵢ`, and one trial.
# - `Yᵢ` follows a categorical distribution with `k` categories and mean `μᵢ`.
#
# Observe that the logistic regression model is a special case of the
# multinomial logistic regression model where `k = 2`.
#
# In this model, the parameters that we want to estimate are the coefficients `β`.
# We need to select prior distributions for these parameters. For each `βᵢ`
# we choose a normal distribution with zero mean and unit variance. Here, `βᵢ`
# denotes the `i`th component of `β`.

# We define this model using the Soss probabilistic programming library:

m = @model X,pool begin
    n = size(X, 1) # number of observations
    p = size(X, 2) # number of features
    k = length(pool.levels) # number of classes
    β ~ Normal(0.0, 1.0) |> iid(p, k) # coefficients
    η = X * β # linear predictor
    μ = NNlib.softmax(η; dims=2) # μ = g⁻¹(η) = softmax(η)
    y_dists = UnivariateFinite(pool.levels, μ; pool=pool) # `UnivariateFinite` is mathematically equivalent to `Categorical`
    y ~ For(j -> y_dists[j], n) # `Yᵢ ~ UnivariateFinite(mean=μᵢ, categories=k)`, which is mathematically equivalent to `Yᵢ ~ Categorical(mean=μᵢ, categories=k)`
end;

# Import the Iris flower data set:

iris = dataset("datasets", "iris");

# Define our feature columns:

feature_columns = [
    :PetalLength,
    :PetalWidth,
    :SepalLength,
    :SepalWidth,
]

# Define our label column:

label_column = :Species

# Convert the Soss model into a `SossMLJModel`:

model = SossMLJModel(;
    model       = m,
    predictor   = MLJBase.UnivariateFinite,
    hyperparams = (pool=iris.Species.pool,),
    infer       = dynamicHMC,
    response    = :y,
);

# Create an MLJ machine for fitting our model:

mach = MLJBase.machine(model, iris[!, feature_columns], iris[!, :Species])

# Fit the machine. This may take several minutes.

MLJBase.fit!(mach)

# Construct the joint posterior:

predictor_joint = MLJBase.predict_joint(mach, iris[!, feature_columns])
typeof(predictor_joint)

# Draw a single sample from the joint posterior:

single_sample = rand(predictor_joint)

# For each row in the dataset, construct the marginal posterior predictive distribution

predictor_marginal = MLJBase.predict(mach, iris[!, feature_columns])

# predictor_marginal is a `UnivariateFiniteVector`

typeof(predictor_marginal)

# `predictor_marginal` has one element for each row in the data set

@show size(predictor_marginal); @show size(iris, 1);

# Use cross-validation to evaluate the model with respect to the Brier score:

evaluate!(mach, resampling=CV(; nfolds = 4, shuffle = true), measure=brier_score, operation=MLJBase.predict)

# Use cross-validation to evaluate the model with respect to accuracy:

evaluate!(mach, resampling=CV(; nfolds = 4, shuffle = true), measure=accuracy, operation=MLJBase.predict_mode)

# The cross-validated accuracy is greater than 90%, which is pretty good!


text:
mutable struct QuadLat{S, T, U} <: AbsLat{S}
  space::QuadSpace{S, T}
  pmat::U
  gram::T

  function QuadLat(K::S, G::T, P::U) where {S, T, U}
    space = QuadSpace(K, G)
    z = new{S, T, U}(space, P)
    return z
  end

  function QuadLat(K::S, G::T) where {S, T}
    n = nrows(G)
    M = pseudo_matrix(identity_matrix(K, n))
    return QuadLat(K, G, M)
  end
end

mutable struct HermLat{S, T, U, V, W} <: AbsLat{S}
  space::HermSpace{S, T, U, W}
  pmat::V
  gram::T

  function HermLat(E::S, G::U, P::V) where {S, U, V}
    @assert degree(E) == 2
    A = automorphisms(E)
    a = gen(E)
    if A[1](a) == a
      involution = A[2]
    else
      involution = A[1]
    end

    K = base_field(E)

    space = HermSpace(E, K, G, involution)

    z = new{S, typeof(K), U, V, typeof(involution)}(space, P)
    return z
  end

  function HermLat(E::S, G::U) where {S, U}
    n = nrows(G)
    M = pseudo_matrix(identity_matrix(E, n))
    return HermLat(E, G, M)
  end
end

function diagonal(L::AbsLat)
  D, _ = _gram_schmidt(gram_matrix_of_basis(L), involution(L))
  return diagonal(D)
end

pseudo_matrix(L::AbsLat) = L.pmat

coefficient_ideals(L::AbsLat) = coefficient_ideals(pseudo_matrix(L))

basis_matrix(L::AbsLat) = matrix(pseudo_matrix(L))

ambient_space(L::AbsLat) = L.space

rank(L::AbsLat) = nrows(L.pmat)

base_algebra(L::QuadLat) = L.space.K

base_algebra(L::HermLat) = L.space.E

involution(L::QuadLat) = identity

involution(L::HermLat) = involution(ambient_space(L))

function gram_matrix_of_basis(L::AbsLat)
  return gram_matrix(ambient_space(L), L.pmat.matrix)
end

function degree(L::AbsLat)
  return dim(ambient_space(L))
end

# Check if one really needs minimal
# Steinitz form is not pretty
function generators(L::AbsLat; minimal::Bool = true)
  St = _steinitz_form(pseudo_matrix(L), Val{false})
  d = nrows(St)
  n = degree(L)
  T = elem_type(base_algebra(L))
  v = Vector{T}[]
  for i in 1:(d - 1)
    @assert isprincipal(coefficient_ideals(St)[i])[1]
    push!(v, T[basis_matrix(L)[i, j] for j in 1:d])
  end

  I = numerator(coefficient_ideals(St)[d])
  den = denominator(coefficient_ideals(St)[d])
  if minimal && base_algebra(L) isa AnticNumberField
    b, a = isprincipal(I)
    if b
      push!(v, T[base_algebra(L)(a)//den * basis_matrix(L)[n, j] for j in 1:d])
    end
    return v
  end

  _assure_weakly_normal_presentation(I)
  push!(v, T[base_algebra(L)(I.gen_one)//den * basis_matrix(L)[n, j] for j in 1:d])
  push!(v, T[base_algebra(L)(I.gen_two)//den * basis_matrix(L)[n, j] for j in 1:d])

  return v
end

function discriminant(L::AbsLat)
  d = det(gram_matrix_of_basis(L))
  v = involution(L)
  C = coefficient_ideals(L)
  I = prod(J for J in C)
  return d * I * v(I)
end

function hasse_invariant(L::QuadLat, p)
  return _hasse_invariant(diagonal(L), p)
end

function hasse_invariant(L::HermLat, p)
  throw(error("The lattice must be quadratic"))
end

function witt_invariant(L::QuadLat, p::NfOrdIdl)
  h = hasse_invariant(L, p)
  F = gram_matrix_of_basis(L)
  dett = det(F)
  K = base_algebra(L)
  ncolsFmod8 = mod(ncols(F), 8)
  if ncolsFmod8 == 3 || ncolsFmod8 == 4
    c = -dett
  elseif ncolsFmod8 == 5 || ncolsFmod8 == 6
    c = K(-1)
  elseif ncolsFmod8 == 7 || ncolsFmod8 == 0
    c = dett
  else
    c = K(1)
  end
  return h * hilbert_symbol(K(-1), c, p)
end

function witt_invariant(L::QuadLat, p::InfPlc)
  if iscomplex(p)
    error("Dont' know what to do. Markus returns false")
  end

  h = hasse_invariant(L, p)
  F = gram_matrix_of_basis(L)
  dett = det(F)
  K = base_algebra(L)
  ncolsFmod8 = mod(ncols(F), 8)
  if ncolsFmod8 == 3 || ncolsFmod8 == 4
    c = -dett
  elseif ncolsFmod8 == 5 || ncolsFmod8 == 6
    c = K(-1)
  elseif ncolsFmod8 == 7 || ncolsFmod8 == 0
    c = dett
  else
    c = K(1)
  end
  @assert !iszero(c)
  if isnegative(c, p)
    return -h
  else
    return h
  end
end

function isrationally_equivalent(L::QuadLat, M::QuadLat, p::NfOrdIdl)
  GL = gram_matrix_of_basis(L)
  GM = gram_matrix_of_basis(M)
  if GL == GM
    return true
  end

  return rank(GL) == rank(GM) && islocal_square(det(GL) * det(GM), p) && hasse_invariant(L, p) == hasse_invariant(M, p)
end

function isequivalent(L::QuadLat, M::QuadLat, p::InfPlc)
  if rank(L) != rank(M)
    return false
  end

  if iscomplex(p)
    return true
  end

  DL = diagonal(L)
  DM = diagonal(M)
  return count(x -> isnegative(x, p), DL) == count(x -> isnegative(x, p), DM)
end


text:
using LinearAlgebra, Arpack
using Superfluids
using Superfluids: argand, cloud, operators, matrices

s = Superfluid{2}(3000, (x,y) -> (x^2+y^2)/2)

# Test values of discretisation and complex

fuzerate(::Type{Discretisation}) = [
    FDDiscretisation{2}(30,0.4),
    FDDiscretisation{2}(30,0.5),
    FDDiscretisation{2}(50,0.4),
    FDDiscretisation{2}(50,0.5)
]

fuzerate(::Type{Complex{Float64}}, n=7) =
    Complex.(randn(n),randn(n))

struct WaveFunction{N} end

function fuzerate(::Type{WaveFunction{2}}) = [
    (x,y) -> exp(-abs2(Complex(x,y))/2)),
    (x,y) -> x*exp(-abs2(Complex(x,y))/2)),
    (x,y) -> y*exp(-abs2(Complex(x,y))/2)),
    (x,y) -> (y+2im*x)*exp(-abs2(Complex(x,y))/2)),
    (x,y) -> (y+2im*x)*exp(-abs2(Complex(x,y))/2)),
    (x,y) -> (y+2im*x)*exp(-abs2(Complex(x,y))/2)),
    (x,y) -> Complex(x,y)*exp(-abs2(Complex(x,y))/2)),
    (x,y) -> Complex(x,-y)*exp(-abs2(Complex(x,y))/2))
]

# # Imprinted offset vortex test
# rv = Complex(1.5, 0.3)
# f1(z) = @. exp(-abs2(z)/18)
# f2(z) = @. f1(z)*(z-rv)/√(abs2(z-rv)+1)
#
# q1 = f2.(argand(d1))
# normalize!(q1)
# q2 = f2.(argand(d2))
# normalize!(q2)
#
# U1 = Superfluids.operators(s,d1,:U) |> only
# U2 = Superfluids.operators(s,d2,:U) |> only

"f(::Discretisation) = operator, g(z) is wave function"
function consistent_operator(f, g)
    zs = fuzerate(Complex{Float64})
    uzs = []
    for d = fuzerate(Discretisation)
        u = f(d)
        q = g.(argand(d))
        normalize!(q)
        uq = Superfluids.interpolate(d, u(q))
        push!(uzs, uq.(zs))
    end
    zs, uzs
end

function sho_error(d)
    s = Superfluid{2}(0, (x,y) -> (x^2+y^2)/2)
    z = argand(d)
    w = normalize(@. exp(-abs2(z)/2))
    q = relaxed_state(sho, d)
    norm(q - w*(w'*q))
end

TODO figure out why this breaks JuliaFormatter
"Expectation values of operators in SHO ground state"
function sho_energies(d)
    z = argand(d)
    w = normalize(@. exp(-abs2(z)/2))
    Δ, U, J = Superfluids.operators(d)
    Dict(
        :T=>dot(w,Δ(w)),
        :U=>dot(w,U(w)),
        :J=>dot(w,J(w))
    )
end

# check energies are (1/2, 0, 1/2) for SHO

# check derivatives of cloud(d)

# φ = ground_state(S,D)
#
# zs, us = consistent_operator(d->only(operators(s, d, :U)), f2)