text:
const term_c = (Tokens.RPAREN, Tokens.RSQUARE, Tokens.RBRACE, Tokens.END, Tokens.ELSE, Tokens.ELSEIF, Tokens.CATCH, Tokens.FINALLY, Tokens.ENDMARKER)

"""
Continue parsing statements until an element of `closers` is hit (usually
`end`). Statements are grouped in a `Block` EXPR.
"""
function parse_block(ps::ParseState, ret::Vector{EXPR} = EXPR[], closers = (Tokens.END,), docable = false)
    while kindof(ps.nt) ∉ closers # loop until an expected closer is hit
        if kindof(ps.nt) ∈ term_c # error handling if an unexpected closer is hit
            if kindof(ps.nt) == Tokens.ENDMARKER
                break
            elseif kindof(ps.nt) == Tokens.RPAREN
                if length(ps.closer.cc) > 1 && :paren == ps.closer.cc[end - 1]
                    break
                else
                    push!(ret, mErrorToken(INSTANCE(next(ps)), UnexpectedToken))
                    ps.errored = true
                end
            elseif kindof(ps.nt) == Tokens.RBRACE
                if length(ps.closer.cc) > 1 && :brace == ps.closer.cc[end - 1]
                    break
                else
                    push!(ret, mErrorToken(INSTANCE(next(ps)), UnexpectedToken))
                    ps.errored = true
                end
            elseif kindof(ps.nt) == Tokens.RSQUARE
                if length(ps.closer.cc) > 1 && :square == ps.closer.cc[end - 1]
                    break
                else
                    push!(ret, mErrorToken(INSTANCE(next(ps)), UnexpectedToken))
                    ps.errored = true
                end
            else
                push!(ret, mErrorToken(INSTANCE(next(ps)), UnexpectedToken))
                ps.errored = true
            end
        else
            if docable
                a = parse_doc(ps)
            else
                a = parse_expression(ps)
            end
            push!(ret, a)
        end
    end
    return ret
end

"""
Parses an iterator, allowing for the preceding keyword `outer`. Returns an
error expression if an invalid expression is parsed (anything other than
`=`, `in`, `∈`).
"""
function parse_iterator(ps::ParseState)
    startbyte = ps.nt.startbyte
    if kindof(ps.nt) == Tokens.OUTER && kindof(ps.nws) != EmptyWS && !Tokens.isoperator(kindof(ps.nnt))
        outer = INSTANCE(next(ps))
        arg = @closer ps :range @closer ps :ws parse_expression(ps)
        if is_range(arg)
            arg.args[1] = setparent!(EXPR(Outer, EXPR[outer, arg.args[1]]), arg)
            arg.fullspan += outer.fullspan
            arg.span = outer.fullspan + arg.span
        else # error handling - incorrect iterator specification
            arg = EXPR(ErrorToken, EXPR[outer, arg])
            arg.kind = InvalidIterator
        end
    else
        arg = @closer ps :range @closer ps :ws parse_expression(ps)
    end
    return arg
end

"""
Parses a group of iterators e.g. used in a `for` loop or generator. Can allow
for a succeeding `Filter` expression.
"""
function parse_iterators(ps::ParseState, allowfilter = false)
    startbyte = ps.nt.startbyte
    arg = parse_iterator(ps)
    if (typof(arg) === Outer && !is_range(arg.args[2])) || !is_range(arg)
        # error handling - incorrect iterator specification
        arg = mErrorToken(arg, InvalidIterator)
        ps.errored = true
    elseif kindof(ps.nt) == Tokens.COMMA # we've hit a comma separated list of iterators.
        arg = EXPR(Block, EXPR[arg])
        while kindof(ps.nt) == Tokens.COMMA
            accept_comma(ps, arg)
            nextarg = parse_iterator(ps)
            if (typof(nextarg) === Outer && !is_range(nextarg.args[2])) || !is_range(nextarg)
                # error handling - incorrect iterator specification
                arg = mErrorToken(arg, InvalidIterator)
                ps.errored = true
            end
            push!(arg, nextarg)
        end
    end

    if allowfilter && kindof(ps.nt) === Tokens.IF # assumes we're inside a generator
        if typof(arg) === Block
            arg = EXPR(Filter, arg.args)
        else
            arg = EXPR(Filter, EXPR[arg])
        end
        push!(arg, mKEYWORD(next(ps)))
        cond = @closer ps :range parse_expression(ps)
        push!(arg, cond)
    end
    return arg
end

function is_range(x::EXPR)
    typof(x) === BinaryOpCall && (is_eq(x.args[2]) || is_in(x.args[2]) || is_elof(x.args[2]))
end

"""
    parse_call(ps, ret)

Parses a function call. Expects to start before the opening parentheses and is passed the expression declaring the function name, `ret`.
"""
function parse_call(ps::ParseState, ret::EXPR, ismacro = false)
    sb = ps.nt.startbyte - ret.fullspan
    # special case for `new` call inside Type definitions, converts `new` to a keyword.
    if typof(ret) === IDENTIFIER && valof(ret) == "new" && :struct in ps.closer.cc
        ret = mKEYWORD(Tokens.NEW, ret.fullspan, ret.span)
    elseif typof(ret) === Curly && valof(ret.args[1]) == "new" && :struct in ps.closer.cc
        ret.args[1] = setparent!(mKEYWORD(Tokens.NEW, ret.args[1].fullspan, ret.args[1].span), ret)
    elseif is_minus(ret) || is_not(ret)
        arg = @closer ps :unary @closer ps :inwhere @precedence ps PowerOp parse_expression(ps)
        if typof(arg) === TupleH
            pushfirst!(arg.args, ret)
            fullspan = ps.nt.startbyte - sb
            ret = EXPR(Call, arg.args, fullspan, fullspan - (last(arg.args).fullspan - last(arg.args).span))
        elseif typof(arg) === WhereOpCall && typof(arg.args[1]) === TupleH
            ret = mWhereOpCall(EXPR(Call, EXPR[ret; arg.args[1].args]), arg.args[2], arg.args[3:end])
        else
            ret = mUnaryOpCall(ret, arg)
        end
    elseif is_and(ret) || is_decl(ret) || is_exor(ret)
        arg = @precedence ps 20 parse_expression(ps)
        if is_exor(ret) && typof(arg) === TupleH && length(arg.args) == 3 && typof(arg.args[2]) === UnaryOpCall && is_dddot(arg.args[2].arg[2])
            arg = EXPR(InvisBrackets, arg.args)
        end
        ret = mUnaryOpCall(ret, arg)
    elseif is_issubt(ret) || is_issupt(ret)
        arg = @precedence ps PowerOp parse_expression(ps)
        ret = EXPR(Call, EXPR[ret; arg.args])
    else
        !ismacro && typof(ret) === MacroName && (ismacro = true)
        args = EXPR[ret, mPUNCTUATION(next(ps))]
        @closeparen ps @default ps parse_comma_sep(ps, args, !ismacro)
        accept_rparen(ps, args)
        fullspan = ps.nt.startbyte - sb
        ret = EXPR(ismacro ? MacroCall : Call, args, fullspan, fullspan - last(args).fullspan + last(args).span)
    end
    return ret
end

_do_kw_convert(ps, a) = !ps.closer.brace && typof(a) === BinaryOpCall && is_eq(a.args[2])
_kw_convert(a) = EXPR(Kw, EXPR[a.args[1], a.args[2], a.args[3]], a.fullspan, a.span)

"""
Parses a comma separated list, optionally allowing for conversion of
assignment (`=`) expressions to `Kw`.
"""
function parse_comma_sep(ps::ParseState, args::Vector{EXPR}, kw = true, block = false, istuple = false)
    @nocloser ps :inwhere @nocloser ps :newline @closer ps :comma while !closer(ps)
        a = parse_expression(ps)
        if kw && _do_kw_convert(ps, a)
            a = _kw_convert(a)
        end
        push!(args, a)
        if kindof(ps.nt) == Tokens.COMMA
            accept_comma(ps, args)
        elseif kindof(ps.ws) == SemiColonWS
            break
        end
    end
    if istuple && length(args) > 2
        block = false
    end

    if kindof(ps.ws) == SemiColonWS
        if @nocloser ps :newline @closer ps :comma @nocloser ps :semicolon closer(ps)
            if block && !(length(args) == 1 && ispunctuation(args[1])) && !(typof(last(args)) === UnaryOpCall && is_dddot(last(args).args[2]))
                push!(args, EXPR(Block, EXPR[pop!(args)]))
            elseif kw && kindof(ps.nt) === Tokens.RPAREN
                push!(args, EXPR(Parameters, EXPR[], 0, 0))
            end
        else
            a = @nocloser ps :newline @closer ps :comma @nocloser ps :inwhere parse_expression(ps)
            if block && !(length(args) == 1 && ispunctuation(args[1])) && !(typof(last(args)) === UnaryOpCall && is_dddot(last(args).args[2])) && !(istuple && kindof(ps.nt) === Tokens.COMMA)
                args1 = EXPR[pop!(args), a]
                @nocloser ps :newline @closer ps :comma while @nocloser ps :semicolon !closer(ps)
                    a = parse_expression(ps)
                    push!(args1, a)
                end
                body = EXPR(Block, args1)
                push!(args, body)
                args = body
            else
                parse_parameters(ps, args, EXPR[a])
            end
        end
    end
    return
end

"""
Parses parameter arguments for a function call (e.g. following a semicolon).
"""
function parse_parameters(ps::ParseState, args::Vector{EXPR}, args1::Vector{EXPR} = EXPR[]; usekw = true)
    if isempty(args1)
        sb = ps.nt.startbyte
        isfirst = true
    else
        sb = ps.nt.startbyte - args1[1].fullspan
        isfirst = false
    end
    @nocloser ps :inwhere @nocloser ps :newline  @closer ps :comma while !isfirst || (@nocloser ps :semicolon !closer(ps))
        if isfirst
            a = parse_expression(ps)
        else
            a = first(args1)
        end
        if usekw && _do_kw_convert(ps, a)
            a = _kw_convert(a)
        end
        if isfirst
            push!(args1, a)
        else
            pop!(args1)
            push!(args1, a)
        end
        if kindof(ps.nt) == Tokens.COMMA
            accept_comma(ps, args1)
        end
        if kindof(ps.ws) == SemiColonWS
            parse_parameters(ps, args1; usekw = usekw)
        end
        isfirst = true
    end
    if !isempty(args1)
        fullspan = ps.nt.startbyte - sb
        paras = EXPR(Parameters, args1, fullspan, fullspan - last(args1).fullspan + last(args1).span)
        push!(args, paras)
    end
    return
end

"""
    parse_macrocall(ps)

Parses a macro call. Expects to start on the `@`.
"""
function parse_macrocall(ps::ParseState)
    sb = ps.t.startbyte
    at = mPUNCTUATION(ps)
    if !isemptyws(ps.ws)
        mname = mErrorToken(INSTANCE(next(ps)), UnexpectedWhiteSpace)
        ps.errored = true
    else
        mname = EXPR(MacroName, EXPR[at, mIDENTIFIER(next(ps))])
    end

    # Handle cases with @ at start of dotted expressions
    if kindof(ps.nt) == Tokens.DOT && isemptyws(ps.ws)
        while kindof(ps.nt) == Tokens.DOT
            op = mOPERATOR(next(ps))
            nextarg = mIDENTIFIER(next(ps))
            mname = mBinaryOpCall(mname, op, EXPR(Quotenode, EXPR[nextarg]))
        end
    end

    if kindof(ps.nt) == Tokens.COMMA
        return EXPR(MacroCall, EXPR[mname], mname.fullspan, mname.span)
    elseif isemptyws(ps.ws) && kindof(ps.nt) == Tokens.LPAREN
        return parse_call(ps, mname, true)
    else
        args = EXPR[mname]
        insquare = ps.closer.insquare
        @default ps while !closer(ps)
            if insquare
                a = @closer ps :insquare @closer ps :inmacro @closer ps :ws @closer ps :wsop parse_expression(ps)
            else
                a = @closer ps :inmacro @closer ps :ws @closer ps :wsop parse_expression(ps)
            end
            push!(args, a)
            if insquare && kindof(ps.nt) == Tokens.FOR
                break
            end
        end
        fullspan = ps.nt.startbyte - sb
        return EXPR(MacroCall, args, fullspan, fullspan - last(args).fullspan + last(args).span)
    end
end

"""
parse_generator(ps)

Having hit `for` not at the beginning of an expression return a generator.
Comprehensions are parsed as SQUAREs containing a generator.
"""
function parse_generator(ps::ParseState, ret::EXPR)
    kw = mKEYWORD(next(ps))
    ret = EXPR(Generator, EXPR[ret, kw])
    ranges = @closesquare ps parse_iterators(ps, true)

    if typof(ranges) === Block
        append!(ret, ranges)
    else
        push!(ret, ranges)
    end

    if typof(ret.args[1]) === Generator || typof(ret.args[1]) === Flatten
        ret = EXPR(Flatten, EXPR[ret])
    end

    return ret
end

"""
Helper function for parsing import/using statements.
"""
function parse_dot_mod(ps::ParseState, is_colon = false)
    args = EXPR[]

    while kindof(ps.nt) == Tokens.DOT || kindof(ps.nt) == Tokens.DDOT || kindof(ps.nt) == Tokens.DDDOT
        d = mOPERATOR(next(ps))
        trailing_ws = d.fullspan - d.span
        if is_dot(d)
            push!(args, mOPERATOR(1 + trailing_ws, 1, Tokens.DOT, false))
        elseif is_ddot(d)
            push!(args, mOPERATOR(1, 1, Tokens.DOT, false))
            push!(args, mOPERATOR(1 + trailing_ws, 1, Tokens.DOT, false))
        elseif is_dddot(d)
            push!(args, mOPERATOR(1, 1, Tokens.DOT, false))
            push!(args, mOPERATOR(1, 1, Tokens.DOT, false))
            push!(args, mOPERATOR(1 + trailing_ws, 1, Tokens.DOT, false))
        end
    end

    # import/export ..
    if kindof(ps.nt) == Tokens.COMMA || kindof(ps.ws) == NewLineWS || kindof(ps.nt) == Tokens.ENDMARKER
        if length(args) == 2
            return EXPR[INSTANCE(ps)]
        end
    end

    while true
        if kindof(ps.nt) == Tokens.AT_SIGN
            at = mPUNCTUATION(next(ps))
            a = INSTANCE(next(ps))
            push!(args, EXPR(MacroName, EXPR[at, a]))
        elseif kindof(ps.nt) == Tokens.LPAREN
            a = EXPR(InvisBrackets, EXPR[mPUNCTUATION(next(ps))])
            push!(a, @closeparen ps parse_expression(ps))
            accept_rparen(ps, a)
            push!(args, a)
        elseif kindof(ps.nt) == Tokens.EX_OR
            a = @closer ps :comma parse_expression(ps)
            push!(args, a)
        elseif !is_colon && isoperator(ps.nt)
            next(ps)
            push!(args, mOPERATOR(ps.nt.startbyte - ps.t.startbyte,  1 + ps.t.endbyte - ps.t.startbyte, kindof(ps.t), false))
        else
            push!(args, INSTANCE(next(ps)))
        end

        if kindof(ps.nt) == Tokens.DOT
            push!(args, mPUNCTUATION(next(ps)))
        elseif isoperator(ps.nt) && (ps.nt.dotop || kindof(ps.nt) == Tokens.DOT)
            push!(args, mPUNCTUATION(Tokens.DOT, 1, 1))
            ps.nt = RawToken(kindof(ps.nt), ps.nt.startpos, ps.nt.endpos, ps.nt.startbyte + 1, ps.nt.endbyte, ps.nt.token_error, false, ps.nt.suffix)
        else
            break
        end
    end
    args
