text:
module MainLoopKernel
using CUDA, Logging,..CUDAGpuUtils, ..ResultListUtils,..WorkQueueUtils,..ScanForDuplicates, Logging,StaticArrays, ..IterationUtils, ..ReductionUtils, ..CUDAAtomicUtils,..MetaDataUtils
using ..BitWiseUtils,..MetadataAnalyzePass, ..ScanForDuplicates, ..ProcessMainDataVerB
export mainKernelLoad,@mainLoopKernelAllocations,getSmallGPUForHousedorff,getBigGPUForHousedorffAfterBoolKernel,@loadDataAtTheBegOfDilatationStep,@mainLoopKernel, @iterateOverWorkQueue,@mainLoop,@mainLoopKernelAllocations,@clearBeforeNextDilatation

"""
clearing data in order to enable their reuse in next iteration
clearIterResShmemLoop - given we trat res shmem as one dimensional array with 32 entries per iteration how many times we need to loop to cover all
clearIterSourceShmemLoop - the same as above but for source shmem
resShmemTotalLength, sourceShmemTotalLength - total (treating as 1 dimensional array) length of the  sourse or result shared memory
"""
macro clearBeforeNextDilatation( clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength, sourceShmemTotalLength,dataBdim)
    return esc(quote
    isMaskFull=true
    # @ifXY 1 1  CUDA.@cuprint "rrrr $(resShmem[10,10,10]) \n" #lll $(length(resShmem)) \n"
    #  for i in 1:length(resShmem) #resShmemTotalLength-1
    #     resShmem[i]=false
    #  end
    # @iterateLinearly clearIterResShmemLoop resShmemTotalLength  resShmem[i]=false
    # @iterateLinearly clearIterSourceShmemLoop sourceShmemTotalLength  sourceShmem[i]=false
    # @iterateLinearly clearIterResShmemLoop 34*20*34  resShmem[i]=false
    # @iterateLinearly clearIterSourceShmemLoop 34*20*34  sourceShmem[i]=false

    @ifY 1 if(threadIdxX()<15) areToBeValidated[threadIdxX()]=false end
    @ifY 2 if(threadIdxX()<8) isAnythingInPadding[threadIdxX()]=false end

end)#quote
end#clearBeforeNextDilatation

"""
allocates memory in the kernel some register and shared memory  (no allocations in global memory here from kernel)
dataBdim - are dimensions of data block - each data block has just one row in the metadata ...

"""
macro mainLoopKernelAllocations(dataBdim)
    return esc(quote
#needed to manage cooperative groups functions
grid_handle = this_grid()

shmemblockData = @cuDynamicSharedMem(UInt32,($dataBdim[1], $dataBdim[2] ,2))
#shmemblockData = @cuStaticSharedMem(Int32,(32, 32 ,2))
# holding values of results
# resShmemblockData = @cuDynamicSharedMem(UInt32,($dataBdim[1], $dataBdim[2] ))

# holding data about result  2)bottom, 3)left 4)right , 5)anterior, 6)posterior , 7) top,  paddings
shmemPaddings = @cuStaticSharedMem(Bool,( 32,32,7 ))

#for storing sums for reductions
shmemSum =  @cuStaticSharedMem(UInt32,(36,16)) # we need this additional spots
#used to load from metadata information are ques to be validated
areToBeValidated =  @cuStaticSharedMem(Bool, (14))
# used to mark wheather there is any true in paddings
isAnythingInPadding =  @cuStaticSharedMem(Bool, (7))
# used to accumulate counts from of fp's and fn's already covered in dilatation steps
alreadyCoveredInQueues =@cuStaticSharedMem(UInt32,(14))

 #coordinates of data in main array
 #we will use this to establish weather we should mark  the data block as empty or full ...
 isMaskFull= true
 #here we will store in registers data uploaded from mask for later verification wheather we should send it or not
#  locArr= UInt32(0)
 offsetIter = UInt32(0)
#  localOffset= UInt32(0)
 #boolean usefull in the iterating over private part of work queue
 isAnyBiggerThanZero =  @cuStaticSharedMem(Bool,1)
 #loading data to shared memory from global about is it even or odd pass
 iterationNumberShmem =  @cuStaticSharedMem(UInt16,1)
 #current spot in tail - usefull to save where we need to access the tail to get proper block
 #shared memory variables needed to marks wheather we are already finished with  any dilatation step
 goldToBeDilatated =  @cuStaticSharedMem(Bool,1)
 segmToBeDilatated =  @cuStaticSharedMem(Bool,1)
 #true when we have more than 0 blocks to analyze in next iteration
 isEvenPass =  @cuStaticSharedMem(Bool,1)
 #keeping information  weather we have not odd or even pass
 workCounterBiggerThan0 =  @cuStaticSharedMem(Bool,1)
 #holding values of  work queue counter ghlobal
 workCounterInshmem= @cuStaticSharedMem(UInt16,1)
 #workCounter of this block - thread local
 workCounterLocalInShmem = @cuStaticSharedMem(UInt16,1)
 # helper variable used to add data to global counter
workCounterHelper = @cuStaticSharedMem(UInt16,1)

 positionInMainWorkQueaue= @cuStaticSharedMem(UInt16,1)
 @ifXY 1 1 iterationNumberShmem[1]= 0
 @ifXY 2 1 isAnyBiggerThanZero[1]= 0
 @ifXY 3 1 workCounterInshmem[1]= 0
 @ifXY 7 1 goldToBeDilatated[1]= 0
 @ifXY 8 1 segmToBeDilatated[1]= 0
 @ifXY 9 1 workCounterBiggerThan0[1]= 0
 @ifXY 19 1 isEvenPass[1]= true# set to true becouse it would be altered befor first dilatation step

 sync_threads()

end)#quote
end

"""
main loop logic
"""
macro mainLoop()
    return esc(quote
    #we already established work queue for first pass
    @clearBeforeNextDilatation( clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength, sourceShmemTotalLength,dataBdim)

    @loadDataAtTheBegOfDilatationStep()

    sync_threads()

    #now we should have all data needed to analyze padding from previous dilatation
    @iterateOverWorkQueue(begin
        @executeIterPadding(dilatationArrs[shmemSum[shmemIndex*4+4]+1]
        ,referenceArrs[shmemSum[shmemIndex*4+4]+1]
        ,shmemSum[shmemIndex*4+1]#xMeta
        ,shmemSum[shmemIndex*4+2]#yMeta
        ,shmemSum[shmemIndex*4+3]#zMeta
        ,shmemSum[shmemIndex*4+4]#isGold
        ,iterationNumberShmem[1]#iterNumb
        )

    end )
   sync_grid(grid_handle)

    #we get dilatation from block its padding will be analyzed later
    @iterateOverWorkQueue(begin
        @executeDataIter(mainArrDims
        ,dilatationArrs[shmemSum[shmemIndex*4+4]+1]
        ,referenceArrs[shmemSum[shmemIndex*4+4]+1]
        ,shmemSum[shmemIndex*4+1]#xMeta
        ,shmemSum[shmemIndex*4+2]#yMeta
        ,shmemSum[shmemIndex*4+3]#zMeta
        ,shmemSum[shmemIndex*4+4]#isGold
        ,iterationNumberShmem[1]#iterNumb
    )
    end )
    sync_grid(grid_handle)
    #every time we are overwriting work queaue
    workQueueCounter[1]=0

    sync_grid(grid_handle)

     MetadataAnalyzePass.@setMEtaDataOtherPasses(locArr,offsetIter,iterationNumberShmem[1])

end) #quote
end#mainLoop

"""
iterating over elements in work queaue  in order to make it work we need  to
    know how many elements are in the work queue - we need workQueueCounter
    workQueaue - matrix with work queue data
    goldToBeDilatated, segmToBeDilatated - shared memory values needed to establish weather we finished dilatations in given pass
    shmemSum will be used to store loaded data from workqueue
    shmemSumLengthMaxDiv4 - number indicating linear length of shmemSum but reduced such that it will be divisible by 4
    ex - actions invoked on the data block when its xMeta,yMeta,zMeta and is gold pass informations are already known
    """
macro iterateOverWorkQueue(ex )
   return esc(quote
    #first part we load data from work queue to shmem sum
    # we will treat shmemSum as 1 dimensional array and write data from work queue
    #mod 1 - xMeta, mod 2 - uMeta, mod 3 - zMeta mod 4 - isGoldPass
    #first we need to  establish how many items in work queue will be analyzed by this block
   numbOfDataBlockPerThreadBlock = cld(workQueueCounter[1],gridDimX() )

    #we need to stuck all of the blocks data into shared memory 4 entries for each block
    @unroll for outerIter in 0: fld((numbOfDataBlockPerThreadBlock*4),shmemSumLengthMaxDiv4)
        workQuueueLinearIndexOffset = ((((numbOfDataBlockPerThreadBlock)*4 ))*(blockIdxX()-1))+ (outerIter*shmemSumLengthMaxDiv4)
        # now we load all needed data into shared memory
        @iterateLinearly cld(shmemSumLengthMaxDiv4,blockDimX()*blockDimY()) shmemSumLengthMaxDiv4 begin
            #checking if we are in range
            if(i<=shmemSumLengthMaxDiv4 )
                if( ((outerIter*shmemSumLengthMaxDiv4)+i)<=((numbOfDataBlockPerThreadBlock*4)) && (workQuueueLinearIndexOffset +i)<=(workQueueCounter[1]*4)  )
                    shmemSum[i] = workQueue[(workQuueueLinearIndexOffset +i)]
                else
                    shmemSum[i] =0
                end
            end

        end
        sync_threads()
    #at this point we had pushed into shared memory as much data as we can fit so we need to start using it
    #second part proper iteration by definition no rounding needed here
    #also we do not make here any attempt of parallelization as this will be done inside the expression we just provide actual metadata for dilatation step
        for shmemIndex in 0:(fld(shmemSumLengthMaxDiv4,4)-1)
                if( ((shmemIndex+1)*4 <=shmemSumLengthMaxDiv4 ) && shmemSum[shmemIndex*4+1]>0  ) #shmemSum[shmemIndex*4+1]>0
                # checking is there any point in futher dilatations of this block
                if((shmemSum[shmemIndex*4+4]==1 && goldToBeDilatated[1]) || (shmemSum[shmemIndex*4+4]==0 && segmToBeDilatated[1]) )
                    #finally all ready for dilatation step to be executed on this particular block
                    #this one below should be removed after problem is diagnosed
                    if(shmemSum[shmemIndex*4+1]< mainArrDims[1]  && shmemSum[shmemIndex*4+2]< mainArrDims[2]&& shmemSum[shmemIndex*4+3]< mainArrDims[2] && shmemSum[shmemIndex*4+4]< 2)
                        $ex
                    end
                end
            end#if in range
        end# main functional loop for dilatation and validation
    end#outer for
    end)#quote
end

"""
main kernel managing
"""
macro mainLoopKernel()
  return esc(quote

    @mainLoopKernelAllocations(dataBdim)

    @clearBeforeNextDilatation(clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength, sourceShmemTotalLength,dataBdim)

    MetadataAnalyzePass.@analyzeMetadataFirstPass()

    sync_grid(grid_handle)

  @loadDataAtTheBegOfDilatationStep()

    sync_grid(grid_handle)
   #we check first wheather next dilatation step should be done or not we also establish some shared memory variables to know wheather both passes should continue or just one
    # checking weather we already finished so we need to check
    # - is amount of results related to gold mask dilatations is equal to false positives or given percent of them
    # - is amount of results related to other  mask dilatations is equal to false negatives or given percent of them
    # - is amount of workQueue that we will want to analyze now is bigger than 0

   #while((goldToBeDilatated[1] || segmToBeDilatated[1]) && workCounterBiggerThan0[1])
    for i in 1:11
        @ifXY 1 1 if(blockIdxX()==1) CUDA.@cuprint " goldToBeDilatated[1] $(goldToBeDilatated[1])  segmToBeDilatated[1] $(segmToBeDilatated[1]) workCounterBiggerThan0[1] $(workCounterBiggerThan0[1]) iterationNumberShmem[1] $(iterationNumberShmem[1]) \n" end
        @mainLoop()
    end
   #end#while we did not yet finished
    #this will basically give the main result
    globalIterationNumb[1]=iterationNumberShmem[1]
end)#quote

end

# """
# after dilatation prepare
# """
# macro prepareForNextDilation()
#     return esc(quote
#     #we update metadata and prepare work queue
#     setMEtaDataOtherPasses()
#     #we clear  and add negation to !isOddPassShmem becouse we want to set the previously updated counter to 0
#     @clearBeforeNextDilatation( clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength, sourceShmemTotalLength,dataBdim)

#     @loadDataAtTheBegOfDilatationStep()
# end)#quote
# end    #prepareForNextDilation

"""
loads data at the begining of each dilatation step
    we need to set some variables in shared memory ro initial values
    fp, fn
"""
macro loadDataAtTheBegOfDilatationStep(  )

    return esc(quote
    #so we know that becouse of sync grid we will have evrywhere the same  iterationNumberShmem and positionInMainWorkQueaue

    @ifXY 1 1 iterationNumberShmem[1]+=1
    #@ifXY 2 2 positionInMainWorkQueaue[1]=0

    @ifXY 2 1 begin
        workCounterInshmem[1]= workQueueCounter[1]
        workCounterBiggerThan0[1]= (workCounterInshmem[1]>0)
                    end
    #we do corection for robustness so we can ignore some of the most distant points - this will reduce the influence of outliers
    @ifXY 3 1 goldToBeDilatated[1]=(globalCurrentFpCount[1] <= ceil(fp[1]*robustnessPercent))
    @ifXY 4 1 segmToBeDilatated[1]=(globalCurrentFnCount[1] <= ceil(fn[1]*robustnessPercent))
        end)#quote
    end

