text:
"""
	triangulate(cmd0::String="", arg1=[]; kwargs...)

Reads one or more ASCII [or binary] files (or standard input) containing x,y[,z] and performs Delaunay
triangulation, i.e., it find how the points should be connected to give the most equilateral
triangulation possible.

Full option list at [`triangulate`](http://gmt.soest.hawaii.edu/doc/latest/triangulate.html)

Parameters
----------
- **C** : **slope_grid** : -- Number --

    Read a slope grid (in degrees) and compute the propagated uncertainty in the
    bathymetry using the CURVE algorithm
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/triangulate.html#c)
- **D** : **derivatives** : -- Str --

    Take either the x- or y-derivatives of surface represented by the planar facets (only used when G is set).
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/triangulate.html#a)
- **E** : **empty** : -- Bool or [] --

    Set the value assigned to empty nodes when G is set [NaN].
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/triangulate.html#e)
- **G** : **grid** : **outgrid** : -- Str or [] --

    Use triangulation to grid the data onto an even grid (specified with R I).
    Append the name of the output grid file.
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/triangulate.html#g)
- **I** : **inc** : -- Str or Number --

    *x_inc* [and optionally *y_inc*] is the grid spacing.
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/triangulate.html#i)
- $(GMT.opt_J)
- **M** : **network** : -- Bool or [] --

    Output triangulation network as multiple line segments separated by a segment header record.
    [`-M`](http://gmt.soest.hawaii.edu/doc/latest/triangulate.html#m)
- **N** : **ids** : -- Bool or [] --

    Used in conjunction with G to also write the triplets of the ids of all the Delaunay vertices
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/triangulate.html#n)
- **Q** : **voronoi** : -- Str or [] --

    Output the edges of the Voronoi cells instead [Default is Delaunay triangle edges]
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/triangulate.html#q)
- $(GMT.opt_R)
- **S** : **triangles** : -- Bool or [] --

    Output triangles as polygon segments separated by a segment header record. Requires Delaunay triangulation.
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/triangulate.html#s)
- **T** : **edges** : -- Bool or [] --

    Output edges or polygons even if gridding has been selected with the G option
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/triangulate.html#t)
- $(GMT.opt_V)
- **Z** : **xyz** : **triplets** : -- Bool or [] --

    [`-Z`](http://gmt.soest.hawaii.edu/doc/latest/triangulate.html#z)
- $(GMT.opt_bi)
- $(GMT.opt_bo)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_r)
- $(GMT.opt_swap_xy)
"""
function triangulate(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && return monolitic("triangulate", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)
	cmd, = parse_R("", d)
	cmd, = parse_J(cmd, d)
	cmd  = parse_V_params(cmd, d)
	cmd, = parse_bi(cmd, d)
	cmd, = parse_bo(cmd, d)
	cmd, = parse_di(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, = parse_i(cmd, d)
	cmd, = parse_r(cmd, d)
	cmd, = parse_swap_xy(cmd, d)

	cmd = add_opt(cmd, 'C', d, [:C :slope_grid])
	cmd = add_opt(cmd, 'D', d, [:D :derivatives])
	cmd = add_opt(cmd, 'E', d, [:E :empty])
    cmd = add_opt(cmd, 'G', d, [:G :grid :outgrid])
	cmd = add_opt(cmd, 'I', d, [:I :inc])
	cmd = add_opt(cmd, 'M', d, [:N :network])
	cmd = add_opt(cmd, 'N', d, [:N :ids])
	cmd = add_opt(cmd, 'Q', d, [:Q :voronoi])
	cmd = add_opt(cmd, 'S', d, [:S :triangles])
	cmd = add_opt(cmd, 'T', d, [:T :edges])
	cmd = add_opt(cmd, 'Z', d, [:Z :xyz :triplets])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "triangulate", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
triangulate(arg1::Array, cmd0::String=""; kw...) = triangulate(cmd0, arg1; kw...)


text:
@testset "{{ title }} Macro Tests" begin

  @test isdefined({{ app }}, Symbol("@{{ name }}")) == true

end

ext: jl


text:
"""
`module  JAC.Statistical`
    ... a submodel of JAC that contains all structs and methods to deal with time-dependent statistical tensors.
"""
module Statistical

    using ..ManyElectron

    """
    `struct  Statistical.ResonanceR`  ... defines a type for a resonance state in the continuum with a well-defined bound-ionic core,
             one or several electrons in the continuum, a widths as well as a loss rate due to additional decay processes.

        + isBound            ::Bool                ... True if it just refers to a bound state with no electron in the continuum.
        + ionLevel           ::Level               ... Level of the bound-state core.
        + widths             ::Float64             ... Widths of the resonance state.
        + lossRate           ::Float64             ... Loss rate due to processes that are not considered explicitly in the time evolution.
    """
    struct ResonanceR
        isBound              ::Bool
        ionLevel             ::Level
        widths               ::Float64
        lossRate             ::Float64
    end

    """
    `Statistical.ResonanceR()`  ... constructor for an `empty` instance of Statistical.ResonanceR().
    """
    function ResonanceR()
        ResonanceR( true, Level(), 0., 0. )
    end

    # `Base.show(io::IO, resonance::Statistical.ResonanceR)`  ... prepares a proper printout of the variable resonance::Statistical.ResonanceR.
    function Base.show(io::IO, resonance::Statistical.ResonanceR)
        println(io, "isBound:              $(resonance.isBound)  ")
        println(io, "ionLevel:             $(resonance.ionLevel)  ")
        println(io, "widths:               $(resonance.widths)  ")
        println(io, "lossRate:             $(resonance.lossRate)  ")
    end

    """
    `struct  Statistical.Tensor`  ... defines a type for a statistical tensor of given rank k, projection q and w.r.t. two resonances.

        + k                  ::AngularJ64          ... Rank of the tensor
        + q                  ::AngularM64          ... Projection of the tensor rank.
        + a                  ::ResonanceR          ... Level/resonance a with total symmetry alpha J.
        + b                  ::ResonanceR          ... Level/resonance a with total symmetry alpha' J'.
        + value              ::Complex{Float64}    ... value of the statistical tensor
    """
    struct Tensor
        k                    ::AngularJ64
        q                    ::AngularM64
        a                    ::ResonanceR
        b                    ::ResonanceR
        value                ::Complex{Float64}
    end

    """
    `Statistical.Tensor()`  ... constructor for an `empty` instance of Statistical.Tensor().
    """
    function Tensor()
        Tensor( AngularJ64(0), AngularM64(0), ResonanceR(), ResonanceR(), Complex(0.) )
    end

    # `Base.show(io::IO, tensor::Statistical.Tensor)`  ... prepares a proper printout of the variable tensor::Statistical.Tensor.
    function Base.show(io::IO, tensor::Statistical.Tensor)
        println(io, "k:              $(tensor.k)  ")
        println(io, "q:              $(tensor.q)  ")
        println(io, "a:              $(tensor.a)  ")
        println(io, "b:              $(tensor.b)  ")
        println(io, "value:          $(tensor.value)  ")
    end

end # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule TtH_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("TtH")
JLLWrappers.@generate_main_file("TtH", UUID("f7b99bbd-2e5d-508a-a6d9-3364a7ade44d"))
end  # module TtH_jll


text:
julia> print_each(X...) = for x in X; println(x); end

julia> print_each(1, "hello", 23.4)
1
hello
23.4

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@testset "Find hblocks" begin
    st = raw"""
        Some text then {{ fill v1 }} and
        {{ if b1 }}
        show stuff here {{ fill v2 }}
        {{ else }}
        show other stuff
        {{ end }}"""

    blocks, tokens = explore_h_steps(st)[:ocblocks]

    @test blocks[1].ss == "{{ fill v1 }}"
    @test blocks[2].ss == "{{ if b1 }}"
    @test blocks[3].ss == "{{ fill v2 }}"
    @test blocks[4].ss == "{{ else }}"
    @test blocks[5].ss == "{{ end }}"
end

@testset "Qual hblocks" begin
    st = raw"""
        Some text then {{ fill v1 }} and
        {{ if b1 }}
        show stuff here {{ fill v2 }}
        {{ else }}
        show other stuff
        {{ end }}
        """

    qblocks, = explore_h_steps(st)[:qblocks]

    @test qblocks[1].fname == "fill"
    @test qblocks[1].params == ["v1"]
    @test qblocks[2].vname == "b1"
    @test qblocks[3].fname == "fill"
    @test qblocks[3].params == ["v2"]
    @test typeof(qblocks[4]) == JuDoc.HElse
    @test typeof(qblocks[5]) == JuDoc.HEnd
end

@testset "Cond block" begin
    st = raw"""
        Some text then {{ fill v1 }} and
        {{ if b1 }}
        show stuff here {{ fill v2 }}
        {{ elseif b2 }}
        other stuff
        {{ else }}
        show other stuff
        {{ end }}
        final text
        """

    cblocks,cdblocks,cpblocks,qblocks = explore_h_steps(st)[:cblocks]

    @test cblocks[1].init_cond == "b1"
    @test cblocks[1].sec_conds == ["b2"]
    @test cblocks[1].actions[1] == "\nshow stuff here {{ fill v2 }}\n"
    @test cblocks[1].actions[2] == "\nother stuff\n"
    @test cblocks[1].actions[3] == "\nshow other stuff\n"
end

@testset "Merge blocks" begin
    st = raw"""
        Some text then {{ fill v1 }} and
        {{ if b1 }}
        show stuff here {{ fill v2 }}
        {{ elseif b2 }}
        other stuff
        {{ else }}
        show other stuff
        {{ end }}
        final text
        """

    hblocks, = explore_h_steps(st)[:hblocks]

    @test hblocks[1].ss == "{{ fill v1 }}"
    @test hblocks[2].ss == "{{ if b1 }}\nshow stuff here {{ fill v2 }}\n{{ elseif b2 }}\nother stuff\n{{ else }}\nshow other stuff\n{{ end }}"
end


text:
module Trajectories

export
    fit_cubic,
    fit_quintic,
    Constant,
    Interpolated,
    Piecewise,
    BezierCurve,
    PointTrajectory,
    FreeVectorTrajectory,
    SE3Trajectory,
    PolyDerivEvaluator

using StaticUnivariatePolynomials
using LinearAlgebra
using StaticArrays
using Rotations
using RigidBodyDynamics

using Base: tail
using StaticUnivariatePolynomials: constant, derivative, exponential_integral, _map

const SUP = StaticUnivariatePolynomials

const BezierCurve = BernsteinPolynomial

@noinline function throw_trajectory_domain_error(x, x0, xf)
    throw(DomainError(x, "Trajectory evaluated outside of range [$x0, $xf]"))
end

zero_tangent(x::Number) = zero(x)
zero_tangent(x::AbstractVector) = zero(x)
zero_tangent(x::Rotation{3, T}) where {T} = zero(SVector{3, T})

include("fit_polynomial.jl")
include("constant.jl")
include("poly_deriv_evaluator.jl")
include("interpolated.jl")
include("piecewise.jl")
include("point_vector.jl")
include("se3.jl")

end


text:
"""
    HMC(n_iters::Int, epsilon::Float64, tau::Int)

Hamiltonian Monte Carlo sampler.

Arguments:

- `n_iters::Int` : The number of samples to pull.
- `epsilon::Float64` : The leapfrog step size to use.
- `tau::Int` : The number of leapfrop steps to use.

Usage:

```julia
HMC(1000, 0.05, 10)
```

Example:

```julia
# Define a simple Normal model with unknown mean and variance.
@model gdemo(x) = begin
    s ~ InverseGamma(2,3)
    m ~ Normal(0, sqrt(s))
    x[1] ~ Normal(m, sqrt(s))
    x[2] ~ Normal(m, sqrt(s))
    return s, m
end

sample(gdemo([1.5, 2]), HMC(1000, 0.05, 10))
```

Tips:

- If you are receiving gradient errors when using `HMC`, try reducing the
`step_size` parameter.

```julia
# Original step_size
sample(gdemo([1.5, 2]), HMC(1000, 0.1, 10))

# Reduced step_size.
sample(gdemo([1.5, 2]), HMC(1000, 0.01, 10))
```
"""
mutable struct HMC{AD, T} <: StaticHamiltonian{AD}
    n_iters   ::  Int       # number of samples
    epsilon   ::  Float64   # leapfrog step size
    tau       ::  Int       # leapfrog step number
    space     ::  Set{T}    # sampling space, emtpy means all
    gid       ::  Int       # group ID
end
HMC(args...) = HMC{ADBackend()}(args...)
function HMC{AD}(epsilon::Float64, tau::Int, space...) where AD
    _space = isa(space, Symbol) ? Set([space]) : Set(space)
    return HMC{AD, eltype(_space)}(1, epsilon, tau, _space, 0)
end
function HMC{AD}(n_iters::Int, epsilon::Float64, tau::Int) where AD
    return HMC{AD, Any}(n_iters, epsilon, tau, Set(), 0)
end
function HMC{AD}(n_iters::Int, epsilon::Float64, tau::Int, space...) where AD
    _space = isa(space, Symbol) ? Set([space]) : Set(space)
    return HMC{AD, eltype(_space)}(n_iters, epsilon, tau, _space, 0)
end
function HMC{AD1}(alg::HMC{AD2, T}, new_gid::Int) where {AD1, AD2, T}
    return HMC{AD1, T}(alg.n_iters, alg.epsilon, alg.tau, alg.space, new_gid)
end
function HMC{AD, T}(alg::HMC, new_gid::Int) where {AD, T}
    return HMC{AD, T}(alg.n_iters, alg.epsilon, alg.tau, alg.space, new_gid)
end

function hmc_step(θ, lj, lj_func, grad_func, H_func, ϵ, alg::HMC, momentum_sampler::Function;
                  rev_func=nothing, log_func=nothing)
  θ_new, lj_new, is_accept, τ_valid, α = _hmc_step(
            θ, lj, lj_func, grad_func, H_func, alg.tau, ϵ, momentum_sampler; rev_func=rev_func, log_func=log_func)
  return θ_new, lj_new, is_accept, α
end

# Below is a trick to remove the dependency of Stan by Requires.jl
# Please see https://github.com/TuringLang/Turing.jl/pull/459 for explanations
DEFAULT_ADAPT_CONF_TYPE = Nothing
STAN_DEFAULT_ADAPT_CONF = nothing

@static if isdefined(Turing, :CmdStan)
    DEFAULT_ADAPT_CONF_TYPE = Union{DEFAULT_ADAPT_CONF_TYPE, CmdStan.Adapt}
    STAN_DEFAULT_ADAPT_CONF = CmdStan.Adapt()
end

Sampler(alg::Hamiltonian) =  Sampler(alg, STAN_DEFAULT_ADAPT_CONF::DEFAULT_ADAPT_CONF_TYPE)
Sampler(alg::Hamiltonian, adapt_conf::DEFAULT_ADAPT_CONF_TYPE) = begin
    info=Dict{Symbol, Any}()

    # For state infomation
    info[:lf_num] = 0
    info[:eval_num] = 0

    # Adapt configuration
    info[:adapt_conf] = adapt_conf

    Sampler(alg, info)
end

function sample(model::Model, alg::Hamiltonian;
                                chunk_size=CHUNKSIZE[],             # set temporary chunk size
                                save_state=false,                   # flag for state saving
                                resume_from=nothing,                # chain to continue
                                reuse_spl_n=0,                      # flag for spl re-using
                                adapt_conf=STAN_DEFAULT_ADAPT_CONF, # adapt configuration
                )
    if ADBACKEND[] == :forward_diff
        default_chunk_size = CHUNKSIZE[]  # record global chunk size
        setchunksize(chunk_size)        # set temp chunk size
    end

    spl = reuse_spl_n > 0 ?
          resume_from.info[:spl] :
          Sampler(alg, adapt_conf)

    @assert isa(spl.alg, Hamiltonian) "[Turing] alg type mismatch; please use resume() to re-use spl"

    alg_str = isa(alg, HMC)   ? "HMC"   :
              isa(alg, HMCDA) ? "HMCDA" :
              isa(alg, SGHMC) ? "SGHMC" :
              isa(alg, SGLD)  ? "SGLD"  :
              isa(alg, NUTS)  ? "NUTS"  : "Hamiltonian"

    # Initialization
    time_total = zero(Float64)
    n = reuse_spl_n > 0 ?
        reuse_spl_n :
        alg.n_iters
    samples = Array{Sample}(undef, n)
    weight = 1 / n
    for i = 1:n
        samples[i] = Sample(weight, Dict{Symbol, Any}())
    end

    vi = if resume_from == nothing
        vi_ = VarInfo()
        model(vi_, HamiltonianRobustInit())
        spl.info[:eval_num] += 1
        vi_
    else
        deepcopy(resume_from.info[:vi])
    end

    if spl.alg.gid == 0
        link!(vi, spl)
        runmodel!(model, vi, spl)
    end

    # HMC steps
    total_lf_num = 0
    total_eval_num = 0
    accept_his = Bool[]
    PROGRESS[] && (spl.info[:progress] = ProgressMeter.Progress(n, 1, "[$alg_str] Sampling...", 0))
    for i = 1:n
        @debug "$alg_str stepping..."

        time_elapsed = @elapsed vi, is_accept = step(model, spl, vi, Val(i == 1))
        time_total += time_elapsed

        if is_accept # accepted => store the new predcits
            samples[i].value = Sample(vi, spl).value
        else         # rejected => store the previous predcits
            samples[i] = samples[i - 1]
        end
        samples[i].value[:elapsed] = time_elapsed
        if haskey(spl.info, :wum)
          samples[i].value[:lf_eps] = getss(spl.info[:wum])
        end

        total_lf_num += spl.info[:lf_num]
        total_eval_num += spl.info[:eval_num]
        push!(accept_his, is_accept)
        PROGRESS[] && ProgressMeter.next!(spl.info[:progress])
    end

    println("[$alg_str] Finished with")
    println("  Running time        = $time_total;")
    if ~isa(alg, NUTS)  # accept rate for NUTS is meaningless - so no printing
        accept_rate = sum(accept_his) / n  # calculate the accept rate
        println("  Accept rate         = $accept_rate;")
    end
    println("  #lf / sample        = $(total_lf_num / n);")
    println("  #evals / sample     = $(total_eval_num / n);")
    if haskey(spl.info, :wum)
      std_str = string(spl.info[:wum].pc)
      std_str = length(std_str) >= 32 ? std_str[1:30]*"..." : std_str   # only show part of pre-cond
      println("  pre-cond. metric    = $(std_str).")
    end

    if ADBACKEND[] == :forward_diff
        setchunksize(default_chunk_size)      # revert global chunk size
    end

    if resume_from != nothing   # concat samples
        pushfirst!(samples, resume_from.value2...)
    end
    c = Chain(0, samples)       # wrap the result by Chain
    if save_state               # save state
        # Convert vi back to X if vi is required to be saved
        if spl.alg.gid == 0 invlink!(vi, spl) end
        save!(c, spl, model, vi)
    end
    return c
end

function step(model, spl::Sampler{<:StaticHamiltonian}, vi::VarInfo, is_first::Val{true})
    spl.info[:wum] = NaiveCompAdapter(UnitPreConditioner(), FixedStepSize(spl.alg.epsilon))
    return vi, true
end

function step(model, spl::Sampler{<:AdaptiveHamiltonian}, vi::VarInfo, is_first::Val{true})
    spl.alg.gid != 0 && link!(vi, spl)
    epsilon = find_good_eps(model, spl, vi) # heuristically find good initial epsilon
    dim = length(vi[spl])
    spl.info[:wum] = ThreePhaseAdapter(spl, epsilon, dim)
    spl.alg.gid != 0 && invlink!(vi, spl)
    return vi, true
end

function step(model, spl::Sampler{<:Hamiltonian}, vi::VarInfo, is_first::Val{false})
    # Get step size
    ϵ = getss(spl.info[:wum])
    @debug "current ϵ: $ϵ"

    # Reset current counters
    spl.info[:lf_num] = 0
    spl.info[:eval_num] = 0

    @debug "X-> R..."
    if spl.alg.gid != 0
        link!(vi, spl)
        runmodel!(model, vi, spl)
    end

    grad_func = gen_grad_func(vi, spl, model)
    lj_func = gen_lj_func(vi, spl, model)
    rev_func = gen_rev_func(vi, spl)
    log_func = gen_log_func(spl)
    momentum_sampler = gen_momentum_sampler(vi, spl, spl.info[:wum].pc)
    H_func = gen_H_func(spl.info[:wum].pc)

    θ, lj = vi[spl], vi.logp

    θ_new, lj_new, is_accept, α = hmc_step(θ, lj, lj_func, grad_func, H_func, ϵ, spl.alg, momentum_sampler;
                                           rev_func=rev_func, log_func=log_func)

    @debug "decide whether to accept..."
    if is_accept
        vi[spl] = θ_new
        setlogp!(vi, lj_new)
    else
        vi[spl] = θ
        setlogp!(vi, lj)
    end

    if PROGRESS[] && spl.alg.gid == 0
        std_str = string(spl.info[:wum].pc)
        std_str = length(std_str) >= 32 ? std_str[1:30]*"..." : std_str
        haskey(spl.info, :progress) && ProgressMeter.update!(
            spl.info[:progress],
            spl.info[:progress].counter;
            showvalues = [(:ϵ, ϵ), (:α, α), (:pre_cond, std_str)],
        )
    end

    if spl.alg isa AdaptiveHamiltonian
        adapt!(spl.info[:wum], α, vi[spl], adapt_M=false, adapt_ϵ=true)
    end

    @debug "R -> X..."
    spl.alg.gid != 0 && invlink!(vi, spl)

    return vi, is_accept
end

function assume(spl::Sampler{<:Hamiltonian}, dist::Distribution, vn::VarName, vi::VarInfo)
    @debug "assuming..."
    updategid!(vi, vn, spl)
    r = vi[vn]
    # acclogp!(vi, logpdf_with_trans(dist, r, istrans(vi, vn)))
    # r
    @debug "dist = $dist"
    @debug "vn = $vn"
    @debug "r = $r" "typeof(r)=$(typeof(r))"
    r, logpdf_with_trans(dist, r, istrans(vi, vn))
end

function assume(spl::Sampler{<:Hamiltonian}, dists::Vector{<:Distribution}, vn::VarName, var::Any, vi::VarInfo)
    @assert length(dists) == 1 "[observe] Turing only support vectorizing i.i.d distribution"
    dist = dists[1]
    n = size(var)[end]

    vns = map(i -> copybyindex(vn, "[$i]"), 1:n)

    rs = vi[vns]  # NOTE: inside Turing the Julia conversion should be sticked to

    # acclogp!(vi, sum(logpdf_with_trans(dist, rs, istrans(vi, vns[1]))))

    if isa(dist, UnivariateDistribution) || isa(dist, MatrixDistribution)
        @assert size(var) == size(rs) "Turing.assume variable and random number dimension unmatched"
        var = rs
    elseif isa(dist, MultivariateDistribution)
        if isa(var, Vector)
            @assert length(var) == size(rs)[2] "Turing.assume variable and random number dimension unmatched"
            for i = 1:n
                var[i] = rs[:,i]
            end
        elseif isa(var, Matrix)
            @assert size(var) == size(rs) "Turing.assume variable and random number dimension unmatched"
            var = rs
        else
            error("[Turing] unsupported variable container")
        end
    end

    var, sum(logpdf_with_trans(dist, rs, istrans(vi, vns[1])))
end

observe(spl::Sampler{<:Hamiltonian}, d::Distribution, value::Any, vi::VarInfo) =
    observe(nothing, d, value, vi)

observe(spl::Sampler{<:Hamiltonian}, ds::Vector{<:Distribution}, value::Any, vi::VarInfo) =
    observe(nothing, ds, value, vi)


text: y = 0; for x = 1:10 println(y); y = 0 end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
struct GNNExplainer <: GeometricInterpretor

end

function explain_node(method :: GNNInterpretation{GNNExplainer}, node :: Int) :: GNNInterpretationResults

end


text:
# This file was generated, do not modify it. # hide
fit!(knn, rows=train)
ŷ = predict(knn, X[test, :]) # or use rows=test
rms(ŷ, y[test])


text:
struct KroneckerMap{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T}
    maps::As
    function KroneckerMap{T, As}(maps::As) where {T, As}
        for A in maps
            promote_type(T, eltype(A)) == T || throw(InexactError())
        end
        return new{T,As}(maps)
    end
end

KroneckerMap{T}(maps::As) where {T, As<:Tuple{Vararg{LinearMap}}} = KroneckerMap{T, As}(maps)

"""
    kron(A::LinearMap, B::LinearMap)
    kron(A, B, Cs...)

Construct a `KroneckerMap <: LinearMap` object, a (lazy) representation of the
Kronecker product of two `LinearMap`s. One of the two factors can be an `AbstractMatrix`,
which is then promoted to a `LinearMap` automatically.

To avoid fallback to the generic [`Base.kron`](@ref) in the multi-map case,
there must be a `LinearMap` object among the first 8 arguments in usage like
`kron(A, B, Cs...)`.

For convenience, one can also use `A ⊗ B` or `⊗(A, B, Cs...)` (typed as `\\otimes+TAB`) to construct the
`KroneckerMap`, even when all arguments are of `AbstractMatrix` type.

If `A`, `B`, `C` and `D` are linear maps of such size that one can form the matrix
products `A*C` and `B*D`, then the mixed-product property `(A⊗B)*(C⊗D) = (A*C)⊗(B*D)`
holds. Upon vector multiplication, this rule is checked for applicability.

# Examples
```jldoctest; setup=(using LinearAlgebra, SparseArrays, LinearMaps)
julia> J = LinearMap(I, 2) # 2×2 identity map
LinearMaps.UniformScalingMap{Bool}(true, 2)

julia> E = spdiagm(-1 => trues(1)); D = E + E' - 2I;

julia> Δ = kron(D, J) + kron(J, D); # discrete 2D-Laplace operator

julia> Matrix(Δ)
4×4 Array{Int64,2}:
 -4   1   1   0
  1  -4   0   1
  1   0  -4   1
  0   1   1  -4
```
"""
Base.kron(A::LinearMap, B::LinearMap) = KroneckerMap{promote_type(eltype(A), eltype(B))}((A, B))
Base.kron(A::LinearMap, B::KroneckerMap) = KroneckerMap{promote_type(eltype(A), eltype(B))}(tuple(A, B.maps...))
Base.kron(A::KroneckerMap, B::LinearMap) = KroneckerMap{promote_type(eltype(A), eltype(B))}(tuple(A.maps..., B))
Base.kron(A::KroneckerMap, B::KroneckerMap) = KroneckerMap{promote_type(eltype(A), eltype(B))}(tuple(A.maps..., B.maps...))
Base.kron(A::LinearMap, B::LinearMap, Cs::LinearMap...) = KroneckerMap{promote_type(eltype(A), eltype(B), map(eltype, Cs)...)}(tuple(A, B, Cs...))
Base.kron(A::AbstractMatrix, B::LinearMap) = kron(LinearMap(A), B)
Base.kron(A::LinearMap, B::AbstractMatrix) = kron(A, LinearMap(B))
# promote AbstractMatrix arguments to LinearMaps, then take LinearMap-Kronecker product
for k in 3:8 # is 8 sufficient?
    Is = ntuple(n->:($(Symbol(:A,n))::AbstractMatrix), Val(k-1))
    # yields (:A1, :A2, :A3, ..., :A(k-1))
    L = :($(Symbol(:A,k))::LinearMap)
    # yields :Ak::LinearMap
    mapargs = ntuple(n -> :(LinearMap($(Symbol(:A,n)))), Val(k-1))
    # yields (:LinearMap(A1), :LinearMap(A2), ..., :LinearMap(A(k-1)))

    @eval Base.kron($(Is...), $L, As::MapOrMatrix...) =
        kron($(mapargs...), $(Symbol(:A,k)), convert_to_lmaps(As...)...)
end

struct KronPower{p}
    function KronPower(p::Integer)
        p > 1 || throw(ArgumentError("the Kronecker power is only defined for exponents larger than 1, got $k"))
        return new{p}()
    end
end

"""
    ⊗(k::Integer)

Construct a lazy representation of the `k`-th Kronecker power
`A^⊗(k) = A ⊗ A ⊗ ... ⊗ A`, where `A` can be an `AbstractMatrix` or a `LinearMap`.
"""
⊗(k::Integer) = KronPower(k)

⊗(A, B, Cs...) = KroneckerMap{promote_type(eltype(A), eltype(B), map(eltype, Cs)...)}(convert_to_lmaps(A, B, Cs...))

Base.:(^)(A::MapOrMatrix, ::KronPower{p}) where {p} =
    ⊗(ntuple(n -> convert_to_lmaps_(A), Val(p))...)

Base.size(A::KroneckerMap) = map(*, size.(A.maps)...)

LinearAlgebra.issymmetric(A::KroneckerMap) = all(issymmetric, A.maps)
LinearAlgebra.ishermitian(A::KroneckerMap{<:Real}) = issymmetric(A)
LinearAlgebra.ishermitian(A::KroneckerMap) = all(ishermitian, A.maps)

LinearAlgebra.adjoint(A::KroneckerMap{T}) where {T} = KroneckerMap{T}(map(adjoint, A.maps))
LinearAlgebra.transpose(A::KroneckerMap{T}) where {T} = KroneckerMap{T}(map(transpose, A.maps))

Base.:(==)(A::KroneckerMap, B::KroneckerMap) = (eltype(A) == eltype(B) && A.maps == B.maps)

#################
# multiplication helper functions
#################

@inline function _kronmul!(y, B, X, At, T)
    na, ma = size(At)
    mb, nb = size(B)
    v = zeros(T, ma)
    temp1 = similar(y, na)
    temp2 = similar(y, nb)
    @views @inbounds for i in 1:ma
        v[i] = one(T)
        A_mul_B!(temp1, At, v)
        A_mul_B!(temp2, X, temp1)
        A_mul_B!(y[((i-1)*mb+1):i*mb], B, temp2)
        v[i] = zero(T)
    end
    return y
end
@inline function _kronmul!(y, B::Union{MatrixMap,UniformScalingMap}, X, At::Union{MatrixMap,UniformScalingMap}, T)
    na, ma = size(At)
    mb, nb = size(B)
    if nb*ma < mb*na
        mul!(reshape(y, (mb, ma)), B, convert(Matrix, X*At))
    else
        mul!(reshape(y, (mb, ma)), convert(Matrix, B*X), At isa MatrixMap ? At.lmap : At.λ)
    end
    return y
end

#################
# multiplication with vectors
#################

Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T,<:NTuple{2,LinearMap}}, x::AbstractVector) where {T}
    require_one_based_indexing(y)
    @boundscheck check_dim_mul(y, L, x)
    A, B = L.maps
    X = LinearMap(reshape(x, (size(B, 2), size(A, 2))); issymmetric=false, ishermitian=false, isposdef=false)
    _kronmul!(y, B, X, transpose(A), T)
    return y
