text:
"""
Representation of a trajectory, ie a Hamiltonian with a discrete integrator that
also checks for divergence.
"""
struct TrajectoryNUTS{Tf,S}
    # "Hamiltonian."
    # H::TH
    "Log density of z (negative log energy) at initial point."
    π₀::Tf
    "Stepsize for leapfrog."
    ϵ::Tf
    "Smallest decrease allowed in the log density."
    min_Δ::Tf
    "Turn statistic configuration."
    turn_statistic_configuration::S
end

function move(tree::Tree{D,T,L}, H::Hamiltonian, trajectory::TrajectoryNUTS, z::PhasePoint{D,T,L}, fwd::Bool) where {D,T,L}
    @unpack ϵ = trajectory
    leapfrog(tree, H, z, fwd ? ϵ : -ϵ)
end

###
### proposals
###

"""
$(SIGNATURES)
Random boolean which is `true` with the given probability `exp(logprob)`, which can be `≥ 1`
in which case no random value is drawn.
"""
function rand_bool_logprob(rng::AbstractRNG, logprob)
    logprob ≥ 0 || (randexp(rng, Float64) > -logprob)
end

function calculate_logprob2(::TrajectoryNUTS, is_doubling, ω₁, ω₂, ω)
    biased_progressive_logprob2(is_doubling, ω₁, ω₂, ω)
end

function combine_proposals(rng, tree::Tree, ::TrajectoryNUTS, z₁, z₂, logprob2::Real, is_forward::Bool)
    z, flag = rand_bool_logprob(rng, logprob2) ? (z₂, z₁.flag) : (z₁, z₂.flag)
    # is_doubling || free_z!(tree, flag)
    free_z!(tree, flag)
    z
end
# function combine_proposals(
    # rng, tree::Tree, ::TrajectoryNUTS, z₁, z₂, logprob2::Real, is_forward, guard₁, guard₂
# )
    # z, flag, guard = rand_bool_logprob(rng, logprob2) ? (z₂, z₁.flag, guard₁) : (z₁, z₂.flag, guard₂)
    # guard || free_z!(tree, flag)
    # z
# end

###
### statistics for visited nodes
###

struct AcceptanceStatistic{T}
    """
    Logarithm of the sum of metropolis acceptances probabilities over the whole trajectory
    (including invalid parts).
    """
    log_sum_α::T
    "Total number of leapfrog steps."
    steps::Int32
end

function combine_acceptance_statistics(A::AcceptanceStatistic, B::AcceptanceStatistic)
    AcceptanceStatistic(logaddexp(A.log_sum_α, B.log_sum_α), A.steps + B.steps)
end

"""
$(SIGNATURES)
Acceptance statistic for a leaf. The initial leaf is considered not to be visited.
"""
function leaf_acceptance_statistic(Δ, is_initial)
    is_initial ? AcceptanceStatistic(oftype(Δ, -Inf), zero(Int32)) : AcceptanceStatistic(Δ < zero(Δ) ? Δ : zero(Δ), one(Int32))
end

"""
$(SIGNATURES)
Return the acceptance rate (a `Real` betwen `0` and `1`).
"""
acceptance_rate(A::AcceptanceStatistic) = (a = exp(A.log_sum_α) / A.steps; a < one(a) ? a : one(a))

combine_visited_statistics(::TrajectoryNUTS, v, w) = combine_acceptance_statistics(v, w)

###
### turn analysis
###

"Statistics for the identification of turning points. See Betancourt (2017, appendix)."
struct GeneralizedTurnStatistic{D,T,L}
    p♯₋::FlaggedVector{D,T,L}
    p♯₊::FlaggedVector{D,T,L}
    ρ::FlaggedVector{D,T,L}
end

"""
A dummy turn statistic to use for convenient type stability; returns null ptrs that should hopefully segfault if accessed
(we'd only access them accidentally in case of a bug!)
"""
function dummy_turn_statistic(::Tree{D,T,L}) where {D,T,L}
    dummy_flag = zero(UInt32)
    dummy_ptr = reinterpret(Ptr{T}, zero(UInt))
    GeneralizedTurnStatistic(
        FlaggedVector{D,T,L}(dummy_ptr, dummy_flag),
        FlaggedVector{D,T,L}(dummy_ptr, dummy_flag),
        FlaggedVector{D,T,L}(dummy_ptr, dummy_flag)
    )
end

function leaf_turn_statistic(tree::Tree{D,T,L}, ::Val{:generalized}, H, z::PhasePoint{D,T,L}) where {D,T,L}
    p♯ = calculate_p♯(tree, H, z)
    GeneralizedTurnStatistic(p♯, p♯, FlaggedVector{D,T,L}(z.p, 0x00000000))
end

function combine_turn_statistics(
    tree::Tree, ::TrajectoryNUTS,
    x::GeneralizedTurnStatistic{D,T,L},
    y::GeneralizedTurnStatistic{D,T,L}
) where {D,T,L}
    ρₓ = x.ρ
    ρʸ = y.ρ
    # @show x.p♯₊.flag, y.p♯₋.flag
    if ρₓ.flag == 0x00000000
        # We are in depth 1; position allocated as part of a PhasePoint, not as part of a turn statistic
        # Therefore we cannot free either of them, and must allocate a new vector to store the results in.
        ρ = undefined_Σρ(tree)
        # Additionally, x.p♯₊ == x.p♯₋ and y.p♯₊ == y.p♯₋, therefore we cannot free them.
    else # we are at a depth of 2 or greater.
        # We can therefore reuse a ρ; arbitrarily, we reuse xρ
        ρ = ρₓ
        # and free yρ as well as the two discarded p♯
        free_Σρ!(tree, ρʸ.flag)
        free_ρ♯!(tree, x.p♯₊.flag)
        free_ρ♯!(tree, y.p♯₋.flag)
    end
    @avx for l in eachindex(ρ)
        ρ[l] = ρₓ[l] + ρʸ[l]
    end
    # x.p♯₊.flag == x.p♯₋.flag || free_ρ♯!(tree, x.p♯₊.flag)
    # y.p♯₊.flag == y.p♯₋.flag || free_ρ♯!(tree, y.p♯₋.flag)
    GeneralizedTurnStatistic(x.p♯₋, y.p♯₊, ρ)
end

function is_turning(::TrajectoryNUTS, τ::GeneralizedTurnStatistic{D,T,L}) where {D,T,L}
    # Uses the generalized NUTS criterion from Betancourt (2017).
    @unpack p♯₋, p♯₊, ρ = τ
    # @argcheck p♯₋ ≢ p♯₊ "internal error: is_turning called on a leaf"
    d♯₋ = zero(T)
    d♯₊ = zero(T)
    @avx for d ∈ eachindex(ρ)
        ρᵈ = ρ[d]
        d♯₋ += ρᵈ * p♯₋[d]
        d♯₊ += ρᵈ * p♯₊[d]
    end
    # Current version always calculates both dot products; alternative:
    # dot(p♯₋, ρ) < 0 || dot(p♯₊, ρ) < 0
    # Calculating both in one loop is faster
    # ( 2D fma + 3D loads vs 2D fma + 4D loads; at most 2/4 ops/cycle
    #   can be loads, so the loads cannot keep up with the 2/4 fma/cycle
    #   and are the bottleneck ),
    # but it disallows conditional checking to skip second dot product.
    #
    # Which is faster on average depends on probability of turning
    # probability of turning will be fairly low for most models.
    (d♯₋ < zero(T)) | (d♯₊ < zero(T))
end

###
### leafs
###

function leaf(tree::Tree, H::Hamiltonian, trajectory::TrajectoryNUTS, z, is_initial)
    @unpack π₀, min_Δ, turn_statistic_configuration = trajectory
    # @show logdensity(H, z), π₀
    Δ = is_initial ? zero(π₀) : logdensity(H, z) - π₀
    isdiv = Δ < min_Δ
    v = leaf_acceptance_statistic(Δ, is_initial)
    if isdiv
        # @show H
        # @show z
        # @show typeof(z)
        (z, Δ, dummy_turn_statistic(tree)), v, true
    else
        τ = leaf_turn_statistic(tree, turn_statistic_configuration, H, z)
        (z, Δ, τ), v, false
    end
end

####
#### NUTS interface
####

"""
$(TYPEDEF)
Implementation of the “No-U-Turn Sampler” of Hoffman and Gelman (2014), including subsequent
developments, as described in Betancourt (2017).
# Fields
$(FIELDS)
"""
struct NUTS{S}
    "Maximum tree depth."
    max_depth::Int
    "Threshold for negative energy relative to starting point that indicates divergence."
    min_Δ::Float64
    """
    Turn statistic configuration. Currently only `Val(:generalized)` (the default) is
    supported.
    """
    turn_statistic_configuration::S
    function NUTS(; max_depth = DEFAULT_MAX_TREE_DEPTH, min_Δ = -1000.0,
                  turn_statistic_configuration::S = Val{:generalized}()) where {S}
        # @argcheck 0 < max_depth ≤ MAX_DIRECTIONS_DEPTH
        # @argcheck min_Δ < 0
        new{S}(Int(max_depth), Float64(min_Δ), turn_statistic_configuration)
    end
end

"""
$(TYPEDEF)
Diagnostic information for a single tree built with the No-U-turn sampler.
# Fields
Accessing fields directly is part of the API.
$(FIELDS)
"""
struct TreeStatisticsNUTS # now 32 bytes; fits in ymm register
    "Log density (negative energy)."
    π::Float64
    "Acceptance rate statistic."
    acceptance_rate::Float64
    "Reason for termination. See [`InvalidTree`](@ref) and [`REACHED_MAX_DEPTH`](@ref)."
    termination::InvalidTree
    "Depth of the tree."
    depth::Int32
    "Number of leapfrog steps evaluated."
    steps::Int32
    # "Directions for tree doubling (useful for debugging)."
    # directions::Directions
end

"""
$(SIGNATURES)
No-U-turn Hamiltonian Monte Carlo transition, using Hamiltonian `H`, starting at evaluated
log density position `Q`, using stepsize `ϵ`. Parameters of `algorithm` govern the details
of tree construction.
Return two values, the new evaluated log density position, and tree statistics.
"""
function sample_tree(rng, tree::Tree, algorithm::NUTS, H::Hamiltonian, z::PhasePoint, ϵ;
                          p = nothing, directions = rand(rng, Directions))
    if p === nothing
        rand_p!(rng, z.p, H.κ)
    else
        @unpack Q, flag = z
        z = PhasePoint(Q, p, flag)
    end
    @unpack max_depth, min_Δ, turn_statistic_configuration = algorithm
    trajectory = TrajectoryNUTS(logdensity(H, z), ϵ, min_Δ, turn_statistic_configuration)
    ζ, v, termination, depth = sample_trajectory(rng, tree, H, trajectory, z, max_depth, directions)
    tree_statistics = TreeStatisticsNUTS(logdensity(H, ζ), acceptance_rate(v), termination, depth, v.steps)#, directions)
    ζ, tree_statistics
end


text:
#
# broadcasting.jl --
#
# Broadcasting of arrays with optional element type conversion.
#

"""
    bcastcopy(A, [T=eltype(A),] dims...)

yields a new array of element type `T` and dimensions `dims` whose values are
given by `A` according to type conversion and broadcasting rules (like for the
`broadcast` method).  Compared to [`bcastlazy`](@ref), it is guaranteed that
the returned array does not share its contents with `A`.

Argument `A` can be a scalar value or an array.

See also [`bcastlazy`](@ref), [`bcastsize`](@ref).

"""
function bcastcopy(A, ::Type{T}, dims::Dims) where {T}
    B = Array{T}(undef, dims)
    B .= A # This expression will clash if dimensions are not compatible.
    return B