"""
allocates memory for small GPU Arrays
"""
function getSmallGPUForHousedorff()
    globalFpResOffsetCounter= CUDA.zeros(UInt32,1)
    globalFnResOffsetCounter= CUDA.zeros(UInt32,1)
    workQueueCounter= CUDA.zeros(UInt32,1)
    globalIterationNumber = CUDA.zeros(UInt32,1)
    globalCurrentFnCount= CUDA.zeros(UInt32,1)
    globalCurrentFpCount= CUDA.zeros(UInt32,1)
    globalIterationNumb= CUDA.zeros(UInt32,1)

return (globalFpResOffsetCounter,globalFnResOffsetCounter,globalIterationNumber,globalCurrentFnCount,globalCurrentFpCount,globalIterationNumb )
end

# """
# allocate memory for bigger arrays
# """
# function getBigGPUForHousedorff()

# end

"""
allocate after prepare bool kernel had finished execution
return metaData,reducedGoldA,reducedSegmA,reducedGoldB,reducedSegmB,resList,workQueue,,workQueuecounter
"""
function getBigGPUForHousedorffAfterBoolKernel(metaData,minxRes,maxxRes,minyRes,maxyRes,minzRes,maxzRes,fn,fp,reducedGoldA,reducedSegmA,dataBdim)
    ###we return only subset of boolean arrays that we are intrested in
    # println( "zzz fp[1] $(fp[1])  fn[1] $(fn[1]) \n")
    # resList = allocateResultLists(fp[1],fn[1])
    resList = allocateResultLists(1000,1000)
    ###we need to return subset of metadata that we are intrested in
    goldArr= reducedGoldA[Int64(minxRes[1]*dataBdim[1]):Int64(maxxRes[1]*dataBdim[1])
            ,Int64(minyRes[1]*dataBdim[2]):Int64(maxyRes[1]*dataBdim[2])
            ,Int64(minzRes[1]):Int64(maxzRes[1])
            ]
    segmArr = reducedSegmA[Int64(minxRes[1]*dataBdim[1]):Int64(maxxRes[1]*dataBdim[1])
                ,Int64(minyRes[1]*dataBdim[2]):Int64(maxyRes[1]*dataBdim[2])
                ,Int64(minzRes[1]):Int64(maxzRes[1])
                ]

    # CUDA.unsafe_free!(reducedGoldA)

    # CUDA.unsafe_free!(reducedSegmA)
    newMeta = metaData[Int64(minxRes[1]):Int64(maxxRes[1]),Int64(minyRes[1]):Int64(maxyRes[1]),Int64(minzRes[1]):Int64(maxzRes[1]),:   ]
    workQueue,workQueueCounter= WorkQueueUtils.allocateWorkQueue( max(length(newMeta),1) )
    metaSize = size(newMeta)

    #here we need to define data structure that will hold all data about paddings
    # it will be composed of 32x32 blocks of UInt8 where first 6 bits will represent paddings
    # number of such blocks will be the same as innersingleDataBlockPass
    paddingStore = CUDA.zeros(UInt8, metaSize[1],metaSize[2],metaSize[3],32,32)

    return(newMeta
            ,goldArr  ,segmArr ,paddingStore
            ,resList,workQueue,workQueueCounter
    )
end

end#MainLoopKernel

# """
# analyzing private part of work queue
# """
# macro privateWorkQueueAnalysis()
#     return esc(quote

#     loadOwnedWorkQueueIntoShmem(mainWorkQueue,mainQuesCounter,toIterWorkQueueShmem,positionInMainWorkQueaue ,numberOfThreadBlocks,tailParts)
#     #at this point if we have anything in the private part of the  work queue we will have it in toIterWorkQueueShmem, in case  we will find some 0 inside it means that queue is exhousted and we need to loo into tail
#     @unroll for i in UInt16(1):32# most outer loop is responsible for z dimension
#         if(toIterWorkQueueShmem[i,1]>0)
#             #we need to check also wheather given dilatation step is already finished - for example it is possible that it do not make sense to dilatate gold mask but still we need dilate other
#             if( (goldToBeDilatated[1]&&toIterWorkQueueShmem[i,4]==1) || (segmToBeDilatated[1]&&toIterWorkQueueShmem[i,4]==0)  )
#                 @innersingleDataBlockPass(toIterWorkQueueShmem[i,4]==1,toIterWorkQueueShmem[i,1],toIterWorkQueueShmem[i,2] ,toIterWorkQueueShmem[i,3] )
#             end
#             sync_threads()
#         else    #at this point we got some 0 in the toIterWorkQueueShmem
#             @ifXY 1 1  isAnyBiggerThanZero[]=false
#             sync_threads()
#         end#if
#     end#for
# end)#quote
# end

# """
# analyzing tail part of work queue
# """
# macro analyzeTail()
#     return esc(quote

#     #below we access tail of working queue in a way that will be atomic
#     @ifXY 1 1 toIterWorkQueueShmem[1]= mainWorkQueueArr[isOddPassShmem[]+1,currentTailPosition[1]]
#     sync_threads()
#     #we need to check also wheather given dilatation step is already finished - for example it is possible that it do not make sense to dilatate gold mask but still we need dilate other
#     if( (goldToBeDilatated[1]&&toIterWorkQueueShmem[1,4]==1) || (segmToBeDilatated[1]&&toIterWorkQueueShmem[1,4]==0)  )
#         @innersingleDataBlockPass(toIterWorkQueueShmem[1,4]==1,toIterWorkQueueShmem[1,1],toIterWorkQueueShmem[1,2] ,toIterWorkQueueShmem[1,3] )
#     end
#     loadDataNeededForTailAnalysisToShmem(currentTailPosition,tailCounter )
#     sync_threads()
# end) #quote
# end


text:
struct OrbitalFieldMatrix <: AbstractEnvironmentFeatureDescriptor
    # stuff
end

# more stuff

ext: jl


text:
GRID_LENGTH = 10
FLASH_THRESHOLD = 9
NEGATIVE_INF_INT = -10000000

function printarr(arr)
  for x in eachrow(arr)
    println(x)
  end
  println()
end

function part1(fname, n_steps)
  octopuses = Array{Int}(undef, (GRID_LENGTH, GRID_LENGTH))
  open(fname, "r") do f
    data = readlines(f)
    for (row_idx, row) in enumerate(data),
        (column_idx, number) in enumerate(row)
        octopuses[row_idx, column_idx] = parse(Int, number)
    end
  end
  n_flashes = 0
  for step in 1:1:n_steps
    octopuses .+= 1
    flashing = octopuses .> FLASH_THRESHOLD
    new_flashes = any(flashing)
    n_flashes += sum(flashing)
    octopuses[flashing] .= NEGATIVE_INF_INT
    while new_flashes
      # This loop is ugly but I did not manage to find an elegant method using,
      # e.g., array slicing.
      for col_idx in 1:1:size(flashing, 1), row_idx in 1:1:size(flashing, 2)
        if !flashing[col_idx, row_idx]
          continue
        end
        for dx in (-1, 0, 1), dy in (-1, 0, 1)
          if (dx == 0) && (dy == 0)
            continue
          end
          try
            octopuses[col_idx+dy, row_idx+dx] += 1
          catch e
            if isa(e, BoundsError)
              continue
            else
              throw(e)
            end
          end
        end
      end
      flashing = octopuses .> FLASH_THRESHOLD
      new_flashes = any(flashing)
      n_flashes += sum(flashing)
      octopuses[flashing] .= NEGATIVE_INF_INT
    end
    if all(octopuses .< 0)
      println("All octopuses flashing at step $step")
      break
    end
    octopuses[octopuses .< 0] .= 0
  end
  println("number of flashes: $n_flashes")
end

part1("test_input.txt", 100)
part1("test_input.txt", 1000)

part1("input.txt", 100)
part1("input.txt", 10000)


text:
module Transport
using LinearAlgebra, StaticArrays
using LightXML, Printf

include("Utils.jl")
include("Constants.jl")
include("IdealGas.jl")
#defne the collision integral Omega11
omega11 = @SMatrix[
    4.008e+000 -1.0380e+00  5.9659e+00 -2.9977e+00  4.9812e-01;
    3.1300e+00 -4.6244e-01  2.2622e+00 -8.4707e-01  1.0961e-01;
    2.6490e+00 -3.6366e-01  1.4543e+00 -4.7309e-01  5.2330e-02;
    2.3140e+00 -2.6379e-01  1.0326e+00 -2.9782e-01  2.7717e-02;
    2.0660e+00 -1.9750e-01  8.0301e-01 -2.1612e-01  1.7991e-02;
    1.8770e+00 -2.5097e-01  9.7452e-01 -4.5890e-01  9.1504e-02;
    1.7290e+00 -1.1931e-01  5.5652e-01 -1.4021e-01  1.0704e-02;
    1.6122e+00 -8.8902e-02  4.6422e-01 -1.0763e-01  6.7026e-03;
    1.5170e+00 -6.6921e-02  3.9820e-01 -8.7569e-02  4.7731e-03;
    1.4400e+00 -5.6468e-02  3.5399e-01 -7.6371e-02  4.0889e-03;
    1.3204e+00 -4.0079e-02  2.8619e-01 -5.9986e-02  3.2544e-03;
    1.2340e+00 -3.4407e-02  2.4399e-01 -5.2429e-02  3.4858e-03;
    1.1680e+00 -1.4465e-02  1.7825e-01 -2.5006e-02 -7.2435e-04;
    1.1166e+00 -1.3643e-02  1.5968e-01 -2.5965e-02  6.4166e-04;
    1.0750e+00 -6.3727e-03  1.2820e-01 -1.4746e-02 -9.7281e-04;
    1.0006e+00 -9.1095e-03  1.0503e-01 -1.6931e-02  6.9808e-04;
    9.5000e-01 -4.6060e-03  7.7117e-02 -1.0881e-02  3.7158e-04;
    9.1310e-01 -1.9581e-03  5.6961e-02 -5.2681e-03 -3.0468e-04;
    8.8450e-01 -2.9286e-03  5.0573e-02 -7.3139e-03  4.8588e-04;
    8.4280e-01 -1.4158e-03  3.3049e-02 -2.7314e-03 -1.2165e-04;
    8.1300e-01 -1.3882e-03  2.4887e-02 -2.0834e-03 -2.3635e-05;
    7.8980e-01 -1.9075e-05  1.7423e-02 -4.1225e-04 -2.2042e-04;
    7.7110e-01  2.0888e-05  1.3897e-02 -2.4013e-04 -1.7556e-04;
    7.5550e-01 -4.5415e-05  1.1903e-02 -6.0156e-04 -3.7943e-05;
    7.4220e-01  5.4903e-04  8.6391e-03  2.7416e-04 -1.5004e-04;
    7.2022e-01 -2.7033e-04  7.7439e-03 -6.1756e-04  4.3613e-05;
    7.0250e-01  8.3428e-04  3.8508e-03  8.2584e-04 -2.0902e-04;
    6.8776e-01  3.9894e-05  3.9815e-03  2.0106e-04 -8.6618e-05;
    6.7510e-01  3.8786e-05  3.3483e-03  5.6978e-05 -3.8795e-05;
    6.6400e-01  3.2925e-04  2.3460e-03  3.2880e-04 -9.2696e-05;
    6.4140e-01 -1.6252e-04  2.1257e-03 -4.1127e-05 -1.6239e-05;
    6.2350e-01  3.8876e-05  1.4220e-03 -2.7050e-05 -1.7042e-06;
    6.0882e-01 -1.2964e-04  1.5607e-03 -4.0253e-04  8.8076e-05;
    5.9640e-01  2.6913e-04  8.8982e-05  6.1738e-04 -1.4247e-04;
    5.7630e-01 -1.8613e-04  9.0887e-04 -2.7653e-04  7.1825e-05;
    5.4150e-01 -5.7840e-05  4.9071e-04 -2.3284e-04  5.6760e-05;
    5.1800e-01  1.5100e-04  7.2383e-04 -5.6918e-04  1.2030e-04
]

# Collision integral Omega22
omega22  = @SMatrix [
    4.1000e+00 -4.9400e-01  5.1705e+00 -2.4579e+00  3.8700e-01;
    3.2630e+00 -4.8580e-01  2.4444e+00 -9.0334e-01  1.1276e-01;
    2.8400e+00 -4.7159e-01  1.5619e+00 -4.5765e-01  4.0936e-02;
    2.5310e+00 -3.6736e-01  1.0905e+00 -2.6097e-01  1.3879e-02;
    2.2840e+00 -2.8365e-01  8.4186e-01 -1.8024e-01  5.8124e-03;
    2.0840e+00 -2.1794e-01  6.9371e-01 -1.4468e-01  4.3636e-03;
    1.9220e+00 -2.0893e-01  6.9031e-01 -1.8193e-01  1.4363e-02;
    1.7902e+00 -1.2932e-01  5.1435e-01 -1.0742e-01  3.8897e-03;
    1.6820e+00 -1.0196e-01  4.5847e-01 -9.7322e-02  4.1072e-03;
    1.5930e+00 -8.5797e-02  4.2360e-01 -9.6560e-02  5.8448e-03;
    1.4550e+00 -4.8607e-02  3.2642e-01 -6.4005e-02  1.9493e-03;
    1.3550e+00 -3.2972e-02  2.7273e-01 -5.1800e-02  1.5361e-03;
    1.2800e+00 -2.4970e-02  2.3581e-01 -4.5617e-02  1.8995e-03;
    1.2220e+00 -1.3736e-02  1.9366e-01 -3.2066e-02  3.8064e-04;
    1.1760e+00 -1.2000e-02  1.6956e-01 -2.7242e-02  3.7995e-04;
    1.0933e+00 -6.3451e-03  1.2403e-01 -1.7897e-02  1.2523e-04;
    1.0390e+00 -2.1846e-03  9.0757e-02 -1.0363e-02 -2.7932e-04;
    9.9960e-01 -1.8453e-03  7.2675e-02 -8.3026e-03  1.1244e-05;
    9.6990e-01 -2.0468e-03  5.9254e-02 -5.6424e-03 -2.8151e-04;
    9.2680e-01 -1.5662e-03  4.1725e-02 -3.9558e-03 -4.8331e-05;
    8.9620e-01 -5.9689e-04  2.9341e-02 -1.8027e-03 -1.5170e-04;
    8.7270e-01  1.2046e-04  2.1443e-02 -6.0037e-04 -2.1344e-04;
    8.5380e-01  2.3847e-04  1.6539e-02 -7.9724e-05 -2.1286e-04;
    8.3790e-01 -2.8887e-04  1.4856e-02 -8.8946e-04 -3.7426e-06;
    8.2430e-01  5.6032e-04  1.0655e-02  3.2629e-04 -1.8930e-04;
    8.0180e-01  5.2094e-04  7.2131e-03  6.8005e-04 -2.1524e-04;
    7.8360e-01 -1.0792e-04  6.8273e-03 -4.4591e-04  5.5406e-05;
    7.6830e-01  5.7900e-04  3.4869e-03  1.1370e-03 -2.8172e-04;
    7.5520e-01  1.0191e-04  3.7993e-03  8.7340e-05 -2.1863e-05;
    7.4360e-01  1.1106e-04  3.1512e-03  1.9598e-04 -6.2322e-05;
    7.1982e-01  5.9318e-05  2.1699e-03  9.4738e-05 -2.8537e-05;
    7.0100e-01  1.2696e-04  1.3406e-03  2.2449e-04 -5.1240e-05;
    6.8545e-01  5.2502e-04 -3.2893e-05  9.4361e-04 -2.1258e-04;
    6.7230e-01  3.0299e-05  1.4509e-03 -5.6073e-04  1.4627e-04;
    6.5100e-01 -3.5478e-05  6.2081e-04  4.7129e-05 -2.3793e-05;
    6.1400e-01 -5.1551e-05  1.6994e-03 -1.2009e-03  2.2986e-04;
    5.8870e-01  4.7810e-04  3.2386e-03 -2.7109e-03  5.3094e-04
]