end


text:
##### Beginning of file

# Automatically generated by:
# PredictMDExtra._print_list_of_package_imports()

##### End of file

ext: jl


text:
INITIAL_FLAG = 1 == 1

export fwhm, grid_resolution, grid_slice

function fwhm(data, b_mean, b_std)
    half_value = b_mean + 2 * b_std
    peak = findall(data .> half_value)
    if (length(peak) <= 5)
        return 0, 0
    end
    return (peak[end][1] - peak[1][1]), (peak[end][1] + peak[1][1]) / 2 #findmid(data)
end

function grid_resolution(data)
    # Input (mean cross section , mean background)
    # Output (mean band tickness, mean space tickness)
    diff_l = data[4:end] - data[1:end-3]
    b_mean = mean(abs.(diff_l))
    b_std = std(abs.(diff_l))
    peak = findall(abs.(diff_l) .> b_mean + b_std)
    flag = diff_l[peak[1]] > 0 ? 1 : 0
    prev = peak[1]
    space, stript, count = 0, 0, 0
    for idx in peak
        if flag == 0 && diff_l[idx] > 0
            space += idx - prev
            flag = 1
            prev = idx
        end
        if flag == 1 && diff_l[idx] < 0
            stript += idx - prev
            flag = 0
            prev = idx
            count += 1
        end
    end
    return stript / count, space / (count - 1)
end

function grid_slice(img)
    img = convert(Array{N0f16}, img)
    b_mean = mean(img)
    flag = Bool.(INITIAL_FLAG)
    width_temp = []
    pos_temp = []
    stript_temp = []
    space_temp = []
    frame = 1
    # loop through every frame
    for i = 5:size(img, 3)
        img1 = img[:, :, i]
        line = mean(img1, dims = 2)
        b_std = std(line)
        # skip frames without signal
        if (b_std > b_mean / 1)
            if (flag == INITIAL_FLAG)
                flag = !flag
                frame = i
            end
            continue
        end
        # take moving average in line signal
        new_line = zeros(size(line, 1), 1)
        for i = 5:size(line, 1)
            new_line[i] = mean(line[i-4:i])
        end
        # calculate image signal range
        (width, mid) = fwhm(new_line, b_mean, b_std)
        # analysis resolution
        #	img_l = mean(img1,dims=1)[1,:,1]
        #	if width != 0
        #		stript_width, space_width = grid_resolution(img_l)
        #	else
        #		stript_width, space_width = (0,0)
        #	end
        width != 0 && push!(width_temp, width)
        mid != 0 && push!(pos_temp, mid)
    end
    img_l = mean(mean(img, dims = 3), dims = 1)[1, :, 1]
    stript_width, space_width = grid_resolution(img_l)
    if isempty(width_temp)
        width_temp = 0
    end
    if isempty(pos_temp)
        pos_temp = 0
    end
    return width_temp, pos_temp, frame, stript_width, space_width
end


text:
function set_state!(vis::Visualizer, model::AbstractGameModel, sce::Scenario{T},
    x::SVx; key::Int=0) where {T,SVx}
    p = model.p
    for i = 1:p
        tr = transRot(model, x, i)
        player_path = "env/player_key$key/player$(sce.player[i].id)"
        settransform!(vis[player_path], tr)
    end
    return nothing
end

function Algames.set_traj!(vis::Visualizer, model::AbstractGameModel, sce::Scenario{T},
    traj; key::Int=0) where {T}
    N = length(traj)
    anim = MeshCat.Animation(10)
    for k = 1:N
        MeshCat.atframe(anim, k) do
            set_state!(vis, model, sce, Algames.state(traj[k]), key=key)
        end
    end
    MeshCat.setanimation!(vis, anim)
    return nothing
end

function set_waypoint!(vis::Visualizer, model::AbstractGameModel, sce::Scenario{T},
    x::SVx, k::Int; key::Int=0) where {T,SVx}
    p = model.p
    for i = 1:p
        tr = transRot(model, x, i)
        waypoint_path = "env/waypoint_key$key/player$(sce.player[i].id)/$k"
        settransform!(vis[waypoint_path], tr)
    end
    return nothing
end

function set_waypoint_traj!(vis::Visualizer, model::AbstractGameModel, sce::Scenario{T},
    traj; key::Int=0) where {T}
    N = length(traj)
    for k = 1:N
        set_waypoint!(vis, model, sce, Algames.state(traj[k]), k, key=key)
    end
    return nothing
end


text:
using ChartmetricScraper
using Test


text:
# ---
# title: 777. Swap Adjacent in LR String
# id: problem777
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Brainteaser
# link: <https://leetcode.com/problems/swap-adjacent-in-lr-string/description/>
# hidden: true
# ---
#
# In a string composed of `'L'`, `'R'`, and `'X'` characters, like
# `"RXXLRXRXL"`, a move consists of either replacing one occurrence of `"XL"`
# with `"LX"`, or replacing one occurrence of `"RX"` with `"XR"`. Given the
# starting string `start` and the ending string `end`, return `True` if and only
# if there exists a sequence of moves to transform one string to the other.
#
#
#
# **Example 1:**
#
#
#
#     Input: start = "RXXLRXRXL", end = "XRLXXRRLX"
#     Output: true
#     Explanation: We can transform start to end following these steps:
#     RXXLRXRXL ->
#     XRXLRXRXL ->
#     XRLXRXRXL ->
#     XRLXXRRXL ->
#     XRLXXRRLX
#
#
# **Example 2:**
#
#
#
#     Input: start = "X", end = "L"
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: start = "LLR", end = "RRL"
#     Output: false
#
#
# **Example 4:**
#
#
#
#     Input: start = "XL", end = "LX"
#     Output: true
#
#
# **Example 5:**
#
#
#
#     Input: start = "XLLR", end = "LXLX"
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `1 <= start.length <= 104`
#   * `start.length == end.length`
#   * Both `start` and `end` will only consist of characters in `'L'`, `'R'`, and `'X'`.
#
#
## @lc code=start
using LeetCode

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


text:
using HDF5

struct ConvexMesh
  problemSpec :: ProblemSpec
  points :: Array{Float64}
  simplIdx :: Array{UInt64}
end

# Loop through simplices
function Base.iterate(m :: ConvexMesh, state=1)
  if state > size(m)[1]
    return nothing
  end

  i = m.simplIdx[state,:]
  (m.points[i,:], state+1)
end

# Length of the mesh is the number of simplices
# Needed for many iteration utilities
Base.length(m :: ConvexMesh) = Base.size(m)[1]

# Size of the mesh is the number of simplices and dimensions
function Base.size(m :: ConvexMesh)
  (size(m.simplIdx)[1], size(m.points)[2])
end

function savehdf5(filename::String, m::ConvexMesh)
  dimnames = keys(m.problemSpec)
  dimranges = values(m.problemSpec)

  h5open(filename, "w") do file
    write(file, "dimnames", dimnames)
    write(file, "support", dimranges)
    write(file, "points", m.points)
    write(file, "simplices", m.simplIdx)
  end
end

function readhdf5(filename::String)
  h5open(filename, "r") do file
    dimnames = read(file, "dimnames")
    support = read(file, "support")
    points = read(file, "points")
    simplIdxs = read(file, "simplices")

    ps = ProblemSpec(x[1] => Tuple(x[2]) for x in zip(dimnames, support))
    ConvexMesh(ps, points, simplIdxs)
  end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    OptimizationLibrary

