text:
# file containing methods for requesting forecast simulations from GEOGloWS
"""
    forecast_stats(parameters::AbstractDict{Symbol,<:Any})

returns statistics calculated from 51 forecast ensemble members.
"""
function forecast_stats(parameters::AbstractDict{Symbol,<:Any})
    endpoint = "ForecastStats"

    r = _make_request(baseurl, endpoint, parameters)

    _read_response(r, parameters[:return_format], with_dates = true)
end

"""
    forecast_stats(reach_id::Int; return_format::Symbol = :csv)

returns statistics calculated from 51 forecast ensemble members.
"""
function forecast_stats(reach_id::Int; return_format::Symbol = :csv)
    parameters = Dict(:reach_id => reach_id, :return_format => return_format)

    forecast_stats(parameters)
end

"""
    forecast_stats(reach_id::Int,date::AbstractString; return_format::Symbol = :csv)

returns statistics calculated from 51 forecast ensemble members.
"""
function forecast_stats(
    reach_id::Int,
    date::AbstractString;
    return_format::Symbol = :csv,
)
    parameters = Dict(
        :reach_id => reach_id,
        :date => date,
        :return_format => return_format,
    )

    forecast_stats(parameters)
end

"""
    forecast_stats(reach_id::Int,date::TimeType; return_format::Symbol = :csv)

returns statistics calculated from 51 forecast ensemble members.
"""
function forecast_stats(
    reach_id::Int,
    date::TimeType;
    return_format::Symbol = :csv,
)
    datestr = try
        Dates.format(date, "YYYY-m-d")
    catch
        Dates.format(date, "yyyymmdd")
    end

    parameters = Dict(
        :reach_id => reach_id,
        :date => datestr,
        :return_format => return_format,
    )

    forecast_stats(parameters)
end

"""
    forecast_stats(lat::Real, lon::Real; return_format::Symbol = :csv)

returns statistics calculated from 51 forecast ensemble members.
"""
function forecast_stats(lat::Real, lon::Real; return_format::Symbol = :csv)
    parameters = Dict(:lat => lat, :lon => lon, :return_format => return_format)

    forecast_stats(parameters)
end

"""
    forecast_stats(lat::Real, lon::Real, date::AbstractString; return_format::Symbol = :csv)

returns statistics calculated from 51 forecast ensemble members.
"""
function forecast_stats(
    lat::Real,
    lon::Real,
    date::AbstractString;
    return_format::Symbol = :csv,
)
    if occursin("-", date)
        date = replace(date, "-" => "")
    end

    parameters = Dict(
        :lat => lat,
        :lon => lon,
        :date => date,
        :return_format => return_format,
    )

    forecast_stats(parameters)
end

"""
    forecast_stats(lat::Real, lon::Real, date::TimeType; return_format::Symbol = :csv)

returns statistics calculated from 51 forecast ensemble members.
"""
function forecast_stats(
    lat::Real,
    lon::Real,
    date::TimeType;
    return_format::Symbol = :csv,
)
    datestr = Dates.format(date, "yyyymmdd")

    parameters = Dict(
        :lat => lat,
        :lon => lon,
        :date => datestr,
        :return_format => return_format,
    )

    forecast_stats(parameters)
end

"""
    forecast_ensembles(parameters::AbstractDict{Symbol,<:Any})

returns a timeseries for each of the 51 normal forecast ensemble members and the 52nd higher resolution forecast
"""
function forecast_ensembles(parameters::AbstractDict{Symbol,<:Any})
    endpoint = "ForecastEnsembles"

    r = _make_request(baseurl, endpoint, parameters)

    _read_response(r, parameters[:return_format], with_dates = true)
end

"""
    forecast_ensembles(reach_id::Int; return_format::Symbol = :csv)

returns a timeseries for each of the 51 normal forecast ensemble members and the 52nd higher resolution forecast
"""
function forecast_ensembles(reach_id::Int; return_format::Symbol = :csv)
    parameters = Dict(:reach_id => reach_id, :return_format => return_format)

    forecast_ensembles(parameters)
end

"""
    forecast_ensembles(reach_id::Int, date::AbstractString; return_format::Symbol = :csv)

returns a timeseries for each of the 51 normal forecast ensemble members and the 52nd higher resolution forecast
"""
function forecast_ensembles(
    reach_id::Int,
    date::AbstractString;
    return_format::Symbol = :csv,
)
    if occursin("-", date)
        date = replace(date, "-" => "")
    end

    parameters = Dict(
        :reach_id => reach_id,
        :date => date,
        :return_format => return_format,
    )

    forecast_ensembles(parameters)
end

"""
    forecast_ensembles(reach_id::Int, date::TimeType; return_format::Symbol = :csv)

returns a timeseries for each of the 51 normal forecast ensemble members and the 52nd higher resolution forecast
"""
function forecast_ensembles(
    reach_id::Int,
    date::TimeType;
    return_format::Symbol = :csv,
)
    datestr = Dates.format(date, "yyyymmdd")

    parameters = Dict(
        :reach_id => reach_id,
        :date => datestr,
        :return_format => return_format,
    )

    forecast_ensembles(parameters)
end

"""
    forecast_ensembles(lat::Real, lon::Real; return_format::Symbol = :csv)

returns a timeseries for each of the 51 normal forecast ensemble members and the 52nd higher resolution forecast
"""
function forecast_ensembles(lat::Real, lon::Real; return_format::Symbol = :csv)
    parameters = Dict(:lat => lat, :lon => lon, :return_format => return_format)

    forecast_ensembles(parameters)
end

"""
    forecast_ensembles(lat::Real, lon::Real, date::AbstractString; return_format::Symbol = :csv)

returns a timeseries for each of the 51 normal forecast ensemble members and the 52nd higher resolution forecast
"""
function forecast_ensembles(
    lat::Real,
    lon::Real,
    date::AbstractString;
    return_format::Symbol = :csv,
)
    if occursin("-", date)
        date = replace(date, "-" => "")
    end

    parameters = Dict(
        :lat => lat,
        :lon => lon,
        :date => date,
        :return_format => return_format,
    )

    forecast_ensembles(parameters)
end

"""
    forecast_ensembles(lat::Real, lon::Real, date::TimeType; return_format::Symbol = :csv)

returns a timeseries for each of the 51 normal forecast ensemble members and the 52nd higher resolution forecast
"""
function forecast_ensembles(
    lat::Real,
    lon::Real,
    date::TimeType;
    return_format::Symbol = :csv,
)
    datestr = Dates.format(date, "yyyymmdd")

    parameters = Dict(
        :lat => lat,
        :lon => lon,
        :date => datestr,
        :return_format => return_format,
    )

    forecast_ensembles(parameters)
end

"""
    forecast_records(parameters::AbstractDict{Symbol,<:Any})

retrieves the rolling record of the mean of the forecasted streamflow during the
first 24 hours of each day's forecast. That is, each day day after the streamflow
forecasts are computed, the average of first 8 of the 3-hour timesteps are recorded to a csv.
"""
function forecast_records(parameters::AbstractDict{Symbol,<:Any})
    endpoint = "ForecastRecords"

    r = _make_request(baseurl, endpoint, parameters)

    _read_response(r, parameters[:return_format], with_dates = true)
end

"""
    forecast_records(reach_id::Int; return_format::Symbol = :csv)

retrieves the rolling record of the mean of the forecasted streamflow during the
first 24 hours of each day's forecast. That is, each day day after the streamflow
forecasts are computed, the average of first 8 of the 3-hour timesteps are recorded to a csv.
"""
function forecast_records(reach_id::Int; return_format::Symbol = :csv)
    parameters = Dict(:reach_id => reach_id, :return_format => return_format)

    forecast_records(parameters)
end

"""
    forecast_records(lat::Real, lon::Real; return_format::Symbol = :csv)

retrieves the rolling record of the mean of the forecasted streamflow during the
first 24 hours of each day's forecast. That is, each day day after the streamflow
forecasts are computed, the average of first 8 of the 3-hour timesteps are recorded to a csv.
"""
function forecast_records(lat::Real, lon::Real; return_format::Symbol = :csv)
    parameters = Dict(:lat => lat, :lon => lon, :return_format => return_format)

    forecast_records(parameters)
end

"""
    forecast_records(reach_id::Int, start_date::TimeType, end_date::TimeType; return_format::Symbol = :csv,)

retrieves the rolling record of the mean of the forecasted streamflow during the
first 24 hours of each day's forecast. That is, each day day after the streamflow
forecasts are computed, the average of first 8 of the 3-hour timesteps are recorded to a csv.
"""
function forecast_records(
    reach_id::Int,
    start_date::TimeType,
    end_date::TimeType;
    return_format::Symbol = :csv,
)
    start_datestr = Dates.format(start_date, "yyyymmdd")
    end_datestr = Dates.format(end_date, "yyyymmdd")

    parameters = Dict(
        :reach_id => reach_id,
        :start_date => start_datestr,
        :end_date => end_datestr,
        :return_format => return_format,
    )

    forecast_records(parameters)
end

"""
    forecast_records(reach_id::Int, start_date::AbstractString, end_date::AbstractString; return_format::Symbol = :csv,)

retrieves the rolling record of the mean of the forecasted streamflow during the
first 24 hours of each day's forecast. That is, each day day after the streamflow
forecasts are computed, the average of first 8 of the 3-hour timesteps are recorded to a csv.
"""
function forecast_records(
    reach_id::Int,
    start_date::AbstractString,
    end_date::AbstractString;
    return_format::Symbol = :csv,
)
    if occursin("-", start_date)
        start_date = replace(start_date, "-" => "")
    end
    if occursin("-", end_date)
        end_date = replace(end_date, "-" => "")
    end

    parameters = Dict(
        :reach_id => reach_id,
        :start_date => start_date,
        :end_date => end_date,
        :return_format => return_format,
    )

    forecast_records(parameters)
end

"""
    forecast_records(lat::Real, lon::Real, start_date::TimeType, end_date::TimeType; return_format::Symbol = :csv,)

retrieves the rolling record of the mean of the forecasted streamflow during the
first 24 hours of each day's forecast. That is, each day day after the streamflow
forecasts are computed, the average of first 8 of the 3-hour timesteps are recorded to a csv.
"""
function forecast_records(
    lat::Real,
    lon::Real,
    start_date::TimeType,
    end_date::TimeType;
    return_format::Symbol = :csv,
)
    start_datestr = Dates.format(start_date, "yyyymmdd")
    end_datestr = Dates.format(end_date, "yyyymmdd")

    parameters = Dict(
        :lat => lat,
        :lon => lon,
        :start_date => start_datestr,
        :end_date => end_datestr,
        :return_format => return_format,
    )

    forecast_records(parameters)
end

"""
    forecast_records(lat::Real, lon::Real, start_date::AbstractString, end_date::AbstractString; return_format::Symbol = :csv,)

retrieves the rolling record of the mean of the forecasted streamflow during the
first 24 hours of each day's forecast. That is, each day day after the streamflow
forecasts are computed, the average of first 8 of the 3-hour timesteps are recorded to a csv.
"""
function forecast_records(
    lat::Real,
    lon::Real,
    start_date::AbstractString,
    end_date::AbstractString;
    return_format::Symbol = :csv,
)
    if occursin("-", start_date)
        start_date = replace(start_date, "-" => "")
    end
    if occursin("-", end_date)
        end_date = replace(end_date, "-" => "")
    end

    parameters = Dict(
        :lat => lat,
        :lon => lon,
        :start_date => start_date,
        :end_date => end_date,
        :return_format => return_format,
    )

    forecast_records(parameters)
end

"""
    forecast_warnings(parameters::AbstractDict{Symbol,<:Any})

returns a csv created to summarize the potential return period level flow events
coming to the reaches in a specified region. The CSV contains a column for the
reach_id, lat of the reach, lon of the reach, maximum forecasted flow in the next
15 day forecast, and a column for each of the return periods (2, 10, 20, 25, 50, 100)
which will contain the date when the forecast is first expected to pass that threshold
"""
function forecast_warnings(parameters::AbstractDict{Symbol,<:Any})
    endpoint = "ForecastWarnings"

    r = _make_request(baseurl, endpoint, parameters)

    _read_response(r, parameters[:return_format])
end

"""
    forecast_warnings(; return_format::Symbol = :csv)

returns a csv created to summarize the potential return period level flow events
coming to the reaches in a specified region. The CSV contains a column for the
reach_id, lat of the reach, lon of the reach, maximum forecasted flow in the next
15 day forecast, and a column for each of the return periods (2, 10, 20, 25, 50, 100)
which will contain the date when the forecast is first expected to pass that threshold
"""
function forecast_warnings(; return_format::Symbol = :csv)
    parameters = Dict(:return_format => return_format)

    forecast_warnings(parameters)
end

"""
    forecast_warnings(date::TimeType; return_format::Symbol = :csv)

returns a csv created to summarize the potential return period level flow events
coming to the reaches in a specified region. The CSV contains a column for the
reach_id, lat of the reach, lon of the reach, maximum forecasted flow in the next
15 day forecast, and a column for each of the return periods (2, 10, 20, 25, 50, 100)
which will contain the date when the forecast is first expected to pass that threshold
"""
function forecast_warnings(date::TimeType; return_format::Symbol = :csv)
    datestr = Dates.format(date, "yyyymmdd")

    parameters = Dict(:date => datestr, :return_format => return_format)

    forecast_warnings(parameters)
end

"""
    forecast_warnings(region::AbstractString; return_format::Symbol = :csv)

returns a csv created to summarize the potential return period level flow events
coming to the reaches in a specified region. The CSV contains a column for the
reach_id, lat of the reach, lon of the reach, maximum forecasted flow in the next
15 day forecast, and a column for each of the return periods (2, 10, 20, 25, 50, 100)
which will contain the date when the forecast is first expected to pass that threshold
"""
function forecast_warnings(region::AbstractString; return_format::Symbol = :csv)
    parameters = Dict(:region => region, :return_format => return_format)

    forecast_warnings(parameters)
end

"""
    forecast_warnings(region::AbstractString, date::TimeType; return_format::Symbol = :csv)

returns a csv created to summarize the potential return period level flow events
coming to the reaches in a specified region. The CSV contains a column for the
reach_id, lat of the reach, lon of the reach, maximum forecasted flow in the next
15 day forecast, and a column for each of the return periods (2, 10, 20, 25, 50, 100)
which will contain the date when the forecast is first expected to pass that threshold
"""
function forecast_warnings(
    region::AbstractString,
    date::TimeType;
    return_format::Symbol = :csv,
)
    datestr = Dates.format(date, "yyyymmdd")

    parameters = Dict(
        :region => region,
        :date => datestr,
        :return_format => return_format,
    )

    forecast_warnings(parameters)
end

"""
    forecast_warnings(region::AbstractString, date::AbstractString; return_format::Symbol = :csv)

returns a csv created to summarize the potential return period level flow events
coming to the reaches in a specified region. The CSV contains a column for the
reach_id, lat of the reach, lon of the reach, maximum forecasted flow in the next
15 day forecast, and a column for each of the return periods (2, 10, 20, 25, 50, 100)
which will contain the date when the forecast is first expected to pass that threshold
"""
function forecast_warnings(
    region::AbstractString,
    date::AbstractString;
    return_format::Symbol = :csv,
)
    if occursin("-", date)
        date = replace(date, "-" => "")
    end

    parameters = Dict(
        region::AbstractString,
        :date => date,
        :return_format => return_format,
    )

    forecast_warnings(parameters)
end


