text:
# tests for curvilinear meshes

function test_curvilinear(mesh, sbp, eqn, opts)

  @assert mesh.sbpface.stencilsize == mesh.numNodesPerElement
  Tdim = mesh.dim
  sbpface = mesh.sbpface

  stencilsize = mesh.sbpface.stencilsize
#  ncontrib = mesh.numTypePerElement[end-1]

  # collect E for all the elements
  E_all = zeros(stencilsize, stencilsize, Tdim, mesh.numEl)
  contrib_count = zeros(Int, mesh.numEl)

  PL = zeros(stencilsize, stencilsize)
  PR = zeros(PL)

  Pnbr = zeros(mesh.numNodesPerFace, mesh.numNodesPerFace)

  println("checking 1^T E 1 = 0")
  # this is a periodic mesh, so we only need to look at interfaces
  for i=1:mesh.numInterfaces
    iface = mesh.interfaces[i]
    contrib_count[iface.elementL] += 1
    contrib_count[iface.elementR] += 1

    # get R_complete, B
    fill!(PL, 0.0)
    fill!(PR, 0.0)
    permMatrix!(sbpface.perm[:, iface.faceL], PL)
    permMatrix!(sbpface.perm[:, iface.faceR], PR)

    RL_complete = sbpface.interp.'*PL
    RR_complete = sbpface.interp.'*PR
    B = diagm(sbpface.wface)

    for dim = 1:Tdim  # loop over parametric directions

      # get N
      NL_dim = diagm(vec(mesh.nrm_face[dim, :, i]))

      # the node order is different for elementR
      permMatrix!(sbpface.nbrperm[:, iface.orient], Pnbr)
      NR_dim = -diagm(Pnbr*vec(mesh.nrm_face[dim, :, i]))

      # compute E in current direction and accumulate into E_all
      E_all[:, :, dim, iface.elementL] += RL_complete.'*NL_dim*B*RL_complete
      E_all[:, :, dim, iface.elementR] += RR_complete.'*NR_dim*B*RR_complete
    end
  end

  Ivec = ones(stencilsize)
  for i=1:mesh.numEl
    @test ( contrib_count[i] )== 3
    for dim=1:Tdim
      E = E_all[:, :, dim, i]
      val = Ivec.'*E*Ivec

      @test isapprox( val[1], 0.0) atol=1e-12
      @test isapprox( norm(E - E.'), 0.0) atol=1e-12
    end
  end

  # test S
  println("checking S and D")
  S = zeros(stencilsize, stencilsize, Tdim)

  for i=1:mesh.numEl
    calcSCurvilinear(sbp, mesh.dxidx[:, :, :, i], S)
    Qxi = sbp.Q[:, :, 1]
    Qeta = sbp.Q[:, :, 2]
    Hinv = inv(diagm(sbp.w./mesh.jac[:, i]))

    for dim=1:Tdim
      @test isapprox( norm(S[:, :, dim] + S[:, :, dim].'), 0.0) atol=1e-12
      D = Hinv*(S[:, :, dim] + 0.5*E_all[:, :, dim, i])
      D1 = D*Ivec
      @test isapprox( norm(D1), 0.0) atol=1e-12
#=
      metrics_xi_x = diagm(vec(mesh.dxidx[1, dim, :, i]))
      metrics_eta_x = diagm(vec(mesh.dxidx[2, dim, :, i]))
      A_xi = metrics_xi_x*Qxi
      A_eta = metrics_eta_x*Qeta

      S2 = 0.5*A_xi - 0.5*A_xi.' + 0.5*A_eta - 0.5*A_eta.'
      println("S2 = \n", S2)
      D = Hinv*(S2 + 0.5*E_all[:, :, dim, i])
      D2 = D*Ivec
      @test isapprox( norm(D2), 0.0) atol=1e-12
=#
    end

  end

  println("checking theorem 6")
  D = zeros(stencilsize, stencilsize, Tdim)
  E = zeros(stencilsize, stencilsize, Tdim)
  for i=1:mesh.numEl
    Hinv = diagm(1./sbp.w)

    calcECurvilinear(sbp, mesh.dxidx[:, :, :, i], E)
    calcDCurvilinear(sbp, mesh.dxidx[:, :, :, i], D)
    for dim=1:Tdim
      lhs = D[:, :, dim]*Ivec
      rhs = Hinv*(E[:, :, dim] - E_all[:, :, dim, i])*Ivec
      @test isapprox( norm(lhs - rhs), 0.0) atol=1e-12
    end
  end
#=
  val = 0.0
  @time for i=1:mesh.numEl
    calcSCurvilinear(sbp, unsafe_view(mesh.dxidx, :, :, :, i), S)

    for dim=1:Tdim
      for i=1:stencilsize
        for j=1:stencilsize
          val += S[j, i, dim]
        end
      end
    end

  end
=#

  return nothing
end

infile = "input_vals_curve.jl"
add_func2!(EulerTests, test_curvilinear, infile, [TAG_CURVILINEAR, TAG_ESS, TAG_SHORTTEST])


text:
#isdefined(Main, :InterpolationKernels) || include("../src/InterpolationKernels.jl")

module InterpolationKernelsTests

#include("../src/interp.jl")
using InterpolationKernels
using InterpolationKernels:
    Bound, Open, Closed, Support,
    LeftAnchoredSupport, RightAnchoredSupport,
    signabs, floating_point_type,
    infimum, supremum, support, brief,
    compute_weights, generic_compute_weights,
    compute_offset_and_weights, generic_compute_offset_and_weights
import InterpolationKernels: support

using Test

shortname(::Nothing) = ""
shortname(m::RegexMatch) = m.captures[1]
shortname(::Type{T}) where {T} = shortname(string(T))
shortname(str::AbstractString) =
    shortname(match(r"([_A-Za-z][_A-Za-z0-9]*)([({]|$)", str))

const FLOATS = (Float32, Float64, BigFloat)

#------------------------------------------------------------------------------
# Non-specialized, non-optimized spline-based kernels to serve as a reference.

abstract type ReferenceSpline{T,S} <: Kernel{T,S} end

struct Spline1{T,M,N1} <: ReferenceSpline{T,1}
    vals::NTuple{M,T}
    c1::NTuple{N1,T}
end

struct Spline1Prime{T,M,N1} <: ReferenceSpline{T,1}
    vals::NTuple{M,T}
    c1::NTuple{N1,T}
end

struct Spline2{T,M,N1} <: ReferenceSpline{T,2}
    vals::NTuple{M,T}
    c1::NTuple{N1,T}
end

struct Spline2Prime{T,M,N1} <: ReferenceSpline{T,2}
    vals::NTuple{M,T}
    c1::NTuple{N1,T}
end

struct Spline3{T,M,N1,N2} <: ReferenceSpline{T,3}
    vals::NTuple{M,T}
    c1::NTuple{N1,T}
    c2::NTuple{N2,T}
end

struct Spline3Prime{T,M,N1,N2} <: ReferenceSpline{T,3}
    vals::NTuple{M,T}
    c1::NTuple{N1,T}
    c2::NTuple{N2,T}
end

struct Spline4{T,M,N1,N2} <: ReferenceSpline{T,4}
    vals::NTuple{M,T}
    c1::NTuple{N1,T}
    c2::NTuple{N2,T}
end

struct Spline4Prime{T,M,N1,N2} <: ReferenceSpline{T,4}
    vals::NTuple{M,T}
    c1::NTuple{N1,T}
    c2::NTuple{N2,T}
end

#Base.convert(::Type{K}, ker::K) where {K<:ReferenceSpline} = ker

Base.values(ker::ReferenceSpline) = ker.vals

@static if VERSION < v"1.4"
    # Method `evalpoly` appeared in Julia 1.4.
    evalpoly(x::Real, c::NTuple{1,Real}) =
        convert(promote_type(typeof(x), typeof(c[1])), c[1])
    evalpoly(x::Real, c::NTuple{2,Real}) = c[1] + x*c[2]
    evalpoly(x::Real, c::NTuple{3,Real}) = c[1] + x*(c[2] + x*c[3])
    evalpoly(x::Real, c::NTuple{4,Real}) = c[1] + x*(c[2] + x*(c[3] + x*c[4]))
end

function (ker::Spline1{T})(_x::Real) where {T}
    x = T(_x)
    half = T(1)/2
    if -half ≤ x < half
        return evalpoly(abs(x), ker.c1)
    else
        return zero(T)
    end
end

function (ker::Spline1Prime{T})(_x::Real) where {T}
    x = T(_x)
    half = T(1)/2
    if -half ≤ x < half
        return sign(x)*evalpoly(abs(x), ker.c1)
    else
        return zero(T)
    end
end

function (ker::Spline2{T})(_x::Real) where {T}
    x = T(_x)
    abs_x = abs(x)
    if abs_x ≥ 1
        return zero(T)
    else
        return evalpoly(abs_x, ker.c1)
    end
end

function (ker::Spline2Prime{T})(_x::Real) where {T}
    x = T(_x)
    abs_x = abs(x)
    if abs_x ≥ 1
        return zero(T)
    else
        return sign(x)*evalpoly(abs_x, ker.c1)
    end
end

function (ker::Spline3{T})(_x::Real) where {T}
    x = T(_x)
    abs_x = abs(x)
    if abs_x ≥ T(3)/2
        return zero(T)
    elseif abs_x ≤ T(1)/2
        return evalpoly(abs_x, ker.c1)
    else
        return evalpoly(abs_x, ker.c2)
    end
end

function (ker::Spline3Prime{T})(_x::Real) where {T}
    x = T(_x)
    abs_x = abs(x)
    if abs_x ≥ T(3)/2
        return zero(T)
    elseif abs_x ≤ T(1)/2
        return sign(x)*evalpoly(abs_x, ker.c1)
    else
        return sign(x)*evalpoly(abs_x, ker.c2)
    end
end

function (ker::Spline4{T})(_x::Real) where {T}
    x = T(_x)
    abs_x = abs(x)
    if abs_x ≥ 2
        return zero(T)
    elseif abs_x ≤ 1
        return evalpoly(abs_x, ker.c1)
    else
        return evalpoly(abs_x, ker.c2)
    end
end

function (ker::Spline4Prime{T})(_x::Real) where {T}
    x = T(_x)
    abs_x = abs(x)
    if abs_x ≥ 2
        return zero(T)
    elseif abs_x ≤ 1
        return sign(x)*evalpoly(abs_x, ker.c1)
    else
        return sign(x)*evalpoly(abs_x, ker.c2)
    end
end

Base.adjoint(ker::Spline1) =
    Spline1Prime(values(ker), derivepoly(ker.c1))

Base.adjoint(ker::Spline2) =
    Spline2Prime(values(ker), derivepoly(ker.c1))

Base.adjoint(ker::Spline3) =
    Spline3Prime(values(ker), derivepoly(ker.c1),  derivepoly(ker.c2))

Base.adjoint(ker::Spline4) =
    Spline4Prime(values(ker), derivepoly(ker.c1),  derivepoly(ker.c2))

# FIXME: can be @generated
derivepoly(c::NTuple{1,T}) where {T<:AbstractFloat} = (zero(T),)
derivepoly(c::NTuple{2,T}) where {T<:AbstractFloat} = (c[2],)
derivepoly(c::NTuple{3,T}) where {T<:AbstractFloat} = (c[2], 2c[3])
derivepoly(c::NTuple{4,T}) where {T<:AbstractFloat} = (c[2], 2c[3], 3c[4])

bspline1(T::Type{<:AbstractFloat} = Float64) =
    Spline1((), (one(T),))

bspline2(T::Type{<:AbstractFloat} = Float64) =
    Spline2((), (one(T), -one(T)))

bspline3(T::Type{<:AbstractFloat} = Float64) =
    Spline3((),
            (T(3//4), T(0), T(-1)),
            (T(9//8), T(-3//2), T(1//2)))

bspline4(T::Type{<:AbstractFloat} = Float64) =
    Spline4((),
            (T(2//3), T(0), T(-1), T(1//2)),
            (T(4//3), T(-2), T( 1), T(-1//6)))

# A family of cubic splines parameterized by a = ker'(1) and b = ker(1).
generic_cubic_spline(a::Real, b::Real) = genericcubic_spline(Float64, a, b)
generic_cubic_spline(T::Type{<:AbstractFloat}, _a::Real, _b::Real) = begin
    a, b = T(_a), T(_b)
    # p1(x) = c0 + c1*x + c2*x^2
    # p2(x) = (c3 + c4*x)*(x - 2)^2
    # c0 = 1 - 2b
    # c1 = 9b - a - 3
    # c2 = 2 + a - 6b
    # c3 = -a - b
    # c4 = a + 2b
    # Spline4((a, b),
    #         (c0, zero(T), c1, c2),
    #         (4c3, 4(c4 - c3), c3 - 4c4, c4))
    Spline4((a, b),
            (1 - 2b, zero(T), 9b - a - 3, 2 + a - 6b),
            (-4a - 4b, 8a + 12b,  -5a - 9b,  a + 2b))
end

# A family of cardinal cubic splines.
keys_spline(a::Real) = keys_spline(Float64, a)
keys_spline(T::Type{<:AbstractFloat}, _a::Real) = begin
    a = T(_a)
    Spline4((a,),
            (one(T), zero(T), -a - 3, a + 2),
            (-4a, 8a, -5a, a))
end

catmull_rom_spline(T::Type{<:AbstractFloat} = Float64) = keys_spline(T, T(-1//2))

mitchell_netravali_spline(b::Real, c::Real) =
    mitchell_netravali_spline(Float64, b, c)

mitchell_netravali_spline(T::Type{<:AbstractFloat}, _b::Real, _c::Real) = begin
    b, c = T(_b), T(_c)
    Spline4((b,c),
            ((6 - 2b)/6, zero(T), (-18 + 12b + 6c)/6, (12 - 9b - 6c)/6),
            ((8b + 24c)/6, (-12b - 48c)/6, (6b + 30c)/6, (-b - 6c)/6))
end

kernel_model(ker::BSpline{1}) = bspline1(eltype(ker))
kernel_model(ker::BSplinePrime{1}) = bspline1(eltype(ker))'

kernel_model(ker::BSpline{2}) = bspline2(eltype(ker))
kernel_model(ker::BSplinePrime{2}) = bspline2(eltype(ker))'

kernel_model(ker::BSpline{3}) = bspline3(eltype(ker))
kernel_model(ker::BSplinePrime{3}) = bspline3(eltype(ker))'

kernel_model(ker::BSpline{4}) = bspline4(eltype(ker))
kernel_model(ker::BSplinePrime{4}) = bspline4(eltype(ker))'

kernel_model(ker::CatmullRomSpline) = catmull_rom_spline(eltype(ker))
kernel_model(ker::CatmullRomSplinePrime) = catmull_rom_spline(eltype(ker))'

kernel_model(ker::CubicSpline) = generic_cubic_spline(eltype(ker), values(ker)...)
#kernel_model(ker::CubicSplinePrime) = cubic_spline(eltype(ker), values(ker)...)'

kernel_model(ker::CardinalCubicSpline) = keys_spline(eltype(ker), values(ker)...)
kernel_model(ker::CardinalCubicSplinePrime) = keys_spline(eltype(ker), values(ker)...)'

kernel_model(::Kernel) = nothing

#------------------------------------------------------------------------------
# Types needed to check other supprot than the symmetric ones.

lower_bound_type(sup::Support) = lower_bound_type(typeof(sup))
lower_bound_type(::Type{<:Support{T,S,L,R}}) where {T,S,L,R} = L

upper_bound_type(sup::Support) = upper_bound_type(typeof(sup))
upper_bound_type(::Type{<:Support{T,S,L,R}}) where {T,S,L,R} = R

struct LeftAnchoredKernel{T,S,L<:Bound,R<:Bound,K<:Kernel{T,S}} <: Kernel{T,S}
    a::T
    ker::K
end

LeftAnchoredKernel(ker::K) where {T,S,K<:Kernel{T,S}} = begin
    sup = support(ker)
    L = lower_bound_type(sup)
    R = upper_bound_type(sup)
    return LeftAnchoredKernel{T,S,L,R,K}(infimum(sup), ker)
end

struct RightAnchoredKernel{T,S,L<:Bound,R<:Bound,K<:Kernel{T,S}} <: Kernel{T,S}
    b::T
    ker::K
end

RightAnchoredKernel(ker::K) where {T,S,K<:Kernel{T,S}} = begin
    sup = support(ker)
    L = lower_bound_type(sup)
    R = upper_bound_type(sup)
    return RightAnchoredKernel{T,S,L,R,K}(supremum(sup), ker)
end

support(ker::LeftAnchoredKernel{T,S,L,R}) where {T,S,L,R} =
    LeftAnchoredSupport{T,S,L,R}(ker.a)

support(ker::RightAnchoredKernel{T,S,L,R}) where {T,S,L,R} =
    RightAnchoredSupport{T,S,L,R}(ker.b)

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

raw_type(ker::Kernel) = raw_type(typeof(ker))
for K in (:BSpline,              :BSplinePrime,
          :CubicSpline,   :CubicSplinePrime,
          :CardinalCubicSpline,  :CardinalCubicSplinePrime,
          :CatmullRomSpline,     :CatmullRomSplinePrime,
          :LanczosKernel,        :LanczosKernelPrime)
    @eval raw_type(::Type{<:$K}) = $K
end
#raw_type(::Type{BSpline{S,T}}) where {S,T} = BSpline{S}
#raw_type(::Type{BSplinePrime{S,T}}) where {S,T} = BSplinePrime{S}

has_size_parameter(ker::Kernel) = has_size_parameter(typeof(ker))
has_size_parameter(::Type{<:Kernel}) = false
has_size_parameter(::Type{<:BSpline}) = true
has_size_parameter(::Type{<:BSplinePrime}) = true
has_size_parameter(::Type{<:LanczosKernel}) = true
has_size_parameter(::Type{<:LanczosKernelPrime}) = true

for (K, Kp) in ((:BSpline,             :BSplinePrime),
                (:CubicSpline,  :CubicSplinePrime),
                (:CardinalCubicSpline, :CardinalCubicSplinePrime),
                (:CatmullRomSpline,    :CatmullRomSplinePrime),
                (:LanczosKernel,       :LanczosKernelPrime))
    if K === :BSpline || K === :LanczosKernel
        @eval begin
            derivative_type(::Type{$K}) = $Kp
            derivative_type(::Type{$K{S}}) where {S} = $Kp{S}
            derivative_type(::Type{$K{S,T}}) where {S,T} = $Kp{S,T}
            primitive_type(::Type{$Kp}) = $K
            primitive_type(::Type{$Kp{S}}) where {S} = $K{S}
            primitive_type(::Type{$Kp{S,T}}) where {S,T} = $K{S,T}
        end
    else
        @eval begin
            derivative_type(::Type{$K}) = $Kp
            derivative_type(::Type{$K{T}}) where {T} = $Kp{T}
            primitive_type(::Type{$Kp}) = $K
            primitive_type(::Type{$Kp{T}}) where {T} = $K{T}
        end
    end
end

is_C¹_continuous(ker::Kernel) = is_C¹_continuous(typeof(ker))
is_C¹_continuous(::Type{<:Kernel}) = true
is_C¹_continuous(::Type{<:BSpline{S}}) where {S} = (S ≥ 3)

derivative_type(ker::Kernel) = derivative_type(typeof(ker))
primitive_type(ker::Kernel) = primitive_type(typeof(ker))

derivative_type(::Type{<:Kernel}) = Any
primitive_type(::Type{<:Kernel}) = Any

function simple_derivative(f::Kernel, x::Real, dx::Real)
    x1 = x - dx
    x2 = x + dx
    return (f(x1) - f(x2))/(x1 - x2)
end

function simple_derivative(f::Kernel{T},
                           x::AbstractArray{T},
                           dx::Real = sqrt(eps(T))) where {T}
    df = similar(x)
    for i in eachindex(x, df)
        df[i] = simple_derivative(f, x[i], dx)
    end
    return df
end

@testset "Interpolation kernels" begin
    @testset "Utilities" begin
        for x in (UInt(3), Int32(-8), -2.1, 0.0f0, 3.7)
            @test signabs(x) === (sign(x), abs(x))
            @test floating_point_type(x) === floating_point_type(typeof(x))
        end
        io = IOBuffer()
        ker = LanczosKernel{4,Float32}()
        show(io, MIME("text/plain"), ker)
        str1 = take!(io)
        show(io, ker)
        str2 = take!(io)
        @test str1 == str2
    end

    @testset "Constructors" begin
        for S in 1:4
            @test eltype(BSpline{S}()) === Float64
            @test eltype(BSpline{S,Float32}()) === Float32
            @test eltype(BSplinePrime{S}()) === Float64
            @test eltype(BSplinePrime{S,Float32}()) === Float32
        end
        @test eltype(CubicSpline(1,0)) === Float64
        @test eltype(CubicSpline(1,0f0)) === Float32
        @test eltype(CubicSpline(1,0.0)) === Float64
        @test eltype(CubicSpline{Float32}(1,0.0)) === Float32

        @test eltype(CubicSplinePrime(1,0)) === Float64
        @test eltype(CubicSplinePrime(1,0f0)) === Float32
        @test eltype(CubicSplinePrime(1,0.0)) === Float64
        @test eltype(CubicSplinePrime{Float32}(1,0.0)) === Float32

        @test eltype(CardinalCubicSpline(0)) === Float64
        @test eltype(CardinalCubicSpline(0f0)) === Float32
        @test eltype(CardinalCubicSpline(0.0)) === Float64
        @test eltype(CardinalCubicSpline{Float32}(0.0)) === Float32

        @test eltype(CardinalCubicSplinePrime(0)) === Float64
        @test eltype(CardinalCubicSplinePrime(0f0)) === Float32
        @test eltype(CardinalCubicSplinePrime(0.0)) === Float64
        @test eltype(CardinalCubicSplinePrime{Float32}(0.0)) === Float32

        @test eltype(CatmullRomSpline()) === Float64
        @test eltype(CatmullRomSpline{Float32}()) === Float32
        @test eltype(CatmullRomSplinePrime()) === Float64
        @test eltype(CatmullRomSplinePrime{Float32}()) === Float32

        for S in (4,6)
            @test eltype(LanczosKernel{S}()) === Float64
            @test eltype(LanczosKernel{S,Float32}()) === Float32
            @test eltype(LanczosKernelPrime{S}()) === Float64
            @test eltype(LanczosKernelPrime{S,Float32}()) === Float32
        end
        @test_throws ArgumentError LanczosKernel{3}()
    end

    @testset "Values" begin
        for ker in (BSpline{1}(),
                    BSplinePrime{1}(), # not C¹ continuous, but the model works
                    BSpline{2}(),
                    BSplinePrime{2}(), # not C¹ continuous, but the model works
                    BSpline{3}(),
                    BSplinePrime{3}(),
                    BSpline{4}(),
                    BSplinePrime{4}(),
                    CubicSpline(-1/3,1/5),
                    CubicSplinePrime(-2/3,1/7),
                    CardinalCubicSpline(-1/3),
                    CardinalCubicSplinePrime(-2/3),
                    CatmullRomSpline(),
                    CatmullRomSplinePrime(),
                    LanczosKernel{4}(),
                    LanczosKernelPrime{4}())
            mdl = kernel_model(ker)
            if mdl === nothing
                continue
            end
            s = 0.123 # step chosen so as to avoid integral values of x
            sup = support(ker)
            x = infimum(sup)-s:s:supremum(sup)+s
            @test ker.(x) ≈ mdl.(x)
        end
    end

    # Check the ability of the generic cubic spline to emulate other kernels.
    @testset "Equivalences" begin
        s = 0.1
        # Emulation of the cubic B-spline.
        let ker = BSpline{4}(),
            mdl = CubicSpline(-1/2,1/6)
            sup = support(ker)
            x = infimum(sup)-s:s:supremum(sup)+s
            @test mdl.(x) ≈ ker.(x)
            @test mdl'.(x) ≈ ker'.(x)
        end
        # Emulation of the family of cardinal cubic splines.
        for a in (-1/3, 0.1, -1/2)
            ker = CardinalCubicSpline(a)
            mdl = CubicSpline(a, 0)
            sup = support(ker)
            x = infimum(sup)-s:s:supremum(sup)+s
            @test mdl.(x) ≈ ker.(x)
            @test mdl'.(x) ≈ ker'.(x)
        end
        # Emulation of the Catmull-Rom kernel.
        let ker = CatmullRomSpline(),
            mdl = CubicSpline(-1/2, 0)
            sup = support(ker)
            x = infimum(sup)-s:s:supremum(sup)+s
            @test mdl.(x) ≈ ker.(x)
            @test mdl'.(x) ≈ ker'.(x)
        end
        # Emulation of the family of Mitchell & Netravali kernels.
        @test MitchellNetravaliSpline() === MitchellNetravaliSpline(1/3,1/3)
        @test MitchellNetravaliSpline{Float32}() === MitchellNetravaliSpline{Float32}(1/3,1/3)
        for (b,c) in ((1/3,1/3), (0.1,0.7), (-0.2,0.9))
            let ker = MitchellNetravaliSpline(b, c),
                mdl = mitchell_netravali_spline(b, c)
                sup = support(ker)
                x = infimum(sup)-s:s:supremum(sup)+s
                @test mdl.(x) ≈ ker.(x)
                @test mdl'.(x) ≈ ker'.(x)
            end
        end
    end

    @testset "Conversion" begin
        for ker in (BSpline{1}(),
                    BSplinePrime{1}(),
                    BSpline{2}(),
                    BSplinePrime{2}(),
                    BSpline{3}(),
                    BSplinePrime{3}(),
                    BSpline{4}(),
                    BSplinePrime{4}(),
                    CatmullRomSpline(),
                    CatmullRomSplinePrime(),
                    CubicSpline(-1/3,1/5),
                    CubicSplinePrime(-2/3,1/7),
                    CardinalCubicSpline(-1/3),
                    CardinalCubicSplinePrime(-2/3),
                    LanczosKernel{4}(),
                    LanczosKernelPrime{4}())
            K = raw_type(ker)
            @test eltype(ker) === Float64
            @test Kernel(ker) === ker
            @test Kernel{eltype(ker)}(ker) === ker
            @test convert(Kernel, ker) === ker
            @test convert(Kernel{eltype(ker)}, ker) === ker
            if has_size_parameter(ker)
                @test Kernel{eltype(ker),length(ker)}(ker) === ker
                @test convert(Kernel{eltype(ker),length(ker)}, ker) === ker
            end
            for T in FLOATS
                if has_size_parameter(ker)
                    S = length(ker)
                    @test typeof(T(ker)) <: K{S,T}
                    if T !== BigFloat
                        @test T(ker) === K{S,T}(ker)
                        @test T(ker) === K{S,T}(values(ker)...)
                    end
                else
                    @test typeof(T(ker)) <: K{T}
                    if T !== BigFloat
                        @test T(ker) === K{T}(ker)
                        @test T(ker) === K{T}(values(ker)...)
                    end
                end
            end
        end
    end

    @testset "Methods" begin
        tol = 1e-14
        for (ker, len, nrml, card) in (
            (BSpline{1}(),                   1, true,  true),
            (BSplinePrime{1}(),              1, false, false),
            (BSpline{2}(),                   2, true,  true),
            (BSplinePrime{2}(),              2, false, false),
            (BSpline{3}(),                   3, true,  false),
            (BSplinePrime{3}(),              3, false, false),
            (BSpline{4}(),                   4, true,  false),
            (BSplinePrime{4}(),              4, false, false),
            (CatmullRomSpline(),             4, true,  true),
            (CatmullRomSplinePrime(),        4, false, false),
            (CubicSpline(-1/3,1/9),          4, true,  false),
            (CubicSplinePrime(-1/3,1/9),     4, false, false),
            (CardinalCubicSpline(-1),        4, true,  true),
            (CardinalCubicSplinePrime(-1/2), 4, false, false),
            (LanczosKernel{2}(),             2, false, true),
            (LanczosKernel{4}(),             4, false, true),
            (LanczosKernel{6}(),             6, false, true),
            (LanczosKernel{8}(),             8, false, true),
            (LanczosKernelPrime{4}(),        4, false, false),
            (LanczosKernelPrime{6}(),        6, false, false))
            name = brief(ker)
            @test isa(name, String)
            @testset "$name" begin
                @test isa(summary(ker), String)
                @test isnormalized(ker) == nrml
                @test iscardinal(ker) == card
                @test length(ker) == len
                @test length(typeof(ker)) == len
                @test shortname(summary(ker)) == shortname(typeof(ker))
                if iscardinal(ker)
                    @test ker(0) == 1
                    @test maximum(abs.(ker([-3,-2,-1,1,2,3]))) ≤ tol
                end
                for κ in (ker,
                          RightAnchoredKernel(ker),
                          LeftAnchoredKernel(ker))
                    sup = support(κ)
                    @test eltype(sup) === eltype(κ)
                    @test length(sup) === length(κ)
                    @test supremum(sup) == infimum(sup) + length(sup)
                end

                # Check compute_weights.  Compare values computed by optimal,
                # generic and reference method.
                k = (length(ker) + 1) >> 1
                err_gen, err_opt = 0.0, 0.0
                for t in (0.0, 0.1, 0.2, 0.3, 0.4)
                    wgt_ref = ntuple(i -> ker(t + k - i), length(ker))
                    wgt_gen = generic_compute_weights(ker, t)
                    wgt_opt = compute_weights(ker, t)
                    err_gen = max(err_gen, maximum(abs.(wgt_gen .- wgt_ref)))
                    err_opt = max(err_opt, maximum(abs.(wgt_opt .- wgt_ref)))
                end
                @test err_gen ≤ tol
                @test err_opt ≤ tol

                # Check compute_offset_and_weights.  Compare values computed by optimal, generic
                # and reference method.
                k = (length(ker) + 1) >> 1
                off_gen_err, wgt_gen_err = 0.0, 0.0
                off_opt_err, wgt_opt_err = 0.0, 0.0
                for x in (50.0, -30.9, 8.2, 11.3, -27.6)
                    # FIXME: These rules are only valid for symmetric kernels.
                    if isodd(length(ker))
                        round_x = round(x)
                        t = x - round_x
                        off_ref = round_x - (length(ker) + 1)/2
                    else
                        floor_x = floor(x)
                        t = x - floor_x
                        off_ref = floor_x - length(ker)/2
                    end
                    wgt_ref = ntuple(i -> ker(t + k - i), length(ker))
                    off_gen, wgt_gen = generic_compute_offset_and_weights(ker, x)
                    off_gen_err = max(off_gen_err, abs(off_gen - off_ref))
                    wgt_gen_err = max(wgt_gen_err, maximum(abs.(wgt_gen .- wgt_ref)))
                    off_opt, wgt_opt = compute_offset_and_weights(ker, x)
                    off_opt_err = max(off_opt_err, abs(off_opt - off_ref))
                    wgt_opt_err = max(wgt_opt_err, maximum(abs.(wgt_opt .- wgt_ref)))
                end
                @test off_gen_err == 0
                @test wgt_gen_err ≤ tol
                @test off_opt_err == 0
                @test wgt_opt_err ≤ tol
            end
        end
        let ker = BSpline{6}()
            @test isa(brief(ker), String)
            @test isa(brief(ker'), String)
        end
    end

    @testset "Derivatives" begin
        for ker in (BSpline{1}(),
                    BSpline{2}(),
                    BSpline{3}(),
                    BSpline{4}(),
                    CatmullRomSpline(),
                    LanczosKernel{4}(),
                    LanczosKernel{6}())
            der = ker'
            @test derivative_type(ker) === typeof(der)
            @test primitive_type(der) === typeof(ker)
            @test values(der) === values(ker)
            @test der === derivative_type(ker)(values(ker)...)
            @test ker === primitive_type(der)(values(der)...)

            s = 0.123 # step chosen so as to avoid integral values of x
            sup = support(ker)
            x = infimum(sup)+s:s:supremum(sup)-s
            @test der.(x) ≈ simple_derivative(ker, x)

        end
    end
end

end # module

ext: jl


text:
"""
Reserved Julia identifiers will be prepended with "_"
"""
const RESERVED_WORDS = ["begin", "while", "if", "for", "try", "return", "break", "continue",
                        "function", "macro", "quote", "let", "local", "global", "const", "do",
                        "struct", "module", "baremodule", "using", "import", "export", "end",
                        "else", "elseif", "catch", "finally", "true", "false"]

"""
    make_name_safe(name::Symbol) -> String
    make_name_safe(name::AbstractString) -> String
Return a valid Julia variable name, prefixed with "_" if the `name` is conflict with Julia's
reserved words.
"""
make_name_safe(name::AbstractString) = (name in RESERVED_WORDS) ? "_" * name : name
make_name_safe(name::Symbol) = make_name_safe(string(name))

"""
    make_symbol_safe(name) -> Symbol
Same as [`make_name_safe`](@ref), but return a Symbol.
"""
make_symbol_safe(x) = Symbol(make_name_safe(x))

"""
    translate(jlty::AbstractJuliaType, options=Dict())
Translate [`AbstractJuliaType`](@ref)s to Julia expressions.
"""
translate(jlty::AbstractJuliaType, options=Dict()) = jlty
function translate(jlty::JuliaUnknown, options=Dict())
    return error("hit a JuliaUnknown($(dumpobj(jlty.x))) when translating the type...")
end
translate(jlty::JuliaCvoid, options=Dict()) = :Cvoid
function translate(jlty::JuliaCbool, options=Dict())
    return get(options, "use_julia_bool", true) ? :Bool : :UInt8
end
translate(jlty::JuliaCuchar, options=Dict()) = :Cuchar
translate(jlty::JuliaCshort, options=Dict()) = :Cshort
translate(jlty::JuliaCushort, options=Dict()) = :Cushort
translate(jlty::JuliaCint, options=Dict()) = :Cint
translate(jlty::JuliaCuint, options=Dict()) = :Cuint
translate(jlty::JuliaClonglong, options=Dict()) = :Clonglong
translate(jlty::JuliaCulonglong, options=Dict()) = :Culonglong
translate(jlty::JuliaCintmax_t, options=Dict()) = :Cintmax_t
translate(jlty::JuliaCuintmax_t, options=Dict()) = :Cuintmax_t
translate(jlty::JuliaCfloat, options=Dict()) = :Cfloat
translate(jlty::JuliaCdouble, options=Dict()) = :Cdouble
translate(jlty::JuliaComplexF32, options=Dict()) = :ComplexF32
translate(jlty::JuliaComplexF64, options=Dict()) = :ComplexF64
translate(jlty::JuliaCptrdiff_t, options=Dict()) = :Cptrdiff_t
translate(jlty::JuliaCssize_t, options=Dict()) = :Cssize_t
translate(jlty::JuliaCsize_t, options=Dict()) = :Csize_t
translate(jlty::JuliaNoReturn, options=Dict()) = :(Union{})
translate(jlty::JuliaPtrCvoid, options=Dict()) = :(Ptr{Cvoid})
translate(jlty::JuliaCstring, options=Dict()) = :Cstring
translate(jlty::JuliaPtrUInt8, options=Dict()) = :(Ptr{UInt8})
translate(jlty::JuliaPtrPtrUInt8, options=Dict()) = :(Ptr{Ptr{UInt8}})
translate(jlty::JuliaAny, options=Dict()) = :Any
translate(jlty::JuliaRefAny, options=Dict()) = :(Ref{Any})
translate(jlty::JuliaCuint128, options=Dict()) = :UInt128
translate(jlty::JuliaCint128, options=Dict()) = :Int128
translate(jlty::JuliaCschar, options=Dict()) = :Int8
translate(jlty::JuliaClongdouble, options=Dict()) = :Float64  # FIXME: the actually type depends on compiler
translate(jlty::JuliaComplex, options=Dict()) = :ComplexF32  # FIXME: correct this
translate(jlty::JuliaChalf, options=Dict()) = :Float16
translate(jlty::JuliaCfloat16, options=Dict()) = :Float16
translate(jlty::JuliaCuint64_t, options=Dict()) = :UInt64
translate(jlty::JuliaCuint32_t, options=Dict()) = :UInt32
translate(jlty::JuliaCuint16_t, options=Dict()) = :UInt16
translate(jlty::JuliaCuint8_t, options=Dict()) = :UInt8
translate(jlty::JuliaCint64_t, options=Dict()) = :Int64
translate(jlty::JuliaCint32_t, options=Dict()) = :Int32
translate(jlty::JuliaCint16_t, options=Dict()) = :Int16
translate(jlty::JuliaCint8_t, options=Dict()) = :Int8
translate(jlty::JuliaCuintptr_t, options=Dict()) = :Csize_t
translate(jlty::JuliaCtm, options=Dict()) = :(Libc.TmStruct)
translate(jlty::JuliaCchar, options=Dict()) = :Cchar
translate(jlty::JuliaClong, options=Dict()) = :Clong
translate(jlty::JuliaCulong, options=Dict()) = :Culong
translate(jlty::JuliaCwchar_t, options=Dict()) = :Cwchar_t
translate(jlty::JuliaCFILE, options=Dict()) = :(Libc.FILE)

function translate(jlty::JuliaCpointer, options=Dict())
    is_jl_funcptr(jlty) && return translate(JuliaPtrCvoid(), options)
    jlptree = tojulia(getPointeeType(jlty.ref))
    if get(options, "always_NUL_terminated_string", false)
        is_jl_char(jlptree) && return :Cstring
        is_jl_wchar(jlptree) && return :Cwstring
    end
    return Expr(:curly, :Ptr, translate(jlptree, options))
end

function translate(jlty::JuliaCconstarray, options=Dict())
    n = getNumElements(jlty.ref)
    elty = translate(tojulia(getElementType(jlty.ref)), options)
    return :(NTuple{$n,$elty})
end

function translate(jlty::JuliaCincompletearray, options=Dict())
    elty = translate(tojulia(getElementType(jlty.ref)), options)
    return Expr(:curly, :Ptr, elty)
end

function translate(jlty::JuliaCtypedef, options=Dict())
    ids = get(options, "DAG_ids", Dict())
    ids_extra = get(options, "DAG_ids_extra", Dict())

    if haskey(ids, jlty.sym)
        return make_symbol_safe(jlty.sym)
    elseif haskey(ids_extra, jlty.sym)
        return translate(ids_extra[jlty.sym], options)
    elseif get(options, "opaque_func_arg_as_PtrCvoid", false)
        return translate(JuliaCvoid(), options)
    else
        return make_symbol_safe(jlty.sym)
    end
end

function translate(jlty::JuliaCenum, options=Dict())
    tags = get(options, "DAG_tags", Dict())
    # for now, we don't distinguish extra tags and ids, this may be improved in the future.
    # tags_extra = get(options, "DAG_tags_extra", Dict())
    ids_extra = get(options, "DAG_tags_extra", Dict())
    if haskey(tags, jlty.sym) || haskey(ids_extra, jlty.sym)
        return make_symbol_safe(jlty.sym)
    else
        # it could a local opaque tag-type
        return translate(JuliaCvoid(), options)
    end
end

function translate(jlty::JuliaCrecord, options=Dict())
    tags = get(options, "DAG_tags", Dict())
    # for now, we don't distinguish extra tags and ids, this may be improved in the future.
    # tags_extra = get(options, "DAG_tags_extra", Dict())
    ids_extra = get(options, "DAG_ids_extra", Dict())
    nested_tags = get(options, "nested_tags", Dict())
    if haskey(tags, jlty.sym) || haskey(ids_extra, jlty.sym)
        return make_symbol_safe(jlty.sym)
    else
        for (id, cursor) in nested_tags
            if is_same(jlty.cursor, cursor)
                @assert isempty(string(jlty.sym))
                return id
            end
        end
        # then it could be a local opaque tag-type
        return translate(JuliaCvoid(), options)
    end
end

translate(jlty::JuliaCfunction, options=Dict()) = translate(JuliaCvoid(), options)


text:
using Base.Test

@testset "linear interpolation base" begin
	@test getindex(LinearInterpolation, one(Float64)) == InterpolationSupport[(1.,1),(0.,2)]
	@test eltype(getindex(LinearInterpolation, one(Float32))) == InterpolationSupport{1, Float32, Int}
	@test getindex(LinearInterpolation, one(Int)) ==  InterpolationSupport[(1,1),(0,2)]
end;

@testset "nearest interpolation base" begin
	@test getindex(NearestInterpolation, one(Float64)) == InterpolationSupport{1, Float64, Int}[(1.,1)]
	@test getindex(NearestInterpolation, one(Float16)) == InterpolationSupport{1, Float16, Int}[(one(Float16),1)]
	@test getindex(NearestInterpolation, one(Int)) == InterpolationSupport{1, Int, Int}[(1,1)]
end;

max_stars_repo_name: JKrehl/Interpol

max_issues_repo_name: JKrehl/Interpol

max_forks_repo_name: JKrehl/Interpol


text:
module LyceumCore

using Base: @pure, @_inline_meta

using StaticNumbers

export @mustimplement, propertytype, argerror
include("util.jl")

export Maybe
export TupleN
export AbsArr, AbsMat, AbsVec
export AbsNestedArr, AbsSimilarNestedArr
export StaticOrVal
include("aliases.jl")

export front, tail, tuple_split, tuple_length, tuple_prod, tuple_minimum
include("tuple.jl")

end # module


text:
"""
    fdweights(z, x, m)

Calculates FD (finite difference) weights for arbitrarily spaced nodes in 1D.

Return an array of size m+1 by length(x) containing, in successive rows, the
weights for derivatives 0, 1, ..., m. The zero (m = 0) derivative corresponds
to interpolation.

# Arguments
- `z`: location at which to approximate the derivative (may but need to be a grid point)
- `x`: vector of x-coordinates (grid points; distinct, otherwise arbitrary)
- `m`: highest derivative for which weights are sought (default value is 0)

# Examples
```julia-repl
julia> fdweights(0, -2:2, 6)
5×5 adjoint(::Matrix{Float64}) with eltype Float64:
 -0.0         0.0        1.0   0.0       -0.0
  0.0833333  -0.666667   0.0   0.666667  -0.0833333
 -0.0833333   1.33333   -2.5   1.33333   -0.0833333
 -0.5         1.0        0.0  -1.0        0.5
  1.0        -4.0        6.0  -4.0        1.0
```
"""
function fdweights(z, x, m = 0)
    n = length(x)
    if m >= n
        m = n - 1 # set to highest possible derivative
    end

    c1 = 1
    c4 = x[1] - z
    C = zeros(n, m+1)
    C[1, 1] = 1

    for i = 1:n-1
        i1 = i + 1
        mn = min(i, m)
        c2 = 1
        c5 = c4
        c4 = x[i1] - z
        for j = 0:i-1
            j1 = j + 1
            c3 = x[i1] - x[j1]
            c2 = c2 * c3
            if j == i-1
                for s = mn:-1:1
                    s1 = s + 1
                    C[i1, s1] = c1 * (s * C[i1-1, s1-1] - c5 * C[i1-1, s1]) / c2
                end
                C[i1, 1] = -c1 * c5 * C[i1-1, 1] / c2
            end
        for s = mn:-1:1
            s1 = s + 1
            C[j1, s1] = (c4 * C[j1, s1] - s*C[j1, s1-1]) / c3
        end
        C[j1, 1] = c4 * C[j1,1] / c3
        end
        c1 = c2
    end
    return C'
end

max_stars_repo_name: williamhunter/fornberg

max_issues_repo_name: williamhunter/fornberg

max_forks_repo_name: williamhunter/fornberg


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -25.420040452488465 <= q <= 21.65452636072683)

                     add_NL_constraint(m, :(softplus(-0.9974376164994152 + 0.7314791518637902*softplus(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4])) + 0.4470938857218645*softplus(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4])) + 0.10010754470216687*softplus(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4])) + -0.8594790184444285*softplus(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4])) + -0.4658863904521331*softplus(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))) + softplus(0.15396139316976853 + -0.4228510094486002*softplus(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4])) + -0.8035341866400141*softplus(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4])) + -0.07617073018159282*softplus(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4])) + 0.3940887684836425*softplus(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4])) + 0.6667343490308526*softplus(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))) + softplus(0.9665409577623398 + -0.449175249339572*softplus(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4])) + -0.2640107596849468*softplus(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4])) + 0.7284419447850534*softplus(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4])) + 0.5662248826250997*softplus(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4])) + -0.21454884110410477*softplus(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))) + softplus(-0.1582635851091898 + 0.5103556986884796*softplus(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4])) + -0.6017569427714062*softplus(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4])) + -0.23749573380464106*softplus(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4])) + -0.2344836196764124*softplus(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4])) + 0.07066758651932892*softplus(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))) + softplus(-0.47638200661148833 + 0.9126631115783219*softplus(-0.8376089108524591 + -0.8472902966255957*$(x[1]) + -0.95763763606758*$(x[2]) + 0.3696542730614176*$(x[3]) + -0.09481436100741814*$(x[4])) + 0.35414400688245484*softplus(0.34708131238661144 + 0.38455737939114654*$(x[1]) + 0.5946120518836455*$(x[2]) + 0.21996792037412316*$(x[3]) + -0.6970126448723701*$(x[4])) + -0.36037763135774403*softplus(-0.38902444039423223 + -0.4012023975891066*$(x[1]) + -0.7732491328671305*$(x[2]) + -0.16777672441934754*$(x[3]) + -0.6363338156602314*$(x[4])) + -0.18657915153735738*softplus(0.16738959928165498 + 0.5119956098239467*$(x[1]) + 0.7821421131636019*$(x[2]) + 0.25850635081101103*$(x[3]) + -0.9052769726690872*$(x[4])) + 0.7872605277677782*softplus(0.1398101826443119 + 0.0438169969151434*$(x[1]) + -0.2051299657215857*$(x[2]) + 0.7229873931159991*$(x[3]) + 0.6169137554983704*$(x[4]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
## UC Model Ref
branches = Dict{Symbol, DeviceModel}(:L => DeviceModel(PSY.Line, PSI.StaticLine),
                                     :T => DeviceModel(PSY.Transformer2W, PSI.StaticTransformer),
                                     :TT => DeviceModel(PSY.TapTransformer, PSI.StaticTransformer),
                                     :dc_line => DeviceModel(PSY.HVDCLine, PSI.HVDCDispatch))

services = Dict{Symbol, PSI.ServiceModel}()

devices = Dict{Symbol, DeviceModel}(:Generators => DeviceModel(PSY.ThermalStandard, PSI.ThermalStandardUnitCommitment),
                                    :Ren => DeviceModel(PSY.RenewableDispatch, PSI.RenewableFullDispatch),
                                    :Loads =>  DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad),
                                    :ILoads =>  DeviceModel(PSY.InterruptibleLoad, PSI.StaticPowerLoad))

template_uc= OperationsTemplate(CopperPlatePowerModel, devices, branches, services);

## ED Model Ref
branches = Dict{Symbol, DeviceModel}(:L => DeviceModel(PSY.Line, PSI.StaticLine),
                                     :T => DeviceModel(PSY.Transformer2W, PSI.StaticTransformer),
                                     :TT => DeviceModel(PSY.TapTransformer, PSI.StaticTransformer),
                                     :dc_line => DeviceModel(PSY.HVDCLine, PSI.HVDCDispatch))

services = Dict{Symbol, PSI.ServiceModel}()

devices = Dict{Symbol, DeviceModel}(:Generators => DeviceModel(PSY.ThermalStandard, PSI.ThermalDispatch, SemiContinuousFF(:P, :ON)),
                                    :Ren => DeviceModel(PSY.RenewableDispatch, PSI.RenewableFullDispatch, RangeFF(:P)),
                                    :Loads =>  DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad),
                                    :ILoads =>  DeviceModel(PSY.InterruptibleLoad, PSI.InterruptiblePowerLoad,))

template_ed= OperationsTemplate(CopperPlatePowerModel, devices, branches, services);


text:
module Dec8Part1
visited_index = Int64[]
function get_new_index_compute_accumulator(op, value, current_index, accumulator)
    if op == "jmp"
        current_index += value
    elseif op == "acc"
        accumulator += value
        current_index += 1
    else
        current_index += 1
    end
    return current_index, accumulator
end

function parse_operation(x::String)::Tuple{String,Int64}
    op, value = split(x, ' ')
    return op, parse(Int64, value)
end

function get_accumulator_run_operations(operations::Array, index::Int64, accumulator::Int64)
    global visited_index
    push!(visited_index, index)
    while true
        op, value = operations[index]
        index, accumulator = get_new_index_compute_accumulator(op, value, index, accumulator)
        if index ∉ visited_index
            push!(visited_index, index)
        else
            return accumulator
        end
    end
end

function main(file_path::String)::Int64
    operations = readlines(file_path)
    current_index = 1
    accumulator = 0
    operations = parse_operation.(operations)
    return get_accumulator_run_operations(operations, current_index, accumulator)
end
end


text:
using Test
include("module-2-exercise-1.jl")

@testset "$(get_module_2_excersise_1_author_name()): odd number tests" begin
     odd_numbers = get_odd_numbers(0)
     @test [] == odd_numbers
     odd_numbers = get_odd_numbers(1)
     @test [1] == odd_numbers
     odd_numbers = get_odd_numbers(5)
     @test [1, 3, 5, 7, 9] == odd_numbers
end


text:
function IRNSSAlmanac(arg0::jint, arg1::jint, arg10::jdouble, arg11::jdouble, arg12::AbsoluteDate, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return IRNSSAlmanac((jint, jint, jdouble, jdouble, AbsoluteDate, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg11, arg12, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function get_af0(obj::IRNSSAlmanac)
    return jcall(obj, "getAf0", jdouble, ())
end

function get_af1(obj::IRNSSAlmanac)
    return jcall(obj, "getAf1", jdouble, ())
end

function get_cic(obj::IRNSSAlmanac)
    return jcall(obj, "getCic", jdouble, ())
end

function get_cis(obj::IRNSSAlmanac)
    return jcall(obj, "getCis", jdouble, ())
end

function get_crc(obj::IRNSSAlmanac)
    return jcall(obj, "getCrc", jdouble, ())
end

function get_crs(obj::IRNSSAlmanac)
    return jcall(obj, "getCrs", jdouble, ())
end

function get_cuc(obj::IRNSSAlmanac)
    return jcall(obj, "getCuc", jdouble, ())
end

function get_cus(obj::IRNSSAlmanac)
    return jcall(obj, "getCus", jdouble, ())
end

function get_date(obj::IRNSSAlmanac)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_e(obj::IRNSSAlmanac)
    return jcall(obj, "getE", jdouble, ())
end

function get_i0(obj::IRNSSAlmanac)
    return jcall(obj, "getI0", jdouble, ())
end

function get_i_dot(obj::IRNSSAlmanac)
    return jcall(obj, "getIDot", jdouble, ())
end

function get_iodec(obj::IRNSSOrbitalElements)
    return jcall(obj, "getIODEC", jint, ())
end

function get_m0(obj::IRNSSAlmanac)
    return jcall(obj, "getM0", jdouble, ())
end

function get_mean_motion(obj::IRNSSAlmanac)
    return jcall(obj, "getMeanMotion", jdouble, ())
end

function get_omega0(obj::IRNSSAlmanac)
    return jcall(obj, "getOmega0", jdouble, ())
end

function get_omega_dot(obj::IRNSSAlmanac)
    return jcall(obj, "getOmegaDot", jdouble, ())
end

function get_pa(obj::IRNSSAlmanac)
    return jcall(obj, "getPa", jdouble, ())
end

function get_prn(obj::IRNSSAlmanac)
    return jcall(obj, "getPRN", jint, ())
end

function get_sma(obj::IRNSSAlmanac)
    return jcall(obj, "getSma", jdouble, ())
end

function get_tgd(obj::IRNSSOrbitalElements)
    return jcall(obj, "getTGD", jdouble, ())
end

function get_time(obj::IRNSSAlmanac)
    return jcall(obj, "getTime", jdouble, ())
end

function get_week(obj::IRNSSAlmanac)
    return jcall(obj, "getWeek", jint, ())
end


text:
type Variance <: ContinuousUnivariateStreamStat
    m::Float64
    sum_sqs::Float64
    v_hat::Float64
    n::Int
end

Variance() = Variance(0.0, 0.0, NaN, 0)

function update!(stat::Variance, x::Real)
    stat.n += 1

    if stat.n == 1
        stat.m = x
        stat.sum_sqs = 0.0
        stat.v_hat = NaN
    else
        m_new = stat.m + (x - stat.m) / stat.n
        stat.sum_sqs += (x - stat.m) * (x - m_new)
        stat.m = m_new
        stat.v_hat = stat.sum_sqs / (stat.n - 1)
    end

    return
end

Base.var(stat::Variance) = stat.v_hat

Base.std(stat::Variance) = sqrt(var(stat))

Base.mean(stat::Variance) = stat.m

state(stat::Variance) = var(stat)

nobs(stat::Variance) = stat.n

function Base.copy(stat::Variance)
    return Variance(stat.m, stat.sum_sqs, stat.v_hat, stat.n)
end

function Base.merge(a::Variance, b::Variance)
    n = a.n + b.n
    m = (a.n / n) * a.m + (b.n / n) * b.m
    sum_sqs = a.sum_sqs + b.sum_sqs
    v_hat = (a.n / n) * a.v_hat + (b.n / n) * b.v_hat
    return Variance(m, sum_sqs, v_hat, n)
end

function Base.empty!(stat::Variance)
    stat.m = 0.0
    stat.sum_sqs = 0.0
    stat.v_hat = NaN
    stat.n = 0
    return
end

function Base.show(io::IO, stat::Variance)
    v = var(stat)
    s = std(stat)
    n = nobs(stat)
    @printf(io, "Online Variance\n")
    @printf(io, " * Variance:  %f\n", v)
    @printf(io, " * Std. Dev.: %f\n", v)
    @printf(io, " * N:         %d\n", n)
    return
end


text:
using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra
using SparseArrays

@testset "MSLP" begin

    target=0;
    nev=3;
    nep=nep_gallery("beam",5);

    # Full matrix test
    λ1,v1=mslp(nep)
    @test norm(compute_Mlincomb(nep,λ1,v1))<1e-10

    Av=[sparse(nep.A[1]),sparse(nep.A[2])];
    nep2 = DEP(nep.n, Av, nep.tauv)
    # Sparse matrix test
    λ2,v2 = mslp(nep2)
    @test norm(compute_Mlincomb(nep2,λ2,v2)) < 1e-10

    @test λ1≈λ2 # They should be exactly the same

    @info "mslp noconv + double"
    @bench @testset  "mslp + double" begin
        nep3 = nep_gallery("dep_double");
        @test_throws NoConvergenceException mslp(nep3, λ=9im, maxit=10)
        λ,v = mslp(nep3, λ=9im, maxit=100,errmeasure=ResidualErrmeasure(nep3))
        @test norm(compute_Mlincomb(nep3, λ, v)) < eps()*100
    end
end


text:
@testset "example files" begin

@testset "bit depths" begin
    for i in 1:32
        b = Bgen(BGEN.datadir("example.$(i)bits.bgen"); idx_path = nothing)
        v = parse_variants(b)
        for (j, (v1, v2)) in enumerate(zip(v, gen_data))
            @test v1.chrom == v2.chrom
            @test v1.varid == v2.varid
            @test v1.rsid == v2.rsid
            @test v1.pos == v2.pos
            @test all(v1.alleles .== v2.alleles)
            @test array_equal(v2.probs, probabilities!(b, v1), i)
        end
    end
end

@testset "zstd" begin
    b = Bgen(BGEN.datadir("example.16bits.zstd.bgen"); idx_path = nothing)
    v = parse_variants(b)
    for (j, (v1, v2)) in enumerate(zip(v, gen_data))
        @test v1.chrom == v2.chrom
        @test v1.varid == v2.varid
        @test v1.rsid == v2.rsid
        @test v1.pos == v2.pos
        @test all(v1.alleles .== v2.alleles)
        @test array_equal(v2.probs, probabilities!(b, v1), 16)
    end
end

@testset "v11" begin
    b = Bgen(BGEN.datadir("example.v11.bgen"); idx_path = nothing)
    v = parse_variants(b)
    for (j, (v1, v2)) in enumerate(zip(v, gen_data))
        @test v1.chrom == v2.chrom
        @test v1.varid == v2.varid
        @test v1.rsid == v2.rsid
        @test v1.pos == v2.pos
        @test all(v1.alleles .== v2.alleles)
        @test array_equal(v2.probs, probabilities!(b, v1), 16)
    end
end

@testset "haplotypes" begin
    b = Bgen(BGEN.datadir("haplotypes.bgen"); idx_path = nothing)
    v = parse_variants(b)
    for (j, (v1, v2)) in enumerate(zip(v, haps_data))
        @test v1.chrom == v2.chrom
        @test v1.varid == v2.varid
        @test v1.rsid == v2.rsid
        @test v1.pos == v2.pos
        @test all(v1.alleles .== v2.alleles)
        @test array_equal(v2.probs, probabilities!(b, v1), 16)
    end
end

@testset "complex" begin
    b = Bgen(BGEN.datadir("complex.bgen"); idx_path = nothing)
    v = parse_variants(b)
    for (j, (v1, v2)) in enumerate(zip(v, vcf_data))
        @test v1.chrom == v2.chrom
        @test v1.varid == v2.varid
        @test v1.rsid == v2.rsid
        @test v1.pos == v2.pos
        @test all(v1.alleles .== v2.alleles)
        @test array_equal(v2.probs, probabilities!(b, v1), 16)
    end
end

@testset "complex bit depths" begin
    for i in 1:32
        b = Bgen(BGEN.datadir("complex.$(i)bits.bgen"); idx_path = nothing)
        v = parse_variants(b)
        for (j, (v1, v2)) in enumerate(zip(v, vcf_data))
            @test v1.chrom == v2.chrom
            @test v1.varid == v2.varid
            @test v1.rsid == v2.rsid
            @test v1.pos == v2.pos
            @test all(v1.alleles .== v2.alleles)
            @test array_equal(v2.probs, probabilities!(b, v1), i)
        end
    end
end

@testset "null" begin
    @test_throws SystemError Bgen(BGEN.datadir("Hello_World.bgen"))
end
end


text:
module Geometry

using StaticArrays, LinearAlgebra, DocStringExtensions
using ..Grids: _ξ1x1, _ξ2x1, _ξ3x1, _ξ1x2, _ξ2x2, _ξ3x2, _ξ1x3, _ξ2x3, _ξ3x3, _M, _MI,
       _x1, _x2, _x3, _JcV
export LocalGeometry, lengthscale, resolutionmetric

"""
    LocalGeometry

The local geometry at a nodal point.

# Constructors

    LocalGeometry(polynomial::Val, vgeo::AbstractArray{T}, n::Integer, e::Integer)

Extracts a `LocalGeometry` object from the `vgeo` array at node `n` in element `e`.

# Fields

$(DocStringExtensions.FIELDS)
"""
struct LocalGeometry{T,P}
  "Polynomial interpolant: currently this is assumed to be `Val{polyorder}`, but this may change in future."
  polynomial::P
  "Cartesian coordinates"
  coord::SVector{3,T}
  "Jacobian from Cartesian to element coordinates: `invJ[i,j]` is ``∂ξ_i/∂x_j``"
  invJ::SMatrix{3,3,T,9}
end

function LocalGeometry(polynomial::Val, vgeo::AbstractArray{T}, n::Integer, e::Integer) where {T}
  coord = @SVector T[vgeo[n, _x1, e], vgeo[n, _x2, e], vgeo[n, _x3, e]]
  invJ = @SMatrix T[vgeo[n, _ξ1x1, e] vgeo[n, _ξ1x2, e] vgeo[n, _ξ1x3, e];
                    vgeo[n, _ξ2x1, e] vgeo[n, _ξ2x2, e] vgeo[n, _ξ2x3, e];
                    vgeo[n, _ξ3x1, e] vgeo[n, _ξ3x2, e] vgeo[n, _ξ3x3, e]]

  LocalGeometry(polynomial, coord, invJ)
end

"""
    resolutionmetric(g::LocalGeometry)

The metric tensor of the discretisation resolution. Given a unit vector `u` in Cartesian
coordinates and `M = resolutionmetric(g)`, `sqrt(u'*M*u)` is the degree-of-freedom density
in the direction of `u`.
"""
function resolutionmetric(g::LocalGeometry{T, Val{polyorder}}) where {T,polyorder}
  S = polyorder*g.invJ/2
  S'*S # TODO: return an eigendecomposition / symmetric object?
end

"""
    lengthscale(g::LocalGeometry)

The effective grid resolution at the point.
"""
function lengthscale(g::LocalGeometry{T, Val{polyorder}}) where {T,polyorder}
  2/(cbrt(det(g.invJ))*polyorder)
end

end # module

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
CONFIG = Dict(
    "MSISDN" => "+xyzxyzxyzxy",
    "USERNAME" => "YOURUSERNAME",
    "PASSWORD" => "YOURPASSWORD"
)


text:
################################################################################
#
#
#
################################################################################

isunit(a::NumFieldElem) = !iszero(a)

canonical_unit(a::NumFieldElem) = a

absolute_degree(A::AnticNumberField) = degree(A)

################################################################################
#
#  Base case for dot products
#
################################################################################

dot(x::fmpz, y::NumFieldElem) = x * y

dot(x::Integer, y::NumFieldElem) = x * y

dot(x::NumFieldElem, y::Integer) = x * y

function dot(a::Array{<: NumFieldElem, 1}, b::Array{fmpz, 1})
  d = zero(parent(a[1]))
  t = zero(d)
  for i=1:length(a)
    mul!(t, a[i], b[i])
    add!(d, d, t)
  end
  return d
end

################################################################################
#
#  Parent
#
################################################################################

@doc Markdown.doc"""
    parent(a::NumFieldElem) -> NumField

Given an element `a` of a number field $K$, this function returns $K$.
"""
parent(a::NumFieldElem)

################################################################################
#
#  Integrality
#
################################################################################

isintegral(a::fmpq) = isone(denominator(a))

@doc Markdown.doc"""
    isintegral(a::NumFieldElem) -> Bool

Returns whether `a` is integral, that is, whether the minimal polynomial of $a$
has integral coefficients.
"""
function isintegral(a::NumFieldElem)
  f = minpoly(a)
  for i = 0:(degree(f) - 1)
    if !isintegral(coeff(f, i))
      return false
    end
  end
  return true
end

################################################################################
#
#  Random elements from arrays of number field elements
#
################################################################################

@doc Markdown.doc"""
    rand(b::Vector{NumFieldElem}, r::UnitRange) -> NumFieldElem

A random linear combination of elements in `b` with coefficients in `r`.
"""
function rand(b::Vector{<: NumFieldElem}, r::UnitRange)
  length(b) == 0 && error("Array must not be empty")
  s = zero(parent(b[1]))
  rand!(s, b, r)
  return s
end

@doc Markdown.doc"""
    rand(b::Vector{NumFieldElem}, r::UnitRange, terms::Int) -> NumFieldElem

A random linear combination (with repetitions) of `terms` elements of `b`
with coefficients in `r`.
"""
function rand(b::Vector{<: NumFieldElem}, r::UnitRange, terms::Int)
  length(b) == 0 && error("Array must not be empty")
  s = zero(parent(b[1]))
  rand!(s, b, r, terms)
  return s
end

@doc Markdown.doc"""
    rand!(c::NumFieldElem, b::Vector{NumFieldElem},
          r::UnitRange, terms::Int) -> NumFieldElem

Sets `c` to a random linear combination (with repetitions) of \code{terms}
elements of `b` with coefficients in `r`. The element `c` is returned.
"""
function rand!(c::T, b::Vector{T}, r::UnitRange,
               terms::Int) where {T <: NumFieldElem}
  length(b) == 0 && error("Array must not be empty")
  (terms <= 0 || terms > length(b)) && error("Number of terms should be at least 1 and cannot exceed the length of the array")

  t = zero(parent(c))

  terms = min(terms, length(b))
  mul!(c, rand(b), rand(r))
  for i = 2:terms
    mul!(t, rand(b), rand(r))
    add!(c, t, c)
  end

  return c
end

@doc Markdown.doc"""
    rand!(c::NumFieldElem, b::Vector{NumFieldElem}, r::UnitRange) -> NumFieldElem

Sets `c` to a random linear combination of elements in `b` with coefficients
in `r`. The element `c` is returned.
"""
function rand!(c::T, b::Vector{T}, r::UnitRange) where {T <: NumFieldElem}
  length(b) == 0 && error("Array must not be empty")

  mul!(c, b[1], rand(r))
  t = zero(parent(c))

  for i = 2:length(b)
    t = mul!(t, b[i], rand(r))
    c = add!(c, t, c)
  end

  return c
end

################################################################################
#
#  Basis matrix
#
################################################################################

@doc Markdown.doc"""
    basis_matrix(v::Vector{NumFieldElem}) -> Mat

Given a vector `v` of `n` elements of a number field `K` of degree `d`, this
function returns an `n x d` matrix with entries in the base field of $K$,
 where row `i` contains the > coefficients of `v[i]` with respect of the
canonical basis of `K`.
"""
basis_matrix(v::Vector{<: NumFieldElem})

################################################################################
#
#  Characteristic polynomial
#
################################################################################

@doc Markdown.doc"""
    charpoly(a::NumFieldElem) -> PolyElem

Given a number field element `a` of a number field `K`, this function returns
the characteristic polynomial of `a` over the base field of `K`.
"""
charpoly(::NumFieldElem)

@doc Markdown.doc"""
    absolute_charpoly(a::NumFieldElem) -> PolyElem

Given a number field element `a` of a number field `K`, this function returns
the characteristic polynomial of `a` over the rationals $\mathbf{Q}$.
"""
absolute_charpoly(::NumFieldElem)

@doc Markdown.doc"""
    minpoly(a::NumFieldElem) -> PolyElem

Given a number field element `a` of a number field `K`, this function returns
the minimal polynomial of `a` over the base field of `K`.
"""
minpoly(::NumFieldElem)

@doc Markdown.doc"""
    absolute_minpoly(a::NumFieldElem) -> PolyElem

Given a number field element `a` of a number field `K`, this function returns
the minimal polynomial of `a` over the rationals $\mathbf{Q}$.
"""
absolute_minpoly(::NumFieldElem)

################################################################################
#
#  Powering with fmpz
#
################################################################################

function ^(x::NumFieldElem, y::fmpz)
  if fits(Int, y)
    return x^Int(y)
  end

  return _power(x, y)
end

# We test once if it fits, otherwise we would have to check for every ^-call
function _power(x::NumFieldElem, y::fmpz)
  res = parent(x)()
  if y < 0
    res = _power(inv(x), -y)
  elseif y == 0
    res = parent(x)(1)
  elseif y == 1
    res = deepcopy(x)
  elseif mod(y, 2) == 0
    z = _power(x, div(y, 2))
    res = z*z
  else
    res = _power(x, y - 1) * x
  end
  return res
end

################################################################################
#
#  SimpleNumFieldElem
#
################################################################################

@doc Markdown.doc"""
    coeff(a::SimpleNumFieldElem, i::Int) -> FieldElem

Given a number field element `a` of a simple number field extension `L/K`, this
function returns the `i`-th coefficient of `a`, when expanded in the canonical
power basis of `L`. The result is an element of `K`.
"""
coeff(::SimpleNumFieldElem, ::Int)

# copy does not do anything (so far), this is only for compatibility with coeffs(::AbsAlgAssElem)

@doc Markdown.doc"""
    coeffs(a::SimpleNumFieldElem, i::Int) -> Vector{FieldElem}

Given a number field element `a` of a simple number field extension `L/K`, this
function returns the coefficients of `a`, when expanded in the canonical
power basis of `L`.
"""
function coeffs(a::SimpleNumFieldElem; copy::Bool = false)
  return [ coeff(a, i) for i = 0:degree(parent(a)) - 1 ]
end

@doc Markdown.doc"""
    (L::SimpleNumField)(c::Vector) -> SimpleNumFieldElem

Given a simple number field extension `L/K` and a vector of elements that are
coercible to `K`, this functions returns the element of `L` with coefficients
`c`.
"""
(K::SimpleNumField)(c::Vector)

function (K::AnticNumberField)(c::Vector{fmpq})
  @assert length(c) == degree(K)
  return K(parent(K.pol)(c))
end

function (L::NfRel{T})(a::Vector{T}) where T
  z = NfRelElem{T}(parent(L.pol)(a))
  z.parent = L
  return z
end

for F in [AnticNumberField, NfRel]
  @eval begin
    function (L::$F)(a::Vector)
      length(a) != degree(L) && throw(error("Vector must have length ($(length(a))) equal to the degree ($(degree(L)))"))
      K = base_field(L)
      aprom = Vector{elem_type(K)}(undef, degree(L))
      for i in 1:degree(L)
        aprom[i] = K(a[i])
      end
      return L(aprom)
    end
  end
end

################################################################################
#
#  Norm and trace
#
################################################################################

@doc Markdown.doc"""
    tr(a::NumFieldElem) -> NumFieldElem

Returns the trace of an element $a$ of a number field extension `L/K`. This
will be an element of `K`.
"""
tr(::NumFieldElem)

@doc Markdown.doc"""
    norm(a::NumFieldElem) -> NumFieldElem

Returns the norm of an element $a$ of a number field extension `L/K`. This
will be an element of `K`.
"""
norm(::NumFieldElem)

function _elem_tr_to(a, k::T) where {T}
  if parent(a) isa T
    @assert parent(a) == k
    return a
  else
    return _elem_tr_to(tr(a), k)
  end
end

function _elem_norm_to(a, k::T) where {T}
  if parent(a) isa T
    @assert parent(a) == k
    return a
  else
    return _elem_norm_to(norm(a), k)
  end
end

@doc Markdown.doc"""
    tr(a::NumFieldElem, k::NumField) -> NumFieldElem

Returns the trace of an element $a$ of a number field `L` with respect to
a subfield $k$ of $L$. This will be an element of `k`.
"""
function tr(a::NumFieldElem, k::NumField)
  _elem_tr_to(a, k)
end

tr(a::NumFieldElem, ::FlintRationalField) = _elem_tr_to(a, FlintQQ)

@doc Markdown.doc"""
    norm(a::NumFieldElem, k::NumField) -> NumFieldElem

Returns the norm of an element $a$ of a number field `L` with respect to
a subfield $k$ of $L$. This will be an element of `k`.
"""
function norm(a::NumFieldElem, k::NumField)
  _elem_norm_to(a, k)
end

norm(a::NumFieldElem, ::FlintRationalField) = _elem_norm_to(a, FlintQQ)

@doc Markdown.doc"""
    absolute_tr(a::NumFieldElem) -> fmpq

Given a number field element $a$, returns the absolute trace of $a$.
"""
function absolute_tr(a::NfRelElem)
  return tr(a, FlintQQ)
end

absolute_tr(a::nf_elem) = tr(a)

absolute_tr(x::fmpq) = x

@doc Markdown.doc"""
    absolute_norm(a::NumFieldElem) -> fmpq

Given a number field element $a$, returns the absolute norm of $a$.
"""
function absolute_norm(a::NfRelElem)
  return norm(a, FlintQQ)
end

absolute_norm(a::nf_elem) = norm(a)

absolute_norm(a::fmpq) = x

################################################################################
#
#  Polynomials
#
################################################################################

@doc Markdown.doc"""
    norm(f::PolyElem{<:NumFieldElem}) -> PolyElem

Returns the norm of $f$, that is, the product of all conjugates of $f$ taken
coefficientwise.
"""
function norm(f::PolyElem{<: NumFieldElem})
  K = base_ring(f)
  P = polynomial_to_power_sums(f, degree(f)*degree(K))
  PQ = elem_type(base_field(K))[tr(x) for x in P]
  return power_sums_to_polynomial(PQ)
end

function isirreducible(f::PolyElem{<: NumFieldElem})
  # TODO (easy): We can do better then this. First do a squarefree factorization
  lf = factor(f)
  return sum(values(lf.fac)) == 1
end

function factor(f::PolyElem{<: NumFieldElem})
  K = base_ring(f)
  Ka, rel_abs, _ = absolute_field(K)

  function map_poly(P::Ring, mp::Map, f::Generic.Poly)
    return P([mp(coeff(f, i)) for i=0:degree(f)])
  end

  fa = map_poly(PolynomialRing(Ka, "T", cached=false)[1], pseudo_inv(rel_abs), f)
  lf = factor(fa)
  res = Fac(map_poly(parent(f), rel_abs, lf.unit), Dict(map_poly(parent(f), rel_abs, k)=>v for (k,v) = lf.fac))

  return res
end

function roots(f::PolyElem{<: NumFieldElem})
  lf = factor(f)
  @assert degree(unit(lf)) == 0
  scale = inv(coeff(unit(lf), 0))
  return elem_type(base_ring(f))[-constant_coefficient(x)*scale for x = keys(lf.fac) if degree(x) == 1]
end

################################################################################
#
#
#
################################################################################

@doc Markdown.doc"""
    representation_matrix(a::NumFieldElem) -> MatElem

Returns the representation matrix of $a$, that is, the matrix representing
multiplication with $a$ with respect to the canonical basis of the parent of $a$.
"""
representation_matrix(a::NumFieldElem)

################################################################################
#
#
#
################################################################################

@doc Markdown.doc"""
    gen(L::SimpleNumField) -> NumFieldElem

Given a simple number field $L = K[x]/(x)$ over $K$, this functions returns the
class of $x$, which is the canonical primitive element of $L$ over $K$.
"""
gen(::SimpleNumField)

@doc Markdown.doc"""
    gens(L::NonSimpleNumField) -> Vector{NumFieldElem}

Given a non-simple number field $L = K[x_1,\dotsc,x_n]/(f_1,\dotsc,f_n)$ over
$K$, this functions returns the list $\bar x_1,\dotsc,\bar x_n$.
"""
gens(::NonSimpleNumField)


text:
using DataFrames, RData, LinearAlgebra, Statistics

years = 2000:2014 # vector with years covered by WIOD Rev. 2016
N = 41 # number of countries
S = 35 # number of industries
dir = "C:/Users/u0148308/Desktop/raw/" # location of raw data

# --------------- Import data ---------------------------------------------------------------------------------------------------------------------------------

function import_excel(dir::String, year::Integer, N::Integer, S::Integer)

    # load data for specified year and transform into a DataFrame (automatically)
    path = dir * "WIOT" * string(year)[end-1:end] * "_ROW_Apr12.xlsx"
    df = DataFrame(XLSX.readxlsx(path)["WIOT_$year"][:], :auto)
    df = df[7:end,5:end]

    # Intermediate demand
    Z = df[1:N*S, 1:N*S] # NS×NS
    Z = Matrix(convert.(Float64, Z))

    # Final demand
    F = df[1:N*S, N*S+1:end-1] # NS×5*N

    # Subset final demand for inventory adjustments
    inventory_columns = 5:5:5*N
    IV = F[:, inventory_columns] # NS×N
    F = F[:, Not(inventory_columns)] # NS×4*N

    IV = Matrix(convert.(Float64, IV)) # NS×N
    F = Matrix(convert.(Float64, F)) # NS×4*N

    return Z, F, IV
end

function antras_trans1(Z::Matrix, F::Matrix, IV::Matrix, N::Integer, S::Integer)

    XX = Z
    FF = F
    NN = IV

    NN = [sum(NN[i,:]) for i in 1:N*S] # NS×1

    YY = [sum(XX[i,:]) + sum(FF[i,:]) + NN[i] for i in 1:N*S] # NS×1
    YY = ifelse.(abs.(YY) .< 1e-6, 0.0, YY)

    adj = YY.-NN
    adj = ifelse.(adj .< 0.0, 0.0, adj)

    XXi = XX .* (repeat(YY,1,N*S) ./ repeat(adj,1,N*S)) # NS×NS, inventory adjustment
    XXi = ifelse.(isnan.(XXi), 0.0, XXi)
    XXi = ifelse.(isinf.(XXi), 0.0, XXi)

    FFi = FF .* (repeat(YY,1,N*4) ./ repeat(adj,1,N*4)) # NS×N*4, inventory adjustment
    FFi = ifelse.(isnan.(FFi), 0.0, FFi)
    FFi = ifelse.(isinf.(FFi), 0.0, FFi)
    FFdi = [sum(FFi[i,j:j+3]) for i in 1:N*S, j in 1:4:N*4] # NS×N

    YYi = [sum(XXi[i,:]) + sum(FFdi[i,:]) for i in 1:N*S] # NS×1

    VVi = [YYi[i] - sum(XXi[:,i]) for i in 1:N*S] # NS×1

    return XXi, FFdi, YYi, VVi
end

function antras_trans2(XXi::Matrix, FFdi::Matrix, YYi::Vector, VVi::Vector)

    gammas = zeros(S+1,N*S)

    for n in 1:N
    for s in 1:S
        inputlist = XXi[:,(n-1)*S+s]
        inputlist = reshape(inputlist, S, N)
        inputlist = [sum(inputlist[i,:]) for i in 1:S]
        inputlist = [inputlist; VVi[(n-1)*S+s]]
        inputratio = inputlist ./ sum(inputlist)
        gammas[:,(n-1)*S+s] = inputratio
    end
    end

    VAcoeff = gammas[S+1,:]
    VAcoeff = ifelse.(isnan.(VAcoeff), 1.0, VAcoeff)
    VAcoeff = reshape(VAcoeff, S, N)

    gammas = gammas[1:S, 1:N*S]
    gammas = ifelse.(isnan.(gammas), 0.0, gammas)

    alphas = zeros(S, N)
    CtyCons = [sum(FFdi[:,i]) for i in 1:N]

    for j in 1:N
        for i in 1:S
            alphas[i, j] = sum(FFdi[i:S:(N-1)*S+i, j]) / CtyCons[j]
        end
    end

    VAn = [sum(VVi[i:i+S-1]) for i in 1:S:N*S] # N×1

    Piagg = [sum(XXi[s:S:(N-1)*S+s,j]) for s in 1:S, j in 1:N*S]
    PiFagg = [sum(FFdi[s:S:(N-1)*S+s,j]) for s in 1:S, j in 1:N]

    Pi = XXi ./ repeat(Piagg, N)
    Pi = ifelse.(isnan.(Pi), 0.0, Pi)

    PiF = FFdi ./ repeat(PiFagg, N)
    PiF = ifelse.(isnan.(PiF), 0.0, PiF)

    M = zeros(N)
    E = zeros(N)

    for j in 1:N
    for s in 1:S
        for i in 1:N
            iiss = (i-1)*S+s
            jjss = (j-1)*S+s
            for r in 1:S
                jjrr = (j-1)*S+r
                iirr = (i-1)*S+r
                M[j] = M[j] + Pi[iiss,jjrr]*gammas[s,jjrr]*YYi[jjrr]
                E[j] = E[j] + Pi[jjss,iirr]*gammas[s,iirr]*YYi[iirr]
            end
            M[j] = M[j] + PiF[iiss,j]*alphas[s,j]*CtyCons[j]
            E[j] = E[j] + PiF[jjss,i]*alphas[s,i]*CtyCons[i]
        end
    end
    end

    Surplusn = E .- M

    return gammas, alphas, Surplusn, VAcoeff, VAn, Pi, PiF, CtyCons
end

Z, F, IV = import_excel(dir, 1995, N, S)
##################
XXi, FFdi, YYi, VVi = antras_trans1(Z, F, IV, N, S)
gammas, alphas, Surplusn, VAcoeff, VAn, Pi, PiF, CtyCons = antras_trans2(XXi, FFdi, YYi, VVi)

Z, F, Y, VA = create_matrices(Z, F, IV, N, S)
π_Z, π_F = create_trade_shares(Z, F, N, S)
γ, α, VA_coeff, TB_ctry, VA_ctry, F_ctry = create_expenditure_shares(Z, F, Y, VA, π_Z, π_F, N, S)
#####################

alphas == α
gammas == γ
PiF == π_F
Pi == π_Z
F == FFdi
F_ctry == CtyCons

VA_ctry.- F_ctry
Surplusn
TB_ctry
Surplusn .- TB_ctry

[sum(α[:,j]) for j in 1:N]
[sum(alphas[:,j]) for j in 1:N]

count(round.(gammas,digits=0) .!= round.(γ,digits=0))

count(abs.(alphas .- α) .> 0.0001)

γ2 = zeros(S+1,N*S)
for n in 1:N
    for s in 1:S
        inputlist = Z[:,(n-1)*S+s]
        inputlist = reshape(inputlist, S, N)
        inputlist = [sum(inputlist[i,:]) for i in 1:S]
        inputlist = [inputlist; VA[(n-1)*S+s]]
        inputratio = inputlist ./ sum(inputlist)
        γ2[:,(n-1)*S+s] = inputratio
    end
end

VA2 = γ2[S+1,:]
VA2 = ifelse.(isnan.(VA2), 1.0, VA2)
VA2 = reshape(VA2, S, N)

γ2 = γ2[1:S, 1:N*S]
γ2 = ifelse.(isnan.(γ2), 0.0, γ2)

count(abs.(gammas .- γ2) .> 5)

count(γ.>1)

findall(γ.>=1)
findall(γ.<0)

findall(VA_coeff.>1)

findall(π_Z.<0)
findall(π_F.<0)

findall(Z.<0)
findall(VA_ctry.<0)

E_A = zeros(N)
M_A = zeros(N)
for j in 1:N
    for s in 1:S
        for i in 1:N
            iiss = (i-1)*S+s
            jjss = (j-1)*S+s
            for r in 1:S
                jjrr = (j-1)*S+r
                iirr = (i-1)*S+r
                M_A[j] += π_Z[iiss,jjrr]*γ[s,jjrr]*Y[jjrr]
                E_A[j] += π_Z[jjss,iirr]*γ[s,iirr]*Y[iirr]
            end
            M_A[j] += π_F[iiss,j]*α[s,j]*F_ctry[j]
            E_A[j] += π_F[jjss,i]*α[s,i]*F_ctry[i]
        end
    end
end

TB_ctry2 = E_A .- M_A # N×1

F_prime_ctry = copy(F_ctry)
Y_prime = copy(Y_prime)

LHS = zeros(N)
RHS = zeros(N)
for j in 1:N
    for s in 1:S
        for i in 1:N
            iiss = (i-1)*S+s
            jjss = (j-1)*S+s
            for r in 1:S
                jjrr = (j-1)*S+r
                iirr = (i-1)*S+r
                LHS[j] += π_Z[iiss,jjrr]*γ[s,jjrr]*Y[jjrr]
                RHS[j] += π_Z[jjss,iirr]*γ[s,iirr]*Y[iirr]
            end
            LHS[j] += π_F[iiss,j]*α[s,j]*F_ctry[j]
            RHS[j] += π_F[jjss,i]*α[s,i]*F_ctry[i]
        end
    end
end

ETB_ctry = RHS .- LHS

findall(abs.(b).>1)

Z_prime = π_Z .* repeat(γ, N) .* repeat(transpose(Y), N*S) # NS×NS
F_prime = π_F .* repeat(α, N) .* repeat(transpose(F_ctry), N*S) # NS×N

E_prime = [sum(Y[i:i+S-1]) - sum(Z_prime[i:i+S-1,i:i+S-1]) - sum(F_prime[i:i+S-1,ceil(Int, i/S)]) for i in 1:S:N*S] # N×1

M_prime_Z = [sum(Z_prime[:,j:j+S-1]) - sum(Z_prime[j:j+S-1,j:j+S-1]) for j in 1:S:N*S]
M_prime_F = [sum(F_prime[:,ceil(Int, j/S)]) - sum(F_prime[j:j+S-1,ceil(Int, j/S)]) for j in 1:S:N*S]
M_prime = M_prime_Z .+ M_prime_F # N×1

ETB_ctry_prime = E_prime .- M_prime

findall(abs.(Z.-Z_prime).>0.1)
findall(abs.(ETB_ctry.-ETB_ctry_prime).>1)

TB_ctry

ETB = ETB_ctry .- TB_ctry

b = ETB./VA_ctry

δ = sign.(b) .* abs.(vfactor.*b)
w_hat = ones(N)
w_hat = w_hat .* (1.0 .+ δ ./ w_hat)

a = TB_ctry./VA_ctry
δ = sign.(a) .* abs.(vfactor.*a)
w_hat = ones(N)
w_hat = w_hat .* (1.0 .+ δ ./ w_hat)

findall(abs.(a).>1)

findall(TB_ctry .> VA_ctry)

max_stars_repo_path: old stuff/transform_WIOD.jl

max_issues_repo_path: old stuff/transform_WIOD.jl

max_forks_repo_path: old stuff/transform_WIOD.jl


text:
using RickTracy

@resetallsnaps
@watch i j k dsum dprod
dsum = 0
for i in 0:2
    for j in 0:9
        for k in 0:9
            dprod = i*j*k
    #         @show i j dsum dprod
            @snapall loc=inner
            dsum += 1
        end
    end
end
#---

tis = RickTracy.traces4symvals(Dict("i"=>2, "j"=>9, "k"=>9))
println("unique lcounts: ", unique(map(ti->ti.lcount, tis)))
println("length(tis): ", length(tis), "\n---------")
foreach(tis[end-4:end]) do ti
    println(ti.exprstr, ": ", ti.val, " (lcount: $(ti.lcount))")
end

#---

#---
push!(Interact.signal(slidees[1]), 10)

#---
#=
Let's say you had a few functions that:

1) do a little setup - initialise vars etc
2) iterate through a large array
3) do some cleanup and return some values

Now, you want to make your program as fast as possible so (I assume) you only want iterate through that array once, since (I assume) it'll be more cache/register friendly that way.

a) do you think my assumptions are correct
b) is there a good way to structure your code to do this?

