text:
@testset "PowerSystems dict parsing" begin
    @info "parsing data from $RTS_GMLC_DIR into ps_dict"
    data_dict = PowerSystems.read_csv_data(RTS_GMLC_DIR, 100.0)
    @test haskey(data_dict, "timeseries_pointers")
end

@testset "CDM parsing" begin
    cdm_dict = nothing
    @info "parsing data from $RTS_GMLC_DIR into ps_dict"
    cdm_dict = PowerSystems.csv2ps_dict(RTS_GMLC_DIR, 100.0)
    @test cdm_dict isa Dict && haskey(cdm_dict, "loadzone")

    @info "assigning time series data for DA"
    cdm_dict = PowerSystems.assign_ts_data(cdm_dict, cdm_dict["forecast"]["DA"])
    @test length(cdm_dict["gen"]["Renewable"]["PV"]["102_PV_1"]["scalingfactor"]) == 24

    @info "making DA System"
    sys_rts_da = System(cdm_dict)
    @test sys_rts_da isa System

    @info "assigning time series data for RT"
    cdm_dict = PowerSystems.assign_ts_data(cdm_dict, cdm_dict["forecast"]["RT"])
    @test length(cdm_dict["gen"]["Renewable"]["PV"]["102_PV_1"]["scalingfactor"]) == 288

    @info "making RT System"
    sys_rts_rt = System(cdm_dict)
    @test sys_rts_rt isa System

    # Verify functionality of the concrete version of System.
    # TODO: Refactor once the ConcreteSystem implementation is finalized.
    sys = ConcreteSystem(sys_rts_da)
    @test length(sys_rts_rt.branches) == length(collect(get_mixed_components(Branch, sys)))
    @test length(sys_rts_rt.loads) == length(collect(get_mixed_components(ElectricLoad, sys)))
    @test length(sys_rts_rt.storage) == length(collect(get_mixed_components(Storage, sys)))
    @test length(sys_rts_rt.generators.thermal) == length(collect(get_mixed_components(ThermalGen, sys)))
    @test length(sys_rts_rt.generators.renewable) == length(collect(get_mixed_components(RenewableGen, sys)))
    @test length(sys_rts_rt.generators.hydro) == length(collect(get_mixed_components(HydroGen, sys)))
    @test length(get_components(Bus, sys)) > 0
    @test length(get_components(ThermalDispatch, sys)) > 0
    for x in (true, false)
        show_component_counts(sys, devnull; show_hierarchy=x)
    end
end

@testset "CDM parsing invalid directory" begin
    baddir = joinpath(RTS_GMLC_DIR, "../../test")
    @info "testing bad directory"
    @test_throws ErrorException PowerSystems.csv2ps_dict(baddir, 100.0)
end

@testset "consistency between CDM and standardfiles" begin
    mp_dict  = parsestandardfiles(joinpath(MATPOWER_DIR, "RTS_GMLC.m"))
    pm_dict = parse_file(joinpath(MATPOWER_DIR, "RTS_GMLC.m"))
    pmmp_dict = PowerSystems.pm2ps_dict(pm_dict)
    mpmmpsys = System(pmmp_dict)

    mpsys = System(mp_dict)

    cdm_dict = PowerSystems.csv2ps_dict(RTS_GMLC_DIR, 100.0)
    cdmsys = System(cdm_dict)

    @test cdmsys.generators.thermal[1].tech.activepowerlimits == mpsys.generators.thermal[1].tech.activepowerlimits
    @test cdmsys.generators.thermal[1].tech.reactivepowerlimits == mpsys.generators.thermal[1].tech.reactivepowerlimits
    @test_skip cdmsys.generators.thermal[1].tech.ramplimits == mpsys.generators.thermal[1].tech.ramplimits

    @test cdmsys.generators.thermal[1].econ.capacity == mpsys.generators.thermal[1].econ.capacity
    @test_skip cdmsys.generators.thermal[1].econ.variablecost == mpsys.generators.thermal[1].econ.variablecost

    @test cdmsys.generators.hydro[1].tech.activepowerlimits == mpsys.generators.hydro[1].tech.activepowerlimits
    @test cdmsys.generators.hydro[1].tech.reactivepowerlimits == mpsys.generators.hydro[1].tech.reactivepowerlimits
    @test cdmsys.generators.hydro[1].tech.installedcapacity == mpsys.generators.hydro[1].tech.installedcapacity
    @test_skip cdmsys.generators.hydro[1].tech.ramplimits == mpsys.generators.hydro[1].tech.ramplimits # this gets adjusted in the pm2ps_dict

    @test cdmsys.generators.hydro[1].econ == mpsys.generators.hydro[1].econ

    @test cdmsys.generators.renewable[1].tech == mpsys.generators.renewable[1].tech

    @test cdmsys.generators.renewable[1].econ == mpsys.generators.renewable[1].econ

    @test cdmsys.branches[1].rate ==
        [b for b in mpsys.branches if
            (b.connectionpoints.from.name == uppercase(cdmsys.branches[1].connectionpoints.from.name))
                & (b.connectionpoints.to.name == uppercase(cdmsys.branches[1].connectionpoints.to.name))][1].rate

    @test cdmsys.branches[6].rate ==
        [b for b in mpsys.branches if
            (b.connectionpoints.from.name == uppercase(cdmsys.branches[6].connectionpoints.from.name))
                & (b.connectionpoints.to.name == uppercase(cdmsys.branches[6].connectionpoints.to.name))][1].rate

    @test cdmsys.branches[120].rate == [b for b in mpsys.branches if
            (b.connectionpoints.from.name == uppercase(cdmsys.branches[120].connectionpoints.from.name))
                & (b.connectionpoints.to.name == uppercase(cdmsys.branches[120].connectionpoints.to.name))][1].rate

end


text:
using BeefBLUP
using DataFrames
using Test

@testset "BeefBLUP.jl" begin
    # Write your tests here.
    correctX = [1 1 0 1; 1 1 0 0; 1 0 1 1; 1 0 1 0; 1 0 1 1; 1 0 1 0; 1 0 0 1]
    fixedfx = DataFrame(year = [1990, 1990, 1991, 1991, 1991, 1991, 1992], sex = ["male", "female", "male", "female", "male", "female", "male"])
    @test BeefBLUP.fixedeffectmatrix(fixedfx)[1] == correctX
    correctA = [1   0   1/2 1/2 1/2 0   0;
                0   1   0   0   1/2 1/2 0;
                1/2 0   1   1/4 1/4 0   0;
                1/2 0   1/4 1   1/4 0   0;
                1/2 1/2 1/4 1/4 1   1/4 0;
                0 1/2   0   0   1/4 1   0;
                0 0     0   0   0   0   1]
    id = collect(1:7)
    dam_id = [missing, missing, missing, missing, 2, 2, missing]
    sire_id = [missing, missing, 1, 1, 1, missing, missing]
    @test BeefBLUP.additiverelationshipmatrix(id, dam_id, sire_id) == correctA
end

ext: jl

max_stars_repo_name: MillironX/beefblup

max_issues_repo_name: MillironX/beefblup

max_forks_repo_name: MillironX/beefblup


text:
import NCDatasets
filename = tempname()

lons = 9.:30.
lats = 54.:66.
times = 0:10

# generate a NetCDF file

NCDatasets.NCDataset(filename,"c") do ds
    NCDatasets.defDim(ds,"lonc",length(lons))
    NCDatasets.defDim(ds,"latc",length(lats))
    NCDatasets.defDim(ds,"time",Inf)

    tempvar = NCDatasets.defVar(ds,"temp",Float32,("lonc","latc","time"))
    timevar = NCDatasets.defVar(ds,"time",Float64,("time",))

    for itime = 1:length(times)
        timevar[itime] = itime * 3600.0
    end
end

# append data

NCDatasets.NCDataset(filename,"a") do ds
    tempvar = ds["temp"]

    for itime = 1:length(times)
        #@show size(tempvar)
        #@show size(tempvar[:,:,itime])
        data = [i+j for i = 1:length(lons), j = 1:length(lats)]
        #tempvar[:,:,itime] = data
        tempvar[:,:,itime] = data[:,:,1:1]
    end
end

rm(filename)


text:
# This file is a part of Julia. Except for the *_kernel functions (see below),
# license is MIT: https://julialang.org/license

struct DoubleFloat64
    hi::Float64
    lo::Float64
end
struct DoubleFloat32
    hi::Float64
end

# sin_kernel and cos_kernel functions are only valid for |x| < pi/4 = 0.7854
# translated from openlibm code: k_sin.c, k_cos.c, k_sinf.c, k_cosf.c.
# acos functions are based on openlibm code: e_acos.c, e_acosf.c.
# asin functions are based on openlibm code: e_asin.c, e_asinf.c. The above
# functions are made available under the following licence:

## Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
##
## Developed at SunPro, a Sun Microsystems, Inc. business.
## Permission to use, copy, modify, and distribute this
## software is freely granted, provided that this notice
## is preserved.

# Trigonometric functions
# sin methods
@noinline sin_domain_error(x) = throw(DomainError(x, "sin(x) is only defined for finite x."))
function sin(x::T) where T<:Union{Float32, Float64}
    absx = abs(x)
    if absx < T(pi)/4 #|x| ~<= pi/4, no need for reduction
        if absx < sqrt(eps(T))
            return x
        end
        return sin_kernel(x)
    elseif isnan(x)
        return T(NaN)
    elseif isinf(x)
        sin_domain_error(x)
    end
    n, y = rem_pio2_kernel(x)
    n = n&3
    if n == 0
        return sin_kernel(y)
    elseif n == 1
        return cos_kernel(y)
    elseif n == 2
        return -sin_kernel(y)
    else
        return -cos_kernel(y)
    end
end
sin(x::Real) = sin(float(x))

# Coefficients in 13th order polynomial approximation on [0; π/4]
#     sin(x) ≈ x + S1*x³ + S2*x⁵ + S3*x⁷ + S4*x⁹ + S5*x¹¹ + S6*x¹³
# D for double, S for sin, number is the order of x-1

const DS1 = -1.66666666666666324348e-01
const DS2 = 8.33333333332248946124e-03
const DS3 = -1.98412698298579493134e-04
const DS4 = 2.75573137070700676789e-06
const DS5 = -2.50507602534068634195e-08
const DS6 = 1.58969099521155010221e-10
"""
    sin_kernel(yhi, ylo)

Computes the sine on the interval [-π/4; π/4].
"""
@inline function sin_kernel(y::DoubleFloat64)
    y² = y.hi*y.hi
    y⁴ = y²*y²
    r  = @horner(y², DS2, DS3, DS4) + y²*y⁴*@horner(y², DS5, DS6)
    y³ = y²*y.hi
    y.hi-((y²*(0.5*y.lo-y³*r)-y.lo)-y³*DS1)
end
@inline function sin_kernel(y::Float64)
    y² =  y*y
    y⁴ =  y²*y²
    r  =  @horner(y², DS2, DS3, DS4) + y²*y⁴*@horner(y², DS5, DS6)
    y³ =  y²*y
    y+y³*(DS1+y²*r)
end

# sin_kernels accepting values from rem_pio2 in the Float32 case
@inline sin_kernel(x::Float32) = sin_kernel(DoubleFloat32(x))
@inline function sin_kernel(y::DoubleFloat32)
    S1 = -0.16666666641626524
    S2 = 0.008333329385889463
    z = y.hi*y.hi
    w = z*z
    r = @horner(z, -0.00019839334836096632, 2.718311493989822e-6)
    s = z*y.hi
    Float32((y.hi + s*@horner(z, S1, S2)) + s*w*r)
end

# cos methods
@noinline cos_domain_error(x) = throw(DomainError(x, "cos(x) is only defined for finite x."))
function cos(x::T) where T<:Union{Float32, Float64}
    absx = abs(x)
    if absx < T(pi)/4
        if absx < sqrt(eps(T)/T(2.0))
            return T(1.0)
        end
        return cos_kernel(x)
    elseif isnan(x)
        return T(NaN)
    elseif isinf(x)
        cos_domain_error(x)
    else
        n, y = rem_pio2_kernel(x)
        n = n&3
        if n == 0
            return cos_kernel(y)
        elseif n == 1
            return -sin_kernel(y)
        elseif n == 2
            return -cos_kernel(y)
        else
            return sin_kernel(y)
        end
    end
end
cos(x::Real) = cos(float(x))

