text:
using Revise
using LinearAlgebra, StaticArrays, GeometryPrimitives, Zygote, ForwardDiff #ChainRules, ReverseDiff
##
Zygote.@adjoint (T::Type{<:SArray})(xs::Number...) = T(xs...), dv -> (nothing, dv...)
Zygote.@adjoint (T::Type{<:SArray})(x::AbstractArray) = T(x), dv -> (nothing, dv)
Zygote.@adjoint (T::Type{<:SMatrix})(xs::Number...) = T(xs...), dv -> (nothing, dv...)
Zygote.@adjoint (T::Type{<:SMatrix})(x::AbstractMatrix) = T(x), dv -> (nothing, dv)
Zygote.@adjoint (T::Type{<:SVector})(xs::Number...) = T(xs...), dv -> (nothing, dv...)
Zygote.@adjoint (T::Type{<:SVector})(x::AbstractVector) = T(x), dv -> (nothing, dv)

Zygote.@adjoint enumerate(xs) = enumerate(xs), diys -> (map(last, diys),)
_ndims(::Base.HasShape{d}) where {d} = d
_ndims(x) = Base.IteratorSize(x) isa Base.HasShape ? _ndims(Base.IteratorSize(x)) : 1
Zygote.@adjoint function Iterators.product(xs...)
                    d = 1
                    Iterators.product(xs...), dy -> ntuple(length(xs)) do n
                        nd = _ndims(xs[n])
                        dims = ntuple(i -> i<d ? i : i+nd, ndims(dy)-nd)
                        d += nd
                        func = sum(y->y[n], dy; dims=dims)
                        ax = axes(xs[n])
                        reshape(func, ax)
                    end
                end

Zygote.refresh()

##
function ε_tensor(n::Float64)
    n² = n^2
    ε =     [	n²      0. 	    0.
                0. 	    n² 	    0.
                0. 	    0. 	    n²  ]
end

# volfrac(vxl::NTuple{2,SVector{2,<:Number}}, nout::SVector{2,<:Real}, r₀::SVector{2,<:Real}) =
#     volfrac((SVector(vxl[N][1],vxl[N][2],0), SVector(vxl[P][1],vxl[P][2],1)),
#             SVector(nout[1], nout[2], 0),
#             SVector(r₀[1], r₀[2], 0))

function surfpt_nearby2(x::Vector{Float64}, s::Sphere{2})
    nout = x==s.c ? SVector(1.0,0.0) : # nout = e₁ for x == s.c
                    normalize(x-s.c)
    return s.c+s.r*nout, nout
end