end
Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerMap{T}, x::AbstractVector) where {T}
    require_one_based_indexing(y)
    @boundscheck check_dim_mul(y, L, x)
    A = first(L.maps)
    B = kron(Base.tail(L.maps)...)
    X = LinearMap(reshape(x, (size(B, 2), size(A, 2))); issymmetric=false, ishermitian=false, isposdef=false)
    _kronmul!(y, B, X, transpose(A), T)
    return y
end
# mixed-product rule, prefer the right if possible
# (A₁ ⊗ A₂ ⊗ ... ⊗ Aᵣ) * (B₁ ⊗ B₂ ⊗ ... ⊗ Bᵣ) = (A₁B₁) ⊗ (A₂B₂) ⊗ ... ⊗ (AᵣBᵣ)
Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{<:Any,<:Tuple{KroneckerMap,KroneckerMap}}, x::AbstractVector)
    B, A = L.maps
    if length(A.maps) == length(B.maps) && all(M -> check_dim_mul(M[1], M[2]), zip(A.maps, B.maps))
        A_mul_B!(y, kron(map(*, A.maps, B.maps)...), x)
    else
        A_mul_B!(y, LinearMap(A)*B, x)
    end
end
# mixed-product rule, prefer the right if possible
# (A₁ ⊗ B₁)*(A₂⊗B₂)*...*(Aᵣ⊗Bᵣ) = (A₁*A₂*...*Aᵣ) ⊗ (B₁*B₂*...*Bᵣ)
Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::CompositeMap{T,<:Tuple{Vararg{KroneckerMap{<:Any,<:Tuple{LinearMap,LinearMap}}}}}, x::AbstractVector) where {T}
    As = map(AB -> AB.maps[1], L.maps)
    Bs = map(AB -> AB.maps[2], L.maps)
    As1, As2 = Base.front(As), Base.tail(As)
    Bs1, Bs2 = Base.front(Bs), Base.tail(Bs)
    apply = all(A -> check_dim_mul(A...), zip(As1, As2)) && all(A -> check_dim_mul(A...), zip(Bs1, Bs2))
    if apply
        A_mul_B!(y, kron(prod(As), prod(Bs)), x)
    else
        A_mul_B!(y, CompositeMap{T}(map(LinearMap, L.maps)), x)
    end