const DC1 = 4.16666666666666019037e-02
const DC2 = -1.38888888888741095749e-03
const DC3 = 2.48015872894767294178e-05
const DC4 = -2.75573143513906633035e-07
const DC5 = 2.08757232129817482790e-09
const DC6 = -1.13596475577881948265e-11
"""
    cos_kernel(y)

Compute the cosine on the interval y∈[-π/4; π/4].
"""
@inline function cos_kernel(y::DoubleFloat64)
    y² = y.hi*y.hi
    y⁴ = y²*y²
    r  = y²*@horner(y², DC1, DC2, DC3) + y⁴*y⁴*@horner(y², DC4, DC5, DC6)
    half_y² = 0.5*y²
    w  = 1.0-half_y²
    w + (((1.0-w)-half_y²) + (y²*r-y.hi*y.lo))
end
@inline function cos_kernel(y::Float64)
    y² = y*y
    y⁴ = y²*y²
    r  = y²*@horner(y², DC1, DC2, DC3) + y⁴*y⁴*@horner(y², DC4, DC5, DC6)
    half_y² = 0.5*y²
    w  = 1.0-half_y²
    w + (((1.0-w)-half_y²) + (y²*r))
end

# cos_kernels accepting values from rem_pio2 in the Float32 case
cos_kernel(x::Float32) = cos_kernel(DoubleFloat32(x))
@inline function cos_kernel(y::DoubleFloat32)
    C0 = -0.499999997251031
    C1 = 0.04166662332373906
    y² = y.hi*y.hi
    y⁴ = y²*y²
    r = @horner(y², -0.001388676377460993, 2.439044879627741e-5)
    Float32(((1.0+y²*C0) + y⁴*C1) + (y⁴*y²)*r)
end

### sincos methods
@noinline sincos_domain_error(x) = throw(DomainError(x, "sincos(x) is only defined for finite x."))

"""
    sincos(x)

Simultaneously compute the sine and cosine of `x`, where the `x` is in radians.
"""
function sincos(x::T) where T<:Union{Float32, Float64}
    if abs(x) < T(pi)/4
        if x == zero(T)
            return x, one(T)
        end
        return sincos_kernel(x)
    elseif isnan(x)
        return T(NaN), T(NaN)
    elseif isinf(x)
        sincos_domain_error(x)
    end
    n, y = rem_pio2_kernel(x)
    n = n&3
    # calculate both kernels at the reduced y...
    si, co = sincos_kernel(y)
    # ... and use the same selection scheme as above: (sin, cos, -sin, -cos) for
    # for sin and (cos, -sin, -cos, sin) for cos
    if n == 0
        return si, co
    elseif n == 1
        return co, -si
    elseif n == 2
        return -si, -co
    else
        return -co, si
    end
end
sincos(x::Real) = sincos(float(x))

# There's no need to write specialized kernels, as inlining takes care of remo-
# ving superfluous calculations.
@inline sincos_kernel(y::Union{Float32, Float64, DoubleFloat32, DoubleFloat64}) = (sin_kernel(y), cos_kernel(y))

# tangent methods
@noinline tan_domain_error(x) = throw(DomainError(x, "tan(x) is only defined for finite x."))
function tan(x::T) where T<:Union{Float32, Float64}
    absx = abs(x)
    if absx < T(pi)/4
        if absx < sqrt(eps(T))/2 # first order dominates, but also allows tan(-0)=-0
            return x
        end
        return tan_kernel(x)
    elseif isnan(x)
        return T(NaN)
    elseif isinf(x)
        tan_domain_error(x)
    end
    n, y = rem_pio2_kernel(x)
    if iseven(n)
        return tan_kernel(y,1)
    else
        return tan_kernel(y,-1)
    end
end
tan(x::Real) = tan(float(x))

@inline tan_kernel(y::Float64) = tan_kernel(DoubleFloat64(y, 0.0), 1)
@inline function tan_kernel(y::DoubleFloat64, k)
    # kernel tan function on ~[-pi/4, pi/4] (except on -0)
    # Input y is assumed to be bounded by ~pi/4 in magnitude.
    # Input k indicates whether tan (if k = 1) or -1/tan (if k = -1) is returned.

    # Algorithm
    #    1. Since tan(-y) = -tan(y), we need only to consider positive y.
    #    2. Callers must return tan(-0) = -0 without calling here since our
    #       odd polynomial is not evaluated in a way that preserves -0.
    #       Callers may do the optimization tan(y) ~ y for tiny y.
    #    3. tan(y) is approximated by a odd polynomial of degree 27 on
    #       [0,0.67434]
    #                       3             27
    #           tan(y) ~ y + T1*y + ... + T13*y ≡ P(y)
    #       where
    #
    #                          |tan(y)         2     4            26   |     -59.2
    #        (tan(y)-P(y))/y = |----- - (1+T1*y +T2*y +.... +T13*y    )| <= 2
    #                          |  y                                    |
    #
    #       Note: tan(y+z) = tan(y) + tan'(y)*z
    #                  ~ tan(y) + (1+y*y)*z
    #       Therefore, for better accuracz in computing tan(y+z), let
    #             3      2      2       2       2
    #        r = y *(T2+y *(T3+y *(...+y *(T12+y *T13))))
    #       then
    #                     3    2
    #        tan(y+z) = y + (T1*y + (y *(r+z)+z))
    #
    #   4. For y in [0.67434,pi/4],  let z = pi/4 - y, then
    #        tan(y) = tan(pi/4-z) = (1-tan(z))/(1+tan(z))
    #               = 1 - 2*(tan(z) - (tan(z)^2)/(1+tan(z)))

    yhi = y.hi
    ylo = y.lo

    if abs(yhi) >= 0.6744
        if yhi < 0.0
            yhi = -yhi
            ylo = -ylo
        end
        # Then, accurately reduce y as "pio4hi"-yhi+"pio4lo"-ylo
        yhi = (pi/4 - yhi) + (3.06161699786838301793e-17 - ylo)
        # yhi is guaranteed to be exact, so ylo is identically zero
        ylo = 0.0
    end
    y² = yhi * yhi
    y⁴ = y² * y²

    # Break P(y)-T1*y³ = y^5*(T[2]+y^2*T[3]+...) into y⁵*r + y⁵*v where
    # r = T[2]+y^4*T[4]+...+y^20*T[12])
    # v = (y^2*(T[3]+y^4*T[5]+...+y^22*[T13]))
    r = @horner(y⁴,
        1.33333333333201242699e-01, # T2
        2.18694882948595424599e-02, # T4
        3.59207910759131235356e-03, # T6
        5.88041240820264096874e-04, # T8
        7.81794442939557092300e-05, # T10
        -1.85586374855275456654e-05) # T12
    v = y² * @horner(y⁴,
        5.39682539762260521377e-02, # T3
        8.86323982359930005737e-03, # T5
        1.45620945432529025516e-03, # T7
        2.46463134818469906812e-04, # T9
        7.14072491382608190305e-05, # T11
        2.59073051863633712884e-05) # T13
    # Precompute y³
    y³ = y² * yhi
    # Calculate  P(y)-y-T1*y³ =  y⁵*r + y⁵*v  = y²(y³*(r+v))
    r = ylo + y² * (y³ * (r + v) + ylo)
    # Calculate P(y)-y = r+T1*y³
    r += 3.33333333333334091986e-01*y³
    # Calculate w = r+y = P(y)
    Px = yhi + r
    if abs(y.hi) >= 0.6744
        # If the original y was above the threshold, then we calculate
        #     tan(y) = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
        #            ≈ 1 - 2*(P(z) - (P(z)^2)/(1+P(z)))
        # where z = y-π/4.
        return (signbit(y.hi) ? -1.0 : 1.0)*(k - 2*(yhi-(Px^2/(k+Px)-r)))
    end
    if k == 1
        # Else, we simply return w = P(y) if k == 1 (integer multiple from argument
        # reduction was even)...
        return Px
    else
        # ...or tan(y) ≈ -1.0/(y+r) if !(k == 1) (integer multiple from argument
        # reduction was odd). If 2ulp error is allowed, simply return the frac-
        # tion directly. Instead, we calculate it accurately.

        # Px0 is w with zeroed out low word
        Px0 = reinterpret(Float64, (reinterpret(UInt64, Px) >> 32) << 32)
        v = r - (Px0 - yhi) # Px0+v = r+y
        t = a = -1.0 / Px
        # zero out low word of t
        t = reinterpret(Float64, (reinterpret(UInt64, t) >> 32) << 32)
        s = 1.0 + t * Px0
        return t + a * (s + t * v)
    end
end

@inline tan_kernel(y::Float32) = tan_kernel(DoubleFloat32(y), 1)
@inline function tan_kernel(y::DoubleFloat32, k)
    # |tan(y)/y - t(y)| < 2**-25.5 (~[-2e-08, 2e-08]). */
    y² = y.hi*y.hi
    r  = @horner(y², 0.00297435743359967304927, 0.00946564784943673166728)
    t  = @horner(y², 0.0533812378445670393523, 0.0245283181166547278873)
    y⁴ = y²*y²
    y³ = y²*y.hi
    u  = @horner(y², 0.333331395030791399758, 0.133392002712976742718)
    Py  = (y.hi+y³*u)+(y³*y⁴)*(t+y⁴*r)
    if k == 1
        return Float32(Py)
    end

    return Float32(-1.0/Py)
end

# fallback methods
sin_kernel(x::Real) = sin(x)
cos_kernel(x::Real) = cos(x)
tan_kernel(x::Real) = tan(x)
sincos_kernel(x::Real) = sincos(x)

# Inverse trigonometric functions
# asin methods
ASIN_X_MIN_THRESHOLD(::Type{Float32}) = 2.0f0^-12
ASIN_X_MIN_THRESHOLD(::Type{Float64}) = sqrt(eps(Float64))

arc_p(t::Float64) =
    t*@horner(t,
    1.66666666666666657415e-01,
    -3.25565818622400915405e-01,
    2.01212532134862925881e-01,
    -4.00555345006794114027e-02,
    7.91534994289814532176e-04,
    3.47933107596021167570e-05)

arc_q(z::Float64) =
    @horner(z,
    1.0,
    -2.40339491173441421878e+00,
    2.02094576023350569471e+00,
    -6.88283971605453293030e-01,
    7.70381505559019352791e-02)

arc_p(t::Float32) =
    t*@horner(t,
    1.6666586697f-01,
    -4.2743422091f-02,
    -8.6563630030f-03)

arc_q(t::Float32) = @horner(t, 1.0f0, -7.0662963390f-01)

@inline arc_tRt(t) = arc_p(t)/arc_q(t)

@inline function asin_kernel(t::Float64, x::Float64)
    # we use that for 1/2 <= x < 1 we have
    #     asin(x) = pi/2-2*asin(sqrt((1-x)/2))
    # Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
    # then for x>0.98
    #     asin(x) = pi/2 - 2*(s+s*z*R(z))
    #         = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
    # For x<=0.98, let pio4_hi = pio2_hi/2, then
    #     f = hi part of s;
    #     c = sqrt(z) - f = (z-f*f)/(s+f)     ...f+c=sqrt(z)
    #  and
    #     asin(x) = pi/2 - 2*(s+s*z*R(z))
    #         = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
    #         = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
    pio2_lo = 6.12323399573676603587e-17
    s = sqrt_llvm(t)
    tRt = arc_tRt(t)
    if abs(x) >= 0.975 # |x| > 0.975
        return flipsign(pi/2 - (2.0*(s + s*tRt) - pio2_lo), x)
    else
        s0 = reinterpret(Float64, (reinterpret(UInt64, s) >> 32) << 32)
        c = (t - s0*s0)/(s + s0)
        p = 2.0*s*tRt - (pio2_lo - 2.0*c)
        q = pi/4 - 2.0*s0
        return flipsign(pi/4 - (p-q), x)
    end
end
@inline function asin_kernel(t::Float32, x::Float32)
    s = sqrt_llvm(Float64(t))
    tRt = arc_tRt(t) # rational approximation
    flipsign(Float32(pi/2 - 2*(s + s*tRt)), x)
end

@noinline asin_domain_error(x) = throw(DomainError(x, "asin(x) is not defined for |x|>1."))
function asin(x::T) where T<:Union{Float32, Float64}
    # Since  asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
    # we approximate asin(x) on [0,0.5] by
    #     asin(x) = x + x*x^2*R(x^2)
    # where
    #     R(x^2) is a rational approximation of (asin(x)-x)/x^3
    # and its remez error is bounded by
    #     |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75)
    absx = abs(x)
    if absx >= T(1.0) # |x|>= 1
        if absx == T(1.0)
            return flipsign(T(pi)/2, x)
        end
        asin_domain_error(x)
    elseif absx < T(1.0)/2
        # if |x| sufficiently small, |x| is a good approximation
        if absx < ASIN_X_MIN_THRESHOLD(T)
            return x
        end
        return muladd(x, arc_tRt(x*x), x)
    end
    # else 1/2 <= |x| < 1
    t = (T(1.0) - absx)/2
    return asin_kernel(t, x)
end
asin(x::Real) = asin(float(x))