# Reduced temperature tstart

tstar = @SVector[0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0,
                2.5, 3.0, 3.5, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0,
                25.0, 30.0, 35.0, 40.0, 50.0, 75.0, 100.0
]

"""
Structure for defining the transport properties of pure species as read from transport.dat file
    While storing all parameters are stored in SI units
"""
struct TransportData
    name::String            # Name of the species
    geom::Int64             # monoatomic, linear, non-linear
    ϵ::Float64              # LJ potential well depth (J)
    σ::Float64              # LJ collision diameter (m)
    μ::Float64              # dipole moment in Debye (m^3/2 J^1/2)
    α::Float64              # polarizability in cubic angstroms (m^3)
    Zrot::Float64           # rotational relaxation collision number
end

"""
Function for the calculation of pure species viscosity in Kg/m-s
#    Usage
     viscosity(sp_trd, T, molwt)
-   sp_trd : Array of species transport data
-   T : Temperature
-   molwt : molecular weights array
"""
function viscosity(sp_trd::Array{TransportData}, T::Float64, molwt::Array{Float64})
    #calculate the reduced temperature
    eps_reciprocal = [1/i.ϵ for i in sp_trd]
    red_Temp = kB*T*eps_reciprocal

    #reduced dipole moment
    δ = [td.μ^2/(2*td.ϵ*td.σ^3) for td in sp_trd]
    Ω22 = get_omega(red_Temp,δ,key="Omega22")

    #Calculate pure species viscosity
    vector_term = [sqrt(molwt[i])/(sp_trd[i].σ^2*Ω22[i]) for i in eachindex(sp_trd)]
    return (5.0*(sqrt(π*kB*T/Na))/(16.0*π))*vector_term
end

"""
Calculation of mixture viscosity in Kg/m-s
#   Usage
    viscosity(sp_trd, T, molwt, mole_fracs)
-   sp_trd: Species transport data
-   T: Temperature
-   molwt : molecular weight vector
-   mole_fracs : Mole fractions vector
"""
function viscosity(sp_trd::Array{TransportData}, T::Float64, molwt::Array{Float64}, mole_fracs::Array{Float64})
    μ = viscosity(sp_trd, T, molwt)
    n = length(mole_fracs)
    Φ_kj = zeros(n)
    X = zeros(n)
    for k in 1:n
        for j in 1:n
            ml_wt_r = molwt[k]/molwt[j]
            visc_r = μ[k]/μ[j]
            # Φ_kj[j] = (1/sqrt(8))*(1+(molwt[k]/molwt[j]))^(-0.5) * (1+ sqrt(μ[k]/μ[j]) * (molwt[j]/molwt[k])^(0.25) )^2
            Φ_kj[j] = (1.0/sqrt(8.0))*(1+ml_wt_r)^(-0.5)*(1+ visc_r^(0.5) *(1.0/ml_wt_r)^(0.25))^2
        end
        denom = sum(mole_fracs .* Φ_kj)
        X[k] = μ[k]*mole_fracs[k]/denom

    end
    return sum(X)
end

"""
Function for the calculation of mixture diffusion coefficients
"""
function D_km!(Dkm::Array{Float64}, sp_trd::Array{TransportData}, T::Float64, p::Float64, molwt::Array{Float64}, molefracs::Array{Float64})
    massfracs = similar(molefracs)
    molefrac_to_massfrac!(massfracs,molefracs,molwt)
    #Get the binary diffusion coefficients
    bdc = D_ij(sp_trd,T,p,molwt)

    for j in 1:length(molefracs)
        denom = 0.0
        for k in 1:length(molefracs)
            if j!=k
                denom += molefracs[k]/bdc[j,k]
            end
        end
        Dkm[j] = (1-massfracs[j])/denom
    end
end

"""
Thermal conductivity of the mixture
"""
function thermal_coductivity(sp_trd::Array{TransportData}, T::Float64, p::Float64, thermo_obj::IdealGas.SpeciesThermoObj,molefracs::Array{Float64})
    molwt = thermo_obj.molwt
    sdc = zeros(length(molwt))
    D_ii!(sdc, sp_trd,T,p,molwt)
    μ = viscosity(sp_trd,T,molwt)
    λ = 0
    cp = IdealGas.cp_all(thermo_obj,T)
    cv = cp .- R

    cv_rot = cv_vib = 0
    #Vector of density
    ρ = (p/R/T)*molwt
    t1 = t2 = 0
    cv_trans = 1.5R
    f_trans = f_rot = f_vib = 0
    for i in eachindex(sp_trd)
        if sp_trd[i].geom == 0
            f_trans = 2.5
        else
            if sp_trd[i].geom == 1
                cv_rot = R
                cv_vib = cv[i] - 2.5R
            else
                cv_rot = 1.5R
                cv_vib = cv[i]-3R
            end

            val = sp_trd[i].ϵ/kB/T
            val0 = sp_trd[i].ϵ/kB/298.0

            ft(x) = 1+ 0.5*π^(1.5)*sqrt(x) + ( (π^2/4.0)+2 )*x + π^(1.5)*x^(1.5)
            Zrot = sp_trd[i].Zrot*ft(val0)/ft(val)

            f_vib = ρ[i]*sdc[i]/μ[i]
            B = Zrot + (2.0/π)*( (5cv_rot/3R) + f_vib )
            A = 2.5-f_vib
            AB = A/B
            f_trans = 2.5*(1- (2cv_rot/π/cv_trans)*AB )
            f_rot = f_vib*(1+2AB/π)
        end
        λ = μ[i]*(f_trans*cv_trans + f_rot*cv_rot + f_vib*cv_vib)/molwt[i]
        t1 += molefracs[i]*λ
        t2 += molefracs[i]/λ
    end
    return 0.5*(t1+1.0/t2)
end

"""
Function for the calculation of binary diffusion coefficients
"""
function D_ij(sp_trd::Array{TransportData}, T::Float64, p::Float64, molwt::Array{Float64})
    n = length(molwt)

    function reduced_T_and_delta(trd_p, trd_np)

        ϵ_jk = σ_jk = μ_jk_sq = 0
        if trd_p.α == 0.0 && trd_np.α == 0.0 || trd_p.α != 0.0 && trd_np.α != 0.0
            ϵ_jk = sqrt(trd_p.ϵ*trd_np.ϵ)
            σ_jk= 0.5*(trd_p.σ+trd_np.σ)
            μ_jk_sq = trd_p.μ * trd_np.μ
        else
            if trd_p.α != 0
                trd_p, trd_np = trd_np, trd_p
            end
            α_n = trd_np.α/trd_np.σ^3
            μ_p = trd_p.μ/sqrt(trd_p.ϵ*trd_p.σ^3)
            ξ = 1+ 0.25 * α_n * μ_p * sqrt(trd_p.ϵ/trd_np.ϵ)
            ϵ_jk = ξ^2 * sqrt(trd_p.ϵ * trd_np.ϵ)
            σ_jk = 0.5*(trd_np.σ + trd_p.σ)*ξ^(-1.0/6.0)
            μ_jk_sq = 0
        end
        rT = kB*T/ϵ_jk
        δ_jk = 0.5*μ_jk_sq/(ϵ_jk*σ_jk^3)
        return rT,δ_jk,σ_jk
    end

    red_Temp = zeros(n,n)
    δ = zeros(n,n)
    σ = zeros(n,n)
    for j in 1:n
        red_Temp[j,j],δ[j,j],σ[j,j] = reduced_T_and_delta(sp_trd[j],sp_trd[j])
        for k in j+1:n
            red_Temp[j,k],δ[j,k],σ[j,k] = reduced_T_and_delta(sp_trd[j],sp_trd[k])
            red_Temp[k,j] = red_Temp[j,k]
            δ[k,j] = δ[j,k]
            σ[k,j] = σ[j,k]
        end
    end

    D = zeros(n,n)
    const_term = 3.0*sqrt(2*π*Na*kB^3*T^3)/(16*p*π)
    for j in 1:n
        Ω11 = get_omega(red_Temp[j,:],δ[j,:],key="Omega11")
        for k in j+1:n
            m_jk = molwt[j]*molwt[k]/(molwt[j]+molwt[k])
            D[j,k] = const_term * sqrt(1/m_jk)/(σ[j,k]^2*Ω11[k])
            D[k,j] = D[j,k]
        end
        # Uncomment the following two lines to print the self diffusion coefficients
        # m_jj = 0.5*molwt[j]
        # D[j,j] = const_term * sqrt(1/m_jj)/(σ[j,j]^2*Ω11[j])
    end

   return D

end

"""
Calculate self diffusion coefficients
"""
function D_ii!(D::Array{Float64},sp_trd::Array{TransportData}, T::Float64, p::Float64, molwt::Array{Float64})
    n = length(molwt)
    #Calculate reduced temperature
    eps_reciprocal = [1/i.ϵ for i in sp_trd]
    red_Temp = kB*T*eps_reciprocal
    #Dipole moment
    δ = [td.μ^2/(2*td.ϵ*td.σ^3) for td in sp_trd]

    Ω11 = get_omega(red_Temp,δ,key="Omega11")
    const_term = 3.0*sqrt(2*π*Na*(kB*T)^3)/(16*p*π)
    for j in 1:n
        D[j]= const_term * sqrt(1/molwt[j])/(sp_trd[j].σ^2*Ω11[j])
    end

end

"""
calculation of collision integral
"""
function get_omega(red_Temp::Array{Float64}, δ::Array{Float64};kwargs...)
    T_coordinate = Array{Int64,1}()
    for tr in red_Temp
        local T_xmin = 0
        for i in eachindex(tstar)
            if tr < tstar[i+1] && tr > tstar[i]
                T_xmin = i
                break
            end
        end
        T_xmax = min(length(tstar),T_xmin+3)
        if T_xmax == length(tstar)
            T_xmin = T_xmax-3
        end
        push!(T_coordinate,T_xmin-1)
    end

    o_xx = Matrix{Float64}(undef,length(T_coordinate),3)
    for i in eachindex(T_coordinate)
        if δ[i] < 1e-2
            if kwargs[Symbol("key")]=="Omega22"
                o_xx[i,:] = omega22[T_coordinate[i]:T_coordinate[i]+2,1]
            else
                o_xx[i,:] = omega11[T_coordinate[i]:T_coordinate[i]+2,1]
            end
        else
            if kwargs[Symbol("key")]=="Omega22"
                o_xx[i,:] = polyfit(δ[i],T_coordinate[i],key="Omega22")
            else
                o_xx[i,:] = polyfit(δ[i],T_coordinate[i],key="Omega11")
            end
        end
    end
    red_Temp = [ t < tstar[1] ? tstar[2] : t > 500 ? 500 : t for t in red_Temp]

    Ω_XX = Array{Float64,1}()
    quad_interpolation!(Ω_XX,red_Temp,T_coordinate,o_xx)
    return Ω_XX
end

function quad_interpolation!(Ω::Array{Float64}, red_Temp::Array{Float64}, tc::Array{Int64} ,o22::Matrix{Float64})
    for i in eachindex(red_Temp)
        if red_Temp[i] < 100
            x12 = 0.5*(tstar[tc[i]]+tstar[tc[i]+1])
            x23 = 0.5*(tstar[tc[i]+1]+tstar[tc[i]+2])
            x2x = 0.5*(tstar[tc[i]+1]+red_Temp[i])
            m12 = (o22[i,2]-o22[i,1])/(tstar[tc[i]+1]-tstar[tc[i]])
            m23 = (o22[i,3]-o22[i,2])/(tstar[tc[i]+2]-tstar[tc[i]+1])
            m2x = (m12*(x23-x2x) + m23*(x2x-x12))/(x23-x12)
            # println(o22[i,2], "\t", m2x, "\t", red_Temp[i], "\t", tstar[tc[i]+1])
            push!(Ω,o22[i,2] + m2x*(red_Temp[i]-tstar[tc[i]+1]))
        else
            push!(Ω,0.703 + red_Temp[i]*(-1.46e-3+ red_Temp[i]*(3.57e-6+ red_Temp[i]*(-3.43e-9))))
        end
    end