end

# Convert dimensions.
bcastcopy(A, ::Type{T}, dims::Tuple{Vararg{Integer}}) where {T} =
    bcastcopy(A, T, to_size(dims))
bcastcopy(A, ::Type{T}, dims::Integer...) where {T} =
    bcastcopy(A, T, dims)

# Guess element type.
bcastcopy(A, dims::Tuple{Vararg{Integer}}) =
    bcastcopy(A, eltype(A), dims)
bcastcopy(A, dims::Integer...) =
    bcastcopy(A, dims)

"""
    bcastlazy(A, [T=eltype(A),] dims...)

yields a *flat* array of type `T` and dimensions `dims` whose values are given
by `A` according to type conversion and broadcasting rules (see `broadcast`
method).  Compared to [`bcastcopy`](@ref), making a copy of `A` is avoided if
it is already an array with the correct type of elements and dimensions or if
it can be reshaped (by the `reshape` method) to the correct type and
dimensions.  This means that the result may share the same contents as `A`.
Argument `A` can be a scalar or an array with 1-based indices.  The result has
1-based indices and contiguous elements which is suitable for fast linear
indexing.

See also [`bcastcopy`](@ref), [`bcastsize`](@ref).

"""
function bcastlazy(A::DenseArray{T},
                   ::Type{T},
                   dims::Dims) where {T}
    has_standard_indexing(A) || throw_non_standard_indexing()
    Adims = size(A)
    Adims == dims && return A
    bcastsize(Adims, dims) == dims ||
        throw(DimensionMismatch("array has incompatible dimensions"))
    return (length(A) == prod(dims) ? reshape(A, dims) : bcastcopy(A, T, dims))
end

# Call `bcastcopy` if argument `A` can certainly not be returned.
bcastlazy(A, ::Type{T}, dims::Dims) where {T} =
    bcastcopy(A, T, dims)

# Convert dimensions.
bcastlazy(A, ::Type{T}, dims::Tuple{Vararg{Integer}}) where {T} =
    bcastlazy(A, T, to_size(dims))
bcastlazy(A, ::Type{T}, dims::Integer...) where {T} =
    bcastlazy(A, T, dims)

# Guess element type.
bcastlazy(A, dims::Tuple{Vararg{Integer}}) =
    bcastlazy(A, eltype(A), dims)
bcastlazy(A, dims::Integer...) =
    bcastlazy(A, dims)

"""
    bcastsize(size(A), size(B), ...) -> siz

yields the size `siz` of the array that would result from applying broadcasting
rules (see `broadcast` method) to arguments `A`, `B`, etc.  The result is a
tuple of integers (of type `Int`).  Call [`check_size`](@ref) if you want to
also make sure that the result is a list of valid dimensions.

The method can also be applied to a single dimension:

    bcastsize(a, b) -> c

to yield the dimension `c` gievn by broadcasting dimensions `a` and `b`
throwing an exception if dimensions are not compatible according to
broadcasting rules.  This is the same as `Base.Broadcasting._bcs1` but it takes
care of converting to `Int`.

See also [`standard_size`](@ref), [`check_size`](@ref), [`bcastcopy`](@ref),
[`bcastlazy`](@ref).

"""
bcastsize(::Tuple{}) = ()
bcastsize(a::Tuple{Vararg{Integer}}) = to_size(a)
bcastsize(a::Tuple{Vararg{Integer}}, b::Tuple{Vararg{Integer}}, args...) =
    bcastsize(bcastsize(a, b), args...)

# Use a recursion to build the dimension list from two lists, if code is
# inlined (for a few number of dimensions), it should be very fast.
bcastsize(::Tuple{}, ::Tuple{}) = ()
bcastsize(::Tuple{}, b::Tuple{Vararg{Integer}}) = bcastsize(b)
bcastsize(a::Tuple{Vararg{Integer}}, ::Tuple{}) = bcastsize(a)
bcastsize(a::Tuple{Vararg{Integer}}, b::Tuple{Vararg{Integer}}) =
    (bcastsize(a[1], b[1]), bcastsize(Base.tail(a), Base.tail(b))...)

bcastsize(a::Integer, b::Integer) = bcastsize(to_int(a), to_int(b))
bcastsize(a::Int, b::Int) =
    (a == b || b == 1 ? a : a == 1 ? b :
    throw(DimensionMismatch("arrays could not be broadcast to a common size")))

ext: jl


text:
using Test
using CLIMA
using CLIMA.HydrostaticBoussinesq
using CLIMA.SimpleBox
using CLIMA.GenericCallbacks
using CLIMA.ODESolvers
using CLIMA.Mesh.Filters
using CLIMA.PlanetParameters
using CLIMA.VariableTemplates
using CLIMA.Mesh.Grids: polynomialorder
import CLIMA.DGmethods: vars_state

function config_simple_box(FT, N, resolution, dimensions)
  prob = HomogeneousBox{FT}(dimensions...)

  cʰ = sqrt(grav * prob.H) # m/s
  model = HydrostaticBoussinesqModel{FT}(prob, cʰ = cʰ)

  config = CLIMA.Ocean_BoxGCM_Configuration("homogeneous_box",
                                            N, resolution, model)

  return config
end

function main(;imex::Bool = false)
  CLIMA.init()

  FT = Float64

  # DG polynomial order
  N = Int(4)

  # Domain resolution and size
  Nˣ = Int(20)
  Nʸ = Int(20)
  Nᶻ = Int(20)
  resolution = (Nˣ, Nʸ, Nᶻ)

  Lˣ = 4e6   # m
  Lʸ = 4e6   # m
  H  = 400   # m
  dimensions = (Lˣ, Lʸ, H)

  timestart = FT(0)    # s
  timeend   = FT(3600) # s

  if imex
    solver_type = CLIMA.IMEXSolverType(linear_model=LinearHBModel)
    Nᶻ = Int(400)
  else
    solver_type = CLIMA.ExplicitSolverType(solver_method=LSRK144NiegemannDiehlBusch)
  end

  driver_config = config_simple_box(FT, N, resolution, dimensions)

  grid = driver_config.grid
  vert_filter = CutoffFilter(grid, polynomialorder(grid)-1)
  exp_filter  = ExponentialFilter(grid, 1, 8)
  modeldata = (vert_filter = vert_filter, exp_filter=exp_filter)

  solver_config = CLIMA.setup_solver(timestart, timeend, driver_config, forcecpu=true,
                                     ode_solver_type=solver_type, modeldata=modeldata)

  result = CLIMA.invoke!(solver_config)

  maxQ = Vars{vars_state(driver_config.bl, FT)}(maximum(solver_config.Q, dims=(1,3)))
  minQ = Vars{vars_state(driver_config.bl, FT)}(minimum(solver_config.Q, dims=(1,3)))

  @test maxQ.θ ≈ minQ.θ
end

main(imex=false)
main(imex=true)

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
using DashBootstrapComponents, DashHtmlComponents