# atan2 methods
ATAN2_PI_LO(::Type{Float32}) = -8.7422776573f-08
ATAN2_RATIO_BIT_SHIFT(::Type{Float32}) = 23
ATAN2_RATIO_THRESHOLD(::Type{Float32}) = 26

ATAN2_PI_LO(::Type{Float64}) = 1.2246467991473531772E-16
ATAN2_RATIO_BIT_SHIFT(::Type{Float64}) = 20
ATAN2_RATIO_THRESHOLD(::Type{Float64}) = 60

function atan2(y::T, x::T) where T<:Union{Float32, Float64}
    # Method :
    #    M1) Reduce y to positive by atan2(y,x)=-atan2(-y,x).
    #    M2) Reduce x to positive by (if x and y are unexceptional):
    #        ARG (x+iy) = arctan(y/x)          ... if x > 0,
    #        ARG (x+iy) = pi - arctan[y/(-x)]   ... if x < 0,
    #
    # Special cases:
    #
    #    S1) ATAN2((anything), NaN ) is NaN;
    #    S2) ATAN2(NAN , (anything) ) is NaN;
    #    S3) ATAN2(+-0, +(anything but NaN)) is +-0  ;
    #    S4) ATAN2(+-0, -(anything but NaN)) is +-pi ;
    #    S5) ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
    #    S6) ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
    #    S7) ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
    #    S8) ATAN2(+-INF,+INF ) is +-pi/4 ;
    #    S9) ATAN2(+-INF,-INF ) is +-3pi/4;
    #    S10) ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
    if isnan(x) || isnan(y) # S1 or S2
        return T(NaN)
    end

    if x == T(1.0) # then y/x = y and x > 0, see M2
        return atan(y)
    end
    # generate an m ∈ {0, 1, 2, 3} to branch off of
    m = 2*signbit(x) + 1*signbit(y)

    if iszero(y)
        if m == 0 || m == 1
            return y # atan(+-0, +anything) = +-0
        elseif m == 2
            return T(pi) # atan(+0, -anything) = pi
        elseif m == 3
            return -T(pi) # atan(-0, -anything) =-pi
        end
    elseif iszero(x)
        return flipsign(T(pi)/2, y)
    end

    if isinf(x)
        if isinf(y)
            if m == 0
                return T(pi)/4  # atan(+Inf), +Inf))
            elseif m == 1
                return -T(pi)/4 # atan(-Inf), +Inf))
            elseif m == 2
                return 3*T(pi)/4 # atan(+Inf, -Inf)
            elseif m == 3
                return -3*T(pi)/4 # atan(-Inf,-Inf)
            end
        else
            if m == 0
                return zero(T)  # atan(+...,+Inf) */
            elseif m == 1
                return -zero(T) # atan(-...,+Inf) */
            elseif m == 2
                return T(pi)    # atan(+...,-Inf) */
            elseif m == 3
                return -T(pi)   # atan(-...,-Inf) */
            end
        end
    end

    # x wasn't Inf, but y is
    isinf(y) && return copysign(T(pi)/2, y)

    ypw = poshighword(y)
    xpw = poshighword(x)
    # compute y/x for Float32
    k = reinterpret(Int32, ypw-xpw)>>ATAN2_RATIO_BIT_SHIFT(T)

    if k > ATAN2_RATIO_THRESHOLD(T) # |y/x| >  threshold
        z=T(pi)/2+T(0.5)*ATAN2_PI_LO(T)
        m&=1;
    elseif x<0 && k < -ATAN2_RATIO_THRESHOLD(T) # 0 > |y|/x > threshold
        z = zero(T)
    else #safe to do y/x
        z = atan(abs(y/x))
    end

    if m == 0
        return z # atan(+,+)
    elseif m == 1
        return -z # atan(-,+)
    elseif m == 2
        return T(pi)-(z-ATAN2_PI_LO(T)) # atan(+,-)
    else # default case m == 3
        return (z-ATAN2_PI_LO(T))-T(pi) # atan(-,-)
    end
end
# acos methods
ACOS_X_MIN_THRESHOLD(::Type{Float32}) = 2.0f0^-26
ACOS_X_MIN_THRESHOLD(::Type{Float64}) = 2.0^-57
PIO2_HI(::Type{Float32}) = 1.5707962513f+00
PIO2_LO(::Type{Float32}) = 7.5497894159f-08
PIO2_HI(::Type{Float64}) = 1.57079632679489655800e+00
PIO2_LO(::Type{Float64}) = 6.12323399573676603587e-17
ACOS_PI(::Type{Float32}) = 3.1415925026f+00
ACOS_PI(::Type{Float64}) = 3.14159265358979311600e+00
@inline ACOS_CORRECT_LOWWORD(::Type{Float32}, x) = reinterpret(Float32, (reinterpret(UInt32, x) & 0xfffff000))
@inline ACOS_CORRECT_LOWWORD(::Type{Float64}, x) = reinterpret(Float64, (reinterpret(UInt64, x) >> 32) << 32)

@noinline acos_domain_error(x) = throw(DomainError(x, "acos(x) not defined for |x| > 1"))
function acos(x::T) where T <: Union{Float32, Float64}
    # Method :
    #    acos(x)  = pi/2 - asin(x)
    #    acos(-x) = pi/2 + asin(x)
    # As a result, we use the same rational approximation (arc_tRt) as in asin.
    # See the comments in asin for more information about this approximation.
    # 1) For |x| <= 0.5
    #    acos(x) = pi/2 - (x + x*x^2*R(x^2))
    # 2) For x < -0.5
    #    acos(x) = pi - 2asin(sqrt((1 - |x|)/2))
    #        = pi - 0.5*(s+s*z*R(z))
    # where z=(1-|x|)/2, s=sqrt(z)
    # 3) For x > 0.5
    #     acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1 - x)/2)))
    #        = 2asin(sqrt((1 - x)/2))
    #        = 2s + 2s*z*R(z)     ...z=(1 - x)/2, s=sqrt(z)
    #        = 2f + (2c + 2s*z*R(z))
    #    where f=hi part of s, and c = (z - f*f)/(s + f) is the correction term
    #    for f so that f + c ~ sqrt(z).

    # Special cases:
    #    4) if x is NaN, return x itself;
    #    5) if |x|>1 throw warning.

    absx = abs(x)
    if absx >= T(1.0)
        # acos(-1) = π, acos(1) = 0
        absx == T(1.0) && return x > T(0.0) ? T(0.0) : T(pi)
        # acos(x) is not defined for |x| > 1
        acos_domain_error(x) # see 5) above
    elseif absx < T(1.0)/2 # see 1) above
        # if |x| sufficiently small, acos(x) ≈ pi/2
        absx < ACOS_X_MIN_THRESHOLD(T) && return T(pi)/2
        # if |x| < 0.5 we have acos(x) = pi/2 - (x + x*x^2*R(x^2))
        return PIO2_HI(T) - (x - (PIO2_LO(T) - x*arc_tRt(x*x)))
    end
    z = (T(1.0) - absx)*T(0.5)
    zRz = arc_tRt(z)
    s = sqrt_llvm(z)
    if x < T(0.0) # see 2) above
        return ACOS_PI(T) - T(2.0)*(s + (zRz*s - PIO2_LO(T)))
    else # see 3) above
        # if x > 0.5 we have
        # acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2)))
        #         = 2asin(sqrt((1-x)/2))
        #         = 2s + 2s*z*R(z)    ...z=(1-x)/2, s=sqrt(z)
        #         = 2f + (2c + 2s*z*R(z))
        # where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term
        # for f so that f+c ~ sqrt(z).
        df = ACOS_CORRECT_LOWWORD(T, s)
        c  = (z - df*df)/(s + df)
        return T(2.0)*(df + (zRz*s + c))
    end
end
acos(x::Real) = acos(float(x))

# multiply in extended precision
function mulpi_ext(x::Float64)
    m = 3.141592653589793
    m_hi = 3.1415926218032837
    m_lo = 3.178650954705639e-8

    x_hi = reinterpret(Float64, reinterpret(UInt64,x) & 0xffff_ffff_f800_0000)
    x_lo = x-x_hi

    y_hi = m*x
    y_lo = x_hi * m_lo + (x_lo* m_hi + ((x_hi*m_hi-y_hi) + x_lo*m_lo))

    DoubleFloat64(y_hi,y_lo)
end
mulpi_ext(x::Float32) = DoubleFloat32(pi*Float64(x))
mulpi_ext(x::Rational) = mulpi_ext(float(x))
mulpi_ext(x::Real) = pi*x # Fallback

"""
    sinpi(x)

Compute ``\\sin(\\pi x)`` more accurately than `sin(pi*x)`, especially for large `x`.
"""
function sinpi(x::T) where T<:AbstractFloat
    if !isfinite(x)
        isnan(x) && return x
        throw(DomainError(x, "`x` cannot be infinite."))
    end

    ax = abs(x)
    s = maxintfloat(T)/2
    ax >= s && return copysign(zero(T),x) # integer-valued

    # reduce to interval [-1,1]
    # assumes RoundNearest rounding mode
    t = 3*s
    rx = x-((x+t)-t) # zeros may be incorrectly signed
    arx = abs(rx)

    if (arx == 0) | (arx == 1)
        copysign(zero(T),x)
    elseif arx < 0.25
        sin_kernel(mulpi_ext(rx))
    elseif arx < 0.75
        y = mulpi_ext(T(0.5) - arx)
        copysign(cos_kernel(y),rx)
    else
        y = mulpi_ext(copysign(one(T),rx) - rx)
        sin_kernel(y)
    end
end

# Integers and Rationals
function sinpi(x::T) where T<:Union{Integer,Rational}
    Tf = float(T)
    if !isfinite(x)
        throw(DomainError(x, "`x` must be finite."))
    end

    # until we get an IEEE remainder function (#9283)
    rx = rem(x,2)
    if rx > 1
        rx -= 2
    elseif rx < -1
        rx += 2
    end
    arx = abs(rx)

    if (arx == 0) | (arx == 1)
        copysign(zero(Tf),x)
    elseif arx < 0.25
        sin_kernel(mulpi_ext(rx))
    elseif arx < 0.75
        y = mulpi_ext(T(0.5) - arx)
        copysign(cos_kernel(y),rx)
    else
        y = mulpi_ext(copysign(one(T),rx) - rx)
        sin_kernel(y)
    end
end

"""
    cospi(x)

Compute ``\\cos(\\pi x)`` more accurately than `cos(pi*x)`, especially for large `x`.
"""
function cospi(x::T) where T<:AbstractFloat
    if !isfinite(x)
        isnan(x) && return x
        throw(DomainError(x, "`x` cannot be infinite."))
    end

    ax = abs(x)
    s = maxintfloat(T)
    ax >= s && return one(T) # even integer-valued

    # reduce to interval [-1,1], then [0,1]
    # assumes RoundNearest rounding mode
    rx = abs(ax-((ax+s)-s))

    if rx <= 0.25
        cos_kernel(mulpi_ext(rx))
    elseif rx < 0.75
        y = mulpi_ext(T(0.5) - rx)
        sin_kernel(y)
    else
        y = mulpi_ext(one(T) - rx)
        -cos_kernel(y)
    end
end

# Integers and Rationals
function cospi(x::T) where T<:Union{Integer,Rational}
    if !isfinite(x)
        throw(DomainError(x, "`x` must be finite."))
    end

    ax = abs(x)
    # until we get an IEEE remainder function (#9283)
    rx = rem(ax,2)
    if rx > 1
        rx = 2-rx
    end

    if rx <= 0.25
        cos_kernel(mulpi_ext(rx))
    elseif rx < 0.75
        y = mulpi_ext(T(0.5) - rx)
        sin_kernel(y)
    else
        y = mulpi_ext(one(T) - rx)
        -cos_kernel(y)
    end
end

sinpi(x::Integer) = x >= 0 ? zero(float(x)) : -zero(float(x))
cospi(x::Integer) = isodd(x) ? -one(float(x)) : one(float(x))
sinpi(x::Real) = sinpi(float(x))
cospi(x::Real) = cospi(float(x))

function sinpi(z::Complex{T}) where T
    F = float(T)
    zr, zi = reim(z)
    if isinteger(zr)
        # zr = ...,-2,-1,0,1,2,...
        # sin(pi*zr) == ±0
        # cos(pi*zr) == ±1
        # cosh(pi*zi) > 0
        s = copysign(zero(F),zr)
        c_pos = isa(zr,Integer) ? iseven(zr) : isinteger(zr/2)
        sh = sinh(pi*zi)
        Complex(s, c_pos ? sh : -sh)
    elseif isinteger(2*zr)
        # zr = ...,-1.5,-0.5,0.5,1.5,2.5,...
        # sin(pi*zr) == ±1
        # cos(pi*zr) == +0
        # sign(sinh(pi*zi)) == sign(zi)
        s_pos = isinteger((2*zr-1)/4)
        ch = cosh(pi*zi)
        Complex(s_pos ? ch : -ch, isnan(zi) ? zero(F) : copysign(zero(F),zi))
    elseif !isfinite(zr)
        if zi == 0 || isinf(zi)
            Complex(F(NaN), F(zi))
        else
            Complex(F(NaN), F(NaN))
        end
    else
        pizi = pi*zi
        Complex(sinpi(zr)*cosh(pizi), cospi(zr)*sinh(pizi))
    end