end

function polyfit(rd::Float64,tc::Int64;kwargs...)
    oxx = Array{Float64,1}()
    n_column = size(omega22,2)
    rdvec = [ rd^i for i in range(start=n_column-1,step=-1,stop=0)]
    reverse!(rdvec)

    if kwargs[Symbol("key")]=="Omega22"
        o22_times_rd = omega22[tc:tc+2,:]' .* rdvec
        oxx = [sum(o22_times_rd[:,i]) for i in 1:size(o22_times_rd,2)]
    elseif kwargs[Symbol("key")]=="Omega11"
        o11_times_rd = omega11[tc:tc+2,:]' .* rdvec
        oxx = [sum(o11_times_rd[:,i]) for i in 1:size(o11_times_rd,2)]
    end
    return oxx
end

"""
A function to read the transport.dat file and create the structure TransportData
#   Usage
    create_transport_data(gasphase, trans_file)
-   gasphase: list of gasphase species
-   trans_file: path to the transport.dat file
"""
function create_transport_data(gasphase::Array{String} , trans_file::AbstractString)
    species_trasport_data = Array{TransportData,1}(undef,length(gasphase))
    open(trans_file) do io
        while !eof(io)
            data_string = readline(io)
            data = split(data_string)
            if in(String(strip(data[1])),gasphase)
                species_trasport_data[get_index(String(strip(data[1])),gasphase)] = get_td_struct(data_string)
            end
        end
    end
    return species_trasport_data
end

function get_td_struct(data_str::String)
    data = split(data_str)
    name = String(strip(data[1]))
    geom = parse(Int64,data[2])
    ljwd = parse(Float64,data[3])*kB   # What is given in the input file is ϵ/kB (Values are stored in J)
    ljcd = parse(Float64,data[4])*Angstroms # Values are stored in m
    dipole = parse(Float64,data[5])*Debye   # Values are stored in m^(3/2) J^(1/2)
    polarizability = parse(Float64,data[6])*Angstroms^3 # Values are in m^3
    zrot = parse(Float64,data[7]) # Unitless
    return TransportData(name,geom,ljwd,ljcd,dipole,polarizability,zrot)
end

"""
This function is for testing the transport code
"""
function run(file_path::AbstractString, lib_dir::AbstractString)
    xmldoc = parse_file(file_path)
    xmlroot = root(xmldoc)

    #get the gasphase species
    gasphase = get_collection_from_xml(xmlroot,"gasphase")
    thermo_file = lib_dir*"therm.dat"
    thermo_all = IdealGas.create_thermo(gasphase,thermo_file)
    mole_fracs = get_molefraction_from_xml(xmlroot,thermo_all.molwt,gasphase)
    local T = get_value_from_xml(xmlroot,"T")
    local p = get_value_from_xml(xmlroot,"p")
    trans_file = lib_dir*"transport.dat"
    sp_tr_data = create_transport_data(gasphase,trans_file )

    #Evaluate viscosity
    μ =  viscosity(sp_tr_data,T,thermo_all.molwt)
    println("\nPure species viscosity:")
    println("     Species \t viscosity(Kg/m-s)")
    println("-----------------------------------")
    for i in eachindex(gasphase)
        @printf("%12s \t      %.4e\n", gasphase[i], μ[i])
    end
    μ = viscosity(sp_tr_data,T,thermo_all.molwt,mole_fracs)
    @printf("Mixture viscosity: %.4e Kg/m-s\n",μ)

    #Evaluate diffusion coefficients
    println("\nMixture diffusion coefficients:")
    Dkm = zeros(length(gasphase))
    D_km!(Dkm,sp_tr_data,T,p,thermo_all.molwt,mole_fracs)
    println("     Species \t Diff.Coeff(m^2/s)")
    println("-----------------------------------")
    for i in eachindex(gasphase)
        @printf("%12s \t      %.4e\n",gasphase[i], Dkm[i])
    end
    bdc =  D_ij(sp_tr_data,T,p,thermo_all.molwt)
    println("\nBinary diffusion coefficients:")
    println("Note: self diffusion coefficients are not printed below:")
    println("-----------------------------------")
    for j in 1:length(gasphase)
        @printf("\t%11s",gasphase[j])
    end
    println()
    for j in 1:length(gasphase)
        @printf("%12s\t",gasphase[j])
        for k in 1:length(gasphase)
            @printf("%.4e\t",bdc[j,k])
        end
        println()
    end

    #Evaluate thermal conductivity
    λ = thermal_coductivity(sp_tr_data,T,p,thermo_all,mole_fracs)
    @printf("\nThermal conductivity of mixture: %.4e (W/m-K)\n", λ)

    return Symbol("Success")
end

end


text:
# New do-block syntax introduced in: https://github.com/JuliaLang/julia/pull/23718
function rewrite_do(expr::Expr)
    expr.head == :do || error("expression is not a do-block")
    call, body = expr.args
    Expr(:call, call.args[1], body, call.args[2:end]...)
end

iskwarg(x::Any) = isa(x, Expr) && (x.head === :parameters || x.head === :kw)

"""
    extract_kwargs(expr::Expr) -> Vector{Union{Expr,Symbol}}

Extract the :parameters and :kw value into an array of :kw expressions
we don't evaluate any expressions for values yet though.
"""
function extract_kwargs(expr::Expr)
    kwargs = Union{Expr,Symbol}[]
    for x in expr.args[2:end]
        if iskwarg(x)
            if x.head === :parameters
                for kw in x.args
                    push!(kwargs, kw)
                end
            else
                push!(kwargs, x)
            end
        end
    end
    return kwargs
end


text:
module Nuklear

include("LibNuklear.jl")
using .LibNuklear

include("backend/GLFW/GLFWBackend.jl")
using .GLFWBackend

end # module


text:
include("../base.jl")

randomInt()=rand()

mutable struct Consequent
    _symbols::Array{String}
    _depth::Number
    _recursive::Bool
    Consequent(symbols::Array{String}) = new(symbols, -Inf, false)
    Consequent(symbols::Array{String}, depth::Number) = new(symbols, depth, false)
end

mutable struct ProductionRule
    _antecedent::String
    _consequent::Array{Consequent}
    _depth::Number

    ProductionRule(antecedent::String, consequent::Array{Consequent}) = new(antecedent, consequent, -Inf)
end

struct Grammar
    _N::Array{String}
    _T::Array{String}
    _P::Array{ProductionRule}
    _S::String
end # struct

getInitialSymbol(gramm::Grammar) = Grammar._S

function hasDepth(productionRule::ProductionRule)
    productionRule._depth != -Inf
end

function getDepth(productionRule::ProductionRule)
    productionRule._depth
end

function getNRules(productionRule::ProductionRule)
    length(productionRule._consequent)
end

function getRule(productionRule::ProductionRule, index::UInt)
    index = index % length(productionRule._consequent) + 1
    productionRule._consequent[index]
end

function getProductionRule(productionRules::Array{ProductionRule}, symbol::String)
    for i=1:length(productionRules)
        if symbol == productionRules[i]._antecedent
            return productionRules[i]
        end
    end
end

function allConsequentAreSet(productionRule::ProductionRule)
    check = true
    for i=1:length(productionRule._consequent)
        if productionRule._consequent[i]._depth == -Inf
            check = false
            break
        end
    end
    return check
end

function isTerminalSymbol(gramm::Grammar, symbol::String)
    check = false
    for i=1:length(gramm._T)
        if gramm._T[i] == symbol
            check = true
            break
        end
    end
    return check
end

function getNRules(gramm::Grammar)
    nRules = 0
    for i=1:length(gramm._P)
        productionRule = gramm._P[i]
        nRules = nRules + getNRules(productionRule)
    end
    return nRules
end

function getNProductionRules(gramm::Grammar)
    return length(gramm._P)
end

function readSymbol(str::String)
    nullChar(element) = findnext(x->x==element, [' ','\n','\r','\t'], 1) != nothing
    count = 0
    len = length(str)
    i = 2
    symbol = ""
    while (str[i]!='>')
        if nullChar(str[i])
            error("There shouldn't be neither whitespaces nor line breaks inside a <symbol>")
        else
            symbol=symbol*str[i]
            i=i+1
        end
        if i > len
            error("While reading a symbol, the string ends prematurely")
        end
    end
    if i == 1
        error("Symbol must have, at least, one character")
    end

    return symbol, i
end

"""
function getPhenotype(GEInd::GEIndividual, Grammar)
    stack = Array{Node}(undef,0)
    pushfirst!(stack, getInitialSymbol(Grammar))

end

"""

println("hola")

N=["expr","op", "pre-op", "var"]
T=["x", "sin", "cos", "tan", "log", "+", "-", "/", "*", "(", ")"]

R=["<expr> -> <expr> <op> <expr>
         | (<expr> <op> <expr>)
         | <pre-op> ( <expr> )
         | <var>",

   "<op> -> +
        | -
        | /
        | *",

    "<pre-op> -> sin
                | cos
                | tan
                | log",

    "<var> -> x"
]

S="expr"

 a=parseNo(N, T, R, S)

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
using ForwardDiff

"Holds the tableau of a fully implicit Runge-Kutta method."
struct TableauFIRK{T} <: AbstractTableauIRK{T}
    @HeaderTableau

    q::CoefficientsRK{T}

    function TableauFIRK{T}(q) where {T}
        if (q.s > 1 && istrilstrict(q.a)) || (q.s==1 && q.a[1,1] == 0)
            @warn "Initializing TableauFIRK with explicit tableau $(q.name).\nYou might want to use TableauERK instead."
        elseif q.s > 1 && istril(q.a)
            @warn "Initializing TableauFIRK with diagonally implicit tableau $(q.name).\nYou might want to use TableauDIRK instead."
        end

        new(q.name, q.o, q.s, q)
    end
end

function TableauFIRK(q::CoefficientsRK{T}) where {T}
    TableauFIRK{T}(q)
end

function TableauFIRK(name::Symbol, order::Int, a::Matrix{T}, b::Vector{T}, c::Vector{T}) where {T}
    TableauFIRK{T}(CoefficientsRK(name, order, a, b, c))
end

# TODO function readTableauFIRKFromFile(dir::AbstractString, name::AbstractString)

"Parameters for right-hand side function of fully implicit Runge-Kutta methods."
mutable struct ParametersFIRK{DT, TT, D, S, ET <: ODE{DT,TT}, FT, JT} <: Parameters{DT,TT}
    equ::ET
    tab::TableauFIRK{TT}
    Δt::TT

    F::FT
    Jconfig::JT

    t::TT
    q::Vector{DT}
end

function ParametersFIRK(equ::ET, tab::TableauFIRK{TT}, Δt::TT) where {DT, TT, ET <: ODE{DT,TT}}
    F = (v,q) -> equ.v(zero(TT), q, v)
    tq = zeros(DT, equ.d)
    tv = zeros(DT, equ.d)
    Jconfig = ForwardDiff.JacobianConfig(F, tv, tq)

    ParametersFIRK{DT, TT, equ.d, tab.s, ET, typeof(F), typeof(Jconfig)}(equ, tab, Δt, F, Jconfig, 0, zeros(DT, equ.d))
end

"""
Fully implicit Runge-Kutta integrator cache.

### Fields

* `q̃`: initial guess of solution
* `ṽ`: initial guess of vector field
* `s̃`: holds shift due to periodicity of solution
* `Q`: internal stages of solution
* `V`: internal stages of vector field
* `Y`: vector field of internal stages
"""
struct IntegratorCacheFIRK{DT,D,S} <: ODEIntegratorCache{DT,D}
    q̃::Vector{DT}
    ṽ::Vector{DT}
    s̃::Vector{DT}

    Q::Vector{Vector{DT}}
    V::Vector{Vector{DT}}
    Y::Vector{Vector{DT}}
    J::Vector{Matrix{DT}}

    function IntegratorCacheFIRK{DT,D,S}() where {DT,D,S}
        Q = create_internal_stage_vector(DT, D, S)
        V = create_internal_stage_vector(DT, D, S)
        Y = create_internal_stage_vector(DT, D, S)
        J = [zeros(DT,D,D) for i in 1:S]
        new(zeros(DT,D), zeros(DT,D), zeros(DT,D), Q, V, Y, J)
    end
end

function IntegratorCache(params::ParametersFIRK{DT, TT, D, S}; kwargs...) where {DT, TT, D, S}
    IntegratorCacheFIRK{DT,D,S}(; kwargs...)
end

"Fully implicit Runge-Kutta integrator."
struct IntegratorFIRK{DT, TT, PT <: ParametersFIRK{DT,TT},
                              ST <: NonlinearSolver{DT},
                              IT <: InitialGuessODE{DT,TT}, N, D, S} <: IntegratorRK{DT,TT}
    params::PT
    solver::ST
    iguess::IT
    cache::IntegratorCacheFIRK{DT,D,S}
end

function IntegratorFIRK(equation::ODE{DT,TT,FT,N}, tableau::TableauFIRK{TT}, Δt::TT; exact_jacobian=false) where {DT,TT,FT,N}
    D = equation.d
    M = equation.n
    S = tableau.s

    # create params
    params = ParametersFIRK(equation, tableau, Δt)

    # create solver
    if exact_jacobian
        solver = create_nonlinear_solver_with_jacobian(DT, D*S, params)
    else
        solver = create_nonlinear_solver(DT, D*S, params)
    end

    # create initial guess
    iguess = InitialGuessODE(get_config(:ig_interpolation), equation, Δt)

    # create cache
    cache = IntegratorCacheFIRK{DT,D,S}()

    # create integrator
    IntegratorFIRK{DT, TT, typeof(params), typeof(solver), typeof(iguess), N, D, S}(
                params, solver, iguess, cache)
end