end

Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) =
    A_mul_B!(y, transpose(A), x)

Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::KroneckerMap, x::AbstractVector) =
    A_mul_B!(y, adjoint(A), x)

###############
# KroneckerSumMap
###############
struct KroneckerSumMap{T, As<:Tuple{LinearMap,LinearMap}} <: LinearMap{T}
    maps::As
    function KroneckerSumMap{T, As}(maps::As) where {T, As}
        for A in maps
            size(A, 1) == size(A, 2) || throw(ArgumentError("operators need to be square in Kronecker sums"))
            promote_type(T, eltype(A)) == T || throw(InexactError())
        end
        return new{T,As}(maps)
    end
end

KroneckerSumMap{T}(maps::As) where {T, As<:Tuple{LinearMap,LinearMap}} = KroneckerSumMap{T, As}(maps)

"""
    kronsum(A, B)
    kronsum(A, B, Cs...)

Construct a `KroneckerSumMap <: LinearMap` object, a (lazy) representation of the
Kronecker sum `A⊕B = A ⊗ Ib + Ia ⊗ B` of two square linear maps of type
`LinearMap` or `AbstractMatrix`. Here, `Ia` and `Ib` are identity operators of
the size of `A` and `B`, respectively. Arguments of type `AbstractMatrix` are
automatically promoted to `LinearMap`s.

For convenience, one can also use `A ⊕ B` or `⊕(A, B, Cs...)` (typed as
`\\oplus+TAB`) to construct the `KroneckerSumMap`.

# Examples
```jldoctest; setup=(using LinearAlgebra, SparseArrays, LinearMaps)
julia> J = LinearMap(I, 2) # 2×2 identity map
LinearMaps.UniformScalingMap{Bool}(true, 2)

julia> E = spdiagm(-1 => trues(1)); D = LinearMap(E + E' - 2I);

julia> Δ₁ = kron(D, J) + kron(J, D); # discrete 2D-Laplace operator, Kronecker sum

julia> Δ₂ = kronsum(D, D);

julia> Δ₃ = D^⊕(2);

julia> Matrix(Δ₁) == Matrix(Δ₂) == Matrix(Δ₃)
true
```
"""
kronsum(A::MapOrMatrix, B::MapOrMatrix) =
    KroneckerSumMap{promote_type(eltype(A), eltype(B))}(convert_to_lmaps(A, B))
kronsum(A::MapOrMatrix, B::MapOrMatrix, C::MapOrMatrix, Ds::MapOrMatrix...) =
    kronsum(A, kronsum(B, C, Ds...))

struct KronSumPower{p}
    function KronSumPower(p::Integer)
        p > 1 || throw(ArgumentError("the Kronecker sum power is only defined for exponents larger than 1, got $k"))
        return new{p}()
    end
end

"""
    ⊕(k::Integer)

Construct a lazy representation of the `k`-th Kronecker sum power `A^⊕(k) = A ⊕ A ⊕ ... ⊕ A`,
where `A` can be a square `AbstractMatrix` or a `LinearMap`.
"""
⊕(k::Integer) = KronSumPower(k)

⊕(a, b, c...) = kronsum(a, b, c...)

Base.:(^)(A::MapOrMatrix, ::KronSumPower{p}) where {p} = kronsum(ntuple(n -> convert_to_lmaps_(A), Val(p))...)

Base.size(A::KroneckerSumMap, i) = prod(size.(A.maps, i))
Base.size(A::KroneckerSumMap) = (size(A, 1), size(A, 2))

LinearAlgebra.issymmetric(A::KroneckerSumMap) = all(issymmetric, A.maps)
LinearAlgebra.ishermitian(A::KroneckerSumMap{<:Real}) = all(issymmetric, A.maps)
LinearAlgebra.ishermitian(A::KroneckerSumMap) = all(ishermitian, A.maps)

LinearAlgebra.adjoint(A::KroneckerSumMap{T}) where {T} = KroneckerSumMap{T}(map(adjoint, A.maps))
LinearAlgebra.transpose(A::KroneckerSumMap{T}) where {T} = KroneckerSumMap{T}(map(transpose, A.maps))

Base.:(==)(A::KroneckerSumMap, B::KroneckerSumMap) = (eltype(A) == eltype(B) && A.maps == B.maps)

Base.@propagate_inbounds function A_mul_B!(y::AbstractVector, L::KroneckerSumMap, x::AbstractVector)
    @boundscheck check_dim_mul(y, L, x)
    A, B = L.maps
    ma, na = size(A)
    mb, nb = size(B)
    X = reshape(x, (nb, na))
    Y = reshape(y, (nb, na))
    mul!(Y, X, convert(AbstractMatrix, transpose(A)))
    mul!(Y, B, X, true, true)
    return y
end

Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) =
    A_mul_B!(y, transpose(A), x)

Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::KroneckerSumMap, x::AbstractVector) =
    A_mul_B!(y, adjoint(A), x)


text:
export SplitExplicitSolverType

"""
# Description
    SplitExplicitSolverType

This solver type constructs an ODE solver using the SplitExplicitLSRK2nSolver.
# Arguments
- `dt_slow` (AbstractFloat): Time step for the slow solver
- `dt_fast` (AbstractFloat): Time step for the fast solver
- `slow_method` (Function): Function defining the explicit
    Runge-Kutta solver for the slow model.
    Default: `LSRK54CarpenterKennedy`
- `fast_method` (Function): Function defining the explicit
    Runge-Kutta solver for the fast model.
    Default: `LSRK54CarpenterKennedy`
"""
struct SplitExplicitSolverType{SM, FM, FT} <: AbstractSolverType
    # Function for an explicit Runge-Kutta method
    slow_method::SM
    # Function for an explicit Runge-Kutta method
    fast_method::FM
    # time step for slow method
    dt_slow::FT
    # time step for fast method
    dt_fast::FT

    function SplitExplicitSolverType{FT}(
        dt_slow,
        dt_fast;
        slow_method = LSRK54CarpenterKennedy,
        fast_method = LSRK54CarpenterKennedy,
    ) where {FT <: AbstractFloat}
        SM = typeof(slow_method)
        FM = typeof(fast_method)

        return new{SM, FM, FT}(slow_method, fast_method, dt_slow, dt_fast)
    end
end

"""
    getdtmodel(ode_solver::AbstractSolverType, bl)

A function which returns a model representing the dynamics
with the most restrictive time-stepping requirements.
"""
function getdtmodel(::SplitExplicitSolverType, bl)
    # For explicit methods, the entire model itself
    # contributes to the total stability of the time-integrator
    return bl
end

"""
# Description
    function solversetup(
        ode_solver::ExplicitSolverType,
        dg,
        Q,
        dt,
        t0,
        diffusion_direction,
    )

Creates an explicit ODE solver.
"""
function solversetup(ode_solver::SplitExplicitSolverType, dg_3D, Q_3D, _, t0, _)
    dg_2D = dg_3D.modeldata.dg_2D
    Q_2D = dg_3D.modeldata.Q_2D

    fast_solver =
        ode_solver.fast_method(dg_2D, Q_2D, dt = ode_solver.dt_fast, t0 = t0)
    slow_solver =
        ode_solver.slow_method(dg_3D, Q_3D, dt = ode_solver.dt_slow, t0 = t0)

    solver = SplitExplicitLSRK2nSolver(slow_solver, fast_solver)

    return solver
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
using OpenBCI
using Base.Test
# Run tests

@time include("test1.jl")


text:
function reserve_slacks(psi_container::PSIContainer, reserve_name::String)
    time_steps = model_time_steps(psi_container)
    var_name = make_variable_name(reserve_name, SLACK_UP)
    variable = add_var_container!(psi_container, var_name, time_steps)

    for jx in time_steps
        variable[jx] = JuMP.@variable(
            psi_container.JuMPmodel,
            base_name = "$(var_name)_{$(jx)}",
            lower_bound = 0.0
        )
        JuMP.add_to_expression!(
            psi_container.cost_function,
            variable[jx] * SERVICES_SLACK_COST,
        )
    end
    return variable
end


text:
"""
```
estimate_bma(m, data, prior; return_output = false, filestring_addl = [])
```

Estimate a Bayesian Model Average.

### Arguments:
- `m::PoolModel`: PoolModel object

### Optional Arguments:
- `data`: well-formed data as `Matrix` or `DataFrame`. If this is not provided, the `load_data` routine will be executed.
- `prior::Float64`: prior probability between the two models

### Keyword Arguments:
- `return_output::Bool`: option to return output. If false, `nothing` is returned.
- `filestring_addl::Vector{String}`: Additional strings to append to output files.
"""
function estimate_bma(m::PoolModel, df::DataFrame, prior::Float64 = 0.5;
                      return_output::Bool = false, save_output::Bool = true,
                      filestring_addl::Vector{String} = Vector{String}(undef,0))
    return estimate_bma(m, df_to_matrix(m, df), prior;
                        return_output = return_output, save_output = save_output,
                        filestring_addl = filestring_addl)
end

function estimate_bma(m::PoolModel, data::Matrix{Float64} = Matrix{Float64}(0,0),
                      prior::Float64 = 0.5; return_output::Bool = false, save_output::Bool = true,
                      filestring_addl::Vector{String} = Vector{String}(undef,0))

    # Compute λ weights
    λ = zeros(size(data,2))
    λ[1] = prior * data[1,1] / (prior * data[1,1] + (1 - prior) * data[2,1])
    for t in 2:size(data,2)
        λ[t] = λ[t-1] * data[1,t] / (λ[t-1] * data[1,t] + (1 - λ[t-1]) * data[2,t])
    end

    # Save weights
    if save_output
        h5open(rawpath(m, "estimate", "bmasave.h5", filestring_addl), "w") do file
            write(file, "bmaparams", λ)
        end
    end

    if return_output
        prob_vec = λ .* vec(data[1,:]) + (1 .- λ) .* vec(data[2,:])
        return λ, prob_vec
    end
end


text:
immutable QuasiHyperbolic{T} <: Model
    a::T
    β::T
    δ::T
end