text:
# # Spatial transformations
#
# DataAugmentation.jl comes with various spatial transformations that you
# can apply to your data. You can apply them to [`Image`](#)s and
# the keypoint-based items [`Keypoints`](#), [`Polygon`](#), and [`BoundingBox`](#).
#
# Let's take this picture of a light house:
# {cell=main style="display:none;" result=false output=false}
using DataAugmentation
using MosaicViews
using Images
using TestImages
using StaticArrays

imagedata = testimage("lighthouse")
imagedata = imresize(imagedata, ratio = 196 / size(imagedata, 1))
#
imagedata = testimage("lighthouse")
# {cell=main style="display:none;"}
imagedata
# To apply a transformation `tfm` to it, wrap it in
# `Image`, apply the transformation and unwrap it using [`itemdata`](#):
# {cell=main }

tfm = Crop((196, 196))
image = Image(imagedata)
apply(tfm, image) |> itemdata

# Now let's say we want to train a light house detector and have a bounding box
# for the light house. We can use the [`BoundingBox`](#) item to represent it.
# It takes the two corners of the bounding rectangle as the first argument. As
# the second argument we have to pass the size of the corresponding image.
#
# {cell=main}
points = SVector{2, Float32}[SVector(23., 120.), SVector(120., 150.)]
bbox = BoundingBox(points, size(imagedata))

# [`showitem`](#) visualizes the two items:
# {cell=main}

showitem((image, bbox))

# If we apply transformations like translation and cropping
# to the image, then the same transformations have to be applied to the bounding
# box. Otherwise, the bounding box will no longer match up with the light house.
#
# Another problem can occur with stochastic transformations like [`RandomResizeCrop`](#).
# If we apply it separately to the image and the bounding box, they will be cropped from
# slightly different locations:
#
# {cell=main}

tfm = RandomResizeCrop((128, 128))
showitem((
    apply(tfm, image),
    apply(tfm, bbox)
))

# Instead, pass a tuple of the items to a single `apply` call so the same
# random state will be used for both image and bounding box:
#
# {cell=main}

apply(tfm, (image, bbox)) |> showitem

#
# ## Gallery
# {cell=main style="display:none;" result=false}

function showtransform(tfm, item, n = 8; ncol = 4)
    return mosaicview(
        [showitem(apply(tfm, item)) for _ in 1:n],
        fillvalue = RGBA(1, 1, 1, 0),
        npad = 8,
        rowmajor = true,
        ncol = ncol)
end

function showtransforms(tfms, item; ncol = length(tfms))
    return mosaicview(
        [showitem(apply(tfm, item)) for tfm in tfms],
        fillvalue = RGBA(1, 1, 1, 0),
        npad = 8,
        rowmajor = true,
        ncol = ncol)
end

# ### [`RandomResizeCrop`](#)`(sz)`
#
# Resizes the sides so that one of them is no longer than `sz` and
# crops a region of size `sz` *from a random location*.
#
# {cell=main result=false}

tfm = RandomResizeCrop((128, 128))

# {cell=main style="display:none;"}

o = showtransform(tfm, (image, bbox), 6, ncol=6)

# ### [`CenterResizeCrop`](#)
#
# Resizes the sides so that one of them is no longer than `sz` and
# crops a region of size `sz` *from the center*.
#
# {cell=main result=false}

tfm = CenterResizeCrop((128, 128))

# {cell=main style="display:none;"}

o = showtransform(tfm, (image, bbox), 1)

# ### [`Crop`](#)`(sz[, from])`
#
# Crops a region of size `sz` from the image, *without resizing* the image first.
#
# {cell=main result=false}

using DataAugmentation: FromOrigin, FromCenter, FromRandom
tfms = [
    Crop((128, 128), FromOrigin()),
    Crop((128, 128), FromCenter()),
    Crop((128, 128), FromRandom()),
    Crop((128, 128), FromRandom()),
    Crop((128, 128), FromRandom()),
    Crop((128, 128), FromRandom()),
]

# {cell=main style="display:none;"}

o = showtransforms(tfms, (image, bbox))

ext: jl


text:
"""
This package provides provides a set of functions to facilitate the field-theoretic
treatment of monolayer graphene using the tight-binding model. The Hamiltonian
employed by this package includes only the nearest-neighbor hopping term with
``t = 2.8``eV. The derivation of the formalism is available
[here](https://arxiv.org/pdf/2007.06984.pdf).
"""
module GrapheneQFT

using Kronecker
using LinearAlgebra
using MLStyle
using QuadGK

export A,
    B,
    SpinUp,
    SpinDown,
    GrapheneCoord,
    GrapheneState,
    crystal_to_cartesian,
    graphene_neighbors,
    graphene_multiple_neighbors,
    Defect,
    ImpurityState,
    LocalSpin,
    Hopping,
    GrapheneSystem,
    mkGrapheneSystem,
    δG_R,
    G_R,
    δΓ,
    Γ,
    δρ_R_graphene,
    peierls_phase

include("computed_quantities.jl")

end # module


text:
#
# LazyAlgebra.jl -
#
# A simple linear algebra system.
#
#-------------------------------------------------------------------------------
#
# This file is part of LazyAlgebra (https://github.com/emmt/LazyAlgebra.jl)
# released under the MIT "Expat" license.
#
# Copyright (c) 2017-2021 Éric Thiébaut.
#

module LazyAlgebra

export
    CirculantConvolution,
    CompressedSparseOperator,
    CroppingOperator,
    Diag,
    Diff,
    FFTOperator,
    GeneralMatrix,
    Gram,
    Id,
    Identity,
    Jacobian,
    LinearMapping,
    Mapping,
    NonuniformScaling,
    RankOneOperator,
    SingularSystem,
    SparseOperator,
    SparseOperatorCOO,
    SparseOperatorCSC,
    SparseOperatorCSR,
    SymbolicLinearMapping,
    SymbolicMapping,
    SymmetricRankOneOperator,
    ZeroPaddingOperator,
    ∇,
    adjoint,
    apply!,
    apply,
    coefficients,
    col_size,
    conjgrad!,
    conjgrad,
    diag,
    gram,
    input_eltype,
    input_ndims,
    input_size,
    input_type,
    is_diagonal,
    is_endomorphism,
    is_linear,
    is_selfadjoint,
    isone,
    iszero,
    jacobian,
    lgemm!,
    lgemm,
    lgemv!,
    lgemv,
    multiplier,
    ncols,
    nnz,
    nonzeros,
    nrows,
    output_eltype,
    output_ndims,
    output_size,
    output_type,
    row_size,
    sparse,
    terms,
    unpack!,
    unscaled,
    unveil,
    vcombine!,
    vcombine,
    vcopy!,
    vcopy,
    vcreate,
    vdot,
    vfill!,
    vmul!,
    vmul,
    vnorm1,
    vnorm2,
    vnorminf,
    vones,
    vproduct!,
    vproduct,
    vscale!,
    vscale,
    vswap!,
    vupdate!,
    vzero!,
    vzeros

using Printf
using ArrayTools

import Base: *, ∘, +, -, \, /, ==
import Base: Tuple, adjoint, inv, axes,
    show, showerror, convert, eltype, ndims, size, length, stride, strides,
    getindex, setindex!, eachindex, first, last, firstindex, lastindex,
    one, zero, isone, iszero, @propagate_inbounds

# Import/using from LinearAlgebra, BLAS and SparseArrays.
using LinearAlgebra
import LinearAlgebra: UniformScaling, diag, ⋅, mul!, rmul!
using LinearAlgebra.BLAS
using LinearAlgebra.BLAS: libblas, @blasfunc,
    BlasInt, BlasReal, BlasFloat, BlasComplex

using SparseArrays: sparse

include("types.jl")
include("traits.jl")
include("utils.jl")
include("methods.jl")
include("vectors.jl")
include("genmult.jl")
import .GenMult: lgemm!, lgemm, lgemv!, lgemv
include("blas.jl")
include("rules.jl")
include("mappings.jl")
include("foundations.jl")

include("sparse.jl")
using .SparseOperators
import .SparseOperators: unpack!

include("cropping.jl")
import .Cropping: CroppingOperator, ZeroPaddingOperator, defaultoffset
include("diff.jl")
import .FiniteDifferences: Diff
include("fft.jl")
import .FFTs: CirculantConvolution, FFTOperator
include("conjgrad.jl")
include("init.jl")

end


text:
function pricer(mcProcess::FinancialMonteCarlo.BaseProcess, StrikeVec::Array{U, 1}, zero_rate::FinancialMonteCarlo.AbstractZeroRateCurve, T::Number, method::CarrMadanMethod, ::FinancialMonteCarlo.CudaMode) where {U <: Number}
    Npow = method.Npow
    N = 2^Npow
    S0 = mcProcess.underlying.S0
    r = FinancialMonteCarlo.integral(zero_rate.r, T) / T
    d = FinancialMonteCarlo.integral(FinancialMonteCarlo.dividend(mcProcess), T) / T
    A = method.A
    #* v-> compute integral as a summation
    dx = A / N
    real_vec = 0:(N-1)
    v = collect(real_vec) * dx            # the final value A is excluded
    v[1] = 1e-312
    cf = FinancialFFT.CharactheristicFunction(mcProcess, T)
    corr = FinancialFFT.CharactheristicExponent(-1im, mcProcess, T)
    CharFunc(v) = cf(v) * exp(-v * 1im * corr)
    integrand_f(v) = exp(1im * (r - d) * v * T) * (CharFunc(v - 1im) - 1) / (1im * v - v^2)
    # Option Price
    weights_ = @. 3 + (-1)^((0:(N-1)) + 1)# Simpson weights
    @views weights_[1] = 1
    @views weights_[end] = 1
    dk = 2 * pi / A
    b = N * dk / 2
    weights_cu = cu(weights_)
    v_cu = cu(v)
    complex_vec_z_k = @. integrand_f(v_cu) * dx * weights_cu / 3 * exp(1im * b * v_cu)
    fft!(complex_vec_z_k)
    z_T = collect(@. real_mod(complex_vec_z_k) / pi)

    ks = -b .+ dk * real_vec
    spline_cub = CubicSplineInterpolation(ks, z_T)
    k_interp = @. log(StrikeVec / S0)
    C = @. S0 * spline_cub(k_interp) + max(S0 - StrikeVec * exp(-(r - d) * T), 0)
    return C * exp(-d * T)
end


text:
using Mamba

## Data
oxford = Dict{Symbol, Any}(
  :r1 =>
    [3, 5, 2, 7, 7, 2, 5, 3, 5, 11, 6, 6, 11, 4, 4, 2, 8, 8, 6, 5, 15, 4, 9, 9,
     4, 12, 8, 8, 6, 8, 12, 4, 7, 16, 12, 9, 4, 7, 8, 11, 5, 12, 8, 17, 9, 3, 2,
     7, 6, 5, 11, 14, 13, 8, 6, 4, 8, 4, 8, 7, 15, 15, 9, 9, 5, 6, 3, 9, 12, 14,
     16, 17, 8, 8, 9, 5, 9, 11, 6, 14, 21, 16, 6, 9, 8, 9, 8, 4, 11, 11, 6, 9,
     4, 4, 9, 9, 10, 14, 6, 3, 4, 6, 10, 4, 3, 3, 10, 4, 10, 5, 4, 3, 13, 1, 7,
     5, 7, 6, 3, 7],
  :n1 =>
    [28, 21, 32, 35, 35, 38, 30, 43, 49, 53, 31, 35, 46, 53, 61, 40, 29, 44, 52,
     55, 61, 31, 48, 44, 42, 53, 56, 71, 43, 43, 43, 40, 44, 70, 75, 71, 37, 31,
     42, 46, 47, 55, 63, 91, 43, 39, 35, 32, 53, 49, 75, 64, 69, 64, 49, 29, 40,
     27, 48, 43, 61, 77, 55, 60, 46, 28, 33, 32, 46, 57, 56, 78, 58, 52, 31, 28,
     46, 42, 45, 63, 71, 69, 43, 50, 31, 34, 54, 46, 58, 62, 52, 41, 34, 52, 63,
     59, 88, 62, 47, 53, 57, 74, 68, 61, 45, 45, 62, 73, 53, 39, 45, 51, 55, 41,
     53, 51, 42, 46, 54, 32],
  :r0 =>
    [0, 2, 2, 1, 2, 0, 1, 1, 1, 2, 4, 4, 2, 1, 7, 4, 3, 5, 3, 2, 4, 1, 4, 5, 2,
     7, 5, 8, 2, 3, 5, 4, 1, 6, 5, 11, 5, 2, 5, 8, 5, 6, 6, 10, 7, 5, 5, 2, 8,
     1, 13, 9, 11, 9, 4, 4, 8, 6, 8, 6, 8, 14, 6, 5, 5, 2, 4, 2, 9, 5, 6, 7, 5,
     10, 3, 2, 1, 7, 9, 13, 9, 11, 4, 8, 2, 3, 7, 4, 7, 5, 6, 6, 5, 6, 9, 7, 7,
     7, 4, 2, 3, 4, 10, 3, 4, 2, 10, 5, 4, 5, 4, 6, 5, 3, 2, 2, 4, 6, 4, 1],
  :n0 =>
    [28, 21, 32, 35, 35, 38, 30, 43, 49, 53, 31, 35, 46, 53, 61, 40, 29, 44, 52,
     55, 61, 31, 48, 44, 42, 53, 56, 71, 43, 43, 43, 40, 44, 70, 75, 71, 37, 31,
     42, 46, 47, 55, 63, 91, 43, 39, 35, 32, 53, 49, 75, 64, 69, 64, 49, 29, 40,
     27, 48, 43, 61, 77, 55, 60, 46, 28, 33, 32, 46, 57, 56, 78, 58, 52, 31, 28,
     46, 42, 45, 63, 71, 69, 43, 50, 31, 34, 54, 46, 58, 62, 52, 41, 34, 52, 63,
     59, 88, 62, 47, 53, 57, 74, 68, 61, 45, 45, 62, 73, 53, 39, 45, 51, 55, 41,
     53, 51, 42, 46, 54, 32],
  :year =>
    [-10, -9, -9, -8, -8, -8, -7, -7, -7, -7, -6, -6, -6, -6, -6, -5, -5, -5,
     -5, -5, -5, -4, -4, -4, -4, -4, -4, -4, -3, -3, -3, -3, -3, -3, -3, -3, -2,
     -2, -2, -2, -2, -2, -2, -2, -2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 0,
     0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2,
     2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6,
     6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 9, 9, 10],
  :K => 120
)
oxford[:K] = length(oxford[:r1])

## Model Specification
model = Model(

  r0 = Stochastic(1,
    (mu, n0, K) ->
      begin
        p = invlogit(mu)
        UnivariateDistribution[Binomial(n0[i], p[i]) for i in 1:K]
      end,
    false
  ),

  r1 = Stochastic(1,
    (mu, alpha, beta1, beta2, year, b, n1, K) ->
      UnivariateDistribution[
        begin
          p = invlogit(mu[i] + alpha + beta1 * year[i] +
                       beta2 * (year[i]^2 - 22.0) + b[i])
          Binomial(n1[i], p)
        end
        for i in 1:K
      ],
    false
  ),

  b = Stochastic(1,
    s2 -> Normal(0, sqrt(s2)),
    false
  ),

  mu = Stochastic(1,
    () -> Normal(0, 1000),
    false
  ),

  alpha = Stochastic(
    () -> Normal(0, 1000)
  ),

  beta1 = Stochastic(
    () -> Normal(0, 1000)
  ),

  beta2 = Stochastic(
    () -> Normal(0, 1000)
  ),

  s2 = Stochastic(
    () -> InverseGamma(0.001, 0.001)
  )

)