function initialize!(int::IntegratorFIRK, sol::AtomicSolutionODE)
    sol.t̅ = sol.t - timestep(int)

    int.params.equ.v(sol.t, sol.q, sol.v)

    initialize!(int.iguess, sol.t, sol.q, sol.v,
                            sol.t̅, sol.q̅, sol.v̅)
end

function update_params!(int::IntegratorFIRK, sol::AtomicSolutionODE)
    # set time for nonlinear solver and copy previous solution
    int.params.t  = sol.t
    int.params.q .= sol.q
end

function initial_guess!(int::IntegratorFIRK, sol::AtomicSolutionODE)
    local offset::Int

    # compute initial guess for internal stages
    for i in eachstage(int)
        evaluate!(int.iguess, sol.q, sol.v, sol.q̅, sol.v̅, int.cache.q̃, int.cache.ṽ, tableau(int).q.c[i])
        for k in eachindex(int.cache.V[i], int.cache.ṽ)
            int.cache.V[i][k] = int.cache.ṽ[k]
        end
    end
    for i in eachstage(int)
        offset = ndims(int)*(i-1)
        for k in eachdim(int)
            int.solver.x[offset+k] = 0
            for j in eachstage(int)
                int.solver.x[offset+k] += tableau(int).q.a[i,j] * int.cache.V[j][k]
            end
        end
    end
end

function compute_stages!(x::Vector{ST}, Q::Vector{Vector{ST}}, V::Vector{Vector{ST}}, Y::Vector{Vector{ST}},
                                    params::ParametersFIRK{DT,TT,D,S}) where {ST,DT,TT,D,S}

    local tᵢ::TT

    @assert S == length(Q) == length(V) == length(Y)

    # copy x to Y and compute Q = q + Δt Y
    for i in 1:S
        @assert D == length(Q[i]) == length(V[i]) == length(Y[i])
        for k in 1:D
            Y[i][k] = x[D*(i-1)+k]
            Q[i][k] = params.q[k] + Y[i][k]
        end
    end

    # compute V = v(Q)
    for i in 1:S
        tᵢ = params.t + params.Δt * params.tab.q.c[i]
        params.equ.v(tᵢ, Q[i], V[i])
    end
end

"Compute stages of fully implicit Runge-Kutta methods."
function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersFIRK{DT,TT,D,S}) where {ST,DT,TT,D,S}
    # temporary variables
    local y1::ST
    local y2::ST

    # create caches for internal stages
    Q = create_internal_stage_vector(ST, D, S)
    V = create_internal_stage_vector(ST, D, S)
    Y = create_internal_stage_vector(ST, D, S)

    # compute stages from nonlinear solver solution x
    compute_stages!(x, Q, V, Y, params)

    # compute b = - (Y-AV)
    for i in 1:S
        for k in 1:D
            y1 = 0
            y2 = 0
            for j in 1:S
                y1 += params.tab.q.a[i,j] * V[j][k]
                y2 += params.tab.q.â[i,j] * V[j][k]
            end
            b[D*(i-1)+k] = - Y[i][k] + params.Δt * (y1 + y2)
        end
    end
end

"Compute stages of fully implicit Runge-Kutta methods."
function jacobian!(x::Vector{DT}, jac::Matrix{DT}, cache::IntegratorCacheFIRK{DT,D,S}, params::ParametersFIRK{DT,TT,D,S}) where {DT,TT,D,S}
    local tᵢ::TT

    for i in 1:S
        for k in 1:D
            cache.Y[i][k] = x[D*(i-1)+k]
            cache.Q[i][k] = params.q[k] + cache.Y[i][k]
        end
        tᵢ = params.t + params.Δt * params.tab.q.c[i]
        F = (v,q) -> params.equ.v(tᵢ, q, v)
        ForwardDiff.jacobian!(cache.J[i], params.F, cache.ṽ, cache.Q[i], params.Jconfig)
    end

    jac .= 0

    for i in 1:S
        for j in 1:S
            for k in 1:D
                for l in 1:D
                    m = D*(i-1)+k
                    n = D*(j-1)+l

                    if i == j && k == l
                        jac[m,n] += -1
                    end

                    jac[m,n] += params.Δt * params.tab.q.a[i,j] * cache.J[j][k,l]
                    jac[m,n] += params.Δt * params.tab.q.â[i,j] * cache.J[j][k,l]
                end
            end
        end
    end
end

"Integrate ODE with fully implicit Runge-Kutta integrator."
function integrate_step!(int::IntegratorFIRK{DT,TT}, sol::AtomicSolutionODE{DT,TT}) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int, sol)

    # compute initial guess
    initial_guess!(int, sol)

    # reset cache
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector field at internal stages
    compute_stages!(int.solver.x, int.cache.Q, int.cache.V, int.cache.Y, int.params)

    # compute final update
    update_solution!(sol.q, sol.q̃, int.cache.V, int.params.tab.q.b, int.params.tab.q.b̂, int.params.Δt)

    # copy solution to initial guess
    update!(int.iguess, sol.t, sol.q, sol.v)
end


text:
using MLJ, CategoricalArrays, PrettyPrinting
import DataFrames: DataFrame, nrow
using UrlDownload

LinearRegressor = @load LinearRegressor pkg=GLM
LinearBinaryClassifier = @load LinearBinaryClassifier pkg=GLM

baseurl = "https://raw.githubusercontent.com/tlienart/DataScienceTutorialsData.jl/master/data/glm/"

dfX = DataFrame(urldownload(baseurl * "X3.csv"))
dfYbinary = DataFrame(urldownload(baseurl * "Y3.csv"))
dfX1 = DataFrame(urldownload(baseurl * "X1.csv"))
dfY1 = DataFrame(urldownload(baseurl * "Y1.csv"));

first(dfX, 3)

first(dfY1, 3)

ms = models() do m
    AbstractVector{Count} <: m.target_scitype
end
foreach(m -> println(m.name), ms)

ms = models() do m
    Vector{Continuous} <: m.target_scitype
end
foreach(m -> println(m.name), ms)

X = copy(dfX1)
y = copy(dfY1)

coerce!(X, autotype(X, :string_to_multiclass))
yv = Vector(y[:, 1])

LinearRegressorPipe = Pipeline(
    Standardizer(),
    OneHotEncoder(drop_last = true),
    LinearRegressor()
)

LinearModel = machine(LinearRegressorPipe, X, yv)
fit!(LinearModel)
fp = fitted_params(LinearModel)

ŷ = MLJ.predict(LinearModel, X)
yhatResponse = [ŷ[i,1].μ for i in 1:nrow(y)]
residuals = y .- yhatResponse
r = report(LinearModel)

k = collect(keys(fp.fitted_params_given_machine))[3]
println("\n Coefficients:  ", fp.fitted_params_given_machine[k].coef)
println("\n y \n ", y[1:5,1])
println("\n ŷ \n ", ŷ[1:5])
println("\n yhatResponse \n ", yhatResponse[1:5])
println("\n Residuals \n ", y[1:5,1] .- yhatResponse[1:5])
println("\n Standard Error per Coefficient \n", r.linear_regressor.stderror[2:end])

round(rms(yhatResponse, y[:,1]), sigdigits=4)

X = copy(dfX)
y = copy(dfYbinary)

coerce!(X, autotype(X, :string_to_multiclass))
yc = CategoricalArray(y[:, 1])
yc = coerce(yc, OrderedFactor)

LinearBinaryClassifierPipe = Pipeline(
    Standardizer(),
    OneHotEncoder(drop_last = true),
    LinearBinaryClassifier()
)

LogisticModel = machine(LinearBinaryClassifierPipe, X, yc)
fit!(LogisticModel)
fp = fitted_params(LogisticModel)

ŷ = MLJ.predict(LogisticModel, X)
residuals = [1 - pdf(ŷ[i], y[i,1]) for i in 1:nrow(y)]
r = report(LogisticModel)

k = collect(keys(fp.fitted_params_given_machine))[3]
println("\n Coefficients:  ", fp.fitted_params_given_machine[k].coef)
println("\n y \n ", y[1:5,1])
println("\n ŷ \n ", ŷ[1:5])
println("\n residuals \n ", residuals[1:5])
println("\n Standard Error per Coefficient \n", r.linear_binary_classifier.stderror[2:end])

yMode = [mode(ŷ[i]) for i in 1:length(ŷ)]
y = coerce(y[:,1], OrderedFactor)
yMode = coerce(yMode, OrderedFactor)
confusion_matrix(yMode, y)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
@testset "DiskStack" begin
ars = [rand(10,5) for i = 1:3]
cars = cat(ars...,dims=3)
dars = map(a->_DiskArray(a,chunksize=(5,2)),ars)
ds = diskstack(dars)
@test size(ds) == (10,5,3)
@test ds isa DiskArrayStack{Float64,3,_DiskArray{Float64,2,Array{Float64,2}},1}
@test isapprox(sum(ds), sum(cars))
@test all(i->getindex_count(i) == 6,dars)
@test ds[2:3,3:4,:] == cars[2:3,3:4,:]
@test all(i->getindex_count(i) == 7,dars)
v = view(ds,:,:,2)
@test v isa DiskArrays.SubDiskArray{Float64,2}
@test v[:,:] == ars[2]
@test getindex_count(dars[2]) == 8
v2 = view(ds,:,:,[3,1])
@test v2[:,:,:] == cars[:,:,[3,1]][:,:,:]
@test all(i->getindex_count(i) == 8,dars)
end


text:
module reprompijl
using MPI
using ArgParse
using Printf
using Pkg
using Statistics

@enum Collective begin
    MPI_Allreduce
    MPI_Alltoall
    MPI_Bcast
    MPI_Scan
    Nanosleep
end

struct Args
    nrep::Int64
    calls::Array{Collective,1}
    sizes::Array{Int64,1}
    operation::MPI.Op
    verbose::Bool
    random::Bool
    check::Bool
    summary::Bool
end

const root = 0

function str_to_collective(str::SubString{String})::Collective
    if str == "MPI_Allreduce"
        return MPI_Allreduce
    elseif str == "MPI_Alltoall"
        return MPI_Alltoall
    elseif str == "MPI_Bcast"
        return MPI_Bcast
    elseif str == "MPI_Scan"
        return MPI_Scan
    else
        return Nanosleep
    end
end

function parse_parameters()::Args
    s = ArgParseSettings()
    @add_arg_table s begin
        "--calls-list", "-c"
            help = "list of comma-separated MPI calls to be benchmarked, e.g., --calls-list=MPI_Bcast,MPI_Allgather"
            required = true
            arg_type = String
        "--msizes-list", "-s"
            help = "list of comma-separated message sizes in Bytes, e.g., --msizes-list=10,1024"
            required = true
            arg_type = String
        "--nrep", "-n"
            help = "set number of experiment repetitions"
            arg_type = Int
            default = 10
        "--verbose", "-v"
            help = "increase verbosity level (print times measured for each process)"
            action = :store_true
        "--random", "-r"
            help = "use random data in collective calls (default is all zeros)"
            action = :store_true
        "--check", "-k"
            help = "check the correctness of the calculated result by recalculating on one core"
            action = :store_true
        "--summary"
            help = "only prints statistical data about results"
            action = :store_true
    end
    args = parse_args(ARGS, s)
    nrep = args["nrep"]
    calls = [str_to_collective(ss) for ss in  split(args["calls-list"], ",")]
    message_sizes = [parse(Int, ss) for ss in split(args["msizes-list"], ",")]
    verbose = args["verbose"]
    random = args["random"]
    check = args["check"]
    summary = args["summary"]
    Args(nrep, calls, message_sizes, MPI.BOR, verbose, random, check, summary)
end

function print_info(args::Args)
    println("#MPI calls:")
    for c in args.calls
        println("#\t", c)
    end
    println("#Message sizes:")
    for s in args.sizes
        println("#\t", s)
    end
    println("#@operation=", args.operation)
    println("#@datatype=", "UInt8")
    println("#@nrep=", args.nrep)
    println("#@root_proc=", 0)
    println("#@nprocs=", MPI.Comm_size(MPI.COMM_WORLD))
    println("#@verbose=", args.verbose)
    println("#@random=", args.random)
    println("#@check=", args.check)
    println("####")
    println("#@Julia Version=", VERSION)
    println("#@MPI.jl Version=", Pkg.installed()["MPI"])
    println("#@MPI Version=", strip(split(read(`ompi_info`, String), "\n")[2]))
end

function print_verbose(times::Array{Float64,1}, args::Args, call::Collective, msize::Int64)
    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)
    size = MPI.Comm_size(comm)

    # alloc space for times
    if rank == root
        all_times = Array{Float64,2}(undef, args.nrep, size)
    else
        all_times = nothing
    end

    # get times from other processes
    MPI.Gather!(times, all_times, args.nrep, root, comm)

    # print times
    if rank == root
        for proc_id in 1:size
            for i in 1:args.nrep
                Printf.@printf("%7d %50s %10d %12ld %14.10f\n", proc_id,
                    call, i, msize, all_times[i, proc_id])
            end
        end
    end
end

function print_simple(times::Array{Float64,1}, args::Args, call::Collective, size::Int64)
    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)

    # calculate max runtimes
    maxRuntimes = MPI.Reduce(times, max, root, comm)

    # print results
    if rank == root
        for i in 1:args.nrep
            Printf.@printf("%50s %10d %12ld %14.10f\n", call, i, size, maxRuntimes[i])
        end
    end
end

function print_summary(times::Array{Float64,1}, args::Args, call::Collective, size::Int64)
    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)
    size = MPI.Comm_size(comm)

    # alloc space for times
    if rank == root
        all_times = Array{Float64,2}(undef, args.nrep, size)
    else
        all_times = nothing
    end

    # get times from other processes
    MPI.Gather!(times, all_times, args.nrep, root, comm)

    # calculate max runtime
    maxRuntimes = MPI.Reduce(times, max, root, comm)

    # calculate min runtime
    minRuntimes = MPI.Reduce(times, min, root, comm)

    # print times
    if rank == root
        maxRuntime = maximum(maxRuntimes)

        minRuntime = minimum(minRuntimes)

        # calculate median runtime
        medianRuntime = median(all_times)

        meanRuntime = mean(all_times)

        Printf.@printf("%50s %10d %14.10f %14.10f %14.10f %14.10f\n",
                    call, args.nrep, meanRuntime, medianRuntime, minRuntime, maxRuntime)
    end