end

function cospi(z::Complex{T}) where T
    F = float(T)
    zr, zi = reim(z)
    if isinteger(zr)
        # zr = ...,-2,-1,0,1,2,...
        # sin(pi*zr) == ±0
        # cos(pi*zr) == ±1
        # sign(sinh(pi*zi)) == sign(zi)
        # cosh(pi*zi) > 0
        s = copysign(zero(F),zr)
        c_pos = isa(zr,Integer) ? iseven(zr) : isinteger(zr/2)
        ch = cosh(pi*zi)
        Complex(c_pos ? ch : -ch, isnan(zi) ? s : -flipsign(s,zi))
    elseif isinteger(2*zr)
        # zr = ...,-1.5,-0.5,0.5,1.5,2.5,...
        # sin(pi*zr) == ±1
        # cos(pi*zr) == +0
        # sign(sinh(pi*zi)) == sign(zi)
        s_pos = isinteger((2*zr-1)/4)
        sh = sinh(pi*zi)
        Complex(zero(F), s_pos ? -sh : sh)
    elseif !isfinite(zr)
        if zi == 0
            Complex(F(NaN), isnan(zr) ? zero(F) : -flipsign(F(zi),zr))
        elseif isinf(zi)
            Complex(F(Inf), F(NaN))
        else
            Complex(F(NaN), F(NaN))
        end
    else
        pizi = pi*zi
        Complex(cospi(zr)*cosh(pizi), -sinpi(zr)*sinh(pizi))
    end
end

"""
    sinc(x)

Compute ``\\sin(\\pi x) / (\\pi x)`` if ``x \\neq 0``, and ``1`` if ``x = 0``.
"""
sinc(x::Number) = x==0 ? one(x)  : oftype(x,sinpi(x)/(pi*x))
sinc(x::Integer) = x==0 ? one(x) : zero(x)
sinc(x::Complex{<:AbstractFloat}) = x==0 ? one(x) : oftype(x, sinpi(x)/(pi*x))
sinc(x::Complex) = sinc(float(x))
sinc(x::Real) = x==0 ? one(x) : isinf(x) ? zero(x) : sinpi(x)/(pi*x)

"""
    cosc(x)

Compute ``\\cos(\\pi x) / x - \\sin(\\pi x) / (\\pi x^2)`` if ``x \\neq 0``, and ``0`` if
``x = 0``. This is the derivative of `sinc(x)`.
"""
cosc(x::Number) = x==0 ? zero(x) : oftype(x,(cospi(x)-sinpi(x)/(pi*x))/x)
cosc(x::Integer) = cosc(float(x))
cosc(x::Complex{<:AbstractFloat}) = x==0 ? zero(x) : oftype(x,(cospi(x)-sinpi(x)/(pi*x))/x)
cosc(x::Complex) = cosc(float(x))
cosc(x::Real) = x==0 || isinf(x) ? zero(x) : (cospi(x)-sinpi(x)/(pi*x))/x

for (finv, f, finvh, fh, finvd, fd, fn) in ((:sec, :cos, :sech, :cosh, :secd, :cosd, "secant"),
                                            (:csc, :sin, :csch, :sinh, :cscd, :sind, "cosecant"),
                                            (:cot, :tan, :coth, :tanh, :cotd, :tand, "cotangent"))
    name = string(finv)
    hname = string(finvh)
    dname = string(finvd)
    @eval begin
        @doc """
            $($name)(x)

        Compute the $($fn) of `x`, where `x` is in radians.
        """ ($finv)(z::T) where {T<:Number} = one(T) / (($f)(z))
        @doc """
            $($hname)(x)

        Compute the hyperbolic $($fn) of `x`.
        """ ($finvh)(z::T) where {T<:Number} = one(T) / (($fh)(z))
        @doc """
            $($dname)(x)

        Compute the $($fn) of `x`, where `x` is in degrees.
        """ ($finvd)(z::T) where {T<:Number} = one(T) / (($fd)(z))
    end
end

for (tfa, tfainv, hfa, hfainv, fn) in ((:asec, :acos, :asech, :acosh, "secant"),
                                       (:acsc, :asin, :acsch, :asinh, "cosecant"),
                                       (:acot, :atan, :acoth, :atanh, "cotangent"))
    tname = string(tfa)
    hname = string(hfa)
    @eval begin
        @doc """
            $($tname)(x)
        Compute the inverse $($fn) of `x`, where the output is in radians. """ ($tfa)(y::T) where {T<:Number} = ($tfainv)(one(T) / y)
        @doc """
            $($hname)(x)
        Compute the inverse hyperbolic $($fn) of `x`. """ ($hfa)(y::T) where {T<:Number} = ($hfainv)(one(T) / y)
    end
end

# multiply in extended precision
function deg2rad_ext(x::Float64)
    m = 0.017453292519943295
    m_hi = 0.01745329238474369
    m_lo = 1.3519960527851425e-10

    u = 134217729.0*x # 0x1p27 + 1
    x_hi = u-(u-x)
    x_lo = x-x_hi

    y_hi = m*x
    y_lo = x_hi * m_lo + (x_lo* m_hi + ((x_hi*m_hi-y_hi) + x_lo*m_lo))

    DoubleFloat64(y_hi,y_lo)
end
deg2rad_ext(x::Float32) = DoubleFloat32(deg2rad(Float64(x)))
deg2rad_ext(x::Real) = deg2rad(x) # Fallback

function sind(x::Real)
    if isinf(x)
        return throw(DomainError(x, "`x` cannot be infinite."))
    elseif isnan(x)
        return oftype(x,NaN)
    end

    rx = copysign(float(rem(x,360)),x)
    arx = abs(rx)

    if rx == zero(rx)
        return rx
    elseif arx < oftype(rx,45)
        return sin_kernel(deg2rad_ext(rx))
    elseif arx <= oftype(rx,135)
        y = deg2rad_ext(oftype(rx,90) - arx)
        return copysign(cos_kernel(y),rx)
    elseif arx == oftype(rx,180)
        return copysign(zero(rx),rx)
    elseif arx < oftype(rx,225)
        y = deg2rad_ext((oftype(rx,180) - arx)*sign(rx))
        return sin_kernel(y)
    elseif arx <= oftype(rx,315)
        y = deg2rad_ext(oftype(rx,270) - arx)
        return -copysign(cos_kernel(y),rx)
    else
        y = deg2rad_ext(rx - copysign(oftype(rx,360),rx))
        return sin_kernel(y)
    end
end

function cosd(x::Real)
    if isinf(x)
        return throw(DomainError(x, "`x` cannot be infinite."))
    elseif isnan(x)
        return oftype(x,NaN)
    end

    rx = abs(float(rem(x,360)))

    if rx <= oftype(rx,45)
        return cos_kernel(deg2rad_ext(rx))
    elseif rx < oftype(rx,135)
        y = deg2rad_ext(oftype(rx,90) - rx)
        return sin_kernel(y)
    elseif rx <= oftype(rx,225)
        y = deg2rad_ext(oftype(rx,180) - rx)
        return -cos_kernel(y)
    elseif rx < oftype(rx,315)
        y = deg2rad_ext(rx - oftype(rx,270))
        return sin_kernel(y)
    else
        y = deg2rad_ext(oftype(rx,360) - rx)
        return cos_kernel(y)
    end
end

tand(x::Real) = sind(x) / cosd(x)

for (fd, f, fn) in ((:sind, :sin, "sine"), (:cosd, :cos, "cosine"), (:tand, :tan, "tangent"))
    name = string(fd)
    @eval begin
        @doc """
            $($name)(x)
        Compute $($fn) of `x`, where `x` is in degrees. """ ($fd)(z) = ($f)(deg2rad(z))
    end
end

for (fd, f, fn) in ((:asind, :asin, "sine"), (:acosd, :acos, "cosine"), (:atand, :atan, "tangent"),
                    (:asecd, :asec, "secant"), (:acscd, :acsc, "cosecant"), (:acotd, :acot, "cotangent"))
    name = string(fd)
    @eval begin
        @doc """
            $($name)(x)

        Compute the inverse $($fn) of `x`, where the output is in degrees. """ ($fd)(y) = rad2deg(($f)(y))
    end
end

max_stars_repo_name: xorJane/julia

max_issues_repo_name: xorJane/julia

max_forks_repo_name: xorJane/julia


text:
module Kosmos

using Dates, Distributed

export Component,
    ComponentSet,
    system!,
    Entity!,
    Kosmo,
    Clock

include("component.jl")
include("clock.jl")
include("kosmo.jl")
include("entity.jl")
include("system.jl")
include("smartlock.jl")

end


text:
x, y = rand(1000), rand(1000)

import HypothesisTests: OneSampleTTest

@test joint_distance_distribution(x, y) isa Vector
@test joint_distance_distribution(OneSampleTTest, x, y) isa OneSampleTTest


text:
import Base: copy, copyto!, ==
import Random: rand, seed!
import RandomNumbers: AbstractRNG, gen_seed, seed_type, unsafe_copyto!, unsafe_compare

"""
```julia
AbstractXorshift1024 <: AbstractRNG{UInt64}
```

The base abstract type for `Xorshift1024`, `Xorshift1024Star` and `Xorshift1024Plus`.
"""
abstract type AbstractXorshift1024 <: AbstractRNG{UInt64} end

for (star, plus) in (
        (false, false),
        (false, true),
        (true, false),
    )
    rng_name = Symbol(string("Xorshift1024", star ? "Star" : plus ? "Plus" : ""))
    @eval begin
        mutable struct $rng_name <: AbstractXorshift1024
            s0::UInt64
            s1::UInt64
            s2::UInt64
            s3::UInt64
            s4::UInt64
            s5::UInt64
            s6::UInt64
            s7::UInt64
            s8::UInt64
            s9::UInt64
            s10::UInt64
            s11::UInt64
            s12::UInt64
            s13::UInt64
            s14::UInt64
            s15::UInt64
            p::Int
        end

        function $rng_name(seed::NTuple{16, UInt64}=gen_seed(UInt64, 16))
            o = zero(UInt64)
            r = $rng_name(
                o, o, o, o, o, o, o, o,
                o, o, o, o, o, o, o, o,
                zero(Int)
            )
            seed!(r, seed)
            r
        end

        function $rng_name(seed::Integer)
            $rng_name(init_seed(seed, UInt64, 16))
        end

        @inline function xorshift_next(r::$rng_name)
            ptr = Ptr{UInt64}(pointer_from_objref(r))
            p = r.p
            s0 = unsafe_load(ptr, p + 1)
            p = (p + 1) % 16
            s1 = unsafe_load(ptr, p + 1)
            s1 ⊻= s1 << 31
            s1 ⊻= s1 >> 11
            s1 ⊻= s0 >> 30
            s1 ⊻= s0
            unsafe_store!(ptr, s1, p + 1)
            r.p = p
            $(star ? :(s1 * 2685821657736338717) :
              plus ? :(s1 + s0) : :s1)
        end
    end
end

@inline seed_type(::Type{T}) where T <: AbstractXorshift1024 = NTuple{16, UInt64}

function copyto!(dest::T, src::T) where T <: AbstractXorshift1024
    unsafe_copyto!(dest, src, UInt64, 16)
    dest.p = src.p
    dest
end

copy(src::T) where T <: AbstractXorshift1024 = copyto!(T(), src)

==(r1::T, r2::T) where T <: AbstractXorshift1024 = unsafe_compare(r1, r2, UInt64, 16) && r1.p == r2.p

function seed!(r::AbstractXorshift1024, seed::NTuple{16, UInt64})
    all(==(0), seed) && error("0 cannot be the seed")
    ptr = Ptr{UInt64}(pointer_from_objref(r))
    @inbounds for i in 1:16
        unsafe_store!(ptr, seed[i], i)
    end
    r.p = 0
    for i in 1:16
        xorshift_next(r)
    end
    r
end

function seed!(r::AbstractXorshift1024, seed::Integer...)
    l = length(seed)
    @assert 0 ≤ l ≤ 16
    if l == 0
        return seed!(r, gen_seed(UInt64, 16))
    end
    seed = [x % UInt64 for x in seed]
    while length(seed) < 16
        push!(seed, splitmix64(seed[end]))
    end
    seed!(r, Tuple(seed))
