text:
using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils, ..PrepareArrtoBool
using ..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates

@testset "getBoolKernel" begin

@testset "getIndexOfQueue" begin
        """
    1)   Left FP
    2)   Left FN
    3)   Right FP
    4)   Right FN
    5)   Posterior FP
    6)   Posterior FN
    7)   Anterior FP
    8)   Anterior FN
    9)   Top FP
    10)   Top FN
    11)   Bottom FP
    12)   Bottom FN
    13)   Total block Fp
    14)   Total block Fn

    """
    dataBdim= (5,5,5)
    @test getIndexOfQueue(1,1,1, (5,5,5),false)==1
    @test getIndexOfQueue(1,1,1, (5,5,5),true)==2

    @test getIndexOfQueue(5,1,1, (5,5,5),false)==3
    @test getIndexOfQueue(5,1,1, (5,5,5),true)==4

    @test getIndexOfQueue(2,1,1, (5,5,5),false)==5
    @test getIndexOfQueue(2,1,1, (5,5,5),true)==6

    @test getIndexOfQueue(2,5,1, (5,5,5),false)==7
    @test getIndexOfQueue(2,5,1, (5,5,5),true)==8

    @test getIndexOfQueue(2,2,1, (5,5,5),false)==9
    @test getIndexOfQueue(2,2,1, (5,5,5),true)==10

    @test getIndexOfQueue(2,2,5, (5,5,5),false)==11
    @test getIndexOfQueue(2,2,5, (5,5,5),true)==12

    @test getIndexOfQueue(5,5,5, (5,5,5),false)==3
    @test getIndexOfQueue(5,5,5, (5,5,5),true)==4

    @test getIndexOfQueue(2,2,2, (5,5,5),false)==13
    @test getIndexOfQueue(2,2,2, (5,5,5),true)==14

end

# metadat= zeros(Int64,100,50,40,15);
# metaDataDims= size(metadat)
# for i in 1:length(metadat)
#     metadat[i]=i
# end
# linIdexMetaPrim= (100*50*3)+303
# linIdexMeta= (-1)+ (100*50*3)+303
# xMeta= mod(linIdexMeta,metaDataDims[1])
# zMeta= fld((linIdexMeta),metaDataDims[1]*metaDataDims[2])
# yMeta= fld((linIdexMeta-((zMeta*metaDataDims[1]*metaDataDims[2] ) + xMeta )),metaDataDims[1])
# qn= 7
# metadat[linIdexMetaPrim+(qn-1)*metaDataDims[1]*metaDataDims[2]*metaDataDims[3] ]
# metadat[xMeta+1,yMeta+1,zMeta+1,qn]

##### iter3dOuter

@testset "iter3dOuter" begin

        using Revise, Parameters, Logging, Test
        using CUDA
        includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
        using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
        using ..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates

    singleVal = CUDA.zeros(14)

    # threads=(32,20)

    threads=(32,8)
    # blocks =2

    # mainArrDims= (7,7,7)
    # dataBdim = (2,2,2)

    blocks =136
    dataBdim = (33,32,32)
    mainArrDims= (512, 512, 500)

    metaDataDims= (cld(mainArrDims[1],dataBdim[1] ),cld(mainArrDims[2],dataBdim[2]),cld(mainArrDims[3],dataBdim[3]))
    #we are iterating here block by block sequentially
    metaDataDims= (16, 16, 16,90)
    metaDataLength= metaDataDims[1]*metaDataDims[2]*metaDataDims[3]
    loopMeta= fld(metaDataLength,blocks )

    # metaDataDims (16, 16, 16, 90)  loopMeta 409 metaDataLength 4096
    function iter3dOuterKernel(mainArrDims,singleVal,metaDataDims,dataBdim, metaDataLength, loopMeta )
        @iter3dOuter(metaDataDims,loopMeta,metaDataLength,
        begin
        @ifXY 1 1 CUDA.@atomic singleVal[]+=1
     @ifXY 1 1    CUDA.@cuprint " linIdexMeta $(linIdexMeta)  zMeta $(zMeta)  yMeta$(yMeta) xMeta$(xMeta)  idX $(blockIdxX()) \n"

    end)

        return
    end
    @cuda threads=threads blocks=blocks iter3dOuterKernel(mainArrDims,singleVal,metaDataDims,dataBdim, metaDataLength, loopMeta )
    @test singleVal[1]==metaDataDims[1]*metaDataDims[2]*metaDataDims[3]

end #"iter3dOuter"

################ iter data block
@testset "iter data block" begin

    using Revise, Parameters, Logging, Test
    using CUDA
    includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
    using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
    using ..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates

    ##### iter data block
    singleVal = CUDA.zeros(Int64,14)
    threads=(32,5)

    blocks =17
    mainArrDims= (67,177,90)
    dataBdim = (43,21,17)
    indices = CUDA.zeros(Bool,mainArrDims)

    metaDataDims= (cld(mainArrDims[1],dataBdim[1] ),cld(mainArrDims[2],dataBdim[2]),cld(mainArrDims[3],dataBdim[3]))

    inBlockLoopX,inBlockLoopY,inBlockLoopZ= (fld(dataBdim[1] ,threads[1]),fld(dataBdim[2] ,threads[2]),dataBdim[3]    )
    #we are iterating here block by block sequentially
    metaDataLength= metaDataDims[1]*metaDataDims[2]*metaDataDims[3]
    loopMeta= fld(metaDataLength,blocks )

    function iterDataBlocksKernel(loopMeta,metaDataLength,mainArrDims,singleVal,metaDataDims,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ ,indices)
        @iter3dOuter(metaDataDims,loopMeta,metaDataLength,
        begin
            @iterDataBlock(mainArrDims,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ,xMeta,yMeta,zMeta,
            begin
           CUDA.@atomic singleVal[1]+=1
            indices[x,y,z]=true
            #CUDA.@cuprint "x $(x) y $(y) z $(z) zMeta $(zMeta)  yMeta $(yMeta) xMeta $(xMeta)  idX $(blockIdxX()) \n"

        end)end)

        return
    end
    @cuda threads=threads blocks=blocks iterDataBlocksKernel(loopMeta,metaDataLength,mainArrDims,singleVal,metaDataDims,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ ,indices)
    @test singleVal[1]==mainArrDims[1]*mainArrDims[2]*mainArrDims[3]
    Int64(singleVal[1])
    7*7*7
    inddd = Array(indices);
    cartt = CartesianIndices(inddd);
    filtered = filter(cc-> !inddd[cc]  ,cartt)
    sum(indices)

end #iter data block

##### @uploadLocalfpFNCounters
@testset "uploadLocalfpFNCounters" begin

    localQuesValues = CUDA.zeros(UInt32,14)
    threads=(32,5)
    blocks =1

    mainArrDims= (99,99,99)
    dataBdim = (17,7,12)

    metaDataDims= (cld(mainArrDims[1],dataBdim[1] ),cld(mainArrDims[2],dataBdim[2]),cld(mainArrDims[3],dataBdim[3]))

    inBlockLoopX,inBlockLoopY,inBlockLoopZ= (fld(dataBdim[1] ,threads[1]),fld(dataBdim[2] ,threads[2]),dataBdim[3]    )
    #we are iterating here block by block sequentially
    metaDataLength= metaDataDims[1]*metaDataDims[2]*metaDataDims[3]
    loopMeta= fld(metaDataLength,blocks )

    function uploadLocalfpFNCountersKernel(mainArrDims,localQuesValues,metaDataDims,loopMeta,metaDataLength,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ)
        # localQuesValues= @cuStaticSharedMem(Float32, 14)

        # @iter3dOuter(metaDataDims,loopMeta,metaDataLength,
        # begin
            #@ifXY 1 1    CUDA.@cuprint "  xMeta $(xMeta) yMeta $(yMeta)  zMeta $(zMeta) \n"
            xMeta= 2
            zMeta= 2
            yMeta= 2

            @iterDataBlock(mainArrDims,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ,xMeta,yMeta,zMeta,
            begin
                        boolGold=true
                        boolSegm=false
                        #coord=PrepareArrtoBool.getIndexOfQueue((xdim * blockDimX())+threadIdxX() ,(ydim * blockDimY())+threadIdxY(),(zdim+1),dataBdim,false)
                        PrepareArrtoBool.@uploadLocalfpFNCounters()
                        #CUDA.@cuprint " spot $(coord)  x $((xdim * blockDimX())+threadIdxX()  )) y $((ydim * blockDimY())+threadIdxY()  )) z $((zdim+1)) dataBdim $(dataBdim[1]),$(dataBdim[2]),$(dataBdim[3]) \n"
                        #CUDA.atomic_add!(pointer(localQuesValues, coord),Float32(1))

                    #@inboundsCUDA.@atomic localQuesValues[coord]+=Float32(1)
                    #CUDAAtomicUtils.atomicallyAddToSpot(Float32,localQuesValues,getIndexOfQueue(x,y,z,dataBdim,boolSegm),1)
                    #PrepareArrtoBool.@uploadLocalfpFNCounters()

                    end)#end)

        return
    end

    @cuda threads=threads blocks=blocks uploadLocalfpFNCountersKernel(mainArrDims,localQuesValues,metaDataDims,loopMeta,metaDataLength,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ)
    @test sum(localQuesValues)>0

    @test Int64(localQuesValues[1])==dataBdim[2]*dataBdim[3]
    @test Int64(localQuesValues[3])==dataBdim[2]*dataBdim[3]

    @test Int64(localQuesValues[5])==dataBdim[1]*dataBdim[3]-2*dataBdim[3]
    @test Int64(localQuesValues[7])==dataBdim[1]*dataBdim[3]-2*dataBdim[3]

    @test Int64(localQuesValues[9])==(dataBdim[1]-2)*(dataBdim[2]-2)
    @test Int64(localQuesValues[11])==(dataBdim[1]-2)*(dataBdim[2]-2)

    @test Int64(localQuesValues[13])==dataBdim[1]*dataBdim[2]*dataBdim[3] - sum(localQuesValues[1:12])

end #"uploadLocalfpFNCounters"

using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
using ..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates

# ########### uploadDataToMetaData
#  @testset "uploadDataToMetaData 1 " begin

#     localQuesValues = CUDA.zeros(Float32,14)
#     threads=(32,5)
#     blocks =8
#     mainArrDims= (634,521,632)
#     dataBdim = (43,21,17)

#     metaDataDims= (cld(mainArrDims[1],dataBdim[1] ),cld(mainArrDims[2],dataBdim[2]),cld(mainArrDims[3],dataBdim[3]))
#     loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength = calculateLoopsIter(dataBdim,threads[1],threads[2],metaDataDims,blocks)

#     #we are iterating here block by block sequentially
#     metaData = MetaDataUtils.allocateMetadata(mainArrDims,dataBdim);
#     metaDataLength= metaDataDims[1]*metaDataDims[2]*metaDataDims[3]
#     loopMeta= fld(metaDataLength,blocks )
#     function uploadDataToMetaDataKernel(loopMeta,metaDataLength,mainArrDims,metaDataDims,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaData)
#         localQuesValues= @cuStaticSharedMem(UInt32, 14)
#         localBool=false
#         anyPositive=false

#         @iter3dOuter(metaDataDims,loopMeta,metaDataLength,
#         begin

#             @iterDataBlock(mainArrDims,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ,
#             begin
#                         boolGold=true
#                         boolSegm=false
#                         localBool=true
#                         if(xMeta==2 && yMeta==2 && zMeta==4)

#                             PrepareArrtoBool.@uploadLocalfpFNCounters()
#                         end

#                     end)
#                     sync_threads()

#                     isAnyPositive= true
#                         PrepareArrtoBool.@uploadDataToMetaData()
#                         sync_threads()

#                         sync_threads()
#                         #threadfence()
#                         @ifY 1 if(threadIdxX()<15)
#                             localQuesValues[threadIdxX()]=0
#                         end
#                 end)

#         return
#     end

#     @cuda threads=threads blocks=blocks uploadDataToMetaDataKernel(loopMeta,metaDataLength,mainArrDims,metaDataDims,loopXMeta,loopYZMeta,yTimesZmeta,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaData)

#     @test metaData[3,3,5,getBeginingOfFpFNcounts()+1]==dataBdim[2]*dataBdim[3]
#     @test metaData[3,3,5,getBeginingOfFpFNcounts()+3]==dataBdim[2]*dataBdim[3]

#     @test metaData[3,3,5,getBeginingOfFpFNcounts()+5]==dataBdim[1]*dataBdim[3]-2*dataBdim[3]
#     @test metaData[3,3,5,getBeginingOfFpFNcounts()+7]==dataBdim[1]*dataBdim[3]-2*dataBdim[3]

#     @test metaData[3,3,5,getBeginingOfFpFNcounts()+9]==(dataBdim[1]-2)*(dataBdim[2]-2)
#     @test metaData[3,3,5,getBeginingOfFpFNcounts()+11]==(dataBdim[1]-2)*(dataBdim[2]-2)

#     @test metaData[3,3,5,getBeginingOfFpFNcounts()+13]==dataBdim[1]*dataBdim[2]*dataBdim[3] - sum(metaData[3,3,5,getBeginingOfFpFNcounts():getBeginingOfFpFNcounts()+12])

#     @test metaData[3,3,5,getBeginingOfFpFNcounts()+15]==dataBdim[1]*dataBdim[2]*dataBdim[3]

#     sum(Array(metaData[3,3,5,getBeginingOfFpFNcounts():getBeginingOfFpFNcounts()+13]))

#     metaData[3,3,5,getBeginingOfFpFNcounts()+1]+ metaData[3,3,5,getBeginingOfFpFNcounts()+3]+ metaData[3,3,5,getBeginingOfFpFNcounts()+5]+ metaData[3,3,5,getBeginingOfFpFNcounts()+7]+ metaData[3,3,5,getBeginingOfFpFNcounts()+9]+ metaData[3,3,5,getBeginingOfFpFNcounts()+11]+ metaData[3,3,5,getBeginingOfFpFNcounts()+13]

# end #test set
# ########### uploadDataToMetaData 2
# @testset "uploadDataToMetaData2" begin

#     # localQuesValues = CUDA.zeros(Float32,14)
#     # threads=(32,5)
#     # blocks =8

#     # mainArrDims= (634,521,632)
#     # dataBdim = (43,21,17)

#     # metaDataDims= (cld(mainArrDims[1],dataBdim[1] ),cld(mainArrDims[2],dataBdim[2]),cld(mainArrDims[3],dataBdim[3]))
#     # metaData = MetaDataUtils.allocateMetadata(mainArrDims,dataBdim);

