text:
"""
Metadata
    - number of fp in main and borders and associated counters (counters will be updated at each dilatation step) and old result counters saving information about counter value in previous dilatation step
    - number of fn in main and borders and associated counters (counters will be updated at each dilatation step) and old result counters saving information about counter value in previous dilatation step
    - offset where is the results for main and borders for both main and borders
    - numbers pointing out begining and end of all rusults associated with this block
    - information is it active  or full or inactive or is it to be activated (this last one set to true if other block set any result through the padding)

1) we iterate through the data of both segm and gold standad array and  proces it - w build up new  3d boolean arrays and we updata data in metadata
    we add the number of false positives and false negatives  to each block with division  for each border so top bottom, left and right anterior posterior
    - ths data will enable skipping some paddings checks in case of empiness and will enable preserving uniqness in case we will want to set new result from both neighbuing blocks ..
    - all of it will be just primary metadata that will be processed later
    - we will also look for in this step for max and min vales of each dimension x,y,z
2) on cpu we create two result queues that will be 1,5 times longer that fp in one case and fn in the other
3) first metadata pass we add to the metadata offset where each block will put its main results and padding results - so all will be stored in result quueue but in diffrent spots
    we need to make ques for paddings longer than number of possible results becouse of possible modifications from neighbouring blocks that can happen simultanously
    we will establish this offsets using atomics- at this pass we will also prepare first work queue with indicies of metadata blocks and booleans indicating is it related to gold pass dilatation step or other pass
4) we do the dilatations in all blocks from work queue and record every time we had stumbled upon newly covered result - in order to later being able to establish precise distance
    we will record in result list not only coordinates but also direction
    - results in main part will be put in the part of result list associated with main block and those from paddings to border lists of surrounding blocks
- block in the end will check is it full if it is it will set such data to metadata and add information that it is inactive
5) we do the metadata pass we analyze only those blocks that are in the borders of intrests - max min x y z was specified in step 1 , we check weather block is set to be activated
 and in case it is not full we will make it active , if the block is set as active we just add it to work queue  that is appropriate to next iteration
  - we need also to scan  the border result ques if there is any duplicate result - if so we set it to zeros and we reduce the border result counter
        of course we check it only in case new counter is bigger than old counter
  concurrently we do some housekeeping so we check weather counters for result are not  indicating that we already finished

Work will have 3 parts - 3 kernels
    1) will prepare data - will return source array in boolean format and will give benchmark
    max and min x,y,z of either mask - to describe smallest possible cube holding all necessary data
    2)first pass through the data  that will mark which data blocks are acteve or full ...  - it will also prepare  the rudamentary worplan for first iteration in phase 3
    3)the actual working part - there will be shceduling block that will pass what work should be done by other blocks - copying da, deciding to terminate etc.
        - rest of the block will iteratively analyze data blocks scheduled by scheduling block

So we will need the block of data representing each wall of the bit cube bit cuve will have edge length of 32 - apart from corner cases

Simplification
Get reduced array of bits from Tpfp..fn modified kernel
Assign blocks of size 32x32x32 for each we will have data structure holding data as pointed below all will be in UInt16
    Is activeOrFullForFirsPassgold - 2 if neither in first mask , 1 if full in first mask; 0 if active in first mask
    Is activeOrFullForFirsPasssegm - 2 if neither in first mask , 1 if full in first mask; 0 if active in first mask
    Is activeOrFullForSecondPassgold - 2 if neither in first mask , 1 if full in first mask; 0 if active in first mask
    Is activeOrFullForSecondPassSegm - 2 if neither in first mask , 1 if full in first mask; 0 if active in first mask
   activeIterNumbMaskA - The number of iteration when block was last active - needed to coordinate which block will be working now
    activeIterNumbMaskB - The number of iteration when block was last active - needed to coordinate which block will be working now
    Data will be stored in UInt16 4 dimensional array in such a way that data blocks block id x,y,z will point to metadata position

ResArray storing all coordinates that were covered and in which iteration it occured  and in pass from which mask it occured

work scheduling structures

workSchedule- matrix that will have  number of rows that is equal to number of thread blocks that will work on a task
    and in each column we will put metadata needed to localize on what data block we are working currently
    - so x,y,z indexes and 1 if we do it from perspective of gold sandard mask and 2 if from second mask
LocalRes -  matrix that will have  number of rows that is equal to number of thread blocks that will work on a task - so blocks will write results to this array
        and then scheduling block will collect those results from all matrix in the ResArray
LocalResLastEntryList - list with length equal to the number of thread blocks - will store the last number of occuppied  spot in Local Res (so block will know where to put next results)
isWorking - bollean list of length equal the number of thread blocks that will point out wheather given thread block is working or not
            - data needed for scheduling block
fp, fn - number of false positive and false negatives - we will get it from algorithm responsible for preparing data
    - and it will be needed  to know how much data needs to be allocated andfor the scheduling block to be able to decide whether we finished

Algo
Check blocks to assign weather they become inactive or full
    We do it by storing in register UInt32 (UInt numb) and then when streaming (each lane is responsible for 32 fields) will modify the bit of stored numb associated with step in a loop we are ; at first pass we do it for both masks at the same time
    Next we can use voting mechanism  and then shared memory reduction to check 2 times weather all numbs are first all 0’s - so numb will be 0 - in this case we will make it inactive second all 1’s so numbs will be max value of UINt32  - in this case we will make it full - we apply the procedure to both masks if neither we will make it active
    We scan in blocks for next active block available and mark using atomics by setting the activeIterNumb of a block to the current iteration number, so no other block can modify it further  !!! we need to check what does the atomic return …
    When working on a block we have 3 jobs
    set to the 1’s all values around each one we find - and write it to global memory immediately or first to shared memory - to be experimented upon
    Add coordinates of covered values to ResArray
    Activate neighbouring blocks if needed
    Activation of neighbour blocks
    When we found the 1 (before dilatation step) in position of max/ min x/y/z we need to save information to the shared memory that this block needs to be activated according to the data present - we will choose the block depending on which dimension was min/max if we have corners or edges we will activate multiple blocks
    Not in all cases the block will get activated it will not get activated if
    Is full
    Neighbours in given direction are all inactive in second mask
    We finish when there are no more active blocks

Consideration for fact that Housdorff needs to be scanned from 2 perspectives
we need to have a copy of original data for both masks so modification of one pass of housdorff will not affect the other - in this case also order is not so important
Each block needs to be aware for which pass it is currently working - hence to know what to dilatate and where to write the results
To maximize occupancy the blocks will concurrently work on both passes at once - it will reduce number of grid sync call at each iterations - and increase occupancy at the end of algorithm

"""
module SimplerHousdorff
using CUDA, ..CUDAGpuUtils, Logging
"""
Metadata
divide array into blocks and assigns the metadata for each
    Is activeOrFullForFirsPassgold - 2 if neither in first mask , 1 if full in first mask; 0 if active in first mask
    Is activeOrFullForFirsPasssegm - 2 if neither in first mask , 1 if full in first mask; 0 if active in first mask
    Is activeOrFullForSecondPassgold - 2 if neither in first mask , 1 if full in first mask; 0 if active in first mask
    Is activeOrFullForSecondPassSegm - 2 if neither in first mask , 1 if full in first mask; 0 if active in first mask
    activeIterNumbMaskA - The number of iteration when block was last active - needed to coordinate which block will be working now
    activeIterNumbMaskB - The number of iteration when block was last active - needed to coordinate which block will be working now
    Data will be stored in n times 34 UInt16 matrix in global memory

So in order to achieve it we need first to get the dimensions of the array we got - as it is reduced array it do not have to have the dimensions that would be simmilar
to the medical image dimensions - yet by construction all dimensions need to be multiple by 32!!!
As we start from place where x,y,z is 0 and we will proceed (concurrently from here )
        we know by construction that min x, min y, and min z will be always some multiply of 32 away from begining
        now we can lounch multiple small blocks where each will calculate the the required metadata - amount of those blocks will be just equal to the number of blocks tha we want to describe
        and will be calculated by ceildiv of given dimension so sth like sizz = size[arrGold]-> will be mapped
        to blocks as sizz-> (celdiv(sizz[1]) ...) by all dimensions ; blocks will be elongated in one dimension and constructed in a way that single warp will be responsible to set metadata for single block

 arrGold - boollean 3 dim Cu array with gold standard (reduced to the smallest block with all true entries of both masks)
 arrSegm - segmentation 3 dim boolean Cu array we want to compare

"""
function getBlockMetaData(arrGold, arrSegm,maxThrPerB::Int64=1024)
    arrGoldDims= size(arrGold)
    metadata = allocateMomory(arrGoldDims)
    metaDataDims = size(metadata)
    blocks,threadNumPerBlock = getBlockNumb(metaDataDims)# for metadata kernel
    @info "blocks " blocks
    @info "threadNumPerBlock " threadNumPerBlock

    @cuda threads=1024 blocks=metaDataDims[1]*metaDataDims[2]*metaDataDims[3] housedorffMetadataKernel(metadata,metaDataDims,arrGoldDims )
    return metadata
end#getBlockMetaData

"""
kernel that output metadata for Housedorff
first all blocks are active - we intend to have all blocks o be active at first pass - and then prograssively  work only on this data that we are intrested in
matadata - 4 dimensional data with block metadata

dataArrs - array of 2 arrays where first is arrGold and second arrSegm
    arrGold - boollean 3 dim Cu array with gold standard (reduced to the smallest block with all true entries of both masks)
    arrSegm - segmentation 3 dim boolean Cu array we want to compare

iterationNumber - variable that we will set in order to mark on what iteration we are curently

arrGoldDims - dimensions of the main array
dimOfThreadBlock - how big is the edge of a cube describing data block - for example for 1024 threads - we will have 32x32x32 size hence dimOfThreadBlock= 32
we will have in basic type 32x32 threads and work on 32x32x32 data block
"""
function housedorffMetadataKernel(matadata
                                ,arrGold
                                ,arrSegm
                                ,arrGoldDims::Tuple{Int64, Int64, Int64}
                                ,iterationNumber::UInt16
                                ,debugginArr
                                ,dimOfThreadBlock::UInt16 )

    #initializing shared array (remember this is autmatically initated to semi random numbers)
    shemm= CuDynamicSharedArray(Bool,(dimOfThreadBlock,dimOfThreadBlock,dimOfThreadBlock) )
    #as we constructed data dimensions o be always multiple of 32 we do not need to do bound checks
    @unroll for k in 1:32
        shemm[blockIdx().x,blockIdx().y,k]=  [ (blockIdx().x-1)*dimOfThreadBlock+threadIdxX(),(blockIdx().y-1)*dimOfThreadBlock+threadIdxY(),(blockIdx().z-1)*dimOfThreadBlock+k]
    end#for

                                # each lane will be responsible for one meta data
    # @unroll for k in 0:metadataDims[1]
    #     matadata[threadIdxX(),k+1, blockIdx().x,1]=(threadIdxX()-1)*32   #min x
    #     matadata[threadIdxX(),k+1, blockIdx().x,2]=min(threadIdxX()*32,arrGoldDims[1])   #max x
    #     matadata[threadIdxX(),k+1, blockIdx().x,3]=k*32   #min y
    #     matadata[threadIdxX(),k+1, blockIdx().x,4]=min((k+1)*32,arrGoldDims[2] )   #max y
    #     matadata[threadIdxX(),k+1, blockIdx().x,5]=(blockIdx().x-1)*32   #min z
    #     matadata[threadIdxX(),k+1, blockIdx().x,6]=min((blockIdx().x)*32,arrGoldDims[3])   #max z
    # end#for
end #housedorffMetadataKernel

"""
so given metadata and some helper constants we need to choose on what data our block of threads will work on
at the begining - first iteration it is not very important as all blocks are active - hence we can schedule work by just evenly dividing all data blocks to thread blocks
Then- some of the blocks will be marked as either inactive or full hence  we need to redistribute the work - so each thread block will have the same (+-1) number of data blocks to work on
    - in order to achieve this -  single block will be scheduled to redistribute work - it will run at the very end of iteration cycle  to redistribute work  for next one
    ; and at the begining - to redistribute blocks activated  in last round of previous iteration
"""
function chooseBlockToWorkOn()

end#chooseBlockToWorkOn

