text:
const _STOP_SIGNAL_KEY = :stop_dignal
_get_stop_signal(gl::GitLink) = get!(gl, _STOP_SIGNAL_KEY, false)
_set_stop_signal!(gl::GitLink, sig::Bool) = set!(gl, _STOP_SIGNAL_KEY, sig)

"""
    run_sync_loop(gl::GitLink; niters = typemax(Int), verbose = true)

Run sucesive `sync_link` operations over the `GilLink`.
It will stop if `niter` or` stop_time` is reached.
It is a lazy method, if no action is require no action will be made (use `force` to avoid it).
A server might call this method asynchronously.
"""
function run_sync_loop(gl::GitLink;
        niters = typemax(Int),
        tout = Inf,
        verbose = true
    )

    t0 = time()
    tout = _LOCK_FORCE_TIME

    for it in 1:niters

        ## ---------------------------------------------------
        # WAIT
        sleep(loop_frec(gl))

        ## ---------------------------------------------------
        # TIME OUT
        (time() - t0) > tout && return gl

        ## ---------------------------------------------------
        # STOP SIGNAL (Dev)
        _get_stop_signal(gl) && return gl

        ## ---------------------------------------------------
        # NEW ITER
        verbose && println("-"^60)
        verbose && @info("Loop iter", it, pid = getpid(), loop_frec = loop_frec(gl))

        ## ---------------------------------------------------
        # FORCE
        force = _do_ping(gl)

        ## ---------------------------------------------------
        # SYNC LINK
        sync_link(gl; verbose, force, tout)

    end

    return gl

end


text:
function deserialize_audio_summary(summary)
    audio = summary.audio

    if audio.content_type == "audio/wav"
        value = load(_format_stream(format"WAV", PipeBuffer(audio.encoded_audio_string)))
    else
        value = nothing
        @warn "unknown audio format $(audio.content_type)"
    end

    return value
end

ext: jl


text:
#=
test_stack_queue:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-19
=#

module TestStackQueue
    # packages
    using Test

    # external modules
    include("my_queue.jl")
    include("my_stack.jl")

    # methods
    function test()
        @testset "StackQueue" begin
            @testset "Stack" begin
                @test_nowarn MyStack.main()
            end
            @testset "Queue" begin
                @test_nowarn MyQueue.main()
            end
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .TestStackQueue
    TestStackQueue.test()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
#=
ProximalMethods.jl

    Provides proximal operator evaluation routines and proximal optimization
	algorithms, such as (accelerated) proximal gradient methods and alternating
	direction method of multipliers (ADMM), for non-smooth/non-differentiable
	objective functions.

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

@date: 2022/01/12
=#

module ProximalMethods

using   LinearAlgebra,
        Optim

export
# Proximal Operators
	soft_thresh,
	block_soft_thresh, block_soft_thresh!,
	shrinkage, shrinkage!,
    smooth!,
# Proximal Gradient Methods
	prox_grad, prox_grad!,
# ADMM
	admm, admm!

# Types
abstract type AbstractAccelScheme end

# Structs
# Line Search
mutable struct BackTrack{Tf, Gf}
	λ::Tf		# stepsize
	λ_prev::Tf 	# previous stepsize
	β::Gf		# contraction factor
end

# Acceleration
# No acceleration
struct NoAccel{Tf} <: AbstractAccelScheme
	ω::Tf	# extrapolation paramater
end
# Simple extrapolation
mutable struct Simple{Tf, Gi} <: AbstractAccelScheme
	ω::Tf	# extrapolation paramater
	k::Gi	# iteration counter
end
# Nesterov momentum
Base.@kwdef mutable struct Nesterov{Tf, G} <: AbstractAccelScheme
	ω::Tf		# extrapolation parameter
	θ::Tf		# momentum parameter
	m::Tf = .0	# convexity parameter
	ls::G		# line search parameters
end

# Proximal Gradient Methods
struct ProxGradState{Tv}
	x::Tv		# current state
	x_prev::Tv	# previous state
	y::Tv		# extrapolated state
	Δ::Tv		# change in state
	∇f::Tv		# gradient of f
end

# ADMM
struct ADMMState{Tv}
	x::Tv		# current state
	z::Tv		# current state
	u::Tv		# current
	z_prev::Tv	# previous state
	r::Tv		# primal residuals
	s::Tv		# dual residuals
end

# Include programs
include("proximal.jl")
include("acceleration.jl")
include("gradient.jl")
include("admm.jl")

end


text:
module test_julia_stream

using Test

path = normpath(@__DIR__, "stream.jl")

io = open(path)
@test position(io) == 0

seek(io, 1)
@test position(io) == 1

seekend(io)

seekstart(io)
@test position(io) == 0

close(io)

end # module test_julia_stream


text:
# Metida

#Simple dataset
df0         = CSV.File(path*"/csv/df0.csv"; types = [String, String, String, String, Float64, Float64]) |> DataFrame

df0m         = CSV.File(path*"/csv/df0miss.csv"; types = [String, String, String, String, Float64, Float64]) |> DataFrame

df1         = CSV.File(path*"/csv/df1.csv"; types = [String, String, String, String, Float64, Float64]) |> DataFrame

ftdf         = CSV.File(path*"/csv/1fptime.csv"; types = [String, String, Float64, Float64]) |> DataFrame

ftdf2        = CSV.File(path*"/csv/1freparma.csv"; types = [String, String, Float64, Float64]) |> DataFrame

ftdf3        = CSV.File(path*"/csv/2f2rand.csv"; types =
[String,  Float64, Float64, String, String, String, String, String]) |> DataFrame


text:
@testset "Test interpolation of outgoing_translation_matrix" begin

    particle_medium = Acoustic(3; ρ=10.0, c=80.0);
    medium = Acoustic(3; ρ=1.0, c=1.0);

    R = 3.0
    r = 1.0

    separation_ratio = 1.02
    ka = 0.1
    k = ka / r
    vol_fraction = 0.3

    basis_order = 2

    ω = k * real(medium.c)

    s1 = Specie(
        particle_medium, Sphere(r);
        number_density = vol_fraction / volume(Sphere(r)),
        exclusion_distance = separation_ratio
    );
    species = [s1]
    material = Material(Sphere(R),species);

    a12 = 2.0 * minimum(outer_radius(s) * s.exclusion_distance for s in material.species)

    tol = 1e-2

    Uinter = outgoing_translation_matrix(ω, medium, material;  basis_order = basis_order, tol = tol)

    rs = LinRange(a12,2R-a12,13);
    φs = LinRange(-pi,pi,13);
    θs = LinRange(0.0,pi,13);

    r2x = radial_to_cartesian_coordinates

    error = [
        begin
            U1 = outgoing_translation_matrix(medium, basis_order, ω, r2x([r,θ,φ]))
            norm(Uinter(r2x([r,θ,φ])) - U1) / norm(U1)
        end
    for r in rs, θ in θs, φ in φs];

    @test maximum(error) < 5 * tol

    @test mean(error) < tol

end

@testset "Legendre polynomial approximation" begin

    P = Legendre{Float64}()

    polynomial_order = 25
    mesh_points = 3*(polynomial_order + 1)
    # mesh_points = (polynomial_order + 1)^2

    rs = LinRange(0,5,mesh_points)
    r2s = LinRange(0,5,4*mesh_points)

    r1_max = maximum(rs);
    rbars = 2.0 .* rs ./ r1_max .- (1.0);
    r2bars = 2.0 .* r2s ./ r1_max .- (1.0);

    data = cos.(rs) .* (1.0 .+ 4 .* exp.(-10 .* rs.^4))
    data2 = cos.(r2s) .* (1.0 .+ 4 .* exp.(-10 .* r2s.^4))
    # plot(rs,data)
    # plot!(r2s,data2)

    Pmat = P[rbars, 1:(polynomial_order + 1)];
    P2mat = P[r2bars, 1:(polynomial_order + 1)];

    # pls_arr * transpose(Pmat) ~ data
    # Pmat * transpose(pls_arr) ~ transpose(data)
    pls_arr = Pmat \ data

    ws = (0:polynomial_order) .+ 1/2
    σs = integration_scheme(rbars; scheme = :trapezoidal)
    tls_arr = ws .* (transpose(Pmat) * (σs .* data))

    λ = 1e-6
    cls_arr = inv(transpose(Pmat) * Pmat + λ*I) * transpose(Pmat) * data

    # norm(Pmat * pls_arr - data) / norm(data)
    # norm(Pmat * tls_arr - data) / norm(data)
    # norm(Pmat * cls_arr - data) / norm(data)
    @test norm(P2mat * pls_arr - data2) / norm(data2) < 5e-3
    @test norm(P2mat * tls_arr - data2) / norm(data2) < 0.6
    @test norm(P2mat * cls_arr - data2) / norm(data2) < 5e-3

# using Plots
#
#     scatter(rs,data)
#     plot!(rs,Pmat * pls_arr, linestyle = :dash)
#     plot!(r2s,P2mat * pls_arr, linestyle = :dash)
#     plot!(r2s,P2mat * tls_arr, linestyle = :dash)
#     plot!(r2s,P2mat * cls_arr, linestyle = :dash)

end

max_stars_repo_name: arturgower/EffectiveWaves


text:
"""
Bridges `CP.SlidingSum` to linear constraints.
"""
struct SlidingSum2LPBridge{T} <: MOIBC.AbstractBridge
    cons::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.Interval{T}}}
end

function MOIBC.bridge_constraint(
    ::Type{SlidingSum2LPBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.SlidingSum{T},
) where {T}
    return MOIBC.bridge_constraint(
        SlidingSum2LPBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{SlidingSum2LPBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.SlidingSum{T},
) where {T}
    f_scalars = MOIU.scalarize(f)
    cons = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.Interval{T}}[
        MOI.add_constraint(
            model,
            sum(f_scalars[i:(i + s.length - 1)]),
            MOI.Interval(s.low, s.high)
        )
        for i in 1:(s.dimension - s.length)
    ]

    return SlidingSum2LPBridge(cons)
end

function MOI.supports_constraint(
    ::Type{SlidingSum2LPBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.SlidingSum{T}},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{SlidingSum2LPBridge{T}}) where {T}
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(::Type{SlidingSum2LPBridge{T}}) where {T}
    return [
        (MOI.ScalarAffineFunction{T}, MOI.Interval{T}),
    ]
end

function MOI.get(
    b::SlidingSum2LPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.Interval{T},
    },
) where {T}
    return length(b.cons)
end

function MOI.get(
    b::SlidingSum2LPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.Interval{T},
    },
) where {T}
    return copy(b.cons)
end


text:
# include("../preamble.jl")
# include((@__DIR__)*"/../../preamble.jl")

T = [-1.1 1.1; 1.0 -1.0]
c = [1.0; -1.0]
b_absorbing = 10.0

T_absorbing = [T [0.0 0.0; 0.0 0.0]; [0.01 0.0 -0.01 0.0; 0.0 0.01 0.0 -0.01]]
c_absorbing = [c; 0.0; 0.0]

P_upr_absorbing = [0.0 0.0 0.0 1.0]
P_lwr_absorbing = [0.0 0.0 1.0 0.0] # both absorbing boundaries

model = BoundedFluidQueue(T_absorbing,c_absorbing,P_lwr_absorbing,P_upr_absorbing,b_absorbing)


text:
module WildBootTests
export BootTestResult, wildboottest, AuxWtType, PType, MAdjType, DistStatType,
       teststat, stattype, p, padj, reps, repsfeas, nbootclust, dof, dof_r, plotpoints, peak, CI, dist, statnumer, statvar, auxweights

using LinearAlgebra, Random, Distributions, SortingAlgorithms, LoopVectorization

include("StrBootTest.jl")
include("utilities.jl")
include("estimators.jl")
include("init.jl")
include("WRE.jl")
include("nonWRE.jl")
include("plot-CI.jl")
include("interface.jl")

# top-level computation routine for OLS/ARubin (and score BS on IV/2SLS); split off to reduce latency when just doing WRE
function boottestOLSARubin!(o::StrBootTest{T}) where T
  if !o.initialized
		Init!(o)
	elseif !o.null
		NoNullUpdate!(o)
		return
  end

  if o.Nw > 1  # if more than one weight group to save memory, make on every call to boottest(), not just once in Init!()
		Random.seed!(o.rng,o.seed)
		MakeWildWeights!(o, last(o.WeightGrp[1])-1, first=true)
	end

  MakeInterpolables!(o)  # make stuff that depends linearly on r, possibly by interpolating, for first weight group

  for w ∈ 1:o.Nw  # do group 1 first because it includes col 1, which is all that might need updating in constructing CI in WCU
		w > 1 && MakeWildWeights!(o, length(o.WeightGrp[w]), first=false)

		MakeNonWREStats!(o, w)

		!o.bootstrapt && UpdateBootstrapcDenom!(o, w)
  end

  o.BFeas = isnan(o.dist[1]) ? 0 : sum(.!(isnan.(o.dist) .| isinf.(o.dist))) - 1
  o.distCDR = zeros(T,0,0)
  nothing
end

# top-level computation routine for non-ARubin WRE; split off to reduce latency when just doing other tests
function boottestWRE!(o::StrBootTest{T}) where T
  if !o.initialized
		Init!(o)
	elseif !o.null
		NoNullUpdate!(o)
		return
  end

  if o.Nw > 1  # if more than one weight group to save memory, make on every call to boottest(), not just once in Init!()
		Random.seed!(o.rng,o.seed)
		MakeWildWeights!(o, last(o.WeightGrp[1])-1, first=true)
	end

  PrepWRE!(o)

  for w ∈ 1:o.Nw  # do group 1 first because it includes col 1, which is all that might need updating in constructing CI in WCU
		w > 1 && MakeWildWeights!(o, length(o.WeightGrp[w]), first=false)

		MakeWREStats!(o, w)

		!o.bootstrapt && UpdateBootstrapcDenom!(o, w)
  end

  o.BFeas = isnan(o.dist[1]) ? 0 : sum(.!(isnan.(o.dist) .| isinf.(o.dist))) - 1
  o.distCDR = zeros(T,0,0)
  nothing