An optimization library (e.g. Optim.jl, Convex.jl, JuMP.jl).
"""
abstract type OptimizationLibrary end

# optimization libraries for dispatch
struct JuliaLib  <: OptimizationLibrary end
struct OptimLib  <: OptimizationLibrary end
struct ConvexLib <: OptimizationLibrary end
struct JuMPLib   <: OptimizationLibrary end


text:
using DataFrames

mutable struct Climate
    climate_data::DataFrame
    rainfall_id::String
    et_id::String
    # t_id::String
end

function subcatchment_data(node::NetworkNode, climate::Climate)
    node_id = node.node_id
    cols = filter(x -> occursin(node_id, string(x)), propertynames(climate))

    return cols
end

function climate_values(node::NetworkNode, climate::Climate, timestep::Union{Nothing, Int}=nothing)
    node_id = node.node_id

    data = climate.climate_data

    rain_col = filter(x -> occursin(node_id, string(x))
                        & occursin(climate.rainfall_id, string(x)),
                        propertynames(data))
    et_col = filter(x -> occursin(node_id, string(x))
                      & occursin(climate.et_id, string(x)),
                      propertynames(data))

    if isempty(rain_col) | isempty(et_col)
        return (missing, missing)
    end

    if isnothing(timestep)
        return select(data, vcat(rain_col, et_col))
    end

    return select(data, vcat(rain_col, et_col))[timestep, :]
end

function sim_length(climate::Climate)
    return nrow(climate.climate_data)
end

max_stars_repo_name: ConnectedSystems/camgraph

max_issues_repo_name: ConnectedSystems/camgraph

max_forks_repo_name: ConnectedSystems/camgraph


text:
# Code for time and frequency response of LTI systems
using LinearAlgebra
# include("Utilities.jl");

# Frequency response at user specified frequencies
function FrequencyResponse(sys::StateSpace,no::Integer,ni::Integer, Om::Array{Float64,1})::Vector{ComplexF64}
    CheckIODimensions(sys,no,ni);

    B = sys.B[:,ni];
    C = reshape(sys.C[no,:],1,length(sys.C[no,:])); # Must ensure row vector
    D = sys.D[no,ni];
    I = im*eye(size(sys.A)[1]);
    TF = Vector{ComplexF64}(undef,length(Om));
    for (i,om) in enumerate(Om)
        G = C*inv(om*I-sys.A)*B; # This is vector of one element
        TF[i] = G[1]+D; # that is why we need G[1], otherwise the addition is across different types.
    end
    return(TF);
end

# Frequency response with no user specified om. This makes an adhoc grid between [1E-3 100*wn.max], with 300 points.
function FrequencyResponse(sys::StateSpace,no::Integer,ni::Integer)::Tuple{Array{Complex{Float64},1},Array{Float64,1}}

    CheckIODimensions(sys,no,ni);

    l,ev = eigen(sys.A);
    lmax = maximum(abs.(l)); # Find highest natural frequency
    Om = logspace(1E-3,lmax*100,300); # Making adhoc decision about 100 times highest natural frequency.

    B = sys.B[:,ni];
    C = reshape(sys.C[no,:],1,length(sys.C[no,:])); # Must ensure row vector
    D = sys.D[no,ni];
    I = im*eye(size(sys.A)[1]);
    TF = Vector{ComplexF64}(undef,length(Om));
    for (i,om) in enumerate(Om)
        G = C*inv(om*I-sys.A)*B; # This is vector of one element ...
        TF[i] = G[1]+D; # ... that is why we need G[1], otherwise the addition is across different types.
    end
    return(TF,Om);
end

function StepResponse(sys::StateSpace,no::Integer,ni::Integer)::Tuple{Array{Float64,1}, Array{Float64,1}}
    l,ev = eigen(sys.A);
    lmin = minimum(abs.(l)); # Find the slowest mode
    lmax = maximum(abs.(l)); # Find the fastest mode
    tmax = 2.5*(2*pi)/lmin; # Quite adhoc -- 10 times the slowest time-constant.
    dt =  (2*pi)/lmax/50; # Quite adhoc -- Sampling is 50 times the fastest mode.
    t = Vector(0:dt:tmax);
    u = ones(length(t));
    y = ForcedResponse(sys,no,ni,t,u);
    return(t,y);
end

function StepResponse(sys::StateSpace,no::Integer,ni::Integer,t::Vector{Float64})::Vector{Float64}
    u = ones(length(t));
    y = ForcedResponse(sys,no,ni,t,u);
    return(y);
end

function ForcedResponse(sys::StateSpace,no::Integer,ni::Integer,t::Vector{Float64},u::Vector{Float64},x0::Vector{Float64}=Array{Float64}(undef,0))::Vector{Float64}
    CheckIODimensions(sys,no,ni);

    # Obtain a discrete-time system.
    dt = t[2]-t[1]; # Assumed to be uniform dt -- probably should put a check in here.
    ns = size(sys.A)[2];

    # We are only looking at SISO time response
    B = sys.B[:,ni];
    C = reshape(sys.C[no,:],1,length(sys.C[no,:])); # Must ensure row vector
    D = sys.D[no,ni];

    M = dt*[sys.A B; zeros(1,ns+1)]
    S = exp(M);
    F = S[1:ns,1:ns];
    G = S[1:ns,ns+1:ns+1];

    # Check initial condition
    if isempty(x0)
        x = zeros(ns);
    else
        x = x0;
    end

    # Time march
    y = similar(t);
    for i = 1:length(t)
        y[i] = (C*x)[1] + D*u[i];
        x = F*x + G*u[i];
    end
    return(y);
end

function ImpulseResponse(sys::StateSpace,no::Integer,ni::Integer)::Tuple{Array{Float64,1}, Array{Float64,1}}
    l,ev = eigen(sys.A);
    lmin = minimum(abs.(l)); # Find the slowest mode
    lmax = maximum(abs.(l)); # Find the fastest mode
    tmax = 2.5*(2*pi)/lmin; # Quite adhoc -- 2.5 times the slowest time-constant.
    dt =  (2*pi)/lmax/50; # Quite adhoc -- Sampling is 50 times the fastest mode.
    t = Vector(0:dt:tmax);
    y = ImpulseResponse(sys,no,ni,t);
    return(t,y);
end

function ImpulseResponse(sys::StateSpace,no::Integer,ni::Integer,t::Vector{Float64})::Vector{Float64}
    y = similar(t);
    B = sys.B[:,ni];
    C = reshape(sys.C[no,:],1,length(sys.C[no,:])); # Must ensure row vector
    for i=1:length(t)
        y[i] = (C*exp(sys.A*t[i])*B)[1];
    end
    return(y);
end

function InitialResponse(sys::StateSpace,x0::Vector{Float64})::Tuple{Vector{Float64}, Matrix{Float64}}
    l,ev = eigen(sys.A);
    lmin = minimum(abs.(l)); # Find the slowest mode
    lmax = maximum(abs.(l)); # Find the fastest mode
    tmax = 2.5*(2*pi)/lmin; # Quite adhoc -- 2.5 times the slowest time-constant.
    dt =  (2*pi)/lmax/50; # Quite adhoc -- Sampling is 50 times the fastest mode.
    t = Vector(0:dt:tmax);
    y = InitialResponse(sys,x0,t);
    return(t,y);
end

function InitialResponse(sys::StateSpace,x0::Vector{Float64},t::Vector{Float64})::Matrix{Float64}
    y = Matrix{Float64}(undef,length(t),length(x0));
    for i=1:length(t)
        y[i,:] = reshape(exp(sys.A*t[i])*x0,1,length(x0));
    end
    return(y);
end

max_stars_repo_name: isrlab/RobustControl

max_issues_repo_name: isrlab/RobustControl

max_forks_repo_name: isrlab/RobustControl


text:
using AcousticFWI,Seismic,PyPlot

function main()

    ns = 10
    ng = 250
    nz = 150
    nx = 250
    nt = 1024
    nf = 1024
    dz = 0.01
    dx = 0.01
    dt = 0.002
    f0 = 10.
    fmin = 1.
    fmax = 30.
    ext = 50
    atten_max = 2.5
    GNiter = 8
    CGiter = 20

    vp = 2.1*ones(nz,nx)
    vp[61:100,:] = 3.0
    vp[101:end,:] = 2.5
    for (i,iz) in enumerate(61:100)
        vp[iz,1:105+i] = 2.7
    end

    ot = zeros(ns)
    isz = 3*ones(Int,ns)
    isx = round(Int,collect(linspace(2,nx-1,ns)))
    igz = 3*ones(Int,ng)
    igx = round(Int,collect(linspace(1,nx,ng)))

    wav = Ricker(f0=f0,dt=dt)
    d = zeros(nt,ng,ns)
    for ishot = 1:ns
        u = Helmholtz(isz[ishot],isx[ishot],ot[ishot],vp,wav,fmin,fmax,nf,nt,dz,dx,dt,ext,atten_max)
        d[:,:,ishot] = u[:,3,igx]
    end

    vp0 = 2.3*ones(vp)
    vp_est,cost,f_inv = FWI(vp0,d,wav,isz,isx,igz,igx,ot,fmin,fmax,nf,dz,dx,dt,ext,atten_max,GNiter,CGiter)

    SeisPlot(d[:,:,round(Int,ns/2)],pclip=98,dx=1,dy=dt,xlabel="Trace",ylabel="t (s)",cmap="gray")
    figure(2) ; vmin = 1000*min(minimum(vp),minimum(vp0),minimum(vp_est))
              ; vmax = 1000*max(maximum(vp),maximum(vp0),maximum(vp_est))
              ; subplot(221) ; imshow(1000*vp,vmin=vmin,vmax=vmax,cmap="YlOrBr",extent=[0.,(nx-1)*1000*dx,(nz-1)*1000*dz,0.])
                             ; xlabel("x (m)")
                             ; ylabel("z (m)")
                             ; colorbar(aspect=15,ticks=[2100,2300,2500,2700,2900,3100])
                             ; ax = gca()
                             ; ax[:text](-600.,-100., "a)", fontsize=15, fontweight="bold")
                             ; ax[:set_xticks]([0,500,1000,1500,2000])
                             ; ax[:set_yticks]([0,300,600,900,1200,1500])
              ; subplot(222) ; imshow(1000*vp0,vmin=vmin,vmax=vmax,cmap="YlOrBr",extent=[0.,(nx-1)*1000*dx,(nz-1)*1000*dz,0.])
                             ; xlabel("x (m)")
                             ; ylabel("z (m)")
                             ; colorbar(aspect=15,ticks=[2100,2300,2500,2700,2900,3100])
                             ; ax = gca()
                             ; ax[:text](-600.,-100., "b)", fontsize=15, fontweight="bold")
                             ; ax[:set_xticks]([0,500,1000,1500,2000])
                             ; ax[:set_yticks]([0,300,600,900,1200,1500])
              ; subplot(223) ; imshow(1000*vp_est,vmin=vmin,vmax=vmax,cmap="YlOrBr",extent=[0.,(nx-1)*1000*dx,(nz-1)*1000*dz,0.])
                             ; xlabel("x (m)")
                             ; ylabel("z (m)")
                             ; colorbar(aspect=15,ticks=[2100,2300,2500,2700,2900,3100])
                             ; ax = gca()
                             ; ax[:text](-600.,-100., "c)", fontsize=15, fontweight="bold")
                             ; ax[:set_xticks]([0,500,1000,1500,2000])
                             ; ax[:set_yticks]([0,300,600,900,1200,1500])
              ; subplot(224) ; plot(1:length(cost[1]),cost[1],color="black",linewidth=2)
                             ; xlabel("Iteration")
                             ; ylabel("Normalized Misfit")
                             ; ax = gca()
                             ; ax[:text](-0.4,1.05, "d)", fontsize=15, fontweight="bold")
                             ; ax[:set_xticks]([1,2,3,4,5,6,7,8])
                             ; ax[:set_yticks]([0.2,0.4,0.6,0.8,1.0])

end

main()


text:
export
    # cudnnActivationMode_t
    CUDNN_ACTIVATION_SIGMOID,
    CUDNN_ACTIVATION_RELU,
    CUDNN_ACTIVATION_TANH,
    CUDNN_ACTIVATION_CLIPPED_RELU

type ActivationDesc
    ptr::Ptr{Void}
end

function ActivationDesc()
    p = Ptr{Void}[0]
    cudnnCreateActivationDescriptor(p)
    desc = ActivationDesc(p[1])
    finalizer(desc, cudnnDestroyActivationDescriptor)
    desc
end

Base.unsafe_convert(::Type{Ptr{Void}}, desc::ActivationDesc) = desc.ptr

function activation(mode, x; relu_nanopt=CUDNN_NOT_PROPAGATE_NAN, relu_ceiling=1.0)
    actdesc = ActivationDesc()
    cudnnSetActivationDescriptor(actdesc, mode, relu_nanopt, relu_ceiling)

    h = handle(x)
    T = eltype(x)
    xdesc = TensorDesc(x)
    y = similar(x)
    cudnnActivationForward(h, actdesc, T[1], xdesc, x, T[0], xdesc, y)

    function backward!(gy, gx)
        isvoid(gx) && return
        cudnnActivationBackward(h, actdesc, T[1], xdesc, y, xdesc, gy, xdesc, x, T[1], xdesc, gx)
    end
    y, backward!
end


text:
# This file was copied from Transducers.jl
# which is available under an MIT license (see LICENSE).
using PkgBenchmark
include("pprinthelper.jl")
group_target =
    PkgBenchmark.readresults(joinpath(@__DIR__, "result-target.json"))
group_baseline =
    PkgBenchmark.readresults(joinpath(@__DIR__, "result-baseline.json"))
judgement = judge(group_target, group_baseline)

displayresult(judgement)

printnewsection("Target result")
displayresult(group_target)

printnewsection("Baseline result")
displayresult(group_baseline)


text:
module Day4

using DelimitedFiles

# Pair of 5x5 markers
mutable struct BingoBoard
    numbers::Matrix{Int}
    markers::Matrix{Bool}
end

struct BingoGame
    boards::Vector{BingoBoard}
    numbers::Vector{Int}
end

parseint(s) = parse(Int, s, base=10)

function parse_input(input::String)
    lines = split(input, "\n")
    # First line is list of numbers called
    numbers = map(parseint, split(lines[1], ","))
    boards = []
    # Followed by groups of six lines
    for i in 2:6:(length(lines)-6)
        board = zeros(5, 5)
        for k in 1:5
            board[k, :] .= map(parseint, split(lines[i+k]))
        end
        push!(boards, board)
    end

    bingoboards = map(boards) do b
        BingoBoard(b, zeros(Bool, 5, 5))
    end
    return BingoGame(bingoboards, numbers)
end

function play(game::BingoGame)
    won_boards = zeros(Bool, length(game.boards))
    first_winner_score = nothing
    for call in game.numbers
        for (i, board) in enumerate(game.boards)
            board.markers += board.numbers .== call

            colsums = sum(board.markers, dims=1)
            rowsums = sum(board.markers, dims=2)
            if max(maximum(colsums), maximum(rowsums)) == 5 && !won_boards[i]
                score = sum(board.numbers[board.markers .== 0]) * call
                if sum(won_boards) == 0
                    first_winner_score = score
                end
                won_boards[i] = 1
                if sum(won_boards) == length(game.boards)
                    return (first_winner_score, score)
                end
            end
        end
    end
end

test = """7,4,9,5,11,17,23,2,0,14,21,24,10,16,13,6,15,25,12,22,18,20,8,19,3,26,1

22 13 17 11  0
 8  2 23  4 24
21  9 14 16  7
 6 10  3 18  5
 1 12 20 15 19

 3 15  0  2 22
 9 18 13 17  5
19  8  7 25 23
20 11 10 24  4
14 21 16 12  6

14 21 17 24  4
10 16 15  9 19
18  8 23 26 20
22 11 13  6  5
 2  0 12  3  7
