text:
"""
    OneHotArray(index => value, shape) -> array::AbstractArray

Create an `array` such that `isequal(array[index], value)`.  The values at
indices other than `index` are undefined.

Currently, this array type is likely only useful as the second argument to
`BangBang.Extras.broadcast_inplace!!` and as the input argument to FLoops.jl's
`@reduce`

# Examples
```jldoctest
julia> using MicroCollections, BangBang.Extras

julia> broadcast_inplace!!(+, ones(Int64, 4), OneHotVector(2 => 3, 4))
4-element Vector{Int64}:
 1
 4
 1
 1

julia> using InitialValues

julia> broadcast_inplace!!(+, InitialValue(+), OneHotVector(2 => 3, 4))
4-element Vector{Int64}:
 0
 3
 0
 0
```
"""
OneHotArray

using BangBang: BangBang, setindex!!
using InitialValues: GenericInitialValue

struct OneHotArray{T,N,Index<:Tuple,Value} <: AbstractArray{T,N}
    index::Index
    value::Value
    shape::NTuple{N,Int}

    global _OneHotArray
    @inline _OneHotArray(
        ::Type{T},
        index::Index,
        value::Value,
        shape::NTuple{N,Int},
    ) where {T,N,Index,Value} = new{T,N,Index,Value}(index, value, shape)
end

const OneHotVector{T,Index<:Tuple,Value} = OneHotArray{T,1,Index,Value}
const OneHotMatrix{T,Index<:Tuple,Value} = OneHotArray{T,2,Index,Value}

indexof(oh::OneHotArray) = oh.index
valueof(oh::OneHotArray) = oh.value

const DimsLike = Union{Tuple{Vararg{Integer}},Integer}

astuple(index) = (index,)
astuple(index::Tuple) = index

@inline function OneHotArray{T}(index_value::Pair, shape::DimsLike) where {T}
    index = astuple(first(index_value))
    value = last(index_value)
    array = _OneHotArray(T, index, value, astuple(shape))
    @boundscheck checkbounds(array, index...)
    return array
end

@inline OneHotArray(index_value::Pair, shape::DimsLike) =
    OneHotArray{typeof(last(index_value))}(index_value, shape)

@inline OneHotArray{T,N}(index_value::Pair, shape::NTuple{N,<:Integer}) where {T,N} =
    OneHotArray{T}(index_value, shape)::OneHotArray{T,N}

@inline OneHotArray{<:Any,N}(index_value::Pair{<:Any,T}, shape::DimsLike) where {N,T} =
    OneHotArray(index_value, astuple(shape))::OneHotArray{<:Any,N}

@inline OneHotVector{T}(index_value::Pair, length::Integer) where {T,N} =
    OneHotArray{T}(index_value, (length,))::OneHotArray{T,N}

Base.size(A::OneHotArray) = A.shape

@inline function Base.getindex(A::OneHotArray{<:Any,N}, I::Vararg{Int,N}) where {N}
    index::NTuple{N,Int} = to_indices(A, A.index)
    if I == index
        return A.value
    else
        throw(UndefRefError())
    end
end

commonsize(a::OneHotArray) = size(a)
Base.@propagate_inbounds function commonsize(a::OneHotArray, b::OneHotArray, rest...)
    @boundscheck if size(a) != size(b)
        error("incompatible size: size(a) = ", size(a), " size(b) = ", size(b))
    end
    commonsize(a, rest...)
    return size(a)
end

Base.@propagate_inbounds function BangBang.Extras.broadcast_inplace!!(
    f,
    inputoutput,
    oh::OneHotArray,
)
    i = indexof(oh)
    v = valueof(oh)
    return setindex!!(inputoutput, f(inputoutput[i...], v), i...)
end

Base.@propagate_inbounds function BangBang.Extras.broadcast_inplace!!(
    op::OP,
    ::GenericInitialValue{OP},
    oh::OneHotArray,
    rest::OneHotArray...,
) where {OP}
    T = mapfoldl(eltype, promote_type, rest; init = eltype(oh))
    output = fill(Base.reduce_empty(op, T), commonsize(oh, rest...))
    output[indexof(oh)...] = valueof(oh)
    for oh in rest  # TODO: unroll?
        output[indexof(oh)...] = valueof(oh)
    end
    return output
end


text:
using Test
import ClimaCore.DataLayouts: VF, IJFH, VIJFH, slab, column

@testset "VIJFH" begin
    Nv = 10 # number of vertical levels
    Nij = 4 # Nij × Nij nodal points per element
    Nh = 10 # number of elements

    TestFloatTypes = (Float32, Float64)
    for FT in TestFloatTypes
        S = Tuple{Complex{FT}, FT}

        # construct a data object with 10 cells in vertical and
        # 10 elements in horizontal with 4 × 4 nodal points per element in horizontal
        array = rand(FT, Nv, Nij, Nij, 3, Nh)

        data = VIJFH{S, Nij}(array)

        @test getfield(data.:1, :array) == @view(array[:, :, :, 1:2, :])
        @test getfield(data.:2, :array) == @view(array[:, :, :, 3:3, :])

        @test size(data) == (Nij, Nij, 1, Nv, Nh)

        # test tuple assignment on columns
        val = (Complex{FT}(-1.0, -2.0), FT(-3.0))

        column(data, 1, 2, 1)[1] = val
        @test array[1, 1, 2, 1, 1] == -1.0
        @test array[1, 1, 2, 2, 1] == -2.0
        @test array[1, 1, 2, 3, 1] == -3.0

        # test value of assing tuple on slab
        sdata = slab(data, 1, 1)
        @test sdata[1, 2] == val

        # sum of all the first field elements
        @test sum(data.:1) ≈
              Complex{FT}(sum(array[:, :, :, 1, :]), sum(array[:, :, :, 2, :]))
        @test sum(x -> x[2], data) ≈ sum(array[:, :, :, 3, :])
    end
end

@testset "VIJFH type safety" begin
    Nv = 1 # number of vertical levels
    Nij = 2 # Nij × Nij nodal points per element
    Nh = 1 # number of elements

    # check that types of the same bitstype throw a conversion error
    SA = (a = 1.0, b = 2.0)
    SB = (c = 1.0, d = 2.0)

    array = zeros(Float64, Nv, Nij, Nij, 2, Nh)
    data = VIJFH{typeof(SA), Nij}(array)

    cdata = column(data, 1, 2, 1)
    cdata[1] = SA
    @test cdata[1] isa typeof(SA)
    @test_throws MethodError cdata[1] = SB

    sdata = slab(data, 1, 1)
    @test sdata[1, 2] isa typeof(SA)
    @test_throws MethodError sdata[1] = SB
end

@testset "broadcasting between VIJFH data object + scalars" begin
    FT = Float64
    array = ones(FT, 2, 2, 2, 2, 2)
    S = Complex{Float64}
    data1 = VIJFH{S, 2}(array)
    res = data1 .+ 1
    @test res isa VIJFH{S}
    @test parent(res) == FT[
        f == 1 ? 2 : 1 for v in 1:2, i in 1:2, j in 1:2, f in 1:2, h in 1:2
    ]
    @test sum(res) == Complex(FT(32.0), FT(16.0))
    @test sum(Base.Broadcast.broadcasted(+, data1, 1)) ==
          Complex(FT(32.0), FT(16.0))
end

@testset "broadcasting between VF + IJFH data object => VIJFH" begin
    FT = Float64
    S = Complex{FT}
    data_vf = VF{S}(ones(FT, 3, 2))
    data_ijfh = IJFH{FT, 2}(ones(FT, 2, 2, 1, 2))
    data_vijfh = data_vf .+ data_ijfh
    @test data_vijfh isa VIJFH{S}
    @test size(data_vijfh) == (2, 2, 1, 3, 2)

    @test parent(data_vijfh) == FT[
        f == 1 ? 2 : 1 for v in 1:3, i in 1:2, j in 1:2, f in 1:2, h in 1:2
    ]
    @test parent(data_vijfh .+ data_vf) == FT[
        f == 1 ? 3 : 2 for v in 1:3, i in 1:2, j in 1:2, f in 1:2, h in 1:2
    ]
    @test parent(data_vijfh .+ data_ijfh) == FT[
        f == 1 ? 3 : 1 for v in 1:3, i in 1:2, j in 1:2, f in 1:2, h in 1:2
    ]
end

ext: jl


text:
using Zygote, StatsFuns
using VIEstimators: chol