end

# if not imposing null and we have returned to boottest!(), then dof=1 or 2; we're plotting or finding CI, and only test stat, not distribution, changes with r
function NoNullUpdate!(o::StrBootTest{T} where T)
  if o.WREnonARubin
		o.numer[:,1] = o.R * o.DGP.Rpar * o.β̂s[1] - o.r
  elseif o.ARubin
		EstimateARubin!(o.DGP, o, o.r)
		o.numer[:,1] = o.v_sd * @view o.DGP.β̈[o.kX₁+1:end,:]  # coefficients on excluded instruments in ARubin OLS
  else
		o.numer[:,1] = o.v_sd * (o.R * (o.ML ? o.β̈ : o.M.β̈) - o.r)  # Analytical Wald numerator; if imposing null then numer[:,1] already equals this. If not, then it's 0 before this
  end
  o.dist[1] = isone(o.dof) ? o.numer[1] / sqrt(o.statDenom[1]) : o.numer[:,1]'invsym(o.statDenom)*o.numer[:,1]
	nothing
end

# compute bootstrap-c denominator from all bootstrap numerators
function UpdateBootstrapcDenom!(o::StrBootTest{T} where T, w::Integer)
  if isone(w)
		tmp = o.numer[:,1]
		o.statDenom = o.numer * o.numer' - tmp * tmp'
		o.numersum = rowsum(o.numer) - tmp
  else
		statDenom .+= o.numer * o.numer'
		numersum .+= rowsum(o.numer)
  end
  if w == o.Nw  # last weight group?
		o.statDenom .= (o.statDenom - o.numersum * o.numersum' / o.B) / o.B
		if o.sqrt
			o.dist .= o.numer ./ sqrtNaN.(o.statDenom)
		else
			negcolquadform!(o.dist, -invsym(o.statDenom), o.numer)  # to reduce latency by minimizing @tturbo instances, work with negative of colquadform in order to fuse code with colquadformminus!
		end
  end
	nothing
end

# include("precompile_WildBootTests.jl")  # source: https://timholy.github.io/SnoopCompile.jl/stable/snoopi_deep_parcel/#SnoopCompile.write
# _precompile_()

end

using StableRNGs, Random, BenchmarkTools
N=1_000_00; G=40; k=2; l=4
Random.seed!(1231)
β=rand(); γ=rand(k); Π=rand(l)
W = rand(N,l)
u₂ = randn(N)
y₂ = W * Π + u₂
u₁ = u₂ + randn(N)
Z = rand(N,k)
y₁ = y₂ * β + Z * γ + u₁
ID = floor.(Int8, collect(0:N-1) / (N/G))
R = [zeros(1,k) 1]; r = [0]
WildBootTests.wildboottest(Float32, R,.46; R1 = [1 zeros(1,k)], r1=[0], resp=y₁, predexog=Z, predendog=y₂, inst=W, clustid=ID, reps=999, issorted=true, rng=StableRNG(1231), LIML=true)

N=1_000_000; G=40; k=12; l=40
Random.seed!(1231)
β=rand(); γ=rand(k); Π=rand(l)
W = rand(N,l)
u₂ = randn(N)
y₂ = W * Π + u₂
u₁ = u₂ + randn(N)
Z = rand(N,k)
y₁ = y₂ * β + Z * γ + u₁
ID = floor.(Int8, collect(0:N-1) / (N/G))
R = [zeros(1,k) 1]; r = [0]
@btime WildBootTests.wildboottest(Float32, R,.4; resp=y₁, predexog=Z, predendog=y₂, inst=W, clustid=ID, reps=9999, issorted=true, getCI=false);

@btime WildBootTests.wildboottest(Float32, R,.4; resp=y₁, predexog=Z, predendog=y₂, inst=W, clustid=ID, reps=9999, issorted=true, getCI=true);


text:
"""
            Functions to write snapshots.

    Author: Ludwig Böss
    Contact: lboess@usm.lmu.de
    Created: 2018-12-12

"""

function write_format2_block_header(f::IOStream, blockname::String)

    # fill blockname to 4 characters
    if length(blockname) < 4
        padding = 4 - length(blockname)
        for i=1:padding
            blockname = blockname * " "
        end
    end

    # transform the name into an array of charcters.
    # Since Julia uses 2 byte chars and C uses 1 byte I have to take
    # a different approach and use Int8s.
    write_name =  Int8.(collect(blockname))

    @info "Writing block: $blockname"

    # write blocksite (8 bytes)
    write(f, Int32(8))

    # write block name
    write(f, write_name)

    return f
end

"""
    write_block(f::IOStream, data,
                blockname::String="";
                snap_format::Int64=2)

Write `data` to a block to an opened file `f`.
"""
function write_block(f::IOStream, data,
                     blockname::String="";
                     snap_format::Int64=2)

    # turn array of arrays d into one large array
    #data = reduce(vcat, d)     # OLD! Keep for later, maybe

    # get total number of particles to write
    N = length(data[:,1])

    # write blocksize
    dtype = typeof(data[1,1])
    dims = length(data[1,:])
    blocksize = Int32(N * sizeof(dtype) * dims)

    if snap_format == 2

        # break if blockname not specified
        if blockname == ""
            error("Please specify blockname!")
        end

        # write header
        f = write_format2_block_header(f, blockname)

        # write end if name header, namely:
        # size to skip from directly after the block name and size of
        # the format 2 block header
        write(f, Int32(blocksize + Int32(8)))
        write(f, Int32(8))
    end

    # write blocksize
    write(f, blocksize)

    # write the block. Since Julia stores the arrays differently in memory
    # they have to be transposed before the can be written.
    write(f, copy(transpose(data)))

    @info "Writing block done."

    write(f, blocksize)
end

"""
    write_header(f::IOStream, h::Header; snap_format::Int64=2)

Writes the header block to an opened file `f`.
"""
function write_header(f::IOStream, h::Header; snap_format::Int64=2)

    if snap_format == 2
        f = write_format2_block_header(f, "HEAD")
        write(f, Int32(264))
        write(f, Int32(8))
    end

    # write blocksize
    write(f, Int32(256))

    # write header to file
    for fields in fieldnames(Header)

        write(f, getfield(h, fields))

    end

    # write blocksize
    write(f, Int32(256))
end

# Format 2 Blocknames
# "POS "
# "VEL "
# "ID  "
# "MASS"
# "U   "
# "RHO "
# "NE  "
# "NH  "
# "HII "
# "HeI "
# "HeII"
# "He3 "
# "H2I "
# "H2II"
# "HM  "
# "HD  "
# "DI  "
# "DII "
# "HeHp"
# "HSML"
# "VALP"
# "SFR "
# "AGE "
# "DETI"
# "HSMS"
# "ACRS"
# "Z   "
# "POT "
# "PHI "
# "PHID"
# "ACCE"
# "ENDT"
# "STRD"
# "STRO"
# "STRB"
# "SHCO"
# "TSTP"
# "BFLD"
# "BFSM"
# "DBDT"
# "VBLK"
# "VRMS"
# "VTAN"
# "VRAD"
# "TNGB"
# "DPP "
# "VDIC"
# "VROT"
# "VORT"
# "DIVB"
# "ACVC"
# "AMDC"
# "VTRB"
# "LTRB"
# "ADYN"
# "EDYN"
# "PHI "
# "XPHI"
# "GPHI"
# "ROTB"
# "SRTB"
# "EULA"
# "EULB"
# "COOR"
# "CONR"
# "DENN"
# "EGYP"
# "EGYC"
# "CRC0"
# "CRQ0"
# "CRP0"
# "CRE0"
# "CRn0"
# "CRco"
# "CRdi"
# "BHMA"
# "ACRB"
# "BHMD"
# "BHPC"
# "BHMB"
# "BHMI"
# "BHMR"
# "MACH"
# "SHSP"
# "SHRH"
# "SHPR"
# "SHVU"
# "SHCP"
# "SHNR"
# "DTEG"
# "PSCS"
# "PSDE"
# "PSEN"
# "PSXC"
# "DJMP"
# "EJMP"
# "CRDE"
# "TIPS"
# "DIPS"
# "CACO"
# "FLDE"
# "STDE"
# "SOMA"
# "PSDE"
# "ANRA"
# "LACA"
# "SHOR"
# "INDE"
# "aTEMP"  # remove a!
# "XNUC"
# "P   "
# "RADG"
# "RADA"
# "ET  "
# "PTSU"
# "DMNB"
# "NTSC"
# "SHSM"
# "SRHO"
# "SVEL"
# "DMHS"
# "DMDE"
# "DMVD"
# "VDMH"
# "VDMD"
# "Zs  "
# "CII "
# "CIa "
# "Cagb"
# "ZAge"
# "ZAlv"
# "iM  "
# "CLDX"
# "HOTT"
# "aTEMP"  # remove a!
# "TRCK"
# "ZsMT"
# "ZSMT"
# "DABN"
# "DABG"
# "DABS"
# "DCFB"
# "DCFA"
# "DCFD"

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
using Test

using POMDPs
using Random

mutable struct A <: POMDP{Int,Bool,Bool} end
@testset "implement" begin

    @test_throws MethodError n_states(A())
    @test_throws MethodError stateindex(A(), 1)

    @test !@implemented discount(::A)
    POMDPs.discount(::A) = 0.95
    @test @implemented discount(::A)

    @test !@implemented reward(::A,::Int,::Bool,::Int)
    @test !@implemented reward(::A,::Int,::Bool)
    POMDPs.reward(::A,::Int,::Bool) = -1.0
    @test @implemented reward(::A,::Int,::Bool,::Int)
    @test @implemented reward(::A,::Int,::Bool)

    @test !@implemented observation(::A,::Int,::Bool,::Int)
    @test !@implemented observation(::A,::Bool,::Int)
    POMDPs.observation(::A,::Bool,::Int) = [true, false]
    @test @implemented observation(::A,::Int,::Bool,::Int)
    @test @implemented observation(::A,::Bool,::Int)
end

mutable struct D end
POMDPs.sampletype(::Type{D}) = Int
struct E end
@testset "sampletype" begin
    @test @implemented sampletype(::D)
    @test sampletype(D()) == Int

    @test_throws MethodError sampletype(E)
    @test_throws MethodError sampletype(E())
end

@testset "infer" begin
    include("test_inferrence.jl")
end

@testset "require" begin
    include("test_requirements.jl")
end

@testset "generative" begin
    include("test_generative.jl")
end

# include("test_tutorials.jl")

@testset "genback" begin
    include("test_generative_backedges.jl")
end

struct CI <: POMDP{Int,Int,Int} end
struct CV <: POMDP{Vector{Float64},Vector{Float64},Vector{Float64}} end

@testset "convert" begin
    @test convert_s(Vector{Float32}, 1, CI()) == Float32[1.0]
    @test convert_s(statetype(CI), Float32[1.0], CI()) == 1
    @test convert_s(statetype(CV), Float32[2.0,3.0], CV()) == [2.0, 3.0]
    @test convert_s(Vector{Float32}, [2.0, 3.0], CV()) == Float32[2.0, 3.0]

    @test convert_a(Vector{Float32}, 1, CI()) == Float32[1.0]
    @test convert_a(statetype(CI), Float32[1.0], CI()) == 1
    @test convert_a(statetype(CV), Float32[2.0,3.0], CV()) == [2.0, 3.0]
    @test convert_a(Vector{Float32}, [2.0, 3.0], CV()) == Float32[2.0, 3.0]

    @test convert_o(Vector{Float32}, 1, CI()) == Float32[1.0]
    @test convert_o(statetype(CI), Float32[1.0], CI()) == 1
    @test convert_o(statetype(CV), Float32[2.0,3.0], CV()) == [2.0, 3.0]
    @test convert_o(Vector{Float32}, [2.0, 3.0], CV()) == Float32[2.0, 3.0]
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# OffsetArrays (arrays with indexing that doesn't start at 1)

# This test file is designed to exercise support for generic indexing,
# even though offset arrays aren't implemented in Base.

# OffsetArrays v1.11.2
# No compat patch and docstrings
module OffsetArrays

using Base: tail, @propagate_inbounds
using Base: IdentityUnitRange

export OffsetArray, OffsetMatrix, OffsetVector

const IIUR = IdentityUnitRange{<:AbstractUnitRange{<:Integer}}

########################################################################################################
# axes.jl
########################################################################################################

struct IdOffsetRange{T<:Integer,I<:AbstractUnitRange{T}} <: AbstractUnitRange{T}
    parent::I
    offset::T

    function IdOffsetRange{T,I}(r::I, offset::T) where {T<:Integer,I<:AbstractUnitRange{T}}
        _bool_check(T, r, offset)
        new{T,I}(r, offset)
    end

    #= This method is necessary to avoid a StackOverflowError in IdOffsetRange{T,I}(r::IdOffsetRange, offset::Integer).
    The type signature in that method is more specific than IdOffsetRange{T,I}(r::I, offset::T),
    so it ends up calling itself if I <: IdOffsetRange.
    =#
    function IdOffsetRange{T,IdOffsetRange{T,I}}(r::IdOffsetRange{T,I}, offset::T) where {T<:Integer,I<:AbstractUnitRange{T}}
        _bool_check(T, r, offset)
        new{T,IdOffsetRange{T,I}}(r, offset)
    end
end

function _bool_check(::Type{Bool}, r, offset)
    # disallow the construction of IdOffsetRange{Bool, UnitRange{Bool}}(true:true, true)
    if offset && (first(r) || last(r))
        throw(ArgumentError("values = $r and offset = $offset can not produce a boolean range"))
    end
    return nothing
end
_bool_check(::Type, r, offset) = nothing

# Construction/coercion from arbitrary AbstractUnitRanges
function IdOffsetRange{T,I}(r::AbstractUnitRange, offset::Integer = 0) where {T<:Integer,I<:AbstractUnitRange{T}}
    rc, o = offset_coerce(I, r)
    return IdOffsetRange{T,I}(rc, convert(T, o+offset)::T)
end
function IdOffsetRange{T}(r::AbstractUnitRange, offset::Integer = 0) where T<:Integer
    rc = convert(AbstractUnitRange{T}, r)::AbstractUnitRange{T}
    return IdOffsetRange{T,typeof(rc)}(rc, convert(T, offset)::T)
end
IdOffsetRange(r::AbstractUnitRange{T}, offset::Integer = 0) where T<:Integer =
    IdOffsetRange{T,typeof(r)}(r, convert(T, offset)::T)

# Coercion from other IdOffsetRanges
IdOffsetRange{T,I}(r::IdOffsetRange{T,I}) where {T<:Integer,I<:AbstractUnitRange{T}} = r
function IdOffsetRange{T,I}(r::IdOffsetRange, offset::Integer = 0) where {T<:Integer,I<:AbstractUnitRange{T}}
    rc, offset_rc = offset_coerce(I, r.parent)
    return IdOffsetRange{T,I}(rc, convert(T, r.offset + offset + offset_rc)::T)
end
IdOffsetRange{T}(r::IdOffsetRange{T}) where {T<:Integer} = r
function IdOffsetRange{T}(r::IdOffsetRange, offset::Integer = 0) where T<:Integer
    return IdOffsetRange{T}(r.parent, r.offset + offset)
end
IdOffsetRange(r::IdOffsetRange) = r

# Constructor to make `show` round-trippable
function IdOffsetRange(; values::AbstractUnitRange{<:Integer}, indices::AbstractUnitRange{<:Integer})
    length(values) == length(indices) || throw(ArgumentError("values and indices must have the same length"))
    offset = first(indices) - 1
    return IdOffsetRange(values .- offset, offset)
end

# Conversions to an AbstractUnitRange{Int} (and to an OrdinalRange{Int,Int} on Julia v"1.6") are necessary
# to evaluate CartesianIndices for BigInt ranges, as their axes are also BigInt ranges
Base.AbstractUnitRange{T}(r::IdOffsetRange) where {T<:Integer} = IdOffsetRange{T}(r)

# TODO: uncomment these when Julia is ready
# # Conversion preserves both the values and the indices, throwing an InexactError if this
# # is not possible.
# Base.convert(::Type{IdOffsetRange{T,I}}, r::IdOffsetRange{T,I}) where {T<:Integer,I<:AbstractUnitRange{T}} = r
# Base.convert(::Type{IdOffsetRange{T,I}}, r::IdOffsetRange) where {T<:Integer,I<:AbstractUnitRange{T}} =
#     IdOffsetRange{T,I}(convert(I, r.parent), r.offset)
# Base.convert(::Type{IdOffsetRange{T,I}}, r::AbstractUnitRange) where {T<:Integer,I<:AbstractUnitRange{T}} =
#     IdOffsetRange{T,I}(convert(I, r), 0)

offset_coerce(::Type{Base.OneTo{T}}, r::Base.OneTo) where T<:Integer = convert(Base.OneTo{T}, r), 0
function offset_coerce(::Type{Base.OneTo{T}}, r::AbstractUnitRange) where T<:Integer
    o = first(r) - 1
    return Base.OneTo{T}(last(r) - o), o
end
# function offset_coerce(::Type{Base.OneTo{T}}, r::IdOffsetRange) where T<:Integer
#     rc, o = offset_coerce(Base.OneTo{T}, r.parent)

# Fallback, specialze this method if `convert(I, r)` doesn't do what you need
offset_coerce(::Type{I}, r::AbstractUnitRange) where I<:AbstractUnitRange =
    convert(I, r)::I, 0

@inline Base.parent(r::IdOffsetRange) = r.parent
@inline Base.axes(r::IdOffsetRange) = (Base.axes1(r),)
@inline Base.axes1(r::IdOffsetRange) = IdOffsetRange(Base.axes1(r.parent), r.offset)
@inline Base.unsafe_indices(r::IdOffsetRange) = (Base.axes1(r),)
@inline Base.length(r::IdOffsetRange) = length(r.parent)
@inline Base.isempty(r::IdOffsetRange) = isempty(r.parent)
Base.reduced_index(i::IdOffsetRange) = typeof(i)(first(i):first(i))
# Workaround for #92 on Julia < 1.4
Base.reduced_index(i::IdentityUnitRange{<:IdOffsetRange}) = typeof(i)(first(i):first(i))
for f in [:firstindex, :lastindex]
    @eval @inline Base.$f(r::IdOffsetRange) = $f(r.parent) + r.offset
end
for f in [:first, :last]
    # coerce the type to deal with values that get promoted on addition (eg. Bool)
    @eval @inline Base.$f(r::IdOffsetRange) = eltype(r)($f(r.parent) + r.offset)
end

# Iteration for an IdOffsetRange
@inline Base.iterate(r::IdOffsetRange, i...) = _iterate(r, i...)
# In general we iterate over the parent term by term and add the offset.
# This might have some performance degradation when coupled with bounds-checking
# See https://github.com/JuliaArrays/OffsetArrays.jl/issues/214
@inline function _iterate(r::IdOffsetRange, i...)
    ret = iterate(r.parent, i...)
    ret === nothing && return nothing
    return (eltype(r)(ret[1] + r.offset), ret[2])
end
# Base.OneTo(n) is known to be exactly equivalent to the range 1:n,
# and has no specialized iteration defined for it,
# so we may add the offset to the range directly and iterate over the result
# This gets around the performance issue described in issue #214
# We use the helper function _addoffset to evaluate the range instead of broadcasting
# just in case this makes it easy for the compiler.
@inline _iterate(r::IdOffsetRange{<:Integer, <:Base.OneTo}, i...) = iterate(_addoffset(r.parent, r.offset), i...)

@inline function Base.getindex(r::IdOffsetRange, i::Integer)
    i isa Bool && throw(ArgumentError("invalid index: $i of type Bool"))
    @boundscheck checkbounds(r, i)
    @inbounds eltype(r)(r.parent[i - r.offset] + r.offset)
end

# Logical indexing following https://github.com/JuliaLang/julia/pull/31829
#= Helper function to perform logical indxeing for boolean ranges
The code implemented is a branch-free version of the following:

    range(first(s) ? first(r) : last(r), length=Int(last(s)))

See https://github.com/JuliaArrays/OffsetArrays.jl/pull/224#discussion_r595635143

Logical indexing does not preserve indices, unlike other forms of vector indexing
=#
@inline function _getindex(r, s::AbstractUnitRange{Bool})
    range(first(r) * first(s) + last(r) * !first(s), length=Int(last(s)))
end
@inline function _getindex(r, s::StepRange{Bool})
    range(first(r) * first(s) + last(r) * !first(s), step = oneunit(step(s)), length=Int(last(s)))
end
@inline function _getindex(r, s::AbstractUnitRange)
    @inbounds rs = r.parent[_subtractoffset(s, r.offset)] .+ r.offset
    _indexedby(rs, axes(s))
end
@inline function _getindex(r, s::StepRange)
    rs = @inbounds r.parent[s .- r.offset] .+ r.offset
    _indexedby(rs, axes(s))
end

for T in [:AbstractUnitRange, :StepRange]
    @eval @inline function Base.getindex(r::IdOffsetRange, s::$T{<:Integer})
        @boundscheck checkbounds(r, s)
        return _getindex(r, s)
    end
end

# These methods are necessary to avoid ambiguity
for R in [:IIUR, :IdOffsetRange]
    @eval @inline function Base.getindex(r::IdOffsetRange, s::$R)
        @boundscheck checkbounds(r, s)
        return _getindex(r, s)
    end
end

# offset-preserve broadcasting
Broadcast.broadcasted(::Base.Broadcast.DefaultArrayStyle{1}, ::typeof(-), r::IdOffsetRange{T}, x::Integer) where T =
    IdOffsetRange{T}(r.parent .- x, r.offset)
Broadcast.broadcasted(::Base.Broadcast.DefaultArrayStyle{1}, ::typeof(+), r::IdOffsetRange{T}, x::Integer) where T =
    IdOffsetRange{T}(r.parent .+ x, r.offset)
Broadcast.broadcasted(::Base.Broadcast.DefaultArrayStyle{1}, ::typeof(+), x::Integer, r::IdOffsetRange{T}) where T =
    IdOffsetRange{T}(x .+ r.parent, r.offset)

Base.show(io::IO, r::IdOffsetRange) = print(io, IdOffsetRange, "(values=",first(r), ':', last(r),", indices=",first(eachindex(r)),':',last(eachindex(r)), ")")

# Optimizations
@inline Base.checkindex(::Type{Bool}, inds::IdOffsetRange, i::Real) = Base.checkindex(Bool, inds.parent, i - inds.offset)

########################################################################################################
# origin.jl
########################################################################################################

struct Origin{T<:Union{Tuple{Vararg{Int}}, Int}}
    index::T
end
Origin(I::Tuple{Vararg{Int}}) = Origin{typeof(I)}(I)
Origin(I::Tuple{Vararg{Number}}) = Origin(map(Int, I))
Origin(I::CartesianIndex) = Origin(Tuple(I))
Origin(I::Number...) = Origin(I)
# Origin(0) != Origin((0, )) but they work the same with broadcasting
Origin(n::Number) = Origin{Int}(Int(n))

Base.Broadcast.broadcastable(o::Origin) = Ref(o)

_showidx(index::Integer) = "(" * string(index) * ")"
_showidx(index::Tuple) = string(index)
Base.show(io::IO, o::Origin) = print(io, "Origin", _showidx(o.index))

########################################################################################################
# utils.jl
########################################################################################################

### 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))

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))