#     # inBlockLoopX,inBlockLoopY,inBlockLoopZ= (fld(dataBdim[1] ,threads[1]),fld(dataBdim[2] ,threads[2]),dataBdim[3]    )
#     # #we are iterating here block by block sequentially
#     # loopXMeta,loopYZMeta= (metaDataDims[1],fld(metaDataDims[2]*metaDataDims[3] ,blocks)  )
#     # yTimesZmeta= metaDataDims[2]*metaDataDims[3]
#     # metaDataLength= metaDataDims[1]*metaDataDims[2]*metaDataDims[3]
#     # loopMeta= fld(metaDataLength,blocks )
#     # function uploadDataToMetaDataKernel(loopMeta,metaDataLength,mainArrDims,metaDataDims,loopXMeta,loopYZMeta,yTimesZmeta,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaData)

#     #     localQuesValues= @cuStaticSharedMem(UInt32, 14)
#     #     localBool=false
#     #     @iter3dOuter(metaDataDims,loopMeta,metaDataLength,
#     #     begin

#     #         @iterDataBlock(mainArrDims,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ,
#     #         begin
#     #                     boolGold=false
#     #                     boolSegm=true
#     #                     localBool=true
#     #                     if(xMeta==2 && yMeta==2 && zMeta==4)

#     #                         PrepareArrtoBool.@uploadLocalfpFNCounters()
#     #                     end

#     #                 end)
#     #                 sync_threads()

#     #                 isAnyPositive= true
#     #                     PrepareArrtoBool.@uploadDataToMetaData()
#     #                     sync_threads()

#     #                     sync_threads()
#     #                     #threadfence()
#     #                     @ifY 1 if(threadIdxX()<15)
#     #                         localQuesValues[threadIdxX()]=0
#     #                     end
#     #             end)

#     #     return
#     # end

#     # @cuda threads=threads blocks=blocks uploadDataToMetaDataKernel(loopMeta,metaDataLength,mainArrDims,metaDataDims,loopXMeta,loopYZMeta,yTimesZmeta,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaData)

#     # @test metaData[3,3,5,getBeginingOfFpFNcounts()+2]==dataBdim[2]*dataBdim[3]
#     # @test metaData[3,3,5,getBeginingOfFpFNcounts()+4]==dataBdim[2]*dataBdim[3]

#     # @test metaData[3,3,5,getBeginingOfFpFNcounts()+6]==dataBdim[1]*dataBdim[3]-2*dataBdim[3]
#     # @test metaData[3,3,5,getBeginingOfFpFNcounts()+8]==dataBdim[1]*dataBdim[3]-2*dataBdim[3]

#     # @test metaData[3,3,5,getBeginingOfFpFNcounts()+10]==(dataBdim[1]-2)*(dataBdim[2]-2)
#     # @test metaData[3,3,5,getBeginingOfFpFNcounts()+12]==(dataBdim[1]-2)*(dataBdim[2]-2)

#     # @test metaData[3,3,5,getBeginingOfFpFNcounts()+14]==dataBdim[1]*dataBdim[2]*dataBdim[3] - sum(metaData[3,3,5,getBeginingOfFpFNcounts():getBeginingOfFpFNcounts()+12])

#     # @test metaData[3,3,5,getBeginingOfFpFNcounts()+16]==dataBdim[1]*dataBdim[2]*dataBdim[3]

#     # sum(Array(metaData[3,3,5,getBeginingOfFpFNcounts():getBeginingOfFpFNcounts()+13]))

# end#testset

#################   uploadMinMaxesToShmem

@testset "uploadMinMaxesToShmem" begin
    localQuesValues = CUDA.zeros(Float32,14)
    threads=(32,5)
    blocks =17

    mainArrDims= (634,521,632)
    dataBdim = (17,7,12)

    metaDataDims= (cld(mainArrDims[1],dataBdim[1] ),cld(mainArrDims[2],dataBdim[2]),cld(mainArrDims[3],dataBdim[3]))

    inBlockLoopX,inBlockLoopY,inBlockLoopZ= (fld(dataBdim[1] ,threads[1]),fld(dataBdim[2] ,threads[2]),dataBdim[3]    )
    #we are iterating here block by block sequentially
    loopXMeta,loopYZMeta= (metaDataDims[1],fld(metaDataDims[2]*metaDataDims[3] ,blocks)  )
    yTimesZmeta= metaDataDims[2]*metaDataDims[3]
    metaData = MetaDataUtils.allocateMetadata(mainArrDims,dataBdim);

    minX= CuArray([Float32(1110.0) ])
    maxX= CuArray([Float32(0.0)])
    minY= CuArray([Float32(1110.0)])
    maxY= CuArray([Float32(0.0)    ])
    minZ= CuArray([Float32(1110.0)])
    maxZ= CuArray([Float32(0.0) ])
    metaDataLength= metaDataDims[1]*metaDataDims[2]*metaDataDims[3]
    loopMeta= fld(metaDataLength,blocks )

    function uploadDataToMetaDataKernel(loopMeta,metaDataLength,minX, maxX, minY,maxY,minZ,maxZ,mainArrDims,metaDataDims,loopXMeta,loopYZMeta,yTimesZmeta,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaData)
        localQuesValues= @cuStaticSharedMem(UInt32, 14)
        isAnyPositive=true
        @iter3dOuter(metaDataDims,loopMeta,metaDataLength,
        begin

            @iterDataBlock(mainArrDims,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ,xMeta,yMeta,zMeta,
            begin

                        if(xMeta>=2 && yMeta>=1 && zMeta>=3 && xMeta<=4 && yMeta<=5 && zMeta<=6  )
                            sync_threads()
                            PrepareArrtoBool.@uploadDataToMetaData()

                            PrepareArrtoBool.@uploadMinMaxesToShmem()
                            sync_threads()

                        end

                    end)

                end)

        return
    end

    @cuda threads=threads blocks=blocks uploadDataToMetaDataKernel(loopMeta,metaDataLength,minX, maxX, minY,maxY,minZ,maxZ,mainArrDims,metaDataDims,loopXMeta,loopYZMeta,yTimesZmeta,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaData)
    @test minX[1]== 2+1
    @test maxX[1]== 4+1
    @test minY[1]== 1+1
    @test maxY[1]== 5+1
    @test minZ[1]== 3+1
    @test maxZ[1]== 6+1
end #test set

###########################   getBoolCubeKernel
@testset "getBoolCubeKernel" begin

    using Revise, Parameters, Logging, Test
    using CUDA
    includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
    using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils, ..PrepareArrtoBool, ..BitWiseUtils
    using ..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates, ..Housdorff

    threads=(32,5)
    blocks =8
    mainArrDims= (512,512,155)
    dataBdim = (2,2,2)
    metaData = MetaDataUtils.allocateMetadata(mainArrDims,dataBdim);
    metaDataDims= size(metaData)
    arrGold = zeros(Int32,mainArrDims);
    arrSegm = zeros(Int32,mainArrDims);

        # so we should get all mins as 2 and all maxes as 5

        arrGold[9,2,5]= 2
        arrGold[3,6,4]= 2
        arrGold[3,3,4]= 2

        arrGold[4,4,4]= 2
        arrSegm[4,4,4]= 2

        arrSegm[3,18,4]= 2
        arrSegm[5,2,11]= 2

        gold3d,segm3d = CuArray(arrGold), CuArray(arrSegm)
        numberToLooFor =Int32(2)
        robustnessPercent = 0.9

        boolKernelArgs, mainKernelArgs,threadsBoolKern,blocksBoolKern ,threadsMainKern,blocksMainKern ,shmemSizeBool,shmemSizeMain=    preparehousedorfKernel(gold3d,segm3d,robustnessPercent,numberToLooFor)
        mainArrDims,dataBdim,metaData,metaDataDims,reducedGoldA,reducedSegmA,loopXinPlane,loopYinPlane,minxRes,maxxRes,minyRes,maxyRes,minzRes,maxzRes,fn,fp ,numberToLooFor,inBlockLoopXZIterWithPadding,shmemblockDataLoop,shmemblockDataLenght,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength=boolKernelArgs

        function locForKernel(goldGPU,segmGPU,mainArrDims,dataBdim,metaData,metaDataDims,reducedGoldA,reducedSegmA,loopXinPlane,loopYinPlane,minxRes,maxxRes,minyRes,maxyRes,minzRes,maxzRes,fn,fp ,numberToLooFor,inBlockLoopXZIterWithPadding,shmemblockDataLoop,shmemblockDataLenght,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength)

        @getBoolCubeKernel()

        return
    end

        @cuda threads=threadsBoolKern blocks=blocksBoolKern locForKernel(gold3d,segm3d,boolKernelArgs...)
            @test  Int64( maxzRes[])==6

        @test Int64(fn[])==3
        @test Int64(fp[])==2

        @test Int64(minxRes[])==2
        @test Int64( maxxRes[])==5
        @test Int64(minyRes[])==1
        @test  Int64(maxyRes[])==9
        @test  Int64( minzRes[])==2
        @test  Int64( maxzRes[])==6

        @test  reducedGoldA[9,2,5]==true
        @test   reducedGoldA[3,6,4]==true
        @test   reducedGoldA[3,3,4]==true
        @test   reducedGoldA[4,4,4]==true
        @test   reducedGoldA[5,5,5]==false
        @test   reducedGoldA[5,2,5]==false
        @test   reducedGoldA[5,2,90]==false

        @test  reducedGoldB[9,2,5]==true
        @test   reducedGoldB[3,6,4]==true
        @test   reducedGoldB[3,3,4]==true
        @test   reducedGoldB[4,4,4]==true
        @test   reducedGoldB[5,5,5]==false
        @test   reducedGoldB[5,2,5]==false
        @test   reducedGoldB[5,2,90]==false

        @test  reducedSegmA[3,18,4]==true
        @test   reducedSegmA[5,2,11]==true
        @test   reducedSegmA[4,4,4]==true
        @test   reducedSegmA[5,5,5]==false
        @test   reducedSegmA[5,2,5]==false
        @test   reducedSegmA[5,2,90]==false

        @test  reducedSegmB[3,18,4]==true
        @test   reducedSegmB[5,2,11]==true
        @test   reducedSegmB[4,4,4]==true
        @test   reducedSegmB[5,5,5]==false
        @test   reducedSegmB[5,2,5]==false
        @test   reducedSegmB[5,2,90]==false
end#test set

###########################   getBoolCubeKernel  2

using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils, ..PrepareArrtoBool
using ..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates

 @testset "getBoolCubeKernel 2 " begin
    localQuesValues = CUDA.zeros(Float32,14)
    threads=(32,5)
    blocks =7

    mainArrDims= (512,512,400)
    dataBdim = (10,10,10)
    metaData = MetaDataUtils.allocateMetadata(mainArrDims,dataBdim);
    metaDataDims= size(metaData)
    arrGold = zeros(Int32,mainArrDims);
    arrSegm = zeros(Int32,mainArrDims);
    loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength = calculateLoopsIter(dataBdim,threads[1],threads[2],metaDataDims,blocks)
    minxRes,maxxRes,minyRes,maxyRes,minzRes,maxzRes,fn,fp  =getSmallForBoolKernel();
    reducedGoldA,reducedSegmA,reducedGoldB,reducedSegmB=  getLargeForBoolKernel(mainArrDims,dataBdim);

        # so we should get all mins as 2 and all maxes as 5

        # 1)   Left FP
        # 2)   Left FN
        arrGold[11,13,13]= 2
        arrSegm[11,13,15]= 2

    # 3)   Right FP
    # 4)   Right FN
        arrGold[20,13,13]= 2
        arrSegm[20,13,15]= 2

        arrGold[20,14,13]= 2
        arrSegm[20,14,15]= 2

    # 5)   Posterior FP

        # 6)   Posterior FN

        arrGold[13,11,14]= 2
        arrSegm[14,11,13]= 2

        arrGold[15,11,13]= 2
        arrSegm[13,11,15]= 2

        arrGold[16,11,13]= 2
        arrSegm[13,11,16]= 2

        # 7)   Anterior FP
        # 8)   Anterior FN
        arrGold[14,20,13]= 2
        arrSegm[13,20,14]= 2

        arrGold[15,20,13]= 2
        arrSegm[14,20,16]= 2

        arrGold[17,20,13]= 2
        arrSegm[15,20,17]= 2

        arrGold[18,20,17]= 2
        arrSegm[16,20,18]= 2

        # 9)   Top FP
        # 10)   Top FN
        arrGold[14,15,11]= 2
        arrSegm[14,16,11]= 2

        arrGold[17,15,11]= 2
        arrSegm[17,16,11]= 2

        arrGold[15,15,11]= 2
        arrSegm[15,16,11]= 2

        arrGold[16,15,11]= 2
        arrSegm[16,16,11]= 2

        arrGold[18,15,11]= 2
        arrSegm[18,16,11]= 2

        # 11)   Bottom FP
        # 12)   Bottom FN

        arrGold[14,15,20]= 2
        arrSegm[14,16,20]= 2

        arrGold[17,15,20]= 2
        arrSegm[17,16,20]= 2

        arrGold[15,15,20]= 2
        arrSegm[15,16,20]= 2

        arrGold[16,15,20]= 2
        arrSegm[16,16,20]= 2

        arrGold[18,15,20]= 2
        arrSegm[18,16,20]= 2

        arrGold[19,15,20]= 2
        arrSegm[19,16,20]= 2

        # 13)   main block Fp
        # 14)   main block Fn

        arrGold[14,15,14]= 2
        arrSegm[14,16,16]= 2

        arrGold[17,15,17]= 2
        arrSegm[17,16,13]= 2

        arrGold[15,15,13]= 2
        arrSegm[15,16,13]= 2

        arrGold[16,15,13]= 2
        arrSegm[16,16,13]= 2

        arrGold[18,15,13]= 2
        arrSegm[18,16,13]= 2

        arrGold[19,15,13]= 2
        arrSegm[19,16,13]= 2

        arrGold[19,12,13]= 2
        arrSegm[19,12,13]= 2

        arrGold[12,12,17]= 2
        arrSegm[12,12,16]= 2

    ########### some noise around

    arrGold[22,15,13]= 2
    arrSegm[18,22,13]= 2

    arrGold[19,22,13]= 2
    arrSegm[22,16,13]= 2

    arrGold[19,22,13]= 2
    arrSegm[19,12,22]= 2

    arrGold[12,12,17]= 2
    arrSegm[12,12,22]= 2

    arrGold[9,15,13]= 2
    arrSegm[18,9,13]= 2

    arrGold[19,22,9]= 2
    arrSegm[9,16,13]= 2

    arrGold[19,9,13]= 2
    arrSegm[19,12,9]= 2

    arrGold[9,12,17]= 2
    arrSegm[12,9,22]= 2

    gold3d,segm3d = CuArray(arrGold), CuArray(arrSegm)
    numberToLooFor =Int32(2)

    function locForKernelB(mainArrDims,dataBdim,metaData,metaDataDims,reducedGoldA,reducedSegmA,reducedGoldB,reducedSegmB,minxRes,maxxRes,minyRes,maxyRes,minzRes,maxzRes,fn,fp ,gold3d,segm3d,numberToLooFor,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength)

    @getBoolCubeKernel()

    return
    end

    @cuda threads=threads blocks=blocks locForKernelB(mainArrDims,dataBdim,metaData,metaDataDims,reducedGoldA,reducedSegmA,reducedGoldB,reducedSegmB,minxRes,maxxRes,minyRes,maxyRes,minzRes,maxzRes,fn,fp ,gold3d,segm3d,numberToLooFor,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength)

            #   @test   metaData[2,2,2,getBeginingOfFpFNcounts()+1]==2
            #   @test   metaData[2,3,2,getBeginingOfFpFNcounts()+2]==3
            #   @test   metaData[2,2,4,getBeginingOfFpFNcounts()+3]==4

            @test Int64(metaData[2,2,2,getBeginingOfFpFNcounts()+1])==1
            @test Int64(metaData[2,2,2,getBeginingOfFpFNcounts()+2])==1
            @test Int64(metaData[2,2,2,getBeginingOfFpFNcounts()+3])==2
            @test Int64(metaData[2,2,2,getBeginingOfFpFNcounts()+4])==2
            @test Int64(metaData[2,2,2,getBeginingOfFpFNcounts()+5])==3
            @test Int64(metaData[2,2,2,getBeginingOfFpFNcounts()+6])==3
            @test Int64(metaData[2,2,2,getBeginingOfFpFNcounts()+7])==4
            @test Int64(metaData[2,2,2,getBeginingOfFpFNcounts()+8])==4
            @test Int64(metaData[2,2,2,getBeginingOfFpFNcounts()+9])==5
            @test Int64(metaData[2,2,2,getBeginingOfFpFNcounts()+10])==5
            @test Int64(metaData[2,2,2,getBeginingOfFpFNcounts()+11])==6
            @test Int64(metaData[2,2,2,getBeginingOfFpFNcounts()+12])==6
            @test Int64(metaData[2,2,2,getBeginingOfFpFNcounts()+13])==7
            @test Int64(metaData[2,2,2,getBeginingOfFpFNcounts()+14])==7
            @test Int64(metaData[2,2,2,getBeginingOfFpFNcounts()+15])== 1+2+3+4+5+6+7
            @test Int64(metaData[2,2,2,getBeginingOfFpFNcounts()+16])== 1+2+3+4+5+6+7