badges = html_span([
    dbc_badge("Primary", href="#", color="primary", className="mr-1"),
    dbc_badge("Secondary", href="#", color="secondary", className="mr-1"),
    dbc_badge("Success", href="#", color="success", className="mr-1"),
    dbc_badge("Warning", href="#", color="warning", className="mr-1"),
    dbc_badge("Danger", href="#", color="danger", className="mr-1"),
    dbc_badge("Info", href="#", color="info", className="mr-1"),
    dbc_badge("Light", href="#", color="light", className="mr-1"),
    dbc_badge("Dark", href="#", color="dark"),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
###############################################################################
#
# Choose text font
#
###############################################################################
"""
    use_sans()

Set default font to sans-serif
"""
function use_sans()
    rc("font"    , family="sans-serif", serif=["sans-serif"]);
    rc("mathtext", fontset="dejavusans");

    return nothing
end

"""
    use_serif()

Set default font to serif
"""
function use_serif()
    rc("font"    , family="serif", serif=["Palatino"]);
    rc("mathtext", fontset="dejavuserif");

    return nothing
end


text:
__precompile__(true)
module DataViewer

#-------------------------------------------------------------------------------
# Dependencies
#-------------------------------------------------------------------------------

using DataFrames

#-------------------------------------------------------------------------------
# Exported methods and types
#-------------------------------------------------------------------------------

export viewdata,
       writehtml

#-------------------------------------------------------------------------------
# Load methods and types
#-------------------------------------------------------------------------------

include("WriteHTML.jl")
include("ViewData.jl")

end


text:
###### This file is part of the MendelImpute.jl package.
###### It contains relevant code to perform chromosome painting and computing
###### computation of poplation admixtures.

"""
    composition(sample_phase::HaplotypeMosaicPair, panelID::Vector{String},
        refID_to_population::Dict{String, String}, [populations::Vector{String}])

Computes a sample's chromosome composition based on phase information. This
function is used for easier plotting a person's admixed proportions.

# Inputs
- `sample_phase`: A `HaplotypeMosaicPair` storing phase information for a
    sample, includes haplotype start position and haplotype label.
- `panelID`: Sample ID's in the reference haplotype panel
- `refID_to_population`: A dictionary mapping each ID in the haplotype
    reference panel to its population origin.

# Optional inputs
- `populations`: A unique list of populations present in `refID_to_population`

# Output
- `composition`: A list of percentages where `composition[i]` equals the
    sample's ancestry (in %) from `populations[i]`
"""
function composition(
    sample_phase::HaplotypeMosaicPair,
    panelID::Vector{String},
    refID_to_population::Dict{String, String};
    populations::Vector{String} = unique_populations(refID_to_population)
    )
    snps = sample_phase.strand1.length
    @assert snps == sample_phase.strand2.length "strands have different length!"

    # return elements
    composition = zeros(length(populations))

    # strand1
    l1 = length(sample_phase.strand1.haplotypelabel)
    for i in 1:l1
        # get complete haplotype index and range of haplotype
        h1 = sample_phase.strand1.haplotypelabel[i]
        cur_range = (i == l1 ?
            (sample_phase.strand1.start[i]:sample_phase.strand1.length) :
            (sample_phase.strand1.start[i]:(sample_phase.strand1.start[i+1]
            - 1)))

        # find current haplotype's population origin
        ref_id = panelID[div(h1, 2, RoundUp)]
        population = refID_to_population[ref_id]

        # update composition
        idx = findfirst(x -> x == population, populations)
        composition[idx] += length(cur_range)
    end

    # strand2
    l2 = length(sample_phase.strand2.haplotypelabel)
    for i in 1:l2
        # get complete haplotype index and range of haplotype
        h2 = sample_phase.strand2.haplotypelabel[i]
        cur_range = (i == l2 ?
            (sample_phase.strand2.start[i]:sample_phase.strand2.length) :
            (sample_phase.strand2.start[i]:(sample_phase.strand2.start[i+1]
            - 1)))

        # find current haplotype's population origin
        ref_id = panelID[div(h2, 2, RoundUp)]
        population = refID_to_population[ref_id]

        # update composition
        idx = findfirst(x -> x == population, populations)
        composition[idx] += length(cur_range)
    end

    sum(composition) == 2snps || error("composition should sum to number of snps")

    return composition ./ 2snps
end

"""
    unique_populations(x::Dict{String, String})

Computes the unique list of populations, preserving order. `x` is a `Dict`
where each sample is a key and populations are values.
"""
function unique_populations(x::Dict{String, String})
    populations = String[]
    for (key, val) in x
        val ∉ populations && push!(populations, val)
    end
    return populations
end

"""
    paint(sample_phase::HaplotypeMosaicPair, panelID::Vector{String},
        refID_to_population::Dict{String, String}, populations::Vector{String})

Converts a person's phased haplotype lengths into segments of percentages. This
function is used for easier plotting a "painted chromosome".

# Inputs
- `sample_phase`: A `HaplotypeMosaicPair` storing phase information for a
    sample, includes haplotype start position and haplotype label.
- `panelID`: Sample ID's in the reference haplotype panel
- `refID_to_population`: A dictionary mapping each ID in the haplotype
    reference panel to its population origin.

# Optional inputs
- `populations`: A unique list of populations present in `refID_to_population`

# Output
- `composition`: A list of percentages where `composition[i]` equals the
    sample's ancestry (in %) from `populations[i]`
"""
function paint(
    sample_phase::HaplotypeMosaicPair,
    panelID::Vector{String},
    refID_to_population::Dict{String, String};
    populations::Vector{String} = unique_populations(refID_to_population)
    )
    snps = sample_phase.strand1.length
    @assert snps == sample_phase.strand2.length "strands have different length!"

    # length of strand 1 and 2
    l1 = length(sample_phase.strand1.haplotypelabel)
    l2 = length(sample_phase.strand2.haplotypelabel)

    # return elements
    s1_composition = (zeros(l1), Vector{String}(undef, l1))
    s2_composition = (zeros(l2), Vector{String}(undef, l2))

    # strand1
    for i in 1:l1
        # get complete haplotype index and range of haplotype
        h1 = sample_phase.strand1.haplotypelabel[i]
        cur_range = (i == l1 ?
            (sample_phase.strand1.start[i]:sample_phase.strand1.length) :
            (sample_phase.strand1.start[i]:(sample_phase.strand1.start[i+1]
            - 1)))

        # find current haplotype's population origin
        ref_id = panelID[div(h1, 2, RoundUp)]
        population = refID_to_population[ref_id]

        # update composition
        s1_composition[1][i] = length(cur_range) / snps
        s1_composition[2][i] = population
    end

    # strand2
    for i in 1:l2
        # get complete haplotype index and range of haplotype
        h2 = sample_phase.strand2.haplotypelabel[i]
        cur_range = (i == l2 ?
            (sample_phase.strand2.start[i]:sample_phase.strand2.length) :
            (sample_phase.strand2.start[i]:(sample_phase.strand2.start[i+1]
            - 1)))

        # find current haplotype's population origin
        ref_id = panelID[div(h2, 2, RoundUp)]
        population = refID_to_population[ref_id]

        # update composition
        s2_composition[1][i] = length(cur_range) / snps
        s2_composition[2][i] = population
    end

    sum(s1_composition[1]) ≈ 1.0 || error("strand 1 compositions should sum to 1")
    sum(s2_composition[1]) ≈ 1.0 || error("strand 1 compositions should sum to 1")

    return s1_composition, s2_composition
end


text:
"""
    abstract type TimeIntegrator

Abstract type for ODE-Integrators.
They should store the update as internal state to
avoid costly reallocations.

# Implementation
For new integrators, the method `step` needs to be
overwritten.
This method should modify the dofs stored in the grid
in-place.
"""
abstract type TimeIntegrator end

"""
    ExplicitEuler(grid::Grid)

Return Euler time-integrator for `grid`.
First order accurate.
"""
struct ExplicitEuler
    dofsupdate::Array{Float64,2}

    function ExplicitEuler(grid::Grid)
        dofsupdate = similar(grid.dofs)
        new(dofsupdate)
    end
end

"""
    step(f, integrator::ExplicitEuler, grid, dt)

Performs an update with the explicit euler method
on `grid` and timestepsize `dt`.
"""
function step(f, integrator::ExplicitEuler, grid, dt)
    integrator.dofsupdate .= 0.0
    f(integrator.dofsupdate, grid.dofs, grid.time)
    grid.dofs .+= dt .* integrator.dofsupdate
    @info "Update has" norm(integrator.dofsupdate)
end

"""
    SSPRK2(grid::Grid)

Return SSPRK2 time-integrator for `grid`.
Two-stage strong-stability preserving Runge Kutta method.
Second order accurate.
"""
struct SSPRK2
    stage1::Array{Float64,2}
    stage2::Array{Float64,2}

    function SSPRK2(grid::Grid)
        stage1 = similar(grid.dofs)
        stage2 = similar(grid.dofs)
        new(stage1, stage2)

    end
end

"""
    step(f, integrator::SSPRK2, grid, dt)

Performs an update with the SSPRK2 method
on `grid` and timestepsize `dt`.
"""

function step(f, integrator::SSPRK2, grid, dt)
    integrator.stage1 .= 0.0
    integrator.stage2 .= 0.0

    f(integrator.stage1, grid.dofs, grid.time)
    integrator.stage1 .= grid.dofs .+ dt .* integrator.stage1

    f(integrator.stage2, integrator.stage1, grid.time + dt)

    grid.dofs .= 0.5 .* (grid.dofs .+ integrator.stage1 .+
        dt .* integrator.stage2)
end

"""
    SSPRK3(grid::Grid)

Return SSPRK3 time-integrator for `grid`.
Three-stage strong-stability preserving Runge Kutta method.
Third order accurate.
"""
struct SSPRK3
    stage1::Array{Float64,2}
    stage2::Array{Float64,2}
    stage3::Array{Float64,2}

    function SSPRK3(grid::Grid)
        stage1 = similar(grid.dofs)
        stage2 = similar(grid.dofs)
        stage3 = similar(grid.dofs)
        new(stage1, stage2, stage3)

    end
end

"""
    step(f, integrator::SSPRK3, grid, dt)

Performs an update with the SSPRK3 method
on `grid` and timestepsize `dt`.
"""
function step(f, integrator::SSPRK3, grid, dt)
    integrator.stage1 .= 0.0
    integrator.stage2 .= 0.0
    integrator.stage3 .= 0.0

    f(integrator.stage1, grid.dofs, grid.time)
    integrator.stage1 .= grid.dofs .+ dt .* integrator.stage1

    f(integrator.stage2, integrator.stage1, grid.time + dt)
    integrator.stage2 .= 0.75 .* grid.dofs .+ 0.25 .* (
        integrator.stage1 + dt *integrator.stage2)

    f(integrator.stage3, integrator.stage2, grid.time + 0.5 * dt)

    grid.dofs .= 1.0/3.0 .* grid.dofs .+ 2.0/3.0 * (
        integrator.stage2 .+ dt .* integrator.stage3)
end

"""
    make_timeintegrator(config::Configuration, grid::Grid)

Returns time integrator from `config` for `grid`.
"""
function make_timeintegrator(config::Configuration, grid::Grid)
    if config.timeintegrator_name == "Euler"
        return ExplicitEuler(grid)
    elseif config.timeintegrator_name == "SSPRK2"
        return SSPRK2(grid)
    elseif config.timeintegrator_name == "SSPRK3"
        return SSPRK3(grid)
    else
        error(string("Unknown timeintegrator name: ", config.timeintegrator_name))
    end
end


text:
using Printf
using LinearAlgebra
using Random

include("../common/PWGrid_v01.jl")
include("../common/ortho_gram_schmidt.jl")
include("../common/wrappers_fft.jl")
include("../common/gen_lattice_pwscf.jl")
include("EnergiesT.jl")
include("PotentialsT.jl")
include("gen_dr.jl")
include("op_K.jl")
include("op_V_ps_loc.jl")
include("op_H.jl")
include("calc_rho.jl")
include("calc_grad.jl")
include("calc_Energies.jl")
include("KS_solve_Emin_sd.jl")
include("KS_solve_Emin_cg.jl")
include("Poisson_solve.jl")
include("LDA_VWN.jl")
include("Kprec.jl")

include("structure_factor.jl")
include("gen_rho.jl")
include("gen_dr.jl")
include("calc_ewald.jl")

function test_main( ns1::Int64,ns2::Int64,ns3::Int64 )

    Ns = [ns1,ns2,ns3]
    LatVecs = gen_lattice_sc(16.0)

    pw = PWGrid( Ns, LatVecs )

    Npoints = pw.Npoints
    Ω = pw.Ω
    r = pw.r
    G = pw.G
    G2 = pw.G2

    Xpos = reshape( [0.0, 0.0, 0.0], (3,1) )

    Sf = structure_factor( Xpos, G )

    E_nn = calc_ewald( pw, Xpos, Sf )

    Vg = zeros(ComplexF64,Npoints)
    prefactor = -4*pi/Ω
    for ig=2:Npoints
        Vg[ig] = prefactor/G2[ig]
    end
    V_ionic = real( G_to_R(Ns, Vg .* Sf) ) * Npoints

    #
    Nstates = 1
    Focc = [1.0]
    psi, Energies, Potentials = KS_solve_Emin_cg( pw, V_ionic, Focc, Nstates, NiterMax=1000 )

    #
    Y = ortho_gram_schmidt(psi)
    mu = Y' * op_H( pw, Potentials, Y )
    evals, evecs = eigen(mu)
    psi = Y*evecs

    for st = 1:Nstates
        @printf("State # %d, Energy = %f\n", st, real(evals[st]))
    end

    @printf("E_nn    = %18.10f\n", E_nn)
    @printf("E total = %18.10f\n", E_nn + Energies.Total)

end

@time test_main( 64, 64, 64 )


text:
using FinancialModels, ARCHModels, Polynomials, JLD
using HypothesisTests: ADFTest, ApproximateTwoSampleKSTest
using Statistics: mean, quantile, std
using StatsBase: kurtosis
using GLM: lm, coef
include("NMTA.jl")
struct Moments # Moments of log-returns
    μ::Float64 # Mean
    σ::Float64 # Standard deviation
    κ::Float64 # Kurtosis
    ks::Float64 # Kolmogorov-Smirnov test statistic for the difference between distributions
    hurst::Float64 # Hurst exponent: hurst < 0.5 => mean reverting; hurst == 0.5 => random walk; hurst > 0.5 => momentum
    gph::Float64 # GPH estimator representing long-range dependence
    adf::Float64 # ADF statistic representing random walk property of returns
    garch::Float64 # GARCH paramaters representing short-range dependence
    hill::Float64 # Hill estimator
    function Moments(log_returns1::Vector{Float64}, log_returns2::Vector{Float64})
        μ = mean(log_returns1); σ = std(log_returns1); κ = kurtosis(log_returns1)
        ks = ApproximateTwoSampleKSTest(log_returns1, log_returns2).δ
        hurst = HurstExponent(log_returns1)
        gph = GPH(abs.(log_returns1))
        adf = ADFTest(log_returns1, :none, 0).stat
        garch = sum(coef(ARCHModels.fit(GARCH{1, 1}, log_returns1))[2:3])
        hill = HillEstimator(log_returns1[findall(x -> (x >= quantile(log_returns1, 0.95)) && (x > 0), log_returns1)], 50)
        new(μ, σ, κ, ks, hurst, gph, adf, garch, hill)
    end
end
function HurstExponent(x, d = 100)
    function Divisors(n, n₀)
        return filter(x -> mod(n, x) == 0, n₀:floor(n/2))
    end
    N = length(x)
    if mod(N, 2) != 0 x = push!(x, (x[N - 1] + x[N]) / 2); N += 1 end
    N₁ = N₀ = min(floor(0.99 * N), N-1); dv = Divisors(N₁, d)
    for i in (N₀ + 1):N
        dw = Divisors(i, d)
        if length(dw) > length(dv) N₁ = i; dv = copy(dw) end
    end
    OptN = Int(N₁); d = dv
    x = x[1:OptN]
    RSempirical = map(i -> RS(x, i), d)
    return coeffs(Polynomials.fit(Polynomial, log10.(d), log10.(RSempirical), 1))[2] # Hurst is slope of log-log linear fit
end
function RS(z, n)
    y = reshape(z, (Int(n), Int(length(z) / n)))
    μ = mean(y, dims = 1)
    σ = std(y, dims = 1)
    temp = cumsum(y .- μ, dims = 1)
    return mean((maximum(temp, dims = 1) - minimum(temp, dims = 1)) / σ)
end
function GPH(x, bandwidthExponent = 0.5)
    n = length(x); g = Int(trunc(n^bandwidthExponent))
    j = 1:g; kk = 1:(n - 1)
    w = 2 .* π .* j ./ n # x .-= mean(x)
    σ = sum(x .^ 2) / n
    Σ = map(k -> sum(x[1:(n - k)] .* x[(1 + k):n]) / n, kk)
    periodogram = map(i -> σ + 2 * sum(Σ .* cos.(w[i] .* kk)), j)
    indeces = j[findall(x -> x > 0, periodogram)]
    x_reg = 2 .* log.(2 .* sin.(w[indeces] ./ 2)); y_reg = log.(periodogram[indeces] ./ (2 * π))
    regression = lm(hcat(ones(length(x_reg)), x_reg), y_reg)
    return abs(coef(regression)[2])
end
function HillEstimator(x, iterations)
    N = length(x)
    logx = log.(x)
    L = minimum(x); R = maximum(x)
    α = 1 / ((sum(logx) / N) - log(L))
    for i in 1:iterations
        C = (log(L) * (L^(-α)) - log(R) * (R^(-α))) / ((L^(-α)) - (R^(-α)))
        D = (R^α * L^α * (log(L) - log(R))^2) / (L^α - R^α)^2
        α = α * (1 + (α * (sum(logx) / N) - α * C - 1) / (α^2 * D - 1))
    end
    return α
end
function MovingBlockBootstrap(mid_price::Vector{Float64}, iterations::Int64 = 1000, window_size::Int64 = 2000)
    log_returns = diff(log.(mid_price))
    bootstrap_moments = Vector{Vector{Float64}}()
    Threads.@threads for i in 1:iterations
        indeces = rand(1:(length(log_returns) - window_size + 1), Int(ceil(length(log_returns) / window_size)))
        bootstrap_returns = Vector{Float64}()
        for index in indeces
            append!(bootstrap_returns, log_returns[index:(index  + window_size - 1)])
        end
        moments = Moments(bootstrap_returns[1:length(log_returns)], log_returns)
        push!(bootstrap_moments, [moments.μ, moments.σ, moments.κ, moments.ks, moments.hurst, moments.gph, moments.adf, moments.garch, moments.hill])
    end
    W = inv(cov(bootstrap_moments))
    save("data/W.jld", "W", W)
end
function WSSE(parameters::JacobLealParameters, replications::Int64, W::Array{Float64, 2}, empirical_moments::Moments, empirical_log_returns::Vector{Float64})
    errormatrix = fill(0.0, (replications, 9))
    Threads.@threads for i in 1:replications
        mid_price = prices(Simulate(parameters, options = Options(time = 10000)))
        log_returns = diff(log.(mid_price))
        simulated_moments = Moments(log_returns, empirical_log_returns)
        errormatrix[i, :] = [simulated_moments.μ-empirical_moments.μ simulated_moments.σ-empirical_moments.σ simulated_moments.κ-empirical_moments.κ simulated_moments.ks-empirical_moments.ks simulated_moments.hurst-empirical_moments.hurst simulated_moments.gph-empirical_moments.gph simulated_moments.adf-empirical_moments.adf simulated_moments.garch-empirical_moments.garch simulated_moments.hill-empirical_moments.hill]
    end
    errors = mean(errormatrix, dims = 1)
    return (errors * W * transpose(errors))[1]
end

mid_price = CSV.File("data/L1LOB.csv", missingstring = "missing", ignoreemptylines = true, select = [:MicroPrice], skipto = 20000, limit = 20000) |> Tables.matrix |> vec |> y -> filter(z -> !ismissing(z), y)
empirical_log_returns = diff(log.(mid_price))
empirical_moments = Moments(empirical_log_returns, empirical_log_returns)
MovingBlockBootstrap(mid_price, 10, 2000)
W = load("W.jld")["W"]
initial_solution = JacobLealParameters()
initial_solution = [initial_solution[1], initial_solution[2], initial_solution[3], initial_solution[4], initial_solution[5], initial_solution[6], initial_solution[7],
initial_solution[8], initial_solution[9], initial_solution[10], initial_solution[11], initial_solution[12],  initial_solution[13]]
ta_rounds = [12, 10, 8, 6]
f_reltol = collect(range(10, 0,  sum(ta_rounds)))
neldermeadstate = nothing
objective = NonDifferentiable(x -> WSSE(x -> JacobLealParameters(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10], x[11], x[12],  x[13]), 10, W, empirical_moments, empirical_log_returns), initial_solution)
optimization_options = Options(show_trace = true, store_trace = true, trace_simplex = true, extended_trace = true, iterations = 30, ξ = 0.15, ta_rounds = ta_rounds, f_reltol = f_reltol)
result = !isnothing(neldermeadstate) ? Optimize(objective, initial_solution, optimization_options, neldermead_state) : Optimize(objective, initial_solution, optimization_options)
save("OptimizationResult.jld", "result", result)
#---------------------------------------------------------------------------------------------------

max_stars_repo_name: IvanJericevich/FinancialModels

max_issues_repo_name: IvanJericevich/FinancialModels

max_forks_repo_name: IvanJericevich/FinancialModels


text:
module MixedAnova
using DataFrames, Statistics, StatsBase, LinearAlgebra, Distributions, Reexport, Printf
@reexport using StatsModels
import StatsBase: fit!, fit, dof, dof_residual, deviance, nobs, coefnames
import StatsModels: TableStatisticalModel, TableRegressionModel, vectorize, kron_insideout,
                    ModelFrame, ModelMatrix, response, columntable, asgn, collect_matrix_terms
import LinearAlgebra.BlasReal
import Tables.istable
import Base: show

export
    # wrappers
    AnovaResult,

    # anova functions
    anova, anova_lm, anova_lme, anova_glm,

    # GoodnessOfFit
    GoodnessOfFit, FTest, LikelihoodRatioTest, LRT, canonicalgoodnessoffit,

    # Others
    lme, glme, lfe, nestedmodels, calcdof, formula, #getterms,
    teststat, pval, anova_test, anova_type,

    # init
    glm_init, mm_init, fem_init

# Test
abstract type GoodnessOfFit end
struct FTest <: GoodnessOfFit end
struct LikelihoodRatioTest <: GoodnessOfFit end
const LRT = LikelihoodRatioTest

# Wrapper for ANOVA
"""
    AnovaResult{M, T, N}

Returned object of `anova`. \n
`M` is a subtype of `Tuple` if multiple models are provided; otherwise, a typeof model. \n
`T` is a subtype of `GoodnessOfFit`; either `FTest` or `LRT`. \n
`N` is the length of parameters.

## Fields

* `model`: full model or tuple of tested models.
* `type`: type of `anova`.
* `dof`: degree of freedoms of models or factors.
* `deviance`: deviance(s) for calculating test statistics. See `deviance` for more details.
* `teststat`: value(s) of test statiscics.
* `pval`: p-value(s) of test statiscics.
* `tests`: `NamedTuple` contained extra statistics.
"""
struct AnovaResult{M, T, N}
    model::M
    type::Int
    dof::NTuple{N, Int}
    deviance::NTuple{N, Float64}
    teststat::NTuple{N, Float64}
    pval::NTuple{N, Float64}
    tests::NamedTuple
end

AnovaResult{T}(model::M,
                type::Int,
                dof::NTuple{N, Int},
                deviance::NTuple{N, Float64},
                teststat::NTuple{N, Float64},
                pval::NTuple{N, Float64},
                tests::NamedTuple) where {M, N, T <: GoodnessOfFit} =
    AnovaResult{M, T, N}(model, type, dof, deviance, teststat, pval, tests)

include("api.jl")
include("termIO.jl")
const path = @__DIR__()

function glm_init()
    include(joinpath(path, "GLM", "anova.jl"))
    include(joinpath(path, "GLM", "io.jl"))
    include(joinpath(path, "GLM", "fit.jl"))
    return
end

function mm_init()
    include(joinpath(path, "MixedModels", "anova.jl"))
    include(joinpath(path, "MixedModels", "io.jl"))
    include(joinpath(path, "MixedModels", "fit.jl"))
    return
end

function fem_init()
    include(joinpath(path, "FixedEffectModels", "anova.jl"))
    include(joinpath(path, "FixedEffectModels", "io.jl"))
    include(joinpath(path, "FixedEffectModels", "fit.jl"))
    return
end

# Appendix
"""
# Appendix I: Type of sum of squares

For exmaple, a two-way ANOVA with factor A and B \n
Type I  : SS(A) -> SS(B | A) -> SS(AB | A,B)  \n
Type II : SS(A | B) -> SS(B | A) -> SS(AB | A,B) \n
Type III: SS(A | A,AB) -> SS(B | AB,B) -> SS(AB | A,B) # equivalent to linear regression  \n

# Appendix II: Examples for linear mixed-effect model
Balanced design: for each level of random effect, there is a fixed number of observations of each within variable
```
16×4 DataFrame
│ Row │ pa   │ time │ BP    │ drug │
│     │ Cat… │ Cat… │ Int64 │ Cat… │
├─────┼──────┼──────┼───────┼──────┤
│ 1   │ 1    │ am   │ 110   │ ACEI │
│ 2   │ 1    │ am   │ 112   │ ACEI │
│ 3   │ 1    │ pm   │ 130   │ ACEI │
│ 4   │ 1    │ pm   │ 135   │ ACEI │
│ 5   │ 2    │ am   │ 120   │ ARB  │
│ 6   │ 2    │ am   │ 127   │ ARB  │
│ 7   │ 2    │ pm   │ 125   │ ARB  │
│ 8   │ 2    │ pm   │ 121   │ ARB  │
│ 9   │ 3    │ am   │ 140   │ ACEI │
│ 10  │ 3    │ am   │ 141   │ ACEI │
│ 11  │ 3    │ pm   │ 165   │ ACEI │
│ 12  │ 3    │ pm   │ 152   │ ACEI │
│ 13  │ 4    │ am   │ 130   │ ARB  │
│ 14  │ 4    │ am   │ 124   │ ARB  │
│ 15  │ 4    │ pm   │ 145   │ ARB  │
│ 16  │ 4    │ pm   │ 151   │ ARB  │
```
BP ~ time * drug + (1|pa), time is within-subjects, drug is between-subjects, pa is random effect.

"""
appendix
end


text:
# This file is part of the TaylorSeries.jl Julia package, MIT license
#
# Luis Benet & David P. Sanders
# UNAM
#
# MIT Expat license
#

## Evaluating ##
"""
    evaluate(a, [dx])

Evaluate a `Taylor1` polynomial using Horner's rule (hand coded). If `dx` is
ommitted, its value is considered as zero. Note that the syntax `a(dx)` is
equivalent to `evaluate(a,dx)`, and `a()` is equivalent to `evaluate(a)`.
"""
function evaluate(a::Taylor1{T}, dx::T) where {T<:Number}
    @inbounds suma = a[end]
    @inbounds for k in a.order-1:-1:0
        suma = suma*dx + a[k]
    end
    suma
end
function evaluate(a::Taylor1{T}, dx::S) where {T<:Number,S<:Number}
    suma = a[end]*one(dx)
    @inbounds for k in a.order-1:-1:0
        suma = suma*dx + a[k]
    end
    suma
end
evaluate(a::Taylor1{T}) where {T<:Number} = a[0]

"""
    evaluate(x, δt)

Evaluates each element of `x::AbstractArray{Taylor1{T}}`,
representing the dependent variables of an ODE, at *time* δt. Note that the
syntax `x(δt)` is equivalent to `evaluate(x, δt)`, and `x()`
is equivalent to `evaluate(x)`.
"""
evaluate(x::AbstractArray{Taylor1{T}}, δt::S) where
    {T<:Number,S<:Number} = evaluate.(x, δt)
evaluate(a::AbstractArray{Taylor1{T}}) where {T<:Number} = evaluate.(a, zero(T))

"""
    evaluate!(x, δt, x0)

Evaluates each element of `x::AbstractArray{Taylor1{T}}`,
representing the Taylor expansion for the dependent variables
of an ODE at *time* `δt`. It updates the vector `x0` with the
computed values.
"""
function evaluate!(x::AbstractArray{Taylor1{T}}, δt::S,
        x0::AbstractArray{T}) where {T<:Number,S<:Number}

    @inbounds for i in eachindex(x, x0)
        x0[i] = evaluate( x[i], δt )
    end
    nothing
end

"""
    evaluate(a, x)

Substitute `x::Taylor1` as independent variable in a `a::Taylor1` polynomial.
Note that the syntax `a(x)` is equivalent to `evaluate(a, x)`.
"""
evaluate(a::Taylor1{T}, x::Taylor1{S}) where {T<:Number,S<:Number} =
    evaluate(promote(a,x)...)

function evaluate(a::Taylor1{T}, x::Taylor1{T}) where {T<:Number}
    if a.order != x.order
        a, x = fixorder(a, x)
    end
    @inbounds suma = a[end]*one(x)
    @inbounds for k = a.order-1:-1:0
        suma = suma*x + a[k]
    end
    suma
end

function evaluate(a::Taylor1{Taylor1{T}}, x::Taylor1{T}) where {T<:Number}
    @inbounds suma = a[end]*one(x)
    @inbounds for k = a.order-1:-1:0
        suma = suma*x + a[k]
    end
    suma
end
function evaluate(a::Taylor1{T}, x::Taylor1{Taylor1{T}}) where {T<:Number}
    @inbounds suma = a[end]*one(x)
    @inbounds for k = a.order-1:-1:0
        suma = suma*x + a[k]
    end
    suma
end

evaluate(p::Taylor1{T}, x::Array{S}) where {T<:Number,S<:Number} =
    evaluate.([p], x)

#function-like behavior for Taylor1
(p::Taylor1)(x) = evaluate(p, x)
(p::Taylor1)() = evaluate(p)

#function-like behavior for Vector{Taylor1}
if VERSION >= v"1.3"
    (p::AbstractArray{Taylor1{T}})(x) where {T<:Number} = evaluate.(p, x)
    (p::AbstractArray{Taylor1{T}})() where {T<:Number} = evaluate.(p)
else
    (p::Array{Taylor1{T}})(x) where {T<:Number} = evaluate.(p, x)
    (p::SubArray{Taylor1{T}})(x) where {T<:Number} = evaluate.(p, x)
    (p::Array{Taylor1{T}})() where {T<:Number} = evaluate.(p)
    (p::SubArray{Taylor1{T}})() where {T<:Number} = evaluate.(p)
end

## Evaluation of multivariable
function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{T,1},
        x0::AbstractArray{T}) where {T<:Number}

    # @assert length(x) == length(x0)
    @inbounds for i in eachindex(x, x0)
        x0[i] = evaluate( x[i], δx )
    end
    nothing