end

function print_results(times::Array{Float64,1}, args::Args, call::Collective, size::Int64)
    if args.verbose
        print_verbose(times, args, call, size)
    elseif args.summary
        print_summary(times, args, call, size)
    else
        print_simple(times, args, call, size)
    end
end

function check_allreduce(msize::Int64, send)
    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)
    size = MPI.Comm_size(comm)

    MPI.Barrier(comm)

    # alloc space for send buffers
    if rank == root
        all_send = zeros(UInt8, msize, size)
    else
        all_send = nothing
    end

    MPI.Gather!(send, all_send, msize, root, comm)

    result = zeros(UInt8, msize)
    # the only supported operation is MPI.BOR
    if rank == root
        for proc_id in 1:size
            for i in 1:msize
                result[i] = result[i] | all_send[i, proc_id]
            end
        end
    end

    return result
end

function check_scan(msize::Int64, send, recv)
    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)
    size = MPI.Comm_size(comm)

    MPI.Barrier(comm)

    # alloc space for send buffers
    if rank == root
        all_send = zeros(UInt8, msize, size)
        all_recv = zeros(UInt8, msize, size)
    else
        all_send = nothing
        all_recv = nothing
    end

    MPI.Gather!(send, all_send, msize, root, comm)
    MPI.Gather!(recv, all_recv, msize, root, comm)

    result = zeros(UInt8, msize, size)
    # the only supported operation is MPI.BOR
    all_correct = true
    if rank == root
        for i in 1:msize
            for proc_id in 1:size
                if proc_id == 1
                    left = all_send[i, 1]
                else
                    left = result[i, proc_id - 1]
                end
                result[i, proc_id] = left | all_send[i, proc_id]
                if result[i, proc_id] != all_recv[i, proc_id]
                    println("Got incorrect result at: [", i, ", ", proc_id, "]")
                    all_correct = false
                end
            end
        end
    end

    return all_correct
end

function bench(args::Args)
    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)
    size = MPI.Comm_size(comm)
    times = zeros(Float64, args.nrep)

    # print header
    if rank == root
        if args.verbose
            Printf.@printf("%7s %50s %10s %12s %14s\n",  "process", "test", "nrep", "count", "time")
        elseif args.summary
            Printf.@printf("%50s %10s %14s %14s %14s %14s\n",
                    "test", "nrep", "mean_sec", "median_sec", "min_sec", "max_sec")
        else
            Printf.@printf("%50s %10s %12s %14s\n", "test", "nrep", "count", "runtime_sec")
        end
    end

    for msize in args.sizes
        for call in args.calls
            # init sync
            MPI.Barrier(comm)

            if call == MPI_Allreduce
                if args.random
                    send = rand(UInt8, msize)
                    recv = rand(UInt8, msize)
                else
                    send = zeros(UInt8, msize)
                    recv = zeros(UInt8, msize)
                end

                if args.check
                    result = check_allreduce(msize, send)
                end

                for i in 1:args.nrep
                    # start sync
                    MPI.Barrier(comm)

                    times[i] = MPI.Wtime()
                    MPI.Allreduce!(send, recv, msize, args.operation, comm)
                    times[i] = MPI.Wtime() - times[i]
                end
            elseif call == MPI_Scan
                if args.random
                    send = rand(UInt8, msize)
                    recv = rand(UInt8, msize)
                else
                    send = zeros(UInt8, msize)
                    recv = zeros(UInt8, msize)
                end

                for i in 1:args.nrep
                    # start sync
                    MPI.Barrier(comm)

                    times[i] = MPI.Wtime()
                    MPI.Scan!(send, recv, msize, args.operation, comm)
                    times[i] = MPI.Wtime() - times[i]
                end
            elseif call == MPI_Alltoall
                if args.random
                    send = rand(UInt8, msize * size)
                    recv = rand(UInt8, msize * size)
                else
                    send = zeros(UInt8, msize * size)
                    recv = zeros(UInt8, msize * size)
                end
                for i in 1:args.nrep
                    # start sync
                    MPI.Barrier(comm)

                    times[i] = MPI.Wtime()
                    MPI.Alltoall!(send, recv, msize, comm)
                    times[i] = MPI.Wtime() - times[i]
                end
            elseif call == MPI_Bcast
                if args.random
                    buf = zeros(UInt8, msize)
                else
                    buf = zeros(UInt8, msize)
                end
                for i in 1:args.nrep
                    # start sync
                    MPI.Barrier(comm)

                    times[i] = MPI.Wtime()
                    MPI.Bcast!(buf, msize, root, comm)
                    times[i] = MPI.Wtime() - times[i]
                end
            end

            # print timing output
            print_results(times, args, call, msize)

            all_correct = true
            if rank == root && args.check
                if call == MPI_Allreduce
                    for i in 1:msize
                        if result[i] != recv[i]
                            println("Got an incorrect result for element ", i)
                            all_correct = false
                        end
                    end
                end
            end

            if args.check && call == MPI_Scan
                if !check_scan(msize, send, recv) && rank == root
                    println("Got an incorrect result")
                    all_correct = false
                end
            end

            if rank == root && args.check && all_correct
                println("Verified every element successfully")
            end
        end
    end
end

function main()
    args = parse_parameters()
    MPI.Init()
    rank = MPI.Comm_rank(MPI.COMM_WORLD)
    if rank == 0
        print_info(args)
    end
    bench(args)
    MPI.Finalize()
end

main()

end # module


text:
# begin enum cusparseStatus_t
typealias cusparseStatus_t Uint32
const CUSPARSE_STATUS_SUCCESS = (uint32)(0)
const CUSPARSE_STATUS_NOT_INITIALIZED = (uint32)(1)
const CUSPARSE_STATUS_ALLOC_FAILED = (uint32)(2)
const CUSPARSE_STATUS_INVALID_VALUE = (uint32)(3)
const CUSPARSE_STATUS_ARCH_MISMATCH = (uint32)(4)
const CUSPARSE_STATUS_MAPPING_ERROR = (uint32)(5)
const CUSPARSE_STATUS_EXECUTION_FAILED = (uint32)(6)
const CUSPARSE_STATUS_INTERNAL_ERROR = (uint32)(7)
const CUSPARSE_STATUS_MATRIX_TYPE_NOT_SUPPORTED = (uint32)(8)
const CUSPARSE_STATUS_ZERO_PIVOT = (uint32)(9)
# end enum cusparseStatus_t

typealias cusparseContext Void
typealias cusparseHandle_t Ptr{cusparseContext}
typealias cusparseMatDescr Void
typealias cusparseMatDescr_t Ptr{cusparseMatDescr}
typealias cusparseSolveAnalysisInfo Void
typealias cusparseSolveAnalysisInfo_t Ptr{cusparseSolveAnalysisInfo}
typealias csrsv2Info Void
typealias csrsv2Info_t Ptr{csrsv2Info}
typealias bsrsv2Info Void
typealias bsrsv2Info_t Ptr{bsrsv2Info}
typealias bsrsm2Info Void
typealias bsrsm2Info_t Ptr{bsrsm2Info}
typealias csric02Info Void
typealias csric02Info_t Ptr{csric02Info}
typealias bsric02Info Void
typealias bsric02Info_t Ptr{bsric02Info}
typealias csrilu02Info Void
typealias csrilu02Info_t Ptr{csrilu02Info}
typealias bsrilu02Info Void
typealias bsrilu02Info_t Ptr{bsrilu02Info}
typealias cusparseHybMat Void
typealias cusparseHybMat_t Ptr{cusparseHybMat}

# begin enum cusparsePointerMode_t
typealias cusparsePointerMode_t Uint32
const CUSPARSE_POINTER_MODE_HOST = (uint32)(0)
const CUSPARSE_POINTER_MODE_DEVICE = (uint32)(1)
# end enum cusparsePointerMode_t

# begin enum cusparseAction_t
typealias cusparseAction_t Uint32
const CUSPARSE_ACTION_SYMBOLIC = (uint32)(0)
const CUSPARSE_ACTION_NUMERIC = (uint32)(1)
# end enum cusparseAction_t

# begin enum cusparseMatrixType_t
typealias cusparseMatrixType_t Uint32
const CUSPARSE_MATRIX_TYPE_GENERAL = (uint32)(0)
const CUSPARSE_MATRIX_TYPE_SYMMETRIC = (uint32)(1)
const CUSPARSE_MATRIX_TYPE_HERMITIAN = (uint32)(2)
const CUSPARSE_MATRIX_TYPE_TRIANGULAR = (uint32)(3)
# end enum cusparseMatrixType_t

# begin enum cusparseFillMode_t
typealias cusparseFillMode_t Uint32
const CUSPARSE_FILL_MODE_LOWER = (uint32)(0)
const CUSPARSE_FILL_MODE_UPPER = (uint32)(1)
# end enum cusparseFillMode_t

# begin enum cusparseDiagType_t
typealias cusparseDiagType_t Uint32
const CUSPARSE_DIAG_TYPE_NON_UNIT = (uint32)(0)
const CUSPARSE_DIAG_TYPE_UNIT = (uint32)(1)
# end enum cusparseDiagType_t

# begin enum cusparseIndexBase_t
typealias cusparseIndexBase_t Uint32
const CUSPARSE_INDEX_BASE_ZERO = (uint32)(0)
const CUSPARSE_INDEX_BASE_ONE = (uint32)(1)
# end enum cusparseIndexBase_t

# begin enum cusparseOperation_t
typealias cusparseOperation_t Uint32
const CUSPARSE_OPERATION_NON_TRANSPOSE = (uint32)(0)
const CUSPARSE_OPERATION_TRANSPOSE = (uint32)(1)
const CUSPARSE_OPERATION_CONJUGATE_TRANSPOSE = (uint32)(2)
# end enum cusparseOperation_t

# begin enum cusparseDirection_t
typealias cusparseDirection_t Uint32
const CUSPARSE_DIRECTION_ROW = (uint32)(0)
const CUSPARSE_DIRECTION_COLUMN = (uint32)(1)
# end enum cusparseDirection_t

# begin enum cusparseHybPartition_t
typealias cusparseHybPartition_t Uint32
const CUSPARSE_HYB_PARTITION_AUTO = (uint32)(0)
const CUSPARSE_HYB_PARTITION_USER = (uint32)(1)
const CUSPARSE_HYB_PARTITION_MAX = (uint32)(2)
# end enum cusparseHybPartition_t

# begin enum cusparseSolvePolicy_t
typealias cusparseSolvePolicy_t Uint32
const CUSPARSE_SOLVE_POLICY_NO_LEVEL = (uint32)(0)
const CUSPARSE_SOLVE_POLICY_USE_LEVEL = (uint32)(1)
# end enum cusparseSolvePolicy_t


text:
cd(@__DIR__)
using Revise, DifferentialEvolutionMCMC, Random, Distributions
import DifferentialEvolutionMCMC: minimize!

Random.seed!(50514)

function sample_prior()
    return [rand(Uniform(-5, 5), 2)]
 end

 function rastrigin(data, x)
     A = 10.0
     n = length(x)
     y = A * n
     for  i in 1:n
         y +=  + x[i]^2 - A * cos(2 * π * x[i])
     end
     return y
 end

bounds = ((-5.0,5.0),)
names = (:x,)

model = DEModel(;
    sample_prior,
    loglike = rastrigin,
    data = nothing,
    names
)

de = DE(;
    sample_prior,
    bounds,
    Np = 6,
    n_groups = 1,
    update_particle! = minimize!,
    evaluate_fitness! = evaluate_fun!
)

n_iter = 10_000
particles = optimize(model, de, n_iter, progress=true);
results = get_optimal(de, model, particles)
println(results)


text:
module BenchSparseIntSet

using DataStructures
using BenchmarkTools
using Random

suite = BenchmarkGroup()

rand_setup =  (
	Random.seed!(1234);
	ids1 = rand(1:30000, 1000);
	ids2 = rand(1:30000, 1000);
)

function create_fill_packed(ids1)
	y = SparseIntSet()
	for i in ids1
		push!(y, i)
	end
end

suite["create_fill"] =
	@benchmarkable create_fill_packed(ids1) setup=rand_setup

suite["in while not in"] =
	@benchmarkable in(23, y) evals=1000 setup=(y = SparseIntSet();)
suite["in while in"] =
	@benchmarkable in(5199, y) evals=1000 setup=(y=SparseIntSet(); push!(y, 5199))

function pop_push(y)
	pop!(y, 5199)
	push!(y, 5199)
end

suite["pop push worst case"] = @benchmarkable pop_push(y) setup=(y=SparseIntSet(); push!(y, 5199))
suite["pop push"] = @benchmarkable pop_push(y) setup=(y=SparseIntSet(); push!(y, 5199); push!(y, 5200))

function iterate_one_bench(x)
	t = 0
	for i in x
		t += i
	end
	return t
end
function iterate_two_bench(x,y)
	t = 0
	for (ix, iy) in zip(x, y)
		t += ix + iy
	end
	return t
end
function iterate_two_exclude_one_bench(x,y,z)
	t = 0
	for (ix, iy) in zip(x, y, exclude=(z,))
		t += ix + iy
	end
	return t
end

x_y_z_setup = (
	Random.seed!(1234);
	x = SparseIntSet(rand(1:30000, 1000));
	y = SparseIntSet(rand(1:30000, 1000));
	z = SparseIntSet(rand(1:30000, 1000));
)

suite["iterate one"] =
	@benchmarkable iterate_one_bench(x) setup=x_y_z_setup

suite["iterate two"] =
	@benchmarkable iterate_two_bench(x,y) setup=x_y_z_setup

