text:
const cudaHostAllocDefault = 0x00
const cudaHostAllocPortable = 0x01
const cudaHostAllocMapped = 0x02
const cudaHostAllocWriteCombined = 0x04
const cudaHostRegisterDefault = 0x00
const cudaHostRegisterPortable = 0x01
const cudaHostRegisterMapped = 0x02
const cudaPeerAccessDefault = 0x00
const cudaStreamDefault = 0x00
const cudaStreamNonBlocking = 0x01
const cudaEventDefault = 0x00
const cudaEventBlockingSync = 0x01
const cudaEventDisableTiming = 0x02
const cudaEventInterprocess = 0x04
const cudaDeviceScheduleAuto = 0x00
const cudaDeviceScheduleSpin = 0x01
const cudaDeviceScheduleYield = 0x02
const cudaDeviceScheduleBlockingSync = 0x04
const cudaDeviceBlockingSync = 0x04
const cudaDeviceScheduleMask = 0x07
const cudaDeviceMapHost = 0x08
const cudaDeviceLmemResizeToMax = 0x10
const cudaDeviceMask = 0x01
const cudaArrayDefault = 0x00
const cudaArrayLayered = 0x01
const cudaArraySurfaceLoadStore = 0x02
const cudaArrayCubemap = 0x04
const cudaArrayTextureGather = 0x08
const cudaIpcMemLazyEnablePeerAccess = 0x01
const cudaMemAttachGlobal = 0x01
const cudaMemAttachHost = 0x02
const cudaMemAttachSingle = 0x04

# Skipping MacroDefinition: cudaDevicePropDontCare { { '\0' } , /* char   name[256];               */ 0 , /* size_t totalGlobalMem;          */ 0 , /* size_t sharedMemPerBlock;       */ 0 , /* int    regsPerBlock;            */ 0 , /* int    warpSize;                */ 0 , /* size_t memPitch;                */ 0 , /* int    maxThreadsPerBlock;      */ { 0 , 0 , 0 } , /* int    maxThreadsDim[3];        */ { 0 , 0 , 0 } , /* int    maxGridSize[3];          */ 0 , /* int    clockRate;               */ 0 , /* size_t totalConstMem;           */ - 1 , /* int    major;                   */ - 1 , /* int    minor;                   */ 0 , /* size_t textureAlignment;        */ 0 , /* size_t texturePitchAlignment    */ - 1 , /* int    deviceOverlap;           */ 0 , /* int    multiProcessorCount;     */ 0 , /* int    kernelExecTimeoutEnabled */ 0 , /* int    integrated               */ 0 , /* int    canMapHostMemory         */ 0 , /* int    computeMode              */ 0 , /* int    maxTexture1D             */ 0 , /* int    maxTexture1DMipmap       */ 0 , /* int    maxTexture1DLinear       */ { 0 , 0 } , /* int    maxTexture2D[2]          */ { 0 , 0 } , /* int    maxTexture2DMipmap[2]    */ { 0 , 0 , 0 } , /* int    maxTexture2DLinear[3]    */ { 0 , 0 } , /* int    maxTexture2DGather[2]    */ { 0 , 0 , 0 } , /* int    maxTexture3D[3]          */ { 0 , 0 , 0 } , /* int    maxTexture3DAlt[3]       */ 0 , /* int    maxTextureCubemap        */ { 0 , 0 } , /* int    maxTexture1DLayered[2]   */ { 0 , 0 , 0 } , /* int    maxTexture2DLayered[3]   */ { 0 , 0 } , /* int    maxTextureCubemapLayered[2] */ 0 , /* int    maxSurface1D             */ { 0 , 0 } , /* int    maxSurface2D[2]          */ { 0 , 0 , 0 } , /* int    maxSurface3D[3]          */ { 0 , 0 } , /* int    maxSurface1DLayered[2]   */ { 0 , 0 , 0 } , /* int    maxSurface2DLayered[3]   */ 0 , /* int    maxSurfaceCubemap        */ { 0 , 0 } , /* int    maxSurfaceCubemapLayered[2] */ 0 , /* size_t surfaceAlignment         */ 0 , /* int    concurrentKernels        */ 0 , /* int    ECCEnabled               */ 0 , /* int    pciBusID                 */ 0 , /* int    pciDeviceID              */ 0 , /* int    pciDomainID              */ 0 , /* int    tccDriver                */ 0 , /* int    asyncEngineCount         */ 0 , /* int    unifiedAddressing        */ 0 , /* int    memoryClockRate          */ 0 , /* int    memoryBusWidth           */ 0 , /* int    l2CacheSize              */ 0 , /* int    maxThreadsPerMultiProcessor */ 0 , /* int    streamPrioritiesSupported */ 0 , /* int    globalL1CacheSupported   */ 0 , /* int    localL1CacheSupported    */ 0 , /* size_t sharedMemPerMultiprocessor; */ 0 , /* int    regsPerMultiprocessor;   */ 0 , /* int    managedMemory            */ 0 , /* int    isMultiGpuBoard          */ 0 , /* int    multiGpuBoardGroupID     */ }

const CUDA_IPC_HANDLE_SIZE = 64

# begin enum cudaError
const cudaError = UInt32
const cudaSuccess = (UInt32)(0)
const cudaErrorMissingConfiguration = (UInt32)(1)
const cudaErrorMemoryAllocation = (UInt32)(2)
const cudaErrorInitializationError = (UInt32)(3)
const cudaErrorLaunchFailure = (UInt32)(4)
const cudaErrorPriorLaunchFailure = (UInt32)(5)
const cudaErrorLaunchTimeout = (UInt32)(6)
const cudaErrorLaunchOutOfResources = (UInt32)(7)
const cudaErrorInvalidDeviceFunction = (UInt32)(8)
const cudaErrorInvalidConfiguration = (UInt32)(9)
const cudaErrorInvalidDevice = (UInt32)(10)
const cudaErrorInvalidValue = (UInt32)(11)
const cudaErrorInvalidPitchValue = (UInt32)(12)
const cudaErrorInvalidSymbol = (UInt32)(13)
const cudaErrorMapBufferObjectFailed = (UInt32)(14)
const cudaErrorUnmapBufferObjectFailed = (UInt32)(15)
const cudaErrorInvalidHostPointer = (UInt32)(16)
const cudaErrorInvalidDevicePointer = (UInt32)(17)
const cudaErrorInvalidTexture = (UInt32)(18)
const cudaErrorInvalidTextureBinding = (UInt32)(19)
const cudaErrorInvalidChannelDescriptor = (UInt32)(20)
const cudaErrorInvalidMemcpyDirection = (UInt32)(21)
const cudaErrorAddressOfConstant = (UInt32)(22)
const cudaErrorTextureFetchFailed = (UInt32)(23)
const cudaErrorTextureNotBound = (UInt32)(24)
const cudaErrorSynchronizationError = (UInt32)(25)
const cudaErrorInvalidFilterSetting = (UInt32)(26)
const cudaErrorInvalidNormSetting = (UInt32)(27)
const cudaErrorMixedDeviceExecution = (UInt32)(28)
const cudaErrorCudartUnloading = (UInt32)(29)
const cudaErrorUnknown = (UInt32)(30)
const cudaErrorNotYetImplemented = (UInt32)(31)
const cudaErrorMemoryValueTooLarge = (UInt32)(32)
const cudaErrorInvalidResourceHandle = (UInt32)(33)
const cudaErrorNotReady = (UInt32)(34)
const cudaErrorInsufficientDriver = (UInt32)(35)
const cudaErrorSetOnActiveProcess = (UInt32)(36)
const cudaErrorInvalidSurface = (UInt32)(37)
const cudaErrorNoDevice = (UInt32)(38)
const cudaErrorECCUncorrectable = (UInt32)(39)
const cudaErrorSharedObjectSymbolNotFound = (UInt32)(40)
const cudaErrorSharedObjectInitFailed = (UInt32)(41)
const cudaErrorUnsupportedLimit = (UInt32)(42)
const cudaErrorDuplicateVariableName = (UInt32)(43)
const cudaErrorDuplicateTextureName = (UInt32)(44)
const cudaErrorDuplicateSurfaceName = (UInt32)(45)
const cudaErrorDevicesUnavailable = (UInt32)(46)
const cudaErrorInvalidKernelImage = (UInt32)(47)
const cudaErrorNoKernelImageForDevice = (UInt32)(48)
const cudaErrorIncompatibleDriverContext = (UInt32)(49)
const cudaErrorPeerAccessAlreadyEnabled = (UInt32)(50)
const cudaErrorPeerAccessNotEnabled = (UInt32)(51)
const cudaErrorDeviceAlreadyInUse = (UInt32)(54)
const cudaErrorProfilerDisabled = (UInt32)(55)
const cudaErrorProfilerNotInitialized = (UInt32)(56)
const cudaErrorProfilerAlreadyStarted = (UInt32)(57)
const cudaErrorProfilerAlreadyStopped = (UInt32)(58)
const cudaErrorAssert = (UInt32)(59)
const cudaErrorTooManyPeers = (UInt32)(60)
const cudaErrorHostMemoryAlreadyRegistered = (UInt32)(61)
const cudaErrorHostMemoryNotRegistered = (UInt32)(62)
const cudaErrorOperatingSystem = (UInt32)(63)
const cudaErrorPeerAccessUnsupported = (UInt32)(64)
const cudaErrorLaunchMaxDepthExceeded = (UInt32)(65)
const cudaErrorLaunchFileScopedTex = (UInt32)(66)
const cudaErrorLaunchFileScopedSurf = (UInt32)(67)
const cudaErrorSyncDepthExceeded = (UInt32)(68)
const cudaErrorLaunchPendingCountExceeded = (UInt32)(69)
const cudaErrorNotPermitted = (UInt32)(70)
const cudaErrorNotSupported = (UInt32)(71)
const cudaErrorHardwareStackError = (UInt32)(72)
const cudaErrorIllegalInstruction = (UInt32)(73)
const cudaErrorMisalignedAddress = (UInt32)(74)
const cudaErrorInvalidAddressSpace = (UInt32)(75)
const cudaErrorInvalidPc = (UInt32)(76)
const cudaErrorIllegalAddress = (UInt32)(77)
const cudaErrorInvalidPtx = (UInt32)(78)
const cudaErrorInvalidGraphicsContext = (UInt32)(79)
const cudaErrorStartupFailure = (UInt32)(127)
const cudaErrorApiFailureBase = (UInt32)(10000)
# end enum cudaError

# begin enum cudaChannelFormatKind
const cudaChannelFormatKind = UInt32
const cudaChannelFormatKindSigned = (UInt32)(0)
const cudaChannelFormatKindUnsigned = (UInt32)(1)
const cudaChannelFormatKindFloat = (UInt32)(2)
const cudaChannelFormatKindNone = (UInt32)(3)
# end enum cudaChannelFormatKind

immutable cudaChannelFormatDesc
    x::Cint
    y::Cint
    z::Cint
    w::Cint
    f::cudaChannelFormatKind
end

const cudaArray = Void
const cudaArray_t = Ptr{cudaArray}
const cudaArray_const_t = Ptr{cudaArray}
const cudaMipmappedArray = Void
const cudaMipmappedArray_t = Ptr{cudaMipmappedArray}
const cudaMipmappedArray_const_t = Ptr{cudaMipmappedArray}

# begin enum cudaMemoryType
const cudaMemoryType = UInt32
const cudaMemoryTypeHost = (UInt32)(1)
const cudaMemoryTypeDevice = (UInt32)(2)
# end enum cudaMemoryType

# begin enum cudaMemcpyKind
const cudaMemcpyKind = UInt32
const cudaMemcpyHostToHost = (UInt32)(0)
const cudaMemcpyHostToDevice = (UInt32)(1)
const cudaMemcpyDeviceToHost = (UInt32)(2)
const cudaMemcpyDeviceToDevice = (UInt32)(3)
const cudaMemcpyDefault = (UInt32)(4)
# end enum cudaMemcpyKind

immutable cudaPitchedPtr
    ptr::Ptr{Void}
    pitch::Csize_t
    xsize::Csize_t
    ysize::Csize_t
end

immutable cudaExtent
    width::Csize_t
    height::Csize_t
    depth::Csize_t
end

immutable cudaPos
    x::Csize_t
    y::Csize_t
    z::Csize_t
end

