text:
function fillGlobalFromShmem(testArrInn,resShmem)
    for z in 1:34
        testArrInn[threadIdxX()+1,threadIdxY()+1,z ]=resShmem[threadIdxX()+1,threadIdxY()+1,z ]
    end

        testArrInn[1,threadIdxX()+1,threadIdxY()+1]=  resShmem[1,threadIdxX()+1,threadIdxY()+1]
        testArrInn[34,threadIdxX()+1,threadIdxY()+1]=  resShmem[34,threadIdxX()+1,threadIdxY()+1]
        testArrInn[threadIdxX()+1,1,threadIdxY()+1]=  resShmem[threadIdxX()+1,1,threadIdxY()+1]
        testArrInn[threadIdxX()+1,34,threadIdxY()+1]=  resShmem[threadIdxX()+1,34,threadIdxY()+1]

end

function getIndiciesWithTrue(arr)
    indicies = CartesianIndices(arr)
    return filter(ind-> arr[ind] ,indicies)

end

# function fillGlobalFromShmem(testArrInn,resShmem)
#     for z in 1:34
#         testArrInn[threadIdxX()+1,threadIdxY()+1,z ]=resShmem[threadIdxX()+1,threadIdxY()+1,z ]
#         sync_threads()
#     end

#     for z in 1:32
#         testArrInn[threadIdxX()+1,threadIdxY()+2,z+1 ]=  resShmem[threadIdxX()+1,threadIdxY()+2,z+1 ]
#         sync_threads()
#         testArrInn[threadIdxX()+1,threadIdxY(),z+1 ]=  resShmem[threadIdxX()+1,threadIdxY()+2,z+1 ]
#         sync_threads()
#         testArrInn[threadIdxX()+2,threadIdxY()+1,z+1 ]=  resShmem[threadIdxX()+2,threadIdxY()+1,z+1 ]
#         sync_threads()
#         testArrInn[threadIdxX(),threadIdxY()+1,z+1 ]=  resShmem[threadIdxX(),threadIdxY()+1,z+1 ]
#         sync_threads()
#     end

# end

# # @testset "processMaskData" begin

#     ####### first it should be zeros
#     testArrIn = CUDA.zeros(Bool,34,34,34);
#     testArrOut = CUDA.zeros(Bool,34,34,34);

#     function testKernprocessMaskData(testArrInn,testArrOut)
#         resShmem =  @cuStaticSharedMem(Bool,(34,34,34))#+2 in order to get the one padding
#         clearMainShmem(resShmem)
#         blockBeginingX=1
#         blockBeginingY=1
#         blockBeginingZ=1
#         isMaskFull= zeros(MVector{1,Bool})
#         isMaskEmpty= ones(MVector{1,Bool})
#         #here we will store in registers data uploaded from mask for later verification wheather we should send it or not
#         locArr= UInt32(0)

#         @unroll for zIter in UInt8(1):UInt8(32)# most outer loop is responsible for z dimension
#            processMaskData( testArrInn[(blockBeginingX+threadIdxX()),(blockBeginingY +threadIdxY()),(blockBeginingZ+zIter)], zIter, resShmem)
#         end
#         fillGlobalFromShmem(testArrOut,resShmem)

#         return
#     end
#     @cuda threads=(32,32) blocks=1 testKernprocessMaskData(testArrIn,testArrOut)
#     @test  sum(testArrOut)==0

#     ########### now  we get one arbitrary point it should lead to 6 in output
#     testArrIn = CUDA.zeros(Bool,32,32,32);
#     testArrOut = CUDA.zeros(Bool,34,34,34);

#     testArrIn[5,5,5]=true
#     function testKernprocessMaskDataB(testArrInn,testArrOut)
#         resShmem =  @cuStaticSharedMem(Bool,(34,34,34))#+2 in order to get the one padding
#         clearMainShmem(resShmem)
#         blockBeginingX=0
#         blockBeginingY=0
#         blockBeginingZ=0
#         isMaskFull= zeros(MVector{1,Bool})
#         isMaskEmpty= ones(MVector{1,Bool})
#         #here we will store in registers data uploaded from mask for later verification wheather we should send it or not
#         locArr= UInt32(0)

#         @unroll for zIter in UInt8(1):UInt8(32)# most outer loop is responsible for z dimension
#            processMaskData( testArrInn[(blockBeginingX+threadIdxX()),(blockBeginingY +threadIdxY()),(blockBeginingZ+zIter)], zIter, resShmem)
#         end
#         sync_threads()
#         fillGlobalFromShmem(testArrOut,resShmem)

#         return
#     end
#     @cuda threads=(32,32) blocks=1 testKernprocessMaskDataB(testArrIn,testArrOut)

#     @test  sum(testArrOut)==6
#     @test  testArrOut[5,5,5]==false
#     @test  testArrOut[7,6,6]==true
#     @test  testArrOut[4+1,5+1,5+1]==true
#     @test  testArrOut[5+1,4+1,5+1]==true
#     @test  testArrOut[5+1,6+1,5+1]==true
#     @test  testArrOut[5+1,5+1,4+1]==true
#     @test  testArrOut[5+1,5+1,6+1]==true

#     getIndiciesWithTrue(testArrOut)

#     ########### checking corner cases
#     testArrIn = CUDA.zeros(Bool,32,32,32);
#     testArrOut = CUDA.zeros(Bool,34,34,34);
#     testLocArr = CUDA.zeros(UInt32,32,32,32);

#     testArrIn[1,1,1]=true
#     function testKernprocessMaskDataBB(testArrInn,testArrOut)
#         resShmem =  @cuStaticSharedMem(Bool,(34,34,34))#+2 in order to get the one padding
#         clearMainShmem(resShmem)
#         blockBeginingX=UInt8(0)
#         blockBeginingY=UInt8(0)
#         blockBeginingZ=UInt8(0)
#         isMaskFull= false
#         isMaskEmpty= false
#         #here we will store in registers data uploaded from mask for later verification wheather we should send it or not
#         locArr= UInt32(0)

#         @unroll for zIter in UInt8(1):UInt8(32)# most outer loop is responsible for z dimension
#             #locArr|= true << zIter
#             processMaskData( testArrInn[(blockBeginingX+threadIdxX()),(blockBeginingY +threadIdxY()),(blockBeginingZ+zIter)], zIter, resShmem)
#             # testLocArr[(blockBeginingX+threadIdxX()),(blockBeginingY +threadIdxY()),(blockBeginingZ+zIter)]=locArr
#         end
#         sync_threads()
#         fillGlobalFromShmem(testArrOut,resShmem)

#         return
#     end
#     @cuda threads=(32,32) blocks=1 testKernprocessMaskDataBB(testArrIn,testArrOut)
#     @test  sum(testArrOut)==6

#     @test  testArrOut[1,1,1]==false
#     @test  testArrOut[2,2,1]==true
#     #locArr|= true << zIter

#     CUDA.reclaim()# just to destroy from gpu our dummy data

# ### stress test with ones data

# testArrIn = CUDA.ones(Bool,34,34,34);
# testArrOut = CUDA.zeros(Bool,34,34,34);

# function testKernprocessMaskDataB(testArrInn,testArrOut)
#     resShmem =  @cuStaticSharedMem(Bool,(34,34,34))#+2 in order to get the one padding
#     clearMainShmem(resShmem)
#     blockBeginingX=0
#     blockBeginingY=0
#     blockBeginingZ=0

#     #here we will store in registers data uploaded from mask for later verification wheather we should send it or not
#     locArr= UInt32(0)

#     @unroll for zIter in UInt8(1):UInt8(32)# most outer loop is responsible for z dimension
#        local locBool = testArrInn[(blockBeginingX+threadIdxX()),(blockBeginingY +threadIdxY()),(blockBeginingZ+zIter)]
#        # locArr|= locBool << zIter
#        processMaskData( locBool, zIter, resShmem)
#        CUDA.unsafe_free!(locBool)
#     end
#     sync_threads()
#     #CUDA.@cuprint "locArr $(Int64(locArr)) \n"

#     fillGlobalFromShmem(testArrOut,resShmem)

#     return
# end
# @cuda threads=(32,32) blocks=1 testKernprocessMaskDataB(testArrIn,testArrOut)
# @test  (length(testArrOut)-(sum(testArrOut)+(4*34)+(8*32)  )) ==0

# using  Test, Revise,CUDA

# includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\utils\\CUDAGpuUtils.jl")
# includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\distanceMetrics\\Housdorff\\mainHouseDorffKernel\\HFUtils.jl")
# includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\distanceMetrics\\Housdorff\\mainHouseDorffKernel\\ProcessMainData.jl")
# includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\GPUtestUtils.jl")

# using ..CUDAGpuUtils, ..HFUtils,Cthulhu

# using ..CUDAGpuUtils, ..HFUtils,Cthulhu

# @testset "executeDataIterFirstPass" begin

#     testArrInCPU= falses(60,60,60)
#     testArrInCPU[CartesianIndex(1,1,1)]= true
#     testArrInCPU[CartesianIndex(5,5,5)]= true
#     testArrIn = CuArray(testArrInCPU);
#     referenceArray= CUDA.ones(Bool,32,32,32);
#     resArray = CUDA.zeros(UInt16,32,32,32);
#     resArraysCounter=CUDA.zeros(Int32,1);
#     blockBeginingX,blockBeginingY,blockBeginingZ =UInt8(0),UInt8(0),UInt8(0)

# function testKernelForExecuteMainPass(analyzedArr, refAray,blockBeginingX,blockBeginingY,blockBeginingZ,resArray,resArraysCounter)
#         resShmem =  @cuStaticSharedMem(Bool,(34,34,34))#+2 in order to get the one padding
#         ProcessMainData.executeDataIterFirstPass(analyzedArr, refAray,blockBeginingX,blockBeginingY,blockBeginingZ,resShmem,resArray,resArraysCounter)
#     return
# end

# @cuda threads=(16,16) blocks=1 testKernelForExecuteMainPass(testArrIn,referenceArray,blockBeginingX,blockBeginingY,blockBeginingZ ,resArray,resArraysCounter)
# @test resArraysCounter[1]==9

# @test resArray[5,5,5]==false
# @test resArray[5,6,5]==true
# @test resArray[6,5,5]==true
# @test resArray[4,5,5]==true
# @test resArray[5,4,5]==true
# @test resArray[5,6,5]==true
# @test resArray[5,5,4]==true
# @test resArray[5,5,6]==true

# # @test resArray[1,1,2]==true
# # @test resArray[1,2,1]==true
# # @test resArray[2,1,1]==true

# end

# # @testset "processMaskData" begin

# #     testArrIn = CUDA.ones(Bool,32,32,32);
# #     referenceArray= CUDA.ones(Bool,32,32,32);
# #     resArray = CUDA.zeros(UInt16,32,32,32);
# #     resArraysCounter=CUDA.zeros(Int32,1);

# #     function testprocessMaskData(testArrInn::CuDeviceArray{Bool, 3, 1},resArray,referenceArray,resArraysCounter)
# #     blockBeginingX,blockBeginingY,blockBeginingZ =UInt8(0),UInt8(0),UInt8(0)

# #     resShmem =  @cuStaticSharedMem(Bool,(34,34,34))#+2 in order to get the one padding
# #     locArr = Int32(0)
# #     isMaskFull= false
# #     isMaskEmpty = false
# #     #locArr.x |= true << UInt8(2)

# #     @unroll for zIter::UInt8 in UInt8(1):UInt8(32)# most outer loop is responsible for z dimension
# #         locBool::Bool = @inbounds testArrInn[(blockBeginingX+threadIdxX()),(blockBeginingY +threadIdxY()),(blockBeginingZ+zIter)]
# #         locArr|= locBool << zIter
# #         ProcessMainData.processMaskData( locBool, zIter, resShmem)
# #     end#for
# #    sync_threads() #we should have in resShmem what we need

# #     @unroll for zIter::UInt8 in UInt8(1):UInt8(32)# most outer loop is responsible for z dimension
# #         locVal::Bool = @inbounds  (locArr>>zIter & Int32(1))==Int32(1)
# #         shmemVal::Bool = @inbounds resShmem[threadIdxX()+1,threadIdxY()+1,zIter+1]
# #         locValOrShmem = (locVal | shmemVal)
# #         isMaskFull= locValOrShmem & isMaskFull
# #         isMaskEmpty = ~locValOrShmem & isMaskEmpty

# #        if(!locVal && shmemVal)
# #          # setting value in global memory
# #          @inbounds  testArrInn[(blockBeginingX+threadIdxX()),(blockBeginingY +threadIdxY()),(blockBeginingZ+zIter)]= true
# #          # if we are here we have some voxel that was false in a primary mask and is becoming now true - if it is additionaly true in reference we need to add it to result
# #          isInReferencaArr::Bool= @inbounds referenceArray[(blockBeginingX+threadIdxX()),(blockBeginingY +threadIdxY()),(blockBeginingZ+zIter)]
# #          if(isInReferencaArr)
# #             @inbounds  resArray[(blockBeginingX+threadIdxX()),(blockBeginingY +threadIdxY()),(blockBeginingZ+zIter)]=UInt16(1)
# #             #CUDA.atomic_inc!(pointer(resArraysCounter), Int32(1))

# #             atomicallyAddOneInt(resArraysCounter)
# #          end#if
# #        end#if

# #          end#for

# #   # IterToValidate()

# #     return

# # end#testprocessMaskData

# # @cuda threads=(32,32) blocks=1 testprocessMaskData(testArrIn,resArray,referenceArray,resArraysCounter)
# # resArraysCounter

# # end

# # # #   nv-nsight-cu-cli --mode=launch julia
# # # using  Test, Revise

# # # includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\aPrfofiling\\profilingProcessMaskData.jl")

# # # CUDA.@profile wrapForProfile()

# # # using  Test, Revise

# # # includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\utils\\CUDAGpuUtils.jl")
# # # includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\distanceMetrics\\Housdorff\\mainHouseDorffKernel\\HFUtils.jl")
# # # includet("C:\\GitHub\\GitHub\\NuclearMedEval\\src\\distanceMetrics\\Housdorff\\mainHouseDorffKernel\\ProcessMainData.jl")
# # # includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\GPUtestUtils.jl")

# # # using ..HFUtils, ..ProcessMainData,CUDA,..CUDAGpuUtils,StaticArrays

# # # function wrapForProfile()

# # #     testArrIn = CUDA.ones(Bool,34,34,34);
# # #     testArrOut = CUDA.zeros(Bool,34,34,34);

