text:
#=
Regex functions for Str strings

Copyright 2018 Gandalf Software, Inc., Scott P. Jones, and other contributors to the Julia language
Licensed under MIT License, see LICENSE.md
Based in part on julia/base/regex.jl
=#

using Base.PCRE

using Base: DEFAULT_COMPILER_OPTS, DEFAULT_MATCH_OPTS
import Base: Regex, match

const Regex_CSEs = Union{ASCIICSE,Latin_CSEs,Binary_CSEs}

# Default is to act as if validated UTF8
cvt_compile(::Type{<:CSE}, co) = co | PCRE.NO_UTF_CHECK | PCRE.UTF
cvt_match(::Type{<:CSE}, co)   = co | PCRE.NO_UTF_CHECK

cvt_compile(::Type{Regex_CSEs}, co) = co & ~PCRE.UTF
cvt_match(::Type{Regex_CSEs}, co)   = co & ~PCRE.NO_UTF_CHECK

cvt_compile(::Type{<:Str{C}}, co) where {C<:CSE} = cvt_compile(C, co)
cvt_match(::Type{<:Str{C}}, co) where {C<:CSE}   = cvt_match(C, co)

# String type is not validate
cvt_compile(::Type{String}, co) = (co & ~PCRE.NO_UTF_CHECK) | PCRE.UTF
cvt_match(::Type{String}, co)   = (co & ~PCRE.NO_UTF_CHECK)

Regex(pattern::MaybeSub{<:Str{C}}, co, mo) where {C<:Byte_CSEs} =
    Regex(String(pattern), cvt_compile(C, co), mo)
Regex(pattern::MaybeSub{<:Str{C}}) where {C<:Byte_CSEs} =
    Regex(String(pattern),
          cvt_compile(C, DEFAULT_COMPILER_OPTS), cvt_match(C, DEFAULT_MATCH_OPTS))

function _update_compiler_opts(flags)
    options = DEFAULT_COMPILER_OPTS
    for f in flags
        options |= f=='i' ? PCRE.CASELESS  :
                   f=='m' ? PCRE.MULTILINE :
                   f=='s' ? PCRE.DOTALL    :
                   f=='x' ? PCRE.EXTENDED  :
                   throw(ArgumentError("unknown regex flag: $f"))
    end
    options
end

Regex(pattern::MaybeSub{<:Str{C}}, flags::AbstractString) where {C<:Byte_CSEs} =
    Regex(pattern, cvt_compile(C, _update_compile_opts(flags)), cvt_match(C, DEFAULT_MATCH_OPTS))

struct RegexMatchStr{T<:AbstractString}
    match::T
    captures::Vector{Union{Nothing,T}}
    offset::Int
    offsets::Vector{Int}
    regex::Regex
end

function show(io::IO, m::RegexMatchStr)
    print(io, "RegexMatchStr(")
    show(io, m.match)
    idx_to_capture_name = PCRE.capture_names(m.regex.regex)
    if !is_empty(m.captures)
        print(io, ", ")
        for i = 1:length(m.captures)
            # If the capture group is named, show the name.
            # Otherwise show its index.
            print(io, get(idx_to_capture_name, i, i), "=")
            show(io, m.captures[i])
            i < length(m.captures) && print(io, ", ")
        end
    end
    print(io, ")")
end

# Capture group extraction
getindex(m::RegexMatchStr, idx::Integer) = m.captures[idx]
function getindex(m::RegexMatchStr, name::Symbol)
    idx = PCRE.substring_number_from_name(m.regex.regex, name)
    idx <= 0 && error("no capture group named $name found in regex")
    m[idx]
end

getindex(m::RegexMatchStr, name::AbstractString) = m[Symbol(name)]

"""
Check if the compile flags match the current search

If not, free up old compilation, and recompile
For better performance, the Regex object should hold spots for 5 compiled regexes,
for 8-bit, UTF-8, 16-bit, UTF-16, and 32-bit code units
"""
function check_compile(::Type{C}, regex::Regex) where {C<:Union{CSE,String}}
    re = regex.regex
    # ASCII is compatible with all (current) types, don't recompile
    C == ASCIICSE && re != C_NULL && return
    oldopt = regex.compile_options
    cvtcomp = cvt_compile(C, oldopt)
    if cvtcomp != oldopt
        regex.compile_options = cvtcomp
        re == C_NULL || (PCRE.free_re(re); regex.regex = re = C_NULL)
        regex.match_data == C_NULL ||
            (PCRE.free_match_data(regex.match_data); regex.match_data = C_NULL)
    end
    if re == C_NULL
        regex.regex = re = PCRE.compile(regex.pattern, cvtcomp)
        PCRE.jit_compile(re)
        regex.match_data = PCRE.create_match_data(re)
        regex.ovec = PCRE.get_ovec(regex.match_data)
    end
    nothing
end

function _match(::Type{C}, re, str, idx, add_opts) where {C<:CSE}
    check_compile(C, re)
    PCRE.exec(re.regex, str, idx-1, cvt_match(C, re.match_options | add_opts), re.match_data) ||
        return nothing
    ovec = re.ovec
    n = div(length(ovec),2) - 1
    mat = SubString(str, ovec[1]+1, prevind(str, ovec[2]+1))
    cap = Union{Nothing,SubString{typeof(str)}}[ovec[2i+1] == PCRE.UNSET ?
                                               nothing :
                                               SubString(str, ovec[2i+1]+1,
                                                         prevind(str, ovec[2i+2]+1)) for i=1:n]
    off = Int[ ovec[2i+1]+1 for i=1:n ]
    RegexMatchStr(mat, cap, Int(ovec[1]+1), off, re)
end

match(re::Regex, str::MaybeSub{<:Str{C}}, idx::Integer,
      add_opts::UInt32=UInt32(0)) where {C<:CSE} =
    error("Regex not supported yet on strings with codeunit == UInt16 or UInt32")
match(re::Regex, str::MaybeSub{<:Str{C}}, idx::Integer,
      add_opts::UInt32=UInt32(0)) where {C<:Regex_CSEs} =
    _match(C, re, str, Int(idx), add_opts)
match(re::Regex, str::MaybeSub{<:Str{UTF8CSE}}, idx::Integer, add_opts::UInt32=UInt32(0)) =
    _match(UTF8CSE, re, str, Int(idx), add_opts)
match(re::Regex, str::MaybeSub{<:Str{_LatinCSE}}, idx::Integer, add_opts::UInt32=UInt32(0)) =
    _match(LatinCSE, re, Str{LatinCSE}(str), Int(idx), add_opts)

function _find(::Type{C}, re, str) where {C}
    isempty(str) && return _not_found
    check_compile(C, re)
    (PCRE.exec(re.regex, str, 0, cvt_match(C, re.match_options), re.match_data)
     ? ((Int(re.ovec[1]) + 1) : prevind(str, Int(re.ovec[2]) + 1)) : _not_found)
end

function _find(::Type{C}, re, str, idx) where {C}
    if idx > ncodeunits(str)+1
        @boundscheck boundserr(str, idx)
        return _not_found
    end
    check_compile(C, re)
    (PCRE.exec(re.regex, str, idx-1, cvt_match(C, re.match_options), re.match_data)
     ? ((Int(re.ovec[1]) + 1) : prevind(str, Int(re.ovec[2]) + 1)) : _not_found)
end

find(::Type{Fwd}, re::Regex, str::MaybeSub{<:AbstractString}, idx::Integer) =
    throw(ArgumentError("regex search is only available for the String or Str types with " *
                        "codeunit == UInt8, or substrings of those types, " *
                        "use UTF8Str to convert"))

find(::Type{Fwd}, re::Regex, str::MaybeSub{<:Str{C}}, idx::Integer) where {C<:Regex_CSEs} =
    _find(C, re, str, idx)
find(::Type{Fwd}, re::Regex, str::MaybeSub{<:Str{UTF8CSE}}, idx::Integer) =
    _find(UTF8CSE, re, str, idx)
find(::Type{Fwd}, re::Regex, str::MaybeSub{<:Str{_LatinCSE}}, idx::Integer) =
    _find(LatinCSE, re, str, idx)
find(::Type{Fwd}, re::Regex, str::MaybeSub{String}, idx::Integer) =
    _find(UTF8CSE, re, str, idx)

find(::Type{First}, re::Regex, str::MaybeSub{<:AbstractString}) = find(Fwd, re, str, 1)
find(::Type{First}, re::Regex, str::MaybeSub{<:Str{C}}) where {C<:Regex_CSEs} = _find(C, re, str)
find(::Type{First}, re::Regex, str::MaybeSub{<:Str{_LatinCSE}}) = _find(LatinCSE, re, str)
find(::Type{First}, re::Regex, str::MaybeSub{String}) = _find(UTF8CSE, re, str)


text:
function BinomialTest()
    return BinomialTest(())
end

function binomial_test(obj::BinomialTest, arg0::jint, arg1::jint, arg2::jdouble, arg3::AlternativeHypothesis)
    return jcall(obj, "binomialTest", jdouble, (jint, jint, jdouble, AlternativeHypothesis), arg0, arg1, arg2, arg3)
end