immutable cudaMemcpy3DParms
    srcArray::cudaArray_t
    srcPos::cudaPos
    srcPtr::cudaPitchedPtr
    dstArray::cudaArray_t
    dstPos::cudaPos
    dstPtr::cudaPitchedPtr
    extent::cudaExtent
    kind::cudaMemcpyKind
end

immutable cudaMemcpy3DPeerParms
    srcArray::cudaArray_t
    srcPos::cudaPos
    srcPtr::cudaPitchedPtr
    srcDevice::Cint
    dstArray::cudaArray_t
    dstPos::cudaPos
    dstPtr::cudaPitchedPtr
    dstDevice::Cint
    extent::cudaExtent
end

const cudaGraphicsResource = Void

# begin enum cudaGraphicsRegisterFlags
const cudaGraphicsRegisterFlags = UInt32
const cudaGraphicsRegisterFlagsNone = (UInt32)(0)
const cudaGraphicsRegisterFlagsReadOnly = (UInt32)(1)
const cudaGraphicsRegisterFlagsWriteDiscard = (UInt32)(2)
const cudaGraphicsRegisterFlagsSurfaceLoadStore = (UInt32)(4)
const cudaGraphicsRegisterFlagsTextureGather = (UInt32)(8)
# end enum cudaGraphicsRegisterFlags

# begin enum cudaGraphicsMapFlags
const cudaGraphicsMapFlags = UInt32
const cudaGraphicsMapFlagsNone = (UInt32)(0)
const cudaGraphicsMapFlagsReadOnly = (UInt32)(1)
const cudaGraphicsMapFlagsWriteDiscard = (UInt32)(2)
# end enum cudaGraphicsMapFlags

# begin enum cudaGraphicsCubeFace
const cudaGraphicsCubeFace = UInt32
const cudaGraphicsCubeFacePositiveX = (UInt32)(0)
const cudaGraphicsCubeFaceNegativeX = (UInt32)(1)
const cudaGraphicsCubeFacePositiveY = (UInt32)(2)
const cudaGraphicsCubeFaceNegativeY = (UInt32)(3)
const cudaGraphicsCubeFacePositiveZ = (UInt32)(4)
const cudaGraphicsCubeFaceNegativeZ = (UInt32)(5)
# end enum cudaGraphicsCubeFace

# begin enum cudaResourceType
const cudaResourceType = UInt32
const cudaResourceTypeArray = (UInt32)(0)
const cudaResourceTypeMipmappedArray = (UInt32)(1)
const cudaResourceTypeLinear = (UInt32)(2)
const cudaResourceTypePitch2D = (UInt32)(3)
# end enum cudaResourceType

# begin enum cudaResourceViewFormat
const cudaResourceViewFormat = UInt32
const cudaResViewFormatNone = (UInt32)(0)
const cudaResViewFormatUnsignedChar1 = (UInt32)(1)
const cudaResViewFormatUnsignedChar2 = (UInt32)(2)
const cudaResViewFormatUnsignedChar4 = (UInt32)(3)
const cudaResViewFormatSignedChar1 = (UInt32)(4)
const cudaResViewFormatSignedChar2 = (UInt32)(5)
const cudaResViewFormatSignedChar4 = (UInt32)(6)
const cudaResViewFormatUnsignedShort1 = (UInt32)(7)
const cudaResViewFormatUnsignedShort2 = (UInt32)(8)
const cudaResViewFormatUnsignedShort4 = (UInt32)(9)
const cudaResViewFormatSignedShort1 = (UInt32)(10)
const cudaResViewFormatSignedShort2 = (UInt32)(11)
const cudaResViewFormatSignedShort4 = (UInt32)(12)
const cudaResViewFormatUnsignedInt1 = (UInt32)(13)
const cudaResViewFormatUnsignedInt2 = (UInt32)(14)
const cudaResViewFormatUnsignedInt4 = (UInt32)(15)
const cudaResViewFormatSignedInt1 = (UInt32)(16)
const cudaResViewFormatSignedInt2 = (UInt32)(17)
const cudaResViewFormatSignedInt4 = (UInt32)(18)
const cudaResViewFormatHalf1 = (UInt32)(19)
const cudaResViewFormatHalf2 = (UInt32)(20)
const cudaResViewFormatHalf4 = (UInt32)(21)
const cudaResViewFormatFloat1 = (UInt32)(22)
const cudaResViewFormatFloat2 = (UInt32)(23)
const cudaResViewFormatFloat4 = (UInt32)(24)
const cudaResViewFormatUnsignedBlockCompressed1 = (UInt32)(25)
const cudaResViewFormatUnsignedBlockCompressed2 = (UInt32)(26)
const cudaResViewFormatUnsignedBlockCompressed3 = (UInt32)(27)
const cudaResViewFormatUnsignedBlockCompressed4 = (UInt32)(28)
const cudaResViewFormatSignedBlockCompressed4 = (UInt32)(29)
const cudaResViewFormatUnsignedBlockCompressed5 = (UInt32)(30)
const cudaResViewFormatSignedBlockCompressed5 = (UInt32)(31)
const cudaResViewFormatUnsignedBlockCompressed6H = (UInt32)(32)
const cudaResViewFormatSignedBlockCompressed6H = (UInt32)(33)
const cudaResViewFormatUnsignedBlockCompressed7 = (UInt32)(34)
# end enum cudaResourceViewFormat

immutable cudaResourceDesc
    resType::cudaResourceType
    res::Void
end

immutable cudaResourceViewDesc
    format::cudaResourceViewFormat
    width::Csize_t
    height::Csize_t
    depth::Csize_t
    firstMipmapLevel::UInt32
    lastMipmapLevel::UInt32
    firstLayer::UInt32
    lastLayer::UInt32
end

immutable cudaPointerAttributes
    memoryType::cudaMemoryType
    device::Cint
    devicePointer::Ptr{Void}
    hostPointer::Ptr{Void}
    isManaged::Cint
end

immutable cudaFuncAttributes
    sharedSizeBytes::Csize_t
    constSizeBytes::Csize_t
    localSizeBytes::Csize_t
    maxThreadsPerBlock::Cint
    numRegs::Cint
    ptxVersion::Cint
    binaryVersion::Cint
    cacheModeCA::Cint
end

# begin enum cudaFuncCache
const cudaFuncCache = UInt32
const cudaFuncCachePreferNone = (UInt32)(0)
const cudaFuncCachePreferShared = (UInt32)(1)
const cudaFuncCachePreferL1 = (UInt32)(2)
const cudaFuncCachePreferEqual = (UInt32)(3)
# end enum cudaFuncCache

# begin enum cudaSharedMemConfig
const cudaSharedMemConfig = UInt32
const cudaSharedMemBankSizeDefault = (UInt32)(0)
const cudaSharedMemBankSizeFourByte = (UInt32)(1)
const cudaSharedMemBankSizeEightByte = (UInt32)(2)
# end enum cudaSharedMemConfig

# begin enum cudaComputeMode
const cudaComputeMode = UInt32
const cudaComputeModeDefault = (UInt32)(0)
const cudaComputeModeExclusive = (UInt32)(1)
const cudaComputeModeProhibited = (UInt32)(2)
const cudaComputeModeExclusiveProcess = (UInt32)(3)
# end enum cudaComputeMode

# begin enum cudaLimit
const cudaLimit = UInt32
const cudaLimitStackSize = (UInt32)(0)
const cudaLimitPrintfFifoSize = (UInt32)(1)
const cudaLimitMallocHeapSize = (UInt32)(2)
const cudaLimitDevRuntimeSyncDepth = (UInt32)(3)
const cudaLimitDevRuntimePendingLaunchCount = (UInt32)(4)
# end enum cudaLimit

# begin enum cudaOutputMode
const cudaOutputMode = UInt32
const cudaKeyValuePair = (UInt32)(0)
const cudaCSV = (UInt32)(1)
# end enum cudaOutputMode

# begin enum cudaDeviceAttr
const cudaDeviceAttr = UInt32
const cudaDevAttrMaxThreadsPerBlock = (UInt32)(1)
const cudaDevAttrMaxBlockDimX = (UInt32)(2)
const cudaDevAttrMaxBlockDimY = (UInt32)(3)
const cudaDevAttrMaxBlockDimZ = (UInt32)(4)
const cudaDevAttrMaxGridDimX = (UInt32)(5)
const cudaDevAttrMaxGridDimY = (UInt32)(6)
const cudaDevAttrMaxGridDimZ = (UInt32)(7)
const cudaDevAttrMaxSharedMemoryPerBlock = (UInt32)(8)
const cudaDevAttrTotalConstantMemory = (UInt32)(9)
const cudaDevAttrWarpSize = (UInt32)(10)
const cudaDevAttrMaxPitch = (UInt32)(11)
const cudaDevAttrMaxRegistersPerBlock = (UInt32)(12)
const cudaDevAttrClockRate = (UInt32)(13)
const cudaDevAttrTextureAlignment = (UInt32)(14)
const cudaDevAttrGpuOverlap = (UInt32)(15)
const cudaDevAttrMultiProcessorCount = (UInt32)(16)
const cudaDevAttrKernelExecTimeout = (UInt32)(17)
const cudaDevAttrIntegrated = (UInt32)(18)
const cudaDevAttrCanMapHostMemory = (UInt32)(19)
const cudaDevAttrComputeMode = (UInt32)(20)
const cudaDevAttrMaxTexture1DWidth = (UInt32)(21)
const cudaDevAttrMaxTexture2DWidth = (UInt32)(22)
const cudaDevAttrMaxTexture2DHeight = (UInt32)(23)
const cudaDevAttrMaxTexture3DWidth = (UInt32)(24)
const cudaDevAttrMaxTexture3DHeight = (UInt32)(25)
const cudaDevAttrMaxTexture3DDepth = (UInt32)(26)
const cudaDevAttrMaxTexture2DLayeredWidth = (UInt32)(27)
const cudaDevAttrMaxTexture2DLayeredHeight = (UInt32)(28)
const cudaDevAttrMaxTexture2DLayeredLayers = (UInt32)(29)
const cudaDevAttrSurfaceAlignment = (UInt32)(30)
const cudaDevAttrConcurrentKernels = (UInt32)(31)
const cudaDevAttrEccEnabled = (UInt32)(32)
const cudaDevAttrPciBusId = (UInt32)(33)
const cudaDevAttrPciDeviceId = (UInt32)(34)
const cudaDevAttrTccDriver = (UInt32)(35)
const cudaDevAttrMemoryClockRate = (UInt32)(36)
const cudaDevAttrGlobalMemoryBusWidth = (UInt32)(37)
const cudaDevAttrL2CacheSize = (UInt32)(38)
const cudaDevAttrMaxThreadsPerMultiProcessor = (UInt32)(39)
const cudaDevAttrAsyncEngineCount = (UInt32)(40)
const cudaDevAttrUnifiedAddressing = (UInt32)(41)
const cudaDevAttrMaxTexture1DLayeredWidth = (UInt32)(42)
const cudaDevAttrMaxTexture1DLayeredLayers = (UInt32)(43)
const cudaDevAttrMaxTexture2DGatherWidth = (UInt32)(45)
const cudaDevAttrMaxTexture2DGatherHeight = (UInt32)(46)
const cudaDevAttrMaxTexture3DWidthAlt = (UInt32)(47)
const cudaDevAttrMaxTexture3DHeightAlt = (UInt32)(48)
const cudaDevAttrMaxTexture3DDepthAlt = (UInt32)(49)
const cudaDevAttrPciDomainId = (UInt32)(50)
const cudaDevAttrTexturePitchAlignment = (UInt32)(51)
const cudaDevAttrMaxTextureCubemapWidth = (UInt32)(52)
const cudaDevAttrMaxTextureCubemapLayeredWidth = (UInt32)(53)
const cudaDevAttrMaxTextureCubemapLayeredLayers = (UInt32)(54)
const cudaDevAttrMaxSurface1DWidth = (UInt32)(55)
const cudaDevAttrMaxSurface2DWidth = (UInt32)(56)
const cudaDevAttrMaxSurface2DHeight = (UInt32)(57)
const cudaDevAttrMaxSurface3DWidth = (UInt32)(58)
const cudaDevAttrMaxSurface3DHeight = (UInt32)(59)
const cudaDevAttrMaxSurface3DDepth = (UInt32)(60)
const cudaDevAttrMaxSurface1DLayeredWidth = (UInt32)(61)
const cudaDevAttrMaxSurface1DLayeredLayers = (UInt32)(62)
const cudaDevAttrMaxSurface2DLayeredWidth = (UInt32)(63)
const cudaDevAttrMaxSurface2DLayeredHeight = (UInt32)(64)
const cudaDevAttrMaxSurface2DLayeredLayers = (UInt32)(65)
const cudaDevAttrMaxSurfaceCubemapWidth = (UInt32)(66)
const cudaDevAttrMaxSurfaceCubemapLayeredWidth = (UInt32)(67)
const cudaDevAttrMaxSurfaceCubemapLayeredLayers = (UInt32)(68)
const cudaDevAttrMaxTexture1DLinearWidth = (UInt32)(69)
const cudaDevAttrMaxTexture2DLinearWidth = (UInt32)(70)
const cudaDevAttrMaxTexture2DLinearHeight = (UInt32)(71)
const cudaDevAttrMaxTexture2DLinearPitch = (UInt32)(72)
const cudaDevAttrMaxTexture2DMipmappedWidth = (UInt32)(73)
const cudaDevAttrMaxTexture2DMipmappedHeight = (UInt32)(74)
const cudaDevAttrComputeCapabilityMajor = (UInt32)(75)
const cudaDevAttrComputeCapabilityMinor = (UInt32)(76)
const cudaDevAttrMaxTexture1DMipmappedWidth = (UInt32)(77)
const cudaDevAttrStreamPrioritiesSupported = (UInt32)(78)
const cudaDevAttrGlobalL1CacheSupported = (UInt32)(79)
const cudaDevAttrLocalL1CacheSupported = (UInt32)(80)
const cudaDevAttrMaxSharedMemoryPerMultiprocessor = (UInt32)(81)
const cudaDevAttrMaxRegistersPerMultiprocessor = (UInt32)(82)
const cudaDevAttrManagedMemory = (UInt32)(83)
const cudaDevAttrIsMultiGpuBoard = (UInt32)(84)
const cudaDevAttrMultiGpuBoardGroupID = (UInt32)(85)
# end enum cudaDeviceAttr