_contiguousindexingtype(r::AbstractUnitRange{<:Integer}) = r

_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

########################################################################################################
# OffsetArrays.jl
########################################################################################################

# Technically we know the length of CartesianIndices but we need to convert it first, so here we
# don't put it in OffsetAxisKnownLength.
const OffsetAxisKnownLength = Union{Integer,AbstractUnitRange}
const OffsetAxis = Union{OffsetAxisKnownLength,Colon}
const ArrayInitializer = Union{UndefInitializer,Missing,Nothing}

## OffsetArray
struct OffsetArray{T,N,AA<:AbstractArray} <: AbstractArray{T,N}
    parent::AA
    offsets::NTuple{N,Int}
    @inline function OffsetArray{T, N, AA}(parent::AA, offsets::NTuple{N, Int}; checkoverflow = true) where {T, N, AA<:AbstractArray{T,N}}
        # allocation of `map` on tuple is optimized away
        checkoverflow && map(overflow_check, axes(parent), offsets)
        new{T, N, AA}(parent, offsets)
    end
end

const OffsetVector{T,AA<:AbstractVector{T}} = OffsetArray{T,1,AA}

const OffsetMatrix{T,AA<:AbstractMatrix{T}} = OffsetArray{T,2,AA}

# checks if the offset may be added to the range without overflowing
function overflow_check(r::AbstractUnitRange, offset::Integer)
    Base.hastypemax(eltype(r)) || return nothing
    # This gives some performance boost https://github.com/JuliaLang/julia/issues/33273
    throw_upper_overflow_error(val) = throw(OverflowError("offset should be <= $(typemax(Int) - val) corresponding to the axis $r, received an offset $offset"))
    throw_lower_overflow_error(val) = throw(OverflowError("offset should be >= $(typemin(Int) - val) corresponding to the axis $r, received an offset $offset"))

    # With ranges in the picture, first(r) might not necessarily be < last(r)
    # we therefore use the min and max of first(r) and last(r) to check for overflow
    firstlast_min, firstlast_max = minmax(first(r), last(r))

    if offset > 0 && firstlast_max > typemax(Int) - offset
        throw_upper_overflow_error(firstlast_max)
    elseif offset < 0 && firstlast_min < typemin(Int) - offset
        throw_lower_overflow_error(firstlast_min)
    end
    return nothing
end

# Tuples of integers are treated as offsets
# Empty Tuples are handled here
@inline function OffsetArray(A::AbstractArray, offsets::Tuple{Vararg{Integer}}; kw...)
    _checkindices(A, offsets, "offsets")
    OffsetArray{eltype(A), ndims(A), typeof(A)}(A, offsets; kw...)
end

# These methods are necessary to disallow incompatible dimensions for
# the OffsetVector and the OffsetMatrix constructors
for (FT, ND) in ((:OffsetVector, :1), (:OffsetMatrix, :2))
    @eval @inline function $FT(A::AbstractArray{<:Any,$ND}, offsets::Tuple{Vararg{Integer}}; kw...)
        _checkindices(A, offsets, "offsets")
        OffsetArray{eltype(A), $ND, typeof(A)}(A, offsets; kw...)
    end
    FTstr = string(FT)
    @eval @inline function $FT(A::AbstractArray, offsets::Tuple{Vararg{Integer}}; kw...)
        throw(ArgumentError($FTstr*" requires a "*string($ND)*"D array"))
    end
end