# # #     function testKernprocessMaskDataB(testArrInn,testArrOut)
# # #         resShmem =  @cuStaticSharedMem(Bool,(34,34,34))#+2 in order to get the one padding
# # #         clearMainShmem(resShmem)
# # #         blockBeginingX=UInt8(0)
# # #         blockBeginingY=UInt8(0)
# # #         blockBeginingZ=UInt8(0)
# # #         isMaskFull= zeros(MVector{1,Bool})
# # #         isMaskEmpty= ones(MVector{1,Bool})
# # #         #here we will store in registers data uploaded from mask for later verification wheather we should send it or not
# # #         locArr= zeros(MVector{32,Bool})

# # #         @unroll for zIter in UInt8(1):UInt8(32)# most outer loop is responsible for z dimension
# # #            processMaskData( testArrInn[(blockBeginingX+threadIdxX()),(blockBeginingY +threadIdxY()),(blockBeginingZ+zIter)], zIter, resShmem,locArr)
# # #         end
# # #         sync_threads()

# # #         fillGlobalFromShmem(testArrOut,resShmem)

# # #         return
# # #     end
# # #     @cuda threads=(32,32) blocks=1 testKernprocessMaskDataB(testArrIn,testArrOut)
# # # end
# # # CUDA.@profile wrapForProfile()

# # #     CUDA.reclaim()# just to destroy from gpu our dummy data

# # # end # processMaskData


text:
using Base.Test

testdir = "tests/"
tests = filter(x -> ismatch(r"^test_.*\.jl$", x), readdir(testdir))

@testset "All tests" begin
    for t in tests
        println("########## ", t, " ##########")
        include(string(testdir, t))
    end
end

ext: jl


text:
@testset "Paginate" begin
    gotd()

    write(joinpath(td, "config.md"), "@def aa = 5")
    mkpath(joinpath(td, "_css"))
    mkpath(joinpath(td, "_layout"))
    write(joinpath(td, "_layout", "head.html"), "HEAD\n")
    write(joinpath(td, "_layout", "foot.html"), "\nFOOT\n")
    write(joinpath(td, "_layout", "page_foot.html"), "\nPG_FOOT\n")
    write(joinpath(td, "index.md"), raw"""
        @def el = ("a", "b", "c", "d")
        ~~~
~~~
        {{paginate el 2}}
        ~~~
~~~
        """)
    write(joinpath(td, "foo.md"), raw"""
        @def a = ["
Item $i
" for i in 1:10]
        Some content
        ~~~
~~~
        {{paginate a 4}}
        ~~~
~~~
        """)
    serve(single=true)
    # expected outputs for index
    @test isfile(joinpath("__site", "index.html"))
    @test isfile(joinpath("__site", "1", "index.html"))
    @test isfile(joinpath("__site", "2", "index.html"))
    @test isfile(joinpath("__site", "foo", "index.html"))
    @test isfile(joinpath("__site", "foo", "1", "index.html"))
    @test isfile(joinpath("__site", "foo", "2", "index.html"))
    @test isfile(joinpath("__site", "foo", "3", "index.html"))
    # expected content
    @test read(joinpath("__site", "index.html"), String) // """
        HEAD

 ab

        PG_FOOT

        FOOT"""
    @test read(joinpath("__site", "2", "index.html"), String) // """
        HEAD

 cd

        PG_FOOT

        FOOT"""
    @test read(joinpath("__site", "foo", "1", "index.html"), String) // """
        HEAD

Some content

Item 1

Item 2

Item 3

Item 4

        PG_FOOT

        FOOT"""
    @test read(joinpath("__site", "foo", "3", "index.html"), String) // """
        HEAD

Some content

Item 9