immutable Array_256_Uint8
    d1::UInt8
    d2::UInt8
    d3::UInt8
    d4::UInt8
    d5::UInt8
    d6::UInt8
    d7::UInt8
    d8::UInt8
    d9::UInt8
    d10::UInt8
    d11::UInt8
    d12::UInt8
    d13::UInt8
    d14::UInt8
    d15::UInt8
    d16::UInt8
    d17::UInt8
    d18::UInt8
    d19::UInt8
    d20::UInt8
    d21::UInt8
    d22::UInt8
    d23::UInt8
    d24::UInt8
    d25::UInt8
    d26::UInt8
    d27::UInt8
    d28::UInt8
    d29::UInt8
    d30::UInt8
    d31::UInt8
    d32::UInt8
    d33::UInt8
    d34::UInt8
    d35::UInt8
    d36::UInt8
    d37::UInt8
    d38::UInt8
    d39::UInt8
    d40::UInt8
    d41::UInt8
    d42::UInt8
    d43::UInt8
    d44::UInt8
    d45::UInt8
    d46::UInt8
    d47::UInt8
    d48::UInt8
    d49::UInt8
    d50::UInt8
    d51::UInt8
    d52::UInt8
    d53::UInt8
    d54::UInt8
    d55::UInt8
    d56::UInt8
    d57::UInt8
    d58::UInt8
    d59::UInt8
    d60::UInt8
    d61::UInt8
    d62::UInt8
    d63::UInt8
    d64::UInt8
    d65::UInt8
    d66::UInt8
    d67::UInt8
    d68::UInt8
    d69::UInt8
    d70::UInt8
    d71::UInt8
    d72::UInt8
    d73::UInt8
    d74::UInt8
    d75::UInt8
    d76::UInt8
    d77::UInt8
    d78::UInt8
    d79::UInt8
    d80::UInt8
    d81::UInt8
    d82::UInt8
    d83::UInt8
    d84::UInt8
    d85::UInt8
    d86::UInt8
    d87::UInt8
    d88::UInt8
    d89::UInt8
    d90::UInt8
    d91::UInt8
    d92::UInt8
    d93::UInt8
    d94::UInt8
    d95::UInt8
    d96::UInt8
    d97::UInt8
    d98::UInt8
    d99::UInt8
    d100::UInt8
    d101::UInt8
    d102::UInt8
    d103::UInt8
    d104::UInt8
    d105::UInt8
    d106::UInt8
    d107::UInt8
    d108::UInt8
    d109::UInt8
    d110::UInt8
    d111::UInt8
    d112::UInt8
    d113::UInt8
    d114::UInt8
    d115::UInt8
    d116::UInt8
    d117::UInt8
    d118::UInt8
    d119::UInt8
    d120::UInt8
    d121::UInt8
    d122::UInt8
    d123::UInt8
    d124::UInt8
    d125::UInt8
    d126::UInt8
    d127::UInt8
    d128::UInt8
    d129::UInt8
    d130::UInt8
    d131::UInt8
    d132::UInt8
    d133::UInt8
    d134::UInt8
    d135::UInt8
    d136::UInt8
    d137::UInt8
    d138::UInt8
    d139::UInt8
    d140::UInt8
    d141::UInt8
    d142::UInt8
    d143::UInt8
    d144::UInt8
    d145::UInt8
    d146::UInt8
    d147::UInt8
    d148::UInt8
    d149::UInt8
    d150::UInt8
    d151::UInt8
    d152::UInt8
    d153::UInt8
    d154::UInt8
    d155::UInt8
    d156::UInt8
    d157::UInt8
    d158::UInt8
    d159::UInt8
    d160::UInt8
    d161::UInt8
    d162::UInt8
    d163::UInt8
    d164::UInt8
    d165::UInt8
    d166::UInt8
    d167::UInt8
    d168::UInt8
    d169::UInt8
    d170::UInt8
    d171::UInt8
    d172::UInt8
    d173::UInt8
    d174::UInt8
    d175::UInt8
    d176::UInt8
    d177::UInt8
    d178::UInt8
    d179::UInt8
    d180::UInt8
    d181::UInt8
    d182::UInt8
    d183::UInt8
    d184::UInt8
    d185::UInt8
    d186::UInt8
    d187::UInt8
    d188::UInt8
    d189::UInt8
    d190::UInt8
    d191::UInt8
    d192::UInt8
    d193::UInt8
    d194::UInt8
    d195::UInt8
    d196::UInt8
    d197::UInt8
    d198::UInt8
    d199::UInt8
    d200::UInt8
    d201::UInt8
    d202::UInt8
    d203::UInt8
    d204::UInt8
    d205::UInt8
    d206::UInt8
    d207::UInt8
    d208::UInt8
    d209::UInt8
    d210::UInt8
    d211::UInt8
    d212::UInt8
    d213::UInt8
    d214::UInt8
    d215::UInt8
    d216::UInt8
    d217::UInt8
    d218::UInt8
    d219::UInt8
    d220::UInt8
    d221::UInt8
    d222::UInt8
    d223::UInt8
    d224::UInt8
    d225::UInt8
    d226::UInt8
    d227::UInt8
    d228::UInt8
    d229::UInt8
    d230::UInt8
    d231::UInt8
    d232::UInt8
    d233::UInt8
    d234::UInt8
    d235::UInt8
    d236::UInt8
    d237::UInt8
    d238::UInt8
    d239::UInt8
    d240::UInt8
    d241::UInt8
    d242::UInt8
    d243::UInt8
    d244::UInt8
    d245::UInt8
    d246::UInt8
    d247::UInt8
    d248::UInt8
    d249::UInt8
    d250::UInt8
    d251::UInt8
    d252::UInt8
    d253::UInt8
    d254::UInt8
    d255::UInt8
    d256::UInt8
end

zero(::Type{Array_256_Uint8}) = Array_256_Uint8(fill(zero(UInt8),256)...)

immutable Array_3_Cint
    d1::Cint
    d2::Cint
    d3::Cint
end

zero(::Type{Array_3_Cint}) = Array_3_Cint(fill(zero(Cint),3)...)

immutable Array_2_Cint
    d1::Cint
    d2::Cint
end

zero(::Type{Array_2_Cint}) = Array_2_Cint(fill(zero(Cint),2)...)

immutable cudaDeviceProp
    name::Array_256_Uint8
    totalGlobalMem::Csize_t
    sharedMemPerBlock::Csize_t
    regsPerBlock::Cint
    warpSize::Cint
    memPitch::Csize_t
    maxThreadsPerBlock::Cint
    maxThreadsDim::Array_3_Cint
    maxGridSize::Array_3_Cint
    clockRate::Cint
    totalConstMem::Csize_t
    major::Cint
    minor::Cint
    textureAlignment::Csize_t
    texturePitchAlignment::Csize_t
    deviceOverlap::Cint
    multiProcessorCount::Cint
    kernelExecTimeoutEnabled::Cint
    integrated::Cint
    canMapHostMemory::Cint
    computeMode::Cint
    maxTexture1D::Cint
    maxTexture1DMipmap::Cint
    maxTexture1DLinear::Cint
    maxTexture2D::Array_2_Cint
    maxTexture2DMipmap::Array_2_Cint
    maxTexture2DLinear::Array_3_Cint
    maxTexture2DGather::Array_2_Cint
    maxTexture3D::Array_3_Cint
    maxTexture3DAlt::Array_3_Cint
    maxTextureCubemap::Cint
    maxTexture1DLayered::Array_2_Cint
    maxTexture2DLayered::Array_3_Cint
    maxTextureCubemapLayered::Array_2_Cint
    maxSurface1D::Cint
    maxSurface2D::Array_2_Cint
    maxSurface3D::Array_3_Cint
    maxSurface1DLayered::Array_2_Cint
    maxSurface2DLayered::Array_3_Cint
    maxSurfaceCubemap::Cint
    maxSurfaceCubemapLayered::Array_2_Cint
    surfaceAlignment::Csize_t
    concurrentKernels::Cint
    ECCEnabled::Cint
    pciBusID::Cint
    pciDeviceID::Cint
    pciDomainID::Cint
    tccDriver::Cint
    asyncEngineCount::Cint
    unifiedAddressing::Cint
    memoryClockRate::Cint
    memoryBusWidth::Cint
    l2CacheSize::Cint
    maxThreadsPerMultiProcessor::Cint
    streamPrioritiesSupported::Cint
    globalL1CacheSupported::Cint
    localL1CacheSupported::Cint
    sharedMemPerMultiprocessor::Csize_t
    regsPerMultiprocessor::Cint
    managedMemory::Cint
    isMultiGpuBoard::Cint
    multiGpuBoardGroupID::Cint
end

immutable Array_64_Uint8
    d1::UInt8
    d2::UInt8
    d3::UInt8
    d4::UInt8
    d5::UInt8
    d6::UInt8
    d7::UInt8
    d8::UInt8
    d9::UInt8
    d10::UInt8
    d11::UInt8
    d12::UInt8
    d13::UInt8
    d14::UInt8
    d15::UInt8
    d16::UInt8
    d17::UInt8
    d18::UInt8
    d19::UInt8
    d20::UInt8
    d21::UInt8
    d22::UInt8
    d23::UInt8
    d24::UInt8
    d25::UInt8
    d26::UInt8
    d27::UInt8
    d28::UInt8
    d29::UInt8
    d30::UInt8
    d31::UInt8
    d32::UInt8
    d33::UInt8
    d34::UInt8
    d35::UInt8
    d36::UInt8
    d37::UInt8
    d38::UInt8
    d39::UInt8
    d40::UInt8
    d41::UInt8
    d42::UInt8
    d43::UInt8
    d44::UInt8
    d45::UInt8
    d46::UInt8
    d47::UInt8
    d48::UInt8
    d49::UInt8
    d50::UInt8
    d51::UInt8
    d52::UInt8
    d53::UInt8
    d54::UInt8
    d55::UInt8
    d56::UInt8
    d57::UInt8
    d58::UInt8
    d59::UInt8
    d60::UInt8
    d61::UInt8
    d62::UInt8
    d63::UInt8
    d64::UInt8
end

zero(::Type{Array_64_Uint8}) = Array_64_Uint8(fill(zero(UInt8),64)...)

immutable cudaIpcEventHandle_st
    reserved::Array_64_Uint8
end

immutable cudaIpcEventHandle_t
    reserved::Array_64_Uint8
end

immutable cudaIpcMemHandle_st
    reserved::Array_64_Uint8
end

immutable cudaIpcMemHandle_t
    reserved::Array_64_Uint8
end