"""
controls allocation of GPU memory - instantiating Cu arrays
 first we alocate the place for  metadata -

   Metadata -  Assign blocks of size 32x32x32 for each we will have data structure holding data as pointed below all will be in Bool - so 4 bytes of metadata per block
      1) isActiveOrFullForPasssegm - true if other  mask is acive for modifications
      2) isActiveOrFullForPassgold -  true if gold standard mask is acive for modifications
      3) isFullSegm - true if other mask is full (only ones)
      4) isFullGold - true if gold standard mask is full (only ones)

      Data will be stored in Bool 4 dimensional array in such a way that data blocks block id x,y,z will point to metadata position

    resArray storing all coordinates that were covered and in which iteration it occured  and in pass from which mask it occured

    work scheduling structures

    mainWorkQueue - basically we need to add to the one dimensional queue all of the  active blocks we find in first pass- crerating basic queue that will be processed by normals passes
    of course during those normal passes some blocks will be added to the queue and some will be removed in order to mage
        it will store 3 indicies (UInt8) of the place of the block in metadata plus  wheater we are referencing main pass or second one (UInt8 wchich will be 1 or 0)
        so entry1 - x position in metadata 2)y pos 3) z pos 4) boolean marking is it gold pass or not
        it we would need additional helper structures
    so we will have:
    mainQuesCounter - telling us how many entries we have in work queue - we will divide this by number of thread blocks + some constant
        - to get some amount of the data blocks to be processed + tail queue that will be accessed in atomic way - but it can be accessed by all blocks
        so if some block will finish work before others it will start processing this tail queue
    mainActiveCounterNow,mainActiveCounterNext - at first main queue will have only active blocks but progressively it will have more and more empty spots
        so we need to get the second counter that will keep track on the ramaining active blocks in curent iteration
            -mainActiveCounterNow - will be reduced every time thread block finish processing block - we will sync grid and start next iteration when it will reach 0
            -mainActiveCounterNext - will be increased every time we activate some block - will become the  mainActiveCounterNow in next iteration if it will reach 0 we will
                call it the end and finish kernel

    fp, fn - number of false positive and false negatives - we will get it from algorithm responsible for preparing data
        - and it will be needed  to know how much data needs to be allocated andfor the scheduling block to be able to decide whether we finished

   arguments
    fpNumb,fnNumb - number of false positives and false negatives
    blocksNum - number of thread blocks that are able to run on the same grid
    arrGoldDims - tuple with dimensions of main data array
    dataBlocksNum- number of dataBlocks
    dimOfThreadBlock - how big is the edge of a cube describing data block - for example for 1024 threads - we will have 32x32x32 size hence dimOfThreadBlock= 32
"""
function allocateMomory(arrGoldDims::Tuple{Int64, Int64, Int64}
                        ,fpNumb::Int64
                        ,fnNumb::Int64
                        ,blocksNum::Int64
                        ,dataBlocksNum::Int64
                        ,dimOfThreadBlock::UInt16  )
    maxresultPoints = fpNumb+fnNumb+1
    #number of data blocks in given dimension
    x=cld(arrGoldDims[1],dimOfThreadBlock)
    y=cld(arrGoldDims[2],dimOfThreadBlock)
    z=cld(arrGoldDims[3],dimOfThreadBlock)

    #we need just some blocks to set in the schedule for the begining - the scheduling block will refine it
    # workNumbPerBlock = cld(x*y*z,dataBlocksNum)
    # workScheduleCPU = zeros(UInt16,blocksNum, workNumbPerBlock  ,4)
    # Threads.@threads for i in 1:blocksNum
    #                         for j in 1:workNumbPerBlock
    #                         workScheduleCPU
    #                         end
    #                     end

    workSchedule= CUDA.zeros(UInt16,blocksNum, cld(x*y*z,dataBlocksNum)  ,4) #in each entry 1) x;2)y;3)z;4)1 if gold standard pass and 2 if segm pass  ; length is set so we will have approximately equal number of blocks to work on

    metaData= CUDA.zeros(UInt16,x,y,z,12)
    resArray=CUDA.zeros(UInt16, blocksNum,maxresultPoints, 4)#in each entry 1) x;2)y;3)z;4)1 if gold standard pass and 2 if segm pass
    localRes= CUDA.zeros(UInt16,maxresultPoints/2,4  ) #maxresultPoints/2 very conservative  we may experiment with far smaller number to  decrese memory usage
    localResLastEntryList= CUDA.zeros(UInt16, blocksNum) # 1 entry per thread block
    innerLoopStep= CUDA.zeros(UInt16, blocksNum) # 1 entry per thread block
    worksScheduleLastStep= CUDA.zeros(UInt16, blocksNum) # 1 entry per thread block
    isWorking = CUDA.zeros(Bool, blocksNum) # 1 entry per thread block
    # for boolean arrays that will store all boolean data about analyzed array#we have two copies as Housdorff is composed of 2 passes
    reducedGoldA= CUDA.zeros(Bool,arrGoldDims)
    reducedSegmA =CUDA.zeros(Bool,arrGoldDims)

    reducedGoldB= CUDA.zeros(Bool,arrGoldDims)
    reducedSegmB =CUDA.zeros(Bool,arrGoldDims)

return (metaData)
end#allocateMomory

end#SimplerHousdorff

    # #first we choose biggest dimension and this will be dimension of elongated block
    # arrGoldDims= size(arrGold)
    # maxDim = argmax(size(arrGold))
    # localArr = [0,0,0]
    # for i in 1:3
    #     if(i!=maxDim)
    #         localArr[i]= cld(arrGoldDims[i],32)
    #     else
    #         #we are now in max dimension
    #         localArr[i]= cld(arrGoldDims[i],maxThrPerB)
    #     end#if
    # end#for

# """
# getting dimensions  needed to run number of blocks that will cover all of data to prepare metadata
# metasataDims - dimensions of metasata array - 4 dims tuple
# maxThrPerB - maximum number of threads per block defoult is 1024
#     return the information to the kernel how to create its blocks of threads (how many); and the number of threads in a block
# """
# function getBlockNumb(metasataDims::Tuple{Int64, Int64, Int64, Int64},maxThrPerB::Int64=1024)

#     return ( metasataDims[3] ,maxThrPerB)
# end#getBlockNumb

# """
# kernel that output metadata for Housedorff
#     matadata - Int16 4 dimensional array where x,y,z dims identical to data blocks  and 4 dim stores 6 numbers
#         Is activeOrFullForFirsPassgold - 2 if neither in first mask , 1 if full in first mask; 0 if active in first mask
#         Is activeOrFullForFirsPasssegm - 2 if neither in first mask , 1 if full in first mask; 0 if active in first mask
#         Is activeOrFullForSecondPassgold - 2 if neither in first mask , 1 if full in first mask; 0 if active in first mask
#         Is activeOrFullForSecondPassSegm - 2 if neither in first mask , 1 if full in first mask; 0 if active in first mask
#         activeIterNumbMaskA - The number of iteration when block was last active - needed to coordinate which block will be working now
#         activeIterNumbMaskB - The number of iteration when block was last active - needed to coordinate which block will be working now
#         metadataDims - dimensions of metadata - 4 dims tuple
#         arrGoldDims - dimensions of the main array
#     IMPORTANT - althought this is 1 dimensional kernel we will output 3 dim data
#         so blockIdx().x - gives info about z - dimension
#         threadIdxX() - tell about position in x dimension
#         in y dimension we will iterate in a loop the number that is equal to the number of the size of the  y dimension (remember we are talking about metadata array)

# """
# function housedorffMetadataKernel(matadata
#                                 ,metadataDims::Tuple{Int64, Int64, Int64, Int64}
#                                 ,arrGoldDims::Tuple{Int64, Int64, Int64} )
#     # each lane will be responsible for one meta data
#     @unroll for k in 0:metadataDims[1]
#         matadata[threadIdxX(),k+1, blockIdx().x,1]=(threadIdxX()-1)*32   #min x
#         matadata[threadIdxX(),k+1, blockIdx().x,2]=min(threadIdxX()*32,arrGoldDims[1])   #max x
#         matadata[threadIdxX(),k+1, blockIdx().x,3]=k*32   #min y
#         matadata[threadIdxX(),k+1, blockIdx().x,4]=min((k+1)*32,arrGoldDims[2] )   #max y
#         matadata[threadIdxX(),k+1, blockIdx().x,5]=(blockIdx().x-1)*32   #min z
#         matadata[threadIdxX(),k+1, blockIdx().x,6]=min((blockIdx().x)*32,arrGoldDims[3])   #max z
#     end#for
# end #housedorffMetadataKernel


text:
# The fraction 49/98 is a curious fraction, as an inexperienced mathematician
# in attempting to simplify it may incorrectly believe that 49/98 = 4/8, which
# is correct, is obtained by cancelling the 9s.
#
# We shall consider fractions like, 30/50 = 3/5, to be trivial examples.
#
# There are exactly four non-trivial examples of this type of fraction, less
# than one in value, and containing two digits in the numerator and denominator.
#
# If the product of these four fractions is given in its lowest common terms,
# find the value of the denominator.

using ProjectEulerSolutions