What do you think about this:
=#

function dosometask1(bigarray, params...)
    #setup stuff1
    #setup stuff2
    function loop_part(i, array_el)
        #set some vars
        #check some stuff
        return all_done #return true if you want to break
    end
    function resultypoo()
        #finish up
        #return good stuff:
        (good, stuff, etc...)
    end

    loop_part, resultypoo
end

#---
using RickTracy
@resetallsnaps
@watch x y z dsum
dsum = 0
# Nxs = rand(1:2)
# Nys = [rand(1:4) for nx in 1:Nxs]
# Nzs = [rand(1:4) for nx in 1:Nxs for ny in 1:Nys[nx]]
Nxs = 3
Nys = [4,3,4]
Nzs = [5,4,5,3,7,4,4,10,8,9,3]
@show Nxs Nys Nzs

i = 0
for x in 1:Nxs
    # @show "----------------"
    for y in 1:Nys[x]
        i += 1
        # @show "*** $i ***"
        for z in 1:Nzs[i]
            # @show "^^^" x y z
            @snapall loc=inner
            dsum += 1
        end
    end
end
# RickTracy.happysnaps

sub2ind_nonsquare(row_lens, k) = begin
    sum = 0
    if length(k) > 1
        for i in 1:(k[1]-1)
            sum += row_lens[i]
        end
    end
    sum += k[end]