end# test set

end#total test set

using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils
using ..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates

##### iter data block
singleVal = CUDA.zeros(Int64,14)
threads=(32,5)

blocks =17
mainArrDims= (67,177,90)
dataBdim = (43,21,17)
indices = CUDA.zeros(Bool,mainArrDims)

metaDataDims= (cld(mainArrDims[1],dataBdim[1] ),cld(mainArrDims[2],dataBdim[2]),cld(mainArrDims[3],dataBdim[3]))

inBlockLoopX,inBlockLoopY,inBlockLoopZ= (fld(dataBdim[1] ,threads[1]),fld(dataBdim[2] ,threads[2]),dataBdim[3]    )
#we are iterating here block by block sequentially
metaDataLength= metaDataDims[1]*metaDataDims[2]*metaDataDims[3]
loopMeta= fld(metaDataLength,blocks )

function iterDataBlocksKernel(loopMeta,metaDataLength,mainArrDims,singleVal,metaDataDims,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ ,indices)
    @iter3dOuter(metaDataDims,loopMeta,metaDataLength,
    begin
        @iterDataBlockZdeepest(mainArrDims,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ,xMeta,yMeta,zMeta,
        begin
       CUDA.@atomic singleVal[1]+=1
        indices[x,y,z]=true
        #CUDA.@cuprint "x $(x) y $(y) z $(z) zMeta $(zMeta)  yMeta $(yMeta) xMeta $(xMeta)  idX $(blockIdxX()) \n"

    end,begin indices[1]=true
                 end )end)

    return
end
@cuda threads=threads blocks=blocks iterDataBlocksKernel(loopMeta,metaDataLength,mainArrDims,singleVal,metaDataDims,dataBdim, inBlockLoopX,inBlockLoopY,inBlockLoopZ ,indices)
@test singleVal[1]==mainArrDims[1]*mainArrDims[2]*mainArrDims[3]
Int64(singleVal[1])
7*7*7
inddd = Array(indices);
cartt = CartesianIndices(inddd);
filtered = filter(cc-> !inddd[cc]  ,cartt)
sum(indices)


text:
using MPIReco

@testset "single- and multi-channel in-memory reconstruction" begin
  bSF = MPIFile(joinpath(datadir, "calibrations", "12.mdf"))
  b = MPIFile(joinpath(datadir, "measurements", "20211226_203916_MultiPatch", "1.mdf"))
  names = (:color, :x, :y, :z, :time)
  values = (1:1,
	    -27.375u"mm":1.25u"mm":11.375u"mm",
	    -11.375u"mm":1.25u"mm":27.375u"mm",
	    0.0u"mm":1.0u"mm":0.0u"mm",
	    0.0u"ms":0.6528u"ms":0.0u"ms")

  # standard reconstruction
  c1 = reconstruction(bSF, b; SNRThresh=5, frames=1, minFreq=80e3,
		      recChannels=1:2, iterations=1, spectralLeakageCorrection=true)
  @test axisnames(c1) == names
  @test axisvalues(c1) == values
  exportImage(joinpath(imgdir, "Reconstruction1.png"), arraydata(c1[1,:,:,1,1]))

  # fused lasso
  c2 = reconstruction(bSF, b; SNRThresh=5, frames=1, minFreq=80e3,
		      recChannels=1:2, iterations=100, solver="fusedlasso",
		      loadasreal=true, λ=[0.01,0.01])
  @test axisnames(c2) == names
  @test axisvalues(c2) == values
  exportImage(joinpath(imgdir, "Reconstruction2.png"), arraydata(c2[1,:,:,1,1]))

  # with interpolation
  c3 = reconstruction(bSF, b; SNRThresh=5, frames=1, minFreq=80e3,
		      recChannels=1:2, gridsize=[100,100,1],iterations=1)
  @test axisnames(c3) == names
  @test axisvalues(c3) == (values[1], -27.8u"mm":0.4u"mm":11.8u"mm", -11.8u"mm":0.4u"mm":27.8u"mm", values[4:5]...)
  exportImage(joinpath(imgdir, "Reconstruction3.png"), arraydata(c3[1,:,:,1,1]))

  # with fov adpation and center shift
  c4 = reconstruction(bSF, b; SNRThresh=5, frames=1, minFreq=80e3,
		      recChannels=1:2, gridsize=[50,50,1],fov=[0.04,0.04,1],
		      center=[0.0,-0.01,0], iterations=1)
  @test axisnames(c4) == names
  # TODO Tobi: does this make sense?
  @test axisvalues(c4) == (values[1], -27.6u"mm":0.8u"mm":11.6u"mm", -11.6u"mm":0.8u"mm":27.6u"mm", 499.5u"mm":1000.0u"mm":499.5u"mm", values[5])
  exportImage(joinpath(imgdir, "Reconstruction4.png"), arraydata(c4[1,:,:,1,1]))

  # multi colored reconstruction (deprecated interface)
  c5 = reconstruction([bSF,bSF], b;
		      SNRThresh=5, frames=1, minFreq=80e3,
		      recChannels=1:2, iterations=1)
  @test axisnames(c5) == names
  @test axisvalues(c5) == (1:2,values[2:end]...)
  exportImage(joinpath(imgdir, "Reconstruction5a.png"), arraydata(c5[1,:,:,1,1]))
  exportImage(joinpath(imgdir, "Reconstruction5b.png"), arraydata(c5[2,:,:,1,1]))

  # multi colored reconstruction (new interface)
  c5 = reconstruction(MultiContrastFile([bSF,bSF]), b;
          SNRThresh=5, frames=1, minFreq=80e3,
          recChannels=1:2, iterations=1)
  @test axisnames(c5) == names
  @test axisvalues(c5) == (1:2,values[2:end]...)
  exportImage(joinpath(imgdir, "Reconstruction5c.png"), arraydata(c5[1,:,:,1,1]))
  exportImage(joinpath(imgdir, "Reconstruction5d.png"), arraydata(c5[2,:,:,1,1]))

  # dict based reco
  r = defaultRecoParams()
  r[:measPath] = filepath(b)
  r[:SFPath] = filepath(bSF)
  r[:frames] =  1
  r[:minFreq] = 80e3
  r[:SNRThresh] = 4
  r[:lambd] = 0.001
  r[:iterations] = 1
  c6 = reconstruction(r)
  @test axisnames(c6) == names
  @test axisvalues(c6) == values
  exportImage(joinpath(imgdir, "Reconstruction6.png"), arraydata(c6[1,:,:,1,1]))

  # channel weighting
  c7a = reconstruction(bSF, b; SNRThresh=5, frames=1, minFreq=80e3,
		      recChannels=1:2, iterations=1, λ=0.1,
          weightType=WeightingType.Channel, channelWeights=[1.0,1.0,1.0])

  exportImage(joinpath(imgdir, "Reconstruction7a.png"), arraydata(c7a[1,:,:,1,1]))

  c7b = reconstruction(bSF, b; SNRThresh=5, frames=1, minFreq=80e3,
		      recChannels=1:2, iterations=1, λ=0.1,
          weightType=WeightingType.Channel, channelWeights=[1.0,0.001,1.0])

  exportImage(joinpath(imgdir, "Reconstruction7b.png"), arraydata(c7b[1,:,:,1,1]))

  c7c = reconstruction(bSF, b; SNRThresh=5, frames=1, minFreq=80e3,
		      recChannels=1:2, iterations=1, λ=0.1,
          weightType=WeightingType.Channel, channelWeights=[0.001,1.0,1.0])

  exportImage(joinpath(imgdir, "Reconstruction7c.png"), arraydata(c7c[1,:,:,1,1]))

end

ext: jl


text:
module CollabUtils2MapSwapTrigger

using ..Ahorn, Maple

@mapdef Trigger "CollabUtils2/MapSwapTrigger" MapSwapTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight,
    map::String="Celeste/1-ForsakenCity", side::String="Normal", room::String="1")

const placements = Ahorn.PlacementDict(
    "Map Swap (Collab Utils 2)" => Ahorn.EntityPlacement(
        MapSwapTrigger,
        "rectangle",
    ),
)

function Ahorn.editingOptions(trigger::MapSwapTrigger)
    return Dict{String, Any}(
        "side" => String["Normal", "BSide", "CSide"]
    )
end

end


text:
module TestScikitLearn

# using Revise
using MLJBase
using Test
using LinearAlgebra
import Random.seed!
seed!(1234)

import MLJModels
import ScikitLearn
using MLJModels.ScikitLearn_
using CategoricalArrays

## CLASSIFIERS
@test_logs (:warn,"kernel parameter is not valid, setting to default=\"rbf\" \n") SVMClassifier(kernel="wrong")
@test_logs (:warn,"penalty parameter is not valid, setting to default=\"l2\" \n") SVMLClassifier(penalty="wrong")
@test_logs (:warn,"loss parameter is not valid, setting to default=\"epsilon_insensitive\" \n") SVMLRegressor(loss="wrong")

plain_classifier = SVMClassifier()
nu_classifier = SVMNuClassifier()
linear_classifier = SVMLClassifier(max_iter=10000)

# test preservation of categorical levels:
task = load_iris();
X, y = X_and_y(task)
train, test = partition(eachindex(y), 0.6) # levels of y are split across split

fitresultC, cacheC, reportC = MLJBase.fit(plain_classifier, 1,
                                          selectrows(X, train), y[train]);
fitresultCnu, cacheCnu, reportCnu = MLJBase.fit(nu_classifier, 1,
                                          selectrows(X, train), y[train]);
fitresultCL, cacheCL, reportCL = MLJBase.fit(linear_classifier, 1,
                                          selectrows(X, train), y[train]);
pcpred = predict(plain_classifier, fitresultC, selectrows(X, test));
nucpred = predict(nu_classifier, fitresultCnu, selectrows(X, test));
lcpred = predict(linear_classifier, fitresultCL, selectrows(X, test));

@test Set(classes(pcpred[1])) == Set(classes(y[1]))
@test Set(classes(nucpred[1])) == Set(classes(y[1]))
@test Set(classes(lcpred[1])) == Set(classes(y[1]))

# test with linear data:
x1 = randn(3000);
x2 = randn(3000);
x3 = randn(3000);
X = (x1=x1, x2=x2, x3=x3);
y = x1 - x2 -2x3;
ycat = map(y) do η
    η > 0 ? "go" : "stop"
end |> categorical;
train, test = partition(eachindex(ycat), 0.8);
fitresultC, cacheC, reportC = MLJBase.fit(plain_classifier, 1,
                                          selectrows(X, train), ycat[train]);
fitresultCnu, cacheCnu, reportCnu = MLJBase.fit(nu_classifier, 1,
                                          selectrows(X, train), ycat[train]);
fitresultCL, cacheCL, reportCL = MLJBase.fit(linear_classifier, 1,
                                          selectrows(X, train), ycat[train]);
pcpred = predict(plain_classifier, fitresultC, selectrows(X, test));
nucpred = predict(nu_classifier, fitresultCnu, selectrows(X, test));
lcpred = predict(linear_classifier, fitresultCL, selectrows(X, test));
@test sum(pcpred .!= ycat[test])/length(ycat) < 0.05
@test sum(nucpred .!= ycat[test])/length(ycat) < 0.05
@test sum(lcpred .!= ycat[test])/length(ycat) < 0.05

## REGRESSORS

plain_regressor = SVMRegressor()
nu_regressor = SVMNuRegressor()
linear_regressor = SVMLRegressor(max_iter=10000)

# test with linear data:
fitresultR, cacheR, reportR = MLJBase.fit(plain_regressor, 1,
                                          selectrows(X, train), y[train]);
fitresultRnu, cacheRnu, reportRnu = MLJBase.fit(nu_regressor, 1,
                                                selectrows(X, train), y[train]);
fitresultRL, cacheRL, reportRL = MLJBase.fit(linear_regressor, 1,
                                             selectrows(X, train), y[train]);

rpred = predict(plain_regressor, fitresultR, selectrows(X, test));
nurpred = predict(nu_regressor, fitresultRnu, selectrows(X, test));
Lrpred = predict(linear_regressor, fitresultRL, selectrows(X, test));

@test norm(rpred - y[test])/sqrt(length(y)) < 0.2
@test norm(nurpred - y[test])/sqrt(length(y)) < 0.2
@test norm(Lrpred - y[test])/sqrt(length(y)) < 0.2

## enet = SCElasticNet()
## enetCV = SCElasticNetCV()

## ELASTIC NET

# generate some synthetic linear data:

x1 = randn(3000);
x2 = randn(3000);
x3 = randn(3000);
X = (x1=x1, x2=x2, x3=x3);
y = x1 - x2 -2x3 + 0.02*randn(3000);

train, test = partition(eachindex(y), 0.7)

# test CV version:
rgsCV = SCElasticNetCV(copy_X=true, cv=5, eps=0.001,
       fit_intercept=true, l1_ratio=[0.5, 0.9], max_iter=1000, n_alphas=100,
       normalize=false, positive=false,
       precompute="auto", selection="cyclic",
       tol=0.0001)

fitresult, cache, report = MLJBase.fit(rgsCV, 0,
                                       selectrows(X, train),
                                       y[train]);

yhat = predict(rgsCV, fitresult, selectrows(X, test));
@test norm(yhat - y[test])/sqrt(length(y)) < 0.2

# test intercept and coefficients close to true:
fitted = fitted_params(rgsCV, fitresult)
@test abs(fitted.intercept) < 0.01
@test abs(fitted.coef[3] .+ 2) < 0.01

alpha, l1_ratio = report.alpha, report.l1_ratio

# test non-CV version:
rgs = SCElasticNet(alpha=alpha, l1_ratio=l1_ratio, fit_intercept=true,
                 normalize=false, precompute=false, max_iter=1000,
                 copy_X=false, tol=0.0001, warm_start=false,
                 positive=false, selection="random")

fitresult, cache, report = MLJBase.fit(rgs, 0,
                                       selectrows(X, train),
                                       y[train]);
yhat = predict(rgs, fitresult, selectrows(X, test));
@test norm(yhat - y[test])/sqrt(length(y)) < 0.2
@test keys(report) == (:n_iter,:dual_gap)
fitted = fitted_params(rgs, fitresult)
@test abs(fitted.intercept) < 0.01
@test abs(fitted.coef[3] .+ 2) < 0.01

info(SVMClassifier)
info(SVMNuClassifier)
info(SVMLClassifier)
info(SVMRegressor)
info(SVMNuRegressor)
info(SVMLRegressor)

info(SCElasticNet)
info(SCElasticNetCV)