# Cycle through numbers using Julia's built-in digits, and store the fractions
# as rationals.  Julia's built-in rationals solve the common denominator
# problem easily.
function p033solution()::Integer
    rationals = Rational[]
    for a in 10:98
        for b in (a+1):99
            if a != b
                val = a / b
                da = digits(a)
                db = digits(b)
                if da[1] == db[1] && da[1] != 0 && da[2] / db[2] == val
                    push!(rationals, a // b)
                elseif da[1] == db[2] && da[2] / db[1] == val
                    push!(rationals, a // b)
                elseif da[2] == db[1] && da[1] / db[2] == val
                    push!(rationals, a // b)
                elseif da[2] == db[2] && da[1] / db[1] == val
                    push!(rationals, a // b)
                end
            end
        end
    end
    return denominator(prod(rationals))
end

p033 = Problems.Problem(p033solution)

Problems.benchmark(p033)

ext: jl

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
using TensorKit, TensorKitManifolds
using Test

spaces = (ℂ^4, ℤ₂Space(2,2), U₁Space(0=>2,1=>1,-1=>1), SU₂Space(0=>2,1/2=>1))
const ϵ = 1e-7
const α = 0.75

@testset "Grassmann with space $V" for V in spaces
    for T in (Float64,)
        W, = leftorth(TensorMap(randn, T, V*V*V, V*V); alg = Polar())
        X = TensorMap(randn, T, space(W))
        Y = TensorMap(randn, T, space(W))
        Δ = @inferred Grassmann.project(X, W)
        Θ = Grassmann.project(Y, W)
        γ = randn(T)
        Ξ = -Δ + γ*Θ
        @test norm(W'*Δ[]) <= sqrt(eps(real(T)))*dim(domain(W))
        @test norm(W'*Θ[]) <= sqrt(eps(real(T)))*dim(domain(W))
        @test norm(W'*Ξ[]) <= sqrt(eps(real(T)))*dim(domain(W))
        @test norm(zero(W)) == 0
        @test (@inferred Grassmann.inner(W, Δ, Θ)) ≈ real(dot(Δ[], Θ[]))
        @test Grassmann.inner(W, Δ, Θ) ≈ real(dot(X, Θ[]))
        @test Grassmann.inner(W, Δ, Θ) ≈ real(dot(Δ[],Y))
        @test Grassmann.inner(W, Δ, Δ) ≈ norm(Δ[])^2

        W2, = @inferred Grassmann.retract(W, Δ, ϵ)
        @test W2 ≈ W + ϵ*Δ[]
        W2, Δ2′ = Grassmann.retract(W, Δ, α)
        @test norm(W2'*Δ2′[]) <= sqrt(eps(real(T)))*dim(domain(W))
        @test Δ2′[] ≈ (first(Grassmann.retract(W, Δ, α + ϵ/2)) -
                        first(Grassmann.retract(W, Δ, α - ϵ/2)))/(ϵ)
        Δ2 = @inferred Grassmann.transport(Δ, W, Δ, α, W2)
        Θ2 = Grassmann.transport(Θ, W, Δ, α, W2)
        Ξ2 = Grassmann.transport(Ξ, W, Δ, α, W2)
        @test Δ2[] ≈ Δ2′[]
        @test norm(W2'*Δ2[]) <= sqrt(eps(real(T)))*dim(domain(W))
        @test norm(W2'*Θ2[]) <= sqrt(eps(real(T)))*dim(domain(W))
        @test norm(W2'*Ξ2[]) <= sqrt(eps(real(T)))*dim(domain(W))
        @test Ξ2[] ≈ -Δ2[] + γ * Θ2[]
        @test Grassmann.inner(W2, Δ2, Θ2) ≈ Grassmann.inner(W, Δ, Θ)
        @test Grassmann.inner(W2, Ξ2, Θ2) ≈ Grassmann.inner(W, Ξ, Θ)

        Wend = TensorMap(randhaar, T, codomain(W), domain(W))
        Δ3, V = Grassmann.invretract(W, Wend)
        @test Wend ≈ retract(W, Δ3, 1)[1] * V
        U = Grassmann.relativegauge(W, Wend)
        V2 = Grassmann.invretract(W, Wend * U)[2]
        @test V2 ≈ one(V2)
    end
end

@testset "Stiefel with space $V" for V in spaces
    for T in (Float64, ComplexF64)
        W = TensorMap(randhaar, T, V*V*V, V*V)
        X = TensorMap(randn, T, space(W))
        Y = TensorMap(randn, T, space(W))
        Δ = @inferred Stiefel.project_euclidean(X, W)
        Θ = Stiefel.project_canonical(Y, W)
        γ = rand()
        Ξ = -Δ + γ*Θ
        @test norm(W'*Δ[] + Δ[]'*W) <= sqrt(eps(real(T)))*dim(domain(W))
        @test norm(W'*Θ[] + Θ[]'*W) <= sqrt(eps(real(T)))*dim(domain(W))
        @test norm(W'*Ξ[] + Ξ[]'*W) <= sqrt(eps(real(T)))*dim(domain(W))
        @test norm(zero(W)) == 0
        @test (@inferred Stiefel.inner_euclidean(W, Δ, Θ)) ≈ real(dot(Δ[], Θ[]))
        @test (@inferred Stiefel.inner_canonical(W, Δ, Θ)) ≈
                                                        real(dot(Δ[], Θ[] - W*(W'*Θ[])/2))
        @test Stiefel.inner_euclidean(W, Δ, Θ) ≈ real(dot(X, Θ[]))
        @test !(Stiefel.inner_euclidean(W, Δ, Θ) ≈ real(dot(Δ[],Y)))
        @test !(Stiefel.inner_canonical(W, Δ, Θ) ≈ real(dot(X, Θ[])))
        @test Stiefel.inner_canonical(W, Δ, Θ) ≈ real(dot(Δ[],Y))
        @test Stiefel.inner_euclidean(W, Δ, Δ) ≈ norm(Δ[])^2
        @test Stiefel.inner_canonical(W, Δ, Δ) ≈ (1//2)*norm(W'*Δ[])^2 + norm(Δ[]-W*(W'Δ[]))^2

        W2, = @inferred Stiefel.retract_exp(W, Δ, ϵ)
        @test W2 ≈ W + ϵ*Δ[]
        W2, Δ2′ = Stiefel.retract_exp(W, Δ, α)
        @test norm(W2'*Δ2′[] + Δ2′[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W))
        @test Δ2′[] ≈ (first(Stiefel.retract_exp(W, Δ, α+ϵ/2)) -
                        first(Stiefel.retract_exp(W, Δ, α-ϵ/2)))/(ϵ)
        Δ2 = @inferred Stiefel.transport_exp(Δ, W, Δ, α, W2)
        Θ2 = Stiefel.transport_exp(Θ, W, Δ, α, W2)
        Ξ2 = Stiefel.transport_exp(Ξ, W, Δ, α, W2)
        @test Δ2′[] ≈ Δ2[]
        @test norm(W2'*Δ2[] + Δ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W))
        @test norm(W2'*Θ2[] + Θ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W))
        @test norm(W2'*Ξ2[] + Ξ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W))
        @test Ξ2[] ≈ -Δ2[] + γ * Θ2[]
        @test Stiefel.inner_euclidean(W2, Δ2, Θ2) ≈ Stiefel.inner_euclidean(W, Δ, Θ)
        @test Stiefel.inner_euclidean(W2, Ξ2, Θ2) ≈ Stiefel.inner_euclidean(W, Ξ, Θ)
        @test Stiefel.inner_canonical(W2, Δ2, Θ2) ≈ Stiefel.inner_canonical(W, Δ, Θ)
        @test Stiefel.inner_canonical(W2, Ξ2, Θ2) ≈ Stiefel.inner_canonical(W, Ξ, Θ)

        W2, = @inferred Stiefel.retract_cayley(W, Δ, ϵ)
        @test W2 ≈ W + ϵ*Δ[]
        W2, Δ2′ = Stiefel.retract_cayley(W, Δ, α)
        @test norm(W2'*Δ2′[] + Δ2′[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W))
        @test Δ2′[] ≈ (first(Stiefel.retract_cayley(W, Δ, α+ϵ/2)) -
                        first(Stiefel.retract_cayley(W, Δ, α-ϵ/2)))/(ϵ)
        @test norm(Δ2′) <= norm(Δ)
        Δ2 = @inferred Stiefel.transport_cayley(Δ, W, Δ, α, W2)
        Θ2 = Stiefel.transport_cayley(Θ, W, Δ, α, W2)
        Ξ2 = Stiefel.transport_cayley(Ξ, W, Δ, α, W2)
        @test !(Δ2′[] ≈ Δ2[])
        @test norm(W2'*Δ2[] + Δ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W))
        @test norm(W2'*Θ2[] + Θ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W))
        @test norm(W2'*Ξ2[] + Ξ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W))
        @test Ξ2[] ≈ -Δ2[] + γ * Θ2[]
        @test Stiefel.inner_euclidean(W2, Δ2, Θ2) ≈ Stiefel.inner_euclidean(W, Δ, Θ)
        @test Stiefel.inner_euclidean(W2, Ξ2, Θ2) ≈ Stiefel.inner_euclidean(W, Ξ, Θ)
        @test Stiefel.inner_canonical(W2, Δ2, Θ2) ≈ Stiefel.inner_canonical(W, Δ, Θ)
        @test Stiefel.inner_canonical(W2, Ξ2, Θ2) ≈ Stiefel.inner_canonical(W, Ξ, Θ)

        W3 = projectisometric!(W + 1e-1 * TensorMap(rand, T, codomain(W), domain(W)))
        Δ3 = Stiefel.invretract(W, W3)
        @test W3 ≈ retract(W, Δ3, 1)[1]
    end
end

@testset "Unitary with space $V" for V in spaces
    for T in (Float64, ComplexF64)
        W, = leftorth(TensorMap(randn, T, V*V*V, V*V); alg = Polar())
        X = TensorMap(randn, T, space(W))
        Y = TensorMap(randn, T, space(W))
        Δ = @inferred Unitary.project(X, W)
        Θ = Unitary.project(Y, W)
        γ = randn()
        Ξ = -Δ + γ*Θ
        @test norm(W'*Δ[] + Δ[]'*W) <= sqrt(eps(real(T)))*dim(domain(W))
        @test norm(W'*Θ[] + Θ[]'*W) <= sqrt(eps(real(T)))*dim(domain(W))
        @test norm(W'*Ξ[] + Ξ[]'*W) <= sqrt(eps(real(T)))*dim(domain(W))
        @test norm(zero(W)) == 0
        @test (@inferred Unitary.inner(W, Δ, Θ)) ≈ real(dot(Δ[], Θ[]))
        @test Unitary.inner(W, Δ, Θ) ≈ real(dot(X, Θ[]))
        @test Unitary.inner(W, Δ, Θ) ≈ real(dot(Δ[],Y))
        @test Unitary.inner(W, Δ, Δ) ≈ norm(Δ[])^2

        W2, = @inferred Unitary.retract(W, Δ, ϵ)
        @test W2 ≈ W + ϵ*Δ[]
        W2, Δ2′ = Unitary.retract(W, Δ, α)
        @test norm(W2'*Δ2′[] + Δ2′[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W))
        @test Δ2′[] ≈ (first(Unitary.retract(W, Δ, α+ϵ/2)) -
                        first(Unitary.retract(W, Δ, α-ϵ/2)))/(ϵ)

        Δ2 = @inferred Unitary.transport_parallel(Δ, W, Δ, α, W2)
        Θ2 = Unitary.transport_parallel(Θ, W, Δ, α, W2)
        Ξ2 = Unitary.transport_parallel(Ξ, W, Δ, α, W2)
        @test Δ2′[] ≈ Δ2[]
        @test norm(W2'*Θ2[] + Θ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W))
        @test norm(W2'*Ξ2[] + Ξ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W))
        @test Ξ2[] ≈ -Δ2[] + γ * Θ2[]
        @test Unitary.inner(W2, Δ2, Θ2) ≈ Unitary.inner(W, Δ, Θ)
        @test Unitary.inner(W2, Ξ2, Θ2) ≈ Unitary.inner(W, Ξ, Θ)

        Δ2 = @inferred Unitary.transport_stiefel(Δ, W, Δ, α, W2)
        Θ2 = Unitary.transport_stiefel(Θ, W, Δ, α, W2)
        Ξ2 = Unitary.transport_stiefel(Ξ, W, Δ, α, W2)
        @test Δ2′[] ≈ Δ2[]
        @test norm(W2'*Δ2[] + Δ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W))
        @test norm(W2'*Θ2[] + Θ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W))
        @test norm(W2'*Ξ2[] + Ξ2[]'*W2) <= sqrt(eps(real(T)))*dim(domain(W))
        @test Ξ2[] ≈ -Δ2[] + γ * Θ2[]
        @test Unitary.inner(W2, Δ2, Θ2) ≈ Unitary.inner(W, Δ, Θ)
        @test Unitary.inner(W2, Ξ2, Θ2) ≈ Unitary.inner(W, Ξ, Θ)
    end
end


text:
# using Distributed
# addprocs(5)

# @everywhere using ES

# info     = Esinfo("", "")
# index = "201807"
# setting = ES.esindexsetting(info, index)[index]["settings"]["index"]
# shards = parse(Int, setting["number_of_shards"])

# # 1
# @distributed (vcat) for i in 0:shards-1
# 	query= @smi(slice= @query(id=i, max= shards) ) |>  df ->  merge(df , @query(size = 100) )
# 	ES.esfsearch(info, index, query)
# end

#2
@fulltext(match = @query(type = "test"))
@fulltext(match = @query(type = "test"), bool = @filter(has("openid")))
@query(size = 1000, @filter(has("openid")))
@query(size = 1000, query = @filter(has("openid")))

#3
# ES.esfsearch(info, index, @query(size=100) )

max_stars_repo_name: frankill/ES

max_issues_repo_name: frankill/ES

max_forks_repo_name: frankill/ES


text:
using ClangCompiler
using Test

const CC = ClangCompiler

@testset "Status" begin
    src = joinpath(@__DIR__, "..", "code", "main.cpp")
    args = get_compiler_args()
    push!(args, "-std=c++11")
    pushfirst!(args, "clang")
    push!(args, "-I$(@__DIR__)")

    instance = CC.CompilerInstance()
    CC.createDiagnostics(instance)
    diag = CC.getDiagnostics(instance)

    invok = CC.createFromCommandLine(src, args, diag)
    CC.setInvocation(instance, invok)
    CC.setTargetAndLangOpts(instance)

    CC.PrintStats(instance, CC.HeaderSearchOptions)
    CC.PrintStats(instance, CC.DiagnosticOptions)
    CC.PrintStats(instance, CC.FrontendOptions)
    CC.PrintStats(instance, CC.CodeGenOptions)
    CC.PrintStats(instance, CC.PreprocessorOptions)
    CC.PrintStats(instance, CC.TargetOptions)

    CC.createFileManager(instance)
    CC.createSourceManager(instance)

    CC.getFileEntry(instance, src)

    CC.createPreprocessor(instance)

    CC.PrintStats(instance, CC.HeaderSearch)

    CC.dispose(instance)
end


text:
function quote_symbols(ex)
    ex
end

function quote_symbols(ex::Symbol)
    QuoteNode(ex)
end

function quote_symbols(ex::Expr)
    to_quote = eachindex(ex.args)

    if ex.head == :call
        if ex.args[1] == :_
            length(ex.args) > 2 ? error("Unqote only a single argument. Right: `_(x)`, Wrong: `_(x, y)`.") : nothing
            return ex.args[2]
        end
        to_quote = to_quote[2:end]
    end
    for i in to_quote
        ex.args[i] = quote_symbols(ex.args[i])
    end
    ex
end

macro quote_symbols(ex)
    esc(quote_symbols(ex))
end


text:
"""
    surfstar(data::Array{<:Real,2}, target::Array{<:Integer,1}, m::Signed=-1,
             dist_func::Any=(e1, e2) -> sum(abs.(e1 .- e2), dims=2);
             f_type::String="continuous")::Array{Float64,1}

Compute feature weights using SURF* algorithm. The f_type argument specifies whether the features are continuous or discrete
and can either have the value of "continuous" or "discrete".

---
# Reference:
- Casey S. Greene, Daniel S. Himmelstein, Jeff Kiralis, and Jason H. Mo-
ore. The informative extremes: Using both nearest and farthest indivi-
duals can improve Relief algorithms in the domain of human genetics.
In Clara Pizzuti, Marylyn D. Ritchie, and Mario Giacobini, editors,
Evolutionary Computation, Machine Learning and Data Mining in Bi-
oinformatics, pages 182–193. Springer, 2010.
"""
function surfstar(data::Array{<:Real,2}, target::Array{<:Integer,1}, m::Signed=-1,
                  dist_func::Any=(e1, e2) -> sum(abs.(e1 .- e2), dims=2);
                  f_type::String="continuous")::Array{Float64,1}

    # Initialize feature weights vector.
    weights = zeros(Float64, 1, size(data, 2))

    # Compute vectors of maximum and minimum feature values.
    max_f_vals = maximum(data, dims=1)
    min_f_vals = minimum(data, dims=1)

    # Sample m samples without replacement. If m has signal value -1, use all samples.
    sample_idxs = StatsBase.sample(1:size(data, 1), if (m==-1) size(data,1) else m end, replace=false)
    if (m == -1) m = size(data, 1) end # If m has signal value -1, set m to total number of examples.

    # Compute pairwise distances between samples (vector form).
    dists = Array{Float64}(undef, Int64((size(data, 1)^2 - size(data, 1))/2 + 1))
    dists[1] = 0  # Set first value of distances vector to 0 - accessed when i == j in square form indices.

    # Construct pairwise distances vector using vectorized distance function.
    top_ptr = 2
    @inbounds for idx = 1:size(data,1)-1
        upper_lim = top_ptr + size(data, 1) - idx - 1
        dists[top_ptr:upper_lim] = dist_func(data[idx:idx, :], data[idx+1:end, :])
        top_ptr = upper_lim + 1
    end

    # Get mean distance between all samples.
    mean_dist = Statistics.mean(dists[2:end])

    # Go over sampled indices.
    @inbounds for idx = sample_idxs

        # Row and column indices for querying pairwise distance vector.
        row_idxs = repeat([idx - 1], size(data, 1))
        col_idxs = collect(0:size(data, 1)-1)

        # Get indices in distance vector (from square form indices).
        neigh_idx = square_to_vec(row_idxs, col_idxs, size(data, 1)) .+ 2

        # Query distances to neighbors and get masks for neighbors that fall within
        # hypersphere with radius mean_dist and those that fall outside hypersphere.
        neigh_mask_near = dists[neigh_idx[neigh_idx .!= 1]] .< mean_dist
        neigh_mask_far = .!neigh_mask_near
        insert!(neigh_mask_near, idx, 0)
        insert!(neigh_mask_far, idx, 0)

        # Get masks for near and far hits/misses.
        hit_neigh_mask_near = neigh_mask_near .& (target .== target[idx])
        miss_neigh_mask_near = neigh_mask_near .& (target .!= target[idx])
        hit_neigh_mask_far = neigh_mask_far .& (target .== target[idx])
        miss_neigh_mask_far = neigh_mask_far .& (target .!= target[idx])

        # Ignore samples that have no hits or misses in radius.
        if sum(hit_neigh_mask_near) == 0 || sum(miss_neigh_mask_near) == 0 ||
            sum(hit_neigh_mask_far) == 0 || sum(miss_neigh_mask_far) == 0
            continue
        end

        # Get class values of near and far misses.
        miss_classes_near = target[miss_neigh_mask_near]
        miss_classes_far = target[miss_neigh_mask_far]

        # Compute weights for near misses and compute weighting vector.
        weights_mult1 = Array{Float64}(undef, length(miss_classes_near))  # Allocate weights multiplier vector.
        cm = countmap(miss_classes_near)                                  # Count unique values.
        u = collect(keys(cm))
        c = collect(values(cm))
        neighbor_weights = c ./ length(miss_classes_near)  # Compute misses' weights.

        @inbounds for (w, val) = zip(neighbor_weights, u)  # Build multiplier vector.
            find_res = findall(miss_classes_near .== val)
            weights_mult1[find_res] .= w
        end

        # Compute weights for far misses and compute weighting vector.
        weights_mult2 = Array{Float64}(undef, length(miss_classes_far))  # Allocate weights multiplier vector.
        cm = countmap(miss_classes_far)                                  # Count unique values.
        u = collect(keys(cm))
        c = collect(values(cm))
        neighbor_weights = c ./ length(miss_classes_far)   # Compute misses' weights.

        @inbounds for (w, val) = zip(neighbor_weights, u)  # Build multiplier vector.
            find_res = findall(miss_classes_far .== val)
            weights_mult2[find_res] .= w
        end

        ### Weights Update ###

        if f_type == "continuous"
            # If features continuous.

            # Penalty term for near neighbors.
            penalty_near = sum(abs.(data[idx:idx, :] .- data[hit_neigh_mask_near, :]) ./ (max_f_vals .- min_f_vals .+ eps(Float64)), dims=1)

            # Reward term for near neighbors.
            reward_near = sum(weights_mult1 .* abs.(data[idx:idx, :] .- data[miss_neigh_mask_near, :]) ./ (max_f_vals .- min_f_vals .+ eps(Float64)), dims=1)

            # Weights values for near neighbors.
            weights_near = weights .- penalty_near ./ (size(data, 1)*sum(hit_neigh_mask_near) + eps(Float64)) .+
                reward_near ./ (size(data, 1)*sum(miss_neigh_mask_near) + eps(Float64))

            # Penalty term for far neighbors.
            penalty_far = sum(abs.(data[idx:idx, :] .- data[hit_neigh_mask_far, :]) ./ (max_f_vals .- min_f_vals .+ eps(Float64)), dims=1)

            # Reward term for far neighbors.
            reward_far = sum(weights_mult2 .* abs.(data[idx:idx, :] .- data[miss_neigh_mask_far, :]) ./ (max_f_vals .- min_f_vals .+ eps(Float64)), dims=1)

            # Weights values for far neighbors.
            weights_far = weights .- penalty_far ./ (size(data, 1)*sum(hit_neigh_mask_far) + eps(Float64)) .+
                reward_far ./ (size(data, 1)*sum(miss_neigh_mask_far) + eps(Float64))

            # Update feature weights.
            weights = weights_near - (weights_far - weights)

        elseif f_type == "discrete"
            # If features discrete.

            # Penalty term for near neighbors.
            penalty_near = sum(Int64.(data[idx:idx, :] .!= data[hit_neigh_mask_near, :]), dims=1)

            # Reward term for near neighbors.
            reward_near = sum(weights_mult1 .* Int64.(data[idx:idx, :] .!= data[miss_neigh_mask_near, :]), dims=1)

            # Weights values for near neighbors.
            weights_near = weights .- penalty_near ./ (size(data, 1)*sum(hit_neigh_mask_near) + eps(Float64)) .+
                reward_near ./ (size(data, 1)*sum(miss_neigh_mask_near) + eps(Float64))

            # Penalty term for far neighbors.
            penalty_far = sum(Int64.(data[idx:idx, :] .!= data[hit_neigh_mask_far, :]), dims=1)

            # Reward term for far neighbors.
            reward_far = sum(weights_mult2 .* Int64.(data[idx:idx, :] .!= data[miss_neigh_mask_far, :]), dims=1)

            # Weights values for far neighbors.
            weights_far = weights .- penalty_far ./ (size(data, 1)*sum(hit_neigh_mask_far) + eps(Float64)) .+
                reward_far ./ (size(data, 1)*sum(miss_neigh_mask_far) + eps(Float64))

            # Update feature weights.
            weights = weights_near - (weights_far - weights)

        else
            throw(DomainError(f_type, "f_type can only be equal to \"continuous\" or \"discrete\"."))
        end

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

    end

    # Return computed feature weights.
    return vec(weights)
end


text:
@testset "Range Utitilities" begin
    @testset "find_first_above_or_equal" begin
        import NonuniformResampling1D.find_first_above_or_equal
        x = 1.0:1.0:10.0
        @test find_first_above_or_equal(.5, x) == 1
        @test find_first_above_or_equal(0, x) == 1
        @test find_first_above_or_equal(-.5, x) == 1
        @test find_first_above_or_equal(1, x) == 1
        @test find_first_above_or_equal(1.1, x) == 2
        @test find_first_above_or_equal(5.4, x) == 6
        @test find_first_above_or_equal(9.9, x) == 10
        @test find_first_above_or_equal(10, x) == 10
        @test_throws Exception find_first_above_or_equal(10.1, x)
    end

    @testset "find_last_below_or_equal" begin
        import NonuniformResampling1D.find_last_below_or_equal
        x = 1.0:1.0:10.0
        @test_throws Exception find_last_below_or_equal(.9, x)
        @test find_last_below_or_equal(1, x) == 1
        @test find_last_below_or_equal(1.1, x) == 1
        @test find_last_below_or_equal(5.4, x) == 5
        @test find_last_below_or_equal(9.9, x) == 9
        @test find_last_below_or_equal(10, x) == 10
        @test find_last_below_or_equal(10.1, x) == 10
        @test find_last_below_or_equal(30, x) == 10
    end

    @testset "find_last_below" begin
        import NonuniformResampling1D.find_last_below
        x = 1.0:1.0:10.0
        @test_throws Exception find_last_below(1, x)
        @test find_last_below(1.1, x) == 1
        @test find_last_below(5.4, x) == 5
        @test find_last_below(9.9, x) == 9
        @test find_last_below(10, x) == 9
        @test find_last_below(10.1, x) == 10
        @test find_last_below(30, x) == 10
    end
end


text:
using LabelledArrays

export rossler_lorenz

@inline @inbounds function eom_rossler_lorenz(u, p, t)
    c_xy, a₁, a₂, a₃, b₁, b₂, b₃ = (p...,)
    x1, x2, x3, y1, y2, y3 = u[1], u[2], u[3], u[4], u[5], u[6]

    dx1 = -a₁*(x2 + x3)
    dx2 = a₁*(x1 + a₂*x2)
    dx3 = a₁*(a₂ + x3*(x1 - a₃))
    dy1 = b₁*(-y1 + y2)
    dy2 = b₂*y1 - y2 - y1*y3 + c_xy*(x2^2)
    dy3 = y1*y2 - b₃*y3

    return SVector{6}(dx1, dx2, dx3, dy1, dy2, dy3)
end

"""
    rossler_lorenz(;u₀ = rand(6), a₁ = -6, a₂ = 6, a₃ = 2.0,
        b₁ = 10, b₂ = 28, b₃ = 8/3, c_xy = 1) → ContinuousDynamicalSystem

Initialise a Rössler-Lorenz system consisting of two independent 3D subsystems:
one Rössler system and one Lorenz system. They are coupled such that the
second component (`x₂`) of the Rössler system unidirectionally forces the
second component (`y₂`) of the Lorenz system.

The parameter `c_xy` controls the coupling strength. The implementation here also
allows for tuning the parameters of each subsystem by introducing the constants
`a₁`, `a₂`, `a₃`, `b₁`, `b₂`, `b₃`. Default values for these parameters are
as in [1].

## Equations of motion

The dynamics is generated by the following vector field

```math
\\begin{aligned}
\\dot x_1 &= a_1(x_2 + x_3) \\\\
\\dot x_2 &= a_2(x_1 + 0.2x_2) \\\\
\\dot x_3 &= a_2(0.2 + x_3(x_1 - a_3)) \\\\
\\dot y_1 &= b_1(y_2 - y_1) \\\\
\\dot y_2 &= y_1(b_2 - y_3) - y_2 +c_{xy}(x_2)^2 \\\\
\\dot y_3 &= y_1 y_2 - b_3y_3
\\end{aligned}
```

with the coupling constant ``c_{xy} \\geq 0``.

## References

1. Krakovská, Anna, et al. "Comparison of six methods for the detection of causality in a
    bivariate time series." Physical Review E 97.4 (2018):042207.
    [https://journals.aps.org/pre/abstract/10.1103/PhysRevE.97.042207](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.97.042207)
"""
function rossler_lorenz(;u₀ = rand(6), a₁ = -6, a₂ = 6, a₃ = 2.0,
    b₁ = 10, b₂ = 28, b₃ = 8/3, c_xy = 1)

    p = @LArray [c_xy, a₁, a₂, a₃, b₁, b₂, b₃] (:c_xy, :a₁, :a₂, :a₃, :b₁, :b₂, :b₃)
    ContinuousDynamicalSystem(eom_rossler_lorenz, u₀, p)
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# `DGMulti` refers to both multiple DG types (polynomial/SBP, simplices/quads/hexes) as well as
# the use of multi-dimensional operators in the solver.
const DGMulti{NDIMS, ElemType, ApproxType, SurfaceIntegral, VolumeIntegral} =
  DG{<:RefElemData{NDIMS, ElemType, ApproxType}, Mortar, SurfaceIntegral, VolumeIntegral} where {Mortar}

# these are necessary for pretty printing
polydeg(dg::DGMulti) = dg.basis.N
Base.summary(io::IO, dg::DG) where {DG <: DGMulti} = print(io, "DGMulti(polydeg=$(polydeg(dg)))")
Base.real(rd::RefElemData{NDIMS, Elem, ApproxType, Nfaces, RealT}) where {NDIMS, Elem, ApproxType, Nfaces, RealT} = RealT

"""
    DGMulti(; polydeg::Integer,
              element_type::AbstractElemShape,
              approximation_type=Polynomial(),
              surface_flux=flux_central,
              surface_integral=SurfaceIntegralWeakForm(surface_flux),
              volume_integral=VolumeIntegralWeakForm(),
              RefElemData_kwargs...)

Create a discontinuous Galerkin method which uses
- approximations of polynomial degree `polydeg`
- element type `element_type` (`Tri()`, `Quad()`, `Tet()`, and `Hex()` currently supported)

Optional:
- `approximation_type` (default is `Polynomial()`; `SBP()` also supported for `Tri()`, `Quad()`,
  and `Hex()` element types).
- `RefElemData_kwargs` are additional keyword arguments for `RefElemData`, such as `quad_rule_vol`.
  For more info, see the [StartUpDG.jl docs](https://jlchan.github.io/StartUpDG.jl/dev/).
"""
function DGMulti(; polydeg::Integer,
                   element_type::AbstractElemShape,
                   approximation_type=Polynomial(),
                   surface_flux=flux_central,
                   surface_integral=SurfaceIntegralWeakForm(surface_flux),
                   volume_integral=VolumeIntegralWeakForm(),
                   kwargs...)
  rd = RefElemData(element_type, approximation_type, polydeg, kwargs...)
  return DG(rd, nothing #= mortar =#, surface_integral, volume_integral)
end

# Type aliases. The first parameter is `ApproxType` since it is more commonly used for dispatch.
const DGMultiWeakForm{ApproxType, ElemType} =
  DGMulti{NDIMS, ElemType, ApproxType, <:SurfaceIntegralWeakForm, <:VolumeIntegralWeakForm} where {NDIMS}

const DGMultiFluxDiff{ApproxType, ElemType} =
  DGMulti{NDIMS, ElemType, ApproxType, <:SurfaceIntegralWeakForm, <:VolumeIntegralFluxDifferencing} where {NDIMS}

# now that DGMulti is defined, we can define constructors for VertexMappedMesh which use dg::DGMulti
"""
    VertexMappedMesh(vertex_coordinates, EToV, dg::DGMulti;
                     is_on_boundary = nothing,
                     is_periodic::NTuple{NDIMS, Bool} = ntuple(_->false, NDIMS)) where {NDIMS, Tv}

Constructor which uses `dg::DGMulti` instead of `rd::RefElemData`.
"""
VertexMappedMesh(vertex_coordinates, EToV, dg::DGMulti; kwargs...) =
  VertexMappedMesh(vertex_coordinates, EToV, dg.basis; kwargs...)

"""
    VertexMappedMesh(triangulateIO, dg::DGMulti, boundary_dict::Dict{Symbol, Int})

Constructor which uses `dg::DGMulti` instead of `rd::RefElemData`.
"""
VertexMappedMesh(triangulateIO, dg::DGMulti, boundary_dict::Dict{Symbol, Int}) =
  VertexMappedMesh(triangulateIO, dg.basis, boundary_dict)

# Todo: simplices. Add traits for dispatch on affine/curved meshes here.

# Matrix type for lazy construction of physical differentiation matrices

# lazy linear combination of B = ∑_i coeffs[i] * A[i]
struct LazyMatrixLinearCombo{Tcoeffs, N, Tv, TA <: AbstractMatrix{Tv}} <: AbstractMatrix{Tv}
  matrices::NTuple{N, TA}
  coeffs::NTuple{N, Tcoeffs}
  function LazyMatrixLinearCombo(matrices, coeffs)
    @assert all(matrix -> size(matrix) == size(first(matrices)), matrices)
    new{typeof(first(coeffs)), length(matrices), eltype(first(matrices)), typeof(first(matrices))}(matrices, coeffs)
  end
end
Base.eltype(A::LazyMatrixLinearCombo) = eltype(first(A.matrices))
Base.IndexStyle(A::LazyMatrixLinearCombo) = IndexCartesian()
Base.size(A::LazyMatrixLinearCombo) = size(first(A.matrices))

@inline function Base.getindex(A::LazyMatrixLinearCombo{<:Real, N}, i, j) where {N}
  val = zero(eltype(A))
  for k in Base.OneTo(N)
    val = val + A.coeffs[k] * getindex(A.matrices[k], i, j)
  end
  return val
end

end # @muladd


text:
using HDF5Logging
using Test
using Logging

@testset "log consistency" begin
    file = tempname() * ".hdf5"
    logger = hdf5_logger(file)
    with_logger(logger) do
        @info "test info" a = 1 b = 2
        @warn "test warn" c = 9
    end
    @test length(logger) == 2
    msg1 = logger[1]
    @test msg1.level == Logging.Info
    @test sort(msg1.data, by = first) == ["a" => 1, "b" => 2]
    msg2 = logger[2]
    @test msg2.level == Logging.Warn
    @test msg2.data == ["c" => 9]
end


text:
import Documenter: deploydocs

deploydocs(
    repo = "github.com/bramtayl/Keys.jl.git",
)


text:
import Observables
import AbstractPlotting

export Whirl
export update

"""
    Represents a whirl.

fields: points, top, bottom, s3rotation, config, segments, color and observable.
"""
mutable struct Whirl <: Sprite
    points::Array{<:S²,1}
    top::S¹
    bottom::S¹
    s3rotation::S³
    config::Biquaternion
    segments::Int
    color::AbstractPlotting.RGBAf0
    observable::Tuple{Observables.Observable{Array{Float64,2}},
                      Observables.Observable{Array{Float64,2}},
                      Observables.Observable{Array{Float64,2}}}
end

"""
    Whirl(scene,
          points,
          [top, [bottom, [s3rotation, [config, [segments, [color, [transparency]]]]]]])

Construct a whirl with the given `scene`, `points` in the base space, `top`, `bottom`,
S³ rotation `s3rotation`, configuration `config`, the number of `segments`, `color` and
`transparency`.
"""
function Whirl(scene::AbstractPlotting.Scene,
               points::Array{<:S²,1};
               top::S¹ = U1(-pi),
               bottom::S¹ = U1(pi),
               s3rotation::S³ = Quaternion(1, 0, 0, 0),
               config::Biquaternion = Biquaternion(ℝ³(0, 0, 0)),
               segments::Int = 36,
               color::AbstractPlotting.RGBAf0 = AbstractPlotting.RGBAf0(1.0,
                                                                        0.2705,
                                                                        0.0,
                                                                        0.5),
               transparency::Bool = false)
    whirl = constructwhirl(points,
                           top = top,
                           bottom = bottom,
                           s3rotation = s3rotation,
                           config = config,
                           segments = segments)
    colorarray = fill(color, segments, segments)
    observable = buildsurface(scene, whirl, colorarray, transparency = transparency)
    Whirl(points, top, bottom, s3rotation, config, segments, color, observable)
end

"""
    update(whirl, points)

Update a Whirl by changing its observable with the given `whirl` and `points` in the base
space.
"""
function update(whirl::Whirl, points::Array{<:S²,1})
    @assert(length(points) == length(whirl.points),
            "The number of the given points must be equal to what it was before.")
    whirl.points = points
    value = constructwhirl(whirl.points,
                           top = whirl.top,
                           bottom = whirl.bottom,
                           s3rotation = whirl.s3rotation,
                           config = whirl.config,
                           segments = whirl.segments)
    updatesurface(value, whirl.observable)
end

"""
    update(whirl, top, bottom)

Update a Whirl by changing its observable with the given `whirl`, `top` and `bottom`.
"""
function update(whirl::Whirl, top::S¹, bottom::S¹)
    whirl.top = top
    whirl.bottom = bottom
    value = constructwhirl(whirl.points,
                           top = whirl.top,
                           bottom = whirl.bottom,
                           s3rotation = whirl.s3rotation,
                           config = whirl.config,
                           segments = whirl.segments)
    updatesurface(value, whirl.observable)
end

"""
    update(whirl, s3rotation)

Update a Whirl by changing its observable with the given `whirl` and S³ rotation
`s3rotation`.
"""
function update(whirl::Whirl, s3rotation::S³)
    whirl.s3rotation = s3rotation
    value = constructwhirl(whirl.points,
                           top = whirl.top,
                           bottom = whirl.bottom,
                           s3rotation = whirl.s3rotation,
                           config = whirl.config,
                           segments = whirl.segments)
    updatesurface(value, whirl.observable)
end

"""
    update(whirl, config)

Update a Whirl by changing its observable with the given `whirl` and configuration `config`.
"""
function update(whirl::Whirl, config::Biquaternion)
    whirl.config = config
    value = constructwhirl(whirl.points,
                           top = whirl.top,
                           bottom = whirl.bottom,
                           s3rotation = whirl.s3rotation,
                           config = whirl.config,
                           segments = whirl.segments)
    updatesurface(value, whirl.observable)
end


text:
# This file is part of JuliaFMI.
# Licensed under MIT: https://github.com/AnHeuermann/JuliaFMI/blob/master/LICENSE.txt

# This file contains the unit testing.

"""
Run all tests for simulating FMU's
"""

using JuliaFMI, Test

include("compare_CSV.jl")
include("tests.jl")
include("fmi-cross-check-tests.jl")

@testset "All Tests" begin
    @testset "Simulating simple FMUs" begin
        @testset "Without events" begin
            testHelloFMI20World()
            if Sys.islinux()
                testCauerLowPassAnalog()
            end
        end
        @testset "With events" begin
            testBouncingBall()
        end
    end

    runFMICrossTests()
end

max_stars_repo_name: AnHeuermann/JuliaFMU

max_issues_repo_name: AnHeuermann/JuliaFMU

max_forks_repo_name: jwscook/JuliaFMI


text:
using LambdaC
using Test
using SoftPredicates
function runtests()
  @test λc"(+ (* 1 2) (+ (+ 4 0) 1))" |> init_interpret == 7
  @test λc"(let x 3 (+ x (* x x)))" |> init_interpret == 12
  @test λc"(! false)" |> init_interpret == true
  @test λc"(& true false)" |> init_interpret == false
  @test λc"((λ x (+ x 1)) 12)" |> init_interpret == 13
  @test λc"(((λ x
            (λ y (+ x y))) 5) 2)" |> init_interpret == 7
  @test λc"(let unif (λ j (λ w (w j)))
           (let x (λ ωq ((unif 1) ωq))
           (let ω (λ i (if (== i 1) then 0.2 else 0.5))
             (x ω))))" |> init_interpret == 0.2

  @test λc"(let x (λ ω (unif 1 ω))
           (let y (λ ω (unif 2 ω))
           (let z (λ ω (>ₛ (x ω) (y ω)))
               ((cond x z)
               (λ i (case i (=> 1 0.123)
                            (=> 2 0.313))))" |> soft_interpret == (0.123, 0.123 >ₛ 0.313)

  @test λc"((λ x (+ x 1)) 12)" |> init_interpret == 13
  # -- make cases
  # -- make apply work
  # -- make "do" work
end

function runtests2()
  @test λc"(+ (* 1 2) (+ (+ 4 0) 1))" |> init_sinterpret == 7
  @test λc"(let x 3 (+ x (* x x)))" |> init_interpret == 12
  @test λc"(! false)" |> init_interpret == true
  @test λc"(& true false)" |> init_interpret == false
  @test λc"((λ x (+ x 1)) 12)" |> init_interpret == 13
  @test λc"(((λ x
            (λ y (+ x y))) 5) 2)" |> init_interpret == 7
  @test λc"(let unif (λ j (λ w (w j)))
           (let x (λ ωq ((unif 1) ωq))
           (let ω (λ i (if (== i 1) then 0.2 else 0.5))
             (x ω))))" |> init_interpret == 0.2

  @test λc"(let x (λ ω (unif 1 ω))
           (let y (λ ω (unif 2 ω))
           (let z (λ ω (>ₛ (x ω) (y ω)))
               ((cond x z)
               (λ i (case i (=> 1 0.123)
                            (=> 2 0.313))))" |> soft_interpret == (0.123, 0.123 >ₛ 0.313)

  @test λc"((λ x (+ x 1)) 12)" |> init_interpret == 13
  # -- make cases
  # -- make apply work
  # -- make "do" work
end

function f()
  λc"(let unif (λ j (λ w (w j)))
          (let x (λ ωq ((unif 1) ωq))
          (let ω (λ i (if (== i 1) then 0.2 else 0.5))
            (x ω))))" |> init_interpret
end

f()

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
fm1 = LMM(ds[:Batch],ds[:Yield])

@test size(fm1.A) == (2,2)
@test size(fm1.trms) == (2,)
@test size(fm1.R) == (2,2)
@test size(fm1.Λ) == (1,)
@test lowerbd(fm1) == zeros(1)
@test fm1[:θ] == ones(1)

fm1[:θ] = [0.713]
@test objective(fm1) ≈ 327.34216280955366

fit(fm1)
@test objective(fm1) ≈ 327.3270598811428

fm2 = LMM(ds2[:Batch],ds2[:Yield])
@test lowerbd(fm2) == zeros(1)
fit(fm2)
@test fm2[:θ] == zeros(1)

X = hcat(ones(size(slp,1)),convert(Array,slp[:Days]))
fm3 = LMM([VectorReMat(slp[:Subject],X')],X,slp[:Reaction])
@test lowerbd(fm3) == [0.,-Inf,0.]
fit(fm3)
@test ReTerms.objective(fm3) ≈ 1751.9393444663153
@test_approx_eq_eps fm3[:θ] [0.9292213074888169,0.01816838485113137,0.22264487095998978] 1.e-6

fm4 = LMM([ReMat(psts[s]) for s in [:Sample,:Batch]],psts[:Strength])
@test lowerbd(fm4) == zeros(2)
@test fm4[:θ] == ones(2)
fit(fm4)
@test ReTerms.objective(fm4) ≈ 247.99446586325791
#@test fm4[:θ] ≈ [3.526885897445589,1.3299228050484744]

fm5 = LMM([ReMat(pen[s]) for s in [:Plate,:Sample]],pen[:Diameter])
@test lowerbd(fm5) == zeros(2)
@test fm5[:θ] == ones(2)
fit(fm5)


text:
############################################################
# judiWavefield ##############################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: June 2017

export judiWavefield, judiWavefieldException, judiDFTwavefield, muteWavefield, dump_wavefield

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

mutable struct judiWavefield{vDT<:Number} <: joAbstractLinearOperator{vDT,vDT}
	name::String
	m::Integer
	n::Integer
	info::Info
    dt
	data
end

mutable struct judiWavefieldException <: Exception
	msg :: String
end

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

## outer constructors

function judiWavefield(info,dt::Real,data::Union{Array, PyCall.PyObject, String}; vDT::DataType=Float32)

	# length of vector
	m = info.n * sum(info.nt)
	n = 1
	dataCell = Array{Any}(undef, info.nsrc)
	for j=1:info.nsrc
		dataCell[j] = data
	end
	return judiWavefield{vDT}("judiWavefield",m,n,info,dt,dataCell)
end

function judiWavefield(info,dt::Real,data::Array{Any,1};vDT::DataType=Float32)
	# length of vector
	m = info.n * sum(info.nt)
	n = 1
	return judiWavefield{vDT}("judiWavefield",m,n,info,dt,data)
end

####################################################################
## overloaded Base functions

# conj(jo)
conj(A::judiWavefield{vDT}) where vDT =
	judiWavefield{vDT}("conj("*A.name*")",A.m,A.n,A.info,A.dt,A.data)

# transpose(jo)
transpose(A::judiWavefield{vDT}) where vDT =
	judiWavefield{vDT}(""*A.name*".'",A.n,A.m,A.info,A.dt,A.data)

# ctranspose(jo)
ctranspose(A::judiWavefield{vDT}) where vDT =
	judiWavefield{vDT}(""*A.name*"'",A.n,A.m,A.info,A.dt,A.data)

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

function vcat(a::judiWavefield{avDT},b::judiWavefield{bvDT}) where {avDT, bvDT}
	m = a.m + b.m
	n = 1
	nsrc = a.info.nsrc + b.info.nsrc
	data = Array{Any}(undef, nsrc)
	nt = Array{Any}(undef, nsrc)
	for j=1:a.info.nsrc
		data[j] = a.data[j]
		nt[j] = a.info.nt[j]
	end
	for j=a.info.nsrc+1:nsrc
		data[j] = b.data[j-a.info.nsrc]
		nt[j] = b.info.nt[j-a.info.nsrc]
	end
	info = Info(a.info.n,nsrc,nt)
	return judiWavefield(info,a.dt,data)
end

# add and subtract, mulitply and divide, norms, dot ...

# DFT operator for wavefields, acts along time dimension
function fft_wavefield(x_in,mode)
	nsrc = x_in.info.nsrc
	if mode==1
		x = judiWavefield(x_in.info,deepcopy(x_in.data); vDT=Complex{Float32})
		for i=1:nsrc
			x.data[i] = convert(Array{Complex{Float32}},x.data[i])
			nx = size(x.data[i],2)
			nz = size(x.data[i],3)
			for j=1:nx
				for k=1:nz
					x.data[i][:,j,k] = fft(x.data[i][:,j,k])
				end
			end
		end
	elseif mode==-1
		x = judiWavefield(x_in.info,deepcopy(x_in.data); vDT=Float32)
		for i=1:nsrc
			nx = size(x.data[i],2)
			nz = size(x.data[i],3)
			for j=1:nx
				for k=1:nz
					x.data[i][:,j,k] = real(ifft(x.data[i][:,j,k]))
				end
			end
			x.data[i] = convert(Array{Float32},real(x.data[i]))
		end
	end
	return x
end

# Sampling mask to extract wavefields from full vector
subsample(u::judiWavefield,srcnum) = judiWavefield(u.info,u.data[srcnum];vDT=eltype(u))

function muteWavefield(u_in::judiWavefield,ts_keep)
	u = deepcopy(u_in)
	for j=1:u.info.nsrc
		idx = ones(size(u.data[j],1));
		idx[ts_keep] = 0
		zero_idx = findall(idx)
		u.data[j][zero_idx,:,:] *= 0.f0
	end
	return u
end

# norm
function norm(a::judiWavefield{avDT}, p::Real=2) where avDT
	np = load_numpy()
    x = 0.f0
    for j=1:a.info.nsrc
        if typeof(a.data[j]) == String
            x += a.dt * sum(np[:abs](np[:load](a.data[j])).^p)
        else
            x += a.dt * sum(np[:abs](a.data[j]["data"]).^p)
        end
    end
    return x^(1.f0/p)
end

# Save wavefield to disk
function dump_wavefield(u::PyObject)
    name = join(["wavefield_", randstring(8), ".dat"])
    u["data"][:dump](name)
    return name
end

function dump_wavefield(u::judiWavefield)
    name = Array{Any}(undef, u.info.nsrc)
    for j=1:u.info.nsrc
        name[j] = join(["wavefield_", randstring(8), ".dat"])
        u.data[j]["data"][:dump](name[j])
    end
    return name
end


text:
using MetabolicFluxAnalysis

x0  = 0.1*rand(2,1000)
tic()
MetabolicFluxAnalysis.hrstep!( [eye(2);-eye(2)] , [ones(2);ones(2)] , 0.8*sqrt(2) , x0 , 2000)
toc()

x0  = 0.1*rand(2,1000)
tic()
MetabolicFluxAnalysis.hrstep_cchr!( [eye(2);-eye(2)] , [ones(2);ones(2)] , 0.8*sqrt(2) , x0 , 2000)
toc()

function benchmark_pencilbox( d::Int , n_samples::Int , n_steps::Int )
    A = 1.0 * [eye(d);-eye(d)]
    b = 1.0 * [1:d;zeros(d)]

    x0  = 0.1*rand(d,n_samples)
    r   = 10.0

    tic()
    MetabolicFluxAnalysis.hrstep_cchr!( A , b , r , x0 , n_steps)
    t1 = toc()
    tic()
    MetabolicFluxAnalysis.hrstep!(      A , b , r , x0 , n_steps)
    t2 = toc()

    return (t1 , t2)
end

benchmark_pencilbox( 100,1000,1000 )

#using Plots
#gr()
#plot()
#scatter!(x0[1,:],x0[2,:])
#plot!()
(lambda_min,lambda_max,x_intersection_min,x_intersection_max) = MetabolicFluxAnalysis.xrayshoot_centerball(1.2,0.0*rand(2,200),rand(2,200))
#using Plots
#gr()
#plot()
#scatter!(x_intersection_max[1,:],x_intersection_max[2,:])
#plot!()

using MetabolicFluxAnalysis
pa1 = [eye(2);-eye(2)]
pb1 = [1.;2;3;4] # size  should be (1+3)*(2+4) = 24
MetabolicFluxAnalysis.vol_mmc_reversed(   pa1 , pb1 , 500,200,10,2  )
MetabolicFluxAnalysis.vol_mmc_reversed(   pa1 * qr(rand(2,2))[1] , pb1 , 500,200,10,2  )
MetabolicFluxAnalysis.vol_mmc_reversed(   pa1 * qr(rand(2,2))[1] , pb1 , 500,200,10,2  )


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#    Tests of reported issues.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Issue #4 - Header with escape sequences" begin

    expected = """
┌───┬───────┬─────┬───┐
│ 1 │   2\\n │   3 │ 4 │
├───┼───────┼─────┼───┤
│ 1 │ false │ 1.0 │ 1 │
│ 2 │  true │ 2.0 │ 2 │
│ 3 │ false │ 3.0 │ 3 │
│ 4 │  true │ 4.0 │ 4 │
│ 5 │ false │ 5.0 │ 5 │
│ 6 │  true │ 6.0 │ 6 │
└───┴───────┴─────┴───┘
"""

    result = pretty_table(String, data, ["1" "2\n" "3" "4"])
    @test result == expected
end

@testset "Issue #9 - Printing table with row number column and sub-headers" begin
    expected = """
┌─────┬───┬───────┬─────┬────┐
│ Row │ 1 │     2 │   3 │  4 │
│     │ 5 │     6 │   7 │  8 │
│     │ 9 │    10 │  11 │ 12 │
├─────┼───┼───────┼─────┼────┤
│   1 │ 1 │ false │ 1.0 │  1 │
│   2 │ 2 │  true │ 2.0 │  2 │
│   3 │ 3 │ false │ 3.0 │  3 │
│   4 │ 4 │  true │ 4.0 │  4 │
│   5 │ 5 │ false │ 5.0 │  5 │
│   6 │ 6 │  true │ 6.0 │  6 │
└─────┴───┴───────┴─────┴────┘
"""
    result = pretty_table(String, data, [1  2  3  4;
                                         5  6  7  8;
                                         9 10 11 12];
                          show_row_number = true)
    @test result == expected
end

@testset "Issue #10 - Size of the row number column without header" begin
    v = 0:1:6

    expected = """
┌───┬───┐
│ 1 │ 0 │
│ 2 │ 1 │
│ 3 │ 2 │
│ 4 │ 3 │
│ 5 │ 4 │
│ 6 │ 5 │
│ 7 │ 6 │
└───┴───┘
"""

    result = pretty_table(String, v; noheader = true, show_row_number = true)
    @test result == expected
end

@testset "Issue #19 - ft_printf with cells that are not numbers" begin
    matrix = [1 1; 2 2; 3 3; "teste" "teste"; 4 4; 5 5; true true; :s :s]

    expected = """
┌────────────┬────────┐
│     Col. 1 │ Col. 2 │
├────────────┼────────┤
│       1.00 │      1 │
│       2.00 │      2 │
│       3.00 │      3 │
│      teste │  teste │
│       4.00 │      4 │
│       5.00 │      5 │
│       1.00 │   true │
│          s │      s │
└────────────┴────────┘
"""

    result = pretty_table(String, matrix; formatters = ft_printf("%10.2f",[1]))
    @test result == expected

end

@testset "Issue #22 - Strings with quotes" begin
    # Without linebreaks.
    matrix = [1 "teste\"teste"
              2 "teste\"\"teste"]

    expected = """
┌────────┬──────────────┐
│ Col. 1 │       Col. 2 │
├────────┼──────────────┤
│      1 │  teste"teste │
│      2 │ teste""teste │
└────────┴──────────────┘
"""

    result = pretty_table(String, matrix)
    @test result == expected

    # With linebreaks.
    matrix = [ 1 """
    function str(str = "one string")
        return str
    end"""
    2 """
    function str(str = "")
        if isempty(str)
            return "one string"
        else
            return str
        end
    end"""]

    expected = """
┌────────┬──────────────────────────────────┐
│ Col. 1 │ Col. 2                           │
├────────┼──────────────────────────────────┤
│ 1      │ function str(str = "one string") │
│        │     return str                   │
│        │ end                              │
│ 2      │ function str(str = "")           │
│        │     if isempty(str)              │
│        │         return "one string"      │
│        │     else                         │
│        │         return str               │
│        │     end                          │
│        │ end                              │
└────────┴──────────────────────────────────┘
"""

    result = pretty_table(String, matrix; alignment = :l, linebreaks = true)

    @test result == expected
end

@testset "Issue #24 - Tables compatibility" begin
    # Named tuple of vectors (satisfies Tables interface).
    ctable = (x = Int64[1, 2, 3, 4], y = ["a", "b", "c", "d"])

    cresult = sprint(pretty_table, ctable)

    expected = """
┌───────┬────────┐
│     x │      y │
│ Int64 │ String │
├───────┼────────┤
│     1 │      a │
│     2 │      b │
│     3 │      c │
│     4 │      d │
└───────┴────────┘
"""

  @test cresult == expected
end

@testset "Issue #28 - Tables.jl API must have priority when printing" begin
    # A NamedTuple is compliant with Tables.jl API.
    table = (x = Int64(1):Int64(3),
             y = 'a':'c',
             z = ["String 1";"String 2";"String 3"]);

    # Thus, the following 3 calls must provide the same results.
    result_1 = pretty_table(String, table)
    result_2 = pretty_table(String, Tables.rowtable(table))
    result_3 = pretty_table(String, Tables.columns(table))

    expected = """
┌───────┬──────┬──────────┐
│     x │    y │        z │
│ Int64 │ Char │   String │
├───────┼──────┼──────────┤
│     1 │    a │ String 1 │
│     2 │    b │ String 2 │
│     3 │    c │ String 3 │
└───────┴──────┴──────────┘
"""

    @test result_1 == result_2 == result_3 == expected
end


text:
module Utils

    function explode_union(T)
        @nospecialize T

        # unpack unionall
        S = T
        vars = []
        while S isa UnionAll
            pushfirst!(vars, S.var)
            S = S.body
        end

        if S isa Union
            Us = Any[explode_union(S.a)..., explode_union(S.b)...]
            Any[foldl((body, var) -> UnionAll(var, body), vars, init=U) for U in Us]
        elseif S == Union{}
            Any[]
        else
            Any[T]
        end
    end

    """
        pointer_from_obj(x)

    Returns `(p, c)` where `Base.pointer_from_objref(p) === x`.

    The pointer remains valid provided the object `c` is not garbage collected.
    """
    function pointer_from_obj(o::T) where {T}
        if T.mutable
            c = o
            p = Base.pointer_from_objref(o)
        else
            c = Ref{Any}(o)
            p = unsafe_load(Ptr{Ptr{Cvoid}}(Base.pointer_from_objref(c)))
        end
        p, c
    end

    """
        ExtraNewline(x)

    An object that displays the same as `x` but with an extra newline in text/plain.
    """
    struct ExtraNewline{T}
        value :: T
    end
    Base.show(io::IO, m::MIME, x::ExtraNewline) = show(io, m, x.value)
    Base.show(io::IO, m::MIME"text/csv", x::ExtraNewline) = show(io, m, x.value)
    Base.show(io::IO, m::MIME"text/tab-separated-values", x::ExtraNewline) = show(io, m, x.value)
    Base.show(io::IO, m::MIME"text/plain", x::ExtraNewline) = (show(io, m, x.value); println(io))
    Base.showable(m::MIME, x::ExtraNewline) = showable(m, x.value)

    const ALL_MIMES = [
        "text/plain",
        "text/html",
        "text/markdown",
        "text/json",
        "text/latex",
        "text/xml",
        "text/csv",
        "application/javascript",
        "application/pdf",
        "application/ogg",
        "image/jpeg",
        "image/png",
        "image/svg+xml",
        "image/gif",
        "image/webp",
        "image/tiff",
        "image/bmp",
        "audio/aac",
        "audio/mpeg",
        "audio/ogg",
        "audio/opus",
        "audio/webm",
        "audio/wav",
        "audio/midi",
        "audio/x-midi",
        "video/mpeg",
        "video/ogg",
        "video/webm",
    ]

    function mimes_for(x)
        @nospecialize x
        # default mimes we always try
        mimes = copy(ALL_MIMES)
        # look for mimes on show methods for this type
        for meth in methods(show, Tuple{IO, MIME, typeof(x)}).ms
            mimetype = meth.sig.parameters[3]
            mimetype isa DataType || continue
            mime = string(mimetype.parameters[1])
            push!(mimes, mime)
        end
        return mimes
    end

    @generated _typeintersect(::Type{T1}, ::Type{T2}) where {T1,T2} = typeintersect(T1, T2)

    @generated _type_ub(::Type{T}) where {T} = begin
        S = T
        while S isa UnionAll
            S = S{S.var.ub}
        end
        S
    end

    @generated _type_lb(::Type{T}) where {T} = begin
        R = T
        while R isa UnionAll
            R = R.body
        end
        if R isa DataType
            S = T
            while S isa UnionAll
                S = S{S.var in R.parameters ? S.var.lb : S.var.ub}
            end
            S
        else
            _type_ub(T)
        end
    end

    @generated _promote_type_bounded(::Type{S}, ::Type{T}, ::Type{B}) where {S,T,B} = begin
        S <: B || error("require S <: B")
        T <: B || error("require T <: B")
        if B isa Union
            return Union{_promote_type_bounded(typeintersect(S, B.a), typeintersect(T, B.a), B.a), _promote_type_bounded(typeintersect(S, B.b), typeintersect(T, B.b), B.b)}
        else
            R = promote_type(S, T)
            if R <: B
                return R
            else
                R = typeintersect(typejoin(S, T), B)
                if R <: B
                    return R
                else
                    return B
                end
            end
        end
    end

    @generated _promote_type_bounded(::Type{T1}, ::Type{T2}, ::Type{T3}, ::Type{B}) where {T1,T2,T3,B} =
        _promote_type_bounded(_promote_type_bounded(T1, T2, B), T3, B)

    # TODO: what is the best way?
    ismutablearray(x::Array) = true
    ismutablearray(x::AbstractArray) = begin
        p = parent(x)
        p === x ? false : ismutablearray(p)
    end

    islittleendian() = Base.ENDIAN_BOM == 0x04030201 ? true : Base.ENDIAN_BOM == 0x01020304 ? false : error()

    isflagset(flags, mask) = (flags & mask) == mask

    size_to_fstrides(elsz::Integer, sz::Tuple{Vararg{Integer}}) =
        isempty(sz) ? () : (elsz, size_to_fstrides(elsz * sz[1], sz[2:end])...)

    size_to_cstrides(elsz::Integer, sz::Tuple{Vararg{Integer}}) =
        isempty(sz) ? () : (size_to_cstrides(elsz * sz[end], sz[1:end-1])..., elsz)

    isfcontiguous(o::AbstractArray) = strides(o) == size_to_fstrides(1, size(o)...)
    isccontiguous(o::AbstractArray) = strides(o) == size_to_cstrides(1, size(o)...)

    struct StaticString{T,N} <: AbstractString
        codeunits :: NTuple{N,T}
    end

    function Base.convert(::Type{String}, x::StaticString)
        cs = collect(x.codeunits)
        i = findfirst(==(0), cs)
        transcode(String, i===nothing ? cs : cs[1:i-1])
    end

    Base.String(x::StaticString) = convert(String, x)

    function Base.convert(::Type{StaticString{T,N}}, x::String) where {T,N}
        cs = transcode(T, x)
        length(cs) > N && throw(InexactError(:convert, StaticString{T,N}, x))
        while length(cs) < N
            push!(cs, 0)
        end
        StaticString{T,N}(NTuple{N,T}(cs))
    end

    Base.convert(::Type{T}, x::AbstractString) where {T<:StaticString} = convert(T, convert(String, x))

    (::Type{T})(x::AbstractString) where {T<:StaticString} = convert(T, x)

    function Base.iterate(x::StaticString, st=nothing)
        if st === nothing
            s = String(x)
            z = iterate(s)
        else
            s, st0 = st
            z = iterate(s, st0)
        end
        if z === nothing
            nothing
        else
            c, newst0 = z
            (c, (s, newst0))
        end
    end
end


text:
# generalspace.jl
#
# Defines the immutible GeneralSpace for a general finite dimensional vector space over an arbitrary
# field, characterized by its field, its dimension and whether or not it is the dual space and/or the
# complex conjugate space. Tensors with GeneralSpace as index spaces make a distinction between covariant
# and contravariant indices, and covariant and contravariant dotted/barred indices.

# GeneralSpace:
#---------------
immutable GeneralSpace{F} <: ElementarySpace{F}
    field::Type{F}
    d::Int
    dual::Bool
    conj::Bool
    GeneralSpace(d::Int,dual::Bool=false,conj::Bool=false) = (d>0 ? new(F, d, dual, !(F<:Real) && conj) : throw(ArgumentError("Dimension of a vector space should be bigger than zero")))
end
GeneralSpace{F}(::Type{F},d::Int,dual::Bool=false,conj::Bool=false) = GeneralSpace{F}(d,dual,conj)
GeneralSpace(d::Int,dual::Bool=false,conj::Bool=false) = GeneralSpace{Real}(d,dual,conj)

# Corresponding methods:
dim(V::GeneralSpace) = V.d
dual(V::GeneralSpace) = GeneralSpace(V.field,V.d, !V.dual, V.conj)
Base.conj(V::GeneralSpace) = GeneralSpace(V.field,V.d, V.dual, !V.conj)
cnumber{F}(V::GeneralSpace{F}) = GeneralSpace(F,1,V.dual,V.conj)
cnumber{F}(::Type{GeneralSpace{F}}) = GeneralSpace(F,1,false,false)
iscnumber(V::ElementarySpace) = dim(V)==1 && V.field <: Number

# Show methods
function Base.show(io::IO, V::GeneralSpace)
    if V.conj
        print(io,"conj(")
    end
    print(io,"$(V.field)($(V.d))")
    if V.dual
        print(io,"*")
    end
    if V.conj
        print(io,")")
    end
end


text:
export optimal
using Optim
using Statistics

# sum of absolute differences (element-wise)
function Δ⃯(𝒄₁, 𝒄₂)
  return sum(abs.(𝒄₁- 𝒄₂))
end

# sum of distances outside bounding box (element-wise)
function Δ͍(𝒄, 𝒄₋, 𝒄₊)
  return sum(max.(𝒄₋-𝒄, 0.0, 𝒄-𝒄₊))
end

# collection 𝒄 with last element dropped
function drop(𝒄)
  return 𝒄[1:end-1]
end

# fill collection with length of 𝒄 with average of 𝒄
function avg(𝒄)
  return fill(mean(𝒄), length(𝒄))
end

# map frequencies and diameters to locations, but drop the flute-length value
⇴ = drop ∘ mapflute
# push a value to front of vector, while dropping last value
⥆ = drop ∘ vcat

# generate error function with scoped constants
function mkerrfn(𝒇, 𝒅₊, 𝒑₋, 𝒑₊)
  𝒍⃯ = 𝒇⇴ 𝒅₊ # positions of holes with maximum allowed diameters
  function errfn(𝒅)
    # locations
    𝒍 = 𝒇⇴ 𝒅      # proposed diameters
    𝒍̲ = 𝒇⇴ avg(𝒅) # average of proposed
    𝒍⃮ = 0.0⥆ 𝒍    # shift proposed locations left (push 0.0)
    𝒍⃭ = 𝒍⃮+𝒑₋      # closest desired position to previous hole
    𝒍⃬ = 𝒍⃮+𝒑₊      # furthest desired position from previous hole
    # error terms
    𝑒̲ = Δ⃯(𝒍, 𝒍̲)       # sum of distances from average
    𝑒⃯ = Δ⃯(𝒍, 𝒍⃯)       # sum of distances from max positions
    𝑒͍ = Δ͍(𝒍, 𝒍⃭, 𝒍⃬)^2  # squared sum of distances outside desired from previous holes
    # sum and weigh error terms
    𝑒 = 𝑒̲ + 2𝑒⃯ + 3𝑒͍
    return 𝑒
  end
  return errfn
end

function optimal(𝒇, 𝒅₋, 𝒅₊, 𝒑₋, 𝒑₊; trace=false)
  # minimize error function
  errfn = mkerrfn(𝒇, 𝒅₊, 𝒑₋, 𝒑₊)
  # box-constrained, initial parameters (bad guess)
  𝒅₁ = [8.8, 8.8, 6.6, 8.8, 9.9, 5.5]
  # simulated annealing (round 1, fast cooldown)
  algo = SAMIN(rt=0.9)
  options = Optim.Options(iterations=Int(1e5), show_trace=trace, show_every=Int(1e4))
  𝒅₂ = Optim.minimizer(optimize(errfn, 𝒅₋, 𝒅₊, 𝒅₁, algo, options))
  # simulated annealing (round 2, slow cooldown)
  algo = SAMIN(rt=0.99)
  options = Optim.Options(iterations=Int(1e6), show_trace=trace, show_every=Int(1e4))
  result = optimize(errfn, 𝒅₋, 𝒅₊, 𝒅₂, algo, options)
  𝒅ₕ = Optim.minimizer(result)
  𝑒 = Optim.minimum(result)
  # proposed diameters
  return 𝒅ₕ, 𝑒
end

max_stars_repo_name: isaacwhanson/scad-flute

max_issues_repo_name: isaacwhanson/scad-flute

max_forks_repo_name: isaacwhanson/scad-flute


text:
# Parallel instance of fwi_objective function # Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

"""
    fwi_objective(model, source, dobs; options=Options())

Evaluate the full-waveform-inversion (reduced state) objective function. Returns a tuple with function value and vectorized \\
gradient. `model` is a `Model` structure with the current velocity model and `source` and `dobs` are the wavelets and \\
observed data of type `judiVector`.

Example
=======

    function_value, gradient = fwi_objective(model, source, dobs)

"""
function fwi_objective(model::Model, source::judiVector, dObs::judiVector; options=Options(), frequencies=[])
# fwi_objective function for multiple sources. The function distributes the sources and the input data amongst the available workers.

    # Process shots from source channel asynchronously
    fwi_objective = retry(TimeModeling.fwi_objective)
    results = Array{Any}(undef, dObs.nsrc)
    @sync begin
        for j=1:dObs.nsrc
            results[j] = @spawn fwi_objective(model, source[j], dObs[j], j; options=options, frequencies=frequencies)
        end
    end

    # Collect and reduce gradients
    gradient = zeros(Float32,prod(model.n)+1)
    for j=1:dObs.nsrc
        gradient += fetch(results[j]); results[j] = []
    end

    # first value corresponds to function value, the rest to the gradient
    return gradient[1], gradient[2:end]
end


text:
using ApproxFun, SingularIntegralEquations, Test

@testset "Fractal" begin
    d = Segment(0.,1.)
    c = Circle()

    for set in (:cantor,:thincantor,:thinnercantor,:thinnestcantor)
        @eval begin
            d = Segment(0.,1.)
            c = Circle()
            @test $set(d,0) == d
            @test $set(d,0,4) == d
            @test $set(d,1) == Segment(0.0,1/3) ∪ Segment(2/3,1.0)
            @test $set(d,1,4) == Segment(0.0,1/4) ∪ Segment(3/4,1.0)
            @test $set(c,0) == c
            @test $set(c,1) == Arc(0.0,1.0,(-2.6179938779914944,-0.523598775598299)) ∪
                                Arc(0.0,1.0,(0.5235987755982987,2.6179938779914944))
        end
    end

    @test cantor(d,2) == Segment(0.0,1/9)∪Segment(2/9,1/3)∪Segment(2/3,7/9)∪Segment(8/9,1.0)
    @test cantor(d,2,4) == Segment(0.0,1/16)∪Segment(3/16,1/4)∪Segment(3/4,13/16)∪Segment(15/16,1.0)

    @test thincantor(d,2) == Segment(0.0,1/12)∪Segment(1/4,1/3)∪Segment(2/3,3/4)∪Segment(11/12,1.0)
    @test thincantor(d,2,4) == Segment(0.0,1/20)∪Segment(1/5,1/4)∪Segment(3/4,4/5)∪Segment(19/20,1.0)

    @test thinnercantor(d,2) == Segment(0.0,1/27)∪Segment(8/27,1/3)∪Segment(2/3,19/27)∪Segment(26/27,1.0)
    @test thinnercantor(d,2,4) == Segment(0.0,1/64)∪Segment(15/64,1/4)∪Segment(3/4,49/64)∪Segment(63/64,1.0)
    @test thinnercantor(d,3,4) == Segment(0.0,1/4096)∪Segment(63/4096,1/64)∪Segment(15/64,961/4096)∪Segment(1023/4096,1/4)∪Segment(3/4,3073/4096)∪Segment(3135/4096,49/64)∪Segment(63/64,4033/4096)∪Segment(4095/4096,1/1)

    @test thinnestcantor(d,2) == Segment(0.0,1/27)∪Segment(8/27,1/3)∪Segment(2/3,19/27)∪Segment(26/27,1.0)
    @test thinnestcantor(d,2,4) == Segment(0.0,1/64)∪Segment(15/64,1/4)∪Segment(3/4,49/64)∪Segment(63/64,1.0)
    @test thinnestcantor(d,3,4) == Segment(0.0,1/16384)∪Segment(255/16384,1/64)∪Segment(15/64,3841/16384)∪Segment(4095/16384,1/4)∪Segment(3/4,12289/16384)∪Segment(12543/16384,49/64)∪Segment(63/64,16129/16384)∪Segment(16383/16384,1/1)

    @test smithvolterracantor(d,0) == d
    @test smithvolterracantor(d,1) == Segment(0.0,3/8)∪Segment(5/8,1.0)
    @test smithvolterracantor(d,2) == Segment(0.0,5/32)∪Segment(7/32,3/8)∪Segment(5/8,25/32)∪Segment(27/32,1.0)
end


text:
@testset "Test Contructor I" begin
    mdl1 = LogisticFunction(2)
    mdl2 = HyperbolicTangent(3)
    mdl3 = LogisticFunction(4)
    mdl4 = Identity(2)
    mdl = FeedforwardNetwork(3, [mdl1, mdl2, mdl3, mdl4])

    @test mdl.ninputs == 3
    @test mdl.noutputs == 2
    @test mdl.mdls[1] == Linear(3, 2)
    @test mdl.mdls[2] == Bias(2)
    @test mdl.mdls[3] == LogisticFunction(2)
    @test mdl.mdls[4] == Linear(2, 3)
    @test mdl.mdls[5] == Bias(3)
    @test mdl.mdls[6] == HyperbolicTangent(3)
    @test mdl.mdls[7] == Linear(3, 4)
    @test mdl.mdls[8] == Bias(4)
    @test mdl.mdls[9] == LogisticFunction(4)
    @test mdl.mdls[10] == Linear(4, 2)
    @test mdl.mdls[11] == Bias(2)
    @test mdl.mdls[12] == Identity(2)
end

@testset "Test Contructor II" begin
    mdl = FeedforwardNetwork(3, 2, [2, 3, 4])

    @test mdl.ninputs == 3
    @test mdl.noutputs == 2
    @test mdl.mdls[1] == Linear(3, 2)
    @test mdl.mdls[2] == Bias(2)
    @test mdl.mdls[3] == HyperbolicTangent(2)
    @test mdl.mdls[4] == Linear(2, 3)
    @test mdl.mdls[5] == Bias(3)
    @test mdl.mdls[6] == HyperbolicTangent(3)
    @test mdl.mdls[7] == Linear(3, 4)
    @test mdl.mdls[8] == Bias(4)
    @test mdl.mdls[9] == HyperbolicTangent(4)
    @test mdl.mdls[10] == Linear(4, 2)
    @test mdl.mdls[11] == Bias(2)
    @test mdl.mdls[12] == Identity(2)
end


text:
function read_VolumeImage(fname::AbstractString)

    @info("Creating volume image from file $fname")

    if occursin(".dat", fname)
        x, y, z, s, t = read_dat(fname)
        method = "CLARA"
        header = Dict()
        units = "nAm/cm^3"
        x = x / 1000
        y = y / 1000
        z = z / 1000
        t = t / 1000
    else
        @warn("Unknown file type")
    end

    header["FileName"] = fname

    coord_system = "?"

    VolumeImage(
        s,
        units,
        collect(x) * u"m",
        collect(y) * u"m",
        collect(z) * u"m",
        collect(t) * u"s",
        method,
        header,
        coord_system,
    )
end


text:
# functions to support transformation of general real or complex matrices to
# Jordan normal form.
#

"""
Transform upper (quasi) triangular input matrix to block-diagonal form.
  The eigenvalues of the matrix must be ordered in a way, that multiple
  eigenvalues are stored in contiguous parts of the diagonal.
  The results stored in place of the input matrix.
"""
function blockdiagonal!(A::StridedMatrix{T}, ilo::Int, ihi::Int) where T <: Number

  while ihi > ilo
    # find lower right part of diagonal with equal eigenvalues
    isep = findlastequal(A, ilo, ihi)
    if isep > ilo
      rA = ilo:isep-1
      rB = isep:ihi
      A[rA, rB] .= sylvester(view(A, rA, rA), -view(A, rB, rB), A[rA, rB])
    end
    ihi = isep - 1
  end
  A
end

"""
Find diagonal position with first eigenvalue equal to last eigenvalue
"""
function findlastequal(A::AbstractMatrix{T}, ilo::Int, ihi::Int) where T <: Complex
  k = ihi
  d = A[k,k]
  while k > ilo && A[k-1,k-1] == d
    k -= 1
  end
  k
end

function findlastequal(A::AbstractMatrix{T}, ilo::Int, ihi::Int) where T <: Real
  k = ihi
  if ihi - ilo > 1
    if A[k,k-1] == 0
      d = A[k,k]
      while k > ilo && A[k-1,k-1] == d
        k -= 1
      end
    else
      k -= 1
      d, di = ev2(A, k)
      di > 0 || error("not expecting real eigenvalue in 2x2-block")
      while k > ilo + 1
        d2, di2 = ev2(A, k-2)
        if di2 > 0 && di2 == di && d2 == d
          k -= 2
        else
          break
        end
      end
    end
  end
  k
end

function ev2(A::AbstractMatrix{T}, k::Int) where T <: Real
  dis = discriminant(A, k, k+1)
  di = sqrt(max(-dis, zero(dis)))
  d = (A[k,k] + A[k+1,k+1]) / 2
  d, di
end

max_stars_repo_name: KlausC/LinAlgBigFloat

max_issues_repo_name: KlausC/LinAlgBigFloat

max_forks_repo_name: KlausC/LinAlgBigFloat


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    signarea(A, B, C)

Compute signed area of triangle formed
by points `A`, `B` and `C`.
"""
function signarea(A::Point{2}, B::Point{2}, C::Point{2})
  ((B - A) × (C - A)) / 2
end

"""
    sideof(point, segment)

Determines on which side of the oriented `segment`
the `point` lies. Possible results are `:LEFT`,
`:RIGHT` or `:ON` the segment.
"""
function sideof(p::Point{2,T}, s::Segment{2,T}) where {T}
  a, b = vertices(s)
  area = signarea(p, a, b)
  ifelse(area > atol(T), :LEFT, ifelse(area < -atol(T), :RIGHT, :ON))
end

"""
    sideof(point, chain)

Determines on which side of the closed `chain` the
`point` lies. Possible results are `:INSIDE` or
`:OUTSIDE` the chain.
"""
function sideof(p::Point{2,T}, c::Chain{2,T}) where {T}
  w = windingnumber(p, c)
  ifelse(isapprox(w, zero(T), atol=atol(T)), :OUTSIDE, :INSIDE)
end


text:
using CImGui
using CImGui.CSyntax
using CImGui.CSyntax.CStatic
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using CImGui.GLFWBackend.GLFW
using CImGui.OpenGLBackend.ModernGL
using Printf

@static if Sys.isapple()
    # OpenGL 3.2 + GLSL 150
    const glsl_version = 150
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
    GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
else
    # OpenGL 3.0 + GLSL 130
    const glsl_version = 130
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
    # GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
end

# setup GLFW error callback
error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
GLFW.SetErrorCallback(error_callback)

# custom callbacks
function CustomUserCallbackMousebutton(window::GLFW.Window, button::GLFW.MouseButton, action::GLFW.Action, mods::Cint)
    println("Trigger user mouse button callback!")
    return nothing
end
SetCustomMouseButtonCallback(CustomUserCallbackMousebutton)

# create window
window = GLFW.CreateWindow(1280, 720, "Demo")
@assert window != C_NULL
GLFW.MakeContextCurrent(window)
GLFW.SwapInterval(1)  # enable vsync

# setup Dear ImGui context
ctx = CImGui.CreateContext()

# setup Dear ImGui style
CImGui.StyleColorsDark()

# setup Platform/Renderer bindings
ImGui_ImplGlfw_InitForOpenGL(window, true)
ImGui_ImplOpenGL3_Init(glsl_version)

show_demo_window = true
clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]
while !GLFW.WindowShouldClose(window)
    # oh my global scope
    global show_demo_window
    global clear_color

    GLFW.PollEvents()
    # start the Dear ImGui frame
    ImGui_ImplOpenGL3_NewFrame()
    ImGui_ImplGlfw_NewFrame()
    CImGui.NewFrame()

    # show the big demo window
    show_demo_window && @c CImGui.ShowDemoWindow(&show_demo_window)

    # rendering
    CImGui.Render()
    GLFW.MakeContextCurrent(window)
    display_w, display_h = GLFW.GetFramebufferSize(window)
    glViewport(0, 0, display_w, display_h)
    glClearColor(clear_color...)
    glClear(GL_COLOR_BUFFER_BIT)
    ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())

    GLFW.MakeContextCurrent(window)
    GLFW.SwapBuffers(window)
end

# cleanup
ImGui_ImplOpenGL3_Shutdown()
ImGui_ImplGlfw_Shutdown()
CImGui.DestroyContext(ctx)

GLFW.DestroyWindow(window)


text:
module Sh

compressors() = Dict{Symbol,AbstractString}(:f=>"cat",:gz=>"gzip")

export viatmp
"viatmp(\"aa.txt\") -> Cmd"
viatmp{S<:AbstractString}(file::S) = viatmp(:f=>file)

export viatmpgz
"viatmpgz(\"aa.gz\") -> Cmd"
viatmpgz{S<:AbstractString}(file::S) = viatmp(:gz=>file)

"""viatmp(:gz=>\"aa.gz\") -> `bash -c \"gzip > aa.gz.TMP && mv aa.gz.TMP aa.gz\"`"""
viatmp{S<:AbstractString}(p::Pair{Symbol,S}) = viatmp(p, compressors())

"""viatmp(:myzip=>aa.zip, Dict(:myzip=>\"zip\"))"""
function viatmp{S1<:AbstractString, S2<:AbstractString}(
		p::Pair{Symbol,S1}, dict::Dict{Symbol,S2} )
 ftype = p[1]
 fname = p[2]
 compr = dict[ftype]
 c("""$compr > \"$fname.TMP\" && mv \"$fname.TMP\" \"$fname\"""")
end

export c
"c(\"echo aa | grep aa\") -> `bash -c \"echo aa | grep aa\"`"
c{S<:AbstractString}(command::S) = `bash -c """$command"""`

export @c_str
""" c\" echo aa | grep aa \" -> `bash -c \"echo aa | grep aa\"` """
macro c_str(command) c(command) end

export @rc_str
""" rc\" command \" # run command"""
macro rc_str(command) run(c(command)) end

end # module

max_stars_repo_name: closescreen/Sh

max_issues_repo_name: closescreen/Sh

max_forks_repo_name: closescreen/Sh


text:
const GIT_USER_NAME = "UnofficialJuliaMirrorBot"
const GIT_USER_EMAIL = "43731649+UnofficialJuliaMirrorBot@users.noreply.github.com"


text:
module PeriodicIndexing

export PeriodicCartesianIndices

using Base: tail

import Base: getindex, length, size, IndexStyle, IndexCartesian, in, iterate, mod,
    first, last

mod_tuple(c::NTuple{N,Int}, size::NTuple{N,Int}) where N  =
    ntuple(k->mod(c[k]-1,size[k])+1, Val(N))

mod(c::CartesianIndex{N}, size::NTuple{N,Int}) where N =
    CartesianIndex(mod_tuple(c.I, size))

add_offset_mod(i::NTuple{N,Int}, start::NTuple{N,Int}, size::NTuple{N,Int}, periodic::NTuple{N,Bool}) where N =
    ntuple(k->( periodic[k]  ? add_offset_mod(i[k], start[k], size[k]) : i[k] ), Val(N))

add_offset_mod(i::Int, start::Int, size::Int) =
    start + mod(i-start, size)

struct PeriodicCartesianIndices{N,K} <:AbstractArray{CartesianIndex{N},N}
    size::NTuple{N,Int}
    iter::CartesianIndices{N,K}
    start::CartesianIndex{N}
    stop::CartesianIndex{N}
    periodic::NTuple{N,Bool}
end

getindex(A::PeriodicCartesianIndices{N}, I::Vararg{Int,N}) where N = mod(getindex(A.iter,I...),A.size)
IndexStyle(A::PeriodicCartesianIndices) = IndexCartesian()
function PeriodicCartesianIndices(size::NTuple{N,Int}, iranges::UnitRange...) where {N}
    start = CartesianIndex(ntuple(k->first(iranges[k]), Val(N)))
    stop  = CartesianIndex(ntuple(k->last(iranges[k]), Val(N)))
    PeriodicCartesianIndices(size, CartesianIndices(iranges), start, stop, ntuple(k->true,Val(N)))
end

function PeriodicCartesianIndices(size::NTuple{N,Int}, start::CartesianIndex{N}, stop::CartesianIndex,  periodic::NTuple{N,Bool}=ntuple(k->true,Val(N))) where {N}
    iranges = ntuple(k-> ((start[k] <= stop[k]) ? (start[k]:stop[k]) : (start[k]:(stop[k]+size[k]))), Val(N))
    PeriodicCartesianIndices(size, CartesianIndices(iranges), start, stop, periodic)
end

size(m::PeriodicCartesianIndices) = size(m.iter)
length(m::PeriodicCartesianIndices) = length(m.iter)

struct PeriodicUnitRange{I} <:AbstractUnitRange{I}
    size::I
    iter::UnitRange{I}
    periodic::Bool

    function PeriodicUnitRange(size::I, iter::UnitRange{I}, periodic=true) where I<:Integer
        if periodic
            e = last(iter)
            l = length(iter)
            e_aligned = mod(e-1,size)+1
            f_aligned = e_aligned - l + 1
            new{I}(size, f_aligned:e_aligned, periodic)
        else
            new{I}(size, max(1,start(iter)):min(size,last(iter)), periodic)
        end
    end
end

first(r::PeriodicUnitRange) = mod(first(r.iter)-1,r.size)+1
last(r::PeriodicUnitRange) = mod(last(r.iter)-1,r.size)+1

function in(i::Integer, r::PeriodicUnitRange{<:Integer})
    if r.periodic
        i_aligned =  mod(i-1,r.size)+1
        return (i_aligned ∈ r.iter) | (i_aligned ∈ r.iter .+ r.size)
    else
        return in(i, r.iter)
    end
end

function iterate(m::PeriodicUnitRange)
    tuple = iterate(m.iter)
    if tuple != nothing
        index, state = tuple
        mod(index-1, m.size)+1, state
    end
end

function iterate(m::PeriodicUnitRange, state)
    tuple = iterate(m.iter, state)
    if tuple != nothing
        index, state = tuple
        mod(index-1, m.size)+1, state
    end
end

length(m::PeriodicUnitRange) = length(m.iter)

nbindexlist(index::Base.CartesianIndex{N}, size::NTuple{N,Int}, periodic=ntuple(k->false, Val(N))) where {N}  =
    PeriodicCartesianIndices(size, index+CartesianIndex(ntuple(k->-1, Val(N))), index+CartesianIndex(ntuple(k->1, Val(N))), periodic)
nbindexlist(index::Int, size::NTuple{1,Int}, periodic=false)  =
    PeriodicUnitRange(size[1], index-1:index+1, periodic)

end # module