end

@inline rand(r::AbstractXorshift1024, ::Type{UInt64}) = xorshift_next(r)


text:
using TuringCallbacks
using Test

@testset "TuringCallbacks.jl" begin
    # Write your tests here.
end


text:
function transform_density!(img::AbstractArray,edges::AbstractArray, cdf::AbstractArray, minval::Union{Real,AbstractGray}, maxval::Union{Real,AbstractGray})
    first_edge = first(edges)
    inv_step_size = 1/step(edges)
    scale = (maxval - minval) / (cdf[end] - first(cdf))
    function transform(val)
        if val >= edges[end]
            newval = cdf[end]
        elseif val < first_edge
            newval = first(cdf)
        else
            index = floor(Int, (val-first_edge)*inv_step_size) + 1
            newval = cdf[index]
        end
        # Scale the new intensity value to so that it lies in the range [minval, maxval].
        newval = minval + (newval - first(cdf)) * scale
    end
    if eltype(img) <: Integer
        map!(val->ceil(transform(val)), img, img)
    else
        map!(transform, img, img)
    end
end

function construct_pdfs(img::AbstractArray, targetimg::AbstractArray, edges::AbstractRange)
    _, histogram = build_histogram(img, edges)
    _, target_histogram = build_histogram(targetimg, edges)
    return edges, histogram / sum(histogram), target_histogram / sum(target_histogram)
end

function construct_pdfs(img::AbstractArray, targetimg::AbstractArray, nbins::Integer = 256)
    if eltype(img) <: AbstractGray
        imin, imax = 0, 1
    else
        imin, imax = min(minfinite(img), minfinite(targetimg)), max(maxfinite(img), maxfinite(targetimg))
    end
    edges, histogram = build_histogram(img, nbins, minval = imin, maxval = imax)
    _, target_histogram = build_histogram(targetimg, edges)
    return edges, histogram / sum(histogram), target_histogram / sum(target_histogram)
end

function lookup_icdf(cdf::AbstractArray, targetcdf::AbstractArray)
    lookup_table = zeros(Int, length(cdf))
    i = 1
    for j = 1:length(cdf)
        p = cdf[j]
        while i < length(targetcdf) && targetcdf[i+1] <= p
            i += 1
        end
        lookup_table[j] = i
    end
    lookup_table
end

function cdf2pdf!(pdf::AbstractArray, cdf::AbstractArray)
    pdf[1] = cdf[1]
    for i = length(cdf)-1:-1:2
        pdf[i] = cdf[i] - cdf[i-1]
    end
end


text:
"""
Train a mixture of probabilistic circuits from data, starting with random example weights.
"""
function train_mixture( pcs::Vector{<:ProbΔ},
                        train_x::XBatches{Bool},
                        pseudocount, num_iters;
                        structure_learner=nothing, learnstruct_step = num_iters + 1, # structure learning
                        logger=nothing, logging_step = 1 # logging or saving results
                      )::AbstractFlatMixture

    # create mixture model with uniform component weights
    mixture_flow = init_mixture_with_flows(FlatMixture(pcs), train_x)

    # reset aggregate statistics
    reset_mixture_aggregate_flows(mixture_flow)

    # do a quick maximization step
    for batch in train_x
        example_weights = random_example_weights(num_examples(batch), num_components(mixture_flow))
        aggregate_flows(mixture_flow, batch, example_weights)
    end
    estimate_parameters(mixture_flow, component_weights(mixture_flow); pseudocount=pseudocount)

    train_mixture(mixture_flow, train_x, pseudocount, num_iters;
                structure_learner=structure_learner, learnstruct_step=learnstruct_step,
                logger=logger, logging_step=logging_step)
end

"""
Train a mixture model from data.
Learning is initialized from the parameters stored in the given mixture.
When a `structure_learner` is given, it will be called between EM steps to update circuit structures.
"""
function train_mixture( mixture::AbstractFlatMixture, # we start from component weights that are already given
                        train_x::XBatches{Bool},
                        pseudocount, num_iters;
                        structure_learner=nothing, learnstruct_step = num_iters + 1, # structure learning
                        logger=nothing, logging_step = 1 # logging or saving results
                      )::AbstractFlatMixture

    @assert feature_type(train_x) == Bool "Can only learn probabilistic circuits on Bool data"

    # initialize data structures
    mixture_flow = init_mixture_with_flows(mixture, train_x)

    if issomething(logger)
        logger(mixture_flow)
    end

    for i in 1:num_iters

        # reset aggregate statistics
        total_component_probability = ones(Float64, num_components(mixture_flow)) .* pseudocount ./ num_components(mixture_flow)
        reset_mixture_aggregate_flows(mixture_flow)

        # are we doing structure learning at the end of this iteration?
        is_learnstruct_iter = issomething(structure_learner) && i % learnstruct_step == 0

        all_example_weights = Vector{Matrix{Float64}}()

        # Expectation step (update example weights given mixture parameters)
        # + collecting aggregate statistics for subsequent maximization step
        for batch in train_x
            log_p_of_x_and_c = log_likelihood_per_instance_component(mixture_flow, batch)
            example_weights = component_weights_per_example(log_p_of_x_and_c)

            # copy the flows already computed by `log_likelihood_per_instance_component` into the underlying aggregate flow circuit
            # this way the maximization step can use them to estimate new parameters
            aggregate_flows_cached(mixture_flow, batch, example_weights)

            # store the aggregated component probabilities such that the maximization step can re-estimate the component weights
            total_component_probability .+= dropdims(sum(example_weights, dims=1), dims=1)

            # cache the example weights for the structure learner at the end of this EM iteration
            is_learnstruct_iter && push!(all_example_weights, example_weights)
        end

        # Maximization step (update mixture parameters given example weights (as stored in aggregate circuits))
        estimate_parameters(mixture_flow, total_component_probability; pseudocount=pseudocount)

        # Structural EM step
        if is_learnstruct_iter
            new_mixture_flow = structure_learner(mixture_flow, train_x, all_example_weights)
            # mixture = replace_prob_circuits(mixture, new_pcs)
            # re-initialize data structures
            mixture_flow = init_mixture_with_flows(new_mixture_flow, train_x)
        end

        if i % logging_step == 0 && issomething(logger)
            logger(mixture_flow)
        end
    end

    return mixture_flow
end

"Ensure we have a FlatMixtureWithFlow where the flow circuits have aggregate flow circuits as origin"
function init_mixture_with_flows(mixture::FlatMixtureWithFlow, ::XBatches{Bool})::FlatMixtureWithFlow
    if ! all(fc -> grand_origin(fc) isa AggregateFlowΔ, mixture.flowcircuits)
        init_mixture_with_flows(origin(mixture))
    else
        mixture
    end
end
function init_mixture_with_flows(mixture::FlatMixture, train_x::XBatches{Bool})::FlatMixtureWithFlow
    aggr_circuits = [AggregateFlowΔ(pc, Float64) for pc in components(mixture)]
    flow_circuits = [FlowΔ(afc, max_batch_size(train_x), Bool, opts_accumulate_flows) for afc in aggr_circuits]
    FlatMixtureWithFlow(mixture, flow_circuits)
end

function reset_mixture_aggregate_flows(mixture_flow::FlatMixtureWithFlow)
    for fc in mixture_flow.flowcircuits
        reset_aggregate_flows(grand_origin(fc))
    end
end

"Compute the component weights for each example from likelihoods"
function component_weights_per_example(log_p_of_x_and_c)
    log_p_of_x = logsumexp(log_p_of_x_and_c, 2) # marginalize out components
    log_p_of_given_x_query_c = mapslices(col -> col .- log_p_of_x, log_p_of_x_and_c, dims=[1])
    p_of_given_x_query_c = exp.(log_p_of_given_x_query_c) # no more risk of underflow, so go to linear space
    @assert sum(p_of_given_x_query_c) ≈ size(log_p_of_x_and_c, 1) # each row has proability 1
    p_of_given_x_query_c
end

"Compute and aggregate flows for mixture components"
function aggregate_flows(mixture_flow, batch, example_weights)
    for i in 1:num_components(mixture_flow)
        fc =  mixture_flow.flowcircuits[i]
        wbatch = weighted_batch_for_component(batch, example_weights,i)
        accumulate_aggr_flows_batch(fc, wbatch)
    end
end

"Aggregate already-computed flows for mixture components"
function aggregate_flows_cached(mixture_flow, batch, example_weights)
    for i in 1:num_components(mixture_flow)
        fc =  mixture_flow.flowcircuits[i]
        wbatch = weighted_batch_for_component(batch, example_weights,i)
        accumulate_aggr_flows_cached(fc, wbatch)
    end
end

function estimate_parameters(mixture_flow, total_component_probability; pseudocount)
    component_weights(mixture_flow) .= total_component_probability ./ sum(total_component_probability)
    for fc in mixture_flow.flowcircuits
        estimate_parameters_cached(grand_origin(fc); pseudocount=pseudocount)
    end
end

"Get a new weighted batch for this component"
weighted_batch_for_component(batch::PlainXData, example_weights, component_i)::WXData =
    WXData(batch, example_weights[:,component_i])

"Create random example weights that sum to one overall components"
function random_example_weights(num_examples::Int, num_components::Int)::Matrix{Float64}
    w = rand(Float64, num_examples, num_components)
    w ./ sum(w;dims=2)
end


text:
#=
    Figure S4
=#
using Optim

xlim = [0.0,0.5]
xgrd = range(xlim...,length=200)

μ = 0.25
σ = [0.01,0.05,0.1]

ω = [-1.5,-0.5,0.0001,0.5,1.5]

σᵢ = σ[1]

plts = [plot() for i = 1:length(σ)]

for (i,σᵢ) in enumerate(σ)

    pltsᵢ = [plot() for i = 1:length(ω)]

    for (j,ωⱼ) in enumerate(ω)

        λ = GammaAlt(μ,σᵢ,ωⱼ)

        # Find mode
        m = optimize(x -> -pdf(λ,x),xlim...).minimizer

        if ωⱼ > 0.0
            xmap = x -> x .> m ? 2m - x : x
        elseif ωⱼ < 0.0
            xmap = x -> x .< m ? 2m - x : x
        else
            xmap = x -> x
        end

        # Plot symmetric distribution
        plot!(pltsᵢ[j],xgrd,pdf(λ,xmap.(xgrd)),lw=2.0,c=:red,α=0.5,ls=:dash)
        plot!(pltsᵢ[j],xgrd,pdf(λ,xgrd),lw=2.0,c=:black)

    end

    plts[i] = plot(pltsᵢ...,layout=grid(1,length(ω)))
    if i < length(σ)
        plot!(plts[i],xticks=(0.0:0.1:0.5,[]))
    else
        plot!(plts[i],xticks=0.0:0.1:0.5)
    end
    plot!(plts[i],subplot=1,ylabel="σ = $σᵢ")

end

figS4 = plot(plts...,layout=grid(length(σ),1),yticks=[],axis=:x,legend=:none,size=(700,400))
[plot!(figS4,subplot=i,title="ω = $ωᵢ") for (i,ωᵢ) in enumerate(ω)]
[plot!(figS4,subplot=i,xlabel="λ") for i = 11:15]
figS4

##############################################################
## Save...
##############################################################

savefig("$(@__DIR__)/FigS4.svg")

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
module Stan

using StanSample, StanOptimize, StanDiagnose, StanVariational

end # module


text:
using ReactionNetworkImporters, Test

@time begin
    @time @testset "BNG Birth-Death Test" begin include("test_nullrxs_odes.jl") end
    @time @testset "BNG Repressilator Test" begin include("test_repressilator_odes.jl") end
    @time @testset "BNG Higher Order Test" begin include("test_higherorder_odes.jl") end
    @time @testset "Matrix Input Test" begin include("test_mats.jl") end
end


text:
"""
    CrystalNets.jl

Module for automatic recognition of crystal net topologies.
To use as an executable, run the source file in a shell:
```bash
julia --project=/path/to/CrystalNets/ /path/to/CrystalNets/src/CrystalNets.jl
```
Otherwise, as a module, to try to recognize the net underlying a crystal given in a
chemical file format called FILE, the entry point is the following execution:
```julia
julia> using CrystalNets

julia> determine_topology(FILE)
```

See the documentation at http://molsim.info/CrystalNets.jl/ for further instructions on
on the use of CrystalNets.jl, and how to do a full installation to reduce latency at
startup.
"""
module CrystalNets

export CrystalNet,
       CrystalNet1D, CrystalNet2D, CrystalNet3D,
       UnderlyingNets,
       TopologicalGenome,
       TopologyResult,
       determine_topology,
       determine_topology_dataset,
       guess_topology,
       guess_topology_dataset,
       parse_chemfile,
       topological_genome,
       recognize_topology,
       StructureType,
       Bonding,
       Clustering,
       ClusterKinds