Base.show(io::IO, ::Type{QuasiHyperbolic}) = print(io, "Classical Quasi-Hyperbolic")

function Base.show(io::IO, m::QuasiHyperbolic)
    @printf(
        io,
        "Classical Quasi-Hyperbolic (a = %.3f, β = %.3f, δ = %.3f)",
        m.a,
        m.β,
        m.δ,
    )
end

function initial_parameters(::Type{QuasiHyperbolic})
    return Float64[log(1.0), logit(0.5), logit(0.5)]
end

canonical_parameters(m::QuasiHyperbolic) = (m.a, m.β, m.δ)

function QuasiHyperbolic(parameters::Vector)
    a = exp(parameters[1])
    β = invlogit(parameters[2])
    δ = invlogit(parameters[3])
    return QuasiHyperbolic(a, β, δ)
end

@inline function predict(
    m::QuasiHyperbolic,
    x2::Real,
    t2::Real,
    x1::Real,
    t1::Real,
)
    a, β, δ = m.a, m.β, m.δ
    u1 = x1 * ifelse(t1 == 0.0, 1.0, β) * δ^t1
    u2 = x2 * ifelse(t2 == 0.0, 1.0, β) * δ^t2
    z = a * (u2 - u1)
    p = invlogit(z)
    return p
end

@inline function gradient_component(
    m::QuasiHyperbolic,
    x2::Real,
    t2::Real,
    x1::Real,
    t1::Real,
)
    a, β, δ = m.a, m.β, m.δ
    βt2 = ifelse(t2 == 0.0, 1.0, β)
    βt1 = ifelse(t1 == 0.0, 1.0, β)
    Δβt2 = ifelse(t2 == 0.0, 0.0, 1.0)
    Δβt1 = ifelse(t1 == 0.0, 0.0, 1.0)
    δt2 = δ^t2
    δt1 = δ^t1
    gr = (
        (x2 * βt2 * δt2 - x1 * βt1 * δt1) * a,
        a * (x2 * Δβt2 * δt2 - x1 * Δβt1 * δt1) * β * (1 - β),
        a * ((x2 * βt2 * t2 * δt2 - x1 * βt1 * t1 * δt1) / δ) * δ * (1 - δ),
    )
    return gr
end


text:
f(x) = [
  -x[2]-x[3]
  x[1] + .1*x[2]
  .1 + x[3]*(x[1]-14)
  ]

f′(x) = [
  0     -1      -1
  1     .1      0
  x[3]   0      (x[1]-14)
  ]

m = 31
N = 2m+1
n = 3
Y = 1000*rand(Complex{Float64}, N*n)

using PyPlot
ion()

function circulant(v)
  N = length(v)
  return [ v[mod(i-j,N)+1] for i in 1:N, j in 1:N ]
end

function Jnum(Y, ϵ=1e-4)
	h₀ = H(Y)
	tmp = cell(length(Y))
	for i in 1:length(Y)
		Y[i] += ϵ
		tmp[i] = vec((H(Y) - h₀) ./ ϵ)
		Y[i] -= ϵ
	end
	return reduce(hcat, tmp)
end

function H(Y)
	Y = reshape(Y,N,n)
	# Y = [Y; conj(Y[end:-1:2,:])]
	tmp = ifft(Y, [1])
	tmp = mapslices(f, tmp, [2])
	tmp = fft(tmp, [1])
	# return vec(tmp)
	return vec(tmp - im .* ifftshift(-m:m) .* Y)
end