f_onbnd(bin,absdin) =  (b = Zygote.@ignore bin; absd = Zygote.@ignore absdin; abs.(b.r.-absd) .≤ Base.rtoldefault(Float64) .* b.r)  # basically b.r .≈ absd but faster
f_isout(b,absd) =  (isout = Zygote.@ignore ((b.r.<absd) .| f_onbnd(b,absd) ); isout)
# isout(bin) =  (b = Zygote.@ignore bin; (b.r.<absd) .| (abs.(b.r.-absd) .≤ Base.rtoldefault(Float64) .* b.r))
f_signs(d) =  (signs = Zygote.@ignore (copysign.(1.0,d)'); signs)
function surfpt_nearby2(x, b::Box{2})
    ax = inv(b.p)
    n0 = b.p ./  [ sqrt(b.p[1,1]^2 + b.p[1,2]^2) sqrt(b.p[2,1]^2 + b.p[2,2]^2)  ]
    d = Array(b.p * (x - b.c))
    cosθ = diag(n0*ax)

    n = n0 .* f_signs(d)
    absd = abs.(d)
    ∆ = (b.r .- absd) .* cosθ
    # onbnd = abs.(dropgrad(b.r).-dropgrad(absd)) .≤ Base.rtoldefault(Float64) .* dropgrad(b.r)  # basically b.r .≈ absd but faster
    # isout = (dropgrad(b.r).<dropgrad(absd)) .| dropgrad(onbnd)
    onbnd = f_onbnd(b,absd)
    isout = f_isout(b,absd)
    projbnd =  all(.!isout .| onbnd)
    onbnd_float =  Float64.(onbnd)
    isout_float =  Float64.(isout)
    # if ( ( abs(b.r[1]-absd[1]) ≤ Base.rtoldefault(Float64) * b.r[1] ) | (b.r[1]<absd[1]) ) | ( ( abs(b.r[2]-absd[2]) ≤ Base.rtoldefault(Float64) * b.r[2] ) | ( b.r[2] < absd[2] ) )
    if count(isout) == 0
        # = not(isout[1]) &  not(isout[2]) case, point is inside box
        l∆x, i = findmin(∆)  # find closest face
        nout = n[i,:]
        ∆x = l∆x * nout
    else
        ∆x = n' * (∆ .* isout_float)
        # = isout[1] | isout[2] case 1: at least one dimension outside or on boundry
        if all(.!isout .| onbnd)
            nout0 = n' * onbnd_float
        else
            nout0 = -∆x
        end
        nout = nout0 / norm(nout0)
    end

    return SVector{2,Float64}(x+∆x), SVector{2,Float64}(nout)
end

# Equation (4) of the paper by Kottke et al.
function τ_trans(ε)
    ε₁₁, ε₂₁, ε₃₁, ε₁₂, ε₂₂, ε₃₂, ε₁₃, ε₂₃, ε₃₃ = ε
    return SMatrix{3,3,Float64,9}(
        -1/ε₁₁, ε₂₁/ε₁₁, ε₃₁/ε₁₁,
        ε₁₂/ε₁₁, ε₂₂ - ε₂₁*ε₁₂/ε₁₁, ε₃₂ - ε₃₁*ε₁₂/ε₁₁,
        ε₁₃/ε₁₁, ε₂₃ - ε₂₁*ε₁₃/ε₁₁, ε₃₃ - ε₃₁*ε₁₃/ε₁₁
    )
end

# Equation (23) of the paper by Kottke et al.
function τ⁻¹_trans(τ)
    τ₁₁, τ₂₁, τ₃₁, τ₁₂, τ₂₂, τ₃₂, τ₁₃, τ₂₃, τ₃₃ = τ
    return SMatrix{3,3,Float64,9}(
        -1/τ₁₁, -τ₂₁/τ₁₁, -τ₃₁/τ₁₁,
        -τ₁₂/τ₁₁, τ₂₂ - τ₂₁*τ₁₂/τ₁₁, τ₃₂ - τ₃₁*τ₁₂/τ₁₁,
        -τ₁₃/τ₁₁, τ₂₃ - τ₂₁*τ₁₃/τ₁₁, τ₃₃ - τ₃₁*τ₁₃/τ₁₁
    )
end

# function kottke_avg_param(param1::SMat3Complex, param2::SMat3Complex, n12::SVec3Float, rvol1::Real)
function kottke_avg_param(param1, param2, n12, rvol1)
    n = n12 / norm(n12) #sqrt(sum2(abs2,n12))

    # Pick a vector that is not along n.
    if any(n .== 0)
    	htemp1 = (n .== 0)
    else
    	htemp1 = SVector(1., 0. , 0.)
    end

    # Create two vectors that are normal to n and normal to each other.
    htemp2 = n × htemp1
    h = htemp2 / norm(htemp2) #sqrt(sum2(abs2,htemp2))
    vtemp = n × h
    v = vtemp / norm(vtemp) #sqrt(sum2(abs2,vtemp))
    # Create a local Cartesian coordinate system.
    S = [n h v]  # unitary

    τ1 = τ_trans(transpose(S) * param1 * S)  # express param1 in S coordinates, and apply τ transform
    τ2 = τ_trans(transpose(S) * param2 * S)  # express param2 in S coordinates, and apply τ transform

    τavg = τ1 .* rvol1 + τ2 .* (1-rvol1)  # volume-weighted average

    return S * τ⁻¹_trans(τavg) * transpose(S)  # apply τ⁻¹ and transform back to global coordinates
end

fign(x) = (y = Zygote.@ignore x; x * y)
fign2(x) = (y = x; x * y)
fign(3.0)
fign2(3.0)
Zygote.gradient(fign,3.0)
Zygote.gradient(fign2,3.0)

function ftest2(a,b,c)
    Zygote.ignore() do
        d = a + 2
        println("d: $d")
    end
    a * b * c * d
end

ftest2(1.,2.,3.)

function foo(p1,p2,p3,p4,p5,p6,p7)
    s = Sphere(					# Instantiate N-D sphere, here N=2 (circle)
        [p1,p2],				# c: center
        p3,						# r: "radii" (half span of each axis)
        ε_tensor(p4),		    # data: any type, data associated with circle shape
        )
    ε_bg = ε_tensor(p5)
    # sum(abs2.(bounds(s)[1])) #*p1
    # bounds(s)[1][1] #*p1
    # r₀,nout = surfpt_nearby2([p6,p7],s)
    r₀,nout = surfpt_nearby2([p6,p7],s)
    vxl1 = SVector(1.3,2.4)
    vxl2 = SVector(-1.2,6.4)
    rvol = volfrac((vxl1,vxl2),nout,r₀)
    kottke_avg_param(s.data, ε_bg, [nout[1];nout[2];0.0], rvol)[1,2]
end

foo(0.13,0.12,1.3,2.24,3.5,0.34,0.41)
@show Zygote.gradient(foo,0.13,0.12,1.3,2.24,3.5,0.34,0.41)

##
@show p8
@show p9
@show p10
@show p11
@show [[1  2];  [3  4]]

ax = [      ax11     ax21
            ax12     ax22      ]
               # d: size of the box in axis directions (cols of ax)
SMatrix{2,2,Float64,4}(ax)

function boo(ax11_v,ax12_v,ax21_v,ax22_v,c1_v,c2_v,d1_v,d2_v,n1_v,n2_v,xx_v,yy_v,dx_v,dy_v)
    ## columns of ax are the axes of the box
    b_ax = [      ax11_v     ax21_v
                ax12_v     ax22_v      ]
    #ax1b,ax2b = GeometryPrimitives.normalize.(([p8,p9], [p10,p11]))
    b = Box(					# Instantiate N-D sphere, here N=2 (circle)
        [c1_v,c2_v],				# c: center
        [d1_v,d2_v],               # d: size of the box in axis directions (cols of ax)
        b_ax,						# r: "radii" (half span of each axis)
        ε_tensor(n1_v),		    # data: any type, data associated with circle shape
        )
    ε_bg = ε_tensor(n2_v)
    # sum(abs2.(bounds(s)[1])) #*p1
    # bounds(s)[1][1] #*p1
    # r₀,nout = surfpt_nearby2([p6,p7],s)
    # r₀,nout = surfpt_nearby2([xx_v; yy_v],b)
    r₀,nout = surfpt_nearby2([xx_v;yy_v],b)
    vxl2 = SVector{2}(xx_v+dx_v/2.,yy_v+dy_v/2.)
    vxl1 = SVector{2}(xx_v-dx_v/2.,yy_v-dy_v/2.)
    rvol = volfrac((vxl1,vxl2),nout,r₀)
    kottke_avg_param(b.data, ε_bg, [nout[1];nout[2];0.0], rvol)[1,1]
end

##
ax11    =   1.1
ax12    =   0.02
ax21    =   0.2
ax22    =   1.05
c1      =   0.1
c2      =   0.
d1      =   1.0
d2      =   1.0
n1      =   2.5
n2      =   3.5
xx      =   0.6
yy      =   0.
dx      =   0.03
dy      =   0.03
##
boo(ax11,ax12,ax21,ax22,c1,c2,d1,d2,n1,n2,xx,yy,dx,dy)
boox(q) = q -> boo(ax11,ax12,ax21,ax22,c1,c2,d1,d2,n1,n2,q,yy,dx,dy)
dboox(q) = q -> Zygote.gradient(boo,ax11,ax12,ax21,ax22,c1,c2,d1,d2,n1,n2,q,yy,dx,dy)[end-3]
using Plots: plot, plot!
x = 0.55:0.001:0.65
plot(x,boox(x))
plot(x,dboox(x))
# 3×3 Array{Float64,2}:
#  9.37196  -0.0           0.0
#  0.0      10.7091        2.66454e-15
#  0.0       2.66454e-15  10.7091
Zygote.gradient(boo,ax11,ax12,ax21,ax22,c1,c2,d1,d2,n1,n2,xx,yy,dx,dy)

ax1b,ax2b = GeometryPrimitives.normalize.(([1.,0.], [0.,1.]))

# SVector(
ntuple(k -> k==1 ? 1.0 : 0.0, 2)

function surfpt_nearby2(x::Vector{Float64}, s::Sphere{2})
    nout = x==s.c ? SVector(1.0,0.0) : # nout = e₁ for x == s.c
                    normalize(x-s.c)
    return s.c+s.r*nout, nout
end

p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12 = 0.13,0.12,1.3,2.24,3.5,0.34,0.41,1.0,0.2,0.0,1.0,1.4
ax1b,ax2b = GeometryPrimitives.normalize.(([p8,p9], [p10,p11]))
b = Box(					# Instantiate N-D sphere, here N=2 (circle)
    [p1,p2],				# c: center
    [p3,p12],
    [ax1b ax2b],						# r: "radii" (half span of each axis)
    ε_tensor(p4),		    # data: any type, data associated with circle shape
    )
for xx=0.4:0.01:5.2,yy=0.4:0.01:5.2

    x = SVector(xx,yy)
    ax = inv(b.p)
    ##
    n = (b.p ./ sqrt.(sum(abs2,b.p,dims=2)[:,1]))  # b.p normalized in row direction.
    cosθ = sum(ax.*n', dims=1)[1,:]  # equivalent to diag(n*ax)
    d = b.p * (x - b.c)
    n = n .* copysign.(1.0,d)'  # operation returns SMatrix (reason for leaving n untransposed)
    absd = abs.(d)
    onbnd = abs.(b.r.-absd) .≤ Base.rtoldefault(Float64) .* b.r  # basically b.r .≈ absd but faster
    isout = (b.r.<absd) .| onbnd
    ∆ = (b.r .- absd) .* cosθ  # entries can be negative
    if count(isout) == 0  # x strictly inside box; ∆ all positive
        case_str_old = "case1: in the box"
        l∆x, i = findmin(∆)  # find closest face
        nout = n[i,:]
        ∆x = l∆x * nout
        xout = x + ∆x
    else  # x outside box or on boundary in one or multiple directions
        ∆x = n' * (∆ .* isout)  # project out .!isout directions
        nout = all(.!isout .| onbnd) ? n'*onbnd : -∆x
        if all(.!isout .| onbnd)
            case_str =  "case2.1: proj onto boundry"
        else
            case_str = "case2.2: just outside "
        end
        nout = normalize(nout)
        xout = x + ∆x
    end

    nout_old = nout
    xout_old = xout

    n0 = b.p ./  [ sqrt(b.p[1,1]^2 + b.p[1,2]^2) sqrt(b.p[2,1]^2 + b.p[2,2]^2)  ]
    d = Array(b.p * (x - b.c))
    cosθ = diag(n0*ax)
    n = n0 .* copysign.(1.0,d)'
    absd = abs.(d)
    ∆ = (b.r .- absd) .* cosθ
    # onbnd = abs.(dropgrad(b).r.-dropgrad(absd)) .≤ Base.rtoldefault(Float64) .* dropgrad(b.r)  # basically b.r .≈ absd but faster
    # isout = (dropgrad(b.r).<dropgrad(absd)) .| dropgrad(onbnd)
    onbnd = Zygote.@ignore abs.(b.r.-absd) .≤ Base.rtoldefault(Float64) .* b.r  # basically b.r .≈ absd but faster
    isout = Zygote.@ignore (b.r.<absd) .| onbnd
    projbnd = Zygote.@ignore all(.!isout .| onbnd)
    onbnd_float = Zygote.@ignore Float64.(onbnd)
    isout_float = Zygote.@ignore Float64.(isout)
    # if ( ( abs(b.r[1]-absd[1]) ≤ Base.rtoldefault(Float64) * b.r[1] ) | (b.r[1]<absd[1]) ) | ( ( abs(b.r[2]-absd[2]) ≤ Base.rtoldefault(Float64) * b.r[2] ) | ( b.r[2] < absd[2] ) )
    if count(isout) == 0
        # = not(isout[1]) &  not(isout[2]) case, point is inside box
        case_str = "case1: in the box"
        l∆x, i = findmin(∆)  # find closest face
        nout = n[i,:]
        ∆x = l∆x * nout
    else
        ∆x = n' * (∆ .* isout_float)
        # = isout[1] | isout[2] case 1: at least one dimension outside or on boundry
        if all(.!isout .| onbnd)
            case_str =  "case2.1: proj onto boundry"
            nout0 = n' * onbnd_float
        else
            case_str = "case2.2: just outside "
            nout0 = -∆x
        end
        nout = nout0 / norm(nout0)
    end
    xout = x + ∆x

    if !((nout_old[1] ≈ nout[1]) & (nout_old[2] ≈ nout[2]) & (xout_old[1] ≈ xout[1]) & (xout_old[2] ≈ xout[2]))
        println("### outputs don't match at x,y: $xx,   $yy")
        println("# old code:")
        println(case_str_old)
        println("nout_old: $nout_old")
        println("xout_old: $xout_old")
        println("# new code:")
        println(case_str_old)
        println("nout: $nout")
        println("xout: $xout")
    end
    # @assert nout_old[1] ≈ nout[1]
    # @assert nout_old[2] ≈ nout[2]
    # @assert xout_old[1] ≈ xout[1]
    # @assert xout_old[2] ≈ xout[2]
end
##
xx1 = 1.3; yy1 = 3.4;
xx2 = 0.2; yy2 = 2.43;
println("### vectors:   x1,y1: $xx1,   $yy1")
println("###            x2,y2: $xx2,   $yy2")
@show v1, v2 =  normalize.(([xx1,yy1], [xx2,yy2]))

dot(v1,v2)

##

sum(abs2,b.p,dims=2)
s = Sphere(					# Instantiate N-D sphere, here N=2 (circle)
    [0.1,0.2],				# c: center
    0.3,						# r: "radii" (half span of each axis)
    ε_tensor(2.2),		    # data: any type, data associated with circle shape
    )
s.c

s = Sphere(					# Instantiate N-D sphere, here N=2 (circle)
    [0.1,0.1],				# c: center
    0.5,						# r: "radii" (half span of each axis)
    ε_tensor(2.2),		    # data: any type, data associated with circle shape
    )
surfpt_nearby([0.2,0.25],s)

surfpt_nearby([0.2,0.1],s)

sum(abs2.(bounds(s)[1]))


text:
using Distributed
const N_CORES = 7
rmprocs(workers())
addprocs(N_CORES)
@everywhere push!(LOAD_PATH,".")
@everywhere using SparseArrays, LinearAlgebra, SharedArrays
include("seif.jl")

seif()

ext: jl


text:
#!/usr/bin/env julia

# Test 3D real-to-complex FFTs.

using PencilFFTs
import MPI

using BenchmarkTools
using LinearAlgebra
using Printf
using Random
using Test

include("include/FourierOperations.jl")
using .FourierOperations

include("include/MPITools.jl")
using .MPITools

const DATA_DIMS_EVEN = (42, 24, 16)
const DATA_DIMS_ODD = DATA_DIMS_EVEN .- 1
const GEOMETRY = ((0.0, 4pi), (0.0, 2pi), (0.0, 2pi))

# Compute and compare ⟨ |u|² ⟩ in physical and spectral space.
function test_global_average(u, uF, plan::PencilFFTPlan,
                             gF::FourierGridIterator)
    comm = get_comm(plan)
    scale = scale_factor(plan)

    sum_u2_local = sqnorm(u)
    sum_uF2_local = sqnorm(uF, gF)

    Ngrid = prod(size_global(pencil(u)))

    avg_u2 = MPI.Allreduce(sum_u2_local, +, comm) / Ngrid

    # To get a physically meaningful quantity, squared values in Fourier space
    # must be normalised by `scale` (and their sum is normalised again by
    # `scale` if one wants the average).
    # Equivalently, uF should be normalised by `sqrt(scale)`.
    @test scale == Ngrid
    avg_uF2 = MPI.Allreduce(sum_uF2_local, +, comm) / (Ngrid * Float64(scale))

    @test avg_u2 ≈ avg_uF2

    nothing
end

# Squared 2-norm of a tuple of arrays using LinearAlgebra.norm.
norm2(x::Tuple) = sum(norm.(x).^2)

function micro_benchmarks(u, uF, gF_global::FourierGrid,
                          gF_local::FourierGridIterator)
    ωF = similar.(uF)

    BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1

    println("Micro-benchmarks:")

    @printf " - %-20s" "divergence global_view..."
    @btime divergence($uF, $gF_global)

    @printf " - %-20s" "divergence local..."
    @btime divergence($uF, $gF_local)

    @printf " - %-20s" "curl! global_view..."
    @btime curl!($ωF, $uF, $gF_global)

    @printf " - %-20s" "curl! local..."
    @btime curl!($ωF, $uF, $gF_local)

    # For these, a generic implementation is used (LinearAlgebra.generic_norm2).
    @printf " - %-20s" "norm2(u)..."
    @btime norm2($u)

    @printf " - %-20s" "norm2(uF)..."
    @btime norm2($uF)

    # These are much faster because parent(u) is a regular Array, and Julia
    # calls BLAS in this case (LinearAlgebra.BLAS.nrm2).
    @printf " - %-20s" "norm2(parent(u))..."
    @btime norm2($(parent.(u)))

    @printf " - %-20s" "norm2(parent(uF))..."
    @btime norm2($(parent.(uF)))

    # Interestingly, this is even faster than BLAS!
    @printf " - %-20s" "sqnorm(u)..."
    @btime sqnorm($u)

    @printf " - %-20s" "sqnorm(parent(u))..."
    @btime sqnorm($(parent.(u)))

    @printf " - %-20s" "sqnorm(uF)..."
    @btime sqnorm($uF, $gF_local)

    nothing
end

function init_random_field!(u::PencilArray{T}, rng) where {T <: Complex}
    fill!(u, zero(T))

    u_g = global_view(u)

    dims_global = size_global(pencil(u))
    ind_space = CartesianIndices(dims_global)
    ind_space_local = CartesianIndices(range_local(pencil(u), LogicalOrder()))
    @assert ndims(ind_space_local) == ndims(ind_space)

    scale = sqrt(2 * prod(dims_global))  # to get order-1 average values

    # Zero-based index of last element of r2c transform (which is set to zero)
    imax = dims_global[1] - 1

    # Loop over global dimensions, so that all processes generate the same
    # random numbers.
    for I in ind_space
        val = scale * randn(rng, T)

        I0 = Tuple(I) .- 1

        # First logical dimension, along which a r2c transform is applied.
        # If zero, Hermitian symmetry must be enforced.
        # (For this, zero-based indices are clearer!)
        i = I0[1]

        # Leave last element of r2c transform as zero.
        # Note: if I don't do this, the norms in physical and Fourier space
        # don't match... This is also the case if I set a real value to these
        # modes.
        i == imax && continue

        # We add in case a previous value was set by the code in the block
        # below (for Hermitian symmetry).
        if I ∈ ind_space_local
            u_g[I] += val
        end

        # If kx != 0, we're done.
        i == 0 || continue

        # Case kx == 0: account for Hermitian symmetry.
        #
        #    u(0, -ky, -kz) = conj(u(0, ky, kz))
        #
        # This also ensures that the zero mode is real.
        J0 = map((i, N) -> i == 0 ? 0 : N - i, I0, dims_global)
        J = CartesianIndex(J0 .+ 1)

        if J ∈ ind_space_local
            u_g[J] += conj(val)
        end
    end

    u
end

function test_rfft(size_in; benchmark=true)
    comm = MPI.COMM_WORLD
    Nproc = MPI.Comm_size(comm)
    rank = MPI.Comm_rank(comm)

    rank == 0 && @info "Input data size: $size_in"

    pdims_2d = let pdims = zeros(Int, 2)
        MPI.Dims_create!(Nproc, pdims)
        pdims[1], pdims[2]
    end

    # Test creating Pencil and PencilArray first, and creating plan afterwards.
    topo = MPITopology(comm, pdims_2d)
    pen = Pencil(topo, size_in, (2, 3))
    u1 = PencilArray{Float64}(undef, pen)

    plan = PencilFFTPlan(u1, Transforms.RFFT())

    # Allocate and initialise vector field in Fourier space.
    uF = allocate_output(plan, Val(3))
    rng = MersenneTwister(42)
    init_random_field!.(uF, (rng, ))

    u = (u1, similar(u1), allocate_input(plan))
    map(u) do v
        @test typeof(v) === typeof(u1)
        @test pencil(v) === pencil(u1)
        @test size(v) == size(u1)
        @test size_local(v) == size_local(u1)
    end
    ldiv!(u, plan, uF)

    gF_global = FourierGrid(GEOMETRY, size_in, permutation(uF))
    gF_local = LocalGridIterator(gF_global, uF)

    # Compare different methods for computing stuff in Fourier space.
    @testset "Fourier operations" begin
        test_global_average(u, uF, plan, gF_local)

        div_global = divergence(uF, gF_global)
        div_local = divergence(uF, gF_local)
        @test div_global ≈ div_local

        ωF_global = similar.(uF)
        ωF_local = similar.(uF)
        curl!(ωF_global, uF, gF_global)
        curl!(ωF_local, uF, gF_local)
        @test all(ωF_global .≈ ωF_local)
    end

    @test sqnorm(u) ≈ norm2(u)

    # These are not the same because `FourierOperations.sqnorm` takes Hermitian
    # symmetry into account, so the result can be roughly twice as large.
    @test 1 < sqnorm(uF, gF_local) / norm2(uF) <= 2 + 1e-8

    rank == 0 && benchmark && micro_benchmarks(u, uF, gF_global, gF_local)

    MPI.Barrier(comm)
end

MPI.Initialized() || MPI.Init()
silence_stdout(MPI.COMM_WORLD)

test_rfft(DATA_DIMS_EVEN)
println()
test_rfft(DATA_DIMS_ODD, benchmark=false)


text:
module StoreCredit

function sc()
    fl = open("A-large-practice.in")
    num_testcase = readline(fl) |> x->strip(x) |> x->parse(Int64,x)
    outfl = open("large_storecredit.out","w")
    i=0
    while (num_testcase>=1)
        i+=1
        c  = readline(fl) |> x->strip(x) |> x->parse(Int64,x)
        l  = readline(fl) |> x->strip(x) |> x->parse(Int64,x)
        ps = readline(fl) |> x->strip(x) |> x->split(x, " ") |> x->map(y->parse(Int64,y),x)
        findinds(outfl,c,l,ps,i)
        num_testcase -= 1
    end
    close(outfl)
    close(fl)
end

function findinds(outfl,c,l,ps,idcase)
    @assert length(ps) == l
    indperm = sortperm(ps)
    ps = ps[indperm]
    ind1,ind2=0,0
    for i in 1:l-1
        for j in i+1:l
            s = ps[i]+ps[j]
            if s == c
                ind1,ind2 = sort([indperm[i],indperm[j]])
            end
        end
    end
    println(outfl, "Case #$(idcase): $ind1 $ind2")
end

end


text:
function getlines(file)
    # Get lines from the file
    f = open(file)
    lines = readlines(f)
    close(f)

    # Clean the elements from whitespaces
    cleanlines = []
    for i in 1:size(lines, 1)
        if lines[i] != " " && lines[i] != ""
            lines[i] = replace(lines[i], "\t" => "")
            lines[i] = strip(lines[i])
            push!(cleanlines, lines[i])
        end
    end
    return cleanlines
end

function getlistnums(lines)
    for i in 1:size(lines, 1)
        if lines[i] == "("
            return parse(Int64, lines[i - 1]), i + 1
        end
    end
end

function getfileitems(file)
    lines = getlines(file)
    listsize, liststart = getlistnums(lines)
    return lines, listsize, liststart
end

function linetolist(line)
    line = replace(line, "(" => "")
    line = replace(line, ")" => "")
    list = split(line, " ")
    return list
end

function linetokeyvalue(line)
    line = replace(line, ";" => "")
    kvlist = split(line, " ")
    return kvlist[1], kvlist[end]
end

function stringtofloatlist(line)
    list = linetolist(line)
    list = parse.(Float64, list)
    return list
end

function stringtointlist(line)
    list = linetolist(line)
    list = parse.(Int64, list)
    return list
end


text:
#Function(Snakesort)

#define First_Core            m->storage[0]
#define Last_Core             m->storage[1]
#define max_weights           m->storage[2]
#define weights_mode          m->storage[3]

# Start_Initialization
# // Verify storage size.
# if (m->storage.size() != 4) {
#   __debugbreak();
# }

function initialize_snakesort2(core)
    weights_ptr_pi      = ag0pi
    weights_ptr         = ag0
    weights_ptr_config  = AG0START
    weights_ptr_config1 = AG1START
    weights_ptr_im      = AG2START

    weights_mode_int = Loc(AsapSim.DMEM, 255)
    ag_config_weights = ((max_weights<<8) + 1)
    ag_config_im = (((max_weights+20)<<8) + 21)

    # Set up address generator strides.
    MOVI(AG0STRIDE, 1)
    MOVI(weights_ptr_config, ag_config_weights)
    MOVI(weights_ptr_config1, ag_config_weights)
    MOVE(weights_mode_int, weights_mode)
    MOVI(weights_ptr_im, ag_config_im)
    MACCL(null, 0, 0)
end

@asap4asm function snakesort2(firstcore, lastcore, max_weights, weights_mode)
    weights_ptr_pi      = :ag0pi
    weights_ptr         = :ag0
    weights_ptr_config  = Loc(AsapSim.AG_START, 0)
    weights_ptr_config1 = Loc(AsapSim.AG_START, 1)
    weights_ptr_im      = Loc(AsapSim.AG_START, 2)

    weights_mode_int = Loc(AsapSim.DMEM, 255)
    ag_config_weights = ((max_weights<<8) + 1)
    ag_config_im = (((max_weights+20)<<8) + 21)

    @label start
    if (weights_mode == 1)
        RPT(max_weights,nop3)
            MOVE(ag0pi, ibuf0,nop3)
        END_RPT()
        MOVI(weights_mode, 0)
    end

    NOP(nop3)
    # if (Last_Core) {
    #     printf("\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d", Dmem[1], Dmem[2], Dmem[3], Dmem[4], Dmem[5], Dmem[6], Dmem[7], Dmem[8], Dmem[9], Dmem[10], Dmem[11], Dmem[12]);
    # }

    MOVE(null, ibuf0)

    @label pass_weights_downstream
    if (!lastcore)
        MOVE(output[0], ibuf0)
        BR(pass_weights_downstream, Z, neg)
        MOVE(null, ibuf0_next)
        BR(pass_weights_downstream, Z, neg)
    end

    if (!lastcore)
        MOVE(output[0], ibuf0)
    end

    if (lastcore)
        MOVE(null, ibuf0)
    end
    RPT(max_weights,nop3)
        MOVE(ag2pi, ibuf0_next, nop3)
        MACL(null, ag1pi,ibuf0,nop3)
    END_RPT()
    NOP(nop3)
    @label pass_im_downstream
    if (!lastcore)
        MOVE(output[0], ibuf0)
        BR(pass_im_downstream, Z, neg)
    end

    if (lastcore)
        MOVI(output[0], acc)
    end

    NOP(nop3)
    # if (Last_Core) {
    #     printf("\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d", Dmem[21], Dmem[22], Dmem[23], Dmem[24], Dmem[25], Dmem[26], Dmem[27], Dmem[28], Dmem[29], Dmem[30], Dmem[31], Dmem[32]);
    # }
    BRL(start)
end


text:
include("pascaltriangle.jl") # define pascaltriangle(m,n)
# include("pascalmax.jl") # degine pascalmax(m,n)

function b2in(vecmb::SparseVector{Int64})

    # this function return the index of many body state at the index of Fock state

    # ignore the last element since it is the coefficient
    Msize = length(vecmb)-1
    Np = sum(vecmb[1:Msize])

    # Lm = mbsize+1
    # Ln = Np + 1
    indfk = 1
    matp = zeros(Int,Msize+1,Np+1);
    pascaltriangle!(Msize,Np,matp) # the size is Msize+1 times Np+1
    # matp = pascaltriangle(Msize,Np) # the size is Msize+1 times Np+1
    # note the indices are m+1 and n+1 for N^m_n

    for jj = 1:Msize
        for ii = 0:vecmb[jj]

            if ii == vecmb[jj]
               continue
            end

            sumparticles = sum(vecmb[1:jj-1]) # numnber of particles taken account so far
            indfk = indfk + matp[Msize-jj+1,Np-ii-sumparticles+1] # (n+m-1)!/n!(m-1)!
            # note the indices are m+1 and n+1 for N^m_n
            # indfk = indfk + pascalmax(Msize-jj,Np-ii-sumparticles) # pascalmax(m,n) = (n+m-1)!/n!(m-1)!

        end
    end

    # unlike Ponomarev's way, I count the index of Fock state in the backwards
    # so that "1" means |Np,0,0,...>.
    indfk = matp[Msize+1,Np+1] + 1 - indfk # the indices are m+1 and n+1 for N^m_n

    return indfk

end

max_stars_repo_name: leios/SOCED

max_issues_repo_name: leios/SOCED

max_forks_repo_name: leios/SOCED


text:
function optimize_route(spatial_data, spatial_data_matrices, cost_data, equipment_data, input_data)

    # weights_ac = snw.weights_ac_cable;
    # snw.meanweight_ac=alg_fac*min(weighy_ac(weighy_ac~=0));
    # weighy_dc=snw.weights_dc_cable;
    # snw.meanweight_dc=alg_fac*min(weighy_dc(weighy_dc~=0));

    # snw.meanweight_dc=alg_fac*min([snw.weights]);
    # snw.meanweight_dc=1;
    # snw.meanweight_ac=snw.meanweight_dc;

    start_position, finish_position = determine_start_and_end_node(spatial_data_matrices, input_data)

    average_weight = sum(spatial_data["segment_weights"]) / size(spatial_data["segment_weights"] ,1)

    minimum_weight = minimum(spatial_data["segment_weights"][spatial_data["segment_weights"] .!= 0])

    from_node = start_position
    to_node = finish_position
    OPEN_COUNT = 1
    path_cost = 0
    offshore_distance = 0
    average_weight = minimum_weight
    target_distance = estimate_initial_distance(start_position, finish_position, spatial_data, average_weight)
    hn = path_cost
    gn = target_distance
    fn = target_distance
    wn = target_distance
    ac_cable = 0
    dc_cable = 0
    ac_dc_transition = 0
    OPEN = my_insert_open(from_node, to_node, hn, gn, fn, wn, ac_dc_transition, ac_cable, dc_cable, offshore_distance, spatial_data)
    OPEN[OPEN_COUNT, 1] = 0

    CLOSED_COUNT = 1
    CLOSED = [from_node]
    NoPath = 1

    #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    # START ALGORITHM
    #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    while (from_node != finish_position) && (NoPath == 1)
        exp_array = my_expand_array(from_node, path_cost, start_position, finish_position, CLOSED, spatial_data, offshore_distance, input_data, average_weight, equipment_data)
        exp_count = size(exp_array, 1)
        # UPDATE LIST OPEN WITH THE SUCCESSOR NODES
        # OPEN LIST FORMAT
        #--------------------------------------------------------------------------
        #IS ON LIST 1/0 |X val |Y val |Parent X val |Parent Y val |h(n) |g(n)|f(n)|
        #--------------------------------------------------------------------------
        #EXPANDED ARRAY FORMAT
        #--------------------------------
        #|X val |Y val ||h(n) |g(n)|f(n)|
        #--------------------------------
        for i = 1 : exp_count
            flag = 0
            if any(OPEN[1 : OPEN_COUNT, 2] .== exp_array[i, 1])
                j = findfirst(OPEN[1 : OPEN_COUNT, 2] .== exp_array[i, 1])
                if exp_array[i, 4] <= OPEN[j, 10]
                    OPEN[j, 3] = from_node
                    OPEN[j, 4] = spatial_data["nodes"][from_node, 2]
                    OPEN[j, 5] = spatial_data["nodes"][from_node, 3]
                    OPEN[j, 8] = exp_array[i, 2]
                    OPEN[j, 9] = exp_array[i, 3]
                    OPEN[j, 10] = exp_array[i, 4]
                    OPEN[j, 11] = exp_array[i, 5]
                    OPEN[j, 12] = exp_array[i, 6]
                    OPEN[j, 13] = exp_array[i, 7]
                    OPEN[j, 14] = exp_array[i, 8]
                    OPEN[j, 15] = exp_array[i, 9]
                end
            else
                OPEN_COUNT = OPEN_COUNT + 1
                new_line = my_insert_open(exp_array[i, 1], from_node, exp_array[i, 2], exp_array[i, 3], exp_array[i,4], exp_array[i,5], exp_array[i,6], exp_array[i,7], exp_array[i,8], exp_array[i,9], spatial_data)
                OPEN = [OPEN; new_line]
            end
        end

        index_min_node = my_min_fn(OPEN, OPEN_COUNT, finish_position)
        if index_min_node != -1
            from_node = Int.(OPEN[index_min_node, 2])
            path_cost = OPEN[index_min_node, 8]
            offshore_distance = OPEN[index_min_node, 15]
            CLOSED_COUNT = CLOSED_COUNT + 1
            CLOSED = [CLOSED; from_node]
            OPEN[index_min_node, 1] = 0
        else
            NoPath = 0
        end
    end
    i = size(CLOSED, 1)
    nodeval = CLOSED[i]
    i = 1
    optimal_path = nodeval
    ac_dc = OPEN[OPEN[:, 2] .== nodeval, 12]
    ac_cab = OPEN[OPEN[:, 2] .== nodeval, 13]
    dc_cab = OPEN[OPEN[:, 2] .== nodeval, 14]
    i = i + 1
    c_tot = 0
    if nodeval == finish_position
        parent_node = Int.(OPEN[OPEN[:, 2] .== nodeval, 3][1])
        c_tot = c_tot + OPEN[OPEN[:,2] .== nodeval, 11][1]
        while parent_node != start_position
            optimal_path = [optimal_path; parent_node]
            ac_dc = [ac_dc; OPEN[OPEN[:, 2] .== parent_node, 12][1]]
            ac_cab = [ac_cab; OPEN[OPEN[:, 2] .== parent_node, 13[1]]]
            dc_cab = [dc_cab; OPEN[OPEN[:, 2] .== parent_node, 14][1]]
            c_tot = c_tot + OPEN[OPEN[:, 2] .== parent_node, 11][1]
            parent_node = Int.(OPEN[OPEN[:, 2] .== parent_node, 3][1])
            i = i + 1
        end
        optimal_path = [optimal_path; start_position]
    end

    return c_tot, optimal_path, ac_dc, ac_cab, dc_cab
end

function determine_start_and_end_node(spatial_data_matrices, input_data)
        idx1 = findall(spatial_data_matrices["ohl"]["nodes"][:, 2] .== input_data["start_node"]["x"])
        idx2 = findall(spatial_data_matrices["ohl"]["nodes"][:, 3] .== input_data["start_node"]["y"])
        start_position = intersect(idx1, idx2)

        idx3 = findall(spatial_data_matrices["ohl"]["nodes"][:, 2] .== input_data["end_node"]["x"])
        idx4 = findall(spatial_data_matrices["ohl"]["nodes"][:, 3] .== input_data["end_node"]["y"])
        finish_position = intersect(idx3, idx4)

        sp = start_position[1]
        fp = finish_position[1]

        return sp, fp
end

function estimate_initial_distance(start_position, finish_position, spatial_data, average_weight)
    distance = sqrt.((spatial_data["nodes"][start_position, 2] - spatial_data["nodes"][finish_position, 2]).^2 + (spatial_data["nodes"][start_position, 3] - spatial_data["nodes"][finish_position, 3]).^2) * average_weight

    d = distance[1]
    return d
end

function my_insert_open(from_node, to_node, hn, gn, fn, wn, ac_dc, ac_cab, dc_cab, offshore_distance, spatial_data)
    new_row = zeros(1, 15)
    new_row[1,1] = 1
    new_row[1,2] = Int.(from_node)
    new_row[1,3] = Int.(to_node)
    new_row[1,4] = spatial_data["nodes"][Int.(from_node) ,2]
    new_row[1,5] = spatial_data["nodes"][Int.(from_node) ,3]
    new_row[1,6] = spatial_data["nodes"][Int.(to_node) ,2]
    new_row[1,7] = spatial_data["nodes"][Int.(to_node) ,3]
    new_row[1,8] = hn
    new_row[1,9] = gn
    new_row[1,10] = fn
    new_row[1,11] = wn
    new_row[1,12] = ac_dc
    new_row[1,13] = ac_cab
    new_row[1,14] = dc_cab
    new_row[1,15] = offshore_distance

    return new_row
end

function my_expand_array(from_node, hn, start_node, finish_node, CLOSED, spatial_data, offshore_distance, input_data, average_weight, equipment_data)
        seg_index1 = findall(spatial_data["segments"][:, 2] .== from_node)
        seg_index2 = findall(spatial_data["segments"][:, 3] .== from_node)

        exp_array = Array{Float64, 2}(undef, 0, 9)
        for idx = 1:length(seg_index1)
            current_node = Int.(spatial_data["segments"][seg_index1[idx], 3])
            if !any(CLOSED .== current_node)
                seg_index = seg_index1[idx]
                new_line = zeros(1, 9)
                new_line[1, 1] = current_node
                distance = sqrt((spatial_data["nodes"][from_node, 2] - spatial_data["nodes"][current_node, 2])^2 + (spatial_data["nodes"][from_node, 3] - spatial_data["nodes"][current_node, 3])^2)
                d_off = input_data["resolution_factor"] * distance * (1 - spatial_data["onshore_flag"][seg_index]) * (1 - spatial_data["onshore_nodes"][from_node]) * (1 - spatial_data["onshore_nodes"][current_node])
                w_off = 1
                if (offshore_distance + d_off) > equipment_data["ugc"]["offshore"]["maximum_distance"]
                    w_off = 1e9
                end
                new_line[1, 2] = hn + distance_w(from_node, current_node, spatial_data, seg_index, w_off)
                if d_off == 0
                    new_line[1, 9] = 0
                else
                    new_line[1, 9] = offshore_distance + d_off
                end
                od = distance_off(finish_node, current_node, spatial_data, seg_index, input_data)
                new_line[1, 3] = distance_est(finish_node, current_node, spatial_data, average_weight)

                new_line[1, 4] = new_line[1, 2] + new_line[1, 3]
                new_line[1, 5] = new_line[1, 2] - hn
                new_line[1, 6] = 0 #nw.ac_dc(seg_index)
                new_line[1, 7] = spatial_data["ohl_ugc_transition"][seg_index1[idx]]
                new_line[1, 8] = spatial_data["ohl_ugc_transition"][seg_index1[idx]]
                exp_array = [exp_array; new_line]
            end
        end

        for idx = 1:length(seg_index2)
            current_node = Int.(spatial_data["segments"][seg_index2[idx], 2])
            if !any(CLOSED .== current_node)
                seg_index = seg_index2[idx]
                new_line = zeros(1, 9)
                new_line[1, 1] = current_node
                distance = sqrt((spatial_data["nodes"][from_node, 2] - spatial_data["nodes"][current_node, 2])^2 + (spatial_data["nodes"][from_node, 3] - spatial_data["nodes"][current_node, 3])^2)
                d_off = input_data["resolution_factor"] * distance * (1 - spatial_data["onshore_flag"][seg_index]) * (1 - spatial_data["onshore_nodes"][from_node]) * (1 - spatial_data["onshore_nodes"][current_node])
                w_off = 1
                if (offshore_distance + d_off ) > equipment_data["ugc"]["offshore"]["maximum_distance"]
                    w_off = 1e9
                end
                new_line[1, 2] = hn + distance_w(from_node, current_node, spatial_data, seg_index, w_off)
                if d_off == 0
                    new_line[1, 9] = 0
                else
                    new_line[1, 9] = offshore_distance + d_off
                end
                od = distance_off(finish_node, current_node, spatial_data, seg_index, input_data)
                new_line[1, 3] = distance_est(finish_node, current_node, spatial_data, average_weight)

                new_line[1, 4] = new_line[1, 2] + new_line[1, 3]
                new_line[1, 5] = new_line[1, 2] - hn
                new_line[1, 6] = 0 #nw.ac_dc(seg_index)
                new_line[1, 7] = spatial_data["ohl_ugc_transition"][seg_index2[idx]]
                new_line[1, 8] = spatial_data["ohl_ugc_transition"][seg_index2[idx]]
                exp_array = [exp_array; new_line]
            end
        end
        return exp_array
end

function distance_w(from_node, to_node, spatial_data, index, w_off)

    distance = sqrt((spatial_data["nodes"][from_node, 2] - spatial_data["nodes"][to_node, 2])^2 + (spatial_data["nodes"][from_node, 3] - spatial_data["nodes"][to_node, 3])^2) * spatial_data["segment_weights"][index] * w_off
    d = distance[1]
   return d
end

function distance_est(from_node, to_node, spatial_data, average_weight)

    distance = sqrt((spatial_data["nodes"][from_node, 2] - spatial_data["nodes"][to_node, 2])^2 + (spatial_data["nodes"][from_node, 3] - spatial_data["nodes"][to_node, 3])^2) * average_weight
    d = distance[1]
    return d
end

function distance_off(from_node, to_node, spatial_data, index, input_data)
    distance = sqrt((spatial_data["nodes"][from_node, 2] - spatial_data["nodes"][to_node, 2])^2 + (spatial_data["nodes"][from_node, 3] - spatial_data["nodes"][to_node, 3])^2) * (1-spatial_data["onshore_flag"][index]) * input_data["resolution_factor"]
    d = distance[1]
    return d
end

function my_min_fn(OPEN, OPEN_COUNT, finish_node)

    index = (1:OPEN_COUNT)'

    temp_array = [OPEN[OPEN[1:OPEN_COUNT, 1] .== 1, :] index[OPEN[1:OPEN_COUNT, 1] .== 1]]
    goal_index = findall(OPEN[:, 1] .== finish_node)

    if !isempty(goal_index)
        i_min = goal_index
    end

    if size(temp_array, 1) != 0
        temp_min = Int.(findfirst(temp_array[:, 10] .== minimum(temp_array[:, 10])))
        i_min = Int.(temp_array[temp_min[1], end])
    else
        i_min = -1
    end
    return i_min
end


text:
using GeoStatsImages
using Base.Test

# no tests


text:
export RandomBatch

"In Memory Data Set"
struct RandomBatch{InputType, OutputType} <: AbstractDataSet{InputType, OutputType}
    data::MemoryDataSet{InputType, OutputType}
    parentSet::AbstractDataSet{InputType, OutputType}

    function RandomBatch(set::AbstractDataSet{I, O}, batch_size::Int) where {I, O}
        newSet = new{I, O}(MemoryDataSet{I, O}(batch_size), set)
        reset(newSet)
        return newSet
    end

    Base.getindex(set::RandomBatch, row) = set.data[row]
    Base.setindex!(set::RandomBatch{I, O}, value::DataEntry{I, O}, row) where {I, O} = set.data[row] = value
    Base.deleteat!(set::RandomBatch, idxes) = deleteat!(set, idxes)
    Base.length(set::RandomBatch)::Int64 = length(set.data)

    function Base.reset(set::RandomBatch)
        for i in 1:length(set)
            set[i] = set.parentSet[rand(1:length(set.parentSet))]
        end
    end

end


text:
"""
    ellipse(; center::Point, xAxis::Real, yAxis::Real, steps::Integer=64, angle::Real=0., pivot::Union{Point, Nothing}=nothing, units::String="kilometers")

Take a Point and calculates the ellipse polygon given two semi-axes expressed in variable units and steps for precision.

# Examples
```jldoctest
julia> center = Point([10, 45])
Point([10.0, 45.0])

julia> ellipse(center=center, xAxis=20, yAxis=7.5, steps=10)
Polygon(Array{Array{Float64,1},1}[[[10.0872, 45.0634], [10.0218, 45.0672], [9.97817, 45.0672], [9.91279, 45.0634], [9.74563, 45.0], [9.91279, 44.9366], [9.97817, 44.9328], [10.0218, 44.9328], [10.0872, 44.9366], [10.2544, 45.0], [10.0872, 45.0634]]])
```
"""
function ellipse(; center::Point, xAxis::Real, yAxis::Real, steps::Integer=64, angle::Real=0., pivot::Union{Point, Nothing}=nothing, units::String="kilometers")
    pivot == nothing && (pivot = center)

    coords = center.coordinates

    if units === "degrees"
        angle = deg2rad(angle)
    else
        xAxis = rhumb_destination(coords, xAxis, 90, units).coordinates
        yAxis = rhumb_destination(coords, yAxis, 0, units).coordinates
        xAxis = xAxis[1] - coords[1]
        yAxis = yAxis[2] - coords[2]
    end

    elCoords = []

    for i in 1:steps
        stepAngle = i * (-360 / steps)
        x = (xAxis * yAxis) / sqrt(yAxis^2 + xAxis^2 * tan(deg2rad(stepAngle))^2)
        y = (xAxis * yAxis) / sqrt(xAxis^2 + yAxis^2 / tan(deg2rad(stepAngle))^2)

        (stepAngle < -90 && stepAngle >= -270) && (x = -x)
        (stepAngle < -180 && stepAngle >= -360) && (y = -y)

        if units === "degrees"
            newx = x * cos(angle) + y * sin(angle)
            newy = y * cos(angle) - x * sin(angle)
            x = newx
            y = newy
        end

        push!(elCoords, [x + coords[1], y + coords[2]])
    end
    push!(elCoords, elCoords[1])

    units === "degrees" && return Polygon([elCoords])

    return transform_rotate(geojson=Polygon([elCoords]), angle=angle, pivot=pivot)

end


text:
using Cite
using Documenter

DocMeta.setdocmeta!(Cite, :DocTestSetup, :(using Cite); recursive=true)

makedocs(;
    modules=[Cite],
    authors="Sebastian Micluța-Câmpeanu <m.c.sebastian95@gmail.com>",
    repo="https://github.com/SebastianM-C/Cite.jl/blob/{commit}{path}#{line}",
    sitename="Cite.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://SebastianM-C.github.io/Cite.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/SebastianM-C/Cite.jl",
)


text:
@computed type A{V <: AbstractVector}
    a::eltype(V)
end

@testset "basic example" begin
    a = A{Vector{Int}}(3.0)
    @test a.a === Int(3)
end

@computed type B{N, M, T}
    a::NTuple{N + M, T}
    B(x::T) = new{N, M, T}(ntuple(i -> x, N + M))
    B{S}(x::S) = B{N, M, T}(convert(T, x))
end

@computed type C{T <: Number}
    a::typeof(one(T) / one(T))
    C() = new(0)
    function C(x)
        return new(x)
    end
end


text:
using PtFEM

data = Dict(
  # Frame(nels, nn, ndim, nst, nip, finite_element(nod, nodof))
  :struc_el => Frame(3, 4, 3, 1, 1, Line(2, 3)),
  :properties => [
    4.0e6 1.0e6 0.3e6 0.3e6],
  :x_coords => [0.0, 5.0, 5.0, 5.0],
  :y_coords => [5.0, 5.0, 5.0, 0.0],
  :z_coords => [5.0, 5.0, 0.0, 0.0],
  :gamma => [0.0, 0.0, 90.0],             # Third element 90 degrees from x-y plane
  :g_num => [
    1 3 4;
    2 2 3],
  :support => [
    (1, [0 0 0 0 0 0]),
    (4, [0 0 0 0 0 0])
    ],
  :loaded_nodes => [
    (2, [0.0 -100.0 0.0 0.0 0.0 0.0])]
)

data |> display
println()

@time fem, dis_df, fm_df = p44(data)
println()

dis_df |> display
println()

println("Actions:")
fm_df |> display
println()


text:
function count_unique(x)
    u = unique(x)
    c = Array{Int64}(undef, length(u))
    for i in 1 : length(u)
        c[i] = length(findall(isequal(u[i]), x))
    end
    return u, c
end

function get_branching_shold(xs)
    u = unique(xs)
    return (u[2:end]+u[1:end-1]) / 2
end

function delta_gini_index(left, right, alpha, branch_count)
    n_left, n_right = length(left), length(right)
    n_total = n_left + n_right
    _, counts = count_unique(left)
    ratio_classes = counts / n_left
    left_gain = (n_left/n_total) * (1-sum(ratio_classes.^2))
    _, counts = count_unique(right)
    ratio_classes = counts / n_right
    right_gain = (n_right/n_total) * (1-sum(ratio_classes.^2))
    return left_gain+right_gain+alpha*branch_count
end

function branch(x, t, alpha, branch_count)
    gains = []
    rules = []
    for feature_id in 1:size(x, 2)
        xs = x[:, feature_id]
        thresholds = get_branching_shold(xs)
        for th in thresholds
            left_y = t[findall(xs.<th)]
            right_y = t[findall(xs.>=th)]
            gain = delta_gini_index(left_y, right_y, alpha, branch_count)
            push!(gains, gain)
            if rules == []
                rules = [feature_id th]
            else
                rules = vcat(rules, [feature_id th])
            end
        end
    end
    best_rule = rules[argmin(gains), :]
    best_gini = minimum(gains)
    feature_id = Int(best_rule[1])
    threshold = best_rule[2]
    left = findall(x[:, feature_id] .< threshold)
    right = findall(x[:, feature_id] .>= threshold)
    return x[left, :], t[left], x[right, :], t[right], feature_id, threshold, best_gini
end

function grow(x, t, classes, gini, alpha, branch_count)
    check_size(x, t)
    uniques, counts = count_unique(t)
    counter = Dict(zip(uniques, counts))
    class_count = [(c in keys(counter)) ? counter[c] : 0 for c in classes]
    this = Dict("class_count" => class_count, "feature_id" => nothing, "threshold" => nothing, "left" => nothing, "right" => nothing)
    if length(t) == 1
        return this
    elseif length(uniques) == 1
        return this
    end
    a = true
    for i in 1:size(x, 2)
        a = a && length(unique(x[:, i])) == 1
    end
    a && return this
    left_x, left_t, right_x, right_t, feature_id, threshold, best_gini = branch(x, t, alpha, branch_count)
    if best_gini < gini
        branch_count += 1
        left = grow(left_x, left_t, classes, best_gini, alpha, branch_count)
        right = grow(right_x, right_t, classes, best_gini, alpha, branch_count)
        return Dict("class_count" => class_count, "feature_id" => feature_id, "threshold" => threshold, "left" => left, "right" => right)
    else
        return this
    end
end

"""
    DecisionTree(; alpha = 0.01)
Normal DecisionTree. `alpha` specify the complexity of the model. If it's small, it's complicated, and if it's big, it's simple.

# Example
```jldoctest classification
julia> tree = DecisionTree()
DecisionTree(0.01, Dict{Any, Any}(), Any[])

julia> fit!(tree, x, t)
Dict{String, Any} with 5 entries:
  "left"        => Dict{String, Any}("left"=>Dict{String, Union{Nothing, Vector…
  "class_count" => [8, 13, 16]
  "threshold"   => 5.7
  "right"       => Dict{String, Any}("left"=>Dict{String, Union{Nothing, Vector…
  "feature_id"  => 1

julia> println(predict(tree, x))
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 1, 2, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
```
"""
mutable struct DecisionTree
    α::Float64
    tree::Dict
    classes::Array{Any, 1}
    DecisionTree(; alpha = 0.01) = new(alpha, Dict(), Array{Any}(undef, 0))
end

function fit!(model::DecisionTree, x, t)
    model.classes = unique(t)
    model.tree = grow(x, t, model.classes, Inf, model.α, 1)
end

function predict(model::DecisionTree, x)
    function predict_one(x, tree)
        while tree["feature_id"] != nothing
            tree = (x[tree["feature_id"]] < tree["threshold"]) ? tree["left"] : tree["right"]
        end
        return argmax(tree["class_count"])
    end
    predicts = [predict_one(x[i, :], model.tree) for i in 1:size(x, 1)]
    d = Dict(zip(collect(1:length(model.classes)), model.classes))
    return map(x->d[x], predicts)
end


text:
"""
    randomize!(root::T, num::Int64=100)::Nothing where T <:AbstractNode

This function randomizes the tree by performing a number of nearest
neighbour interchange (NNI) moves and a randomization of the branch lengths.
The number of NNI moves is specified in the parameter num.

* `root` : root node of tree to be edited.

* `num` : amount of NNI moves to perform.
"""
function randomize!(root::T, num::Int64=100)::Nothing where T <:AbstractNode
    n_nodes = size(root)[1]
    i = 0
    while i < num
        n = rand(1:n_nodes)
        NNI!(root, n)
        i += 1
    end
    blv = rand(n_nodes)
    set_branchlength_vector!(root, blv)
end

"""
    randomize(root::T, num::Int64=100)::T where T <:AbstractNode

This function returns a randomized copy of the tree by performing a number of nearest
neighbour interchange (NNI) moves and a randomization of the branch lengths.
The number of NNI moves is specified in the parameter num.

* `root` : root node of tree to be edited.

* `num` : amount of NNI moves to perform.
"""
function randomize(root::T, num::Int64=100)::T where T <:AbstractNode
    R = deepcopy(root)
    randomize!(R, num)
    R
end


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

# Test conversion to and from unix
@test Dates.unix2datetime(Dates.datetime2unix(DateTime(2000,1,1))) == DateTime(2000,1,1)
@test Dates.value(Dates.DateTime(1970)) == Dates.UNIXEPOCH

# Tests from here: https://en.wikipedia.org/wiki/Unix_time
@test string(Dates.unix2datetime(1095379198.75)) == string("2004-09-16T23:59:58.75")
@test string(Dates.unix2datetime(1095379199.00)) == string("2004-09-16T23:59:59")
@test string(Dates.unix2datetime(1095379199.25)) == string("2004-09-16T23:59:59.25")
@test string(Dates.unix2datetime(1095379199.50)) == string("2004-09-16T23:59:59.5")
@test string(Dates.unix2datetime(1095379199.75)) == string("2004-09-16T23:59:59.75")
@test string(Dates.unix2datetime(1095379200.00)) == string("2004-09-17T00:00:00")
@test string(Dates.unix2datetime(1095379200.25)) == string("2004-09-17T00:00:00.25")
@test string(Dates.unix2datetime(1095379200.50)) == string("2004-09-17T00:00:00.5")
@test string(Dates.unix2datetime(1095379200.75)) == string("2004-09-17T00:00:00.75")
@test string(Dates.unix2datetime(1095379201.00)) == string("2004-09-17T00:00:01")
@test string(Dates.unix2datetime(1095379201.25)) == string("2004-09-17T00:00:01.25")
@test string(Dates.unix2datetime(915148798.75)) == string("1998-12-31T23:59:58.75")
@test string(Dates.unix2datetime(915148799.00)) == string("1998-12-31T23:59:59")
@test string(Dates.unix2datetime(915148799.25)) == string("1998-12-31T23:59:59.25")
@test string(Dates.unix2datetime(915148799.50)) == string("1998-12-31T23:59:59.5")
@test string(Dates.unix2datetime(915148799.75)) == string("1998-12-31T23:59:59.75")
@test string(Dates.unix2datetime(915148800.00)) == string("1999-01-01T00:00:00")
@test string(Dates.unix2datetime(915148800.25)) == string("1999-01-01T00:00:00.25")
@test string(Dates.unix2datetime(915148800.50)) == string("1999-01-01T00:00:00.5")
@test string(Dates.unix2datetime(915148800.75)) == string("1999-01-01T00:00:00.75")
@test string(Dates.unix2datetime(915148801.00)) == string("1999-01-01T00:00:01")
@test string(Dates.unix2datetime(915148801.25)) == string("1999-01-01T00:00:01.25")

# Test conversion to and from Rata Die
@test Date(Dates.rata2datetime(734869)) == Dates.Date(2013,1,1)
@test Dates.datetime2rata(Dates.rata2datetime(734869)) == 734869

# Tests from here: http://mysite.verizon.net/aesir_research/date/back.htm#JDN
@test Dates.julian2datetime(1721119.5) == Dates.DateTime(0,3,1)
@test Dates.julian2datetime(1721424.5) == Dates.DateTime(0,12,31)
@test Dates.julian2datetime(1721425.5) == Dates.DateTime(1,1,1)
@test Dates.julian2datetime(2299149.5) == Dates.DateTime(1582,10,4)
@test Dates.julian2datetime(2415020.5) == Dates.DateTime(1900,1,1)
@test Dates.julian2datetime(2415385.5) == Dates.DateTime(1901,1,1)
@test Dates.julian2datetime(2440587.5) == Dates.DateTime(1970,1,1)
@test Dates.julian2datetime(2444239.5) == Dates.DateTime(1980,1,1)
@test Dates.julian2datetime(2452695.625) == Dates.DateTime(2003,2,25,3)
@test Dates.datetime2julian(Dates.DateTime(2013,12,3,21)) == 2456630.375

@test typeof(Dates.now()) <: Dates.DateTime
@test typeof(Dates.today()) <: Dates.Date
@test typeof(Dates.now(Dates.UTC)) <: Dates.DateTime

if is_apple()
    withenv("TZ" => "UTC") do
        @test abs(Dates.now() - now(Dates.UTC)) < Dates.Second(1)
    end
end
@test abs(Dates.now() - now(Dates.UTC)) < Dates.Hour(16)

# Issue #9171, #9169
let t = Dates.Period[Dates.Week(2), Dates.Day(14), Dates.Hour(14*24), Dates.Minute(14*24*60), Dates.Second(14*24*60*60), Dates.Millisecond(14*24*60*60*1000)]
    for i = 1:length(t)
        Pi = typeof(t[i])
        for j = 1:length(t)
            @test t[i] == t[j]
            @test Int(convert(Pi,t[j])) == Int(t[i])
        end
        for j = i+1:length(t)
            Pj = typeof(t[j])
            tj1 = t[j] + one(Pj)
            @test t[i] < tj1
            @test_throws InexactError Pi(tj1)
            @test_throws InexactError Pj(Pi(typemax(Int64)))
            @test_throws InexactError Pj(Pi(typemin(Int64)))
        end
    end
end
@test Dates.Year(3) == Dates.Month(36)
@test Int(convert(Dates.Month, Dates.Year(3))) == 36
@test Int(convert(Dates.Year, Dates.Month(36))) == 3
@test Dates.Year(3) < Dates.Month(37)
@test_throws InexactError convert(Dates.Year, Dates.Month(37))
@test_throws InexactError Dates.Month(Dates.Year(typemax(Int64)))

# Ensure that conversion of 32-bit integers work
let dt = DateTime(1915,1,1,12)
    unix = Int32(Dates.datetime2unix(dt))
    julian = Int32(Dates.datetime2julian(dt))

    @test Dates.unix2datetime(unix) == dt
    @test Dates.julian2datetime(julian) == dt
end

# Conversions to/from numbers
a = Dates.DateTime(2000)
b = Dates.Date(2000)
@test convert(Real,b) == 730120
@test convert(Float64,b) == 730120.0
@test convert(Int32,b) == 730120
@test convert(Real,a) == 63082368000000
@test convert(Float64,a) == 63082368000000.0
@test convert(Int64,a) == 63082368000000
@test convert(DateTime,63082368000000) == a
@test convert(DateTime,63082368000000.0) == a
@test convert(Date,730120) == b
@test convert(Date,730120.0) == b
@test convert(Date,Int32(730120)) == b

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
using SnoopCompile

tinf = @snoopi_deep include("runtests.jl")
itrigs = inference_triggers(tinf)

fitrigs = filter(itrig -> itrig.callerframes[end].linfo.def.module === FranklinParser, itrigs)

# Jan 30, 3 triggers, can probably ignore.
# July 16: itrigs has 4 elements for Base._* (probably irrelevant)
#          fitrigs has 0 elements


text:
// count function parameters as initialized

int main() {
  printInt(fac(5));
  return 0 ;
}

int fac (int a) {
  int r;
  int n;
  r = 1;
  n = a;
  while (n > 0)
  {
    r = r * n;
    n = n - 1;
  }
  return r;
}

max_stars_repo_name: jwilk/jtc

max_issues_repo_name: jwilk/jtc

max_forks_repo_name: jwilk/jtc


text:
using Xlib

function x11demo()
    # Open connection to the server.
    dpy = XOpenDisplay(C_NULL)
    dpy == C_NULL && error("unable to open display")
    scr = DefaultScreen(dpy)

    # Create a window.
    win = XCreateSimpleWindow(dpy, RootWindow(dpy, scr), 10, 10, 300, 100, 1,
                              BlackPixel(dpy, scr), WhitePixel(dpy, scr))

    # Select the kind of events we are interested in.
    XSelectInput(dpy, win, ExposureMask | KeyPressMask)

    # Show or in x11 terms map window.
    XMapWindow(dpy, win)

    # Run event loop.
    evt = Ref(XEvent())
    while true
        XNextEvent(dpy, evt)

        # Draw or redraw the window.
        if EventType(evt) == Expose
            XFillRectangle(dpy, win, DefaultGC(dpy, scr), 24, 24, 16, 16)
            XDrawString(dpy, win, DefaultGC(dpy, scr), 50, 50, "Hello, World! Press any key to exit.")
        end

        # Exit whenever a key is pressed.
        if EventType(evt) == KeyPress
            break
        end
    end

    # Shutdown server connection
    XCloseDisplay(dpy)
end

x11demo()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@testset "1103.distribute-candies-to-people.jl" begin
    @test distribute_candies(7, 4) == [1, 2, 3, 1]
    @test distribute_candies(10, 3) == [5, 2, 3]
    @test distribute_candies(1958512, 40) == [
        49050,
        49100,
        49150,
        49200,
        49250,
        49300,
        49350,
        49400,
        49450,
        49500,
        49550,
        49600,
        49650,
        49700,
        49750,
        49800,
        49850,
        49900,
        49252,
        48020,
        48069,
        48118,
        48167,
        48216,
        48265,
        48314,
        48363,
        48412,
        48461,
        48510,
        48559,
        48608,
        48657,
        48706,
        48755,
        48804,
        48853,
        48902,
        48951,
        49000,
    ]
end


text:
# Code used to generate SDP problems with a quadratic objective and a system of equality constraints save in Testproblems/SDPandEqConstraints with an accurante
# solution from MOSEK
# The problem data is generated based on the Optimality conditions:
# Ax*+s*=b, s in K, Px*+q+A'y*=0, y in K*

# steps:
# 1. Create random matrix A in R^mxn with 50% nonzeros
# 2. Half the elements in s belong to equality constraints, other half to one sdp cone, s=[s1;s2], with s1 in {0}-cone and s2 in psd-cone
# 3. Create random solution vector xtrue and strue with s1 = 0s, and mat(s2) in psd-cone
# 4. Determine right side of Ax+s=b
# 5. Create random pos def matrix P, and random dual variable ytrue=[y1;y2] with y1 in R+ and y2 in psd-cone
# 6. calculate q based on P,xtrue and ytrue

# workspace()
# include("./Helper.jl")
# include("../../OSSDP/Code/src/Solver.jl")

# using JuMP, Mosek, JLD, HelperFunctions, OSSDP, Base.Test

rng = MersenneTwister(12345)
dirPath = "../DataFiles/Julia/SDPandEqConstraints/"
!ispath(dirPath) && mkdir(dirPath)
nn = 30

   # @testset "Random SDP Mosek vs. QOCS" begin

 for iii =1:1:nn
  # choose size of problem
  n = rand(rng,10:100)
  r = rand(rng,2:10)
  # dimension of psd-cone (vector form)
  m2 = r^2
  # number of eq constraints
  m1 = m2
  m = m1+m2

  A = sprand(rng,m,n,0.5)
  xtrue = randn(rng,n,1)

  # create primal feasible point -> b
  s1 = zeros(m1,1)
  S2 = generatePosDefMatrix(r,rng)
  s2 = vec(S2)
  strue = [s1;s2]
  b = A*xtrue + strue

  # create dual feasible point -> P,q
   P = generatePosDefMatrix(n,rng)
   Phalf = P^(1/2)

  y1 = rand(rng,m1,1)
  Y2 = generatePosDefMatrix(r,rng)
  y2 = vec(Y2)
  ytrue = [y1;y2]
  q = (-P*xtrue -  A'*ytrue)[:]
  ra = 0.
  Kf = m1
  Kl = 0
  Kq = []
  Ks = [r^2]

  PhalfInvQ=Phalf\q

  # K = Cone(m1,0,[],[r^2])
  # setOFF = OSSDPSettings(rho=0.1,sigma=1e-6,alpha=1.6,max_iter=1500,verbose=false,checkTermination=1,scaling = 0,eps_abs = 1e-5,eps_rel=1e-5,adaptive_rho=false)
  # setON = OSSDPSettings(rho=0.1,sigma=1e-6,alpha=1.6,max_iter=1500,verbose=false,checkTermination=1,scaling = 10,eps_abs = 1e-5,eps_rel=1e-5,adaptive_rho=false)
  # setAdaptiveOFF = OSSDPSettings(rho=0.1,sigma=1e-6,alpha=1.6,max_iter=1500,verbose=false,checkTermination=1,scaling = 0,eps_abs = 1e-5,eps_rel=1e-5,adaptive_rho=true)
  # setAdaptiveON = OSSDPSettings(rho=0.1,sigma=1e-6,alpha=1.6,max_iter=1500,verbose=false,checkTermination=1,scaling = 10,eps_abs = 1e-5,eps_rel=1e-5,adaptive_rho=true)
  # res1,nothing = OSSDP.solve(P,q,A,b,K,setOFF);
  # print("\n.")
  # res2,nothing = OSSDP.solve(P,q,A,b,K,setON);
  # print(".")
  # res3,nothing = OSSDP.solve(P,q,A,b,K,setAdaptiveOFF);
  # print(".")
  # res4,nothing = OSSDP.solve(P,q,A,b,K,setAdaptiveON);
  # print(".")

  # println(">>>QOCS: Problem $(iii), res1: cost: $(res1.cost), status: $(res1.status)")
  # println(">>>QOCS: Problem $(iii), res2: cost: $(res2.cost), status: $(res2.status)")
  # println(">>>QOCS: Problem $(iii), res3: cost: $(res3.cost), status: $(res3.status)")
  # println(">>>QOCS: Problem $(iii), res4: cost: $(res4.cost), status: $(res4.status)")

  # solve accurately once with mosek s
  model = Model(solver=MosekSolver())
  @variable(model, x[1:n])
  @variable(model, s1[1:m1])
  @variable(model, S2[1:r,1:r],SDP)

  # in case of quadratic objective the following problem is solved:
  # min x'Px + 2q'x
  # s.t. Ax+s=b, s in S+
  # instead solve equivalent problem
  # min t
  # s.t. ||Phalf*x+inv(Phalf)*q|| <= t
  # Ax+s=b, s in S+
  @variable(model, t)
  @objective(model, Min, t)

  @constraint(model,norm(Phalf*x+PhalfInvQ) <= t)
  s = [s1;vec(S2)]
  @constraint(model,A*x+s.==b)
  @constraint(model,s1.== 0.)
  status = JuMP.solve(model)
  objTrue = getobjectivevalue(model)

  objTrue = 0.5*( objTrue^2-q'*(P\q) )
  solTrue = getvalue(x)

 # @test abs(objTrue-res1.cost) < 1e-2
 # @test norm(solTrue-res1.x,Inf) < 1e-2

  nr = "$(iii)"
  if iii < 10
    nr = "0$(iii)"
  end

  fn = "SDPandEqConstraints"*nr*".jld"
  problemType = "Random SDP with equality constraints"
  problemName = "SDPandEqConstraints"*nr

   JLD.save(dirPath*fn,"n",n,"m",m,"A",A,"b",b,"P",P,"q",q,"r",ra,"objTrue",objTrue,"solTrue",solTrue,"problemType",problemType,"problemName",problemName,"Kf",Kf,"Kl",Kl,"Kq",Kq,"Ks",Ks)
  println("$(iii)/$(nn) completed!")
end
   # end


text:
#
include("Yfind.jl")
include("is_it_exit2.jl")
include("is_it.jl")

@everywhere function escape_exit_function_parallel(mesh_list,t_end,Energy)
    Q=mesh_list[1]; P=mesh_list[2];
    H=(2*Energy)^2
    tol_dist=1e-5
    Y=Yfind(Q,P,H)
    if ~isempty(Y)
        u0=zeros(5)
        u0[1]=0 #X
        u0[2]=P #P
        u0[3]=Q #Q
        u0[4]=Y[1] #Y
        u0[5]=0 #Y

        prob = ODEProblem(Eq_of_M,u0,(0., t_end))
        # sol=solve(prob,RK4(),maxiters=1e20, reltol=1e-6,abstol=1e-8,callback=cb)
        sol=solve(prob, Vern9(),maxiters=1e10, reltol=1e-8,abstol=1e-10,callback=cb)
        uf=zeros(5)
        uf[1]=sol[1,end] #X
        uf[2]=sol[2,end] #P
        uf[3]=sol[3,end] #Q
        uf[4]=sol[4,end] #Y
        uf[5]=0
        dH=abs(H_test(uf)-H)
        if dH<1e-5
            if is_it(uf,H,tol_dist)
                return sol.u[end][5]
            else
                return -1
            end
        end
    end
return -2
end


text:
using sparsePolyChaos
using PolyChaos

maxDeg = 10
interactions = 1

## Define the model
# μ = 0
# σ = 1
# model(x) = exp(μ + σ * x)
# model(x) = x^7 - 21 * x^5 + 105 * x^3 - 105 * x + x^4 - 6 * x^2 + 3 # He7 + He4
# model(x) = x^4 - 3* x^2
model(x) = x^10 - 45 * x^8 + 630 * x^6 - 3150 * x^4 + 4725 * x^2 - 945 # He10. Bad conditioning: Restart needed until sample size = 1800
# model(x) = 0

## Build the full model
op = GaussOrthoPoly(maxDeg)

## Compute sparse PCE solution
pce = sparsePCE(op, model; pMax = maxDeg, jMax = interactions)

println("Basis polynomials: ", pce[:Ap])
println("PCE coefficients: ", pce[:coeffs])
# println("Max degree:: ", p)
println("R² error: ", pce[:R²])
println("Q² error: ", pce[:Q²])

# ## Sparse PCE with custom sampling method

# println("\n\n== Sparse PCE with custom sampling method ==")
# function sampleFromGaussianMixture(n::Int,μ::Vector{},σ::Vector{},w::Vector{})
#     X = Float64[]
#     for i in 1:n
#         k = findfirst(x -> x > rand(), cumsum(w))
#         push!(X, μ[k] + σ[k]*randn())
#     end
#     return X
# end

# function ρ_gauss(x,μ,σ)
#     1 / sqrt(2*π*σ^2) * exp(-(x - μ)^2 / (2σ^2))
# end

# μ, σ, w = [2.1, 3.2], [0.3, 0.4], [0.3, 0.7]
# ρ(x) = sum( w[i]*ρ_gauss(x,μ[i],σ[i]) for i in 1:length(w) )
# meas = Measure("my_GaussMixture", ρ, (-Inf,Inf), false, Dict(:μ=>μ,:σ=>σ,:w=>w)) # build measure
# op_gm = OrthoPoly("my_op",maxDeg,meas;Nquad=150,Nrec = 5*maxDeg, discretization=stieltjes) # construct orthogonal polynomial

# sampleFun(sampleSize) = sampleFromGaussianMixture(sampleSize,μ,σ,w)
# pce = sparsePCE(op_gm, model, sampleFun; pMax = maxDeg)

# println("Basis polynomials: ", pce[:Ap])
# println("PCE coefficients: ", pce[:coeffs])
# # println("Max degree:: ", pce[:p])
# println("R² error: ", pce[:R²])
# println("Q² error: ", pce[:Q²])

#TODO: Comparison to monte carlo
# # ------- Monte-Carlo -------
# # Evaluate model function on same set of previously drawn samples X
# y_mc = model.(X)

max_stars_repo_name: adriangrupp/sparsePolyChaos

max_issues_repo_name: adriangrupp/sparsePolyChaos

max_forks_repo_name: adriangrupp/sparsePolyChaos


text:
export FisherTest, FisherTestResult

struct FisherTest <: AbstractHLAModel end

struct FisherTestResult
    replacements::Vector{Replacement}
    counts::Vector{Dict{HLAAllele, Array{Int, 2}}}
    p_values::Vector{Dict{HLAAllele, Float64}}
    log_odds::Vector{Dict{HLAAllele, Float64}}
end

function Escape.run(model::FisherTest, data::AbstractHLAData; fdr::Bool = false)

    r = replacements(data)
    m = hla_matrix(data.hla_types; allele_depth = allele_depth)
    alleles = sort(unique_alleles(data.hla_types, allele_depth = allele_depth))

    counts = [Dict{HLAAllele, Array{Int, 2}}() for _ in 1:length(r)]
    p_values = [Dict{HLAAllele, Float64}() for _ in 1:length(r)]
    log_odds = [Dict{HLAAllele, Float64}() for _ in 1:length(r)]

    @showprogress "Fisher's exact test " for (i, replacement) in enumerate(r)
        y = targets(replacement, data)

        for (j, allele) in enumerate(alleles)
            mutation_present = Bool.(collect(skipmissing(y)))
            allele_present = Bool.(min.(1, m[.!ismissing.(y), j]))

            c, p, beta = fisher_exact_test(mutation_present, allele_present)
            counts[i][allele] = c
            p_values[i][allele] = p
            log_odds[i][allele] = beta
        end
    end

    if fdr
        adjust_p_values!(p_values)
    end

    return FisherTestResult(r, counts, p_values, log_odds)
end

function adjust_p_values!(x::Vector{Dict{HLAAllele, Float64}})
    for d in x
        p_values = collect(values(d))
        adjusted_p_values = adjust(p_values, BenjaminiHochberg())
        a = Dict(p_values .=> adjusted_p_values)

        for (k, v) in d
            d[k] = a[v]
        end
    end

    return x
end

function Escape.run(model::FisherTest, data::AbstractHLAData, replacement::Replacement;
                    wp = WorkerPool(workers()), allele_depth::Int = 1)
    y = targets(replacement, data)
    m = hla_matrix(data.hla_types; allele_depth = allele_depth)

    alleles = sort(unique_alleles(data.hla_types, allele_depth = allele_depth))

    counts = Dict{HLAAllele, Array{Int, 2}}()
    p_values = Dict{HLAAllele, Float64}()
    log_odds = Dict{HLAAllele, Float64}()

    for (i, allele) in enumerate(alleles)
        mutation_present = Bool.(collect(skipmissing(y)))
        allele_present = Bool.(min.(1, m[.!ismissing.(y), i]))

        c, p, beta = fisher_exact_test(mutation_present, allele_present)
        counts[allele] = c
        p_values[allele] = p
        log_odds[allele] = beta
    end

    return FisherTestResult([replacement], [counts], [p_values], [log_odds])
end

function fisher_exact_test(a::AbstractVector{Bool}, b::AbstractVector{Bool})
    c = zeros(Int, 2, 2)

    for i in eachindex(a)
        c[2 - a[i], 2 - b[i]] += 1
    end

    if (c[:, 1] == [0, 0]) | (c[:, 2] == [0, 0])
        return (counts = c, p = 1, beta = Inf)
    end

    f = FisherExactTest(c[1,1], c[1,2], c[2,1], c[2,2])

    return (counts = c, p = pvalue(f, method = :minlike), beta = log(f.ω))
end


text:
test = (4, 8)
puzzle = (1, 10)

rot(v::Integer, m::Integer = 10) = mod(v - 1, m) + 1

function q1(p1::Integer, p2::Integer)::Int
    roll(i) = sum(rot.([i, i + 1, i + 2], 100))

    rolls = [[roll(i), roll(i + 3)] for i in 1:6:1000]
    fields = [rot.(fs) for fs in cumsum([[p1, p2], rolls...])]
    scores = cumsum(fields[2:end])
    last_round = findfirst(v -> any(>=(1000), v), scores)
    winner = argmax(scores[last_round])
    loser_score = winner == 1 ? scores[last_round - 1][2] : scores[last_round][1]
    total_rolls = 3 * (2 * (last_round - 1) + winner)
    return total_rolls * loser_score
end

function q2(p1::Integer, p2::Integer)::Int
    (s1, s2) = quantum_scores.((p1, p2))
    return foldl(2:size(s1, 1); init = (0, 0)) do (w1, w2), r
        # p1 reaches 21 exactly after r rounds this many times
        p1w = sum(s1[r, 21:end]) - 27 * sum(s1[r - 1, 21:end])
        # p2 stays below 21 after r rounds this many times
        p2l = sum(s2[r - 1, 1:20])
        # vice versa
        p2w = sum(s2[r, 21:end]) - 27 * sum(s2[r - 1, 21:end])
        p1l = sum(s1[r, 1:20])
        # count wins
        return (w1 + p1w * p2l, w2 + p2w * p1l)
    end |> maximum
end

function quantum_scores(start::Integer; rounds = 10, max_score = 90)
    rolls = [3 => 1, 4 => 3, 5 => 6, 6 => 7, 7 => 6, 8 => 3, 9 => 1]

    # states[r, f, s] = number of parallel states after a solo game of `r - 1`
    #   rounds ((r - 1) * 3 rolls) end on field `f`, with a score of `s`.
    states = fill(0, rounds + 1, 10, max_score + 1)
    states[1, start, 1] = 1
    for round in 1:rounds, score in 1:max_score, field in 1:min(score, 10)
        states[round + 1, field, score + 1] =
            sum([count * states[round, rot(field - roll), score + 1 - field]
                 for (roll, count) in rolls])
    end

    # clean, check and trim
    r = reshape(sum(states[2:end, :, 2:end]; dims = 2), rounds, max_score)
    @assert all(iszero, r[end, 1:20]) "Round count too low, not all games end!"
    @assert sum(r[end, :]) == 27 ^ rounds "Max score too low!"
    return r[CartesianIndex(1, 1):findlast(!iszero, r)]
end

@assert q1(test...) == 739785
@assert q2(test...) == 444356092776315

println("Q1: ", q1(puzzle...))
println("Q2: ", q2(puzzle...))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
# A simple async-safe logger implemented with a macro.

type LogLevel
  name::String
  esc::String
  lvl::Int8
end

log_levels = Dict{Symbol, LogLevel}({
  :error => LogLevel("error", "\e[1;31m", 1),
  :info => LogLevel("info", "\e[1;32m", 2),
  :request => LogLevel("req", "\e[1;34m", 3),
  :debug => LogLevel("debug", "\e[1;35m", 4),
  :default => LogLevel("log", "\e[1;30m", 10)
})

log_level = log_levels[:default]

log_dir = "/tmp/jlogs"

function noop() end

macro log(lev::Symbol, msg...)
  local lvl = get(log_levels, lev, log_levels[:default])
  if lvl.name == "log"
    msg = [lev]
  end
  if lvl.lvl > log_level.lvl
    return :(noop())
  end
  return :( print("\e[1;30m[\e[0m "*$lvl.esc*$lvl.name*" \e[1;30m]\e[0m "*
                  (join($(Expr(:tuple,map(i->:(string($i)),msg)...)), " "))*"\n") )
end

try
  log_level = lowercase(ENV["LOG_LEVEL"])
catch end

if !isdir(log_dir)
  try
    mkdir(log_dir)
  catch ex
   @log error "Could not create log directory" log_dir ":" ex
  end
end

log_buffer_size = 0

function get_log_file()
  log_file_path = log_dir * "/out.log"
  global log_file = open(log_file_path, "w")
end

get_log_file()

function log_to_file(str::String)
  log_buffer_size += sizeof(str)
  write(log_file, str)
end


text:
"""
sprob
"""

function _getstates(states::Vector{Symbol}, x::Vector{Symbol})::Vector{Int}
    h = Dict(s=>i for (i,s)=enumerate(states))
    [h[s] for s = x]
end

function sprob(Q::Matrix{Tv}; maxiter=5000, steps=20, rtol=1.0e-6) where Tv
    _sprob(Q, maxiter, steps, rtol)
end

function sprob(Q::SparseCSC{Tv,Ti}; maxiter=5000, steps=20, rtol=1.0e-6) where {Tv,Ti}
    _sprob(Q, maxiter, steps, rtol)
end

function sprob(Q::SparseMatrixCSC{Tv,Ti}; maxiter=5000, steps=20, rtol=1.0e-6) where {Tv,Ti}
    _sprob(Q, maxiter, steps, rtol)
end

function sprob(m::CTMCModel{Tv}; states = nothing, maxiter=5000, steps=20, rtol=1.0e-6) where Tv
    if states == nothing
        _sprob(m.Q, maxiter, steps, rtol)
    else
        s = _getstates(m.states, states)
        _sprob(m.Q, maxiter, steps, rtol)[s]
    end
end

function _sprob(Q::Matrix{Tv}, maxiter, steps, rtol)::Vector{Tv} where {Tv<:Number}
    gth(Q)
end

function _sprob(Q::SparseCSC{Tv,Ti}, maxiter, steps, rtol)::Vector{Tv} where {Tv<:Number,Ti}
    x, conv, iter, rerror = stgs(Q, maxiter=maxiter, steps=steps, rtol=rtol)
    if conv == false
        @warn "GS did not converge", iter, rerror
    end
    x
end

function _sprob(Q::SparseMatrixCSC{Tv,Ti}, maxiter, steps, rtol)::Vector{Tv} where {Tv<:Number,Ti}
    x, conv, iter, rerror = stgs(Q, maxiter=maxiter, steps=steps, rtol=rtol)
    if conv == false
        @warn "GS did not converge", iter, rerror
    end
    x
end

##

function _sprobsen(Q::Matrix{Tv}, pis::Vector{Tv}, b::Vector{Tv}, maxiter, steps, rtol)::Vector{Tv} where {Tv<:Number}
    stsen(Q, pis, b)
end

function _sprobsen(Q::SparseCSC{Tv,Ti}, pis::Vector{Tv}, b::Vector{Tv}, maxiter, steps, rtol)::Vector{Tv} where {Tv<:Number,Ti}
    s, conv, iter, rerror = stsengs(Q, pis, b, maxiter=maxiter, steps=steps, rtol=rtol)
    if conv == false
        @warn "GSsen did not converge", iter, rerror
    end
    s
end

function _sprobsen(Q::SparseMatrixCSC{Tv,Ti}, pis::Vector{Tv}, b::Vector{Tv}, maxiter, steps, rtol)::Vector{Tv} where {Tv<:Number,Ti}
    s, conv, iter, rerror = stsengs(Q, pis, b, maxiter=maxiter, steps=steps, rtol=rtol)
    if conv == false
        @warn "GSsen did not converge", iter, rerror
    end
    s
end

##

mutable struct SymbolicCTMCProbExpression{Tv} <: AbstractVectorSymbolic{Tv}
    params::Set{Symbol}
    op::Symbol
    Q::AbstractMatrixSymbolic{Tv}
    options::Dict{Symbol,Any}
    dim::Int
end

function _toexpr(x::SymbolicCTMCProbExpression)
    Expr(:call, x.op, _toexpr(x.Q))
end

function Base.show(io::IO, x::SymbolicCTMCProbExpression{Tv}) where Tv
    Base.show(io, "SymbolicCTMCProbExpression $(objectid(x))")
end

##

function sprob(Q::AbstractMatrixSymbolic{Tv}; maxiter=5000, steps=20, rtol=1.0e-6) where Tv
    _sprob(Q, maxiter, steps, rtol)
end

function _sprob(Q::Matrix{<:AbstractSymbolic{Tv}}, maxiter, steps, rtol) where {Tv<:Number}
    _sprob(convert(AbstractMatrixSymbolic{Tv}, Q), maxiter, steps, rtol)
end

function _sprob(Q::SparseCSC{<:AbstractSymbolic{Tv},Ti}, maxiter, steps, rtol) where {Tv<:Number,Ti}
    _sprob(convert(AbstractMatrixSymbolic{Tv}, Q), maxiter, steps, rtol)
end

function _sprob(Q::SparseMatrixCSC{<:AbstractSymbolic{Tv},Ti}, maxiter, steps, rtol) where {Tv<:Number,Ti}
    _sprob(convert(AbstractMatrixSymbolic{Tv}, Q), maxiter, steps, rtol)
end

function _sprob(Q::AbstractMatrixSymbolic{Tv}, maxiter, steps, rtol) where {Tv<:Number}
    SymbolicCTMCProbExpression{Tv}(Q.params, :sprob, Q, Dict(:maxiter=>maxiter, :steps=>steps, :rtol=>rtol), Q.dim[1])
end

"""
seval
"""

function _eval(::Val{:sprob}, f::SymbolicCTMCProbExpression{Tv}, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where {Tv<:Number}
    Q = seval(f.Q, env, cache)
    _sprob(Q, f.options[:maxiter], f.options[:steps], f.options[:rtol])
end

function _eval(::Val{:sprob}, f::SymbolicCTMCProbExpression{Tv}, dvar::Symbol, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where {Tv<:Number}
    pis = seval(f, env, cache)
    Q = seval(f.Q, env, cache)
    dQ = seval(f.Q, dvar, env, cache)
    _sprobsen(Q, pis, dQ' * pis, f.options[:maxiter], f.options[:steps], f.options[:rtol])
end

function _eval(::Val{:sprob}, f::SymbolicCTMCProbExpression{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Vector{Tv} where {Tv<:Number}
    pis = seval(f, env, cache)
    dpis_a = seval(f, dvar[1], env, cache)
    dpis_b = seval(f, dvar[1], env, cache)

    Q = seval(f.Q, env, cache)
    dQ_a = seval(f.Q, dvar[1], env, cache)
    dQ_b = seval(f.Q, dvar[2], env, cache)
    dQ_ab = seval(f.Q, dvar, env, cache)

    _sprobsen(Q, pis, dQ_ab' * pis + dQ_a' * dpis_b + dQ_b' * dpis_a, f.options[:maxiter], f.options[:steps], f.options[:rtol])
end


text:
module TestShow

using Test
using Mimi
import Mimi:
    compdefs, compdef, ComponentId, MimiStruct, ParameterDef

@defcomp X begin
    x = Parameter(index = [time])
    y = Variable(index = [time])
    function run_timestep(p, v, d, t)
        v.y[t] = 1
    end
end

function test_show(obj, expected::AbstractString)
    buf = IOBuffer()
    show(buf, obj)
    @test String(take!(buf)) == expected
end

function showstr(obj)
    buf = IOBuffer()
    show(buf, obj)
    String(take!(buf))
end

compid = ComponentId(TestShow, :something)
test_show(compid, "<ComponentId TestShow.something>")

struct Foo <: MimiStruct
    a::Dict
    b::Int
    c::Float64
end

foo = Foo(Dict((:x=>10, :y=>20)), 4, 44.4)
# Use typeof(Foo) so it works in test mode and using include(), which results in Main.TestShow.Foo.
# test_show(foo, "$(typeof(foo))\n  a: Dict{Symbol,Int64}\n    y => 20\n    x => 10\n  b: 4\n  c: 44.4")

# (:name, :datatype, :dim_names, :description, :unit, :default)
p = ParameterDef(:v1, Float64, [:time], "description string", "Mg C", 101)
# test_show(p, "ParameterDef\n  name: :v1\n  datatype: Float64\n  dim_names: Symbol[:time]\n  description: \"description string\"\n  unit: \"Mg C\"\n  default: 101")

p = ParameterDef(:v1, Float64, [:time], "", "", nothing)
# test_show(p, "ParameterDef\n  name: :v1\n  datatype: Float64\n  dim_names: Symbol[:time]\n  default: nothing")

m = Model()
set_dimension!(m, :time, 2000:2005)
add_comp!(m, X)                         # Original component X
set_param!(m, :X, :x, zeros(6))

expected = """
Model
  md: ModelDef(##anonymous#)
    comp_id: <ComponentId Mimi.##anonymous#>
    variables: OrderedCollections.OrderedDict{Symbol,VariableDef}
    parameters: OrderedCollections.OrderedDict{Symbol,ParameterDef}
    dim_dict: OrderedCollections.OrderedDict{Symbol,Union{Nothing, Dimension}}
      time => [2000, 2001, 2002, 2003, 2004, 2005]
    first: nothing
    last: nothing
    is_uniform: true
    comps_dict: OrderedCollections.OrderedDict{Symbol,AbstractComponentDef}
      X => ComponentDef(X)
        comp_id: <ComponentId TestShow.X>
        variables: OrderedCollections.OrderedDict{Symbol,VariableDef}
          y => VariableDef(y::Number)
              1: time
        parameters: OrderedCollections.OrderedDict{Symbol,ParameterDef}
          x => ParameterDef(x::Number)
              1: time
            default: nothing
        dim_dict: OrderedCollections.OrderedDict{Symbol,Union{Nothing, Dimension}}
          time => nothing
        first: nothing
        last: nothing
        is_uniform: true
      1: ExternalParameterConnection
        comp_name: :X
        param_name: :x
        external_param: :x
    external_params: Dict{Symbol,ModelParameter}
      x => ArrayModelParameter{TimestepArray{FixedTimestep{2000,1,2005},Float64,1}}
        values: TimestepArray{FixedTimestep{2000,1,2005},Float64,1}
            1: 0.0
            2: 0.0
            3: 0.0
            4: 0.0
            ...
            6: 0.0
          1: time
    number_type: Float64
  mi: nothing"""   # ignore (most) whitespace

# Quote regex special characters
# Modified from https://github.com/JuliaLang/julia/issues/6124
function quotemeta(s::AbstractString)
    res = replace(s, r"([()[\]{}?*\$.&~#=!<>|:])" => s"\\\1")
    replace(res, "\0" => "\\0")
end

re = quotemeta(expected)

# remove the number in the gensym since it changes each run
output = replace(showstr(m), r"(##anonymous#)\d+" => s"\1")

@test match(Regex(re), output) !== nothing

end # module


text:
mutable struct FITSDataSet
    # metadata
    datasetid::String
    x::Float64

    function FITSDataSet()
        new("DUMMY", -1.0)
    end

    function FITSDataSet(datasetid)
        new(datasetid, 7.0)
    end
end

a = FITSDataSet()
b = FITSDataSet("ALMA")

finale(x) = @async println("Finalizing $(x.datasetid).")
finalizer(finale, a)

#finalizer(a) do x
#    @async println("Finalizing $x.")
#end

println(a)
println(b)

a = b
println(a)


text:
using Test
using ParallelStencil
import ParallelStencil: @reset_parallel_stencil, @is_initialized, @get_package,  @get_numbertype, @get_ndims, SUPPORTED_PACKAGES, PKG_CUDA, PKG_NONE, NUMBERTYPE_NONE, NDIMS_NONE
import ParallelStencil: @require, @symbols, longnameof
import ParallelStencil: checkargs_init, check_already_initialized, set_initialized, is_initialized, check_initialized, set_package, set_numbertype, set_ndims
using ParallelStencil.Exceptions
TEST_PACKAGES = SUPPORTED_PACKAGES
@static if PKG_CUDA in TEST_PACKAGES
    import CUDA
    if !CUDA.functional() TEST_PACKAGES = filter!(x->x≠PKG_CUDA, TEST_PACKAGES) end
end

@static for package in TEST_PACKAGES  eval(:(
    @testset "$(basename(@__FILE__)) (package: $(nameof($package)))" begin
        @testset "1. initialization of ParallelStencil" begin
            @require !@is_initialized()
            @init_parallel_stencil($package, ComplexF32, 3)
            @testset "initialized" begin
                @test @is_initialized()
                @test @get_package() == $package
                @test @get_numbertype() == ComplexF32
                @test @get_ndims() == 3
            end;
            @testset "Data" begin
                @test @isdefined(Data)
                @test length(@symbols($(@__MODULE__), Data)) > 1
                @test Symbol("Number") in @symbols($(@__MODULE__), Data)
                @test Symbol("Array") in @symbols($(@__MODULE__), Data)
                @test Symbol("DeviceArray") in @symbols($(@__MODULE__), Data)
            end;
            @reset_parallel_stencil()
        end;
        @testset "2. Exceptions" begin
            @testset "already initialized" begin
                set_initialized(true)
                set_package(:CUDA)
                set_numbertype(Float64)
                set_ndims(3)
                @require is_initialized()
                @test_throws IncoherentCallError check_already_initialized(:Threads, Float64, 3)
                @test_throws IncoherentCallError check_already_initialized(:CUDA, Float32, 3)
                @test_throws IncoherentCallError check_already_initialized(:CUDA, Float64, 2)
                set_initialized(false)
                set_package(PKG_NONE)
                set_numbertype(NUMBERTYPE_NONE)
                set_ndims(NDIMS_NONE)
            end;
            @testset "arguments" begin
                MyType = Float32
                @test_throws ArgumentError checkargs_init(5, Float64, 3)
                @test_throws ArgumentError checkargs_init($package, MyType, 3)
                @test_throws ArgumentError checkargs_init($package, Float64, Float64)
            end;
            @testset "check_initialized" begin
                @require !is_initialized()
                @test_throws NotInitializedError check_initialized()
            end;
            @reset_parallel_stencil()
         end;
    end;
)) end == nothing || true;


text:
"""
    MEstimation_results(results::Union{NLsolve.SolverResults, Optim.MultivariateOptimizationResults, Optim.UnivariateOptimizationResults},
                        theta::Vector,
                        data::Any,
                        template::Union{objective_function_template, estimating_function_template},
                        regularizer::Function,
                        br::Bool,
                        has_objective::Bool,
                        has_regularizer::Bool,
                        br_method::String)

[Composite type](https://docs.julialang.org/en/v1/manual/types/#Composite-Types-1) for the output of [`fit`](@ref) for an [`objective_function_template`](@ref) or an [`estimating_function_template`](@ref).

Arguments
===
+ `results`: a `Union{NLsolve.SolverResults, Optim.MultivariateOptimizationResults, Optim.UnivariateOptimizationResults}` object holding the optimization results as returned from `Optim.optimize` or `NLsolve.nlsolve`.
+ `theta`: a `Vector{Float64}` with the M-estimates or their reduced-bias version.
+ `data`: an object of [composite type](https://docs.julialang.org/en/v1/manual/types/#Composite-Types-1) with all the data required to compute the objective function or the estimating functions; see [`objective_function`](@ref)` and [`estimating_function`](@ref).
+ `template`: an [`objective_function_template`](@ref) or [`estimating_function_template`](@ref) object.
+ `regularizer`: a function of `theta` and `data` returning either a `Vector{Float64}` of dimension equal to the number of the estimating functions (for [`estimating_function_template`](@ref)) or a `Float64` (for [`objective_function_template`](@ref)). See [`fit`](@ref) for details.
+ `br`: a `Bool`; if `false` then `results` are from the computation of M-estimates, otherwise from the computation of the RBM-estimates.
+ `has_objective`: a `Bool`; if `true` then `template` is an [`objective_function_template`](@ref).
+ `has_regularizer`: a `Bool`; if `true` then a regularizer function has been used during optimization; see [`fit`](@ref).
+ `br_method`: either "implicit_trace" (default) or "explicit_trace"; see [`fit`](@ref).
"""
struct MEstimation_results
    results::Union{NLsolve.SolverResults, Optim.MultivariateOptimizationResults, Optim.UnivariateOptimizationResults}
    theta::Vector
    data::Any
    template::Union{objective_function_template, estimating_function_template}
    regularizer::Function
    br::Bool
    has_objective::Bool
    has_regularizer::Bool
    br_method::String
end

"""
    vcov(results::MEstimation_results)

Compute an estimate of the variance-covariance matrix of the `M`-estimator or its reduced-bias version at `results.theta`, from a [`MEstimation_results`](@ref) object.

Arguments
===
+ `results`: a [`MEstimation_results`](@ref) object.

Result
===
The `length(coef(results))` times `length(coef(results))` estimated variance covariance matrix for the parameters. This matrix is the empirical sandwich variance covariance matrix for M- and RBM-estimators. See, for example, [Stefanski and Boos (2002, expression 10)](https://www.jstor.org/stable/3087324).

"""
function vcov(results::MEstimation_results)
    if (results.has_objective)
        obj_quantities(results.theta, results.data, results.template, false)[1]
    else
        ef_quantities(results.theta, results.data, results.template, false)[1]
    end
end

"""
    tic(results::MEstimation_results)

Compute the Takeuchi Information Criterion at `results.theta`, from a [`MEstimation_results`](@ref) object.

Arguments
===
+ `results`: a [`MEstimation_results`](@ref) object.

Details
===
`nothing` is returned if `results.template` is an [`estimating_function_template`](@ref).
"""
function tic(results::MEstimation_results)
    if (results.has_objective)
        obj = objective_function(results.theta, results.data, results.template, false)
        quants = obj_quantities(results.theta, results.data, results.template, true)
        -2 * (obj + 2 * quants[1])
    end
end

"""
    aic(results::MEstimation_results)

Compute the Akaike Information Criterion at `results.theta`, from a [`MEstimation_results`](@ref) object with an [`objective_function_template`](@ref).

Arguments
===
+ `results`: a [`MEstimation_results`](@ref) object.

Details
===
`nothing` is returned if `results.template` is an [`estimating_function_template`](@ref).
"""
function aic(results::MEstimation_results)
    if (results.has_objective)
        obj = objective_function(results.theta, results.data, results.template, false)
        p = length(results.theta)
        -2 * (obj - p)
    end
end

"""
    coef(results::MEstimation_results)

Extract the parameter estimates from a [`MEstimation_results`](@ref) object.

Arguments
===
+ `results`: a [`MEstimation_results`](@ref) object.

Details
===
`coef(results)` returns `results.theta`
"""
function coef(results::MEstimation_results)
    results.theta
end

"""
    show(io::IO,
         results::MEstimation_results;
         digits::Int64 = 4)

`show` method for `MEstimation_results` objects.

Arguments
===
+ `io`: an `IO` object; see [`show`](@ref) for details.
+ `results`: a [`MEstimation_results`](@ref) object.

[Keyword arguments](https://docs.julialang.org/en/v1/manual/functions/#Keyword-Arguments-1)
===
+ `digits`: an `Int64` indicating the number of digits to display for the various summaries. Default is `4`.

Details
===
If `MEstimation_results.has_objective == true`, then the result of `aic(results)` and `tic(results)` are also printed.
"""
function Base.show(io::IO, results::MEstimation_results;
                   digits::Real = 4)
    theta = results.theta
    p = length(theta)
    v = vcov(results)
    if results.has_objective
        if results.has_regularizer
            println(io,
                    (results.br ? "RBM" : "M") * "-estimation with objective contributions `",
                    results.template.obj_contribution, "` and user-supplied regularizer")
        else
            println(io,
                    (results.br ? "RBM" : "M") * "-estimation with objective contributions `",
                    results.template.obj_contribution, "`")
        end
    else
        if results.has_regularizer
            println(io,
                    (results.br ? "RBM" : "M") * "-estimation with estimating function contributions `",
                    results.template.ef_contribution, "` and user-supplied regularizer")
        else
            println(io,
                    (results.br ? "RBM" : "M") * "-estimation with estimating function contributions `",
                    results.template.ef_contribution, "`")
        end
    end
    if (results.br)
        println(io, "Bias reduction method: ", results.br_method)
    end
    println(io)
    # println("Parameter\tEstimate\tS.E")
    # for i in 1:p
    #     est = theta[i]
    #     std = sqrt(v[i, i])
    #     println("theta[$(i)]", "\t", round(est, digits = digits), "\t\t", round(std, digits = digits))
    # end
    show(io, coeftable(results))
    if results.has_objective
        objfun = objective_function(results.theta, results.data, results.template, results.br)
        if results.br
            print(io, "\nPenalized objetive:\t\t", round(objfun, digits = digits))
        else
            print(io, "\nObjective:\t\t\t", round(objfun, digits = digits))
        end
        print(io, "\nTakeuchi information criterion:\t", round(tic(results), digits = digits))
        print(io, "\nAkaike information criterion:\t", round(aic(results), digits = digits))
        print(io, "\nConverged: ", Optim.converged(results.results))
    else
        estfun = estimating_function(results.theta, results.data, results.template, results.br)
        if results.br
            print(io, "\nAdjusted estimating functions:\t", estfun)
        else
            print(io, "\nEstimating functions:\t", estfun)
        end
        print(io, "\nConverged: ", NLsolve.converged(results.results))
    end
end

"""
    stderror(results::MEstimation_results)

Compute estimated standard errors from a from a [`MEstimation_results`](@ref) object.

Arguments
===
+ `results`: a [`MEstimation_results`](@ref) object.

Details
===
The estimated standard errors are computed as `sqrt.(diag(vcov(results)))`.
"""
function stderror(results::MEstimation_results)
    sqrt.(diag(vcov(results)))
end

"""
    coeftable(results::MEstimation_results;
              level::Real=0.95)

Return a `StatsBase.CoefTable` from a [`MEstimation_results`](@ref) object.

Arguments
===
+ `results`: a [`MEstimation_results`](@ref) object.

[Keyword arguments](https://docs.julialang.org/en/v1/manual/functions/#Keyword-Arguments-1)
===
+ `level`: a `Real` that determines the level of the reported confidence intervals; default is `0.95`; see Details.

Details
===
The reported confidence intervals are Wald-type of nominal level `level`, using quantiles of the standard normal distribution.
"""
function coeftable(results::MEstimation_results; level::Real=0.95)
    cc = coef(results)
    se = stderror(results)
    zz = cc ./ se
    p = 2 * ccdf.(Ref(Normal()), abs.(zz))
    ci = se * quantile(Normal(), (1-level)/2)
    levstr = isinteger(level*100) ? string(Integer(level*100)) : string(level*100)
    CoefTable(hcat(cc, se, zz, p, cc + ci, cc - ci),
              ["Estimate","Std. Error","z value","Pr(>|z|)","Lower $levstr%","Upper $levstr%"],
              ["theta[$i]" for i = 1:length(cc)], 4)
end


text:
abstract type AbstractPenaltyEvaluator <: AbstractNLPEvaluator end

n_variables(ag::AbstractPenaltyEvaluator) = n_variables(ag.inner)
# All constraints moved inside the objective with penalties!
n_constraints(ag::AbstractPenaltyEvaluator) = 0

constraints_type(ag::AbstractPenaltyEvaluator) = :bound

# Penalty evaluators don't have constraints so Lagrangian is objective
has_hessian_lagrangian(nlp::AbstractPenaltyEvaluator) = has_hessian(nlp)

# Getters
get(ag::AbstractPenaltyEvaluator, attr::AbstractNLPAttribute) = get(ag.inner, attr)
get(ag::AbstractPenaltyEvaluator, attr::AbstractVariable) = get(ag.inner, attr)
get(ag::AbstractPenaltyEvaluator, attr::PS.AbstractNetworkAttribute) = get(ag.inner, attr)

# Setters
function setvalues!(ag::AbstractPenaltyEvaluator, attr::PS.AbstractNetworkValues, values)
    setvalues!(ag.inner, attr, values)
end

# Initial position
initial(ag::AbstractPenaltyEvaluator) = initial(ag.inner)

# Bounds
bounds(ag::AbstractPenaltyEvaluator, ::Variables) = bounds(ag.inner, Variables())
bounds(ag::AbstractPenaltyEvaluator, ::Constraints) = Float64[], Float64[]

# based objective
objective_original(ag::AbstractPenaltyEvaluator, u) = objective(ag.inner, u)

function constraint!(ag::AbstractPenaltyEvaluator, cons, u)
    @assert length(cons) == 0
    return
end

function jacobian!(ag::AbstractPenaltyEvaluator, jac, u)
    @assert length(jac) == 0
    return
end

jacobian_structure(ag::AbstractPenaltyEvaluator) = (Int[], Int[])
function jacobian_structure!(ag::AbstractPenaltyEvaluator, rows, cols)
    @assert length(rows) == length(cols) == 0
end

function hessian_structure(ag::AbstractPenaltyEvaluator)
    return hessian_structure(ag.inner)
end

primal_infeasibility!(ag::AbstractPenaltyEvaluator, cons, u) = primal_infeasibility!(ag.inner, cons, u)
primal_infeasibility(ag::AbstractPenaltyEvaluator, u) = primal_infeasibility(ag.inner, u)


text:
@testset verbose=true "Test Signatures" begin
    using JavaCall.Signatures
    using JavaCall.JNI

    OBJECT_SYM = Symbol("java.lang.Object")
    STRING_SYM = Symbol("java.lang.String")

    @testset "Primitive from Types" begin
        @test signature(jboolean) == "Z"
        @test signature(jbyte) == "B"
        @test signature(jchar) == "C"
        @test signature(jshort) == "S"
        @test signature(jint) == "I"
        @test signature(jlong) == "J"
        @test signature(jfloat) == "F"
        @test signature(jdouble) == "D"
        @test signature(jvoid) == "V"
    end

    @testset "Primitive Arrays from Types" begin
        @test signature(Vector{jboolean}) == "[Z"
        @test signature(Matrix{jbyte}) == "[[B"
        @test signature(Array{jvoid, 4}) == "[[[[V"
        @test signature(Array{jint, 0}) == "I"
    end

    @testset "Objects from Symbols" begin
        @test signature(OBJECT_SYM) == "Ljava/lang/Object;"
        @test signature(Vector{OBJECT_SYM}) == "[Ljava/lang/Object;"
        @test signature(Array{OBJECT_SYM, 4}) == "[[[[Ljava/lang/Object;"
    end

    @testset "Primitive Method Signatures" begin
        @test signature(MethodSignature(jboolean, [])) == "()Z"
        @test signature(MethodSignature(jchar, [jint, jdouble])) == "(ID)C"
        @test signature(MethodSignature(jshort, [jlong, jfloat, Vector{jdouble}])) == "(JF[D)S"
    end

    @testset "Object Method Signatures" begin
        @test signature(MethodSignature(jvoid, [OBJECT_SYM])) == "(Ljava/lang/Object;)V"
        @test signature(MethodSignature(jvoid, [Vector{OBJECT_SYM}])) == "([Ljava/lang/Object;)V"
        @test signature(MethodSignature(STRING_SYM, [OBJECT_SYM, Vector{STRING_SYM}])) ==
            "(Ljava/lang/Object;[Ljava/lang/String;)Ljava/lang/String;"
    end
end


text:
struct XCBButtonCode{T}
    XCBButtonCode(v::Val{T}) where {T} = new{v}()
end

# xcb_button_release_event_t is a type alias for xcb_button_press_event_t
# so there is no need to define methods for xcb_button_release_event_t

XCBButtonCode(mouse_event::xcb_button_press_event_t) = XCBButtonCode(Val(Int(mouse_event.detail)))

button(::XCBButtonCode{Val(1)}) = ButtonLeft()
button(::XCBButtonCode{Val(2)}) = ButtonMiddle()
button(::XCBButtonCode{Val(3)}) = ButtonRight()
button(::XCBButtonCode{Val(4)}) = ButtonScrollUp()
button(::XCBButtonCode{Val(5)}) = ButtonScrollDown()

MouseState(mouse_event::xcb_button_press_event_t) = MouseState((mouse_event.state .| [XCB_BUTTON_MASK_1, XCB_BUTTON_MASK_2, XCB_BUTTON_MASK_3, XCB_BUTTON_MASK_4, XCB_BUTTON_MASK_5, XCB_BUTTON_MASK_ANY] .== mouse_event.state)...)

MouseEvent(mouse_event::xcb_button_press_event_t) = MouseEvent(button(XCBButtonCode(mouse_event)), MouseState(mouse_event), response_type(mouse_event) == XCB_BUTTON_PRESS ? ButtonPressed() : ButtonReleased())

function KeyContext(key_event::xcb_key_press_event_t)
    state = key_event.state
    KeyContext((state .| [XCB_MOD_MASK_2, XCB_MOD_MASK_LOCK] .== state)...)
end

function KeyModifierState(key_event::xcb_key_press_event_t)
    state = key_event.state
    KeyModifierState((state .| [XCB_MOD_MASK_SHIFT, XCB_MOD_MASK_CONTROL, XCB_MOD_MASK_1, XCB_MOD_MASK_4] .== state)...)
end


text:
function cumulative_probability(obj::RealDistribution, arg0::jdouble)
    return jcall(obj, "cumulativeProbability", jdouble, (jdouble,), arg0)
end

function density(obj::RealDistribution, arg0::jdouble)
    return jcall(obj, "density", jdouble, (jdouble,), arg0)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_numerical_mean(obj::RealDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::RealDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

function get_support_lower_bound(obj::RealDistribution)
    return jcall(obj, "getSupportLowerBound", jdouble, ())
end

function get_support_upper_bound(obj::RealDistribution)
    return jcall(obj, "getSupportUpperBound", jdouble, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function inverse_cumulative_probability(obj::AbstractRealDistribution, arg0::jdouble)
    return jcall(obj, "inverseCumulativeProbability", jdouble, (jdouble,), arg0)
end

function is_support_connected(obj::RealDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end

function log_density(obj::AbstractRealDistribution, arg0::jdouble)
    return jcall(obj, "logDensity", jdouble, (jdouble,), arg0)
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function probability(obj::AbstractRealDistribution, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "probability", jdouble, (jdouble, jdouble), arg0, arg1)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end