Item 10

        PG_FOOT

        FOOT"""

    # WARNINGS
    # (see warnings)

    # ERRORS
    write(joinpath(td, "foo.md"), raw"""
        @def a = ["
Item $i
" for i in 1:10]
        Some content
        ~~~
~~~
        {{paginate a}}
        ~~~
~~~
        """)
    @test_throws Franklin.HTMLFunctionError serve(single=true)

end


text:
module mocylpull14
using FinEtools
using Compat.Test
function test()
    E1=1.0;
    nu23=0.19;
    rin=1.;
    rex =1.2;
    Length = 1*rex
    tolerance=rin/1000.

    MR = DeforModelRed2DAxisymm

    fens,fes = Q4block(rex-rin,Length,5,20);
    fens.xyz[:, 1] = fens.xyz[:, 1] .+ rin
    bdryfes = meshboundary(fes);

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),2))
    l1 =selectnode(fens; box=[0 rex 0 0], inflate = tolerance)
    setebc!(u,l1,true, 2, 0.0)
    l1 =selectnode(fens; box=[0 rex Length Length], inflate = tolerance)
    setebc!(u,l1,true, 2, 0.0)
    applyebc!(u)
    numberdofs!(u)
    @test u.nfreedofs == 240

    material=MatDeforElastIso(MR, 00.0, E1, nu23, 0.0)
    println("success? ")
    femm = FEMMDeforLinear(MR, IntegData(fes, GaussRule(2, 2), true), material, true)
    println("failure? ")
    femm = FEMMDeforLinear(MR, IntegData(fes, GaussRule(2, 2), true), material)

    true
end
end
using .mocylpull14
mocylpull14.test()


text:
"""
    b̃p(p; k=1)

Closure for the circular wave function in polar coordinates.

The closure captures the values of the mode number ``p`` as well as the
wavenumber ``k`` and does not include any normalization.
The function can then be evaluated at any ``(r,θ)`` point.

The circular waves in polar coordinates are defined by
```math
b̃_p := (r,θ) ↦ J_{p}(kr) e^{ı p θ}, \\qquad ∀ p ∈ \\mathbb{Z}.
```
"""
b̃p(p; k=1) = (r,θ) -> besselj(p, k*r) * exp(im * p * θ)

"""
    βp(p; k=1)

Computation of the normalization constant of the circular waves.

The circular wave is defined by [`b̃p`](@ref) and the normalization is done
using the ``k``-weighted ``H^1`` norm.

By definition
```math
    \\| b̃_p \\|_{H^{1}}^2 = \\| b̃_p \\|_{L^{2}}^2 + k^{-2} \\| \\nabla b̃_p \\|_{L^{2}}^2.
```
Using integration by parts,
```math
    \\| \\nabla b̃_p \\|_{L^{2}}^2 = k^2 \\| b̃_p \\|_{L^{2}}^2 + (\\partial_n b̃_p , b̃_p).
```
On the one hand
```math
    \\| b̃_p \\|_{L^{2}}^2 = \\pi (J_p^2(k) - J_{p-1}(k)J_{p+1}(k)).
```
On the other hand
```math
    (\\partial_n b̃_p , b̃_p) = \\pi k (J_{p-1}(k) - J_{p+1}(k))J_p(k).
```
Hence
```math
    \\| b̃_p \\|_{H^{1}}^2 = 2\\pi (J_p^2(k) - J_{p-1}(k)J_{p+1}(k))
    + \\pi k^{-1} (J_{p-1}(k) - J_{p+1}(k))J_p(k).
```
"""
function βp(p; k=1)
    L2nrm2 = π * (besselj(p,k)^2 - besselj(p-1,k) * besselj(p+1,k))
    H1nrm2 = 2 * L2nrm2 + (π/k) * (besselj(p-1,k) - besselj(p+1,k)) * besselj(p,k)
    return 1 / sqrt(H1nrm2)
end

"""
    bp(p; k=1)

Closure for the normalized circular waves.

The closure captures the values of the mode number ``p`` and the wavenumber
``k``.
The normalization is the default normalization defined in the function
[`βp`](@ref) and is precomputed once for all.
The function can then be evaluated at any ``(r,θ)`` point.

The normalized circular waves are defined by
```math
b_p := β_p b̃_p, \\qquad ∀ p ∈ \\mathbb{Z}.
```
"""
function bp(p; k=1)
    β = βp(p; k=k)
    return (r,θ) -> β * b̃p(p; k=k)(r,θ)
end

"""
    solution_surrogate(U; k=1)

Computes a solution surrogate from a set of coefficients `U`.

The parameter `U` is a vector of coefficients ``u_p`` for ``p`` ranging from
``-P`` to ``P`` (hence of size ``2P+1``).
The solution surrogate is then
```math
\\mathbf{x} ↦ \\sum_{|p| \\leq P} u_p b_{p}(\\mathbf{x}).
```
"""
function solution_surrogate(U; k=1)
    P = (length(U) - 1) // 2
    return (r, θ) -> sum([u * bp(p; k=k)(r, θ) for (p, u) in zip(-P:P, U)])
end

max_stars_repo_name: EmileParolin/STrAW

max_issues_repo_name: EmileParolin/STrAW

max_forks_repo_name: EmileParolin/STrAW


text:
module mushr_pf

include("MotionModels/MotionModels.jl")
include("SensorModels/SensorModels.jl")
include("ParticleFilter/ParticleFilter.jl")

end # module


text:
using CUDNN: CD, cudnnConvolutionForward, CUDNN_CONVOLUTION, CUDNN_CROSS_CORRELATION, cudnnGetConvolutionNdDescriptor

function cudnnConvolutionForward(x,w; cd=nothing)
    nd, pads, strides, upscale, mode, dtype = cd != nothing ? cudnnGetConvolutionNdDescriptor(cd) : (2, (0,0), (1,1), (1,1), CUDNN_CONVOLUTION, Float32)
    pad, stride = pads[1], strides[1]
end

# outputDim = 1 + (inputDim + 2*pad - filterDim)/convolutionStride
###
# cudnnConvolutionForward(x, w, y;
#    padding=c.padding, stride=c.stride, upscale=c.upscale, mode=c.mode,
#                            algorithm=c.algorithm, workSpace=c.workSpace, workSpaceSizeInBytes=c.workSpaceSizeInBytes,
#                                                        alpha=c.alpha, beta=c.beta)
###
function kudnnConvolutionForward(x,w; cd=nothing)
    # x: (W,H,C,N)
    # w: (W,H,C,K)
    # y: (W,H,K,N)
    nd, pads, strides, upscale, mode, dtype = cd != nothing ? cudnnGetConvolutionNdDescriptor(cd) : (2, (0,0), (1,1), (1,1), CUDNN_CONVOLUTION, Float32)
    pad, stride = pads[1], strides[1]
    Wx,Hx,Cx,N = size(x)
    Ww,Hw,Cw,K = size(w)
    @assert Cx==Cw
    @assert mode == CUDNN_CONVOLUTION
    Wy,Hy = floor(Int, 1 + (Int[Wx,Hx] + 2*pad - Int[Ww,Hw]) / stride)
    # Wy,Hy = 1 + div(Int[Wx,Hx] + 2*pad - Int[Ww,Hw], stride)
    y = zeros(dtype, (Wy,Hy,K,N))
    for n in 1:N, k in 1:K, c in 1:Cx
        y[:,:,k,n] += _conv2(x[:,:,c,n], w[:,:,c,k]; pad=pad, stride=stride, xcorr=CUDNN_CONVOLUTION!=mode)
    end
    return y
end

# dw = rot180(xcorr(x,dy))
function kudnnConvolutionBackwardFilter(x,dy,dw; cd=nothing)
    # x:    (Wx,Hx,Cx,N)
    # dy:   (Wy,Hy,K,N)
    # dw:    (Ww,Hw,Cw,K)
    nd, pads, strides, upscale, mode, dtype = cd != nothing ? cudnnGetConvolutionNdDescriptor(cd) : (2, (0,0), (1,1), (1,1), CUDNN_CONVOLUTION, Float32)
    pad, stride = pads[1], strides[1]
    @assert pad==0 && stride==1
    @assert mode == CUDNN_CONVOLUTION
    Wx,Hx,C,Nx = size(x)
    Wy,Hy,K,Ny = size(dy)
    dw2 = zeros(dw)
    for c in 1:C, k in 1:K, n in 1:Ny
        dw2[:,:,c,k] += rot180(_conv2(x[:,:,c,n], dy[:,:,k,n]; pad=pad, stride=stride, xcorr=true))
    end
    return dw2
end

# dx = xcorr(dy, w, 'full')
function kudnnConvolutionBackwardData(dy,w,dx; cd=nothing) # TODO change the order of parameters
    nd, pads, strides, upscale, mode, dtype = cd != nothing ? cudnnGetConvolutionNdDescriptor(cd) : (2, (0,0), (1,1), (1,1), CUDNN_CONVOLUTION, Float32)
    pad, stride = pads[1], strides[1]
    @assert pad==0 && stride==1
    @assert mode == CUDNN_CONVOLUTION
    Wy,Hy,Ky,N = size(dy)
    Ww,Hw,C,Kw = size(w)
    @assert Ky==Kw
    dx2 = zeros(dx)
    for n in 1:N, c in 1:C, k in 1:Kw
        t = conv2(dy[:,:,k,n], rot180(w[:,:,c,k]))
        # t = _conv2(dy[:,:,k,n], w[:,:,c,k]; pad=Ww-1, stride=stride, xcorr=true)
        dx2[:,:,c,n] += t
    end
    return dx2
end


text:
# This file contains code that was formerly part of Julia. License is MIT: https://julialang.org/license

using AbstractFFTs
using AbstractFFTs: Plan
using ChainRulesTestUtils

using LinearAlgebra
using Random
using Test

import Unitful

Random.seed!(1234)

include("testplans.jl")

@testset "rfft sizes" begin
    A = rand(11, 10)
    @test @inferred(AbstractFFTs.rfft_output_size(A, 1)) == (6, 10)
    @test @inferred(AbstractFFTs.rfft_output_size(A, 2)) == (11, 6)
    A1 = rand(6, 10); A2 = rand(11, 6)
    @test @inferred(AbstractFFTs.brfft_output_size(A1, 11, 1)) == (11, 10)
    @test @inferred(AbstractFFTs.brfft_output_size(A2, 10, 2)) == (11, 10)
    @test_throws AssertionError AbstractFFTs.brfft_output_size(A1, 10, 2)
end

@testset "Custom Plan" begin
    # DFT along last dimension, results computed using FFTW
    for (x, fftw_fft) in (
        (collect(1:7),
         [28.0 + 0.0im,
          -3.5 + 7.267824888003178im,
          -3.5 + 2.7911568610884143im,
          -3.5 + 0.7988521603655248im,
          -3.5 - 0.7988521603655248im,
          -3.5 - 2.7911568610884143im,
          -3.5 - 7.267824888003178im]),
        (collect(1:8),
         [36.0 + 0.0im,
          -4.0 + 9.65685424949238im,
          -4.0 + 4.0im,
          -4.0 + 1.6568542494923806im,
          -4.0 + 0.0im,
          -4.0 - 1.6568542494923806im,
          -4.0 - 4.0im,
          -4.0 - 9.65685424949238im]),
        (collect(reshape(1:8, 2, 4)),
         [16.0+0.0im  -4.0+4.0im  -4.0+0.0im  -4.0-4.0im;
          20.0+0.0im  -4.0+4.0im  -4.0+0.0im  -4.0-4.0im]),
        (collect(reshape(1:9, 3, 3)),
         [12.0+0.0im  -4.5+2.598076211353316im  -4.5-2.598076211353316im;
          15.0+0.0im  -4.5+2.598076211353316im  -4.5-2.598076211353316im;
          18.0+0.0im  -4.5+2.598076211353316im  -4.5-2.598076211353316im]),
    )
        # FFT
        dims = ndims(x)
        y = AbstractFFTs.fft(x, dims)
        @test y ≈ fftw_fft
        P = plan_fft(x, dims)
        @test eltype(P) === ComplexF64
        @test P * x ≈ fftw_fft
        @test P \ (P * x) ≈ x

        fftw_bfft = complex.(size(x, dims) .* x)
        @test AbstractFFTs.bfft(y, dims) ≈ fftw_bfft
        P = plan_bfft(x, dims)
        @test P * y ≈ fftw_bfft
        @test P \ (P * y) ≈ y

        fftw_ifft = complex.(x)
        @test AbstractFFTs.ifft(y, dims) ≈ fftw_ifft
        P = plan_ifft(x, dims)
        @test P * y ≈ fftw_ifft
        @test P \ (P * y) ≈ y

        # real FFT
        fftw_rfft = fftw_fft[
            (Colon() for _ in 1:(ndims(fftw_fft) - 1))...,
            1:(size(fftw_fft, ndims(fftw_fft)) ÷ 2 + 1)
        ]
        ry = AbstractFFTs.rfft(x, dims)
        @test ry ≈ fftw_rfft
        P = plan_rfft(x, dims)
        @test eltype(P) === Int
        @test P * x ≈ fftw_rfft
        @test P \ (P * x) ≈ x

        fftw_brfft = complex.(size(x, dims) .* x)
        @test AbstractFFTs.brfft(ry, size(x, dims), dims) ≈ fftw_brfft
        P = plan_brfft(ry, size(x, dims), dims)
        @test P * ry ≈ fftw_brfft
        @test P \ (P * ry) ≈ ry

        fftw_irfft = complex.(x)
        @test AbstractFFTs.irfft(ry, size(x, dims), dims) ≈ fftw_irfft
        P = plan_irfft(ry, size(x, dims), dims)
        @test P * ry ≈ fftw_irfft
        @test P \ (P * ry) ≈ ry
    end
end

@testset "Shift functions" begin
    @test @inferred(AbstractFFTs.fftshift([1 2 3])) == [3 1 2]
    @test @inferred(AbstractFFTs.fftshift([1, 2, 3])) == [3, 1, 2]
    @test @inferred(AbstractFFTs.fftshift([1 2 3; 4 5 6])) == [6 4 5; 3 1 2]

    @test @inferred(AbstractFFTs.fftshift([1 2 3; 4 5 6], 1)) == [4 5 6; 1 2 3]
    @test @inferred(AbstractFFTs.fftshift([1 2 3; 4 5 6], ())) == [1 2 3; 4 5 6]
    @test @inferred(AbstractFFTs.fftshift([1 2 3; 4 5 6], (1,2))) == [6 4 5; 3 1 2]
    @test @inferred(AbstractFFTs.fftshift([1 2 3; 4 5 6], 1:2)) == [6 4 5; 3 1 2]

    @test @inferred(AbstractFFTs.ifftshift([1 2 3])) == [2 3 1]
    @test @inferred(AbstractFFTs.ifftshift([1, 2, 3])) == [2, 3, 1]
    @test @inferred(AbstractFFTs.ifftshift([1 2 3; 4 5 6])) == [5 6 4; 2 3 1]

    @test @inferred(AbstractFFTs.ifftshift([1 2 3; 4 5 6], 1)) == [4 5 6; 1 2 3]
    @test @inferred(AbstractFFTs.ifftshift([1 2 3; 4 5 6], ())) == [1 2 3; 4 5 6]
    @test @inferred(AbstractFFTs.ifftshift([1 2 3; 4 5 6], (1,2))) == [5 6 4; 2 3 1]
    @test @inferred(AbstractFFTs.ifftshift([1 2 3; 4 5 6], 1:2)) == [5 6 4; 2 3 1]
end

@testset "FFT Frequencies" begin
    @test fftfreq(8) isa Frequencies
    @test copy(fftfreq(8)) isa Frequencies

    # N even
    @test fftfreq(8) == [0.0, 0.125, 0.25, 0.375, -0.5, -0.375, -0.25, -0.125]
    @test rfftfreq(8) == [0.0, 0.125, 0.25, 0.375, 0.5]
    @test fftshift(fftfreq(8)) == -0.5:0.125:0.375

    # N odd
    @test fftfreq(5) == [0.0, 0.2, 0.4, -0.4, -0.2]
    @test rfftfreq(5) == [0.0, 0.2, 0.4]
    @test fftshift(fftfreq(5)) == -0.4:0.2:0.4

    # Sampling Frequency
    @test fftfreq(5, 2) == [0.0, 0.4, 0.8, -0.8, -0.4]
    # <:Number type compatibility
    @test eltype(fftfreq(5, ComplexF64(2))) == ComplexF64

    @test_throws ArgumentError Frequencies(12, 10, 1)

    @testset "scaling" begin
        @test fftfreq(4, 1) * 2 === fftfreq(4, 2)
        @test fftfreq(4, 1) .* 2 === fftfreq(4, 2)
        @test 2 * fftfreq(4, 1) === fftfreq(4, 2)
        @test 2 .* fftfreq(4, 1) === fftfreq(4, 2)

        @test fftfreq(4, 1) / 2 === fftfreq(4, 1/2)
        @test fftfreq(4, 1) ./ 2 === fftfreq(4, 1/2)

        @test 2 \ fftfreq(4, 1) === fftfreq(4, 1/2)
        @test 2 .\ fftfreq(4, 1) === fftfreq(4, 1/2)
    end

    @testset "extrema" begin
        function check_extrema(freqs)
            for f in [minimum, maximum, extrema]
                @test f(freqs) == f(collect(freqs)) == f(fftshift(freqs))
            end
        end
        for f in (fftfreq, rfftfreq), n in (8, 9), multiplier in (2, 1/3, -1/7, 1.0*Unitful.mm)
            freqs = f(n, multiplier)
            check_extrema(freqs)
        end
    end
end

@testset "normalization" begin
    # normalization should be inferable even if region is only inferred as ::Any,
    # need to wrap in another function to test this (note that p.region::Any for
    # p::TestPlan)
    f9(p::Plan{T}, sz) where {T} = AbstractFFTs.normalization(real(T), sz, p.region)
    @test @inferred(f9(plan_fft(zeros(10), 1), 10)) == 1/10
end

@testset "ChainRules" begin
    @testset "shift functions" begin
        for x in (randn(3), randn(3, 4), randn(3, 4, 5))
            for dims in ((), 1, 2, (1,2), 1:2)
                any(d > ndims(x) for d in dims) && continue

                # type inference checks of `rrule` fail on old Julia versions
                # for higher-dimensional arrays:
                # https://github.com/JuliaMath/AbstractFFTs.jl/pull/58#issuecomment-916530016
                check_inferred = ndims(x) < 3 || VERSION >= v"1.6"

                test_frule(AbstractFFTs.fftshift, x, dims)
                test_rrule(AbstractFFTs.fftshift, x, dims; check_inferred=check_inferred)

                test_frule(AbstractFFTs.ifftshift, x, dims)
                test_rrule(AbstractFFTs.ifftshift, x, dims; check_inferred=check_inferred)
            end
        end
    end

    @testset "fft" begin
        for x in (randn(3), randn(3, 4), randn(3, 4, 5))
            N = ndims(x)
            complex_x = complex.(x)
            for dims in unique((1, 1:N, N))
                for f in (fft, ifft, bfft)
                    test_frule(f, x, dims)
                    test_rrule(f, x, dims)
                    test_frule(f, complex_x, dims)
                    test_rrule(f, complex_x, dims)
                end

                test_frule(rfft, x, dims)
                test_rrule(rfft, x, dims)

                for f in (irfft, brfft)
                    for d in (2 * size(x, first(dims)) - 1, 2 * size(x, first(dims)) - 2)
                        test_frule(f, x, d, dims)
                        test_rrule(f, x, d, dims)
                        test_frule(f, complex_x, d, dims)
                        test_rrule(f, complex_x, d, dims)
                    end
                end
            end
        end
    end
end


text:
"""
    function copyright()
---
Show copyright information.
"""
function copyright()
    w = 60
    year = Dates.year(now())
    msg = ["",
           "GEAS: Genome Editing Assisted Selection",
           "===========",
           "Copyright © Xijiang Yu, $year",
           "MIT license"]
    @info join(lpad.(msg, w), "\n")
end

function isLaterBeagle(a, b)
    mon = Dict("Jan" => 1, "Feb" => 2, "Mar" => 3,
               "Apr" => 4, "May" => 5, "Jun" => 6,
               "Jul" => 7, "Aug" => 8, "Sep" => 9,
               "Oct" => 10, "Nov" => 11, "Dec" => 12, "mmm" => 0)
    da, db = a[8:9], b[8:9]
    ma, mb = mon[a[10:12]], mon[b[10:12]]
    ya, yb = a[13:14], b[13:14]
    if ya == yb
        if ma == mb
            return da > db ? a : b
        else
            return ma > mb ? a : b
        end
    else
        return ya > yb ? a : b
    end
end

"""
    function update_beagle()
---
Return the latest beagle URL
"""
function update_beagle()
    BeagleURL = "https://faculty.washington.edu/browning/beagle/"
    # Note Beagle is not version controlled by git.
    # It seems that names for new versions are following the pattern like
    # beagle.ddMmmyy.xxx.jar
    # so I just grab all the file names in $BeagleURL
    # There also must be better way to get the URL, like XML?
    beagle = "beagle.00mmm00.xxx.jar"
    web = HTTP.request("GET", BeagleURL)
    txt = split(String(web.body), "\n")
    for line in txt
        r = match(r"beagle\.\d{2}\w{3}\d{2}.{5}jar", line)
        r == nothing || (beagle = isLaterBeagle(beagle, r.match))
    end
    @info "Updating Beagle to $beagle."
    local_beagle = joinpath(bin_dir, "beagle.jar")
    isfile(local_beagle) && rm(local_beagle, force=true)
    download(joinpath(BeagleURL, beagle), local_beagle)

    beagle2vcfURL = "https://faculty.washington.edu/browning/beagle_utilities/beagle2vcf.jar"
    beagle2vcf = joinpath(bin_dir, "beagle2vcf.jar")

    if !isfile(joinpath(bin_dir, beagle2vcf))
        @info "Downloading beagle2vcf.jar"
        download(beagle2vcfURL, beagle2vcf)
    end
end

"""
    function update_plink()
---
Return URL of the latest plink version 1 for Linux x86_64.
"""
function update_plink()
    @info "Updating plink to the latest"
    plinkURL = "http://s3.amazonaws.com/plink1-assets/plink_linux_x86_64_latest.zip"
    cd(joinpath(dat_dir, "test"))
    download(plinkURL, "plink.zip")
    run(`unzip -f plink.zip`)
    cp("plink", joinpath(bin_dir, "plink"), force = true)
    cd(prj_dir)
end

"""
    function update_macs()
---
Clone and compile marcs into `bin`.
"""
function update_macs()
    @info "Updating macs"
    marcs = "https://github.com/gchen98/macs"
    cd(dat_dir)
    if isdir("macs")
        cd("macs")
        run(`git pull --no-rebase`)
        cd("..")
    else
        run(`git clone $marcs`)
    end
    cd("macs")
    run(`g++ -o macs -O3 -Wall simulator.cpp algorithm.cpp datastructures.cpp`)
    run(`g++ -o msformatter -O3 -Wall msformat.cpp`)
    cp("macs", joinpath(bin_dir, "macs"), force = true)
    cp("msformatter", joinpath(bin_dir, "msformatter"), force = true)
    cd(prj_dir)
end

"""
    function get_QMSim()
---
Pull the latest QMSim to package bin directory
"""
function get_qmsim(;force=false)
    url = "http://www.aps.uoguelph.ca/~msargol/qmsim/QMSim_Linux.zip"
    if !isfile("$bin_dir/qmsim") || force
        @info "Fetch `QMSim`"
        tst = joinpath(dat_dir, "test")
        isdir(tst) || mkpath(tst)
        cd(tst)
        download(url, "qmsim.zip")
        run(`unzip -f qmsim.zip`)
        cp("QMSim_Linux/QMSim", joinpath(bin_dir, "qmsim"), force = true)
        cd(prj_dir)
    else
        @info join(["",
                    "QMSim was already downloaded",
                    "Call GEAS.get_qmsim(force=true) to force download"],
                   "\n")
    end
end

"""
    Update()
---
Every time one starts this package, the package will automatically check if `plink` and `beagle.jar` exist.
If not, the package will download the latest version of these two files.
The package will update, or download them on 17th every month anyway.
"""
function Update()
    @info "Updating binaries"
    isdir(bin_dir) || mkdir(bin_dir)
    update_beagle()
    update_plink()
    update_macs()
    get_qmsim()
end


text:
export LogHomotopy

"""
    LogHomotopy(H::AbstractHomotopy)

Reparameterizes the homotopy ``H(x,t)`` as ``Ĥ(x,s) := H(x,e^{-s})``.
"""
struct LogHomotopy{H<:AbstractHomotopy} <: AbstractHomotopy
    homotopy::H
end

Base.size(H::LogHomotopy) = size(H.homotopy)

struct LogHomotopyCache{HC} <: AbstractHomotopyCache
    cache::HC
end

function cache(H::LogHomotopy, x, s)
    t = exp(-s)
    c = cache(H.homotopy, x, t)
    LogHomotopyCache(c)
end

function evaluate!(u, H::LogHomotopy, x, s, c::LogHomotopyCache)
    evaluate!(u, H.homotopy, x, exp(-s), c.cache)
end
function jacobian!(U, H::LogHomotopy, x, s, c::LogHomotopyCache)
    jacobian!(U, H.homotopy, x, exp(-s), c.cache)
end
function dt!(u, H::LogHomotopy, x, s, c::LogHomotopyCache)
    t = exp(-s)
    dt!(u, H.homotopy, x, t, c.cache)
    LinearAlgebra.rmul!(u, -t)
    u
end

function evaluate_and_jacobian!(u, U, H::LogHomotopy, x, s, c::LogHomotopyCache)
    evaluate_and_jacobian!(u, U, H.homotopy, x, exp(-s), c.cache)
end

function jacobian_and_dt!(U, u, H::LogHomotopy, x, s, c::LogHomotopyCache)
    t = exp(-s)
    jacobian_and_dt!(U, u, H.homotopy, x, t, c.cache)
    LinearAlgebra.rmul!(u, -t)

    nothing
end

function evaluate(H::LogHomotopy, x, s, c::LogHomotopyCache)
    evaluate(H.homotopy, x, exp(-s), c.cache)
end
function jacobian(H::LogHomotopy, x, s, c::LogHomotopyCache)
    jacobian(H.homotopy, x, exp(-s), c.cache)
end
function dt(H::LogHomotopy, x, s, c::LogHomotopyCache)
    t = exp(-s)
    u = dt(H.homotopy, x, t, c.cache)
    LinearAlgebra.rmul!(u, -t)
    u
end

basehomotopy(H::LogHomotopy) = H.homotopy


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "LAMMPS"
version = v"1.0.0" # Equivalent to 2020-10-29

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/lammps/lammps.git", "88fd96ec52f86dba4b222623f3a06632a32e42f1")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/lammps/
mkdir build && cd build/
cmake ../cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DBUILD_SHARED_LIBS=ON \
    -DLAMMPS_EXCEPTIONS=ON \
    -DPKG_SNAP=ON
make -j${nproc}
make install

if [[ "${target}" == *mingw* ]]; then
    cp *.dll ${prefix}/bin/
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental=true))
filter!(p -> libc(p) != "musl", platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("liblammps", :liblammps),
    ExecutableProduct("lmp", :lmp),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"8")

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
println("==================================================================================")
println("Beginning execution of script GMHExamples.jl/ode/springmass/scripts/SpringMass1.jl")
println("==================================================================================")

import GeneralizedMetropolisHastings #pre-include to avoid race condition when running in parallel
import GMHExamples

println("Number of parallel processes running: ",nprocs())
println("Use addprocs() in the REPL if you want to run on more processes")

@everywhere using GeneralizedMetropolisHastings
@everywhere using GMHExamples
using PyPlot

#Standard M-H for nproposals == 1
#Generalized M-H for nproposals > 1
nproposals = 20

#MCMC iteration specifications
nburnin = 200
niterations = 1000
ntunerperiod = 20

#Time points to simulate the spring-mass ODE
timepoints = 0.0:0.1:10.0

###Initial conditions for the spring-mass ODE (position and speed)
initialposition = -1.0 #in meters
initialvelocity = 1.0 #in meters/second

###Values of the model parameters (spring stiffness K and mass M)
K = 50.0 #in Newton/meter
M = 10.0 #in kg
lows = [K-K/5,M-M/5]
highs = [K+K/5,M+M/5]

###The variance of the normal noise on the input data
variance = [0.01,0.09]

println("==========================================")
println("Simulation parameters defined successfully")
println("==========================================")

###Create a Spring-Mass model with measurement data and ODE function
m = springmassmodel(timepoints,[initialposition,initialvelocity],[K,M],variance,lows,highs)

###Show the model
println("==========================")
println("Model defined successfully")
println("==========================")
show(m)

###Plot the measurement data (simmulated data + noise)
figure(string("SpringMass1-",nprocs())) ; clf()
subplot(221)
plot(dataindex(m),measurements(m)[:,1];label="location")
plot(dataindex(m),measurements(m)[:,2];label="velocity")
xlabel("Time")
ylabel("Amplitude")
title("Spring-Mass measurement data")
grid("on")
legend(loc="upper right",fancybox="true")

###Create a Metropolis sampler with a Normal proposal density
s = sampler(:mh,:normal,0.1,ones(2))
println("============================")
println("Sampler defined successfully")
println("============================")
show(s)

###Create a tuner that scales the proposal density
t = tuner(:scale,ntunerperiod,0.5,:erf)
println("==========================")
println("Tuner defined successfully")
println("==========================")
show(t)

###Create a Generalized Metropolis-Hastings runner (which will default to Standard MH when nproposals=1)
p = policy(:mh,nproposals;initialize=:prior)
r = runner(p,niterations,nproposals;numburnin = nburnin)
println("===========================")
println("Runner defined successfully")
println("===========================")
show(r)

###Run the MCMC (can take quite a bit of time)
println("=======================")
println("Run the MCMC simulation")
println("=======================")
@time c = run!(r,m,s,t)
println("=========================")
println("Completed MCMC simulation")
println("=========================")

###Show the result of the simulations
show(c)

nparas = numparas(m)
meanparamvals = mean(samples(c),2)
stdparamvals = std(samples(c),2)

println("Results of the MCMC simulation:")
println(" mean K: ",meanparamvals[1])
println(" mean M: ",meanparamvals[2])
println(" mean K/M: ",meanparamvals[1]/meanparamvals[2])
println("Mean K/M should be close to $(K/M)")

println("================")
println("Plotting results")
println("================")

###Plot the loglikelihood values across samples
###After an initial few low values, this should remain relatively high
subplot(222)
plot(1:numsamples(c),logposterior(c,m))
title("Log-Posterior values across samples")
xlabel("Samples")
ylabel("Log-Posterior")

###Plot a scatter plot of K vs M values
###These should be spread around the K/M == 10.0 line (the diagonal in the figure)
ax3 = subplot(223)
ax3[:set_xlim]([lows[1],highs[1]])
ax3[:set_ylim]([lows[2],highs[2]])
scatter(vec(getindex(samples(c),1,:)),vec(getindex(samples(c),2,:)),marker=".",color="blue")
ax3[:set_aspect](abs(highs[1]-lows[1])/abs(highs[2]-lows[2]))
title("MCMC samples of Spring-Mass ODE parameters")
xlabel("Stiffness K (N/m)")
ylabel("Mass M (kg)")
grid("on")

###Plot a histogram of K/M values, which should peak around the true ratio of K/M
ax4 = subplot(224)
kml,kmu = K/M-K/M/20.0,K/M+K/M/20.0
ax4[:set_xlim]([kml,kmu])
nbins = linspace(kml,kmu,100)
h = PyPlot.plt[:hist]((vec(getindex(samples(c),1,:)./getindex(samples(c),2,:))),nbins)
grid("on")
xlabel("K/M")
ylabel("Number of Samples")
title("Histogram of K/M values")

###Finally, plot the average model results in the data window
subplot(221)
modeldata = evaluate!(m,vec(meanparamvals))
plot(dataindex(m),modeldata[:,1];label="model location")
plot(dataindex(m),modeldata[:,2];label="model velocity")


text:
export domains
export PiecewisePoly
export update!
export lengthpoints

mutable struct PiecewisePoly{T<:Number} <: AbstractPolynomial{T}
    polys::Vector{AbstractPolynomial{T}}#array of polynomials
    domains::Vector{Interval{T}}
    var::Symbol
    function PiecewisePoly{T}(subP::AbstractPolynomial{T}, pieDomain::Interval{T}, var::Symbol) where {T<:Number}
        checkends(pieDomain,subP,true)
        checkinfs(pieDomain, domain(subP))
        return new{T}([subP], [pieDomain], var) # new initializing globalPoly must be created with an initial poly to avoid checks further down the line
    end
end

function PiecewisePoly{T}(subP::N, lower, upper; var::SymbolLike=:x,  containlower::Bool=containslower(domain(subP)), containupper::Bool=containsupper(domain(subP))) where {T,N<:AbstractPolynomial}
    lower = convert(T, lower)
    upper = convert(T, upper)
    subP = convert(Base.typename(N).wrapper{T},subP)
    pieDomain = Interval(lower, upper, containlower, containupper)
    return PiecewisePoly{T}(subP, pieDomain, var)
end
function PiecewisePoly{T}(subP::AbstractPolynomial; kwargs...) where {T}
    lower=first(domain(subP))
    upper=last(domain(subP))
    if T<:Int && (lower==Inf || upper==Inf)
        error("PiecewisePoly of type $T does not support infinite endpoints, consider setting upper and lower bounds or converting to type <:AbstractFloat")
    end
    PiecewisePoly{T}(subP, lower, upper; kwargs...)
end
function PiecewisePoly(subP::AbstractPolynomial{T}, lower, upper;kwargs...) where {T}
    return PiecewisePoly{T}(subP, lower, upper;kwargs...)
end
function PiecewisePoly(subP::AbstractPolynomial{T};kwargs...) where {T}
    return PiecewisePoly{T}(subP;kwargs...)
end

function showterm(io::IO, ::Type{<:PiecewisePoly{T}}, pj::Any, var, j, first::Bool, mimetype) where {T}
    !first &&  print(io, " ")
    print(io, "$(pj)")
    return true
end

length(pieP::PiecewisePoly) = length(pieP.polys)

function lengthpoints(pieP::PiecewisePoly)
    lengthp = 0
    pieDomains = domains(pieP)
    lengthp += length(pieP[1])
    for i in 2:length(pieP)
        if iscontinuous(pieP) && last(pieDomains[i-1]) == first(pieDomains[i])
            lengthp -= 1
        end
        @inbounds lengthp += length(pieP[i])
    end
    return lengthp
end

Base.eachindex(pieP::PiecewisePoly) = 1:length(pieP)
# indexing global polys
Base.getindex(pieP::PiecewisePoly, i::Int) = pieP.polys[i]
#Base.setindex!(pieP::PiecewisePoly, p::AbstractPolynomial, i) = #todo
Base.firstindex(pieP::PiecewisePoly) = 1
Base.lastindex(pieP::PiecewisePoly) = length(pieP.polys)

domains(pieP::PiecewisePoly) = pieP.domains

iscontinuous(pieP::PiecewisePoly) = true # update this in future

function push!(pieP::PiecewisePoly, subP::AbstractPolynomial, pieDomain::Interval)
    pieDomains = domains(pieP)
    checkends(pieDomain,subP,iscontinuous(pieP))
    subDomain = checkinfs(pieDomain, domain(subP)) # check this mapping for infinite local domains mapping to non infinite global domians, see internal comment
    if last(pieDomains[end]) <= first(pieDomain)
        if last(pieDomains[end]) == first(pieDomain) && iscontinuous(pieP) && pieP(last(pieDomains[end])) != subP(first(subDomain))
            error("new poly \"N\" lower endpoint y value is not coincidence with PiecewisePoly \"G\" (G=$(pieP(last(pieDomains[end]))) vs N=$(subP(first(domain(subP))))). Use a non continuous PiecewisePoly if this is what you want.")
        end
        push!(pieP.polys, subP)
        push!(pieP.domains, pieDomain)
        return pieP
    end
    if first(pieDomains[1]) >= last(pieDomain)
        if first(pieDomains[1]) == last(pieDomain) && iscontinuous(pieP) && pieP(first(pieDomains[1])) != subP(last(subDomain))
            error("new poly \"N\" upper endpoint y value is not coincidence with PiecewisePoly \"G\" (G=$(pieP(first(pieDomains[1]))) vs N=$(subP(last(domain(subP))))). Use a non continuous PiecewisePoly if this is what you want.")
        end
        pushfirst!(pieP.polys, subP)
        return pushfirst!(pieP.domains, pieDomain)
    end
    for i in 1:length(pieP)-1
        if last(pieDomains[i+1]) > first(pieDomain)
            if  last(pieDomains[i]) > first(pieDomain)
                conflict = intersect(pieDomain, pieDomains[i])
                error("the domain intersects with an existing polynomial at location [$(i)]. Intersection is $conflict .")
            end
            if first(pieDomains[i+1]) < last(pieDomain)
                conflict = intersect(pieDomain, pieDomains[i+1])
                error("the domain intersects with an existing polynomial at location [$(i+1)]. Intersection is $conflict .")
            end
            if last(pieDomains[i]) == first(pieDomain) && iscontinuous(pieP) && pieP(last(pieDomains[i])) != subP(first(subDomain))
                error("new poly \"N\" lower endpoint y value is not coincidence with PiecewisePoly \"G\" (G=$(pieP(last(pieDomains[i]))) vs N=$(subP(first(domain(subP))))). Use a non continuous PiecewisePoly if this is what you want.")
            end
            if first(pieDomains[i+1]) == last(pieDomain) && iscontinuous(pieP) && pieP(first(pieDomains[i+1])) != subP(last(subDomain))
                error("new poly \"N\" upper endpoint y value is not coincidence with PiecewisePoly \"G\" (G=$(pieP(first(pieDomains[i+1]))) vs N=$(subP(last(domain(subP))))). Use a non continuous PiecewisePoly if this is what you want.")
            end
            insert!(pieP.polys, i+1, subP)
            insert!(pieP.domains, i+1, pieDomain)
            return pieP
        end
    end

    conflict = intersect(pieDomain, pieDomains[end])
    error("new domain endpoints intersect existing polynomial. Intersection is $conflict")
end
function push!(pieP::PiecewisePoly{T}, subP::AbstractPolynomial, gloLower, gloUpper; containlower=containslower(domain(subP)), containupper=containsupper(domain(subP))) where {T}
    if T<:Int && (gloLower==Inf || gloUupper==Inf)
        error("PiecewisePoly of type $T does not support infinite endpoints, consider setting upper and lower bounds or converting to type <:AbstractFloat")
    end
    pieDomain = Interval(gloLower, gloUpper, containlower, containupper)
    return push!(pieP, subP, pieDomain)
end
function push!(pieP::PiecewisePoly, subP::AbstractPolynomial; kwargs...)
    lower=first(domain(subP))
    upper=last(domain(subP))
    return push!(pieP, subP, lower, upper; kwargs...)
end

function domain(pieP::PiecewisePoly)
    low,up = domains(pieP)[1], domains(pieP)[end]
    return Interval(first(low), last(up), containslower(low), containsupper(up))
end

function (pieP::PiecewisePoly{T})(x::S) where {T,S<:Number}
    pieDomain = domains(pieP)
    for i in 1:length(pieP)
        if x in pieDomain[i]
            return pieP[i](global2local(domain(pieP[i]), pieDomain[i], x))
        end
    end
    if x ∉ domain(pieP)
        pDomain = domain(pieP)
        throw(DomainError(x,"outside of domain $pDomain"))
    end
    closestMatch = argmin([min(abs.(x .- [first(domains(pieP)[i]), last(domains(pieP)[i])])...) for i in 1:length(pieP)])
    throw(DomainError(x,"not found in PiecewisePoly domains. Closest polynomial has domain $(domains(pieP)[closestMatch]) at index [$closestMatch]."))
end

function update!(pieP::PiecewisePoly{T}, y::AbstractVector{T}) where {T}
    pieDomains = domains(pieP)
    continuouspLengths = [ iscontinuous(pieP) && last(pieDomains[i-1]) == first(pieDomains[i]) ? length(pieP[i]) - 1 : length(pieP[i]) for i in 2:length(pieP)]
    pushfirst!(continuouspLengths,length(pieP[1]))
    sum(continuouspLengths) !=length(y) && throw(DimensionMismatch("tried to assign $(length(y)) elements to $(sum(continuouspLengths)) destinations"))
    Threads.@threads for i in 1:length(pieP)
        if continuouspLengths[i] == length(pieP[i]) # case where endpoints are not coincident
            pieP[i][:] = y[sum(continuouspLengths[1:i-1])+1:sum(continuouspLengths[1:i])]
        else
            pieP[i][:] = y[sum(continuouspLengths[1:i-1]):sum(continuouspLengths[1:i])]
        end
    end
    return pieP
end

function checkinfs(pieDomain::Interval, subDomain::Interval)
    subLower = first(subDomain)
    subUpper = last(subDomain)
    pieLower = first(pieDomain)
    pieUpper = last(pieDomain)
    if isinf(pieLower) || isinf(pieUpper) || isinf(subLower) || isinf(subUpper)
        if isinf(pieLower) && isinf(pieUpper)
            if xor(isinf(subLower),isinf(subUpper))
                throw(DomainError(pieDomain,"this piecewise domain is not compatible with sub domain $subDomain"))
            else
                return subDomain
            end
        elseif isinf(subLower)&&isinf(subUpper)
            return pieDomain # return global domain so domain mapping doesn't get confused, check if this needs to be done elsewhere too
        elseif !issubset(pieDomain,subDomain)
            if !isinf(pieLower)&&!isinf(pieUpper)
                throw(DomainError(pieDomain,"this piecewise domain is not a subset of sub domain $subDomain"))
            end
            if issubset(Interval(-Inf,Inf, false, false),Interval(min(subLower,pieLower), max(subUpper,pieUpper)))
                throw(DomainError(pieDomain,"this piecewise domain is not compatible with sub domain $subDomain"))
            end
        end
    end
    return subDomain
end

function checkends(pieDomain::Interval,subP::AbstractPolynomial, iscontinuous::Bool)
    subDomain = domain(subP)
    if iscontinuous
        ( !containslower(subDomain) || !containsupper(subDomain) ) && throw(DomainError(subP),"continuous PiecewisePoly can only contain sub polynomial that have closed endpoints")
    else
        !containslower(subDomain) && containslower(pieDomain) && throw(DomainError(pieDomain,"piecewise lower endpoints incompatible with sub polynomial domain $(subDomain). Consider either opening the piecewise endpoint or closing the sub endpoint."))#check first that local and global bound endpoints are compatible
        !containsupper(subDomain) && containsupper(pieDomain) && throw(DomainError(pieDomain,"piecewise upper endpoints incompatible with sub polynomial domain $(subDomain). Consider either opening the piecewise endpoint or closing the sub endpoint."))
    end
end

# following Julia convert() convention where thing being converted to comes first
mapendpoints(from::Interval, to::Interval, x::Number) = first(to) + (last(to) - first(to)) * (x - first(from)) / ( last(from) - first(from) )
function global2local(subDomain::Interval, pieDomain::Interval, x::Number)
    subLower = first(subDomain)
    subUpper = last(subDomain)
    pieLower = first(pieDomain)
    pieUpper = last(pieDomain)
    if isinf(pieLower) || isinf(pieUpper) || isinf(subLower) || isinf(subUpper)
        if (!isinf(pieLower)&&!isinf(pieUpper)) || (isinf(subLower) && isinf(subUpper))
            return x
        elseif (!isinf(subLower)&&!isinf(subUpper))
            if (isinf(pieLower) && isinf(pieUpper))
                return mapendpoints(Interval(-1,1),subDomain,2*atan(x) / pi)
            elseif isinf(pieUpper)
                return mapendpoints(Interval(-1,1),subDomain,(x-1)/(1+x)) # allow user to specify custom function here?
            else
                return mapendpoints(Interval(-1,1),subDomain,(-x-1)/(x-1))
            end
        elseif isinf(pieUpper)
            return x - pieLower + subLower
        else
            return x - pieUpper + subUpper
        end

    end
    return subLower + (subUpper - subLower) * (x - pieLower) / ( pieUpper - pieLower )
end


text:
using Revise
using gtpmip

graph = optiplan_solve("test/pddl/domain.pddl", "test/pddl/problem.pddl")
# solve_optiplan!(model)
# 1


text:
##

using ACE, StaticArrays
using Random, Printf, Test, LinearAlgebra, ACE.Testing
using ACE: evaluate, evaluate_d, SymmetricBasis, PIBasis, rand_vec3
using ACE.Random: rand_rot, rand_refl
using ACEbase.Testing: fdtest

##

# construct the 1p-basis
maxdeg = 6
ord = 3
Bsel = SimpleSparseBasis(ord, maxdeg)

B1p = ACE.Utils.RnYlm_1pbasis(; maxdeg=maxdeg)

# generate a configuration
nX = 10
_randX() = State(rr = rand_vec3(B1p["Rn"]))
Xs = [ _randX() for _=1:nX ]
cfg = ACEConfig(Xs)

##

@info("SymmetricBasis construction and evaluation: EuclideanVector")

φ = ACE.EuclideanVector(Float64)
pibasis = PIBasis(B1p, Bsel; property = φ)
basis = SymmetricBasis(φ, pibasis)
@time SymmetricBasis(φ, pibasis)

BB = evaluate(basis, cfg)

Iz = findall(iszero, sum(norm, basis.A2Bmap, dims=1)[:])
if !isempty(Iz)
   @warn("The A2B map for EuclideanVector has $(length(Iz))/$(length(basis.pibasis)) zero-columns!!!!")
end

##

@info("Test FIO")
using ACEbase.Testing: test_fio
@warn("turned off failing FIO test")
# println_slim(@test(all(test_fio(basis; warntype = false))))

##

@info("Test equivariance properties for real version")

tol = 1e-10

@info("check for rotation, permutation and inversion equivariance")
for ntest = 1:30
   local Xs, BB
   Xs = [ _randX() for _=1:nX ]
   BB = evaluate(basis, ACEConfig(Xs))
   Q = rand([-1,1]) * ACE.Random.rand_rot()
   Xs_rot = Ref(Q) .* shuffle(Xs)
   BB_rot = evaluate(basis, ACEConfig(Xs_rot))
   print_tf(@test all([ norm(Q' * b1 - b2) < tol
                        for (b1, b2) in zip(BB_rot, BB)  ]))
end
println()

# @info("Test equivariance properties for complex version")
#
# basis = SymmetricBasis(φ, pibasis; isreal=false)
# # a stupid but necessary test
# BB = evaluate(basis, cfg)
# BB1 = basis.A2Bmap * evaluate(basis.pibasis, cfg)
# println(@test isapprox(BB, BB1, rtol=1e-10)) # MS: This test will fail for isreal=true
#
#
# @info("check for rotation, permutation and inversion equivariance")
# for ntest = 1:30
#    local Xs, BB
#    Xs = rand(PositionState{Float64}, B1p["Rn"].basis, nX)
#    BB = evaluate(basis, ACEConfig(Xs))
#    Q = rand([-1,1]) * ACE.Random.rand_rot()
#    Xs_rot = Ref(Q) .* shuffle(Xs)
#    BB_rot = evaluate(basis, ACEConfig(Xs_rot))
#    print_tf(@test all([ norm(Q' * b1 - b2) < tol
#                         for (b1, b2) in zip(BB_rot, BB)  ]))
# end
# println()

# ## keep for further profiling
#
# φ = ACE.EuclideanVector(Complex{Float64})
# pibasis = PIBasis(B1p, ord, maxdeg; property = φ, isreal = false)
# basis = SymmetricBasis(pibasis, φ)
# @time SymmetricBasis(pibasis, φ);
#
# Profile.clear(); # Profile.init(; delay = 0.0001)
# @profile SymmetricBasis(pibasis, φ);
# ProfileView.view()

##

@info(" ... derivatives")
_rrval(x::ACE.XState) = x.rr
for ntest = 1:30
   Us = randn(SVector{3, Float64}, length(Xs))
   C = randn(typeof(φ.val), length(basis))
   F = t -> sum( sum(c .* b.val)
                 for (c, b) in zip(C, ACE.evaluate(basis, ACEConfig(Xs + t[1] * Us))) )
   dF = t -> [ sum( sum(c .* db)
                    for (c, db) in zip(C, _rrval.(ACE.evaluate_d(basis, ACEConfig(Xs + t[1] * Us))) * Us) ) ]
   print_tf(@test fdtest(F, dF, [0.0], verbose=false))
end
println()

##


text:
"""
  AbstractFTPlan{T,D,R}

Abstract type for any NFFT-like plan (NFFT, NFFT, NFCT, NFST).
* T is the element type (Float32/Float64)
* D is the number of dimensions of the input array.
* R is the number of dimensions of the output array.
"""
abstract type AbstractFTPlan{T,D,R} end

"""
  AbstractRealFTPlan{T,D,R}

Abstract type for either an NFCT plan or an NFST plan.
* T is the element type (Float32/Float64)
* D is the number of dimensions of the input array.
* R is the number of dimensions of the output array.
"""
abstract type AbstractRealFTPlan{T,D,R} <: AbstractFTPlan{T,D,R} end

"""
  AbstractComplexFTPlan{T,D,R}

Abstract type for either an NFFT plan or an NNFFT plan.
* T is the element type (Float32/Float64)
* D is the number of dimensions of the input array.
* R is the number of dimensions of the output array.
"""
abstract type AbstractComplexFTPlan{T,D,R} <: AbstractFTPlan{T,D,R} end

"""
  AbstractNFFTPlan{T,D,R}

Abstract type for an NFFT plan.
* T is the element type (Float32/Float64)
* D is the number of dimensions of the input array.
* R is the number of dimensions of the output array.
"""
abstract type AbstractNFFTPlan{T,D,R} <: AbstractComplexFTPlan{T,D,R} end

"""
  AbstractNFCTPlan{T,D,R}

Abstract type for an NFCT plan.
* T is the element type (Float32/Float64)
* D is the number of dimensions of the input array.
* R is the number of dimensions of the output array.
"""
abstract type AbstractNFCTPlan{T,D,R} <: AbstractRealFTPlan{T,D,R} end

"""
  AbstractNFSTPlan{T,D,R}

Abstract type for an NFST plan.
* T is the element type (Float32/Float64)
* D is the number of dimensions of the input array.
* R is the number of dimensions of the output array.
"""
abstract type AbstractNFSTPlan{T,D,R} <: AbstractRealFTPlan{T,D,R} end

"""
  AbstractNNFFTPlan{T,D,R}

Abstract type for an NNFFT plan.
* T is the element type (Float32/Float64)
* D is the number of dimensions of the input array.
* R is the number of dimensions of the output array.
"""
abstract type AbstractNNFFTPlan{T,D,R} <: AbstractComplexFTPlan{T,D,R} end

#####################
# Function needed to make AbstractFTPlan an operator
#####################

Base.eltype(p::AbstractComplexFTPlan{T}) where T = Complex{T}
Base.eltype(p::AbstractRealFTPlan{T}) where T = T

Base.size(p::AbstractFTPlan) = (prod(size_out(p)), prod(size_in(p)))
Base.length(p::AbstractFTPlan) = prod(size(p))

Base.size(p::Adjoint{Complex{T}, U}) where {T, U<:AbstractComplexFTPlan{T}} =
  (prod(size_in(p.parent)), prod(size_out(p.parent)))

Base.size(p::Transpose{T, U}) where {T, U<:AbstractRealFTPlan{T}} =
  (prod(size_in(p.parent)), prod(size_out(p.parent)))

LinearAlgebra.adjoint(p::AbstractComplexFTPlan{T,D,R}) where {T,D,R} =
  Adjoint{Complex{T}, typeof(p)}(p)

LinearAlgebra.transpose(p::AbstractRealFTPlan{T,D,R}) where {T,D,R} =
  Transpose{T, typeof(p)}(p)

size_in(p::Adjoint{Complex{T},<:AbstractComplexFTPlan{T,D,R}}) where {T,D,R} = size_out(p.parent)
size_out(p::Adjoint{Complex{T},<:AbstractComplexFTPlan{T,D,R}}) where {T,D,R} = size_in(p.parent)

size_in(p::Transpose{T,<:AbstractRealFTPlan{T,D,R}}) where {T,D,R} = size_out(p.parent)
size_out(p::Transpose{T,<:AbstractRealFTPlan{T,D,R}}) where {T,D,R} = size_in(p.parent)

# The following are required both, because there seem to be special array handling in base
function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, p::Transpose{T,<:AbstractRealFTPlan{T,D,R}}) where {T,D,R}
  print(io, "Transposed ")
  print(io, p.parent)
end

function Base.show(io::IO, p::Transpose{T,<:AbstractRealFTPlan{T,D,R}}) where {T,D,R}
  print(io, "Transposed ")
  print(io, p.parent)
end

function Base.show(io::IO, p::Adjoint{Complex{T},<:AbstractComplexFTPlan{T,D,R}}) where {T,D,R}
  print(io, "Adjoint ")
  print(io, p.parent)
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, p::Adjoint{Complex{T},<:AbstractComplexFTPlan{T,D,R}}) where {T,D,R}
  print(io, "Adjoint ")
  print(io, p.parent)
end

function Base.copy(p::Transpose{T,<:AbstractRealFTPlan{T,D,R}}) where {T,D,R}
  transpose(copy(p.parent))
end

function Base.copy(p::Adjoint{Complex{T},<:AbstractComplexFTPlan{T,D,R}}) where {T,D,R}
  adjoint(copy(p.parent))
end

#####################
# define interface
#####################

"""
    mul!(fHat, p, f) -> fHat

Inplace NFFT/NFCT/NFST/NNFFT transforming the `D` dimensional array `f` to the `R` dimensional array `fHat`.
The transformation is applied along `D-R+1` dimensions specified in the plan `p`.
"""
@mustimplement LinearAlgebra.mul!( fHat::AbstractArray, p::AbstractFTPlan{T}, f::AbstractArray) where {T}

"""
    mul!(f, p, fHat) -> f

Inplace adjoint NFFT/NNFFT transforming the `R` dimensional array `fHat` to the `D` dimensional array `f`.
The transformation is applied along `D-R+1` dimensions specified in the plan `p`.
"""
@mustimplement LinearAlgebra.mul!(f::AbstractArray, p::Adjoint{Complex{T},<:AbstractComplexFTPlan{T}}, fHat::AbstractArray) where {T}

"""
    mul!(f, p, fHat) -> f

Inplace transposed NFCT/NFST transforming the `R` dimensional array `fHat` to the `D` dimensional array `f`.
The transformation is applied along `D-R+1` dimensions specified in the plan `p`.
"""
@mustimplement LinearAlgebra.mul!(f::AbstractArray, p::Transpose{T,<:AbstractRealFTPlan{T}}, fHat::AbstractArray) where {T}

"""
    size_in(p)

Size of the input array for the plan p (NFFT/NFCT/NFST/NNFFT).
The returned tuple has `R` entries.
Note that this will be the output size for the transposed / adjoint operator.
"""
@mustimplement size_in(p::AbstractFTPlan{T,D,R}) where {T,D,R}

"""
    size_out(p)

Size of the output array for the plan p (NFFT/NFCT/NFST/NNFFT).
The returned tuple has `R` entries.
Note that this will be the input size for the transposed / adjoint operator.
"""
@mustimplement size_out(p::AbstractFTPlan{T,D,R}) where {T,D,R}

"""
    nodes!(p, x) -> p

Change nodes `x` in the plan `p` operation and return the plan.
"""
@mustimplement nodes!(p::AbstractFTPlan{T}, x::Matrix{T}) where {T}

## Optional Interface ##
# The following methods can but don't need to be implemented

@mustimplement deconvolve!(p::AbstractNFFTPlan, f::AbstractArray, g::AbstractArray)
@mustimplement deconvolve_transpose!(p::AbstractNFFTPlan, g::AbstractArray, f::AbstractArray)
@mustimplement convolve!(p::AbstractNFFTPlan, g::AbstractArray, fHat::AbstractArray)
@mustimplement convolve_transpose!(p::AbstractNFFTPlan, fHat::AbstractArray, g::AbstractArray)


text:
using LibPQ, Tables, Dates
using SHA
using JSON

function creatConnection()
    return LibPQ.Connection("dbname=klingelbeutel")
end

function createuser(_login::String, _vorname::String, _nachname::String, _email::String, password::String )
    conn = creatConnection()
    execute(conn, "BEGIN;")
    p = bytes2hex(sha256(password))
    sql = "INSERT INTO app_user(login, vorname, nachname, email, password)values(\$1, \$2, \$3, \$4, \$5)"
    LibPQ.load!((login=[_login], vorname=[_vorname], nachname=[_nachname], email=[_email], password=[p]),conn, sql)
    execute(conn, "COMMIT;")
end

function getNews()
    conn = creatConnection()
    sql = "select url, news_short, news_text, datum from news order by datum desc limit 10;"
    result = execute(conn, sql)
    data = columntable(result)
    ergo_array = []
    for i in 1:length(data[2])
        row = Dict("news_short" => data[2][i],"news_text" => data[3][i], "datum" => data[4][i] )
        push!(ergo_array,row)

    end
    println(JSON.json(ergo_array))

    println(data[2][1])
    println(length(data[2]))
end

getNews()
#createuser("tom", "Thomas", "Schmidt", "t_schmidt@web.de", "Geheim#1")

max_stars_repo_name: Kompendix/Klingelbeutel

max_issues_repo_name: Kompendix/Klingelbeutel

max_forks_repo_name: Kompendix/Klingelbeutel


text:
module StopWatchesTests

import StopWatches
using Test

if !@isdefined(isnothing)
    using Compat: isnothing
end

include("utils.jl")

function include_tests(m = @__MODULE__, dir = @__DIR__)
    for file in readdir(dir)
        if match(r"^test_.*\.jl$", file) !== nothing
            Base.include(m, joinpath(dir, file))
        end
    end
end

include_tests()

function collect_modules(root::Module)
    modules = Module[]
    for n in names(root, all = true)
        m = getproperty(root, n)
        m isa Module || continue
        m === root && continue
        startswith(string(nameof(m)), "Test") || continue
        push!(modules, m)
    end
    # Hack: move TestAnalysis at the end
    let i = findfirst(==(TestAnalysis), modules)
        if i !== nothing
            deleteat!(modules, i)
            push!(modules, TestAnalysis)
        end
    end
    return modules
end

collect_modules() = collect_modules(@__MODULE__)

this_project() = joinpath(dirname(@__DIR__), "Project.toml")

function is_in_path()
    project = this_project()
    paths = Base.load_path()
    project in paths && return true
    realproject = realpath(project)
    realproject in paths && return true
    matches(path) = path == project || path == realproject
    return any(paths) do path
        matches(path) || matches(realpath(path))
    end
end

function with_project(f)
    is_in_path() && return f()
    load_path = copy(LOAD_PATH)
    push!(LOAD_PATH, this_project())
    try
        f()
    finally
        append!(empty!(LOAD_PATH), load_path)
    end
end

function runtests_impl(modules = collect_modules())
    @testset "$(nameof(m))" for m in modules
        tests = map(names(m, all = true)) do n
            n == :test || startswith(string(n), "test_") || return nothing
            f = getproperty(m, n)
            f !== m || return nothing
            parentmodule(f) === m || return nothing
            applicable(f) || return nothing  # removed by Revise?
            return f
        end
        filter!(!isnothing, tests)
        @testset "$f" for f in tests
            @debug "Testing $m.$f"
            f()
        end
    end
end

runtests() = with_project(runtests_impl)

end # module


text:
using DSGE, ModelConstructors, Dates, Plots
fn = dirname(@__FILE__)

# This script shows how to calculate state space systems
# with time-varying information sets and shows how they
# affect forecasts. See the documentation
# for mathematical and theoretical clarifications on
# what we mean precisely by time-varying information sets.

# What do you want to do?
calculate_system = true
modal_forecast   = true
create_plots     = true

## Set up

# Initialize model w/desired dates
m = Model1002("ss10")
usual_model_settings!(m, "200930"; cdvt = "200930", fcast_date = quartertodate("2020-Q3"))
m <= Setting(:use_population_forecast, false)
m <= Setting(:saveroot, joinpath("$(fn)", "..", "save"))
m <= Setting(:dataroot, joinpath("$(fn)", "..", "save", "input_data"))
m <= Setting(:date_forecast_end, quartertodate("2029-Q4"))

if calculate_system
    # Set up two different sets of equilibrium conditions,
    # one of which corresponds to an expected temporary ZLB
    # from 2020:Q3 to 2023:Q4 and another which
    # corresponds to an expected temporary ZLB from
    # 2020:Q3 to 2024:Q4 (one additional year).
    m <= Setting(:replace_eqcond, true) # turn on replace_eqcond and gensys2
    m <= Setting(:gensys2, true)        # to allow temporary alternative policies (namely the ZLB)
    regime_eqcond_info_shortzlb = Dict{Int, EqcondEntry}()
    replace_eqcond_func_dict_shortzlb = Dict()
    for i in 2:(length(DSGE.quarter_range(Date(2020, 9, 30), Date(2023, 12, 31))) + 1)
        # replace_eqcond_func_dict_shortzlb[i] = DSGE.zero_rate_replace_eq_entries
        regime_eqcond_info_shortzlb[i] = EqcondEntry(zero_rate())
    end
    regime_eqcond_info_longzlb = Dict{Int, EqcondEntry}()
    replace_eqcond_func_dict_longzlb = Dict()
    for i in 2:(length(DSGE.quarter_range(Date(2020, 9, 30), Date(2024, 12, 31))) + 1)
        # replace_eqcond_func_dict_longzlb[i] = DSGE.zero_rate_replace_eq_entries
        regime_eqcond_info_longzlb[i] = EqcondEntry(zero_rate())
    end
    m <= Setting(:zero_rate_zlb_value, 0.)
    m <= Setting(:temporary_altpolicy_names, [:zero_rate])

    # Set up regime-switching dates and indices
    regime_dates = Dict{Int, Date}(1 => date_presample_start(m))
    for (i, d) in enumerate(DSGE.quarter_range(Date(2020, 9, 30), Date(2025, 3, 31)))
        # Add regimes for each quarter up to 2025:Q1 b/c the same :regime_dates dictionary
        # will be used when calculating the state-space system for the shorter ZLB policy
        # and the longer ZLB policy. If we added regimes only up to 2024:Q1, then
        # the longer ZLB policy will run into an error.
        regime_dates[i + 1] = d
    end
    m <= Setting(:regime_dates, regime_dates)
    m <= Setting(:regime_switching, true)
    setup_regime_switching_inds!(m)

    ## Calculate state-space system with time-varying information sets

    # First, let's calculate an example where people learn about the short ZLB
    # starting in 2020:Q3. In 2020:Q2 (and before), people do NOT expect a ZLB.
    reg_2024Q1 = get_setting(m, :n_regimes) - 4 # This is the regime corresponding to 2024:Q1
    @assert get_setting(m, :regime_dates)[reg_2024Q1] == Date(2024, 3, 31)

    # How the time-varying information set is constructed:
    # In regime 1, agents only need to use the regime 1 transition equation
    # From 2020:Q3 to 2024:Q1, people need to use the transition equations for 2020:Q3 to 2024:Q1
    # to construct the measurement equation. In 2024:Q2 onward, people again only need to use
    # the current regime's transition equation. Alternatively, since the transition equation
    # in 2024:Q1 is the same as the one in 2024:Q2 - 2025:Q1, we could have also assumed that,
    # during 2020:Q3 to 2024:Q1, people use the transition equation from 2020:Q3 to 2025:Q1.
    # This alternative approach is mathematically equivalent.
    m <= Setting(:tvis_information_set,
                 vcat([1:1],
                      [i:reg_2024Q1 for i in 2:reg_2024Q1],
                      [[i:i] for i in (reg_2024Q1 + 1):get_setting(m, :n_regimes)]...))
    # m <= Setting(:replace_eqcond_func_dict, replace_eqcond_func_dict_shortzlb)
    m <= Setting(:regime_eqcond_info, regime_eqcond_info_shortzlb)
    system_shortzlb = compute_system(m; tvis = true)

    # We could use also these settings to generate the same results.
    # m <= Setting(:tvis_replace_eqcond_func_dict, [replace_eqcond_func_dict_shortzlb])
    # m <= Setting(:tvis_select_system, ones(Int, 7))

    # Note that if :tvis_replace_eqcond_func_dict is already a setting in m,
    # then we use this setting to determine the dictionary of functions for
    # :replace_eqcond_func_dict rather than what is already stored
    # in the setting :replace_eqcond_func_dict.

    # Now, let's calculate an example where people learn about the short ZLB
    # starting in 2020:Q3, but in 2021:Q2, people learn that the long ZLB
    # will go into action.

    # First, we add info about the equilibrium conditions
    reg_2021Q1 = DSGE.subtract_quarters(Date(2021, 3, 31), Date(2020, 9, 30)) + 2
    # m <= Setting(:tvis_replace_eqcond_func_dict, [replace_eqcond_func_dict_shortzlb, replace_eqcond_func_dict_longzlb])
    m <= Setting(:tvis_regime_eqcond_info, [regime_eqcond_info_shortzlb, regime_eqcond_info_longzlb])
    m <= Setting(:tvis_select_system, vcat(ones(Int, reg_2021Q1), # Use set 1 of conditions for regimes up to (and incl.) 2021:Q1
                                           fill(2, get_setting(m, :n_regimes) - reg_2021Q1))) # Use set 2 for remaining regimes

    # Second, we add info about the information sets
    m <= Setting(:tvis_information_set,
                 vcat([1:1], # Regime 1 doesn't change
                      [i:reg_2024Q1 for i in 2:reg_2021Q1], # This is the same as before, but the info set changes after 2021:Q1
                      [i:get_setting(m, :n_regimes) for i in (reg_2021Q1 + 1):get_setting(m, :n_regimes)])) # Now we use regimes to 2025:Q1 for the info set
    system_shortlong_zlb = compute_system(m; tvis = true)

    # Some tests that should pass
    for i in 1:reg_2021Q1
        @assert system_shortzlb[i, :ZZ] == system_shortlong_zlb[i, :ZZ]
        @assert system_shortzlb[i, :TTT] == system_shortlong_zlb[i, :TTT]
        @assert system_shortzlb[i, :ZZ] != system_shortlong_zlb[(reg_2021Q1 + 1), :ZZ]
        @assert system_shortzlb[i, :TTT] != system_shortlong_zlb[(reg_2021Q1 + 1), :TTT]
    end
    for i in (reg_2021Q1 + 1):(get_setting(m, :n_regimes) - 2)
        @assert system_shortzlb[i, :ZZ] != system_shortlong_zlb[i, :ZZ]
        @assert system_shortzlb[i, :TTT] != system_shortlong_zlb[i, :TTT]
    end
    @assert system_shortzlb[get_setting(m, :n_regimes) - 1, :ZZ] ≈ # b/c transition equation matters only for
        system_shortlong_zlb[get_setting(m, :n_regimes) - 1, :ZZ]  # forward expectations in measurement equation
    @assert system_shortzlb[get_setting(m, :n_regimes), :ZZ] ==
        system_shortlong_zlb[get_setting(m, :n_regimes), :ZZ]
    @assert system_shortzlb[get_setting(m, :n_regimes), :TTT] ==
        system_shortlong_zlb[get_setting(m, :n_regimes), :TTT]
end

if modal_forecast
    # Set up in case calculate_system is not true
    if !calculate_system
        m <= Setting(:replace_eqcond, true)
        m <= Setting(:gensys2, true)
        regime_eqcond_info_shortzlb = Dict{Int, EqcondEntry}()
        replace_eqcond_func_dict_shortzlb = Dict()
        for i in 2:(length(DSGE.quarter_range(Date(2020, 9, 30), Date(2023, 12, 31))) + 1)
            # replace_eqcond_func_dict_shortzlb[i] = DSGE.zero_rate_replace_eq_entries
            regime_eqcond_info_shortzlb[i] = EqcondEntry(zero_rate())
        end
        regime_eqcond_info_longzlb = Dict{Int, EqcondEntry}()
        replace_eqcond_func_dict_longzlb = Dict()
        for i in 2:(length(DSGE.quarter_range(Date(2020, 9, 30), Date(2024, 12, 31))) + 1)
            # replace_eqcond_func_dict_longzlb[i] = DSGE.zero_rate_replace_eq_entries
            regime_eqcond_info_longzlb[i] = EqcondEntry(zero_rate())
        end

        regime_dates = Dict{Int, Date}(1 => date_presample_start(m))
        for (i, d) in enumerate(DSGE.quarter_range(Date(2020, 9, 30), Date(2025, 3, 31)))
            regime_dates[i + 1] = d
        end
        m <= Setting(:zero_rate_zlb_value, 0.)
        m <= Setting(:temporary_altpolicy_names, [:zero_rate])

        m <= Setting(:regime_dates, regime_dates)
        m <= Setting(:regime_switching, true)
        setup_regime_switching_inds!(m)
    end

    # Set up two modal forecasts, to show the effect of extending the ZLB
    reg_2024Q1 = get_setting(m, :n_regimes) - 4
    m <= Setting(:tvis_information_set, vcat([1:1], [i:reg_2024Q1 for i in 2:reg_2024Q1],
                                             [[i:i] for i in (reg_2024Q1 + 1):get_setting(m, :n_regimes)]...))
    # m <= Setting(:tvis_replace_eqcond_func_dict, [replace_eqcond_func_dict_shortzlb])
    m <= Setting(:tvis_regime_eqcond_info, [regime_eqcond_info_shortzlb])
    # m <= Setting(:replace_eqcond_func_dict, replace_eqcond_func_dict_shortzlb)
    m <= Setting(:regime_eqcond_info, regime_eqcond_info_shortzlb)

    usual_model_forecast(m, :mode, :none, [:histobs, :forecastobs, :histpseudo, :forecastpseudo];
                         forecast_string = "short", params = [x.value for x in m.parameters],
                         check_empty_columns = false)
    mbobs_short = read_mb(m, :mode, :none, :forecastobs; forecast_string = "short")
    mbpse_short = read_mb(m, :mode, :none, :forecastpseudo; forecast_string = "short")

    reg_2021Q1 = DSGE.subtract_quarters(Date(2021, 3, 31), Date(2020, 9, 30)) + 2
    # m <= Setting(:tvis_replace_eqcond_func_dict, [replace_eqcond_func_dict_shortzlb, replace_eqcond_func_dict_longzlb])
    m <= Setting(:tvis_regime_eqcond_info, [regime_eqcond_info_shortzlb, regime_eqcond_info_longzlb])
    m <= Setting(:tvis_select_system, vcat(ones(Int, reg_2021Q1), fill(2, get_setting(m, :n_regimes) - reg_2021Q1)))
    m <= Setting(:tvis_information_set, vcat([1:1], [i:reg_2024Q1 for i in 2:reg_2021Q1],
                                             [i:get_setting(m, :n_regimes) for i in (reg_2021Q1 + 1):get_setting(m, :n_regimes)]))
    usual_model_forecast(m, :mode, :none, [:histobs, :forecastobs, :histpseudo, :forecastpseudo];
                         forecast_string = "shortlong", params = [x.value for x in m.parameters],
                         check_empty_columns = false)
    mbobs_shortlong = read_mb(m, :mode, :none, :forecastobs; forecast_string = "shortlong")
    mbpse_shortlong = read_mb(m, :mode, :none, :forecastpseudo; forecast_string = "shortlong")

    if create_plots
        plot_dicts = Dict()
        dates = [DSGE.iterate_quarters(d, -1) for d in mbobs_short.means[!, :date]]

        plot_dicts[:obs_gdp] = plot()
        plot!(dates, mbobs_shortlong.means[!, :obs_gdp], color = :black, label = "2024:Q4 ZLB", linewidth = 3)
        plot!(dates, mbobs_short.means[!, :obs_gdp], color = :blue, label = "2023:Q4 ZLB", linewidth = 3)

        plot_dicts[:obs_corepce] = plot()
        plot!(dates, mbobs_shortlong.means[!, :obs_corepce], color = :black, label = "2024:Q4 ZLB", linewidth = 3)
        plot!(dates, mbobs_short.means[!, :obs_corepce], color = :blue, label = "2023:Q4 ZLB", linewidth = 3)

        plot_dicts[:obs_nominalrate] = plot()
        plot!(dates, mbobs_shortlong.means[!, :obs_nominalrate], color = :black, label = "2024:Q4 ZLB", linewidth = 3)
        plot!(dates, mbobs_short.means[!, :obs_nominalrate], color = :blue, label = "2023:Q4 ZLB", linewidth = 3)

        plot_dicts[:obs_longinflation] = plot()
        plot!(dates, mbobs_shortlong.means[!, :obs_longinflation], color = :black, label = "2024:Q4 ZLB", linewidth = 3)
        plot!(dates, mbobs_short.means[!, :obs_longinflation], color = :blue, label = "2023:Q4 ZLB", linewidth = 3)

        plot_dicts[:obs_spread] = plot()
        plot!(dates, mbobs_shortlong.means[!, :obs_spread], color = :black, label = "2024:Q4 ZLB", linewidth = 3)
        plot!(dates, mbobs_short.means[!, :obs_spread], color = :blue, label = "2023:Q4 ZLB", linewidth = 3)

        plot_dicts[:NaturalRate] = plot()
        plot!(dates, mbpse_shortlong.means[!, :NaturalRate], color = :black, label = "2024:Q4 ZLB", linewidth = 3)
        plot!(dates, mbpse_short.means[!, :NaturalRate], color = :blue, label = "2023:Q4 ZLB", linewidth = 3)
    end
end


text:
using AllTools
using Test

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


text:
using DataFrames
using CSV
using Statistics
using Base.Threads: @threads
using StatsBase: sample
using StaticArrays
using Revise
using BenchmarkTools
using EvoTrees
using EvoTrees: get_gain, get_edges, binarize, get_max_gain, update_grads!, grow_tree, grow_gbtree, SplitInfo, Tree, TrainNode, TreeNode, EvoTreeRegressor, predict, predict!, sigmoid
using EvoTrees: find_bags, update_bags!, find_split_static!, pred_leaf, sigmoid, logit
using Plots
using Distributions

# prepare a dataset
features = rand(10_000) .* 5
X = reshape(features, (size(features)[1], 1))
Y = sin.(features) .* 0.5 .+ 0.5
Y = logit(Y) + randn(size(Y))
Y = sigmoid(Y)
𝑖 = collect(1:size(X,1))
𝑗 = collect(1:size(X,2))

# train-eval split
𝑖_sample = sample(𝑖, size(𝑖, 1), replace = false)
train_size = 0.8
𝑖_train = 𝑖_sample[1:floor(Int, train_size * size(𝑖, 1))]
𝑖_eval = 𝑖_sample[floor(Int, train_size * size(𝑖, 1))+1:end]

X_train, X_eval = X[𝑖_train, :], X[𝑖_eval, :]
Y_train, Y_eval = Y[𝑖_train], Y[𝑖_eval]
𝑖 = collect(1:size(X_train,1))

# set parameters
params1 = EvoTreeRegressor(
    loss=:gaussian, metric=:gaussian,
    nrounds=100, nbins=100,
    λ = 0.0, γ=0.0, η=0.1,
    max_depth = 6, min_weight = 1.0,
    rowsample=0.5, colsample=1.0, seed=123)

mean(Y_train.^2)

# initial info
@time δ, δ² = zeros(SVector{params1.K, Float64}, size(X_train, 1)), zeros(SVector{params1.K, Float64}, size(X_train, 1))
𝑤 = zeros(SVector{1, Float64}, size(X_train, 1)) .+ 1
pred = zeros(SVector{params1.K,Float64}, size(X_train,1))
@time update_grads!(params1.loss, params1.α, pred, Y_train, δ, δ², 𝑤)
∑δ, ∑δ², ∑𝑤 = sum(δ[𝑖]), sum(δ²[𝑖]), sum(𝑤[𝑖])
@time gain = get_gain(params1.loss, ∑δ, ∑δ², ∑𝑤, params1.λ)
# @btime gain = get_gain($params1.loss, $∑δ, $∑δ², $∑𝑤, $params1.λ)

# initialize train_nodes
train_nodes = Vector{TrainNode{params1.K, Float64, BitSet, Array{Int64, 1}, Int}}(undef, 2^params1.max_depth-1)
for node in 1:2^params1.max_depth-1
    train_nodes[node] = TrainNode(0, SVector{params1.K, Float64}(fill(-Inf, params1.K)), SVector{params1.K, Float64}(fill(-Inf, params1.K)), SVector{1, Float64}(fill(-Inf, 1)), -Inf, BitSet([0]), [0])
    # train_nodes[feat] = TrainNode(0, fill(-Inf, params1.K), fill(-Inf, params1.K), -Inf, -Inf, BitSet([0]), [0])
end

# initializde node splits info and tracks - colsample size (𝑗)
splits = Vector{SplitInfo{params1.K, Float64, Int}}(undef, size(𝑗, 1))
for feat in 1:size(𝑗, 1)
    splits[feat] = SplitInfo{params1.K, Float64, Int}(gain, SVector{params1.K, Float64}(zeros(params1.K)), SVector{params1.K, Float64}(zeros(params1.K)), SVector{1, Float64}(zeros(1)), SVector{params1.K, Float64}(zeros(params1.K)), SVector{params1.K, Float64}(zeros(params1.K)), SVector{1, Float64}(zeros(1)), -Inf, -Inf, 0, feat, 0.0)
end

# binarize data and create bags
@time edges = get_edges(X_train, params1.nbins)
@time X_bin = binarize(X_train, edges)
@time bags = Vector{Vector{BitSet}}(undef, size(𝑗, 1))
function prep(X_bin, bags)
    @threads for feat in 1:size(𝑗, 1)
         bags[feat] = find_bags(X_bin[:,feat])
    end
    return bags
end
@time bags = prep(X_bin, bags)

# initialize histograms
feat=1
hist_δ = Vector{Vector{SVector{params1.K, Float64}}}(undef, size(𝑗, 1))
hist_δ² = Vector{Vector{SVector{params1.K, Float64}}}(undef, size(𝑗, 1))
hist_𝑤 = Vector{Vector{SVector{1, Float64}}}(undef, size(𝑗, 1))
for feat in 1:size(𝑗, 1)
    hist_δ[feat] = zeros(SVector{params1.K, Float64}, length(bags[feat]))
    hist_δ²[feat] = zeros(SVector{params1.K, Float64}, length(bags[feat]))
    hist_𝑤[feat] = zeros(SVector{1, Float64}, length(bags[feat]))
end

# grow single tree
#  0.135954 seconds (717.54 k allocations: 15.219 MiB)
@time train_nodes[1] = TrainNode(1, SVector(∑δ), SVector(∑δ²), SVector(∑𝑤), gain, BitSet(𝑖), 𝑗)
@time tree = grow_tree(bags, δ, δ², 𝑤, hist_δ, hist_δ², hist_𝑤, params1, train_nodes, splits, edges, X_bin)
# @btime tree = grow_tree($bags, $δ, $δ², $𝑤, $hist_δ, $hist_δ², $hist_𝑤, $params1, $train_nodes, $splits, $tracks, $edges, $X_bin)
@time pred_train = predict(tree, X_train, params1.K)
# 705.901 μs (18 allocations: 626.08 KiB)
# @btime pred_train = predict($tree, $X_train, $params1.K)
@time pred_leaf_ = pred_leaf(params1.loss, train_nodes[1], params1, δ²)
# 1.899 ns (0 allocations: 0 bytes)
# @btime pred_leaf_ = pred_leaf($params1.loss, $train_nodes[1], $params1, $δ²)
# @btime pred_train = predict($tree, $X_train, params1.K)

@time model = grow_gbtree(X_train, Y_train, params1, print_every_n = 20)
# @btime model = grow_gbtree($X_train, $Y_train, $params1, print_every_n = 1)
@time pred_train = predict(model, X_train)
# @btime pred_train = predict($model, $X_train)

x_perm = sortperm(X_train[:,1])
plot(X_train, Y_train, ms = 1, mcolor = "gray", mscolor = "lightgray", background_color = RGB(1, 1, 1), seriestype=:scatter, xaxis = ("feature"), yaxis = ("target"), legend = true, label = "")
plot!(X_train[:,1][x_perm], pred_train[:,1][x_perm], color = "navy", linewidth = 1.5, label = "Median")
# σ²
plot!(X_train[:,1][x_perm], sqrt.(pred_train[:,2][x_perm]), color = "red", linewidth = 1.5, label = "sigma")
# q20
dist = Normal.(pred_train[:,1][x_perm], sqrt.(pred_train[:,2][x_perm]))
pred_train_q20 = quantile.(dist, 0.2)
pred_train_q80 = quantile.(dist, 0.8)
sum(pred_train_q20 .< Y_train[x_perm]) / size(Y_train,1)
sum(pred_train_q80 .< Y_train[x_perm]) / size(pred_train,1)
plot!(X_train[:,1][x_perm], pred_train_q20, color = "green", linewidth = 1.5, label = "Q20")
plot!(X_train[:,1][x_perm], pred_train_q80, color = "green", linewidth = 1.5, label = "Q80")
savefig("gaussian_likelihood.png")


text:
function absorption_line(x::T; mid=zero(T), width=one(T), depth=one(T)) where T<:AF
    return one(T) - depth * exp(-((x-mid)/width)^2.0/2.0)
end

function gaussian_line(x::T; mid=zero(T), width=one(T), depth=one(T)) where T<:AF
    return absorption_line(x, mid=mid, width=width, depth=depth)
end

function lorentzian_line(x::T; mid=zero(T), width=one(T), depth=one(T)) where T<:AF
    y = (x - mid)/(width/2)
    return one(T) - (depth / (one(T) + y^2))
end

# calculate doppler factor
function doppler_factor(vel::T) where T<:AF
    num = one(T) + vel/c_ms
    den = one(T) - vel/c_ms
    return sqrt(num/den)
end

# Width of line due to thermal broadening
function width_thermal(; λ::T1=1.0, M::T1=1.0, T::T1=5778.0, v_turb::T1=0.0) where T1<:AF
    return sqrt((2.0*kB/mH) * (T/M) + v_turb^2) * (λ/c)
end

# Quadratic limb darkening law.
# Takes μ = cos(heliocentric angle) and LD parameters, u1 and u2.
function quad_limb_darkening(μ::T, u1::T, u2::T) where T<:AF
    return !iszero(μ) * (one(T) - u1*(one(T)-μ) - u2*(one(T)-μ)^2)
end

function quad_limb_darkening(x::T, y::T, u1::T, u2::T) where T<:AF
    return quad_limb_darkening(calc_mu(x,y), u1, u2)
end

"""
    rotation_period(sin_lat; A=14.713, B=2.396, C=1.787)

Calculate stellar rotation period in days at given sine of latitude.

# Arguments
- `sin_lat::Float64`: Sine of stellar latitude.
- `A::Float64=14.713`: Coefficient from Snodgrass & Ulrich (1990)
- `B::Float64=2.396`: Coefficient from Snodgrass & Ulrich (1990)
- `C::Float64=1.787`: Coefficient from Snodgrass & Ulrich (1990)
"""
function rotation_period(sin_lat::T; A::T=14.713, B::T=-2.396, C::T=-1.787) where T<:AF
    @assert -1.0 <= sin_lat <= 1.0
    return 360.0/(A + B * sin_lat^2 + C * sin_lat^4)
end

"""
    patch_velocity_los(x, y; rstar, pole)

Compute line of sight velocity of a patch of stellar surface given by x,y (assumed in [-1,1]).
Return value is in (Rsol/day)/speed of light (i.e., dimensionless like z = v/c)

# Arguments
- `rstar::Float64=1.0`: in R_sol (affects return velocity, but not x,y)
- `pole = (0,1,0)`: unit vector for stellar rotation axis (default is equator-on)
"""
function patch_velocity_los(x::T,y::T; rstar=one(T), pole=(zero(T), one(T), zero(T))) where T<:AF
    polex, poley, polez = pole
    v0 = 0.000168710673 # in (Rsol/day)/speed_of_light
    z = sqrt(rstar - calc_r2(x,y))
    sin_lat = (x*polex) + (y*poley) + (z*polez)
    vmax = v0*rstar/rotation_period(sin_lat)
    return vmax * (polex*y - poley*x)
end

function patch_velocity_los(t::Tuple{T,T}; rstar=one(T), pole=(zero(T), one(T), zero(T))) where T<:AF
    return patch_velocity_los(t[1], t[2], rstar=rstar, pole=pole)
end

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
@testset "LibGEOS regressions" begin

    # https://github.com/JuliaGeo/LibGEOS.jl/issues/29
    pts = [[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0]]
    @test_throws LibGEOS.GEOSError Polygon([pts])

    # https://github.com/JuliaGeo/LibGEOS.jl/issues/24
    point = Point(2, 3)
    @test GeoInterface.geotype(point) == :Point

    # https://github.com/JuliaGeo/LibGEOS.jl/issues/25
    a = LibGEOS.createCoordSeq([0.0, 2.0])
    @test LibGEOS.getCoordinates(a) == Array{Float64,1}[[0.0, 2.0]]

    # https://github.com/JuliaGeo/LibGEOS.jl/issues/12
    mp = LibGEOS.MultiPoint(Vector{Float64}[[0, 0], [10, 0], [10, 10], [11, 10]])
    @test GeoInterface.geotype(mp) == :MultiPoint

    # https://github.com/JuliaGeo/LibGEOS.jl/issues/20
    # LibGEOS doesn't support Extended WKT
    ewkt = "SRID=32756; POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10))"
    # ParseException: Unknown type: 'SRID=32756;'
    @test_throws LibGEOS.GEOSError readgeom(ewkt)

    # https://github.com/JuliaGeo/LibGEOS.jl/issues/40
    # IllegalArgumentException: Points of LinearRing do not form a closed linestring
    @test_throws LibGEOS.GEOSError readgeom("POLYGON((-1. 1., 1. 3., 0. 2., -1. 1.5))")

end


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e 'include(popfirst!(ARGS))' \
    "${BASH_SOURCE[0]}" "$@"
    =#

# setprecision(100000)
import Base.π # needed in order to redefine it

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

# Dot minus (monus):
# i.e., truncated negation
@doc raw"""
Truncated negation: an operator that acts like subtraction without going to negatives.

This can be written in the REPL as `\dotminus<tab>`.

Examples:
```julia
julia> 6 ∸ 3
3

julia> 3 ∸ 1
2

julia> 10 ∸ 12
0

julia> 1 ∸ 2
0
```
"""
function ∸(x::Number, y::Number)
    z = x - y
    return z > 0 ? z : 0
end

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

pairntuple_error = "This function is only defined for natural numbers.  Use cℤ."

@doc raw"""
```julia
pair_tuple(x::Integer, y::Integer)
pair_tuple(x::Integer, y::Integer, z::Integer...)
pair_tuple(t::Tuple)
```

The pairing function is a function that uniquely maps any number of integer inputs to a single integer:

```math
\left(k_1,k_2\right) \mapsto \frac{1}{2}\left(k_1+k_2\right)\left(k_1+k_2+1\right)+k_1\\
\mathbb{N}^2 \to \mathbb{N}
```

This concept can also be generalised to take integers ``x_0`` to ``x_{n-1}`` so that

```math
(x_0, x_1, \ldots, x_{n - 1}) \mapsto \left\langle\left\langle x_0, \ldots, x_{n - 2}\right\rangle, x_{n - 1}\right\rangle\\
\mathbb{N}^n \to \mathbb{N}
```

and returns their "pair".

As this maps every combination of integers to a single output, this number gets massive very fast.  As such, `pair_tuple` will return a `BigInt`.

To depair the tuple, see ``\pi``.

For more information, see [https://en.m.wikipedia.org/wiki/Pairing_function](https://en.m.wikipedia.org/wiki/Pairing_function).

!!! note

    There are two variations to this process **which will change the output**:
      - As per the visualisation in the link above, the first number given is found along the ``x`` axis, and the second along the ``y`` axis.  **In this implementation, it is the other way around.**  This can easily by modified, however, by adding ``k_1`` at the end instead of ``k_2``.
      - The other variation is which direction of nesting we decide on for more than 2 inputs.  The alternative (which **we do not use in this implementation**) is
        ```math
        \left(x_0, x_1, \ldots, x_{n - 1}\right) \mapsto \left\langle x_0, \left\langle x_1, \ldots, x_{n - 1}\right\rangle\right\rangle
        ```

---

### Examples

```julia
julia> pair_tuple(5,7) # code pair of natural numbers as a natural number
83
```
"""
pair_tuple(x::Integer, y::Integer) = x < 0 || y < 0 ? throw(error("$pairntuple_error")) : big(x) + binomial(big(x)+big(y)+1, 2)
pair_tuple(x::Integer, y::Integer, z::Integer...) = pair_tuple(pair_tuple(x, y), z...)
pair_tuple(t::Tuple) = pair_tuple(t...)

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

# Singleton type with a descriptive name, and use that as the second argument
# In the algebraic unpairing function
# struct Algebra end
# const algebraic = Algebra()

# The brute-force unpairing function:
# Takes in integer m (that is, a natural number)
# and integer n (that is, m ⟼ <x1, ..., xn>)
# returns x1, ..., xn.
# Alternatively, given
@generated function _π_brute_force(::Val{n}, m::Integer) where {n}
    quote
        iszero(n) && return nothing
        isone(n) && return m

        # promote type to BigInt
        # we probably do not need n to be a big int
        m = big(m)

        @inbounds @fastmath Base.Cartesian.@nloops $n i d -> 0:m begin
            if isequal(pair_tuple((Base.Cartesian.@ntuple $n i)...), m)
                return Base.Cartesian.@ntuple $n i
            end
        end
    end
end
# Ensuring the function is more readible
# (i.e., switch the inputs) such that
# π(m, n) ⟼ <x1, ..., xn> = m
_π_brute_force(m::Integer, n::Integer) = _π_brute_force(Val(n), m)
# default to n=2
_π_brute_force(m::Integer) = _π_brute_force(m, 2)
# Defining a selection function that obtains
# the kth element in the tuple obtained using π;
# e.g., π(83, 2, 0) = 5 = \pi_2^0(83)
# note: by convention this is indexed from zero,
# which is why we need to offset it by one.
_π_brute_force(m::Integer, n::Integer, k::Integer) = _π_brute_force(m, n)[k+1]

# Algebraic unpairing function:
# Unpairing if n=2 (base case):
function _π(m::Integer)
    m = big(m)

    w = (Base.isqrt(8*m + 1) - 1) ÷ 2
    t = (w^2 + w) ÷ 2
    x = m - t
    y = w - x

    if ! isequal(pair_tuple(x, y), m)
        throw(error("The provided m = $m is not equal to ⟨ $x, $y ⟩, and so there has been an error in the calculation."))
    end

    return BigInt(x), BigInt(y)
end

# Generalised case:
function _π(m::Integer, n::Integer)
    iszero(n) && return nothing
    isone(n) && return m

    appended_tuple = _π(m)

    while n != 2
        appended_tuple = (_π(appended_tuple[1])..., appended_tuple[2:end]...)
        n -= 1
    end

    return appended_tuple
end
# Selection function taking input k and outputting the
# kth element in the tuple obtained by the algebraic π
_π(m::Integer, n::Integer, k::Integer) = _π(m, n)[k+1]

@doc raw"""
```julia
π(m::Integer..., algebraic::Bool = true)
```

The depairing function.  Takes in an integer and find its unique pair of integers such that those integers paired is exactly the function input:

```math
\pi(m) \mapsto \left\langle m_0, m_1\right\rangle\\
\mathbb{N} \to \mathbb{N}^2
```

Setting boolean `algebraic` flag to `false` calls a brute force search method which is extremely slow for large ``n`` (see the `n` parameter below) but exists for completeness.  That is,

```math
\pi(m) := \left(\mu x\leq m\right)\left(\exists y\leq m\right)m=\left\langle x,y\right\rangle\\
\pi_0\left(\left\langle x_0,x_1\right\rangle\right) = x_0\\
\pi_1\left(\left\langle x_0,x_1\right\rangle\right) = x_1
```

This function also takes in a parameter ``n``:

```math
\pi(m, n) \mapsto \left\langle m_1, \ldots, m_n\right\rangle\\
\mathbb{N} \to \mathbb{N}^n
```

This function may also take in a third parameter ``k`` so that

```math
\pi(m, n, k) \equiv \pi_k^n \mapsto \left\langle m_1, \ldots, m_n\right\rangle \mapsto m_k\quad\text{where }1\leq k\leq n
```

!!! note

    See the note from documentation on `pair_tuple`.  This function may not do exactly what you expect.

---

### Examples

```julia
julia> π(83, algebraic = false) # code a natural number into a 2-tuple
(5, 7)

julia> π(83, 2, algebraic = true) # use algebraic method of depairing rather than search (much faster)
(5, 7)

julia> π(83, 2, 1, algebraic = false) # code a natural number into a 2-tuple and get the the number in the tuple indexed by 1 (index starting from zero)
7
```
"""
π(m::Integer...; algebraic::Bool = true) =
    algebraic ? _π(m...) : _π_brute_force(m...)
depair(m::Integer...; algebraic::Bool = true) = π(m..., algebraic = algebraic)

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

@doc raw"""
```julia
cℤ(z::Integer)
cℤ(z::Integer, w::Integer...)
```

This function uniquely maps an integer to a natural number.  The output is a BigInt.  This function can be written in the REPL by `c\bbZ<tab>`.

```math
z \mapsto c\mathbb{Z}(z)\\
z, w, \ldots \mapsto c\mathbb{Z}(z), c\mathbb{Z}(w), \ldots\\
\mathbb{Z} \to \mathbb{N}
```

---

```julia
cℤ(zs::AbstractArray{<:Integer})
```

Given an array of integers, `cℤ` will convert every element of the array into natrual numbers.

---

```julia
cℤ(r::NTuple{N, Integer})
```

Given a tuple, `cℤ` will convert ever element of the tuple into natural numbers, and pair the result:

```math
\left(z, w\right) \mapsto c\mathbb{Z}\left(\left\langle z, w\right\rangle\right)\\
\mathbb{Z}^2 \to \mathbb{N}
```
```math
\left(z_1, \ldots, z_n\right) \mapsto c\mathbb{Z}\left(\left\langle z_1\right\rangle\right), \ldots, c\mathbb{Z}\left(\left\langle z_n\right\rangle\right)\\
\mathbb{Z}^n \to \mathbb{N}
```

---

### Examples

```julia
julia> cℤ(-10) # code integer as a natural number
19

julia> cℤ((-1,2)) # cℤ pairs the code of each
16
```
"""
cℤ(z::Integer) = z >= 0 ? 2 * big(z) : (2 * abs(big(z))) - 1
cℤ(z::Integer, w::Integer...) = cℤ(z), cℤ(w...)...
cℤ(zs::AbstractArray{<:Integer}) = cℤ.(zs)
cℤ(r::NTuple{N, Integer}) where {N} = pair_tuple(cℤ.(r)...)

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

cℤ⁻¹_error = "Invalid input. We have only defined this function for natural numbers.  Why are you even using it?"

@doc raw"""
```julia
cℤ⁻¹(n::Integer)
```

The inverse of `cℤ⁻¹`.  Uniquely maps natural numbers to integers.  This function can be written in the REPL by `c\bbZ<tab>\^1<tab>`.

---

```julia
cℤ⁻¹(zs::AbstractArray{<:Integer})
```

Given an array of integers, `cℤ⁻¹` will convert every element of the array into integers.

---

### Examples

```julia
julia> cℤ⁻¹(19)
-10
```
"""
cℤ⁻¹(n::Integer) = n < 0 ? throw(error("$cℤ⁻¹_error")) : (iseven(n) ? Int(big(n) / 2) : -Int(floor(big(n) / 2) + 1))
cℤ⁻¹(ns::AbstractArray{<:Integer}) = cℤ⁻¹.(ns)


text:
# Autogenerated wrapper script for Arb_jll for i686-w64-mingw32
export libarb

using FLINT_jll
using GMP_jll
using MPFR_jll
JLLWrappers.@generate_wrapper_header("Arb")
JLLWrappers.@declare_library_product(libarb, "libarb-2.dll")
function __init__()
    JLLWrappers.@generate_init_header(FLINT_jll, GMP_jll, MPFR_jll)
    JLLWrappers.@init_library_product(
        libarb,
        "bin\\libarb-2.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()