import LinearAlgebra: det, norm, rank, cross
using Base.Threads
import Serialization

using PeriodicGraphEmbeddings
using PeriodicGraphs
using StaticArrays
using Graphs

import Logging

const DOWARN = Base.RefValue{Bool}(false)
const DOEXPORT = Base.RefValue{Bool}(false)

"""
    toggle_warning(to=nothing)

Toggle warnings on (if `to == true`) or off (if `to == false`).
Without an argument, toggle on and off repeatedly at each call.
"""
function toggle_warning(to=nothing)
    global DOWARN
    DOWARN[] = to isa Nothing ? !DOWARN[] : to
end

"""
    toggle_export(to=nothing)

Toggle exports on (if `to == true`) or off (if `to == false`).
Without an argument, toggle on and off repeatedly at each call.
"""
function toggle_export(to=nothing)
    global DOEXPORT
    DOEXPORT[] = to isa Nothing ? !DOEXPORT[] : to
end

function __init__()
    toggle_warning("--no-warn" ∉ ARGS)
    toggle_export("--no-export" ∉ ARGS)
    @static if abspath(PROGRAM_FILE) == @__FILE__
        global DOWARN
        if DOWARN[]
            @info "Proceed to a full installation (see README.md) for better performance."
        end
        @show "EXECUTING JULIA MAIN"
        ret_code = CrystalNets.julia_main()
        exit(ret_code)
    end
    nothing
end

include("utils.jl")
__precompile__(true)
include("specialsolver.jl") # Exact sparse symmetric integer matrix factorization
include("options.jl") # Computation options
include("types.jl") # Main internal type definitions used to represent topologies
include("input.jl") # Crystal file parsing and conversion to an internal type
include("archive.jl") # Manipulation of the topological archive
include("output.jl") # Crystal file exports
include("arithmetics.jl") # Handling of (possibly sparse) integer/rational matrices
include("stability.jl") # Functions related to unstable nets
include("topology.jl") # Main functions of the algorithm
include("query.jl") # Entry point for the user-facing functions
include("executable.jl") # Entry point for the argument parsing of the executable
include("precompile.jl")

end # module CrystalNets

nothing


text:
mutable struct TestPlan{T,N} <: Plan{T}
    region
    sz::NTuple{N,Int}
    pinv::Plan{T}
    function TestPlan{T}(region, sz::NTuple{N,Int}) where {T,N}
        return new{T,N}(region, sz)
    end
end

mutable struct InverseTestPlan{T,N} <: Plan{T}
    region
    sz::NTuple{N,Int}
    pinv::Plan{T}
    function InverseTestPlan{T}(region, sz::NTuple{N,Int}) where {T,N}
        return new{T,N}(region, sz)
    end
end

Base.size(p::TestPlan) = p.sz
Base.ndims(::TestPlan{T,N}) where {T,N} = N
Base.size(p::InverseTestPlan) = p.sz
Base.ndims(::InverseTestPlan{T,N}) where {T,N} = N

function AbstractFFTs.plan_fft(x::AbstractArray{T}, region; kwargs...) where {T}
    return TestPlan{T}(region, size(x))
end
function AbstractFFTs.plan_bfft(x::AbstractArray{T}, region; kwargs...) where {T}
    return InverseTestPlan{T}(region, size(x))
end
function AbstractFFTs.plan_inv(p::TestPlan{T}) where {T}
    unscaled_pinv = InverseTestPlan{T}(p.region, p.sz)
    unscaled_pinv.pinv = p
    pinv = AbstractFFTs.ScaledPlan(
        unscaled_pinv, AbstractFFTs.normalization(T, p.sz, p.region),
    )
    return pinv
end
function AbstractFFTs.plan_inv(p::InverseTestPlan{T}) where {T}
    unscaled_pinv = TestPlan{T}(p.region, p.sz)
    unscaled_pinv.pinv = p
    pinv = AbstractFFTs.ScaledPlan(
        unscaled_pinv, AbstractFFTs.normalization(T, p.sz, p.region),
    )
    return pinv
end

# Just a helper function since forward and backward are nearly identical
# The function does not check if the size of `y` and `x` are compatible, this
# is done in the function where `dft!` is called since the check differs for FFTs
# with complex and real-valued signals
function dft!(
    y::AbstractArray{<:Complex,N},
    x::AbstractArray{<:Union{Complex,Real},N},
    dims,
    sign::Int
) where {N}
    # check that dimensions that are transformed are unique
    allunique(dims) || error("dimensions have to be unique")

    T = eltype(y)
    # we use `size(x, d)` since for real-valued signals
    # `size(y, first(dims)) = size(x, first(dims)) ÷ 2 + 1`
    cs = map(d -> T(sign * 2π / size(x, d)), dims)
    fill!(y, zero(T))
    for yidx in CartesianIndices(y)
        # set of indices of `x` on which `y[yidx]` depends
        xindices = CartesianIndices(
            ntuple(i -> i in dims ? axes(x, i) : yidx[i]:yidx[i], Val(N))
        )
        for xidx in xindices
            y[yidx] += x[xidx] * cis(sum(c * (yidx[d] - 1) * (xidx[d] - 1) for (c, d) in zip(cs, dims)))
        end
    end
    return y
end

function mul!(
    y::AbstractArray{<:Complex,N}, p::TestPlan, x::AbstractArray{<:Union{Complex,Real},N}
) where {N}
    size(y) == size(p) == size(x) || throw(DimensionMismatch())
    dft!(y, x, p.region, -1)
end
function mul!(
    y::AbstractArray{<:Complex,N}, p::InverseTestPlan, x::AbstractArray{<:Union{Complex,Real},N}
) where {N}
    size(y) == size(p) == size(x) || throw(DimensionMismatch())
    dft!(y, x, p.region, 1)
end

Base.:*(p::TestPlan, x::AbstractArray) = mul!(similar(x, complex(float(eltype(x)))), p, x)
Base.:*(p::InverseTestPlan, x::AbstractArray) = mul!(similar(x, complex(float(eltype(x)))), p, x)

mutable struct TestRPlan{T,N} <: Plan{T}
    region
    sz::NTuple{N,Int}
    pinv::Plan{T}
    TestRPlan{T}(region, sz::NTuple{N,Int}) where {T,N} = new{T,N}(region, sz)
end

mutable struct InverseTestRPlan{T,N} <: Plan{T}
    d::Int
    region
    sz::NTuple{N,Int}
    pinv::Plan{T}
    function InverseTestRPlan{T}(d::Int, region, sz::NTuple{N,Int}) where {T,N}
        sz[first(region)::Int] == d ÷ 2 + 1 || error("incompatible dimensions")
        return new{T,N}(d, region, sz)
    end
end

function AbstractFFTs.plan_rfft(x::AbstractArray{T}, region; kwargs...) where {T}
    return TestRPlan{T}(region, size(x))
end
function AbstractFFTs.plan_brfft(x::AbstractArray{T}, d, region; kwargs...) where {T}
    return InverseTestRPlan{T}(d, region, size(x))
end
function AbstractFFTs.plan_inv(p::TestRPlan{T,N}) where {T,N}
    firstdim = first(p.region)::Int
    d = p.sz[firstdim]
    sz = ntuple(i -> i == firstdim ? d ÷ 2 + 1 : p.sz[i], Val(N))
    unscaled_pinv = InverseTestRPlan{T}(d, p.region, sz)
    unscaled_pinv.pinv = p
    pinv = AbstractFFTs.ScaledPlan(
        unscaled_pinv, AbstractFFTs.normalization(T, p.sz, p.region),
    )
    return pinv
end
function AbstractFFTs.plan_inv(p::InverseTestRPlan{T,N}) where {T,N}
    firstdim = first(p.region)::Int
    sz = ntuple(i -> i == firstdim ? p.d : p.sz[i], Val(N))
    unscaled_pinv = TestRPlan{T}(p.region, sz)
    unscaled_pinv.pinv = p
    pinv = AbstractFFTs.ScaledPlan(
        unscaled_pinv, AbstractFFTs.normalization(T, sz, p.region),
    )
    return pinv
end

Base.size(p::TestRPlan) = p.sz
Base.ndims(::TestRPlan{T,N}) where {T,N} = N
Base.size(p::InverseTestRPlan) = p.sz
Base.ndims(::InverseTestRPlan{T,N}) where {T,N} = N

function real_invdft!(
    y::AbstractArray{<:Real,N},
    x::AbstractArray{<:Union{Complex,Real},N},
    dims,
) where {N}
    # check that dimensions that are transformed are unique
    allunique(dims) || error("dimensions have to be unique")

    firstdim = first(dims)
    size_x_firstdim = size(x, firstdim)
    iseven_firstdim = iseven(size(y, firstdim))
    # we do not check that the input corresponds to a real-valued signal
    # (i.e., that the first and, if `iseven_firstdim`, the last value in dimension
    # `haldim` of `x` are real values) due to numerical inaccuracies
    # instead we just use the real part of these entries

    T = eltype(y)
    # we use `size(y, d)` since `size(x, first(dims)) = size(y, first(dims)) ÷ 2 + 1`
    cs = map(d -> T(2π / size(y, d)), dims)
    fill!(y, zero(T))
    for yidx in CartesianIndices(y)
        # set of indices of `x` on which `y[yidx]` depends
        xindices = CartesianIndices(
            ntuple(i -> i in dims ? axes(x, i) : yidx[i]:yidx[i], Val(N))
        )
        for xidx in xindices
            coeffimag, coeffreal = sincos(
                sum(c * (yidx[d] - 1) * (xidx[d] - 1) for (c, d) in zip(cs, dims))
            )

            # the first and, if `iseven_firstdim`, the last term of the DFT are scaled
            # with 1 instead of 2 and only the real part is used (see note above)
            xidx_firstdim = xidx[firstdim]
            if xidx_firstdim == 1 || (iseven_firstdim && xidx_firstdim == size_x_firstdim)
                y[yidx] += coeffreal * real(x[xidx])
            else
                xreal, ximag = reim(x[xidx])
                y[yidx] += 2 * (coeffreal * xreal - coeffimag * ximag)
            end
        end
    end

    return y
end

to_real!(x::AbstractArray) = map!(real, x, x)

function Base.:*(p::TestRPlan, x::AbstractArray)
    size(p) == size(x) || error("array and plan are not consistent")

    # create output array
    firstdim = first(p.region)::Int
    d = size(x, firstdim)
    firstdim_size = d ÷ 2 + 1
    T = complex(float(eltype(x)))
    sz = ntuple(i -> i == firstdim ? firstdim_size : size(x, i), Val(ndims(x)))
    y = similar(x, T, sz)

    # compute DFT
    dft!(y, x, p.region, -1)

    # we clean the output a bit to make sure that we return real values
    # whenever the output is mathematically guaranteed to be a real number
    to_real!(selectdim(y, firstdim, 1))
    if iseven(d)
        to_real!(selectdim(y, firstdim, firstdim_size))
    end

    return y
end

function Base.:*(p::InverseTestRPlan, x::AbstractArray)
    size(p) == size(x) || error("array and plan are not consistent")

    # create output array
    firstdim = first(p.region)::Int
    d = p.d
    sz = ntuple(i -> i == firstdim ? d : size(x, i), Val(ndims(x)))
    y = similar(x, real(float(eltype(x))), sz)

    # compute DFT
    real_invdft!(y, x, p.region)

    return y
end


text:
module TimeMachine

include("Dtw.jl")

end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =================================================================================================
"""
    tess_text(
        inst::TessInst
    )::Union{String, Nothing}

Extract the text from the image. If there is an error `nothing` will be returned.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | inst |         | The instance to grab the text from.

__Details:__

This method will call `tess_recognize()` if it has not been called yet for the image.

__Example:__

```jldoctest
using Tesseract

download_languages()

instance = TessInst()
pix = sample_pix()

tess_image(instance, pix)
tess_resolution(instance, 72)

text = tess_text(instance)

for line in split(text, '\\n'; keepempty = false)[1:5]
    println(line)
end

# output

No one would have believed in the last years of the
the nineteenth century that this world was being watched
watched keenly and closely by intelligences greater than
than man’s and yet as mortal as his own; that as men busied
busied themselves about their various concerns they were
```

See also: [`tess_hocr`](@ref), [`tess_alto`](@ref), [`tess_tsv`](@ref),
          [`tess_parsed_tsv`](@ref), [`tess_confidences`](@ref)
"""
function tess_text(
            inst::TessInst
        )::Union{String, Nothing}
    local retval = nothing

    if is_valid(inst) == true
        local text = ccall(
            (:TessBaseAPIGetUTF8Text, TESSERACT),
            Cstring,
            (Ptr{Cvoid},),
            inst
        )

        if text != C_NULL
            retval = unsafe_string(text)
            delete_text(text)
        end
    else
        @error "Instance has been freed."
    end

    return retval