const cudaError_t = cudaError
const CUstream_st = Void
const cudaStream_t = Ptr{CUstream_st}
const CUevent_st = Void
const cudaEvent_t = Ptr{CUevent_st}
const cudaGraphicsResource_t = Ptr{cudaGraphicsResource}
const CUuuid_st = Void
const cudaUUID_t = Void
const cudaOutputMode_t = cudaOutputMode

immutable char1
    x::UInt8
end

immutable uchar1
    x::Cuchar
end

immutable char3
    x::UInt8
    y::UInt8
    z::UInt8
end

immutable uchar3
    x::Cuchar
    y::Cuchar
    z::Cuchar
end

immutable short1
    x::Int16
end

immutable ushort1
    x::UInt16
end

immutable short3
    x::Int16
    y::Int16
    z::Int16
end

immutable ushort3
    x::UInt16
    y::UInt16
    z::UInt16
end

immutable int1
    x::Cint
end

immutable uint1
    x::UInt32
end

immutable int3
    x::Cint
    y::Cint
    z::Cint
end

immutable uint3
    x::UInt32
    y::UInt32
    z::UInt32
end

immutable long1
    x::Clong
end

immutable ulong1
    x::Culong
end

immutable long3
    x::Clong
    y::Clong
    z::Clong
end

immutable ulong3
    x::Culong
    y::Culong
    z::Culong
end

immutable float1
    x::Cfloat
end

immutable float3
    x::Cfloat
    y::Cfloat
    z::Cfloat
end

immutable longlong1
    x::Clonglong
end

immutable ulonglong1
    x::Culonglong
end

immutable longlong3
    x::Clonglong
    y::Clonglong
    z::Clonglong
end

immutable ulonglong3
    x::Culonglong
    y::Culonglong
    z::Culonglong
end

immutable double1
    x::Cdouble
end

immutable double3
    x::Cdouble
    y::Cdouble
    z::Cdouble
end

immutable dim3
    x::UInt32
    y::UInt32
    z::UInt32
end

const CUDART_VERSION = 6050

()

const cudaStreamCallback_t = Ptr{Void}


text:
function apply_boundary_conditions_on_z_axis!(  rbpot::Array{T, 4}, ix::Int, iy::Int, rbi::Int, ax::DiscreteAxis{T, :infinite, :fixed}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[ix, iy, 1, rbi] = grid_boundary_factors[1] * rbpot[ix, iy, 2, rbi] # infinity boundaries in r
    nothing
end
function apply_boundary_conditions_on_z_axis!(  rbpot::Array{T, 4}, ix::Int, iy::Int, rbi::Int, ax::DiscreteAxis{T, :infinite, :infinite}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[ix, iy, 1, rbi] = grid_boundary_factors[1] * rbpot[ix, iy, 2, rbi]
    rbpot[ix, iy, end, rbi] = grid_boundary_factors[2] * rbpot[ix, iy, end - 1, rbi]
    nothing
end
function apply_boundary_conditions_on_z_axis!(  rbpot::Array{T, 4}, ix::Int, iy::Int, rbi::Int, ax::DiscreteAxis{T, :reflecting, :fixed}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[ix, iy, 1, rbi] = rbpot[ix, iy, 3, rbi] # infinity boundaries in r
    nothing
end
function apply_boundary_conditions_on_z_axis!(  rbpot::Array{T, 4}, ix::Int, iy::Int, rbi::Int, ax::DiscreteAxis{T, :fixed, :infinite}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[ix, iy, end, rbi] = grid_boundary_factors[2] * rbpot[ix, iy, end - 1, rbi] # infinity boundaries in r
    nothing
end
function apply_boundary_conditions_on_z_axis!(  rbpot::Array{T, 4}, ix::Int, iy::Int, rbi::Int, ax::DiscreteAxis{T, :fixed, :reflecting}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[ix, iy, end, rbi] = rbpot[ix, iy, end - 2, rbi] # infinity boundaries in r
    nothing
end
function apply_boundary_conditions_on_z_axis!(  rbpot::Array{T, 4}, ix::Int, iy::Int, rbi::Int, ax::DiscreteAxis{T, :reflecting, :reflecting}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[ix, iy, 1, rbi] = rbpot[ix, iy, 3, rbi]
    rbpot[ix, iy, end, rbi] = rbpot[ix, iy, end - 2, rbi]
    nothing
end
function apply_boundary_conditions_on_z_axis!(  rbpot::Array{T, 4}, ix::Int, iy::Int, rbi::Int, ax::DiscreteAxis{T, :fixed, :fixed}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    nothing
end
function apply_boundary_conditions_on_z_axis!(  rbpot::Array{T, 4}, ix::Int, iy::Int, rbi::Int, ax::DiscreteAxis{T, :periodic, :periodic}, int::Interval{:closed, :open, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[ix, iy,   1, rbi] = rbpot[ ix, iy, end - 1, rbi]
    rbpot[ix, iy, end, rbi] = rbpot[ ix, iy,       2, rbi]
    nothing
end

function apply_boundary_conditions_on_x_axis!(  rbpot::Array{T, 4}, iy::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :fixed, :fixed}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    nothing
end
function apply_boundary_conditions_on_x_axis!(  rbpot::Array{T, 4}, iy::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :infinite, :infinite}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[  1, iy, iz, rbi] = grid_boundary_factors[1] * rbpot[ 2, iy, iz, rbi]  # hmm, this is probably not fully correct since this is the red-black dimension
    rbpot[end, iy, iz, rbi] = grid_boundary_factors[2] * rbpot[ end - 1, iy, iz, rbi]  # anyhow its just an approximation
    nothing
end
function apply_boundary_conditions_on_x_axis!(  rbpot::Array{T, 4}, iy::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :fixed, :infinite}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[end, iy, iz, rbi] = grid_boundary_factors[2] * rbpot[ end - 1, iy, iz, rbi]    # infinity boundaries in x (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_x_axis!(  rbpot::Array{T, 4}, iy::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :infinite, :fixed}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[  1, iy, iz, rbi] = grid_boundary_factors[1] * rbpot[ 2, iy, iz, rbi]  # infinity boundaries in x (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_x_axis!(  rbpot::Array{T, 4}, iy::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :infinite, :reflecting}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[  1, iy, iz, rbi] = grid_boundary_factors[1] * rbpot[ 2, iy, iz, rbi]  # infinity boundaries in x (only ±1 because this is the compressed (redblack) dimension)
    rbpot[end, iy, iz, rbi] = rbpot[end - 1, iy, iz, rbi] # (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_x_axis!(  rbpot::Array{T, 4}, iy::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :reflecting, :infinite}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[1, iy, iz, rbi]   = rbpot[2, iy, iz, rbi] #  (only ±1 because this is the compressed (redblack) dimension)
    rbpot[end, iy, iz, rbi] = grid_boundary_factors[2] * rbpot[ end - 1, iy, iz, rbi]    # infinity boundaries in x (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_x_axis!(  rbpot::Array{T, 4}, iy::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :reflecting, :reflecting}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[1, iy, iz, rbi]   = rbpot[2, iy, iz, rbi] #  (only ±1 because this is the compressed (redblack) dimension)
    rbpot[end, iy, iz, rbi] = rbpot[end - 1, iy, iz, rbi] # (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_x_axis!(  rbpot::Array{T, 4}, iy::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :fixed, :reflecting}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[end, iy, iz, rbi] = rbpot[end - 1, iy, iz, rbi] # (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_x_axis!(  rbpot::Array{T, 4}, iy::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :reflecting, :fixed}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[1, iy, iz, rbi]   = rbpot[2, iy, iz, rbi] #  (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_x_axis!(  rbpot::Array{T, 4}, iy::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :periodic, :periodic}, int::Interval{:closed, :open, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[  1, iy, iz, rbi] = rbpot[end, iy, iz, rbi]
    rbpot[end, iy, iz, rbi] = rbpot[  1, iy, iz, rbi]
    nothing
end

function apply_boundary_conditions_on_y_axis!(  rbpot::Array{T, 4}, ix::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :fixed, :fixed}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T}, only2d::Val{only_2d})::Nothing where {T, only_2d}
    nothing
end
function apply_boundary_conditions_on_y_axis!(  rbpot::Array{T, 4}, ix::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :infinite, :infinite}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T}, only2d::Val{only_2d})::Nothing where {T, only_2d}
    rbpot[ ix, 1, iz, rbi]   = grid_boundary_factors[1] * rbpot[ ix, 2, iz, rbi]
    rbpot[ ix, end, iz, rbi] = grid_boundary_factors[2] * rbpot[ ix, end - 1, iz, rbi]
    nothing
end
function apply_boundary_conditions_on_y_axis!(  rbpot::Array{T, 4}, ix::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :fixed, :infinite}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T}, only2d::Val{only_2d})::Nothing where {T, only_2d}
    rbpot[ ix, end, iz, rbi] = grid_boundary_factors[2] * rbpot[ ix, end - 1, iz, rbi]
    nothing
end
function apply_boundary_conditions_on_y_axis!(  rbpot::Array{T, 4}, ix::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :infinite, :fixed}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T}, only2d::Val{only_2d})::Nothing where {T, only_2d}
    rbpot[ ix, 1, iz, rbi]   = grid_boundary_factors[1] * rbpot[ ix, 2, iz, rbi]
    nothing
end
function apply_boundary_conditions_on_y_axis!(  rbpot::Array{T, 4}, ix::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :reflecting, :reflecting}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T}, only2d::Val{only_2d})::Nothing where {T, only_2d}
    rbpot[ ix, 1, iz, rbi]   = rbpot[ ix, 3, iz, rbi]
    rbpot[ ix, end, iz, rbi] = rbpot[ ix, end - 2, iz, rbi]
    nothing
end
function apply_boundary_conditions_on_y_axis!(  rbpot::Array{T, 4}, ix::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :fixed, :reflecting}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T}, only2d::Val{only_2d})::Nothing where {T, only_2d}
    rbpot[ ix, end, iz, rbi] = rbpot[ ix, end - 2, iz, rbi]
    nothing
end
function apply_boundary_conditions_on_y_axis!(  rbpot::Array{T, 4}, ix::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :reflecting, :fixed}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T}, only2d::Val{only_2d})::Nothing where {T, only_2d}
    rbpot[ ix, 1, iz, rbi]   = rbpot[ ix, 3, iz, rbi]
    nothing
end
function apply_boundary_conditions_on_y_axis!(  rbpot::Array{T, 4}, ix::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :infinite, :reflecting}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T}, only2d::Val{only_2d})::Nothing where {T, only_2d}
    rbpot[ ix, 1, iz, rbi]   = grid_boundary_factors[1] * rbpot[ ix, 2, iz, rbi]
    rbpot[ ix, end, iz, rbi] = rbpot[ ix, end - 2, iz, rbi]
    nothing
end
function apply_boundary_conditions_on_y_axis!(  rbpot::Array{T, 4}, ix::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :reflecting, :infinite}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T}, only2d::Val{only_2d})::Nothing where {T, only_2d}
    rbpot[ ix, 1, iz, rbi]   = rbpot[ ix, 3, iz, rbi]
    rbpot[ ix, end, iz, rbi] = grid_boundary_factors[2] * rbpot[ ix, end - 1, iz, rbi]
    nothing
end
function apply_boundary_conditions_on_y_axis!(  rbpot::Array{T, 4}, ix::Int, iz::Int, rbi::Int, ax::DiscreteAxis{T, :periodic, :periodic}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T}, only2d::Val{only_2d})::Nothing where {T, only_2d}
    rbpot[ ix, 1, iz, rbi]   = rbpot[ ix, end - 1, iz, rbi]
    rbpot[ ix, end, iz, rbi] = rbpot[ ix, 2, iz, rbi]
    nothing
end

function apply_boundary_conditions!(fssrb::PotentialSimulationSetupRB{T, N1, N2, :cartesian}, update_even_points::Val{even_points}, only2d::Val{only_2d}) where {T, N1, N2, even_points, only_2d}
    rbi::Int = even_points ? rb_even::Int : rb_odd::Int

    if only_2d
        iy::Int = 2
        error("Boundary handling for 2D simulation (x&z) in cartesian coordinated not yet implemented.")
    else
        @inbounds for ix in axes(fssrb.potential, 1)
            for iy in axes(fssrb.potential, 2)
                apply_boundary_conditions_on_z_axis!( fssrb.potential, ix, iy, rbi, fssrb.grid.axes[3], fssrb.grid.axes[3].interval, fssrb.grid_boundary_factors[3])
            end
            for iz in axes(fssrb.potential, 3)
                apply_boundary_conditions_on_y_axis!( fssrb.potential, ix, iz, rbi, fssrb.grid.axes[2], fssrb.grid.axes[2].interval, fssrb.grid_boundary_factors[2], only2d)
            end
        end
        @inbounds for iy in axes(fssrb.potential, 2) # z boundaries
            for iz in axes(fssrb.potential, 3)
                apply_boundary_conditions_on_x_axis!( fssrb.potential, iy, iz, rbi, fssrb.grid.axes[1], fssrb.grid.axes[1].interval, fssrb.grid_boundary_factors[1])
            end
        end
    end
    nothing