function binomial_test(obj::BinomialTest, arg0::jint, arg1::jint, arg2::jdouble, arg3::AlternativeHypothesis, arg4::jdouble)
    return jcall(obj, "binomialTest", jboolean, (jint, jint, jdouble, AlternativeHypothesis, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

ext: jl


text:
ds = Table("wikipedia")
intervals = [Interval("a", "b")]
granularity = SimpleGranularity("hour")
dimspec = DefaultDS("dimension")
tpnmetric = Dimension()
query = Contains("a")

@testset "Timeseries" begin
    @test_throws TypeError(:typeassert, "", Filter, 1) Timeseries(
        ds, intervals, granularity, filter=1)
    @test_throws TypeError(:typeassert, "", Vector{<:Aggregator}, 1) Timeseries(
        ds, intervals, granularity, aggregations=1)
    @test_throws TypeError(:typeassert, "", Vector{<:PostAggregator}, 1) Timeseries(
        ds, intervals, granularity, postAggregations=1)
    @test_throws TypeError(:typeassert, "", Vector{<:VirtualColumn}, 1) Timeseries(
        ds, intervals, granularity, virtualColumns=1)
    @test_throws TypeError(:typeassert, "", Bool, 1) Timeseries(
        ds, intervals, granularity, descending=1)
    @test_throws ErrorException("limit must be a non-negative integer") Timeseries(
        ds, intervals, granularity, limit=-1)
    @test_throws TypeError(:typeassert, "", Dict, 1) Timeseries(
        ds, intervals, granularity, context=1)
end

@testset "Dimension" begin
    @test_throws ErrorException("Invalid ordering") Dimension(ordering="a")
    @test_throws TypeError(:typeassert, "", String, 1) Dimension(previousStop=1)
end

@testset "TopN" begin
    @test_throws TypeError(:typeassert, "", Filter, 1) TopN(
        ds, intervals, granularity, dimspec, 1, tpnmetric, filter=1)
    @test_throws TypeError(:typeassert, "", Vector{<:Aggregator}, 1) TopN(
        ds, intervals, granularity, dimspec, 1, tpnmetric, aggregations=1)
    @test_throws TypeError(:typeassert, "", Vector{<:PostAggregator}, 1) TopN(
        ds, intervals, granularity, dimspec, 1, tpnmetric, postAggregations=1)
    @test_throws TypeError(:typeassert, "", Vector{<:VirtualColumn}, 1) TopN(
        ds, intervals, granularity, dimspec, 1, tpnmetric, virtualColumns=1)
    @test_throws TypeError(:typeassert, "", Dict, 1) TopN(
        ds, intervals, granularity, dimspec, 1, tpnmetric, context=1)
end

@testset "OrderByColumn" begin
    @test_throws ErrorException("Invalid direction") OrderByColumn("dimension", "a")
    @test_throws ErrorException("Invalid dimensionOrder") OrderByColumn("dimension", "ascending", dimensionOrder="a")
end

@testset "DefaultLS" begin
    @test_throws ErrorException("limit must be a non-negative integer") DefaultLS(limit=-1)
    @test_throws ErrorException("offset must be a non-negative integer") DefaultLS(offset=-1)
end

@testset "GroupBy" begin
    @test_throws TypeError(:typeassert, "", Filter, 1) GroupBy(
        ds, [dimspec], intervals, granularity, filter=1)
    @test_throws TypeError(:typeassert, "", Vector{<:Aggregator}, 1) GroupBy(
        ds, [dimspec], intervals, granularity, aggregations=1)
    @test_throws TypeError(:typeassert, "", Vector{<:PostAggregator}, 1) GroupBy(
        ds, [dimspec], intervals, granularity, postAggregations=1)
    @test_throws TypeError(:typeassert, "", Vector{<:VirtualColumn}, 1) GroupBy(
        ds, [dimspec], intervals, granularity, virtualColumns=1)
    @test_throws TypeError(:typeassert, "", LimitSpec, 1) GroupBy(
        ds, [dimspec], intervals, granularity, limitSpec=1)
    @test_throws TypeError(:typeassert, "", Dict, 1) GroupBy(
        ds, [dimspec], intervals, granularity, context=1)
    @test_throws TypeError(:typeassert, "", Union{HavingSpec, Filter}, 1) GroupBy(
        ds, [dimspec], intervals, granularity, having=1)
    @test_throws TypeError(:typeassert, "", Vector{Vector{String}}, 1) GroupBy(
        ds, [dimspec], intervals, granularity, subtotalsSpec=1)
end

@testset "Scan" begin
    @test_throws TypeError(:typeassert, "", Filter, 1) Scan(
        ds, intervals, filter=1)
    @test_throws TypeError(:typeassert, "", Vector{String}, 1) Scan(
        ds, intervals, columns=1)
    @test_throws ErrorException("Invalid order") Scan(
        ds, intervals, order="a")
    @test_throws TypeError(:typeassert, "", Vector{<:VirtualColumn}, 1) Scan(
        ds, intervals, virtualColumns=1)
    @test_throws ErrorException("limit must be a non-negative integer") Scan(
        ds, intervals, limit=-1)
    @test_throws ErrorException("offset must be a non-negative integer") Scan(
        ds, intervals, offset=-1)
    @test_throws ErrorException("Invalid resultFormat") Scan(
        ds, intervals, resultFormat=1)
    @test_throws ErrorException("batchSize must be a non-negative integer") Scan(
        ds, intervals, batchSize=-1)
    @test_throws TypeError(:typeassert, "", Dict, 1) Scan(
        ds, intervals, context=1)
    @test_throws TypeError(:typeassert, "", Bool, 1) Scan(
        ds, intervals, legacy=1)
end

@testset "Search" begin
    @test_throws TypeError(:typeassert, "", Granularity, 1) Search(
        ds, intervals, query, granularity=1)
    @test_throws TypeError(:typeassert, "", Filter, 1) Search(
        ds, intervals, query, filter=1)
    @test_throws TypeError(:typeassert, "", Vector{<:VirtualColumn}, 1) Search(
        ds, intervals, query, virtualColumns=1)
    @test_throws ErrorException("Invalid sort value") Search(
        ds, intervals, query, sort="a")
    @test_throws ErrorException("limit must be a non-negative integer") Search(
        ds, intervals, query, limit=-1)
    @test_throws TypeError(:typeassert, "", Dict, 1) Search(
        ds, intervals, query, context=1)
end

@testset "TimeBoundary" begin
    @test_throws ErrorException("Invalid bound") TimeBoundary(ds, bound="a")
    @test_throws TypeError(:typeassert, "", Filter, 1) TimeBoundary(ds, filter=1)
    @test_throws TypeError(:typeassert, "", Vector{<:VirtualColumn}, 1) TimeBoundary(ds, virtualColumns=1)
    @test_throws TypeError(:typeassert, "", Dict, 1) TimeBoundary(ds, context=1)
end

@testset "SegmentMetadata" begin
    @test_throws TypeError(:typeassert, "", Vector{<:Interval}, 1) SegmentMetadata(
        ds, intervals=1)
    @test_throws TypeError(:typeassert, "", Union{String, Vector{String}}, 1) SegmentMetadata(
        ds, toInclude=1)
    @test_throws TypeError(:typeassert, "", Bool, 1) SegmentMetadata(
        ds, merge=1)
    @test_throws ErrorException("Invalid analysisTypes") SegmentMetadata(
        ds, analysisTypes=["a"])
    @test_throws TypeError(:typeassert, "", Bool, 1) SegmentMetadata(
        ds, lenientAggregatorMerge=1)
    @test_throws TypeError(:typeassert, "", Vector{<:VirtualColumn}, 1) SegmentMetadata(
        ds, virtualColumns=1)
    @test_throws TypeError(:typeassert, "", Dict, 1) SegmentMetadata(
        ds, context=1)
end

@testset "DatasourceMetadata" begin
    @test_throws TypeError(:typeassert, "", Dict, 1) DatasourceMetadata(ds, context=1)
end

@testset "Parameter" begin
    @test_throws ErrorException("Invalid 'type'") Parameter("a", 1)
end

@testset "SQL" begin
    @test_throws TypeError(:typeassert, "", Vector{Parameter}, 1) Sql("query", parameters=1)
    @test_throws ErrorException("Invalid resultFormat") Sql("query", resultFormat="a")
    @test_throws TypeError(:typeassert, "", Bool, 1) Sql("query", header=1)
    @test_throws TypeError(:typeassert, "", Dict, 1) Sql("query", context=1)
end


text:
# ---
# title: 846. Hand of Straights
# id: problem846
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Ordered Map
# link: <https://leetcode.com/problems/hand-of-straights/description/>
# hidden: true
# ---
#
# Alice has a `hand` of cards, given as an array of integers.
#
# Now she wants to rearrange the cards into groups so that each group is size
# `W`, and consists of `W` consecutive cards.
#
# Return `true` if and only if she can.
#
# **Note:** This question is the same as 1296:
# <https://leetcode.comhttps://leetcode.com/problems/divide-array-in-sets-of-k-consecutive-numbers/>
#
#
#
# **Example 1:**
#
#
#
#     Input: hand = [1,2,3,6,2,3,4,7,8], W = 3
#     Output: true
#     Explanation: Alice's hand can be rearranged as [1,2,3],[2,3,4],[6,7,8]
#
#
# **Example 2:**
#
#
#
#     Input: hand = [1,2,3,4,5], W = 4
#     Output: false
#     Explanation: Alice's hand can't be rearranged into groups of 4.
#
#
#
#
#
# **Constraints:**
#
#   * `1 <= hand.length <= 10000`
#   * `0 <= hand[i] <= 10^9`
#   * `1 <= W <= hand.length`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module test_julia_wasm

using Test

s = """
0000000: 0061 736d                                 ; WASM_BINARY_MAGIC
0000004: 0100 0000                                 ; WASM_BINARY_VERSION
; section "Type" (1)
0000008: 01                                        ; section code
0000009: 00                                        ; section size (guess)
000000a: 01                                        ; num types
; type 0
000000b: 60                                        ; func
000000c: 02                                        ; num params
000000d: 7e                                        ; i64
000000e: 7e                                        ; i64
000000f: 01                                        ; num results
0000010: 7e                                        ; i64
0000009: 07                                        ; FIXUP section size
; section "Function" (3)
0000011: 03                                        ; section code
0000012: 00                                        ; section size (guess)
0000013: 01                                        ; num functions
0000014: 00                                        ; function 0 signature index
0000012: 02                                        ; FIXUP section size
; section "Code" (10)
0000015: 0a                                        ; section code
0000016: 00                                        ; section size (guess)
0000017: 01                                        ; num functions
; function body 0
0000018: 00                                        ; func body size (guess)
0000019: 00                                        ; local decl count
000001a: 20                                        ; get_local
000001b: 00                                        ; local index
000001c: 20                                        ; get_local
000001d: 01                                        ; local index
000001e: 7c                                        ; i64.add
000001f: 0f                                        ; return
0000020: 0b                                        ; end
0000018: 08                                        ; FIXUP func body size
0000016: 0a                                        ; FIXUP section size
"""

function wasm_bytes(s)
    bytes = UInt8[]
    pat = r"[0-9a-f]{2}"
    for line in split(s, '\n')
        isempty(line) && continue
        line[1] == ';' && continue
        hex = line[9:51]
        bs = parse.(UInt8, base=16, [m.match for m in eachmatch(pat, hex)])
        push!(bytes, bs...)
    end
    bytes
end

@test wasm_bytes(s)[1:4] == UInt8[0x00, 0x61, 0x73, 0x6d]

end # module test_julia_wasm

using Jive
@If VERSION >= v"1.2-DEV" module test_julia_wasm_isjsvm

using Test

@test !Sys.isjsvm()

end # @If VERSION >= v"1.2-DEV" module test_julia_wasm_isjsvm


text:
export
    C1, C2, C3,
    C2Shared, C3Shared,
    AbstractBasisType,
    BasisType,
    twoPhotons,
    retrieveSolution,
    alpharange,
    calculate_basis,
    FullCorrelation_parallized,
    energy,
    flab,
    calculate_triple_products,
    calculate_triple_products_fast,
    FullC3,
    c3_slice,
    integrate_c3_shell,
    renormalize_correlation,
    symmetrize_correlation,
    add_Gaussian_filter,
    postprocess_correlations,
    c2_difference

const C1        = Vector{Float64}
const C1Shared  = SharedArray{Float64, 1}
const C2        = Array{Float64, 3}
const C2Shared  = SharedArray{Float64, 3}
const C3        = Array{Float64, 5}
const C3Shared  = SharedArray{Float64, 5}

abstract type AbstractBasisType end

"""Datatype for precalculated three photon correlation basis function"""
type BasisType <: AbstractBasisType
    wignerlist::Array{Float64,1}
    indices::Array{Int64,2}
    PAcombos::Array{Int64,2}
    B::Array{Float32,2}
    h_P::Array{Float32,2}
    basislen::Int64
    N::Int64
    L::Int64
    LMAX::Int64
    lrange::StepRange
    ctr::Dict
    rtc::Dict
    K::Int64
    lambda::Float64
    dq::Float64

    function BasisType( wignerlist::Array{Float64,1}, indices::Array{Int64,2}, PAcombos::Array{Int64,2}, B::Array{Float32,2}, h_P::Array{Float32,2}, basislen::Int64, N::Int64, L::Int64, LMAX::Int64, lrange::StepRange, ctr::Dict, rtc::Dict, K::Int64, lambda::Float64, dq::Float64)
        new(wignerlist, indices, PAcombos, B, h_P, basislen, N, L, LMAX, lrange, ctr, rtc, K, lambda, dq)
    end

    function BasisType(N::Int64,L::Int64,LMAX::Int64,K::Int64,lambda::Float64,dq::Float64)
        lrange = (0:2:L)
        ctr = Dict(l => Umat(l) for l=lrange)
        rtc = Dict(l => ctranspose(ctr[l]) for l=lrange)

        new(zeros(Float64,1), zeros(Int64,2,2), zeros(Int64,2,2), zeros(Float64,2,2), Array{Float32}(2,2), 0, N, L, LMAX, lrange, ctr, rtc, K, lambda, dq)
    end
end

"""Helper function to calculate basis"""
function calculate_combolist(K::Int64, mcombolist)
    PAcombos = Int64[]
    ki::Int64 = 1
    num::Int64 = 0
    for k1 = 1:K
        for k2 = 1:k1
            for k3 = 1:k2
                for (l1,l2,l3,mcombos, jstart) in mcombolist
                    num+=1
                    append!(PAcombos, [k1,k2,k3,ki,l1,l2,l3,jstart,mcombos])
                end
                ki+=1
            end
        end
    end
    return reshape(PAcombos, 9, num)
end

"""Precalculated values for calculating the three-photon and two-photon correlation"""
function calculate_basis(L::Int64, LMAX::Int64, N::Int64, K::Int64, lambda::Float64=0.0, dq::Float64=0.0, forIntensity=true)

    println("Calculating complex basis with N=$N L=$L K=$K (LMAX=$LMAX, lambda=$lambda, dq=$dq).")

    lrange = forIntensity ? (0:2:L) : (0:L)
    basislen = sum(abs(wigner(l1,m1, l2,m2, l3,-m3)) > 100*eps() ? 1 : 0 for l1=lrange for l2=lrange for l3=lrange for m1=-l1:l1 for m2=-l2:l2 for m3=-l3:l3)
    klength = Integer(K*(K+1)*(K+2)/6)
    qlist = Float64[acos(k*lambda*dq/(4*pi)) for k=1:K]

    wignerlist = Array{Float64}(basislen)
    indiceslist = Array{Int64}(9, basislen)
    # B = SharedArray(Float32, 2*N^2, basislen)
    # P = SharedArray(Float32, klength, basislen)
    B = zeros(Float32, 2*N^2, basislen)
    P = zeros(Float32, klength, basislen)
    mcombolist = Vector{Int64}[]

    ctr = Dict(l => Umat(l) for l=lrange)
    rtc = Dict(l => ctranspose(ctr[l]) for l=lrange)

    i = 1
    for l1 = lrange
        for l2 = lrange
            for l3 = lrange
                mcombos = 0
                mcombo_istart = i
                for m1 = -l1:l1
                    for m2 = -l2:l2
                        for m3 = -l3:l3
                            w = wigner(l1,m1, l2,m2, l3,-m3)
                            if abs(w) > 100*eps()
                                wignerlist[i] = w
                                indiceslist[:,i] = Int64[seanindex(l1,m1,LMAX), seanindex(l2,m2,LMAX),seanindex(l3,-m3,LMAX), l1,m1, l2,m2, l3,-m3]
                                i+=1
                                mcombos += 1
                            end
                        end
                    end
                end
                if mcombos>0 push!(mcombolist, Int64[l1,l2,l3,mcombos,mcombo_istart]) end
            end
        end
    end

    PAcombos = calculate_combolist(K,mcombolist)

    #Commented out parallelism below because of BUS error on cluster, has to to with SharedArrays, may go away with Julia 0.6 upgrade
    #Maybe try @everywhere gc() every 100 iterations later
    # @time begin
        # @sync @parallel for i=1:basislen
        for i=1:basislen
            l1,m1,l2,m2,l3,m3 = indiceslist[4:9,i]
            w = wignerlist[i]
            #Here, only the real part of the complex exponential plays a role
            B[:,i] = reshape(Float64[cos(m2*a + m3*b) for a in alpharange(N), b in alpharange_2pi(2*N)], 2*N^2)
            P[:,i] = reshape(Float64[w*sphPlm(l1,m1,qlist[k1]) * sphPlm(l2,m2,qlist[k2]) * sphPlm(l3,m3,qlist[k3]) for k1=1:K for k2=1:k1 for k3=1:k2], klength)
            # if mod(i,1000) == 0
            #     @everywhere gc()
            # end
        end
    # end
    h_P = transpose(sdata(P))

    println("Calculation complete ($basislen basislen).")

    return BasisType(wignerlist, indiceslist, PAcombos, sdata(B), h_P, basislen, N, L, LMAX, lrange, ctr, rtc, K, lambda, dq)
end

complexBasis_choice = calculate_basis

function calculate_triple_products_fast(intensity::SphericalHarmonicsVolume, basis::BasisType)
    coeff = intensity.coeff
    PAcombos = basis.PAcombos
    indices = basis.indices
    wignerlist = basis.wignerlist
    P = basis.h_P
    PA = SharedArray{Float64}(Base.size(P))

    @sync @parallel for i = 1:Base.size(PAcombos)[2]
        k1,k2,k3,ki,l1,l2,l3,jstart,mcombos = PAcombos[:,i]
        ck1,ck2,ck3 = coeff[k1],coeff[k2],coeff[k3]
        As::Float64=0.0
        for n = 1:mcombos
            j = jstart + n - 1
            @inbounds @fastmath As += wignerlist[j]*real(ck1[indices[1, j]]*ck2[indices[2, j]]*ck3[indices[3, j]])
        end
        @inbounds @fastmath PA[jstart:jstart+mcombos-1,ki] = As * P[jstart:jstart+mcombos-1,ki]
    end
    return sdata(PA)
end

"Calculates the three photon correlation from spherical harmonics coefficients in a paralllized way."
function FullCorrelation_parallized(intensity::SphericalHarmonicsVolume, basis::BasisType, minimal::Bool=true, normalize::Bool=false, return_raw::Bool=false)

    PA = calculate_triple_products_fast(intensity, basis)
    res = basis.B*PA

    if return_raw return res end

    t = zeros(Float64, basis.N, 2*basis.N, basis.K, basis.K, basis.K)
    i = 1
    for k1=1:basis.K for k2=1:k1 for k3=1:k2
        t[:,:, k3, k2, k1] = reshape(res[:, i], basis.N, 2*basis.N)
        i += 1
    end end end
    return normalize ? t / sum(abs,t) : t
end

"""Central energy calculation function"""
function energy(intensity::SphericalHarmonicsVolume, basis::AbstractBasisType, c3ref::C3, measure::String="Bayes", negativity_factor::Float64=0.0)
    return energy(intensity, basis, c3ref, 1:basis.K, measure, negativity_factor)
end

"""Central energy calculation function"""
function energy(intensity::SphericalHarmonicsVolume, basis::AbstractBasisType, c3ref::C3, K3_range::UnitRange{Int64}, measure::String="Bayes", negativity_factor::Float64=0.0)
    c3 = FullCorrelation_parallized(intensity, basis, true, true, true)
    c3 = max.(c3, 1e-30) #Filter out results where a negative c3 value is expected. This happens in particular when starting structures are derived from sparse (histogrammed) photon correlations.

    # average_negativity = negativityCheck(deleteTerms(intensity,basis.K, basis.L))/basis.K
    # neg = exp(negativity_factor * average_negativity)
    neg = 1.0

    res = 0.0
    if measure == "Bayes"
        i = 1
        for k1 in K3_range
            for k2 = minimum(K3_range):k1
                for k3 = minimum(K3_range):k2
                    p = reshape(c3[:, i], basis.N, 2*basis.N)
                    i += 1
                    q = c3ref[:,:,k3,k2,k1]
                    p = p / sum(abs,p)
                    q = q / sum(abs,q)
                    @fastmath res += -1.0*sum(q.*log.(p))*tripletFactor(k1,k2,k3)
                end
            end
        end
    end
    return neg*res
end

"""Calculates the two photon correlation from spherical harmonics coefficients"""
function twoPhotons(volume::SphericalHarmonicsVolume, basis::BasisType, K2::Int64, minimal::Bool=false, normalize::Bool=false)
    c = zeros(Float64,basis.N,K2,K2)
    mdq = dq(volume)
    for k1=1:K2
        for k2=1:(minimal ? k1 : K2)
            slice = zeros(Float64, basis.N)
            for l = basis.lrange
                #NOTE: No need for conj for complex dot products (it's included in the definition)
                fac = Base.dot(cvec_get(volume,k1,l), cvec_get(volume,k2,l))
                slice += fac * Float64[ Plm(l,0,alpha_star(alpha, k1, k2, mdq, basis.lambda)) for alpha = alpharange(basis.N)]
            end
            c[:,k2,k1] = real(slice)
        end
    end
    return normalize ? c/sum(abs,c) : 1/(4*pi)*c
end

"""Alpharange of theory"""
function alpharange(N)
    da = pi/N
    linspace(da/2,pi-da/2,N)
end

"""Alpharange of histograms"""
function alpharange_2pi(N)
    da = 2*pi/N
    linspace(da/2,2*pi-da/2,N)
end

"""Retrieves a set of spherical harmonics coefficeints
Modified version to allow for ranges
"""
function retrieveSolution(c2::C2, L::Int64, LMAX::Int64, K2_range::UnitRange{Int64}, qmax::Float64, lambda::Float64)
    N,_,_ = size(c2)
    K2 = length(K2_range)
    K2_high = maximum(K2_range)
    K2_low = minimum(K2_range)
    @assert L < K2/2 #2*L+1 = K2

    #Create empty Spherical Harmonics volume
    intensity = SphericalHarmonicsVolume(LMAX, K2_high, qmax)
    println("Extracting solution over K2_range=$(K2_range) with K2=$(K2) and L=$L and qmax=$qmax.")
    mdq = dq(intensity)
    eigenvecs = Dict()
    eigenvals = Dict()
    Gmatrices = Dict()

    for l = 0:2:L
        G = zeros(Float64, K2, K2)
        for k1 in K2_range
            for k2 = K2_low:k1
                slice = c2[:,k2,k1]

                #symmetrize 2 photon correlation if lambda = 0.0
                if lambda == 0.0
                    slice = 0.5*(slice + reverse(slice))
                end
                #TODO: This might change in the future, with other scaling depending on the structure
                #Only takes low l orders into account for lower shells, otherwise fitting is unstable
                upperL = min(L, min(k1,k2)*2-2)
                A = Float64[ (1/(4*pi))*Plm(l,0,alpha_star(alpha, k1, k2, mdq, lambda)) for alpha = alpharange(N), l = 0:2:upperL]
                #This inversion is stable
                fac = A \ slice

                if l <= upperL
                    val = fac[round(Int64,l/2)+1]
                    G[k2-K2_low+1,k1-K2_low+1] = val
                    G[k1-K2_low+1,k2-K2_low+1] = val
                end
            end
        end
        #Diagonalize the matrix
        F = eigfact(Symmetric(G),K2-(2*l+1)+1:K2)
        # F = eigfact(G, permute=true, scale=false)

        eigenval, eigenvectors = F[:values], F[:vectors]

        #Calculate the vectors
        eigenvalmatrix = diagm(sqrt.(max.(0.0, eigenval)))
        eigenvecs[l] = eigenvectors'
        eigenvals[l] = eigenvalmatrix
        Gmatrices[l] = G
        m = eigenvalmatrix*eigenvectors'

        for k in K2_range
            cvec_set(intensity,k,l,(m[:,k-K2_low+1]))
            # vec = m[K2-(2*l):K2,k]
            # cvec_set(intensity,k,l, vec)
        end
    end

    intensity = real_to_comp(intensity)

    if negativityCheck(intensity, K2_range) > 0.33
        for k = 1:intensity.KMAX intensity.coeff[k] *= -1.0 end
    end

    return intensity#,eigenvecs,eigenvals,Gmatrices
end

"""Retrieves a set of spherical harmonics coefficeints"""
function retrieveSolution(c2::C2, L::Int64, LMAX::Int64, K2::Int64, qmax::Float64, lambda::Float64)
    return retrieveSolution(c2, L, LMAX, 1:K2, qmax, lambda)
end

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

"Part of the three photon correlation basis functions"
function flab(k1::Int64, k2::Int64, k3::Int64, dq::Float64, lambda::Float64, l1::Int64, l2::Int64, l3::Int64, a::Float64, b::Float64)
    fl = zero(Complex{Float64})
    theta1 = acos(k1*lambda*dq/(4*pi))
    theta2 = acos(k2*lambda*dq/(4*pi))
    theta3 = acos(k3*lambda*dq/(4*pi))
    if !(abs(l1-l2) <= l3 <= (l1+l2)) return 0.0 end

    for m1p=-l1:l1
        for m2p=-l2:l2
            for m3p=-l3:l3

                fac = wigner(l1,m1p, l2,m2p, l3,-m3p)

                if abs(fac) > 100*eps()
                    cont = sphPlm(l1,m1p,theta1) * sphPlm(l2,m2p,theta2) * sphPlm(l3,-m3p,theta3) * exp(1im*(m2p*a - m3p*b))

                    #The following is wrong, the curvature is already in the thetas, no need to adjust the alpha/betas (phew)
                    # cont = sphPlm(l1,m1p,theta1) * sphPlm(l2,m2p,theta2) * sphPlm(l3,-m3p,theta3) * exp(1im*(m2p*alpha_star(a, k1, k2, dq, lambda) - m3p*alpha_star(b, k1, k3, dq, lambda)))

                    fl += fac* cont
                end
            end
        end
    end
    return Float64(real(fl))
end

"Part of the three photon correlation basis functions"
function flab(l1::Int64, l2::Int64, l3::Int64, a::Float64, b::Float64)
    return flab(0, 0, 0, 0.0, 0.0, l1, l2, l3, a, b)
end

function c3_slice(ck1::Vector{Complex{Float64}}, ck2::Vector{Complex{Float64}}, ck3::Vector{Complex{Float64}}, N::Int64, L::Int64, LMAX::Int64, k1::Int64, k2::Int64, k3::Int64, dq::Float64, lambda::Float64)
    slice = zeros(Complex{Float64}, N, 2*N)
    for l1 = 0:2:L
        for l2 = 0:2:L
            for l3 = 0:2:L
                m = Float64[ flab(k1,k2,k3,dq,lambda, l1,l2,l3, a,b) for a=alpharange(N), b=alpharange_2pi(2*N)]
                if norm(m) > 3*eps()
                    for m1 = -l1:l1
                        for m2 = -l2:l2
                            for m3 = -l3:l3
                                w = wigner(l1,m1, l2,m2, l3,-m3)
                                if abs(w) > 3*eps()
                                    slice += m*w * ck1[seanindex(l1,m1,LMAX)] * ck2[seanindex(l2,m2,LMAX)] * ck3[seanindex(l3,-m3,LMAX)]
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    return real(slice)
end

"""Calculate complete three-photon correlation
This is just for testing because it's a slow calculation."""
function FullC3(volume::SphericalHarmonicsVolume, L::Int64, K::Int64, N::Int64, LMAX::Int64)

    coeff = volume.coeff
    c3 = zeros(Float64, N, 2*N, K, K, K)

    for k1 = 1:K
        ck1 = coeff[k1]
        for k2 = 1:k1
            ck2 = coeff[k2]
            for k3 = 1:k2
                ck3 = coeff[k3]
                c3[:,:, k3,k2,k1] = c3_slice(ck1, ck2, ck3, N, L)
            end
        end
    end
    return c3 #/ sum(abs, c3)
end

# """From a two-photon correlation with k1>k2 restriction, calculates the full version."""
# function complete_two_photon_correlation(c2::C2)
#     N,K,_= Base.size(c2)
#     Float64[k1 >= k2 ? c2[a,k2,k1] : c2[a,k1,k2] for a = 1:N, k1=1:K, k2=1:K]
# end
#
# """From a three-photon correlation with k1>k2>k3 restriction, calculates the full version.
# This is used to complete the histogrammed three-photon correlation before integration."""
# function complete_three_photon_correlation(c3::C3)
#     N,_,K,_,_ = Base.size(c3)
#     c3new = deepcopy(c3)
#     for k1 = 1:K
#         for k2 = 1:K
#             for k3 = 1:K
#                 if     k1 >= k3 && k3 >= k2 c3new[:,:,k3, k2, k1] = c3[:,:,k2,k3,k1]'
#
#                 elseif k3 >= k1 && k1 >= k2
#                     b = c3[:,:,k2,k1,k3]
#                     c3new[:,:,k3, k2, k1] = [b[j,mod(j-i,N)+1] for i = 1:N, j=1:N]
#                 elseif k3 >= k2 && k2 >= k1
#                     b = c3[:,:,k1,k2,k3]
#                     c3new[:,:,k3, k2, k1] = [b[mod(j-i,N)+1,j] for i = 1:N, j=1:N]
#
#                 elseif k2 >= k3 && k3 >= k1
#                     b = c3[:,:,k1,k3,k2]
#                     c3new[:,:,k3, k2, k1] = [b[mod(i-j,N)+1,i] for i = 1:N, j=1:N]
#                 elseif k2 >= k1 && k1 >= k3
#                     b = c3[:,:,k3,k1,k2]
#                     c3new[:,:,k3, k2, k1] = [b[i,mod(i-j,N)+1] for i = 1:N, j=1:N]
#                 end
#             end
#         end
#     end
#     return c3new
# end

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

function integrateShell_2pc_alt(intensity::SphericalHarmonicsVolume, N::Int64, K::Int64, iterations::Int64=Integer(1e7))
    da = pi/N
    c2 = zeros(N, K, K)
    c2counts = copy(c2)

    rotations = Matrix{Float64}[random_rotation(3) for k = 1: 10000]
    dq = dr(intensity)
    surf = getSurfaceVolume(intensity)

    for i = 1:iterations
        k1 = Int64(rand(1:K))
        k2 = Int64(rand(1:K))

        p1 = k1*dq*rotations[rand(1:length(rotations))]*[0,1,0]
        p2 = k2*dq*rotations[rand(1:length(rotations))]*[0,1,0]

        a = angle_between(p1,p2)
        ai = Int64(mod(floor(Int64, a/da),N)+1)

        @inbounds c2[ai,k2,k1] += real(getVolumeInterpolated(surf, p1))* real(getVolumeInterpolated(surf, p2))
        @inbounds c2counts[ai,k2,k1] += 1.0

    end
    return c2./c2counts
end

function integrateShell_3pc(intensity::SphericalHarmonicsVolume, N::Int64, K::Int64, iterations::Int64=Integer(1e7))
    da = pi/N
    c3 = zeros(N,N, K, K,K)
    c3counts = copy(c3)

    rotations = Matrix{Float64}[random_rotation(3) for k = 1: 10000]
    dq = dr(intensity)
    surf = getSurfaceVolume(intensity)

    for i = 1:iterations
        k1 = Int64(rand(1:K))
        k2 = Int64(rand(1:K))
        k3 = Int64(rand(1:K))

        a = rand()*2.0*pi
        b = rand()*2.0*pi

        rot = rotations[rand(1:length(rotations))]
        p1 = k1*dq*rot*[0,1,0]
        p2 = k2*dq*rot*[cos(a) sin(a) 0; -sin(a) cos(a) 0; 0 0 0]*[0,1,0]
        p3 = k3*dq*rot*[cos(b) sin(b) 0; -sin(b) cos(b) 0; 0 0 0]*[0,1,0]

        alpha,beta = mod(angle_between(p1,p2), pi),mod(angle_between(p1,p3), pi)
        ai,bi = Int64(mod(floor(Int64, a/da),N)+1),Int64(mod(floor(Int64, b/da),N)+1)

        @inbounds c3[ai,bi,k3,k2,k1] += real(getVolumeInterpolated(surf, p1))* real(getVolumeInterpolated(surf, p2))* real(getVolumeInterpolated(surf, p3))
        @inbounds c3counts[ai,bi,k3,k2,k1] += 1.0
    end
    #Removing the cout normalization doesn't change anything
    return Float64[c3counts[a,b,k3,k2,k1] > 0 ? c3[a,b,k3,k2,k1] /c3counts[a,b,k3,k2,k1] : 0.0  for a=1:N,b=1:N,k3=1:K,k2=1:K,k1=1:K]
end

function integrateShell_2pc(intensity::SphericalHarmonicsVolume, N::Int64, K::Int64, iterations::Int64=Integer(1e7))
    da = pi/N
    c2 = zeros(N, K, K)
    c2counts = copy(c2)

    rotations = Matrix{Float64}[random_rotation(3) for k = 1: 10000]
    dq = dr(intensity)
    surf = getSurfaceVolume(intensity)

    for i = 1:iterations
        k1 = Int64(rand(1:K))
        k2 = Int64(rand(1:K))

        a = rand()*2.0*pi
        rot = rotations[rand(1:length(rotations))]
        p1 = k1*dq*rot*[0,1,0]
        p2 = k2*dq*rot*[cos(a) sin(a) 0; -sin(a) cos(a) 0; 0 0 0]*[0,1,0]

        ai = Int64(mod(floor(Int64, a/da),N)+1)

        @inbounds c2[ai,k2,k1] += real(getVolumeInterpolated(surf, p1))* real(getVolumeInterpolated(surf, p2))
        @inbounds c2counts[ai,k2,k1] += 1.0
    end
    return c2./c2counts
end

function integrateShell_3pc_alt(intensity::SphericalHarmonicsVolume, N::Int64, K::Int64, iterations::Int64=Integer(1e7))
    rotations = Matrix{Float64}[random_rotation(3) for k = 1: 300]
    dq = dr(intensity)
    surf = getSurfaceVolume(intensity)
    da = pi/N
    c3 = zeros(N,N, K, K, K)
    #     c3counts = copy(c3)

    up = Float64[0.0, 1.0, 0.0]

    for k1::Int64=1:1K
        p1 = k1*dq*up
        for k2::Int64=1:K
            for k3::Int64=1:K

                for ai::Int64 = 1:N
                    #                     ai = Int64(mod(floor(Int64, a/da),N)+1)
                    a = Float64(da*ai-da/2)
                    p2 = k2*dq*[cos(a) sin(a) 0; -sin(a) cos(a) 0; 0 0 0]*up

                    for bi::Int64 = 1:N
                        #                         bi = Int64(mod(floor(Int64, b/da),N)+1)
                        b = Float64(da*bi-da/2)
                        p3 = k3*dq*[cos(b) sin(b) 0; -sin(b) cos(b) 0; 0 0 0]*up

                        sum = 0.0
                        for k = 1:length(rotations)
                            rot = rotations[k]
                            sum += real(getVolumeInterpolated(surf, rot*p1))* real(getVolumeInterpolated(surf, rot*p2))* real(getVolumeInterpolated(surf, rot*p3))
                        end
                        c3[ai, bi, k3, k2, k1] = sum
                        #                         c3counts[ai, bi, k3, k2, k1]
                    end
                end
            end
        end
    end
    return c3 / sum(abs,c3)
end

"""Integrate out a single shell for comparison with theory"""
function integrate_c3_shell(intensity::SphericalHarmonicsVolume, k1::Int64, k2::Int64, k3::Int64, N::Int64, lambda::Float64, iterations::Int64=Integer(1e6))
    da = pi/N
    mdq = dq(intensity)
    surf = getSurfaceVolume(intensity)

    c3,c3counts = @sync @parallel ( (a,b) -> (a[1]+b[1], a[2]+b[2])) for i = 1:nworkers()
        c3_local = zeros(Float64,N,2*N)
        c3counts_local = zeros(Float64,N,2*N)
        rotations = Matrix{Float64}[random_rotation(3) for k = 1: 10000]

        for j = 1:iterations

            a = rand()*2.0*pi
            b = rand()*2.0*pi

            rot = rotations[rand(1:length(rotations))]
            p1 = k1*mdq*[0,1,0]
            p2 = k2*mdq*[cos(a) sin(a) 0; -sin(a) cos(a) 0; 0 0 0]*[0,1,0]
            p3 = k3*mdq*[cos(b) sin(b) 0; -sin(b) cos(b) 0; 0 0 0]*[0,1,0]

            alpha,beta = angle_between(p1,p2),angle_between(p1,p3)
            if alpha > pi
                alpha = 2*pi - alpha
                beta = 2*pi - beta
            end

            ai,bi = Int64(mod(floor(Int64, alpha/da),N)+1),Int64(mod(floor(Int64, beta/da),2*N)+1)

            p1e = rot*detector_to_Ewald_sphere(p1, lambda)
            p2e = rot*detector_to_Ewald_sphere(p2, lambda)
            p3e = rot*detector_to_Ewald_sphere(p3, lambda)

            @inbounds c3_local[ai,bi] += real(getVolumeInterpolated(surf, p1e))* real(getVolumeInterpolated(surf, p2e))* real(getVolumeInterpolated(surf, p3e))
            @inbounds c3counts_local[ai,bi] += 1.0

        end
        (c3_local,c3counts_local)
    end
    #Removing the cout normalization doesn't change anything
    return sdata(c3)./sdata(c3counts)
end

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

function c2_difference(a::C2, b::C2, K_range::UnitRange{Int64})
    sum(abs, sum(a[:,k2,k1]/sum(abs, a[:,k2,k1]) - b[:,k2,k1]/sum(abs, b[:,k2,k1]) for k1=K_range for k2=minimum(K_range):k1))/length(K_range)
end

function renormalize_correlation(c2::C2)
    N,K2,_ = Base.size(c2)
    c2_sym = deepcopy(c2)
    for k1=1:K2
        for k2=1:K2
            c2_sym[:,k2,k1] = sum(abs, c2[:,k2,k1]) > eps() ? c2[:,k2,k1]/sum(abs, c2[:,k2,k1]) : ones(N)
        end
    end
    return c2_sym
end
function renormalize_correlation(c3::C3)
    _,_,K3,_,_ = Base.size(c3)
   c3_sym = deepcopy(c3)
    for k1=1:K3
        for k2=1:K3
            for k3=1:K3
                c3_sym[:,:,k3,k2,k1] = c3[:,:,k3,k2,k1]/sum(abs, c3[:,:,k3,k2,k1])# + 1.0
            end
        end
    end
    return c3_sym
end

function symmetrize_correlation(c2::C2)
    _,K2,_ = Base.size(c2)
    c2_sym = deepcopy(c2)
    for k1=1:K2
        for k2=1:K2
            c2_sym[:,k2,k1] = 0.5*(c2[:,k2,k1] + reverse(c2[:,k2,k1]))
        end
    end
    return c2_sym
end

function symmetrize_correlation(c3::C3)
    _,_,K3,_,_ = Base.size(c3)
    c3_sym = deepcopy(c3)
    for k1=1:K3
        for k2=1:K3
            for k3=1:K3
                c3_sym[:,:,k3,k2,k1] = 0.5*(c3[:,:,k3,k2,k1] + rot180(c3[:,:,k3,k2,k1]))
            end
        end
    end
    return c3_sym
end

function add_Gaussian_filter(c2::C2, sigma::Float64=1.0)
    filtered_c2 = deepcopy(c2)
    _,K2,_ = Base.size(c2)
    for k1=1:K2
        for k2=1:K2
            filtered_c2[:,k2,k1] = gaussian_filter(c2[:,k2,k1], sigma)
        end
    end
    return filtered_c2
end

function add_Gaussian_filter(c3::C3, sigma::Float64=1.0)
    c3_filtered = deepcopy(c3)
    _,_,K3,_,_ = Base.size(c3)
    for k1=1:K3
        for k2=1:k1
            for k3=1:k2
                c3_filtered[:,:,k3,k2,k1] = gaussian_filter(c3_filtered[:,:,k3,k2,k1], 1.0)
            end
        end
    end
    return c3_filtered
end

"""Corrects the two-photon and three-photon correlation with the beamstop and adds a Gauss filter to compensate for pixel noise"""
function correct_beamstop(c2::C2, c3::C3, c2_beamstop::C2, c3_beamstop::C3)
    #Renormalize c2
    c2_corrected = c2 ./ renormalize_correlation(c2_beamstop)

    #Renormalize c3
    c3_corrected = c3 ./ renormalize_correlation(c3_beamstop)

    return c2_corrected, c3_corrected
end


text:
"""
    ref_int(a=1, b=1; tmax=10000, rtol=1E-12)

Calculates Coulomb repulsion between two Gaussian charge distributions
that have same origin.

# Arguments
- `a=1`  :   with of first Gaussian exp(-ax^2)
- `b=1`  :   with of second Gaussian exp(-bx^2)

# Keywords
- `tmax=10000`   :  maximum t-value for integration
- `rtol=1E-12`   :  relative tolerance for t-integration
"""
function ref_int(a=1, b=1; tmax=10000, rtol=1E-12)
    _f(t) = 1/sqrt( (a+b)*t^2 + a*b)
    norm = a*b/π
    return 2*√π * norm * quadgk(_f, 0, tmax; rtol=rtol)[1]
end

function test_eri_accuracy(d, ne, np; a=1, b=1)
    basis = Basis(-0.5*d, 0.5*d, ne, np)
    ρ1 = exp.(-a.*basis.^2) .* sqrt(a/π)
    ρ2 = exp.(-b.*basis.^2) .* sqrt(b/π)
    ω = get_weight(basis)
    w = diagm(ω)
    l = length(basis)

    J = zeros(l,l)
    for i in 1:l, j in 1:l
        J[i,j] = eri(basis, i, j)
    end
    cal = ρ2'*w*J*ρ1

    ref = ref_int(a, b)

    return Dict("ref"=>ref, "cal"=>cal)
end

max_stars_repo_name: tjjarvinen/OneDimensionalCorrelation

max_issues_repo_name: tjjarvinen/OneDimensionalCorrelation

max_forks_repo_name: tjjarvinen/OneDimensionalCorrelation


text:
abstract type AbstractGRUESN <: AbstractEchoStateNetwork end

struct GRUESN{T, S<:AbstractArray{T}, I, B, F, N, G} <: AbstractGRUESN
    res_size::I
    in_size::I
    out_size::I
    train_data::S
    alpha::T
    nla_type::N
    activation::F
    W::S
    W_in::S
    gates::G
    states::S
    extended_states::B
end

"""
    GRUESN(W::AbstractArray{T}, train_data::AbstractArray{T}, W_in::AbstractArray{T} [, gates_weight::T, activation::Any, alpha::T, nla_type::NonLinearAlgorithm, extended_states::Bool])

Return a Gated Recurrent Unit [1] ESN struct.

[1] Cho, Kyunghyun, et al. “Learning phrase representations using RNN encoder-decoder for statistical machine translation.” arXiv preprint arXiv:1406.1078 (2014).
"""
function GRUESN(W::AbstractArray{T},
        train_data::AbstractArray{T},
        W_in::AbstractArray{T};
        gates_weight::T = 0.9,
        activation::Any = tanh,
        alpha::T = 1.0,
        nla_type::NonLinearAlgorithm = NLADefault(),
        extended_states::Bool = false) where T<:AbstractFloat

    in_size = size(train_data, 1)
    out_size = size(train_data, 1)
    res_size = size(W, 1)

    if size(W_in, 1) != res_size
        throw(DimensionMismatch("size(W_in, 1) must be equal to size(W, 1)"))
    elseif size(W_in, 2) != in_size
        throw(DimensionMismatch("size(W_in, 2) must be equal to in_size"))
    end

    gates = GRUGates(res_size, in_size, gates_weight)
    states = gru_states(W, W_in, train_data, alpha, activation, extended_states, gates_weight, gates)

    return GRUESN{T, typeof(train_data),
        typeof(res_size),
        typeof(extended_states),
        typeof(activation),
        typeof(nla_type),
        typeof(gates)}(res_size, in_size, out_size, train_data,
    alpha, nla_type, activation, W, W_in, gates, states, extended_states)
end

"""
    GRUESNpredict(esn::AbstractGRUESN, predict_len::Int, W_out::AbstractArray{Float64})

Return the prediction for a given length of the constructed GRUESN.
"""
function GRUESNpredict(esn::AbstractGRUESN,
    predict_len::Int,
    W_out::AbstractArray{Float64})

    output = zeros(Float64, esn.in_size, predict_len)
    x = esn.states[:, end]

    if esn.extended_states == false
        for i=1:predict_len
            x_new = nla(esn.nla_type, x)
            out = (W_out*x_new)
            output[:, i] = out
            x = gru(esn.gates, esn.activation, esn.W, esn.W_in, x, out)
        end
    elseif esn.extended_states == true
        for i=1:predict_len
            x_new = nla(esn.nla_type, x)
            out = (W_out*x_new)
            output[:, i] = out
            x = vcat(gru(esn.gates, esn.activation, esn.W, esn.W_in, x[1:esn.res_size], out), out)
        end
    end
    return output
end

abstract type AbstractGates end

struct GRUGates{T, S<:AbstractArray{T}} <: AbstractGates
    U_r::S
    W_r::S
    b_r::S
    U_z::S
    W_z::S
    b_z::S
end

function GRUGates(res_size, in_size, gates_weight)
    U_r = irrational_sign_input(res_size, in_size, gates_weight, start = res_size*in_size)
    W_r = irrational_sign_input(res_size, res_size, gates_weight, start = 2*res_size*in_size)
    b_r = irrational_sign_input(res_size, 1, gates_weight, start = res_size*(2*in_size+res_size))
    U_z = irrational_sign_input(res_size, in_size, gates_weight, start = res_size*(2*in_size+res_size+1))
    W_z = irrational_sign_input(res_size, res_size, gates_weight, start = res_size*(3*in_size+res_size+1))
    b_z = irrational_sign_input(res_size, 1, gates_weight, start = res_size*(2*in_size+2*res_size+1))

    return GRUGates{typeof(gates_weight), typeof(U_r)}(U_r, W_r, b_r, U_z, W_z, b_z)
end

function gru(gates, activation, W, W_in, x, y)
    reset_gate = gates.U_r*y + gates.W_r*x+gates.b_r
    reset_gate = 1 ./(1 .+exp.(-reset_gate))
    update = activation.(W_in*y+W*(reset_gate.*x))
    update_gate = gates.U_z*y + gates.W_z*x + gates.b_z
    update_gate = 1 ./(1 .+exp.(-update_gate))
    return update_gate.*update + (1 .- update_gate) .*x
end

function gru_states(W::AbstractArray{Float64},
        W_in::AbstractArray{Float64},
        train_data::AbstractArray{Float64},
        alpha::Float64,
        activation::Function,
        extended_states::Bool,
        gates_weight::Float64,
        gates::AbstractGates)

    train_len = size(train_data, 2)
    res_size = size(W, 1)
    in_size = size(train_data, 1)
    states = zeros(Float64, res_size, train_len)

    for i=1:train_len-1
        states[:, i+1] = gru(gates, activation, W, W_in, states[:, i], train_data[:, i])
    end

    if extended_states == true
        ext_states = vcat(states, hcat(zeros(Float64, in_size), train_data[:,1:end-1]))
        return ext_states
    else
        return states
    end
end


text:
using Base.Test
include("splinesParallel.jl")

@testset "splines" begin
	point = [0.5, 0.7, 1.2]
	V = [1.0,3.4,5.5,6.9,2.4,8.0];
	FV = [[1,2,4],[2,3,5],[3,5,6],[4,5,7],[5,7,8],[6,8,9]];
	@testset "sel" begin
 		@test typeof(sel(1)(V))==Float64
		@test typeof(sel(1)(FV))==Array{Int64,1}
		@test sel(4)(V)==6.9
		@test sel(5)(FV)==[5,7,8]
	end
	@testset "larBernsteinBasis" begin
		@test typeof(larBernsteinBasis(sel(1),3)([0.5,2,7.3]))==Array{Float64,1}
		@test typeof(larBernsteinBasis(sel(5), 6)(V))==
		Array{Float64,1}
		H=[1.2,3.4,5.6,7.8,9.0]
		F=larBernsteinBasis(sel(4),5)(H)
		@test length(F)==6
	end

	@testset "plarTensorProdSurface" begin
		controlpoints = [[[0,0,0],[2,-4,2]],[[0,3,1],[4,0,0]]]
		basis = larBernsteinBasis(sel(1), 1)
		point = [0.5, 2, 7.3]
		T=plarTensorProdSurface([basis, basis], controlpoints)(point)
		@test typeof(T)==SharedArray{Float64,1}
		@test length(T)==3
	end

	@testset "plarBilinearSurface" begin
		controlpoints = [[[0,0,0],[2,-4,2]],[[0,3,1],[4,0,0]]]
		Bil=plarBilinearSurface(controlpoints)([0.5,0.7])
		@test typeof(Bil)==SharedArray{Float64,1}
		@test length(Bil)==3
	end

	@testset "plarBiquadraticSurface" begin
		controlpoints = [[[0,0,0],[2,0,1],[3,1,1]],[[1,3,-1],[2,2,0],[3,2,0]],[[-2,4,0],[2,5,1],[1,3,2]]]
		mapping = plarBiquadraticSurface(controlpoints)([0.7, 0.2])
		@test typeof(mapping)==SharedArray{Float64,1}
		@test length(mapping)==3
	end

	@testset "plarBicubicSurface" begin
		controlpoints = [[[0,0,0], [0,3,4], [0,6,3], [0,10,0]],[[3,0,2], [2,2.5,5], [3,6,5], [4,8,2]],[[6,0,2], [8,3,5], [7,6,4.5], [6,10,2.5]],[[10,0,0], [11,3,4], [11,6,3], [10,9,0]]]
		mapping = plarBicubicSurface(controlpoints)([0.7,0.2])
		@test typeof(mapping)==SharedArray{Float64,1}
		@test length(mapping)==3
	end

	@testset "plarBezier"begin
		controldata = [[0,1],[0,0],[1,1],[1,0]]
		U = sel(2)
		@test typeof(plarBezier(U, controldata)(point))==
		SharedArray{Float64,1}
		G=plarBezier(U, controldata)(point)
		@test length(G)==2
		@test typeof(plarBezier(sel(3), controldata)(point))==
		SharedArray{Float64,1}
		F=plarBezier(sel(3), controldata)(point)
		@test length(F)==2
	end

	@testset "plarBezierCurve" begin
		controldata = [[1,2,-5],[2,-4,2]]
		@test typeof(plarBezierCurve(controldata)(point))==SharedArray{Float64,1}
		@test length(plarBezierCurve(controldata)(point))==3
	end

	@testset "bsplineBasis" begin
		knots = [0,0,0,1,1,2,2,3,3,4,4,4]
		ncontrols = 9
		degree = 2
		H=bsplineBasis(degree, knots, ncontrols)([0.7,2.5,3.0])
		@test typeof(bsplineBasis(degree, knots, ncontrols)([0.7,2.5,3.0]))==Array{Float64,1}
		@test length(bsplineBasis(degree, knots, ncontrols)([0.7,2.5,3.0]))==9
	end

	@testset "plarCoonsPatch" begin
		Su0 = plarBezier(sel(1), [[0,0,0],[10,0,0]])
		Su1 = plarBezier(sel(1), [[0,10,0],[2.5,10,3],[5,10,-3],[7.5,10,3],[10,10,0]])
		Sv0 = plarBezier(sel(1), [[0,0,0],[0,0,3],[0,10,3],[0,10,0]])
		Sv1 = plarBezier(sel(2), [[10,0,0],[10,5,3],[10,10,0]])
		point = [0.7, 0.2]
		out = plarCoonsPatch([Su0,Su1,Sv0,Sv1])(point)
		@test typeof(out)==Array{Float64,1}
	end

	@testset "pt_bspline"begin
		b1 = plarBezier(S1, [[0,1,0],[0,1,5]])
		b2 = plarBezier(S1, [[0,0,0],[0,0,5]])
		b3 = plarBezier(S1, [[1,0,0],[2,-1,2.5],[1,0,5]])
		b4 = plarBezier(S1, [[1,1,0],[1,1,5]])
		b5 = plarBezier(S1, [[0,1,0],[0,1,5]])
		controls = [b1, b2, b3, b4, b5]
		knots = [0,1,2,3,4,5,6,7]           # periodic B-spline
		knots = [0,0,0,1,2,3,3,3]           # non-periodic B-spline
		ptbspline = pt_bspline(S2, 2, knots, controls)
		B=ptbspline([0.7, 0.2])
		@test typeof(B)==SharedArray{Float64,1}
		@test length(B)==3
	end

	@testset "pt_nurbs" begin
		knots = [0,0,0,1,1,2,2,3,3,4,4,4]
		p=sqrt(2)/2.0
		controls = [[-1,0,1], [-p,p,p], [0,1,1], [p,p,p],[1,0,1], [p,-p,p],[0,-1,1], [-p,-p,p], [-1,0,1]]
		c1 = plarBezier(S1, [[-1,0,0,1],[-1,0,1,1]])
		c2 = plarBezier(S1, [[-p,p,0,p],[-p,p,p,p]])
		c3 = plarBezier(S1, [[0,1,0,1],[0,1,1,1]])
		c4 = plarBezier(S1, [[p,p,0,p],[p,p,p,p]])
		c5 = plarBezier(S1, [[1,0,0,1],[1,0,1,1]])
		c6 = plarBezier(S1, [[p,-p,0,p],[p,-p,p,p]])
		c7 = plarBezier(S1, [[0,-1,0,1],[0,-1,1,1]])
		c8 = plarBezier(S1, [[-p,-p,0,p],[-p,-p,p,p]])
		c9 = plarBezier(S1, [[-1,0,0,1],[-1,0,1,1]])
		controls = [c1,c2,c3,c4,c5,c6,c7,c8,c9]
		tnurbs = pt_nurbs(S2,2,knots,controls)
		N=tnurbs([0.7, 0.2])
		@test typeof(N)==Array{Float64,1}
		@test length(N)==3
	end

end


text:
export AbstractTrajectory, get_trace, RTSA, SARTSA

"""
    AbstractTrajectory{names,types} <: AbstractArray{NamedTuple{names,types},1}

A trajectory is used to record some useful information
during the interactions between agents and environments.

# Parameters
- `names`::`NTuple{Symbol}`, indicate what fields to be recorded.
- `types`::`Tuple{DataType...}`, the datatypes of `names`.

The length of `names` and `types` must match.

Required Methods:

- [`get_trace`](@ref)
- `Base.push!(t::AbstractTrajectory, kv::Pair{Symbol})`
- `Base.pop!(t::AbstractTrajectory, s::Symbol)`

Optional Methods:

- `Base.length`
- `Base.size`
- `Base.lastindex`
- `Base.isempty`
- `Base.empty!`
"""
abstract type AbstractTrajectory{names,types} <: AbstractArray{NamedTuple{names,types},1} end

# some typical trace names
"An alias of `(:reward, :terminal, :state, :action)`"
const RTSA = (:reward, :terminal, :state, :action)

"An alias of `(:state, :action, :reward, :terminal, :next_state, :next_action)`"
const SARTSA = (:state, :action, :reward, :terminal, :next_state, :next_action)

"""
    get_trace(t::AbstractTrajectory, s::NTuple{N,Symbol}) where {N}
"""
get_trace(t::AbstractTrajectory, s::NTuple{N,Symbol}) where {N} =
    NamedTuple{s}(get_trace(t, x) for x in s)

"""
    get_trace(t::AbstractTrajectory, s::Symbol...)
"""
get_trace(t::AbstractTrajectory, s::Symbol...) = get_trace(t, s)

"""
    get_trace(t::AbstractTrajectory{names}) where {names}
"""
get_trace(t::AbstractTrajectory{names}) where {names} =
    NamedTuple{names}(get_trace(t, x) for x in names)

Base.length(t::AbstractTrajectory) = maximum(length(x) for x in get_trace(t))
Base.size(t::AbstractTrajectory) = (length(t),)
Base.lastindex(t::AbstractTrajectory) = length(t)
Base.getindex(t::AbstractTrajectory{names,types}, i::Int) where {names,types} =
    NamedTuple{names,types}(Tuple(x[i] for x in get_trace(t)))

Base.isempty(t::AbstractTrajectory) = all(isempty(t) for t in get_trace(t))

function Base.empty!(t::AbstractTrajectory)
    for x in get_trace(t)
        empty!(x)
    end
end

"""
    Base.push!(t::AbstractTrajectory; kwargs...)
"""
function Base.push!(t::AbstractTrajectory; kwargs...)
    for kv in kwargs
        push!(t, kv)
    end
end

"""
    Base.pop!(t::AbstractTrajectory{names}) where {names}
`pop!` out one element of each trace in `t`
"""
function Base.pop!(t::AbstractTrajectory{names}) where {names}
    pop!(t, names...)
end

"""
    Base.pop!(t::AbstractTrajectory, s::Symbol...)
`pop!` out one element of the traces specified in `s`
"""
function Base.pop!(t::AbstractTrajectory, s::Symbol...)
    NamedTuple{s}(pop!(t, x) for x in s)
end


text:
function load(io::IO, ::Type{IAEAHeader})
    read_header(io)
end

const MANDATORY_KEYWORDS = [
:IAEA_INDEX,
:TITLE,
:FILE_TYPE,
:CHECKSUM,
:RECORD_CONTENTS,
:RECORD_CONSTANT,
:RECORD_LENGTH,
:BYTE_ORDER,
:TRANSPORT_PARAMETERS,
:MACHINE_TYPE,
:MONTE_CARLO_CODE_VERSION,
:GLOBAL_PHOTON_ENERGY_CUTOFF,
:GLOBAL_PARTICLE_ENERGY_CUTOFF,
:COORDINATE_SYSTEM_DESCRIPTION,
]
# $ORIG_HISTORIES:

function validate_header_dict(d)
    for key in MANDATORY_KEYWORDS
        @argcheck key ∈ keys(d)
    end
    ft = parse(Int, stripcomments(d[:FILE_TYPE]))
    if ft == 0 # phase space
        @argcheck :ORIG_HISTORIES ∈ keys(d)
        @argcheck :PARTICLES ∈ keys(d)
    elseif ft == 1 # event generator
        @argcheck :INPUT_FILE_FOR_EVENT_GENERATOR ∈ keys(d)
    end
end

const RKEY = r"\$(.*):"

iskeyline(s) = occursin(RKEY, s)
isemptyline(s) = strip(s) == ""
parsekey(line) = Symbol(match(RKEY, line)[1])
stripcomments(s) = strip(first(split(s, "//")))

function read_header_dict(io::IO; validate::Bool=false)
    d = OrderedDict{Symbol, String}()
    val_lines = String[]
    line = readline(io)
    key = parsekey(line)
    while !eof(io)
        line = readline(io)
        if iskeyline(line)
            val = join(val_lines, '\n')
            d[key] = val
            key = parsekey(line)
            empty!(val_lines)
        else
            push!(val_lines, line)
        end
    end
    d[key] = join(val_lines, '\n')
    @assert eof(io)
    if validate
        validate_header_dict(d)
    end
    d
end

function cleanup_record(s)
    lines = split(s, '\n')
    lines = map(stripcomments, lines)
    filter(l -> !isemptyline(l), lines)
end

function read_header(io::IO)
    d = read_header_dict(io)
    IAEAHeader(d)
end

function IAEAHeader(d::AbstractDict)

    contents  = cleanup_record(d[:RECORD_CONTENTS])
    constants = cleanup_record(d[:RECORD_CONSTANT])

    read_next!(T, c) = parse(T, popfirst!(c))
    read_next!(::Type{Bool}, c) = Bool(parse(Int, popfirst!(c)))
    function add_next_default!(d,key,contents, constants)
        stored_in_phsp = read_next!(Bool, contents)
        stored_in_header = !stored_in_phsp
        T = Float32
        if stored_in_header
            d[key] = read_next!(T, constants)
        end
    end

    rcd = OrderedDict()
    add_next_default!(rcd,:x,contents, constants)
    add_next_default!(rcd,:y,contents, constants)
    add_next_default!(rcd,:z,contents, constants)
    add_next_default!(rcd,:u,contents, constants)
    add_next_default!(rcd,:v,contents, constants)
    add_next_default!(rcd,:w,contents, constants)
    @assert !haskey(rcd, :w)
    add_next_default!(rcd, :weight, contents, constants)
    Nf = read_next!(Int, contents)
    Ni = read_next!(Int, contents)
    @assert isempty(constants)
    rc = (;rcd...)
    raw = OrderedDict{Symbol, String}(
        (k=>chomp(v)) for (k,v) in pairs(d) if
        !(k in [:RECORD_CONSTANT, :RECORD_CONSTANT]))

    # NT = typeof(rc)
    IAEAHeader{Nf, Ni}(rc, raw)
end


text:
Random.seed!(1234)
N, M = 50, 200
G = erdos_renyi(N, M)
# println("Is the ER graph connected: ", is_connected(G))

L = Matrix(laplacian_matrix(G))
𝛌, 𝚽 = eigen(L); standardize_eigenvectors!(𝚽)
Q = incidence_matrix(G; oriented = true)
∇𝚽 = Q' * 𝚽


text:
"""
    cartesian_sample_from_logits(logits; dims=1)

Returns an array `X` of indices sampled from `CartesianIndices(logits)` along
dimensions `dims`, with probabilities `P = softmax(logits; dims)`.
In particular, dimensions `dims` of `X` are singleton.
"""
function cartesian_sample_from_logits(logits::AbstractArray; dims=1)
    return cartesian_sample_from_logits_gumbel(logits; dims)
end

# cartesian_sample_from_logits using Gumbel trick
function cartesian_sample_from_logits_gumbel(logits::AbstractArray; dims = 1)
    z = logits .+ randgumbel.(float(eltype(logits)))
    return argmax(z; dims)
end

"""
    randgumbel(T = Float64)

Generates a random Gumbel variate.
"""
randgumbel(::Type{T} = Float64) where {T} = -log(Random.randexp(T))


text:
struct JuMPProblem{E, M, V}
    evaluator::E
    model::M
    vars::V
end

struct JuMPEvaluator{XL, XU, X, CL, CU, O, C, G, J, JS, H, HS} <: MOI.AbstractNLPEvaluator
    nvars::Int
    xlb::XL
    xub::XU
    x0::X
    nconstr::Int
    clb::CL
    cub::CU
    njacvalues::Int
    nhessvalues::Int
    obj::O
    eval_g::C
    eval_grad_f::G
    eval_jac_g::J
    jac_structure::JS
    eval_h::H
    lag_hess_structure::HS
end

function get_jump_problem(
    model::VecModel, x0 = getinit(model);
    optimizer, first_order, kwargs...,
)
    integers = model.integer
    eq = if length(model.eq_constraints.fs) == 0
        nothing
    else
        model.eq_constraints
    end
    ineq = if length(model.ineq_constraints.fs) == 0
        nothing
    else
        model.ineq_constraints
    end
    obj = CountingFunction(getobjective(model))
    return get_jump_problem(
        obj, ineq, eq, x0, integers, getmin(model),
        getmax(model), first_order, optimizer,
    ), obj.counter
end
function get_jump_problem(
    obj, ineq_constr, eq_constr, x0, integers,
    xlb, xub, first_order, optimizer,
)
    nvars = 0
    if ineq_constr !== nothing
        ineqJ0 = Zygote.jacobian(ineq_constr, x0)[1]
        ineq_nconstr, nvars = size(ineqJ0)
        Joffset = nvalues(ineqJ0)
    else
        ineqJ0 = nothing
        ineq_nconstr = 0
        Joffset = 0
    end
    if eq_constr !== nothing
        eqJ0 = Zygote.jacobian(eq_constr, x0)[1]
        eq_nconstr, nvars = size(eqJ0)
    else
        eqJ0 = nothing
        eq_nconstr = 0
    end
    njacvals = nvalues(ineqJ0) + nvalues(eqJ0)
    @assert nvars > 0
    lag(factor, y) = x -> begin
        factor * obj(x) +
            _dot(ineq_constr, x, @view(y[1:ineq_nconstr])) +
            _dot(eq_constr, x, @view(y[ineq_nconstr+1:end]))
    end
    clb = [fill(-Inf, ineq_nconstr); zeros(eq_nconstr)]
    cub = zeros(ineq_nconstr + eq_nconstr)

    function eval_g(x::AbstractVector{Float64}, g::AbstractVector{Float64})
        if ineq_constr !== nothing
            g[1:ineq_nconstr] .= ineq_constr(x)
        end
        if eq_constr !== nothing
            g[ineq_nconstr+1:end] .= eq_constr(x)
        end
        return g
    end
    function eval_grad_f(x::AbstractVector{Float64}, grad_f::AbstractVector{Float64})
        grad_f .= Zygote.gradient(obj, x)[1]
    end
    function jac_structure()
        rows = fill(0, njacvals)
        cols = fill(0, njacvals)
        ineqJ0 === nothing || fill_indices!(rows, cols, ineqJ0)
        eqJ0 === nothing || fill_indices!(rows, cols, eqJ0, offset = Joffset, row_offset = ineq_nconstr)
        return tuple.(rows, cols)
    end
    function eval_jac_g(x::AbstractVector{Float64}, values::AbstractVector{Float64})
        values .= 0
        if ineq_constr !== nothing
            ineqJ = Zygote.jacobian(ineq_constr, x)[1]
            add_values!(values, ineqJ)
        end
        if eq_constr !== nothing
            eqJ = Zygote.jacobian(eq_constr, x)[1]
            add_values!(values, eqJ, offset = Joffset)
        end
        return values
    end

    if first_order
        eval_h = (x...) -> 0.0
        lag_hess_structure = () -> Tuple{Int,Int}[]
        Hnvalues = 0
    else
        HL0 = LowerTriangular(
            Zygote.hessian(
                lag(1.0, ones(ineq_nconstr + eq_nconstr)),
                x0,
            ),
        )
        Hnvalues = nvalues(HL0)
        lag_hess_structure = function ()
            rows = fill(0, Hnvalues)
            cols = fill(0, Hnvalues)
            fill_indices!(rows, cols, HL0)
            return tuple.(rows, cols)
        end
        eval_h = function (x::AbstractVector{Float64}, obj_factor::Float64, lambda::AbstractVector{Float64}, values::AbstractVector{Float64})
            HL = LowerTriangular(
                Zygote.hessian(lag(obj_factor, lambda), x),
            )
            values .= 0
            add_values!(values, HL)
            return values
        end
    end

    evaluator = JuMPEvaluator(
        nvars, xlb, xub, x0, ineq_nconstr + eq_nconstr, clb, cub,
        njacvals, Hnvalues, obj, eval_g, eval_grad_f, eval_jac_g,
        jac_structure, eval_h, lag_hess_structure,
    )
    jump_model = JuMP.Model(optimizer)
    moi_model = JuMP.backend(jump_model)
    MOI.empty!(moi_model)
    vars = map(1:nvars) do i
        v = MOI.add_variable(moi_model)
        if integers[i]
            if xlb[i] > -1 && xub[i] < 2
                MOI.add_constraint(moi_model, MOI.SingleVariable(v), MOI.ZeroOne())
            else
                MOI.add_constraint(moi_model, MOI.SingleVariable(v), MOI.Integer())
            end
        end
        if xub[i] != Inf
            MOI.add_constraint(
                moi_model,
                MOI.SingleVariable(v),
                MOI.LessThan(xub[i]),
            )
        end
        if xlb[i] != Inf
            MOI.add_constraint(
                moi_model,
                MOI.SingleVariable(v),
                MOI.GreaterThan(xlb[i]),
            )
        end
        MOI.set(moi_model, MOI.VariablePrimalStart(), v, x0[i])
        return v
    end
    block_data = MOI.NLPBlockData(MOI.NLPBoundsPair.(clb, cub), evaluator, true)
    MOI.set(moi_model, MOI.NLPBlock(), block_data)
    MOI.set(moi_model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    return JuMPProblem(evaluator, jump_model, vars)
end

function MOI.initialize(d::JuMPEvaluator, requested_features::Vector{Symbol})
    for feat in requested_features
        if !(feat in MOI.features_available(d))
            error("Unsupported feature $feat")
            # TODO: implement Jac-vec products for solvers that need them
        end
    end
end

function MOI.features_available(::JuMPEvaluator)
    return [:Grad, :Jac, :Hess]
end

MOI.eval_objective(d::JuMPEvaluator, x) = d.obj(x)

function MOI.eval_constraint(d::JuMPEvaluator, g, x)
    d.eval_g(x, g)
    return
end

function MOI.eval_objective_gradient(d::JuMPEvaluator, grad_f, x)
    d.eval_grad_f(x, grad_f)
    return
end

function MOI.jacobian_structure(d::JuMPEvaluator)
    return d.jac_structure()
end
# lower triangle only
function MOI.hessian_lagrangian_structure(d::JuMPEvaluator)
    return d.lag_hess_structure()
end

function MOI.eval_constraint_jacobian(d::JuMPEvaluator, J, x)
    d.eval_jac_g(x, J)
    return
end

function MOI.eval_hessian_lagrangian(d::JuMPEvaluator, H, x, σ, μ)
    d.eval_h(x, σ, μ, H)
   return
end

nvalues(::Nothing) = 0
nvalues(J::Matrix) = length(J)
nvalues(H::LowerTriangular{<:Real, <:Matrix}) = (size(H, 1) + 1) * size(H, 1) ÷ 2
nvalues(H::SparseMatrixCSC) = length(H.nzval)

# Implement these for sparse matrices
function fill_indices!(rows, cols, J0::Matrix; offset = 0, row_offset = 0)
    nconstr, nvars = size(J0)
    for j in 1:nvars
        cols[offset + 1 : offset + nconstr] .= j
        rows[offset + 1 : offset + nconstr] .= row_offset+1:row_offset+nconstr
        offset += nconstr
    end
    return rows, cols
end
function fill_indices!(rows, cols, HL::LowerTriangular{<:Real, <:Matrix}; offset = 0, row_offset = 0)
    nvars = size(HL, 1)
    for j in 1:nvars
        cols[offset + 1 : offset + nvars - j + 1] .= j
        rows[offset + 1 : offset + nvars - j + 1] .= row_offset + j : row_offset + nvars
        offset += nvars - j + 1
    end
    return rows, cols
end
function fill_indices!(rows, cols, HL::SparseMatrixCSC; offset = 0, row_offset = 0)
    for col in 1:length(HL.colptr)-1
        indices = HL.colptr[col]:HL.colptr[col+1]-1
        nvars = length(indices)
        cols[offset + 1 : offset + nvars] .= col
        rows[offset + 1 : offset + nvars] = row_offset .+ HL.rowval[indices]
        offset += nvars
    end
    return rows, cols
end

function add_values!(values, J::Matrix; offset = 0)
    nvars = length(J)
    values[offset+1:offset+nvars] .+= vec(J)
    return values
end
function add_values!(values, HL::LowerTriangular{<:Real, <:Matrix}; factor = 1, offset = 0)
    nvars = size(HL, 1)
    for j in 1:nvars
        values[offset + 1 : offset + nvars - j + 1] .+= HL[j:nvars, j] .* factor
        offset += nvars - j + 1
    end
    return values
end
function add_values!(values, HL::SparseMatrixCSC; factor = 1, offset = 0)
    nvars = length(HL.nzval)
    values[offset+1:offset+nvars] .= HL.nzval .* factor
    return values
end

_dot(f, x, y) = dot(f(x), y)
_dot(::Nothing, ::Any, ::Any) = 0.0


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -48.987102771054815 <= q <= 45.974621434893855)

                     add_NL_constraint(m, :(softplus(-0.368688824745 + 0.553104331767496*softplus(-0.34679271547764046 + -0.03387429483407178*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.41639023650850415*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.7737560883943391*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.09930846098451918*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.9337943134926245*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + -0.17925952214788987*softplus(0.6682158984675199 + 0.6233864552377666*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.35483759096347534*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4760912088664373*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + -0.23820709887370928*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9242931300611921*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.4698400936211189*softplus(0.8534714115378583 + 0.6703763846356168*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.08490767686166478*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.5656859981927576*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.4980009846607585*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.014109452266772848*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + -0.7276938776175426*softplus(-0.09327714623589056 + -0.24359337663872616*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.2759577392076742*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4166615549845427*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.17410413013320136*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9439406543948436*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.4202042464647322*softplus(-0.09609361556421225 + -0.8712983668552394*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.8013603269787453*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.3524955117104245*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.28183131704908826*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.6019873452940772*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3])))) + softplus(-0.8761821794177092 + 0.4562914248397121*softplus(-0.34679271547764046 + -0.03387429483407178*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.41639023650850415*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.7737560883943391*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.09930846098451918*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.9337943134926245*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.3827259206071312*softplus(0.6682158984675199 + 0.6233864552377666*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.35483759096347534*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4760912088664373*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + -0.23820709887370928*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9242931300611921*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.20062384412841983*softplus(0.8534714115378583 + 0.6703763846356168*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.08490767686166478*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.5656859981927576*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.4980009846607585*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.014109452266772848*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + -0.28064116456943733*softplus(-0.09327714623589056 + -0.24359337663872616*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.2759577392076742*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4166615549845427*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.17410413013320136*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9439406543948436*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.4949171063694977*softplus(-0.09609361556421225 + -0.8712983668552394*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.8013603269787453*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.3524955117104245*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.28183131704908826*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.6019873452940772*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3])))) + softplus(0.7716665307334409 + 0.20199447224764766*softplus(-0.34679271547764046 + -0.03387429483407178*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.41639023650850415*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.7737560883943391*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.09930846098451918*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.9337943134926245*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.3172295658101083*softplus(0.6682158984675199 + 0.6233864552377666*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.35483759096347534*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4760912088664373*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + -0.23820709887370928*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9242931300611921*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.10629001518557368*softplus(0.8534714115378583 + 0.6703763846356168*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.08490767686166478*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.5656859981927576*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.4980009846607585*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.014109452266772848*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.4930771030063288*softplus(-0.09327714623589056 + -0.24359337663872616*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.2759577392076742*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4166615549845427*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.17410413013320136*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9439406543948436*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.5476037221052898*softplus(-0.09609361556421225 + -0.8712983668552394*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.8013603269787453*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.3524955117104245*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.28183131704908826*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.6019873452940772*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3])))) + softplus(0.6820199557891846 + 0.4693680154187745*softplus(-0.34679271547764046 + -0.03387429483407178*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.41639023650850415*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.7737560883943391*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.09930846098451918*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.9337943134926245*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + -0.6552396777765495*softplus(0.6682158984675199 + 0.6233864552377666*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.35483759096347534*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4760912088664373*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + -0.23820709887370928*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9242931300611921*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + -0.535415419107649*softplus(0.8534714115378583 + 0.6703763846356168*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.08490767686166478*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.5656859981927576*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.4980009846607585*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.014109452266772848*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + -0.022291893207070146*softplus(-0.09327714623589056 + -0.24359337663872616*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.2759577392076742*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4166615549845427*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.17410413013320136*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9439406543948436*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + -0.04141935716867051*softplus(-0.09609361556421225 + -0.8712983668552394*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.8013603269787453*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.3524955117104245*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.28183131704908826*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.6019873452940772*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3])))) + softplus(-0.779255880286021 + 0.8816016628889711*softplus(-0.34679271547764046 + -0.03387429483407178*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.41639023650850415*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.7737560883943391*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.09930846098451918*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.9337943134926245*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.3409151924818685*softplus(0.6682158984675199 + 0.6233864552377666*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.35483759096347534*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4760912088664373*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + -0.23820709887370928*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9242931300611921*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + -0.5165976261581497*softplus(0.8534714115378583 + 0.6703763846356168*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.08490767686166478*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.5656859981927576*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.4980009846607585*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.014109452266772848*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + -0.06878582640645536*softplus(-0.09327714623589056 + -0.24359337663872616*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + 0.2759577392076742*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + -0.4166615549845427*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.17410413013320136*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + 0.9439406543948436*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3]))) + 0.8469129987678201*softplus(-0.09609361556421225 + -0.8712983668552394*softplus(-0.3632300212500721 + -0.3374701832267393*$(x[1]) + 0.8756028312962347*$(x[2]) + 0.4387543479647671*$(x[3])) + -0.8013603269787453*softplus(0.09153485555592855 + 0.23042216426293383*$(x[1]) + 0.7064087974334008*$(x[2]) + -0.7947656255493469*$(x[3])) + 0.3524955117104245*softplus(-0.9772041331486054 + -0.49882047002052676*$(x[1]) + 0.9295268810329342*$(x[2]) + -0.9980546314915881*$(x[3])) + 0.28183131704908826*softplus(-0.6572599840304187 + 0.9275280812259474*$(x[1]) + 0.48175699367938973*$(x[2]) + 0.6558120152589955*$(x[3])) + -0.6019873452940772*softplus(-0.39686236988964074 + -0.053441112110775624*$(x[1]) + -0.8651914231914781*$(x[2]) + 0.0572006213281977*$(x[3])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using Adapt: Adapt, adapt
using Base: index_shape, index_dimsum, index_ndims, to_indices, front, tail
using Base: mightalias, unalias, dataids, unsafe_convert
using BenchmarkTools
using ElasticArrays
using MacroTools
using UnPack
using Random
using Test
using UnsafeArrays

using SpecialArrays
using SpecialArrays: TypedBool, True, False
using SpecialArrays: tuple_getindex, tuple_setindex, tuple_map, tuple_split, invert

include("util.jl")


text:
module rrtTypes
    export Node, Edge, Obstacle, Room, Point

    struct Point
        x::Int
        y::Int
    end

    struct Node
        id::Int
        iPrev::Int #parent
        state::Tuple{Point} #hrm, don't know how to restrict to two points.  i guess it's immutable so use tuple
    end
    Base.start(n::Node) = 1
    Base.done(n::Node, i) = nfields(n) < i
    Base.next(n::Node, i) = getfield(n, i), i+1

    struct Edge
        startnode::Int
        endnode::Int
    end

    struct Obstacle
        SW
        NE
    end

    ### Some types?
    struct Room #corners of the room
        x1::Int
        y1::Int #bottom left corner?
        x2::Int
        y2::Int
        obstacleList::Array{Obstacle}
    end

    struct robot
    end

end


text:
using Random: GLOBAL_RNG

"""
    AbstractNoiseSource

# Interface

To define your own source, simply provide the function for how signal is manipulated given an `AbstractRNG`

* `(::MySource)(rng, input)`
* `noise(::MySource, input)`
* `limits(::MySource)` (optional)

See [`ShotNoise`](@ref) for an example of an implementation of the interface
"""
abstract type AbstractNoiseSource <: Function end

"""
    (::AbstractNoiseSource)([rng], input)

Process a signal following the transfer function defined by the noise source. A random number generator can be passed if the process is stochastic.
"""
(source::AbstractNoiseSource)(input) = source(GLOBAL_RNG, input)

"""
    noise(::AbstractNoiseSource, input)

Return the RMS noise given the input signal
"""
noise(source::AbstractNoiseSource, input)

"""
    limimts(::AbstractNoiseSource)

Returns the limits for signal input. By Default, is `(0, Inf)`
"""
limits(::AbstractNoiseSource) = 0, Inf

"""
    ShotNoise()

A simple Poisson process for modifying incoming signal according to counting statistics.

# Examples

Simulate a Poisson process with average value of 1 three times.
```julia
julia> proc = ShotNoise()

julia> [proc(1) for _ in 1:3]
3-element Vector{Int64}:
 0
 3
 2

julia> noise(proc, 100)
10.0
```
"""
Base.@kwdef struct ShotNoise <: AbstractNoiseSource end

(::ShotNoise)(rng, input) = rand(rng, Poisson(input))

noise(::ShotNoise, input) = sqrt(input)

"""
    ReadNoise(;noise)

Gaussian read noise, with the given RMS noise.
"""
Base.@kwdef struct ReadNoise{T} <: AbstractNoiseSource
    noise::T
end

(tf::ReadNoise)(rng, input)= tf.noise * randn(rng) + input

noise(tf::ReadNoise, input) = tf.noise


text:
# This file includes methods for data generation.

export getdata, getpoint

"""
    getdata(f, vtx, npts)

Returns a vector of (d + m)-dimensional points that are dispersed randomly over d-dimensioanl domain with vertex points
`vtx`. `f` is the vector-valued function used to map each point in the domain. `npts` is the number of points to returned.

    getdata(vtx, npts)

Returns a vector of d-dimensional points that are dispersed randomly over d-dimensioanl domain with vertex points
`vtx`. `npts` is the number of points to be returned`
"""
function getdata end

getdata(f, vtx::AbstractVector{<:AbstractVector}, args::Int...) = getdata(f, ngon(vtx), args...)
getdata(vtx::AbstractVector{<:AbstractVector}, args::Int...) = getdata(ngon(vtx), args...)

getdata(domain::Line, nedgepoints::Int) = boundarypoints(domain, nedgepoints)
getdata(f, domain::Line, nedgepoints::Int) = [Point(pnt..., f(pnt...)...) for pnt in boundarypoints(domain, nedgepoints)]

# nedgepoints=10 is for backwards compatibility
getdata(domain::Ngon, ninterpoints::Int, nedgepoints::Int=10) = disperse(domain, ninterpoints, nedgepoints)
getdata(f, domain::Ngon, ninterpoints::Int, nedgepoints::Int=10) =
    [Point(pnt..., f(pnt...)...) for pnt in disperse(domain, ninterpoints, nedgepoints)]

ngon(pts::AbstractVector{<:AbstractVector}) = Ngon(SVector{length(pts)}(map(item -> Point(item...), pts)))

"""
    getpoint(ngon::Ngon, maxiter::Int=100_000)

Returns a valid that is inside of `ngon`.

    getpoint(line::Line)

Returns a valid point that is inside `line`.
"""
function getpoint end

function getpoint(ngon::Ngon{Dim, T, N, P}; maxiter::Int=100_000) where {Dim, T, N, P}
    tess = spt.Delaunay(coordinates(ngon))
    A, b = boundboxtransforms(ngon)
    iter = 1
    while iter ≤ maxiter
        val = A * rand(T, Dim) + b
        pnt = Point(val...)
        isvalidpoint(pnt, tess) && return pnt
        iter += 1
    end
    return Point(NaN, NaN)
end

function getpoint(line::Line{1, T}) where {T}
    p0, p1 = only(line[1]), only(line[2])
    pnt = p0 + (p1 - p0) * rand(T)
    Point(pnt)
end

isvalidpoint(pnt::AbstractPoint, tess) = tess.find_simplex(pnt)[1] ≥ 0 && pnt !== Point(NaN, NaN)

function disperse(ngon::Ngon, ninterpoints::Int, nedgepoints::Int)
    interpoints = interiorpoints(ngon, ninterpoints)
    boundpoints = boundarypoints(ngon, nedgepoints)
    vcat(boundpoints, interpoints) |> unique
end

function interiorpoints(ngon::Ngon, ninterpoints::Int)
    allpnts = [getpoint(ngon) for i in 1 : 10 * ninterpoints]
    [Point(val[1], val[2]) for val in eachcol(kmeans(hcat(collect.(allpnts)...), ninterpoints).centers)]
end

boundarypoints(line::Line, nedgepoints::Int=10) = edgepoints(line[1], line[2], nedgepoints)
boundarypoints(ngon::Ngon, nedgepoints::Int=10) = vcat(map(((pnt1, pnt2),) -> edgepoints(pnt1, pnt2, nedgepoints),
    TupleView{2, 1}(SVector([ngon.points; [ngon.points[1]]]...)))...)

function edgepoints(p1::AbstractPoint{2, T}, p2::AbstractPoint{2, T}, nedgepoints::Int=10) where {T}
    x = collect(LinRange(p1[1], p2[1], nedgepoints))
    y = collect(LinRange(p1[2], p2[2], nedgepoints))
    [Point(xi, yi) for (xi, yi) in zip(x, y)]
end
edgepoints(p1::AbstractPoint{1, T}, p2::AbstractPoint{1, T}, nedgepoints::Int=10) where {T} =
    [Point(xi) for xi in collect(LinRange(p1[1], p2[1], nedgepoints))]

function boundboxtransforms(ngon::Ngon)
    coords = coordinates(ngon)
    x, y = getindex.(coords, 1), getindex.(coords, 2)
    xmin, xmax, ymin, ymax = minimum(x), maximum(x), minimum(y), maximum(y)
    xwidth, ywidth = xmax - xmin, ymax - ymin
    A = [xwidth 0; 0 ywidth]
    b = [xmin, ymin]
    A, b
end

# TODO: Complete function
function filtertriangle(quality, args...; kwargs...) end


text:
immutable LowRankPertOperator{OO,LR,T} <: Operator{T}
    op::OO
    pert::LR

    function LowRankPertOperator(a::OO,b::LR)
        @assert domainspace(a)==domainspace(b)
        @assert rangespace(a)==rangespace(b)
        new(a,b)
    end
end

function LowRankPertOperator(Bin::Operator,Lin::LowRankOperator)
    B,L2=promotedomainspace([Bin,Lin])
    rsp=rangespace(B)  # use rangespace of B because LowRankOperator only
                        # needs convert, and its unlikely that the rangespaces
                        # will be inferred from L
    L=promoterangespace(L2,rsp)

    LowRankPertOperator{typeof(B),typeof(L),promote_type(eltype(L),eltype(B))}(B,L)
end

Base.convert{OT<:Operator}(::Type{Operator},V::Vector{OT})=LowRankPertOperator(V)

Base.getindex(L::LowRankPertOperator,k::Integer,j::Integer)=L.op[k,j]+L.pert[k,j]

domainspace(L::LowRankPertOperator)=domainspace(L.op)
rangespace(L::LowRankPertOperator)=rangespace(L.op)
datasize(L::LowRankPertOperator,k...)=datasize(L.pert,k...)

for OP in (:promotedomainspace,:promoterangespace)
    @eval $OP(L::LowRankPertOperator,sp::Space)=LowRankPertOperator($OP(L.op,sp),$OP(L.pert,sp))
end

## algebra

for OP in (:+,:-)
    @eval $OP(A::LowRankPertOperator,B::LowRankPertOperator) =
        LowRankPertOperator($OP(A.op,B.op),$OP(A.pert,B.pert))
end

*(L::LowRankPertOperator,f::Fun)=L.op*f+L.pert*f

*(L::LowRankPertOperator,B::LowRankOperator)=L.op*B+L.pert*B
*(B::LowRankOperator,L::LowRankPertOperator)=B*L.op+B*L.pert

*(A::LowRankPertOperator,B::LowRankPertOperator)=A.op*B + A.pert*B

# ambiguity
for TYP in (:TimesOperator,:ZeroOperator,:PlusOperator,:Conversion,:Operator)
    @eval begin
        +(L::LowRankOperator,B::$TYP) = LowRankPertOperator(B,L)
        +(B::$TYP,L::LowRankOperator) = LowRankPertOperator(B,L)

        -(L::LowRankOperator,B::$TYP)=LowRankPertOperator(-B,L)
        -(B::$TYP,L::LowRankOperator)=LowRankPertOperator(B,-L)

        *(L::LowRankPertOperator,B::$TYP)=LowRankPertOperator(L.op*B,L.pert*B)
        *(B::$TYP,L::LowRankPertOperator)=LowRankPertOperator(B*L.op,B*L.pert)
    end
end


text:
@testset "543.diameter-of-binary-tree.jl" begin
    t = convert(TreeNode{Int}, [1, 2, 3, 4, 5])
    @test diameter_of_binary_tree(t) == 3
end