"""

test_game = parse_input(test)
first, last = play(test_game)
@assert first == 4512
@assert last == 1924

input = read("resources/day4_input.txt", String)
game = parse_input(input)
play(game) |> display

end


text:
# Copyright 2014 Mitchell Kember. Subject to the MIT License.
# Google Code Jam Qualification Round 2014
# Problem A: Magic Trick

# Attempts to determine the volunteer's card given two card arrangements ca and
# cb, and the row that the card was on for each arrangement, a and b. Returns a
# vector of the cards that the volunteer could have chosen.
function magic(a, b, ca, cb)
	ra = ca[a,:]
	rb = cb[b,:]
	intersect(ra, rb)
end

# Given the possible cards, returns the appropriate message to print.
function message(cards)
	len = length(cards)
	if len == 1
		cards[1]
	elseif len == 0
		"Volunteer cheated!"
	else
		"Bad magician!"
	end
end

# Parses the volunteer's answer and the cards.
function parse()
	a = int(readline())
	cards = Array(Int, 4, 4)
	for i = 1:4
		cards[i,:] = int(split(readline()))
	end
	(a, cards)
end

function main()
	n = int(readline())
	for i = 1:n
		a, ca = parse()
		b, cb = parse()
		cards = magic(a, b, ca, cb)
		str = message(cards)
		println("Case #$i: ", str)
	end
end

main()


text:
@testset "Event store" begin
    store = Sialia.Events.get_event_store(;auto_snapshot=2)

    # Clean slate for testing
    Sialia.Events.reset!(store)

    @testset "Start" begin
        # Create connection, indexes, etc.
        Sialia.Events.run!(store)
        @test !isnothing(store.client)

        # Already started
        @test_throws Sialia.Events.EventException Sialia.Events.run!(store)
    end

    id1 = UUIDs.uuid4()
    id2 = UUIDs.uuid4()

    @testset "Saves events / snapshots" begin
        # Store an event
        event1 = Sialia.Events.Event(
            "testcontext",
            "testaggregate",
            "testevent";
            aggregate_id=id1,
            correlation_id=UUIDs.uuid4(),
            causation_id=UUIDs.uuid4(),
            data=Dict("testkey" => "testdata"),
            revision=Int(1),
            state=Dict("id" => id1, "revision" => Int(1)),
        )

        @test Sialia.Events.last_sequence(store, "events") == 0
        stored = Sialia.Events.save_events(store, [event1])
        @test Sialia.Events.last_sequence(store, "events") == 1
        @test length(stored) == 1
        @test stored[1].metadata.position == 1
        @test isnothing(Sialia.Events.get_snapshot(store, id1))

        event1dup = Sialia.Events.Event(
            "testcontext",
            "testaggregate",
            "testevent";
            aggregate_id=id1,
            correlation_id=UUIDs.uuid4(),
            causation_id=UUIDs.uuid4(),
            data=Dict("testkey" => "testdata"),
            revision=Int(1),
            state=Dict("id" => id1, "revision" => Int(1)),
        )

        @test_throws Sialia.Events.EventException Sialia.Events.save_events(store, [event1dup])

        # Sequence does not update when save errors
        @test Sialia.Events.last_sequence(store, "events") == 1

        # Store another event
        event2 = Sialia.Events.Event(
            "testcontext",
            "testaggregate",
            "testevent";
            aggregate_id=id1,
            correlation_id=UUIDs.uuid4(),
            causation_id=UUIDs.uuid4(),
            data=Dict("testkey" => "testdata2"),
            revision=Int(2),
            state=Dict("id" => id1, "revision" => Int(2)),
        )

        stored = Sialia.Events.save_events(store, [event2])
        @test Sialia.Events.last_sequence(store, "events") == 2
        @test length(stored) == 1
        @test stored[1].metadata.position == 2

        # Should create a snapshot
        @test Sialia.Events.get_snapshot(store, id1)["state"] == event2.state

        # Make some more events
        events = Vector{Sialia.Events.Event}()
        for i in 3:6
            # Store another event
            event = Sialia.Events.Event(
                "testcontext",
                "testaggregate",
                "testevent";
                aggregate_id=id1,
                correlation_id=UUIDs.uuid4(),
                causation_id=UUIDs.uuid4(),
                data=Dict("testkey" => "testdata$i"),
                revision=Int(i),
                state=Dict("id" => id1, "revision" => Int(i)),
            )

            push!(events, event)
        end

        stored = Sialia.Events.save_events(store, events)

        @test length(stored) == 4
        @test stored[4].metadata.position == 6

        @test Sialia.Events.get_snapshot(store, id1)["state"]["revision"] == 6

        # Store an event for a different aggregate
        event3 = Sialia.Events.Event(
            "testcontext",
            "testaggregate",
            "testevent";
            aggregate_id=id2,
            correlation_id=UUIDs.uuid4(),
            causation_id=UUIDs.uuid4(),
            data=Dict("testkey" => "testdata"),
            revision=Int(1),
            state=Dict("id" => id2, "revision" => Int(1)),
        )

        stored = Sialia.Events.save_events(store, [event3])
        @test length(stored) == 1
        @test stored[1].metadata.position == 7
        @test isnothing(Sialia.Events.get_snapshot(store, id2))
    end

    @testset "Aggregate events" begin
        iter = Sialia.Events.aggregate_events_iterator(store, id1)
        last_revision = 0
        for event in iter
            @test event.aggregate.revision - 1 == last_revision
            @test event.aggregate.id == id1
            last_revision = event.aggregate.revision
        end

        @test last_revision == 6

        iter = Sialia.Events.aggregate_events_iterator(store, id1, from_revision=2, to_revision=5)
        last_revision = 1
        for event in iter
            @test event.aggregate.revision - 1 == last_revision
            @test event.aggregate.id == id1
            last_revision = event.aggregate.revision
        end

        @test last_revision == 5
    end

    @testset "Unpublished events / Mark as published" begin
        iter = Sialia.Events.unpublished_events_iterator(store)
        unpublished = Base.collect(iter)

        @test length(unpublished) == 7

        Sialia.Events.mark_published_events(store, unpublished)

        iter = Sialia.Events.unpublished_events_iterator(store)
        unpublished = Base.collect(iter)

        @test length(unpublished) == 0
    end

    @testset "Last event" begin
        event = Sialia.Events.aggregate_last_event(store, id1)
        @test event.aggregate.revision == 6
        @test event.metadata.position == 6

        event = Sialia.Events.aggregate_last_event(store, id2)
        @test event.aggregate.revision == 1
        @test event.metadata.position == 7

        event = Sialia.Events.aggregate_last_event(store, UUIDs.uuid4())
        @test isnothing(event)
    end

    @testset "Replay" begin
        iter = Sialia.Events.replay_events(store)
        last_position = 0
        for event in iter
            @test event.metadata.position - 1 == last_position
            last_position = event.metadata.position
        end

        @test last_position == 7

        iter = Sialia.Events.replay_events(store, from=2, to=5)
        last_position = 1
        for event in iter
            @test event.metadata.position - 1 == last_position
            last_position = event.metadata.position
        end

        @test last_position == 5
    end

    @testset "Stop" begin
        Sialia.Events.stop!(store)
        @test isnothing(store.client)

        # Already stopped
        @test_throws Sialia.Events.EventException Sialia.Events.stop!(store)
    end

    @testset "Restart" begin
        # Can restart from existing data
        Sialia.Events.run!(store)
        Sialia.Events.stop!(store)
    end
end


text:
## Compatability for julia0.3
## ==========================

immutable Logic
  val::Int8
end

 const qf_uf = Logic(0)         # Uninterpreted Functions
 const qf_nra = Logic(1)        # Non-Linear Real Arithmetic
 const qf_nra_ode = Logic(2)    # Non-Linear Real Arithmetic
 const qf_bv = Logic(3)         # BitVectors
 const qf_rdl = Logic(4)        # Real difference logics
 const qf_idl = Logic(5)        # Integer difference logics
 const qf_lra = Logic(6)        # Real linear arithmetic
 const qf_lia = Logic(7)        # Integer linear arithmetic
 const qf_ufidl = Logic(8)      # UF + IDL
 const qf_uflra = Logic(9)      # UF + LRA
 const qf_bool = Logic(10)      # Only booleans
 const qf_ct = Logic(11)        # Costqf


text:
function _gaussian(E::T, p::T, E0::T, p0::T; dE = 5.0, dthetap = 0.01) where T<:Number
    thetap = acos(p)/(2pi)
    thetap0 = acos(p0)/(2pi)
    p = clamp(p,nextfloat(-1.0),prevfloat(1.0))
    a = ((E - E0)/dE)^2 + ((thetap - thetap0)/dthetap)^2
    J = inv(sqrt(1.0 - p^2))/2pi
    return J*exp(-0.5*a)/(2pi*dE*dthetap)
end

function gaussian(energy::AbstractVector, pitch::AbstractVector, E0, p0; kwargs...)
    nenergy = length(energy)
    npitch = length(pitch)

    f = zeros(nenergy,npitch)
    for j=1:npitch, i=1:nenergy
        f[i,j] = _gaussian(energy[i],pitch[j], E0, p0; kwargs...)
    end
    return f
end

function lnΔ_ee(ne, Te)
    # Thermal electron-electron collisions
    # NRL Plasma Formulary ne in cm^-3, Te in eV
    return 23.5 - log(sqrt(ne)*(Te^(-5/4))) - sqrt(1e-5 + ((log(Te) - 2)^2)/16)
end

function lnΔ_ei(ne, Te, ni, Ti; mu = H2_amu, Z = 1)
    # Electron-ion collisions
    # NRL Plasma Formulary ne/i in cm^-3, Te/i in eV, mu in amu
    mr = e_amu/mu
    if Ti*mr < Te < 10Z^2
        return 23 - log(sqrt(ne)*Z*Te^(-3/2))
    elseif Ti*mr < 10Z^2 < Te
        return 24 - log(sqrt(ne)/Te)
    elseif Te < Ti*mr
        return 16 - log(mu*sqrt(ni)*(Ti^(-3/2))*Z^2)
    end
    return 16.0
end

function lnΔ_ii(ni1, Ti1, ni2, Ti2; mu1=H2_amu, mu2=H2_amu, Z1=1,Z2=1)
    # Mixed ion-ion collisions
    # NRL Plasma Formulary ni in cm^-3 Ti in eV, mu in amu

    return 23 - log(((Z1*Z2*(mu1+mu2))/(mu1*Ti2 + mu2*Ti1))*sqrt((ni1*Z1^2)/Ti1 + (ni2*Z2^2)/Ti2))
end

function lnΔ_ii(ni1, Ti1, ni2, Ti2, ne, Te, beta_D; mu1=H2_amu, mu2=H2_amu, Z1=1, Z2=1)
    # Counter-streaming ions (relative velocity v_D = beta_D*c) in the presence of warm electrons
    # NRL Plasma Formulary ne in cm^-3, Te in eV

    m_i1 = mass_u*mu1
    m_i2 = mass_u*mu2
    m_e = mass_u*e_amu
    L1 = Ti1*e0/m_i1
    L2 = Ti2*e0/m_i2
    U = Te*e0/m_e

    if max(L1, L2) < (beta_D*c0)^2 < U
        return 43 - log((Z1*Z2*(mu1 + mu2)/(mu1*mu2*beta_D^2))*sqrt(ne/Te))
    else
        return lnΔ_ii(ni1, Ti1, ni2, Ti2; mu1=mu1, mu2=mu2, Z1=Z1, Z2=Z2)
    end
end

function legendre(x::T, nterms) where T

    a = one(T)
    b = x

    p = zeros(T, nterms)

    for l=0:min(1,nterms-1)
        p[l+1] = l == 0 ? a : b
    end

    for l in 2:(nterms-1)
        a, p[l+1] = p[l], ((2l-1)*x*p[l] - (l-1)*a)/l
    end

    return p
end

function heaviside(t::T) where T
    t < 0 ? zero(T) : one(T)
end

function _slowing_down_legendre_expansion(u::T, P::Vector, P0::Vector) where T
    L = 0:(length(P)-1)
    S = 0.5*sum((2*l + 1) * P[l+1] * P0[l+1] * u^(l*(l+1)) for l in L)
    return S
end

function _slowing_down_legendre_expansion(xi::T, xi0::T, u::T, nterms) where T
    P = legendre(xi, nterms)
    P0 = legendre(xi0, nterms)
    _slowing_down_legendre_expansion(u, P, P0)
end

function slowing_down_time(ne, Te, Ti, Zeff;
                           Ai = H2_amu, Ab = H2_amu, Aimp=C6_amu,
                           Zi = 1, Zimp = 6, Zb=1)
    # Calculate Slowing down time on electrons
    # Heating of toroidal plasmas by neutral injection: T H Stix 1972 Plasma Phys. 14 367
    # Te in keV, ne in cm^-3

    nimp = Zeff > 1 ? ne*(Zeff - 1)/(Zimp*(Zimp-1)) : zero(ne)
    ni = max(ne - Zimp*nimp,zero(ne))

    lnΔ_e = lnΔ_ei(ne, Te*1e3, ni, Ti*1e3; mu=Ab, Z=Zb)

    tau_s = 6.27e8 * Ab*((Te*1e3)^1.5)/(ne*lnΔ_e*Zb^2)

    return tau_s
end

function _electron_ion_drag_difference(Eb, ne, Te, Ti, Zeff; Ab=H2_amu, Ai=H2_amu, Aimp=C6_amu, Zb=1,Zi=1,Zimp=6)
    # Calculates difference of electron drag - ion drag
    # Energetic ion distribution resulting from neutral beam injection in tokamaks: J Gaffey 1976 J. Plasma Phys. 16 149
    # E0/Te in keV, ne in cm^-3

    m_e = e_amu*mass_u
    m_b = Ab*mass_u
    m_imp = Aimp*mass_u
    m_i = Ai*mass_u

    v_b = sqrt(2*Eb*e0*1e3/m_b)
    v_e = sqrt(2*Te*e0*1e3/m_e)

    nimp = Zeff > 1 ? ne*(Zeff-1)/((Zimp*(Zimp-1))) : zero(ne)
    ni = max(ne - Zimp*nimp,zero(ne))

    lnΔ_be = lnΔ_ei(ne, Te*1e3, ni, Ti*1e3; mu=Ab, Z=Zb)
    Γ_be = (2*pi*ne*(e0^4)*Zb^2*lnΔ_be)/(m_b^2)

    electron_drag = ((8*Γ_be*m_b)/(3*sqrt(pi)*m_e*v_e^3))*v_b^3

    lnΔ_bi = lnΔ_ii(ni, Ti*1e3, ni, Ti*1e3, ne, Te*1e3, v_b/c0; mu1=Ab, Z1=Zb, mu2=Ai, Z2=Zi)
    Γ_bi = (2*pi*ni*(e0^4)*Zi^2*Zb^2*lnΔ_bi)/(m_b^2)

    lnΔ_bimp = lnΔ_ii(ni, Ti*1e3, nimp, Ti*1e3, ne, Te*1e3, v_b/c0; mu1=Ab, Z1=Zb, mu2=Aimp, Z2=Zimp)
    Γ_bimp = (2*pi*nimp*(e0^4)*(Zimp^2)*(Zb^2)*lnΔ_bimp)/(m_b^2)

    ion_drag = 2*m_b*(Γ_bi/m_i + Γ_bimp/m_imp)

    return electron_drag - ion_drag
end

function critical_energy(ne, Te, Ti, Zeff; Emax = 300.0,
                         Ai = H2_amu, Aimp = C6_amu, Ab=H2_amu,
                         Zi = 1, Zimp=6, Zb=1)
    # Calculates critical energy using root finding on drags
    # Energetic ion distribution resulting from neutral beam injection in tokamaks: J Gaffey 1976 J. Plasma Phys. 16 149
    # E0/Te in keV, ne in cm^-3

    Ec = find_zero(x -> _electron_ion_drag_difference(x, ne, Te, Ti, Zeff;
                   Ai=Ai, Aimp=Aimp, Ab=Ab, Zi=Zi,Zimp=Zimp, Zb=Zb),
                  (0, Emax))
    return Ec
end

function approx_critical_energy(ne, Te, Zeff; correction_factor = 1.0,
                                Ai = H2_amu, Aimp = C6_amu, Ab=H2_amu,
                                Zi = 1, Zimp=6, Zb=1)
    # Calculate critical energy assuming lnΔ_be == lnΔ_bi
    # For DIII-D a correction factor of (lnΔ_bi/lnΔ_be)^(2/3) ≈ 1.2 can be used
    # Heating of toroidal plasmas by neutral injection: T H Stix 1972 Plasma Phys. 14 367
    # E0/Te in keV, ne in cm^-3

    nimp = Zeff > 1 ? ne*(Zeff - 1)/(Zimp*(Zimp-1)) : zero(ne)
    ni = max(ne - Zimp*nimp,zero(ne))

    # average charge to mass ratio (avg_cmr) is where approximation comes in
    avg_cmr = (ni*(Zi^2/Ai) + nimp*(Zimp^2/Aimp))/ne
    Ec = 14.8 * Ab * Te * correction_factor*avg_cmr^(2/3)

    return Ec
end

function _thermalization_time(v_b, v_c, tau_s)
    vb3 = v_b^3
    vc3 = v_c^3
    tau_s*log((vb3 + vc3)/vc3)/3
end

function thermalization_time(E0; ne = 1.0, Zeff = 1.0, Te = 1.0, Ti = Te,
                             Ai = H2_amu, Ab = Ai, Aimp = C6_amu,
                             Zi=1, Zimp=6, Zb=Zi, kwargs...) where T<:Number

    m_b = mass_u*Ab
    v_b = sqrt(2*(e0*E0*1e3)/(m_b))
    tau_s = slowing_down_time(ne, Te, Ti, Zeff; Ai = Ai, Aimp=Aimp, Ab = Ab, Zi=Zi, Zimp=Zimp, Zb=Zb)
    Ec = critical_energy(ne, Te, Ti, Zeff; Ai = Ai, Aimp=Aimp, Ab = Ab, Zi=Zi, Zimp=Zimp, Zb=Zb)
    v_c = sqrt(2*(e0*Ec*1e3)/m_b)

    return _thermalization_time(v_b, v_c, tau_s)
end

function _slowing_down(v::T, v_b::T, v_c::T, tau_s::T, Zeff::T, P::Vector, P0::Vector;
                       tau_on = 0.0, tau_off = 1.0, tau = 1.0) where T<:Number

    v >= v_b && return zero(T)

    v3 = v^3
    vb3 = v_b^3
    vc3 = v_c^3

    inv_v3vc3 = inv(v3 + vc3)
    u1 = v3/vb3
    u2 = (vb3 + vc3)*inv_v3vc3
    u = (u1*u2)^(Zeff/6)

    if tau_on > tau_off
        tau_off = Inf
    end
    t_b = _thermalization_time(v_b, v_c, tau_s)
    t_0 = tau_on*t_b
    t_th = tau_s*log(u2)/3 # threshold time
    t_1 = tau_off*t_b - t_0
    t = tau*t_b - t_0

    S = min(t,t_b,t_1)/t_b
    U = heaviside(t - t_th) - heaviside(t - t_1 - t_th)
    if U == 0
        return zero(T)
    end
    lex = _slowing_down_legendre_expansion(u, P, P0)
    g = max(S*tau_s*inv_v3vc3*lex,zero(T))/t_b
end

function _slowing_down(E::T, p::T, E0::T, p0::T, nterms; kwargs...) where T<:Number
    P = legendre(p, nterms)
    P0 = legendre(p0, nterms)
    _slowing_down(E, E0, P, P0; kwargs...)
end

function slowing_down(energy::AbstractVector, pitch::AbstractVector, E0, p0; nterms=50,
                      ne = 1.0, Zeff = 1.0, Te = 1.0, Ti = Te,
                      Ai = H2_amu, Ab = Ai, Aimp = C6_amu,
                      Zi=1, Zimp=6, Zb=Zi, kwargs...) where T<:Number

    nenergy = length(energy)
    npitch = length(pitch)

    # Calculate quantities that do not depend on energy/pitch
    m_b = mass_u*Ab
    v_b = sqrt(2*(e0*E0*1e3)/(m_b))
    tau_s = slowing_down_time(ne, Te, Ti, Zeff; Ai = Ai, Aimp=Aimp, Ab = Ab, Zi=Zi, Zimp=Zimp, Zb=Zb)
    Ec = critical_energy(ne, Te, Ti, Zeff; Ai = Ai, Aimp=Aimp, Ab = Ab, Zi=Zi, Zimp=Zimp, Zb=Zb)
    v_c = sqrt(2*(e0*Ec*1e3)/m_b)

    # Legendre expansion at birth pitch
    P0 = legendre(p0, nterms)

    f_slow = zeros(nenergy,npitch)
    for j=1:npitch
        P = legendre(pitch[j], nterms)
        for i=1:nenergy
            v = sqrt(2*(e0*energy[i]*1e3)/(m_b))
            J = v*((e0*1e3)/m_b)
            f_slow[i,j] = J*_slowing_down(v, v_b, v_c, tau_s, Zeff, P, P0; kwargs...)
        end
    end

    return f_slow
end

function _approx_slowing_down(v, v_b, v_c, p::T, p0::T, Zeff::T, tau_s::T, mi_mb::T;
                             tau = 1.0, tau_off=1.0, tau_on = 0.0) where T<:Number

    v >= v_b && return 0.0

    v3 = v^3
    vb3 = v_b^3
    vc3 = v_c^3

    inv_v3vc3 = inv(v3 + vc3)
    u1 = v3/vb3
    u2 = (vb3 + vc3)*inv_v3vc3
    u = (u1*u2)^(Zeff/6)

    if tau_on > tau_off
        tau_off = Inf
    end
    t_b = _thermalization_time(v_b, v_c, tau_s)
    t_0 = tau_on*t_b
    t_th = tau_s*log(u2)/3 # threshold time
    t_1 = tau_off*t_b - t_0
    t = tau*t_b - t_0

    S = min(t,t_b,t_1)/t_b
    U = heaviside(t - t_th) - heaviside(t - t_1 - t_th)
    if U == 0
        return zero(T)
    end

    beta = mi_mb*Zeff/2
    alpha = (beta/3)*(1 - p0^2)*log((1 + vc3/v3)/(1+vc3/vb3))
    g = tau_s*inv_v3vc3*inv(sqrt(4*pi*alpha))*exp(-((p - p0)^2)/(4*alpha))
    return g*S*U
end

function approx_slowing_down(energy::AbstractVector, pitch::AbstractVector, E0, p0;
                      ne = 1.0, Zeff = 1.0, Te = 1.0, Ti = Te,
                      Ai = H2_amu, Ab = Ai, Aimp = C6_amu,
                      Zi=1, Zimp=6, Zb=Zi, kwargs...) where T<:Number

    nenergy = length(energy)
    npitch = length(pitch)

    # Calculate quantities that do not depend on energy/pitch
    m_b = mass_u*Ab
    v_b = sqrt(2*(e0*E0*1e3)/(m_b))
    tau_s = slowing_down_time(ne, Te, Ti, Zeff; Ai = Ai, Aimp=Aimp, Ab = Ab, Zi=Zi, Zimp=Zimp, Zb=Zb)
    Ec = critical_energy(ne, Te, Ti, Zeff; Ai = Ai, Aimp=Aimp, Ab = Ab, Zi=Zi, Zimp=Zimp, Zb=Zb)
    v_c = sqrt(2*(e0*Ec*1e3)/m_b)

    nenergy = length(energy)
    npitch = length(pitch)

    f_slow = zeros(nenergy,npitch)
    for j=1:npitch, i=1:nenergy
        v = sqrt(2*(e0*energy[i]*1e3)/(m_b))
        J = v*((e0*1e3)/m_b)
        f_slow[i,j] = J*_approx_slowing_down(v, v_b, v_c, pitch[j], p0, Zeff, tau_s, Ai/Ab; kwargs...)
    end

    return f_slow
end

function _bimaxwellian(energy::S, pitch::S, T_perp::S, T_para::S;
                      vd_para=0.0, sigma=-1, Ai = H2_amu) where S<:Number

    m_i = mass_u*Ai

    T_perp = T_perp*1e3*e0
    T_para = T_para*1e3*e0
    Ed_para = (0.5*m_i*vd_para^2)

    E = energy*e0*1e3
    xi = clamp(sigma*pitch,nextfloat(-1.0),prevfloat(1.0))
    J = e0*1e3
    A = sqrt(E/(pi*T_para*T_perp^2))
    T1 = -(E*xi^2 + Ed_para - vd_para*xi*sqrt(2*m_i*E))/T_para
    T2 = -(1-xi^2)*E/T_perp
    return A*exp(T1 + T2)*J
end

function bimaxwellian(energy::AbstractVector, pitch::AbstractVector, T_perp, T_para; kwargs...)
    nenergy = length(energy)
    npitch = length(pitch)
    f_ep = zeros(nenergy, npitch)
    for j=1:npitch, i=1:nenergy
        f_ep[i,j] = _bimaxwellian(energy[i], pitch[j], T_perp, T_para; kwargs...)
    end
    T = (T_para + 2*T_perp)/3
    return f_ep
end

function maxwellian(energy, pitch, T; kwargs...)
    bimaxwellian(energy, pitch, T, T; kwargs...)
end


text:
function shifted_fft(x)
    return fftshift(fft(fftshift(x)))
end

function shifted_ifft(x)
    return fftshift(ifft(fftshift(x)))
end

function trapz(x, y; dim=1)
    perm = [dim:max(ndims(y),dim); 1:dim-1];
    y = permutedims(y, perm);
    if ndims(y) == 1
        m = 1;
    else
        m = size(y,1);
    end

    if m == 1
        M = length(y);
        out = sum(diff(x).*(y[1:M-1] + y[2:M])/2.0);
    else
        out = transpose(diff(x)) * (reshape(y,Val(2))[1:m-1,:] + reshape(y,Val(2))[2:m,:])/2.0;
        siz = size(y); siz = collect(siz); siz[1] = 1;
        out = reshape(out, tuple(siz...));
        out = permutedims(out, invperm(perm));
        ind = findall(collect(size(out)).==1);
        out = dropdims(out,dims=ind[1]);
        if length(out) == 1;
            out = out[1];
        end
    end

    return out
end

function notfinite(x)
    return ~isfinite(x)
end

function f0(q)
    #=
    same as
    W[~sp.isfinite(W)] = 0.0
    example
    WE = f0(W)
    =#
    # for (i,j) in enumerate(q)
    #     if isfinite(j) == false
    #         q[i]::Float64 = 0
    #     end
    # end
    # q
    if ~isfinite(q)
        return 0
    else
        return q
    end
end

function f1(q)
    #=
    same as
    Na[(~sp.isfinite(Na)) | (Na < 0.0)] = 1.0

    example:
    Ne = [f1(i,0) for i in Na]
    =#
    if isfinite(q) == false || q < 0
        q = 1.0
    end
    q
end

function f10(q)
    #=
    same as
    Na[(~sp.isfinite(Na)) | (Na < 0.0)] = 1.0

    example:
    Ne = [f1(i,0) for i in Na]
    =#
    if isfinite(q) == false || q < 0
        q = 0.0
    end
    q
end

function f2(para1, para2)
    #=
    same as RT[τ<0] = 0
    example:

    raman = sdo(τ, 62.5e-15, 77e-15)
    f2(raman, τ)
    plot(τ, raman)

    =#
    for i in 1:length(para1)
        if para2[i] < 0
            para1[i] = 0
        end
    end
    para1
end

function f3(x)
    return 2pi*3e8/x
end

function isPositive(para)
    if para < 0
        return 0
    end
end

function getElectricField(pulse, fiber)
    ω0 = 2pi*c/pulse.λ0
    Et = pulse.data .* exp.(-1im.*ω0 .* pulse.τ)
    Et .= real(Et)
    Et = sqrt(2/(ϵ0*c*effectiveArea(fiber.diameter))).*Et
    return Et
end

function flip(arr)
    arr1 = zeros(Complex{Float64}, (size(arr)[2], size(arr)[1]))
    for i in 1:size(arr)[2]
        arr1[i,:] = arr[:,i]
    end
    arr1
end

function bound(tmin, tmax, x)
    if x>=tmin && x<=tmax return true
    else return false
    end
end

function spectrogram(τ, U; window=1e-15)
    g(t) = exp.(-1.0*((τ.-t)./(2*window)).^2)
    points = floor(Int, (τ[end] - τ[1] )/window)+1
    out = zeros(Complex{Float64}, (length(τ), points*500))
    counter = minimum(τ)
    for i in 1:points*500
        out[:,i] = shifted_ifft((U .* g(counter)))
        counter += window/500
    end
    #out1 = flip(out)
    abs2.(out)
end

function sorder(γ, P0, τ0, β2)
    return sqrt(γ*P0*τ0/abs(β2))
end

function dlength(τ0, β2)
    return τ0^2/abs(β2)
end

function nlength(γ, P0)
    return 1/(γ*P0)
end

function getElectricFieldE(pulse, fiber)
    ω0 = 2pi*c/pulse.λ0
    Et = pulse.data .* cos.(ω0 .* pulse.τ)
    Et = sqrt(2/(ϵ0*c*effectiveArea(fiber.diameter))).*Et
    return Et
end

function get_env(ef)
    ew = fft(ef)
    ew[Int(size(ew)[1]/2+1):end] .= 0.0
    ep = 2.0.*ifft(ew)
    ep
end

function mfft(x)
    y = fft(x)
    if isreal(x[1])
        return y[1:Int(length(y)/2 + 1)]
    else
        return y
    end
end

function mifft(x)
    if length(x) % 2 == 0
        x = cat(x,conj(x[end-1:-1:2]),dims=1)
        return real.(ifft(x))
    else
        return ifft(x)
    end
end

function smoothstep(left, right, loss, x)
    x = clamp((x-left)/(right-left), 0, 1)
    return x*x*(3-2*x)*loss
end

function flipsmoothstep(left, right, loss,x)
    if x<left
        return 1.0
    elseif x>right
        return loss
    else
        x = clamp1(1-(x-left)/(right-left), 0, 1.0, loss)
        y = x*x*(3.0-2.0*x)
        return (y*(1-loss)+loss)
    end
end

function clamp1(x, lo, hi, loss)
    if x<lo
        return x = 1.0
    elseif x>hi
        return x = loss
    end
    return x
end

function ml_0!(X, Y, Z, A, z, q)
  k = size(X)[1]
  # @inbounds @simd for j in 1:900
    @inbounds @simd for i in eachindex(X)
      # Ew0[i] = Ew0[i]+1
      X[i] = Y[i]*A[i]*exp(Z[i]*z)*q
    # end
  end
  nothing
end

function ml!(X, Y, z)
  # @inbounds @simd for j in 1:900
    @inbounds @simd for i in eachindex(X)
      # x = X[i]
      X[i] = X[i]*exp(-Y[i]*z)
    # end
  end
  nothing
end

function ml_1!(X, Y, z)
  # @inbounds @simd for j in 1:900
    @inbounds @simd for i in eachindex(X)
      # Ew0[i] = Ew0[i]+1
      X[i] = Y[i]*z
    # end
  end
  nothing
end

function ml_2!(X, Y, Z, z)
  # @inbounds @simd for j in 1:900
    @inbounds @simd for i in eachindex(X)
      # Ew0[i] = Ew0[i]+1
      X[i] = Y[i]*Z[i]*z
    # end
  end
  nothing
end

function ift!(Et, ift, convw)
  # @inbounds @simd for j in 1:900
    @inbounds @simd for i in eachindex(Et)
      # Ew0[i] = Ew0[i]+1
      Et[i] = ift[i]*convw
    # end
  end
  nothing
end

δ(k::Integer,j::Integer) = k == j ? 1 : 0

function cumtrapz(x::AbstractVector, y::AbstractVector, retarr::AbstractVector)
    @assert length(x) == length(y) "x and y vectors must be of the same length!"
    for i in 2 : length(y)
        retarr[i] = retarr[i-1] + (x[i] - x[i-1]) * (y[i] + y[i-1])
    end
    for i in eachindex(retarr)
        retarr[i] = retarr[i]*1//2
    end
end

function trapz(x::AbstractVector, y::AbstractVector)
    @assert length(x) == length(y) "x and y vectors must be of the same length!"
    retval = 0.0
    for i in 1 : length(y)-1
        retval += (x[i+1] - x[i]) * (y[i] + y[i+1])
    end
    return 1//2 * retval
end

max_stars_repo_name: SabbahMohammed/Dev

max_issues_repo_name: SabbahMohammed/Dev

max_forks_repo_name: SabbahMohammed/Dev


text:
"""
$(TYPEDEF)