end

ext: jl


text:
#derrida_map calculates derrida map for a net up to a  H hamming distance
function derrida_map!(net::Net2{T},nrand_init::Int=prod([length(sigma.state_range) for sigma in net.nodes]);tau::Int=1,update_mode::Symbol=:sy,init_h::Symbol=:rand,Hmax::Int=sum(length(sigma.state_range)>1 for sigma in net.nodes)*(1 + (update_mode==:sy ? 0 : maximum(net.in_taus)))) where {T<:state_type}
    index_nodes_to_flip = find(sigma->length(sigma.state_range)>1,net.nodes)
    net_hist_size = 1 + (update_mode==:sy ? 0 : maximum(net.in_taus))
    Nmax = length(index_nodes_to_flip)*net_hist_size
    net_size = length(net.nodes)
    1<=Hmax<=Nmax || error("Hmax must be less or equal than network size")
    state_space_size = prod([length(sigma.state_range) for sigma in net.nodes])
    if nrand_init>state_space_size
        println("Oversampling error: Random initial conditions bigger than possible number of configurations")
        nrand_init = state_space_size
    end
    mut_indices = collect(product(index_nodes_to_flip,1:net_hist_size))
    hXa = Array{Int}(net_size,net_hist_size)
    hXatp = Array{Int}(net_size,net_hist_size)
    hXb = Array{Int}(net_size,net_hist_size)
    hXbtp = Array{Int}(net_size,net_hist_size)
    # sim_buffer = Array{T}(net_size,tau + net_hist_size)
    Hmap = Tuple{Float64,Float64}[]
    i = 0
    while i < nrand_init
        rand_init_net!(net;update_mode=update_mode,init_h=init_h)
        t = rand(1:net_hist_size)-1
        evolve_net!(net,net_hist_size-1,keep=false,update_mode=update_mode,t_init=t)
        copy!(hXa,net.sim_buffer[:,end-net_hist_size+1:end])
        evolve_net!(net,tau,keep=false,update_mode=update_mode,t_init=t)
        copy!(hXatp,net.sim_buffer[:,end-net_hist_size+1:end])
        for H in 1:Hmax
            copy!(hXb,hXa)
            for (n,l) in sample(mut_indices, H, replace = false)
                hXb[n,l] = rand(filter(x -> x!=hXa[n,l],net.nodes[n].state_range))
            end
            # print("H=$H :",sum(hXa.!=hXb)," -> ")
            # init_net!(net;update_mode=update_mode,init_array=vec(hXb))
            for ni in 1:net_size
                net.sim_buffer[ni,end-net_hist_size+1:end] = collect(hXb[ni,:])
            end
            evolve_net!(net,tau,keep=false,update_mode=update_mode,t_init=t)
            copy!(hXbtp,net.sim_buffer[:,end-net_hist_size+1:end])
            # println(sum(net.sim_buffer[:,end-net_hist_size+1:end].!=hXatp))
            push!(Hmap,(H/Nmax,sum(hXatp.!=hXbtp)/Nmax))
            if mod(i,10000)==1
                @printf("%0.2f%% progress, %d out of %d random initial conditions explored at distance H(0)=%d\n", 100*i/nrand_init, i, nrand_init,H)
                println("H(0)=$H:")
                # println(join(map(string,hXa),","))
                # println(join(map(string,hXb),","))
                println("H($tau)=",sum(hXatp.!=hXbtp),":")
                # println(join(map(string,hXatp),","))
                # println(join(map(string,hXbtp),","))
            end
        end
        i += 1
    end
    return reshape(Hmap,Hmax,nrand_init)
end

function derrida_map_mod!(netWT::Net2{T},nrand_init::Int=prod([length(sigma.state_range) for sigma in netWT.nodes]);mutants::Vector{Pair{String,Symbol}} = Vector{Pair{String,Symbol}}(),tau::Int=1,fix_targets::Bool=true,update_mode::Symbol=:sy,init_h::Symbol=:rand,Hmax::Int=sum([length(sigma.state_range)>1 for sigma in netWT.nodes])) where {T<:state_type}
    netMT = deepcopy(netWT)
    for (id,pert) in mutants
        if pert == :KO
            knock_out!(netMT,id)
        elseif pert == :KD
            knock_down!(netMT,id)
        else
            error("Invalid method")
        end

    end
    Mutant_state_space_size = prod([length(sigma.state_range) for sigma in netMT.nodes])
    Mutant_state_space_size = length(combinations(1:Mutant_state_space_size,2))
    if nrand_init>Mutant_state_space_size
        nrand_init = Mutant_state_space_size
        println("Oversampling error: Random initial conditions bigger than possible number of configurations")
    end
    net_hist_size = 1 + (update_mode==:sy ? 0 : maximum(netWT.in_taus))
    index_nodes_to_fix = Integer[]
    if fix_targets
        for (id,_) in mutants
            i = findfirst(x->x.id==id,netWT.nodes)
            j = findfirst(x->x.id==id,netMT.nodes)
            (i != 0 && j != 0) || error("Trying to set unexisting $id node")
            #i == j || error("Comparative derrida map requires networks to have the same size and same node order")
            push!(index_nodes_to_fix,i)
        end
    end
    Hmin = length(index_nodes_to_fix)
    index_nodes_to_flip = find(x->length(x.state_range)>1, netWT.nodes)
    deleteat!(index_nodes_to_flip,findin(index_nodes_to_flip,index_nodes_to_fix))
    NWT = length(netWT.nodes)
    NMT = length(netMT.nodes)
    net_size = NWT
    #Nmax = length(index_nodes_to_flip)*net_hist_size
    Nmax = NWT*net_hist_size
    Hmin<=Hmax<=Nmax || error("Hmax must be greater or equal than size of fixed targets, and equal or less than network size")
    pairs = Dict(i=>Dict{AbstractString,AbstractString}() for i in Hmin:Hmax)
    mut_indices = shuffle(collect(product(index_nodes_to_flip,1:net_hist_size)))
    hXWT = Array{Int}(NWT,net_hist_size)
    hXWTtptau = collect(Array{Int}(NWT,net_hist_size) for _ in 1:tau)
    hXMT = Array{Int}(NMT,net_hist_size)
    hXMTtptau = collect(Array{Int}(NMT,net_hist_size) for _ in 1:tau)
    Hmap = Tuple{Int64,Vararg{Float64,2}}[]
    i = 0
    while i < nrand_init
        rand_init_net!(netMT;update_mode=update_mode,init_h=init_h)
        copy!(hXMT,netMT.sim_buffer[:,end-net_hist_size+1:end])
        t = rand(1:net_hist_size)-1
        netMTtimeseries = evolve_net!(netMT,tau,keep=true,update_mode=update_mode,t_init=t)
        for ti in 1:tau
            copy!(hXMTtptau[ti],netMTtimeseries[:,1+tau-net_hist_size+1:1+tau])
        end
        for H in Hmin:Hmax
            copy!(hXWT,hXMT)
            for (n,l) in sample(mut_indices, H-Hmin, replace = false)
                hXWT[n,l] = rand(filter(x -> x!=hXMT[n,l],netWT.nodes[n].state_range))
            end
            for n in index_nodes_to_fix
                hXWT[n,end] = rand(filter(x -> x!=hXMT[n,end],netWT.nodes[n].state_range))
            end # those nodes to be shut off in the knock out network are set in the wild type network to any other state different than that
            for ni in 1:NWT
                netWT.sim_buffer[ni,end-net_hist_size+1:end] = collect(hXWT[ni,:])
            end
            init_net!(netWT;t0_mode=:custom,update_mode=update_mode,init_h=:last)
            netWTtimeseries = evolve_net!(netWT,tau,keep=true,update_mode=update_mode,t_init=t)
            for ti in 1:tau
                copy!(hXWTtptau[ti],netWTtimeseries[:,1+tau-net_hist_size+1:1+tau])
                push!(Hmap,(ti,sum(hXWT.!=hXMT)/Nmax,sum(hXWTtptau[ti].!=hXMTtptau[ti])/Nmax))
            end
        end
        i += 1
        if mod(i,10000)==0
            @printf("%0.2f%% progress, %d out of %d random initial conditions explored\n", 100*i/nrand_init, i, nrand_init)
        end
    end
    return reshape(Hmap,tau*(Hmax-Hmin+1),nrand_init)
end

# get_basins explores all possible configurations and returns a transition graph, with a dictionary storing labels for
# each configuration