end

# =================================================================================================
"""
    tess_hocr(
        inst::TessInst,
        page::Integer = Int32(1)
    )::Union{String, Nothing}

Extract the text in hOCR format from the image.  Returns `nothing` if there is an error.

__Arguments:__

| T | Name | Default    | Description
|:--| :--- | :--------- | :----------
| R | inst |            | The instance to grab the text from.
| O | page | `Int32(1)` | The page to extract the hOCR text for.

__Details:__

This method will call `tess_recognize()` if it has not been called yet for the image.  The current
hOCR spec can be accessed at http://kba.cloud/hocr-spec/1.2/.

__Example:__

```jldoctest
using Tesseract

download_languages()

instance = TessInst()
pix = sample_pix()

tess_image(instance, pix)
tess_resolution(instance, 72)

hocr = tess_hocr(instance)

for line in split(hocr, '\\n'; keepempty = false)[1:5]
    println(strip(line))
end

# output

<div class='ocr_page' id='page_1' title='image ""; bbox 0 0 500 600; ppageno 0'>
<div class='ocr_carea' id='block_1_1' title="bbox 10 9 489 523">
<p class='ocr_par' id='par_1_1' lang='eng' title="bbox 11 9 417 23">
<span class='ocr_line' id='line_1_1' title="bbox 11 9 417 23; baseline 0 -3; x_size 22.717392; x_descenders 5.5; x_ascenders 5.7391305">
<span class='ocrx_word' id='word_1_1' title='bbox 11 9 25 20; x_wconf 95'>No</span>
```

See also: [`tess_text`](@ref), [`tess_alto`](@ref), [`tess_tsv`](@ref),
          [`tess_parsed_tsv`](@ref), [`tess_confidences`](@ref)
"""
function tess_hocr(
            inst::TessInst,
            page::Integer = Int32(1)
        )::Union{String, Nothing}
    local retval = nothing

    if is_valid(inst) == true
        local text = ccall(
            (:TessBaseAPIGetHOCRText, TESSERACT),
            Cstring,
            (Ptr{Cvoid},Cint),
            inst,
            page-1

        )
        if text != C_NULL
            retval = unsafe_string(text)
            delete_text(text)
        end
    else
        @error "Instance has been freed."
    end

    return retval
end

# =================================================================================================
"""
    tess_alto(
        inst::TessInst,
        page::Integer = Int32(1)
    )::Union{String, Nothing}

Extract the text in ALTO format from the image.  Returns `nothing` if there is an error.

__Arguments:__

| T | Name | Default    | Description
|:--| :--- | :--------- | :----------
| R | inst |            | The instance to grab the text from.
| O | page | `Int32(1)` | The page to extract the ALTO text for.

__Details:__

This method will call `tess_recognize()` if it has not been called yet for the image.  The current
ALTO spec can be accessed at https://github.com/altoxml/documentation/wiki/Versions.

__Example:__

```jldoctest
using Tesseract

download_languages()

instance = TessInst()
pix = sample_pix()

tess_image(instance, pix)
tess_resolution(instance, 72)

alto = tess_alto(instance)

for line in split(alto, '\\n'; keepempty = false)[1:5]
    println(strip(line))
end

# output

<Page WIDTH="500" HEIGHT="600" PHYSICAL_IMG_NR="0" ID="page_0">
<PrintSpace HPOS="0" VPOS="0" WIDTH="500" HEIGHT="600">
<ComposedBlock ID="cblock_0" HPOS="10" VPOS="9" WIDTH="479" HEIGHT="514">
<TextBlock ID="block_0" HPOS="11" VPOS="9" WIDTH="406" HEIGHT="14">
<TextLine ID="line_0" HPOS="11" VPOS="9" WIDTH="406" HEIGHT="14">
```

See also: [`tess_text`](@ref), [`tess_hocr`](@ref), [`tess_tsv`](@ref), [`tess_parsed_tsv`](@ref),
[`tess_confidences`](@ref)
"""
function tess_alto(
            inst::TessInst,
            page::Integer = Int32(1)
        )::Union{String, Nothing}
    local retval = nothing

    if is_valid(inst) == true
        local text = ccall(
            (:TessBaseAPIGetAltoText, TESSERACT),
            Cstring,
            (Ptr{Cvoid},Cint),
            inst,
            page - 1
        )
        if text != C_NULL
            retval = unsafe_string(text)
            delete_text(text)
        end
    else
        @error "Instance has been freed."
    end

    return retval
end

# =================================================================================================
"""
    tess_tsv(
        inst::TessInst,
        page::Integer = Int32(1)
    )::Union{String, Nothing}

Retrieve the TSV results from an recognition as a string with tabbed separated values.

__Arguments:__

| T | Name | Default    | Description
|:--| :--- | :--------- | :----------
| R | inst |            | The Tesseract instance to call.
| O | page | `Int32(1)` | The page to get the data for.

__Details:__

This method will call `tess_recognize()` if it has not been called yet for the image.

If you want to read the data take a look at `tess_parsed_tsv()` which parses it into a form that
is easier to use.  If you just want to display the data to the user or write it to a file then
this method will probably be fine.

This call basically returns the results of the recognition process.  Each line in the string
identifies an "object" found by Tesseract.  The 12 values per line are as follows:

  * level - Identifies what the line describes.
  * page - This is the page number passed into the `tess_tsv()` method.
  * block - Identifies the block on the page.
  * paragraph - The paragraph number in the block.
  * line - The line in the paragraph.
  * word - The word in the line.
  * left - Left edge of the item in pixels.
  * top - Top edge of the item in pixels.
  * width - Width of the item in pixels.
  * height - Height of the item in pixels.
  * conf - How confident the OCR engine is of the word (`0` - `100`). `-1` if level is not 5.
  * text - The word that was decoded from the image.

Level identifies what information the line is providing:

* `1` - Page information, added at the start of the page.
* `2` - Block information, added at the start of a block.
* `3` - Paragraph information, added at the start of a paragraph.
* `4` - Line information, added at the start of a line.
* `5` - Word information, identifies a word that was read from the page.

The left, top, width, and height values define a box in pixels that encompases the item.  So if the
level is 1, the box describes the whole image.  If the level is `1`, then the box encloses the block
that was extracted, and so on down to the word that was extracted.

__Example:__

```jldoctest; filter = r"(\\s+)"
using Tesseract

download_languages()

instance = TessInst()
pix = sample_pix()

tess_image(instance, pix)
tess_resolution(instance, 72)

tsv = tess_tsv(instance)

for line in split(tsv, '\\n'; keepempty = false)[2:6]
    println(strip(line))
end

# output

2 1   1   0   0   0   10  9   479 514 -1
3 1   1   1   0   0   11  9   406 14  -1
4 1   1   1   1   0   11  9   406 14  -1
5 1   1   1   1   1   11  9   14  11  95  No
5 1   1   1   1   2   35  12  22  8   95  one
```

See also: [`tess_tsv`](@ref), [`tess_hocr`](@ref), [`tess_alto`](@ref), [`tess_parsed_tsv`](@ref),
[`tess_confidences`](@ref)
"""
function tess_tsv(
            inst::TessInst,
            page::Integer = Int32(1)
        )::Union{String, Nothing}
    local retval = nothing

    if is_valid(inst) == true
        local text = ccall(
            (:TessBaseAPIGetTsvText, TESSERACT),
            Cstring,
            (Ptr{Cvoid},Cint),
            inst,
            page - 1
        )
        if text != C_NULL
            retval = unsafe_string(text)
            delete_text(text)
        end
    else
        @error "Instance has been freed."
    end

    return retval
end

# =================================================================================================
"""
    tess_text_box(
        inst::TessInst,
        page::Integer = Int32(1)
    )::Union{String, Nothing}

Retrieve the boxes for the identified characters on the page.  If there is an error 'nothing' is
returned.

__Arguments:__

| T | Name | Default    | Description
|:--| :--- | :--------- | :----------
| R | inst |            | The Tesseract instance to call.
| O | page | `Int32(1)` | The page to get the data for.

__Details:__

This method will call `tess_recognize()` if it has not been called yet for the image.  The results
are primarly useful for training purposes.

Each line in the result is a character identified with 6 values:

  * The character that was recognized.
  * The left edge of the character measured in pixels from the left edge.
  * The bottom edge of the character measured in pixels from the bottom of the image.
  * The right edge of the character measured in pixels from the left edge.
  * The top edge of the character measured in pixels from the bottom of the image.
  * The page used in the recogniztion 0 based.

__Example:__

```jldoctest"
using Tesseract

download_languages()

instance = TessInst()
pix = sample_pix()

tess_image(instance, pix)
tess_resolution(instance, 72)

box = tess_text_box(instance)

for line in split(box, '\\n'; keepempty = false)[1:5]
  println(line)
end

# output

N 11 580 17 591 0
o 11 580 25 591 0
o 35 580 42 588 0
n 35 580 41 588 0
e 43 580 57 588 0
```

See also: [`tess_unlv`](@ref), [`tess_lstm_box`](@ref), [`tess_word_box`](@ref)
"""
function tess_text_box(
            inst::TessInst,
            page::Integer = Int32(1)
        )::Union{String, Nothing}
    local retval = nothing

    if is_valid(inst) == true
        local text = ccall(
            (:TessBaseAPIGetBoxText, TESSERACT),
            Cstring,
            (Ptr{Cvoid}, Cint),
            inst,
            page - 1
        )
        if text != C_NULL
            retval = unsafe_string(text)
            delete_text(text)
        end
    else
        @error "Instance has been freed."
    end

    return retval
end

# =================================================================================================
"""
    tess_word_box(
        inst::TessInst,
        page::Integer = Int32(1)
    )::Union{String, Nothing}

Create a UTF8 box file with WordStr strings,  If there is an error 'nothing' is returned.

__Arguments:__

| T | Name | Default    | Description
|:--| :--- | :--------- | :----------
| R | inst |            | The Tesseract instance to call.
| O | page | `Int32(1)` | The page to get the data for.

__Details:__

This method will call `tess_recognize()` if it has not been called yet for the image.  The results
are probably used for training.

__Example:__

```jldoctest; filter = r"(\\s+)"
using Tesseract

download_languages()

instance = TessInst()
pix = sample_pix()

tess_image(instance, pix)
tess_resolution(instance, 72)

box = tess_word_box(instance)

for line in split(box, '\\n'; keepempty = false)[1:5]
    println(line)
end

# output

WordStr 11 577 417 591 0 #No one would have believed in the last years of the
    418 577 422 591 0
WordStr 11 557 457 571 0 #the nineteenth century that this world was being watched
    458 557 462 571 0
WordStr 10 537 457 551 0 #watched keenly and closely by intelligences greater than
```

See also: [`tess_unlv`](@ref), [`tess_lstm_box`](@ref), [`tess_text_box`](@ref)
"""
function tess_word_box(
            inst::TessInst,
            page::Integer = Int32(1)
        )::Union{String, Nothing}
    local retval = nothing

    if is_valid(inst) == true
        local text = ccall(
            (:TessBaseAPIGetWordStrBoxText, TESSERACT),
            Cstring,
            (Ptr{Cvoid}, Cint),
            inst,
            page - 1
        )
        if text != C_NULL
            retval = unsafe_string(text)
            delete_text(text)
        end
    else
        @error "Instance has been freed."
    end

    return retval
end

# =================================================================================================
"""
    tess_lstm_box(
        inst::TessInst,
        page::Integer = Int32(1)
    )::Union{String, Nothing}

Return the UTF-8 box file for LSTM training.  If there is an error 'nothing' is returned.

__Arguments:__

| T | Name | Default    | Description
|:--| :--- | :--------- | :----------
| R | inst |            | The Tesseract instance to call.
| O | page | `Int32(1)` | The page to get the data for.

__Details:__

This method will call `tess_recognize()` if it has not been called yet for the image.

__Example:__

```jldoctest; filter = r"(\\s+)"
using Tesseract

download_languages()

instance = TessInst()
pix = sample_pix()

tess_image(instance, pix)
tess_resolution(instance, 72)

box = tess_lstm_box(instance)

for line in split(box, '\\n'; keepempty = false)[1:5]
    println(line)
end

# output

N 11 577 422 591 0
o 11 577 422 591 0
  11 577 422 591 0
o 11 577 422 591 0
n 11 577 422 591 0
```

See also: [`tess_unlv`](@ref), [`tess_word_box`](@ref), [`tess_text_box`](@ref)
"""
function tess_lstm_box(
            inst::TessInst,
            page::Integer = Int32(1)
        )::Union{String, Nothing}
    local retval = nothing

    if is_valid(inst) == true
        local text = ccall(
            (:TessBaseAPIGetLSTMBoxText, TESSERACT),
            Cstring,
            (Ptr{Cvoid}, Cint),
            inst,
            page - 1
        )
        if text != C_NULL
            retval = unsafe_string(text)
            delete_text(text)
        end
    else
        @error "Instance has been freed."
    end

    return retval
