text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""
Module to enclose QType polynomial specific functionality. For reference see:
1. [_Robust, efficient computational methods for axially symmetric optical aspheres_ - G. W. Forbes, 2010](https://www.osapublishing.org/viewmedia.cfm?uri=oe-18-19-19700&seq=0)
2. [_Characterizing the shape of freeform optics_ - G. W. Forbes, 2012](https://www.osapublishing.org/viewmedia.cfm?uri=oe-20-3-2483&seq=0)

"""
module QType
using StaticArrays
using Plots
using ..OpticSim: QTYPE_PRECOMP

function F(m::Int, n::Int)::Float64
    # 2 eq A.13
    @assert m > 0
    if n === 0
        return Float64((m^2 * factorial2(big(2m - 3))) / (2^(m + 1) * factorial(big(m - 1))))
    elseif n > 0 && m === 1
        return Float64((4(n - 1)^2 * n^2 + 1) / (8 * (2n - 1)^2) + 11 * Int(n === 1) / 32)
    elseif n > 0 && m > 1
        χ = m + n - 2
        return Float64((2n * χ * (3 - 5m + 4n * χ) + m^2 * (3 - m + 4n * χ)) / ((m + 2n - 3) * (m + 2n - 2) * (m + 2n - 1) * (2n - 1)) * γ(m, n))
    else
        throw(ErrorException("Invalid n and m"))
    end
end

# 2 eq A.14
γ(a::Int, b::Int) = (factorial(big(b)) * factorial2(big(2a + 2b - 3))) / (2^(a + 1) * factorial(big(a + b - 3)) * factorial2(big(2b - 1)))

function G(m::Int, n::Int)::Float64
    # 2 eq A.15
    @assert m > 0
    if n === 0
        return Float64(factorial2(big(2m - 1)) / (2^(m + 1) * factorial(big(m - 1))))
    elseif n > 0 && m === 1
        return Float64(-((2 * n^2 - 1) * (n^2 - 1)) / (8 * (4 * n^2 - 1)) - Int(n === 1) / 24)
    elseif n > 0 && m > 0
        return Float64(-((2n * (m + n - 1) - m) * (n + 1) * (2m + 2n - 1)) / ((m + 2n - 2) * (m + 2n - 1) * (m + 2n) * (2n + 1)) * γ(m, n))
    else
        throw(ErrorException("Invalid n and m"))
    end
end

function factorial2(n::I)::I where {I<:Signed}
    # Defined in 2 appendix A
    @assert n <= 21 || n isa BigInt # not sure what number is the limit but this should do it...
    if n <= 0
        return I(1)
    else
        return n * factorial2(n - 2)
    end
end

function f(m::Int, n::Int, force::Bool = false)::Float64
    # 2 eq A.18b
    @assert m > 0
    if m <= QTYPE_PRECOMP && n < QTYPE_PRECOMP && !force
        return @inbounds PRECOMP_f[m, n + 1]
    else
        if n === 0
            return sqrt(F(m, 0))
        else
            return sqrt(F(m, n) - g(m, n - 1, force)^2)
        end
    end
end

function g(m::Int, n::Int, force::Bool = false)::Float64
    # 2 eq A.18a
    @assert m > 0
    if m <= QTYPE_PRECOMP && n < QTYPE_PRECOMP && !force
        return @inbounds PRECOMP_g[m, n + 1]
    else
        return G(m, n) / f(m, n, force)
    end
end

@inline function A(m::Int, n::Int, force::Bool = false)::Float64
    # 2 eq A.3a
    @assert m > 0
    if m <= QTYPE_PRECOMP && n < QTYPE_PRECOMP && !force
        return @inbounds PRECOMP_A[m, n + 1]
    else
        if m === 1 && n === 0
            return 2.0
        elseif m === 1 && n === 1
            return -4.0 / 3.0
        elseif n === 0 && m > 1
            return 2m - 1.0
        else
            return (2n - 1) * (m + 2n - 2) * (4n * (m + n - 2) + (m - 3) * (2m - 1)) / D(m, n)
        end
    end
end

@inline function B(m::Int, n::Int, force::Bool = false)::Float64
    # 2 eq A.3b
    @assert m > 0
    if m <= QTYPE_PRECOMP && n < QTYPE_PRECOMP && !force
        return @inbounds PRECOMP_B[m, n + 1]
    else
        if n === 1 && m === 1
            return -8.0 / 3.0
        elseif n === 0
            if m === 1
                return -1.0
            else
                return 2.0 * (1.0 - m)
            end
        else
            return -2.0 * (2n - 1) * (m + 2n - 3) * (m + 2n - 2) * (m + 2n - 1) / D(m, n)
        end
    end
end

@inline function C(m::Int, n::Int, force::Bool = false)::Float64
    # 2 eq A.3c
    @assert m > 0
    if n === 0
        return NaN
    elseif m <= QTYPE_PRECOMP && n < QTYPE_PRECOMP && !force
        return @inbounds PRECOMP_C[m, n + 1]
    else
        if m === 1 && n === 1
            return -11.0 / 3.0
        elseif m === 1 && n === 2
            return 0.0
        else
            return n * (2n - 3) * (m + 2n - 1) * (2m + 2n - 3) / D(m, n)
        end
    end
end

@inline function D(m::Int, n::Int)::Float64
    # 2 eq A.3d
    @assert m > 0
    return Float64((4n^2 - 1) * (m + n - 2) * (m + 2n - 3))
end

const PRECOMP_g = Matrix{Float64}(reshape(collect(g(m, n, true) for n in 0:(QTYPE_PRECOMP - 1) for m in 1:QTYPE_PRECOMP), (QTYPE_PRECOMP, QTYPE_PRECOMP)))
const PRECOMP_f = Matrix{Float64}(reshape(collect(f(m, n, true) for n in 0:(QTYPE_PRECOMP - 1) for m in 1:QTYPE_PRECOMP), (QTYPE_PRECOMP, QTYPE_PRECOMP)))
const PRECOMP_A = Matrix{Float64}(reshape(collect(A(m, n, true) for n in 0:(QTYPE_PRECOMP - 1) for m in 1:QTYPE_PRECOMP), (QTYPE_PRECOMP, QTYPE_PRECOMP)))
const PRECOMP_B = Matrix{Float64}(reshape(collect(B(m, n, true) for n in 0:(QTYPE_PRECOMP - 1) for m in 1:QTYPE_PRECOMP), (QTYPE_PRECOMP, QTYPE_PRECOMP)))
const PRECOMP_C = Matrix{Float64}(reshape(collect(C(m, n, true) for n in 0:(QTYPE_PRECOMP - 1) for m in 1:QTYPE_PRECOMP), (QTYPE_PRECOMP, QTYPE_PRECOMP)))

"""
    S(coeffs::SVector{NP1,T}, m::Int x::T) -> T

Evaluates ``\\sum_{n=0}^{N}c_n^mQ_n^m(x)`` where ``c_n^m`` is either an ``\\alpha`` or ``\\beta`` QType coefficient and ``m \\gt 0``.
"""
function S(coeffs::SVector{NP1,T}, m::Int, x::T)::T where {T<:Real,NP1}
    @assert m > 0
    if all(iszero.(coeffs))
        return zero(T)
    end
    N = NP1 - 1 # offset for indexing
    if N === 0
        return coeffs[1]
    end
    αₙ₊₂ = zero(T)
    αₙ₊₁ = zero(T)
    αₙ = zero(T)
    @inbounds for n in N:-1:0
        # 2 eq B.6
        αₙ = coeffs[n + 1] + (A(m, n) + B(m, n) * x) * αₙ₊₁ - C(m, n + 1) * αₙ₊₂
        if n > 0
            αₙ₊₂ = αₙ₊₁
            αₙ₊₁ = αₙ
        end
    end
    # 2 eq B.9
    if m === 1 && N > 2
        return αₙ / 2 - 2 / 5 * αₙ₊₂
    else
        return αₙ / 2
    end
end

"""
    dSdx(coeffs::SVector{NP1,T}, x::T) -> T

Evaluates ``\\frac{\\partial}{\\partial x}\\sum_{n=0}^{N}c_n^mQ_n^m(x)`` where ``c_n^m`` is either an ``\\alpha`` or ``\\beta`` QType coefficient and ``m \\gt 0``.
"""
function dSdx(coeffs::SVector{NP1,T}, m::Int, x::T)::T where {T<:Real,NP1}
    @assert m > 0
    if all(iszero.(coeffs))
        return zero(T)
    end
    N = NP1 - 1 # offset for indexing
    if N === 0
        return zero(T)
    end
    αₙ₊₂ = zero(T)
    αₙ₊₁ = zero(T)
    dαₙ₊₂ = zero(T)
    dαₙ₊₁ = zero(T)
    dαₙ = zero(T)
    @inbounds for n in N:-1:0
        # calcualte deriv
        Bmn = B(m, n)
        k = (A(m, n) + Bmn * x)
        Cmnp1 = C(m, n + 1)
        if n < N # derivative calcualted from N-1 to 0, so ignore first iteration
            # 2 eq B.11
            dαₙ = Bmn * αₙ₊₁ + k * dαₙ₊₁ - Cmnp1 * dαₙ₊₂
            if n > 0
                dαₙ₊₂ = dαₙ₊₁
                dαₙ₊₁ = dαₙ
            end
        end
        # calculate height for use in next iter deriv
        # 2 eq B.6
        αₙ = coeffs[n + 1] + k * αₙ₊₁ - Cmnp1 * αₙ₊₂
        αₙ₊₂ = αₙ₊₁
        αₙ₊₁ = αₙ
    end
    # 2 eq B.10
    if m === 1 && N > 2
        return dαₙ / 2 - 2 / 5 * dαₙ₊₂
    else
        return dαₙ / 2
    end
end

# for cases where m == 0 below:

function f0(n::Int, force::Bool = false)::Float64
    if n < QTYPE_PRECOMP && !force
        return @inbounds PRECOMP_f0[n + 1]
    elseif n === 0
        return 2.0
    elseif n === 1
        return sqrt(19) / 2
    else
        # 1 eq A.16
        return sqrt(n * (n + 1) + 3.0 - g0(n - 1, force)^2 - h0(n - 2, force)^2)
    end
end

function g0(n::Int, force::Bool = false)::Float64
    if n < QTYPE_PRECOMP && !force
        return @inbounds PRECOMP_g0[n + 1]
    elseif n === 0
        return -0.5
    else
        # 1 eq A.15
        return -(1 + g0(n - 1, force) * h0(n - 1, force)) / f0(n, force)
    end
end

function h0(n::Int, force::Bool = false)::Float64
    if n < QTYPE_PRECOMP && !force
        return @inbounds PRECOMP_h0[n + 1]
    else
        # 1 eq A.14
        return -(n + 2) * (n + 1) / (2 * f0(n, force))
    end
end

const PRECOMP_g0 = Vector{Float64}(collect(g0(n, true) for n in 0:(QTYPE_PRECOMP - 1)))
const PRECOMP_f0 = Vector{Float64}(collect(f0(n, true) for n in 0:(QTYPE_PRECOMP - 1)))
const PRECOMP_h0 = Vector{Float64}(collect(h0(n, true) for n in 0:(QTYPE_PRECOMP - 1)))

"""
    S0(coeffs::SVector{NP1,T}, x::T) -> T

Evaluates ``\\sum_{n=0}^{N}\\alpha_n^0Q_n^0(x)``.
"""
function S0(coeffs::SVector{NP1,T}, x::T)::T where {T<:Real,NP1}
    if all(iszero.(coeffs))
        return zero(T)
    end
    N = NP1 - 1 # offset for indexing
    if N === 0
        return coeffs[1]
    end
    k = T(2) - 4 * x
    # 1 eq 3.7a
    αₙ₊₂ = coeffs[N + 1]
    # 1 eq 3.7b
    αₙ₊₁ = coeffs[N] + k * αₙ₊₂
    for n in (N - 2):-1:0
        # 1 eq 3.8
        αₙ = coeffs[n + 1] + k * αₙ₊₁ - αₙ₊₂
        αₙ₊₂ = αₙ₊₁
        αₙ₊₁ = αₙ
    end
    # 1 eq 3.9
    return 2 * (αₙ₊₁ + αₙ₊₂) # a0 and a1
end

"""
    dS0dx(coeffs::SVector{NP1,T}, x::T) -> T

Evaluates ``\\frac{\\partial}{\\partial x}\\sum_{n=0}^{N}\\alpha_n^0Q_n^0(x)``.
"""
function dS0dx(coeffs::SVector{NP1,T}, x::T)::T where {T<:Real,NP1}
    if all(iszero.(coeffs))
        return zero(T)
    end
    N = NP1 - 1 # offset for indexing
    if N == 0
        return zero(T)
    end
    k = T(2) - 4 * x
    αₙ₊₂ = coeffs[N + 1]
    αₙ₊₁ = coeffs[N] + k * αₙ₊₂
    # 1 eq 3.10a
    dαₙ₊₂ = zero(T)
    # 1 eq 3.10b
    dαₙ₊₁ = -4 * αₙ₊₂
    for n in (N - 2):-1:0
        # 1 eq 3.11
        dαₙ = k * dαₙ₊₁ - dαₙ₊₂ - 4 * αₙ₊₁
        dαₙ₊₂ = dαₙ₊₁
        dαₙ₊₁ = dαₙ
        αₙ = coeffs[n + 1] + k * αₙ₊₁ - αₙ₊₂
        αₙ₊₂ = αₙ₊₁
        αₙ₊₁ = αₙ
    end
    # 1 eq 3.12
    return 2 * (dαₙ₊₁ + dαₙ₊₂) # a0 and a1
end

########################################
# for testing, not used in actual calculation because there is a more efficient method to find the sum directly (as in S and S0)

function P(m::Int, n::Int, x::T)::T where {T<:Real}
    if m === 0
        if n === 0
            return T(2)
        elseif n === 1
            return T(6) - 8 * x
        else
            # 1 eq 2.6
            return (T(2) - 4 * x) * P(0, n - 1, x) - P(0, n - 2, x)
        end
    else
        if n === 0
            return one(T) / 2
        elseif n === 1
            if m === 1
                # 2 eq A.5a
                return one(T) - x / 2
            else
                # 2 eq A.5b
                return m - 1 / 2 - (m - 1) * x
            end
        else
            # 2 eq A.2
            return (A(m, n - 1) + B(m, n - 1) * x) * P(m, n - 1, x) - C(m, n - 1) * P(m, n - 2, x)
        end
    end
end

function Q(m::Int, n::Int, x::T)::T where {T<:Real}
    if m === 0
        if n === 0
            return one(T)
        elseif n === 1
            return (T(13) - 16x) / sqrt(19)
        else
            # 1 eq 2.7
            return (P(0, n, x) - g0(n - 1) * Q(0, n - 1, x) - h0(n - 2) * Q(0, n - 2, x)) / f0(n)
        end
    else
        if n === 0
            return 1 / (2 * f(m, 0))
        else
            # 2 eq A.22
            return (P(m, n, x) - g(m, n - 1) * Q(m, n - 1, x)) / f(m, n)
        end
    end
end

function plot!(m::Int, n::Int)
    @assert m >= 0 && n >= 0
    us = 0:0.01:1
    vs = []
    for u in us
        if m === 0
            push!(vs, u^2 * (1 - u^2) * Q(m, n, u^2))
        else
            push!(vs, u^m * Q(m, n, u^2))
        end
    end
    Plots.plot!(us, vs, label = m === 0 ? "\$u^2(1-u^2)Q^{$m}_{$n}(u^2)\$" : "\$u^mQ^{$m}_{$n}(u^2)\$")
end

end # module QType

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

"""
    QTypeSurface{T,D,M,N} <: ParametricSurface{T,D}

Surface incorporating the QType polynomials - radius and conic are defined relative to absolute semi-diameter, QType terms are normalized according to the `normradius` parameter.
`T` is the datatype, `D` is the dimensionality, `M` and `N` are the maximum QType terms used.

The surface is centered at the origin and treated as being the cap of an infinite cylinder, thus creating a true half-space.
Outside of 0 <= ρ <= 1 the height of the surface is not necessarily well defined, so NaN may be returned.

```julia
QTypeSurface(semidiameter; radius = Inf, conic = 0.0, αcoeffs = nothing, βcoeffs = nothing, normradius = semidiameter)
```

`αcoeffs` and `βcoeffs` should be a vector of tuples of the form `(m, n, v)` where `v` is the value of the coefficient ``α_n^m`` or ``β_n^m`` respectively.

The sag is defined by the equation

```math
\\begin{aligned}
z(r,\\phi) = & \\frac{cr^2}{1 + \\sqrt{1 - (1+k)c^2r^2}} + \\frac{\\sqrt{1 + kc^2r^2}}{\\sqrt{1-(1+k)c^2r^2}} \\cdot \\\\
             & \\left\\{ \\rho^2(1-\\rho^2)\\sum_{n=0}^{N}\\alpha_n^0 Q_n^0 (\\rho^2) + \\sum_{m=1}^{M}\\rho^m\\sum_{n=0}^N \\left[ \\alpha_n^m\\cos{m\\phi} +\\beta_n^m\\sin{m\\phi}\\right]Q_n^m(\\rho^2) \\right\\}
\\end{aligned}
```

where ``\\rho = \\frac{r}{\\texttt{normradius}}``, ``c = \\frac{1}{\\texttt{radius}}``, ``k = \\texttt{conic}`` and ``Q_n^m`` is the QType polynomial index ``m``, ``n``.
"""
struct QTypeSurface{T,D,M,N} <: ParametricSurface{T,D}
    semidiameter::T
    curvature::T
    conic::T
    boundingcylinder::Cylinder{T,D}
    b0coeffs::SVector{N,T}
    dαcoeffs::SVector{M,SVector{N,T}} # m goes from 1:M while n goes from 0:N-1
    dβcoeffs::SVector{M,SVector{N,T}}
    normradius::T
    maxheight::T

    function QTypeSurface(semidiameter::T; radius::T = typemax(T), conic::T = zero(T), αcoeffs::Union{Nothing,Vector{Tuple{Int,Int,T}}} = nothing, βcoeffs::Union{Nothing,Vector{Tuple{Int,Int,T}}} = nothing, normradius::T = semidiameter) where {T<:Real}
        @assert !isnan(semidiameter) && !isnan(radius) && !isnan(conic)
        @assert semidiameter > zero(T)
        @assert one(T) - (1 / radius)^2 * (conic + one(T)) * semidiameter^2 > 0 "Invalid surface (conic/radius combination: $radius, $conic)"
        # work out maximum coefficient value
        N = -1
        M = 0
        if αcoeffs !== nothing
            αcoeffs = αcoeffs::Vector{Tuple{Int,Int,T}}
            for α in αcoeffs
                m, n, v = α
                if n > N
                    N = n
                end
                if m > M
                    M = m
                end
            end
        end
        if βcoeffs !== nothing
            βcoeffs = βcoeffs::Vector{Tuple{Int,Int,T}}
            for β in βcoeffs
                m, n, v = β
                if n > N
                    N = n
                end
                if m > M
                    M = m
                end
            end
        end
        # process the inputs to get parameter matrices
        α0coeffs = zeros(MVector{N + 1,T})
        αcoeffsproc = zeros(MMatrix{M,N + 1,T})
        βcoeffsproc = zeros(MMatrix{M,N + 1,T})
        if αcoeffs !== nothing
            for α in αcoeffs
                m, n, v = α
                @assert m >= 0 && n >= 0
                if m == 0
                    α0coeffs[n + 1] = v
                else
                    αcoeffsproc[m, n + 1] = v
                end
            end
        end
        if βcoeffs !== nothing
            for β in βcoeffs
                m, n, v = β
                @assert m >= 0 && n >= 0
                βcoeffsproc[m, n + 1] = v
            end
        end
        # calculate the α term coefficients for m = 0
        b0coeffs = zeros(MVector{N + 1,T})
        if N >= 0
            # 1 eq 3.4a
            bₙp2 = α0coeffs[N + 1] / QType.f0(N)
            b0coeffs[N + 1] = bₙp2
            if N > 0
                #1 eq 3.4b
                bₙp1 = (α0coeffs[N] - QType.g0(N - 1) * bₙp2) / QType.f0(N - 1)
                b0coeffs[N] = bₙp1
                if N > 1
                    for n in (N - 2):-1:0
                        # 1 eq 3.5
                        bₙ = (α0coeffs[n + 1] - QType.g0(n) * bₙp1 - QType.h0(n) * bₙp2) / QType.f0(n)
                        b0coeffs[n + 1] = bₙ
                        bₙp2 = bₙp1
                        bₙp1 = bₙ
                    end
                end
            end
        end
        # m==0 max value is <0.4 for any n, otherwise max is 1 for any m and n, for simplicity just sum everything though it may not be the tightest
        maxheight = zero(T)
        maxheight += 0.4 * sum(abs.(α0coeffs))
        maxheight += 0.4 * sum(abs.(b0coeffs))
        @inbounds @simd for m in 1:M
            maxheight += sum(abs.(αcoeffsproc[m, :]))
            maxheight += sum(abs.(βcoeffsproc[m, :]))
        end
        # calculate the α and β term coefficients
        dαcoeffs = zeros(MVector{M,SVector{N + 1,T}})
        dβcoeffs = zeros(MVector{M,SVector{N + 1,T}})
        for m in 1:M
            thisα = zeros(MVector{N + 1,T})
            thisβ = zeros(MVector{N + 1,T})
            if N >= 0
                fmN = QType.f(m, N)
                lastdαₙ = αcoeffsproc[m, N + 1] / fmN
                thisα[N + 1] = lastdαₙ
                lastdβₙ = βcoeffsproc[m, N + 1] / fmN
                thisβ[N + 1] = lastdβₙ
                for n in (N - 1):-1:0
                    gmn = QType.g(m, n)
                    fmn = QType.f(m, n)
                    # 2 eq B.4
                    dαₙ = (αcoeffsproc[m, n + 1] - gmn * lastdαₙ) / fmn
                    dβₙ = (βcoeffsproc[m, n + 1] - gmn * lastdβₙ) / fmn
                    lastdαₙ = dαₙ
                    lastdβₙ = dβₙ
                    thisα[n + 1] = dαₙ
                    thisβ[n + 1] = dβₙ
                end
            end
            dαcoeffs[m] = SVector{N + 1,T}(thisα)
            dβcoeffs[m] = SVector{N + 1,T}(thisβ)
        end
        NP1 = N + 1
        new{T,3,M,NP1}(semidiameter, 1 / radius, conic, Cylinder(semidiameter, interface = opaqueinterface(T)), SVector{NP1,T}(b0coeffs), SVector{M,SVector{NP1,T}}(dαcoeffs), SVector{M,SVector{NP1,T}}(dβcoeffs), normradius, maxheight) # TODO!! incorrect interface on cylinder
    end
end
export QTypeSurface

uvrange(::Type{QTypeSurface{T,D,M,N}}) where {T<:Real,D,M,N} = ((zero(T), one(T)), (-T(π), T(π))) # ρ and θ

semidiameter(z::QTypeSurface{T}) where {T<:Real} = z.semidiameter

function point(z::QTypeSurface{T,3,M,N}, ρ::T, θ::T)::SVector{3,T} where {T<:Real,M,N}
    rad = z.semidiameter
    r = ρ * rad
    # ρ is normalised [0, 1]
    # r is absolute
    r2 = r^2
    ρ2 = ρ^2
    c2 = z.curvature^2
    q = (one(T) + z.conic) * c2 * r2
    if q > one(T)
        return SVector{3,T}(NaN, NaN, NaN)
    end
    sqrtq = sqrt(one(T) - q)
    h = z.curvature * r2 / (one(T) + sqrtq)
    if N > 0
        # 2 eq B.1
        p = one(T) + z.conic * c2 * r2
        if p < zero(T)
            return SVector{3,T}(NaN, NaN, NaN)
        end
        u = ρ * rad / z.normradius
        u2 = u^2
        tot = u2 * (one(T) - u2) * QType.S0(z.b0coeffs, u2)
        @inbounds @simd for m in 1:M
            tot += u^m * (cos(m * θ) * QType.S(z.dαcoeffs[m], m, u2)::T + sin(m * θ) * QType.S(z.dβcoeffs[m], m, u2)::T)
        end
        h += (sqrt(p) / sqrtq) * tot
    end
    return SVector{3,T}(r * cos(θ), r * sin(θ), h)
end

function partials(z::QTypeSurface{T,3,M,N}, ρ::T, θ::T)::Tuple{SVector{3,T},SVector{3,T}} where {T<:Real,M,N}
    rad = z.semidiameter
    r = ρ * rad
    # ρ is normalised [0, 1]
    # r is absolute
    r2 = r^2
    ρ2 = ρ^2
    k = z.conic
    c = z.curvature
    c2 = c^2
    q = (one(T) + k) * c2 * r2
    if q > one(T)
        return SVector{3,T}(NaN, NaN, NaN), SVector{3,T}(NaN, NaN, NaN)
    end
    sqrtq = sqrt(one(T) - q)
    dhdρ = -((rad * c * r * sqrtq) / (q - one(T)))
    dhdθ = zero(T)
    if N > 0
        p = one(T) + k * c2 * r2
        if p < zero(T)
            return SVector{3,T}(NaN, NaN, NaN), SVector{3,T}(NaN, NaN, NaN)
        end
        sqrtp = sqrt(p)
        a = sqrtp / sqrtq
        dadρ = (c2 * (2k + one(T)) * rad * r) / (sqrtp * sqrtq^3)
        n = rad / z.normradius
        u = ρ * n
        u2 = u^2
        S0v = QType.S0(z.b0coeffs, u2)::T
        b = u2 * (one(T) - u2) * S0v
        dbdu = 2 * u * (one(T) - 2 * u2) * S0v + 2 * u^3 * (one(T) - u2) * QType.dS0dx(z.b0coeffs, u2)::T
        c = zero(T)
        dcdu = zero(T)
        if M > 0
            @inbounds @simd for m in 1:M
                Sα = QType.S(z.dαcoeffs[m], m, u2)::T
                Sβ = QType.S(z.dβcoeffs[m], m, u2)::T
                sinmθ = sin(m * θ)
                cosmθ = cos(m * θ)
                dhdθ += u^m * m * (-sinmθ * Sα + cosmθ * Sβ)
                c += u^m * (cosmθ * Sα + sinmθ * Sβ)
                dcdu += u^(m - 1) * (m * (cosmθ * Sα + sinmθ * Sβ) + 2 * u2 * (cosmθ * QType.dSdx(z.dαcoeffs[m], m, u2)::T + sinmθ * QType.dSdx(z.dβcoeffs[m], m, u2)::T))
            end
        end
        dhdθ *= a
        dhdρ += dadρ * (b + c) + a * n * (dbdu + dcdu) # want the derivative wrt ρ, not u
    end
    cosθ = cos(θ)
    sinθ = sin(θ)
    pu = SVector{3,T}(rad * cosθ, rad * sinθ, dhdρ)
    pv = SVector{3,T}(r * -sinθ, r * cosθ, dhdθ)
    return pu, pv
end

function normal(z::QTypeSurface{T,D,M,N}, ρ::T, θ::T)::SVector{3,T} where {T<:Real,D,M,N}
    du, dv = partials(z, ρ, θ)
    if ρ == zero(T) && norm(dv) == zero(T)
        # in cases where there is no δθ at ρ = 0 (i.e. anything which is rotationally symetric)
        # then we get some big problems, hardcoding this case solves the problems
        return SVector{3,T}(0, 0, 1)
    end
    return normalize(cross(du, dv))
end

function uv(z::QTypeSurface{T,3,M,N}, p::SVector{3,T})::SVector{2,T} where {T<:Real,M,N}
    # avoid divide by zero for ForwardDiff
    ϕ = NaNsafeatan(p[2], p[1])
    if p[1] == zero(T) && p[2] == zero(T)
        ρ = zero(T)
    else
        ρ = sqrt(p[1]^2 + p[2]^2) / semidiameter(z)
    end
    return SVector{2,T}(ρ, ϕ)
end

function onsurface(surf::QTypeSurface{T,3,M,N}, p::SVector{3,T}) where {T<:Real,M,N}
    ρ, θ = uv(surf, p)
    if ρ > one(T)
        return false
    else
        surfpoint = point(surf, ρ, θ)
        return samepoint(p[3], surfpoint[3])
    end
end

function inside(surf::QTypeSurface{T,3,M,N}, p::SVector{3,T}) where {T<:Real,M,N}
    ρ, θ = uv(surf, p)
    if ρ > one(T)
        return false
    else
        surfpoint = point(surf, ρ, θ)
        return p[3] < surfpoint[3]
    end
end

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

# Assumes the ray has been transformed into the canonical coordinate frame which has the vertical axis passing through (0,0,0) and aligned with the z axis.
function surfaceintersection(surf::AcceleratedParametricSurface{T,3,QTypeSurface{T,3,M,N}}, r::AbstractRay{T,3}) where {T<:Real,M,N}
    cylint = surfaceintersection(surf.surface.boundingcylinder, r)
    if cylint isa EmptyInterval{T}
        return EmptyInterval(T)
    else
        if doesintersect(surf.triangles_bbox, r) || inside(surf.triangles_bbox, origin(r))
            surfint = triangulatedintersection(surf, r)
            if !(surfint isa EmptyInterval{T})
                return intervalintersection(cylint, surfint)
            end
        end
        # hasn't hit the surface
        if lower(cylint) isa RayOrigin{T} && upper(cylint) isa Infinity{T}
            if inside(surf.surface, origin(r))
                return Interval(RayOrigin(T), Infinity(T))
            else
                return EmptyInterval(T)
            end
            # otherwise check that the intersection is underneath the surface
        else
            p = point(closestintersection(cylint, false))
            ρ, ϕ = uv(surf, p)
            surfpoint = point(surf.surface, ρ, ϕ)
            if p[3] < surfpoint[3]
                return cylint # TODO!! UV (and interface) issues?
            else
                return EmptyInterval(T)
            end
        end
    end
end

function AcceleratedParametricSurface(surf::S, numsamples::Int = 17; interface::NullOrFresnel{T} = NullInterface(T)) where {T<:Real,N,S<:QTypeSurface{T,N}}
    # Zernike users ρ, ϕ uv space so need to modify extension of triangulation
    a = AcceleratedParametricSurface(surf, triangulate(surf, numsamples, true, false, true, false), interface = interface)
    emptytrianglepool!(T)
    return a
end

function BoundingBox(surf::QTypeSurface{T,3}) where {T<:Real}
    xmin = -semidiameter(surf)
    xmax = semidiameter(surf)
    ymin = -semidiameter(surf)
    ymax = semidiameter(surf)
    # curvature only goes one way
    q = one(T) - (one(T) + surf.conic) * surf.curvature^2 * surf.semidiameter^2
    if q < zero(T)
        throw(ErrorException("The surface is invalid, no bounding box can be constructed"))
    end
    hmax = surf.curvature * surf.semidiameter^2 / (one(T) + sqrt(q))
    if hmax > zero(T)
        zmax = hmax + surf.maxheight
        zmin = -surf.maxheight
    else
        zmax = surf.maxheight
        zmin = hmax - surf.maxheight
    end
    return BoundingBox(xmin, xmax, ymin, ymax, zmin, zmax)
end


text:
module SampleMultiplePackage1

using Test

test() = @test true

end  # module SampleMultiplePackage1

ext: jl


text:
#= Matrix in ModelingToolkit
=#

using ModelingToolkit

@variables x[1:2, 1:2]


text:
"""
`rowcoord(df::DataFrame)`

Returns the :latitude and :longitude columns of a dataframe
as a tuple of (latitude, longitude). If the dataframe has more
than 1 row a vector of these tuples is returned.
"""
function rowcoord(df::DataFrame)
  if nrow(df) == 1
    return singlerowcoord(df)
  else
    return [singlerowcoord(df[i, :]) for i=1:nrow(df)]
  end
end

# Returns Coordinate from a single dataframe row
function singlerowcoord(df::DataFrame)::Coordinate
  df[:latitude][1], df[:longitude][1]
end

# Convert distance from default meters
function dist_convert(distance::AbstractFloat,
                      newunit::Unitful.FreeUnits,
                      returntype::Type=AbstractFloat)
   v = 1.0u"m" * distance |> newunit
   return v.val |> returntype
end


text:
include("two_link.jl")
using TwoLink
using DynamicMovementPrimitives
using Requires

cpoints = [0.5 -0.5;
0.5 0;
1.5 0;
1.5 -0.5]

y,yd,ydd = connect_points(cpoints,20)
q = inverse_kin(y,:up)

# Setup ============================================
Nbasis = 20
αz = 25.
αx = 1.

opts = DMPopts(Nbasis,αx,αz, :canonical)

T = size(y,1)
t = LinRange(0,T,T)
τ = t[end]/3
h = t[2]-t[1]
qd = centraldiff(q) / h
qdd = centraldiff(qd) / h
g=q[end,:][:]
dmp = fit(q,qd,qdd,t,opts)

to,qo,zo,xo = solve(dmp)
yo = forward_kin(qo)

@require Plots begin
fig1 = plot(q, xlabel="Time", title="Joint space reference and output", lab=["Reference_"].*["x" "y"])
plot!(qo, l=:dash, lab=["Fit_"].*["x" "y"])

fig2 = plot(y, xlabel="Time", title="Cartesian space reference and output", lab=["Reference_"].*["x" "y"])
plot!(yo, l=:dash, lab=["Fit_"].*["x" "y"])

fig3 = plot(q[:,1],q[:,2], xlabel="x", ylabel="y", title="Joint space reference and output", lab="Reference")
plot!(qo[:,1],qo[:,2], l=:dash, lab="Fit")

fig4 = plot(y[:,1],y[:,2], xlabel="x", ylabel="y", title="Cartesian space reference and output", lab="Reference")
plot!(yo[:,1],yo[:,2], l=:dash, lab="Fit")

subplot(fig1,fig2,fig3,fig4); gui()
end


text:
# Tests for Von-Mises distribution

using Distributions
using Test

function test_vonmises(μ::Float64, κ::Float64)
    d = VonMises(μ, κ)
    @test length(d) == 1
    @test mean(d) == μ
    @test median(d) == μ
    @test mode(d) == μ
    @test d == typeof(d)(params(d)...,d.I0κx)
    @test d == deepcopy(d)
    @test partype(d) == Float64
    # println(d)

    # conversions
    @test typeof(convert(VonMises{Float32}, d)) == VonMises{Float32}

    # Support
    @test support(d) == RealInterval(d.μ-π,d.μ+π)
    @test pdf(d, d.μ-2π) == 0.0
    @test pdf(d, d.μ+2π) == 0.0

end

## General testing

for (μ, κ) in [(2.0, 1.0),
               (2.0, 5.0),
               (3.0, 1.0),
               (3.0, 5.0),
               (5.0, 2.0)]

    test_vonmises(μ, κ)
end


text:
using Test

function sum(arr)
    total = 0
    for x in arr
        total += x
    end
    return total
end

@test sum([1,2,3,4]) == 10


text:
using Statistics
using StatsModels
using GLM
using DataFrames
using CSV
using Lathe

data = DataFrame(CSV.File("automobili.csv"))
data = Lathe.preprocess.predict(data, :model)
dataTrain, dataTest = Lathe.preprocess.TrainTestSplit(data,.80)
fm = @formula(cena ~ kubikaza+godiste+model)
linearRegressor = lm(fm, dataTrain)

predvidjenaCena = predict(linearRegressor, dataTest)
println("Procenjene cene vozila: ")
for i in 1:length(predvidjenaCena)
   println("X$i: $(dataTest.proizvodjac[i])  $(dataTest.model[i]) $(dataTest.kubikaza[i]) $(dataTest.godiste[i]) Cena: $(dataTest.cena[i]) Procena: $(predvidjenaCena[i])")
end

# Procena greske
# Racunanje i ispis gresaka
errorsTest = dataTest.cena-predvidjenaCena
println()
println("Spisak svih gresaka pri testiranju je: $(round.(abs.(errorsTest); digits = 2))")
# Racunanje i ispis proseka greske
absMeanErrorTest = mean(abs.(errorsTest))
println("Prosecna absolutna greska pri testiranju je: $absMeanErrorTest evra")
# Racunanje prosecne relativne greske u procentima MAPE
mapeTest = mean(abs.(errorsTest./dataTest.cena))
println("Prosecna relativna greska pri testiranju je: $(mapeTest*100) procenata")

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using CSV
using DataFrames
using Turing
using Dates
using StatsBase
using Distributions
using CovidSurvey
using Underscores
# using Gadfly
using StatsPlots
using LaTeXStrings
using NamedTupleTools
using StatsFuns
using Random
using TransformVariables, LogDensityProblems
using Optim
using ReverseDiff
using BSON
import TimeSeriesDecompositions as TSD
plotlyjs()

# using DataFramesMeta
# plotlyjs()
const dateform = "yyyymmdd"
include( projectdir("scripts/survey/rawsurvey_utils.jl") )

## ============================================================================
function plot_contacts_per_day(df; region=0)
    gp = groupby(df[:,[:date]], :date, sort=true)
    cb = combine(gp, nrow)
    p = plot(cb.date, cb.nrow, ylab="# contacts", xticks=:native)
    savefig(p, projectdir("scripts/survey/","contacts_per_day_region=$region.html"))
    return p
end

function main()
    plotlyjs()
    for region in 1:5
        @info "load data region $region"
        rawcontacts = readcontacts( "/home/and/data/covidsurvey/rawcontacts.csv", select_region=region )
        plot_contacts_per_day(rawcontacts; region)
        threshold = (family=50, colleagues=100, friends=100, strangers=1000)
        contacts = filteroutliers(rawcontacts, threshold) #NOTE National: p = 99.9

        @info "zero inflated regression region $region"
        surveyproblem = SurveyInferenceProblem(
            ZINegBinomialProblem,
            ctype,
            contacts
        )

        results = optimize(surveyproblem, LBFGS(),
            Optim.Options(show_trace=true, iterations=2000);
        )

        contactrates = DataFrame(results)
        save(projectdir("data","raw_contact_rates_region=$region.csv"), contactrates)

        let
            fname = projectdir("scripts/survey/","survey_inference_results_method=$(results.problemtype)_region=$region.bson")
            save(fname, θ = results)
        end

        plot_confidence_timeseries(results)

        # pyplot()
        # p = plot_survey_histograms(results, day=nothing)
        # savefig(p, projectdir("figures/histogram_survey_data_vs_fake_data_region=$region.png"))
        # plotlyjs()

        @info "decomposition region $region"
        slowtrends = TSD.decompose(contactrates, Optim.Options(show_trace=true, iterations=2000); fdir=projectdir("figures/"))
        save(projectdir("data","smoothed_contact_rates_region=$region.csv"), slowtrends)
    end
end

# main()

## =============================================================================
# apply simple average

function rollingmean(xs, window)
    n  = length(xs)-window+1
    ys = Vector{Float64}(undef, n)
    for i in 1:n
        ys[i] = mean(xs[i:i+window-1])
    end
    return ys
end

function smoothing(df, window)
    @assert isodd(window)
    Δ = window ÷ 2
    smoothed = DataFrame( :date => df.date[1+Δ:end-Δ] )
    cols = filter(x->nonmissingtype( eltype( df[!,x] ) ) <: Number, names(df))
    for col in cols
        smoothed[!,col] = rollingmean(df[:,col], window)
    end
    smoothed
end

ps = (
    threshold = (family=50, colleagues=100, friends=100, strangers=1000),
    window = 7,
    fname = normpath( homedir(), "data/covidsurvey/rawcontacts.csv"),
    ctypes = [:family, :colleagues, :friends, :strangers],
)

function main(ps)
    contacts = @_ ps.fname |>
        readcontacts |>
        filteroutliers(__, ps.threshold) |>
        groupby(__, :date, sort = true) |>
        combine(__, ps.ctypes .=> sum, nrow, renamecols = false) |>
        DataFrames.select(__, All(), AsTable(ps.ctypes) => ByRow(sum) => :total) |>
        transform(__, [ps.ctypes..., :total] .=> ( x -> x ./ __.nrow ), renamecols = false) |>
        disallowmissing |>
        smoothing(__, ps.window)

    save(projectdir("data/contacts/dk","smoothed_contact_rates_window=$(ps.window).csv"), contacts)
    return contacts
end

contacts = main(ps)

## =========================================================================================================
@info "load data"
rawcontacts = readcontacts( normpath( homedir(), "data/covidsurvey/rawcontacts.csv"), select_region=0 )
p = plot_contacts_per_day(rawcontacts)
threshold = 99.9 #NOTE National
contacts = filteroutliers(rawcontacts, threshold)
contacts[!,:total] .= 0
for c in ctype
    contacts[!,:total] += contacts[:,c]
end
contacts = contacts[:, [:total, :strangers, :date]]
## =========================================================================================================
@info "zero inflated regression"
surveyproblem = SurveyInferenceProblem(
    ZINegBinomialProblem,
    datacols(contacts),
    contacts
)

results = optimize(surveyproblem, LBFGS(),
    Optim.Options(show_trace=true, iterations=2000);
)

contactrates = DataFrame(results)
save(projectdir("data","raw_total_contact_rates.csv"), contactrates)

let
    fname = projectdir("scripts/survey/","survey_inference_results_method=$(results.problemtype)_TOTAL.bson")
    save(fname, θ = results)
end

plot_confidence_timeseries(results)

pyplot()
p = plot_survey_histograms(results, day=nothing)
# savefig(p, projectdir("figures/histogram_survey_data_vs_fake_data_region=$region.png"))
plotlyjs()

@info "decomposition region $region"
slowtrends = TSD.decompose(contactrates, Optim.Options(show_trace=true, iterations=2000); fdir=projectdir("figures/"))
save(projectdir("data","smoothed_total_contact_rates.csv"), slowtrends)

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
using Brusselator

using Plots, Printf

using FFTW: irfft

# ## Let's prescibe parameter values and solve the PDE
#
# We are now ready to write up a program that sets up parameter values, constructs
# the problem `prob`, # time steps the solutions `prob.sol` and plots it.

# ## Choosing a device: CPU or GPU

dev = CPU()    # Device (CPU/GPU)
nothing # hide

# ## Numerical parameters and time-stepping parameters

     nx = 512            # grid resolution
stepper = "RK4"          # timestepper
     dt = 0.01           # timestep
 nsteps = 2500           # total number of time-steps
 nsubs  = 25             # number of time-steps for intermediate logging/plotting (nsteps must be multiple of nsubs)
nothing # hide

# ## Physical parameters

E = 1.4

L = 137.37     # Domain length

# Parameters used by Tzou et al. (2013)
ε = 0.1
μ = 25
ρ = 0.178

D_c = ((sqrt(1 + E^2) - 1) / E)^2
B_H = (1 + E * sqrt(D_c))^2

B = B_H + ε^2 * μ
D = D_c + ε^2 * ρ

kc = sqrt(E / sqrt(D))

nothing # hide

# ## Construct the `struct`s and you are ready to go!
# Create a `grid` and also `params`, `vars`, and the `equation` structs. Then
# give them all as input to the `FourierFlows.Problem()` constructor to get a
# problem struct, `prob`, that contains all of the above.

grid = OneDGrid(dev, nx, L)
params = Brusselator.Params(B, D, E)
vars = Brusselator.Vars(dev, grid)
equation = Brusselator.Equation(dev, params, grid)

prob = FourierFlows.Problem(equation, stepper, dt, grid, vars, params, dev)
nothing #hide

get_u(prob) = irfft(prob.sol[:, 1], prob.grid.nx)
get_v(prob) = irfft(prob.sol[:, 2], prob.grid.nx)

u_solution = Diagnostic(get_u, prob; nsteps=nsteps+1)
diags = [u_solution]

# ## Visualizing the simulation

# We define a function that plots $u$ and $v$.

function plot_output(prob)
  plot_u = plot(grid.x, vars.u,
                 color = :blue,
                legend = false,
             linewidth = 2,
                 alpha = 0.7,
                 xlims = (-L/2, L/2),
                 ylims = (E-1.5, E+1.5),
                xlabel = "x",
                ylabel = "u")

  plot_v = plot(grid.x, vars.v,
                 color = :red,
                legend = false,
             linewidth = 2,
                 alpha = 0.7,
                 xlims = (-L/2, L/2),
                 ylims = (B/E-1.5, B/E+1.5),
                xlabel = "x",
                ylabel = "v")

  title = plot(title = "Brusselator", grid = false, showaxis = false, bottom_margin = -20Plots.px)

  return plot(title, plot_u, plot_v,
           layout = @layout([A{0.01h}; [B; C]]),
             size = (600, 400),
              dpi = 150)
end
nothing # hide

# ## Setting initial conditions

# We reproduce here (to our best knowledge possible) the initial condition used
# by Tzou et al. (2013).

l = 28;
θ = @. (grid.x > -l/2) & (grid.x < l/2)

au, av = -(E^2 + kc^2) / B, -1
cu, cv = -E * (E + im) / B, 1

u0 = @.  E  + ε * real( au * exp(im * kc * grid.x) * θ + cu * (1 - θ) )
v0 = @. B/E + ε * real( av * exp(im * kc * grid.x) * θ + cv * (1 - θ) )

set_uv!(prob, u0, v0)

# Let's plot the initial condition.

plot_output(prob)

# ## Time-stepping the `Problem` forward

# We time-step the `Problem` forward in time. We update variables by calling
# `updatevars!()` and we also update the plot. We enclose the `for` loop in
# an `@animate` macro to produce an animation of the solution.

startwalltime = time()

p = plot_output(prob)

anim = @animate for j=0:Int(nsteps/nsubs)

  log = @sprintf("step: %04d, t: %.1f, walltime: %.2f min",
  prob.clock.step, prob.clock.t, (time()-startwalltime)/60)

  if j%(500/nsubs)==0; println(log) end

  updatevars!(prob)

  p[2][1][:y] = vars.u    # updates the plot for u
  p[2][:title] = "t = "*@sprintf("%.1f", prob.clock.t) # updates time in the title
  p[3][1][:y] = vars.v    # updates the plot for v

  stepforward!(prob, diags, nsubs)
end

mp4(anim, "brusselator.mp4", fps=18)

# Now let's plot a space-time diagram of the $u$ solution.

t = u_solution.t[1:nsteps+1]
U_xt = zeros(grid.nx, nsteps+1)
[ U_xt[:, j] = u_solution.data[j] for j=1:nsteps+1 ]

heatmap(grid.x, t, U_xt',
             c = :grayC,
         clims = (1, 3))


text:
module Queryverse

using Reexport

@reexport using DataValues
import IterableTables
@reexport using Query
@reexport using QueryTables
@reexport using DataFrames
@reexport using FileIO
@reexport using ExcelFiles
@reexport using StatFiles
@reexport using CSVFiles
@reexport using FeatherFiles
@reexport using ParquetFiles
@reexport using VegaLite
@reexport using DataVoyager

export DV

export @tee

macro tee(ex)
    quote
        x -> begin
            x |> $(esc(ex))
            return x
        end
    end
end

const DV = DataValues.DataValue

end # module


text:
import Base: +, -, *, /, .*, ./, ^, !=, <, >, ==, .<, .>, .>=, .<=, norm

export Point, O, randompoint, randompointarray
"""
The Point type holds two coordinates. Currently it's immutable, so remember not try to
change the values of the x and y values directly.
"""
immutable Point
   x::Float64
   y::Float64
end

"""
O is a shortcut for the current origin, `0/0`
"""
const O = Point(0, 0)

+(z::Number, p1::Point)              = Point(p1.x + z,    p1.y + z)
+(p1::Point, z::Number)              = Point(p1.x + z,    p1.y + z)
+(p1::Point, p2::Point)              = Point(p1.x + p2.x, p1.y + p2.y)
-(p1::Point, p2::Point)              = Point(p1.x - p2.x, p1.y - p2.y)
-(p::Point)                          = Point(-p.x,        -p.y)
-(z::Number, p1::Point)              = Point(p1.x - z,    p1.y - z)
-(p1::Point, z::Number)              = Point(p1.x - z,    p1.y - z)
*(k::Number, p2::Point)              = Point(k * p2.x,    k * p2.y)
*(p2::Point, k::Number)              = Point(k * p2.x,    k * p2.y)
/(p2::Point, k::Number)              = Point(p2.x/k,      p2.y/k)
.*(k::Number, p2::Point)             = Point(k * p2.x,    k * p2.y)
.*(p2::Point, k::Number)             = Point(k * p2.x,    k * p2.y)
.*(p1::Point, p2::Point)             = Point(p1.x * p2.x, p1.y * p2.y)
./(p2::Point, k::Number)             = Point(p2.x/k,      p2.y/k)
^(p::Point, e::Integer)              = Point(p.x^e,       p.y^e)
^(p::Point, e::Float64)              = Point(p.x^e,       p.y^e)
inner(p1::Point, p2::Point)          = p1.x * p2.x + p1.y * p2.y

# comparisons (perform lexicographically)

isequal(p1::Point, p2::Point)         = (p1.x == p2.x && p1.y == p2.y )
isless(p1::Point, p2::Point)          = (p1.x < p2.x || (p1.x == p2.x && p1.y < p2.y) )
!=(p1::Point, p2::Point)              = !isequal(p1, p2)
<(p1::Point, p2::Point)               = isless(p1,p2)
>(p1::Point, p2::Point)               = p2 < p1
==(p1::Point, p2::Point)              = isequal(p1, p2)
.<(p1::Point, p2::Point)              = p1 < p2
.>(p1::Point, p2::Point)              = p2 < p1
.>=(p1::Point, p2::Point)             = p1 <= p2
.<=(p1::Point, p2::Point)             = p2 <= p1
cmp(p1::Point, p2::Point)             = p1<p2 ? -1 : p2  <p1 ? 1 : 0

function randomordinate(low, high)
    low + rand() * abs(high - low)
end

"""
    randompoint(lowpt, highpt)

Return a random point somewhere inside the rectangle defined by the two points.
"""
function randompoint(lowpt, highpt)
  Point(randomordinate(lowpt.x, highpt.x), randomordinate(lowpt.y, highpt.y))
end

"""
    randompoint(lowx, lowy, highx, highy)

Return a random point somewhere inside a rectangle defined by the four values.
"""
function randompoint(lowx, lowy, highx, highy)
    Point(randomordinate(lowx, highx), randomordinate(lowy, highy))
end

"""
    randompointarray(lowpt, highpt, n)

Return an array of `n` random points somewhere inside the rectangle defined by two points.
"""
function randompointarray(lowpt, highpt, n)
  array = Point[]
  for i in 1:n
    push!(array, randompoint(lowpt, highpt))
  end
  array
end

"""
    randompointarray(lowx, lowy, highx, highy, n)

Return an array of `n` random points somewhere inside the rectangle defined by the four coordinates.

"""
function randompointarray(lowx, lowy, highx, highy, n)
    array = Point[]
    for i in 1:n
        push!(array, randompoint(lowx, lowy, highx, highy))
    end
    array
end

"""
Find the norm of two points (two argument form).

    norm(p1::Point, p2::Point)
"""
function norm(p1::Point, p2::Point)
    norm([p1.x, p1.y] - [p2.x, p2.y])
end

"""
Find the distance between a point `p` and a line between two points `a` and `b`.

    point_line_distance(p::Point, a::Point, b::Point)

"""
function point_line_distance(p::Point, a::Point, b::Point)
    # area of triangle
    area = abs(0.5 * (a.x * b.y + b.x * p.y + p.x * a.y - b.x * a.y - p.x * b.y - a.x * p.y))
    # length of the bottom edge
    dx = a.x - b.x
    dy = a.y - b.y
    bottom = sqrt(dx * dx + dy * dy)
    return area / bottom
end

"""
    midpoint(p1, p2)

Find the midpoint between two points.
"""
midpoint(p1::Point, p2::Point) = Point((p1.x + p2.x) / 2., (p1.y + p2.y) / 2.)

"""
Find midpoint between the first two elements of an array of points.

    midpoint(a)
"""
midpoint(pt::Array) = midpoint(pt[1], pt[2])

"""
Find intersection of two lines `p1`-`p2` and `p3`-`p4`

    intersection(p1, p2, p3, p4)

This returns a tuple: `(false, 0)` or `(true, Point)`.
"""
function intersection(p1, p2, p3, p4)
    flag = false
    ip = 0
    s1 = p2 - p1
    s2 = p4 - p3
    u = p1 - p3
    ip = 1 / (-s2.x * s1.y + s1.x * s2.y)
    s = (-s1.y * u.x + s1.x * u.y) * ip
    t = ( s2.x * u.y - s2.y * u.x) * ip
    if (s >= 0) && (s <= 1) && (t >= 0) && (t <= 1)
        if isapprox(ip, 0, atol=0.1)
            ip = p1 + (s1 * t)
            flag = true
        end
    end
    return (flag, ip)
end

# end


text:
module HoJBot

export start_bot

using CSV
using DataFrames
using Dates
using Discord
using Downloads
using ExpiringCaches
using Formatting
using JSON
using OrderedCollections
using Plots
using Pretend
using PrettyTables
using TimeZones
using UUIDs
using Pkg

import HTTP
import JSON3
import StructTypes

const COMMAND_PREFIX = get(ENV, "HOJBOT_COMMAND_PREFIX", ",")

include("constants.jl")
include("dispatcher.jl")
include("util.jl")
include("discord.jl")
include("main.jl")

include("command/tz.jl")
include("command/j.jl")
include("command/gm.jl")
include("command/react.jl")
include("command/src.jl")

include("handler/reaction.jl")
include("handler/whistle.jl")
include("handler/mod.jl")

include("type/discourse.jl")
include("command/discourse.jl")
include("handler/discourse.jl")

include("type/ig.jl")
include("command/ig.jl")

end # module


text:
using MakeFshared
using Test

@testset "MakeFshared" begin

println("# printout function check")
@test printout() == nothing

println("# calc function in eins module check")
@test calc(Int32[4], Int32[5]) == 20
@test calc(Int32[8], Int32[8]) == 64
@test calc(Float32[9.0], Float32[2.0]) == 0
@test_throws MethodError calc("c", "a")

println("# foo function in zwei module check")
a = Int32[3]
b = Int32[5]
@test foo(a, b) == nothing
@test b[1] == 8

println("# moo function in zwei module check")
a = Int32[5, 7]
b = Int32[8, 9]
@test moo(a, b) == nothing
@test b[1] == 13
@test b[2] == 16

end


text:
using UnfoldMakie
using Test

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


text:
using MAGMA
using Test, Random, LinearAlgebra, CuArrays, CUDAnative, CUDAapi, CUDAdrv

using Test, LinearAlgebra, CuArrays, CUDAnative, CUDAapi, CUDAdrv

@testset "dense linear algebra" begin
    include("dense/svd.jl")
    include("dense/linearsystemsolver.jl")
end


text:
using DocStringExtensions,
      FillArrays,
      ForwardDiff,
      LinearAlgebra,
      RecipesBase,
      Reexport,
      UnPack

@reexport using MathematicalSystems

import CommonSolve: solve, init, solve!


text:
using AdaptiveFilters
using Test, Statistics

@testset "AdaptiveFilters.jl" begin
    y = sin.(0:0.1:100)

    yh = adaptive_filter(y, lr=0.01, order=4)
    @test mean(abs2, y-yh) < 1e-3

    yh = adaptive_filter(y, lr = 0.99, order=4)
    @test 1e-2 < mean(abs2, y-yh) < 1e-1

    yh = adaptive_filter(y, ExponentialWeight, order=2, lr = 0.01)
    @test mean(abs2, y-yh) < 1e-4

    @test length(y) == length(yh)

    y = sin.(0:0.1:100)
    yh = focused_adaptive_filter(y, (0.01,2), (2pi)/0.1, lr=0.01, order=4)
    @test mean(abs2, y-yh) < 1e-2

end


text:
abstract type AbstractQuadratureRule end

struct QuadratureRule{shape,N,T} <: AbstractQuadratureRule
    weights::Vector{T}
    points::Vector{Tensors.Vec{N,T}}
end

getweights(quad::QuadratureRule) = quad.weights
getpoints(quad::QuadratureRule) = quad.points

"""
function integrate(f::function,qr::QuadratureRule)
integrate function f on reference shape of quadrature qr
"""
function integrate(f::Function,qr::QuadratureRule)
    p = getpoints(qr); w = getweights(qr)
    int_val = zero(typeof(f(p[1])))
    for (i,x) in enumerate(p)
        int_val += f(x)*w[i]
    end
    int_val
end

#quadratures
struct GrundmannMoeller <: AbstractQuadratureRule end
struct Strang <: AbstractQuadratureRule end
struct DefaultQuad <: AbstractQuadratureRule end
struct GaussLegendre <: AbstractQuadratureRule end
struct Gauss <: AbstractQuadratureRule end

function (::Type{QuadratureRule{Triangle}})(quad_type::DefaultQuad, order::Int)
    if order <= 6
        return QuadratureRule{Triangle}(Strang(),order)
    elseif order > 6 && isodd(order)
        s = Int((order-1)/2)
        return QuadratureRule{Triangle}(GrundmannMoeller(),s)
    else
        throw(ArgumentError("Quadrature rule of order $order not available"))
    end
end

function (::Type{QuadratureRule{Tetrahedron}})(quad_type::DefaultQuad, order::Int)
    if order <= 4
        return QuadratureRule{Tetrahedron}(Gauss(),order)
    elseif order > 6 && isodd(order)
        s = Int((order-1)/2)
        return QuadratureRule{Tetrahedron}(GrundmannMoeller(),s)
    else
        throw(ArgumentError("Quadrature rule of order $order not available"))
    end
end

function (::Type{QuadratureRule{HyperCube{dim}}})(quad_type::DefaultQuad, order::Int) where {dim}
    return QuadratureRule{HyperCube{dim}}(GaussLegendre(),order)
end

# Get GaussLegendre weigths and points from FastGaussQuadrature
function (::Type{QuadratureRule{Segment}})(quad_type::GaussLegendre, order::Int)
    points, weights = FastGaussQuadrature.gausslegendre(order)
    # Shift interval from (-1,1) to (0,1)
    weights *= 0.5
    points = points .+ 1.0; points /= 2.0
    return QuadratureRule{Segment,Float64}(weights, [Tensors.Tensor{1,1}([x]) for x in points])
end

function (::Type{QuadratureRule{Segment}})(quad_type::DefaultQuad, order::Int)
    return QuadratureRule{Segment}(GaussLegendre(), order)
end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
# Stubs for conditionally defined functions

"""
Wannerize the obtained bands using wannier90. By default all converged
bands from the `scfres` are employed (change with `n_bands` kwargs)
and `n_wannier = n_bands` wannier functions are computed using
random Gaussians as guesses. All keyword arguments supported by
Wannier90 for the disentanglement may be added as keyword arguments.
The function returns the `fileprefix`.

!!! warning "Experimental feature"
    Currently this is an experimental feature, which has not yet been tested
    to full depth. The interface is considered unstable and may change
    incompatibly in the future. Use at your own risk and please report bugs
    in case you encounter any.
"""
function run_wannier90 end


text:
module StringArrayEditor

export
    load,
    Line,
    Lines,
    Collection,
    insert!,
    delete!,
    append!,
    replace!,
    move!,
    moveafter!,
    Range,
    Ranges,
    value,
    save,
    map!,
    map,
    prev,
    next,
    +,-,
    peek

import Base:show,delete!,insert!,copy,append!,==,length,iterate,getindex,lastindex,in,replace!,findfirst,findnext,map!,map,+,-

using PrePostCall

using InteractiveUtils
import InteractiveUtils.edit

"""
    Reference

Every new reference i.e. line or range must be a subtype of Reference

"""
abstract type Reference end

struct File
    data::Vector{String}
    references::Vector{Reference}
    function File(data::Vector{T}) where T<: AbstractString
        new(data,Vector{Reference}())
    end
end

"""
    load(path::AbstractString)

Load the file from `path` into memory and return a File object.
"""
function load(path::AbstractString)
    return File(readlines(path))
end

"""
    save(f::File,path::AbstractString)

Save file `f` to `path`.
"""
function save(f::File,path::AbstractString)
    open(path,"w") do d
        for l in f.data
            write(d,l,"\n")
        end
    end
end

function show(io::IO,f::File)
    print(io,"File(<Lines $(length(f.data))>,<Ref $(length(f.references))>)")
end

writeToFile(io::IOStream,l::String) = write(io,l)
function writeToFile(io::IOStream,lines::Vector{String})
    for l in lines
        write(io,l,"\n")
    end
end

function deletefromreferences!(f::File,del::Reference)
    for (idx,ref) in enumerate(f.references)
        if ref === del
            deleteat!(f.references,idx)
            return nothing
        end
    end
    error("Reference $del not found in $f")
end

function notify_delete(f::File,del::Reference)
    for ref in f.references|>copy
        ref === del && continue # Skip del
        notify_delete(ref,del)
    end
end

function notify_insert(f::File,ins::Reference,lines::Int)
    for ref in f.references
        ref === ins && continue # Skip ins
        notify_insert(ref,ins,lines)
    end
end

function notify_append(f::File,app::Reference,lines::Int)
    for ref in f.references
        ref === app && continue # Skip ins
        notify_append(ref,app,lines)
    end
end

function notify_replace(f::File,rep::Reference,lines::Int)
    for ref in f.references|>copy
        ref === rep && continue # Skip ins
        notify_replace(ref,rep,lines)
    end
end

@pre alive(self::T) where T<:Reference = isdestroyed(self) && error("Can't perform operation on deleted $T.")

include("./Line.jl")
include("./Range.jl")
include("./Collection.jl")
include("./Search.jl")
include("./misc.jl")

end


text:
module Xyce
using CxxWrap
using XyceWrapper_jll
@wrapmodule(XyceWrapper_jll.xycelib)

function __init__()
    @initcxx
end

function simulate(netlist; type=TRAN, name="handler", vectors=["TIME", "V(*)", "I(*)"])
    cxxvec = StdVector(StdString.(vectors))
    oh = OutputHandler(name, type, cxxvec)
    x = GenCouplingSimulator()
    # this avoids aborting when Xyce encounters an error
    set_report_handler()
    arg = ["xyce", netlist]
    argch = Base.cconvert.(CxxPtr{CxxChar}, arg)
    argv = Base.cconvert(CxxPtr{CxxPtr{CxxChar}}, argch)
    initializeEarly(x, 2, argv)
    addOutputInterface(x, CxxPtr(oh))
    initializeLate(x)
    runSimulation(x)
    names = String.(getFieldnames(oh))
    res = Dict()
    for (i, n) in enumerate(names)
        res[n] = Array(getRealData(oh, i-1))
    end
    return res
end

end # module


text:
# `rebuild` and `dims` are key methods to add for a new type

"""
    rebuild(x, args...)
    rebuild(x; kw...)

Rebuild an object struct with updated field values.

`x` can be a `AbstractDimArray`, a `Dimension`, `LookupArray` or other custom types.

This is an abstraction that alows inbuilt and custom types to be rebuilt
to update their fields, as most objects in DimensionalData.jl are immutable.

The arguments version can be concise but depends on a fixed order defined for some
DimensionalData objects. It should be defined based on the object type in DimensionalData,
adding the fields specific to your object.

The keyword version ignores order, and is mostly automated
using `ConstructionBase.setproperties`. It should only be defined if your object has
missing fields or fields with different names to DimensionalData objects.

The arguments required are defined for the abstract type that has a `rebuild` method.
"""
function rebuild end

"""
    dims(x, [dims::Tuple]) => Tuple{Vararg{<:Dimension}}
    dims(x, dim) => Dimension

Return a tuple of `Dimension`s for an object, in the order that matches
the axes or columns of the underlying data.

`dims` can be `Dimension`, `Dimension` types, or `Symbols` for `Dim{Symbol}`.

The default is to return `nothing`.
"""
function dims end

"""
    refdims(x, [dims::Tuple]) => Tuple{Vararg{<:Dimension}}
    refdims(x, dim) => Dimension

Reference dimensions for an array that is a slice or view of another
array with more dimensions.

`slicedims(a, dims)` returns a tuple containing the current new dimensions
and the new reference dimensions. Refdims can be stored in a field or disgarded,
as it is mostly to give context to plots. Ignoring refdims will simply leave some
captions empty.

The default is to return an empty `Tuple` `()`.
"""
function refdims end
refdims(x, lookup) = dims(refdims(x), lookup)

"""
    val(x)
    val(dims::Tuple) => Tuple

Return the contained value of a wrapper object.

`dims` can be `Dimension`, `Dimension` types, or `Symbols` for `Dim{Symbol}`.

Objects that don't define a `val` method are returned unaltered.
"""
function val end

# Methods to retreive Object/Dimension/LookupArray properties
#
# These work on AbstractDimStack, AbstractDimArray, Dimension
# LookupArray, and tuples of Dimension/LookupArray. A `dims` argument
# can be supplied to select a subset of dimensions or a single
# Dimension.

"""
    index(x) => Tuple{Vararg{<:AbstractArray}}
    index(x, dims::Tuple) => Tuple{Vararg{<:AbstractArray}}
    index(dims::Tuple) => Tuple{Vararg{<:AbstractArray}}}
    index(x, dim) => AbstractArray
    index(dim::Dimension) => AbstractArray

Return the contained index of a `Dimension`.

Only valid when a `Dimension` contains an `AbstractArray`
or a Val tuple like `Val{(:a, :b)}()`. The `Val` is unwrapped
to return just the `Tuple`

`dims` can be a `Dimension`, or a tuple of `Dimension`.
"""
function index end

"""
    lookup(x::Dimension) => LookupArray
    lookup(x, [dims::Tuple]) => Tuple{Vararg{<:LookupArray}}
    lookup(x::Tuple) => Tuple{Vararg{<:LookupArray}}
    lookup(x, dim) => LookupArray

Returns the [`LookupArray`](@ref) of a dimension. This dictates
properties of the dimension such as array axis and index order,
and sampling properties.

`dims` can be a `Dimension`, a dimension type, or a tuple of either.
"""
function lookup end

"""
    metadata(x) => (object metadata)
    metadata(x, dims::Tuple)  => Tuple (Dimension metadata)
    metadata(xs::Tuple) => Tuple

Returns the metadata for an object or for the specified dimension(s)

Second argument `dims` can be `Dimension`s, `Dimension` types,
or `Symbols` for `Dim{Symbol}`.
"""
function metadata end

"""
    name(x) => String
    name(xs:Tuple) => NTuple{N,String}
    name(x, dims::Tuple) => NTuple{N,String}
    name(x, dim) => String

Get the name of an array or Dimension, or a tuple of of either.

Second argument `dims` can be `Dimension`s, `Dimension` types,
or `Symbols` for `Dim{Symbol}`.
"""
function name end

"""
    units(x) => Union{Nothing,Any}
    units(xs:Tuple) => Tuple
    unit(A::AbstractDimArray, dims::Tuple) => Tuple
    unit(A::AbstractDimArray, dim) => Union{Nothing,Any}

Get the units of an array or `Dimension`, or a tuple of of either.

Units do not have a set field, and may or may not be included in `metadata`.
This method is to facilitate use in labels and plots when units are available,
not a guarantee that they will be. If not available, `nothing` is returned.

Second argument `dims` can be `Dimension`s, `Dimension` types,
or `Symbols` for `Dim{Symbol}`.
"""
function units end

"""
    label(x) => String
    label(x, dims::Tuple) => NTuple{N,String}
    label(x, dim) => String
    label(xs::Tuple) => NTuple{N,String}

Get a plot label for data or a dimension. This will include the name and units
if they exist, and anything else that should be shown on a plot.

Second argument `dims` can be `Dimension`s, `Dimension` types,
or `Symbols` for `Dim{Symbol}`.
"""
function label end

"""
    bounds(xs, [dims::Tuple]) => Tuple{Vararg{<:Tuple{T,T}}}
    bounds(xs::Tuple) => Tuple{Vararg{<:Tuple{T,T}}}
    bounds(x, dim) => Tuple{T,T}
    bounds(dim::Union{Dimension,LookupArray}) => Tuple{T,T}

Return the bounds of all dimensions of an object, of a specific dimension,
or of a tuple of dimensions.

`dims` can be a `Dimension`, a dimension type, or a tuple of either.
"""
function bounds end

"""
    order(x, [dims::Tuple]) => Tuple
    order(xs::Tuple) => Tuple
    order(x::Union{Dimension,LookupArray}) => Order

Return the `Ordering` of the dimension index for each dimension:
`ForwardOrdered`, `ReverseOrdered`, or [`Unordered`](@ref)

Second argument `dims` can be `Dimension`s, `Dimension` types,
or `Symbols` for `Dim{Symbol}`.
"""
function order end

"""
    sampling(x, [dims::Tuple]) => Tuple
    sampling(x, dim) => Sampling
    sampling(xs::Tuple) => Tuple{Vararg{<:Sampling}}
    sampling(x:Union{Dimension,LookupArray}) => Sampling

Return the [`Sampling`](@ref) for each dimension.

Second argument `dims` can be `Dimension`s, `Dimension` types,
or `Symbols` for `Dim{Symbol}`.
"""
function sampling end

"""
    span(x, [dims::Tuple]) => Tuple
    span(x, dim) => Span
    span(xs::Tuple) => Tuple{Vararg{<:Span,N}}
    span(x::Union{Dimension,LookupArray}) => Span

Return the [`Span`](@ref) for each dimension.

Second argument `dims` can be `Dimension`s, `Dimension` types,
or `Symbols` for `Dim{Symbol}`.
"""
function span end

"""
    locus(x, [dims::Tuple]) => Tuple
    locus(x, dim) => Locus
    locus(xs::Tuple) => Tuple{Vararg{<:Locus,N}}
    locus(x::Union{Dimension,LookupArray}) => Locus

Return the [`Locus`](@ref) for each dimension.

Second argument `dims` can be `Dimension`s, `Dimension` types,
or `Symbols` for `Dim{Symbol}`.
"""
function locus end

"""
    hasselection(x, selector) => Bool
    hasselection(x, selectors::Tuple) => Bool

Check if indexing into x with `selectors` can be performed, where
x is some object with a `dims` method, and `selectors` is a `Selector`
or `Dimension` or a tuple of either.
"""
function hasselection end


text:
using Revise, Jive
watch(@__DIR__, sources=[]) do path
    @info :changed path
    revise()
    runtests(@__DIR__, skip=["revise.jl"])
end
# Jive.stop(watch)


text:
# Autogenerated wrapper script for SBML_jll for aarch64-linux-gnu-cxx03
export libsbml

using XML2_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("SBML")
JLLWrappers.@declare_library_product(libsbml, "libsbml.so.5")
function __init__()
    JLLWrappers.@generate_init_header(XML2_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libsbml,
        "lib/libsbml.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export
    dnnAlgorithmConvolutionGemm,
    dnnAlgorithmConvolutionDirect,
    dnnAlgorithmConvolutionFFT,
    dnnBorderZeros, dnnBorderExtrapolation # dnnBorder_t

function convolution_desc{N}(x::Array, w::Array,
    padding::NTuple{N,Int}, strides::NTuple{N,Int};
    algo=dnnAlgorithmConvolutionDirect, border=dnnBorderZeros)

    csize_x = Csize_t[size(x)...]
    dims = ntuple(i -> (size(x,i)+2padding[i]-size(w,i)) ÷ strides[i] + 1, N)
    csize_y = Csize_t[dims..., size(w,ndims(w)), size(x,ndims(x))]
    csize_w = Csize_t[size(w)...]
    cstrides = Csize_t[strides...]
    cpadding = Cint[padding...]
    (algo, ndims(x), csize_x, csize_y, csize_w, cstrides, cpadding, border)
end

function convolution{T}(x::Array{T}, w::Array{T}, padding, strides)
    desc = convolution_desc(x, w, padding, strides)
    y = Array{T}(desc[4]...)
    p = Ptr{Cvoid}[0]
    attr = dnnPrimitiveAttributesCreate()
    dnnConvolutionCreateForward_F32(p, attr, desc...)
    p_conv = p[1]

    resources = Ptr{Cvoid}[pointer(x),pointer(y),pointer(w)]
    dnnExecute(p_conv, resources)

    dnnPrimitiveAttributesDestroy(attr)
    dnnDelete(p_conv)
    y
end

function ∇convolution_data{T}(x::Array{T}, w::Array{T}, dy::Array{T}, padding, strides)
    p = Ptr{Cvoid}[0]
    attr = dnnPrimitiveAttributesCreate()
    desc = convolution_desc(x, w, padding, strides)
    dnnConvolutionCreateBackwardData_F32(p, attr, desc...)
    p_conv = p[1]

    dx = similar(x)
    resources = Ptr{Cvoid}[0,0,pointer(w),0,pointer(dx),0,0,pointer(dy)]
    dnnExecute(p_conv, resources)

    dnnPrimitiveAttributesDestroy(attr)
    dnnDelete(p_conv)
    dx
end

function ∇convolution_filter{T}(x::Array{T}, w::Array{T}, dy::Array{T}, padding, strides)
    p = Ptr{Cvoid}[0]
    attr = dnnPrimitiveAttributesCreate()
    desc = convolution_desc(x, w, padding, strides)
    dnnConvolutionCreateBackwardFilter_F32(p, attr, desc...)
    p_conv = p[1]

    dw = similar(w)
    resources = Ptr{Cvoid}[pointer(x),0,0,0,0,pointer(dw),0,pointer(dy)]
    dnnExecute(p_conv, resources)

    dnnPrimitiveAttributesDestroy(attr)
    dnnDelete(p_conv)
    dw
end

function ∇convolution_bias{T}(x::Array{T}, w::Array{T}, b::Array{T}, dy::Array{T}, padding, strides)
    p = Ptr{Cvoid}[0]
    attr = dnnPrimitiveAttributesCreate()
    algorithm,dimension,dstSize
    desc = convolution_desc(x, w, padding, strides)
    dnnConvolutionCreateBackwardBias_F32(p, attr, desc...)
    p_conv = p[1]

    db = similar(b)
    resources = Ptr{Cvoid}[0,0,0,0,0,0,pointer(db),pointer(dy)]
    dnnExecute(p_conv, resources)

    dnnPrimitiveAttributesDestroy(attr)
    dnnDelete(p_conv)
    db
end


text:
module VulkanGen

using Transducers, DataStructures, JuliaFormatter, Parameters, EzXML, AbstractTrees, LightGraphs
using VulkanCore.vk
using VulkanCore:vk
using VulkanCore

include("spec.jl")

include("utils.jl")
include("naming_conventions.jl")
include("signatures.jl")
include("types.jl")
include("codegen.jl")
include("vulkan.jl")
include("parse.jl")

api = API("", OrderedDict.([(), (), (), ()])..., nothing) # this variable is replaced by the api provided in the function wrap

include("dependency_resolution.jl")
include("wrap_api.jl")

export vk,
CamelCaseLower,
CamelCaseUpper,
FDefinition,
SDefinition,
Signature,
SnakeCaseLower,
SnakeCaseUpper,
Statement,
VulkanCore,
argnames,
argtypes,
detect_convention,
enforce_convention,
generate,
prefix_vk,
remove_parts,
remove_prefix,
statements,
vk_prefix,
@sym_str,
remove_prefix,
wrap,
API,
parse_api,
CDefinition,
EDefinition,
PositionalArgument,
KeywordArgument,
StructWrapper,
ConstWrapper,
FuncWrapper,
patterns,
vulkan_to_julia,
nc_convert,
Generated

end # module VulkanGen


text:
# Option Pricing Formulas
using DataFrames
using Distributions
using Roots
using Dates
using RobustModels
using MixedModels
using GLM
using ShiftedArrays

using BLPData

"""
    bsopm_call(S, K, r, σ, t; q=0.0, u=0.0)
    Black-Scholes Pricing Model for Call Options

    returns call value, delta and vega as a named tuple
"""
function bsopm_call(S, K, r, σ, t; q=0.0, u=0.0)
    n = Normal()
    d₁ = (log(S/K) + ((r-q-u) + 0.5 * σ^2) * t) /(σ * √t)
    d₂ = d₁ - σ * √t
    value = S * exp(-(q+u)*t) * cdf(n, d₁) - K * exp(-r*t) * cdf(n, d₂)

    δ = exp(-(q+u)*t) * cdf(n, d₁)
    vega = S * exp(-(q+u)*t) * √t * pdf(n, d₁)

    return (value=value, δ=δ, vega=vega)
end  # function bsomp_call

function bsopm_impvol(opt_val, S, K, r, t; q=0.0, u=0.0)

    f(x) = opt_val - bsopm_call(S, K, r, x, t, q=q, u=u).value
    r = try
        find_zero(f, (0.0, 500.0)) # , Bisection())
    catch e
        println(e)
        return 0.0
    end
    return r
end  # function bsopm_impvol

"""
    calc_term(valuedate::Date, id::String, ipodate::Date,
        enddate::Date;
        dadate::Union{Missing,Date}=missing,
        mergerdate::Union{Missing,Date}=missing,
        closingdate::Union{Missing,Date}=missing,
        da_to_merger_days::Int64=120,
        merger_to_effective_days::Int64=30,
        ann_factor::Int64=252)

calculates exercise date, maturity date and option
term in years.
"""
function calc_term(valuedate::Date, id::String, ipodate::Date;
    enddate::Union{Missing,Date}=missing,
    dadate::Union{Missing,Date}=missing,
    mergerdate::Union{Missing,Date}=missing,
    closingdate::Union{Missing,Date}=missing,
    da_to_merger_days::Int64=120,
    merger_to_effective_days::Int64=30,
    ann_factor::Int64=252)

    spac_term = (ismissing(enddate) ? closingdate : enddate) - ipodate

    if ismissing(dadate) && ismissing(mergerdate)
        maturity_date = ipodate + Month(60)
        first_exercise_date = ipodate + Day(spac_term)
    elseif ismissing(mergerdate)
        maturity_date = dadate + Day(da_to_merger_days) + Month(60)
        first_exercise_date = dadate + Day(da_to_merger_days) +
            Day(merger_to_effective_days)
    else
        maturity_date = mergerdate + Month(60)
        first_exercise_date =  max(ipodate + Month(12),
            mergerdate + Day(merger_to_effective_days))
    end

    t = max(1, (maturity_date - valuedate).value)/365.25
    Nstep = max(30, convert(Int, round(t * ann_factor)))

    e = max(0, (first_exercise_date - valuedate).value)/365.25
    Estep = convert(Int, round(e * ann_factor))

    return (id=id, date=valuedate, t=t, Nstep=Nstep, Estep=Estep)
end

"""
    calc_hr(S, K, r, σ, t, a, b, c)

calculates warrant HR based on HW market model
"""
function calc_hr(S, K, r, σ, t, a, b, c)
    value, δ, vega =  bsopm_call(S, K, r, σ, t)
    δₘᵥ = clamp(δ + (vega / (S * sqrt(t))) * (a + b * δ + c * δ^2), 0, 1)

    return (S=S, value=value, δ=δ, vega=vega, δₘᵥ=δₘᵥ)
end  # function calc_hr

"""
    fortify_mv_df(dyn_df::DataFrame; K::Float64=11.5, r::Float64=0.0008)

generates price data and computes min-variance deltas for a stock/warrant pair.
"""
function fortify_mv_df!(dyn_df::DataFrame; K::Float64=11.5, r::Float64=0.0008)

    # create phase variables

    # transform(z_u, [:date, :dadate, :closingdate] => ByRow( (x, y, z) ->
    #     !(ismissing(y) || ismissing(z)) ? ifelse(x >= y && x < z, 2,  3) : 1) => :phase)

    dyn_df[!, :ls] = lag(dyn_df.s)
    dyn_df[!, :ds] = dyn_df.s .- dyn_df.ls

    dyn_df[!, :lw] = lag(dyn_df.w)
    dyn_df[!, :dw] = dyn_df.w .- dyn_df.lw

    z = DataFrame(map(row -> calc_term(row.date, row.id, row.ipodate, enddate=row.enddate,
            dadate=row.dadate, mergerdate=row.mergerdate, closingdate=row.closingdate),
            eachrow(dyn_df)))

    leftjoin!(dyn_df, z, on=[:id, :date])
    select!(dyn_df, Not([:id, :ipodate, :dadate, :closingdate, :enddate, :mergerdate, :Nstep, :Estep]))

    # fortify data with bsopm delta and vega
    transform!(dyn_df, [:w, :s, :t] => ByRow( (x, y, z) -> bsopm_impvol(x, y, K, r, z)) => :σ)
    transform!(dyn_df, [:s, :σ, :t] => ByRow( (a, b, c) -> !isnothing(b) ?
        bsopm_call(a, K, r, b, c) : (value=missing, δ=missing, vega=missing)) => [:value, :δ, :vega])

    return dyn_df
end  # function fortify_mv_df!

"""
    compute_mv_delta(dyn_df::DataFrame;
        window::Union{Missing,Int64}=missing, robust::Bool=true)

computes min-variance deltas for a stock/warrant pair.
"""
function compute_mv_delta(dyn_df::DataFrame;
        window::Union{Missing,Int64}=missing, robust::Bool=true)

    if ismissing(window)
        window = nrow(dyn_df)
    end

    sub_df = last(dyn_df, window)
    Y = sub_df.dw .- sub_df.δ .* sub_df.ds
    μ = (sub_df.vega ./ sqrt.(sub_df.t)) .* (sub_df.ds./sub_df.s)
    x₁ = μ
    x₂ = μ .* sub_df.δ
    x₃ = μ .* sub_df.δ.^2

    if robust
        hw_reg = rlm(@formula(Y ~ -1 + x₁ + x₂ + x₃),
            DataFrame(Y=Y, x₁=x₁, x₂=x₂, x₃=x₃),
            MEstimator{HuberLoss}(); method=:chol,
            initial_scale=15.0)
    else
        hw_reg = lm(@formula(Y ~ -1 + x₁ + x₂ + x₃),
            DataFrame(Y=Y, x₁=x₁, x₂=x₂, x₃=x₃))
    end

    a, b, c = coef(hw_reg)
    sub_df[!, :δₘᵥ] = max.(0.0, sub_df.δ .+ (sub_df.vega ./
        (sub_df.s .* sqrt.(sub_df.t))) .* (a .+ b .* sub_df.δ + c .* sub_df.δ.^2))

    return (sub_df=sub_df, hw_reg=hw_reg, params=(a, b, c))

end  # function compute_mv_delta

"""
    lme_mv_delta(dyn_df::DataFrame;
        start_date::Union{Missing,Date}=missing,
        end_date::Union{Missing,Date}=missing)

computes min-variance deltas for a stock/warrant pair.
"""
function lme_mv_delta(dyn_df::DataFrame;
        start_date::Union{Missing,Date}=missing,
        end_date::Union{Missing,Date}=missing)

    if !ismissing(end_date)
        sub_df = subset(dyn_df, :date => x -> x .<= end_date, skipmissing=true)
    elseif !ismissing(start_date)
        sub_df = subset(dyn_df, :date => x -> x .>= start_date, skipmissing=true)
    elseif !(ismissing(start_date) && ismissing(end_date))
        sub_df = subset(dyn_df, :date => x -> x .>= start_date && x .<= end_date,
            skipmissing=true)
    else
        sub_df = dyn_df
    end

    Y = sub_df.dw .- sub_df.δ .* sub_df.ds
    μ = (sub_df.vega ./ sqrt.(sub_df.t)) .* (sub_df.ds./sub_df.s)
    x₁ = μ
    x₂ = μ .* sub_df.δ
    x₃ = μ .* sub_df.δ.^2

    estim_df = DataFrame(Y=Y, x₁=x₁, x₂=x₂, x₃=x₃, symbol=sub_df.symbol)

    fm = @formula(Y ~ 0 + x₁ + x₂ + x₃ + (0 + x₁ + x₂ + x₃|symbol))
    fm1 = fit(MixedModel, fm, estim_df)

    a, b, c = fixef(fm1)
    fixef_df = DataFrame(symbol="AVG", x₁=a, x₂=b, x₃=c)
    ranef_df = DataFrame(only(raneftables(fm1)))
    param_df = reduce(vcat, [fixef_df, ranef_df])
    param_df[!, :obsdate] .= maximum(dyn_df.date)

    dyn_df[!, :δₘᵥ] = max.(0.0, dyn_df.δ .+ (dyn_df.vega ./
        (dyn_df.s .* sqrt.(dyn_df.t))) .* (a .+ b .* dyn_df.δ + c .* dyn_df.δ.^2))

    return (estim=estim_df, fm1=fm1, param_df=param_df)
end

max_stars_repo_name: GodolphinCapitalManagement/Iex

max_issues_repo_name: GodolphinCapitalManagement/Iex

max_forks_repo_name: GodolphinCapitalManagement/Iex


text:
using BenchmarkTools, SparseArrays, LinearAlgebra, LinearMaps
import ExpmV, Expokit

#BenchmarkTools.DEFAULT_PARAMETERS.seconds = 120

const p = 1e-3
const m = 8
const dt = 0.5 + 0.1*rand() # time

SUITE = BenchmarkGroup()

SUITE["Expm"] = BenchmarkGroup()
SUITE["ExpmV"] = BenchmarkGroup()
SUITE["Expokit"] = BenchmarkGroup()
SUITE["ExpmV_LinMap"] = BenchmarkGroup()
SUITE["Expokit_LinMap"] = BenchmarkGroup()
SUITE["ExpmV_Precomp"] = BenchmarkGroup()
SUITE["Expokit_Precomp"] = BenchmarkGroup()
SUITE["Expokit_LinMap_Precomp"] = BenchmarkGroup()

dimensions = [d for d in 2 .^ (5:10)]

for d in dimensions
    A = sprandn(ComplexF64,d,d,m/d) + m*I
    L = LinearMap(A)
    b = normalize(randn(ComplexF64,d))

    SUITE["Expm"][d]            = @benchmarkable exp(dt_A)*$b setup=(dt_A = Matrix($dt*$A))
    SUITE["ExpmV"][d]           = @benchmarkable ExpmV.expmv($dt,$A,$b;shift=true)
    SUITE["Expokit"][d]         = @benchmarkable Expokit.expmv($dt,$A,$b)
    SUITE["ExpmV_LinMap"][d]    = @benchmarkable ExpmV.expmv($dt,$L,$b;shift=false)
    SUITE["Expokit_LinMap"][d]  = @benchmarkable Expokit.expmv($dt,$L,$b)

    f, b1, b2 = similar(b), similar(b), similar(b)
    M, _ = ExpmV.select_taylor_degree(A,1)
    SUITE["ExpmV_Precomp"][d]   = @benchmarkable ExpmV.expmv!($f,$dt,$A,$b,$M,$norm,$opnorm,$b1,$b2)
    SUITE["Expokit_Precomp"][d] = @benchmarkable Expokit.expmv!($f,$dt,$A,$b)
    SUITE["Expokit_LinMap_Precomp"][d] = @benchmarkable Expokit.expmv!($f,$dt,$A,$b)
end

results = run(SUITE, verbose = true, seconds = 5)

for d in dimensions
    results_d = results[@tagged d]
    names = collect(keys(results_d))
    len = maximum(length, names)
    sort!(names; by = k -> minimum(results[k][d]).time)

    @show d
    for name in names
        print(rpad(name * ": ", len+2))
        show(results[name][d])
        println("")
    end
    println("")
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/MortarContact2DAD.jl/blob/master/LICENSE

using MortarContact2DAD, Test

# Matching mesh, undeformed, gap = 1.0 between surfaces

X = Dict(
    1 => [0.0, 0.0],
    2 => [0.0, 2.0],
    3 => [1.0, 2.0],
    4 => [1.0, 0.0])

u = Dict(
    1 => [0.0, 0.0],
    2 => [0.0, 0.0],
    3 => [0.0, 0.0],
    4 => [0.0, 0.0])

slave = Element(Seg2, [1, 2])
master = Element(Seg2, [3, 4])
update!([slave, master], "geometry", X)
update!([slave, master], "displacement", u)
problem = Problem(Contact2DAD, "test problem", 2, "displacement")
push!(problem.properties.always_in_contact, 1, 2)
problem.properties.rotate_normals = true
add_slave_elements!(problem, [slave])
add_master_elements!(problem, [master])
problem.assembly.u = zeros(8)
problem.assembly.la = zeros(8)
assemble!(problem, 0.0)

n = slave("normal", 0.0)
@test isapprox(n[1], n[2])
@test isapprox(n[1], [1.0, 0.0])

C1 = Matrix(sparse(problem.assembly.C1, 4, 8))
C2 = Matrix(sparse(problem.assembly.C2, 4, 8))
K = Matrix(sparse(problem.assembly.K, 4, 8))
D = Matrix(sparse(problem.assembly.D, 4, 8))
f = Vector(sparse(problem.assembly.f, 4, 1)[:])
g = Vector(sparse(problem.assembly.g, 4, 1)[:])

C1_expected = [
 1.0  0.0  0.0  0.0   0.0   0.0  -1.0   0.0
 0.0  1.0  0.0  0.0   0.0   0.0   0.0  -1.0
 0.0  0.0  1.0  0.0  -1.0   0.0   0.0   0.0
 0.0  0.0  0.0  1.0   0.0  -1.0   0.0   0.0]
C2_expected = [
 1.0   0.5  0.0 -0.5   0.0  0.0  -1.0  0.0
 0.0   0.0  0.0  0.0   0.0  0.0   0.0  0.0
 0.0   0.5  1.0 -0.5  -1.0  0.0   0.0  0.0
 0.0   0.0  0.0  0.0   0.0  0.0   0.0  0.0]
D_expected = [
 0.0   0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  -1.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  -1.0  0.0  0.0  0.0  0.0]
@test isapprox(C1, C1_expected)
@test isapprox(K, zeros(4, 8))
@test isapprox(f, zeros(4))
@test isapprox(g, [1.0, 0.0, 1.0, 0.0])
@test isapprox(C2, C2_expected)
@test isapprox(D, D_expected)