end

function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{Taylor1{T},1},
        x0::AbstractArray{Taylor1{T}}) where {T<:NumberNotSeriesN}

    # @assert length(x) == length(x0)
    @inbounds for i in eachindex(x, x0)
        x0[i] = evaluate( x[i], δx )
    end
    nothing
end

function evaluate!(x::AbstractArray{TaylorN{T}}, δx::Array{TaylorN{T},1},
        x0::AbstractArray{TaylorN{T}}; sorting::Bool=true) where {T<:NumberNotSeriesN}

    @inbounds for i in eachindex(x, x0)
        x0[i] = _evaluate( x[i], δx, Val(sorting) )
    end
    nothing
end

function evaluate!(x::AbstractArray{TaylorN{T}}, δt::T,
        x0::AbstractArray{TaylorN{T}}) where {T<:Number}

    # @assert length(x) == length(x0)
    @inbounds for i in eachindex(x, x0)
        x0[i] = evaluate( x[i], δt )
    end
    nothing
end

"""
    evaluate(a, [vals])

Evaluate a `HomogeneousPolynomial` polynomial at `vals`. If `vals` is ommitted,
it's evaluated at zero. Note that the syntax `a(vals)` is equivalent to
`evaluate(a, vals)`; and `a()` is equivalent to `evaluate(a)`.
"""
function evaluate(a::HomogeneousPolynomial{T}, vals::NTuple{N,S} ) where
        {T<:Number,S<:Number,N}

    @assert N == get_numvars()

    return _evaluate(a, vals)