end

ies = RickTracy.get_iter_extrema([:x, :y, :z])
@show "------ end ------"
# @show Nxs Nys Nzs
maxNs = [[Nxs], Nys, Nzs]
for i in 1:3
    for (ien,(k,v)) in enumerate(ies[i])
        j = i==1 ? 1 : sub2ind_nonsquare(maxNs[i-1], k)
        maxv = maxNs[i][j]
        if v != (1,maxv)
            @show "----" i k v (1,maxv) ien
        end
    end
end

#---

#---
[(k,ies[3][k]) for k in sort([keys(ies[3])...;])]
#---

NTuple{2, String}(("1","2")) |> typeof
Tuple{Vararg{String}}(("1","2", "#")) |> typeof
NTuple{1, String}() |> typeof
#---
(1:0...)
#---
N = 3
iter_extrema = [Dict{NTuple{i-1, Int}, Tuple{Int,Int}}() for i in 1:N]
i = 1
tplkey = ()
extreme_default = (typemax(Int), typemin(Int)) #min, max worst case scenarios
iex = get!(iter_extrema[i], tplkey, extreme_default)
#---
isimmutable(iex)

#---

max_stars_repo_name: JobJob/RickTracy


text:
math_namespace_modify = Namespace(
    "sin" => AetherlangObject((line::Ref{UInt}, ns::Namespace, arg::AetherlangObject)->AetherlangObject(sin(arg.dataref)), Ref(AETH_BUILTIN_TYPES.callable)),
    "cos" => AetherlangObject((line::Ref{UInt}, ns::Namespace, arg::AetherlangObject)->AetherlangObject(cos(arg.dataref)), Ref(AETH_BUILTIN_TYPES.callable)),
    "tan" => AetherlangObject((line::Ref{UInt}, ns::Namespace, arg::AetherlangObject)->AetherlangObject(tan(arg.dataref)), Ref(AETH_BUILTIN_TYPES.callable)),
    "cot" => AetherlangObject((line::Ref{UInt}, ns::Namespace, arg::AetherlangObject)->AetherlangObject(cot(arg.dataref)), Ref(AETH_BUILTIN_TYPES.callable)),
    "mean" => AetherlangObject((line::Ref{UInt}, ns::Namespace, args...)->AetherlangObject(sum([arg.dataref for arg in args])/length(args)), Ref(AETH_BUILTIN_TYPES.callable)),
    "pow" => AetherlangObject((line::Ref{UInt}, ns::Namespace, arg1::AetherlangObject, arg2::AetherlangObject)->AetherlangObject(arg1.dataref^arg2.dataref), Ref(AETH_BUILTIN_TYPES.callable)),
    "sqrt" => AetherlangObject((line::Ref{UInt}, ns::Namespace, arg1::AetherlangObject)->AetherlangObject(sqrt(arg.dataref)), Ref(AETH_BUILTIN_TYPES.callable)),
    "div" => AetherlangObject((line::Ref{UInt}, ns::Namespace, arg1::AetherlangObject, arg2::AetherlangObject)->AetherlangObject(arg1.dataref//arg2.dataref), Ref(AETH_BUILTIN_TYPES.callable)),
    "mod" => AetherlangObject((line::Ref{UInt}, ns::Namespace, arg1::AetherlangObject, arg2::AetherlangObject)->AetherlangObject(arg1.dataref%arg2.dataref), Ref(AETH_BUILTIN_TYPES.callable)),
    "pi" => AetherlangObject(pi),
    "eul" => AetherlangObject(exp(1)),
    "exp" => AetherlangObject((line::Ref{UInt}, ns::Namespace, arg::AetherlangObject)->AetherlangObject(exp(arg.dataref)), Ref(AETH_BUILTIN_TYPES.callable)),
    "ln" => AetherlangObject((line::Ref{UInt}, ns::Namespace, arg::AetherlangObject)->AetherlangObject(log(arg.dataref)), Ref(AETH_BUILTIN_TYPES.callable)),
    "log" => AetherlangObject((line::Ref{UInt}, ns::Namespace, arg1::AetherlangObject, arg2::AetherlangObject)->AetherlangObject(log(arg1.dataref, arg2.dataref)), Ref(AETH_BUILTIN_TYPES.callable)),
    "lg" => AetherlangObject((line::Ref{UInt}, ns::Namespace, arg::AetherlangObject)->AetherlangObject(log10(arg.dataref)), Ref(AETH_BUILTIN_TYPES.callable)),
    "min" => AetherlangObject((line::Ref{UInt}, ns::Namespace, args...)->AetherlangObject(min([arg.dataref for arg in args]...)), Ref(AETH_BUILTIN_TYPES.callable)),
    "max" => AetherlangObject((line::Ref{UInt}, ns::Namespace, args...)->AetherlangObject(max([arg.dataref for arg in args]...)), Ref(AETH_BUILTIN_TYPES.callable))
)
math_namespace_modify["%"] = math_namespace_modify["mod"]

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
"""
The `MultiZAF` structure holds the information necessary to perform matrix correction on a collection of
characteristic X-rays that were measured simultaneously from the same element.
Use `zafcorrection(...)` to construct these rather than the internal constructor.
"""
struct MultiZAF
    xrays::Vector{CharXRay}
    zafs::Dict{AtomicSubShell,ZAFCorrection}
    function MultiZAF(xrays, zafs)
        elm = element(xrays[1])
        mat = material(first(values(zafs)))
        e0 = beamEnergy(first(values(zafs)))
        @assert all(f -> isequal(element(f), elm), xrays)
        "MultiZAF constructor: All the characteristic X-rays must be from the same element."
        @assert all(f -> isequal(element(f), elm), keys(zafs))
        "MultiZAF constructor: All the shells must be from the same element as the X-rays."
        @assert all(f -> haskey(zafs, inner(f)), xrays)
        "MultiZAF constructor: There must be a ZAF correction for each characteristic X-ray.",
        @assert all(f -> isequal(material(f), mat), values(zafs))
        "MultiZAF constructor: All the materials must match."
        @assert all(f -> isequal(beamEnergy(f), e0), values(zafs))
        "MultiZAF constructor: All the beam energies must match."
        return new(xrays, zafs)
    end
end

"""
    zafcorrection(
      mctype::Type{<:MatrixCorrection},
      fctype::Type{<:FluorescenceCorrection},
      cctype::Type{<:CoatingCorrection},
      mat::Material,
      cxrs,
      e0,
      coating=missing
    )

Constructs a MultiZAF around the mctype and fctype algorithms for a collection of CharXRay `cxrs`.
"""
function zafcorrection(
    mctype::Type{<:MatrixCorrection},
    fctype::Type{<:FluorescenceCorrection},
    cctype::Type{<:CoatingCorrection},
    mat::Material,
    cxrs,
    e0::Real,
    coating::Union{Film,AbstractVector{Film},Missing} = missing,
)
    mat = asnormalized(mat)
    shells = union(filter(sh->energy(sh)<e0, inner.(cxrs)))
    zafs = Dict((sh, zafcorrection(mctype, fctype, cctype, mat, sh, e0, coating)) for sh in shells)
    return MultiZAF(cxrs, zafs)
end

"""
    zafcorrection(
      mctype::Type{<:MatrixCorrection},
      fctype::Type{<:FluorescenceCorrection},
      cctype::Type{<:CoatingCorrection},
      unk::Material,
      std::Material,
      cxrs,
      e0;
      unkCoating::Union{Film,AbstractVector{Film},Missing} = missing,
      stdCoating::Union{Film,AbstractVector{Film},Missing} = missing,
    )

Constructs a tuple of MultiZAF around the mctype and fctype correction algorithms for the unknown and standard for a
collection of CharXRay `cxrs`.
"""
function zafcorrection(
    mctype::Type{<:MatrixCorrection},
    fctype::Type{<:FluorescenceCorrection},
    cctype::Type{<:CoatingCorrection},
    unk::Material,
    std::Material,
    cxrs,
    e0::Real;
    unkCoating::Union{Film,AbstractVector{Film},Missing} = missing,
    stdCoating::Union{Film,AbstractVector{Film},Missing} = missing,
)
    xrays = filter(cxr->energy(inner(cxr))<e0, cxrs)
    return (
        zafcorrection(mctype, fctype, cctype, unk, xrays, e0, unkCoating),
        zafcorrection(mctype, fctype, cctype, std, xrays, e0, stdCoating),
    )
end

"""
    shells(mz::MultiZAF)

A set of all sub-shells supported by this MultiZAF
"""
NeXLCore.atomicsubshells(mz::MultiZAF) = keys(mz.zafs)
"""
    NeXLCore.element(mz::MultiZAF)

The element associate with this `MultiZAF`
"""
NeXLCore.element(mz::MultiZAF) = element(mz.xrays[1])

"""
    NeXLCore.characteristic(mz::MultiZAF)

The X-rays associated with this `MultiZAF`.
"""
NeXLCore.characteristic(mz::MultiZAF) = mz.xrays

NeXLCore.material(mz::MultiZAF) = material(first(values(mz.zafs)))

beamEnergy(mz::MultiZAF) = beamEnergy(first(values(mz.zafs)))

NeXLCore.name(mz::MultiZAF) = repr(brightest(mz.xrays)) * "+" * string(length(mz.xrays) - 1) * " others"

commonXrays(unk::MultiZAF, std::MultiZAF) = union(characteristic(unk), characteristic(std))

"""
    Z(unk::MultiZAF, std::MultiZAF)

The Z (atomic number) correction for `unk` relative to `std`.
"""
function Z(unk::MultiZAF, std::MultiZAF)
    z, n = 0.0, 0.0
    for (sh, cxrs2) in splitbyshell(commonXrays(unk, std))
        norm = sum(weight.(NormalizeByShell, cxrs2))
        n += norm
        z += Z(unk.zafs[sh], std.zafs[sh]) * norm
    end
    return z / n
end

"""
    A(unk::MultiZAF, std::MultiZAF, θunk::AbstractFloat, θstd::AbstractFloat)

The A (absorption) correction for `unk` relative to `std`.
"""
function A(unk::MultiZAF, std::MultiZAF, θunk::AbstractFloat, θstd::AbstractFloat)
    a, n = 0.0, 0.0
    for (sh, cxrs2) in splitbyshell(commonXrays(unk, std))
        zafU, zafS = unk.zafs[sh], std.zafs[sh]
        for cxr in cxrs2
            w = weight(NormalizeByShell, cxr)
            a += w * A(zafU, zafS, cxr, θunk, θstd)
            n += w
        end
    end
    return a / n
end

"""
    F(unk::MultiZAF, std::MultiZAF, θunk::AbstractFloat, θstd::AbstractFloat)

The F (fluoresence) correction for `unk` relative to `std`.
"""
function F(unk::MultiZAF, std::MultiZAF, θunk::AbstractFloat, θstd::AbstractFloat)
    f, n = 0.0, 0.0
    for (sh, cxrs2) in splitbyshell(commonXrays(unk, std))
        zafU, zafS = unk.zafs[sh], std.zafs[sh]
        for cxr in cxrs2
            w = weight(NormalizeByShell, cxr)
            f += w * F(zafU, zafS, cxr, θunk, θstd)
            n += w
        end
    end
    return f / n
end

"""
    generation(unk::MultiZAF, std::MultiZAF)

The generation correction for `unk` relative to `std`.  Usually, 1.0 unless the standard and unknown were collected
at different beam energies.
"""
function generation(unk::MultiZAF, std::MultiZAF)
    g, n = 0.0, 0.0
    for (sh, cxrs2) in splitbyshell(commonXrays(unk, std))
        zafU, zafS = unk.zafs[sh], std.zafs[sh]
        icxU = ionizationcrosssection(sh, beamEnergy(zafU))
        icxS = ionizationcrosssection(sh, beamEnergy(zafS))
        for cxr in cxrs2
            w = weight(NormalizeByShell, cxr)
            g += w * (icxU / icxS)
            n += w
        end
    end
    return g / n
end

"""
    coating(unk::MultiZAF, std::MultiZAF, θunk::AbstractFloat, θstd::AbstractFloat)

The conductive (or other) coating correction factor.
"""
function coating(unk::MultiZAF, std::MultiZAF, θunk::AbstractFloat, θstd::AbstractFloat)
    c, n = 0.0, 0.0
    for (sh, cxrs2) in splitbyshell(commonXrays(unk, std))
        zafU, zafS = unk.zafs[sh], std.zafs[sh]
        for cxr in cxrs2
            c += weight(NormalizeByShell, cxr) * coating(zafU, zafS, cxr, θunk, θstd)
            n += weight(NormalizeByShell, cxr)
        end
    end
    return c / n
end
"""
    gZAFc(unk::MultiZAF, std::MultiZAF, θunk::AbstractFloat, θstd::AbstractFloat)
    gZAFc(kr::KRatio, unkComp::Material; mc::Type{<:MatrixCorrection} = XPP, fc::Type{<:FluorescenceCorrection} = ReedFluorescence, cc::Type{<:CoatingCorrection} = Coating)

The combined generation, atomic number, absorption and generation corrections.
"""
function gZAFc(unk::MultiZAF, std::MultiZAF, θunk::AbstractFloat, θstd::AbstractFloat)
    a, n = 0.0, 0.0
    for (sh, cxrs2) in splitbyshell(commonXrays(unk, std))
        zafU, zafS = unk.zafs[sh], std.zafs[sh]
        icxU = ionizationcrosssection(sh, beamEnergy(zafU))
        icxS = ionizationcrosssection(sh, beamEnergy(zafS))
        for cxr in cxrs2
            w = weight(NormalizeByShell, cxr)
            a +=
                w *
                (icxU / icxS) *
                ZA(zafU, zafS, cxr, θunk, θstd) *
                F(zafU, zafS, cxr, θunk, θstd) *
                coating(zafU, zafS, cxr, θunk, θstd)
            n += w
        end
    end
    return a / n
end

function gZAFc(
    kr::KRatio,
    unkComp::Material;
    mc::Type{<:MatrixCorrection} = XPP,
    fc::Type{<:FluorescenceCorrection} = ReedFluorescence,
    cc::Type{<:CoatingCorrection} = Coating,
)
    elm = kr.element
    xrays = filter(cxr->energy(inner(cxr))<min(kr.unkProps[:BeamEnergy],kr.stdProps[:BeamEnergy]), kr.xrays)
    zu = zafcorrection(mc, fc, cc, unkComp, xrays, kr.unkProps[:BeamEnergy])
    zs = zafcorrection(mc, fc, cc, kr.standard, xrays, kr.stdProps[:BeamEnergy])
    return gZAFc(zu, zs, kr.unkProps[:TakeOffAngle], kr.stdProps[:TakeOffAngle])
end

"""
    k(unk::MultiZAF, std::MultiZAF, θunk::AbstractFloat, θstd::AbstractFloat)

The computed k-ratio for the unknown relative to standard.
"""
function k(unk::MultiZAF, std::MultiZAF, θunk::AbstractFloat, θstd::AbstractFloat)
    elm = element(unk)
    return (nonneg(material(unk), elm) / nonneg(material(std), elm)) * gZAFc(unk, std, θunk, θstd)
end

"""
    asa(::Type{DataFrame}, unk::MultiZAF, std::MultiZAF, θunk::AbstractFloat, θstd::AbstractFloat)::DataFrame

Tabulate a matrix correction relative to the specified unknown and standard in a DataFrame.
"""
function NeXLUncertainties.asa(
    ::Type{DataFrame},
    unk::MultiZAF,
    std::MultiZAF,
    θunk::AbstractFloat,
    θstd::AbstractFloat,
)::DataFrame
    tot = gZAFc(unk, std, θunk, θstd)
    @assert isequal(element(unk), element(std)) "The unknown and standard's elements must match."
    return DataFrame(
        Unknown = [name(material(unk))],
        E₀ᵤ = [beamEnergy(unk)],
        Standard = [name(material(std))],
        E₀ₛ = [beamEnergy(std)],
        Xrays = [name(union(characteristic(unk), characteristic(std)))],
        Generation = [generation(unk, std)],
        Z = [Z(unk, std)],
        A = [A(unk, std, θunk, θstd)],
        F = [F(unk, std, θunk, θstd)],
        coating = [coating(unk, std, θunk, θstd)],
        gZAFc = [tot],
        k = [tot * material(unk)[element(unk)] / material(std)[element(unk)]],
    )
end

"""
    detail(::Type{DataFrame}, unk::MultiZAF, std::MultiZAF)::DataFrame

Tabulate each term in the MultiZAF matrix correction in a DataFrame.
"""
function detail(::Type{DataFrame}, unk::MultiZAF, std::MultiZAF, θunk::AbstractFloat, θstd::AbstractFloat)::DataFrame
    stds, stdE0, unks = Vector{String}(), Vector{Float64}(), Vector{String}()
    unkE0, xray, g = Vector{Float64}(), Vector{CharXRay}(), Vector{Float64}()
    z, a, f = Vector{Float64}(), Vector{Float64}(), Vector{Float64}()
    c, wgt, zaf, k = Vector{Float64}(), Vector{Float64}(), Vector{Float64}(), Vector{Float64}()
    for (sh, cxrs2) in splitbyshell(commonXrays(unk, std))
        zafU, zafS = unk.zafs[sh], std.zafs[sh]
        matU, matS = material(zafU), material(zafS)
        elm = element(sh)
        for cxr in cxrs2
            push!(unks, name(material(zafU)))
            push!(unkE0, beamEnergy(zafU))
            push!(stds, name(material(zafS)))
            push!(stdE0, beamEnergy(zafS))
            push!(xray, cxr)
            push!(wgt, weight(NormalizeByShell, cxr))
            push!(g, generation(zafU, zafS, inner(cxr)))
            push!(z, Z(zafU, zafS))
            push!(a, A(zafU, zafS, cxr, θunk, θstd))
            push!(f, F(zafU, zafS, cxr, θunk, θstd))
            push!(c, coating(zafU, zafS, cxr, θunk, θstd))
            tot = ZAFc(zafU, zafS, cxr, θunk, θstd)
            push!(zaf, tot)
            push!(k, tot * matU[elm] / matS[elm])
        end
    end
    return DataFrame(
        Unknown = unks,
        E0unk = 0.001 * unkE0,
        Standard = stds,
        E0std = 0.001 * stdE0,
        Xray = xray,
        Weight = wgt,
        Generation = g,
        Z = z,
        A = a,
        F = f,
        c = c,
        ZAF = zaf,
        k = k,
    )
end

"""
    detail(::Type{DataFrame}, mzs::AbstractArray{Tuple{MultiZAF, MultiZAF}})::DataFrame

Tabulate the details of a matrix correction relative to the specified unknown and standard in a DataFrame.
"""
detail(::Type{DataFrame}, mzs::AbstractArray{Tuple{MultiZAF,MultiZAF}}, θunk::AbstractFloat, θstd::AbstractFloat) =
    mapreduce(tmm -> detail(tmm[1], tmm[2], θunk, θstd), append!, mzs)

"""
    asa(::Type{DataFrame}, mzs::AbstractArray{Tuple{MultiZAF,MultiZAF}}, θunk::AbstractFloat, θstd::AbstractFloat)::DataFrame

Tabulate a matrix correction relative to a specified Dict of unknowns and standards in a DataFrame.
"""
NeXLUncertainties.asa(
    ::Type{DataFrame},
    mzs::AbstractArray{Tuple{MultiZAF,MultiZAF}},
    θunk::AbstractFloat,
    θstd::AbstractFloat,
) = mapreduce(tmm -> asa(DataFrame, tmm[1], tmm[2], θunk, θstd), append!, mzs)

function NeXLUncertainties.asa(
    ::Type{DataFrame},
    unk::Material,
    std::Material,
    xrays::Vector{CharXRay},
    e0::Float64,
    toa::Float64;
    mc::Type{<:MatrixCorrection} = XPP,
    fc::Type{<:FluorescenceCorrection} = ReedFluorescence,
    coating::Type{<:CoatingCorrection} = Coating,
)
    flines = filter(cxr->energy(inner(cxr)) < e0, xrays)
    zafs = zafcorrection(mc, fc, coating, unk, std, flines, e0)
    asa(DataFrame, zafs..., toa, toa)
end

function aspure(unk::Material, cxr::CharXRay, e0::Float64, toa::Float64)
    zs = zafcorrection(XPP, ReedFluorescence, NullCoating, unk, pure(element(cxr)), [cxr], e0)
    return k(zs..., toa, toa)
end

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
"""
    abstract type Flow

Supertype for NormalisingFlow and InverseNormalisingFlow.
"""
abstract type Flow end

"""
    naive_init(T::Type{<:Flow}, rng::AbstractRNG, D::Int, invertibles::Vector{<:Invertible})

Get parameters to create a valid `Flow`, doesn't matter how performant it is.
"""
naive_init(T::Type{<:Flow}, rng::AbstractRNG, D::Int, invertibles::AbstractVector) =
    naive_init!(Vector{Float64}(nparams(T, D, invertibles)), T, rng, D, invertibles)

"""
    identity_init(T::Type{<:Flow}, D::Int, invertibles::Vector{<:Invertible})

Get parameters to create a `Flow`, whose transform is the identity function.
"""
identity_init(T::Type{<:Flow}, rng::AbstractRNG, D::Int, invertibles::AbstractVector) =
    identity_init!(Vector{Float64}(nparams(T, D, invertibles)), T, rng, D, invertibles)

# include("flows/normalising_flow.jl")
include("flows/inverse_normalising_flow.jl")

######################## Functions generic to all flows. ########################

dim(d::Flow) = dim(d.p0)
params(d::Flow) = vcat(params.(d.transforms)...)

# Determine number of parameters before instantiating a Flow.
nparams(::Type{<:Flow}, D::Int, ctors::AbstractVector) = sum(nparams.(ctors, D))

# Determine number of parameters using an instantiation of a Flow.
nparams(d::Flow) = sum(nparams.(d.transforms))

for flow in [:NormalisingFlow, :InverseNormalisingFlow]
    @eval @unionise function $flow(p0, θ::Vector{<:Real}, ctors::AbstractVector)
        @assert all(issubtype.(ctors, Invertible))
        @assert length(θ) == sum(nparams.(ctors, dim(p0)))
        pos, transforms = 1, Vector{Invertible}(length(ctors))
        for t in eachindex(ctors)
            Δ = nparams(ctors[t], dim(p0))
            transforms[t] = ctors[t](θ[pos:pos+Δ-1], dim(p0))
            pos += Δ
        end
        return $flow(p0, transforms)
    end
end