## Initial Values
inits = [
  Dict(:r0 => oxford[:r0], :r1 => oxford[:r1], :alpha => 0, :beta1 => 0,
       :beta2 => 0, :s2 => 1, :b => zeros(oxford[:K]),
       :mu => zeros(oxford[:K])),
  Dict(:r0 => oxford[:r0], :r1 => oxford[:r1], :alpha => 1, :beta1 => 1,
       :beta2 => 1, :s2 => 10, :b => zeros(oxford[:K]),
       :mu => zeros(oxford[:K]))
]

## Sampling Scheme
scheme = [AMWG([:alpha, :beta1, :beta2], 1.0),
          Slice(:s2, 1.0),
          Slice(:mu, 1.0),
          Slice(:b, 1.0)]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, oxford, inits, 12500, burnin=2500, thin=2, chains=2)
describe(sim)


text:
using Yao, Yao.YaoBlocks.ConstGate, BenchmarkTools
using JSON
using LinearAlgebra, Pkg
using CuYao, CuArrays

BLAS.set_num_threads(1)

const nqubits=4:25
const benchmarks = Dict()

layer(nbit::Int, x::Symbol) = layer(nbit, Val(x))
layer(nbit::Int, ::Val{:first}) = chain(nbit, put(i=>chain(Rx(0), Rz(0))) for i = 1:nbit);
layer(nbit::Int, ::Val{:last}) = chain(nbit, put(i=>chain(Rz(0), Rx(0))) for i = 1:nbit)
layer(nbit::Int, ::Val{:mid}) = chain(nbit, put(i=>chain(Rz(0), Rx(0), Rz(0))) for i = 1:nbit);
entangler(pairs) = chain(control(ctrl, target=>X) for (ctrl, target) in pairs);
function build_circuit(n, nlayers, pairs)
    circuit = chain(n)
    push!(circuit, layer(n, :first))
    for i in 2:nlayers
        push!(circuit, entangler(pairs))
        push!(circuit, layer(n, :mid))
    end
    push!(circuit, entangler(pairs))
    push!(circuit, layer(n, :last))
    return circuit
end

qft_rotation(i, j) = control(i, j=>shift(2π/(1<<(i-j+1))))
qft_layer(n, k) = chain(n, j==k ? put(k=>H) : qft_rotation(j, k) for j in k:n)
qft_reverse(n) = chain(qft_layer(n, k) for k in 1:n)

function qft(n)
	circuit = chain(n)
	push!(circuit, qft_reverse(n))

	for i in 2:n/2
		push!(circuit, swap(i, n-i-1))
	end
	return circuit
end

qft(n) = chain(qft_reverse(n))

macro task(name::String, nqubits_ex, body)
    nqubits = nqubits_ex.args[2]
    msg = "benchmarking $name"
    quote
        @info $msg
        benchmarks[$(name)] = Dict()
        benchmarks[$(name)]["nqubits"] = $(esc(nqubits))
        benchmarks[$(name)]["times"] = $(esc(body))
    end
end

@task "X" nqubits=nqubits begin
    map(nqubits) do k
        t = @benchmark apply!(st, $(put(k, 2=>X))) setup=(st=rand_state($k))
        minimum(t).time
    end
end

@task "H" nqubits=nqubits begin
    map(nqubits) do k
        t = @benchmark apply!(st, $(put(k, 2=>H))) setup=(st=rand_state($k))
        minimum(t).time
    end
end

@task "T" nqubits=nqubits begin
    map(nqubits) do k
        t = @benchmark apply!(st, $(put(k, 2=>T))) setup=(st=rand_state($k))
        minimum(t).time
    end
end

@task "CNOT" nqubits=nqubits begin
    map(nqubits) do k
        t = @benchmark apply!(st, $(control(k, 2, 3=>X))) setup=(st=rand_state($k))
        minimum(t).time
    end
end

@task "Rx(0.5)" nqubits=nqubits begin
    map(nqubits) do k
        t = @benchmark apply!(st, $(put(k, 2=>Rx(0.5)))) setup=(st=rand_state($k))
        minimum(t).time
    end
end

@task "Rz(0.5)" nqubits=nqubits begin
    map(nqubits) do k
        t = @benchmark apply!(st, $(put(k, 2=>Rz(0.5)))) setup=(st=rand_state($k))
        minimum(t).time
    end
end

@task "CRx(0.5)" nqubits=nqubits begin
    map(nqubits) do k
        t = @benchmark apply!(st, $(control(k, 2, 3=>Rx(0.5)))) setup=(st=rand_state($k))
        minimum(t).time
    end
end

@task "Toffoli" nqubits=nqubits begin
    map(nqubits) do k
        t = @benchmark apply!(st, $(control(k, (2, 3), 1=>X))) setup=(st=rand_state($k))
        minimum(t).time
    end
end

const qft_nqubits = 4:5

@task "QFT" nqubits=qft_nqubits begin
    map(qft_nqubits) do k
        t = @benchmark apply!(st, $(qft(k))) setup=(st=zero_state($k))
        minimum(t).time
    end
end

@task "QFT (batch)" nqubits=qft_nqubits begin
    map(qft_nqubits) do k
        t = @benchmark apply!(st, $(qft(k))) setup=(st=zero_state($k, nbatch=1000))
        minimum(t).time
    end
end

@static if CuArrays.functional()

    @task "QFT (cuda)" nqubits=qft_nqubits begin
        map(qft_nqubits) do k
            t = @benchmark CuArrays.@sync(apply!(st, $(qft(k)))) setup=(st=cu(zero_state($k)))
            minimum(t).time
        end
    end

    @task "QFT (batch) (cuda)" nqubits=4:15 begin
        map(4:15) do k
            t = @benchmark CuArrays.@sync(apply!(st, $(qft(k)))) setup=(st=cu(zero_state($k, nbatch=1000)))
            minimum(t).time
        end
    end

end

const qcbm_nqubits = 4:25

@task "QCBM" nqubits=qcbm_nqubits begin
    map(qcbm_nqubits) do k
        t = @benchmark apply!(st, $(build_circuit(k, 9, [(i, mod1(i+1, k)) for i in 1:k]))) setup=(st=zero_state($k))
        minimum(t).time
    end
end

@task "QCBM (batch)" nqubits=4:15 begin
    map(4:15) do k
        t = @benchmark apply!(st, $(build_circuit(k, 9, [(i, mod1(i+1, k)) for i in 1:k]))) setup=(st=zero_state($k, nbatch=1000))
        minimum(t).time
    end
end

@static if CuArrays.functional()

    @task "QCBM (cuda)" nqubits=qcbm_nqubits begin
        map(qcbm_nqubits) do k
            t = @benchmark CuArrays.@sync(apply!(st, $(build_circuit(k, 9, [(i, mod1(i+1, k)) for i in 1:k])))) setup=(st=cu(zero_state($k)))
            minimum(t).time
        end
    end

    @task "QCBM (batch) (cuda)" nqubits=4:15 begin
        map(4:15) do k
            t = @benchmark CuArrays.@sync(apply!(st, $(build_circuit(k, 9, [(i, mod1(i+1, k)) for i in 1:k])))) setup=(st=cu(zero_state($k, nbatch=1000)))
            minimum(t).time
        end
    end

end

if !ispath("data")
    mkpath("data")
end

write("data/data.json", JSON.json(benchmarks))


text:
# Dunavant, D. (1985). High degree efficient symmetrical Gaussian quadrature rules
# for the triangle. International journal for numerical methods in engineering, 21(6),
# 1129-1148.
using LinearAlgebra
using Optim
using StaticArrays
using BenchmarkTools

setprecision(BigFloat, 512)
const T = BigFloat

# Polynomial degree
const p = parse(Int64, ARGS[1])

# m values from Table I.
const m_all = SVector(1, 2, 3, 4, 5, 7, 8, 10, 12, 14, 16, 19, 21, 24, 27, 30, 33, 37, 40, 44)
const m = m_all[p]