## OffsetArray constructors
for FT in (:OffsetArray, :OffsetVector, :OffsetMatrix)
    # Nested OffsetArrays may strip off the wrapper and collate the offsets
    # empty tuples are handled here
    @eval @inline function $FT(A::OffsetArray, offsets::Tuple{Vararg{Int}}; checkoverflow = true)
        _checkindices(A, offsets, "offsets")
        # ensure that the offsets may be added together without an overflow
        checkoverflow && map(overflow_check, axes(A), offsets)
        I = map(+, _offsets(A, parent(A)), offsets)
        $FT(parent(A), I, checkoverflow = false)
    end
    @eval @inline function $FT(A::OffsetArray, offsets::Tuple{Integer,Vararg{Integer}}; kw...)
        $FT(A, map(Int, offsets); kw...)
    end

    # In general, indices get converted to AbstractUnitRanges.
    # CartesianIndices{N} get converted to N ranges
    @eval @inline function $FT(A::AbstractArray, inds::Tuple{Any,Vararg{Any}}; kw...)
        $FT(A, _toAbstractUnitRanges(to_indices(A, axes(A), inds)); kw...)
    end

    # convert ranges to offsets
    @eval @inline function $FT(A::AbstractArray, inds::Tuple{AbstractUnitRange,Vararg{AbstractUnitRange}}; kw...)
        _checkindices(A, inds, "indices")
        # Performance gain by wrapping the error in a function: see https://github.com/JuliaLang/julia/issues/37558
        throw_dimerr(lA, lI) = throw(DimensionMismatch("supplied axes do not agree with the size of the array (got size $lA for the array and $lI for the indices"))
        lA = size(A)
        lI = map(length, inds)
        lA == lI || throw_dimerr(lA, lI)
        $FT(A, map(_offset, axes(A), inds); kw...)
    end

    @eval @inline $FT(A::AbstractArray, inds...; kw...) = $FT(A, inds; kw...)
    @eval @inline $FT(A::AbstractArray; checkoverflow = false) = $FT(A, ntuple(zero, Val(ndims(A))), checkoverflow = checkoverflow)

    @eval @inline $FT(A::AbstractArray, origin::Origin; checkoverflow = true) = $FT(A, origin.index .- first.(axes(A)); checkoverflow = checkoverflow)
end

(o::Origin)(A::AbstractArray) = OffsetArray(no_offset_view(A), o)
Origin(A::AbstractArray) = Origin(first.(axes(A)))

# conversion-related methods
@inline OffsetArray{T}(M::AbstractArray, I...; kw...) where {T} = OffsetArray{T,ndims(M)}(M, I...; kw...)

@inline function OffsetArray{T,N}(M::AbstractArray{<:Any,N}, I...; kw...) where {T,N}
    M2 = _of_eltype(T, M)
    OffsetArray{T,N}(M2, I...; kw...)
end
@inline OffsetArray{T,N}(M::OffsetArray{T,N}, I...; kw...) where {T,N} = OffsetArray(M, I...; kw...)
@inline OffsetArray{T,N}(M::AbstractArray{T,N}, I...; kw...) where {T,N} = OffsetArray{T,N,typeof(M)}(M, I...; kw...)

@inline OffsetArray{T,N,A}(M::AbstractArray{<:Any,N}, I...; kw...) where {T,N,A<:AbstractArray{T,N}} = OffsetArray{T,N,A}(M, I; kw...)
@inline function OffsetArray{T,N,A}(M::AbstractArray{<:Any,N}, I::NTuple{N,Int}; checkoverflow = true) where {T,N,A<:AbstractArray{T,N}}
    checkoverflow && map(overflow_check, axes(M), I)
    Mv = no_offset_view(M)
    MvA = convert(A, Mv)::A
    Iof = map(+, _offsets(M), I)
    OffsetArray{T,N,A}(MvA, Iof, checkoverflow = false)
end
@inline function OffsetArray{T, N, AA}(parent::AbstractArray{<:Any,N}, offsets::NTuple{N, Integer}; kw...) where {T, N, AA<:AbstractArray{T,N}}
    OffsetArray{T, N, AA}(parent, map(Int, offsets)::NTuple{N,Int}; kw...)
end
@inline function OffsetArray{T,N,A}(M::AbstractArray{<:Any,N}, I::Tuple{AbstractUnitRange,Vararg{AbstractUnitRange}}; kw...) where {T,N,A<:AbstractArray{T,N}}
    _checkindices(M, I, "indices")
    # Performance gain by wrapping the error in a function: see https://github.com/JuliaLang/julia/issues/37558
    throw_dimerr(lA, lI) = throw(DimensionMismatch("supplied axes do not agree with the size of the array (got size $lA for the array and $lI for the indices"))
    lM = size(M)
    lI = map(length, I)
    lM == lI || throw_dimerr(lM, lI)
    OffsetArray{T,N,A}(M, map(_offset, axes(M), I); kw...)
end
@inline function OffsetArray{T,N,A}(M::AbstractArray{<:Any,N}, I::Tuple; kw...) where {T,N,A<:AbstractArray{T,N}}
    OffsetArray{T,N,A}(M, _toAbstractUnitRanges(to_indices(M, axes(M), I)); kw...)
end
@inline function OffsetArray{T,N,A}(M::AbstractArray{<:Any,N}; kw...) where {T,N,A<:AbstractArray{T,N}}
    Mv = no_offset_view(M)
    MvA = convert(A, Mv)::A
    OffsetArray{T,N,A}(MvA, _offsets(M); kw...)
end
@inline OffsetArray{T,N,A}(M::A; checkoverflow = false) where {T,N,A<:AbstractArray{T,N}} = OffsetArray{T,N,A}(M, ntuple(zero, Val(N)); checkoverflow = checkoverflow)

Base.convert(::Type{T}, M::AbstractArray) where {T<:OffsetArray} = M isa T ? M : T(M)

@inline AbstractArray{T,N}(M::OffsetArray{S,N}) where {T,S,N} = OffsetArray{T}(M)

# array initialization
@inline function OffsetArray{T,N}(init::ArrayInitializer, inds::Tuple{Vararg{OffsetAxisKnownLength}}; kw...) where {T,N}
    _checkindices(N, inds, "indices")
    AA = Array{T,N}(init, map(_indexlength, inds))
    OffsetArray{T, N, typeof(AA)}(AA, map(_indexoffset, inds); kw...)
end
@inline function OffsetArray{T, N}(init::ArrayInitializer, inds::Tuple; kw...) where {T, N}
    OffsetArray{T, N}(init, _toAbstractUnitRanges(inds); kw...)
end
@inline OffsetArray{T,N}(init::ArrayInitializer, inds...; kw...) where {T,N} = OffsetArray{T,N}(init, inds; kw...)

@inline OffsetArray{T}(init::ArrayInitializer, inds::NTuple{N, OffsetAxisKnownLength}; kw...) where {T,N} = OffsetArray{T,N}(init, inds; kw...)
@inline function OffsetArray{T}(init::ArrayInitializer, inds::Tuple; kw...) where {T}
    OffsetArray{T}(init, _toAbstractUnitRanges(inds); kw...)
end
@inline OffsetArray{T}(init::ArrayInitializer, inds...; kw...) where {T} = OffsetArray{T}(init, inds; kw...)

Base.IndexStyle(::Type{OA}) where {OA<:OffsetArray} = IndexStyle(parenttype(OA))
parenttype(::Type{OffsetArray{T,N,AA}}) where {T,N,AA} = AA
parenttype(A::OffsetArray) = parenttype(typeof(A))

Base.parent(A::OffsetArray) = A.parent

# TODO: Ideally we would delegate to the parent's broadcasting implementation, but that
#       is currently broken in sufficiently many implementation, namely RecursiveArrayTools, DistributedArrays
#       and StaticArrays, that it will take concentrated effort to get this working across the ecosystem.
#       The goal would be to have `OffsetArray(CuArray) .+ 1 == OffsetArray{CuArray}`.
# Base.Broadcast.BroadcastStyle(::Type{<:OffsetArray{<:Any, <:Any, AA}}) where AA = Base.Broadcast.BroadcastStyle(AA)

@inline Base.size(A::OffsetArray) = size(parent(A))

@inline Base.axes(A::OffsetArray) = map(IdOffsetRange, axes(parent(A)), A.offsets)
@inline Base.axes(A::OffsetArray, d) = d <= ndims(A) ? IdOffsetRange(axes(parent(A), d), A.offsets[d]) : IdOffsetRange(axes(parent(A), d))
@inline Base.axes1(A::OffsetArray{T,0}) where {T} = IdOffsetRange(axes(parent(A), 1))  # we only need to specialize this one

# Utils to translate a function to the parent while preserving offsets
unwrap(x) = x, identity
unwrap(x::OffsetArray) = parent(x), data -> OffsetArray(data, x.offsets, checkoverflow = false)
function parent_call(f, x)
    parent, wrap_offset = unwrap(x)
    wrap_offset(f(parent))
end

Base.similar(A::OffsetArray, ::Type{T}, dims::Dims) where T =
    similar(parent(A), T, dims)
function Base.similar(A::AbstractArray, ::Type{T}, shape::Tuple{OffsetAxisKnownLength,Vararg{OffsetAxisKnownLength}}) where T
    # strip IdOffsetRanges to extract the parent range and use it to generate the array
    new_shape = map(_strip_IdOffsetRange, shape)
    # route through _similar_axes_or_length to avoid a stack overflow if map(_strip_IdOffsetRange, shape) === shape
    # This tries to use new_shape directly in similar if similar(A, T, ::typeof(new_shape)) is defined
    # If this fails, it calls similar(A, T, map(_indexlength, new_shape)) to use the size along each axis
    # to generate the new array
    P = _similar_axes_or_length(A, T, new_shape, shape)
    return OffsetArray(P, map(_offset, axes(P), shape))
end
Base.similar(::Type{A}, sz::Tuple{Vararg{Int}}) where {A<:OffsetArray} = similar(Array{eltype(A)}, sz)
function Base.similar(::Type{T}, shape::Tuple{OffsetAxisKnownLength,Vararg{OffsetAxisKnownLength}}) where {T<:AbstractArray}
    new_shape = map(_strip_IdOffsetRange, shape)
    P = _similar_axes_or_length(T, new_shape, shape)
    OffsetArray(P, map(_offset, axes(P), shape))
end
# Try to use the axes to generate the parent array type
# This is useful if the axes have special meanings, such as with static arrays
# This method is hit if at least one axis provided to similar(A, T, axes) is an IdOffsetRange
# For example this is hit when similar(A::OffsetArray) is called,
# which expands to similar(A, eltype(A), axes(A))
_similar_axes_or_length(A, T, ax, ::Any) = similar(A, T, ax)
_similar_axes_or_length(AT, ax, ::Any) = similar(AT, ax)
# Handle the general case by resorting to lengths along each axis
# This is hit if none of the axes provided to similar(A, T, axes) are IdOffsetRanges,
# and if similar(A, T, axes::AX) is not defined for the type AX.
# In this case the best that we can do is to create a mutable array of the correct size
_similar_axes_or_length(A, T, ax::I, ::I) where {I} = similar(A, T, map(_indexlength, ax))
_similar_axes_or_length(AT, ax::I, ::I) where {I} = similar(AT, map(_indexlength, ax))

# reshape accepts a single colon
Base.reshape(A::AbstractArray, inds::OffsetAxis...) = reshape(A, inds)
function Base.reshape(A::AbstractArray, inds::Tuple{OffsetAxis,Vararg{OffsetAxis}})
    AR = reshape(no_offset_view(A), map(_indexlength, inds))
    O = OffsetArray(AR, map(_offset, axes(AR), inds))
    return _popreshape(O, axes(AR), _filterreshapeinds(inds))
end

# Reshaping OffsetArrays can "pop" the original OffsetArray wrapper and return
# an OffsetArray(reshape(...)) instead of an OffsetArray(reshape(OffsetArray(...)))
# Short-circuit for AbstractVectors if the axes are compatible to get around the Base restriction
# to 1-based vectors
function _reshape(A::AbstractVector, inds::Tuple{OffsetAxis})
    @noinline throw_dimerr(ind::Integer) = throw(
        DimensionMismatch("parent has $(size(A,1)) elements, which is incompatible with length $ind"))
    @noinline throw_dimerr(ind) = throw(
        DimensionMismatch("parent has $(size(A,1)) elements, which is incompatible with indices $ind"))
    _checksize(first(inds), size(A,1)) || throw_dimerr(first(inds))
    A
end
_reshape(A, inds) = _reshape2(A, inds)
_reshape2(A, inds) = reshape(A, inds)
# avoid a stackoverflow by relegating to the parent if no_offset_view returns an offsetarray
_reshape2(A::OffsetArray, inds) = reshape(parent(A), inds)
_reshape_nov(A, inds) = _reshape(no_offset_view(A), inds)

Base.reshape(A::OffsetArray, inds::Tuple{OffsetAxis,Vararg{OffsetAxis}}) =
    OffsetArray(_reshape(parent(A), inds), map(_toaxis, inds))
# And for non-offset axes, we can just return a reshape of the parent directly
Base.reshape(A::OffsetArray, inds::Tuple{Union{Integer,Base.OneTo},Vararg{Union{Integer,Base.OneTo}}}) = _reshape_nov(A, inds)
Base.reshape(A::OffsetArray, inds::Dims) = _reshape_nov(A, inds)
Base.reshape(A::OffsetVector, ::Colon) = A
Base.reshape(A::OffsetVector, ::Tuple{Colon}) = A
Base.reshape(A::OffsetArray, ::Colon) = reshape(A, (Colon(),))
Base.reshape(A::OffsetArray, inds::Union{Int,Colon}...) = reshape(A, inds)
Base.reshape(A::OffsetArray, inds::Tuple{Vararg{Union{Int,Colon}}}) = _reshape_nov(A, inds)
# The following two additional methods for Colon are added to resolve method ambiguities to
# Base: https://github.com/JuliaLang/julia/pull/45387#issuecomment-1132859663
Base.reshape(A::OffsetArray, inds::Colon) = _reshape_nov(A, inds)
Base.reshape(A::OffsetArray, inds::Tuple{Colon}) = _reshape_nov(A, inds)