end
true


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct AdminInviteRequestsDeniedApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_adminInviteRequestsDeniedList(_api::AdminInviteRequestsDeniedApi; team_id=nothing, cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/admin.inviteRequests.denied.list", ["slackAuth"])
    Swagger.set_param(_ctx.query, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.query, "cursor", cursor)  # type String
    Swagger.set_param(_ctx.query, "limit", limit)  # type Int32
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

List all denied workspace invite requests.
Param: team_id::String
Param: cursor::String
Param: limit::Int32
Return: Dict{String, Any}
"""
function adminInviteRequestsDeniedList(_api::AdminInviteRequestsDeniedApi; team_id=nothing, cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminInviteRequestsDeniedList(_api; team_id=team_id, cursor=cursor, limit=limit, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminInviteRequestsDeniedList(_api::AdminInviteRequestsDeniedApi, response_stream::Channel; team_id=nothing, cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminInviteRequestsDeniedList(_api; team_id=team_id, cursor=cursor, limit=limit, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export adminInviteRequestsDeniedList


text: Dict{String,Any}("table" => Dict{String,Any}())

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
module MacroExpandJL

export macroexpand_jl

macroexpand_jl(io::IO, ex::Expr) = (show_unquoted(io, macroexpand(ex)); print(io, "\n"))

function macroexpand_jl(outfilename, infilename)
    str = readall(infilename)
    open(outfilename, "w") do file
        println(file, "# Auto-generated by MacroExpandJL")
        println(file, "# DO NOT EDIT---instead, edit the source file $infilename")
        indx = 1
        while indx < length(str)
            ex, indx = parse(str, indx)
            macroexpand_jl(file, ex)
        end
    end
end

function macroexpand_jl(infilename)
    base, ext = splitext(infilename)
    outfilename = string(base, "_expanded", ext)
    macroexpand_jl(outfilename, infilename)
end

### The code below is a lightly-edited version of what's in base/show.jl

## Abstract Syntax Tree (AST) printing ##

# Summary:
#   print(io, ex) defers to show_unquoted(io, ex)
#   show(io, ex) defers to show_unquoted(io, QuoteNode(ex))
#   show_unquoted(io, ex) does the heavy lifting
#
# AST printing should follow two rules:
#   1. parse(string(ex)) == ex
#   2. eval(parse(repr(ex))) == ex
#
# Rule 1 means that printing an expression should generate Julia code which
# could be reparsed to obtain the original expression. This code should be
# unambiguous and as readable as possible.
#
# Rule 2 means that showing an expression should generate a quoted version of
# print’s output. Parsing and then evaling this output should return the
# original expression.
#
# This is consistent with many other show methods, i.e.:
#   show(Set([1,2,3]))                # ==> "Set{Int64}([2,3,1])"
#   eval(parse("Set{Int64}([2,3,1])”) # ==> An actual set
# While this isn’t true of ALL show methods, it is of all ASTs.

typealias ExprNode Union{Expr, QuoteNode, SymbolNode, LineNumberNode,
                         LabelNode, GotoNode, TopNode}
show_unquoted(io::IO, ex)              = show_unquoted(io, ex, 0, 0)
show_unquoted(io::IO, ex, indent::Int) = show_unquoted(io, ex, indent, 0)
show_unquoted(io::IO, ex, ::Int,::Int) = show(io, ex)

## AST printing constants ##

const indent_width = 4
const quoted_syms = Set{Symbol}([:(:),:(::),:(:=),:(=),:(==),:(===),:(=>)])
const uni_ops = Set{Symbol}([:(+), :(-), :(!), :(~), :(<:), :(>:)])
const bin_ops_by_prec = [
    "= := += -= *= /= //= .//= .*= ./= \\= .\\= ^= .^= %= .%= |= &= \$= => <<= >>= >>>= ~ .+= .-=",
    "?",
    "||",
    "&&",
    "-- -->",
    "> < >= <= == === != !== .> .< .>= .<= .== .!= .= .! <: >:",
    "|> <|",
    ": ..",
    "+ - .+ .- | \$",
    "<< >> >>> .<< .>> .>>>",
    "* / ./ % .% & .* \\ .\\",
    "// .//",
    "^ .^",
    "::",
    "."
]
const bin_op_precs = Dict{Symbol,Int}(merge([Dict{Any, Any}([symbol(op)=>i for op=split(bin_ops_by_prec[i])]) for i=1:length(bin_ops_by_prec)]...))
const bin_ops = Set{Symbol}(keys(bin_op_precs))
const expr_infix_wide = Set([:(=), :(+=), :(-=), :(*=), :(/=), :(\=), :(&=),
    :(|=), :($=), :(>>>=), :(>>=), :(<<=), :(&&), :(||)])
const expr_infix = Set([:(:), :(<:), :(->), :(=>), symbol("::")])
const expr_calls  = Dict(:call =>('(',')'), :ref =>('[',']'), :curly =>('{','}'))
const expr_parens = Dict(:tuple=>('(',')'), :vcat=>('[',']'), :cell1d=>('{','}'),
                     :hcat =>('[',']'), :row =>('[',']'))

## AST decoding helpers ##

is_id_start_char(c::Char) = ccall(:jl_id_start_char, Cint, (UInt32,), c) != 0
is_id_char(c::Char) = ccall(:jl_id_char, Cint, (UInt32,), c) != 0
function isidentifier(s::AbstractString)
    i = start(s)
    done(s, i) && return false
    (c, i) = next(s, i)
    is_id_start_char(c) || return false
    while !done(s, i)
        (c, i) = next(s, i)
        is_id_char(c) || return false
    end
    return true
end
isoperator(s::ByteString) = ccall(:jl_is_operator, Cint, (Ptr{UInt8},), s) != 0
isoperator(s::AbstractString) = isoperator(bytestring(s))

is_expr(ex, head::Symbol)         = (isa(ex, Expr) && (ex.head == head))
is_expr(ex, head::Symbol, n::Int) = is_expr(ex, head) && length(ex.args) == n

is_linenumber(ex::LineNumberNode) = true
is_linenumber(ex::Expr)           = is(ex.head, :line)
is_linenumber(ex)                 = false

is_quoted(ex)            = false
is_quoted(ex::QuoteNode) = true
is_quoted(ex::Expr)      = is_expr(ex, :quote, 1)

unquoted(ex::QuoteNode)  = ex.value
unquoted(ex::Expr)       = ex.args[1]

## AST printing helpers ##

const indent_width = 4

function show_expr_type(io::IO, ty)
    if !is(ty, Any)
        if is(ty, Function)
            print(io, "::F")
        elseif is(ty, IntrinsicFunction)
            print(io, "::I")
        else
            print(io, "::$ty")
        end
    end
end

show_linenumber(io::IO, line)       = print(io," # line ",line,':')
show_linenumber(io::IO, line, file) = print(io," # ",file,", line ",line,':')

# show a block, e g if/for/etc
function show_block(io::IO, head, args::Vector, body, indent::Int)
    print(io, head, ' ')
    show_list(io, args, ", ", indent)

    ind = is(head, :module) ? indent : indent + indent_width
    local exs
    if is_expr(body, :block) || is_expr(body, :body)
        if !isempty(body.args) && isa(body.args[1], Expr) && (body.args[1]::Expr).head == :boundscheck
            # Print as @inbounds ...
            exblock = copy(body.args[2])
            pop!(exblock.args)  # remove the Expr(:boundscheck, :(Base.pop))
            if (body.args[1]::Expr).args[1] == false
                print(io, '\n', " "^(indent+indent_width), "@inbounds ")
            end
            show_block(io, "begin", exblock, indent+indent_width)
            print(io, "end", '\n', " "^indent)
            return
        end
        exs = body.args
    else
        exs = Any[body]
    end
    for ex in exs
        if !is_linenumber(ex); print(io, '\n', " "^ind); end
        show_unquoted(io, ex, ind)
    end
    print(io, '\n', " "^indent)
end
show_block(io::IO,head,block,i::Int) = show_block(io,head,Any[],block,i)
function show_block(io::IO, head, arg, block, i::Int)
    if is_expr(arg, :block)
        show_block(io, head, arg.args, block, i)
    else
        show_block(io, head, Any[arg], block, i)
    end
end

# show an indented list
function show_list(io::IO, items, sep, indent::Int, prec::Int=0)
    n = length(items)
    if n == 0; return end
    indent += indent_width
    show_unquoted(io, items[1], indent, prec)
    for item in items[2:end]
        print(io, sep)
        show_unquoted(io, item, indent, prec)
    end
end
# show an indented list inside the parens (op, cl)
function show_enclosed_list(io::IO, op, items, sep, cl, indent, prec=0)
    print(io, op); show_list(io, items, sep, indent, prec); print(io, cl)
end

## AST printing ##

show_unquoted(io::IO, sym::Symbol, ::Int, ::Int)        = print(io, sym)
show_unquoted(io::IO, ex::LineNumberNode, ::Int, ::Int) = show_linenumber(io, ex.line)
show_unquoted(io::IO, ex::LabelNode, ::Int, ::Int)      = print(io, ex.label, ": ")
show_unquoted(io::IO, ex::GotoNode, ::Int, ::Int)       = print(io, "goto ", ex.label)
show_unquoted(io::IO, ex::TopNode, ::Int, ::Int)        = print(io,"top(",ex.name,')')

function show_unquoted(io::IO, ex::SymbolNode, ::Int, ::Int)
    print(io, ex.name)
    show_expr_type(io, ex.typ)
end

show_unquoted(io::IO, ex::QuoteNode, indent::Int, prec::Int) =
    show_unquoted_quote_expr(io, ex.value, indent, prec)

function show_unquoted_quote_expr(io::IO, value, indent::Int, prec::Int)
    if isa(value, Symbol) && !(value in quoted_syms)
        s = string(value)
        if (isidentifier(s) || isoperator(s)) && s != "end"
            print(io, ":")
            print(io, value)
        else
            print(io, "symbol(\"", escape_string(s), "\")")
        end
    else
        if isa(value,Expr) && value.head === :block
            show_block(io, "quote", value, indent)
            print(io, "end")
        else
            print(io, ":(")
            show_unquoted(io, value, indent+indent_width, 0)
            print(io, ")")
        end
    end
end

# TODO: implement interpolated strings
function show_unquoted(io::IO, ex::Expr, indent::Int, prec::Int)
    head, args, nargs = ex.head, ex.args, length(ex.args)

    # dot (i.e. "x.y")
    if is(head, :(.))
        show_unquoted(io, args[1], indent + indent_width)
        print(io, '.')
        if is_quoted(args[2])
            show_unquoted(io, unquoted(args[2]), indent + indent_width)
        else
            print(io, '(')
            show_unquoted(io, args[2], indent + indent_width)
            print(io, ')')
        end

    # infix (i.e. "x<:y" or "x = y")
    elseif (head in expr_infix && nargs==2) || (is(head,:(:)) && nargs==3)
        show_list(io, args, head, indent)
    elseif head in expr_infix_wide && nargs == 2
        func_prec = get(bin_op_precs, head, 0)
        if func_prec < prec
            show_enclosed_list(io, '(', args, " $head ", ')', indent, func_prec)
        else
            show_list(io, args, " $head ", indent, func_prec)
        end

    # list (i.e. "(1,2,3)" or "[1,2,3]")
    elseif haskey(expr_parens, head)               # :tuple/:vcat/:cell1d
        op, cl = expr_parens[head]
        if head === :vcat && !isempty(args) && is_expr(args[1], :row)
            sep = ";"
        elseif head === :hcat || head === :row
            sep = " "
        else
            sep = ","
        end
        head !== :row && print(io, op)
        show_list(io, args, sep, indent)
        if is(head, :tuple) && nargs == 1; print(io, ','); end
        head !== :row && print(io, cl)

    # function call
    elseif haskey(expr_calls, head) && nargs >= 1  # :call/:ref/:curly
        func = args[1]
        func_prec = get(bin_op_precs, func, 0)
        func_args = args[2:end]

        # scalar multiplication (i.e. "100x")
        if func == :(*) && length(func_args)==2 && isa(func_args[1], Real) && isa(func_args[2], Symbol)
            if func_prec <= prec
                show_enclosed_list(io, '(', func_args, "", ')', indent, func_prec)
            else
                show_list(io, func_args, "", indent, func_prec)
            end

        # unary operator (i.e. "!z")
        elseif func in uni_ops && length(func_args) == 1
            show_unquoted(io, func, indent)
            if isa(func_args[1], Expr) || length(func_args) > 1
                show_enclosed_list(io, '(', func_args, ",", ')', indent, func_prec)
            else
                show_unquoted(io, func_args[1])
            end

        # binary operator (i.e. "x + y")
        elseif func in bin_ops
            if length(func_args) > 1
                sep = func_prec >= bin_op_precs[:(^)] ? "$func" : " $func "
                if func_prec <= prec
                    show_enclosed_list(io, '(', func_args, sep, ')', indent, func_prec)
                else
                    show_list(io, func_args, sep, indent, func_prec)
                end
            else
                # 1-argument call to normally-binary operator
                op, cl = expr_calls[head]
                print(io, "(")
                show_unquoted(io, func, indent)
                print(io, ")")
                show_enclosed_list(io, op, func_args, ",", cl, indent)
            end

        # normal function (i.e. "f(x,y)")
        else
            op, cl = expr_calls[head]
            show_unquoted(io, func, indent)
            show_enclosed_list(io, op, func_args, ",", cl, indent)
        end
    elseif is(head, :ccall)
        show_unquoted(io, :ccall, indent)
        show_enclosed_list(io, '(', args, ",", ')', indent)

    # comparison (i.e. "x < y < z")
    elseif is(head, :comparison) && nargs >= 3 && (nargs&1==1)
        comp_prec = minimum([get(bin_op_precs, comp, 0) for comp=args[2:2:end]])
        if comp_prec <= prec
            show_enclosed_list(io, '(', args, " ", ')', indent, comp_prec)
        else
            show_list(io, args, " ", indent, comp_prec)
        end

    # block with argument
    elseif head in (:for,:while,:function,:if) && nargs==2
        show_block(io, head, args[1], args[2], indent); print(io, "end")

    elseif head == :module
        modulesym = args[2]
        args[1] ? println(io, "module ", string(args[2])) : println(io, "baremodule ", string(modulesym))
        modulebody = args[3]
        if length(modulebody.args) > 2
            # Skip the initial `eval(x) = top(Core).eval($modulesym,x)` and the scoped version
            ex1 = modulebody.args[1]
            if isa(ex1, Expr) && (ex1::Expr).head == :(=)
                ex2 = (ex1::Expr).args[1]
                if isa(ex2, Expr) && (ex2::Expr).head == :call && (ex2::Expr).args[1] == :eval
                    # Starts with `eval(x) = ...`, good enough
                    modulebody = copy(modulebody)
                    deleteat!(modulebody.args, (1,2))
                end
            end
        end
        show_unquoted(io, modulebody)
        println(io, "end")

    # type declaration
    elseif is(head, :type) && nargs==3
        show_block(io, args[1] ? :type : :immutable, args[2], args[3], indent)
        print(io, "end")

    # empty return (i.e. "function f() return end")
    elseif is(head, :return) && nargs == 1 && is(args[1], nothing)
        print(io, head)

    # type annotation (i.e. "::Int")
    elseif is(head, symbol("::")) && nargs == 1
        print(io, "::")
        show_unquoted(io, args[1], indent)

    # var-arg declaration or expansion
    # (i.e. "function f(L...) end" or "f(B...)")
    elseif is(head, :(...)) && nargs == 1
        show_unquoted(io, args[1], indent)
        print(io, "...")

    elseif (nargs == 1 && head in (:return, :abstract, :const)) ||
                          head in (:local,  :global)
        print(io, head, ' ')
        show_list(io, args, ", ", indent)
    elseif is(head, :macrocall) && nargs >= 1
        show_list(io, args, ' ', indent)
    elseif is(head, :typealias) && nargs == 2
        print(io, "typealias ")
        show_list(io, args, ' ', indent)
    elseif is(head, :line) && 1 <= nargs <= 2
        show_linenumber(io, args...)
    elseif is(head, :if) && nargs == 3           # if/else
        show_block(io, "if",   args[1], args[2], indent)
        show_block(io, "else", args[3], indent)
        print(io, "end")
    elseif is(head, :try) && 3 <= nargs <= 4
        show_block(io, "try", args[1], indent)
        if is_expr(args[3], :block)
            show_block(io, "catch", is(args[2], false) ? [] : args[2], args[3], indent)
        end
        if nargs >= 4 && is_expr(args[4], :block)
            show_block(io, "finally", [], args[4], indent)
        end
        print(io, "end")
    elseif is(head, :let) && nargs >= 1
        show_block(io, "let", args[2:end], args[1], indent); print(io, "end")
    elseif is(head, :block) || is(head, :body)
        show_block(io, "begin", ex, indent); print(io, "end")
    elseif is(head, :quote) && nargs == 1
        show_unquoted_quote_expr(io, args[1], indent, 0)
    elseif is(head, :gotoifnot) && nargs == 2
        print(io, "unless ")
        show_list(io, args, " goto ", indent)
    elseif is(head, :string) && nargs == 1 && isa(args[1], AbstractString)
        show(io, args[1])
    elseif is(head, :null)
        print(io, "nothing")
    elseif is(head, :kw) && length(args)==2
        show_unquoted(io, args[1], indent+indent_width)
        print(io, '=')
        show_unquoted(io, args[2], indent+indent_width)
    elseif is(head, :string)
        a = map(args) do x
            if !isa(x,AbstractString)
                if isa(x,Symbol) && !(x in quoted_syms)
                    string("\$", x)
                else
                    string("\$(", sprint(show_unquoted,x), ")")
                end
            else
                sprint(print_escaped, x, "\"\$")
            end
        end
        print(io, '"', a..., '"')

    elseif is(head, :&) && length(args) == 1
        print(io, '&')
        show_unquoted(io, args[1])

    # print anything else as "Expr(head, args...)"
    else
        print(io, "\$(Expr(")
        show(io, ex.head)
        for arg in args
            print(io, ", ")
            if isa(arg, ExprNode)
                show_unquoted(io, arg)
            else
                show(io, arg)
            end
        end
        print(io, "))")
    end

    show_expr_type(io, ex.typ)
end
end


text:
module ContractionSequenceOptimization

  using NDTensors

  import NDTensors:
    dim

  export optimized_contraction_sequence

  include("utils.jl")
  include("three_tensors.jl")
  include("breadth_first_constructive.jl")

  # The integer type of the dimensions and costs.
  # Needs to be large to avoid overflow.
  const DimT = UInt128

  """
      optimized_contraction_sequence(T)

  Returns a contraction sequence for contracting the tensors `T`. The sequence is generally optimal (currently, outer product contractions are skipped, but some optimal sequences require outer product contractions).
  """
  function optimized_contraction_sequence(T)
    if length(T) == 1
      return Any[1]
    elseif length(T) == 2
      return Any[1, 2]
    elseif length(T) == 3
      return optimized_contraction_sequence(T[1], T[2], T[3])
    end
    return breadth_first_constructive(T)
  end

end # module ContractionSequenceOptimization


text:
@inline norm2(x) = norm(x)
@inline norm2(x::Dual) = hypot(x.value, x.epsilon)

@inline errcheck(x, y, tol) = isfinite(x) && isfinite(y) && (norm2(x-y) > max(norm2(x), norm2(y))*tol)

# Same as in FastTransforms.jl
"""
Pochhammer symbol ``(x)_n = \\frac{\\Gamma(x+n)}{\\Gamma(x)}`` for the rising factorial.
"""
function pochhammer(x::T, n::Integer) where T
    S = typeof(zero(T)/one(T))
    ret = one(S)
    if n ≥ 0
        for i = 0:n-1
            ret *= x+i
        end
    else
        ret /= pochhammer(x+n, -n)
    end
    ret
end

pochhammer(x::T, n::Number) where T = isinteger(n) ? pochhammer(x, Int(n)) : ogamma(x)/ogamma(x+n)

function pochhammer(x::Number, n::UnitRange{T}) where T <: Real
    ret = Vector{promote_type(typeof(x),T)}(undef, length(n))
    ret[1] = pochhammer(x, first(n))
    for i = 2:length(n)
        ret[i] = (x+n[i]-1)*ret[i-1]
    end
    ret
end

ogamma(x::Number) = (isinteger(x) && x<0) ? zero(float(x)) : inv(gamma(x))

"""
    @clenshaw(x, c...)

Compute `∑_{i=1}^N cᵢ Tᵢ₋₁(x)`.
"""
macro clenshaw(x, c...)
    a, b = :(zero(t)), :(zero(t))
    as = []
    N = length(c)
    for k = N:-1:2
        ak = Symbol("a", k)
        push!(as, :($ak = $a))
        a = :(muladd(t, $a, $(esc(c[k]))-$b))
        b = :($ak)
    end
    ex = Expr(:block, as..., :(muladd(t/2, $a, $(esc(c[1]))-$b)))
    Expr(:block, :(t = $(esc(2))*$(esc(x))), ex)
end

const ρ = 0.72
const ρϵ = 0.71

struct ℕ end

Base.in(n::Integer, ::Type{ℕ}) = n > 0
Base.in(n::Real, ::Type{ℕ}) = (ν = round(Int, n); n == ν && ν ∈ ℕ)
Base.in(n::Complex, ::Type{ℕ}) = imag(n) == 0 && real(n) ∈ ℕ
Base.in(n::Dual, ::Type{ℕ}) = dualpart(n) == 0 && realpart(n) ∈ ℕ

struct ℕ₀ end

Base.in(n::Integer, ::Type{ℕ₀}) = n ≥ 0
Base.in(n::Real, ::Type{ℕ₀}) = (ν = round(Int, n); n == ν && ν ∈ ℕ₀)
Base.in(n::Complex, ::Type{ℕ₀}) = imag(n) == 0 && real(n) ∈ ℕ₀
Base.in(n::Dual, ::Type{ℕ₀}) = dualpart(n) == 0 && realpart(n) ∈ ℕ₀

struct ℤ end

Base.in(n::Integer, ::Type{ℤ}) = true
Base.in(n::Real, ::Type{ℤ}) = n == round(Int, n)
Base.in(n::Complex, ::Type{ℤ}) = imag(n) == 0 && real(n) ∈ ℤ
Base.in(n::Dual, ::Type{ℤ}) = dualpart(n) == 0 && realpart(n) ∈ ℤ

abeqcd(a, b, cd) = isequal(a, b) && isequal(b, cd)
abeqcd(a, b, c, d) = isequal(a, c) && isequal(b, d)

absarg(z) = abs(angle(z))

sqrtatanhsqrt(x) = x == 0 ? one(x) : (s = sqrt(-x); atan(s)/s)
sqrtasinsqrt(x) = x == 0 ? one(x) : (s = sqrt(x); asin(s)/s)
sinnasinsqrt(n, x) = x == 0 ? one(x) : (s = sqrt(x); sin(n*asin(s))/(n*s))
cosnasinsqrt(n, x) = cos(n*asin(sqrt(x)))
expnlog1pcoshatanhsqrt(n, x) = x == 0 ? one(x) : (s = sqrt(x); (exp(n*log1p(s))+exp(n*log1p(-s)))/2)
expnlog1psinhatanhsqrt(n, x) = x == 0 ? one(x) : (s = sqrt(x); (exp(n*log1p(s))-exp(n*log1p(-s)))/(2n*s))

sqrtatanhsqrt(x::Union{T, Dual{T}}) where {T<:Real} = x == 0 ? one(x) : x > 0 ? (s = sqrt(x); atanh(s)/s) : (s = sqrt(-x); atan(s)/s)
sqrtasinsqrt(x::Union{T, Dual{T}}) where {T<:Real} = x == 0 ? one(x) : x > 0 ? (s = sqrt(x); asin(s)/s) : (s = sqrt(-x); asinh(s)/s)
sinnasinsqrt(n, x::Union{T, Dual{T}}) where {T<:Real} = x == 0 ? one(x) : x > 0 ? (s = sqrt(x); sin(n*asin(s))/(n*s)) : (s = sqrt(-x); sinh(n*asinh(s))/(n*s))
cosnasinsqrt(n, x::Union{T, Dual{T}}) where {T<:Real} = x > 0 ? cos(n*asin(sqrt(x))) : cosh(n*asinh(sqrt(-x)))
expnlog1pcoshatanhsqrt(n, x::Union{T, Dual{T}}) where {T<:Real} = x == 0 ? one(x) : x > 0 ? exp(n/2*log1p(-x))*cosh(n*atanh(sqrt(x))) : exp(n/2*log1p(-x))*cos(n*atan(sqrt(-x)))
expnlog1psinhatanhsqrt(n, x::Union{T, Dual{T}}) where {T<:Real} = x == 0 ? one(x) : x > 0 ? (s = sqrt(x); exp(n/2*log1p(-x))*sinh(n*atanh(s))/(n*s)) : (s = sqrt(-x); exp(n/2*log1p(-x))*sin(n*atan(s))/(n*s))

expm1nlog1p(n, x) = x == 0 ? one(x) : expm1(n*log1p(x))/(n*x)

log1pover(s) = log1p(s)/s
logandpoly(x) = x == 0 ? one(x) : 6*(-2x+(x-2)*log1p(-x))/x^3
function logandpoly(x::Union{Float64, ComplexF64})
    if abs(x) > 0.2
        6*(-2x+(x-2)*log1p(-x))/x^3
    else
        logandpolyseries(x)
    end
end

logandpolyseries(x::Union{Float64, Dual128, ComplexF64, DualComplex256}) = @evalpoly(x, 1.0, 1.0, 0.9, 0.8, 0.7142857142857143, 0.6428571428571429, 0.5833333333333334, 0.5333333333333333, 0.4909090909090909, 0.45454545454545453, 0.4230769230769231, 0.3956043956043956, 0.37142857142857144, 0.35, 0.33088235294117646, 0.3137254901960784, 0.2982456140350877, 0.28421052631578947, 0.2714285714285714, 0.2597402597402597)

speciallog(x::Real) = iszero(x) ? one(x) : (x > 0 ? (s = sqrt(x); 3(atanh(s)-s)/s^3) : (s = sqrt(-x); 3(s-atan(s))/s^3))
speciallog(x) = iszero(x) ? one(x) : (s = sqrt(-x); 3(s-atan(s))/s^3)
function speciallog(x::Float64)
    if x > 0.2
        s = sqrt(x)
        3(atanh(s)-s)/s^3
    elseif x < -0.2
        s = sqrt(-x)
        3(s-atan(s))/s^3
    else
        speciallogseries(x)
    end
end
function speciallog(x::ComplexF64)
    if abs(x) > 0.2
        s = sqrt(-x)
        3(s-atan(s))/s^3
    else
        speciallogseries(x)
    end
end
# The Maclaurin series fails to be accurate to 1e-15 near x ≈ ±0.2. So we use a highly accurate Chebyshev expansion.
speciallogseries(x::Union{Float64, Dual128}) = @clenshaw(5.0x, 1.0087391788544393911192, 1.220474262857857637288e-01, 8.7957928919918696061703e-03, 6.9050958578444820505037e-04, 5.7037120050065804396306e-05, 4.8731405131379353370205e-06, 4.2648797509486828820613e-07, 3.800372208946157617901e-08, 3.434168059359993493634e-09, 3.1381484326392473547608e-10, 2.8939845618385022798906e-11, 2.6892186934806386106143e-12, 2.5150879096374730760324e-13, 2.3652490233687788117887e-14, 2.2349973917002118259929e-15, 2.120769988408948118084e-16)
speciallogseries(x::Union{ComplexF64, DualComplex256}) = @evalpoly(x, 1.0000000000000000000000, 5.9999999999999999999966e-01, 4.2857142857142857142869e-01, 3.3333333333333333333347e-01, 2.7272727272727272727292e-01, 2.3076923076923076923072e-01, 1.9999999999999999999996e-01, 1.7647058823529411764702e-01, 1.5789473684210526315786e-01, 1.4285714285714285714283e-01, 1.3043478260869565217384e-01, 1.2000000000000000000000e-01, 1.1111111111111111111109e-01, 1.0344827586206896551722e-01, 9.6774193548387096774217e-02, 9.0909090909090909090938e-02, 8.5714285714285714285696e-02, 8.1081081081081081081064e-02, 7.6923076923076923076907e-02, 7.3170731707317073170688e-02)

tanpi(z) = sinpi(z)/cospi(z)

const libm = Base.libm_name

unsafe_gamma(x::Float64) = ccall((:tgamma, libm),  Float64, (Float64, ), x)
unsafe_gamma(x::Float32) = ccall((:tgammaf, libm),  Float32, (Float32, ), x)
function unsafe_gamma(x::BigFloat)
    z = BigFloat()
    ccall((:mpfr_gamma, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, Base.MPFR.ROUNDING_MODE[])
    return z
end
unsafe_gamma(z::Dual) = (r = realpart(z);w = unsafe_gamma(r); dual(w, w*digamma(r)*dualpart(z)))
unsafe_gamma(z) = gamma(z)
"""
    @lanczosratio(z, ϵ, c₀, c...)

Compute `∑_{i=1}^N cᵢ/(z-1+i)/(z-1+i+ϵ) / ( c₀ + ∑_{i=1}^N cᵢ/(z-1+i) )`.
"""
macro lanczosratio(z, ϵ, c₀, c...)
    ex_num = :(zero(zm1))
    ex_den = esc(c₀)
    for i = 0:length(c)-1
        temp = :(inv(zm1+$i))
        ex_num = :(muladd($(esc(c[i+1])), $temp*inv(zm1pϵ+$i), $ex_num))
        ex_den = :(muladd($(esc(c[i+1])), $temp, $ex_den))
    end
    ex = :($ex_num/$ex_den)
    Expr(:block, :(zm1 = $(esc(z))), :(zm1pϵ = $(esc(z))+$(esc(ϵ))), ex)
end

lanczosratio(z::Union{Float64, ComplexF64, Dual128, DualComplex256}, ϵ::Union{Float64, ComplexF64, Dual128, DualComplex256}) = @lanczosratio(z, ϵ, 0.99999999999999709182, 57.156235665862923517, -59.597960355475491248, 14.136097974741747174, -0.49191381609762019978, 0.33994649984811888699E-4, 0.46523628927048575665E-4, -0.98374475304879564677E-4, 0.15808870322491248884E-3, -0.21026444172410488319E-3, 0.21743961811521264320E-3, -0.16431810653676389022E-3, 0.84418223983852743293E-4, -0.26190838401581408670E-4, 0.36899182659531622704E-5)

function lanczosapprox(z::Union{Float64, ComplexF64, Dual128, DualComplex256}, ϵ::Union{Float64, ComplexF64, Dual128, DualComplex256})
    zm0p5 = z-0.5
    zpgm0p5 = zm0p5+4.7421875
    return zm0p5*log1p(ϵ/zpgm0p5) + ϵ*log(zpgm0p5+ϵ) - ϵ + log1p(-ϵ*lanczosratio(z, ϵ))
end

function H(z::Union{Float64, ComplexF64, Dual128, DualComplex256}, ϵ::Union{Float64, ComplexF64, Dual128, DualComplex256})
    zm0p5 = z-0.5
    zpgm0p5 = zm0p5+4.7421875
    if real(z) ≥ 1/2
        if z == z+ϵ # ϵ is numerical 0
            zm0p5/zpgm0p5 + log(zpgm0p5) - 1 - lanczosratio(z, ϵ)
        else
            expm1( zm0p5*log1p(ϵ/zpgm0p5) + ϵ*log(zpgm0p5+ϵ) - ϵ + log1p(-ϵ*lanczosratio(z, ϵ)) )/ϵ
        end
    else
        tpz = tanpi(z)
        if z == z+ϵ # ϵ is numerical 0
            H(1-z, ϵ) - π/tpz
        else
            temp = (cospi(ϵ) + sinpi(ϵ)/tpz)*H(1-z, -ϵ) + .5ϵ*(π*sinc(.5ϵ))^2 - π*sinc(ϵ)/tpz
            temp/(1-ϵ*temp)
        end
    end
end

"""
Compute the function (1/Γ(z)-1/Γ(z+ϵ))/ϵ
"""
function G(z::Union{Float64, ComplexF64, Dual128, DualComplex256}, ϵ::Union{Float64, ComplexF64, Dual128, DualComplex256})
    n, zpϵ = round(Int, real(z)), z+ϵ
    if abs(ϵ) > 0.1
        (inv(unsafe_gamma(z))-inv(unsafe_gamma(zpϵ)))/ϵ
    elseif z ≠ zpϵ
        m = round(Int, real(zpϵ))
        if z == n && n ≤ 0
            -inv(ϵ*unsafe_gamma(zpϵ))
        elseif zpϵ == m && m ≤ 0
            inv(ϵ*unsafe_gamma(z))
        elseif abs(z+abs(n)) < abs(zpϵ+abs(m))
            H(z, ϵ)/unsafe_gamma(zpϵ)
        else
            H(zpϵ, -ϵ)/unsafe_gamma(z)
        end
    else # ϵ is numerical 0
        if z == n && n ≤ 0
            (-1)^(n+1)*unsafe_gamma(1-n)
        else
            digamma(z)/unsafe_gamma(z)
        end
    end
end

G(z::T, ϵ::T) where {T<:Number} = ϵ == 0 ? digamma(z)/unsafe_gamma(z) : (inv(unsafe_gamma(z))-inv(unsafe_gamma(z+ϵ)))/ϵ
G(z::Number, ϵ::Number) = G(promote(z, ϵ)...)

"""
Compute the function ((z+ϵ)ₘ-(z)ₘ)/ϵ
"""
function P(z::Number, ϵ::Number, m::Int)
    n₀ = -round(Int, real(z))
    if ϵ == 0
        if 0 ≤ n₀ < m
            ret1, ret2, n = one(z), zero(z), 0
            while n < m
                n == n₀ && (n+=1; continue)
                ret1 *= z+n
                ret2 += inv(z+n)
                n += 1
            end
            ret1 + pochhammer(z, m)*ret2
        else
            ret = zero(z)
            for n=0:m-1
                ret += inv(z+n)
            end
            pochhammer(z, m)*ret
        end
    else
        if 0 ≤ n₀ < m
            zpϵ, ret1, ret2, n = z+ϵ, one(z), zero(z), 0
            while n < m
                n == n₀ && (n+=1; continue)
                ret1 *= zpϵ+n
                ret2 += log1p(ϵ/(z+n))
                n += 1
            end
            ret1 + pochhammer(z, m)*expm1(ret2)/ϵ
        else
            ret = zero(z)
            for n=0:m-1
                ret += log1p(ϵ/(z+n))
            end
            pochhammer(z, m)*expm1(ret)/ϵ
        end
    end
end

E(z::Number, ϵ::Number) = ϵ == 0 ? z : expm1(ϵ*z)/ϵ

G(z::AbstractVector{BigFloat}, ϵ::BigFloat) = BigFloat[G(zi, ϵ) for zi in z]

# Transformation formula w = 1-z

function reconeα₀(a, b, c, m::Int, ϵ)
    _a, _b, _c, _ϵ = promote(a, b, c, ϵ)
    return _reconeα₀(_a, _b, _c, m, _ϵ)
end
function _reconeα₀(a::T, b::T, c::T, m::Int, ϵ::T) where {T}
    if ϵ == 0
        return (-1)^m*gamma(real(T)(m))*gamma(c)/(gamma(a+m)*gamma(b+m))
    else
        return gamma(c)/(ϵ*gamma(1-m-ϵ)*gamma(a+m+ϵ)*gamma(b+m+ϵ))
    end
end
function reconeβ₀(a, b, c, w, m::Int, ϵ)
    _a, _b, _c, _, _ϵ = promote(a, b, c, real(w), ϵ)
    _w, _ = promote(w, zero(_a))
    return _reconeβ₀(_a, _b, _c, _w, m, _ϵ)
end
function _reconeβ₀(a::T, b::T, c::T, w::Number, m::Int, ϵ::T) where {T}
    if abs(ϵ) > 0.1
        return ( pochhammer(a, m)*pochhammer(b, m)/(gamma(1-ϵ)*gamma(a+m+ϵ)*gamma(b+m+ϵ)*gamma(real(T)(m)+1)) - w^ϵ/(gamma(a)*gamma(b)*gamma(m+1+ϵ)) )*gamma(c)*w^m/ϵ
    else
        return ( (G(1, -ϵ)/gamma(real(T)(m)+1)+G(m+1, ϵ))/(gamma(a+m+ϵ)*gamma(b+m+ϵ)) - (G(a+m, ϵ)/gamma(b+m+ϵ)+G(float(b)+m, ϵ)/gamma(a+m))/gamma(m+1+ϵ) - E(log(w), ϵ)/(gamma(a+m)*gamma(b+m)*gamma(m+1+ϵ)) )*gamma(c)*pochhammer(a, m)*pochhammer(b, m)*w^m
    end
end
function reconeγ₀(a, b, c, w, m::Int, ϵ)
    _a, _b, _c, _, _ϵ = promote(a, b, c, real(w), ϵ)
    _w, _ = promote(w, zero(_a))
    return _reconeγ₀(_a, _b, _c, _w, m, _ϵ)
end
_reconeγ₀(a::T, b::T, c::T, w::Number, m::Int, ϵ::T) where {T} = gamma(c)*pochhammer(a, m)*pochhammer(b, m)*w^m/(gamma(a+m+ϵ)*gamma(b+m+ϵ)*gamma(real(T)(m)+1)*gamma(1-ϵ))

function Aone(a, b, c, w, m::Int, ϵ)
    αₙ = reconeα₀(a, b, c, m, ϵ)*one(w)
    ret = m ≤ 0 ? zero(w) : αₙ
    for n = 0:m-2
        αₙ *= (a+n)*(b+n)/((n+1)*(1-m-ϵ+n))*w
        ret += αₙ
    end
    ret
end

function Bone(a, b, c, w, m::Int, ϵ)
    βₙ, γₙ = reconeβ₀(a, b, c, w, m, ϵ)*one(w), reconeγ₀(a, b, c, w, m, ϵ)*w
    ret, n = βₙ, 0
    while abs(βₙ) > 10abs(ret)*eps() || n ≤ 0
        βₙ = (a+m+n+ϵ)*(b+m+n+ϵ)/((m+n+1+ϵ)*(n+1))*w*βₙ + ( (a+m+n)*(b+m+n)/(m+n+1) - (a+m+n) - (b+m+n) - ϵ + (a+m+n+ϵ)*(b+m+n+ϵ)/(n+1) )*γₙ/((m+n+1+ϵ)*(n+1-ϵ))
        ret += βₙ
        γₙ *= (a+m+n)*(b+m+n)/((m+n+1)*(n+1-ϵ))*w
        n += 1
    end
    ret
end

function _₂F₁one(a, b, c, z)
    m = round(Int, real(c-(a+b)))
    ϵ = c-(a+b)-m
    w = 1-z
    (-1)^m/sinc(ϵ)*(Aone(a, b, c, w, m, ϵ) + Bone(a, b, c, w, m, ϵ))
end

# Transformation formula w = 1/z

function recInfα₀(a, b, c, m::Int, ϵ)
    _a, _b, _c, _ϵ = promote(a, b, c, ϵ)
    return _recInfα₀(_a, _b, _c, m, _ϵ)
end
function _recInfα₀(a::T, b::T, c::T, m::Int, ϵ::T) where {T}
    if ϵ == 0
        return (-1)^m*gamma(real(T)(m))*gamma(c)/(gamma(a+m)*gamma(c-a))
    else
        return gamma(c)/(ϵ*gamma(1-m-ϵ)*gamma(a+m+ϵ)*gamma(c-a))
    end
end
function recInfβ₀(a, b, c, w, m::Int, ϵ)
    _a, _b, _c, _, _ϵ = promote(a, b, c, real(w), ϵ)
    _w, _ = promote(w, zero(_a))
    return _recInfβ₀(_a, _b, _c, _w, m, _ϵ)
end
function _recInfβ₀(a::T, b::T, c::T, w::Number, m::Int, ϵ::T) where {T}
    if abs(ϵ) > 0.1
        return ( pochhammer(a, m)*pochhammer(1-c+a, m)/(gamma(1-ϵ)*gamma(a+m+ϵ)*gamma(c-a)*gamma(real(T)(m)+1)) -
            (-w)^ϵ*pochhammer(1-c+a+ϵ, m)/(gamma(a)*gamma(c-a-ϵ)*gamma(m+1+ϵ)) )*gamma(c)*w^m/ϵ
    else
        return ( (pochhammer(1-c+a+ϵ, m)*G(1, -ϵ)-P(1-c+a, ϵ, m)/gamma(1-ϵ))/(gamma(c-a)*gamma(a+m+ϵ)*gamma(real(T)(m)+1)) +
            pochhammer(1-c+a+ϵ, m)*( (G(m+1, ϵ)/gamma(a+m+ϵ) - G(a+m, ϵ)/gamma(m+1+ϵ))/gamma(c-a) -
            (G(c-a, -ϵ) - E(-log(-w), -ϵ)/gamma(c-a-ϵ))/(gamma(m+1+ϵ)*gamma(a+m)) ) )*gamma(c)*pochhammer(a, m)*w^m
    end
end
function recInfγ₀(a, b, c, w, m::Int, ϵ)
    _a, _b, _c, _, _ϵ = promote(a, b, c, real(w), ϵ)
    _w, _ = promote(w, zero(_a))
    return _recInfγ₀(_a, _b, _c, _w, m, _ϵ)
end
_recInfγ₀(a::T, b::T, c::T, w::Number, m::Int, ϵ::T) where {T} = gamma(c)*pochhammer(a, m)*pochhammer(1-c+a, m)*w^m/(gamma(a+m+ϵ)*gamma(c-a)*gamma(real(T)(m)+1)*gamma(1-ϵ))

function AInf(a, b, c, w, m::Int, ϵ)
    αₙ = recInfα₀(a, b, c, m, ϵ)*one(w)
    ret = m ≤ 0 ? zero(w) : αₙ
    for n = 0:m-2
        αₙ *= (a+n)*(1-c+a+n)/((n+1)*(1-m-ϵ+n))*w
        ret += αₙ
    end
    ret
end

function BInf(a, b, c, win, m::Int, ϵ)
    w = win
    βₙ, γₙ = recInfβ₀(a, b, c, win, m, ϵ)*one(w), recInfγ₀(a, b, c, win, m, ϵ)*w
    ret, n = βₙ, 0
    while abs(βₙ) > 10abs(ret)*eps() || n ≤ 0
        βₙ = (a+m+n+ϵ)*(1-c+a+m+n+ϵ)/((m+n+1+ϵ)*(n+1))*w*βₙ + ( (a+m+n)*(1-c+a+m+n)/(m+n+1) - (a+m+n) - (1-c+a+m+n) - ϵ + (a+m+n+ϵ)*(1-c+a+m+n+ϵ)/(n+1) )*γₙ/((m+n+1+ϵ)*(n+1-ϵ))
        ret += βₙ
        γₙ *= (a+m+n)*(1-c+a+m+n)/((m+n+1)*(n+1-ϵ))*w
        n += 1
    end
    ret
end

function _₂F₁Inf(a, b, c, z)
    m = round(Int, real(b-a))
    ϵ = b-a-m
    w = inv(z)
    (-1)^m*(-w)^a/sinc(ϵ)*(AInf(a, b, c, w, m, ϵ) + BInf(a, b, c, w, m, ϵ))
end

function _₂F₁maclaurin(a::Number, b::Number, c::Number, z::Number)
    T = float(promote_type(typeof(a), typeof(b), typeof(c), typeof(z)))
    S₀, S₁, j = one(T), one(T)+a*b*z/c, 1
    while errcheck(S₀, S₁, 10eps(real(T))) || j ≤ 1
        rⱼ = (a+j)/(j+1)*(b+j)/(c+j)
        S₀, S₁ = S₁, S₁+(S₁-S₀)*rⱼ*z
        j += 1
    end
    return S₁
end

function _₂F₁maclaurinalt(a::Number, b::Number, c::Number, z::Number)
    T = float(promote_type(typeof(a), typeof(b), typeof(c), typeof(z)))
    C, S, j = one(T), one(T), 0
    while abs(C) > 10abs(S)*eps(real(T)) || j ≤ 1
        C *= (a+j)/(j+1)*(b+j)/(c+j)*z
        S += C
        j += 1
    end
    return S
end

function _₂F₁continuation(s::Number, t::Number, c::Number, z₀::Number, z::Number)
    T = float(promote_type(typeof(s), typeof(t), typeof(c), typeof(z₀), typeof(z)))
    izz₀, d0, d1 = inv(z-z₀), one(T), s/(2s-t+one(T))*((s+1)*(1-2z₀)+(t+1)*z₀-c)
    S₀, S₁, izz₀j, j = one(T), one(T)+d1*izz₀, izz₀, 2
    while errcheck(S₀, S₁, 10eps(real(T))) || j ≤ 2
        d0, d1, izz₀j = d1, (j+s-one(T))/j/(j+2s-t)*(((j+s)*(1-2z₀)+(t+1)*z₀-c)*d1 + z₀*(1-z₀)*(j+s-2)*d0), izz₀j*izz₀
        S₀, S₁ = S₁, S₁+d1*izz₀j
        j += 1
    end
    return S₁
end

function _₂F₁continuationalt(a::Number, c::Number, z₀::Number, z::Number)
    T = float(promote_type(typeof(a), typeof(c), typeof(z₀), typeof(z)))
    izz₀ = inv(z-z₀)
    e0, e1 = one(T), (a+one(T))*(one(T)-2z₀)+(2a+one(T))*z₀-c
    f0, f1 = zero(T), one(T)-2z₀
    cⱼ = log(z₀-z)+2digamma(one(T))-digamma(a)-digamma(c-a)
    S₀ = cⱼ
    cⱼ += 2/one(T)-one(T)/a
    C = a*izz₀
    S₁, j = S₀+(e1*cⱼ-f1)*C, 2
    while errcheck(S₀, S₁, 10eps(real(T))) || j ≤ 2
        f0, f1 = f1, (((j+a)*(1-2z₀)+(2a+1)*z₀-c)*f1+z₀*(1-z₀)*(j-1)*f0+(1-2z₀)*e1+2z₀*(1-z₀)*e0)/j
        e0, e1 = e1, (((j+a)*(1-2z₀)+(2a+1)*z₀-c)*e1+z₀*(1-z₀)*(j-1)*e0)/j
        C *= (a+j-1)*izz₀/j
        cⱼ += 2/T(j)-one(T)/(a+j-one(T))
        S₀, S₁ = S₁, S₁+(e1*cⱼ-f1)*C
        j += 1
    end
    return S₁
end

function _₂F₁logsum(a::Number, b::Number, z::Number, w::Number, s::Int)
    T = float(promote_type(typeof(a), typeof(b), typeof(z), typeof(w)))
    cⱼ = 2digamma(one(T))-digamma(a)-digamma(b)+s*log1p(-z)
    C, S, j = one(T), cⱼ, 0
    while abs(C) > 10abs(S)*eps(real(T)) || j ≤ 1
        C *= (a+j)/(j+1)^2*(b+j)*w
        cⱼ += 2/(j+one(T))-one(T)/(a+j)-one(T)/(b+j)
        S += C*cⱼ
        j += 1
    end
    return S
end

function _₂F₁logsumalt(a::Number, b::Number, z::Number, w::Number)
    T = float(promote_type(typeof(a), typeof(b), typeof(z), typeof(w)))
    d, cⱼ = one(T)-b, 2digamma(one(T))-digamma(a)-digamma(b)-log(-w)
    C, S, j = one(T), cⱼ, 0
    while abs(C) > 10abs(S)*eps(real(T)) || j ≤ 1
        C *= (a+j)/(j+1)^2*(d+j)*w
        cⱼ += 2/(j+one(T))-one(T)/(a+j)+one(T)/(b-(j+one(T)))
        S += C*cⱼ
        j += 1
    end
    return S
end

function _₂F₁taylor(a::Number, b::Number, c::Number, z::Number)
    T = float(promote_type(typeof(a), typeof(b), typeof(c), typeof(z)))
    z₀ = abs(z) < 1 ? ρϵ*sign(z) : sign(z)/ρϵ
    q₀, q₁ = _₂F₁(a, b, c, z₀), a*b/c*_₂F₁(a+1, b+1, c+1, z₀)
    S₀, zz₀ = q₀, z-z₀
    S₁, zz₀j, j = S₀+q₁*zz₀, zz₀, 0
    while errcheck(S₀, S₁, 10eps(real(T))) || j ≤ 1
        q₀, q₁ = q₁, ((j*(2z₀-one(T))-c+(a+b+one(T))*z₀)*q₁ + (a+j)*(b+j)/(j+one(T))*q₀)/(z₀*(one(T)-z₀)*(j+2))
        zz₀j *= zz₀
        S₀, S₁ = S₁, S₁+q₁*zz₀j
        j += 1
    end
    return S₁
end

function _₁F₁maclaurin(a::Number, b::Number, z::Number)
    T = float(promote_type(typeof(a), typeof(b), typeof(z)))
    S₀, S₁, j = one(T), one(T)+a*z/b, 1
    while errcheck(S₀, S₁, 10eps(real(T))) || j ≤ 1
        rⱼ = (a+j)/((b+j)*(j+1))
        S₀, S₁ = S₁, S₁+(S₁-S₀)*rⱼ*z
        j += 1
    end
    return S₁
end

function _₃F₂maclaurin(a₁, a₂, a₃, b₁, b₂, z)
    T = float(promote_type(typeof(a₁), typeof(a₂), typeof(a₃), typeof(b₁), typeof(b₂), typeof(z)))
    S₀, S₁, j = one(T), one(T)+(a₁*a₂*a₃*z)/(b₁*b₂), 1
    while errcheck(S₀, S₁, 10eps(real(T))) || j ≤ 1
        rⱼ = ((a₁+j)*(a₂+j)*(a₃+j))/((b₁+j)*(b₂+j)*(j+1))
        S₀, S₁ = S₁, S₁+(S₁-S₀)*rⱼ*z
        j += 1
    end
    return S₁
end

function pFqmaclaurin(a::AbstractVector{S}, b::AbstractVector{U}, z::V) where {S, U, V}
    T = float(promote_type(S, U, V))
    S₀, S₁, j = one(T), one(T)+prod(a)*z/prod(b), 1
    while errcheck(S₀, S₁, 10eps(real(T))) || j ≤ 1
        rⱼ = inv(j+one(T))
        for i=1:length(a) rⱼ *= a[i]+j end
        for i=1:length(b) rⱼ /= b[i]+j end
        S₀, S₁ = S₁, S₁+(S₁-S₀)*rⱼ*z
        j += 1
    end
    return S₁
end

function continuedfraction(v::V, u::U, rtol::T) where {V<:Function, U<:Function, T<:Number}
    n = 4
    a, b = continuedfraction(v, u, n, 2n)
    n *= 2
    while n <= 2^16
        isapprox(a, b, rtol=rtol) && return b
        n *= 2
        a = deepcopy(b)
        b = continuedfraction(v, u, n)
    end
    isapprox(a, b, rtol=rtol) && return b
    error("Convergence failure for continued fraction for hypergeometric function")
end

function continuedfraction(v::V, u::U, n::Int, m::Int) where {V<:Function, U<:Function}
    @assert m > n
    output_m = u(m) / v(m)
    for i ∈ reverse(n+1:m-1)
        output_m = u(i) / (v(i) + output_m)
    end
    ui, vi = u(n), v(n)
    output_m = ui / (vi + output_m)
    output_n = ui / vi
    for i ∈ reverse(1:n-1)
        ui, vi = u(i), v(i)
        output_m = ui / (vi + output_m)
        output_n = ui / (vi + output_n)
    end
    vi = v(0)
    return vi + output_n, vi + output_m
end

function continuedfraction(v::V, u::U, n::Int) where {V<:Function, U<:Function}
    output = u(n) / v(n)
    for i ∈ reverse(1:n-1)
        output = u(i) / (v(i) + output)
    end
    return v(0) + output
end


text:
"""
Monte Carlo Planning solver with DPW and Bayesian Optimization Action Selection

Fields:

    depth::Int64:
        Maximum rollout horizon and tree depth.
        default: 10

    exploration_constant::Float64:
        Specified how much the solver should explore.
        In the UCB equation, Q + c*sqrt(log(t/N)), c is the exploration constant.
        default: 1.0

    n_iterations::Int64
        Number of iterations during each action() call.
        default: 100

    max_time::Float64
        Maximum amount of CPU time spent iterating through simulations.
        default: Inf

    k_action::Float64
    alpha_action::Float64
    k_state::Float64
    alpha_state::Float64
        These constants control the double progressive widening. A new state
        or action will be added if the number of children is less than or equal to kN^alpha.
        defaults: k:10, alpha:0.5

    check_repeat_state::Bool
    check_repeat_action::Bool
        When constructing the tree, check whether a state or action has been seen before (there is a computational cost to maintaining the dictionaries necessary for this)
        default: true

    tree_in_info::Bool:
        If true, return the tree in the info dict when action_info is called. False by default because it can use a lot of memory if histories are being saved.
        default: false

    rng::AbstractRNG:
        Random number generator

    estimate_value::Any (rollout policy)
        Function, object, or number used to estimate the value at the leaf nodes.
        If this is a function `f`, `f(mdp, s, depth)` will be called to estimate the value.
        If this is an object `o`, `estimate_value(o, mdp, s, depth)` will be called.
        If this is a number, the value will be set to that number.
        default: RolloutEstimator(RandomSolver(rng))

    init_Q::Any
        Function, object, or number used to set the initial Q(s,a) value at a new node.
        If this is a function `f`, `f(mdp, s, a)` will be called to set the value.
        If this is an object `o`, `init_Q(o, mdp, s, a)` will be called.
        If this is a number, Q will always be set to that number.
        default: 0.0

    init_N::Any
        Function, object, or number used to set the initial N(s,a) value at a new node.
        If this is a function `f`, `f(mdp, s, a)` will be called to set the value.
        If this is an object `o`, `init_N(o, mdp, s, a)` will be called.
        If this is a number, N will always be set to that number.
        default: 0

    default_action::Any
        Function, action, or Policy used to determine the action if POMCP fails with exception `ex`.
        If this is a Function `f`, `f(pomdp, belief, ex)` will be called.
        If this is a Policy `p`, `action(p, belief)` will be called.
        If it is an object `a`, `default_action(a, pomdp, belief, ex)` will be called, and if this method is not implemented, `a` will be returned directly.
        default: `ExceptionRethrow()`
"""
mutable struct BOMCPSolver <: Solver
    action_selector::ActionSelector
    belief_updater::Updater
    depth::Int64
    exploration_constant::Float64
    n_iterations::Int64
    max_time::Float64
    k_action::Float64
    alpha_action::Float64
    k_belief::Float64
    alpha_belief::Float64
    tree_in_info::Bool
    rng::AbstractRNG
    estimate_value::Any
    init_Q::Any
    init_N::Any
    default_action::Any
end

"""
    BOMCPSolver()

Use keyword arguments to specify values for the fields
"""
function BOMCPSolver(
                    action_selector::ActionSelector,
                    belief_updater::Updater
                    ;depth::Int=10,
                    exploration_constant::Float64=1.0,
                    n_iterations::Int=100,
                    max_time::Float64=Inf,
                    k_action::Float64=10.0,
                    alpha_action::Float64=0.5,
                    k_belief::Float64=10.0,
                    alpha_belief::Float64=0.5,
                    tree_in_info::Bool=false,
                    rng::AbstractRNG=Random.GLOBAL_RNG,
                    estimate_value::Any = RolloutEstimator(RandomSolver(rng)),
                    init_Q::Any = 0.0,
                    init_N::Any = 0,
                    default_action::Any = ExceptionRethrow()
                   )
    BOMCPSolver(action_selector, belief_updater, depth, exploration_constant, n_iterations, max_time, k_action, alpha_action, k_belief, alpha_belief, tree_in_info, rng, estimate_value, init_Q, init_N, default_action)
end

mutable struct BOMCPPlanner{P<:Union{MDP,POMDP}, B, A, SE, RNG} <: Policy
    solver::BOMCPSolver
    p::P
    tree::Union{Nothing, BOMCPTree{B,A}}
    solved_estimate::SE
    rng::RNG
    gp::Union{Dict, GPLA}
end

function BOMCPPlanner(solver::BOMCPSolver, p::P) where P<:Union{POMDP,MDP}
    se = convert_estimator(solver.estimate_value, solver, p)
    gp = solver.action_selector.gp
    return BOMCPPlanner{P,
                      BOMCP.belief_type(solver.belief_updater),
                      actiontype(P),
                      typeof(se),
                      typeof(solver.rng)}(solver,
                                          p,
                                          nothing,
                                          se,
                                          solver.rng,
                                          gp
                     )
end

Random.seed!(p::BOMCPPlanner, seed) = Random.seed!(p.rng, seed)
POMDPs.solve(solver::BOMCPSolver, p::Union{POMDP,MDP}) = BOMCPPlanner(solver, p)


text:
using LinearAlgebra
using JuLIP, JuLIP.Potentials, SKTB, SparseArrays
using JuLIP.Potentials: site_energy, site_energy_d
TB = SKTB

@testset "Virial" begin

# gamma point only
beta = 10.0        # temperature / smearing paramter: 10 to 50 for room temperature
DIM = (2,2,2)
# define the model
tbm = TB.NRLTB.NRLTBModel(:Si, TB.FermiDiracSmearing(beta), bzquad = TB.GammaPoint() )
at = DIM * bulk(:Si, pbc=(true, true, true), cubic = true)
# JuLIP.rattle!(at, 0.02)
F = cell(at)'
pert = zeros(Float64, 3, 3)
@show length(at)
println("Finite difference test for virial, with Gamma point only")
E = energy(tbm, at)
vir = virial(tbm, at)
errors = Float64[]
for p = 2:9
    h = 0.1^p
    dEh = zero(vir)
    for n = 1:length(F)
        pert[n] += h
        Defm = I + pert * inv(F)
        apply_defm!(at, Defm)
        # set_cell!(at, (F+pert)'; updatepositions=true)
        Eh = energy(tbm, at)
        apply_defm!(at, inv(Defm))
        dEh[n] = (Eh - E) / h
        pert[n] -= h
    end
    push!(errors, norm(vir + dEh, Inf))
    println( " ", p, " | ", errors[end])
end
println(@test minimum(errors) < 1e-3 * maximum(errors))

# multiple k-points
nquad = 2
bzquad = TB.MPGrid(at, (nquad, nquad, nquad))
tbm = TB.NRLTB.NRLTBModel(:Si, TB.FermiDiracSmearing(beta), bzquad = bzquad)
DIM = (2,2,2)
at = DIM * bulk(:Si, pbc=(true, true, true), cubic = true)
# JuLIP.rattle!(at, 0.02)
F = cell(at)'
pert = zeros(Float64, 3, 3)
@show length(at)

println("Finite difference test for virial, with ", nquad, " points in each direction")
E = energy(tbm, at)
vir = virial(tbm, at)
errors = Float64[]
for p = 2:9
    h = 0.1^p
    dEh = zero(vir)
    for n = 1:length(F)
        pert[n] += h
        Defm  = I + pert * inv(F)
        apply_defm!(at, Defm)
        # set_cell!(at, (F+pert)'; updatepositions=true)
        Eh = energy(tbm, at)
        dEh[n] = (Eh - E) / h
        apply_defm!(at, inv(Defm))
        pert[n] -= h
    end
    push!(errors, norm(vir + dEh, Inf))
    println( " ", p, " | ", errors[end])
end
println(@test minimum(errors) < 1e-3 * maximum(errors))

end


text:
function compute_reservation_wage(mcm; return_values = false)
    V, U = solve_mccall_model(mcm)
    w_idx = searchsortedfirst(V .- U, 0)

    if w_idx == length(V)
        w_bar = Inf
    else
        w_bar = mcm.w_vec[w_idx]
    end

    if return_values == false
        return w_bar
    else
        return w_bar, V, U
    end
end

max_stars_repo_name: yywingliang/lecture-source-jl

max_issues_repo_name: yywingliang/lecture-source-jl

max_forks_repo_name: yywingliang/lecture-source-jl


text:
module CollidingCylindersSimulation

using  Unitful
using  ModiaLang
import Modia3D
using  Modia3D.ModiaInterface

vmatRed   = Modia3D.VisualMaterial(color="Red")
vmatGreen = Modia3D.VisualMaterial(color="Green")
vmatBlue  = Modia3D.VisualMaterial(color="Blue")
vmatGrey  = Modia3D.VisualMaterial(color="Grey", transparency=0.5)

LxTable = 2.0
LyTable = 2.0
LzTable = 0.02
diameterCyl      = 0.1
lengthCyl        = 0.5
innerDiameterCyl = 0.08

collidingCylinders = Model(
    world = Object3D(feature=Scene(gravityField=UniformGravityField(g=9.81, n=[0, 0, -1]),
                                   enableContactDetection=true)),

    cylinderX = Object3D(feature=Solid(shape=Cylinder(axis=1, diameter=diameterCyl, length=lengthCyl, innerDiameter=innerDiameterCyl),
                                       solidMaterial="BilliardBall",
                                       visualMaterial=vmatRed,
                                       collision=true)),
    jointX = FreeMotion(obj1=:world, obj2=:cylinderX,
                        r=Var(init=[-0.3, -0.1, 0.5]),
                        rot=Var(init=[0.0, -90*u"°", -80*u"°"])),

    cylinderY = Object3D(feature=Solid(shape=Cylinder(axis=2, diameter=diameterCyl, length=lengthCyl, innerDiameter=innerDiameterCyl),
                                       solidMaterial="BilliardBall",
                                       visualMaterial=vmatGreen,
                                       collision=true)),
    jointY = FreeMotion(obj1=:world, obj2=:cylinderY,
                        r=Var(init=[0.0, -0.1, 0.5]),
                        rot=Var(init=[-10*u"°", 0.0, 0.0])),

    cylinderZ = Object3D(feature=Solid(shape=Cylinder(axis=3, diameter=diameterCyl, length=lengthCyl, innerDiameter=innerDiameterCyl),
                                       solidMaterial="BilliardBall",
                                       visualMaterial=vmatBlue,
                                       collision=true)),
    jointZ = FreeMotion(obj1=:world, obj2=:cylinderZ,
                        r=Var(init=[0.3, -0.1, 0.5]),
                        rot=Var(init=[80*u"°", 0.0, 0.0])),

    cylinder2 = Object3D(parent=:world,
                         translation=[0.0, 0.1, 0.2],
                         feature=Solid(shape=Cylinder(axis=1, diameter=diameterCyl, length=1.0),
                                       solidMaterial="BilliardBall",
                                       visualMaterial=vmatGrey,
                                       collision=true)),
    cylinder3 = Object3D(parent=:world,
                         translation=[0.0, -0.1, 0.0],
                         feature=Solid(shape=Cylinder(axis=1, diameter=diameterCyl, length=1.0),
                                       solidMaterial="BilliardBall",
                                       visualMaterial=vmatGrey,
                                       collision=true))
)

CollidingCylinders = @instantiateModel(buildModia3D(collidingCylinders), unitless=true, log=false, logStateSelection=false, logCode=false)

stopTime = 0.7
interval = 0.001
tolerance = 1e-8
finalStates = [-0.3000014767733773, -1.082193896696148, -0.538807311042302, -3.7592282983652485e-6, -2.5117143103911763, -3.747655431035406, 3.434780596077875, -1.557228662514374, 1.617103948525645, 0.03212608204313396, 0.005364383932787988, 0.5421088990973728, -3.042579138039232e-7, -1.0813130417689532, -0.5361375150624056, -7.007175995563201e-7, -2.5091653603011435, -3.740308350308572, 0.27706268760792585, 0.011027339253852072, -0.0002832776767408359, -0.4592873394400626, 0.021172006600497985, -0.005240358556866595, 0.3000003035109844, -1.0823403796115827, -0.5378963004593084, 8.132044608411955e-7, -2.5121486704670146, -3.7451155263449594, 1.8287413219661708, 0.0005856583646075422, -0.019142491087689967, -0.511322469433515, -0.007676420293499906, -0.04351478584418659]
simulate!(CollidingCylinders, stopTime=stopTime, interval=interval, tolerance=tolerance, log=true, logStates=true, logEvents=true, requiredFinalStates=finalStates)

@usingModiaPlot
plot(CollidingCylinders, [("jointX.r", "jointY.r", "jointZ.r") ("jointX.rot", "jointY.rot", "jointZ.rot")
                                    ("jointX.v", "jointY.v", "jointZ.v") ("jointX.w", "jointY.w", "jointZ.w")], figure=1)

end


text:
# SINGLE NODE/CORE VERSION

import Pkg; Pkg.instantiate()

include("/users/yh31/scratch/projects/gigaword_64k/src/gigaword_64k.jl")

const path_apw_dir = "/users/yh31/scratch/datasets/entity_linking/raw_data/gigaword/giga5/data/apw_eng"

const path_intermed_data = "/users/yh31/scratch/datasets/entity_linking/raw_data/gigaword/giga_topk/intermediate_files"
const path_intermed_apw = joinpath(path_intermed_data, "apw")

process_part_of_tree(path_apw_dir, path_intermed_apw, 4)


text:
using SimultaneousPerturbationStochasticApproximation
using Test, Random, Statistics, LinearAlgebra

using SimultaneousPerturbationStochasticApproximation: estimate_grad, init!, LearningRateUpdate, update!, ConvergenceWelchTest, converged!

@testset "estimate gradient" begin
f(x) = sum(abs2, x) + .2*randn()
@test mean(estimate_grad(Random.GLOBAL_RNG, [2, 1], .01, f)[1] for _ in 1:10^5) ≈ [4, 2] atol = .5
end # estimate gradient tests

@testset "heuristic" begin
f(x) = sum(abs2, x) + .2*randn()
h = SimpleHeuristic()
init!(Random.GLOBAL_RNG, h, f, [2, 1], 10^4, [.2, .2])
@test h.c ≈ .2 atol = .15

h = SimpleHeuristic()
init!(Random.GLOBAL_RNG, h, x -> sum(abs2, x), [2, 1], 10^4, [.2, .2])
@test h.c == h.c_min

h = SimpleHeuristic(c = .1)
init!(Random.GLOBAL_RNG, h, f, [2, 1], 10^4, [.2, .2])
@test h.c == .1

h = SimpleHeuristic(elementwise = true)
init!(Random.GLOBAL_RNG, h, f, [2, 1], 10^4, [.2, .2])
@test length(h.a) == 2

# Update heuristics
f(x) = dot([100, 1], x) + .1 * randn()
lb = fill(-10, 2)
ub = fill(10, 2)
h = SimpleHeuristic(elementwise = true,
                    update = LearningRateUpdate(lb, ub, N = 10))
theta = [2., 1]
init!(Random.GLOBAL_RNG, h, f, theta, 10^4, [.2, .2])
olda = copy(h.a)
for _ in 1:10
    theta .*= .8 # trend
    update!(h, theta)
end
@test (&)((h.a .> olda)...)

Random.seed!(2);
olda = copy(h.a)
for _ in 1:10
    theta  = 100*randn(2) # span
    update!(h, theta)
end
@test (&)((h.a .< olda)...)
end # heuristic tests

@testset "convergence" begin
f(x) = sum(abs2, x) + .1*randn()
Random.seed!(2);
c = ConvergenceWelchTest(N = 1, K = 500)
global ret = false
for _ in 1:10
    global ret = converged!(c, f, [0, 0]) # converged
end
@test ret == true

c = ConvergenceWelchTest(N = 10)
global ret = false
theta = [5., 1]
for _ in 1:41
    theta .*= .98 # not converged
    global ret = converged!(c, f, theta)
end
@test ret == false

spsa = SPSA(lower = fill(-10, 2), upper = fill(10, 2),
            convergencetest = ConvergenceWelchTest())
res = minimize!(spsa, x -> sum(abs2, x), maxfevals = 10^8, verbose = false)
@test res[1].evals < 10^8

spsa = SPSA(lower = fill(-10, 2), upper = fill(10, 2),
            convergencetest = ConvergenceWelchTest())
res = minimize!(spsa, f, maxfevals = 10^8, verbose = false)

end # convergence

@testset "minimize" begin
spsa = SPSA(lower = fill(-10, 2), upper = fill(10, 2))
res = minimize!(spsa, x -> sum(abs2, x), maxfevals = 10^3)
@test res[1].x ≈ zeros(2) atol = 1e-12

res = minimize(x -> sum(abs2, x), lower = fill(-10, 2), upper = fill(10, 2),
               maxfevals = 10^3)
@test res[1].x ≈ zeros(2) atol = 1e-12
end


text:
using YaoExtensions, YaoPlots, Yao
using Compose, Cairo

_save(str) = PNG(joinpath(@__DIR__, str))

YaoPlots.CircuitStyles.linecolor[] = "#000000"
# qft circuit
vizcircuit(qft_circuit(5)) |> _save("qft.png")

# variational circuit
vizcircuit(variational_circuit(5)) |> _save("variational.png")
# vizcircuit(variational_circuit(5; mode=:Merged))

# general U4 gate
vizcircuit(general_U4()) |> _save("u4.png")

# quantum supremacy circuit
vizcircuit(rand_supremacy2d(2, 2, 8)) |> _save("supremacy2d.png")

# google 52 qubit
vizcircuit(rand_google53(10)) |> _save("google53.png")

# control blocks
vizcircuit(chain(control(5, (2,-3), 4=>X), control(5, (-4, -2), 1=>Z))) |> _save("controls.png")

# controlled kron
control(4, 2, (1, 3)=>kron(X, X)) |> vizcircuit |> _save("cxx.png")

control(4, -collect(1:4-1), 4=>-Z) |> vizcircuit |> _save("reflect.png")

chain(5, [put(5, 2=>X), Yao.Measure(5; locs=(2,3)), Yao.Measure(5;locs=(2,)), Yao.Measure(5; resetto=bit"00110")]) |> vizcircuit |> _save("measure.png")

chain(5, [put(5, 2=>ConstGate.Sdag), put(5, 3=>ConstGate.Tdag),
    put(5, (2,3)=>ConstGate.CNOT), put(5, (1,4)=>ConstGate.CZ), put(5, (1,2,5)=>ConstGate.Toffoli),
    put(5, (2,3)=>ConstGate.SWAP), put(5, (1,)=>ConstGate.P0), put(5, (1,)=>ConstGate.I2),
    put(5, (2,)=>ConstGate.P1), put(5, (1,)=>ConstGate.Pu), put(5, (1,)=>ConstGate.Pd),
    put(5, (2,)=>ConstGate.T),
    ]) |> vizcircuit |> _save("constgates.png")

chain(5, [put(5, (2,3)=>label(SWAP, "SWAP")'), put(5, 2=>label(I2, "id")), put(5, 2=>label(X, "X")), control(5, (5,3), (2,4,1)=>put(3, (1,3)=>label(SWAP, "SWAP")))]) |> vizcircuit |>  _save("multiqubit.png")

YaoPlots.CircuitStyles.linecolor[] = "#FFFFFF"
# qft circuit
vizcircuit(qft_circuit(5)) |> _save("qft-white.png")