function get_basins!(net::Net2,out_dir::AbstractString,root_name::AbstractString;trav::Symbol=:rw,update_mode::Symbol=:sy,max_traj::Integer=1 + (update_mode==:sy ? 0 : maximum(net.in_taus)))
    if isfile(joinpath(out_dir,"$(root_name)_att_stats.txt"))
        println("Previous analysis already existing. Please delete files to proceed.:",joinpath(out_dir,"$(root_name)_att_stats.txt"))
        return nothing
    end
    trav == :rw || trav == :l || error("Invalid option for traversal method")
    net_size = nv(net.graph)
    state_space_size = prod(length(sigma.state_range) for sigma in net.nodes)
    if state_space_size==0 && all(x->isinteger(x)&&x>0,length(sigma.state_range) for sigma in net.nodes)
        error("Full exploration of state space exceeds memory limit")
    end
    itstspace = product((sigma.state_range for sigma in net.nodes)...)
    state_space_size = length(itstspace)
    transition_graph = DiGraph()
    tags_nid = update_mode == :sy ? Dict{Int,Int}() : Dict{Vector{Int},Int}()
    nid_tags = update_mode == :sy ? Vector{Int}() : Vector{Vector{Int}}()
    tag_base = maximum(length(sigma.state_range) for sigma in net.nodes)
    println("Calculating attraction landscape for $root_name")
    tau_max = maximum(net.in_taus)
    net_hist_size = tau_max + 1
    if update_mode == :sy
        i = 0
        for leave in itstspace
            init_net!(net;update_mode=update_mode,init_array=round.(Int,collect(leave)))
            leave_tag = parse(Int,string(map(x->round.(Int,x),leave)...),tag_base)
            if !haskey(tags_nid,leave_tag)
                add_vertex!(transition_graph)
                i = nv(transition_graph)
                push!(nid_tags,leave_tag)
                tags_nid[leave_tag] = i
            end
            shoot = evolve_net!(net,1;keep=false,update_mode=update_mode)
            shoot_tag = parse(Int,string(round.(Int,shoot)...),tag_base)
            if !haskey(tags_nid,shoot_tag)
                add_vertex!(transition_graph)
                i = nv(transition_graph)
                push!(nid_tags,shoot_tag)
                tags_nid[shoot_tag] = i
            end
            add_edge!(transition_graph,tags_nid[leave_tag],tags_nid[shoot_tag])
            if mod(i,100000)==0
                @printf "%0.3f%% configurations explored\n" 100*i/state_space_size
            end
        end
    elseif update_mode == :fd || update_mode == :md
        hX = Array{Int8}(net_size,net_hist_size)
        hX_tag = Vector{Int}(net_hist_size)
        hXtp1 = Array{Int8}(net_size,net_hist_size)
        hXtp1_tag = Vector{Int}(net_hist_size)
        time_steps = 1
        sim_buffer = Array{Int8}(net_size,time_steps + net_hist_size)
        for (i,X) in enumerate(itstspace)
            init_net!(net;init_h=:default,update_mode=update_mode,init_array=round.(Int,collect(X)))
            hX[:] = net.sim_buffer
            # for (nsigma,sigma) in enumerate(netstate)
            #     init_node!(sigma,X[nsigma])
            #     hX[nsigma,:] = sigma.state
            # end
            for t in 1:net_hist_size
                hX_tag[t] = parse(Int,string(round.(Int,hX[:,t])...),tag_base)
            end
            if !haskey(tags_nid,hX_tag)
                add_vertex!(transition_graph)
                nvertex = nv(transition_graph)
                push!(nid_tags, collect(hX_tag))
                tags_nid[collect(hX_tag)] = nvertex
                # println(hX_tag, " tag added, vertex ",tags_nid[hX_tag])
            end
            # hXtp1_tag[:] = hX_tag[:]
            # shift!(hXtp1_tag)
            for t in 1:max_traj
                hXtp1[:] = evolve_net!(net,1,keep=false,update_mode=update_mode,sim=sim_buffer,t_init=t-1)
                for ti in 1:net_hist_size
                    hXtp1_tag[ti] = parse(Int,string(round.(Int,hXtp1[:,ti])...),tag_base)
                end
                # Xtp1_tag = parse(Int,string(round(Int,Xtp1)...),tag_base)
                # push!(hXtp1_tag,Xtp1_tag)
                if !haskey(tags_nid,hXtp1_tag)
                    add_vertex!(transition_graph)
                    nvertex = nv(transition_graph)
                    push!(nid_tags,collect(hXtp1_tag))
                    tags_nid[collect(hXtp1_tag)] = nvertex
                    # println(hXtp1_tag, " tag added, vertex ",tags_nid[hXtp1_tag])
                end
                # foreach(x->println(x," ",typeof(x)," ",hX_tag==x),keys(tags_nid))
                add_edge!(transition_graph,tags_nid[hX_tag],tags_nid[hXtp1_tag]) || break
                copy!(hX_tag,hXtp1_tag)
                # shift!(hXtp1_tag)
            end
            if mod(i,100000)==0
                @printf "%0.2f%% initial conditions explored\n" 100*i/state_space_size
            end
        end
    elseif update_mode == :st || update_mode == :mst
        hX = Array{Int8}(net_size,net_hist_size)
        hX_tag = Vector{Int}(net_hist_size+1)
        hXtp1 = Array{Int8}(net_size,net_hist_size)
        hXtp1_tag = Vector{Int}(net_hist_size+1)
        time_steps = 1
        sim_buffer = Array{Int8}(net_size,time_steps + net_hist_size)
        for (i,X) in enumerate(itstspace)
            init_net!(net;init_h=:default,update_mode=update_mode,init_array=round.(Int,collect(X)))
            hX[:] = net.sim_buffer
            # for (nsigma,sigma) in enumerate(netstate)
            #     init_node!(sigma,X[nsigma])
            #     hX[nsigma,:] = sigma.state
            # end
            for t in 1:net_hist_size
                hX_tag[1+t] = parse(Int,string(round.(Int,hX[:,t])...),tag_base)
            end
            hX_tag[1] = 0
            if !haskey(tags_nid,hX_tag)
                add_vertex!(transition_graph)
                nvertex = nv(transition_graph)
                push!(nid_tags, collect(hX_tag))
                tags_nid[collect(hX_tag)] = nvertex
                # println(hX_tag, " tag added, vertex ",tags_nid[hX_tag])
            end
            # hXtp1_tag[:] = hX_tag[:]
            # shift!(hXtp1_tag)
            # tmax = 0
            for t in 1:max_traj
                hXtp1[:] = evolve_net!(net,1,keep=false,update_mode=update_mode,sim=sim_buffer,t_init=t-1)
                for ti in 1:net_hist_size
                    hXtp1_tag[1+ti] = parse(Int,string(round.(Int,hXtp1[:,ti])...),tag_base)
                end
                hXtp1_tag[1] = mod(t,net_hist_size)
                # Xtp1_tag = parse(Int,string(round(Int,Xtp1)...),tag_base)
                # push!(hXtp1_tag,Xtp1_tag)
                # println(t, " ",keys(tags_nid))
                if !haskey(tags_nid,hXtp1_tag)
                    add_vertex!(transition_graph)
                    nvertex = nv(transition_graph)
                    push!(nid_tags, collect(hXtp1_tag))
                    tags_nid[collect(hXtp1_tag)] = nvertex
                    # println(hXtp1_tag, " tag added, vertex ",tags_nid[hXtp1_tag])
                end
                # println(t, " ",keys(tags_nid))
                # tmax = t
                add_edge!(transition_graph,tags_nid[hX_tag],tags_nid[hXtp1_tag]) || break
                copy!(hX_tag,hXtp1_tag)
                # shift!(hXtp1_tag)
            end
            # println(tmax)
            if mod(i,100000)==0
                @printf "%0.2f%% initial conditions explored\n" 100*i/state_space_size
            end
        end
    end
    state_space_size = nv(transition_graph)
    println("Decomposing transition graph into attraction basins")
    basin_node_list = weakly_connected_components(transition_graph)
    Natt = length(basin_node_list)
    println(Natt," attractors found")
    sort!(basin_node_list,by=length,rev=true)
    println("Saving individual basins")
    if !isdir(out_dir)
        mkpath(out_dir)
    end
    basin_tree_list = Vector{LightGraphs.DiGraph}(Natt)
    basin_tags = Vector{Vector{Integer}}(Natt)
    att_list = Vector{Vector{Integer}}(Natt)
    periods = Vector{Integer}(Natt)
    branch_sizes = Vector{Vector{Integer}}(Natt)
    patt_basins = Vector{Vector{Integer}}(Natt)
    branch_tree_list = Vector{Vector{LightGraphs.DiGraph}}(Natt)
    branch_patt = Vector{Vector{Integer}}(Natt)
    g = GZip.open(joinpath(out_dir,"$(root_name)_att_tagging.dat.gz"),"w")
    h = GZip.open(joinpath(out_dir,"$(root_name)_trans_dist.dat.gz"),"w")
    try

        for (bi,basin) in enumerate(basin_node_list)
            basin_tree_list[bi] = induced_subgraph(transition_graph,basin)[1]
            println("Basin $bi has ",nv(basin_tree_list[bi])," space state size")
            savegraph(joinpath(out_dir,"$(root_name)_att_$(bi).gml.gz"),basin_tree_list[bi],join([root_name, bi],"_"),GMLFormat(),compress = true)
            println("$(root_name)_att_$(bi) basin saved")
            if update_mode == :sy
                writedlm(g,hcat(fill(bi,nv(basin_tree_list[bi])),nid_tags[basin]))
            else
                writedlm(g,hcat(fill(bi,nv(basin_tree_list[bi])),reshape([nid_tags[v][d] for d in eachindex(nid_tags[1]) for v in basin],length(basin),length(nid_tags[1]))))
            end
            println("Computing attractor of basin $bi")
            probe = basin[1]
            path_to_att = saw(transition_graph,probe,length(basin))
            path_in_att = saw(transition_graph,path_to_att[end],length(basin))
            att_list[bi] = path_in_att
            periods[bi] = length(path_in_att)
            println("Attractor $bi has period ",periods[bi])
            patt_basin = Vector{Integer}()
            for point_in_att in path_in_att
                push!(patt_basin, findfirst(basin,point_in_att))
            end
            patt_basins[bi] = patt_basin
        end
        transition_graph = DiGraph()
        for (bi,basin_tree) in enumerate(basin_tree_list)
            println("Calculating branches associated to attractor $bi")
            edges_att = hcat(circshift(patt_basins[bi],1),patt_basins[bi])
            for pi in 1:length(patt_basins[bi])
                rem_edge!(basin_tree,edges_att[pi,:]...)
            end
            branch_node_list = weakly_connected_components(basin_tree)
            branch_tree_list[bi] = Vector{LightGraphs.DiGraph}(periods[bi])
            patt_to_look = deepcopy(patt_basins[bi])
            branch_patt[bi] = Vector{Integer}(periods[bi])
            for (bri,branch) in enumerate(branch_node_list)
                for point_in_att in patt_to_look
                    if point_in_att in branch
                        branch_patt[bi][bri] = findfirst(patt_basins[bi],point_in_att)
                        branch_tree_list[bi][bri] = induced_subgraph(basin_tree,branch)[1]
                        # println("Branch $bi:$bri has $(nv(branch_tree_list[bi][bri])) nodes")
                        deleteat!(patt_to_look,findfirst(patt_to_look,point_in_att))
                        break
                    end
                end
            end
        end
        basin_tree_list = Vector{LightGraphs.DiGraph}()
        for (bi,basin_tree) in enumerate(branch_tree_list)
            println("Computing transients distribution associated to attractor $bi")
            dists = Integer[]
            brsize = Vector{Integer}(periods[bi])
            for (bri,branch_tree) in enumerate(basin_tree)
                println("Measuring transient distances for branch $bi:$bri")
                dists_br = Integer[]
                if (trav==:rw)
                    maxbr_size = nv(branch_tree)
                    for pbi in vertices(branch_tree)
                      # push!(dists_br,length(a_star(basin_tree,pbranch,branch_patt[bri])))
                        push!(dists_br,length(non_backtracking_randomwalk(branch_tree,pbi,maxbr_size)-1))
                        if mod(pbi,100000)==0
                            @printf("%0.2f%% of branch %d:%d explored\n",100*pbi/maxbr_size, bi, bri)
                        end
                    end
                    filter!(x->x!=0,dists_br)
                else
                    m = mapreduce(x->in_neighbors(branch_tree,x),vcat,vertices(branch_tree))
                    d = sparsevec(m,ones(Integer,length(m)),nv(branch_tree))
                    l = 1
                    println("Level $l")
                    while !isempty(m)
                        m = mapreduce(x->in_neighbors(branch_tree,x),vcat,m)
                        d[m] += 1
                        l+=1
                        println("Level $l")
                    end
                    dists_br = nonzeros(d)
                end
                append!(dists,dists_br)
                brsize[branch_patt[bi][bri]] = length(dists_br)
            end
            branch_sizes[bi] = brsize
            writedlm(h,dists',',')
            println("Maximum transient distance:$(maximum(dists))")
        end
    finally
        close(g)
        close(h)
    end
    f = open(joinpath(out_dir,"$(root_name)_att_stats.txt"),"w")
    sum_att = sum(periods)
    sum_basin = state_space_size - sum_att
    try
        println("Writing statistics")
        for bi in eachindex(att_list)
            if bi > 1
                write(f,"\n")
            end
            write(f,"Attractor:$bi\t")
            period = periods[bi]
            write(f,"Period:$period\t")
            basin_size = mapreduce(x->nv(x),+,branch_tree_list[bi])-period
            write(f,"Basin size:$basin_size/$sum_basin ($(format(basin_size/sum_basin*100.0,precision=2,zeropadding=true))%)\n")
            if update_mode==:sy
                for (pi,conf_id) in enumerate(att_list[bi])
                    conf = reverse(digits(nid_tags[conf_id],tag_base,net_size))
                    write(f,"$(join(conf,","))\t$(branch_sizes[bi][pi])\n")
                end
            else
                for (pi,conf_id) in enumerate(att_list[bi])
                    conf = reverse(digits(nid_tags[conf_id][end],tag_base,net_size))
                    write(f,"$(join(conf,","))\t$(branch_sizes[bi][pi])\n")
                end
            end
        end
    finally
        close(f)
    end
    return basin_tags,basin_tree_list,att_list
    # return transition_graph
end

# get_attractor finds attractors exploring multiple random initial conditions (nrand_init) within a simulation
# time range limited by max_steps, after a transient given by transient_size
function find_attractor!(net::Net2{T},transient_size::Integer,max_steps::Integer,nrand_init::Integer;update_mode::Symbol=:sy,forcing_rhythms::Dict{Symbol,Tuple{Vector{Bool},Vector{Int}}}=Dict{Symbol,Tuple{Vector{Bool},Vector{Int}}}(),init_constraints...) where {T<:state_type}
    i = 0
    # tag_base = maximum([length(sigma.state_range) for sigma in net.nodes])
    net_hist_size = 1 + (update_mode==:sy ? 0 : maximum(net.in_taus))
    net_size = length(net.nodes)
    att_list = Vector{Vector{Tuple{Vararg{Tuple{Vararg{Int8,net_size}},net_hist_size}}}}()
    hX = Array{T}(net_size,net_hist_size)
    sim_trans_buffer = Array{T}(net_size,transient_size + net_hist_size)
    time_steps = 1
    sim_att_buffer = Array{T}(net_size,time_steps + net_hist_size)
    unsuccess = 0
    basin_sizes = Int[]
    sigma_ext_forcing = collect(sigma_symbol=>rhythm_values[2] for (sigma_symbol,rhythm_values) in forcing_rhythms)
	init_constraints = vcat(init_constraints,sigma_ext_forcing)
    if update_mode in [:roa,:ga]
        error("Not yet implemented for random asyncrhonous updating")
    end
    while i < nrand_init
        rand_init_net!(net;update_mode=update_mode,init_h=:rand,init_constraints...)
        init_probe = tuple(map(x->tuple(x...),(net.sim_buffer[:,ti] for ti in 1:net_hist_size))...)
        p_att = Vector{Tuple{Vararg{Tuple{Vararg{Int8,net_size}},net_hist_size}}}()
        hX[:] = evolve_net!(net,transient_size,keep=false,update_mode=update_mode,sim=sim_trans_buffer,forcing_rhythms=forcing_rhythms)
        probe = tuple(map(x->tuple(x...),(hX[:,ti] for ti in 1:net_hist_size))...)
        push!(p_att,probe)
        add_attractor = false
        already_found = false
        attractor_found = false
        tmax = 1
        for t in transient_size:transient_size+max_steps-1
            hX[:] = evolve_net!(net,1,keep=false,update_mode=update_mode,sim=sim_att_buffer,t_init=t,forcing_rhythms=forcing_rhythms)
            ntp1 = tuple(map(x->tuple(x...),(hX[:,ti] for ti in 1:net_hist_size))...)
            tmax = t
            if probe[end] == ntp1[end]
                # println("Likely attractor found")
                for (pi,att) in enumerate(att_list)
                    for p_att_i in p_att
                        if p_att_i in att
                            already_found = true
                            break
                        end
                    end
                    if already_found
                        if !(init_probe in p_att)
                            basin_sizes[pi]+=1
                        end
                        break
                    end
                end
                attractor_found = true
                if !already_found
                    push!(basin_sizes,1)
                    push!(att_list,p_att)
                    println("Attractor ",size(att_list,1)," has period ",size(p_att,1))
                    for att_i in p_att
                        println(join(att_i[end],","))
                    end
                end
                break
            end
            push!(p_att,ntp1)
        end
        i += 1
        if !(already_found||attractor_found)
            println("No attractor found After $tmax steps")
            print("probe:")
            println(join(probe[end],","))
            print("ntp1 :")
            println(join(hX[:,end],","))
            #  println("No attractors after reaching maxtimum time steps")
            unsuccess += 1
            # println("No attractor found for initial condition $i:")
            # for probe_i in probe
            #     conf = reverse(digits(probe_i,tag_base,net_size))
            #     println(join(conf,","))
            # end
        end
    end
    if length(att_list)!=0
        println("Estimated basin sizes:")
        for (pi,att) in enumerate(att_list)
            println("Attractor $pi:",format(100.0*basin_sizes[pi]/sum(basin_sizes),precision=2,zeropadding=true),"%")
        end
    end
    println(unsuccess," unsuccessful tries after reaching maximmum time steps set")
    return att_list
end

function find_attractor2!(net::Net2{T},max_steps::Integer,nrand_init::Integer;update_mode::Symbol=:sy,forcing_rhythms::Dict{Symbol,Tuple{Vector{Bool},Vector{Int}}}=Dict{Symbol,Tuple{Vector{Bool},Vector{Int}}}(),init_constraints...) where {T<:state_type}
    i = 0
    # tag_base = maximum([length(sigma.state_range) for sigma in net.nodes])
    net_hist_size = 1 + (update_mode==:sy ? 0 : maximum(net.in_taus))
    net_size = length(net.nodes)
    att_list = Vector{Vector{Tuple{Vararg{Tuple{Vararg{Int8,net_size}},net_hist_size}}}}()
    hX = Array{T}(net_size,net_hist_size)
    sim_att_buffer = Array{T}(net_size,1 + net_hist_size)
    unsuccess = 0
    basin_sizes = Int[]
    transients = Vector{Vector{Int}}()
    sigma_ext_forcing = collect(sigma_symbol=>rhythm_values[2] for (sigma_symbol,rhythm_values) in forcing_rhythms)
	init_constraints = vcat(init_constraints,sigma_ext_forcing)
    if update_mode in [:roa,:ga]
        error("Not yet implemented for random asyncrhonous updating")
    end
    while i < nrand_init
        rand_init_net!(net;update_mode=update_mode,init_h=:rand,init_constraints...)
        p_att = Vector{Tuple{Vararg{Tuple{Vararg{Int8,net_size}},net_hist_size}}}()
        init_probe = tuple(map(x->tuple(x...),(net.sim_buffer[:,ti] for ti in 1:net_hist_size))...)
        push!(p_att,init_probe)
        add_attractor = false
        already_found = false
        attractor_found = false
        tmax = 0
        for t in 0:max_steps-1
            hX[:] = evolve_net!(net,1,keep=false,update_mode=update_mode,sim=sim_att_buffer,t_init=t,forcing_rhythms=forcing_rhythms)
            ntp1 = tuple(map(x->tuple(x...),(hX[:,ti] for ti in 1:net_hist_size))...)
            tmax = t
            for nti in collect(size(p_att,1):-1:1)
                if ntp1[end] == p_att[nti][end]
                    #println("Likely attractor found:$i at t=$t")
                    # Check if it was previously found
                    for (pi,att) in enumerate(att_list)
                        for p_att_i in p_att
                            if p_att_i in att
                                already_found = true
                                if !(init_probe in att)
                                    basin_sizes[pi]+=1
                                    push!(transients[pi],nti-1)
                                end
                                break
                            end
                        end
                        if already_found
                            break
                        end
                    end
                    if !already_found
                        deleteat!(p_att,1:nti-1)
                        push!(att_list,p_att)
                        push!(basin_sizes,init_probe in p_att ? 0 : 1)
                        push!(transients,[nti-1])
                        println("Attractor ",size(att_list,1)," has period ",size(p_att,1))
                        for att_i in p_att
                            println(join(att_i[end],","))
                        end
                    end
                    attractor_found = true
                    break
                end
            end
            if attractor_found
                #println(i," : ",nti-1)
                break
            end
            push!(p_att,ntp1)
        end
        i += 1
        if !(already_found||attractor_found)
            println("No attractor found After $tmax steps")
            print("probe:")
            println(join(init_probe[end],","))
            #  println("No attractors after reaching maxtimum time steps")
            unsuccess += 1
            # println("No attractor found for initial condition $i:")
            # for probe_i in probe
            #     conf = reverse(digits(probe_i,tag_base,net_size))
            #     println(join(conf,","))
            # end
        end
    end
    if length(att_list)!=0
        println("Estimated basin sizes:")
        for (pi,att) in enumerate(att_list)
            println("Attractor $pi:",format(100.0*basin_sizes[pi]/sum(basin_sizes),precision=2,zeropadding=true),"%")
        end
    end
    println(unsuccess," unsuccessful tries after reaching maximmum time steps set")
    return att_list, transients
end

function read_metadata(filepath::AbstractString,pars::Vector{Symbol})
    param_found = falses(length(pars))
    if isfile(filepath)
		println("\nLoading metafile ",filepath)
        par_dict = Dict(map(x->x=>0,pars))
		open(filepath,"r") do f
			for l in eachline(f)
                for (pi,p) in enumerate(pars)
				    if ismatch(Regex("^$p = \\d+"),l)
                        m = match(r"\d+",l)
                        par_dict[p] = parse(Int,m.match)
                        param_found[pi] = true
                    end
                    !(all(param_found)) || break
                end
            end
        end
    end
    if all(param_found)
        println("\nSucessfully read parameters ",join(map(string,pars),","))
        return par_dict
    else
        return nothing
    end
end


text:
module ReferenceTests

using Test
using Images
using FileIO
using ImageInTerminal
using ColorTypes
using SHA
using DeepDiffs
using Random

export
    @withcolor,
    @io2str,
    @test_reference

include("utils.jl")
include("test_reference.jl")
include("core.jl")
include("handlers.jl")

end # module


text:
#=
This file is part of the replication code for: Hasenzagl, T., Pellegrino, F., Reichlin, L., & Ricco, G. (2020). A Model of the Fed's View on Inflation.
Please cite the paper if you are using any part of the code for academic work (including, but not limited to, conference and peer-reviewed papers).
=#

function kalman_diffuse!(par::ParSsm, do_loglik=0::Int64, do_smoother=0::Int64, do_sim_smoother=0::Int64, F_tol=1e-8::Float64)

# ----------------------------------------------------------------------------------------------------------------------
# Diffuse Kalman Filter, Smoother and Simulation Smoother (univariate approach for multivariate models).
# The function update the loglikelihood in par if do_loglik is 1.
#
# The model is specified as follows:
#
# y_{t} = d + Z*a_{t}   + e_{t}, where e_{t} ~ N(0, R)
# α_{t} = c + T*α_{t-1} + u_{t}, where u_{t} ~ N(0, Q)
#
# ----------------------------------------------------------------------------------------------------------------------
#
# Kalman routines for non-stationary systems with time-varying system matrices and missing data. The Durbin and
# Koopman's simulation smoother includes the Jarocinski (2015) modifications. The univariate approach for multivariate
# models introduced by Anderson and Moore (1979) and extended by Durbin and Koopman (2012) is used to improve the
# estimation speed.
#
# See: Harvey Koopman and Shephard (2004) or Durbin and Koopman (2012)
#
# ----------------------------------------------------------------------------------------------------------------------
# Input arguments:
#
# y       observed data                                              (n x m)
# d       constant of the observation equation                       (n x 1)
# Z       coefficients of the observation equation                   (n x k)
# R       variance-covariance matrix of the observation equation     (n x n)
# c       constant of the transition equation                        (k x 1)
# T       coefficients of the transition equation                    (k x k)
# Q       variance-covariance matrix of the transition equation      (k x k)
# α¹      initial value of the states                                (k x 1)
# P¹      initial variance-covariance matrix of the states           (k x k)
# P̄¹      diagonal matrix: 1 for non-stationary states (0 otherwise) (k x k)
# F_tol   small value used to determine whether F is zero or not     (1 x 1)
#
# ----------------------------------------------------------------------------------------------------------------------
#
# Optional input arguments:
#
# loglik_bool    enable the estimation of the loglikelihood        (1 x 1)
# sim_smoother   enable the simulation smoother                    (1 x 1)
#
# ----------------------------------------------------------------------------------------------------------------------
#
# Output
# α̃              smoothed state vector                             (k x m)
# Pˢ             smoothed variance-covariance matrix of the states (k x k)
# loglik         loglikelihood (if enabled)                        (1 x 1)
# ----------------------------------------------------------------------------------------------------------------------

     # -----------------------------------------------------------------------------------------------------------------
     # Initialisation
     # -----------------------------------------------------------------------------------------------------------------

     # Dimensions
     k    = size(par.T)[1]
     n, m = size(par.y);

     # Kalman filter output
     ν  = zeros(m, n);
     K  = zeros(k, n, m);
     K̄  = zeros(k, n, m);
     F  = zeros(m, n);
     F̄  = zeros(m, n);
     α̂ᶠ = zeros(k, m);
     Pᶠ = zeros(k, k, m);
     P̄ᶠ = zeros(k, k, m);

     # Kalman smoother output
     α̂ˢ = zeros(k, m);
     Pˢ = zeros(k, k, m);
     P̄ˢ = zeros(k, k, m);

     # Loglikelihood
     if do_loglik == 1
          par.loglik = 0;
     end

     if do_sim_smoother == 1
          do_smoother = 1;
     end

     # Initial variance
     par.P¹[par.P̄¹ .== 1] .= 0;

     # -----------------------------------------------------------------------------------------------------------------
     # Simulation smoother (if enabled)
     # -----------------------------------------------------------------------------------------------------------------

     if do_sim_smoother == 1

          # Draw: y⁺ and α⁺

          y⁺ = zeros(n, m);
          α⁺ = zeros(k, m);

          # ------------------------------------------------------------------------------------------------------------
          # Durbin, Koopman 2002, Algorithm 2.
          #
          # Generate y⁺ and α⁺ - y and α are drawn from their unconditional distributions using the demeaned model,
          # i.e. zero initial state and zero constant terms.
          #
          # Note: Jarocinski (2015) modification
          # ------------------------------------------------------------------------------------------------------------

          for t=1:m
               if t == 1
                    # draw draw the first state with α1=0
                    α⁺[:, 1] = get_random_disturbance(k, par.P¹);
               else
                    α⁺[:, t] = par.T*α⁺[:, t-1] + get_random_disturbance(k, par.Q);
               end

               y⁺[:, t] = par.Z*α⁺[:, t] + get_random_disturbance(n, par.R);
          end

          # Define antithetic variable yᵃ
          yᵃ = par.y - y⁺;

     else
          yᵃ = par.y;
     end

     # -----------------------------------------------------------------------------------------------------------------
     # In order to use the univariate approach we must deal with correlated measurement errors (if any)
     # -----------------------------------------------------------------------------------------------------------------

     # 1. Check if there are correlated measurement errors
     R_ind               = par.R .!= 0;
     R_ind[Matrix(I,n,n) .== 1] .= false;

     # 2. If there are correlated measurement errors, apply the following Shur decomposition R = S⋅H⋅S', where H is
     #    diagonal and holds the eigenvalues of R, while S is orthogonal. Pre-multiply the observation equation by S',
     #    so that the variance of the errors will be diagonal and equal to H
     if sum(sum(R_ind)) > 0
          H, S, _ = schur(copy(par.R));
          yᵃ = S'*yᵃ;
          d̃  = S'*copy(par.d);
          Z̃  = S'*copy(par.Z);
          R̃  = H;
     else
          d̃ = par.d;
          Z̃ = par.Z;
          R̃ = par.R;
     end

     # -----------------------------------------------------------------------------------------------------------------
     # Kalman Filter - Univariate (or sequential) filtering
     # -----------------------------------------------------------------------------------------------------------------

     α̂ᶠᵢ = par.α¹;
     Pᶠᵢ = par.P¹;
     P̄ᶠᵢ = par.P̄¹;

     α̂ᶠ[:, 1]    = α̂ᶠᵢ;
     Pᶠ[:, :, 1] = Pᶠᵢ;
     P̄ᶠ[:, :, 1] = P̄ᶠᵢ;

     # Loop over t
     for t=1:m

          # Get last estimates for α and P
          α̂ᶠᵢ = α̂ᶠ[:, t];
          Pᶠᵢ = Pᶠ[:, :, t];
          P̄ᶠᵢ = P̄ᶠ[:, :, t];

          # Handle missing data
          no_na     = ismissing.(yᵃ[:, t]) .== false;
          ind_no_na = findall(no_na);
          yᵃᵗ       = convert(Array{Float64, 1}, yᵃ[no_na, t]);
          d̃ᵗ        = d̃[no_na, :];
          Z̃ᵗ        = Z̃[no_na, :];
          R̃ᵗ        = R̃[no_na, no_na];

          # Loop over not missing data
          for i=1:sum(no_na)

               # Kalman filter i-th output
               Fᵗᵢ = (Z̃ᵗ[i, :]'*Pᶠᵢ*Z̃ᵗ[i, :] + R̃ᵗ[i, i])[1];
               F̄ᵗᵢ = (Z̃ᵗ[i, :]'*P̄ᶠᵢ*Z̃ᵗ[i, :])[1];

               if F̄ᵗᵢ > F_tol || Fᵗᵢ > F_tol
                    Kᵗᵢ = Pᶠᵢ*Z̃ᵗ[i, :];
                    K̄ᵗᵢ = P̄ᶠᵢ*Z̃ᵗ[i, :];
                    νᵗᵢ = (yᵃᵗ[i, :] .- d̃ᵗ[i, :] .- Z̃ᵗ[i, :]'*α̂ᶠᵢ)[1];

                    # --------------------------------------------------------------------------------------------------
                    # Case 1: F̄ᵗᵢ greater than zero
                    # --------------------------------------------------------------------------------------------------

                    if F̄ᵗᵢ > F_tol

                         if t < m
                              # Update α̂ᶠᵢ, Pᶠᵢ and P̄ᶠᵢ
                              α̂ᶠᵢ = α̂ᶠᵢ + K̄ᵗᵢ/F̄ᵗᵢ*νᵗᵢ;
                              Pᶠᵢ = Pᶠᵢ + K̄ᵗᵢ*K̄ᵗᵢ'*Fᵗᵢ/(F̄ᵗᵢ^2) - (Kᵗᵢ*K̄ᵗᵢ'+K̄ᵗᵢ*Kᵗᵢ')/F̄ᵗᵢ;
                              P̄ᶠᵢ = P̄ᶠᵢ - K̄ᵗᵢ/F̄ᵗᵢ*K̄ᵗᵢ';
                         end

                         # Estimate loglikelihood  (if enabled)
                         if do_loglik == 1
                              par.loglik = par.loglik -0.5*(log(2*pi) + log(F̄ᵗᵢ));
                         end

                    # --------------------------------------------------------------------------------------------------
                    # Case 2: F̄ᵗᵢ is equal zero
                    # --------------------------------------------------------------------------------------------------

                    else

                         if t < m
                              # Update α̂ᶠᵢ, Pᶠᵢ and P̄ᶠᵢ (the latter doesn't change in this case)
                              α̂ᶠᵢ = α̂ᶠᵢ + Kᵗᵢ/Fᵗᵢ*νᵗᵢ;
                              Pᶠᵢ = Pᶠᵢ - Kᵗᵢ/Fᵗᵢ*Kᵗᵢ';
                         end

                         # Estimate loglikelihood  (if enabled)
                         if do_loglik == 1
                              par.loglik = par.loglik -0.5*(log(2*pi) + log(Fᵗᵢ) + (νᵗᵢ^2)/Fᵗᵢ);
                         end
                    end

                    # Store Kalman filter output (1/2): i-th values
                    if do_smoother == 1
                         ν[t, ind_no_na[i]]    = νᵗᵢ;
                         K[:, ind_no_na[i], t] = Kᵗᵢ;
                         K̄[:, ind_no_na[i], t] = K̄ᵗᵢ;
                         F[t, ind_no_na[i]]    = Fᵗᵢ;
                         F̄[t, ind_no_na[i]]    = F̄ᵗᵢ;
                    end
               end

               if t < m
                    # Store Kalman filter output (2/2): t-th values
                    α̂ᶠ[:, t+1]    = par.c + par.T*α̂ᶠᵢ;
                    Pᶠ[:, :, t+1] = par.T*Pᶠᵢ*par.T' + par.Q;
                    P̄ᶠ[:, :, t+1] = par.T*P̄ᶠᵢ*par.T';
               end
          end

          if sum(no_na) == 0 && t < m

               # Store Kalman filter predictions
               α̂ᶠ[:, t+1]    = par.c + par.T*α̂ᶠᵢ;
               Pᶠ[:, :, t+1] = par.T*Pᶠᵢ*par.T' + par.Q;
               P̄ᶠ[:, :, t+1] = par.T*P̄ᶠᵢ*par.T';
          end
     end

     # -----------------------------------------------------------------------------------------------------------------
     # Kalman Smoother - Univariate (or sequential) smoother
     # -----------------------------------------------------------------------------------------------------------------

     if do_smoother == 1

          rᵢ = zeros(2*k);
          Nᵢ = zeros(2*k, 2*k);
          TT = kron(Matrix(I,2,2), par.T);

          # Loop over t
          for t=m:-1:1

               if t < m
                    rᵢ = TT'*rᵢ;
                    Nᵢ = TT'*Nᵢ*TT;
               end

               # Handle missing data
               no_na     = ismissing.(yᵃ[:, t]) .== false;
               ind_no_na = findall(no_na);
               yᵃᵗ       = convert(Array{Float64, 1}, yᵃ[no_na, t]);
               Z̃ᵗ        = Z̃[no_na, :];

               # Loop over not missing data
               for i=sum(no_na):-1:1

                    # Retrieve Kalman Filter output (1/2)
                    Fᵗᵢ = F[t, ind_no_na[i]];
                    F̄ᵗᵢ = F̄[t, ind_no_na[i]];

                    if F̄ᵗᵢ > F_tol || Fᵗᵢ > F_tol

                         # Retrieve Kalman Filter output (2/2)
                         νᵗᵢ = ν[t, ind_no_na[i]];
                         Kᵗᵢ = K[:, ind_no_na[i], t];
                         K̄ᵗᵢ = K̄[:, ind_no_na[i], t];

                         # ---------------------------------------------------------------------------------------------
                         # Case 1: F̄ᵗᵢ greater than zero
                         # ---------------------------------------------------------------------------------------------

                         if F̄ᵗᵢ > F_tol

                              # Estimate Lᵗᵢ, L̄ᵗᵢ and Mᵗᵢ using the Kalman Filter output
                              Lᵗᵢ = (K̄ᵗᵢ*Fᵗᵢ/F̄ᵗᵢ-Kᵗᵢ)*Z̃ᵗ[i, :]'/F̄ᵗᵢ;
                              L̄ᵗᵢ = Matrix(I,k,k) - K̄ᵗᵢ*Z̃ᵗ[i, :]'/F̄ᵗᵢ;
                              Mᵗᵢ = [L̄ᵗᵢ Lᵗᵢ; zeros(size(L̄ᵗᵢ)) L̄ᵗᵢ];

                              # Estimate rᵢ and Nᵢ
                              temp¹ = Z̃ᵗ[i, :]/F̄ᵗᵢ;
                              temp² = temp¹*Z̃ᵗ[i, :]';
                              rᵢ    = [zeros(k); temp¹*νᵗᵢ] + Mᵗᵢ'*rᵢ;
                              Nᵢ    = [zeros(k, k) temp²; temp² temp²*Fᵗᵢ/F̄ᵗᵢ] + Mᵗᵢ'*Nᵢ*Mᵗᵢ;

                         # ---------------------------------------------------------------------------------------------
                         # Case 2: F̄ᵗᵢ is equal zero
                         # ---------------------------------------------------------------------------------------------

                         else
                              # Estimate Lᵗᵢ, L̄ᵗᵢ and LMᵗᵢ using the Kalman Filter output
                              Lᵗᵢ = Matrix(I,k,k) - Kᵗᵢ*Z̃ᵗ[i, :]'/Fᵗᵢ;
                              Mᵗᵢ = [Lᵗᵢ zeros(size(Lᵗᵢ)); zeros(size(Lᵗᵢ)) Lᵗᵢ];

                              # Estimate rᵢ and Nᵢ
                              temp¹ = Z̃ᵗ[i, :]/Fᵗᵢ;
                              rᵢ    = [temp¹*νᵗᵢ; zeros(k)] + Mᵗᵢ'*rᵢ;
                              Nᵢ    = [temp¹*Z̃ᵗ[i, :]' zeros(k, k); zeros(k, 2*k)] + Mᵗᵢ'*Nᵢ*Mᵗᵢ;
                         end
                    end
               end

               # Smoothed states and covariance
               P̃ᶠ          = [Pᶠ[:, :, t] P̄ᶠ[:, :, t]];
               α̂ˢ[:, t]    = α̂ᶠ[:, t] + P̃ᶠ*rᵢ;
               Pˢ[:, :, t] = Pᶠ[:, :, t] - P̃ᶠ*Nᵢ*P̃ᶠ';
          end
     end

     # -----------------------------------------------------------------------------------------------------------------
     # Return output
     # -----------------------------------------------------------------------------------------------------------------

     if do_sim_smoother == 1
          α = α̂ˢ + α⁺;
          P = Pˢ;

     elseif do_smoother == 1
          α = α̂ˢ;
          P = Pˢ;

     else
          α = α̂ᶠ;
          P = Pᶠ;
     end

     return α, P;
end


text:
using epilepsy: sample_entropy, multiscale_entropy
using NPZ: npzread

eeg = npzread("test_data/eeg_processed.npy")
eeg_segment = eeg[1, 1:10000]

epoched_eeg = [[1, 2, 3, 4, 5, 4], [1, 1, 1, 1, 1, 1], [1, 2, 1, 2, 1, 2]]
sine = sin.(collect(1:0.1:1000))

m = 2
r = 0.2
δ = 1
τ = 3

@test sample_entropy(ones(10), m, r) == (0.0, 28, 28)
@test sample_entropy(sine, m, r)[1] <= 0.5
@test multiscale_entropy(sine, 1, m, r) == sample_entropy(sine, m, r)


text:
using MKtest
using Test

# Running MKtest approximation with default parameters
@testset "MKtest.jl" begin

	adap  = MKtest.parameters(N=1000,n=661,al_low=0.2,B_bins=[0.2,0.999])
	binom = MKtest.binom_op!(adap)
	MKtest.set_θ!(adap)
	MKtest.set_ppos!(adap)

	@test adap.ppos_h == 0.0002072562337323172
	@test adap.θ_noncoding == 1.6433217979109007e-6

end