# permutedims in Base does not preserve axes, and can not be fixed in a non-breaking way
# This is a stopgap solution
Base.permutedims(v::OffsetVector) = reshape(v, (1, axes(v, 1)))

Base.fill(v, inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {N} =
    fill!(similar(Array{typeof(v)}, inds), v)
Base.zeros(::Type{T}, inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {T, N} =
    fill!(similar(Array{T}, inds), zero(T))
Base.ones(::Type{T}, inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {T, N} =
    fill!(similar(Array{T}, inds), one(T))
Base.trues(inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {N} =
    fill!(similar(BitArray, inds), true)
Base.falses(inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {N} =
    fill!(similar(BitArray, inds), false)

Base.zero(A::OffsetArray) = parent_call(zero, A)
Base.fill!(A::OffsetArray, x) = parent_call(Ap -> fill!(Ap, x), A)

## Indexing

# Note this gets the index of the parent *array*, not the index of the parent *range*
# Here's how one can think about this:
#   Δi = i - first(r)
#   i′ = first(r.parent) + Δi
# and one obtains the result below.
parentindex(r::IdOffsetRange, i) = i - r.offset

@propagate_inbounds Base.getindex(A::OffsetArray{<:Any,0})  = A.parent[]

@inline function Base.getindex(A::OffsetArray{<:Any,N}, I::Vararg{Int,N}) where N
    @boundscheck checkbounds(A, I...)
    J = map(parentindex, axes(A), I)
    @inbounds parent(A)[J...]
end

@propagate_inbounds Base.getindex(A::OffsetArray{<:Any,N}, c::Vararg{Colon,N}) where N =
    parent_call(x -> getindex(x, c...), A)

# With one Colon we use linear indexing.
# In this case we may forward the index to the parent, as the information about the axes is lost
# The exception to this is with OffsetVectors where the axis information is preserved,
# but that case is handled by getindex(::OffsetArray{<:Any,N}, ::Vararg{Colon,N})
@propagate_inbounds Base.getindex(A::OffsetArray, c::Colon) = A.parent[:]

@inline function Base.getindex(A::OffsetVector, i::Int)
    @boundscheck checkbounds(A, i)
    @inbounds parent(A)[parentindex(Base.axes1(A), i)]
end
@propagate_inbounds Base.getindex(A::OffsetArray, i::Int)  = parent(A)[i]

@inline function Base.setindex!(A::OffsetArray{T,N}, val, I::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A, I...)
    J = map(parentindex, axes(A), I)
    @inbounds parent(A)[J...] = val
    A
end

@inline function Base.setindex!(A::OffsetVector, val, i::Int)
    @boundscheck checkbounds(A, i)
    @inbounds parent(A)[parentindex(Base.axes1(A), i)] = val
    A
end
@propagate_inbounds function Base.setindex!(A::OffsetArray, val, i::Int)
    parent(A)[i] = val
    A
end

@inline Base.iterate(a::OffsetArray, i...) = iterate(parent(a), i...)

Base.in(x, A::OffsetArray) = in(x, parent(A))
Base.copy(A::OffsetArray) = parent_call(copy, A)

Base.strides(A::OffsetArray) = strides(parent(A))
Base.elsize(::Type{OffsetArray{T,N,A}}) where {T,N,A} = Base.elsize(A)
@inline Base.unsafe_convert(::Type{Ptr{T}}, A::OffsetArray{T}) where {T} = Base.unsafe_convert(Ptr{T}, parent(A))

# For fast broadcasting: ref https://discourse.julialang.org/t/why-is-there-a-performance-hit-on-broadcasting-with-offsetarrays/32194
Base.dataids(A::OffsetArray) = Base.dataids(parent(A))
Broadcast.broadcast_unalias(dest::OffsetArray, src::OffsetArray) = parent(dest) === parent(src) ? src : Broadcast.unalias(dest, src)

### Special handling for AbstractRange
const OffsetRange{T} = OffsetVector{T,<:AbstractRange{T}}
const OffsetUnitRange{T} = OffsetVector{T,<:AbstractUnitRange{T}}

Base.step(a::OffsetRange) = step(parent(a))

Base.checkindex(::Type{Bool}, inds::AbstractUnitRange, or::OffsetRange) = Base.checkindex(Bool, inds, parent(or))

# Certain special methods for linear indexing with integer ranges (or OffsetRanges)
# These may bypass the default getindex(A, I...) pathway if the parent types permit this
# For example AbstractUnitRanges and Arrays have special linear indexing behavior defined

# If both the arguments are offset, we may unwrap the indices to call (::OffsetArray)[::AbstractRange{Int}]
@propagate_inbounds function Base.getindex(A::OffsetArray, r::OffsetRange{Int})
    _indexedby(A[parent(r)], axes(r))
end
# If the indices are offset, we may unwrap them and pass the parent to getindex
@propagate_inbounds function Base.getindex(A::AbstractRange, r::OffsetRange{Int})
    _indexedby(A[parent(r)], axes(r))
end

# An OffsetUnitRange might use the rapid getindex(::Array, ::AbstractUnitRange{Int}) for contiguous indexing
@propagate_inbounds function Base.getindex(A::Array, r::OffsetUnitRange{Int})
    B = A[_contiguousindexingtype(parent(r))]
    OffsetArray(B, axes(r), checkoverflow = false)
end

# Linear Indexing of OffsetArrays with AbstractUnitRanges may use the faster contiguous indexing methods
@inline function Base.getindex(A::OffsetArray, r::AbstractUnitRange{Int})
    @boundscheck checkbounds(A, r)
    # nD OffsetArrays do not have their linear indices shifted, so we may forward the indices provided to the parent
    @inbounds B = parent(A)[_contiguousindexingtype(r)]
    _indexedby(B, axes(r))
end
@inline function Base.getindex(A::OffsetVector, r::AbstractUnitRange{Int})
    @boundscheck checkbounds(A, r)
    # OffsetVectors may have their linear indices shifted, so we subtract the offset from the indices provided
    @inbounds B = parent(A)[_subtractoffset(r, A.offsets[1])]
    _indexedby(B, axes(r))
end

# This method added mainly to index an OffsetRange with another range
@inline function Base.getindex(A::OffsetVector, r::AbstractRange{Int})
    @boundscheck checkbounds(A, r)
    @inbounds B = parent(A)[_subtractoffset(r, A.offsets[1])]
    _indexedby(B, axes(r))
end

# In general we would pass through getindex(A, I...) which calls to_indices(A, I) and finally to_index(I)
# An OffsetUnitRange{Int} has an equivalent IdOffsetRange with the same values and axes,
# something similar also holds for OffsetUnitRange{BigInt}
# We may replace the former with the latter in an indexing operation to obtain a performance boost
@inline function Base.to_index(r::OffsetUnitRange{<:Union{Int,BigInt}})
    of = first(axes(r,1)) - 1
    IdOffsetRange(_subtractoffset(parent(r), of), of)
end

@inline function _boundscheck_index_retaining_axes(r, s)
    @boundscheck checkbounds(r, s)
    @inbounds pr = r[UnitRange(s)]
    _indexedby(pr, axes(s))
end
@inline _boundscheck_return(r, s) = (@boundscheck checkbounds(r, s); s)

for OR in [:IIUR, :IdOffsetRange]
    for R in [:StepRange, :StepRangeLen, :LinRange, :UnitRange]
        @eval @inline Base.getindex(r::$R, s::$OR) = _boundscheck_index_retaining_axes(r, s)
    end

    # this method is needed for ambiguity resolution
    @eval @inline function Base.getindex(r::StepRangeLen{T,<:Base.TwicePrecision,<:Base.TwicePrecision}, s::$OR) where T
        _boundscheck_index_retaining_axes(r, s)
    end
end
Base.getindex(r::Base.OneTo, s::IdOffsetRange) = _boundscheck_index_retaining_axes(r, s)

# These methods are added to avoid ambiguities with Base.
# The ones involving Base types should be ported to Base and version-limited here
@inline Base.getindex(r::IdentityUnitRange, s::IIUR) = _boundscheck_return(r, s)
@inline Base.getindex(r::IdentityUnitRange, s::IdOffsetRange) = _boundscheck_return(r, s)
if IdentityUnitRange !== Base.Slice
    @inline Base.getindex(r::Base.Slice, s::IIUR) = _boundscheck_return(r, s)
    @inline Base.getindex(r::Base.Slice, s::IdOffsetRange) = _boundscheck_return(r, s)
end

# eltype conversion
# This may use specialized map methods for the parent
Base.map(::Type{T}, O::OffsetArray) where {T} = parent_call(x -> map(T, x), O)
Base.map(::Type{T}, r::IdOffsetRange) where {T<:Real} = _indexedby(map(T, UnitRange(r)), axes(r))
if eltype(IIUR) === Int
    # This is type-piracy, but there is no way to convert an IdentityUnitRange to a non-Int type in Base
    Base.map(::Type{T}, r::IdentityUnitRange) where {T<:Real} = _indexedby(map(T, UnitRange(r)), axes(r))
end

# mapreduce is faster with an IdOffsetRange than with an OffsetUnitRange
# We therefore convert OffsetUnitRanges to IdOffsetRanges with the same values and axes
function Base.mapreduce(f, op, A1::OffsetUnitRange{<:Integer}, As::OffsetUnitRange{<:Integer}...; kw...)
    As = (A1, As...)
    ofs = map(A -> first(axes(A,1)) - 1, As)
    AIds = map((A, of) -> IdOffsetRange(_subtractoffset(parent(A), of), of), As, ofs)
    mapreduce(f, op, AIds...; kw...)
end

# Optimize certain reductions that treat an OffsetVector as a list
for f in [:minimum, :maximum, :extrema, :sum]
    @eval Base.$f(r::OffsetRange) = $f(parent(r))
end

function Base.show(io::IO, r::OffsetRange)
    show(io, r.parent)
    print(io, " with indices ", UnitRange(axes(r, 1)))
end
Base.show(io::IO, ::MIME"text/plain", r::OffsetRange) = show(io, r)

### Some mutating functions defined only for OffsetVector ###

Base.resize!(A::OffsetVector, nl::Integer) = (resize!(A.parent, nl); A)
Base.push!(A::OffsetVector, x...) = (push!(A.parent, x...); A)
Base.pop!(A::OffsetVector) = pop!(A.parent)
Base.append!(A::OffsetVector, items) = (append!(A.parent, items); A)
Base.empty!(A::OffsetVector) = (empty!(A.parent); A)

# These functions keep the summary compact
function Base.inds2string(inds::Tuple{Vararg{Union{IdOffsetRange, IdentityUnitRange{<:IdOffsetRange}}}})
    Base.inds2string(map(UnitRange, inds))
end
Base.showindices(io::IO, ind1::IdOffsetRange, inds::IdOffsetRange...) = Base.showindices(io, map(UnitRange, (ind1, inds...))...)

function Base.showarg(io::IO, @nospecialize(a::OffsetArray), toplevel)
    print(io, "OffsetArray(")
    Base.showarg(io, parent(a), false)
    Base.showindices(io, axes(a)...)
    print(io, ')')
    if toplevel
        print(io, " with eltype ", eltype(a))
    end
end

function Base.replace_in_print_matrix(A::OffsetArray{<:Any,2}, i::Integer, j::Integer, s::AbstractString)
    J = map(parentindex, axes(A), (i,j))
    Base.replace_in_print_matrix(parent(A), J..., s)
end
function Base.replace_in_print_matrix(A::OffsetArray{<:Any,1}, i::Integer, j::Integer, s::AbstractString)
    ip = parentindex(axes(A,1), i)
    Base.replace_in_print_matrix(parent(A), ip, j, s)
end

no_offset_view(A::OffsetArray) = no_offset_view(parent(A))
no_offset_view(a::Array) = a
no_offset_view(i::Number) = i
no_offset_view(A::AbstractArray) = _no_offset_view(axes(A), A)
_no_offset_view(::Tuple{}, A::AbstractArray{T,0}) where T = A
_no_offset_view(::Tuple{Base.OneTo, Vararg{Base.OneTo}}, A::AbstractArray) = A
# the following method is needed for ambiguity resolution
_no_offset_view(::Tuple{Base.OneTo, Vararg{Base.OneTo}}, A::AbstractUnitRange) = A
_no_offset_view(::Any, A::AbstractArray) = OffsetArray(A, Origin(1))
_no_offset_view(::Any, A::AbstractUnitRange) = UnitRange(A)

#####
# center/centered
# These two helpers are deliberately not exported; their meaning can be very different in
# other scenarios and will be very likely to cause name conflicts if exported.
#####
function center(A::AbstractArray, r::RoundingMode=RoundDown)
    map(axes(A)) do inds
        round(Int, (length(inds)-1)/2, r) + first(inds)
    end
end

centered(A::AbstractArray, cp::Dims=center(A)) = OffsetArray(A, .-cp)

centered(A::AbstractArray, i::CartesianIndex) = centered(A, Tuple(i))

# we may pass the searchsorted* functions to the parent, and wrap the offset
for f in [:searchsortedfirst, :searchsortedlast, :searchsorted]
    _safe_f = Symbol("_safe_" * String(f))
    @eval function $_safe_f(v::OffsetArray, x, ilo, ihi, o::Base.Ordering)
        offset = firstindex(v) - firstindex(parent(v))
        $f(parent(v), x, ilo - offset, ihi - offset, o) .+ offset
    end
    @eval Base.$f(v::OffsetVector, x, ilo::T, ihi::T, o::Base.Ordering) where T<:Integer =
        $_safe_f(v, x, ilo, ihi, o)
end

##
# Deprecations
##

# This is a bad API design as it introduces counter intuitive results (#250)
@deprecate centered(A::AbstractArray, r::RoundingMode) OffsetArray(A, .-center(A, r)) false

end # module

max_stars_repo_name: TechPenguineer/julia

max_issues_repo_name: TechPenguineer/julia

max_forks_repo_name: Nexuscompute/julia


text:
@defcomp climatech4cycle begin
    # Global CH4 emissions in Mt of CH4
    globch4 = Parameter(index=[time])

    # Atmospheric CH4 concentration
    acch4 = Variable(index=[time])

    # CH4 decay
    ch4decay = Variable()

    lifech4 = Parameter(default = 12.0)

    #  CH4 pre industrial
    ch4pre = Parameter(default = 790.0)

    function run_timestep(p, v, d, t)

        if is_first(t)
            v.ch4decay = 1.0 / p.lifech4

            v.acch4[1] = 1222.0
        else
            # Calculate CH4 concentrations
            v.acch4[t] = v.acch4[t - 1] + 0.3597 * p.globch4[t] - v.ch4decay * (v.acch4[t - 1] - p.ch4pre)

            if v.acch4[t] < 0
                error("ch4 atmospheric concentration out of range in $t")
            end
        end
    end
end


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license

export mass, set_mass!, charge, set_charge!, name, set_name!, fullname,
vdw_radius, covalent_radius, atomic_number, atom_type, set_atom_type!,
set_property!, property, AtomType

"""
Create an atom with the given ``name`` and set the atom ``type`` to be the same
as ``name``.
"""
function Atom(name::String)
    return Atom(lib.chfl_atom(pointer(name)))
end

"""
Get a copy of the ``atom`` at the given ``index`` from a ``frame``
"""
function Atom(frame::Frame, index::Integer)
    handle = lib.chfl_atom_from_frame(frame.handle, UInt64(index))
    return Atom(handle)
end

"""
Get a copy of the ``atom`` at a given ``index`` from a ``topology``
"""
function Atom(topology::Topology, index::Integer)
    handle = lib.chfl_atom_from_topology(topology.handle, UInt64(index))
    return Atom(handle)
end

"""
Free the allocated memory for the ``Atom`` object.
"""
function _free(atom::Atom)
    lib.chfl_atom_free(atom.handle)
end

"""
Get the mass of an ``atom``.

The mass is given in atomic mass units.
"""
function mass(atom::Atom)
    m = Ref{Float64}(0)
    _check(
        lib.chfl_atom_mass(atom.handle, m)
    )
    return m[]
end

"""
Set the mass of an ``atom`` to ``mass``.

The mass must be in atomic mass units.
"""
function set_mass!(atom::Atom, mass)
    _check(
        lib.chfl_atom_set_mass(atom.handle, Float64(mass))
    )
    return nothing
end

"""
Get the charge of an ``atom``.

The charge is in number of the electron charge *e*.
"""
function charge(atom::Atom)
    c = Ref{Float64}(0)
    _check(
        lib.chfl_atom_charge(atom.handle, c)
    )
    return c[]
end

"""
Set the charge of an ``atom`` to ``charge``.

The charge must be in number of the electron charge *e*.
"""
function set_charge!(atom::Atom, charge)
    _check(
        lib.chfl_atom_set_charge(atom.handle, Float64(charge))
    )
    return nothing
end

"""
Get the name of an ``atom``
"""
function name(atom::Atom)
    return _call_with_growing_buffer(
        (buffer, size) -> _check(lib.chfl_atom_name(atom.handle, buffer, size))
    )
end

"""
Set the name of an ``atom`` to ``name``.
"""
function set_name!(atom::Atom, name::String)
    _check(
        lib.chfl_atom_set_name(atom.handle, pointer(name))
    )
    return nothing
end

"""
Get the type of an ``atom``.
"""
function atom_type(atom::Atom)
    return _call_with_growing_buffer(
        (buffer, size) -> _check(lib.chfl_atom_type(atom.handle, buffer, size))
    )
end

"""
Set the type of an ``atom`` to ``type``.
"""
function set_atom_type!(atom::Atom, atom_type::String)
    _check(
        lib.chfl_atom_set_type(atom.handle, pointer(atom_type))
    )
    return nothing
end

"""
Get the full name of an ``atom`` from the atom type.

For example, the full name of an atom with type "He" is "Helium".
"""
function Base.fullname(atom::Atom)
    return _call_with_growing_buffer(
        (buffer, size) -> _check(lib.chfl_atom_full_name(atom.handle, buffer, size))
    )
end

"""
Get the Van der Waals radius of an ``atom`` from the atom type.

If the radius can not be found, this function returns -1.
"""
function vdw_radius(atom::Atom)
    radius = Ref{Float64}(0)
    _check(
        lib.chfl_atom_vdw_radius(atom.handle, radius)
    )
    return radius[]
end

"""
Get the covalent radius of an ``atom`` from the atom type.

If the radius can not be found, returns -1.
"""
function covalent_radius(atom::Atom)
    radius = Ref{Float64}(0)
    _check(
        lib.chfl_atom_covalent_radius(atom.handle, radius)
    )
    return radius[]
end

"""
Get the atomic number of an ``atom`` from the atom type.

If the atomic number can not be found, returns 0.
"""
function atomic_number(atom::Atom)
    number = Ref{UInt64}(0)
    _check(
        lib.chfl_atom_atomic_number(atom.handle, number)
    )
    return number[]
end

"""
Set a named property for the given atom.
"""
function set_property!(atom::Atom, name::String, property)

    prop = Property(property)

    _check(
        lib.chfl_atom_set_property(atom.handle, pointer(name), prop.handle)
    )
    return nothing
end

"""
Get a named property for the given atom.
"""
function property(atom::Atom, name::String)
    get(Property(lib.chfl_atom_get_property(atom.handle, pointer(name))))
end

"""
Make a deep copy of an ``atom``.
"""
function Base.deepcopy(atom::Atom)
    handle = lib.chfl_atom_copy(atom.handle)
    return Atom(handle)
end


text:
include("../src/autodiff_module.jl")
include("../../__lib__/math/opti/src/opti_module.jl")

module optitest
using ..autodiff
using ..opti

function quadratic(x)
    y = zero(eltype(x))
    for ii in eachindex(x)
        y += x[ii]^2
    end
    return y
end

function quadratic(x, x0)
    y = zero(eltype(x))
    for ii in eachindex(x)
        y += (x[ii] - x0[ii])^2
    end
    return y
end

# x1a = autodiff.DualNumber(3.0)
# x2a = autodiff.DualNumber(4.0)

# x_vec = [x1a, x2a]
# g = autodiff.gradient(quadratic, x_vec)

x0 = [3.0, 4.0]

c = [5.0, 6.0]
obj_fcn(x)              = quadratic(x, c)
grad_fcn(objfcn, x)     = autodiff.gradient(autodiff.ForwardMode(), objfcn, x)
dir_diff_fcn(fcn, x, s) = autodiff.directional_derivative(fcn, x, s)

tol_options      = opti.ToleranceOptions(
    step_size_tol = 1e-15,
    f_abs_tol = 1e-15,
    f_rel_tol = 1e-15,
    g_abs_tol = 1e-15,
    max_iter = 100)

# (x_sol, y_sol, iter, stopping_crit) = opti.quasi_newton(
#     obj_fcn,
#     x0;
#     tol_options          = tol_options,
#     gradient_fcn         = grad_fcn,
#     directional_diff_fcn = dir_diff_fcn)

stat = opti.cg_pr(
    obj_fcn,
    x0;
    tol_options          = tol_options,
    gradient_fcn         = grad_fcn,
    directional_diff_fcn = dir_diff_fcn)

# directional diff
# lim[h -> 0] (fcn(x + h * s) - fcn(x)) / h

# x0 = [3.0, 4.0]
# obj_fcn(x) = quadratic(x)

# # s = [-4.0, 3.0]
# # s = -[-3.0, 4.0]
# s = [0.0, 5.0]
# s = [5.0, 0.0]

# dg = autodiff.directional_derivative(obj_fcn, x0, s)

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
function testbed_chain(
    device,
    S,
    U = complex(1.0, 0.0),
    Y = complex(270.0, -380.0),
    P_ref = t -> -2,
    Q_ref = t -> 0.0,
)
    busses = [SlackAlgebraic(; U = U), device, PQAlgebraic(; S = S)]
    lines = [
        StaticLine(; from = 1, to = 2, Y = Y),
        StaticLine(; from = 2, to = 3, Y = Y),
    ]
    pg = PowerGrid(busses, lines)
    ADN(pg, typeof(device), P_ref, Q_ref)
end

function testbed(
    device,
    U = complex(1.0, 0.0);
    P_ref = t -> -2,
    Q_ref = t -> 0.0,
)
    T = OLTC(
        from = 1,
        to = 2,
        Uos = 110E3 / base_voltage_HV, # Bemessungsspannung Oberspannungsseite in kV
        Uus = 20E3 / base_voltage,# Bemessungsspannung Unterspannungsseite in kV
        k = 0, # Kennzahl der Schaltgruppe, Bsp: YD5-> 5
        ssp = 0.0, # Stufenschalterposition
        stufe_us = 0.625, # Stufung pro Stufenschalterstellung in %
        Sr = 25E6 / base_power, #Bemessungsscheinleistung in MVA
        uk = 12.0, # Kurzschlussspannung in %
        Pvk = 25E3 / base_power, # Kupferverluste in MW
        Pvl = 0.0, # Eisenverluste in kW
        iLeer = 0.0, # Leerlaufstrom in %
    )
    busses = [SlackAlgebraic(; U = U), device]
    lines = [T,]
    pg = PowerGrid(busses, lines)
    ADN(pg, typeof(device), P_ref, Q_ref)
end

function testbed_dynamicline(
    device,
    U = complex(1.0, 0.0),
    R = 0.5,
    L = 0.7;
    P_ref = t -> -2,
    Q_ref = t -> 0.0,
)
    busses = [SlackAlgebraic(; U = U), device, device]
    lines = [
        RLLine(; from = 1, to = 2, R = R, L = L, ω0 = ω),
        RLLine(; from = 1, to = 2, R = R, L = L, ω0 = ω),
    ]
    pg = PowerGrid(busses, lines)
    ADN(pg, typeof(device), P_ref, Q_ref)
end

function testbed_line(
    device,
    U = complex(1.0, 0.0),
    Y = complex(270.0, -380.0);
    P_ref = t -> -2,
    Q_ref = t -> 0.0,
)
    busses = [SlackAlgebraic(; U = U), device, device]
    lines = [
        StaticLine(; from = 1, to = 2, Y = Y),
        StaticLine(; from = 2, to = 3, Y = Y),
    ]
    pg = PowerGrid(busses, lines)
    ADN(pg, typeof(device), P_ref, Q_ref)
end

# # all values in p.u.
# device = DGUnitPLLPQTracking(;
#     K_pll = 0.1, #Hz/pu
#     K_PT1 = 1.0, # unit: [y]/[u] = [P]/[P] = 1
#     T_PT1 = 10.0, # unit: s
#     S_pq = complex(-1.0, -0.25), # pu
#     Y_n = 0.0,
# )

function testbed_CIGRE(
    DG_unit_type;
    idx = 1:12,
    P_ref = t -> -2,
    Q_ref = t -> 0.0,
)
    busses_static, lines, T, elist, Zs, Yshs = CIGRE_static()
    dir = @__DIR__
    power_flow = pf_sol(
        busses_static,
        lines;
        matfile = "$dir/LFCigre.mat",
        outfile = nothing,
    )
    busses = copy(busses_static)
    DG_locs = 2:12
    for i in DG_locs
        S_bkgrnd = zero(im)
        try
            S_bkgrnd = busses_static[i].S
        catch
            S_bkgrnd = complex(busses_static[i].P, busses_static[i].Q)
        end
        busses[i] = eval(Symbol(DG_unit_type))(;
            K_pll = 1632.993, #Hz/pu
            K_PT1 = 1.0,
            T_PT1 = 0.1,
            K_FRT = 2.0,
            I_max = 1.0, # pu
            Vref = 1.0, #power_flow.V[i], #1.0, # pu
            Vdead = 0.1, #pu
            S_pq = S_bkgrnd, # TODO: check sign of PQ node
            Y_n = 0.0,
        )
    end
    pg = PowerGrid(busses[idx], lines[1:length(idx)-1])
    ADN(pg, DG_unit_type, P_ref, Q_ref), power_flow
end

function testbed_CIGRE_dynamic(
    DG_unit_type;
    idx = 1:12,
    P_ref = t -> -2,
    Q_ref = t -> 0.0,
)
    busses_static, lines_static, T, elist, Zs, Yshs = CIGRE_static()
    busses = copy(busses_static)
    DG_locs = 2:12
    for i in DG_locs
        S_bkgrnd = zero(im)
        try
            S_bkgrnd = busses_static[i].S
        catch
            S_bkgrnd = complex(busses_static[i].P, busses_static[i].Q)
        end
        busses[i] = eval(Symbol(DG_unit_type))(;
            K_pll = 1632.993, #Hz/pu
            K_PT1 = 1.0,
            T_PT1 = 0.1,
            K_FRT = 2.0,
            I_max = 1.0, # pu
            Vref = 1.0, #power_flow.V[i], #1.0, # pu
            Vdead = 0.1, #pu
            S_pq = S_bkgrnd, # TODO: check sign of PQ node
            Y_n = 0.0,
        )
    end

    lines = Array{AbstractLine,1}([])
    # # For ssp=0, replace it with PiModelLine
    # # assume t_km = 1, t_mk real
    Y = PiModel(T) #./ admittance_m
    ue = sqrt(-Y[2, 2] / Y[1, 1])
    Z = inv(Y[1, 2] / ue)
    TD = RLLine(
        from = 1,
        to = 2,
        R = real(Z),
        L = (imag(Z) / ω),
        ω0 = ω
    )
    push!(lines, TD)
    for (e, Z, Ysh) in zip(elist, Zs, Yshs)
        push!(
            lines,
            RLLine(
                from = first(e),
                to = last(e),
                R = real(Z) * base_admittance,
                L = (imag(Z) / ω)  * base_admittance,
                ω0 = ω
            ),
        )
    end
    pg = PowerGrid(busses[idx], lines[1:length(idx)-1])
    ADN(pg, DG_unit_type, P_ref, Q_ref)
end

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


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

"""
    pointgradients(img, points; [method])

Compute the gradients along all dimensions of N-dimensional `img`
at `points` using `method`. Default method is `:ando3`.
"""
function pointgradients(img::AbstractArray{T,N},
                        points::AbstractVector;
                        method=:ando3) where {N,T}
  dims = size(img)
  npts = length(points)

  # smoothing weights
  weights = (method == :sobel ? [1,2,1] :
             method == :ando3 ? [.112737,.274526,.112737] :
             error("Unknown gradient method: $method"))

  # pad input image
  padimg = padarray(img, Pad(:replicate, 1, 1))

  # gradient matrix
  G = zeros(npts, N)

  # compute gradient for all directions at specified points
  for i in 1:N
    # kernel = centered difference + perpendicular smoothing
    if dims[i] > 1
      # centered difference
      kern = reshape([-1,0,1], ntuple(k -> k == i ? 3 : 1, N))
      # perpendicular smoothing
      for j in setdiff(1:N, i)
        if dims[j] > 1
          wvec = reshape(weights, ntuple(k -> k == j ? 3 : 1, N))
          kern = broadcast(*, kern, wvec)
        end
      end

      A = zeros(size(kern))
      for (k, icenter) in enumerate(points)
        i1 = CartesianIndex(ntuple(i->1, N))
        for ii in CartesianIndices(A)
          A[ii] = padimg[ii + icenter - i1]
        end

        G[k,i] = sum(kern .* A)
      end
    end
  end

  G
end


text:
#=ЗАДАЧА:На прямоугольном поле произвольных размеров расставить маркеры в виде "шахматных" клеток,
 начиная с юго-западного угла поля, когда каждая отдельная "шахматная" клетка имеет размер n x n клеток поля (n - это параметр функции).
 Начальное положение Робота - произвольное, конечное - совпадает с начальным.
 Клетки на севере и востоке могут получаться"обрезанными" - зависит от соотношения размеров поля и "шахматных" клеток.
 (Подсказка: здесь могут быть полезными две глобальных переменных,
  в которых будут содержаться текущие декартовы координаты Робота относительно начала координат в левом нижнем углу поля, например)
=#
function chess(r::Robot)
    to_Sud=moves!(r,Sud)
    to_West=moves!(r,West)
    putmarker!(r)#маркер в углу

    num_steps_1=checking(r,Nord)
    back(r,Sud)
    num_steps_2=checking(r,Ost)
    back(r,West)
    if num_steps_1<num_steps_2
        steps=num_steps_1
    else
        steps=num_steps_2
    end
        putmarkers!(r,Nord,steps)
        back(r,Sud)
        putmarkers!(r,Ost,steps)
        back(r,West)
    for i in 1:steps
        move!(r,Nord)
        putmarkers!(r,Ost,steps)
        back(r,West)
    end
    back(r,Sud)

    moves!(r,Ost,to_West)
    moves!(r,Nord,to_Sud)
end
function check_of_marker(r::Robot,side::HorizonSide)#если робот на клетке с маркером-следующая-без(и наоборот)
     if ismarker(r)==true
        move!(r,side)
     else
        move!(r,side)
        putmarker!(r)
     end
end
function putmarkers!(r::Robot, side::HorizonSide,steps::Int)
    for j in 1:steps
        check_of_marker(r,side)
    end
end
function checking(r::Robot,side::HorizonSide)
     num_steps=0
     while  isborder(r,side)==false
        move!(r,side)
        num_steps+=1
     end
    return num_steps
end
function back(r::Robot,side::HorizonSide)
    while isborder(r,side)==false
        move!(r,side)
    end
end
function moves!(r::Robot, side::HorizonSide)
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps+=1
    end
    return num_steps
end

function moves!(r::Robot,side::HorizonSide,num_steps::Int)
    for _ in 1:num_steps
        move!(r,side)
    end
end

max_stars_repo_name: ilya-design/mirea-progs

max_issues_repo_name: ilya-design/mirea-progs

max_forks_repo_name: ilya-design/mirea-progs


text:
function log_network(
    state::Tuple{AbstractGraph, AbstractArray}, tick_nr::Int64
)
    graph, agent_list = state
    agent_opinion = [a.opinion for a in agent_list]
    agent_perceiv_publ_opinion = [a.perceiv_publ_opinion for a in agent_list]
    agent_indegree = indegree(graph)
    agent_outdegree = outdegree(graph)
    agent_centrality = closeness_centrality(graph)
    agent_cc = local_clustering_coefficient(graph)
    agent_component = (
        [i for i in keys(connected_components(graph)), v in 1:nv(graph)
        if v in connected_components(graph)[i]]
    )

    return DataFrame(
        TickNr = tick_nr,
        AgentID = 1:length(agent_list),
        Opinion = agent_opinion,
        PerceivPublOpinion = agent_perceiv_publ_opinion,
        Indegree = agent_indegree,
        Outdegree = agent_outdegree,
        Centrality = agent_centrality,
        CC = agent_cc,
        Component = agent_component
    )
end

# suppress output of include()
;


text:
import Base: @deprecate

## deprecated in v0.12

@deprecate(Cart_get(comm, maxdims), Cart_get(comm), false)

@deprecate(Cart_coords!(comm, rank, maxdims::Integer, coords),
           Cart_coords!(comm, rank, coords), false)

@deprecate(Cart_coords(comm, maxdims::Integer), Cart_coords(comm), false)

## deprecated in v0.11

@deprecate Reduce_in_place!(buf, count::Integer, op, root::Integer, comm::Comm) begin
    @assert count == length(buf)
    Reduce!(buf, op, root, comm)
end false

@deprecate Scatter_in_place!(buf, count::Integer, root::Integer, comm::Comm) begin
    if root == MPI.Comm_rank(comm)
        Scatter!(buf, nothing, count, root, comm)
    else
        Scatter!(nothing, buf, count, root, comm)
    end
end false

@deprecate Scatterv_in_place!(buf, counts::Vector, root::Integer, comm::Comm) begin
    if root == MPI.Comm_rank(comm)
        Scatterv!(buf, nothing, counts, root, comm)
    else
        Scatterv!(nothing, buf, counts, root, comm)
    end
end false

@deprecate Gather_in_place!(buf, count::Integer, root::Integer, comm::Comm) begin
    if root == MPI.Comm_rank(comm)
        Gather!(nothing, buf, count, root, comm)
    else
        Gather!(buf, nothing, count, root, comm)
    end
end false

@deprecate Gatherv_in_place!(buf, counts::Vector{Cint}, root::Integer, comm::Comm) begin
    if root == MPI.Comm_rank(comm)
        Gatherv!(nothing, buf, counts, root, comm)
    else
        Gatherv!(buf, nothing, counts, root, comm)
    end
end false

@deprecate(mpitype(T), Datatype(T), false)

@deprecate(Type_Create_Subarray(ndims::Integer, sizes::MPIBuffertype{Cint}, subsizes::MPIBuffertype{Cint},
                                starts::MPIBuffertype{Cint}, order::Integer, oldtype),
           Types.create_subarray(sizes, subsizes, starts, Datatype(oldtype); rowmajor = order == MPI_ORDER_C), false)
@deprecate(Type_Create_Struct(nfields::Integer, blocklengths::MPIBuffertype{Cint},
                              displacements::MPIBuffertype{Cptrdiff_t}, types::MPIBuffertype{MPI_Datatype}),
           Types.create_struct(blocklengths, displacements, types), false)
@deprecate(Type_Commit!(datatype), Types.commit!(datatype), false)

@deprecate(Send(buf, count::Integer, datatype::Datatype, dest::Integer, tag::Integer, comm::Comm),
           Send(Buffer(buf, count, datatype), dest, tag, comm), false)
@deprecate(Send(buf::AbstractArray, count::Integer, dest::Integer, tag::Integer, comm::Comm),
           Send(view(buf, 1:count), dest, tag, comm), false)
@deprecate(Send(buf::Ref, count::Integer, dest::Integer, tag::Integer, comm::Comm),
           Send(buf, dest, tag, comm), false)

@deprecate(Isend(buf, count::Integer, datatype::Datatype, dest::Integer, tag::Integer, comm::Comm),
           Isend(Buffer(buf,count,datatype), dest, tag, comm), false)
@deprecate(Isend(buf::AbstractArray, count::Integer, dest::Integer, tag::Integer, comm::Comm),
           Isend(view(buf,1:count), dest, tag, comm), false)
@deprecate(Isend(buf::Ref, count::Integer, dest::Integer, tag::Integer, comm::Comm),
           Isend(buf, dest, tag, comm), false)

@deprecate(Recv!(buf, count::Integer, datatype::Datatype, src::Integer, tag::Integer, comm::Comm),
           Recv!(Buffer(buf, count, datatype), src, tag, comm), false)
@deprecate(Recv!(buf::AbstractArray, count::Integer, src::Integer, tag::Integer, comm::Comm),
           Recv!(view(buf, 1:count), src, tag, comm), false)
@deprecate(Recv!(buf::Ref, count::Integer, src::Integer, tag::Integer, comm::Comm),
           Recv!(buf, src, tag, comm), false)

@deprecate(Irecv!(buf, count::Integer, datatype::Datatype, src::Integer, tag::Integer, comm::Comm),
           Irecv!(Buffer(buf,count,datatype), src, tag, comm), false)
@deprecate(Irecv!(buf::AbstractArray, count::Integer, src::Integer, tag::Integer, comm::Comm),
           Irecv!(view(buf,1:count), src, tag, comm), false)
@deprecate(Irecv!(buf::Ref, count::Integer, src::Integer, tag::Integer, comm::Comm),
           Irecv!(buf, src, tag, comm), false)

@deprecate(Sendrecv!(sendbuf, sendcount::Integer, sendtype,   dest::Integer, sendtag::Integer,
                     recvbuf, recvcount::Integer, recvtype, source::Integer, recvtag::Integer,
                     comm::Comm),
           Sendrecv!(Buffer(sendbuf, sendcount, sendtype), dest, sendtag,
                     Buffer(recvbuf, recvcount, recvtype), source, recvtag, comm), false)
@deprecate(Sendrecv!(sendbuf, sendcount::Integer,   dest::Integer, sendtag::Integer,
                     recvbuf, recvcount::Integer, source::Integer, recvtag::Integer,
                     comm::Comm),
           Sendrecv!(view(sendbuf, 1:sendcount), dest, sendtag,
                     view(recvbuf, 1:recvcount), source, recvtag, comm), false)


text:
#! /usr/bin/julia

# Rosetta Code, Date format

ts = time()

println("Today's date is:")
println("    ", strftime("%F", ts))
println("    ", strftime("%A, %B %d, %Y", ts))

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
export iauPmpx
"""
Proper motion and parallax.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   rc,dc  double     ICRS RA,Dec at catalog epoch (radians)
   pr     double     RA proper motion (radians/year; Note 1)
   pd     double     Dec proper motion (radians/year)
   px     double     parallax (arcsec)
   rv     double     radial velocity (km/s, +ve if receding)
   pmt    double     proper motion time interval (SSB, Julian years)
   pob    double[3]  SSB to observer vector (au)

Returned:
   pco    double[3]  coordinate direction (BCRS unit vector)

Notes:

   1. The proper motion in RA is dRA/dt rather than cos(Dec)*dRA/dt.

   2. The proper motion time interval is for when the starlight
   reaches the solar system barycenter.

   3). To avoid the need for iteration, the Roemer effect (i.e. the
   small annual modulation of the proper motion coming from the
   changing light time) is applied approximately, using the
   direction of the star at the catalog epoch.

References:

   1984 Astronomical Almanac, pp B39-B41.

   Urban, S. & Seidelmann, P. K. (eds), Explanatory Supplement to
   the Astronomical Almanac, 3rd ed., University Science Books
   (2013), Section 7.2.

Called:
   iauPdp       scalar product of two p-vectors
   iauPn        decompose p-vector into modulus and direction

This revision:   2013 October 9

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauPmpx(double rc, double dc, double pr, double pd,
#              double px, double rv, double pmt, double pob[3],
#              double pco[3])

function iauPmpx(rc::Real, dc::Real, pr::Real, pd::Real,
                 px::Real, rv::Real, pmt::Real, pob::AbstractVector{<:Real})
   pco = zeros(Float64, 3)

   ccall((:iauPmpx, libsofa_c), Cvoid,
       (Cdouble, Cdouble, Cdouble, Cdouble,
        Cdouble, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}),
       convert(Float64, rc), convert(Float64, dc),
       convert(Float64, pr), convert(Float64, pd),
       convert(Float64, px), convert(Float64, rv),
       convert(Float64, pmt),
       convert(Array{Float64, 1}, pob),
       pco)

   return SVector{3}(pco)
end


text:
struct WrappedLearner <: Learner
    learner
    stop
    askhook
    tellhook
end
WrappedLearner(learner,stop) = WrappedLearner(learner, stop, (x,y)->nothing, (x,y)->nothing)

function ask!(learner::WrappedLearner,input)
    if learner.stop(learner.learner)
        val = nothing
    else
        val = ask!(learner.learner,input)
    end
    learner.askhook(learner.learner,val)
    return val
end

function tell!(learner::WrappedLearner,message)
    tell!(learner.learner,message)
    learner.tellhook(learner.learner,message)
end

### The interface probabily needs to be a little different
mutable struct IgnorantLearner <: Learner
    iter
    state
end
IgnorantLearner(iter) = IgnorantLearner(iter,nothing)

function ask!(learner::IgnorantLearner,input)

    if input==nothing
        return nothing
    end

    temp = learner.state==nothing ? iterate(learner.iter) : iterate(learner.iter,learner.state)
    if temp==nothing
        return nothing
    else
        val,learner.state = temp
        return val
    end
end

tell!(learner::IgnorantLearner,m) = nothing

### Something which would help to reproduce the state of the Learner.
mutable struct HistoryMaster <: Master
    tasks
    results
    np
end

nparalel(master::HistoryMaster) = master.np

function HistoryMaster(output,np)
    N = length(output)
    tasks = Channel{Any}(N)
    results = Channel{Tuple{Any,Any}}(N)

    for oi in output
        put!(results,oi)
    end

    HistoryMaster(tasks,results,np)
end

###

struct PlayRes
    input::Union{Nothing,Vector}
    output::Vector
    nprocs::Int
end

function play!(learner, f::Function, stop::Function)#,cond::Function) ### I could also pass WorkerPool here
    master = WorkMaster(f)
    nprocs = length(master.slaves)
    loop = Loop(master, learner)
    nodes = evaluate!(loop, stop)
    releaseall!(master)
    return PlayRes(nothing, nodes, nprocs)
end

function replay!(learner,nodes::Array{Any}) ### I could make a st function as an optional argument
    for n in nodes
        tell!(learner,n)
    end
end

function replay!(learner, ldata::PlayRes, n::Int)
    @assert n <= length(ldata.output)
    master = HistoryMaster(ldata.output,ldata.nprocs)
    loop = Loop(master,learner)
    iter = 1:n
    evaluate!(loop,iter)
end

replay!(learner, ldata::PlayRes) = replay!(learner, ldata, length(ldata.output))


text:
# Counts of discrete values

#################################################
#
#  counts on given levels
#
#################################################

@compat IntUnitRange{T<:Integer} = UnitRange{T}

#### functions for counting a single list of integers (1D)
"""
    addcounts!(r, x, levels, [wv])

Add the number of occurrences in `x` of each value in `levels` to an existing
array `r`. If a weighting vector `wv` is specified, the sum of weights is used
rather than the raw counts.
"""
function addcounts!(r::AbstractArray, x::IntegerArray, levels::IntUnitRange)
    # add counts of integers from x to r

    k = length(levels)
    length(r) == k || throw(DimensionMismatch())

    m0 = levels[1]
    m1 = levels[end]
    b = m0 - 1

    @inbounds for i in 1 : length(x)
        xi = x[i]
        if m0 <= xi <= m1
            r[xi - b] += 1
        end
    end
    return r
end

function addcounts!(r::AbstractArray, x::IntegerArray, levels::IntUnitRange, wv::WeightVec)
    k = length(levels)
    length(r) == k || throw(DimensionMismatch())

    m0 = levels[1]
    m1 = levels[end]
    b = m0 - 1
    w = values(wv)

    @inbounds for i in 1 : length(x)
        xi = x[i]
        if m0 <= xi <= m1
            r[xi - b] += w[i]
        end
    end
    return r
end

"""
    counts(x, levels=span(x), [wv::WeightVec])

Count the number of times that values in the range `levels` occur in
`x`. The output is a vector of length `length(levels)`. If a weighting
vector `wv` is specified, the sum of the weights is used rather than the
raw counts.
"""
counts(x::IntegerArray, levels::IntUnitRange) =
    addcounts!(zeros(Int, length(levels)), x, levels)
counts(x::IntegerArray, levels::IntUnitRange, wv::WeightVec) =
    addcounts!(zeros(eltype(wv), length(levels)), x, levels, wv)

"""
    counts(x, k::Integer, [wv::WeightVec])

Count the number of times integers in the range 1 to `k` occur in `x`.
"""
counts(x::IntegerArray, k::Integer) = counts(x, 1:k)
counts(x::IntegerArray, k::Integer, wv::WeightVec) = counts(x, 1:k, wv)
counts(x::IntegerArray) = counts(x, span(x))
counts(x::IntegerArray, wv::WeightVec) = counts(x, span(x), wv)

"""
    proportions(x, levels=span(x), [wv::WeightVec])

Return the proportion of values in the range `levels` that occur in `x`.
Equivalent to `counts(x, levels) / length(x)`. If a weighting vector `wv`
is specified, the sum of the weights is used rather than the raw counts.
"""
proportions(x::IntegerArray, levels::IntUnitRange) = counts(x, levels) .* inv(length(x))
proportions(x::IntegerArray, levels::IntUnitRange, wv::WeightVec) =
    counts(x, levels, wv) .* inv(sum(wv))

"""
    proportions(x, k::Integer, [wv::WeightVec])

Return the proportion of integers in 1 to `k` that occur in `x`.
"""
proportions(x::IntegerArray, k::Integer) = proportions(x, 1:k)
proportions(x::IntegerArray, k::Integer, wv::WeightVec) = proportions(x, 1:k, wv)
proportions(x::IntegerArray) = proportions(x, span(x))
proportions(x::IntegerArray, wv::WeightVec) = proportions(x, span(x), wv)

#### functions for counting a single list of integers (2D)

function addcounts!(r::AbstractArray, x::IntegerArray, y::IntegerArray, levels::NTuple{2,IntUnitRange})
    # add counts of integers from x to r

    n = length(x)
    length(y) == n || throw(DimensionMismatch())

    xlevels, ylevels = levels

    kx = length(xlevels)
    ky = length(ylevels)
    size(r) == (kx, ky) || throw(DimensionMismatch())

    mx0 = xlevels[1]
    mx1 = xlevels[end]
    my0 = ylevels[1]
    my1 = ylevels[end]

    bx = mx0 - 1
    by = my0 - 1

    for i = 1:n
        xi = x[i]
        yi = y[i]
        if (mx0 <= xi <= mx1) && (my0 <= yi <= my1)
            r[xi - bx, yi - by] += 1
        end
    end
    return r
end

function addcounts!(r::AbstractArray, x::IntegerArray, y::IntegerArray,
                    levels::NTuple{2,IntUnitRange}, wv::WeightVec)
    # add counts of integers from x to r

    n = length(x)
    length(y) == length(wv) == n || throw(DimensionMismatch())

    xlevels, ylevels = levels

    kx = length(xlevels)
    ky = length(ylevels)
    size(r) == (kx, ky) || throw(DimensionMismatch())

    mx0 = xlevels[1]
    mx1 = xlevels[end]
    my0 = ylevels[1]
    my1 = ylevels[end]

    bx = mx0 - 1
    by = my0 - 1
    w = values(wv)

    for i = 1:n
        xi = x[i]
        yi = y[i]
        if (mx0 <= xi <= mx1) && (my0 <= yi <= my1)
            r[xi - bx, yi - by] += w[i]
        end
    end
    return r
end

# facet functions

function counts(x::IntegerArray, y::IntegerArray, levels::NTuple{2,IntUnitRange})
    addcounts!(zeros(Int, length(levels[1]), length(levels[2])), x, y, levels)
end

function counts(x::IntegerArray, y::IntegerArray, levels::NTuple{2,IntUnitRange}, wv::WeightVec)
    addcounts!(zeros(eltype(wv), length(levels[1]), length(levels[2])), x, y, levels, wv)
end

counts(x::IntegerArray, y::IntegerArray, levels::IntUnitRange) =
    counts(x, y, (levels, levels))
counts(x::IntegerArray, y::IntegerArray, levels::IntUnitRange, wv::WeightVec) =
    counts(x, y, (levels, levels), wv)

counts(x::IntegerArray, y::IntegerArray, ks::NTuple{2,Integer}) =
    counts(x, y, (1:ks[1], 1:ks[2]))
counts(x::IntegerArray, y::IntegerArray, ks::NTuple{2,Integer}, wv::WeightVec) =
    counts(x, y, (1:ks[1], 1:ks[2]), wv)
counts(x::IntegerArray, y::IntegerArray, k::Integer) = counts(x, y, (1:k, 1:k))
counts(x::IntegerArray, y::IntegerArray, k::Integer, wv::WeightVec) =
    counts(x, y, (1:k, 1:k), wv)
counts(x::IntegerArray, y::IntegerArray) = counts(x, y, (span(x), span(y)))
counts(x::IntegerArray, y::IntegerArray, wv::WeightVec) = counts(x, y, (span(x), span(y)), wv)

proportions(x::IntegerArray, y::IntegerArray, levels::NTuple{2,IntUnitRange}) =
    counts(x, y, levels) .* inv(length(x))
proportions(x::IntegerArray, y::IntegerArray, levels::NTuple{2,IntUnitRange}, wv::WeightVec) =
    counts(x, y, levels, wv) .* inv(sum(wv))

proportions(x::IntegerArray, y::IntegerArray, ks::NTuple{2,Integer}) =
    proportions(x, y, (1:ks[1], 1:ks[2]))
proportions(x::IntegerArray, y::IntegerArray, ks::NTuple{2,Integer}, wv::WeightVec) =
    proportions(x, y, (1:ks[1], 1:ks[2]), wv)
proportions(x::IntegerArray, y::IntegerArray, k::Integer) = proportions(x, y, (1:k, 1:k))
proportions(x::IntegerArray, y::IntegerArray, k::Integer, wv::WeightVec) =
    proportions(x, y, (1:k, 1:k), wv)
proportions(x::IntegerArray, y::IntegerArray) = proportions(x, y, (span(x), span(y)))
proportions(x::IntegerArray, y::IntegerArray, wv::WeightVec) =
    proportions(x, y, (span(x), span(y)), wv)

#################################################
#
#  countmap on unknown levels
#
#  These methods are based on dictionaries, and
#  can be used on any kind of hashable values.
#
#################################################

## auxiliary functions

function _normalize_countmap{T}(cm::Dict{T}, s::Real)
    r = Dict{T,Float64}()
    for (k, c) in cm
        r[k] = c / s
    end
    return r
end

## 1D

function addcounts!{T}(cm::Dict{T}, x::AbstractArray{T})
    for v in x
        cm[v] = get(cm, v, 0) + 1
    end
    return cm
end

function addcounts!{T,W}(cm::Dict{T}, x::AbstractArray{T}, wv::WeightVec{W})
    n = length(x)
    length(wv) == n || throw(DimensionMismatch())
    w = values(wv)
    z = zero(W)

    for i = 1 : n
        @inbounds xi = x[i]
        @inbounds wi = w[i]
        cm[xi] = get(cm, xi, z) + wi
    end
    return cm
end

"""
    countmap(x)

Return a dictionary mapping each unique value in `x` to its number
of occurrences.
"""
countmap{T}(x::AbstractArray{T}) = addcounts!(Dict{T,Int}(), x)
countmap{T,W}(x::AbstractArray{T}, wv::WeightVec{W}) = addcounts!(Dict{T,W}(), x, wv)

"""
    proportionmap(x)

Return a dictionary mapping each unique value in `x` to its
proportion in `x`.
"""
proportionmap(x::AbstractArray) = _normalize_countmap(countmap(x), length(x))
proportionmap(x::AbstractArray, wv::WeightVec) = _normalize_countmap(countmap(x, wv), sum(wv))


text:
# Load Julia packages (libraries)

using Pkg, DrWatson

#@quickactivate "SR2StanPluto"
using StanSample
using StatisticalRethinking

# Define the Stan language model

stan8_4s = "
data{
  int N;
  vector[N] y;
}
parameters{
  real sigma;
  real alpha;
  real a1;
  real a2;
}
model{
  real mu;
  mu = a1 + a2;
  y ~ normal( mu , sigma );
}
";

# Define the Stanmodel and set the output format to :mcmcchains.

m8_4s = SampleModel("m8.4s", stan8_4s);

# Input data for cmdstan

m8_4_data = Dict("N" => 100, "y" => rand(Normal(0, 1), 100));

# Sample using cmdstan

rc8_4s = stan_sample(m8_4s, data=m8_4_data);

rethinking = "
        mean   sd  5.5% 94.5% n_eff Rhat
alpha 0.06 1.90 -2.22  2.49  1321    1
sigma 2.15 2.32  0.70  5.21   461    1
";

# Describe the draws

if success(rc)
  part8_4s = read_samples(m8_4s, :particles)
  part8_4s |> display
end


text:
export Point

struct Point
    latitude
    longitude
end

latitude(p::Point) = p.latitude
longitude(p::Point) = p.longitude


text:
# Particle in a box
using Plots
using LinearAlgebra: eigen
using HelFEM: HelFEM

a, b = -33.3, 7.1
pb = HelFEM.PolynomialBasis(:lip, 10)
B = HelFEM.FEMBasis(pb, range(a, b, length=21))
DD = HelFEM.radial_integral(B, q->1.0, lderivative=true, rderivative=true)
S = HelFEM.radial_integral(B, q->1.0)
e = eigen(DD, S)
# Compare energies with analytical energies (En = (πn/L)^2)
N=length(B)
es = (1:N) .^ 2 .* (π / (b - a))^2
e.values .- es
# Plot the wavefunctions -- they should be sinewaves:
xs = range(a, b, length=1001)
Bri = B(xs)
ns = 1:10
plot(xs, Bri*e.vectors[:,ns], alpha=(ns .^ -0.6)', legend=false)