end

function _evaluate(a::HomogeneousPolynomial{T}, vals::NTuple{N,S} ) where
        {T<:Number,S<:Number,N}

    ct = coeff_table[a.order+1]
    suma = zero(a[1])*vals[1]

    for (i,a_coeff) in enumerate(a.coeffs)
        iszero(a_coeff) && continue
        @inbounds tmp = prod( vals .^ ct[i] )
        suma += a_coeff * tmp
    end

    return suma
end

evaluate(a::HomogeneousPolynomial{T}, vals::Array{S,1} ) where
        {T<:Number,S<:NumberNotSeriesN} = _evaluate(a, (vals...,))

evaluate(a::HomogeneousPolynomial{T}, v, vals::Vararg{S,N}) where {T,S,N} = evaluate(a, (v, vals...,))

evaluate(a::HomogeneousPolynomial{T}, v) where {T<:Number} = evaluate(a, [v...])

function evaluate(a::HomogeneousPolynomial)
    a.order == 0 && return a[1]
    zero(a[1])
end

#function-like behavior for HomogeneousPolynomial
(p::HomogeneousPolynomial)(x) = evaluate(p, x)

(p::HomogeneousPolynomial)(x, v::Vararg{T, N}) where {T,N} = evaluate(p, (x, v...,))

(p::HomogeneousPolynomial)() = evaluate(p)

"""
    evaluate(a, [vals]; sorting::Bool=true)

Evaluate the `TaylorN` polynomial `a` at `vals`.
If `vals` is ommitted, it's evaluated at zero. The
keyword parameter `sorting` can be used to avoid
sorting (in increasing order by `abs2`) the
terms that are added.

Note that the syntax `a(vals)` is equivalent to
`evaluate(a, vals)`; and `a()` is equivalent to
`evaluate(a)`. No extension exists that incorporates
`sorting`.
"""
evaluate(a::TaylorN{T}, vals::NTuple; sorting::Bool=true) where {T<:Number} =
    _evaluate(a, vals, Val(sorting))

evaluate(a::TaylorN, vals; sorting::Bool=true) = _evaluate(a, (vals...,), Val(sorting))

evaluate(a::TaylorN, v, vals::Vararg{T, N}; sorting::Bool=true) where {T,N} =
    _evaluate(a, (v, vals...,), Val(sorting))

function _evaluate(a::TaylorN{T}, vals) where {T<:Number}
    @assert get_numvars() == length(vals)
    R = promote_type(T,typeof(vals[1]))
    a_length = length(a)
    suma = zeros(R, a_length)
    @inbounds for homPol in length(a):-1:1
        suma[homPol] = evaluate(a.coeffs[homPol], vals)
    end
    return suma
end

function _evaluate(a::TaylorN{T}, vals::NTuple, ::Val{true}) where {T<:Number}
    suma = _evaluate(a, vals)
    return sum( sort!(suma, by=abs2) )
end
function _evaluate(a::TaylorN{T}, vals::NTuple, ::Val{false}) where {T<:Number}
    suma = _evaluate(a, vals)
    return sum( suma )
end

function evaluate(a::TaylorN{T}, vals::NTuple{N,Taylor1{S}}) where
        {T<:Number,S<:NumberNotSeries,N}

    @assert N == get_numvars()

    R = promote_type(T,S)
    ord = maximum( get_order.(vals) )
    suma = Taylor1(zeros(R, ord))

    @inbounds for homPol in length(a):-1:1
        suma += evaluate(a.coeffs[homPol], vals)
    end

    return suma
end

evaluate(a::TaylorN{T}, vals::Array{Taylor1{S},1}) where
    {T<:Number,S<:NumberNotSeriesN} = evaluate(a, (vals...,))

function evaluate(a::TaylorN{Taylor1{T}}, vals::NTuple{N, Taylor1{T}}) where
        {T<:NumberNotSeries,N}

    @assert N == get_numvars()

    ord = maximum( get_order.(vals) )
    suma = Taylor1(zeros(T, ord))

    for homPol in length(a):-1:1
        suma += evaluate(a.coeffs[homPol], vals)
    end

    return suma
end

evaluate(a::TaylorN{Taylor1{T}}, vals::Array{Taylor1{T},1}) where
    {T<:NumberNotSeries} = evaluate(a, (vals...,))

function evaluate(a::TaylorN{T}, vals::NTuple{N, TaylorN{S}}) where
        {T<:Number,S<:NumberNotSeries,N}

    @assert length(vals) == get_numvars()

    suma = zero(constant_term(a))*vals[1]

    for homPol in length(a):-1:1
        suma += evaluate(a.coeffs[homPol], vals)
    end

    return suma
end

evaluate(a::TaylorN{T}, vals::Array{TaylorN{S},1}) where
    {T<:Number,S<:NumberNotSeries} = evaluate(a, (vals...,))

function evaluate(a::TaylorN{T}, s::Symbol, val::S) where
        {T<:Number, S<:NumberNotSeriesN}
    vars = get_variables(T)
    ind = lookupvar(s)
    vars[ind] = val
    evaluate(a, vars)
end

evaluate(a::TaylorN{T}, x::Pair{Symbol,S}) where {T<:Number,S<:NumberNotSeriesN} =
    evaluate(a, first(x), last(x))

evaluate(a::TaylorN{T}) where {T<:Number} = a[0][1]

#High-dim array evaluation
function evaluate(A::AbstractArray{TaylorN{T},N}, δx::Vector{S}) where {T<:Number,S<:Number,N}
    R = promote_type(T,S)
    return evaluate(convert(Array{TaylorN{R},N},A), convert(Vector{R},δx))
end
function evaluate(A::Array{TaylorN{T}}, δx::Vector{T}) where {T<:Number}
    Anew = Array{T}(undef, size(A)...)
    evaluate!(A, δx, Anew)
    return Anew
end
evaluate(A::AbstractArray{TaylorN{T}}) where {T<:Number} = evaluate.(A)

#function-like behavior for TaylorN
(p::TaylorN)(x) = evaluate(p, x)
(p::TaylorN)() = evaluate(p)
(p::TaylorN)(s::Symbol, x) = evaluate(p, s, x)
(p::TaylorN)(x::Pair) = evaluate(p, first(x), last(x))
(p::TaylorN)(x, v::Vararg{T, N}) where {T,N} = evaluate(p, (x, v...,))
(p::TaylorN)(b::Bool, x) = evaluate(p, x, sorting=b)
(p::TaylorN)(b::Bool, x, v::Vararg{T, N}) where {T,N} =
    evaluate(p, (x, v...,), sorting=b)

#function-like behavior for AbstractArray{TaylorN{T}}
if VERSION > v"1.1"
    (p::AbstractArray{TaylorN{T}})(x) where {T<:Number} = evaluate(p, x)
    (p::AbstractArray{TaylorN{T}})() where {T<:Number} = evaluate(p)