A system of partial differential equations.

# Fields
$(FIELDS)

# Example

```julia
using ModelingToolkit

@parameters t x
@variables u(..)
Dxx = Differential(x)^2
Dtt = Differential(t)^2
Dt = Differential(t)

#2D PDE
C=1
eq  = Dtt(u(t,x)) ~ C^2*Dxx(u(t,x))

# Initial and boundary conditions
bcs = [u(t,0) ~ 0.,# for all t > 0
       u(t,1) ~ 0.,# for all t > 0
       u(0,x) ~ x*(1. - x), #for all 0 < x < 1
       Dt(u(0,x)) ~ 0. ] #for all  0 < x < 1]

# Space and time domains
domains = [t ∈ IntervalDomain(0.0,1.0),
           x ∈ IntervalDomain(0.0,1.0)]

pde_system = PDESystem(eq,bcs,domains,[t,x],[u])
```
"""
struct PDESystem <: ModelingToolkit.AbstractSystem
    "The equations which define the PDE"
    eqs
    "The boundary conditions"
    bcs
    "The domain for the independent variables."
    domain
    "The independent variables"
    indvars
    "The dependent variables"
    depvars
    "The parameters"
    ps
    """
    defaults: The default values to use when initial conditions and/or
    parameters are not supplied in `ODEProblem`.
    """
    defaults::Dict
    @add_kwonly function PDESystem(eqs, bcs, domain, indvars, depvars, ps = SciMLBase.NullParameters(), defaults = Dict())
        new(eqs, bcs, domain, indvars, depvars, ps, defaults)
    end
end

Base.getproperty(x::PDESystem, sym::Symbol) = getfield(x, sym)

Base.summary(prob::PDESystem) = string(nameof(typeof(prob)))
function Base.show(io::IO, ::MIME"text/plain", sys::PDESystem)
    println(io,summary(sys))
    println(io,"Equations: ", get_eqs(sys))
    println(io,"Boundary Conditions: ", get_bcs(sys))
    println(io,"Domain: ", get_domain(sys))
    println(io,"Dependent Variables: ", get_depvars(sys))
    println(io,"Independent Variables: ", get_indvars(sys))
    println(io,"Parameters: ", get_ps(sys))
    print(io,"Default Parameter Values", get_defaults(sys))
    return nothing
end


text:
function gaussian_pdf(n, μ, σ)
    x = collect(1:1<<n)
    pl = @. 1 / sqrt(2pi * σ^2) * exp(-(x - μ)^2 / (2 * σ^2))
    pl / sum(pl)
end

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
"""
    DiagnosticVar