end

# =========================================================================================
"""
    tess_unlv(
        inst::TessInst
    )::Union{String, Nothing}

Extract the text in UNLV format UTF-8 with reject and suspect codes.  If there is an error
`nothing` is returned.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | inst |         | The Tesseract instance to call.

__Details:__

This method will call `tess_recognize()` if it has not been called yet for the image.

This method is more used to test the OCR results than anything else.  If you want the
original Latin1 encoding use [`tess_unlv_latin1`](@ref) method.

__Example:__

```jldoctest
using Tesseract

download_languages()

instance = TessInst()
pix = sample_pix()

tess_image(instance, pix)
tess_resolution(instance, 72)

unlv = tess_unlv(instance)

for line in split(unlv, '\\n'; keepempty = false)[1:5]
    println(line)
end

# output

No one would have believed in the last years of the
the nineteenth century that this world was being watched
watched keenly and closely by intelligences greater than
than man's and yet as mortal as his own; that as men busied
busied themselves about their various concerns they were
```

See also: [`tess_lstm_box`](@ref), [`tess_word_box`](@ref), [`tess_text_box`](@ref),
          [`tess_unlv_latin1`](@ref)
"""
function tess_unlv(
            inst::TessInst
        )::Union{String, Nothing}
    local retval = nothing

    if is_valid(inst) == true
        local bytes = ccall(
            (:TessBaseAPIGetUNLVText, TESSERACT),
            Ptr{UInt8},
            (Ptr{Cvoid},),
            inst
        )

        if bytes == C_NULL
            return nothing
        end

        local buffer = Vector{UInt8}()

        local offset = 1
        while unsafe_load(bytes, offset) != 0x00
            local by = unsafe_load(bytes, offset)
            if by < 128
                push!(buffer, by)
            else
                push!(buffer, 0xc0 | (by >> 6))
                push!(buffer, 0x80 | (by & 0x3f))
            end
            offset += 1
        end
        delete_text(bytes)

        retval = String(buffer)
    else
        @error "Instance has been freed."
    end

    return retval
end

# =========================================================================================
"""
    tess_unlv_latin1(
        inst::TessInst
    )::Union{String, Nothing}

Extract the text in UNLV format Latin-1 with reject and suspect codes.  If there is an error
`nothing` is returned.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | inst |         | The Tesseract instance to call.

__Details:__

This method will call `tess_recognize()` if it has not been called yet for the image.

This method is more used to test the OCR results than anything else.  This method returns
the OCR data in Latin1 encoding.  If you want to use it as a string in Julia use the
[`tess_unlv`](@ref) method which will convert it to UTF-8 for you.

__Example:__

```jldoctest
using Tesseract

download_languages()

instance = TessInst()
pix = sample_pix()

tess_image(instance, pix)
tess_resolution(instance, 72)

unlv = tess_unlv_latin1(instance)

# output

1470-element Array{UInt8,1}:
 0x4e
 0x6f
 0x20
 0x6f
 0x6e
 0x65
 0x20
 0x77
 0x6f
 0x75
    ⋮
 0x6f
 0x6e
 0x6d
 0x65
 0x6e
 0x74
 0x2e
 0x0a
 0x0a
```

See also: [`tess_lstm_box`](@ref), [`tess_word_box`](@ref), [`tess_text_box`](@ref),
          [`tess_unlv`](@ref)
"""
function tess_unlv_latin1(
            inst::TessInst
        )::Union{Vector{UInt8}, Nothing}
    local retval = nothing

    if is_valid(inst) == true
        local bytes = ccall(
            (:TessBaseAPIGetUNLVText, TESSERACT),
            Ptr{UInt8},
            (Ptr{Cvoid},),
            inst
        )

        if bytes == C_NULL
            return nothing
        end

        local size = 0
        while unsafe_load(bytes, size + 1) != 0x00
            size += 1
        end

        retval = Vector{UInt8}(undef, size)

        unsafe_copyto!(pointer(retval), bytes, size)

        delete_text(bytes)
    else
        @error "Instance has been freed."
    end

    return retval
end

# =================================================================================================
"""
    tess_confidences(
        inst::TessInst
    )::Union{Vector{Int}, Nothing}

Extract the confidences in the words that where extracted from the image.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | inst |         | The Tesseract instance to call.

__Details:__

This method will call `tess_recognize()` if it has not been called yet for the image.

__Example:__

```jldoctest
using Tesseract

download_languages()

instance = TessInst()
pix = sample_pix()

tess_image(instance, pix)
tess_resolution(instance, 72)

confidence = tess_confidences(instance)

# output

256-element Array{Int64,1}:
 95
 95
 93
 93
 96
 96
 96
 96
 96
 96
  ⋮
 96
 96
 96
 96
 96
 96
 96
 87
 87
```

See also: [`tess_tsv`](@ref), [`tess_parsed_tsv`](@ref)
"""
function tess_confidences(
            inst::TessInst
        )::Union{Vector{Int}, Nothing}

    # ---------------------------------------------------------------------------------------------
    # Make sure the object is still valid.
    if is_valid(inst) == false
        @error "Instance has been freed."
        return nothing
    end

    # ---------------------------------------------------------------------------------------------
    # Make the C call.
    local list = ccall(
        (:TessBaseAPIAllWordConfidences, TESSERACT),
        Ptr{Cint},
        (Ptr{Cvoid},),
        inst
    )

    if list == C_NULL
        return nothing
    end

    # ---------------------------------------------------------------------------------------------
    # Walk the pointer until we find a -1.
    local retval = Vector{Cint}()
    local offset = 1
    while unsafe_load(list, offset) != -1
        push!(retval, unsafe_load(list, offset))
        offset += 1
    end

    # ---------------------------------------------------------------------------------------------
    # Free the memory.
    delete_array(list)

    return retval
end


text:
function time_stepping_works(arch, ft)
    Nx, Ny, Nz = 16, 16, 16
    Lx, Ly, Lz = 1, 2, 3
    Δt = 1

    model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), arch=arch, float_type=ft)
    time_step!(model, 1, Δt)

    # Just testing that no errors/crashes happen when time stepping.
    return true
end

function run_first_AB2_time_step_tests(arch, ft)
    Nx, Ny, Nz = 16, 16, 16
    Lx, Ly, Lz = 1, 2, 3
    Δt = 1

    model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), arch=arch, float_type=ft)

    add_ones(args...) = 1.0
    model.forcing = Forcing(nothing, nothing, nothing, add_ones, nothing)

    time_step!(model, 1, Δt)

    # Test that GT = 1 after first time step and that AB2 actually reduced to forward Euler.
    @test all(model.G.Gu.data .≈ 0)
    @test all(model.G.Gv.data .≈ 0)
    @test all(model.G.Gw.data .≈ 0)
    @test all(model.G.GT.data .≈ 1.0)
    @test all(model.G.GS.data .≈ 0)

    return nothing
end


text:
module IDS

export append, tupleid, base, combine, append, increment, defID, singletonid

# # IDs
# IDs are used primarily in Omega to give identities to random variables
# Conceptually, an ID is a tuple of integers, but other data structures may be used

# # Generic Interface
"`base(::Type{T}, i)` singleton (`i`,) of collection type `T` "
function base end

"`combine(a, b)` Uniquely combine ids `a` and `b`"
function combine end

"`append(a, b)` append `b` to the end of `a`, like `append!` but functional"
function append end

"`Increment(id::ID)` the id"
function increment end

include("uintid.jl")
include("tupleid.jl")
include("vectorid.jl")
include("conversions.jl")

"default ID type"
defID(args...) = VectorID

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Runs through the entire ACAS X pipeline from jsons to test output. Includes:
processing data from json, creating a filtered dataset, running acasx_mc_tree.
Intended use is to duplicate a small test dataset (e.g., DASC) and perform a test run.

Usage:

ptest = pipelinetest() #produces data under GrammarExpts/test/PipelineTest

#inspect Dataset and processing, if desired

cleanup(ptest) #remove created artifacts, except results
"""
module PipelineTest

export pipelinetest, cleanup

import Compat.ASCIIString
using GrammarExpts, CASJson2DataFrame, FilterRecords
using DataFrameSets
using DataFrames
using Datasets
using Base.Test

using ACASX_MC, ACASX_MC_Tree
using FilterNMACInfo

const DIR = dirname(@__FILE__)
const DATADIR = joinpath(DIR, "..", "..", "..", "data")
const DASC_JSON = joinpath(DATADIR, "dasc", "json")
const TESTDIR = joinpath(DIR, "..", "..", "..", "test", "PipelineTest")

type PTest
  dirs::Vector{ASCIIString}
end

function pipelinetest(orig_dataname::AbstractString="dasc",
    orig_jsondir::AbstractString=DASC_JSON,
    test_dataname::AbstractString="exampledata")
  test_csvdir = joinpath(DATADIR, test_dataname, "csv")
  test_datadir = Pkg.dir("Datasets", "data", test_dataname)

  orig_datafilt = orig_dataname * "filt"
  test_datafilt = test_dataname * "filt"
  test_datafiltdir = Pkg.dir("Datasets", "data", test_datafilt)

  #process jsons
  script_base(orig_jsondir, test_csvdir, test_datadir;
    fromjson=true, correct_coc=true)

  #check data
  dset1 = dataset(orig_dataname)
  dset2 = dataset(test_dataname)
  for i in length(dset1)
    @test getrecords(dset1, i) == getrecords(dset2, i)
  end

  #test run
  outdir = joinpath(TESTDIR, test_dataname)
  acasx_mc(outdir=outdir, data=test_dataname)

  #generate filtered data
  remove_cpa(test_dataname, test_datafiltdir; t_min=35, n_before=5)

  #check filt data
  dset1 = dataset(orig_datafilt)
  dset2 = dataset(test_datafilt)
  for i in length(dset1)
    @test getrecords(dset1, i) == getrecords(dset2, i)
  end

  #test run on filt
  outdir = joinpath(TESTDIR, test_datafilt)
  acasx_mc_tree(outdir=outdir, data=test_datafilt)

  #for cleanup
  intermediatedir = joinpath(DATADIR, test_dataname)
  PTest([intermediatedir, test_datadir, test_datafiltdir])
end

function cleanup(ptest::PTest)
  for d in ptest.dirs
    rm(d; recursive=true)
  end
end

end #module


text:
using LinearAlgebra, ForwardDiff
using Convex, SCS, ECOS

n = 3
m = 3
k = n + m + 1

c = zeros(n)
A = Diagonal(ones(n))
b = zeros(n)

"Convex.jl"
x = Variable(n)
prob = minimize(c' * x)
prob.constraints += norm(b - A * x) <= 0.0
@time solve!(prob, ECOS.Optimizer)

@show prob.status
@show x.value
@show prob.constraints[1].dual
prob.optval

Q = Array([zeros(n, n) A' c;
     -A zeros(m, m) b;
     -c' -b'      0.0])
Q_vec = vec(Q)

function F(z, Q_vec)
    ũ = z[1:k]
    u = z[k .+ (1:k)]
    v = z[2 * k .+ (1:k)]

    [(I + reshape(Q_vec, k, k)) * ũ - (u + v);
     u - P_soc(ũ - v);
     ũ - u]
end

z = rand(3k)
F(z, Q_vec)
Fz(x) = F(x, Q_vec)
FQ(x) = F(z, x)
Jz = ForwardDiff.jacobian(Fz, z)
JQ = ForwardDiff.jacobian(FQ, Q_vec)

norm((Jz' * Jz) \ (Jz' * JQ), Inf)


text:
"""
    similarity(sender::Agent, receiver::Agent)

Compute the ordinal similarity of two agents.
Argument names are chosen to match assimilate! function.
"""
function similarity(sender::Agent, receiver::Agent)
    absolute_difference = sum(abs.(sender.opinions .- receiver.opinions))
    highest_possible_difference = 2 * length(sender.opinions)
    return 1 - (absolute_difference / highest_possible_difference)
end

"""
    assimilate!(sender::Agent, receiver::Agent)

Change one of the receiver's opinions to sender's opinion.
"""
function assimilate!(sender::Agent, receiver::Agent)
    i = Random.rand(1:length(sender.opinions))
    receiver.opinions[i] = StatsBase.mean([sender.opinions[i], receiver.opinions[i]])
    return receiver
end