else
    (p::Array{TaylorN{T}})(x) where {T<:Number} = evaluate(p, x)
    (p::SubArray{TaylorN{T}})(x) where {T<:Number} = evaluate(p, x)
    (p::Array{TaylorN{T}})() where {T<:Number} = evaluate(p)
    (p::SubArray{TaylorN{T}})() where {T<:Number} = evaluate(p)
end


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

@doc raw"""ResourcePolicyRule is a predicate that matches some resource requests, testing the request's verb and the target resource. A ResourcePolicyRule matches a resource request if and only if: (a) at least one member of verbs matches the request, (b) at least one member of apiGroups matches the request, (c) at least one member of resources matches the request, and (d) least one member of namespaces matches the request.

    IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule(;
        apiGroups=nothing,
        clusterScope=nothing,
        namespaces=nothing,
        resources=nothing,
        verbs=nothing,
    )

    - apiGroups::Vector{String} : `apiGroups` is a list of matching API groups and may not be empty. \"*\" matches all API groups and, if present, must be the only entry. Required.
    - clusterScope::Bool : `clusterScope` indicates whether to match requests that do not specify a namespace (which happens either because the resource is not namespaced or the request targets all namespaces). If this field is omitted or false then the `namespaces` field must contain a non-empty list.
    - namespaces::Vector{String} : `namespaces` is a list of target namespaces that restricts matches.  A request that specifies a target namespace matches only if either (a) this list contains that target namespace or (b) this list contains \"*\".  Note that \"*\" matches any specified namespace but does not match a request that _does not specify_ a namespace (see the `clusterScope` field for that). This list may be empty, but only if `clusterScope` is true.
    - resources::Vector{String} : `resources` is a list of matching resources (i.e., lowercase and plural) with, if desired, subresource.  For example, [ \"services\", \"nodes/status\" ].  This list may not be empty. \"*\" matches all resources and, if present, must be the only entry. Required.
    - verbs::Vector{String} : `verbs` is a list of matching verbs and may not be empty. \"*\" matches all verbs and, if present, must be the only entry. Required.
"""
mutable struct IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule <: SwaggerModel
    apiGroups::Any # spec type: Union{ Nothing, Vector{String} } # spec name: apiGroups
    clusterScope::Any # spec type: Union{ Nothing, Bool } # spec name: clusterScope
    namespaces::Any # spec type: Union{ Nothing, Vector{String} } # spec name: namespaces
    resources::Any # spec type: Union{ Nothing, Vector{String} } # spec name: resources
    verbs::Any # spec type: Union{ Nothing, Vector{String} } # spec name: verbs

    function IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule(;apiGroups=nothing, clusterScope=nothing, namespaces=nothing, resources=nothing, verbs=nothing)
        o = new()
        validate_property(IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule, Symbol("apiGroups"), apiGroups)
        setfield!(o, Symbol("apiGroups"), apiGroups)
        validate_property(IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule, Symbol("clusterScope"), clusterScope)
        setfield!(o, Symbol("clusterScope"), clusterScope)
        validate_property(IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule, Symbol("namespaces"), namespaces)
        setfield!(o, Symbol("namespaces"), namespaces)
        validate_property(IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule, Symbol("resources"), resources)
        setfield!(o, Symbol("resources"), resources)
        validate_property(IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule, Symbol("verbs"), verbs)
        setfield!(o, Symbol("verbs"), verbs)
        o
    end
end # type IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule

const _property_map_IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule = Dict{Symbol,Symbol}(Symbol("apiGroups")=>Symbol("apiGroups"), Symbol("clusterScope")=>Symbol("clusterScope"), Symbol("namespaces")=>Symbol("namespaces"), Symbol("resources")=>Symbol("resources"), Symbol("verbs")=>Symbol("verbs"))
const _property_types_IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule = Dict{Symbol,String}(Symbol("apiGroups")=>"Vector{String}", Symbol("clusterScope")=>"Bool", Symbol("namespaces")=>"Vector{String}", Symbol("resources")=>"Vector{String}", Symbol("verbs")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule }) = collect(keys(_property_map_IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule))
Swagger.property_type(::Type{ IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule[name]))}
Swagger.field_name(::Type{ IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule }, property_name::Symbol) =  _property_map_IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule[property_name]

function check_required(o::IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule)
    (getproperty(o, Symbol("apiGroups")) === nothing) && (return false)
    (getproperty(o, Symbol("resources")) === nothing) && (return false)
    (getproperty(o, Symbol("verbs")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule }, name::Symbol, val)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license
# givensAlgorithm functions are derived from LAPACK, see below

abstract type AbstractRotation{T} end

transpose(R::AbstractRotation) = error("transpose not implemented for $(typeof(R)). Consider using adjoint instead of transpose.")

function *(R::AbstractRotation{T}, A::AbstractVecOrMat{S}) where {T,S}
    TS = typeof(zero(T)*zero(S) + zero(T)*zero(S))
    lmul!(convert(AbstractRotation{TS}, R), TS == S ? copy(A) : convert(AbstractArray{TS}, A))
end
*(A::AbstractVector, adjR::Adjoint{<:Any,<:AbstractRotation}) = _absvecormat_mul_adjrot(A, adjR)
*(A::AbstractMatrix, adjR::Adjoint{<:Any,<:AbstractRotation}) = _absvecormat_mul_adjrot(A, adjR)
function _absvecormat_mul_adjrot(A::AbstractVecOrMat{T}, adjR::Adjoint{<:Any,<:AbstractRotation{S}}) where {T,S}
    R = adjR.parent
    TS = typeof(zero(T)*zero(S) + zero(T)*zero(S))
    rmul!(TS == T ? copy(A) : convert(AbstractArray{TS}, A), adjoint(convert(AbstractRotation{TS}, R)))
end
"""
    LinearAlgebra.Givens(i1,i2,c,s) -> G

A Givens rotation linear operator. The fields `c` and `s` represent the cosine and sine of
the rotation angle, respectively. The `Givens` type supports left multiplication `G*A` and
conjugated transpose right multiplication `A*G'`. The type doesn't have a `size` and can
therefore be multiplied with matrices of arbitrary size as long as `i2<=size(A,2)` for
`G*A` or `i2<=size(A,1)` for `A*G'`.

See also: [`givens`](@ref)
"""
struct Givens{T} <: AbstractRotation{T}
    i1::Int
    i2::Int
    c::T
    s::T
end
mutable struct Rotation{T} <: AbstractRotation{T}
    rotations::Vector{Givens{T}}
end

convert(::Type{T}, r::T) where {T<:AbstractRotation} = r
convert(::Type{T}, r::AbstractRotation) where {T<:AbstractRotation} = T(r)

Givens{T}(G::Givens{T}) where {T} = G
Givens{T}(G::Givens) where {T} = Givens(G.i1, G.i2, convert(T, G.c), convert(T, G.s))
Rotation{T}(R::Rotation{T}) where {T} = R
Rotation{T}(R::Rotation) where {T} = Rotation{T}([Givens{T}(g) for g in R.rotations])
AbstractRotation{T}(G::Givens) where {T} = Givens{T}(G)
AbstractRotation{T}(R::Rotation) where {T} = Rotation{T}(R)