The base type for all diagnostic variables.

Various kinds of diagnostic variables (such as `HorizontalAverage`)
are defined as abstract sub-types of this type and implement the group
of functions listed at the end for use by the diagnostics group code
generator.

A particular diagnostic variable is defined as a concrete sub-type of
one of the kinds of diagnostic variables. The type itself is generated as
are the group of functions below.

A diagnostic variable is always associated with a `ClimateMachine`
configuration type, so as to allow the same name to be used by different
configurations.
"""
abstract type DiagnosticVar end

"""
    dv_dg_points_length(::CT, ::Type{DVT}) where {
        CT <: ClimateMachineConfigType,
        DVT <: DiagnosticVar,
    }

Returns an expression that evaluates to the length of the points
dimension of a DG array of the diagnostic variable type.
"""
function dv_dg_points_length end

"""
    dv_dg_points_index(::CT, ::Type{DVT}) where {
        CT <: ClimateMachineConfigType,
        DVT <: DiagnosticVar,
    }

Returns an expression that, within a `@traverse_dg_grid`, evaluates
to an index into the points dimension of a DG array of the diagnostic
variable type.
"""
function dv_dg_points_index end

"""
    dv_dg_elems_length(::CT, ::Type{DVT}) where {
        CT <: ClimateMachineConfigType,
        DVT <: DiagnosticVar,
    }

Returns an expression that evaluates to the length of the elements
dimension of a DG array of the diagnostic variable type.
"""
function dv_dg_elems_length end

"""
    dv_dg_elems_index(::CT, ::Type{DVT}) where {
        CT <: ClimateMachineConfigType,
        DVT <: DiagnosticVar,
    }

Returns an expression that, within a `@traverse_dg_grid`, evaluates to
an index into the elements dimension of a DG array of the diagnostic
variable type.
"""
function dv_dg_elems_index end

"""
    dv_dg_dimnames(::CT, ::Type{DVT}) where {
        CT <: ClimateMachineConfigType,
        DVT <: DiagnosticVar,
    }

Returns a tuple of the names of the dimensions of the diagnostic
variable type for when interpolation is _not_ used.
"""
function dv_dg_dimnames end

"""
    dv_dg_dimranges(::CT, ::Type{DVT}) where {
        CT <: ClimateMachineConfigType,
        DVT <: DiagnosticVar,
    }

Returns a tuple of expressions that evaluate to the range of the
dimensions for the diagnostic variable type for when interpolation is
_not_ used.
"""
function dv_dg_dimranges end

"""
    dv_i_dimnames(::CT, ::Type{DVT}) where {
        CT <: ClimateMachineConfigType,
        DVT <: DiagnosticVar,
    }

Returns a tuple of the names of the dimensions of the diagnostic
variable type or of an expression, for when interpolation is used.
"""
function dv_i_dimnames end

"""
    dv_op(::CT, ::Type{DVT}, lhs, rhs) where {
        CT <: ClimateMachineConfigType,
        DVT <: DiagnosticVar,
    }

Returns an expression that, within a `@traverse_dg_grid`, evaluates to
an assignment of `rhs` to `lhs`, where `rhs` is the implementation of
the diagnostic variable and `lhs` is the appropriate location in the
array containing the computed diagnostic variable values.
"""
function dv_op end

"""
    dv_reduce(::CT, ::Type{DVT}, array_name) where {
        CT <: ClimateMachineConfigType,
        DVT <: DiagnosticVar,
    }

"""
function dv_reduce end

"""
    dv_name(::CT, ::DVT) where {
        CT <: ClimateMachineConfigType,
        DVT <: DiagnosticVar,
    }

Returns the name of the diagnostic variable as a `String`. Generated.
"""
function dv_name end

"""
    dv_attrib(::CT, ::DVT) where {
        CT <: ClimateMachineConfigType,
        DVT <: DiagnosticVar,
    }

Returns a `Dict` of diagnostic variable attributes, primarily for NetCDF.
Generated.
"""
function dv_attrib end

# Default method for variable attributes.
dv_attrib(::ClimateMachineConfigType, ::DiagnosticVar) = Dict()

"""
    dv_args(::CT, ::DVT) where {
        CT <: ClimateMachineConfigType,
        DVT <: DiagnosticVar,
    }