suite["iterate two exclude one"] =
	@benchmarkable iterate_two_exclude_one_bench(x,y,z) setup=x_y_z_setup

end  # module
BenchSparseIntSet.suite

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

module Delays

using Actors
import Actors: newLink

"""
```
@delayed expr
@delayed expr timeout
```
Test an expression `expr` and return it after it becomes true.
This has a default timeout of 1 second if unspecified and
`timeout` seconds if specified. It does polling every 0.1 second.
"""
macro delayed(expr, timeout)
    return quote
        timedwait($timeout) do
            try
                return $(esc(expr))
            catch
                return false
            end
        end
        $(esc(expr))
    end
end
macro delayed(expr)
    return :(@delayed $(esc(expr)) 1)
end

function Base.:(==)(l1::L, l2::L) where L <:Link
    return l1.chn == l2.chn && l1.pid == l2.pid && l1.mode == l2.mode
end
Base.:(==)(t1::R, t2::R) where R<:Base.RefValue{Task} = t1[] == t2[]

Base.copy(rt::Base.RefValue{Task}) = Ref(rt[])
function Base.copy(lk::Link)
    lk1 = newLink()
    lk1.chn  = lk.chn
    lk1.pid  = lk.pid
    lk1.mode = lk.mode
    return lk1
end

"""
    changed(var, timeout=1)

Return a mutable variable `var` after it has been changed. Wait for
the change at maximum `timeout` seconds. Poll every 0.1 seconds.
"""
function changed(var, timeout=1)
    oldvar = copy(var)
    timedwait(timeout) do
        oldvar != var
    end
    return var
end

export @delayed, changed

end # Delays


text:
test_orthogonality_orthonormality(dict::Dictionary, measure::BasisFunctions.Measure=measure(dict)) =
    test_orthogonality_orthonormality(dict, isorthogonal(dict, measure), isorthonormal(dict, measure), measure)

function test_orthogonality_orthonormality(dict::Dictionary, isorthog::Bool, isorthonorm::Bool, measure::BasisFunctions.Measure=measure(dict); atol=1e-6,rtol=1e-6, opts...)
    isorthonormal(dict, measure) && (@test isorthogonal(dict, measure))
    @test isorthogonal(dict, measure) == isorthog
    @test isorthonormal(dict, measure) == isorthonorm

    G = gram(dict, measure)
    @test eltype(G) == coefficienttype(dict)
    M = Matrix(BasisFunctions.default_gram(dict, measure; warnslow=false, atol=atol,rtol=rtol,opts...))

    if isorthog
        @test isdiag(G)
    else
        @test !isdiag(G)
        @test norm(M - Diagonal(diag(M))) > 1e-3 || any(abs.(diag(M)).< 1e-5)
    end
    if isorthonorm
        @test BasisFunctions.isidentity(G)
    else
        @test !BasisFunctions.isidentity(G)
        @test !(diag(M) ≈ Ones(length(dict)))
    end

    @test norm(M - Matrix(G)) < 1e-5
end


text:
# create a new "build_series_args" which converts all inputs into xs = Any[xitems], ys = Any[yitems].
# Special handling for: no args, xmin/xmax, parametric, dataframes
# Then once inputs have been converted, build the series args, map functions, etc.
# This should cut down on boilerplate code and allow more focused dispatch on type
# note: returns meta information... mainly for use with automatic labeling from DataFrames for now

const FuncOrFuncs{F} = Union{F, Vector{F}, Matrix{F}}
const MaybeNumber = Union{Number, Missing}
const MaybeString = Union{AbstractString, Missing}
const DataPoint = Union{MaybeNumber, MaybeString}

prepareSeriesData(x) = error("Cannot convert $(typeof(x)) to series data for plotting")
prepareSeriesData(::Nothing) = nothing
prepareSeriesData(f::Function) = f
prepareSeriesData(a::AbstractArray{<:MaybeNumber}) = replace!(
                                    x -> ismissing(x) || isinf(x) ? NaN : x,
                                    map(float,a))
prepareSeriesData(a::AbstractArray{<:MaybeString}) = replace(x -> ismissing(x) ? "" : x, a)
prepareSeriesData(s::Surface{<:AMat{<:MaybeNumber}}) = Surface(prepareSeriesData(s.surf))
prepareSeriesData(s::Surface) = s  # non-numeric Surface, such as an image
prepareSeriesData(v::Volume) = Volume(prepareSeriesData(v.v), v.x_extents, v.y_extents, v.z_extents)

# default: assume x represents a single series
convertToAnyVector(x, plotattributes) = Any[prepareSeriesData(x)]

# fixed number of blank series
convertToAnyVector(n::Integer, plotattributes) = Any[zeros(0) for i in 1:n]

# vector of data points is a single series
convertToAnyVector(v::AVec{<:DataPoint}, plotattributes) = Any[prepareSeriesData(v)]

# list of things (maybe other vectors, functions, or something else)
function convertToAnyVector(v::AVec, plotattributes)
    if all(x -> x isa MaybeNumber, v)
        convertToAnyVector(Vector{MaybeNumber}(v), plotattributes)
    elseif all(x -> x isa MaybeString, v)
        convertToAnyVector(Vector{MaybeString}(v), plotattributes)
    else
        vcat((convertToAnyVector(vi, plotattributes) for vi in v)...)
    end
end

# Matrix is split into columns
function convertToAnyVector(v::AMat{<:DataPoint}, plotattributes)
    if all3D(plotattributes)
        Any[prepareSeriesData(Surface(v))]
    else
        Any[prepareSeriesData(v[:, i]) for i in 1:size(v, 2)]
    end
end

# --------------------------------------------------------------------
# Fillranges & ribbons

process_fillrange(range::Number, plotattributes) = [range]
process_fillrange(range, plotattributes) = convertToAnyVector(range, plotattributes)

process_ribbon(ribbon::Number, plotattributes) = [ribbon]
process_ribbon(ribbon, plotattributes) = convertToAnyVector(ribbon, plotattributes)
# ribbon as a tuple: (lower_ribbons, upper_ribbons)
process_ribbon(ribbon::Tuple{Any,Any}, plotattributes) = collect(zip(convertToAnyVector(ribbon[1], plotattributes),
                                                     convertToAnyVector(ribbon[2], plotattributes)))

# --------------------------------------------------------------------

# TODO: can we avoid the copy here?  one error that crops up is that mapping functions over the same array
#       result in that array being shared.  push!, etc will add too many items to that array

compute_x(x::Nothing, y::Nothing, z)      = 1:size(z,1)
compute_x(x::Nothing, y, z)            = 1:size(y,1)
compute_x(x::Function, y, z)        = map(x, y)
compute_x(x, y, z)                  = copy(x)

# compute_y(x::Void, y::Function, z)  = error()
compute_y(x::Nothing, y::Nothing, z)      = 1:size(z,2)
compute_y(x, y::Function, z)        = map(y, x)
compute_y(x, y, z)                  = copy(y)

compute_z(x, y, z::Function)        = map(z, x, y)
compute_z(x, y, z::AbstractMatrix)  = Surface(z)
compute_z(x, y, z::Nothing)            = nothing
compute_z(x, y, z)                  = copy(z)

nobigs(v::AVec{BigFloat}) = map(Float64, v)
nobigs(v::AVec{BigInt}) = map(Int64, v)
nobigs(v) = v

@noinline function compute_xyz(x, y, z)
    x = compute_x(x,y,z)
    y = compute_y(x,y,z)
    z = compute_z(x,y,z)
    nobigs(x), nobigs(y), nobigs(z)
end

# not allowed
compute_xyz(x::Nothing, y::FuncOrFuncs{F}, z) where {F<:Function}       = error("If you want to plot the function `$y`, you need to define the x values!")
compute_xyz(x::Nothing, y::Nothing, z::FuncOrFuncs{F}) where {F<:Function} = error("If you want to plot the function `$z`, you need to define x and y values!")
compute_xyz(x::Nothing, y::Nothing, z::Nothing)        = error("x/y/z are all nothing!")

# --------------------------------------------------------------------

# we are going to build recipes to do the processing and splitting of the args

# ensure we dispatch to the slicer
struct SliceIt end

# the catch-all recipes
@recipe function f(::Type{SliceIt}, x, y, z)

    # handle data with formatting attached
    if typeof(x) <: Formatted
        xformatter := x.formatter
        x = x.data
    end
    if typeof(y) <: Formatted
        yformatter := y.formatter
        y = y.data
    end
    if typeof(z) <: Formatted
        zformatter := z.formatter
        z = z.data
    end

    xs = convertToAnyVector(x, plotattributes)
    ys = convertToAnyVector(y, plotattributes)
    zs = convertToAnyVector(z, plotattributes)

    fr = pop!(plotattributes, :fillrange, nothing)
    fillranges = process_fillrange(fr, plotattributes)
    mf = length(fillranges)

    rib = pop!(plotattributes, :ribbon, nothing)
    ribbons = process_ribbon(rib, plotattributes)
    mr = length(ribbons)

    # @show zs

    mx = length(xs)
    my = length(ys)
    mz = length(zs)
    if mx > 0 && my > 0 && mz > 0
        for i in 1:max(mx, my, mz)
            # add a new series
            di = copy(plotattributes)
            xi, yi, zi = xs[mod1(i,mx)], ys[mod1(i,my)], zs[mod1(i,mz)]
            di[:x], di[:y], di[:z] = compute_xyz(xi, yi, zi)

            # handle fillrange
            fr = fillranges[mod1(i,mf)]
            di[:fillrange] = isa(fr, Function) ? map(fr, di[:x]) : fr

            # handle ribbons
            rib = ribbons[mod1(i,mr)]
            di[:ribbon] = isa(rib, Function) ? map(rib, di[:x]) : rib

            push!(series_list, RecipeData(di, ()))
        end
    end
    nothing  # don't add a series for the main block
end

# this is the default "type recipe"... just pass the object through
@recipe f(::Type{T}, v::T) where {T<:Any} = v

# this should catch unhandled "series recipes" and error with a nice message
@recipe f(::Type{V}, x, y, z) where {V<:Val} = error("The backend must not support the series type $V, and there isn't a series recipe defined.")

_apply_type_recipe(plotattributes, v) = RecipesBase.apply_recipe(plotattributes, typeof(v), v)[1].args[1]

# Handle type recipes when the recipe is defined on the elements.
# This sort of recipe should return a pair of functions... one to convert to number,
# and one to format tick values.
function _apply_type_recipe(plotattributes, v::AbstractArray)
    isempty(skipmissing(v)) && return Float64[]
    x = first(skipmissing(v))
    args = RecipesBase.apply_recipe(plotattributes, typeof(x), x)[1].args
    if length(args) == 2 && typeof(args[1]) <: Function && typeof(args[2]) <: Function
        numfunc, formatter = args
        Formatted(map(numfunc, v), formatter)
    else
        v
    end
end

# # special handling for Surface... need to properly unwrap and re-wrap
# function _apply_type_recipe(plotattributes, v::Surface)
#     T = eltype(v.surf)
#     @show T
#     if T <: Integer || T <: AbstractFloat
#         v
#     else
#         ret = _apply_type_recipe(plotattributes, v.surf)
#         if typeof(ret) <: Formatted
#             Formatted(Surface(ret.data), ret.formatter)
#         else
#             v
#         end
#     end
# end

# don't do anything for ints or floats
_apply_type_recipe(plotattributes, v::AbstractArray{T}) where {T<:Union{Integer,AbstractFloat}} = v

# handle "type recipes" by converting inputs, and then either re-calling or slicing
@recipe function f(x, y, z)
    did_replace = false
    newx = _apply_type_recipe(plotattributes, x)
    x === newx || (did_replace = true)
    newy = _apply_type_recipe(plotattributes, y)
    y === newy || (did_replace = true)
    newz = _apply_type_recipe(plotattributes, z)
    z === newz || (did_replace = true)
    if did_replace
        newx, newy, newz
    else
        SliceIt, x, y, z
    end
end
@recipe function f(x, y)
    did_replace = false
    newx = _apply_type_recipe(plotattributes, x)
    x === newx || (did_replace = true)
    newy = _apply_type_recipe(plotattributes, y)
    y === newy || (did_replace = true)
    if did_replace
        newx, newy
    else
        SliceIt, x, y, nothing
    end
end
@recipe function f(y)
    newy = _apply_type_recipe(plotattributes, y)
    if y !== newy
        newy
    else
        SliceIt, nothing, y, nothing
    end
end

# if there's more than 3 inputs, it can't be passed directly to SliceIt
# so we'll apply_type_recipe to all of them
@recipe function f(v1, v2, v3, v4, vrest...)
    did_replace = false
    newargs = map(v -> begin
        newv = _apply_type_recipe(plotattributes, v)
        if newv !== v
            did_replace = true
        end
        newv
    end, (v1, v2, v3, v4, vrest...))
    if !did_replace
        error("Couldn't process recipe args: $(map(typeof, (v1, v2, v3, v4, vrest...)))")
    end
    newargs
end

# # --------------------------------------------------------------------
# # 1 argument
# # --------------------------------------------------------------------

# helper function to ensure relevant attributes are wrapped by Surface
function wrap_surfaces(plotattributes::KW)
    if haskey(plotattributes, :fill_z)
        v = plotattributes[:fill_z]
        if !isa(v, Surface)
            plotattributes[:fill_z] = Surface(v)
        end
    end
end

@recipe f(n::Integer) = is3d(get(plotattributes,:seriestype,:path)) ? (SliceIt, n, n, n) : (SliceIt, n, n, nothing)

all3D(plotattributes::KW) = trueOrAllTrue(st -> st in (:contour, :contourf, :heatmap, :surface, :wireframe, :contour3d, :image, :plots_heatmap), get(plotattributes, :seriestype, :none))