@testset "util" begin
    to_test = [
        to_positive_definite,
        to_positive_definite_softplus,
        to_positive_definite_square,
    ]
    @testset "$to_psd" for to_psd in to_test
        rng = MersenneTwister(123456)
        D = 11
        A = randn(rng, D, D)
        U = UpperTriangular(A)
        V = to_psd(U)

        @test all(diag(V) .> 0)

        adjoint_test(
            A->to_psd(UpperTriangular(A)),
            UpperTriangular(randn(rng, D, D)),
            A,
        )

        if to_psd === to_positive_definite_softplus
            @test to_positive_definite_softplus_inv(V) ≈ U

            ΔV = randn(rng, D, D)
            ΔV[diagind(ΔV)] .= softplus.(ΔV[diagind(ΔV)])
            frule_test(to_positive_definite_softplus_inv, (V, ΔV); fdm=forward_fdm(5, 1))
        end
    end
    @testset "upper_triangular" begin
        rng = MersenneTwister(123456)
        D = 13
        A = randn(rng, D, D)
        U = upper_triangular(A)

        @test U == UpperTriangular(A)
        @test first(Zygote.gradient(A->sum(upper_triangular(A)), A)) isa Matrix
        adjoint_test(upper_triangular, randn(rng, D, D), A)
    end
    @testset "chol" begin
        rng = MersenneTwister(123456)
        D = 3

        # Check correctness of chol. Should be trivial.
        U_A = to_positive_definite_softplus(UpperTriangular(randn(rng, D, D)))
        A = Symmetric(U_A'U_A)
        @test chol(A) ≈ cholesky(A).U

        # Check correctness of frule for chol.
        ΔU_A = to_positive_definite_softplus(UpperTriangular(randn(rng, D, D)))
        ΔA = Symmetric(ΔU_A'ΔU_A)
        frule_test(chol, (A, ΔA); fdm=forward_fdm(5, 1))
    end
    @testset "inv(::Cholesky)" begin
        rng = MersenneTwister(123456)
        D = 3

        A = randn(rng, D, D)
        adjoint_test(A->inv(cholesky(Symmetric(A'A + I))), randn(rng, D, D), A)
    end
    # @testset "chol_inv(::AbstractMatrix{<:Real})" begin
    #     rng = MersenneTwister(123456)
    #     D = 3

    #     A = randn(rng, D, D)
    #     S = Symmetric(A'A + I)
    #     ΔS = randn(rng, D, D)

    #     S_inv, back_inv = Zygote.pullback(inv, S)
    #     Δ_inv = first(back_inv(ΔS))

    #     S_chol_inv, back_chol_inv = Zygote.pullback(VIEstimators.chol_inv, S)
    #     Δ_chol_inv = first(back_chol_inv(ΔS))

    #     @test S_inv ≈ S_chol_inv
    #     @test Δ_inv ≈ Δ_chol_inv
    # end
end


text:
# A quadratic quadrilateral in 3D space.

struct Quadrilateral8_3D{T <: AbstractFloat}
    # The points are assumed to be ordered  in counter clockwise order as follows
    # p₁ = vertex A
    # p₂ = vertex B
    # p₃ = vertex C
    # p₄ = vertex D
    # p₅ = point on the quadratic segment from A to B
    # p₆ = point on the quadratic segment from B to C
    # p₇ = point on the quadratic segment from C to D
    # p₈ = point on the quadratic segment from D to A
    points::NTuple{8, Point_3D{T}}
end

# Constructors
# -------------------------------------------------------------------------------------------------
Quadrilateral8_3D(p₁::Point_3D{T}, p₂::Point_3D{T}, p₃::Point_3D{T}, p₄::Point_3D{T},
                  p₅::Point_3D{T}, p₆::Point_3D{T}, p₇::Point_3D{T}, p₈::Point_3D{T}
        ) where {T <: AbstractFloat} = Quadrilateral8_3D((p₁, p₂, p₃, p₄, p₅, p₆, p₇, p₈))

# Methods
# -------------------------------------------------------------------------------------------------
function (quad8::Quadrilateral8_3D{T})(r::R, s::S) where {T <: AbstractFloat,
                                                          R <: Real,
                                                          S <: Real}
    # See The Visualization Toolkit: An Object-Oriented Approach to 3D Graphics, 4th Edition
    # Chapter 8, Advanced Data Representation, in the interpolation functions section
    ξ = 2T(r) - 1; η = 2T(s) - 1
    return (1 - ξ)*(1 - η)*(-ξ - η - 1)/4*quad8.points[1] +
           (1 + ξ)*(1 - η)*( ξ - η - 1)/4*quad8.points[2] +
           (1 + ξ)*(1 + η)*( ξ + η - 1)/4*quad8.points[3] +
           (1 - ξ)*(1 + η)*(-ξ + η - 1)/4*quad8.points[4] +
                      (1 - ξ^2)*(1 - η)/2*quad8.points[5] +
                      (1 - η^2)*(1 + ξ)/2*quad8.points[6] +
                      (1 - ξ^2)*(1 + η)/2*quad8.points[7] +
                      (1 - η^2)*(1 - ξ)/2*quad8.points[8]
end

function derivatives(quad8::Quadrilateral8_3D{T}, r::R, s::S) where {T <: AbstractFloat,
                                                                     R <: Real,
                                                                     S <: Real}
    # Chain rule
    # ∂Q   ∂Q ∂ξ  ∂Q   ∂Q ∂η
    # -- = -- --, -- = -- --
    # ∂r   ∂ξ ∂r  ∂s   ∂η ∂s
    ξ = 2T(r) - 1; η = 2T(s) - 1
    ∂Q_∂ξ = (1 - η)*(2ξ + η)/4*quad8.points[1] +
            (1 - η)*(2ξ - η)/4*quad8.points[2] +
            (1 + η)*(2ξ + η)/4*quad8.points[3] +
            (1 + η)*(2ξ - η)/4*quad8.points[4] +
                    -ξ*(1 - η)*quad8.points[5] +
                   (1 - η^2)/2*quad8.points[6] +
                    -ξ*(1 + η)*quad8.points[7] +
                  -(1 - η^2)/2*quad8.points[8]

    ∂Q_∂η = (1 - ξ)*( ξ + 2η)/4*quad8.points[1] +
            (1 + ξ)*(-ξ + 2η)/4*quad8.points[2] +
            (1 + ξ)*( ξ + 2η)/4*quad8.points[3] +
            (1 - ξ)*(-ξ + 2η)/4*quad8.points[4] +
                   -(1 - ξ^2)/2*quad8.points[5] +
                     -η*(1 + ξ)*quad8.points[6] +
                    (1 - ξ^2)/2*quad8.points[7] +
                     -η*(1 - ξ)*quad8.points[8]

    return 2*∂Q_∂ξ, 2*∂Q_∂η
end

function area(quad8::Quadrilateral8_3D{T}; N::Int64=15) where {T <: AbstractFloat}
    # Numerical integration required. Gauss-Legendre quadrature over a quadrilateral is used.
    # Let Q(r,s) be the interpolation function for quad8,
    #                             1  1
    # A = ∬ ||∂Q/∂r × ∂Q/∂s||dA = ∫  ∫ ||∂Q/∂r × ∂Q/∂s|| ds dr
    #      D                      0  0
    #
    #       N   N
    #   =   ∑   ∑  wᵢwⱼ||∂Q/∂r(rᵢ,sⱼ) × ∂Q/∂s(rᵢ,sⱼ)||
    #      i=1 j=1
    # N is the square root of the number of points used in the quadrature.
    # See tuning/Quadrilateral8_3D_area.jl for more info on how N = 15 was chosen.
    w, r = gauss_legendre_quadrature(T, N)
    a = T(0)
    for i = 1:N, j = 1:N
        ∂Q_∂r, ∂Q_∂s = derivatives(quad8, r[i], r[j])
        a += w[i]*w[j]*norm(∂Q_∂r × ∂Q_∂s)
    end
    return a
end

function triangulate(quad8::Quadrilateral8_3D{T}, N::Int64) where {T <: AbstractFloat}
    # N is the number of divisions of each edge
    triangles = Vector{Triangle_3D{T}}(undef, 2*(N+1)*(N+1))
    if N === 0
        triangles[1] = Triangle_3D(quad8.points[1], quad8.points[2], quad8.points[3])
        triangles[2] = Triangle_3D(quad8.points[3], quad8.points[4], quad8.points[1])
    else
        for j = 0:N, i = 0:N
            triangles[2*(N+1)*j + 2i + 1] = Triangle_3D(quad8(    i/(N+1),     j/(N+1)),
                                                        quad8((i+1)/(N+1),     j/(N+1)),
                                                        quad8(    i/(N+1), (j+1)/(N+1)))
            triangles[2*(N+1)*j + 2i + 2] = Triangle_3D(quad8(    i/(N+1), (j+1)/(N+1)),
                                                        quad8((i+1)/(N+1),     j/(N+1)),
                                                        quad8((i+1)/(N+1), (j+1)/(N+1)))
        end
    end
    return triangles
end

function intersect(l::LineSegment_3D{T}, quad8::Quadrilateral8_3D{T}; N::Int64 = 20) where {T <: AbstractFloat}
    # Triangulate and intersect each triangle
    triangles = triangulate(quad8, N)
    npoints = 0
    p₁ = Point_3D(T, 0)
    p₂ = Point_3D(T, 0)
    intersections = l .∩ triangles
    bools = map(x->x[1], intersections)
    points = map(x->x[2], intersections)
    npoints = count(bools)
    if npoints === 0
        return false, 0, p₁, p₂
    elseif npoints === 1
        p₁ = points[argmax(bools)]
        return true, 1, p₁, p₂
    elseif npoints === 2
        indices = findall(bools)
        p₁ = points[indices[1]]
        p₂ = points[indices[2]]
        # Check uniqueness
        if p₁ ≈ p₂
            return true, 1, p₁, p₂
        else
            return true, 2, p₁, p₂
        end
    else
        # Account for 3 points and 4 points?
        # If intersection is on edge shared by two triangles on entrance and/or exit 3/4 intersections
        # can be detected
        return true, -1, p₁, p₂
    end
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
function run(world::VLGameWorld)::VLMarketSimulationResult

    try

        # get data about the world -
        agent_array = world.agent_array                                         # get the agent array
        asset_price_array = world.asset_price_array                             # get the asset price array
        number_of_iterations = world.number_of_iterations                       # how iterations are we going to do -
        number_of_agents = length(agent_array)
        (number_of_iterations, number_of_assets) = size(asset_price_array)

        # initialize simulation data structures -
        agent_desired_orders_array = Array{VLOrderModel,1}()

        # main game loop -
        for iteration_index = 2:number_of_iterations

            # for each agent, we need to pass the current price array to the logic function -
            for agent_index = 1:number_of_agents

                # get the agent model -
                agent_model = agent_array[agent_index]

                # what is the desired position size in the next time step?
                # desired_position_size_array: number of assets x 1 array
                tmp_desired_order_array = agent_model.agent_trade_logic(agent_model, asset_price_array, iteration_index)
                for asset_index = 1:number_of_assets
                    tmp_order = tmp_desired_order_array[asset_index]
                    push!(agent_desired_orders_array, tmp_order)
                end
            end

            # ok, so let's post the desired trades to the exchange -> and get the exachange output
            asset_price_array = world.exchange_logic(iteration_index, world, agent_desired_orders_array, asset_price_array)

            @show asset_price_array

            # clean-up for next iteration -
            empty!(agent_desired_orders_array)
        end

        # build result object -
        result_object = VLMarketSimulationResult(agent_array, asset_price_array)

        # return -
        return result_object
    catch error
        rethrow(error)
    end
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: mwaa
using AWS.Compat
using AWS.UUIDs

"""
    CreateCliToken()

Create a CLI token to use Airflow CLI.

# Required Parameters
- `Name`: Create a CLI token request for a MWAA environment.

"""
create_cli_token(Name; aws_config::AWSConfig=global_aws_config()) = mwaa("POST", "/clitoken/$(Name)"; aws_config=aws_config)
create_cli_token(Name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mwaa("POST", "/clitoken/$(Name)", args; aws_config=aws_config)

"""
    CreateEnvironment()

JSON blob that describes the environment to create.

# Required Parameters
- `DagS3Path`: The relative path to the DAG folder on your Amazon S3 storage bucket. For example, dags. For more information, see Importing DAGs on Amazon MWAA.
- `ExecutionRoleArn`: The Amazon Resource Name (ARN) of the execution role for your environment. An execution role is an AWS Identity and Access Management (IAM) role that grants MWAA permission to access AWS services and resources used by your environment. For example, arn:aws:iam::123456789:role/my-execution-role. For more information, see Managing access to Amazon Managed Workflows for Apache Airflow.
- `Name`: The name of your MWAA environment.
- `NetworkConfiguration`: The VPC networking components you want to use for your environment. At least two private subnet identifiers and one VPC security group identifier are required to create an environment. For more information, see Creating the VPC network for a MWAA environment.
- `SourceBucketArn`: The Amazon Resource Name (ARN) of your Amazon S3 storage bucket. For example, arn:aws:s3:::airflow-mybucketname.

# Optional Parameters
- `AirflowConfigurationOptions`: The Apache Airflow configuration setting you want to override in your environment. For more information, see Environment configuration.
- `AirflowVersion`: The Apache Airflow version you want to use for your environment.
- `EnvironmentClass`: The environment class you want to use for your environment. The environment class determines the size of the containers and database used for your Apache Airflow services.
- `KmsKey`: The AWS Key Management Service (KMS) key to encrypt and decrypt the data in your environment. You can use an AWS KMS key managed by MWAA, or a custom KMS key (advanced). For more information, see Customer master keys (CMKs) in the AWS KMS developer guide.
- `LoggingConfiguration`: The Apache Airflow logs you want to send to Amazon CloudWatch Logs.
- `MaxWorkers`: The maximum number of workers that you want to run in your environment. MWAA scales the number of Apache Airflow workers and the Fargate containers that run your tasks up to the number you specify in this field. When there are no more tasks running, and no more in the queue, MWAA disposes of the extra containers leaving the one worker that is included with your environment.
- `PluginsS3ObjectVersion`: The plugins.zip file version you want to use.
- `PluginsS3Path`: The relative path to the plugins.zip file on your Amazon S3 storage bucket. For example, plugins.zip. If a relative path is provided in the request, then PluginsS3ObjectVersion is required. For more information, see Importing DAGs on Amazon MWAA.
- `RequirementsS3ObjectVersion`: The requirements.txt file version you want to use.
- `RequirementsS3Path`: The relative path to the requirements.txt file on your Amazon S3 storage bucket. For example, requirements.txt. If a relative path is provided in the request, then RequirementsS3ObjectVersion is required. For more information, see Importing DAGs on Amazon MWAA.
- `Tags`: The metadata tags you want to attach to your environment. For more information, see Tagging AWS resources.
- `WebserverAccessMode`: The networking access of your Apache Airflow web server. A public network allows your Airflow UI to be accessed over the Internet by users granted access in your IAM policy. A private network limits access of your Airflow UI to users within your VPC. For more information, see Creating the VPC network for a MWAA environment.
- `WeeklyMaintenanceWindowStart`: The day and time you want MWAA to start weekly maintenance updates on your environment.
"""
create_environment(DagS3Path, ExecutionRoleArn, Name, NetworkConfiguration, SourceBucketArn; aws_config::AWSConfig=global_aws_config()) = mwaa("PUT", "/environments/$(Name)", Dict{String, Any}("DagS3Path"=>DagS3Path, "ExecutionRoleArn"=>ExecutionRoleArn, "NetworkConfiguration"=>NetworkConfiguration, "SourceBucketArn"=>SourceBucketArn); aws_config=aws_config)
create_environment(DagS3Path, ExecutionRoleArn, Name, NetworkConfiguration, SourceBucketArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mwaa("PUT", "/environments/$(Name)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("DagS3Path"=>DagS3Path, "ExecutionRoleArn"=>ExecutionRoleArn, "NetworkConfiguration"=>NetworkConfiguration, "SourceBucketArn"=>SourceBucketArn), args)); aws_config=aws_config)

"""
    CreateWebLoginToken()

Create a JWT token to be used to login to Airflow Web UI with claims based Authentication.

# Required Parameters
- `Name`: Create an Airflow Web UI login token request for a MWAA environment.

"""
create_web_login_token(Name; aws_config::AWSConfig=global_aws_config()) = mwaa("POST", "/webtoken/$(Name)"; aws_config=aws_config)
create_web_login_token(Name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mwaa("POST", "/webtoken/$(Name)", args; aws_config=aws_config)

"""
    DeleteEnvironment()

Delete an existing environment.

# Required Parameters
- `Name`: The name of the environment to delete.

"""
delete_environment(Name; aws_config::AWSConfig=global_aws_config()) = mwaa("DELETE", "/environments/$(Name)"; aws_config=aws_config)
delete_environment(Name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mwaa("DELETE", "/environments/$(Name)", args; aws_config=aws_config)

"""
    GetEnvironment()

Get details of an existing environment.

# Required Parameters
- `Name`: The name of the environment to retrieve.

"""
get_environment(Name; aws_config::AWSConfig=global_aws_config()) = mwaa("GET", "/environments/$(Name)"; aws_config=aws_config)
get_environment(Name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mwaa("GET", "/environments/$(Name)", args; aws_config=aws_config)

"""
    ListEnvironments()

List Amazon MWAA Environments.

# Optional Parameters
- `MaxResults`: The maximum results when listing MWAA environments.
- `NextToken`: The Next Token when listing MWAA environments.
"""
list_environments(; aws_config::AWSConfig=global_aws_config()) = mwaa("GET", "/environments"; aws_config=aws_config)
list_environments(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = mwaa("GET", "/environments", args; aws_config=aws_config)

"""
    ListTagsForResource()

List the tags for MWAA environments.

# Required Parameters
- `ResourceArn`: The ARN of the MWAA environment.

"""
list_tags_for_resource(ResourceArn; aws_config::AWSConfig=global_aws_config()) = mwaa("GET", "/tags/$(ResourceArn)"; aws_config=aws_config)
list_tags_for_resource(ResourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mwaa("GET", "/tags/$(ResourceArn)", args; aws_config=aws_config)

"""
    PublishMetrics()

An operation for publishing metrics from the customers to the Ops plane.

# Required Parameters
- `EnvironmentName`: Publishes environment metric data to Amazon CloudWatch.
- `MetricData`: Publishes metric data points to Amazon CloudWatch. CloudWatch associates the data points with the specified metrica.

"""
publish_metrics(EnvironmentName, MetricData; aws_config::AWSConfig=global_aws_config()) = mwaa("POST", "/metrics/environments/$(EnvironmentName)", Dict{String, Any}("MetricData"=>MetricData); aws_config=aws_config)
publish_metrics(EnvironmentName, MetricData, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mwaa("POST", "/metrics/environments/$(EnvironmentName)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("MetricData"=>MetricData), args)); aws_config=aws_config)

"""
    TagResource()

Add tag to the MWAA environments.

# Required Parameters
- `ResourceArn`: The tag resource ARN of the MWAA environments.
- `Tags`: The tag resource tag of the MWAA environments.

"""
tag_resource(ResourceArn, Tags; aws_config::AWSConfig=global_aws_config()) = mwaa("POST", "/tags/$(ResourceArn)", Dict{String, Any}("Tags"=>Tags); aws_config=aws_config)
tag_resource(ResourceArn, Tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mwaa("POST", "/tags/$(ResourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Tags"=>Tags), args)); aws_config=aws_config)

"""
    UntagResource()

Remove a tag from the MWAA environments.

# Required Parameters
- `ResourceArn`: The tag resource ARN of the MWAA environments.
- `tagKeys`: The tag resource key of the MWAA environments.

"""
untag_resource(ResourceArn, tagKeys; aws_config::AWSConfig=global_aws_config()) = mwaa("DELETE", "/tags/$(ResourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(ResourceArn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mwaa("DELETE", "/tags/$(ResourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), args)); aws_config=aws_config)

"""
    UpdateEnvironment()

Update an MWAA environment.

# Required Parameters
- `Name`: The name of your Amazon MWAA environment that you wish to update.

# Optional Parameters
- `AirflowConfigurationOptions`: The Airflow Configuration Options to update of your Amazon MWAA environment.
- `AirflowVersion`: The Airflow Version to update of your Amazon MWAA environment.
- `DagS3Path`: The Dags folder S3 Path to update of your Amazon MWAA environment.
- `EnvironmentClass`: The Environment Class to update of your Amazon MWAA environment.
- `ExecutionRoleArn`: The Executio Role ARN to update of your Amazon MWAA environment.
- `LoggingConfiguration`: The Logging Configuration to update of your Amazon MWAA environment.
- `MaxWorkers`: The Maximum Workers to update of your Amazon MWAA environment.
- `NetworkConfiguration`: The Network Configuration to update of your Amazon MWAA environment.
- `PluginsS3ObjectVersion`: The Plugins.zip S3 Object Version to update of your Amazon MWAA environment.
- `PluginsS3Path`: The Plugins.zip S3 Path to update of your Amazon MWAA environment.
- `RequirementsS3ObjectVersion`: The Requirements.txt S3 ObjectV ersion to update of your Amazon MWAA environment.
- `RequirementsS3Path`: The Requirements.txt S3 Path to update of your Amazon MWAA environment.
- `SourceBucketArn`: The S3 Source Bucket ARN to update of your Amazon MWAA environment.
- `WebserverAccessMode`: The Webserver Access Mode to update of your Amazon MWAA environment.
- `WeeklyMaintenanceWindowStart`: The Weekly Maintenance Window Start to update of your Amazon MWAA environment.
"""
update_environment(Name; aws_config::AWSConfig=global_aws_config()) = mwaa("PATCH", "/environments/$(Name)"; aws_config=aws_config)
update_environment(Name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mwaa("PATCH", "/environments/$(Name)", args; aws_config=aws_config)


text:
import NLsolve
import HDF5
import QuadGK

# --- Influence Coefficients ---

"""
applies for both Ay and Rx depending on which function ifunc(x, y, phi)
is passed in
"""
function panelIntegration(xvec, yvec, thetavec, ifunc)

    # initialize
    nx = length(xvec)
    ntheta = length(thetavec)
    dtheta = thetavec[2] - thetavec[1]  # assumes equally spaced
    A = zeros(nx, ntheta)

    for i in eachindex(xvec)
        # redefine function so it has one parameter for use in quadgk
        integrand(phi) = ifunc(xvec[i], yvec[i], phi)

        for j in eachindex(thetavec)
            # an Adaptive Gauss-Kronrod quadrature integration.  Tried trapz but this was faster.
            A[i, j], error = QuadGK.quadgk(integrand, thetavec[j]-dtheta/2.0, thetavec[j]+dtheta/2.0, atol=1e-10)
        end

    end

    return A
end

"""
integrand used for computing Dx
"""
function Dxintegrand(x, y, phi)
    v1 = x + sin(phi)
    v2 = y - cos(phi)
    # v1 and v2 should never both be zero b.c. we never integrate self.  RxII handles that case.
    return (v1*sin(phi) - v2*cos(phi))/(2*pi*(v1*v1 + v2*v2))
end

"""
integrand used for computing Ay
"""
function Ayintegrand(x, y, phi)
    v1 = x + sin(phi)
    v2 = y - cos(phi)
    if abs(v1) < 1e-12 && abs(v2) < 1e-12  # occurs when integrating self, function symmetric around singularity, should integrate to zero
        return 0.0
    end
    return (v1*cos(phi) + v2*sin(phi))/(2*pi*(v1*v1 + v2*v2))
end

function AyIJ(xvec, yvec, thetavec)
    return panelIntegration(xvec, yvec, thetavec, Ayintegrand)
end

function DxIJ(xvec, yvec, thetavec)
    return panelIntegration(xvec, yvec, thetavec, Dxintegrand)
end

function WxIJ(xvec, yvec, thetavec)

    # initialize
    nx = length(xvec)
    ntheta = length(thetavec)
    dtheta = thetavec[2] - thetavec[1]  # assumes equally spaced
    Wx = zeros(nx, ntheta)

    for i in eachindex(xvec)
        if yvec[i] >= -1.0 && yvec[i] <= 1.0 && xvec[i] >= 0.0 && xvec[i]^2 + yvec[i]^2 >= 1.0
        # if yvec[i] >= -1.0 && yvec[i] <= 1.0 && (xvec[i] >= 0.0 || (xvec[i] >= -1 && xvec[i]^2 + yvec[i]^2 <= 1.0))
            thetak = acos(yvec[i])
            k = findfirst(thetavec + dtheta/2 .> thetak)  # index of intersection
            Wx[i, k] = -1.0
            Wx[i, ntheta-k+1] = 1.0
        end
    end

    return Wx
end

function DxII(thetavec)

    # initialize
    ntheta = length(thetavec)
    dtheta = thetavec[2] - thetavec[1]  # assumes equally spaced
    Rx = dtheta/(4*pi)*ones(ntheta, ntheta)

    for i in eachindex(thetavec)
        if i <= ntheta/2
            Rx[i, i] = (-1 + 1.0/ntheta)/2.0
        else
            Rx[i, i] = (1 + 1.0/ntheta)/2.0
        end
    end

    return Rx
end

function WxII(thetavec)

    # initialize
    ntheta = length(thetavec)
    Wx = zeros(ntheta, ntheta)

    for i = div(ntheta,2)+1:ntheta
        Wx[i, ntheta+1-i] = -1
    end

    return Wx
end

function precomputeMatrices(ntheta)

    # precompute self influence matrices

    # setup discretization (all the same, and uniformly spaced in theta)
    dtheta = 2*pi/ntheta
    theta = collect(dtheta/2:dtheta:2*pi)

    Dxself = DxII(theta)
    Wxself = WxII(theta)
    Ayself = AyIJ(-sin.(theta), cos.(theta), theta)

    # write to file
    HDF5.h5open("$modulepath/theta-$ntheta.h5", "w") do file
        HDF5.write(file, "theta", theta)
        HDF5.write(file, "Dx", Dxself)
        HDF5.write(file, "Wx", Wxself)
        HDF5.write(file, "Ay", Ayself)
    end

end

function matrixAssemble(centerX, centerY, radii, ntheta)
    """
    centerX, centerY: array of x,y coordinates for centers of the VAWTs in the farm
    radii: corresponding array of their radii
    """

    file = "$modulepath/theta-$ntheta.h5"
    if !isfile(file)
        precomputeMatrices(ntheta)
    end

    theta = HDF5.h5read(file, "theta")
    Dxself = HDF5.h5read(file, "Dx")
    Wxself = HDF5.h5read(file, "Wx")
    Ayself = HDF5.h5read(file, "Ay")

    # initialize global matrices
    nturbines = length(radii)
    Dx = zeros(nturbines*ntheta, nturbines*ntheta)
    Wx = zeros(nturbines*ntheta, nturbines*ntheta)
    Ay = zeros(nturbines*ntheta, nturbines*ntheta)

    # iterate through turbines
    for I in eachindex(radii)
        for J in eachindex(radii)

            # find normalized i locations relative to center of turbine J
            x = (centerX[I].-radii[I]*sin.(theta) .- centerX[J])/radii[J]
            y = (centerY[I].+radii[I]*cos.(theta) .- centerY[J])/radii[J]

            # self-influence is precomputed
            if I == J
                Dxsub = Dxself
                Wxsub = Wxself
                Aysub = Ayself

            # pairs can be mapped for same radius
            elseif J < I && radii[I] == radii[J]

                # grab cross-diagonal I,J -> J,I matrix
                Dxsub = Dx[(J-1)*ntheta+1:J*ntheta, (I-1)*ntheta+1:I*ntheta]
                Aysub = Ay[(J-1)*ntheta+1:J*ntheta, (I-1)*ntheta+1:I*ntheta]

                # mapping index for coefficients that are the same
                idx = [div(ntheta,2)+1:ntheta; 1:div(ntheta,2)]

                # directly map over
                Dxsub = Dxsub[idx, idx]
                Aysub = Aysub[idx, idx]

                # wake term must be recomptued
                Wxsub = WxIJ(x, y, theta)

            # # if VAWTs are very far apart we can approximate some of the influence coefficients
            # elseif approxfar && sqrt((centerX[I]-centerX[J])^2 + (centerY[I]-centerY[J])^2) > 10*radii[I]
            #     println("far apart")
            #     xc = (centerX[I] - centerX[J])/radii[J]
            #     yc = (centerY[I] - centerY[J])/radii[J]

            #     Rxsub = RxIJFar(xc, yc, theta)
            #     Wxsub = zeros(ntheta, ntheta)  # should have negligible wake impact
            #     Aysub = AyIJFar(xc, yc, theta)

            else
                Dxsub = DxIJ(x, y, theta)
                Wxsub = WxIJ(x, y, theta)
                Aysub = AyIJ(x, y, theta)
            end

            # assemble into global matrix
            Dx[(I-1)*ntheta+1:I*ntheta, (J-1)*ntheta+1:J*ntheta] = Dxsub
            Wx[(I-1)*ntheta+1:I*ntheta, (J-1)*ntheta+1:J*ntheta] = Wxsub
            Ay[(I-1)*ntheta+1:I*ntheta, (J-1)*ntheta+1:J*ntheta] = Aysub

        end
    end
    Ax = Dx + Wx

    return Ax, Ay, theta
end
# ---------------------------

# ------- Force Coefficients ---------

mutable struct Turbine{TF,TI,TFN}
    r::TF
    chord::TF
    twist::TF
    delta::TF
    B::TI
    af::TFN
    Omega::TF
    centerX::TF
    centerY::TF
end

mutable struct Environment{TF}
    Vinf::TF
    rho::TF
    mu::TF
end

function radialforce(uvec, vvec, thetavec,
    turbine::Turbine, env)
    # u, v, theta - arrays of size ntheta
    # r, chord, twist, Vinf, Omega, rho, mu - scalars

    # unpack
    r = turbine.r
    chord = turbine.chord
    twist = turbine.twist
    delta = turbine.delta
    B = turbine.B
    Omega = turbine.Omega

    Vinf = env.Vinf
    rho = env.rho

    # set the rotation direction
    rotation = sign(Omega)

    # velocity components and angles
    Vn = Vinf*(1.0 .+ uvec).*sin.(thetavec) .- Vinf*vvec.*cos.(thetavec)
    Vt = rotation*(Vinf*(1.0 .+ uvec).*cos.(thetavec) .+ Vinf*vvec.*sin.(thetavec)) .+ abs(Omega)*r
    W = sqrt.(Vn.^2 .+ Vt.^2)
    phi = atan.(Vn, Vt)
    alpha = phi .- twist
    # Re = rho*W*chord/mu  # currently no Re dependence

    # airfoil
    # cl = turbine.af.cl[alpha]
    # cd = turbine.af.cd[alpha]
    cl, cd = turbine.af(alpha)

    # rotate force coefficients
    cn = cl.*cos.(phi) .+ cd.*sin.(phi)
    ct = cl.*sin.(phi) .- cd.*cos.(phi)

    # radial force
    sigma = B*chord/r
    q = sigma/(4*pi)*cn.*(W/Vinf).^2

    # instantaneous forces
    qdyn = 0.5*rho*W.^2
    Rp = -cn.*qdyn*chord
    Tp = ct.*qdyn*chord/cos(delta)
    Zp = -cn.*qdyn*chord*tan(delta)

    # nonlinear correction factor
    integrand = (W/Vinf).^2 .* (cn.*sin.(thetavec) .- rotation*ct.*cos.(thetavec)/cos(delta))
    CT = sigma/(4*pi) * pInt(thetavec, integrand)
    if CT > 2.0
        a = 0.5*(1.0 + sqrt(1.0 + CT))
        ka = 1.0 / (a-1)

    elseif CT > 0.96
        a = 1.0/7*(1 + 3.0*sqrt(7.0/2*CT - 3))
        ka = 18.0*a / (7*a^2 - 2*a + 4)

    else
        a = 0.5*(1 - sqrt(1.0 - CT))
        ka = 1.0 / (1-a)
    end

    # power coefficient
    H = 1.0  # per unit height
    Sref = 2*r*H
    Q = r*Tp
    P = abs(Omega)*B/(2*pi)*pInt(thetavec, Q)
    CP = P / (0.5*rho*Vinf^3 * Sref)

    return q, ka, CT, CP, Rp, Tp, Zp
end

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

# ------ Solve System --------------

function residual(w, A, theta,
    k, turbines, env)

    # setup
    ntheta = length(theta)
    nturbines = Int(length(w)/2/ntheta)
    q = zeros(ntheta*nturbines)
    ka = 0.0

    for i = 1:nturbines
        idx = collect((i-1)*ntheta+1:i*ntheta)

        u = w[idx]
        v = w[ntheta*nturbines .+ idx]

        q[idx], ka, _, _, _, _, _ = radialforce(u, v, theta, turbines[i], env)
    end

    if nturbines == 1  # if only one turbine use the k from the analysis
        k = [ka]
    end  # otherwise, use k that was input to this function

    # reformat to multiply in correct locations
    kmult = repeat(k, inner=[ntheta])
    kmult = [kmult; kmult]

    return (A*q).*kmult - w
end

function actuatorcylinder(turbines, env, ntheta)

    # list comprehensions
    centerX = [turbine.centerX for turbine in turbines]
    centerY = [turbine.centerY for turbine in turbines]
    radii = [turbine.r for turbine in turbines]

    # assemble global matrices
    Ax, Ay, theta = matrixAssemble(centerX, centerY, radii, ntheta)

    # setup
    ntheta = length(theta)
    nturbines = length(turbines)
    tol = 1e-6
    CT = zeros(nturbines)
    CP = zeros(nturbines)
    Rp = zeros(ntheta, nturbines)
    Tp = zeros(ntheta, nturbines)
    Zp = zeros(ntheta, nturbines)

    q = zeros(ntheta)

    # compute nonlinear correction factors (each turbine individaully)
    k = zeros(nturbines)

    for i = 1:nturbines
        w0 = zeros(ntheta*2)

        idx = collect((i-1)*ntheta+1:i*ntheta)

        resid_single(x) = residual(x,[Ax[idx, idx]; Ay[idx, idx]], theta, [1.0], turbines, env)
        result = NLsolve.nlsolve(resid_single, w0, ftol=tol)
        w = result.zero
        if !NLsolve.converged(result)
            println("NLsolve terminated prematurely. info = ", info)
        end

        idx = collect(1:ntheta)
        u = w[idx]
        v = w[ntheta .+ idx]
        q, k[i], CT[i], CP[i], Rp[:, i], Tp[:, i], Zp[:, i] = radialforce(u, v, theta, turbines[i], env)

    end

    if nturbines == 1
        return CT, CP, Rp, Tp, Zp, theta
    end

    # Solve coupled system
    w0 = zeros(nturbines*ntheta*2)

    resid_multiple(x) = residual(x,[Ax; Ay], theta, k, turbines, env)
    result = NLsolve.nlsolve(resid_multiple, w0, ftol=tol)
    w = result.zero
    if !NLsolve.converged(result)
        println("NLsolve terminated prematurely. info = ", info)
    end

    for i = 1:nturbines
        idx = collect((i-1)*ntheta+1:i*ntheta)

        u = w[idx]
        v = w[ntheta*nturbines .+ idx]
        _, _, CT[i], CP[i], Rp[:, i], Tp[:, i], Zp[:, i] = radialforce(u, v, theta, turbines[i], env)

    end

    return CT, CP, Rp, Tp, Zp, theta
end

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

# ---------- helper methods --------------

# trapezoidal integration
function trapz(x, y)  # integrate y w.r.t. x

    integral = 0.0
    for i = 1:length(x)-1
        integral += (x[i+1]-x[i])*0.5*(y[i] + y[i+1])
    end
    return integral
end

# integration for a periodic function where end points don't reach ends (uses trapezoidal method)
function pInt(theta, f)

    integral = trapz(theta, f)

    # add end points
    dtheta = 2*theta[1]  # assumes equally spaced, starts at 0
    integral += dtheta * 0.5*(f[1] + f[end])

    return integral
end

max_stars_repo_name: kevmoor/vawt-ac

max_issues_repo_name: kevmoor/vawt-ac

max_forks_repo_name: kevmoor/vawt-ac


text:
__precompile__()

module HermiteQuadratureRules

using #Reduce,
      Compat,
      Base.Cartesian,
     # Distributions,
      SpecialFunctions

export  generate_weights,
        eval_on_nodes,
        herm_design,
        calc_weights,
        rule_lengths,
        gk_nodes,
        weighted_gk_weights

include("rules_grad_free.jl")
#include("rules.jl")
#include("distribution_derivatives.jl")

function __init__()
    setprecision(512)
#    Reduce.Rational(false)
    nothing
end

@generated function normal_moment_general(::Val{N}) where N
    expr = rcall(:(int(x^$N * exp(-x^2/2) / y, x, -infinity, infinity)))
    quote
        y = sqrt(big(2)π)
        $expr
    end
end
function normal_moment_general(N::Int)
    y = big(2)π
    rcall(:(int(x^N * exp(-x^2/2) / y, x, -infinity, infinity)))
end
@generated function normal_moment(::Val{N}) where N
    N%2 == 1 ? 0.0 : normal_moment_general(N)
end
@generated function poly_integral_tuple(::Val{N}) where N
    ntuple( i -> normal_moment_general(i) , Val(N))
end

function eval_poly(x::T, β) where T
    @boundscheck length(β) > 1
    @inbounds out = β[1]
    @inbounds for i ∈ 2:length(β)
        out = fma(out, x, β[i])
    end
    out
end

"""
Order is function val, and then derivatives.
"""
function herm_design(nodes, ::Type{T} = BigFloat) where T
    n = length(nodes)
    X = Matrix{T}(2n, 2n)
    buffer = T.(nodes)
    ∂buffer = similar(buffer)
    for iₙ ∈ 1:n
        X[2iₙ-1,1] = one(T)
        X[2iₙ,1] = zero(T)
    end
    for iₙ ∈ 1:n
        X[2iₙ-1,2] = buffer[iₙ]
        X[2iₙ,2] = one(T)
    end
    for j ∈ 2:2n-1
        ∂buffer .= buffer .* j
        buffer .*= nodes
        for iₙ ∈ 1:n
            X[2iₙ-1,j+1] = buffer[iₙ]
            X[2iₙ,j+1] = ∂buffer[iₙ]
        end
    end
    X
end
function herm_design_sym(nodes::AbstractVector{T}) where T
    n = length(nodes)
    X = Matrix{T}(4n+2, 4n+2)
    buffer = T.(nodes)
    ∂buffer = similar(buffer)
    for iₙ ∈ 1:2n+1
        X[2iₙ-1,1] = one(T)
        X[2iₙ,1] = zero(T)
    end
    for iₙ ∈ 1:n
        X[2iₙ-1,2] = -buffer[n-iₙ+1]
        X[2iₙ,2] = one(T)
    end
    X[2n+1,2] = zero(T)
    X[2n+2,2] = one(T)
    for iₙ ∈ 1:n
        X[2iₙ+1+2n,2] = buffer[iₙ]
        X[2iₙ+2+2n,2] = one(T)
    end
    for j ∈ 2:4n+1
        ∂buffer .= buffer .* j
        buffer .*= nodes
        for iₙ ∈ 1:n
            X[2iₙ-1,j+1] = (-1)^j* buffer[n-iₙ+1]
            X[2iₙ,j+1]  = (-1)^(j-1)*∂buffer[n-iₙ+1]
        end
        X[2n+1,j+1] = zero(T)
        X[2n+2,j+1] = one(T)
        for iₙ ∈ 1:n
            X[2iₙ+1+2n,j+1] = buffer[iₙ]
            X[2iₙ+2+2n,j+1] = ∂buffer[iₙ]
        end
    end
    X
end
function herm_design_sym_reversed_order(nodes::AbstractVector{T}) where T
    n = length(nodes)
    N = 4n+2
    X = Matrix{T}(N,N)
    buffer = Array{T}(uninitialized, n)
    ∂buffer = Array{T}(uninitialized, n)
    for j ∈ 1:N-2
        exponent = N - j
        @. ∂buffer = exponent * nodes ^ (exponent - 1)
        @. buffer = nodes ^ exponent
        for iₙ ∈ 1:n
            X[2iₙ-1,j] = (-1)^exponent * buffer[n-iₙ+1]
            X[2iₙ,j]  = (-1)^(exponent-1) * ∂buffer[n-iₙ+1]
        end
        X[2n+1,j] = zero(T)
        X[2n+2,j] = one(T)
        for iₙ ∈ 1:n
            X[2iₙ+1+2n,j] = buffer[iₙ]
            X[2iₙ+2+2n,j] = ∂buffer[iₙ]
        end
    end
    @. buffer = T(nodes)
    for iₙ ∈ 1:n
        X[2iₙ-1,N-1] = -buffer[n-iₙ+1]
        X[2iₙ,N-1] = one(T)
    end
    X[2n+1,N-1] = zero(T)
    X[2n+2,N-1] = one(T)
    for iₙ ∈ 1:n
        X[2iₙ+1+2n,N-1] = buffer[iₙ]
        X[2iₙ+2+2n,N-1] = one(T)
    end
    for iₙ ∈ 1:2n+1
        X[2iₙ-1,N] = one(T)
        X[2iₙ,N] = zero(T)
    end
    X
end

@generated function weight_mat(Xⁱ::AbstractMatrix{T}, ::Val{N}) where {N,T}
    quote
        weights = Matrix{BigFloat}($N, size(Xⁱ,2))
        for j ∈ 1:size(Xⁱ,2)
            @nexprs $N i -> begin
                weights[i,j] = normal_moment_general(Val{2(i-1)}()) * Xⁱ[2(i-1)+1,j]
            end
        end
        weights
    end
end
@generated function weight_vec(Xⁱ::AbstractMatrix{T}, ::Val{N}) where {N,T}
    quote
        weights = fill(zero(T), size(Xⁱ,2))
        for j ∈ 1:size(Xⁱ,2)
            @nexprs $N i -> begin
                weights[j] += normal_moment_general(Val{2(i-1)}()) * Xⁱ[2(i-1)+1,j]
            end
        end
        weights
    end
end
function weight_vec(Xⁱ::AbstractMatrix{T}) where T#Dynamic dispatch
    weight_vec(Xⁱ, Val{size(Xⁱ,1) ÷ 2}())::Vector{T}
end

function generate_weights(nodes)
    X = herm_design(nodes)
    weight_vec(inv(X))
end

function generate_weights_sym(nodes)
    X = herm_design_sym(nodes)
    weight_vec(inv(X))
end

function eval_on_nodes(f, nodes::AbstractVector{T}) where T
    out = Matrix{T}(2, length(nodes))
    for i ∈ eachindex(nodes)
        fn = f(nodes[i] + eps()*im )
        out[1,i] = imag(fn) / eps()
        out[2,i] = real(fn)
    end
    out
end
# function eval_on_nodes(f::Distributions.UnivariateDistribution, nodes::AbstractVector{T}) where T
#     out = Matrix{T}(2, length(nodes))
#     for i ∈ eachindex(nodes)
#         out[:,i] .= fg(f, nodes[i])
#     end
#     out
# end
# function eval_weighted_nodes(f::Distributions.UnivariateDistribution, nodes::AbstractVector{T}) where T
#     out = Matrix{T}(2, length(nodes))
#     for i ∈ eachindex(nodes)
#         out[:,i] .= fgw(f, nodes[i])
#     end
#     out
# end

# using Optim, ForwardDiff
# initial_n = BigFloat.(sort(rand(4)));
# node_weight_magnitude(nodes) = sum(x -> (x^8-1//(4length(nodes)+2) ), generate_weights_sym(nodes))
# td = TwiceDifferentiable(node_weight_magnitude, initial_n; autodiff = :forward)

# opt = optimize(td, initial_n, Newton())
end # module


text:
import ..Resampling.resample
import ..UncertainValues.AbstractUncertainScalarKDE
import Statistics

#########################################
# Statistics on `AbstractUncertainValueDataset`s
#########################################

"""
    mean(d::AbstractUncertainValueDataset, n::Int)

Computes the element-wise mean of a dataset of uncertain values. Takes
the mean of an `n`-draw sample for each element.
"""
Statistics.mean(d::AbstractUncertainValueDataset, n::Int; kwargs...) =
	[mean(resample(d[i], n), kwargs...) for i = 1:length(d)]

"""
    median(d::AbstractUncertainValueDataset, n::Int)

Computes the element-wise median of a dataset of uncertain values. Takes
the median of an `n`-draw sample for each element.
"""
Statistics.median(d::AbstractUncertainValueDataset, n::Int; kwargs...) =
	[median(resample(d[i], n), kwargs...) for i = 1:length(d)]

"""
    middle(d::AbstractUncertainValueDataset, n::Int)

Compute the middle of `n` realisations of an `AbstractUncertainValueDataset`.
"""
Statistics.middle(d::AbstractUncertainValueDataset, n::Int; kwargs...) =
	[middle(resample(d[i], n), kwargs...) for i = 1:length(d)]

"""
    quantile(d::AbstractUncertainValueDataset, p, n::Int; kwargs...)

Compute element-wise quantile(s) `p `of a dataset consisting of uncertain
values. Takes the quantiles of an `n`-draw sample for each element.
"""
Statistics.quantile(d::AbstractUncertainValueDataset, p, n::Int; kwargs...) =
    [quantile(resample(d[i], n), p, kwargs...) for i = 1:length(d)]

"""
    std(d::AbstractUncertainValueDataset, n::Int; kwargs...)

Computes the element-wise standard deviation of a dataset of uncertain values.
Takes the standard deviation of an `n`-draw sample for each element.
"""
Statistics.std(d::AbstractUncertainValueDataset, n::Int; kwargs...) =
	[std(resample(d[i], n), kwargs...) for i = 1:length(d)]

"""
    var(d::AbstractUncertainValueDataset, n::Int; kwargs...)

Computes the element-wise sample variance of a dataset of uncertain values.
Takes the sample variance of an `n`-draw sample for each element.
"""
Statistics.var(d::AbstractUncertainValueDataset, n::Int; kwargs...) =
	[var(resample(d[i], n), kwargs...) for i = 1:length(d)]

"""
    cor(d1::AbstractUncertainValueDataset, d2::AbstractUncertainValueDataset, n::Int; kwargs...)

Obtain a distribution for the Pearson correlation between two uncertain datasets `d1` and `d2`.
This is done by resampling both datasets multiple times and compute the correlation between
those draws. This yields a distribution of correlation estimates.

The procedure is as follows.

1. First, draw a realisation of `d1` according to the distributions furnishing its uncertain values.
2. Then, draw a realisation `d2` according to its furnishing distributions.
3. Compute the correlation between those two draws/realisations, both of which are vectors of
    length `L = length(d1) = length(d2)`.
4. Repeat the procedure `n` times, drawing `n` separate pairs of realisations of `d1` and `d2`.
    This yields `n` estimates of the correlation between `d1` and `d2`, which is returned as
    a vector.
"""
Statistics.cor(d1::DT, d2::DT, n::Int; kwargs...) where {DT <: AbstractUncertainValueDataset} =
    [cor(resample(d1), resample(d2), kwargs...) for i = 1:n]

"""
    cor(d1::AbstractUncertainValueDataset, d2::AbstractUncertainValueDataset; kwargs...)

Obtain a single estimate of the Pearson correlation between two uncertain datasets `d1` and `d2`.
This is done by resampling both datasets independently, that is: First draw a realisation of
`d1` according to the distributions furnishing its uncertain values. Then, draw a realisation
`d2` according to its furnishing distributions. These those two draws/realisations are now
both vectors of length `L = length(d1) = length(d2)`. Finally, compute the correlation between
those draws. This yields a single estimates of the correlation between `d1` and `d2`.
"""
Statistics.cor(d1::DT, d2::DT, kwargs...) where {DT <: AbstractUncertainValueDataset} =
    cor(resample(d1), resample(d2), kwargs...)

"""
    cov(d1::AbstractUncertainValueDataset, d2::AbstractUncertainValueDataset, n::Int; kwargs...)

Obtain a distribution for the covariance between two uncertain datasets `d1` and `d2`.
This is done by resampling both datasets multiple times and compute the covariance between
those draws. This yields a distribution of covariance estimates.

The procedure is as follows.

1. First, draw a realisation of `d1` according to the distributions furnishing its uncertain values.
2. Then, draw a realisation `d2` according to its furnishing distributions.
3. Compute the covariance between those two draws/realisations, both of which are vectors of
    length `L = length(d1) = length(d2)`.
4. Repeat the procedure `n` times, drawing `n` separate pairs of realisations of `d1` and `d2`.
    This yields `n` estimates of the covariance between `d1` and `d2`, which is returned as
    a vector.
"""
Statistics.cov(d1::DT, d2::DT, n::Int; kwargs...) where DT <: AbstractUncertainValueDataset =
    [cov(resample(d1), resample(d2), kwargs...) for i = 1:n]

"""
    cov(d1::AbstractUncertainValueDataset, d2::AbstractUncertainValueDataset; kwargs...)

Obtain a single estimate for the covariance between two uncertain datasets `d1` and `d2`.
This is done by resampling both datasets independently, that is: First draw a realisation of
`d1` according to the distributions furnishing its uncertain values. Then, draw a realisation
`d2` according to its furnishing distributions. These those two draws/realisations are now
both vectors of length `L = length(d1) = length(d2)`. Finally, compute the covariance between
those draws. This yields a single estimates of the covariance between `d1` and `d2`.
"""
Statistics.cov(d1::DT, d2::DT, kwargs...) where {DT <: AbstractUncertainValueDataset} =
    cov(resample(d1), resample(d2), kwargs...)

export cor, cov, var, std, quantile, mean, median


text:
include("../../utils.jl")
using Flux, Printf, Statistics

# Simple model
function gen_model()
    return Chain(
        Dense(8, 8, σ),
        Dense(8, 8, σ),
        Dense(8, 8, σ),
        Dense(8, 1, σ),
    )
end

# Static dataset
function gen_dataset(batch_size = 128, num_minibatches = 256)
    # Return an array of (x, y) tuples
    return [(
        randn(8, batch_size),
        randn(1, batch_size),
    ) for _ in 1:num_minibatches]
end

# Helper function to print some `@timed` stats
function info_stats(msg, stats, num_epochs)
    @info(
        msg,
        time=stats.time,
        time_per_epoch=stats.time/num_epochs,
        gc=@sprintf("%.1f%%", stats.gctime*100.0/stats.time),
        allocated=Base.format_bytes(stats.bytes),
    )
end

function test()
   # Warm up codegen
#   @info("Warming up Flux.train!()")
   begin
     model = gen_model()
     dataset = gen_dataset()
     opt = Flux.Optimise.ADAM(1e-4)
     num_epochs = 20
     stats = @timed begin
        for idx in 1:num_epochs
            Flux.train!(
                (x, y) -> Flux.Losses.mse(model(x), y),
                Flux.params(model),
                dataset,
                opt,
            )
        end
    end
#    info_stats("First warm-up completed", stats, num_epochs)

    stats = @timed begin
        for idx in 1:4
            Flux.train!(
                (x, y) -> Flux.Losses.mse(model(x), y),
                Flux.params(model),
                dataset,
                opt,
            )
        end
    end
#    info_stats("Second warm-up completed", stats, num_epochs)
   end

   # Train `num_models` in parallel, each on its own thread
   num_models = 32
   datasets = [gen_dataset() for _ in 1:num_models]
   models = [gen_model() for _ in 1:num_models]
   training_stats = [Any[] for _ in 1:num_models]

 #  @warn("Beginning training with $(Threads.nthreads()) threads")
   Threads.@threads for model_idx in 1:num_models
      model = models[model_idx]
      dataset = datasets[model_idx]
      opt = Flux.Optimise.ADAM(1e-4)

    # Threads will die off one by one, so we can see how reducing the number of threads
    # eases off GC pressure over the whole cohort
    for idx in 1:(model_idx*num_epochs)
        push!(training_stats[model_idx], @timed begin
            Flux.train!(
                (x, y) -> Flux.Losses.mse(model(x), y),
                Flux.params(model),
                dataset,
                opt,
            )
        end)
    end

    # Calculate mean statistics over the tail end of this thread's lifetime
    tail_stats = training_stats[model_idx][end-num_epochs+1:end]
    mean_stats = (;
        time = mean(s.time for s in tail_stats),
        gctime = mean(s.gctime for s in tail_stats),
        bytes = mean(s.bytes for s in tail_stats),
    )
#    info_stats("Finished model $(model_idx)", mean_stats, model_idx*num_epochs)
    end
end

@gctime test()

max_stars_repo_name: vtjnash/GCBenchmarks

max_issues_repo_name: vtjnash/GCBenchmarks

max_forks_repo_name: vtjnash/GCBenchmarks


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

include("construction.jl")
include("consistency.jl")
include("numeric.jl")
include("trig.jl")
include("hyperbolic.jl")
include("non_BigFloat.jl")
include("linear_algebra.jl")
include("loops.jl")
include("parsing.jl")
include("rounding_macros.jl")
include("rounding.jl")

# include("complex.jl")


text:
using Pkg

Pkg.activate(pwd())
Pkg.Registry.add(RegistrySpec(url = "https://github.com/JuliaRegistries/General"))
Pkg.update()
Pkg.resolve()

max_stars_repo_name: mattiapenati/web-frameworks

max_issues_repo_name: mattiapenati/web-frameworks

max_forks_repo_name: mattiapenati/web-frameworks


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

"""
    AbstractArrayOfSimilarArrays{T,M,N} <: AbstractArray{<:AbstractArray{T,M},N}

An array that contains arrays that have the same size/axes. The array is
internally stored in flattened form as some kind of array of dimension
`M + N`. The flattened form can be accessed via `flatview(A)`.

Subtypes must implement (in addition to typical array operations):

    flatview(A::SomeArrayOfSimilarArrays)::AbstractArray{T,M+N}

The following type aliases are defined:

* `AbstractVectorOfSimilarArrays{T,M} = AbstractArrayOfSimilarArrays{T,M,1}`
* `AbstractArrayOfSimilarVectors{T,N} = AbstractArrayOfSimilarArrays{T,1,N}`
* `AbstractVectorOfSimilarVectors{T} = AbstractArrayOfSimilarArrays{T,1,1}`
"""
abstract type AbstractArrayOfSimilarArrays{T,M,N} <: AbstractArray{Array{T,M},N} end
export AbstractArrayOfSimilarArrays

const AbstractVectorOfSimilarArrays{T,M} = AbstractArrayOfSimilarArrays{T,M,1}
export AbstractVectorOfSimilarArrays

const AbstractArrayOfSimilarVectors{T,N} = AbstractArrayOfSimilarArrays{T,1,N}
export AbstractArrayOfSimilarVectors

const AbstractVectorOfSimilarVectors{T} = AbstractArrayOfSimilarArrays{T,1,1}
export AbstractVectorOfSimilarVectors

"""
    ArrayOfSimilarArrays{T,M,N,L,P} <: AbstractArrayOfSimilarArrays{T,M,N}

Represents a view of an array of dimension `L = M + N` as an array of
dimension M with elements that are arrays with dimension N. All element arrays
implicitly have equal size/axes.

Constructors:

    ArrayOfSimilarArrays{T,M,N}(flat_data::AbstractArray)
    ArrayOfSimilarArrays{T,M}(flat_data::AbstractArray)

The following type aliases are defined:

* `VectorOfSimilarArrays{T,M} = AbstractArrayOfSimilarArrays{T,M,1}`
* `ArrayOfSimilarVectors{T,N} = AbstractArrayOfSimilarArrays{T,1,N}`
* `VectorOfSimilarVectors{T} = AbstractArrayOfSimilarArrays{T,1,1}`

`VectorOfSimilarArrays` supports `push!()`, etc., provided the underlying
array supports resizing of it's last dimension (e.g. an `ElasticArray`).

The nested array can also be created using the function [`nestedview`](@ref)
and the wrapped flat array can be accessed using [`flatview`](@ref)
afterwards:

```julia
A_flat = rand(2,3,4,5,6)
A_nested = nestedview(A_flat, 2)
A_nested isa AbstractArray{<:AbstractArray{T,2},3} where T
flatview(A_nested) === A_flat
```
"""
struct ArrayOfSimilarArrays{
    T, M, N, L,
    P<:AbstractArray{T,L}
} <: AbstractArrayOfSimilarArrays{T,M,N}
    data::P

    function ArrayOfSimilarArrays{T,M,N}(flat_data::AbstractArray{U,L}) where {T,M,N,L,U}
        size_inner, size_outer = split_tuple(size(flat_data), Val{M}())
        require_ndims(flat_data, _add_vals(Val{M}(), Val{N}()))
        conv_parent = _convert_elype(T, flat_data)
        P = typeof(conv_parent)
        new{T,M,N,L,P}(conv_parent)
    end

    function ArrayOfSimilarArrays{T,M}(flat_data::AbstractArray{U,L}) where {T,M,L,U}
        size_inner, size_outer = split_tuple(size(flat_data), Val{M}())
        N = length(size_outer)
        conv_parent = _convert_elype(T, flat_data)
        P = typeof(conv_parent)
        new{T,M,N,L,P}(conv_parent)
    end
end

export ArrayOfSimilarArrays

function ArrayOfSimilarArrays{T,M,N}(A::AbstractArray{<:AbstractArray{U,M},N}) where {T,M,N,U}
    B = ArrayOfSimilarArrays{T,M,N}(Array{T}(undef, innersize(A)..., size(A)...))
    copyto!(B, A)
end

ArrayOfSimilarArrays{T}(A::AbstractArray{<:AbstractArray{U,M},N}) where {T,M,N,U} =
    ArrayOfSimilarArrays{T,M,N}(A)

ArrayOfSimilarArrays(A::AbstractArray{<:AbstractArray{T,M},N}) where {T,M,N} =
    ArrayOfSimilarArrays{T,M,N}(A)

Base.convert(R::Type{ArrayOfSimilarArrays{T,M,N}}, flat_data::AbstractArray{U,L}) where {T,M,N,L,U} = R(flat_data)
Base.convert(R::Type{ArrayOfSimilarArrays{T,M}}, flat_data::AbstractArray{U,L}) where {T,M,L,U} = R(flat_data)

Base.convert(R::Type{ArrayOfSimilarArrays{T,M,N}}, A::AbstractArray{<:AbstractArray{U,M},N}) where {T,M,N,U} = R(A)
Base.convert(R::Type{ArrayOfSimilarArrays{T}}, A::AbstractArray{<:AbstractArray{U,M},N}) where {T,M,N,U} = R(A)
Base.convert(R::Type{ArrayOfSimilarArrays}, A::AbstractArray{<:AbstractArray{T,M},N}) where {T,M,N} = R(A)

@inline function innersize(A::ArrayOfSimilarArrays{T,M,N}) where {T,M,N}
    front_tuple(size(A.data), Val{M}())
end

@inline function _innerlength(A::AbstractArray{<:AbstractArray{T,M},N}) where {T,M,N}
    prod(innersize(A))
end

import Base.==
(==)(A::ArrayOfSimilarArrays{T,M,N}, B::ArrayOfSimilarArrays{T,M,N}) where {T,M,N} =
    (A.data == B.data)

"""
    flatview(A::ArrayOfSimilarArrays{T,M,N,L,P})::P

Returns the array of dimensionality `L = M + N` wrapped by `A`. The shape of
the result may be freely changed without breaking the inner consistency of
`A`.
"""
flatview(A::ArrayOfSimilarArrays) = A.data

Base.size(A::ArrayOfSimilarArrays{T,M,N}) where {T,M,N} = split_tuple(size(A.data), Val{M}())[2]

Base.@propagate_inbounds Base.getindex(A::ArrayOfSimilarArrays{T,M,N}, idxs::Vararg{Integer,N}) where {T,M,N} =
    view(A.data, _ncolons(Val{M}())..., idxs...)

Base.@propagate_inbounds Base.setindex!(A::ArrayOfSimilarArrays{T,M,N}, x::AbstractArray{U,M}, idxs::Vararg{Integer,N}) where {T,M,N,U} =
    setindex!(A.data, x, _ncolons(Val{M}())..., idxs...)

Base.@propagate_inbounds function Base.unsafe_view(A::ArrayOfSimilarArrays{T,M,N}, idxs::Vararg{Union{Real, AbstractArray},N}) where {T,M,N}
    dataview = view(A.data, _ncolons(Val{M}())..., idxs...)
    L = length(size(dataview))
    N_view = L - M
    ArrayOfSimilarArrays{T,M,N_view}(dataview)
end

@inline function Base.resize!(A::ArrayOfSimilarArrays{T,M,N}, dims::Vararg{Integer,N}) where {T,M,N}
    resize!(A.data, innersize(A)..., dims...)
    A
end

function Base.similar(A::ArrayOfSimilarArrays{T,M,N}, ::Type{<:AbstractArray{U}}, dims::Dims) where {T,M,N,U}
    data = A.data
    size_inner, size_outer = split_tuple(size(data), Val{M}())
    # ToDo: Don't use similar if data is an ElasticArray?
    ArrayOfSimilarArrays{T,M,N}(similar(data, U, size_inner..., dims...))
end

function Base.deepcopy(A::ArrayOfSimilarArrays{T,M,N}) where {T,M,N}
    ArrayOfSimilarArrays{T,M,N}(deepcopy(A.data))
end

function Base.copyto!(dest::ArrayOfSimilarArrays{T,M,N}, src::ArrayOfSimilarArrays{U,M,N}) where {T,M,N,U}
    copyto!(dest.data, src.data)
    dest
end

function Base.append!(dest::ArrayOfSimilarArrays{T,M,N}, src::ArrayOfSimilarArrays{U,M,N}) where {T,M,N,U}
    innersize(dest) != innersize(src) && throw(DimensionMismatch("Can't append, shape of element arrays of source and dest are not equal"))
    append!(dest.data, src.data)
    dest
end

Base.append!(dest::ArrayOfSimilarArrays{T,M,N}, src::AbstractArray{<:AbstractArray{U,M},N}) where {T,M,N,U} =
    append!(dest, ArrayOfSimilarArrays(src))

function Base.prepend!(dest::ArrayOfSimilarArrays{T,M,N}, src::ArrayOfSimilarArrays{U,M,N}) where {T,M,N,U}
    innersize(dest) != innersize(src) && throw(DimensionMismatch("Can't prepend, shape of element arrays of source and dest are not equal"))
    prepend!(dest.data, src.data)
    dest
end

Base.prepend!(dest::ArrayOfSimilarArrays{T,M,N}, src::AbstractArray{<:AbstractArray{U,M},N}) where {T,M,N,U} =
    prepend!(dest, ArrayOfSimilarArrays(src))

UnsafeArrays.unsafe_uview(A::ArrayOfSimilarArrays{T,M,N}) where {T,M,N} =
    ArrayOfSimilarArrays{T,M,N}(uview(A.data))

function innermap(f::Base.Callable, A::ArrayOfSimilarArrays{T,M,N}) where {T,M,N}
    new_data = map(f, A.data)
    U = eltype(new_data)
    ArrayOfSimilarArrays{U,M,N}(new_data)
end

function deepmap(f::Base.Callable, A::ArrayOfSimilarArrays{T,M,N}) where {T,M,N}
    new_data = deepmap(f, A.data)
    U = eltype(new_data)
    ArrayOfSimilarArrays{U,M,N}(new_data)
end

Base.@pure _result_is_nested(idxs_outer::Tuple, idxs_inner::Tuple) =
    Val{!(Base.index_dimsum(idxs_outer...) isa Tuple{}) && !(Base.index_dimsum(idxs_inner...) isa Tuple{})}()

Base.@pure ndims_after_getindex(idxs::Tuple) = Val{length(Base.index_dimsum(idxs...))}()

Base.@propagate_inbounds function deepgetindex(A::ArrayOfSimilarArrays{T,M,N,L}, idxs::Vararg{Any,L}) where {T,M,N,L}
    idxs_outer, idxs_inner = split_tuple(idxs, Val{N}())
    nested = _result_is_nested(idxs_outer, idxs_inner)
    _deepgetindex_impl_aosa(A, idxs_outer, idxs_inner, nested)
end

Base.@propagate_inbounds _deepgetindex_impl_aosa(A::ArrayOfSimilarArrays, idxs_outer::Tuple, idxs_inner::Tuple, nested::Val{false}) =
    getindex(A.data, idxs_inner..., idxs_outer...)

Base.@propagate_inbounds function _deepgetindex_impl_aosa(A::ArrayOfSimilarArrays, idxs_outer::Tuple, idxs_inner::Tuple, nested::Val{true})
    new_data = getindex(A.data, idxs_inner..., idxs_outer...)
    nestedview(new_data, ndims_after_getindex(idxs_inner))
end

Base.@propagate_inbounds function deepsetindex!(A::ArrayOfSimilarArrays{T,M,N,L}, x, idxs::Vararg{Any,L}) where {T,M,N,L}
    idxs_outer, idxs_inner = split_tuple(idxs, Val{N}())
    _deepsetindex_impl_aosa!(A, x, idxs_outer, idxs_inner, _result_is_nested(idxs_outer, idxs_inner))
    A
end

Base.@propagate_inbounds _deepsetindex_impl_aosa!(A::ArrayOfSimilarArrays, x, idxs_outer::Tuple, idxs_inner::Tuple, nested::Val{false}) =
    setindex!(A.data, x, idxs_inner..., idxs_outer...)

Base.@propagate_inbounds _deepsetindex_impl_aosa!(A::ArrayOfSimilarArrays, x, idxs_outer::Tuple, idxs_inner::Tuple, nested::Val{true}) =
    _deepsetindex_impl!(A, x, idxs_outer, idxs_inner)

# TODO: Specialized method _deepsetindex_impl_aosa!(A::ArrayOfSimilarArrays, x::ArrayOfSimilarArrays, idxs_outer::Tuple, idxs_inner::Tuple, nested::Val{true}) =

Base.@propagate_inbounds function deepview(A::ArrayOfSimilarArrays{T,M,N,L}, idxs::Vararg{Any,L}) where {T,M,N,L}
    idxs_outer, idxs_inner = split_tuple(idxs, Val{N}())
    nested = _result_is_nested(idxs_outer, idxs_inner)
    _deepview_impl_aosa(A, idxs_outer, idxs_inner, nested)
end

Base.@propagate_inbounds _deepview_impl_aosa(A::ArrayOfSimilarArrays, idxs_outer::Tuple, idxs_inner::Tuple, nested::Val{false}) =
    view(A.data, idxs_inner..., idxs_outer...)

Base.@propagate_inbounds function _deepview_impl_aosa(A::ArrayOfSimilarArrays, idxs_outer::Tuple, idxs_inner::Tuple, nested::Val{true})
    new_data = view(A.data, idxs_inner..., idxs_outer...)
    nestedview(new_data, ndims_after_getindex(idxs_inner))
end

const VectorOfSimilarArrays{
    T, M, L,
    P<:AbstractArray{T,L}
} = ArrayOfSimilarArrays{T,M,1,L,P}

export VectorOfSimilarArrays

VectorOfSimilarArrays{T}(flat_data::AbstractArray{U,L}) where {T,U,L} =
    ArrayOfSimilarArrays{T,length(Base.front(size(flat_data))),1}(flat_data)

VectorOfSimilarArrays(flat_data::AbstractArray{T,L}) where {T,L} =
    ArrayOfSimilarArrays{T,length(Base.front(size(flat_data))),1}(flat_data)

VectorOfSimilarArrays{T}(A::AbstractVector{<:AbstractArray{U,M}}) where {T,M,U} =
    VectorOfSimilarArrays{T,M}(A)

VectorOfSimilarArrays(A::AbstractVector{<:AbstractArray{T,M}}) where {T,M} =
    VectorOfSimilarArrays{T,M}(A)

Base.convert(R::Type{VectorOfSimilarArrays{T}}, flat_data::AbstractArray{U,L}) where {T,U,L} = R(flat_data)
Base.convert(R::Type{VectorOfSimilarArrays}, flat_data::AbstractArray{T,L}) where {T,L} = R(flat_data)
Base.convert(R::Type{VectorOfSimilarArrays{T}}, A::AbstractVector{<:AbstractArray{U,M}}) where {T,M,U} = R(A)
Base.convert(R::Type{VectorOfSimilarArrays}, A::AbstractVector{<:AbstractArray{T,M}}) where {T,M} = R(A)

@inline Base.IndexStyle(A::VectorOfSimilarArrays) = IndexLinear()

function Base.push!(V::VectorOfSimilarArrays{T,M}, x::AbstractArray{U,M}) where {T,M,U}
    size(x) != Base.front(size(V.data)) && throw(DimensionMismatch("Can't push, shape of source and elements of target is incompatible"))
    append!(V.data, x)
    V
end

function Base.pop!(V::VectorOfSimilarArrays)
    isempty(V) && throw(ArgumentError("array must be non-empty"))
    x = V[end]
    resize!(V, size(V, 1) - 1)
    x
end

function Base.pushfirst!(V::VectorOfSimilarArrays{T,M}, x::AbstractArray{U,M}) where {T,M,U}
    size(x) != Base.front(size(V.data)) && throw(DimensionMismatch("Can't push, shape of source and elements of target is incompatible"))
    prepend!(V.data, x)
    V
end

# Will need equivalent of resize! that resizes in front of data instead of in back:
# popfirst!(V::ArrayOfSimilarArrays) = ...

function Base.empty(A::VectorOfSimilarArrays{T,M}, ::Type{<:AbstractArray{U}}) where {T,M,U}
    data = A.data
    size_inner, size_outer = split_tuple(size(data), Val{1}())
    empty_size_outer = map(x -> zero(x), size_outer)

    # ToDo: Don't use similar if data is an ElasticArray?
    VectorOfSimilarArrays{T,M}(similar(data, U, size_inner..., empty_size_outer...))
end

Base.empty(A::VectorOfSimilarArrays{T,M}) where {T,M} = empty(A, Array{T})

const ArrayOfSimilarVectors{
    T, N, L,
    P<:AbstractArray{T,L}
} = ArrayOfSimilarArrays{T,1,N,L,P}

export ArrayOfSimilarVectors

ArrayOfSimilarVectors{T}(flat_data::AbstractArray{U,L}) where {T,U,L} =
    ArrayOfSimilarArrays{T,1,length(Base.front(size(flat_data)))}(flat_data)

ArrayOfSimilarVectors(flat_data::AbstractArray{T,L}) where {T,L} =
    ArrayOfSimilarArrays{T,1,length(Base.front(size(flat_data)))}(flat_data)

ArrayOfSimilarVectors{T}(A::AbstractArray{<:AbstractVector{U},N}) where {T,N,U} =
    ArrayOfSimilarVectors{T,N}(A)

ArrayOfSimilarVectors(A::AbstractArray{<:AbstractVector{T},N}) where {T,N} =
    ArrayOfSimilarVectors{T,N}(A)

Base.convert(R::Type{ArrayOfSimilarVectors{T}}, flat_data::AbstractArray{U,L}) where {T,U,L} = R(flat_data)
Base.convert(R::Type{ArrayOfSimilarVectors}, flat_data::AbstractArray{T,L}) where {T,L} = R(flat_data)
Base.convert(R::Type{ArrayOfSimilarVectors{T}}, A::AbstractArray{<:AbstractVector{U},N}) where {T,N,U} = R(A)
Base.convert(R::Type{ArrayOfSimilarVectors}, A::AbstractArray{<:AbstractVector{T},N}) where {T,N} = R(A)

const VectorOfSimilarVectors{
    T,
    P<:AbstractArray{T,2}
} = ArrayOfSimilarArrays{T,1,1,2,P}

export VectorOfSimilarVectors

VectorOfSimilarVectors{T}(flat_data::AbstractArray{U,2}) where {T,U} =
    ArrayOfSimilarArrays{T,1,1}(flat_data)

VectorOfSimilarVectors(flat_data::AbstractArray{T,2}) where {T} =
    VectorOfSimilarVectors{T}(flat_data)

VectorOfSimilarVectors{T}(A::AbstractVector{<:AbstractVector{U}}) where {T,U} =
    ArrayOfSimilarArrays{T,1}(A)

VectorOfSimilarVectors(A::AbstractVector{<:AbstractVector{T}}) where {T} =
    VectorOfSimilarVectors{T}(A)

Base.convert(R::Type{VectorOfSimilarVectors{T}}, flat_data::AbstractArray{U,2}) where {T,U} = R(flat_data)
Base.convert(R::Type{VectorOfSimilarVectors}, flat_data::AbstractArray{T,2}) where {T} = R(flat_data)
Base.convert(R::Type{VectorOfSimilarVectors{T}}, A::AbstractVector{<:AbstractVector{U}}) where {T,U} = R(A)
Base.convert(R::Type{VectorOfSimilarVectors}, A::AbstractVector{<:AbstractVector{T}}) where {T} = R(A)

@inline Base.IndexStyle(A::VectorOfSimilarVectors) = IndexLinear()

"""
    sum(X::AbstractVectorOfSimilarArrays)
    sum(X::AbstractVectorOfSimilarArrays, w::StatsBase.AbstractWeights)

Compute the sum of the elements vectors of `X`. Equivalent to `sum` of
`flatview(X)` along the last dimension.
"""
Base.sum(X::AbstractVectorOfSimilarArrays{T,M}) where {T,M} = sum(flatview(X); dims = M + 1)

"""
    mean(X::AbstractVectorOfSimilarArrays)
    mean(X::AbstractVectorOfSimilarArrays, w::StatsBase.AbstractWeights)

Compute the mean of the elements vectors of `X`. Equivalent to `mean` of
`flatview(X)` along the last dimension.
"""
Statistics.mean(X::AbstractVectorOfSimilarArrays{T,M}) where {T,M} =
    mean(flatview(X); dims = M + 1)

"""
    var(X::AbstractVectorOfSimilarArrays; corrected::Bool = true)
    var(X::AbstractVectorOfSimilarArrays, w::StatsBase.AbstractWeights; corrected::Bool = true)

Compute the sample variance of the elements vectors of `X`. Equivalent to
`var` of `flatview(X)` along the last dimension.
"""
Statistics.var(X::AbstractVectorOfSimilarArrays{T,M}; corrected::Bool = true) where {T,M} =
    var(flatview(X); dims = M + 1, corrected = corrected)

"""
    cov(X::AbstractVectorOfSimilarVectors; corrected::Bool = true)
    cov(X::AbstractVectorOfSimilarVectors, w::StatsBase.AbstractWeights; corrected::Bool = true)

Compute the covariance matrix between the elements of the elements of `X`
along `X`. Equivalent to `cov` of `flatview(X)` along dimension 2.
"""
Statistics.cov(X::AbstractVectorOfSimilarVectors; corrected::Bool = true) =
    cov(flatview(X); dims = 2, corrected = corrected)

"""
    cor(X::AbstractVectorOfSimilarVectors)
    cor(X::AbstractVectorOfSimilarVectors, w::StatsBase.AbstractWeights)

Compute the Pearson correlation matrix between the elements of the elements of
 `X` along `X`. Equivalent to `cor` of `flatview(X)` along dimension 2.
"""
Statistics.cor(X::AbstractVectorOfSimilarVectors) =
    cor(flatview(X); dims = 2)


text:
using NamedDims
@testset "named" begin
    A = rand(Int8, 2,2)
    B = rand(Int8, 3,3)
    ndA = NamedDimsArray(A, (:ia, :ja))
    ndB = NamedDimsArray(B, (:ib, :jb))

    @test NamedDims.names(kronecker(ndA, ndB)) == (Symbol("iaᵡib"), Symbol("jaᵡjb"))
    @test NamedDims.names(kronecker(ndA, B)) == (Symbol("iaᵡ_"), Symbol("jaᵡ_"))
    @test NamedDims.names(kronecker(A, ndB)) == (Symbol("_ᵡib"), Symbol("_ᵡjb"))

    sum(abs2, kronecker(ndA, ndB) .- kron(A, B)) == 0

    @test 0 == @allocated Kronecker._join(:i, :j)
    @test 0 == @allocated Kronecker.kron_names((:ia, :ja), (:ib, :jb))

    @test NamedDims.names(kronecker(ndA, 3)) == (Symbol("iaᵡiaᵡiaᵡia"), Symbol("jaᵡjaᵡjaᵡja"))
    @test 0 == @allocated Kronecker.kron_names((:i, :j), Val(3))
end


text:
function _dft_prop_changed_old_state(e::CustomEvent, prop::Symbol)
    get!(e, prop, nothing)
end

function _dft_prop_changed_new_state(e::CustomEvent, prop::Symbol)
    obj = getindex(extras(e), "OBJ")
    return getproperty(obj, prop)
end

function _dft_prop_changed_update!(e::CustomEvent, prop::Symbol)
    new_mtime = new_state(e, prop)
    setindex!(e, deepcopy(new_mtime), prop)
    return e
end

function _dft_prop_changed_trigger(::CustomEvent, old_val, new_val)
    return !isequal(old_val, new_val)
end

function PropertyChangedEvent(obj)
    e = CustomEvent{Symbol, Any}(
        old_state = _dft_prop_changed_old_state,
        new_state = _dft_prop_changed_new_state,
        update! = _dft_prop_changed_update!,
        trigger = _dft_prop_changed_trigger,
    )
    setindex!(extras(e), obj, "OBJ")
    return e
end


text:
mag_struct = Nemo.mag_struct
acb_calc_integrate_opts = Nemo.acb_calc_integrate_opts

# Code taken from Nemo with minor changes to allow using the analytic
# flag
function acb_calc_func_wrap(res::Ptr{acb}, x::Ptr{acb}, param::Ptr{Nothing}, order::Int, prec::Int)
    xx = unsafe_load(x)
    xx.parent = AcbField(prec)
    F = unsafe_pointer_to_objref(param)
    w = F(xx)
    ccall((:acb_set, Nemo.libarb), Cvoid, (Ptr{acb}, Ref{acb}), res, w)
    return zero(Cint)
end

acb_calc_func_wrap_c() = @cfunction(acb_calc_func_wrap, Cint,
                                    (Ptr{acb}, Ptr{acb}, Ptr{Nothing}, Int, Int))

function acb_calc_func_wrap_analytic(res::Ptr{acb}, x::Ptr{acb}, param::Ptr{Nothing}, order::Int, prec::Int)
    @assert order <= 1

    xx = unsafe_load(x)
    xx.parent = AcbField(prec)

    F = unsafe_pointer_to_objref(param)

    w = F(xx, analytic = order != 0)
    ccall((:acb_set, Nemo.libarb), Cvoid, (Ptr{acb}, Ref{acb}), res, w)
    return zero(Cint)
end

acb_calc_func_wrap_analytic_c() = @cfunction(acb_calc_func_wrap_analytic, Cint,
                                             (Ptr{acb}, Ptr{acb}, Ptr{Nothing}, Int, Int))

const ARB_CALC_SUCCESS = UInt(0)
const ARB_CALC_NO_CONVERGENCE = UInt(2)

function integrate(C::AcbField, F, a, b;
                   rel_tol = -1.0,
                   abs_tol = -1.0,
                   deg_limit::Int = 0,
                   eval_limit::Int = 0,
                   depth_limit::Int = 0,
                   use_heap::Int = 0,
                   verbose::Int = 0,
                   checkanalytic = false)
    opts = acb_calc_integrate_opts(deg_limit, eval_limit, depth_limit,
                                   Cint(use_heap), Cint(verbose))

    lower = C(a)
    upper = C(b)

    cgoal = 0

    if rel_tol === -1.0
        cgoal = precision(C)
    else
        t = BigFloat(rel_tol, RoundDown)
        cgoal_clong = Ref{Clong}()
        ccall((:mpfr_get_d_2exp, :libmpfr), Float64,
              (Ref{Clong}, Ref{BigFloat}, Cint),
              cgoal_clong, t, Base.MPFR.MPFRRoundDown)
        cgoal = -Int(cgoal_clong[]) + 1
    end

    ctol = mag_struct(0, 0)
    ccall((:mag_init, Nemo.libarb), Nothing, (Ref{mag_struct},), ctol)

    if abs_tol === -1.0
        ccall((:mag_set_ui_2exp_si, Nemo.libarb), Nothing, (Ref{mag_struct}, UInt, Int), ctol, 1, -precision(C))
    else
        t = BigFloat(abs_tol, RoundDown)
        expo = Ref{Clong}()
        d = ccall((:mpfr_get_d_2exp, :libmpfr), Float64,
                  (Ref{Clong}, Ref{BigFloat}, Cint),
                  expo, t, Base.MPFR.MPFRRoundDown)
        ccall((:mag_set_d, Nemo.libarb), Nothing, (Ref{mag_struct}, Float64), ctol, d)
        ccall((:mag_mul_2exp_si, Nemo.libarb), Nothing,
              (Ref{mag_struct}, Ref{mag_struct}, Int), ctol, ctol, Int(expo[]))
    end

    res = C()

    status = ccall((:acb_calc_integrate, Nemo.libarb), UInt,
                   (Ref{acb},                       #res
                    Ptr{Nothing},                      #func
                    Any,                            #params
                    Ref{acb},                       #a
                    Ref{acb},                       #b
                    Int,                            #rel_goal
                    Ref{mag_struct},                #abs_tol
                    Ref{acb_calc_integrate_opts},   #opts
                    Int),
                   res,
                   ifelse(checkanalytic, acb_calc_func_wrap_analytic_c(), acb_calc_func_wrap_c()),
                   F, lower, upper, cgoal, ctol, opts, precision(C))

    ccall((:mag_clear, Nemo.libarb), Nothing, (Ref{mag_struct},), ctol)

    if status == ARB_CALC_SUCCESS
        nothing
    elseif status == ARB_CALC_NO_CONVERGENCE
        @warn "Integration did converge"
    end
    return res
end


text:
#==============================================================================
    Code for solving the Hamiltonian Jacboi Bellman for
	   an basic model: a neoclassical growth model: ρV(k) = max_{c} U(c) + V'(k)[F(k)-δk - c]
	   Where s(k) = F(k) - δk - c(k) and c(k) = U'^{-1}(V'(k))

	Translated Julia code from Matlab code by Ben Moll:
        http://www.princeton.edu/~moll/HACTproject.htm

    Updated to Julia 1.0.0
==============================================================================#

using Distributions, Plots, SparseArrays, LinearAlgebra

σ= 2.0 #
ρ = 0.05 #the discount rate
δ = 0.05 # the depreciation rate
A = 1.0
α= 0.3

k_ss = (α*A/(ρ+δ))^(1/(1-α))

H = 10000
k_min = 0.001*k_ss
k_max = 2.0*k_ss

k = LinRange(k_min, k_max, H)
k = convert(Array, k) # create grid for a values
dk = (k_max-k_min)/(H-1)

Δ = 1000

maxit = 10
ε = 10e-6

dVf, dVb = [zeros(H,1) for i in 1:6]

#initial guess for V
v0 = (A.*k.^α).^(1-σ)/(1-σ)/ρ
v= v0

dist=[]

for n=1:maxit
	V=v

    # forward difference
	dVf[1:H-1] = (V[2:H]-V[1:H-1])/dk
	dVf[H] = (A.*k_max.^α - δ.*k_max)^(-σ)

	# backward difference
	dVb[2:H] = (V[2:H]-V[1:H-1])/dk
	dVb[1] = (A.*k_min.^α - δ.*k_min)^(-σ) # the boundary condition

	I_concave = dVb .> dVf

    # consumption and savings with forward difference
    cf = dVf.^(-1/σ)
    μ_f = A.*k.^α - δ.*k -cf

    # consumption and savings with backward difference
    cb = dVb.^(-1/σ)
    μ_b = A.*k.^α - δ.*k -cb

	c0 = A.*k.^α - δ.*k
    dV0 = c0.^(-σ)

    # Now to make a choice between forward and backward difference
    If = μ_f .> 0
    Ib = μ_b .< 0
    I0 = (1.0.-If-Ib)

    global dV_Upwind= dVf.*If + dVb.*Ib + dV0.*I0

    global c = dV_Upwind.^(-1/σ)
    u = c.^(1-σ)/(1-σ)

    # create the transition matrix
    X = -min.(μ_b,0)/dk
    Y = -max.(μ_f,0)/dk + min.(μ_b,0)/dk
    Z = max.(μ_f,0)/dk

    AA = sparse(Diagonal((Y[:]))) + [zeros(1,H); sparse(Diagonal((X[2:H]))) zeros(H-1,1)] + [zeros(H-1,1) sparse(Diagonal((Z[1:H-1]))); zeros(1,H)]
    B = (ρ + 1/Δ)*sparse(I,H,H) - AA
    b = u + V./Δ

    V = B\b
    V_change = V-v
    global v= V

	push!(dist,findmax(abs.(V_change))[1])
	if dist[n] .< ε
		println("Value Function Converged Iteration=")
		println(n)
		break
	end
end

plot(dist, grid=false,
		xlabel="Iteration", ylabel="||V^{n+1} - V^n||",
		ylims=(-0.001,0.030),
		legend=false, title="")
png("Convergence")

v_err = c.^(1-σ)/(1-σ) + dV_Upwind.*(A.*k.^α - δ.*k -c) - ρ.*v

plot(k, v_err, grid=false,
		xlabel="k", ylabel="Error in the HJB equation",
		xlims=(k_min,k_max),
		legend=false, title="")
png("HJB_error")

plot(k, v, grid=false,
		xlabel="k", ylabel="V(k)",
		xlims=(k_min,k_max),
		legend=false, title="")
png("Value_function_vs_k")

plot(k, c, grid=false,
		xlabel="k", ylabel="c(k)",
		xlims=(k_min,k_max),
		legend=false, title="")
png("c(k)_vs_k")

# approximation at the borrowing constraint
k_dot = (A.*k.^α - δ.*k -c)

plot(k, k_dot, grid=false,
		xlabel="k", ylabel="s(k)",
		xlims=(k_min,k_max), title="", label="s(k)", legend=false)
plot!(k, zeros(H,1), label="", line=:dash)
png("stateconstraint")


text:
using LinearAlgebra, ManifoldMeasures, Manifolds, MeasureTheory, Random, StatsBase
using Test
using Manifolds: ℝ

function test_vMF(
    d, μ; check_logdensity=true, estimate_mode=true, max_norm=0.1, ndraws=10_000
)
    M = base_manifold(d)
    x = rand(normalize(basemeasure(d)))

    # check mode is μ
    @test mode(d) ≈ μ

    if check_logdensity
        v = project(M, μ, x)
        @test logdensity(d, μ) ≥ logdensity(d, exp(M, μ, v, 1e-5))
        @test logdensity(d, μ) ≥ logdensity(d, exp(M, μ, v, -1e-5))
    end

    # check rand returns the right type
    @test is_point(M, rand(d))
    @test real(eltype(rand(d))) <: Float64
    @test real(eltype(rand(Float32, d))) <: Float32

    if estimate_mode
        xs = [rand(d) for _ in 1:ndraws]
        @test all(x -> is_point(M, x), xs)
        if d isa VonMises{ℝ}
            z̄ = mean(cis, xs)
            μ̂ = angle(z̄)
        else
            d isa VonMises{ℂ}
            z̄ = mean(xs)
            μ̂ = project(M, z̄)
        end
        @test norm(μ̂ - μ) ≤ max_norm
    end
end

@testset "von Mises-Fisher" begin
    # TODO: use frequentist tests from Mardia & Jupp to check samplers
    @testset "VonMises" begin
        @testset "real (μ, κ)" begin
            @test VonMises(; μ=0.5, κ=1.0) isa VonMisesFisher{Circle{ℝ},(:μ, :κ)}
            @testset for κ in (0.1, 50, 1e5)
                μ = rand() * 2π - π
                test_vMF(VonMises(; μ=μ, κ=κ), μ; estimate_mode=κ == 50)
            end
        end
        @testset "complex (μ, κ)" begin
            @test VonMises(ℂ; μ=im, κ=1.0) isa VonMisesFisher{Circle{ℂ},(:μ, :κ)}
            @testset for κ in (0.1, 50, 1e5)
                μ = sign(randn(ComplexF64))
                test_vMF(VonMises(ℂ; μ=μ, κ=κ), μ; estimate_mode=κ == 50)
            end
        end
        @testset "complex (c,)" begin
            @test VonMises(ℂ; c=3 + 4im) isa VonMisesFisher{Circle{ℂ},(:c,)}
            @testset for κ in (0.1, 50, 1e5)
                c = randn(ComplexF64)
                test_vMF(VonMises(ℂ; c=c), sign(c); estimate_mode=κ == 50)
            end
        end
    end

    @testset "VonMisesFisher" begin
        @testset "show" begin
            c = [1, 2, 3]
            d = VonMisesFisher(Sphere(2, ℂ); c=c)
            @test sprint(show, "text/plain", d) ==
                  "VonMisesFisher(Sphere(2, ℂ); c = [1, 2, 3])"
        end
        @testset "Sphere" begin
            @testset "Sphere($n, $𝔽)" for 𝔽 in (ℝ, ℂ), n in (0, 1, 4)
                T = 𝔽 === ℂ ? ComplexF64 : Float64
                @test VonMisesFisher(n + 1, 𝔽) === VonMisesFisher(Sphere(n, 𝔽))
                @testset "(μ, κ)" begin
                    @testset for κ in (0.1, 50, 1e5)
                        μ = normalize(randn(T, n + 1))
                        test_vMF(
                            VonMisesFisher(n + 1, 𝔽; μ=μ, κ=κ), μ; estimate_mode=κ == 50
                        )
                    end
                end
                @testset "(c,)" begin
                    @testset for κ in (0.1, 50, 1e5)
                        c = randn(T, n + 1)
                        test_vMF(
                            VonMisesFisher(n + 1, 𝔽; c=c),
                            normalize(c);
                            estimate_mode=κ == 50,
                        )
                    end
                end
            end
        end

        @testset "Stiefel" begin
            @testset "Stiefel($n, $k, $𝔽)" for 𝔽 in (ℝ, ℂ), (n, k) in ((3, 1), (4, 3))
                T = 𝔽 === ℂ ? ComplexF64 : Float64
                @test VonMisesFisher(n, k, 𝔽) === VonMisesFisher(Stiefel(n, k, 𝔽))
                @testset "(F,)" begin
                    F = 50 * randn(T, n, k)
                    U, D, V = svd(F)
                    test_vMF(
                        VonMisesFisher(n, k, 𝔽; F=F),
                        U * V';
                        check_logdensity=false,
                        estimate_mode=false,
                    )
                end
                @testset "(U,D,V)" begin
                    U = rand(normalize(Hausdorff(Stiefel(n, k, 𝔽))))
                    V = rand(normalize(Hausdorff(Stiefel(k, k, 𝔽))))
                    D = sort(50 .* abs.(rand(k)))
                    test_vMF(
                        VonMisesFisher(n, k, 𝔽; U=U, D=D, V=V),
                        U * V';
                        check_logdensity=false,
                    )
                end
                @testset "(H,P)" begin
                    H = rand(normalize(Hausdorff(Stiefel(n, k, 𝔽))))
                    V = rand(normalize(Hausdorff(Stiefel(k, k, 𝔽))))
                    D = sort(50 .* abs.(randn(k)))
                    P = V * Diagonal(D) * V'
                    test_vMF(
                        VonMisesFisher(n, k, 𝔽; H=H, P=P),
                        H;
                        check_logdensity=false,
                        estimate_mode=false,
                    )
                end
            end
        end
    end
end


text:
## Memory handling for memory

#
struct CodeMem
    slots::Vector{Symbol}
    slot_names::Vector{String}
    special_names::Dict{Symbol,String}
end

# Creates a memory object
# k is the number of slots
# slot_name_fun(i) should return a string
#    which is the name of the node
function CodeMem(k, slot_name_fun)
    slots = map(i -> :Free, 1:k) # Special keyword
    #slot_names=Vector{String}(undef,k);
    slot_names = slot_name_fun.(1:k)
    return CodeMem(slots, slot_names, Dict{Symbol,String}())
end

function get_slot_number(mem::CodeMem, node::Symbol)
    return findfirst(mem.slots .== node)
end

function set_slot_number!(mem, slot_number, node)
    return mem.slots[slot_number] = node
end

function get_slot_name(mem::CodeMem, node::Symbol)
    if (haskey(mem.special_names, node))
        return mem.special_names[node]
    else
        if (node in mem.slots)
            return mem.slot_names[findfirst(mem.slots .== node)]
        else
            @show mem.slots
            error("Unable to find memory slot for $node")
        end
    end
end

function is_allocated(mem::CodeMem, node::Symbol)
    return node in mem.slots
end

# Returns a reference to a memory slot
# which is either a new slot or
# one of the elements in the Vector allow_inplace
#
# Returns a Tuple with a reference to the slot
# and the name
#
# This does not actually allocate
function get_free_slot(mem::CodeMem, allow_inplace = [])
    for i = 1:size(mem.slots, 1)
        if mem.slots[i] == :Free || mem.slots[i] in allow_inplace
            return (i, mem.slot_names[i])
        end
    end
    return error("No free memory slot found!")
end

# Allocates a slot that has been found using
# get_free_slot
function alloc_slot!(mem::CodeMem, i, node::Symbol)
    return mem.slots[i] = node
end

function free!(mem::CodeMem, i)
    return mem.slots[i] = :Free
end


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

JLLWrappers.@generate_main_file_header("TensorRT")
JLLWrappers.@generate_main_file("TensorRT", UUID("2eaff018-c6f9-5ddb-97b6-8b0d18888fbb"))
end  # module TensorRT_jll


text:
module BEASTTreeUtils

export rtree,
       standardize_height!,
       writeTopology,
       vcv

# Internal library for working with trees. functions rely heavily on PhyloNetworks.jl

using PhyloNetworks, Distributions, LinearAlgebra
pn = PhyloNetworks

# for backward compatibility with earlier code
function parse_newick(newick::String)
    return readTopology(newick)
end

@deprecate parse_newick(newick) PhyloNetworks.readTopology(newick)

function rnewick(taxa::Vector{String}; args...)
    return writeTopology(rtree(taxa; args...))
end

### Generate HybridTree with random topologies / branch lenghts

# recursively add nodes to HybridNetwork
function random_tips!(net::HybridNetwork,
                        parent::pn.Node,
                        rng::UnitRange{Int}
                        )

    split_ind = rand(first(rng):(last(rng) - 1))

    l_rng = first(rng):split_ind
    r_rng = (split_ind + 1):last(rng)

    num_internal = parent.number

    for new_rng in (l_rng, r_rng)
        if length(new_rng) > 1
            new_parent = pn.Node(num_internal - 1, false)

            pn.parseTreeNode!(new_parent, parent, net)

            random_tips!(net, new_parent, new_rng)

            num_internal -= 2 * length(new_rng) - 1
        else
            leaf_node = pn.Node(new_rng[1], true)
            pn.parseTreeNode!(leaf_node, parent, net)
        end
    end

end

function rtree(n::Int;
                keep_order::Bool = false, # keep (true) or scramble (false) initial label ordering
                ultrametric::Bool = false, # all leaf nodes at same time (true) or not (false)
                edge_dist::ContinuousUnivariateDistribution = Exponential(1.0)
                )
    labels = ["t$i" for i = 1:n]
    return rtree(labels, keep_order = keep_order, ultrametric = ultrametric,
                    edge_dist = edge_dist)
end

# generates a random HybridNetwork
function rtree(labels::Array{<:AbstractString}; # leaf labels
                keep_order::Bool = false, # keep (true) or scramble (false) initial label ordering
                ultrametric::Bool = false, # all leaf nodes at same time (true) or not (false)
                edge_dist::ContinuousUnivariateDistribution = Exponential(1.0)) # distribution for generating random edge lengths
                                                                                # (only applies with ultrametric = false)

    n = length(labels)

    if !ultrametric

        net = HybridNetwork()

        root = pn.Node(-2, false)

        random_tips!(net, root, 1:n)

        pn.pushNode!(net, root)
        net.root = 2 * n - 1

        for edge in net.edge
            edge.length = rand(edge_dist)
        end

    else

        # code largely copied from nj.jl in PhyloNetworks.jl

        available_nodes = [pn.Node(i, true) for i = 1:n]

        node_heights = zeros(n)
        t = 0.0

        net = HybridNetwork(copy(available_nodes), pn.Edge[])
        for i = 1:(n - 1)
            m = n + 1 - i
            i1, i2 = sample(1:m, 2, replace = false)

            n1 = available_nodes[i1]
            n2 = available_nodes[i2]
            t += rand(Exponential(1/m)) #new maximum height

            t1 = t - node_heights[i1]
            t2 = t - node_heights[i2]

            e1 = pn.Edge(net.numEdges + 1, t1)
            e2 = pn.Edge(net.numEdges + 2, t2)

            new_node = pn.Node(net.numNodes + 1, false, false, [e1, e2])

            pn.setNode!(e1, pn.Node[n1, new_node])
            pn.setNode!(e2, pn.Node[n2, new_node])

            pn.setEdge!(n1, e1)
            pn.setEdge!(n2, e2)

            pn.pushEdge!(net, e1)
            pn.pushEdge!(net, e2)

            pn.pushNode!(net, new_node)

            available_nodes[i1] = new_node
            node_heights[i1] = t

            deleteat!(available_nodes, i2) #TODO: make more memory efficient with @view
            deleteat!(node_heights, i2)
        end

        net.root = net.numNodes

    end

    directEdges!(net) # not sure this is necessary
    preorder!(net) # not sure this is necessary

    tip_order = collect(1:n)

    if !keep_order
        tip_order = sortperm(rand(n))
    end

    for i = 1:n
        nm = labels[tip_order[i]]
        net.leaf[i].name = nm
        push!(net.names, nm)
    end

    return readTopology(writeTopology(net)) # hacky, but should catch some problems
end

function traverse_node_distances!(dists::Vector{Float64},
                                net::HybridNetwork,
                                node::PhyloNetworks.Node,
                                next_leaf::Int
                                )

    start_ind = next_leaf
    for edge in node.edge
        child = edge.isChild1 ? edge.node[1] : edge.node[2]
        if child !== node
            if child.leaf
                # @show net.leaf[next_leaf].name
                # @show child.name
                @assert net.leaf[next_leaf] === child
                dists[next_leaf] = edge.length
                next_leaf += 1
            else
                rng = traverse_node_distances!(dists, net, child, next_leaf)
                for i in rng
                    dists[i] += edge.length
                end

                next_leaf = last(rng) + 1
            end
        end
    end
    return start_ind:(next_leaf - 1)
end

function leaf_distances(net::HybridNetwork)
    root = net.node[net.root]
    dists = zeros(length(net.leaf))

    traverse_node_distances!(dists, net, root, 1)
    return dists
end

import PhyloNetworks.vcv

function vcv(net::HybridNetwork, taxa::Vector{String})
    v = PhyloNetworks.vcv(net)
    sim_taxa = string.(names(v))

    perm = indexin(taxa, sim_taxa)
    V = Matrix(v)[perm, perm]
end

function vcv(nwk::String, taxa::Vector{String})
    return vcv(readTopology(nwk), taxa)
end

function scale_to!(net::HybridNetwork, height::Float64)
    preorder!(net)
    original_height = maximum(getNodeAges(net))
    mult = height / original_height
    for edge in net.edge
        edge.length *= mult
    end
end

function standardize_height!(net::HybridNetwork)
    scale_to!(net, 1.0)
end

end

max_stars_repo_name: Julia-BEAST-utils/BEASTTreeUtils

max_issues_repo_name: Julia-BEAST-utils/BEASTTreeUtils

max_forks_repo_name: Julia-BEAST-utils/BEASTTreeUtils


text:
function compute_zi(snap::SNAPParams, runtime_arrays::RuntimeArrays)
    idouble = 0;
    num_elements = snap.prebuilt_arrays.num_elements
    # n_elements = length(snap.elements)
    for elem1 = 0:(num_elements-1)
        for elem2 = 0:(num_elements-1)

            for jjz = 1:snap.prebuilt_arrays.idxz_max
                j1 = snap.prebuilt_arrays.idxz[jjz].j1;
                j2 = snap.prebuilt_arrays.idxz[jjz].j2;
                j = snap.prebuilt_arrays.idxz[jjz].j;
                ma1min = snap.prebuilt_arrays.idxz[jjz].ma1min;
                ma2max = snap.prebuilt_arrays.idxz[jjz].ma2max;
                na = snap.prebuilt_arrays.idxz[jjz].na;
                mb1min = snap.prebuilt_arrays.idxz[jjz].mb1min;
                mb2max = snap.prebuilt_arrays.idxz[jjz].mb2max;
                nb = snap.prebuilt_arrays.idxz[jjz].nb;

                cgblock = snap.prebuilt_arrays.cglist[snap.prebuilt_arrays.idxcg_block[j1+1, j2+1, j+1]+1:end];
                # println("cgblock $cgblock")
                runtime_arrays.zlist_r[idouble*snap.prebuilt_arrays.idxz_max + jjz] = 0.0;
                runtime_arrays.zlist_i[idouble*snap.prebuilt_arrays.idxz_max + jjz] = 0.0;

                jju1 = snap.prebuilt_arrays.idxu_block[j1+1] + (j1 + 1) * mb1min;
                jju2 = snap.prebuilt_arrays.idxu_block[j2+1] + (j2 + 1) * mb2max;
                icgb = mb1min * (j2 + 1) + mb2max;
                for ib = 0:(nb-1)

                    suma1_r = 0.0;
                    suma1_i = 0.0;
                    u1_r = runtime_arrays.ulisttot_r[(elem1*snap.prebuilt_arrays.idxu_max+jju1+1):end];
                    u1_i = runtime_arrays.ulisttot_i[(elem1*snap.prebuilt_arrays.idxu_max+jju1+1):end];
                    u2_r = runtime_arrays.ulisttot_r[(elem2*snap.prebuilt_arrays.idxu_max+jju2+1):end];
                    u2_i = runtime_arrays.ulisttot_i[(elem2*snap.prebuilt_arrays.idxu_max+jju2+1):end];

                    ma1 = ma1min;
                    ma2 = ma2max;
                    icga = ma1min * (j2 + 1) + ma2max;

                    for ia = 0:(na-1)
                        suma1_r += cgblock[icga+1] * (u1_r[ma1+1] * u2_r[ma2+1] - u1_i[ma1+1] * u2_i[ma2+1]);
                        suma1_i += cgblock[icga+1] * (u1_r[ma1+1] * u2_i[ma2+1] + u1_i[ma1+1] * u2_r[ma2+1]);
                        ma1+=1;
                        ma2-=1;
                        icga += j2;
                    end # loop over ia
                    # println("suma1_r $suma1_r, suma1_i $suma1_i, cgblock[icgb] $(cgblock[icgb+1])")
                    runtime_arrays.zlist_r[jjz] += cgblock[icgb+1] * suma1_r;
                    runtime_arrays.zlist_i[jjz] += cgblock[icgb+1] * suma1_i;

                    jju1 += j1 + 1;
                    jju2 -= j2 + 1;
                    icgb += j2;
                end # loop over ib
                if (snap.bnorm_flag)
                    runtime_arrays.zlist_r[jjz] /= (j+1);
                    runtime_arrays.zlist_i[jjz] /= (j+1);
                end

            end # loop over jjz
            idouble+=1;
        end
    end
end