Returns a tuple of `(arg_name, arg_type, slurp, default)` (as returned
by `MacroTools.splitarg` for the arguments specified by the
implementation of the diagnostic variable.
"""
function dv_args end

"""
    dv_scale(::CT, ::DVT) where {
        CT <: ClimateMachineConfigType,
        DVT <: DiagnosticVar,
    }

If scaling was specified for the diagnostic variable, return the
diagnostic variable with which the scaling should be done, otherwise
return `nothing`.
"""
function dv_scale end

# Default method for diagnostic variable scaling.
dv_scale(::ClimateMachineConfigType, ::DiagnosticVar) = nothing

"""
    dv_project(::CT, ::DVT) where {
        CT <: ClimateMachineConfigType,
        DVT <: DiagnosticVar,
    }

Return `true` if the specified diagnostic variable should be
projected after interpolation.
"""
function dv_project end

# Default method for diagnostic variable scaling.
dv_project(::ClimateMachineConfigType, ::DiagnosticVar) = false

####

# Generate a standardized type name from:
# - the configuration type,
# - the diagnostic variable kind, and
# - the diagnostic variable name.
function dv_type_name(dvtype, config_type, name)
    let uppers_in(s) =
            foldl((f, c) -> isuppercase(c) ? f * c : f, String(s), init = "")
        return uppers_in(config_type) *
               "_" *
               uppers_in(dvtype) *
               "_" *
               String(name)
    end
end

"""
    generate_dv_interface(
        dvtype,
        config_type,
        name,
        units,
        long_name,
        standard_name,
    )

Generate the type for a diagnostic variable, add an instance of the
type into `AllDiagnosticVars`, and generate `dv_name` and `dv_attrib`.
"""
function generate_dv_interface(
    dvtype,
    config_type,
    name,
    units = "",
    long_name = "",
    standard_name = "",
)
    dvtypname = Symbol(dv_type_name(dvtype, config_type, name))
    attrib_ex = quote end
    if any(a -> a != "", [units, long_name, standard_name])
        attrib_ex = quote
            dv_attrib(::$config_type, ::$dvtypname) = OrderedDict(
                "units" => $units,
                "long_name" => $long_name,
                "standard_name" => $standard_name,
            )
        end
    end
    quote
        struct $dvtypname <: $dvtype end
        DiagnosticsMachine.AllDiagnosticVars[$config_type][$(String(name))] =
            $dvtypname()
        dv_name(::$config_type, ::$dvtypname) = $(String(name))
        $(attrib_ex)
    end
end

"""
    generate_dv_function(dvtype, config_type, name, impl)

Generate `dv_args` for a diagnostic variable as well as the
implementation function: `dv_<dvtype>`, adding the configuration
type and the diagnostic variable type as the first two parameters
for dispatch.

The implementation _must_ be defined as:
```
f(
    [<component-name>::<component-type>,]
    bl::<balance-law-type>,
    states::States,
    curr_time::Float64,
    cache::Dict{Symbol, Any},
)
```
Where `<component-name>` is the name of the property within the balance
law type, `<component-type>` and `<balance-law-type>` are used for
dispatch, and `cache` may be used to store intermediate computations.
"""
function generate_dv_function(dvtype, config_type, name, impl)
    dvfun = Symbol("dv_", dvtype)
    dvtypname = Symbol(dv_type_name(dvtype, config_type, name))
    @capture(impl, ((args__,),) -> (body_)) ||
        @capture(impl, (args_) -> (body_)) ||
        error("Bad implementation for $(esc(names[1]))")
    split_fun_args = map(splitarg, args)
    fun_args = map(a -> :($(a[1])::$(a[2])), split_fun_args)
    quote
        function dv_args(::$config_type, ::$dvtypname)
            $split_fun_args
        end
        function $dvfun(::$config_type, ::$dvtypname, $(fun_args...))
            $body
        end
    end
end

"""
    generate_dv_scale(dvtype, config_type, name, scale)

Generate `dv_scale` for a diagnostic variable, returning the diagnostic
variable with which it is to be scaled.
"""
function generate_dv_scale(dvtype, config_type, name, scale)
    dvtypname = Symbol(dv_type_name(dvtype, config_type, name))
    sdv = nothing
    if !isnothing(scale)
        cfg_type_name = getfield(ConfigTypes, config_type)
        sdv = DiagnosticsMachine.AllDiagnosticVars[cfg_type_name][String(scale)]
    end
    quote
        dv_scale(::$config_type, ::$dvtypname) = $sdv
    end
end

"""
    generate_dv_project(dvtype, config_type, name, scale)

Generate `dv_project` for a diagnostic variable, returning `true` if
the diagnostic variable is to be projected.
"""
function generate_dv_project(dvtype, config_type, name, project)
    dvtypname = Symbol(dv_type_name(dvtype, config_type, name))
    quote
        dv_project(::$config_type, ::$dvtypname) = $project
    end
end

"""
    States

Composite of the various states, used as a parameter to diagnostic
variable implementations.
"""
struct States{PS, GFS, AS}
    prognostic::PS
    gradient_flux::GFS
    auxiliary::AS
end

# The various kinds of diagnostic variables as well as interfaces to
# create variables of these kinds.
include("pointwise.jl")
include("horizontal_average.jl")


text:
#####
##### 3D modes on pyramids
#####

function basis(elem::Wedge,N,r,s,t)
    V_tri,Vr_tri,Vs_tri = basis(Tri(),N,r,s)
    V_1D,Vt_1D = basis(Line(),N,t)
    Np = (N+1)*(N+2)*(N+1)÷2
    V,Vr,Vs,Vt = ntuple(x->zeros(length(r),Np),4)
    id = 1
    for i = 1:size(V_tri,2)
        for j = 1:size(V_1D,2)
            @. V[:,id]  = V_tri[:,i]  * V_1D[:,j]
            @. Vr[:,id] = Vr_tri[:,i] * V_1D[:,j]
            @. Vs[:,id] = Vs_tri[:,i] * V_1D[:,j]
            @. Vt[:,id] = V_tri[:,i]  * Vt_1D[:,j]
            id += 1
        end
    end
    return V,Vr,Vs,Vt
end

function equi_nodes(elem::Wedge,N)
    t1D = LinRange(-1,1,N+1)
    r_tri,s_tri = equi_nodes(Tri(),N)
    r,t = vec.(meshgrid(r_tri,t1D))
    s,_ = vec.(meshgrid(s_tri,t1D))
    return r,s,t
end

function nodes(elem::Wedge,N)
    t1D,_ = gauss_lobatto_quad(0,0,N)
    r_tri,s_tri = nodes(Tri(),N)
    r,t = vec.(meshgrid(r_tri,t1D))
    s,_ = vec.(meshgrid(s_tri,t1D))
    return r,s,t
end

function quad_nodes(elem::Wedge,N)
    r_tri,s_tri,w_tri = quad_nodes(Tri(),N)
    t1D,wt1D = quad_nodes(Line(),N)
    wrs,wt = vec.(meshgrid(w_tri,wt1D))
    r,t = vec.(meshgrid(r_tri,t1D))
    s,_ = vec.(meshgrid(s_tri,t1D))
    return r,s,t,wrs.*wt
end


text:
###Expressions
ex=:x+:y
@test isa(ex,Expression)
@test has(ex,:x)
ex*=ex
@test has(ex,:x)

#addparse
ex=:x+:y
ap=terms(ex)
@test length(ap)==2
@test expression(ap)==ex #maybe have isequivalent instead to give addparse some room, maybe have == = isequivalent, equivalent=isequivalent?
i=1
for term in ex
	@test term==ap[i]
	i+=1
end

#componify
ex=:x+1-(:x+1)
ex=componify(ex)
@test has(ex,:x)

#simplify
ex=(:x+:x)/2
@test simplify(ex)==:x
@test simplify(:x+:x)==2*:x

#evaluate

#original unsorted tests
ex1=simplify(1+:x*2-1)
ex2=push!(2,:x)
@test ex1==ex2

ex1=:x*:y
ex2=Equations.sumsym(-:x+:x*:y+:x)
@test ex1==ex2

@test Equations.sumsym(1-:x+:x)==1
@test Equations.sumnum(1+:x-1)==:x

@test evaluate(:x+:z,Dict(:x=>3,:z=>4))==7

#Units
l=Physical(:l,:m);t=Physical(:t,:s);v=l/t;
@test v==Physical(:l/:t,:m/:s)


text:
module HHNSGA
    using Pkg
    try
        using Plots
        using UnicodePlots
        using OrderedCollections
        using LinearAlgebra
        using UnicodePlots
        using PyCall
        using NSGAII
        export NSGAII
        #export nS

    catch
        include("install.jl")
    end
    export SNN

    include("../src/SpikingNeuralNetworks.jl")
    include("../src/units.jl")
    include("../src/plot.jl")
    SNN = SpikingNeuralNetworks.SNN
    using Random
    py"""
    import copy
    from neuronunit.optimisation import algorithms
    from neuronunit.optimisation import optimisations
    from neuronunit.optimisation.optimization_management import TSD
    from neuronunit.optimisation.optimization_management import OptMan
    import joblib

    import matplotlib
    import matplotlib.pyplot as plt
    import numpy as np
    matplotlib.get_backend()
    import sys
    import os
    sys.path.append(os.getcwd())
    """
    py"""
    from neuronunit import tests
    from neuronunit.tests import fi
    from simple_with_injection import SimpleModel
    from neo import AnalogSignal
    import quantities as pq
    try:
       from julia import Main
    except:
       import julia
       julia.install()
       from julia import Main

    """
    simple = py"SimpleModel()"
    py"""
    import pickle
    import numpy as np
    import random
    from sciunit.scores.collections import ScoreArray
    from sciunit import TestSuite
    import sciunit

    from neuronunit.tests.fi import RheobaseTestP
    from collections import OrderedDict
    cell_tests = pickle.load(open('multicellular_constraints.p','rb'))
    for test in cell_tests.values():
        if "Rheobase test" in test.keys():
            temp_test = {k:v for k,v in test.items()}
            break
    rt = temp_test["Rheobase test"]
    #rtp = RheobaseTestP(rt.observation)
    JHH = {
    'Vr': -68.9346,
    'Cm': 0.0002,
    'gl': 1.0 * 1e-5,
    'El': -65.0,
    'EK': -90.0,
    'ENa': 50.0,
    'gNa': 0.02,
    'gK': 0.006,
    'Vt': -63.0
    }

    ranges = OrderedDict({k:[v-0.5*np.abs(v),v+0.5*np.abs(v)] for k,v in copy.copy(JHH).items()})
    N = len(JHH)
    model = SimpleModel()
    """
    N = py"N"
    ranges = py"ranges"
    H1=[values(ranges)]
    current_params = py"rt.params"
    simple.attrs = py"JHH"
    #using Py2Jl
    py"""

    def evaluate(test,god):
        model = SimpleModel()
        model.attrs = god
        rt = test['Rheobase test']
        #rtp = RheobaseTestP(rt.observation)

        rheo = rt.generate_prediction(model)
        #    errors = tuple([100.0,100.0,100.0,100.0,100.0,100.0,100.0,100.0])
        #    return errors
        scores_ = []
        #assert 'value' in rheo.keys()

        # for temp_test in test.values():
        #     #print(temp_test.name,temp_test.params)
        #     #print(temp_test.params['injected_square_current']['amplitude'])
        #     #if 'InjectedCurrentAPWidthTest' not in temp_test.params.keys():
        #     #   temp_test.params['injected_square_current'] = {}
        #     if temp_test.name in "InjectedCurrentAPWidthTest":
        #        temp_test.params['injected_square_current']['amplitude'] = rheo['value']
        #     if temp_test.name in "InjectedCurrentAPThresholdTest":
        #        temp_test.params['injected_square_current']['amplitude'] = rheo['value']
        #     if temp_test.name in "InjectedCurrentAPThresholdTest":
        #        temp_test.params['injected_square_current']['amplitude'] = rheo['value']
        #     if temp_test.name in "RheobaseTest":
        #        temp_test.params['injected_square_current']['amplitude'] = rheo['value']

        for key, temp_test in test.items():
            #if key == 'Rheobase test':
            # def to_map(key_val):
            temp_test.score_type = sciunit.scores.ZScore
            if key in ['Injected current AP width test', 'Injected current AP threshold test', 'Injected current AP amplitude test', 'Rheobase test']:
                temp_test.params['amplitude'] = rheo['value']

            try:
                score = temp_test.judge(model)
                score = np.abs(score.log_norm_score)

            except:
                try:
                    score = np.abs(float(score.raw_score))
                except:
                    score = 100#np.inf
            if isinstance(score, sciunit.scores.incomplete.InsufficientDataScore):
                score = 100#np.inf
                #if score == 100:
                print(temp_test.name)
                print(temp_test.params)
                if not rheo['value'] is None:

                    import pdb
                    pdb.set_trace()
            scores_.append(score)

        tt = TestSuite(list(test.values()))
        SA = ScoreArray(tt, scores_)
        errors = tuple(SA.values,)

        return errors
    """

    py"""
    def z_py(evaluate,god):
        cell_tests = pickle.load(open('multicellular_constraints.p','rb'))
        tests = cell_tests[list(cell_tests.keys())[0]]
        SA = evaluate(tests,god)
        return SA
    """

    function references_for_H()
        py"""
        grid = ParameterGrid(ranges)
        genes = []
        for g in grid:
            genes.append(list(g.values()))
        """
        genes = py"genes"
        return genes
    end

    export init_function2
    function init_function2(x)
        py"""
        def one(x):
            if str('x') in locals():
                genes_out = x
                genes_out_dic = OrderedDict({k:genes_out[i] for i,k in enumerate(ranges.keys()) })
            else:
                lower_list = [v[0] for k,v in ranges.items()]
                upper_list = [v[1] for k,v in ranges.items()]
                genes_out = [random.uniform(lower, upper) for lower, upper in zip(lower_list, upper_list)]
                genes_out_dic = OrderedDict({k:genes_out[i] for i,k in enumerate(ranges.keys()) })

            return genes_out, genes_out_dic
        """
        decoded = x

        try
           decoded = NSGAII.decode(x, HHNSGA.bc)
        catch
           decoded = x
           println("no")
        end
        genes_out, genes_out_dic = py"one"(decoded)
        #@show(genes_out)

        bincoded = NSGAII.encode(genes_out, bc)
        #@show(bincoded)
        return bincoded, genes_out_dic#py"genes_out",py"genes_out_dic"
    end
    py"""
    lower_list = [v[0] for k,v in ranges.items()]
    upper_list = [v[1] for k,v in ranges.items()]
    #print(len(lower_list))
    """

    const bc = BinaryCoding(9, [:Int,:Int,:Int,:Int,:Int,:Int,:Int,:Int,:Int], py"lower_list", py"upper_list")
    export init_function
    function init_function()
        py"""
        def one():
            import random
            lower_list = [v[0] for k,v in ranges.items()]
            upper_list = [v[1] for k,v in ranges.items()]
            genes_out = [random.uniform(lower, upper) for lower, upper in zip(lower_list, upper_list)]
            genes_out_dic = OrderedDict({k:genes_out[i] for i,k in enumerate(ranges.keys()) })
            return genes_out, genes_out_dic
        """
        genes_out, genes_out_dic = py"one"()
        #@show(genes_out)
        bincoded = NSGAII.encode(genes_out, bc)
        #@show(bincoded)
        return bincoded
    end

    export z
    function z(x)
        genes_out,god = init_function2(x)
        #x = genes_out
        # The slow part can it be done in parallel.
        contents = py"z_py"(py"evaluate",god)
        @show(contents)
        println("gene fitness calculated")
        return contents
    end

end

max_stars_repo_name: russelljjarvis/juliaSNN

max_issues_repo_name: russelljjarvis/juliaSNN

max_forks_repo_name: russelljjarvis/juliaSNN


text:
"""
Materialize a `Selection` into a `DataFrame`.
"""
function materialize(node::Selection)
    mask = generate_column(node.source, node.predicate)
    @assert eltype(mask) <: Union{Missing, Bool}
    materialize(node.source)[coalesce.(mask, false), :]
end


text:
"""
```julia
AGeom=Vector{AGeomss,1}
```
An acquisition geometry, bundled into the mutable `AGeom` type,
has to be specified in order to
create any `Expt` variable.
Acquisition has supersources, sources and receivers.
For `AGeom` is a mutable type to define source-receiver geometry.
Each supersource has `ns` sources that are
injected (or active) simultaneously, and
a set of `nr` receivers to
record.
`AGeomss` is a subtype for each supersource, therefore:
This package provides tools to easily create commonly-used `AGeom` instances,
after deciding on a spatial grid `mgrid` for the experiment.

## Indexing and Fields
Let `ageom` be an instance of this type, then fields
can be accessed using:

* `ageom[i]` : acquisition for ith supersource
* `ageom[i].sx` : x positions of sources
* `ageom[i].sz` : z positions of sources
* `ageom[i].rx` : x positions of receivers
* `ageom[i].rz` : z positions of receivers
* `ageom[i].ns` : number of sources
* `ageom[i].nr` : number of receivers

```julia
ageom=AGeom(mgrid, attrib, SSrcs(10), Recs(20))
```
Return pre-defined acquisition geometries (with 10 sources and 20 receivers), based on `attrib`, on an input mesh `mgrid`.
The sources and receivers positions are same for all supersources, and are not placed on the edges of the mesh.
Choose `attrib::Symbol`
* `=:xwell` cross-well acquisition;
* `=:surf` surface acquisition;
* `=:vsp` vertical seismic profiling;
* `=:rvsp`  reverse vertical seismic profiling;
* `=:downhole` downhole sources and receivers.
"""
AGeom=Array{AGeomss,1}

function Array{AGeomss,1}(mgrid::AbstractArray{T}, attrib::Symbol, ss::SSrcs=SSrcs(1), r::Recs=Recs(10)) where {T<:StepRangeLen}
	otx=(0.9*mgrid[2][1]+0.1*mgrid[2][end]); ntx=(0.1*mgrid[2][1]+0.9*mgrid[2][end]);
	otwx=(0.95*mgrid[2][1]+0.05*mgrid[2][end]); ntwx=(0.05*mgrid[2][1]+0.95*mgrid[2][end]);
	otz=(0.9*mgrid[1][1]+0.1*mgrid[1][end]); ntz=(0.1*mgrid[1][1]+0.9*mgrid[1][end]);
	otwz=(0.95*mgrid[1][1]+0.05*mgrid[1][end]); ntwz=(0.05*mgrid[1][1]+0.95*mgrid[1][end]);
	quatx = (0.75*mgrid[2][1]+0.25*mgrid[2][end]); quatz = (0.75*mgrid[1][1]+0.25*mgrid[1][end])
	tquatx = (0.25*mgrid[2][1]+0.75*mgrid[2][end]); tquatz = (0.25*mgrid[1][1]+0.75*mgrid[1][end])
	halfx = 0.5*(mgrid[2][1]+mgrid[2][end]);	halfz = 0.5*(mgrid[1][1]+mgrid[1][end]);
	ageom=AGeom(mgrid, ss, Srcs(1), r)
	if(attrib == :xwell)
		update!(ageom, SSrcs(), [otz,otx], [ntz,otx])
		update!(ageom, Recs(), [otz,ntx], [ntz,ntx])
	elseif(attrib == :surf)
		update!(ageom, SSrcs(), [otwz,otx], [otwz,ntx])
		update!(ageom, Recs(), [otwz,otx], [otwz,ntx])
	elseif(attrib == :vsp)
		update!(ageom, SSrcs(), [otwz,otx], [otwz,ntx])
		update!(ageom, Recs(), [quatz,otx], [ntz,ntx])
	elseif(attrib == :rvsp)
		update!(ageom, SSrcs(), [quatz,otx], [ntz,otx])
		update!(ageom, Recs(),  [otwz,otx], [otwzz,ntx])
	elseif(attrib == :downhole)
		update!(ageom, SSrc(), [quatz,quatx], [otz,quatx])
		update!(ageom, Recs(),  [quatz,quatx], [otz,quatx])
	else
		error("invalid attrib")
	end
	return ageom
end


text:
using LighthouseFlux
using Documenter

makedocs(; modules=[LighthouseFlux], sitename="LighthouseFlux",
         authors="Beacon Biosignals and other contributors",
         pages=["API Documentation" => "index.md"])

deploydocs(; repo="github.com/beacon-biosignals/LighthouseFlux.jl.git", devbranch="main")


text:
"""
    MixedModelPermutation{T<:AbstractFloat}

Object returned by the various permutation methods with fields
- `fits`: the parameter estimates from the permutation replicates as a vector of named tuples.
- `λ`: `Vector{LowerTriangular{T,Matrix{T}}}` containing copies of the λ field from `ReMat` model terms
- `inds`: `Vector{Vector{Int}}` containing copies of the `inds` field from `ReMat` model terms
- `lowerbd`: `Vector{T}` containing the vector of lower bounds (corresponds to the identically named field of [`OptSummary`](@ref))
- `fcnames`: NamedTuple whose keys are the grouping factor names and whose values are the column names

The schema of `fits` is, by default,
```
Tables.Schema:
 :objective  T
 :σ          T
 :β          NamedTuple{β_names}{NTuple{p,T}}
 :se         StaticArrays.SArray{Tuple{p},T,1,p}
 :θ          StaticArrays.SArray{Tuple{k},T,1,k}
```
where the sizes, `p` and `k`, of the `β` and `θ` elements are determined by the model.

Characteristics of the permutation replicates can be extracted as properties.  The `σs` and
`σρs` properties unravel the `σ` and `θ` estimates into estimates of the standard deviations
and correlations of the random-effects terms.
"""
struct MixedModelPermutation{T<:AbstractFloat} <: MixedModels.MixedModelFitCollection{T}
    fits::Vector
    λ::Vector{<:Union{LowerTriangular{T},Diagonal{T}}}
    inds::Vector{Vector{Int}} # XXX should this be explicitly Int64 or <:Integer?
    lowerbd::Vector{T}
    fcnames::NamedTuple
end

"""
_perm2boot(p::MixedModelPermutation)

Non copying conversion to MixedModels.MixedModelBootstrap.
"""
_perm2boot(p::MixedModelPermutation) = MixedModelBootstrap(p.fits, p.λ, p.inds, p.lowerbd, p.fcnames)


text:
struct MCMCResult{T} #<: Distributions.ContinuousUnivariateDistribution
    Probs::Matrix{T}
    RevCholWisharts::ScatteredMatrix{T,2,RevCholWishart{T}}
    CholInvWisharts::ScatteredMatrix{T,2,CholInvWishart{T}}
end
function MCMCResult(NG, iter, ::Type{T} = Float32) where T
    probs = Array{T}(undef, NG, iter)
    revcholwisharts = ScatteredMatrix{T,2,RevCholWishart{T}}(undef, NG, iter)
    cholinvwisharts = ScatteredMatrix{T,2,CholInvWishart{T}}(undef, NG, iter)
    MCMCResult(probs, revcholwisharts, cholinvwisharts)
end

function MCMCResults(NG, iter, ::Type{T} = Float32, nthread = nthreads()) where T
    mcmc_results = Vector{MCMCResult{T}}(undef, nthread)
    @threads for n ∈ 1:nthread
        mcmc_results[n] = MCMCResult(NG, iter, T)
    end
    mcmc_results
end

function Base.cat(a::MCMCResult{T}, b::MCMCResult{T}) where T
    MCMCResult(
        hcat(a.Probs,           b.Probs),
        ScatteredMatrix{T,2,RevCholWishart{T}}(cat(a.RevCholWisharts.data, b.RevCholWisharts.data, dims = 2)),
        ScatteredMatrix{T,2,CholInvWishart{T}}(cat(a.CholInvWisharts.data, b.CholInvWisharts.data, dims = 2))
    )
end

function Base.cat(a::MCMCResult{T}, b::MCMCResult{T}, iter) where T
    MCMCResult(
        hcat(@view(a.Probs[:,1:iter]),           @view(b.Probs[:,1:iter])),
        ScatteredMatrix{T,2,RevCholWishart{T}}(cat(@view(a.RevCholWisharts[:,1:iter,:]), @view(b.RevCholWisharts[:,1:iter,:]), dims = 2)),
        ScatteredMatrix{T,2,CholInvWishart{T}}(cat(@view(a.CholInvWisharts[:,1:iter,:]), @view(b.CholInvWisharts[:,1:iter,:]), dims = 2))
    )
end

@generated function Distributions.logpdf(res::MCMCResult{T}, x::SVector{3}) where {T}
    # W, Wshift = VectorizationBase.pick_vector_width_shift(T)
    # V = SVec{W,T}
    quote
        probs = res.Probs
        rcw = res.RevCholWisharts
        N = length(probs)
        # vrcw = vectorizable(rcw)
        # vprobs = vectorizable(probs)
        p = zero(Float64)
        @vectorize Float64 for n ∈ 1:N
            vRCW = rcw[n]
            vRCWx = vRCW * x
            lkernel = SLEEFwrap.log(1.0 + (vRCWx' * vRCWx))
            ld = logdetw(vRCW)
            ν = extract_ν(vRCW)
            nexponent = 0.5 * ν + 1.5
            base = SLEEFwrap.log(probs[n]) + ld + SLEEFwrap.lgamma(nexponent) - SLEEFwrap.lgamma(0.5ν) - 1.5*SLEEFwrap.log(ν)

            p += SLEEFwrap.exp(base - nexponent * lkernel)
        end
        # log(p) - log($T(size(probs,2))) - $(T(1.5log(π)))
        log(p) - log(size(probs,2)) - $(1.5log(π))
    end
end
@generated function logpdf2(res::MCMCResult{T}, x::SVector{3}) where {T}
    # W, Wshift = VectorizationBase.pick_vector_width_shift(T)
    # V = SVec{W,T}
    quote
        probs = res.Probs
        rcw = res.RevCholWisharts
        N = length(probs)
        # vrcw = vectorizable(rcw)
        # vprobs = vectorizable(probs)
        p = zero(Float64)
        for n ∈ 1:N
            rcw32 = rcw[n]
            vRCW = RevCholWishart{Float64}(Float64.(rcw32.data))
            vRCWx = vRCW * x
            lkernel = log(1.0 + (vRCWx' * vRCWx))
            ld = logdetw(vRCW)
            ν = extract_ν(vRCW)
            nexponent = 0.5 * ν + 1.5
            base = log(Float64(probs[n])) + ld + lgamma(nexponent) - lgamma(0.5ν) - 1.5*log(ν)

            pn = exp(base - nexponent * lkernel)
            # @show pn
            p += pn

        end
        # log(p) - log($T(size(probs,2))) - $(T(1.5log(π)))
        log(p) - log(size(probs,2)) - $(1.5log(π))
    end
end

@generated function Random.randn(rng::AbstractRNG, ::Type{SVec{W,T}}) where {W,T}
    :(SVec(@ntuple $W w -> VE(randn(rng, $T))))
end

@generated function sample_distances!(rng::AbstractRNG, wdistances::WeightedSamples{T}, res::MCMCResult{T}) where {T}
    W = 2*VectorizationBase.pick_vector_width(T)
    V = SVec{W,T}
    quote
        # @show $V
        # ciw = res.CholInvWisharts
        copyto!(wdistances.weights, res.Probs)
        ciw = res.CholInvWisharts
        distances = wdistances.distances
        # ciw = res.CholInvWisharts
        N = length(distances)
        vciw = vectorizable(ciw)
        # @show typeof(ciw), typeof(vciw)
        ptr_distances = vectorizable(distances)
        for n ∈ 1:($W):N+$(1-W)
            vZ = SVector{3}(randn(rng, $V),randn(rng, $V),randn(rng, $V))
            vCIW = vload($V, vciw, n)
            vx = vCIW * vZ
            # @show extract_ν(vCIW)
            u = randchisq(rng, extract_ν(vCIW))
            vstore!(ptr_distances + n - 1, sqrt( vx' * vx / u ))
        end
        for n ∈ N+1-(N % $W):N
            z = SVector{3}(randn(rng, T),randn(rng, T),randn(rng, T))
            CIW = ciw[n]
            x = CIW * z
            u = randchisq(rng, extract_ν(CIW))
            distances[n] = sqrt( x' * x / u )
        end
        distances
    end
end

@generated function sample_misses!(table::Array{SVec{W,T},3}, revcholwisharts, probs, r1, v1, r2, v2, c1, c2, N) where {T,W}
    # W = VectorizationBase.pick_vector_width(T)
    quote
        rng = VectorizedRNG.GLOBAL_vPCG
        # Transpose U; broadcast the individual elements.
        chol_c2 = chol(c2)
        for i ∈ eachindex(probs)
            p = vbroadcast(SVec{$W, $T}, probs[i])
            c2_temp = xtx(UpperTriangle3(revcholwisharts[i]) * chol_c2)
            Σ, δr₀ = CovarianceRealism.RIC_to_2D(r1, v1, c1, r2, v2, c2_temp)
            U = chol(Σ)

            L11 = vbroadcast(SVec{$W,$T}, U[1,1])
            L21 = vbroadcast(SVec{$W,$T}, U[1,2])
            L22 = vbroadcast(SVec{$W,$T}, U[2,2])
            vδr₀ = vbroadcast(SVec{$W,$T}, δr₀)
            ν = extract_ν(revcholwisharts[i])
            νinv = 1/ν
            vν = vbroadcast(SVec{$(2W),$T}, (ν < one($T) ? muladd(T(0.5),ν,one($T)) : T(0.5)*ν ))
            @inbounds for n ∈ 1:N
                # Sampling random normals in batches of 4W is fastest on the tested architecture.
                z1_and_2 = randn(rng, Vec{$(4W),$T})

                # We want coordinates z1 and z2. Because we sampled 4W, we have 2W of each z1 and z2.
                # In some versions of Julia there is a bug that causes a more typical version of constructing these tuples:
                # z1_1 = ntuple(Val($W)) do i x1_and_2[i] end
                # or equivalently
                # z1_1 = ntuple(i -> x1_and_2[i], Val($W))
                # So we construct the tuple explicitly here. This is fast in all versions of Julia > 1.
                z1_1 = SVec($(Expr(:tuple, [:(z1_and_2[$w]) for w ∈    1: W]...)))
                z1_2 = SVec($(Expr(:tuple, [:(z1_and_2[$w]) for w ∈  W+1:2W]...)))
                z2_1 = SVec($(Expr(:tuple, [:(z1_and_2[$w]) for w ∈ 2W+1:3W]...)))
                z2_2 = SVec($(Expr(:tuple, [:(z1_and_2[$w]) for w ∈ 3W+1:4W]...)))

                u1_and_2 = vbroadcast(SVec{$(2W),$T}, T(0.5)) /
                    (
                        ν < one($T) ?
                            SIMDPirates.evmul(SLEEF.exp(-SVec(randexp(rng, Vec{$(2W),$T}))*νinv), randgamma_g1(rng, vν)) :
                            randgamma_g1(rng, vν)
                    ) |> SIMDPirates.extract_data

                u1 = sqrt(SVec($(Expr(:tuple, [:(u1_and_2[$w]) for w ∈   1: W]...))))
                u2 = sqrt(SVec($(Expr(:tuple, [:(u1_and_2[$w]) for w ∈ W+1:2W]...))))
                # Multiply Calculate L * z = x.
                # Additrionally, scaling by u, the random chi square, to sample from the multivariate t distribution.
                # Additionally, subtract the distance δr₀, which we treat as [δr₀, 0]
                x1_1 = u1 * L11 * z1_1 - vδr₀
                x1_2 = u2 * L11 * z1_2 - vδr₀

                x2_1 = u1 * (L21 * z1_1 + L22 * z2_1)
                x2_2 = u2 * (L21 * z1_2 + L22 * z2_2)

                # Calculate distance, multiply 1000 to translate RevCholWishartsfrom kilometers to meters.
                δ_1 = SIMDPirates.evmul(sqrt( x1_1*x1_1 + x2_1*x2_1 ), T(10^3))
                δ_2 = SIMDPirates.evmul(sqrt( x1_2*x1_2 + x2_2*x2_2 ), T(10^3))
                vi = vbroadcast(SVec{$W,$T}, $T(100))
                l_1_5 = δ_1 < vi
                l_2_5 = δ_2 < vi
                # Every 5 times, we check if we can break early.
                for i ∈ 50:-5:1
                    ((l_1_5 === SVec($(Expr(:tuple, [:(Core.VecElement(false)) for w ∈ 1:W]...)))) && (l_2_5 === SVec($(Expr(:tuple, [:(Core.VecElement(false)) for w ∈ 1:W]...))))) && break
                    # (any(l_1_5) || any(l_2_5)) || break
                    Base.Cartesian.@nexprs 5 j -> begin
                        vi_j = vbroadcast(SVec{$W,$T}, $T(i+1-j))
                        ti_1_j = table[1,j,i+1-j]
                        l_1_j = δ_1 < vi_j
                        ti_1_j = vifelse(l_1_j, ti_1_j + p, ti_1_j)
                        table[1,j,i+1-j] = ti_1_j
                        ti_2_j = table[2,j,i+1-j]
                        l_2_j = δ_2 < vi_j
                        ti_2_j = vifelse(l_2_j, ti_2_j + p, ti_2_j)
                        table[2,j,i+1-j] = ti_2_j
                    end
                end
            end

        end

    end
end

function sample_Pc!(rng::AbstractRNG, wpc_array::WeightedSamples{T1}, res::MCMCResult{T2}, r1, v1, c1, r2, v2, c2, HBR) where {T1,T2}
    sample_Pc!(rng, wpc_array, res.RevCholWisharts, res.Probs, r1, v1, c1, r2, v2, c2, HBR)
end

function sample_Pc!(rng::AbstractRNG, wpc_array::WeightedSamples{T1}, rcws::ScatteredMatrix{T2,2,RevCholWishart{T2}}, probs, r1, v1, c1, r2, v2, c2, HBR) where {T1,T2}
    chol_c2 = chol(c2)
    pc_array = wpc_array.distances
    copyto!(wpc_array.weights, probs)
    for i ∈ eachindex(pc_array)
        rcws64 = RevCholWishart{Float64}(Float64.(rcws[i].data))
        c2_temp = xtx(randinvwishartfactor(rng, rcws64) * chol_c2)
        pc_array[i] = pc2dfoster_RIC(r1, v1, c1, r2, v2, c2_temp, HBR)
    end
    wpc_array
end
function sample_Pc!(rng::AbstractRNG, pc_array::AbstractVector{T1}, rcws::ScatteredMatrix{T2,2,RevCholWishart{T2}}, r1, v1, c1, r2, v2, c2, HBR) where {T1,T2}
    chol_c2 = chol(c2)
    for i ∈ eachindex(pc_array)
        rcws64 = RevCholWishart{Float64}(Float64.(rcws[i].data))
        c2_temp = xtx(randinvwishartfactor(rng, rcws64) * chol_c2)
        pc_array[i] = pc2dfoster_RIC(r1, v1, c1, r2, v2, c2_temp, HBR)
    end
    pc_array
end
# This method also calculates covariance ratios
function sample_Pc!(rng::AbstractRNG, pc_array::AbstractMatrix{T1}, rcws::ScatteredMatrix{T2,2,RevCholWishart{T2}}, r1, v1, c1, r2, v2, c2, HBR) where {T1,T2}
    chol_c2 = chol(c2)
    ldc = logdet(chol_c2)
#    @show c2
    for i ∈ 1:size(pc_array,1)
        # RevCholWishart is precision factor
        rcws64 = RevCholWishart{Float64}(Float64.(rcws[i].data))
        rw = randinvwishartfactor(rng, rcws64)
        cov_factor = rw * chol_c2
        c11 = cov_factor[1,1]
        c22 = cov_factor[2,2]
        c33 = cov_factor[3,3]
        if !((c11 > 0) && (c22 > 0) && (c33 > 0))
            @show c11, c22, c33
            @show rcws64
            @show rw
            @show chol_c2
            @show cov_factor
        end

        pc_array[i,2] = log(cov_factor[1,1]) + log(cov_factor[2,2]) + log(cov_factor[3,3]) - ldc
#        pc_array[i,2] = logdet(cov_factor) - ldc
        c2_temp = xtx(cov_factor)
#        if i <= 4
#            @show rcws64.data
#            @show rw
#            @show c2_temp
#        end
        pc_array[i,1] = pc2dfoster_RIC(r1, v1, c1, r2, v2, c2_temp, HBR)
    end
    pc_array
end