# return a surface if this is a 3d plot, otherwise let it be sliced up
@recipe function f(mat::AMat{T}) where T<:Union{Integer,AbstractFloat,Missing}
    if all3D(plotattributes)
        n,m = size(mat)
        wrap_surfaces(plotattributes)
        SliceIt, 1:m, 1:n, Surface(mat)
    else
        SliceIt, nothing, mat, nothing
    end
end

# if a matrix is wrapped by Formatted, do similar logic, but wrap data with Surface
@recipe function f(fmt::Formatted{T}) where T<:AbstractMatrix
    if all3D(plotattributes)
        mat = fmt.data
        n,m = size(mat)
        wrap_surfaces(plotattributes)
        SliceIt, 1:m, 1:n, Formatted(Surface(mat), fmt.formatter)
    else
        SliceIt, nothing, fmt, nothing
    end
end

# assume this is a Volume, so construct one
@recipe function f(vol::AbstractArray{T,3}, args...) where T<:Union{Number,Missing}
    seriestype := :volume
    SliceIt, nothing, Volume(vol, args...), nothing
end

# # images - grays
function clamp_greys!(mat::AMat{T}) where T<:Gray
    for i in eachindex(mat)
        mat[i].val < 0 && (mat[i] = Gray(0))
        mat[i].val > 1 && (mat[i] = Gray(1))
    end
    mat
end

@recipe function f(mat::AMat{T}) where T<:Gray
    n, m = size(mat)
    if is_seriestype_supported(:image)
        seriestype := :image
        yflip --> true
        SliceIt, 1:m, 1:n, Surface(clamp_greys!(mat))
    else
        seriestype := :heatmap
        yflip --> true
        cbar --> false
        fillcolor --> ColorGradient([:black, :white])
        SliceIt, 1:m, 1:n, Surface(clamp!(convert(Matrix{Float64}, mat), 0., 1.))
    end
end

# # images - colors

@recipe function f(mat::AMat{T}) where T<:Colorant
	n, m = size(mat)

    if is_seriestype_supported(:image)
        seriestype := :image
        yflip --> true
        SliceIt, 1:m, 1:n, Surface(mat)
    else
        seriestype := :heatmap
        yflip --> true
        cbar --> false
        z, plotattributes[:fillcolor] = replace_image_with_heatmap(mat)
        SliceIt, 1:m, 1:n, Surface(z)
    end
end

#
# # plotting arbitrary shapes/polygons

@recipe function f(shape::Shape)
    seriestype --> :shape
    coords(shape)
end

@recipe function f(shapes::AVec{Shape})
    seriestype --> :shape
    coords(shapes)
end

@recipe function f(shapes::AMat{Shape})
    seriestype --> :shape
    for j in 1:size(shapes,2)
        @series coords(vec(shapes[:,j]))
    end
end

# Dicts: each entry is a data point (x,y)=(key,value)

@recipe f(d::AbstractDict) = collect(keys(d)), collect(values(d))

# function without range... use the current range of the x-axis

@recipe function f(f::FuncOrFuncs{F}) where F<:Function
    plt = plotattributes[:plot_object]
    xmin, xmax = if haskey(plotattributes, :xlims)
        plotattributes[:xlims]
    else
        try
            axis_limits(plt[1], :x)
        catch
            xinv = invscalefunc(get(plotattributes, :xscale, :identity))
            xm = tryrange(f, xinv.([-5,-1,0,0.01]))
            xm, tryrange(f, filter(x->x>xm, xinv.([5,1,0.99, 0, -0.01])))
        end
    end

    f, xmin, xmax
end

# try some intervals over which the function may be defined
function tryrange(F::AbstractArray, vec)
    rets = [tryrange(f, vec) for f in F] # get the preferred for each
    maxind = maximum(indexin(rets, vec)) # get the last attempt that succeeded (most likely to fit all)
    rets .= [tryrange(f, vec[maxind:maxind]) for f in F] # ensure that all functions compute there
    rets[1]
end

function tryrange(F, vec)
    for v in vec
        try
            tmp = F(v)
            return v
        catch
        end
    end
    error("$F is not a Function, or is not defined at any of the values $vec")
end
#
# # --------------------------------------------------------------------
# # 2 arguments
# # --------------------------------------------------------------------
#
#
# # if functions come first, just swap the order (not to be confused with parametric functions...
# # as there would be more than one function passed in)

@recipe function f(f::FuncOrFuncs{F}, x) where F<:Function
    F2 = typeof(x)
    @assert !(F2 <: Function || (F2 <: AbstractArray && F2.parameters[1] <: Function))  # otherwise we'd hit infinite recursion here
    x, f
end

#
# # --------------------------------------------------------------------
# # 3 arguments
# # --------------------------------------------------------------------
#
#
# # 3d line or scatter

@recipe function f(x::AVec, y::AVec, z::AVec)
    # st = get(plotattributes, :seriestype, :none)
    # if st == :scatter
    #     plotattributes[:seriestype] = :scatter3d
    # elseif !is3d(st)
    #     plotattributes[:seriestype] = :path3d
    # end
    SliceIt, x, y, z
end

@recipe function f(x::AMat, y::AMat, z::AMat)
    # st = get(plotattributes, :seriestype, :none)
    # if size(x) == size(y) == size(z)
    #     if !is3d(st)
    #         seriestype := :path3d
    #     end
    # end
    wrap_surfaces(plotattributes)
    SliceIt, x, y, z
end

#
# # surface-like... function

@recipe function f(x::AVec, y::AVec, zf::Function)
    # x = X <: Number ? sort(x) : x
    # y = Y <: Number ? sort(y) : y
    wrap_surfaces(plotattributes)
    SliceIt, x, y, Surface(zf, x, y)  # TODO: replace with SurfaceFunction when supported
end

#
# # surface-like... matrix grid

@recipe function f(x::AVec, y::AVec, z::AMat)
    if !like_surface(get(plotattributes, :seriestype, :none))
        plotattributes[:seriestype] = :contour
    end
    wrap_surfaces(plotattributes)
    SliceIt, x, y, Surface(z)
end

# # images - grays

@recipe function f(x::AVec, y::AVec, mat::AMat{T}) where T<:Gray
    if is_seriestype_supported(:image)
        seriestype := :image
        yflip --> true
        SliceIt, x, y, Surface(mat)
    else
        seriestype := :heatmap
        yflip --> true
        cbar --> false
        fillcolor --> ColorGradient([:black, :white])
        SliceIt, x, y, Surface(convert(Matrix{Float64}, mat))
    end
end

# # images - colors

@recipe function f(x::AVec, y::AVec, mat::AMat{T}) where T<:Colorant
    if is_seriestype_supported(:image)
        seriestype := :image
        yflip --> true
        SliceIt, x, y, Surface(mat)
    else
        seriestype := :heatmap
        yflip --> true
        cbar --> false
        z, plotattributes[:fillcolor] = replace_image_with_heatmap(mat)
        SliceIt, x, y, Surface(z)
    end
end

#
#
# # --------------------------------------------------------------------
# # Parametric functions
# # --------------------------------------------------------------------

#
# # special handling... xmin/xmax with parametric function(s)
@recipe function f(f::Function, xmin::Number, xmax::Number)
    xscale, yscale = [get(plotattributes, sym, :identity) for sym=(:xscale,:yscale)]
    _scaled_adapted_grid(f, xscale, yscale, xmin, xmax)
end
@recipe function f(fs::AbstractArray{F}, xmin::Number, xmax::Number) where F<:Function
    xscale, yscale = [get(plotattributes, sym, :identity) for sym=(:xscale,:yscale)]
    xs = Array{Any}(undef, length(fs))
    ys = Array{Any}(undef, length(fs))
    for (i, (x, y)) in enumerate(_scaled_adapted_grid(f, xscale, yscale, xmin, xmax) for f in fs)
        xs[i] = x
    	ys[i] = y
    end
    xs, ys
end
@recipe f(fx::FuncOrFuncs{F}, fy::FuncOrFuncs{G}, u::AVec) where {F<:Function,G<:Function}  = mapFuncOrFuncs(fx, u), mapFuncOrFuncs(fy, u)
@recipe f(fx::FuncOrFuncs{F}, fy::FuncOrFuncs{G}, umin::Number, umax::Number, n = 200) where {F<:Function,G<:Function} = fx, fy, range(umin, stop = umax, length = n)

function _scaled_adapted_grid(f, xscale, yscale, xmin, xmax)
    (xf, xinv), (yf, yinv) =  ((scalefunc(s),invscalefunc(s)) for s in (xscale,yscale))
    xs, ys = adapted_grid(yf∘f∘xinv, xf.((xmin, xmax)))
    xinv.(xs), yinv.(ys)
end

#
# # special handling... 3D parametric function(s)
@recipe function f(fx::FuncOrFuncs{F}, fy::FuncOrFuncs{G}, fz::FuncOrFuncs{H}, u::AVec) where {F<:Function,G<:Function,H<:Function}
    mapFuncOrFuncs(fx, u), mapFuncOrFuncs(fy, u), mapFuncOrFuncs(fz, u)
end
@recipe function f(fx::FuncOrFuncs{F}, fy::FuncOrFuncs{G}, fz::FuncOrFuncs{H}, umin::Number, umax::Number, numPoints = 200) where {F<:Function,G<:Function,H<:Function}
    fx, fy, fz, range(umin, stop = umax, length = numPoints)
end

#
#
# # --------------------------------------------------------------------
# # Lists of tuples and GeometryTypes.Points
# # --------------------------------------------------------------------
#

@recipe f(v::AVec{<:Tuple})               = unzip(v)
@recipe f(v::AVec{<:GeometryTypes.Point}) = unzip(v)
@recipe f(tup::Tuple)             = [tup]
@recipe f(p::GeometryTypes.Point) = [p]

# Special case for 4-tuples in :ohlc series
@recipe f(xyuv::AVec{<:Tuple{R1,R2,R3,R4}}) where {R1,R2,R3,R4} = get(plotattributes,:seriestype,:path)==:ohlc ? OHLC[OHLC(t...) for t in xyuv] : unzip(xyuv)

#
# # --------------------------------------------------------------------
# # handle grouping
# # --------------------------------------------------------------------

# @recipe function f(groupby::GroupBy, args...)
#     for (i,glab) in enumerate(groupby.groupLabels)
#         # create a new series, with the label of the group, and an idxfilter (to be applied in slice_and_dice)
#         # TODO: use @series instead
#         @show i, glab, groupby.groupIds[i]
#         di = copy(plotattributes)
#         get!(di, :label, string(glab))
#         get!(di, :idxfilter, groupby.groupIds[i])
#         push!(series_list, RecipeData(di, args))
#     end
#     nothing
# end

splittable_kw(key, val, lengthGroup) = false
splittable_kw(key, val::AbstractArray, lengthGroup) = !(key in (:group, :color_palette)) && size(val,1) == lengthGroup
splittable_kw(key, val::Tuple, lengthGroup) = all(splittable_kw.(key, val, lengthGroup))
splittable_kw(key, val::SeriesAnnotations, lengthGroup) = splittable_kw(key, val.strs, lengthGroup)

split_kw(key, val::AbstractArray, indices) = val[indices, fill(Colon(), ndims(val)-1)...]
split_kw(key, val::Tuple, indices) = Tuple(split_kw(key, v, indices) for v in val)
function split_kw(key, val::SeriesAnnotations, indices)
    split_strs = split_kw(key, val.strs, indices)
    return SeriesAnnotations(split_strs, val.font, val.baseshape, val.scalefactor)
end

function groupedvec2mat(x_ind, x, y::AbstractArray, groupby, def_val = y[1])
    y_mat = Array{promote_type(eltype(y), typeof(def_val))}(undef, length(keys(x_ind)), length(groupby.groupLabels))
    fill!(y_mat, def_val)
    for i in 1:length(groupby.groupLabels)
        xi = x[groupby.groupIds[i]]
        yi = y[groupby.groupIds[i]]
        y_mat[getindex.(Ref(x_ind), xi), i] = yi
    end
    return y_mat
end

groupedvec2mat(x_ind, x, y::Tuple, groupby) = Tuple(groupedvec2mat(x_ind, x, v, groupby) for v in y)

group_as_matrix(t) = false

# split the group into 1 series per group, and set the label and idxfilter for each
@recipe function f(groupby::GroupBy, args...)
    lengthGroup = maximum(union(groupby.groupIds...))
    if !(group_as_matrix(args[1]))
        for (i,glab) in enumerate(groupby.groupLabels)
            @series begin
                label     --> string(glab)
                idxfilter --> groupby.groupIds[i]
                for (key,val) in plotattributes
                    if splittable_kw(key, val, lengthGroup)
                        :($key) := split_kw(key, val, groupby.groupIds[i])
                    end
                end
                args
            end
        end
    else
        g = args[1]
        if length(g.args) == 1
            x = zeros(Int, lengthGroup)
            for indexes in groupby.groupIds
                x[indexes] = 1:length(indexes)
            end
            last_args = g.args
        else
            x = g.args[1]
            last_args = g.args[2:end]
        end
        x_u = unique(sort(x))
        x_ind = Dict(zip(x_u, 1:length(x_u)))
        for (key,val) in plotattributes
            if splittable_kw(key, val, lengthGroup)
                :($key) := groupedvec2mat(x_ind, x, val, groupby)
            end
        end
        label --> reshape(groupby.groupLabels, 1, :)
        typeof(g)((x_u, (groupedvec2mat(x_ind, x, arg, groupby, NaN) for arg in last_args)...))
    end
end


text:
##################################################################################################
function nodal_expression(canonical_model::CanonicalModel,
                         devices,
                         system_formulation::Type{S}) where {S<:PM.AbstractPowerModel}

    if model_has_parameters(canonical_model)
        _nodal_expression_param(canonical_model, devices, system_formulation)
    else
        _nodal_expression_fixed(canonical_model, devices, system_formulation)
    end

    return

end


text:
function ConcurrentCollections.maybepopfirst!(ch::Channel)
    y = iterate(ch)
    y === nothing && return nothing
    return Some(first(y))
end