# Polar moments from Table I.
const ν_all = SVector(
        +1//1,
        +1//4,
        -1//10,
        +1//10,
        -2//35,
       +29//560,
        +1//28,
        -1//28,
       +11//350,
        +1//40,
       -37//1540,
        -1//55,
       +13//616,
        +1//55,
       -49//2860,
        -2//143,
      +425//28028,
      +137//10010,
        +1//91,
       -64//5005,
        -1//91,
      +523//45760,
       +85//8008,
        +1//112,
     -6733//680680,
      -109//12376,
        -1//136,
      +217//24310,
      +209//24752,
        +1//136,
     -2909//369512,
       -65//9044,
        -2//323,
    +66197//9237800,
     +8069//1175720,
      +317//51680,
        +1//190,
     -3769//587860,
       -77//12920,
        -1//190,
    +83651//14226212,
    +11303//1989680,
       +92//17765,
        +1//220)
const ν = SVector{m, T}(ν_all[1:m])

# Polar moment indices from Table I.
const iν_all = SMatrix{44, 2, Int64, 88}(
    0, 2, 3, 4, 5, 6, 6, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 12, 13, 13, 14, 14, 14, 15, 15, 15, 16, 16, 16, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19, 20, 20, 20, 20,
    0, 0, 3, 0, 3, 0, 6, 3, 0, 6, 3, 9,  0,  6,  3,  9,  0,  6, 12,  3,  9,  0,  6, 12,  3,  9, 15,  0,  6, 12,  3,  9, 15,  0,  6, 12, 18,  3,  9, 15,  0,  6, 12, 18
   )
const iν = SMatrix{m, 2, Int64, 2m}(iν_all[1:m, :])

# n₀, n₁, n₂ values provided in Table IV,
# the number of groups of multiplicity 1, 3, and 6
const nᵢ= SVector((1, 0, 0),
                  (0, 1, 0),
                  (1, 1, 0),
                  (0, 2, 0),
                  (1, 2, 0),
                  (0, 2, 1),
                  (1, 2, 1),
                  (1, 3, 1),
                  (1, 4, 1),
                  (1, 2, 3),
                  (0, 5, 2),
                  (0, 5, 3),
                  (1, 6, 3),
                  (0, 6, 4),
                  (0, 6, 5),
                  (1, 7, 5),
                  (1, 8, 6),
                  (1, 9, 7),
                  (1, 8, 8),
                  (1, 10, 8))
const n₀, n₁, n₂ = nᵢ[p]
# calculate the number of unknowns
const n = n₀ + 2n₁ + 3n₂
# calculate the number of gaussian points and weights
const ng = n₀ + 3n₁ + 6n₂

# Given α, β, and weights in Appendix II
if p === 1
    const αβw_Float64 = @SVector [
      (0.333333333333333, 0.333333333333333, 0.500000000000000)]
elseif p === 2
    const αβw_Float64 = @SVector [
      (0.666666666666667, 0.166666666666667, 0.166666666666667)]
elseif p === 3
    const αβw_Float64 = @SVector [
      (0.333333333333333, 0.333333333333333, -0.281250000000000),
      (0.600000000000000, 0.200000000000000, 0.260416666666667)]
elseif p === 4
    const αβw_Float64 = @SVector [
      (0.108103018168070, 0.445948490915965, 0.111690794839005),
      (0.816847572980459, 0.091576213509771, 0.054975871827661)]
elseif p === 5
    const αβw_Float64 = @SVector [
      (0.333333333333333, 0.333333333333333, 0.112500000000000),
      (0.059715871789770, 0.470142064105115, 0.066197076394253),
      (0.797426985353087, 0.101286507323456, 0.062969590272414)]
elseif p === 6
    const αβw_Float64 = @SVector [
      (0.501426509658179, 0.249286745170910, 0.058393137863189),
      (0.873821971016996, 0.063089014491502, 0.025422453185104),
      (0.053145049844816, 0.636502499121399, 0.041425537809187)]
elseif p === 7
    const αβw_Float64 = @SVector [
      (0.333333333333333, 0.333333333333333, -0.074785022233841),
      (0.479308067841920, 0.260345966079040, 0.087807628716604),
      (0.869739794195568, 0.065130102902216, 0.026673617804419),
      (0.048690315425316, 0.312865496004874, 0.038556880445128)]
elseif p === 8
    const αβw_Float64 = @SVector [
      (0.333333333333333, 0.333333333333333, 0.072157803838894),
      (0.081414823414554, 0.459292588292723, 0.047545817133643),
      (0.658861384496480, 0.170569307751760, 0.051608685267359),
      (0.898905543365938, 0.050547228317031, 0.016229248811599),
      (0.008394777409958, 0.263112829634638, 0.013615157087217)]
elseif p === 9
    const αβw_Float64 = @SVector [
      (0.3333333333330, 0.3333333333333, 0.04856789814140),
      (0.0206349616025, 0.4896825191990, 0.01566735011355),
      (0.1258208170140, 0.4370895914930, 0.03891377050240),
      (0.6235929287620, 0.1882035356190, 0.03982386946360),
      (0.9105409732110, 0.0447295133945, 0.01278883782935),
      (0.0368384120547, 0.2219629891613, 0.02164176968865)]
elseif p === 10
    const αβw_Float64 = @SVector [
      (0.3333333333330, 0.3333333333330, 0.04540899519140),
      (0.0288447332327, 0.4855776333840, 0.01836297887825),
      (0.7810368490300, 0.1094815754850, 0.02266052971775),
      (0.1417072194150, 0.5503529418210, 0.03637895842270),
      (0.0250035347627, 0.7283239045970, 0.01416362126555),
      (0.0095408154003, 0.9236559335870, 0.00471083348185)]
elseif p === 11
    const αβw_Float64 = @SVector [
      (-0.0692220965415, +0.5346110482710, 0.00046350316448),
      (+0.2020613940680, +0.3989693029660, 0.03857476745740),
      (+0.5933801991370, +0.2033099004310, 0.02966148869040),
      (+0.7612981754350, +0.1193509122830, 0.01809227025170),
      (+0.9352701037770, +0.0323649481113, 0.00682986550135),
      (+0.0501781383105, +0.3566206482610, 0.02616855598110),
      (+0.0210220165362, +0.1714889803040, 0.01035382981955)]
elseif p === 12
    const αβw_Float64 = @SVector [
      (0.0235652204524, 0.4882173897740, 0.01286553322025),
      (0.1205512154110, 0.4397243922940, 0.02184627226900),
      (0.4575792299760, 0.2712103850120, 0.03142911210895),
      (0.7448477089170, 0.1275761455420, 0.01739805646535),
      (0.9573652990940, 0.0213173504532, 0.00308313052578),
      (0.1153434945350, 0.2757132696860, 0.02018577888320),
      (0.0228383322223, 0.6958360867880, 0.01117838660115),
      (0.0257340505483, 0.1162519159080, 0.00865811555435)]
elseif p === 13
    const αβw_Float64 = @SVector [
      (0.33333333333330, 0.33333333333330, 0.026260461700400),
      (0.00990363012059, 0.49504818494000, 0.005640072604650),
      (0.06256672978090, 0.46871663511000, 0.015711759181250),
      (0.17095732639700, 0.41452133680100, 0.023536251252100),
      (0.54120085591400, 0.22939957204300, 0.023681793268200),
      (0.77115100960700, 0.11442449519600, 0.015583764522900),
      (0.95037721727300, 0.02481139136350, 0.003987885732535),
      (0.09485382837960, 0.26879499705900, 0.018424201364350),
      (0.01810077327880, 0.29173006673400, 0.008700731651900),
      (0.02223307667410, 0.12635738549200, 0.007760893419500)]
elseif p === 14
    const αβw_Float64 = @SVector [
      (0.02207217927560, 0.48896391036200, 0.01094179068470),
      (0.16471056131900, 0.41764471934000, 0.01639417677205),
      (0.45304494338200, 0.27347752830900, 0.02588705225365),
      (0.64558893517500, 0.17720553241300, 0.02108129436850),
      (0.87640023381800, 0.06179988309090, 0.00721684983490),
      (0.96121807750300, 0.01939096124870, 0.00246170180120),
      (0.05712475740360, 0.17226668782100, 0.01233287660630),
      (0.09291624935700, 0.33686145979600, 0.01928575539355),
      (0.01464695005570, 0.29837288213600, 0.00721815405675),
      (0.00126833093287, 0.11897449769700, 0.00250511441925)]
elseif p === 15
    const αβw_Float64 = @SVector [
      (-0.0139458337165, +0.5069729168580, 0.000958437821425),
      (+0.1371872914340, +0.4314063542830, 0.022124513635550),
      (+0.4446127103060, +0.2776936448470, 0.025593274359450),
      (+0.7470702179170, +0.1264648910410, 0.011843867935350),
      (+0.8583832280510, +0.0708083859747, 0.006644887845000),
      (+0.9620696595180, +0.0189651702411, 0.002374458304095),
      (+0.1337341619670, +0.2613113711400, 0.019275036299800),
      (+0.0363666773969, +0.3880467670900, 0.013607907160300),
      (-0.0101748831266, +0.2857122200500, 0.001091038683400),
      (+0.0368438698759, +0.2155996640720, 0.010752659923850),
      (+0.0124598093312, +0.1035756165760, 0.003836971315525)]
elseif p === 16
    const αβw_Float64 = @SVector [
      (+0.33333333333330, +0.33333333333330, 0.023437848713800),
      (+0.00523891610312, +0.49738054194800, 0.003202939289290),
      (+0.17306112290100, +0.41346943854900, 0.020855148369700),
      (+0.05908280186600, +0.47045859906700, 0.013445742125050),
      (+0.51889250006100, +0.24055374997000, 0.021066261380800),
      (+0.70406841155500, +0.14796579422300, 0.015000133421400),
      (+0.84906962468500, +0.07546518765750, 0.007100049462500),
      (+0.96680719475400, +0.01659640262300, 0.001791231175635),
      (+0.10357569224500, +0.29655559658000, 0.016386573730300),
      (+0.02008341165540, +0.33772306340300, 0.007649153124200),
      (-0.00434100261414, +0.20474828164300, 0.001193122096420),
      (+0.04194178646800, +0.18935849213100, 0.009542396377950),
      (+0.01431732023070, +0.08528361568270, 0.003425027273270)]
elseif p === 17
    const αβw_Float64 = @SVector [
      (0.33333333333330, 0.33333333333330, 0.016718599645400),
      (0.00565891888645, 0.49717054055700, 0.002546707720255),
      (0.03564735475080, 0.48217632262500, 0.007335432263800),
      (0.09952006195840, 0.45023996902100, 0.012175439176850),
      (0.19946752124500, 0.40026623937700, 0.015553775434500),
      (0.49571746405800, 0.25214126797100, 0.015628555609300),
      (0.67590599068300, 0.16204700465800, 0.012407827169850),
      (0.84824823547900, 0.07587588226070, 0.007028036535300),
      (0.96869054606400, 0.01565472696780, 0.001597338086890),
      (0.01018692882690, 0.33431986736400, 0.004059827659495),
      (0.13544087167100, 0.29222153779700, 0.013402871141600),
      (0.05442392429060, 0.31957488542300, 0.009229996605400),
      (0.01286856083360, 0.19070422419200, 0.004238434267165),
      (0.06716578241350, 0.18048321164900, 0.009146398385000),
      (0.01466318222480, 0.08071131367960, 0.003332816002085)]
elseif p === 18
    const αβw_Float64 = @SVector [
      (+0.33333333333330, +0.33333333333330, +0.015404969968800),
      (+0.01331038273820, +0.49334480863100, +0.004536218339700),
      (+0.06157881151610, +0.46921059424200, +0.009380658469800),
      (+0.12743720822600, +0.43628139588700, +0.009720548992750),
      (+0.21030765865300, +0.39484617067300, +0.013876974305400),
      (+0.50041086239400, +0.24979456880300, +0.016128112675750),
      (+0.67713561251200, +0.16143219374400, +0.012537016308450),
      (+0.84680354502900, +0.07659822748540, +0.007635963985900),
      (+0.95149512129300, +0.02425243935350, +0.003396961011480),
      (+0.91370726556600, +0.04314636721700, -0.001111549364960),
      (+0.00843053620242, +0.35891149494100, +0.003165957038205),
      (+0.13118655173700, +0.29440247675200, +0.013628769024550),
      (+0.05020315156570, +0.32501780164200, +0.008838392824750),
      (+0.06632926381090, +0.18473755966600, +0.009189742319050),
      (+0.01199619456620, +0.21879680001300, +0.004052366404095),
      (+0.01485810059010, +0.10117959713600, +0.003817064535365),
      (-0.03522201528790, +0.02087475528260, +2.3093830397e-05)]
elseif p === 19
    const αβw_Float64 = @SVector [
      (0.33333333333330, 0.33333333333330, 0.016453165694450),
      (0.02078002585400, 0.48960998707300, 0.005165365945650),
      (0.09092621460420, 0.45453689269800, 0.011193623631500),
      (0.19716663870100, 0.40141668064900, 0.015133062934750),
      (0.48889669119400, 0.25555165440300, 0.015245483901100),
      (0.64584411569600, 0.17707794215200, 0.012079606370800),
      (0.77987789354400, 0.11006105322800, 0.008025401793400),
      (0.88894275149600, 0.05552862425180, 0.004042290130890),
      (0.97475627244600, 0.01262186377720, 0.001039681013740),
      (0.00361141784841, 0.39575478735700, 0.001942438452490),
      (0.13446675453100, 0.30792998388000, 0.012787080306000),
      (0.01444602577610, 0.26456694840700, 0.004440451786670),
      (0.04693357883820, 0.35853935220600, 0.008062273380850),
      (0.00286112035057, 0.15780740596900, 0.001245970908745),
      (0.22386142409800, 0.07505059697590, 0.009121420059500),
      (0.03464707481680, 0.14242160111300, 0.005129281868100),
      (0.01016111929630, 0.06549462808290, 0.001899964427650)]
elseif p === 20
    const αβw_Float64 = @SVector [
      (+0.3333333333333, +0.3333333333333, +0.016528527770800),
      (-0.0019009287044, +0.5009504643520, +0.000433509592831),
      (+0.0235740841305, +0.4882129579350, +0.005830026358200),
      (+0.0897266360994, +0.4551366869500, +0.011438468178200),
      (+0.1960074813630, +0.4019962593180, +0.015224491336950),
      (+0.4882141804810, +0.2558929097590, +0.015312445862700),
      (+0.6470234880100, +0.1764882559950, +0.012184028838400),
      (+0.7916582893260, +0.1041708553370, +0.007998716016000),
      (+0.8938620723180, +0.0530689638409, +0.003849150907800),
      (+0.9167625696080, +0.0416187151960, -0.000316030248744),
      (+0.9768361571860, +0.0115819214068, +0.000875567150595),
      (+0.0487415836648, +0.3448557702290, +0.008232919594800),
      (+0.0063141159486, +0.3778432695950, +0.002419516770245),
      (+0.1343165205470, +0.3066354790620, +0.012902453267350),
      (+0.0139738939624, +0.2494193627750, +0.004235545527220),
      (+0.0755491329098, +0.2127757248030, +0.009177457053150),
      (-0.0083681532082, +0.1469654360530, +0.000352202338954),
      (+0.0266860632587, +0.1377269788290, +0.005056342463750),
      (+0.0105477192941, +0.0596961091490, +0.001786954692975)]
else
    error("Unsupported polynomial order")
end

# Given α and β, find r and a in the area function.
# x = [r, a]
function αβ_to_ra_error(x, α, β)
    r = x[1]
    a = x[2]
    α_err = (1 - 2r*cos(a))/3 - α
    β_err = (1 + r*cos(a) - sqrt(T(3))*r*sin(a))/3 - β
    return α_err^2 + β_err^2
end

function αβw_to_initial_guess()
    initial_guess = SizedVector{n, T}(zeros(T, n))
    icount = 1
    if n₀ === 1
        α,β,w = αβw_Float64[1]
        initial_guess[1] = 2w
        icount += 1
    end
    for i = n₀ + 1 : 2 : n₀ + 2n₁
        α,β,w = αβw_Float64[icount]
        initial_guess[i]   = 6w
        initial_guess[i+1] = (1 - 3α)/2
        icount += 1
    end
    for i = n₀ + 2n₁ + 1: 3 : n
        α,β,w = αβw_Float64[icount]
        initial_guess[i] = 12w
        res = optimize(x -> αβ_to_ra_error(x, α, β),
                       [0.0, 0.0],
                       Optim.Options(g_tol=1e-30,
                                     iterations=Int64(1e4),
                                    )
                      )
        if !Optim.converged(res)
            println("r, a from α, β")
            println(res)
        end
        x = Optim.minimizer(res)
        initial_guess[i+1] = x[1]
        initial_guess[i+2] = x[2]
        icount += 1
    end
    return initial_guess
end

# The error in equation 22b for moment iM given
# x = {w₀ if n₀ = 1} ∪ {wᵢ, rᵢ for i = 1:n₁} ∪ {wᵢ, rᵢ, 3αᵢ for i = n₁+1:n₂}
function moment_error(x::SizedVector{n, T, Vector{T}}, iM::Int64)
    err = -ν[iM]
    j = iν[iM, 1]
    k = iν[iM, 2]
    if n₀ === 1 && iM === 1
        err += x[1]
    end
    for i = n₀ + 1 : 2 : n₀ + 2n₁
        err += x[i]*x[i+1]^j
    end
    for i = n₀ + 2n₁ + 1: 3 : n
        err += x[i]*x[i+1]^j * cos(k*x[i+2])
    end
    return err
end

# The sum or the squared error in each unknown
function objective(x::SizedVector{n, T, Vector{T}})
    sum_err = zero(T)
    for iM = 1:m
        sum_err += moment_error(x, iM)^2
    end
    return sum_err
end

# Convert to barycentric coordinates
function area(x)
    nsum = n₀ + n₁ + n₂
    w = zeros(T, nsum)
    α = zeros(T, nsum)
    β = zeros(T, nsum)
    γ = zeros(T, nsum)
    icount = 1
    if !(n₀ === 0)
        w[icount] = x[1]
        α[icount] = T(1//3)
        β[icount] = T(1//3)
        γ[icount] = T(1//3)
        icount += 1
    end
    for i = n₀ + 1 : 2 : n₀ + 2n₁
        w[icount] = x[i]/3
        α[icount] = (1 - 2x[i+1])/3
        β[icount] = (1 - α[icount])/2
        γ[icount] = β[icount]
        icount += 1
    end
    for i = n₀ + 2n₁ + 1: 3 : n
        r = x[i+1]
        a = x[i+2]
        w[icount] = x[i]/6
        α[icount] = (1 - 2r*cos(a))/3
        β[icount] = (1 + r*cos(a) - sqrt(T(3))*r*sin(a))/3
        γ[icount] = 1 - α[icount] - β[icount]
        icount += 1
    end
    return [ (w[i], α[i], β[i], γ[i]) for i ∈ 1:nsum ]
end

# Integration error function from Appendix II
function poly(wαβγ)
    exact = 2/(T(p+1)*T(p+2))
    sum = zero(T)
    if n₀ !== 0
        w, α, β, γ = wαβγ[1]
        sum += w*α^p
    end
    for i = n₀ + 1:n₀ + n₁
        w, α, β, γ = wαβγ[i]
        sum += w*(α^p + 2*β^p)
    end
    for i = n₀ + n₁ + 1: n₀ + n₁ + n₂
        w, α, β, γ = wαβγ[i]
        sum += w*2*(α^p + β^p + γ^p)
    end
    approx = sum
    err = (exact - approx)/approx
    return err
end

initial_guess = αβw_to_initial_guess()
println("Polynomial degree: $p")
println("nᵢ = $n₀, $n₁, $n₂")
@time res = optimize(objective,
                     initial_guess,
                     Optim.Options(g_tol=1e-200,
                                   iterations=Int64(1e6),
                                   show_trace=true,
                                   show_every=10000
                                  )
                    )
x = Optim.minimizer(res)
if !Optim.converged(res)
    println(res)
end
wαβγ = area(x)
for t in wαβγ
    println(Float64.(t))
end
println("w, α, β, γ errors")
for i = 1:length(wαβγ)
    α_ref,β_ref,w_ref = αβw_Float64[i]
    γ_ref = 1 - α_ref - β_ref
    w, α, β, γ = wαβγ[i]
    w_err = Float64(w)-2w_ref
    α_err = Float64(α)-α_ref
    β_err = Float64(β)-β_ref
    γ_err = Float64(γ)-γ_ref
    errs = (w_err, α_err, β_err, γ_err)
    for err in errs
        if abs(err) < 1e-14
            print("OK ")
        else
            print(err, " ")
        end
    end
    println()
end
println("Integration error")
println(poly(wαβγ))

weights = T[]
points = NTuple{2, T}[]
if n₀ !== 0
    push!(weights, wαβγ[1][1]/2)
    push!(points, (wαβγ[1][2], wαβγ[1][3]))
end
for i = n₀ + 1:n₀ + n₁
    w, α, β, γ = wαβγ[i]
    for j = 1:3
        push!(weights, w/2)
    end
    push!(points, (β, β))
    push!(points, (β, α))
    push!(points, (α, β))
end
for i = n₀ + n₁ + 1: n₀ + n₁ + n₂
    w, α, β, γ = wαβγ[i]
    for j = 1:6
        push!(weights, w/2)
    end
    push!(points, (γ, β))
    push!(points, (β, α))
    push!(points, (α, γ))
    push!(points, (γ, α))
    push!(points, (β, γ))
    push!(points, (α, β))
end

io = open("legendre_triangle.jl", "a")
println(io, "@generated function gauss_quadrature(form::Val{:legendre},")
println(io, "                                     shape::RefTriangle,")
println(io, "                                     degree::Val{$p},")
println(io, "                                     type::Type{T}) where {T}")
println(io, """    weights = SVector(:(\$(T(big"$(weights[1])"))),""")
for i = 2:ng
println(io, """                      :(\$(T(big"$(weights[i])"))),""")
end
println(io, "                       )")
println(io, """    points = SVector(:(\$(tuple(T(big"$(points[1][1])"), T(big"$(points[1][2])")))),""")
for i in 2:ng
println(io, """                     :(\$(tuple(T(big"$(points[i][1])"), T(big"$(points[i][2])")))),""")
end
println(io, "                      )")
println(io, "    return weights, points")
println(io, "end")

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
# based on the contents of align.jl
# consider changing terminology from align to calibrate.

function alignregionκ(y_cost::Vector{Complex{T}},
    U_cost,
    P_cost,
    As,
    Bs,
    Es,
    fs,
    SW,
    Δsys_cs,
    a_setp, b_setp, κs_β_DOFs, κs_β_orderings,
    shift_lb::Vector{T},
    shift_ub::Vector{T};
    w = ones(T, length(As)),
    LS_inds = 1:length(U_cost),
    N_starts = 100,
    local_optim_algorithm = NLopt.LN_BOBYQA,
    xtol_rel = 1e-3,
    maxeval = 50,
    maxtime = Inf,
    β_optim_algorithm = :GN_DIRECT_L,
    κ_lb_default = 1e-2,
    κ_ub_default = 1e2,
    β_max_iters = 500,
    β_xtol_rel = 1e-9,
    β_ftol_rel = 1e-9,
    β_maxtime = Inf) where T <: Real

    # prepare.
    N_d = sum( getNd(Bs[n]) for n = 1:length(Bs) )
    @assert length(shift_ub) == length(shift_lb) == N_d

    U_rad_cost = U_cost .* (2*π)

    # setup inner optim over β.
    q, updatedfunc, getshiftfunc, N_vars_set,
    run_optim, obj_func_β, E_BLS, κ_BLS, b_BLS, updateβfunc, updateκfunc,
    q_β = setupcostnesteddwarp(Es, Bs, As, fs, SW, LS_inds, U_rad_cost,
        y_cost, Δsys_cs, a_setp, b_setp, κs_β_DOFs, κs_β_orderings;
        w = w,
        optim_algorithm = β_optim_algorithm,
        κ_lb_default = κ_lb_default,
        κ_ub_default = κ_ub_default,
        max_iters = β_max_iters,
        xtol_rel = β_xtol_rel,
        ftol_rel = β_ftol_rel,
        maxtime = β_maxtime)

    # set up outer optim over shifts.
    #N_β = sum( getNβ(κs_β_DOFs[n], Bs[n]) for n = 1:length(Bs) )
    N_β = sum( getNβ(Bs[n]) for n = 1:length(Bs) )
    p_β = zeros(T, N_β) # persistant buffer.

    obj_func = pp->costnestedd(U_rad_cost, y_cost, updatedfunc, updateκfunc, pp,
    #Es, Bs, κs_β_orderings, κs_β_DOFs, q, run_optim, E_BLS, κ_BLS, b_BLS, p_β)
    Es, Bs, run_optim, E_BLS, κ_BLS, b_BLS, p_β)

    # optim.
    prob = MultistartOptimization.MinimizationProblem(obj_func, shift_lb, shift_ub)

    local_method = MultistartOptimization.NLoptLocalMethod(; algorithm = local_optim_algorithm,
    xtol_rel = xtol_rel,
    maxeval = maxeval,
    maxtime = maxtime)

    multistart_method = MultistartOptimization.TikTak(N_starts)
    ret_mo = MultistartOptimization.multistart_minimization(multistart_method,
        local_method, prob)
    #
    println("ret_mo = ", ret_mo)
    minf = ret_mo.value
    minx = ret_mo.location
    ret = :None
    if hasproperty(ret_mo, :ret)
        ret = ret_mo.ret
    end

    return obj_func, minf, minx, ret
end

"""
Requires user to supply a surrogate.
"""
function aligncompoundκ(y::Vector{Complex{T}}, U_y, P_y, As, Bs, Es, fs, SW,
    Δsys_cs, a_setp, b_setp, κs_β_DOFs, κs_β_orderings::Vector{Vector{Vector{Int}}},
    shift_lb::Vector{T},
    shift_ub::Vector{T},
    cost_inds_set::Vector{Vector{Int}};
    loop_range = 1:length(cost_inds_set),
    w = ones(T, length(As)),
    N_starts = 100,
    local_optim_algorithm = NLopt.LN_BOBYQA,
    xtol_rel = 1e-3,
    maxeval = 50,
    maxtime = Inf,
    β_optim_algorithm = :GN_DIRECT_L,
    κ_lb_default = 1e-2,
    κ_ub_default = 1e2,
    β_max_iters = 500,
    β_xtol_rel = 1e-9,
    β_ftol_rel = 1e-9,
    β_maxtime = Inf) where T <: Real

    N_regions = length(loop_range)
    minxs = Vector{Vector{T}}(undef, N_regions)
    rets = Vector{Symbol}(undef, N_regions)
    minfs = Vector{T}(undef, N_regions)
    obj_funcs = Vector{Function}(undef, N_regions)

    for r in loop_range
        println("Working on region $(r)")
        obj_funcs[r], minfs[r], minxs[r], rets[r] = alignregionκ(y[cost_inds_set[r]],
            U_y[cost_inds_set[r]],
            P_y[cost_inds_set[r]],
            As,
            Bs,
            Es,
            fs,
            SW,
            Δsys_cs,
            a_setp, b_setp, κs_β_DOFs, κs_β_orderings,
            shift_lb,
            shift_ub;
            w = w,
            N_starts = N_starts,
            local_optim_algorithm = local_optim_algorithm,
            xtol_rel = xtol_rel,
            maxeval = maxeval,
            maxtime = maxtime,
            β_optim_algorithm = β_optim_algorithm,
            κ_lb_default = κ_lb_default,
            κ_ub_default = κ_ub_default,
            β_max_iters = β_max_iters,
            β_xtol_rel = β_xtol_rel,
            β_ftol_rel = β_ftol_rel,
            β_maxtime = β_maxtime)
    end

    return obj_funcs, minfs, minxs, rets
end

function aligncompoundsingleregion(y::Vector{Complex{T}}, U_y, P_y, As, Bs, Es, fs, SW,
    Δsys_cs, a_setp, b_setp, #κs_β_DOFs, κs_β_orderings::Vector{Vector{Vector{Int}}},
    shift_lb::Vector{T},
    shift_ub::Vector{T},
    cost_inds_set::Vector{Vector{Int}};
    loop_range = 1:length(cost_inds_set),
    w = ones(T, length(As)),
    N_starts = 100,
    local_optim_algorithm = NLopt.LN_BOBYQA,
    xtol_rel = 1e-3,
    maxeval = 50,
    maxtime = Inf,
    β_optim_algorithm = :GN_DIRECT_L,
    κ_lb_default = 1e-2,
    κ_ub_default = 1e2,
    β_max_iters = 500,
    β_xtol_rel = 1e-9,
    β_ftol_rel = 1e-9,
    β_maxtime = Inf) where T <: Real

    tmp = collect( cost_inds_set[r] for r = 1:length(cost_inds_set) )
    inds_redundant = NMRSpectraSimulator.combinevectors(tmp)
    inds = inds_redundant[indexin( unique(inds_redundant), inds_redundant)]

    U_cost = U_y[inds]
    P_cost = P_y[inds]
    y_cost = y[inds]

    obj_func, minf, minx, ret = alignregion(y_cost,
        U_cost,
        P_cost,
        As,
        Bs,
        Es,
        fs,
        SW,
        Δsys_cs,
        a_setp, b_setp, #κs_β_DOFs, κs_β_orderings,
        shift_lb,
        shift_ub;
        w = w,
        N_starts = 100,
        local_optim_algorithm = NLopt.LN_BOBYQA,
        xtol_rel = 1e-3,
        maxeval = 50,
        maxtime = Inf,
        β_optim_algorithm = :GN_DIRECT_L,
        κ_lb_default = 1e-2,
        κ_ub_default = 1e2,
        β_max_iters = 500,
        β_xtol_rel = 1e-9,
        β_ftol_rel = 1e-9,
        β_maxtime = Inf)

    return obj_func, minf, minx, ret, P_cost, y_cost
end


text:
@testset "AbstractLatent" begin
    kernel = SqExponentialKernel()
    T = Float64
    D = 5
    m0 = ZeroMean()
    opt = ADAM()
    X = rand(1, D)
    K = kernelmatrix(kernel, X; obsdim=1)
    @testset "LatentGP" begin
        m = AGP.LatentGP(T, D, kernel, m0, opt)
        @test m isa AGP.LatentGP
        @test AGP.prior(m) isa AGP.GPPrior
        @test AGP.posterior(m) isa AGP.Posterior
        # Get functions
        @test AGP.kernel(m) == kernel
        @test AGP.pr_mean(m) == m0
        @test AGP.pr_mean(m, vec(X)) == zeros(D)
        @test AGP.pr_cov(m) == PDMat(Matrix{T}(I(D)))
        @test mean(m) == zeros(D)
        @test cov(m) == PDMat(Matrix{T}(I(D)))
        AGP.pr_cov!(m, PDMat(K))
        @test AGP.pr_cov(m) == K
        @test dim(m) == D
    end
    @testset "VarLatent" begin
        m = AGP.VarLatent(T, D, kernel, m0, opt)
        @test m isa AGP.VarLatent
        @test AGP.prior(m) isa AGP.GPPrior
        @test AGP.posterior(m) isa AGP.VarPosterior
        # Get functions
        @test AGP.kernel(m) == kernel
        @test AGP.pr_mean(m) == m0
        @test AGP.pr_mean(m, vec(X)) == zeros(D)
        @test AGP.pr_cov(m) == PDMat(Matrix{T}(I(D)))
        @test mean(m) == zeros(D)
        @test cov(m) == PDMat(Matrix{T}(I(D)))
        AGP.pr_cov!(m, PDMat(K))
        @test AGP.pr_cov(m) == K
        @test dim(m) == D
    end
end


text:
using Colors
using Images
using Boltzmann
using MNIST
using ImageView
using Base.Test

function run_mnist()
    X, y = testdata()  # test data is smaller, no need to downsample
    binarize!(X;threhold=0.01)

    HiddenUnits = 256;
    Epochs = 10;
    LearningRate = 0.1;

    m = BernoulliRBM(28*28, HiddenUnits, (28,28); momentum=0.0, dataset=X)
    mM = BernoulliRBM(28*28, HiddenUnits, (28,28); momentum=0.5, dataset=X)

    # Attempt without Momentum
    info("Running Without Momentum")
    fit(m, X;n_iter=Epochs,lr=LearningRate)

    # Attempt with Momentum
    info("Running With Momentum")
    fit(mM, X;n_iter=Epochs,lr=LearningRate)

    return m
end

run_mnist()


text:
export nodeFluxPlot

function nodeFluxPlot( mpSim::MPsim, timeGrid::Vector{T}, fluxType::Symbol,
    nodes::AbstractString...; plotType::Symbol=:normal, showPlot::Bool=true,
    savePlot::Bool=false, filename::AbstractString="",
    timeFactor::Real=12.0 ) where T <: Real

    savePlot |= filename != ""
    plotType = haskey( plotTypes, plotType ) ? plotType : :normal
    plotFunction = plotTypes[plotType][1]

    # Don't do anything if the plot mustn't be shown or saved.
    if !( showPlot || savePlot )
        return
    end  # if !( showPlot || savePlot )

    # Generate the report.
    fluxReport = nodeFluxReport( mpSim, timeGrid, fluxType, nodes... )

    if isempty( fluxReport )
        return
    end  # if isempty( fluxReport )

    generateNodeFluxPlot( fluxReport, nodes, fluxType, timeFactor, showPlot,
        savePlot, filename, plotFunction )

end  # nodeFluxPlot( mpSim, timeGrid, fluxType, nodes, plotType, showPlot,
     #   savePlot, filename, timeFactor )

nodeFluxPlot( mpSim::MPsim, timeRes::Real, fluxType::Symbol,
    nodes::AbstractString...; plotType::Symbol=:normal, showPlot::Bool=true,
    savePlot::Bool=false, filename::AbstractString="", timeFactor::Real=12.0 ) =
    nodeFluxPlot( mpSim, MP.generateTimeGrid( mpSim, timeRes ), fluxType,
    nodes..., plotType=plotType, showPlot=showPlot, savePlot=savePlot,
    filename=filename, timeFactor=timeFactor )

include( "private/nodeflux.jl" )


text:
using Distributions
using LogExpFunctions
using PSIS
using Test

@testset "effective sample size" begin
    w = fill(0.01, 100)
    reff = rand()
    @test ess_is(w) ≈ 100
    @test ess_is(w; reff=reff) ≈ 100 .* reff

    w = zeros(100)
    i = rand(1:100)
    w[i] = 1
    @test ess_is(w) ≈ 1
    @test ess_is(w; reff=reff) ≈ 1 .* reff

    logw = randn(100)
    logw_norm = logsumexp(logw)
    result = PSISResult(logw, logw_norm, 1.5, 20, PSIS.GeneralizedPareto(0.0, 1.0, 0.6))
    @test ess_is(result) ≈ ess_is(exp.(logw .- logw_norm); reff=1.5)

    result = PSISResult(logw, logw_norm, 1.5, 20, PSIS.GeneralizedPareto(0.0, 1.0, 0.71))
    @test ismissing(ess_is(result))
    @test ess_is(result; bad_shape_missing=false) ≈
        ess_is(exp.(logw .- logw_norm); reff=1.5)

    logw = randn(3, 100, 4)
    logw_norm = dropdims(logsumexp(logw; dims=(2, 3)); dims=(2, 3))
    tail_dists = [
        PSIS.GeneralizedPareto(0.0, 1.0, 0.69),
        PSIS.GeneralizedPareto(0.0, 1.0, 0.71),
        missing,
    ]
    reff = [1.5, 0.8, 1.0]
    result = PSISResult(logw, logw_norm, reff, [20, 20, 20], tail_dists)
    ess = ess_is(result)
    @test ess isa Vector
    @test length(ess) == 3
    @test ess[1] ≈ ess_is(result.weights; reff=reff)[1]
    @test ismissing(ess[2])
    @test ismissing(ess[3])
    ess = ess_is(result; bad_shape_missing=false)
    @test ess[1:2] ≈ ess_is(result.weights; reff=reff)[1:2]
    @test ismissing(ess[3])
end


text:
#Official site: https://julialang.org/

#Simple syntax, no need to include any libraries for terminal output
println("Hello world!")

#From Julia REPL
# julia> include("/path/to/hello_world.jl")

#From terminal
# $ /path/to/julia hello_world.jl

#Output: Hello world!

#Here is some more information: https://juliabyexample.helpmanual.io/


text:
const s0 = (
	0xd1310ba6, 0x98dfb5ac, 0x2ffd72db, 0xd01adfb7, 0xb8e1afed, 0x6a267e96,
	0xba7c9045, 0xf12c7f99, 0x24a19947, 0xb3916cf7, 0x0801f2e2, 0x858efc16,
	0x636920d8, 0x71574e69, 0xa458fea3, 0xf4933d7e, 0x0d95748f, 0x728eb658,
	0x718bcd58, 0x82154aee, 0x7b54a41d, 0xc25a59b5, 0x9c30d539, 0x2af26013,
	0xc5d1b023, 0x286085f0, 0xca417918, 0xb8db38ef, 0x8e79dcb0, 0x603a180e,
	0x6c9e0e8b, 0xb01e8a3e, 0xd71577c1, 0xbd314b27, 0x78af2fda, 0x55605c60,
	0xe65525f3, 0xaa55ab94, 0x57489862, 0x63e81440, 0x55ca396a, 0x2aab10b6,
	0xb4cc5c34, 0x1141e8ce, 0xa15486af, 0x7c72e993, 0xb3ee1411, 0x636fbc2a,
	0x2ba9c55d, 0x741831f6, 0xce5c3e16, 0x9b87931e, 0xafd6ba33, 0x6c24cf5c,
	0x7a325381, 0x28958677, 0x3b8f4898, 0x6b4bb9af, 0xc4bfe81b, 0x66282193,
	0x61d809cc, 0xfb21a991, 0x487cac60, 0x5dec8032, 0xef845d5d, 0xe98575b1,
	0xdc262302, 0xeb651b88, 0x23893e81, 0xd396acc5, 0x0f6d6ff3, 0x83f44239,
	0x2e0b4482, 0xa4842004, 0x69c8f04a, 0x9e1f9b5e, 0x21c66842, 0xf6e96c9a,
	0x670c9c61, 0xabd388f0, 0x6a51a0d2, 0xd8542f68, 0x960fa728, 0xab5133a3,
	0x6eef0b6c, 0x137a3be4, 0xba3bf050, 0x7efb2a98, 0xa1f1651d, 0x39af0176,
	0x66ca593e, 0x82430e88, 0x8cee8619, 0x456f9fb4, 0x7d84a5c3, 0x3b8b5ebe,
	0xe06f75d8, 0x85c12073, 0x401a449f, 0x56c16aa6, 0x4ed3aa62, 0x363f7706,
	0x1bfedf72, 0x429b023d, 0x37d0d724, 0xd00a1248, 0xdb0fead3, 0x49f1c09b,
	0x075372c9, 0x80991b7b, 0x25d479d8, 0xf6e8def7, 0xe3fe501a, 0xb6794c3b,
	0x976ce0bd, 0x04c006ba, 0xc1a94fb6, 0x409f60c4, 0x5e5c9ec2, 0x196a2463,
	0x68fb6faf, 0x3e6c53b5, 0x1339b2eb, 0x3b52ec6f, 0x6dfc511f, 0x9b30952c,
	0xcc814544, 0xaf5ebd09, 0xbee3d004, 0xde334afd, 0x660f2807, 0x192e4bb3,
	0xc0cba857, 0x45c8740f, 0xd20b5f39, 0xb9d3fbdb, 0x5579c0bd, 0x1a60320a,
	0xd6a100c6, 0x402c7279, 0x679f25fe, 0xfb1fa3cc, 0x8ea5e9f8, 0xdb3222f8,
	0x3c7516df, 0xfd616b15, 0x2f501ec8, 0xad0552ab, 0x323db5fa, 0xfd238760,
	0x53317b48, 0x3e00df82, 0x9e5c57bb, 0xca6f8ca0, 0x1a87562e, 0xdf1769db,
	0xd542a8f6, 0x287effc3, 0xac6732c6, 0x8c4f5573, 0x695b27b0, 0xbbca58c8,
	0xe1ffa35d, 0xb8f011a0, 0x10fa3d98, 0xfd2183b8, 0x4afcb56c, 0x2dd1d35b,
	0x9a53e479, 0xb6f84565, 0xd28e49bc, 0x4bfb9790, 0xe1ddf2da, 0xa4cb7e33,
	0x62fb1341, 0xcee4c6e8, 0xef20cada, 0x36774c01, 0xd07e9efe, 0x2bf11fb4,
	0x95dbda4d, 0xae909198, 0xeaad8e71, 0x6b93d5a0, 0xd08ed1d0, 0xafc725e0,
	0x8e3c5b2f, 0x8e7594b7, 0x8ff6e2fb, 0xf2122b64, 0x8888b812, 0x900df01c,
	0x4fad5ea0, 0x688fc31c, 0xd1cff191, 0xb3a8c1ad, 0x2f2f2218, 0xbe0e1777,
	0xea752dfe, 0x8b021fa1, 0xe5a0cc0f, 0xb56f74e8, 0x18acf3d6, 0xce89e299,
	0xb4a84fe0, 0xfd13e0b7, 0x7cc43b81, 0xd2ada8d9, 0x165fa266, 0x80957705,
	0x93cc7314, 0x211a1477, 0xe6ad2065, 0x77b5fa86, 0xc75442f5, 0xfb9d35cf,
	0xebcdaf0c, 0x7b3e89a0, 0xd6411bd3, 0xae1e7e49, 0x00250e2d, 0x2071b35e,
	0x226800bb, 0x57b8e0af, 0x2464369b, 0xf009b91e, 0x5563911d, 0x59dfa6aa,
	0x78c14389, 0xd95a537f, 0x207d5ba2, 0x02e5b9c5, 0x83260376, 0x6295cfa9,
	0x11c81968, 0x4e734a41, 0xb3472dca, 0x7b14a94a, 0x1b510052, 0x9a532915,
	0xd60f573f, 0xbc9bc6e4, 0x2b60a476, 0x81e67400, 0x08ba6fb5, 0x571be91f,
	0xf296ec6b, 0x2a0dd915, 0xb6636521, 0xe7b9f9b6, 0xff34052e, 0xc5855664,
	0x53b02d5d, 0xa99f8fa1, 0x08ba4799, 0x6e85076a,
)::NTuple{256, UInt32}

const s1 = (
	0x4b7a70e9, 0xb5b32944, 0xdb75092e, 0xc4192623, 0xad6ea6b0, 0x49a7df7d,
	0x9cee60b8, 0x8fedb266, 0xecaa8c71, 0x699a17ff, 0x5664526c, 0xc2b19ee1,
	0x193602a5, 0x75094c29, 0xa0591340, 0xe4183a3e, 0x3f54989a, 0x5b429d65,
	0x6b8fe4d6, 0x99f73fd6, 0xa1d29c07, 0xefe830f5, 0x4d2d38e6, 0xf0255dc1,
	0x4cdd2086, 0x8470eb26, 0x6382e9c6, 0x021ecc5e, 0x09686b3f, 0x3ebaefc9,
	0x3c971814, 0x6b6a70a1, 0x687f3584, 0x52a0e286, 0xb79c5305, 0xaa500737,
	0x3e07841c, 0x7fdeae5c, 0x8e7d44ec, 0x5716f2b8, 0xb03ada37, 0xf0500c0d,
	0xf01c1f04, 0x0200b3ff, 0xae0cf51a, 0x3cb574b2, 0x25837a58, 0xdc0921bd,
	0xd19113f9, 0x7ca92ff6, 0x94324773, 0x22f54701, 0x3ae5e581, 0x37c2dadc,
	0xc8b57634, 0x9af3dda7, 0xa9446146, 0x0fd0030e, 0xecc8c73e, 0xa4751e41,
	0xe238cd99, 0x3bea0e2f, 0x3280bba1, 0x183eb331, 0x4e548b38, 0x4f6db908,
	0x6f420d03, 0xf60a04bf, 0x2cb81290, 0x24977c79, 0x5679b072, 0xbcaf89af,
	0xde9a771f, 0xd9930810, 0xb38bae12, 0xdccf3f2e, 0x5512721f, 0x2e6b7124,
	0x501adde6, 0x9f84cd87, 0x7a584718, 0x7408da17, 0xbc9f9abc, 0xe94b7d8c,
	0xec7aec3a, 0xdb851dfa, 0x63094366, 0xc464c3d2, 0xef1c1847, 0x3215d908,
	0xdd433b37, 0x24c2ba16, 0x12a14d43, 0x2a65c451, 0x50940002, 0x133ae4dd,
	0x71dff89e, 0x10314e55, 0x81ac77d6, 0x5f11199b, 0x043556f1, 0xd7a3c76b,
	0x3c11183b, 0x5924a509, 0xf28fe6ed, 0x97f1fbfa, 0x9ebabf2c, 0x1e153c6e,
	0x86e34570, 0xeae96fb1, 0x860e5e0a, 0x5a3e2ab3, 0x771fe71c, 0x4e3d06fa,
	0x2965dcb9, 0x99e71d0f, 0x803e89d6, 0x5266c825, 0x2e4cc978, 0x9c10b36a,
	0xc6150eba, 0x94e2ea78, 0xa5fc3c53, 0x1e0a2df4, 0xf2f74ea7, 0x361d2b3d,
	0x1939260f, 0x19c27960, 0x5223a708, 0xf71312b6, 0xebadfe6e, 0xeac31f66,
	0xe3bc4595, 0xa67bc883, 0xb17f37d1, 0x018cff28, 0xc332ddef, 0xbe6c5aa5,
	0x65582185, 0x68ab9802, 0xeecea50f, 0xdb2f953b, 0x2aef7dad, 0x5b6e2f84,
	0x1521b628, 0x29076170, 0xecdd4775, 0x619f1510, 0x13cca830, 0xeb61bd96,
	0x0334fe1e, 0xaa0363cf, 0xb5735c90, 0x4c70a239, 0xd59e9e0b, 0xcbaade14,
	0xeecc86bc, 0x60622ca7, 0x9cab5cab, 0xb2f3846e, 0x648b1eaf, 0x19bdf0ca,
	0xa02369b9, 0x655abb50, 0x40685a32, 0x3c2ab4b3, 0x319ee9d5, 0xc021b8f7,
	0x9b540b19, 0x875fa099, 0x95f7997e, 0x623d7da8, 0xf837889a, 0x97e32d77,
	0x11ed935f, 0x16681281, 0x0e358829, 0xc7e61fd6, 0x96dedfa1, 0x7858ba99,
	0x57f584a5, 0x1b227263, 0x9b83c3ff, 0x1ac24696, 0xcdb30aeb, 0x532e3054,
	0x8fd948e4, 0x6dbc3128, 0x58ebf2ef, 0x34c6ffea, 0xfe28ed61, 0xee7c3c73,
	0x5d4a14d9, 0xe864b7e3, 0x42105d14, 0x203e13e0, 0x45eee2b6, 0xa3aaabea,
	0xdb6c4f15, 0xfacb4fd0, 0xc742f442, 0xef6abbb5, 0x654f3b1d, 0x41cd2105,
	0xd81e799e, 0x86854dc7, 0xe44b476a, 0x3d816250, 0xcf62a1f2, 0x5b8d2646,
	0xfc8883a0, 0xc1c7b6a3, 0x7f1524c3, 0x69cb7492, 0x47848a0b, 0x5692b285,
	0x095bbf00, 0xad19489d, 0x1462b174, 0x23820e00, 0x58428d2a, 0x0c55f5ea,
	0x1dadf43e, 0x233f7061, 0x3372f092, 0x8d937e41, 0xd65fecf1, 0x6c223bdb,
	0x7cde3759, 0xcbee7460, 0x4085f2a7, 0xce77326e, 0xa6078084, 0x19f8509e,
	0xe8efd855, 0x61d99735, 0xa969a7aa, 0xc50c06c2, 0x5a04abfc, 0x800bcadc,
	0x9e447a2e, 0xc3453484, 0xfdd56705, 0x0e1e9ec9, 0xdb73dbd3, 0x105588cd,
	0x675fda79, 0xe3674340, 0xc5c43465, 0x713e38d8, 0x3d28f89e, 0xf16dff20,
	0x153e21e7, 0x8fb03d4a, 0xe6e39f2b, 0xdb83adf7,
)::NTuple{256, UInt32}

const s2 = (
	0xe93d5a68, 0x948140f7, 0xf64c261c, 0x94692934, 0x411520f7, 0x7602d4f7,
	0xbcf46b2e, 0xd4a20068, 0xd4082471, 0x3320f46a, 0x43b7d4b7, 0x500061af,
	0x1e39f62e, 0x97244546, 0x14214f74, 0xbf8b8840, 0x4d95fc1d, 0x96b591af,
	0x70f4ddd3, 0x66a02f45, 0xbfbc09ec, 0x03bd9785, 0x7fac6dd0, 0x31cb8504,
	0x96eb27b3, 0x55fd3941, 0xda2547e6, 0xabca0a9a, 0x28507825, 0x530429f4,
	0x0a2c86da, 0xe9b66dfb, 0x68dc1462, 0xd7486900, 0x680ec0a4, 0x27a18dee,
	0x4f3ffea2, 0xe887ad8c, 0xb58ce006, 0x7af4d6b6, 0xaace1e7c, 0xd3375fec,
	0xce78a399, 0x406b2a42, 0x20fe9e35, 0xd9f385b9, 0xee39d7ab, 0x3b124e8b,
	0x1dc9faf7, 0x4b6d1856, 0x26a36631, 0xeae397b2, 0x3a6efa74, 0xdd5b4332,
	0x6841e7f7, 0xca7820fb, 0xfb0af54e, 0xd8feb397, 0x454056ac, 0xba489527,
	0x55533a3a, 0x20838d87, 0xfe6ba9b7, 0xd096954b, 0x55a867bc, 0xa1159a58,
	0xcca92963, 0x99e1db33, 0xa62a4a56, 0x3f3125f9, 0x5ef47e1c, 0x9029317c,
	0xfdf8e802, 0x04272f70, 0x80bb155c, 0x05282ce3, 0x95c11548, 0xe4c66d22,
	0x48c1133f, 0xc70f86dc, 0x07f9c9ee, 0x41041f0f, 0x404779a4, 0x5d886e17,
	0x325f51eb, 0xd59bc0d1, 0xf2bcc18f, 0x41113564, 0x257b7834, 0x602a9c60,
	0xdff8e8a3, 0x1f636c1b, 0x0e12b4c2, 0x02e1329e, 0xaf664fd1, 0xcad18115,
	0x6b2395e0, 0x333e92e1, 0x3b240b62, 0xeebeb922, 0x85b2a20e, 0xe6ba0d99,
	0xde720c8c, 0x2da2f728, 0xd0127845, 0x95b794fd, 0x647d0862, 0xe7ccf5f0,
	0x5449a36f, 0x877d48fa, 0xc39dfd27, 0xf33e8d1e, 0x0a476341, 0x992eff74,
	0x3a6f6eab, 0xf4f8fd37, 0xa812dc60, 0xa1ebddf8, 0x991be14c, 0xdb6e6b0d,
	0xc67b5510, 0x6d672c37, 0x2765d43b, 0xdcd0e804, 0xf1290dc7, 0xcc00ffa3,
	0xb5390f92, 0x690fed0b, 0x667b9ffb, 0xcedb7d9c, 0xa091cf0b, 0xd9155ea3,
	0xbb132f88, 0x515bad24, 0x7b9479bf, 0x763bd6eb, 0x37392eb3, 0xcc115979,
	0x8026e297, 0xf42e312d, 0x6842ada7, 0xc66a2b3b, 0x12754ccc, 0x782ef11c,
	0x6a124237, 0xb79251e7, 0x06a1bbe6, 0x4bfb6350, 0x1a6b1018, 0x11caedfa,
	0x3d25bdd8, 0xe2e1c3c9, 0x44421659, 0x0a121386, 0xd90cec6e, 0xd5abea2a,
	0x64af674e, 0xda86a85f, 0xbebfe988, 0x64e4c3fe, 0x9dbc8057, 0xf0f7c086,
	0x60787bf8, 0x6003604d, 0xd1fd8346, 0xf6381fb0, 0x7745ae04, 0xd736fccc,
	0x83426b33, 0xf01eab71, 0xb0804187, 0x3c005e5f, 0x77a057be, 0xbde8ae24,
	0x55464299, 0xbf582e61, 0x4e58f48f, 0xf2ddfda2, 0xf474ef38, 0x8789bdc2,
	0x5366f9c3, 0xc8b38e74, 0xb475f255, 0x46fcd9b9, 0x7aeb2661, 0x8b1ddf84,
	0x846a0e79, 0x915f95e2, 0x466e598e, 0x20b45770, 0x8cd55591, 0xc902de4c,
	0xb90bace1, 0xbb8205d0, 0x11a86248, 0x7574a99e, 0xb77f19b6, 0xe0a9dc09,
	0x662d09a1, 0xc4324633, 0xe85a1f02, 0x09f0be8c, 0x4a99a025, 0x1d6efe10,
	0x1ab93d1d, 0x0ba5a4df, 0xa186f20f, 0x2868f169, 0xdcb7da83, 0x573906fe,
	0xa1e2ce9b, 0x4fcd7f52, 0x50115e01, 0xa70683fa, 0xa002b5c4, 0x0de6d027,
	0x9af88c27, 0x773f8641, 0xc3604c06, 0x61a806b5, 0xf0177a28, 0xc0f586e0,
	0x006058aa, 0x30dc7d62, 0x11e69ed7, 0x2338ea63, 0x53c2dd94, 0xc2c21634,
	0xbbcbee56, 0x90bcb6de, 0xebfc7da1, 0xce591d76, 0x6f05e409, 0x4b7c0188,
	0x39720a3d, 0x7c927c24, 0x86e3725f, 0x724d9db9, 0x1ac15bb4, 0xd39eb8fc,
	0xed545578, 0x08fca5b5, 0xd83d7cd3, 0x4dad0fc4, 0x1e50ef5e, 0xb161e6f8,
	0xa28514d9, 0x6c51133c, 0x6fd5c7e7, 0x56e14ec4, 0x362abfce, 0xddc6c837,
	0xd79a3234, 0x92638212, 0x670efa8e, 0x406000e0,
)::NTuple{256, UInt32}

const s3 = (
	0x3a39ce37, 0xd3faf5cf, 0xabc27737, 0x5ac52d1b, 0x5cb0679e, 0x4fa33742,
	0xd3822740, 0x99bc9bbe, 0xd5118e9d, 0xbf0f7315, 0xd62d1c7e, 0xc700c47b,
	0xb78c1b6b, 0x21a19045, 0xb26eb1be, 0x6a366eb4, 0x5748ab2f, 0xbc946e79,
	0xc6a376d2, 0x6549c2c8, 0x530ff8ee, 0x468dde7d, 0xd5730a1d, 0x4cd04dc6,
	0x2939bbdb, 0xa9ba4650, 0xac9526e8, 0xbe5ee304, 0xa1fad5f0, 0x6a2d519a,
	0x63ef8ce2, 0x9a86ee22, 0xc089c2b8, 0x43242ef6, 0xa51e03aa, 0x9cf2d0a4,
	0x83c061ba, 0x9be96a4d, 0x8fe51550, 0xba645bd6, 0x2826a2f9, 0xa73a3ae1,
	0x4ba99586, 0xef5562e9, 0xc72fefd3, 0xf752f7da, 0x3f046f69, 0x77fa0a59,
	0x80e4a915, 0x87b08601, 0x9b09e6ad, 0x3b3ee593, 0xe990fd5a, 0x9e34d797,
	0x2cf0b7d9, 0x022b8b51, 0x96d5ac3a, 0x017da67d, 0xd1cf3ed6, 0x7c7d2d28,
	0x1f9f25cf, 0xadf2b89b, 0x5ad6b472, 0x5a88f54c, 0xe029ac71, 0xe019a5e6,
	0x47b0acfd, 0xed93fa9b, 0xe8d3c48d, 0x283b57cc, 0xf8d56629, 0x79132e28,
	0x785f0191, 0xed756055, 0xf7960e44, 0xe3d35e8c, 0x15056dd4, 0x88f46dba,
	0x03a16125, 0x0564f0bd, 0xc3eb9e15, 0x3c9057a2, 0x97271aec, 0xa93a072a,
	0x1b3f6d9b, 0x1e6321f5, 0xf59c66fb, 0x26dcf319, 0x7533d928, 0xb155fdf5,
	0x03563482, 0x8aba3cbb, 0x28517711, 0xc20ad9f8, 0xabcc5167, 0xccad925f,
	0x4de81751, 0x3830dc8e, 0x379d5862, 0x9320f991, 0xea7a90c2, 0xfb3e7bce,
	0x5121ce64, 0x774fbe32, 0xa8b6e37e, 0xc3293d46, 0x48de5369, 0x6413e680,
	0xa2ae0810, 0xdd6db224, 0x69852dfd, 0x09072166, 0xb39a460a, 0x6445c0dd,
	0x586cdecf, 0x1c20c8ae, 0x5bbef7dd, 0x1b588d40, 0xccd2017f, 0x6bb4e3bb,
	0xdda26a7e, 0x3a59ff45, 0x3e350a44, 0xbcb4cdd5, 0x72eacea8, 0xfa6484bb,
	0x8d6612ae, 0xbf3c6f47, 0xd29be463, 0x542f5d9e, 0xaec2771b, 0xf64e6370,
	0x740e0d8d, 0xe75b1357, 0xf8721671, 0xaf537d5d, 0x4040cb08, 0x4eb4e2cc,
	0x34d2466a, 0x0115af84, 0xe1b00428, 0x95983a1d, 0x06b89fb4, 0xce6ea048,
	0x6f3f3b82, 0x3520ab82, 0x011a1d4b, 0x277227f8, 0x611560b1, 0xe7933fdc,
	0xbb3a792b, 0x344525bd, 0xa08839e1, 0x51ce794b, 0x2f32c9b7, 0xa01fbac9,
	0xe01cc87e, 0xbcc7d1f6, 0xcf0111c3, 0xa1e8aac7, 0x1a908749, 0xd44fbd9a,
	0xd0dadecb, 0xd50ada38, 0x0339c32a, 0xc6913667, 0x8df9317c, 0xe0b12b4f,
	0xf79e59b7, 0x43f5bb3a, 0xf2d519ff, 0x27d9459c, 0xbf97222c, 0x15e6fc2a,
	0x0f91fc71, 0x9b941525, 0xfae59361, 0xceb69ceb, 0xc2a86459, 0x12baa8d1,
	0xb6c1075e, 0xe3056a0c, 0x10d25065, 0xcb03a442, 0xe0ec6e0e, 0x1698db3b,
	0x4c98a0be, 0x3278e964, 0x9f1f9532, 0xe0d392df, 0xd3a0342b, 0x8971f21e,
	0x1b0a7441, 0x4ba3348c, 0xc5be7120, 0xc37632d8, 0xdf359f8d, 0x9b992f2e,
	0xe60b6f47, 0x0fe3f11d, 0xe54cda54, 0x1edad891, 0xce6279cf, 0xcd3e7e6f,
	0x1618b166, 0xfd2c1d05, 0x848fd2c5, 0xf6fb2299, 0xf523f357, 0xa6327623,
	0x93a83531, 0x56cccd02, 0xacf08162, 0x5a75ebb5, 0x6e163697, 0x88d273cc,
	0xde966292, 0x81b949d0, 0x4c50901b, 0x71c65614, 0xe6c6c7bd, 0x327a140a,
	0x45e1d006, 0xc3f27b9a, 0xc9aa53fd, 0x62a80f00, 0xbb25bfe2, 0x35bdd2f6,
	0x71126905, 0xb2040222, 0xb6cbcf7c, 0xcd769c2b, 0x53113ec0, 0x1640e3d3,
	0x38abbd60, 0x2547adf0, 0xba38209c, 0xf746ce76, 0x77afa1c5, 0x20756060,
	0x85cbfe4e, 0x8ae88dd8, 0x7aaaf9b0, 0x4cf9aa7e, 0x1948c25c, 0x02fb8a8c,
	0x01c36ae4, 0xd6ebe1f9, 0x90d4f869, 0xa65cdea0, 0x3f09252d, 0xc208e69f,
	0xb74e6132, 0xce77e25b, 0x578fdfe3, 0x3ac372e6,
)::NTuple{256, UInt32}

const p = (
	0x243f6a88, 0x85a308d3, 0x13198a2e, 0x03707344, 0xa4093822, 0x299f31d0,
	0x082efa98, 0xec4e6c89, 0x452821e6, 0x38d01377, 0xbe5466cf, 0x34e90c6c,
	0xc0ac29b7, 0xc97c50dd, 0x3f84d5b5, 0xb5470917, 0x9216d5d9, 0x8979fb1b,
)::NTuple{18, UInt32}


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Miscellaneous functions for the functions of IAU-2006 theory.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Compute the sum as used in many computations of IAU-2006 theory. The first
# parameter is a tuple with the coefficient matrices, the second is the Julian
# century [TT], and the rest are the fundamental arguments.
@inline function _iau2006_sum(
    coefs::Tuple, T_TT::Number, M_s::Number,
    M_m::Number,
    u_Mm::Number,
    D_s::Number,
    Ω_m::Number,
    λ_M☿::Number,
    λ_M♀::Number,
    λ_Me::Number,
    λ_M♂::Number,
    λ_M♃::Number,
    λ_M♄::Number,
    λ_M⛢::Number,
    λ_M♆::Number,
    p_λ::Number
)
    # Result of the sum.
    r = 0.0

    # Auxiliary variable to compute the powers of T_TT.
    T_TT_power = one(T_TT)

    # Number of sums.
    num_sums = length(coefs)

    @inbounds for i = 1:num_sums
        ci = coefs[i]

        # Result of this sum.
        rp = zero(r)

        # Number of coefficients in this sum.
        num_coefs = size(ci, 2)

        # Notice that the matrices were transposed when created to improve the
        # performance due to memory alignment.
        for j = 1:num_coefs
            As     = ci[ 2,j]
            Ac     = ci[ 3,j]
            ap     = ci[ 4,j] * M_m  + ci[ 5,j] * M_s  + ci[ 6,j] * u_Mm +
                     ci[ 7,j] * D_s  + ci[ 8,j] * Ω_m  + ci[ 9,j] * λ_M☿ +
                     ci[10,j] * λ_M♀ + ci[11,j] * λ_Me + ci[12,j] * λ_M♂ +
                     ci[13,j] * λ_M♃ + ci[14,j] * λ_M♄ + ci[15,j] * λ_M⛢ +
                     ci[16,j] * λ_M♆ + ci[17,j] * p_λ
            sj, cj = sincos(ap)
            rp    += (As * sj + Ac * cj) / 1e6
        end

        # Accumulate in the output variable.
        r += rp * T_TT_power

        # Update the T_TT power for the next pass.
        T_TT_power *= T_TT
    end

    return r
end


text:
module TestEspy
using EspyInsideFunction,Test#,StaticArrays

## Test request

struct Element  end
requestable(el::Element)= (X=(3,2),  gp=forloop(2, (F=(2,2), material=(ε=(2,2),Σ=(2,2)))) )
el = Element()
r8 = @request X,gp[].(F,material.(Σ,ε))

ra    = requestable(el)

@testset "requestable" begin
    @test ra.X == (3,2)
    @test ra.gp.body == (F = (2, 2), material = (ε = (2, 2), Σ = (2, 2)))
end

@testset "Makekey" begin
    @test EspyInsideFunction.makekey_symbol(2,(2,3)) == ([3 5 7; 4 6 8], 8)
    @test EspyInsideFunction.makekey_tuple(2,:((a,b)),(a=(2,3),b=(2,2))) == ((a=[3 5 7; 4 6 8], b=[9 11; 10 12]), 12)
    @test EspyInsideFunction.makekey_tuple(2,:((a,b)),(a=(2,3),b=scalar)) == ((a=[3 5 7; 4 6 8], b=9), 9)
    @test EspyInsideFunction.makekey_loop(2,:((a,b)),forloop(3,(a=(2,3),b=(2,2)))) == ([(a=[3 5 7; 4 6 8], b=[9 11; 10 12]), (a=[13 15 17; 14 16 18], b=[19 21; 20 22]), (a=[23 25 27; 24 26 28], b=[29 31; 30 32])], 32)
    @test EspyInsideFunction.makekey_tuple(2,:((a[].(b,)),),(a=forloop(2,(b=(2,3),)),)) == (( a=[(b=[3 5 7; 4 6 8],), (b=[9 11 13; 10 12 14],)] ,), 14)
    @test EspyInsideFunction.makekey_tuple(2,:(a[].b),(a=forloop(2,(b=(2,3),)),)) == ((a=[(b=[3 5 7; 4 6 8],), (b=[9 11 13; 10 12 14],)],), 14)
    @test makekey(r8,ra) == ((X = [1 4; 2 5; 3 6], gp = [(F = [7 9; 8 10], material = (Σ = [11 13; 12 14], ε = [15 17; 16 18])), (F = [19 21; 20 22], material = (Σ = [23 25; 24 26], ε = [27 29; 28 30]))]), 30)
end

## annotated code (written by user)

struct El  end
requestable(el::El)= (gp=forloop(2, (z=scalar,s=scalar, material=(a=scalar,b=scalar))),)
el = El()

@espy function residual(x,y)
    ngp=2
    r = 0
    for igp=1:ngp
        :z = x[igp]+y[igp]
        :s,dum  = :material(z)
        r += s
    end
    return r
end
@espy function material(z)
    :a = z+1
    :b = a*z
    return b,3.
end

## macro'ed code (returned by @espy)
# function residual(out,key,x,y)
#     ngp = 2
#     r   = 0
#     for igp = 1:ngp
#         @espy_loop key gp                     # key_gp = key.gp[igp]
#         z = x[igp]+y[igp]
#         @espy_record out key_gp z             # out[key_gp.z] = z
#         s = @espy_call out key_gp material(z) # s = material(out,key_gp.material,z)
#         @espy_record out key_gp s             # out[key_gp.s] = s
#         r += s
#     end
#     return r
# end
# function material(out,key,z)
#     a = z+1
#     @espy_record out key a                    # out[key.a] = a
#     b = a*z
#     @espy_record out key b                    # out[key.b] = b
#     return b
# end
## final code (after @espy_call, @espy_loop and @espy_record are run)
# function residual(out,key,x,y) #
#     ngp = 2
#     r = 0
#     for igp = 1:ngp
#         key_gp = key.gp[igp] #
#         z = x[igp]+y[igp]
#         out[key_gp.z] = z #
#         s  = material(out,key_gp.material,z) #
#         out[key_gp.s] = s #
#         r += s
#     end
#     return r
# end
# function material(out,key,z) #
#     a = z+1
#     out[key.a] = a #
#     b = a*z
#     out[key.b] = b #
#     return b
# end

## Test result extraction
req       = @request gp[].(s,z,material.(a,b))  # generates an expression, not a variable
key,nkey  = makekey(req,requestable(el))
nstep,nel = 2,3
out       = Array{Float64,3}(undef,nkey,nel,nstep)
istep,iel = 1,2
x,y       = [1.,2.],[.5,.2]
r         = residual(@view(out[:,iel,istep]),key,x,y)

@testset "Espy" begin
    @test key == (gp=[(s=1, z=2, material=(a=3, b=4)), (s=5, z=6, material=(a=7, b=8))],)
    @test out[:,iel,istep] ≈ [3.75, 1.5, 2.5, 3.75, 7.04, 2.2, 3.2, 7.04]
end

out = Vector{Float64}(undef,5)
out[[1,2]] .= [1,2]
out[[3,4]] .= (3,4)
# out[[5,6]] .= @SVector [5,6] works, but requires StaticArrays in the test environment
out[[5]]   .= 5
@testset "Assign to out" begin
    @test out == [1.,2.,3.,4.,5.]
end

end

max_stars_repo_name: PhilippeMaincon/Espy


text:
using DataDeps

register(DataDep(
  "skin-segmentation",
  "http://archive.ics.uci.edu/ml/datasets/Skin+Segmentation",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/00229/Skin_NonSkin.txt",
  "e30c0a845385dcc95a45c45ed263465674a49638e98ef740afd520769c7714a4",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=4,
        feature_indices=1:3,
        delim='\t',
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
module Multifractals

###exports first to show API

###type declarations

###functions

end # module


text:
using ChemistryFeaturization

# Batching Utilities
"""
    batch_graph_data(laplacians, encoded_features)

Takes vectors of laplacians and encoded features and joins them
into a single graph of disjoint subgraphs. The resulting graph
is massive and hence the return types are sparse. Few of the layers
don't work with Sparse Arrays (specifically on GPUs), so it would
make sense to convert them to dense.
"""
struct BatchedAtomicGraph{T1,T2,S}
    laplacians::T1
    encoded_features::T2
    sizes::S
end

Flux.@functor BatchedAtomicGraph (laplacians, encoded_features)

batch_graph_data(t::Tuple) = batch_graph_data(t[1], t[2])

function batch_graph_data(laplacians, encoded_features)
    # Ideally should be sparse arrays but currently doesn't work well on GPUs
    _sizes = map(x -> size(x, 1), laplacians)
    total_nodes = sum(_sizes)
    batched_laplacian = zeros(eltype(laplacians[1]), total_nodes, total_nodes)
    idx = 1
    for i in 1:length(laplacians)
        batched_laplacian[idx:(idx + _sizes[i] - 1), idx:(idx + _sizes[i] - 1)] .= laplacians[i]
        idx += _sizes[i]
    end
    _sizes = vcat(0, cumsum(_sizes))
    enc_feats = hcat(encoded_features...)
    return BatchedAtomicGraph(batched_laplacian, enc_feats, _sizes)
end

function BatchedAtomicGraph(batch_size::Int, atoms::Vector{FeaturizedAtoms})
    return BatchedAtomicGraph(batch_size, map(x -> x.atoms.laplacian, atoms), map(x -> x.encoded_features, atoms))
end

function BatchedAtomicGraph(batch_size::Int, laplacians, encoded_features)
    return batch_graph_data.(zip(Iterators.partition(laplacians, batch_size),
                                 Iterators.partition(encoded_features, batch_size)))
end


text:
include("../src/mcmc.jl")
include("../src/truth.jl")
using Plots
using Plots.PlotMeasures
using LinearAlgebra #required to load MCMC object
using .MCMC
using JLD2

#backend
#pyplot()
"""
Plotting script of the trace for the Markov Chain Monte Carlo (MCMC),

Requires:
--------
- MCMC object

produces one plot over the MCMC iterations per parameter in the MCMC learning.

"""

function main()
    #create an Array of the MCMC stored
    homedir=split(pwd(),"/test")[1]*"/"
    outdir=homedir*"output/"
    datdir=outdir*"exp_name"
    mcmcfile = datdir*"MCMC.jld2"

    @load mcmcfile mcmc
    posterior=mcmc.posterior #without burnin

    for i = 1:size(posterior)[2]
        plot(posterior[:,i],
             legend=false,
             dpi=300,
             left_margin=50px,
             bottom_margin=50px)
        xlabel!("iteration")
        ylabel!("param_"*string(i) )
        savefig(outdir*"mcmctrace_param_"*string(i)*".pdf")
    end

end

main()


text:
module ChangeBasisTests

using InplaceArrays.Fields
using InplaceArrays.Polynomials

xi = Point(0.5,0.5)
np = 5
x = fill(xi,np)

nodes = Point{2,Int}[(0,0),(1,0),(0,1),(1,1)]

order = 1
V = Float64
G = gradient_type(V,xi)
H = gradient_type(G,xi)
f = MonomialBasis{2}(V,order)

change = inv(evaluate(f,nodes))

g = change_basis(f,change)

gx = V[1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0]
∇gx = G[
  (-1.0, -1.0) (1.0, 0.0) (0.0, 1.0) (0.0, 0.0);
  (-1.0, 0.0) (1.0, -1.0) (0.0, 0.0) (0.0, 1.0);
  (0.0, -1.0) (0.0, 0.0) (-1.0, 1.0) (1.0, 0.0);
  (0.0, 0.0) (0.0, -1.0) (-1.0, 0.0) (1.0, 1.0)]
∇∇gx = H[
  (0.0, 1.0, 1.0, 0.0) (0.0, -1.0, -1.0, 0.0) (0.0, -1.0, -1.0, 0.0) (0.0, 1.0, 1.0, 0.0);
  (0.0, 1.0, 1.0, 0.0) (0.0, -1.0, -1.0, 0.0) (0.0, -1.0, -1.0, 0.0) (0.0, 1.0, 1.0, 0.0);
  (0.0, 1.0, 1.0, 0.0) (0.0, -1.0, -1.0, 0.0) (0.0, -1.0, -1.0, 0.0) (0.0, 1.0, 1.0, 0.0);
  (0.0, 1.0, 1.0, 0.0) (0.0, -1.0, -1.0, 0.0) (0.0, -1.0, -1.0, 0.0) (0.0, 1.0, 1.0, 0.0)]

test_field(g,nodes,gx,grad=∇gx,hessian=∇∇gx)

end # module


text:
#if lowercase(get(ENV, "CI", "false")) == "true"
#    include("install_pycall.jl")
#end
using SafeTestsets, Pkg, Distributions

# Now we need to grab all the subpackages for testing
const GROUP = get(ENV, "GROUP", "ALL")
#using PyCall

# add develop to subpkg
function dev_subpkg(subpkg)
    subpkg_path = joinpath(dirname(@__DIR__), "lib", subpkg)
    Pkg.develop(PackageSpec(path=subpkg_path))
end

# activate subpkg and instantiate it
function activate_subpkg_env(subpkg)
    subpkg_path = joinpath(dirname(@__DIR__), "lib", subpkg)
    Pkg.activate(subpkg_path)
    Pkg.develop(PackageSpec(path=subpkg_path))
    Pkg.instantiate()
end

# Now split depending on what kind of test we are doing
if GROUP == "ALL" || GROUP == "Core"
    dev_subpkg("ComradeGalactic")
    @safetestset "CORE Comrade.jl" begin
        include(joinpath(@__DIR__, "Core/core.jl"))
    end
else
    dev_subpkg(GROUP)
    subpkg_path = joinpath(dirname(@__DIR__), "lib", GROUP)
    Pkg.test(PackageSpec(name=GROUP, path=subpkg_path))
end


text:
include("intervalTreeNode.jl")
include("intervalSearchTree.jl")
include("treeTrimming.jl")


text:
# This file is a part of SimilaritySearch.jl

"""
    find_neighborhood(index::SearchGraph{T}, item, res, vstate)

Searches for `item` neighborhood in the index, i.e., if `item` were in the index whose items should be its neighbors (intenal function).
`res` is always reused since `reduce` creates a new KnnResult from it (a copy if `reduce` in its simpler terms)
"""
function find_neighborhood(index::SearchGraph, item, res, vstate)
    n = length(index)

    if n > 0
        empty!(res, index.neighborhood.ksearch)
        reduce(index.neighborhood.reduce, search(index, item, res; vstate), index)
    else
        KnnResult(index.neighborhood.ksearch)
    end
end

"""
    push_neighborhood!(index::SearchGraph, item, neighbors::KnnResult; apply_callbacks=true)

Inserts the object `item` into the index, i.e., creates an edge from items listed in L and the
vertex created for ìtem` (internal function)
"""
function push_neighborhood!(index::SearchGraph, item, neighbors::KnnResult; apply_callbacks=true)
    push!(index.db, item)
    push!(index.links, neighbors)
    push!(index.locks, Threads.SpinLock())
    n = length(index)
    k = index.neighborhood.k

    @inbounds for (id, dist) in neighbors
        vertex = index.links[id]
        vertex.k = max(maxlength(vertex), k) # adjusting maximum size to the current allowed neighborhood size
        push!(vertex, n, dist)
    end

    apply_callbacks && callbacks(index)

    if index.verbose && length(index) % 100_000 == 0
        println(stderr, "added n=$(length(index)), neighborhood=$(length(neighbors)), $(string(index.search_algo)), $(Dates.now())")
    end
end

const GlobalSatKnnResult = [KnnResult(1)]
function __init__neighborhood()
    for i in 2:Threads.nthreads()
        push!(GlobalSatKnnResult, KnnResult(1))
    end
end

"""
    SatNeighborhood()

New items are connected with a small set of items computed with a SAT like scheme (**cite**).
It starts with `k` near items that are reduced to a small neighborhood due to the SAT partitioning stage.
"""
struct SatNeighborhood <: NeighborhoodReduction end

Base.copy(::SatNeighborhood) = SatNeighborhood()

"""
    reduce(sat::SatNeighborhood, res::KnnResult, index::SearchGraph)

Reduces `res` using the SAT strategy.
"""
function Base.reduce(sat::SatNeighborhood, res::KnnResult, index::SearchGraph)
    near = GlobalSatKnnResult[Threads.threadid()]
    N = KnnResult(ceil(Int, maxlength(res)/2))

    #@inbounds for i in lastindex(res):-1:firstindex(res)  # DistSat => works a little better but also produces a bit larger neighborhoods
    #id, dist = res[i]
    @inbounds for (id, dist) in res
        pobj = index[id]
        empty!(near)
        push!(near, 0, dist)
        for nearID in keys(N)
            d = evaluate(index.dist, index[nearID], pobj)
            push!(near, nearID, d)
        end

        argmin(near) == 0 && push!(N, id, dist)
    end

    N
end

"""
    struct IdentityNeighborhood

It does not modifies the given neighborhood
"""
struct IdentityNeighborhood <: NeighborhoodReduction end

Base.copy(::IdentityNeighborhood) = IdentityNeighborhood()
Base.reduce(sat::IdentityNeighborhood, res::KnnResult, index::SearchGraph) = copy(res)


text:
### Low-level utilities ###

_indexoffset(r::AbstractRange) = first(r) - 1
_indexoffset(i::Integer) = 0
_indexlength(r::AbstractRange) = length(r)
_indexlength(i::Integer) = Int(i)
_indexlength(i::Colon) = Colon()

# utility methods used in reshape
# we don't use _indexlength in this to avoid converting the arguments to Int
_checksize(ind::Integer, s) = ind == s
_checksize(ind::AbstractUnitRange, s) = length(ind) == s

_toaxis(i::Integer) = Base.OneTo(i)
_toaxis(i) = i

_strip_IdOffsetRange(r::IdOffsetRange) = parent(r)
_strip_IdOffsetRange(r) = r

_offset(axparent::AbstractUnitRange, ax::AbstractUnitRange) = first(ax) - first(axparent)
_offset(axparent::AbstractUnitRange, ::Union{Integer, Colon}) = 1 - first(axparent)

_offsets(A::AbstractArray) = map(ax -> first(ax) - 1, axes(A))
_offsets(A::AbstractArray, B::AbstractArray) = map(_offset, axes(B), axes(A))

"""
    OffsetArrays.AxisConversionStyle(typeof(indices))

`AxisConversionStyle` declares if `indices` should be converted to a single `AbstractUnitRange{Int}`
or to a `Tuple{Vararg{AbstractUnitRange{Int}}}` while flattening custom types into indices.
This method is called after `to_indices(A::Array, axes(A), indices)` to provide
further information in case `to_indices` does not return a `Tuple` of `AbstractUnitRange{Int}`.

Custom index types should extend `AxisConversionStyle` and return either `OffsetArray.SingleRange()`,
which is the default, or `OffsetArray.TupleOfRanges()`. In the former case, the type `T` should
define `Base.convert(::Type{AbstractUnitRange{Int}}, ::T)`, whereas in the latter it should define
`Base.convert(::Type{Tuple{Vararg{AbstractUnitRange{Int}}}}, ::T)`.

An example of the latter is `CartesianIndices`, which is converted to a `Tuple` of
`AbstractUnitRange{Int}` while flattening the indices.

# Example
```jldoctest; setup=:(using OffsetArrays)
julia> struct NTupleOfUnitRanges{N}
           x ::NTuple{N, UnitRange{Int}}
       end

julia> Base.to_indices(A, inds, t::Tuple{NTupleOfUnitRanges{N}}) where {N} = t;

julia> OffsetArrays.AxisConversionStyle(::Type{NTupleOfUnitRanges{N}}) where {N} = OffsetArrays.TupleOfRanges();

julia> Base.convert(::Type{Tuple{Vararg{AbstractUnitRange{Int}}}}, t::NTupleOfUnitRanges) = t.x;

julia> a = zeros(3, 3);

julia> inds = NTupleOfUnitRanges((3:5, 2:4));

julia> oa = OffsetArray(a, inds);

julia> axes(oa, 1) == 3:5
true

julia> axes(oa, 2) == 2:4
true
```
"""
abstract type AxisConversionStyle end
struct SingleRange <: AxisConversionStyle end
struct TupleOfRanges <: AxisConversionStyle end

AxisConversionStyle(::Type) = SingleRange()
AxisConversionStyle(::Type{<:CartesianIndices}) = TupleOfRanges()

_convertTupleAbstractUnitRange(x) = _convertTupleAbstractUnitRange(AxisConversionStyle(typeof(x)), x)
_convertTupleAbstractUnitRange(::SingleRange, x) = (convert(AbstractUnitRange{Int}, x),)
_convertTupleAbstractUnitRange(::TupleOfRanges, x) = convert(Tuple{Vararg{AbstractUnitRange{Int}}}, x)

_toAbstractUnitRanges(t::Tuple) = (_convertTupleAbstractUnitRange(first(t))..., _toAbstractUnitRanges(tail(t))...)
_toAbstractUnitRanges(::Tuple{}) = ()

# ensure that the indices are consistent in the constructor
_checkindices(A::AbstractArray, indices, label) = _checkindices(ndims(A), indices, label)
function _checkindices(N::Integer, indices, label)
    throw_argumenterror(N, indices, label) = throw(ArgumentError(label*" $indices are not compatible with a $(N)D array"))
    N == length(indices) || throw_argumenterror(N, indices, label)
end

@inline _indexedby(r::AbstractVector, ax::Tuple{Any}) = _indexedby(r, ax[1])
@inline _indexedby(r::AbstractUnitRange{<:Integer}, ::Base.OneTo) = no_offset_view(r)
@inline _indexedby(r::AbstractUnitRange{Bool}, ::Base.OneTo) = no_offset_view(r)
@inline _indexedby(r::AbstractVector, ::Base.OneTo) = no_offset_view(r)
@inline function _indexedby(r::AbstractUnitRange{<:Integer}, ax::AbstractUnitRange)
	of = convert(eltype(r), first(ax) - 1)
	IdOffsetRange(_subtractoffset(r, of), of)
end
@inline _indexedby(r::AbstractUnitRange{Bool}, ax::AbstractUnitRange) = OffsetArray(r, ax)
@inline _indexedby(r::AbstractVector, ax::AbstractUnitRange) = OffsetArray(r, ax)

# These functions are equivalent to the broadcasted operation r .- of
# However these ensure that the result is an AbstractRange even if a specific
# broadcasting behavior is not defined for a custom type
@inline _subtractoffset(r::AbstractUnitRange, of) = UnitRange(first(r) - of, last(r) - of)
@inline _subtractoffset(r::AbstractRange, of) = range(first(r) - of, stop = last(r) - of, step = step(r))

# similar to _subtractoffset, except these evaluate r .+ of
@inline _addoffset(r::AbstractUnitRange, of) = UnitRange(first(r) + of, last(r) + of)
@inline _addoffset(r::AbstractRange, of) = range(first(r) + of, stop = last(r) + of, step = step(r))

if VERSION <= v"1.7.0-DEV.1039"
    _contiguousindexingtype(r::AbstractUnitRange{<:Integer}) = UnitRange{Int}(r)
else
    _contiguousindexingtype(r::AbstractUnitRange{<:Integer}) = r
end

_of_eltype(::Type{T}, M::AbstractArray{T}) where {T} = M
_of_eltype(T, M::AbstractArray) = map(T, M)

# filter the arguments to reshape to check if there are any ranges
# If not, we may pop the parent array
_filterreshapeinds(t::Tuple{AbstractUnitRange, Vararg{Any}}) = t
_filterreshapeinds(t::Tuple) = _filterreshapeinds(tail(t))
_filterreshapeinds(t::Tuple{}) = t
_popreshape(A::AbstractArray, ax::Tuple{Vararg{Base.OneTo}}, inds::Tuple{}) = no_offset_view(A)
_popreshape(A::AbstractArray, ax, inds) = A