adjoint(G::Givens) = Adjoint(G)
adjoint(R::Rotation) = Adjoint(R)
Base.copy(aG::Adjoint{<:Any,<:Givens}) = (G = aG.parent; Givens(G.i1, G.i2, conj(G.c), -G.s))
Base.copy(aR::Adjoint{<:Any,Rotation{T}}) where {T} = Rotation{T}(reverse!([copy(r') for r in aR.parent.rotations]))

realmin2(::Type{Float32}) = reinterpret(Float32, 0x26000000)
realmin2(::Type{Float64}) = reinterpret(Float64, 0x21a0000000000000)
realmin2(::Type{T}) where {T} = (twopar = 2one(T); twopar^trunc(Integer,log(realmin(T)/eps(T))/log(twopar)/twopar))

# derived from LAPACK's dlartg
# Copyright:
# Univ. of Tennessee
# Univ. of California Berkeley
# Univ. of Colorado Denver
# NAG Ltd.
function givensAlgorithm(f::T, g::T) where T<:AbstractFloat
    onepar = one(T)
    twopar = 2one(T)
    T0 = typeof(onepar) # dimensionless
    zeropar = T0(zero(T)) # must be dimensionless

    # need both dimensionful and dimensionless versions of these:
    safmn2 = realmin2(T0)
    safmn2u = realmin2(T)
    safmx2 = one(T)/safmn2
    safmx2u = oneunit(T)/safmn2

    if g == 0
        cs = onepar
        sn = zeropar
        r = f
    elseif f == 0
        cs = zeropar
        sn = onepar
        r = g
    else
        f1 = f
        g1 = g
        scalepar = max(abs(f1), abs(g1))
        if scalepar >= safmx2u
            count = 0
            while true
                count += 1
                f1 *= safmn2
                g1 *= safmn2
                scalepar = max(abs(f1), abs(g1))
                if scalepar < safmx2u break end
            end
            r = sqrt(f1*f1 + g1*g1)
            cs = f1/r
            sn = g1/r
            for i = 1:count
                r *= safmx2
            end
        elseif scalepar <= safmn2u
            count = 0
            while true
                count += 1
                f1 *= safmx2
                g1 *= safmx2
                scalepar = max(abs(f1), abs(g1))
                if scalepar > safmn2u break end
            end
            r = sqrt(f1*f1 + g1*g1)
            cs = f1/r
            sn = g1/r
            for i = 1:count
                r *= safmn2
            end
        else
            r = sqrt(f1*f1 + g1*g1)
            cs = f1/r
            sn = g1/r
        end
        if abs(f) > abs(g) && cs < 0
            cs = -cs
            sn = -sn
            r = -r
        end
    end
    return cs, sn, r
end

# derived from LAPACK's zlartg
# Copyright:
# Univ. of Tennessee
# Univ. of California Berkeley
# Univ. of Colorado Denver
# NAG Ltd.
function givensAlgorithm(f::Complex{T}, g::Complex{T}) where T<:AbstractFloat
    twopar, onepar = 2one(T), one(T)
    T0 = typeof(onepar) # dimensionless
    zeropar = T0(zero(T)) # must be dimensionless
    czero = complex(zeropar)

    abs1(ff) = max(abs(real(ff)), abs(imag(ff)))
    safmin = realmin(T0)
    safmn2 = realmin2(T0)
    safmn2u = realmin2(T)
    safmx2 = one(T)/safmn2
    safmx2u = oneunit(T)/safmn2
    scalepar = max(abs1(f), abs1(g))
    fs = f
    gs = g
    count = 0
    if scalepar >= safmx2u
        while true
            count += 1
            fs *= safmn2
            gs *= safmn2
            scalepar *= safmn2
            if scalepar < safmx2u break end
        end
    elseif scalepar <= safmn2u
        if g == 0
            cs = onepar
            sn = czero
            r = f
            return cs, sn, r
        end
        while true
            count -= 1
            fs *= safmx2
            gs *= safmx2
            scalepar *= safmx2
            if scalepar > safmn2u break end
        end
    end
    f2 = abs2(fs)
    g2 = abs2(gs)
    if f2 <= max(g2, oneunit(T))*safmin
        # This is a rare case: F is very small.
        if f == 0
            cs = zero(T)
            r = complex(hypot(real(g), imag(g)))
            # do complex/real division explicitly with two real divisions
            d = hypot(real(gs), imag(gs))
            sn = complex(real(gs)/d, -imag(gs)/d)
            return cs, sn, r
        end
        f2s = hypot(real(fs), imag(fs))
        # g2 and g2s are accurate
        # g2 is at least safmin, and g2s is at least safmn2
        g2s = sqrt(g2)
        # error in cs from underflow in f2s is at most
        # unfl / safmn2 .lt. sqrt(unfl*eps) .lt. eps
        # if max(g2,one)=g2, then f2 .lt. g2*safmin,
        # and so cs .lt. sqrt(safmin)
        # if max(g2,one)=one, then f2 .lt. safmin
        # and so cs .lt. sqrt(safmin)/safmn2 = sqrt(eps)
        # therefore, cs = f2s/g2s / sqrt( 1 + (f2s/g2s)**2 ) = f2s/g2s
        cs = f2s/g2s
        # make sure abs(ff) = 1
        # do complex/real division explicitly with 2 real divisions
        if abs1(f) > 1
            d = hypot(real(f), imag(f))
            ff = complex(real(f)/d, imag(f)/d)
        else
            dr = safmx2*real(f)
            di = safmx2*imag(f)
            d = hypot(dr, di)
            ff = complex(dr/d, di/d)
        end
        sn = ff*complex(real(gs)/g2s, -imag(gs)/g2s)
        r = cs*f + sn*g
    else
        # This is the most common case.
        # Neither F2 nor F2/G2 are less than SAFMIN
        # F2S cannot overflow, and it is accurate
        f2s = sqrt(onepar + g2/f2)
        # do the f2s(real)*fs(complex) multiply with two real multiplies
        r = complex(f2s*real(fs), f2s*imag(fs))
        cs = onepar/f2s
        d = f2 + g2
        # do complex/real division explicitly with two real divisions
        sn = complex(real(r)/d, imag(r)/d)
        sn *= conj(gs)
        if count != 0
            if count > 0
                for i = 1:count
                    r *= safmx2
                end
            else
                for i = 1:-count
                    r *= safmn2
                end
            end
        end
    end
    return cs, sn, r
end

"""

    givens(f::T, g::T, i1::Integer, i2::Integer) where {T} -> (G::Givens, r::T)

Computes the Givens rotation `G` and scalar `r` such that for any vector `x` where
```
x[i1] = f
x[i2] = g
```
the result of the multiplication
```
y = G*x
```
has the property that
```
y[i1] = r
y[i2] = 0
```

See also: [`LinearAlgebra.Givens`](@ref)
"""
function givens(f::T, g::T, i1::Integer, i2::Integer) where T
    if i1 == i2
        throw(ArgumentError("Indices must be distinct."))
    end
    c, s, r = givensAlgorithm(f, g)
    if i1 > i2
        s = -conj(s)
        i1,i2 = i2,i1
    end
    Givens(i1, i2, convert(T, c), convert(T, s)), r
end
"""
    givens(A::AbstractArray, i1::Integer, i2::Integer, j::Integer) -> (G::Givens, r)

Computes the Givens rotation `G` and scalar `r` such that the result of the multiplication
```
B = G*A
```
has the property that
```
B[i1,j] = r
B[i2,j] = 0
```

See also: [`LinearAlgebra.Givens`](@ref)
"""
givens(A::AbstractMatrix, i1::Integer, i2::Integer, j::Integer) =
    givens(A[i1,j], A[i2,j],i1,i2)

"""
    givens(x::AbstractVector, i1::Integer, i2::Integer) -> (G::Givens, r)

Computes the Givens rotation `G` and scalar `r` such that the result of the multiplication
```
B = G*x
```
has the property that
```
B[i1] = r
B[i2] = 0
```

See also: [`LinearAlgebra.Givens`](@ref)
"""
givens(x::AbstractVector, i1::Integer, i2::Integer) =
    givens(x[i1], x[i2], i1, i2)

function getindex(G::Givens, i::Integer, j::Integer)
    if i == j
        if i == G.i1 || i == G.i2
            G.c
        else
            oneunit(G.c)
        end
    elseif i == G.i1 && j == G.i2
        G.s
    elseif i == G.i2 && j == G.i1
        -conj(G.s)
    else
        zero(G.s)
    end
end

function lmul!(G::Givens, A::AbstractVecOrMat)
    m, n = size(A, 1), size(A, 2)
    if G.i2 > m
        throw(DimensionMismatch("column indices for rotation are outside the matrix"))
    end
    @inbounds @simd for i = 1:n
        a1, a2 = A[G.i1,i], A[G.i2,i]
        A[G.i1,i] =       G.c *a1 + G.s*a2
        A[G.i2,i] = -conj(G.s)*a1 + G.c*a2
    end
    return A
end
function rmul!(A::AbstractMatrix, adjG::Adjoint{<:Any,<:Givens})
    G = adjG.parent
    m, n = size(A, 1), size(A, 2)
    if G.i2 > n
        throw(DimensionMismatch("column indices for rotation are outside the matrix"))
    end
    @inbounds @simd for i = 1:m
        a1, a2 = A[i,G.i1], A[i,G.i2]
        A[i,G.i1] =  a1*G.c + a2*conj(G.s)
        A[i,G.i2] = -a1*G.s + a2*G.c
    end
    return A
end

function lmul!(G::Givens, R::Rotation)
    push!(R.rotations, G)
    return R
end
function lmul!(R::Rotation, A::AbstractMatrix)
    @inbounds for i = 1:length(R.rotations)
        lmul!(R.rotations[i], A)
    end
    return A
end
function rmul!(A::AbstractMatrix, adjR::Adjoint{<:Any,<:Rotation})
    R = adjR.parent
    @inbounds for i = 1:length(R.rotations)
        rmul!(A, adjoint(R.rotations[i]))
    end
    return A
end
*(G1::Givens{T}, G2::Givens{T}) where {T} = Rotation(push!(push!(Givens{T}[], G2), G1))

# TODO: None of the following disambiguation methods are great. They should perhaps
# instead be MethodErrors, or revised.
#
# dismabiguation methods: *(Adj/Trans of AbsVec or AbsMat, Adj of AbstractRotation)
*(A::Adjoint{<:Any,<:AbstractVector}, B::Adjoint{<:Any,<:AbstractRotation}) = copy(A) * B
*(A::Adjoint{<:Any,<:AbstractMatrix}, B::Adjoint{<:Any,<:AbstractRotation}) = copy(A) * B
*(A::Transpose{<:Any,<:AbstractVector}, B::Adjoint{<:Any,<:AbstractRotation}) = copy(A) * B
*(A::Transpose{<:Any,<:AbstractMatrix}, B::Adjoint{<:Any,<:AbstractRotation}) = copy(A) * B
# dismabiguation methods: *(Adj/Trans of AbsTri or RealHermSymComplex{Herm|Sym}, Adj of AbstractRotation)
*(A::Adjoint{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:AbstractRotation}) = copy(A) * B
*(A::Transpose{<:Any,<:AbstractTriangular}, B::Adjoint{<:Any,<:AbstractRotation}) = copy(A) * B
*(A::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::Adjoint{<:Any,<:AbstractRotation}) = copy(A) * B
*(A::Transpose{<:Any,<:RealHermSymComplexSym}, B::Adjoint{<:Any,<:AbstractRotation}) = copy(A) * B
# dismabiguation methods: *(Diag/AbsTri, Adj of AbstractRotation)
*(A::Diagonal, B::Adjoint{<:Any,<:AbstractRotation}) = A * copy(B)
*(A::AbstractTriangular, B::Adjoint{<:Any,<:AbstractRotation}) = A * copy(B)

max_stars_repo_name: nreigl/julia

max_issues_repo_name: BenChung/Julia-BTI

max_forks_repo_name: BenChung/Julia-BTI


text:
using SWIM
using Test

@testset "SWIM.jl" begin
    c = Canvas(10, 10, 2.83)
    initialize(c)

    nodes = generate_nodes(c, 10, 100)
    simulate(c,nodes,100,2.0)

end


text:
#!/usr/bin/env julia
#
# This is an adapted version of format.jl from JuliaFormatter with the
# following license:
#
#    MIT License
#    Copyright (c) 2019 Dominique Luna
#
#    Permission is hereby granted, free of charge, to any person obtaining a
#    copy of this software and associated documentation files (the
#    "Software"), to deal in the Software without restriction, including
#    without limitation the rights to use, copy, modify, merge, publish,
#    distribute, sublicense, and/or sell copies of the Software, and to permit
#    persons to whom the Software is furnished to do so, subject to the
#    following conditions:
#
#    The above copyright notice and this permission notice shall be included
#    in all copies or substantial portions of the Software.
#
#    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
#    OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
#    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
#    NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
#    DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
#    OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
#    USE OR OTHER DEALINGS IN THE SOFTWARE.
#
# and from climaformat.jl from ClimateMachine.jl with the following license:
#
#    Copyright 2019 Climate Modeling Alliance
#
#    Licensed under the Apache License, Version 2.0 (the "License");
#    you may not use this file except in compliance with the License.
#    You may obtain a copy of the License at
#
#        http://www.apache.org/licenses/LICENSE-2.0
#
#    Unless required by applicable law or agreed to in writing, software
#    distributed under the License is distributed on an "AS IS" BASIS,
#    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#    See the License for the specific language governing permissions and
#    limitations under the License.

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using JuliaFormatter

include("formatter_options.jl")

help = """
Usage: format.jl [flags] [FILE/PATH]...

Formats the given julia files using the CLIMA formatting options.  If paths
are given it will format the julia files in the paths. Otherwise, it will
format all changed julia files.

    -v, --verbose
        Print the name of the files being formatted with relevant details.

    -h, --help
        Print this message.
"""

function parse_opts!(args::Vector{String})
  i = 1
  opts = Dict{Symbol, Union{Int, Bool}}()
  while i ≤ length(args)
    arg = args[i]
    if arg[1] != '-'
      i += 1
      continue
    end
    if arg == "-v" || arg == "--verbose"
      opt = :verbose
    elseif arg == "-h" || arg == "--help"
      opt = :help
    else
      error("invalid option $arg")
    end
    if opt in (:verbose, :help)
      opts[opt] = true
      deleteat!(args, i)
    end
  end
  return opts
end

opts = parse_opts!(ARGS)
if haskey(opts, :help)
  write(stdout, help)
  exit(0)
end
filenames = ARGS

format(filenames; formatter_options..., opts...)


text:
export planeconnect

"""
    planeconnect(nodes::Array{Node,1})
    planeconnect(node::Node)

Provides objects needed to connect to a plane.

Returns a tuple of a `deepcopy` of the `nodes` passed and a group of `Equiv` objects
 connecting the nodes passed to the nodes returned.

`(planenodearray, equivgroup) = planeconnect(nodearray)`
"""
function planeconnect(nodes::Array{Node,1})
  planenodearray = deepcopy(nodes)
  equivarray = Array{Equiv}(undef,length(nodes))
  for i in eachindex(nodes)
    equivarray[i] = Equiv([nodes[i],planenodearray[i]])
  end
  equivgroup = Group(equivarray)
  return (planenodearray, equivgroup)
end

planeconnect(node::Node) = planeconnect([node])


text:
module CodecZstd

export
    ZstdCompressor,
    ZstdCompressorStream,
    ZstdDecompressor,
    ZstdDecompressorStream

import TranscodingStreams:
    TranscodingStreams,
    TranscodingStream,
    Memory,
    Error,
    initialize,
    finalize,
    splitkwargs

using Zstd_jll

# Module containing directly wrapped ccalls generated by CLang.jl
include("libzstd/LibZstd.jl")
# Library functions that have a Julian interface. This file originally preceded the above module
include("libzstd.jl")
include("compression.jl")
include("decompression.jl")

end # module


text:
__precompile__(true)
module MacroTools

using Compat
using Compat.Markdown
export @match, @capture

include("match.jl")
include("types.jl")
include("union.jl")
include("macro.jl")
include("utils.jl")

include("examples/destruct.jl")
include("examples/threading.jl")
include("examples/forward.jl")

const animals = Symbol[]

function __init__()
  animals_file = joinpath(dirname(@__FILE__), "..", "animals.txt")
  _animals = split(read(animals_file, String))
  resize!(animals, length(_animals))
  animals .= Symbol.(lowercase.(_animals))
  Compat.Random.shuffle!(animals)
end

end # module


text:
module Zm_par

using Printf, JuLIP, Distributions, LinearAlgebra, ASE

export Stationary, VelocityVerlet, Zmethod, MaxwellBoltzmann_scale, MaxwellBoltzmann, kB, fs

kB = 8.617330337217213e-05 #units taken from Python ASE
fs = 0.09822694788464063

function MaxwellBoltzmann(at, temp)
    d = Normal()
    M = reshape(collect(rand(d, 3*length(at.M))), (length(at.M),3))
    M2 = M .* sqrt.(at.M .* temp)
    set_momenta!(at, collect(transpose(M2)))

    return at
end

function Stationary(at)
    p0 = sum(at.P)
    mtot = sum(at.M)
    v0 = p0 / mtot

    momenta = [at.M[i] * v0 for i in 1:length(at.M)]
    set_momenta!(at, collect(at.P - momenta))

    return at
end

function VelocityVerlet(IP, at, dt)
    V = at.P ./ at.M
    A = forces(IP, at) ./ at.M

    set_positions!(at, at.X + (V .* dt) + (.5 * A * dt^2))

    nA = forces(IP, at) ./ at.M
    nV = V + (.5 * (A + nA) * dt)
    set_momenta!(at, collect(nV .* at.M))

    return at
end

function MaxwellBoltzmann_scale(at, temp)
    d = Normal()
    M = reshape(collect(rand(d, 3*length(at.M))), (length(at.M),3))

    N = []

    for j in 1:1000
        d = Normal()
        M = reshape(collect(rand(d, 3*length(at.M))), (length(at.M),3))
        n = 0
        for i in 1:1:length(M[:,1])
            n += norm(M[i,:])
        end
        push!(N,n)
    end

    d = Normal()
    M = reshape(collect(rand(d, 3*length(at.M))), (length(at.M),3))
    Mnm = mean(N)

    for i in 1:10000
        n = 0
        for i in 1:length(M[:,1])
            n += norm(M[i,:])
        end
        if n < Mnm
            M = M*1.0001
        else
            M = M*0.9999
        end
    end

    M2 = M .* sqrt.(at.M .* temp)
    set_momenta!(at, transpose(M2))

    return at #, Ml, Mnm, Sl, M
end

function Zmethod(IP, at, nsteps, R, dt, A, N, save_config)
    E0 = energy(IP, at)

    m = at.M

    E_tot = zeros(nsteps*R)
    E_pot = zeros(nsteps*R)
    E_kin = zeros(nsteps*R)
    P = zeros(nsteps*R)
    T = zeros(nsteps*R)

    al = []

    for j in 0:R-1
        for i in 1:nsteps
            k = (j*nsteps)+i
            at = VelocityVerlet(IP, at, dt * fs)
            Ek = ((0.5 * sum(at.M) * norm(at.P ./ at.M)^2)/length(at.M)) / length(at.M)
            Ep = (energy(IP, at) - E0) / length(at.M)
            E_tot[k] = Ek + Ep
            E_pot[k] = Ep
            E_kin[k] = Ek
            T[k] = Ek / (1.5 * kB)
            P[k] = -tr(stress(IP, at))/3.0

            if k % 100 == 0
                Temp =  (Ek / (1.5 * kB))
                println("iteration: ($k), temp: ($Temp)")
            end

            if i % save_config == 0
                push!(al, deepcopy(at))
            end
        end

        v = at.P ./ m
        C = A/norm(v)

        set_momenta!(at, collect((v + C*v) .* m))
    end

    return E_tot, E_pot, E_kin, P, T, al
end

end


text:
if !@isdefined(TabLevel)
    TabLevel = ""
end
println(TabLevel*"Open OpIsLess.jl")
TabLevel=TabLevel*"    "

import Base.isless

###Inequality used for sorting
isless(x::MajElem, y::MajElem) = OpIsLess(x,y)
isless(x::MajoranaOp,y ::MajoranaOp) = OpIsLess(x,y)

function OpIsLess(a::MajoranaOp,b::MajoranaOp) ###Check if operators are ordered
    if a.OpList < b.OpList ##Smaller
        return true
    else ###Check for the reverse
        if a.OpList > b.OpList ##larger
            return false
        else ###Check if they are the same
            if b.OpList == a.OpList ##larger
                return true
            else
                error("Elements $a and $b cannot be ordered!")
            end
        end
    end
end

function OpIsLess(ElemA::MajElem,ElemB::MajElem)
    ###Check if majorana elements are ordered
    if ElemA.Type == ElemB.Type ##Same Operator
        return ElemA.Index < ElemB.Index ###Order by the second index
    else  ###Order by the first index
        return do_OpOrdered(ElemA.Type,ElemB.Type) ##If OpA < OpB?
    end
end

function OpOrdered(a::MajoranaOp,b::MajoranaOp) ###Check if operators are ordered
    return a.OpList <= b.OpList
end

function do_OpOrdered(OpA::String,OpB::String)
    ####Check if OpA<=OpB? (If they are the same we assume they are ordered)
    ###The order is bx < c < by < bz such that by and bz get pushed to the right and can be converted if necessary...
    ###We assume that OpA and OpB is in the desired range
    if OpA == "bx"
        return true ###bx is always first
    elseif OpA == "c"
        if OpB == "bx"
            return false
        else
            return true
        end
    elseif OpA == "by"
        if OpB == "bx" || OpB == "c"
            return false
        else
            return true
        end
    elseif OpA == "bz"
        if OpB == "bz"
            return true
        else
            return false
        end
    else
        error("Elements $OpA and $OpB cannot be ordered!")
    end
end

TabLevel=TabLevel[1:end-4]
println(TabLevel*"Close OpIsLess.jl")


text:
import Distances: pairwise
import Base: eachindex

"""
    CondCache

Cache for use by `ConditionalMean`s, `ConditionalKernel`s and `ConditionalCrossKernel`s.
Avoids recomputing the covariance `Σff` and the Kriging vector `α`.
"""
struct CondCache
    Σff::LazyPDMat
    α::AV{<:Real}
    function CondCache(mf::AV{<:Real}, Σff::AM, f::AV{<:Real})
        return new(Σff, Σff \ (f - mf))
    end
end
function CondCache(kff::Kernel, μf::MeanFunction, X::AV, f::AV{<:Real})
    return CondCache(map(μf, X), pairwise(kff, X), f)
end
show(io::IO, ::CondCache) = show(io, "CondCache")

"""
    ConditionalMean <: MeanFunction

Computes the mean of a GP `g` conditioned on observations of another GP `f`.
"""
struct ConditionalMean <: MeanFunction
    c::CondCache
    μg::MeanFunction
    kfg::CrossKernel
end
length(μ::ConditionalMean) = length(μ.μg)
eachindex(μ::ConditionalMean) = eachindex(μ.μg)
(μ::ConditionalMean)(x::Number) = map(μ, [x])[1]
(μ::ConditionalMean)(x::AbstractVector) = map(μ, ColsAreObs(reshape(x, length(x), 1)))[1]
function _map(μ::ConditionalMean, Xg::AV)
    return map(μ.μg, Xg) + pairwise(μ.kfg, :, Xg)' * μ.c.α
end

"""
    ConditionalKernel <: Kernel

Computes the (co)variance of a GP `g` conditioned on observations of another GP `f`.
"""
struct ConditionalKernel <: Kernel
    c::CondCache
    kfg::CrossKernel
    kgg::Kernel
end
(k::ConditionalKernel)(x::Number, x′::Number) = map(k, [x], [x′])[1]
(k::ConditionalKernel)(x::AV, x′::AV) =
    map(k, ColsAreObs(reshape(x, length(x), 1)), ColsAreObs(reshape(x′, length(x′), 1)))[1]
length(k::ConditionalKernel) = length(k.kgg)
eachindex(k::ConditionalKernel) = eachindex(k.kgg)

function _map(k::ConditionalKernel, X::AV)
    σgg = map(k.kgg, X)
    Σfg_X = pairwise(k.kfg, :, X)
    σ′gg = diag_Xᵀ_invA_X(k.c.Σff, Σfg_X)
    return (σgg .- σ′gg) .* .!(σgg .≈ σ′gg)
end
function _map(k::ConditionalKernel, X::AV, X′::AV)
    σgg = map(k.kgg, X, X′)
    Σfg_X = pairwise(k.kfg, :, X)
    Σfg_X′ = pairwise(k.kfg, :, X′)
    σ′gg = diag_Xᵀ_invA_Y(Σfg_X, k.c.Σff, Σfg_X′)
    return (σgg .- σ′gg) .* .!(σgg .≈ σ′gg)
end

function _pairwise(k::ConditionalKernel, X::AV)
    Σgg = AbstractMatrix(pairwise(k.kgg, X))
    Σfg_X = pairwise(k.kfg, :, X)
    Σ′gg = AbstractMatrix(Xt_invA_X(k.c.Σff, Σfg_X))
    return (Σgg .- Σ′gg) .* .!(Σgg .≈ Σ′gg)
end
function _pairwise(k::ConditionalKernel, X::AV, X′::AV)
    Σgg = pairwise(k.kgg, X, X′)
    Σfg_X = pairwise(k.kfg, :, X)
    Σfg_X′ = pairwise(k.kfg, :, X′)
    Σ′gg = Xt_invA_Y(Σfg_X, k.c.Σff, Σfg_X′)
    return (Σgg .- Σ′gg) .* .!(Σgg .≈ Σ′gg)
end

"""
    ConditionalCrossKernel <: CrossKernel

Computes the covariance between `g` and `h` conditioned on observations of a process `f`.
"""
struct ConditionalCrossKernel <: CrossKernel
    c::CondCache
    kfg::CrossKernel
    kfh::CrossKernel
    kgh::CrossKernel
end
(k::ConditionalCrossKernel)(x::Number, x′::Number) = map(k, [x], [x′])[1]
(k::ConditionalCrossKernel)(x::AV, x′::AV) =
    map(k, ColsAreObs(reshape(x, length(x), 1)), ColsAreObs(reshape(x′, length(x′), 1)))[1]
size(k::ConditionalCrossKernel, N::Int) = size(k.kgh, N)
eachindex(k::ConditionalCrossKernel, N::Int) = eachindex(k.kgh, N)

function _map(k::ConditionalCrossKernel, X::AV, X′::AV)
    σgh = map(k.kgh, X, X′)
    Σfg_X = pairwise(k.kfg, :, X)
    Σfh_X′ = pairwise(k.kfh, :, X′)
    σ′gh = diag_Xᵀ_invA_Y(Σfg_X, k.c.Σff, Σfh_X′)
    return (σgh .- σ′gh) .* .!(σgh .≈ σ′gh)
end

function _pairwise(k::ConditionalCrossKernel, Xg::AV, Xh::AV)
    Σgh = pairwise(k.kgh, Xg, Xh)
    Σfg_Xg = pairwise(k.kfg, :, Xg)
    Σfh_Xh = pairwise(k.kfh, :, Xh)
    Σ′gh = Xt_invA_Y(Σfg_Xg, k.c.Σff, Σfh_Xh)
    return (Σgh .- Σ′gh) .* .!(Σgh .≈ Σ′gh)
end