@noinline function J(Y)
	Y = reshape(Y,N,n)
	tmp = ifft(Y, [1])
	M = Array{Complex{Float64}}(N,n,n)
	for i in 1:N
		M[i,:,:] = f′(tmp[i,:]')
	end
	ifft!(M,[1])
	K = spdiagm(im*ifftshift(-m:m))
	tmp = [ circulant(M[:,k,l]) - (k==l ? K : 0) for k in 1:n, l in 1:n ]
	tmp = reducedim(hcat, tmp, [2], Array{Complex}(N,0))
	tmp = reducedim(vcat, tmp, [1], Array{Complex}(0,n*N))[1]
end

matshow(real(Jnum(Y) - J(Y))); colorbar()
matshow(imag(Jnum(Y) - J(Y))); colorbar()
matshow(abs(Jnum(Y)) ./ (abs(J(Y)))+1); colorbar()
matshow(real(J(Y))); colorbar()
matshow(real(Jnum(Y))); colorbar()
matshow(imag(J(Y))); colorbar()
matshow(imag(Jnum(Y))); colorbar()

close()

newton(f,f′,x) = x - inv(f′(x))*f(x)

Ytmpa = 5*rand(Complex{Float64}, N*n)
Ytmp = Ytmpa
Ytmp = newton(H,J,Ytmp); println(norm(H(Ytmp)))
Ytmp = newton(H,Jnum,Ytmp); println(norm(H(Ytmp)))
H(Ytmp)
J(Ytmp)

@noinline function J2(Y)
	Y = reshape(Y,N,n)
	tmp = ifft(Y, [1])
	M = Array{Complex{Float64}}(N,n,n)
	for i in 1:N
		M[i,:,:] = f′(tmp[i,:]')
	end
	ifft!(M,[1])
	K = spdiagm(im*ifftshift(-m:m))
	tmp = [ circulant(M[:,k,l]) - (k==l ? K : 0) for k in 1:n, l in 1:n ]
end

function realJ(Y)
	Y = reshape(Y,m+1+m,n)
	Y = Y[1:m+1,:] + im*[zeros(1,n); Y[m+2:end,:]]
	Y = [Y; conj(Y[end:-1:2,:])]
	tmp = J2(Y)
	tmp = map(x->x[1:m+1,1:m+1], tmp)
	tmp = map(x->[real(x) -imag(x)[:,2:end]; imag(x)[2:end,:] real(x)[2:end,2:end]], tmp)
	tmp = reducedim(hcat, tmp, [2], Array{Float64}(N,0))
	tmp = reducedim(vcat, tmp, [1], Array{Float64}(0,N*n))[1]
end

Y = 1000000000*rand(N*n)
tmp1 = realJ(Y)
tmp2 = Jgen(Y)
tmp3 = Jroessler(Y)

matshow(tmp2-tmp3); colorbar()
matshow(abs(tmp1)./(abs(tmp3)+.01)); colorbar()
matshow(tmp3); colorbar()

function Jroessler(V)
	local α=.1, β=.1, ℵ=14
	local X₀,Xᵣ,Xᵢ,Y₀,Yᵣ,Yᵢ,Z₀,Zᵣ,Zᵢ
	V = reshape(V,N,n)

	X₀,Y₀,Z₀ = V[1,1], V[1,2], V[1,3]
	Xᵣ,Yᵣ,Zᵣ = V[2:2+m-1,1], V[2:2+m-1,2], V[2:2+m-1,3]
	Xᵢ,Yᵢ,Zᵢ = V[2+m:end,1], V[2+m:end,2], V[2+m:end,3]

	local D = (1:m)

	# Jacobian of circular convolution of coefficients of real functions in appropriate format.
	# ∂/∂X cconv(X,Y) = ∂/∂X F(x⋅y) = rCCD(Y)
	rCCD(C) = rCCD(real(C[1]), real(C[2:end]), imag(C[2:end]))
	function rCCD(V₀,Vᵣ,Vᵢ)
		local I1,I2,Wᵣ,Wᵢ
		I1 = [ mod(i-j, 2m+1)+1 for i in 0:m, j in 0:m ]
		I2 = [ mod(i+j, 2m+1)+1 for i in 0:m, j in 0:m ]
		Wᵣ = [V₀; Vᵣ; Vᵣ[end:-1:1]]
		Wᵢ = [.0; Vᵢ; -Vᵢ[end:-1:1]]

		local rtn =  [
			(Wᵣ[I1]+Wᵣ[I2])					(-Wᵢ[I1]+Wᵢ[I2])[:,2:end]
			(Wᵢ[I1]+Wᵢ[I2])[2:end,:]		(Wᵣ[I1]-Wᵣ[I2])[2:end,2:end]
		] / (2m+1)
		rtn[:,1] /= 2.0
		return rtn
	end

	local AbyX, AbyY, AbyZ, Abyω, Abyℵ,
		BbyX, BbyY, BbyZ, Bbyω, Bbyℵ,
		CbyX, CbyY, CbyZ, Cbyω, Cbyℵ,
		LbyX, LbyY, LbyZ, Lbyω, Lbyℵ

	AbyX = [
		.0					zeros(1,m)			zeros(1,m)
		zeros(m)			zeros(m,m)			diagm(D)
		zeros(m)			diagm(-D)			zeros(m,m)
	]

	AbyY = AbyZ = -eye(2m+1)

	BbyX = eye(2m+1)

	BbyY = [
		α					zeros(1,m)			zeros(1,m)
		zeros(m)			α*eye(m)			diagm(D)
		zeros(m)			diagm(-D)			α*eye(m)
	]

	BbyZ = zeros(2m+1,2m+1)

	CbyX = rCCD(Z₀,Zᵣ,Zᵢ)

	CbyY = zeros(2m+1, 2m+1)

	CbyZ = rCCD(X₀,Xᵣ,Xᵢ) + [
		-ℵ					zeros(1,m)			zeros(1,m)
		zeros(m)			-ℵ*eye(m,m)			diagm(D)
		zeros(m)			diagm(-D)			-ℵ*eye(m,m)
	]

	LbyX = [ 1.0			ones(1,m)			zeros(1,m) ]
	LbyY = zeros(1,2m+1)
	LbyZ = zeros(1,2m+1)

	Abyω = [
		.0
		(1:m).*Xᵢ
		-(1:m).*Xᵣ
	]

	Bbyω = [
		.0
		(1:m).*Yᵢ
		-(1:m).*Yᵣ
	]

	Cbyω = [
		.0
		(1:m).*Zᵢ
		-(1:m).*Zᵣ
	]

	Lbyω = .0

	Abyℵ = Bbyℵ = zeros(2m+1)
	Cbyℵ = [-Z₀; -Zᵣ; -Zᵢ]
	Lbyℵ = .0

	return [
		AbyX AbyY AbyZ
		BbyX BbyY BbyZ
		CbyX CbyY CbyZ
	]
end

function fToH(f)
	return function H(V::Vector{Float64})
		V = reshape(V,N,n)
		anchor = sum(V[1:m+1]) #TODO hm...

		ω,ρ = V[end-1:end]
		V = reshape(V[1:end-2],2m+1,3)
		V = complex(V[1:m+1,:], [zeros(3)'; V[m+2:end,:]])

		defect = rfft(mapslices(x->f(0,[x;ρ]), irfft(V, 2m, [1]), [2]), [1]) - im*ω*(0:m) .* V
		defect = vec([ real(defect); imag(defect)[2:end,:] ])

		return [defect; anchor]
	end
end

function f′ToJ(f′)
	return function J(V::Vector{Float64})
		V = reshape(V,N,n)
		V = complex(V[1:m+1,:], [zeros(1,n);V[m+2:end,:]])
		T = irfft(V, N, [1])

		M = Array{Float64}(size(T,1), size(T,2), size(T,2)) #n, comp, deriv
		for i in 1:N
			M[i,:,:] = f′(T[i,:]')
		end

		M′ = [ rCCD(rfft(M[:,i,j])) for i in 1:size(T,2), j in 1:size(T,2) ]

		for i in 1:size(T,2)
			D = diagm(0:m)
			M′[i,i][m+2:end,1:m+1] -= D[2:end,:]
			M′[i,i][1:m+1, m+2:end] += D[:,2:end]
		end

		M′ = reducedim(vcat, M′, [1], Array{Float64}(0,size(M′[1,1],2)))
		M′ = reducedim(hcat, M′, [2], Array{Float64}(size(M′[1,1],1),0))[1]

		return M′
	end
end

Jgen = f′ToJ(f′)

# Jacobian of circular convolution of coefficients of real functions in appropriate format.
# ∂/∂X cconv(X,Y) = ∂/∂X F(x⋅y) = rCCD(Y)
rCCD(C) = rCCD(real(C[1]), real(C[2:end]), imag(C[2:end]))
function rCCD(V₀,Vᵣ,Vᵢ)
	m = length(Vᵣ)
	I1 = [ mod(i-j, 2m+1)+1 for i in 0:m, j in 0:m ]
	I2 = [ mod(i+j, 2m+1)+1 for i in 0:m, j in 0:m ]
	Wᵣ = [V₀; Vᵣ; Vᵣ[end:-1:1]]
	Wᵢ = [.0; Vᵢ; -Vᵢ[end:-1:1]]

	local rtn =  [
		(Wᵣ[I1]+Wᵣ[I2])					(-Wᵢ[I1]+Wᵢ[I2])[:,2:end]
		(Wᵢ[I1]+Wᵢ[I2])[2:end,:]		(Wᵣ[I1]-Wᵣ[I2])[2:end,2:end]
	] / (2m)
	rtn[:,1] /= 2
	return rtn
end


text:
export Interval

struct Interval{T<:Real}
    a::T
    b::T
end

Base.minimum(it::Interval) = it.a
Base.maximum(it::Interval) = it.b
Base.in(x, it::Interval) = (x <= maximum(it)) && (x >= minimum(it))
Base.isempty(it::Interval) = ~(maximum(it) > minimum(it))

function Base.intersect(A::Interval, B::Interval)
    if isempty(A) || isempty(B)
        return Interval(0.0, 0.0)
    end
    Amin = minimum(A)
    Amax = maximum(A)
    Bmin = minimum(B)
    Bmax = maximum(B)
    if (Bmin > Amax) || (Amin > Bmax)
        return Interval(0.0, 0.0)
    else
        Cmin = max(Amin, Bmin)
        Cmax = min(Amax, Bmax)
        return Interval(Cmin, Cmax)
    end
end


text:
function holtslag_psim(xi)
if xi>0
    y=1+5*xi;
else
    y=(1-16*xi)^(-1/4);
end
return y
end


text:
################# load data
using JSON
using HDF5
function load_data(f)
    nstripes = length(f)
    coeffs_1 = f[1]["legendre_coeffs"]
    num_pos = length(coeffs_1)
    num_coeffs = length(coeffs_1[1])

    log10_τ = zeros(nstripes)
    position = zeros(num_pos, nstripes)
    # optical coefficients: coefficients x positions on stripe x stripe id
    optical = zeros(num_pos, num_coeffs, nstripes)
    helper(x) = tuple(x...)
    rescaling_parameters = []
    T_max = zeros(nstripes)
    log10_τ = zeros(nstripes)
    for i in 1:nstripes
        # stripe conditions
        T_max[i] = f[i]["meta"]["temp"]
        log10_τ[i] = log10(f[i]["meta"]["dwell"])

        # centered position on the wafer
        x = float.(f[i]["pos_wafer"]) # position
        c = float(f[i]["center_real"]) # center
        @. position[:, i] = x - c

        # optical coefficients
        coeffs_i = f[i]["legendre_coeffs"]
        for j = 1:num_pos
            optical[j, :, i] = float.(coeffs_i[j])
        end

        # rescaling parameters
        push!(rescaling_parameters, helper.(f[i]["rescaling_parameters"]))
    end
    position = [x for x in eachcol(position)] # convert to vector for each stripe index
    optical = [[x for x in eachcol(optical[:, :, i])] for i in 1:size(optical, 3)]
    return position, optical, rescaling_parameters, T_max, log10_τ
end


text:
"""
distance(tree::Tree,
         node1::Int64,
         node2::Int64)

Distance between two `Node`s on a `Tree`
"""
function distance(tree::Tree,
                  node1::Int64,
                  node2::Int64)
  path = branchpath(tree, node1, node2)
  dist = 0.
  for i in path
    dist += tree.branches[i].length
  end
  return dist
end

"""
distance(tree::Tree)

Pairwise distances between all leaf `Node`s on a `Tree`
"""
function distance(tree::Tree)
  leaves = findleaves(tree)
  return [distance(tree, i, j) for i in leaves, j in leaves]
end

"""
distance(tree::Tree,
         node::Int64)

Distance between a `Node` and it's associated root
"""
function distance(tree::Tree,
                  node::Int64)
  path = branchpath(tree, node)
  dist = 0.
  for i in path
    dist += tree.branches[i].length
  end
  return dist
end


text:
#=
forecast.jl

    Forecasting routines for linear Gaussian state space models

@author: Quint Wiersma <q.wiersma@vu.nl>

@date: 2021/12/06
=#
# Struct
struct Forecast{Ty, TF, Ta, TP}
    y_hat::Ty   # forecasts
    F_hat::TF   # forecast error variance
    a_hat::Ta   # forecasted states
    P_hat::TP   # forecasted state variances
end
# Constructor
function Forecast(n::Integer, p::Integer, h::Integer, T::Type)
    # forecast components
    y_hat= Matrix{T}(undef, n, h)
    F_hat= Array{T,3}(undef, n, n, h)
    a_hat= Matrix{T}(undef, p, h)
    P_hat= Array{T,3}(undef, p, p, h)

    return Forecast(y_hat, F_hat, a_hat, P_hat)
end

"""
	forecast(sys, h)

Compute ``h``-step ahead forecasts for observations and states with
corresponding forecast error variances and forecasted state variances for a
state space model with system matrices `sys`.

#### Arguments
  - `sys::StateSpaceSystem`	: state space system matrices
  - `h::Integer`			: forecast horizon

#### Returns
  - `f::Forecast`   : forecasts and forecast error variances
"""
function forecast(sys::LinearTimeInvariant, h::Integer)
	# get dims
	(n,T_len)= size(sys.y)
    p= length(sys.a1)

    # Type
    T= eltype(sys.y)

    # Kalman filter
    filter= MultivariateFilter(n, p, T_len, T)
    kalman_filter!(filter, sys)

    # forecasts
    f= Forecast(n, p, h, T)
    @inbounds @fastmath for i in 1:h
        # Store filter forecast output
        f.F_hat[:,:,i]= view(filter.F,:,:,T_len-h+i)
        f.a_hat[:,i]= view(filter.a,:,T_len-h+i)
        f.P_hat[:,:,i]= view(filter.P,:,:,T_len-h+i)

        # forecast
        y_hat= view(f.y_hat, :, i)
        a_hat= view(f.a_hat, :, i)
        y_hat.= sys.d
        mul!(y_hat, sys.Z, a_hat, 1., 1.)
    end

	return f
end

function forecast(sys::LinearTimeVariant, h::Integer)
	# get dims
	(n,T_len)= size(sys.y)
    p= length(sys.a1)

    # Type
    T= eltype(sys.y)

    # Kalman filter
    filter= MultivariateFilter(n, p, T_len, T)
    kalman_filter!(filter, sys)

    # forecasts
    f= Forecast(n, p, h, T)
    @inbounds @fastmath for i in 1:h
        # Store filter forecast output
        f.F_hat[:,:,i]= view(filter.F,:,:,T_len-h+i)
        f.a_hat[:,i]= view(filter.a,:,T_len-h+i)
        f.P_hat[:,:,i]= view(filter.P,:,:,T_len-h+i)

        # forecast
        y_hat= view(f.y_hat, :, i)
        a_hat= view(f.a_hat, :, i)
        y_hat.= view(sys.d, :, T_len-h+i)
        mul!(y_hat, sys.Z[T_len-h+i], a_hat, 1., 1.)
    end

	return f
end


text:
import Base.+

function get_matrices(
  ar::AbstractArray{T, 1},
  cr::AbstractArray{T, 1},
  ac::AbstractArray{T, 1},
  bc::AbstractArray{T, 1},
  cc::AbstractArray{T, 1},
  dc::AbstractArray{T, 1},
  x::AbstractArray{T, 1},
  diag::AbstractArray{T, 1}
) where T
  N = size(x, 1)
  Jr = size(ar, 1)
  Jc = size(ac, 1)

  a = diag .+ (sum(ar) + sum(ac))
  dx = x[2:N] - x[1:N-1]

  arg = x * transpose(dc)
  ca = cos.(arg)
  sa = sin.(arg)

  U = cat(
    ones(N) * transpose(ar),
    ca .* transpose(ac) + sa .* transpose(bc),
    sa .* transpose(ac) - ca .* transpose(bc),
    dims=2
  )
  V = cat(ones(N, Jr), ca, sa, dims=2)
  Pc = exp.(-dx * transpose(cc))
  P = cat(exp.(-dx * transpose(cr)), Pc, Pc, dims=2)

  return a, U, V, P
end

abstract type Term end

function get_terms(term::Term)
  return tuple(term)
end

function get_coefficients(term::Term)
  return zeros(0), zeros(0), zeros(0), zeros(0), zeros(0), zeros(0)
end

function get_matrices(term::Term, x, diag)
  return get_matrices(get_coefficients(term)..., x, diag)
end

struct TermSum <: Term
  terms
end

function get_terms(term::TermSum)
  return term.terms
end

function get_coefficients(term_sum::TermSum)
  ar = zeros(0)
  cr = zeros(0)
  ac = zeros(0)
  bc = zeros(0)
  cc = zeros(0)
  dc = zeros(0)
  for term in term_sum.terms
      coeffs = get_coefficients(term)
      ar = vcat(ar, coeffs[1])
      cr = vcat(cr, coeffs[2])
      ac = vcat(ac, coeffs[3])
      bc = vcat(bc, coeffs[4])
      cc = vcat(cc, coeffs[5])
      dc = vcat(dc, coeffs[6])
  end
  return ar, cr, ac, bc, cc, dc
end

function +(t1::Term, t2::Term)
  return TermSum((get_terms(t1)..., get_terms(t2)...))
end

struct RealTerm <: Term
  a::DenseArray{Float64, 1}
  c::DenseArray{Float64, 1}
  RealTerm(a::Float64, c::Float64) = new(a .+ zeros(1), c .+ zeros(1))
end

function get_coefficients(term::RealTerm)
  return term.a, term.c, zeros(0), zeros(0), zeros(0), zeros(0)
end

struct ComplexTerm <: Term
  a::DenseArray{Float64, 1}
  b::DenseArray{Float64, 1}
  c::DenseArray{Float64, 1}
  d::DenseArray{Float64, 1}
  ComplexTerm(a::Float64, b::Float64, c::Float64, d::Float64) = new(a .+ zeros(1), b .+ zeros(1), c .+ zeros(1), d .+ zeros(1))
end

function get_coefficients(term::ComplexTerm)
  return zeros(0), zeros(0), term.a, term.b, term.c, term.d
end

struct SHOTerm <: Term
  S0::Float64
  w0::Float64
  Q::Float64
end

function get_coefficients(term::SHOTerm)
  eps = 1e-5
  if term.Q < 0.5
    f = sqrt(max(1 - 4 * term.Q^2, eps))
    a = 0.5 * term.S0 * term.w0 * term.Q
    c = 0.5 * term.w0 / term.Q
    return (
      [a * (1 + 1 / f), a * (1 - 1 / f)],
      [c * (1 - f), c * (1 + f)],
      zeros(0), zeros(0), zeros(0), zeros(0)
    )
  end
  f = sqrt(max(4 * term.Q^2 - 1, eps))
  a = term.S0 * term.w0 * term.Q
  c = 0.5 * term.w0 / term.Q
  return zeros(0), zeros(0), [a], [a/f], [c], [c*f]
end


text:
"""
`module  JAC.GreenFunction`
    ... a submodel of JAC that contains all methods for computing approximate many-electron Green functions.
"""
module GreenFunction

    using Printf, ..AngularMomentum, ..Basics, ..Defaults, ..ManyElectron, ..Nuclear, ..Radial, ..TableStrings

    """
    `abstract type GreenFunction.AbstractGreenApproach`
        ... defines an abstract and a number of singleton types for approximating a many-electron Green
            function for calculating second-order processes.

      + struct SingleCSFwithoutCI
        ... to approximate the many-electron multiplets (gMultiplet) for every chosen level symmetry by dealing with each CSF
            independently, and without any configuration interaction. This is a fast but also very rough approximation.

      + struct CoreSpaceCI
        ... to approximate the many-electron multiplets (gMultiplet) by taking the electron-electron interaction between the
            bound-state orbitals into account.

      + struct DampedSpaceCI
        ... to approximate the many-electron multiplets (gMultiplet) by taking the electron-electron interaction for all, the
            bound and free-electron, orbitals into account but by including a damping factor e^{tau*r}, tau > 0 into the
            electron densities rho_ab (r) --> rho_ab (r) * e^{tau*r}
    """
    abstract type  AbstractGreenApproach                           end
    struct         SingleCSFwithoutCI  <:  AbstractGreenApproach   end
    struct         CoreSpaceCI         <:  AbstractGreenApproach   end
    struct         DampedSpaceCI       <:  AbstractGreenApproach   end

    """
    `struct  GreenFunction.Settings`
        ... defines a type for defining the details and parameters of the approximated Green functions.

        + approach                 ::AbstractGreenApproach            ... Approach used to approximate the representation.
        + excitationScheme         ::Basics.AbstractExcitationScheme  ... Applied excitation scheme with regard to the
                                                                          given bound configuration.
        + nMax                     ::Int64            ... maximum principal quantum numbers of (single-electron)
                                                          excitations that are to be included into the representation.
        + lValues                  ::Array{Int64,1}   ... List of (non-relativistic) orbital angular momenta for which
                                                          (single-electron) excitations are to be included.
        + levelSymmetries          ::Array{LevelSymmetry,1}           ... Total level symmetries J^P to be included into
                                                                          the  Green function representation.
        + printBeforeComputation   ::Bool             ... True if a short overview is to be printed before.
        + selectLevels             ::Bool             ... True if individual levels are selected for the computation.
        + selectedLevels           ::Array{Int64,1}   ... List of selected levels.
    """
    struct Settings
        approach                   ::AbstractGreenApproach
        excitationScheme           ::Basics.AbstractExcitationScheme
        nMax                       ::Int64
        lValues                    ::Array{Int64,1}
        levelSymmetries            ::Array{LevelSymmetry,1}
        printBeforeComputation     ::Bool
        selectLevels               ::Bool
        selectedLevels             ::Array{Int64,1}
    end

    """
    `GreenFunction.Settings()`  ... constructor for an `empty` instance of GreenFunction.
    """
    function Settings()
        Settings( GreenFunction.SingleCSFwithoutCI(), Basics.NoExcitationScheme(), 0, Int64[], LevelSymmetry[], false, false, Int64[])
    end

    # `Base.show(io::IO, settings::GreenFunction.Settings)`  ... prepares a proper printout of the variable settings::GreenFunction.Settings.
    function Base.show(io::IO, settings::GreenFunction.Settings)
        println(io, "approach:                 $(settings.approach)  ")
        println(io, "excitationScheme:         $(settings.excitationScheme)  ")
        println(io, "nMax:                     $(settings.nMax)  ")
        println(io, "lValues:                  $(settings.lValues)  ")
        println(io, "levelSymmetries:          $(settings.levelSymmetries)  ")
        println(io, "printBeforeComputation:   $(settings.printBeforeComputation)  ")
        println(io, "selectLevels:             $(settings.selectLevels)  ")
        println(io, "selectedLevels:           $(settings.selectedLevels)  ")
    end

    """
    `struct  GreenFunction.Channel`
        ... defines a type for a single symmetry channel of an (approximate) Green function.

        + symmetry          ::LevelSymmetry    ... Level symmetry of this part of the representation.
        + gMultiplet        ::Multiplet        ... Multiplet of (scattering) levels of this symmetry.
    """
    struct Channel
        symmetry            ::LevelSymmetry
        gMultiplet          ::Multiplet
    end

    """
    `GreenFunction.Channel()`  ... constructor for an `empty` instance of GreenFunction.Channel.
    """
    function Channel()
        Channel( LevelSymmetry(0, Basics.plus), ManyElectron.Multiplet)
    end

    # `Base.show(io::IO, channel::GreenFunction.Channel)`  ... prepares a proper printout of the variable channel::GreenFunction.Channel.
    function Base.show(io::IO, channel::GreenFunction.Channel)
        println(io, "symmetry:                $(channel.symmetry)  ")
        println(io, "gMultiplet:              $(channel.gMultiplet)  ")
    end

    """
    `struct  GreenFunction.Representation`
        ... defines a type to keep a representation of an (approximate) Green function that is associated
            with a given set of bound configurations.

        + approach          ::AbstractGreenApproach           ... Approach used to approximate the representation.
        + excitationScheme  ::Basics.AbstractExcitationScheme ... Applied excitation scheme with regard to the
                                                                  given bound configuration.
        + boundConfigs      ::Array{Configuration,1}          ... List of bound configurations to which the Green
                                                                  function refers to.
        + NoElectrons       ::Int64                           ... Number of electrons.
        + channels          ::Array{GreenFunction.Channel,1}  ... List of channels with different level symmetry.
    """
    struct Representation
        approach            ::AbstractGreenApproach
        excitationScheme    ::Basics.AbstractExcitationScheme
        boundConfigs        ::Array{Configuration,1}
        NoElectrons         ::Int64
        channels            ::Array{GreenFunction.Channel,1}
    end

    """
    `GreenFunction.Representation()`  ... constructor for an `empty` instance of GreenFunction.Representation.
    """
    function Representation()
        Representation( GreenFunction.SingleCSFwithoutCI(), Basics.NoExcitationScheme(), Configuration[], 0, GreenFunction.Channel[])
    end

    # `Base.show(io::IO, representation::GreenFunction.Representation)`  ... prepares a proper printout of the
    #       variable representation::GreenFunction.Representation.
    function Base.show(io::IO, representation::GreenFunction.Representation)
        println(io, "approach:                $(representation.approach)  ")
        println(io, "excitationScheme:        $(representation.excitationScheme)  ")
        println(io, "boundConfigs:            $(representation.boundConfigs)  ")
        println(io, "NoElectrons:             $(representation.NoElectrons)  ")
        println(io, "channels:                $(representation.channels)  ")
    end

    """
    `GreenFunction.computeRepresentation(configs::Array{Configuration,1}, multiplet::Multiplet, nm::Nuclear.Model,
                                         grid::Radial.Grid, settings::GreenFunction.Settings; output=true)`
        ... to compute an approximate Green function representation for the levels of the given bound configurations and for the
            given excitation scheme.
    """
    function computeRepresentation(configs::Array{Configuration,1}, multiplet::Multiplet, nm::Nuclear.Model,
                                   grid::Radial.Grid, settings::GreenFunction.Settings; output=true)
        println("")
        printstyled("GreenFunction.computeRepresentation(): The computation of approximate Green functions starts now ... \n", color=:light_green)
        printstyled("---------------------------------------------------------------------------------------------------- \n", color=:light_green)

        gChannels = GreenFunction.Channel[]

        # Generate all (non-relativistic) configurations from the bound configurations due to the given excitation scheme
        confList = Basics.generateConfigurationsForExcitationScheme(configs, settings.excitationScheme, settings.nMax, settings.lValues)
        # Print (if required) information about the generated configuration list
        if  settings.printBeforeComputation    GreenFunction.displayConfigurationList(confList)    end
        # Determine all (selected) levels and generate an averaged potential for these levels; use a ... potential for the selected levels
        pot = GreenFunction.generateMeanPotential(multiplet, nm, grid, settings)
        # Generate a full single-electron spectrum for this potential
        spectrum = Bspline.generateSpectrum()
        # Cycle over all selected level symmetries to generate in turn all channels of the requires Green function representation
        for  levelSymmetry  in  settings.levelSymmetries
            # Generate all CSFs for this symmetry
            csfs  = Basics.generateCsfList(levelSymmetry, confList)
            # Compute a multiplet representation due to the given approach
            gMultiplet = GreenFunction.computeMultipletForApproach(settings.approach, csfs, spectrum)
            channel    = GreenFunction.Channel( levelSymmetry, gMultiplet )
        end

        gfRep = GreenFunction.Representation( settings.approach, settings.excitationScheme, configs, confList[1].NoElectrons , gChannels)
        #
        # Print all results to screen
        GreenFunction.displayRepresentation(stdout, gfRep)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary    GreenFunction.displayRepresentation(iostream, gfRep)   end
        #
        if    output    return( gfRep )
        else            return( nothing )
        end
    end

    """
    `GreenFunction.displayConfigurationList(configs::Array{Configuration,1})`
        ... displays the gnerated list of configurations in a compact form; nothing is returned in this case.
    """
    function  displayConfigurationList(configs::Array{Configuration,1})
        println(" ")
        println("  Generated list of (non-relativistic) configurations that contribute to the Green functions representation:")
        println(" ")
        println("  ", TableStrings.hLine(60))
        for  i = 1:length(configs)
            sa = string("(", i, ")" );    sa = TableStrings.flushright(9, sa, na=3);   sa = sa * string(configs[i])
            println(sa)
        end
        println("  ", TableStrings.hLine(60))

        return( nothing )
    end

    """
    `GreenFunction.displayRepresentation(representation::GreenFunction.Representation)`
        ... to display
    """
    function  displayRepresentation(representation::GreenFunction.Representation)
        println(" ")
        println("  Selected GreenFunction levels:  ... not yet implemented")
        println(" ")
        println("  ", TableStrings.hLine(43))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(10, "Level"; na=2);                             sb = sb * TableStrings.hBlank(12)
        sa = sa * TableStrings.center(10, "J^P";   na=4);                             sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(14, "Energy"; na=4);
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=4)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(43))
        #
        for  outcome in outcomes
            sa  = "  ";    sym = LevelSymmetry( outcome.level.J, outcome.level.parity)
            sa = sa * TableStrings.center(10, TableStrings.level(outcome.level.index); na=2)
            sa = sa * TableStrings.center(10, string(sym); na=4)
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", outcome.level.energy)) * "    "
            println( sa )
        end
        println("  ", TableStrings.hLine(43))
        #
        return( nothing )
    end

    """
    `GreenFunction.generateMeanPotential(multiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid, settings::GreenFunction.Settings)`
        ... generates a mean radial potential for the selected levels of the multiplet that is used to compute the
            single-electron spectrum for the representation of the Green function; a potential::Radial.Potential
            is returned.
    """
    function  generateMeanPotential(multiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid, settings::GreenFunction.Settings)
        nuclearPotential  = Nuclear.nuclearPotential(nm, grid)
        ## wp1 = compute("radial potential: core-Hartree", grid, wLevel)
        ## wp2 = compute("radial potential: Hartree-Slater", grid, wLevel)
        ## wp3 = compute("radial potential: Kohn-Sham", grid, wLevel)
        ## wp4 = compute("radial potential: Dirac-Fock-Slater", grid, wLevel)

        if  settings.selectLevels
            wpx = zeros[grid.NoPoints];   nx = 0
            for  i = 1:length(multiplet.levels)
                if haskey(settings.selectedLevels, i)
                    wp = compute("radial potential: Dirac-Fock-Slater", grid, multiplet.level[i]);   np = length(wp.Zr)
                    wpx[1:np] = wpx[1:np] + wp.Zr[1:np];    nx = nx + 1
                end
            end
            wpx = wpx / nx
            wp = Radial.Potential(wp.name, wpx, grid)
        else
            wp = compute("radial potential: Dirac-Fock-Slater", grid, multiplet.levels[1].basis)
        end

        pot = Basics.add(nuclearPotential, wp)
        println("Mean potential for the generation of the (quasi-complete) single-electron spectrum:  $(pot)")

        return( pot )
    end

end # module


text:
using Test

using Roots
using Unitful
using UnitfulAtomic

using EquationsOfState
using EquationsOfState.Collections
using EquationsOfState.Find

# Data in the following tests are from
# https://github.com/materialsproject/pymatgen/blob/1f0957b8525ddc7d12ea348a19caecebe6c7ff34/pymatgen/analysis/tests/test_eos.py
@testset "Test data from Pymatgen" begin
    results = [
        25.987454833,
        26.9045702104,
        27.8430241908,
        28.8029649591,
        29.7848370694,
        30.7887887064,
        31.814968055,
        32.8638196693,
        33.9353435494,
        35.0299842495,
        36.1477417695,
        37.2892088485,
        38.4543854865,
        39.6437162376,
        40.857201102,
        42.095136449,
        43.3579668329,
        44.6456922537,
        45.9587572656,
        47.2973100535,
        48.6614988019,
        50.0517680652,
        51.4682660281,
        52.9112890601,
        54.3808371612,
        55.8775030703,
        57.4014349722,
        58.9526328669,
    ] .* u"angstrom^3"
    energies = [
        -7.63622156576,
        -8.16831294894,
        -8.63871612686,
        -9.05181213218,
        -9.41170988374,
        -9.72238224345,
        -9.98744832526,
        -10.210309552,
        -10.3943401353,
        -10.5427238068,
        -10.6584266073,
        -10.7442240979,
        -10.8027285713,
        -10.8363890521,
        -10.8474912964,
        -10.838157792,
        -10.8103477586,
        -10.7659387815,
        -10.7066179666,
        -10.6339907853,
        -10.5495538639,
        -10.4546677714,
        -10.3506386542,
        -10.2386366017,
        -10.1197772808,
        -9.99504030111,
        -9.86535084973,
        -9.73155247952,
    ] .* u"eV"
    isapprox(
        map(
            e -> findvolume(
                EnergyForm(),
                BirchMurnaghan3rd(
                    40.98926572528106 * u"angstrom^3",
                    0.5369258245417454 * u"eV/angstrom^3",
                    4.178644235500821 * u"1000mm/m",
                    -10.842803908240892 * u"eV",
                ),
                e,
                (eps(), 100) .* u"angstrom^3",
            ),
            energies,
        ),
        results,
    )
end

@testset "Test `findvolume` with random unit" begin
    pressures = collect(0:20:200) .* u"GPa"
    eos = BirchMurnaghan3rd(167 * u"angstrom^3", 2600 * u"kbar", 4.0 * u"1000mm/m")
    volumes = map(
        p -> findvolume(PressureForm(), eos, p, (eps() * u"bohr^3", eos.v0 * 1.3)),
        pressures,
    )
    @test isapprox(
        ustrip.(map(apply(PressureForm(), eos), volumes) - pressures),
        zeros(11),
        atol = 1e-5,
    )
end # testset


text:
module Yields

import BSplineKit
import ForwardDiff
using LinearAlgebra
using UnicodePlots

# don't export type, as the API of Yields.Zero is nicer and
# less polluting than Zero and less/equally verbose as ZeroYieldCurve or ZeroCurve
export rate, discount, accumulation, forward,
    LinearSpline, CubicSpline

include("Rate.jl")
include("AbstractYield.jl")
include("utils.jl")
include("bootstrap.jl")
include("SmithWilson.jl")
include("generics.jl")
include("RateCombination.jl")

end


text:
#RSDeltaSigmaPort: Generate frequency spectrum plots
#-------------------------------------------------------------------------------

#==Frequency-spectrum plots
===============================================================================#
function plotSpec()
	plot = cons(:plot, linlin, title = "Signal Spectrum", legend=true,
		xyaxes=set(xmin=0, xmax=0.5, ymin=-120, ymax=0),
		labels=set(xaxis="Normalized Frequency", yaxis="dBFS"),
	)
	return plot
end

function plotSpec!(plot, sig; M::Int=1, id::String="signal", color=:blue)
	sig = conv2seriesmatrix2D(sig)
	M, N = size(sig)
	𝑓pts = div(N,2)+1 #Only need half the spectrum
	𝑓 = range(0, stop=0.5, length=𝑓pts)
	𝑓 = ones(M)*𝑓' #2D frequency matrix for downstream algorithms
	spec = fft(sig)/(M*N/4)
	specdB = dbv.(spec[:,1:𝑓pts])

	#Convert to waveforms:
	specdB = waveform(𝑓, specdB)

	push!(plot,
		cons(:wfrm, specdB, line=set(style=:solid, color=color, width=2), label=id),
	)
	return plot
end

function plotSpec(sig; M::Int=1, id::String="signal", color=:blue)
	plot = plotSpec()
	return plotSpec!(plot, sig, id=id, color=color, M=M)
end

#==Modulator output spectrum
===============================================================================#
function plotModSpectrum()
	plot = plotSpec()
	plot.title = "Modulator Output Spectrum"
	return plot
end

function plotModSpectrum!(plot, specinfo; id::String="Simulation", color=:blue)
	𝑓 = specinfo.freq

	#Convert to waveforms:
	ePSD = waveform(𝑓, specinfo.ePSD)
	specdB = waveform(𝑓, specinfo.specdB)

	snr_str = @sprintf("SNR = %4.1fdB", specinfo.SNR) #orig. coords: (0.05,-10)
	nbw_str = @sprintf("NBW = %4.1E × 𝑓s", specinfo.NBW) #orig. coords: (0.5, -90)
	ann_str = string(snr_str, "\n", nbw_str)
	push!(plot,
		cons(:wfrm, specdB, line=set(style=:solid, color=color, width=2), label=id),
		cons(:wfrm, ePSD, line=set(style=:solid, color=:red, width=2), label="Expected PSD"),
		cons(:atext, ann_str, y=-90, reloffset=set(x=0.5), align=:cc),
	)

	return plot
end

"""`plotModSpectrum(sig, NTF, iband, f; M=1, title, id)`

Generate plot of modulator output spectrum:
 - `sig`: Time domain ΔΣ signal.
 - `NTF`:
 - `f`: Bin location of input signal.
 - `M`: Number of modulator steps (= nlev-1).
"""
function plotModSpectrum(specinfo; id::String="Simulation", color=:blue)
	plot = plotModSpectrum()
	plotModSpectrum!(plot, specinfo, id=id, color=color)
	return plot
end

#==plotSpectrum()
===============================================================================#
function plotSpectrum()
	plot = cons(:plot, loglin, title = "Smoothed Spectrum", legend=true,
		xyaxes=set(xmin=1e-3, xmax=0.5, ymin=-120, ymax=0),
		labels=set(xaxis="Normalized Frequency", yaxis=""),
	)
	return plot
end

function plotSpectrum!(plot, X, fin; nbin::Int=8, n::Int=3, lw=1, color=:blue, id="")
	(f, p) = logsmooth(X, fin, nbin=nbin, n=n)
	pw = waveform(f, p)
	push!(plot,
		cons(:wfrm, pw, line=set(style=:solid, color=color, width=lw), label=id),
	)
	return plot
end

"""`plotSpectrum(X, fin, nbin=8, n=3, lw=1)`

Plot a smoothed spectrum.
"""
function plotSpectrum(X, fin; nbin::Int=8, n::Int=3, lw=1, color=:blue, id="")
	plot = plotSpectrum()
	plotSpectrum!(plot, X, fin, nbin=nbin, n=n, lw=lw, color=color, id=id)
	return plot
end

#==plotExampleSpectrum
===============================================================================#
"""`plotExampleSpectrum(NTF, OSR=64, M=1, f0=0, quadrature=false, ampdB=-3, ftest=nothing, N=2^12, sizes=nothing)`

#Inputs
 - `M`: Number of modulator steps (= nlev-1).
 - `ampdB`: Test tone amplitude, relative to full-scale.
"""
function plotExampleSpectrum(NTF; OSR::Int=64, M::Int=1, f0::Float64=0, quadrature::Bool=false,
		ampdB::Float64=-3.0, ftest=nothing, N::Int=2^12, sizes=nothing)
	#Computed defaults
	if isnothing(ftest)
		ftest = default_ftest(OSR, f0=f0, quadrature=quadrature)
	end
	fband = default_fband(OSR, f0=f0, quadrature=quadrature)
	delta = 2

	usrsizes = sizes
	sizes = Dict{Symbol, Any}( #Defaults
		:lw => 1, #LineWidth
		:ms => 5, #MarkerSize
		:fs => 12, #FontSize
		:fw => "normal", #FontWeight
	)
	if isa(usrsizes, Dict) #Overwrite with user supplied values
		push!(sizes, usrsizes...)
	end

	#Plot an example spectrum
#	fin = round(Int, ftest*N)
	local u, simresult
	if !quadrature
		phi0=π/2 #Switch from sin to cos
		(u, iftest) = genTestTone_sin(ampdB, ftest, M=M, phi0=phi0, N=N)
		simresult = simulateDSM(u, NTF, nlev=M+1)
	else
		(u, iftest) = genTestTone_quad(ampdB, ftest, M=M, N=N)
		simresult = simulateQDSM(u, NTF, nlev=M+1)
	end

	plot = cons(:plot, linlin, title = "Modulator Output Spectrum", legend=true,
		labels=set(xaxis="Normalized Frequency", yaxis="PSD [dBFS/NBW]"),
		xyaxes=set(ymin=-140, ymax=0),
	)

	specinfo = calcSpecInfo(simresult.v, NTF, fband, ftest, M=M, quadrature=quadrature)
	if !quadrature
		if 0==f0
#			set(plot, loglin, xyaxes=set(xmin=.001))
		end

		#circular smoothing; not re-scaled
		spec_smoothed = circ_smooth(abs.(specinfo.spec0).^2, 16)
		A = dbv(specinfo.spec0[specinfo.itest+1])

		#Convert to waveforms & add to plot:
		𝑓 = specinfo.freq; 𝑓pts = length(𝑓)
		specw = waveform(𝑓, specinfo.specdB)
		spec_sw = waveform(𝑓, dbp.(spec_smoothed[1:𝑓pts]))
		ePSDw = waveform(𝑓, specinfo.ePSD)

		push!(plot,
			cons(:wfrm, specw, line=set(style=:solid, color=:cyan, width=2), label=""),
			cons(:wfrm, spec_sw, line=set(style=:solid, color=:blue, width=2), label=""),
			cons(:wfrm, ePSDw, line=set(style=:solid, color=:magenta, width=1), label=""),
		)
		alignleft = f0<0.25
		align = alignleft ? :tl : :tr
		annx  = alignleft ? f0+0.5/OSR : f0-0.5/OSR
		tstr = @sprintf(" SQNR = %.1fdB\n @ A = %.1fdBFS\n  & OSR = %.0f\n",
			specinfo.SNR, A, OSR
		)
		nbwstr = @sprintf("NBW=%.1e ", specinfo.NBW)
		#TODO: use `sizes` atext:
		push!(plot,
			cons(:atext, tstr, x=annx, y=-5, align=align),
			cons(:atext, nbwstr, x=0.5, y=-135, align=:cr),
		)
	else
		error("plotExampleSpectrum not implemented for quadrature modulators.")
#=
		spec0 = fftshift(spec0/2)
		freq = range(-0.5, stop=0.5, length=N+1)
		freq = deleteat!(collect(freq), N+1)
		plot(freq,dbv(spec0),'c','Linewidth',1);
		spec_smoothed = circ_smooth(abs(spec0).^2, 16);
		plot(freq, dbp(spec_smoothed), 'b', 'Linewidth', 3);
		Snn = abs(evalTF(NTF,exp(2i*pi*freq))).^2 * 2/12 * (delta/M)^2;
		plot(freq, dbp(Snn*NBW), 'm', 'Linewidth', 1);
		SNR = calculateSNR(spec0(N/2+1+[f1_bin:f2_bin]),fin-f1_bin);
		msg = sprintf('SQNR = %.1fdB\n @ A=%.1fdBFS & OSR=%.0f\n', ...
		SNR, dbv(spec0(N/2+fin+1)), OSR );
		if f0>=0
			text(f0-0.05, -15, msg,'Hor','Right');
		else
			text(f0+0.05, -15, msg,'Hor','Left');
		end
		text(-0.5,-135,sprintf(' NBW=%.1e',NBW),'hor','left');
		figureMagic([-0.5 0.5],0.125,2, [-140 0],10,2);
=#
	end

	return plot
end

"""`plotExampleSpectrum(dsm, NTF=nothing; ampdB=-3, ftest=nothing, N=2^12, sizes=sizes)`
"""
function plotExampleSpectrum(dsm::AbstractDSM, NTF=nothing; ampdB=-3.0, ftest=nothing, N=2^12, sizes=nothing)
	isnothing(NTF) && (NTF = synthesizeNTF(dsm))
	return plotExampleSpectrum(NTF, OSR=dsm.OSR, M=dsm.M, f0=dsm.f0, quadrature=isquadrature(dsm),
		ampdB=ampdB, ftest=ftest, N=N, sizes=sizes
	)
end

#Last line


text:
println("===============================================")
println("Testing elements.jl")
println("===============================================")

function springpot_spring_reduce_G()
    t = Vector{RheoFloat}(0.0:0.1:10.0)

    springpot = relaxmod(RheoModel(Springpot, cᵦ = 1.5, β=0),t)
    spring = relaxmod(Spring,t, k=1.5)

    all(i -> isapprox(springpot[i], spring[i]), eachindex(t))
end
@test springpot_spring_reduce_G()

function springpot_spring_reduce_J()
    t = Vector{RheoFloat}(0.0:0.1:10.0)

    springpot_J = creepcomp(Springpot, cᵦ = 1.5, β=0)
    springpot = springpot_J(t)
    spring = creepcomp(Spring,t, [1.5])

    all(i -> isapprox(springpot[i], spring[i]), eachindex(t))
end
@test springpot_spring_reduce_J()

function springpot_spring_reduce_Gp()
    ω = Vector{RheoFloat}(0.0:0.1:10.0)

    springpot = storagemod(Springpot, ω, cᵦ = 1.5, β=0)
    spring = storagemod(Spring, ω, [1.5])

    all(i -> isapprox(springpot[i], spring[i]), eachindex(ω))
end
@test springpot_spring_reduce_Gp()

function springpot_spring_reduce_Gpp()
    ω = Vector{RheoFloat}(0.0:0.1:10.0)

    springpot_Gpp = lossmod(RheoModel(Springpot, cᵦ = 1.5, β=0))
    springpot = springpot_Gpp(ω)
    spring = lossmod(Spring, ω, [1.5])

    all(i -> isapprox(springpot[i], spring[i]), eachindex(ω))
end
@test springpot_spring_reduce_Gpp()

function springpot_dashpot_reduce_G()
    # This gives NaN as defined presently,
    # no test written
    return true
end
@test springpot_dashpot_reduce_G()

function springpot_dashpot_reduce_J()
    t = Vector{RheoFloat}(0.0:0.1:10.0)

    springpot = creepcomp(Springpot,t, [1.5, 1.0])
    dashpot = creepcomp(Dashpot,t, [1.5])

    all(i -> isapprox(springpot[i], dashpot[i]), eachindex(t))
end
@test springpot_dashpot_reduce_J()

function springpot_dashpot_reduce_Gp(tol)
    ω = Vector{RheoFloat}(0.0:0.1:10.0)

    springpot = storagemod(Springpot, ω, [1.5, 1.0])
    dashpot = real(dynamicmod(Dashpot, ω, [1.5]))

    # lower tolerance due to cos(0) innacuracy in Julia
    all(i -> isapprox(springpot[i], dashpot[i], atol=tol), eachindex(ω))
end
@test springpot_dashpot_reduce_Gp(tol)

function springpot_dashpot_reduce_Gpp()
    ω = Vector{RheoFloat}(0.0:0.1:10.0)

    springpot = lossmod(Springpot, ω, [1.5, 1.0])
    dashpot = imag(dynamicmod(Dashpot, ω, [1.5]))

    all(i -> isapprox(springpot[i], dashpot[i]), eachindex(ω))
end
@test springpot_dashpot_reduce_Gpp()


text:
const Days = convert(Vector{Float64},slp[:Days])
const vf = VectorReMat(slp[:Subject],hcat(ones(length(Days)),Days)',:Subject,["(Intercept)","Days"])
const Reaction = convert(Array,slp[:Reaction])

@test size(vf) == (180,36)
const vrp = vf'vf
@test (vf'ones(size(vf,1)))[1:4] == [10.,45,10,45]
@test isa(vrp,MixedModels.HBlkDiag{Float64})
@test eltype(vrp) == Float64
@test size(vrp) == (36,36)
const rhs1 = ones(36,2)
const x = similar(rhs1)
const b1 = copy(vrp.arr[:,:,1]) + I
@test sub(MixedModels.inflate!(vrp).arr,:,:,1) == b1
const cf = cholfact(b1)
@test triu!(sub(cholfact!(vrp).arr,:,:,1)) == cf[:U]

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text: Dict{String,Any}("point" => Dict{String,Any}("x" => Dict{String,Any}("value" => "1","type" => "integer"),"y" => Dict{String,Any}("value" => "2","type" => "integer")),"name" => Dict{String,Any}("first" => Dict{String,Any}("value" => "Tom","type" => "string"),"last" => Dict{String,Any}("value" => "Preston-Werner","type" => "string")),"str-key" => Dict{String,Any}("a" => Dict{String,Any}("value" => "1","type" => "integer")),"simple" => Dict{String,Any}("a" => Dict{String,Any}("value" => "1","type" => "integer")),"table-array" => Any[Dict{String,Any}("a" => Dict{String,Any}("value" => "1","type" => "integer")), Dict{String,Any}("b" => Dict{String,Any}("value" => "2","type" => "integer"))])

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
#utilities
#some of these are candidates for Finometrics

@inline function finiteormissing!(df::DataFrame, s::Symbol)::Nothing
  df[!,s] .= (f::MFloat64->coalesce(isfinite(f),false) ? f : missing).(df[!,s])

  return nothing
end

@inline function finiteormissing!(df::AbstractDataFrame, s::Symbol)::Nothing
  @mpar for r ∈ eachrow(df)
    r[s] = coalesce(isfinite(r[s]),false) ? r[s] : missing
  end

  return nothing
end

#@inline Base.isfinite(::Missing)=missing

#NOTE: The below adjusts for logged vs unlogged returns
#const RETAINED_FACT_COLS = (( RETURN_FIELD == :lret) ?
#  [RETURN_FIELD; FACTOR_FIELDS; :lrfr, :lmkt] :
#  [RETURN_FIELD; FACTOR_FIELDS; :rfr, :mkt]

#consider promiting to the below to finometrics
@inline obliterate!(df::DataFrame)::DataFrame = select!(df, Not(names(df)))

#checks if a column or multiple columns are a key
@inline function checkdfkey(df::AbstractDataFrame, keys::Union{Symbol, Vector{Symbol}})
  N::Int = size(df,1)
  Nunique::Int = length(groupby(df, keys))

  (N==Nunique) || error("$keys is not a key!!")

  return nothing
end

#converts gzip to lz4
#=function gz2lz4(file::String; filetype::String = "csv", location::String = "",
  gzextension::String = "gz", lz4extension::String = "lz4",
  deleteold::Bool = false)

  infile::String = "$location\\$file.$filetype.$gzextension"
  outfile::String = "$location\\$file.$filetype.$lz4extension"

  isfile(infile) || error("file not found. Full path: $location\\$file.$filetype.$gzextension")

  #Below doesnt work
  indata = open(infile)
  outdata = open(outfile, "w")
  s = GzipDecompressorStream(LZ4HCCompressorStream(outdata))
  write(s, indata)
  close(s)
end=#

#converts gzip to lz4
function gz2zstd(file::String; filetype::String = "csv", location::Union{String,Nothing} = nothing,
  gzextension::String = "gz", zstdextension::String = "zstd",
  deleteold::Bool = false, level=1)

  local prefix::String

  prefix = isnothing(location) ?  "" : "$location\\"
  infile::String = "$prefix$file.$filetype.$gzextension"
  outfile::String = "$prefix$file.$filetype.$zstdextension"

  isfile(infile) || error("file not found. Full path: $prefix$file.$filetype.$gzextension")

  #Below doesnt work
  indata = open(infile)
  outdata = open(outfile, "w")
  s = GzipDecompressorStream(ZstdCompressorStream(outdata, level=level))

  try
    write(s, indata)
  finally
    close(s)
  end
end

function recompressfiles(files = FILES_TO_RECOMPRESS; datapath=DATA_PATH)::Nothing

  for f ∈ files
    gz2zstd(f, location=datapath)
  end

  return nothing
end

recompressfile(file; datapath=DATA_PATH)::Nothing = recompressfiles(
  [file]; datapath=datapath)

#check if all columns in a df or subdf are equal
#extra logic is to handle missing values
@inline function eachrowequal(df::AbstractDataFrame; cols::Vector{Symbol}=propertynames(df))::Bool
  for c ∈ eachcol(view(df, :, cols))
    completec = collect(skipmissing(c))
    (length(completec) == 0) && continue #all missing case

    #check for equality
    if  (length(completec) == length(c)) && (all(c[1] .== c))
     continue
    end

    return false
  end
  return true
end

#some useful functions for conforming numbers
@inline floororx(::Missing, ::Float64 = DENOM_FLOOR) = missing
@inline floororx(x::Real, floor::Float64 = DENOM_FLOOR) = max(x, floor)

@inline boundunit(::Missing, ::Float64=UNIT_TOL) = missing
@inline boundunit(x::Real, tol::Float64=UNIT_TOL) = max(min(x, 1.0 - tol), tol)
@inline boundunit1(x::Union{Real,Missing}) = boundunit(x, 0.0)

@inline boundabsunit(::Missing, ::Float64=UNIT_TOL) = missing
@inline boundabsunit(x::Real, tol::Float64=UNIT_TOL) = max(min(x, 1.0 -tol), -1.0 + tol)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module BalancedIncompleteBlockDesigns
using ChooseOptimizer
using JuMP, MathProgBase, Primes

export BIBD, projective, BIBD_check, blocks

"""
`BIBD(b,v,r,k,λ)` finds a `(b,v,r,k,λ)`-balanced incomplete block design.
The return is a `v`-by-`b` matrix of 0s and 1s whose i,j-entry  is 1
exactly when vertex i is in block j. Errors erupt if no such design exists.

The shorter form `BIBD(b,k,λ)` returns a symmetric design and is equivalent
to `BIBD(b,b,k,k,λ)`.
"""
function BIBD(b::Int,v::Int,r::Int,k::Int,l::Int)
    errmsg = "No ($b,$v,$r,$k,$l)-BIBD can be found"
    #set_solver(Gurobi)
    #set_solver_verbose(true)

    @assert b>0 && v>0 && r>0 && k>0 && l>0 && k<v "Invalid parameters: "*errmsg
    @assert b*k==v*r && r*(k-1)==l*(v-1) errmsg

    M = Model(get_solver())

    @variable(M, x[1:v,1:b], Bin)       # 1{vtx i in block B}
    @variable(M, y[1:v,1:v,1:b], Bin)   # 1{vtcs i,j in block B} where i≠j

    # Basics for y
    for i=1:v
        for B=1:b
            @constraint(M, y[i,i,B]==0)  # i≠j
        end
    end

    for i=1:v
        for j=1:v
            for B=1:b
                @constraint(M, y[i,j,B]==y[j,i,B]) # symmetric in i,j
            end
        end
    end

    # Each vertex is in r blocks
    for i=1:v
        @constraint(M, sum(x[i,B] for B=1:b) == r)
    end

    # Each block has k vertices
    for B=1:b
        @constraint(M, sum(x[i,B] for i=1:v) == k)
    end

    # every pair i≠j are in l blocks together
    for i=1:v
        for j=1:v
            if i!=j
                @constraint(M, sum(y[i,j,B] for B=1:b) == l)
            end
        end
    end

    # if i,j in B, then x[i,B] and x[j,B] must be 1
    for i=1:v
        for j=1:v
            for B=1:b
                @constraint(M, x[i,B] >= y[i,j,B])
            end
        end
    end

    # if i,j in B, force y[i,j,B] to 1
    for i=1:v
        for j=1:v
            if i!=j
                for B=1:b
                    @constraint(M, y[i,j,B] >= x[i,B] + x[j,B] - 1)
                end
            end
        end
    end

    optimize!(M)
    status = Int(termination_status(M))

    if status != 1
        error(errmsg)
    end

    X = Int.(value.(x))
end

BIBD(b::Int,k::Int,l::Int) = BIBD(b,b,k,k,l)

"""
`projective(n)` returns a projective plane of order `n` as a 0,1-matrix.
The rows of the matrix correspond to the points and the columns to the lines.
The i,j-entry is 1 iff point i is on line j.
"""
function projective(n::Int)
    @assert n>1 "Order of the plane must be at least 2"
    return BIBD(n*n+n+1,n+1,1)
end

function _all_same(vec::Array{T,1}) where T
    n = length(vec)
    if n < 2
        return true
    end
    for k=1:n-1
        if vec[k] != vec[k+1]
            return false
        end
    end
    return true
end

"""
`BIBD_check(A)` checks if the matrix `A` represents a
balanced incomplete block design. If so, it returns
its parameters as a 5-tuple `(b,v,r,k,λ)`.
"""
function BIBD_check(A::Array{Int,2})
    if maximum(A)>1 || minimum(A)<0
        error("Not a 0,1-matrix")
    end

    v,b = size(A)

    rlist = sum(A;dims=2)[:]
    if !_all_same(rlist)
        error("Replication number fail")
    end
    r = rlist[1]

    klist = sum(A;dims=1)[:]
    if !_all_same(klist)
        error("Block size fail")
    end

    k = klist[1]

    if k==v
        error("Incompleteness fail: v==k==$v")
    end

    B = A*A'
    l = B[1,2]
    for i=1:v
        for j=1:v
            if i!=j && B[i,j]!=l
                error("Lambda fail")
            end
        end
    end

    return (b,v,r,k,l)
end

"""
`blocks(A)` gives a list of the blocks of the balanced incomplete block design
specified by the matrix `A` (e.g., the output of `BIBD`)
"""
function blocks(A::Array{Int,2})
    b,v,r,k,l = BIBD_check(A)
    [ Set(findall(A[:,j] .> 0)) for j=1:b]
end

include("prime_plane.jl")

end  # end of module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct OauthV2Api <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_oauthV2Access(_api::OauthV2Api, code::String; client_id=nothing, client_secret=nothing, redirect_uri=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/oauth.v2.access", ["slackAuth"])
    Swagger.set_param(_ctx.query, "client_id", client_id)  # type String
    Swagger.set_param(_ctx.query, "client_secret", client_secret)  # type String
    Swagger.set_param(_ctx.query, "code", code)  # type String
    Swagger.set_param(_ctx.query, "redirect_uri", redirect_uri)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Exchanges a temporary OAuth verifier code for an access token.
Param: code::String (required)
Param: client_id::String
Param: client_secret::String
Param: redirect_uri::String
Return: Dict{String, Any}
"""
function oauthV2Access(_api::OauthV2Api, code::String; client_id=nothing, client_secret=nothing, redirect_uri=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_oauthV2Access(_api, code; client_id=client_id, client_secret=client_secret, redirect_uri=redirect_uri, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function oauthV2Access(_api::OauthV2Api, response_stream::Channel, code::String; client_id=nothing, client_secret=nothing, redirect_uri=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_oauthV2Access(_api, code; client_id=client_id, client_secret=client_secret, redirect_uri=redirect_uri, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export oauthV2Access


text:
"""
    module Simulate_Billiards_OneBall

Simulates a billiard ball on a (non-standard) square billiard table
that is hitting the cushing several times. This model is used to
test the collision with the cushing.
"""
module Simulate_Billiards_OneBall

using  Modia3D
import Modia3D.ModiaMath

vmatBalls = Material(color="Red"      , transparency=0.0)
vmatTable = Material(color="DarkGreen", transparency=0.4)

diameter = 0.06
radius   = diameter/2

TableX = 1.24
TableY = 1.12
TableZ = 0.03
LyBox  = 0.03
hz     = 7*radius/5
hzz    = hz/3

@assembly Cushion(world) begin
  cushionPart11 = Solid(SolidBox(LyBox, TableY, hz) , "BilliardCushion", vmatTable)
  cushionPart12 = Solid(SolidBox(1.5*LyBox, TableY, hzz,rsmall=0.0) , "BilliardCushion", vmatTable; contactMaterial="BilliardCushion")
  cushionPart21 = Solid(SolidBox(TableX, LyBox , hz) , "BilliardCushion", vmatTable)
  cushionPart22 = Solid(SolidBox(TableX, 1.5*LyBox, hzz,rsmall=0.0) , "BilliardCushion", vmatTable; contactMaterial="BilliardCushion")

  cushion11 = Object3D(world    , cushionPart11, fixed=true, r=[TableX/2 - LyBox/2 , 0.0                 , hz/2])
  cushion12 = Object3D(cushion11, cushionPart12, fixed=true, r=[-LyBox/4           , 0.0                 , hz/2+hzz/2])
  cushion31 = Object3D(world    , cushionPart11, fixed=true, r=[-TableX/2 + LyBox/2, 0.0                 , hz/2])
  cushion32 = Object3D(cushion31, cushionPart12, fixed=true, r=[LyBox/4            , 0.0                 , hz/2+hzz/2])
  cushion21 = Object3D(world    , cushionPart21, fixed=true, r=[0.0                , TableY/2 - LyBox/2  , hz/2])
  cushion22 = Object3D(cushion21, cushionPart22, fixed=true, r=[0.0                , -LyBox/4            , hz/2+hzz/2])
  cushion41 = Object3D(world    , cushionPart21, fixed=true, r=[0.0                , -TableY/2 + LyBox/2 , hz/2])
  cushion42 = Object3D(cushion21, cushionPart22, fixed=true, r=[0.0                , LyBox/4             , hz/2+hzz/2])
end

distance_balls = sqrt(3)/2*diameter
dist           = 0.0001
solidBall      = Solid(SolidSphere(diameter), "BilliardBall", vmatBalls ; contactMaterial="BilliardBall")
solidTable     = Solid(SolidBox(TableX, TableY, TableZ), "BilliardTable", vmatTable; contactMaterial="BilliardTable")

@assembly Billiards_OneBall() begin
  world   = Object3D()
  cushion = Cushion(world)
  table   = Object3D(world, solidTable, fixed=true , r=[0.0, 0.0, -TableZ/2])
  ball    = Object3D(world, solidBall , fixed=false, r=[0.0, 0.0, diameter/2], v_start=[6.0, 3.0, 0.0], visualizeFrame=true)
end

gravField         = UniformGravityField(g=9.81, n=[0,0,-1])
billiards_OneBall = Billiards_OneBall(sceneOptions=SceneOptions(gravityField=gravField, visualizeFrames=false,
                                                                defaultFrameLength=0.1, enableContactDetection=true))

# Modia3D.visualizeAssembly!( billiards_OneBall )
model = SimulationModel( billiards_OneBall, maxNumberOfSteps=1000 )   # default maxNumberOfSteps gives an error
# ModiaMath.print_ModelVariables(model)
result = ModiaMath.simulate!(model; stopTime=0.5, tolerance=1e-8, log=false)

ModiaMath.plot(result, [("ball.r[1]", "ball.r[2]"),
                        ("ball.v[1]", "ball.v[2]"),
                        ("ball.w[1]", "ball.w[2]")])

println("... success of examples/collisions/Simulate_Billiards_OneBall!")
end


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

__precompile__(true)

module ValueShapes

using ArraysOfArrays
using Distributions
using ElasticArrays
using FillArrays
using Random
using Statistics
using StatsBase

import EmpiricalDistributions
import IntervalSets
import Tables
import TypedTables

include("value_shape.jl")
include("value_accessor.jl")
include("scalar_shape.jl")
include("array_shape.jl")
include("const_value_shape.jl")
include("named_tuple_shape.jl")
include("functions.jl")
include("distributions.jl")
include("const_value_dist.jl")
include("named_tuple_dist.jl")

end # module


text:
module Diffractor

using StructArrays

export ∂⃖, gradient

include("runtime.jl")
include("interface.jl")
include("utils.jl")
include("tangent.jl")
include("jet.jl")

include("stage1/generated.jl")
include("stage1/termination.jl")
include("stage1/forward.jl")
include("stage1/recurse_fwd.jl")
include("stage1/mixed.jl")
include("stage1/broadcast.jl")

include("extra_rules.jl")

include("higher_fwd_rules.jl")

include("debugutils.jl")

end


text:
@with_kw mutable struct Sobol <: GSAMethod
    order::Array{Int}=[0,1]
    nboot::Int=0
    conf_int::Float64=0.95
end

mutable struct SobolResult{T1,T2,T3,T4}
    S1::T1
    S1_Conf_Int::T2
    S2::T3
    S2_Conf_Int::T4
    ST::T1
    ST_Conf_Int::T2
end

function fuse_designs(A,B)
    d = size(A,1)
    Aᵦ = [copy(A) for i in 1:d]
    for i in 1:d
        Aᵦ[i][i,:] = B[i,:]
    end
    Aᵦ
end

function gsa(f,method::Sobol,A::AbstractMatrix,B::AbstractMatrix;
             batch=false,Ei_estimator = :Jansen1999)
    Aᵦ = fuse_designs(A,B)
    d,n = size(A)
    multioutput = false
    all_points = hcat(A,B,reduce(hcat,Aᵦ))

    if batch
        all_y = f(all_points)
        multioutput = all_y isa AbstractMatrix
    else
        _y = [f(all_points[:,i]) for i in 1:size(all_points,2)]
        multioutput = !(eltype(_y) <: Number)
        all_y = multioutput ? reduce(hcat,_y) : _y
    end

    if !multioutput

        Ey = mean(all_y[1:2n])
        Vary = var(all_y[1:2n])

        fA = all_y[1:n]
        fB = all_y[(n+1):(2n)]
        fAⁱ= [all_y[(j*n+1):((j+1)*n)] for j in 2:(d+1)]
        if 2 in method.order
            fAⁱʲ= all_y[Int(end - (d*(d-1))/2):end]
        end
        Vᵢ = [sum(fB.*(fAⁱ[i].-fA)) for i in 1:d]./n
        if 2 in method.order
            Vᵢⱼ = [sum(abs2, fAⁱ[i] - fAⁱ[j]) for i in 1:d for j in i+1:d]./(2n)
        end
        if Ei_estimator == :Homma1996
            Eᵢ = [Vary .- sum(fA .* fAⁱ[i])./(n) + Ey.^2 for i in 1:d]
        elseif Ei_estimator == :Sobol2007
            Eᵢ = [sum(abs2,fA-fAⁱ[i]) for i in 1:d]./(2n)
        elseif Ei_estimator == :Jansen1999
            Eᵢ = [sum(fA.*(fA.-fAⁱ[i])) for i in 1:d]./(n)
        end

    else

        Ey = mean(all_y[:,1:2n],dims=2)
        Vary = var(all_y[:,1:2n],dims=2)

        fA = all_y[:,1:n]
        fB = all_y[:,(n+1):(2n)]
        fAⁱ= [all_y[:,(j*n+1):((j+1)*n)] for j in 2:(d+1)]

        Vᵢ = reduce(hcat,[sum(fB.*(fAⁱ[i].-fA),dims=2) for i in 1:d]./n)

        if 2 in method.order
            Vᵢⱼ = reduce(hcat, [sum(abs2, fAⁱ[i] - fAⁱ[j],dims=2) for i in 1:d for j in i+1:d]./(2n))
        end
        if Ei_estimator == :Homma1996
            Eᵢ = reduce(hcat,[Vary .- sum(fA .* fAⁱ[i],dims=2)./(n) + Ey.^2 for i in 1:d])
        elseif Ei_estimator == :Sobol2007
            Eᵢ = reduce(hcat,[sum(abs2,fA-fAⁱ[i],dims=2) for i in 1:d]./(2n))
        elseif Ei_estimator == :Jansen1999
            Eᵢ = reduce(hcat,[sum(fA.*(fA.-fAⁱ[i]),dims=2) for i in 1:d]./(n))
        end

    end
    if 2 in method.order
        Sᵢⱼ= (Vᵢⱼ)./Vary
    end

    Sᵢ = Vᵢ ./Vary
    Tᵢ = Eᵢ ./Vary
    SobolResult(Sᵢ,nothing,2 in method.order ? Sᵢⱼ : nothing ,nothing,Tᵢ,nothing)
end