text:
"""
     FciqmcRunStrategy{T}
Abstract type representing the strategy for running and terminating
[`fciqmc!()`](@ref). The type parameter `T` is relevant for reporting the shift
and the norm.

Implemented strategies:

   * [`RunTillLastStep`](@ref)
"""
abstract type FciqmcRunStrategy{T} end

"""
    RunTillLastStep(step::Int = 0 # number of current/starting timestep
                 laststep::Int = 100 # number of final timestep
                 shiftMode::Bool = false # whether to adjust shift
                 shift = 0.0 # starting/current value of shift
                 dτ::Float64 = 0.01 # current value of time step
    ) <: FciqmcRunStrategy
Parameters for running [`fciqmc!()`](@ref) for a fixed number of time steps.
For alternative strategies, see [`FciqmcRunStrategy`](@ref).
"""
@with_kw mutable struct RunTillLastStep{T} <: FciqmcRunStrategy{T}
    step::Int = 0 # number of current/starting timestep
    laststep::Int = 100 # number of final timestep
    shiftMode::Bool = false # whether to adjust shift
    shift::T = 0.0 # starting/current value of shift
    dτ::Float64 = 0.01 # time step
end


text:
@cenum ELFSegment::UInt32 begin
   SEGMENT_GLOBAL_PROGRAM = 0
   SEGMENT_GLOBAL_AGENT = 1
   SEGMENT_READONLY_AGENT = 2
   SEGMENT_CODE_AGENT = 3
   SEGMENT_LAST = 4
end

@cenum SymbolAllocation::UInt32 begin
   SYMBOL_ALLOCATION_DEFAULT = 0
   SYMBOL_ALLOCATION_GLOBAL_PROGRAM = 1
   SYMBOL_ALLOCATION_GLOBAL_AGENT = 2
   SYMBOL_ALLOCATION_READONLY_AGENT = 3
   SYMBOL_ALLOCATION_LAST = 4
end

@cenum SymbolFlag::UInt32 begin
   SYMBOL_FLAG_DEFAULT = 0
   SYMBOL_FLAG_CONST = 1
   SYMBOL_FLAG_LAST = 2
end

@cenum MetadataKind::UInt16 begin
   METADATA_KIND_NONE = 0
   METADATA_KIND_INIT_SAMP = 1
   METADATA_KIND_INIT_ROIMG = 2
   METADATA_KIND_INIT_WOIMG = 3
   METADATA_KIND_INIT_RWIMG = 4
end

@cenum SamplerCoord::UInt8 begin
   SAMPLER_COORD_UNNORMALIZED = 0
   SAMPLER_COORD_NORMALIZED = 1
end

@cenum SamplerFilter::UInt8 begin
   SAMPLER_FILTER_NEAREST = 0
   SAMPLER_FILTER_LINEAR = 1
end

@cenum SamplerAddressing::UInt8 begin
   SAMPLER_ADDRESSING_UNDEFINED = 0
   SAMPLER_ADDRESSING_CLAMP_TO_EDGE = 1
   SAMPLER_ADDRESSING_CLAMP_TO_BORDER = 2
   SAMPLER_ADDRESSING_REPEAT = 3
   SAMPLER_ADDRESSING_MIRRORED_REPEAT = 4
end

@cenum ImageGeometry::UInt8 begin
   IMAGE_GEOMETRY_1D = 0
   IMAGE_GEOMETRY_2D = 1
   IMAGE_GEOMETRY_3D = 2
   IMAGE_GEOMETRY_1DA = 3
   IMAGE_GEOMETRY_2DA = 4
   IMAGE_GEOMETRY_1DB = 5
   IMAGE_GEOMETRY_2DDEPTH = 6
   IMAGE_GEOMETRY_2DADEPTH = 7
end

@cenum ImageChannelOrder::UInt8 begin
   IMAGE_CHANNEL_ORDER_A = 0
   IMAGE_CHANNEL_ORDER_R = 1
   IMAGE_CHANNEL_ORDER_RX = 2
   IMAGE_CHANNEL_ORDER_RG = 3
   IMAGE_CHANNEL_ORDER_RGX = 4
   IMAGE_CHANNEL_ORDER_RA = 5
   IMAGE_CHANNEL_ORDER_RGB = 6
   IMAGE_CHANNEL_ORDER_RGBX = 7
   IMAGE_CHANNEL_ORDER_RGBA = 8
   IMAGE_CHANNEL_ORDER_BGRA = 9
   IMAGE_CHANNEL_ORDER_ARGB = 10
   IMAGE_CHANNEL_ORDER_ABGR = 11
   IMAGE_CHANNEL_ORDER_SRGB = 12
   IMAGE_CHANNEL_ORDER_SRGBX = 13
   IMAGE_CHANNEL_ORDER_SRGBA = 14
   IMAGE_CHANNEL_ORDER_SBGRA = 15
   IMAGE_CHANNEL_ORDER_INTENSITY = 16
   IMAGE_CHANNEL_ORDER_LUMINANCE = 17
   IMAGE_CHANNEL_ORDER_DEPTH = 18
   IMAGE_CHANNEL_ORDER_DEPTH_STENCIL = 19
end

@cenum ELFSection::UInt32 begin
   RODATA_GLOBAL_PROGRAM = 0
   RODATA_GLOBAL_AGENT = 1
   RODATA_READONLY_AGENT = 2
   DATA_GLOBAL_PROGRAM = 3
   DATA_GLOBAL_AGENT = 4
   DATA_READONLY_AGENT = 5
   BSS_GLOBAL_PROGRAM = 6
   BSS_GLOBAL_AGENT = 7
   BSS_READONLY_AGENT = 8
   SECTION_LAST = 9
end

@cenum Status::UInt32 begin
    STATUS_SUCCESS = 0
    STATUS_INFO_BREAK = 1
    STATUS_ERROR = 4096
    STATUS_ERROR_INVALID_ARGUMENT = 4097
    STATUS_ERROR_INVALID_QUEUE_CREATION = 4098
    STATUS_ERROR_INVALID_ALLOCATION = 4099
    STATUS_ERROR_INVALID_AGENT = 4100
    STATUS_ERROR_INVALID_REGION = 4101
    STATUS_ERROR_INVALID_SIGNAL = 4102
    STATUS_ERROR_INVALID_QUEUE = 4103
    STATUS_ERROR_OUT_OF_RESOURCES = 4104
    STATUS_ERROR_INVALID_PACKET_FORMAT = 4105
    STATUS_ERROR_RESOURCE_FREE = 4106
    STATUS_ERROR_NOT_INITIALIZED = 4107
    STATUS_ERROR_REFCOUNT_OVERFLOW = 4108
    STATUS_ERROR_INCOMPATIBLE_ARGUMENTS = 4109
    STATUS_ERROR_INVALID_INDEX = 4110
    STATUS_ERROR_INVALID_ISA = 4111
    STATUS_ERROR_INVALID_ISA_NAME = 4119
    STATUS_ERROR_INVALID_CODE_OBJECT = 4112
    STATUS_ERROR_INVALID_EXECUTABLE = 4113
    STATUS_ERROR_FROZEN_EXECUTABLE = 4114
    STATUS_ERROR_INVALID_SYMBOL_NAME = 4115
    STATUS_ERROR_VARIABLE_ALREADY_DEFINED = 4116
    STATUS_ERROR_VARIABLE_UNDEFINED = 4117
    STATUS_ERROR_EXCEPTION = 4118
    STATUS_ERROR_INVALID_CODE_SYMBOL = 4120
    STATUS_ERROR_INVALID_EXECUTABLE_SYMBOL = 4121
    STATUS_ERROR_INVALID_FILE = 4128
    STATUS_ERROR_INVALID_CODE_OBJECT_READER = 4129
    STATUS_ERROR_INVALID_CACHE = 4130
    STATUS_ERROR_INVALID_WAVEFRONT = 4131
    STATUS_ERROR_INVALID_SIGNAL_GROUP = 4132
    STATUS_ERROR_INVALID_RUNTIME_STATE = 4133
    STATUS_ERROR_FATAL = 4134
end

@cenum AccessPermission::UInt32 begin
    ACCESS_PERMISSION_RO = 1
    ACCESS_PERMISSION_WO = 2
    ACCESS_PERMISSION_RW = 3
end

@cenum Endianness::UInt32 begin
    ENDIANNESS_LITTLE = 0
    ENDIANNESS_BIG = 1
end

@cenum MachineModel::UInt32 begin
    MACHINE_MODEL_SMALL = 0
    MACHINE_MODEL_LARGE = 1
end

@cenum Profile::UInt32 begin
    PROFILE_BASE = 0
    PROFILE_FULL = 1
end

@cenum SystemInfo::UInt32 begin
    SYSTEM_INFO_VERSION_MAJOR = 0
    SYSTEM_INFO_VERSION_MINOR = 1
    SYSTEM_INFO_TIMESTAMP = 2
    SYSTEM_INFO_TIMESTAMP_FREQUENCY = 3
    SYSTEM_INFO_SIGNAL_MAX_WAIT = 4
    SYSTEM_INFO_ENDIANNESS = 5
    SYSTEM_INFO_MACHINE_MODEL = 6
    SYSTEM_INFO_EXTENSIONS = 7
    AMD_SYSTEM_INFO_BUILD_VERSION = 512
end

@cenum Extension::UInt32 begin
    EXTENSION_FINALIZER = 0
    EXTENSION_IMAGES = 1
    EXTENSION_PERFORMANCE_COUNTERS = 2
    EXTENSION_PROFILING_EVENTS = 3
    EXTENSION_STD_LAST = 3
    AMD_FIRST_EXTENSION = 512
    EXTENSION_AMD_PROFILER = 512
    EXTENSION_AMD_LOADER = 513
    EXTENSION_AMD_AQLPROFILE = 514
    AMD_LAST_EXTENSION = 514
end

@cenum AgentFeature::UInt32 begin
    AGENT_FEATURE_KERNEL_DISPATCH = 1
    AGENT_FEATURE_AGENT_DISPATCH = 2
end

@cenum DeviceType::UInt32 begin
    DEVICE_TYPE_CPU = 0
    DEVICE_TYPE_GPU = 1
    DEVICE_TYPE_DSP = 2
end

@cenum DefaultFloatRoundingMode::UInt32 begin
    DEFAULT_FLOAT_ROUNDING_MODE_DEFAULT = 0
    DEFAULT_FLOAT_ROUNDING_MODE_ZERO = 1
    DEFAULT_FLOAT_ROUNDING_MODE_NEAR = 2
end

@cenum AgentInfo::UInt32 begin
    AGENT_INFO_NAME = 0
    AGENT_INFO_VENDOR_NAME = 1
    AGENT_INFO_FEATURE = 2
    AGENT_INFO_MACHINE_MODEL = 3
    AGENT_INFO_PROFILE = 4
    AGENT_INFO_DEFAULT_FLOAT_ROUNDING_MODE = 5
    AGENT_INFO_BASE_PROFILE_DEFAULT_FLOAT_ROUNDING_MODES = 23
    AGENT_INFO_FAST_F16_OPERATION = 24
    AGENT_INFO_WAVEFRONT_SIZE = 6
    AGENT_INFO_WORKGROUP_MAX_DIM = 7
    AGENT_INFO_WORKGROUP_MAX_SIZE = 8
    AGENT_INFO_GRID_MAX_DIM = 9
    AGENT_INFO_GRID_MAX_SIZE = 10
    AGENT_INFO_FBARRIER_MAX_SIZE = 11
    AGENT_INFO_QUEUES_MAX = 12
    AGENT_INFO_QUEUE_MIN_SIZE = 13
    AGENT_INFO_QUEUE_MAX_SIZE = 14
    AGENT_INFO_QUEUE_TYPE = 15
    AGENT_INFO_NODE = 16
    AGENT_INFO_DEVICE = 17
    AGENT_INFO_CACHE_SIZE = 18
    AGENT_INFO_ISA = 19
    AGENT_INFO_EXTENSIONS = 20
    AGENT_INFO_VERSION_MAJOR = 21
    AGENT_INFO_VERSION_MINOR = 22
end

@cenum ExceptionPolicy::UInt32 begin
    EXCEPTION_POLICY_BREAK = 1
    EXCEPTION_POLICY_DETECT = 2
end

@cenum CacheInfo::UInt32 begin
    CACHE_INFO_NAME_LENGTH = 0
    CACHE_INFO_NAME = 1
    CACHE_INFO_LEVEL = 2
    CACHE_INFO_SIZE = 3
end

@cenum SignalCondition::UInt32 begin
    SIGNAL_CONDITION_EQ = 0
    SIGNAL_CONDITION_NE = 1
    SIGNAL_CONDITION_LT = 2
    SIGNAL_CONDITION_GTE = 3
end

@cenum WaitState::UInt32 begin
    WAIT_STATE_BLOCKED = 0
    WAIT_STATE_ACTIVE = 1
end

@cenum QueueType::UInt32 begin
    QUEUE_TYPE_MULTI = 0
    QUEUE_TYPE_SINGLE = 1
end

@cenum QueueFeature::UInt32 begin
    QUEUE_FEATURE_KERNEL_DISPATCH = 1
    QUEUE_FEATURE_AGENT_DISPATCH = 2
end

@cenum PacketType::UInt32 begin
    PACKET_TYPE_VENDOR_SPECIFIC = 0
    PACKET_TYPE_INVALID = 1
    PACKET_TYPE_KERNEL_DISPATCH = 2
    PACKET_TYPE_BARRIER_AND = 3
    PACKET_TYPE_AGENT_DISPATCH = 4
    PACKET_TYPE_BARRIER_OR = 5
end

@cenum FenceScope::UInt32 begin
    FENCE_SCOPE_NONE = 0
    FENCE_SCOPE_AGENT = 1
    FENCE_SCOPE_SYSTEM = 2
end

@cenum PacketHeader::UInt32 begin
    PACKET_HEADER_TYPE = 0
    PACKET_HEADER_BARRIER = 8
    PACKET_HEADER_SCACQUIRE_FENCE_SCOPE = 9
    PACKET_HEADER_ACQUIRE_FENCE_SCOPE = 9
    PACKET_HEADER_SCRELEASE_FENCE_SCOPE = 11
    PACKET_HEADER_RELEASE_FENCE_SCOPE = 11
end

@cenum PacketHeaderWidth::UInt32 begin
    PACKET_HEADER_WIDTH_TYPE = 8
    PACKET_HEADER_WIDTH_BARRIER = 1
    PACKET_HEADER_WIDTH_SCACQUIRE_FENCE_SCOPE = 2
    PACKET_HEADER_WIDTH_ACQUIRE_FENCE_SCOPE = 2
    PACKET_HEADER_WIDTH_SCRELEASE_FENCE_SCOPE = 2
    PACKET_HEADER_WIDTH_RELEASE_FENCE_SCOPE = 2
end

@cenum KernelDispatchPacketSetup::UInt32 begin
    KERNEL_DISPATCH_PACKET_SETUP_DIMENSIONS = 0
end

@cenum KernelDispatchPacketSetupWidth::UInt32 begin
    KERNEL_DISPATCH_PACKET_SETUP_WIDTH_DIMENSIONS = 2
end

@cenum RegionSegment::UInt32 begin
    REGION_SEGMENT_GLOBAL = 0
    REGION_SEGMENT_READONLY = 1
    REGION_SEGMENT_PRIVATE = 2
    REGION_SEGMENT_GROUP = 3
    REGION_SEGMENT_KERNARG = 4
end

@cenum RegionGlobalFlag::UInt32 begin
    REGION_GLOBAL_FLAG_KERNARG = 1
    REGION_GLOBAL_FLAG_FINE_GRAINED = 2
    REGION_GLOBAL_FLAG_COARSE_GRAINED = 4
end

@cenum RegionInfo::UInt32 begin
    REGION_INFO_SEGMENT = 0
    REGION_INFO_GLOBAL_FLAGS = 1
    REGION_INFO_SIZE = 2
    REGION_INFO_ALLOC_MAX_SIZE = 4
    REGION_INFO_ALLOC_MAX_PRIVATE_WORKGROUP_SIZE = 8
    REGION_INFO_RUNTIME_ALLOC_ALLOWED = 5
    REGION_INFO_RUNTIME_ALLOC_GRANULE = 6
    REGION_INFO_RUNTIME_ALLOC_ALIGNMENT = 7
end

@cenum ISAInfo::UInt32 begin
    ISA_INFO_NAME_LENGTH = 0
    ISA_INFO_NAME = 1
    ISA_INFO_CALL_CONVENTION_COUNT = 2
    ISA_INFO_CALL_CONVENTION_INFO_WAVEFRONT_SIZE = 3
    ISA_INFO_CALL_CONVENTION_INFO_WAVEFRONTS_PER_COMPUTE_UNIT = 4
    ISA_INFO_MACHINE_MODELS = 5
    ISA_INFO_PROFILES = 6
    ISA_INFO_DEFAULT_FLOAT_ROUNDING_MODES = 7
    ISA_INFO_BASE_PROFILE_DEFAULT_FLOAT_ROUNDING_MODES = 8
    ISA_INFO_FAST_F16_OPERATION = 9
    ISA_INFO_WORKGROUP_MAX_DIM = 12
    ISA_INFO_WORKGROUP_MAX_SIZE = 13
    ISA_INFO_GRID_MAX_DIM = 14
    ISA_INFO_GRID_MAX_SIZE = 16
    ISA_INFO_FBARRIER_MAX_SIZE = 17
end

@cenum FPType::UInt32 begin
    FP_TYPE_16 = 1
    FP_TYPE_32 = 2
    FP_TYPE_64 = 4
end

@cenum FlushMode::UInt32 begin
    FLUSH_MODE_FTZ = 1
    FLUSH_MODE_NON_FTZ = 2
end

@cenum RoundMethod::UInt32 begin
    ROUND_METHOD_SINGLE = 1
    ROUND_METHOD_DOUBLE = 2
end

@cenum WavefrontInfo::UInt32 begin
    WAVEFRONT_INFO_SIZE = 0
end

@cenum ExecutableState::UInt32 begin
    EXECUTABLE_STATE_UNFROZEN = 0
    EXECUTABLE_STATE_FROZEN = 1
end

@cenum ExecutableInfo::UInt32 begin
    EXECUTABLE_INFO_PROFILE = 1
    EXECUTABLE_INFO_STATE = 2
    EXECUTABLE_INFO_DEFAULT_FLOAT_ROUNDING_MODE = 3
end

@cenum SymbolKind::UInt32 begin
    SYMBOL_KIND_VARIABLE = 0
    SYMBOL_KIND_KERNEL = 1
    SYMBOL_KIND_INDIRECT_FUNCTION = 2
end

@cenum SymbolLinkage::UInt32 begin
    SYMBOL_LINKAGE_MODULE = 0
    SYMBOL_LINKAGE_PROGRAM = 1
end

@cenum VariableAllocation::UInt32 begin
    VARIABLE_ALLOCATION_AGENT = 0
    VARIABLE_ALLOCATION_PROGRAM = 1
end

@cenum VariableSegment::UInt32 begin
    VARIABLE_SEGMENT_GLOBAL = 0
    VARIABLE_SEGMENT_READONLY = 1
end

@cenum ExecutableSymbolInfo::UInt32 begin
    EXECUTABLE_SYMBOL_INFO_TYPE = 0
    EXECUTABLE_SYMBOL_INFO_NAME_LENGTH = 1
    EXECUTABLE_SYMBOL_INFO_NAME = 2
    EXECUTABLE_SYMBOL_INFO_MODULE_NAME_LENGTH = 3
    EXECUTABLE_SYMBOL_INFO_MODULE_NAME = 4
    EXECUTABLE_SYMBOL_INFO_AGENT = 20
    EXECUTABLE_SYMBOL_INFO_VARIABLE_ADDRESS = 21
    EXECUTABLE_SYMBOL_INFO_LINKAGE = 5
    EXECUTABLE_SYMBOL_INFO_IS_DEFINITION = 17
    EXECUTABLE_SYMBOL_INFO_VARIABLE_ALLOCATION = 6
    EXECUTABLE_SYMBOL_INFO_VARIABLE_SEGMENT = 7
    EXECUTABLE_SYMBOL_INFO_VARIABLE_ALIGNMENT = 8
    EXECUTABLE_SYMBOL_INFO_VARIABLE_SIZE = 9
    EXECUTABLE_SYMBOL_INFO_VARIABLE_IS_CONST = 10
    EXECUTABLE_SYMBOL_INFO_KERNEL_OBJECT = 22
    EXECUTABLE_SYMBOL_INFO_KERNEL_KERNARG_SEGMENT_SIZE = 11
    EXECUTABLE_SYMBOL_INFO_KERNEL_KERNARG_SEGMENT_ALIGNMENT = 12
    EXECUTABLE_SYMBOL_INFO_KERNEL_GROUP_SEGMENT_SIZE = 13
    EXECUTABLE_SYMBOL_INFO_KERNEL_PRIVATE_SEGMENT_SIZE = 14
    EXECUTABLE_SYMBOL_INFO_KERNEL_DYNAMIC_CALLSTACK = 15
    EXECUTABLE_SYMBOL_INFO_KERNEL_CALL_CONVENTION = 18
    EXECUTABLE_SYMBOL_INFO_INDIRECT_FUNCTION_OBJECT = 23
    EXECUTABLE_SYMBOL_INFO_INDIRECT_FUNCTION_CALL_CONVENTION = 16
end

@cenum CodeObjectType::UInt32 begin
    CODE_OBJECT_TYPE_PROGRAM = 0
end

@cenum CodeObjectInfo::UInt32 begin
    CODE_OBJECT_INFO_VERSION = 0
    CODE_OBJECT_INFO_TYPE = 1
    CODE_OBJECT_INFO_ISA = 2
    CODE_OBJECT_INFO_MACHINE_MODEL = 3
    CODE_OBJECT_INFO_PROFILE = 4
    CODE_OBJECT_INFO_DEFAULT_FLOAT_ROUNDING_MODE = 5
end

@cenum CodeSymbolInfo::UInt32 begin
    CODE_SYMBOL_INFO_TYPE = 0
    CODE_SYMBOL_INFO_NAME_LENGTH = 1
    CODE_SYMBOL_INFO_NAME = 2
    CODE_SYMBOL_INFO_MODULE_NAME_LENGTH = 3
    CODE_SYMBOL_INFO_MODULE_NAME = 4
    CODE_SYMBOL_INFO_LINKAGE = 5
    CODE_SYMBOL_INFO_IS_DEFINITION = 17
    CODE_SYMBOL_INFO_VARIABLE_ALLOCATION = 6
    CODE_SYMBOL_INFO_VARIABLE_SEGMENT = 7
    CODE_SYMBOL_INFO_VARIABLE_ALIGNMENT = 8
    CODE_SYMBOL_INFO_VARIABLE_SIZE = 9
    CODE_SYMBOL_INFO_VARIABLE_IS_CONST = 10
    CODE_SYMBOL_INFO_KERNEL_KERNARG_SEGMENT_SIZE = 11
    CODE_SYMBOL_INFO_KERNEL_KERNARG_SEGMENT_ALIGNMENT = 12
    CODE_SYMBOL_INFO_KERNEL_GROUP_SEGMENT_SIZE = 13
    CODE_SYMBOL_INFO_KERNEL_PRIVATE_SEGMENT_SIZE = 14
    CODE_SYMBOL_INFO_KERNEL_DYNAMIC_CALLSTACK = 15
    CODE_SYMBOL_INFO_KERNEL_CALL_CONVENTION = 18
    CODE_SYMBOL_INFO_INDIRECT_FUNCTION_CALL_CONVENTION = 16
end

@cenum KernelCodeVersion::UInt32 begin
    KERNEL_CODE_VERSION_MAJOR = 1
    KERNEL_CODE_VERSION_MINOR = 1
end

@cenum MachineKind::UInt16 begin
    MACHINE_KIND_UNDEFINED = 0
    MACHINE_KIND_AMDGPU = 1
end

@cenum FloatRoundMode::UInt16 begin
    FLOAT_ROUND_MODE_NEAREST_EVEN = 0
    FLOAT_ROUND_MODE_PLUS_INFINITY = 1
    FLOAT_ROUND_MODE_MINUS_INFINITY = 2
    FLOAT_ROUND_MODE_ZERO = 3
end

@cenum FloatDenormMode::UInt32 begin
    FLOAT_DENORM_MODE_FLUSH_SOURCE_OUTPUT = 0
    FLOAT_DENORM_MODE_FLUSH_OUTPUT = 1
    FLOAT_DENORM_MODE_FLUSH_SOURCE = 2
    FLOAT_DENORM_MODE_NO_FLUSH = 3
end

@cenum ComputePgmRsrcOne::UInt32 begin
    COMPUTE_PGM_RSRC_ONE_GRANULATED_WORKITEM_VGPR_COUNT_SHIFT = 0
    COMPUTE_PGM_RSRC_ONE_GRANULATED_WORKITEM_VGPR_COUNT_WIDTH = 6
    COMPUTE_PGM_RSRC_ONE_GRANULATED_WORKITEM_VGPR_COUNT = 63
    COMPUTE_PGM_RSRC_ONE_GRANULATED_WAVEFRONT_SGPR_COUNT_SHIFT = 6
    COMPUTE_PGM_RSRC_ONE_GRANULATED_WAVEFRONT_SGPR_COUNT_WIDTH = 4
    COMPUTE_PGM_RSRC_ONE_GRANULATED_WAVEFRONT_SGPR_COUNT = 960
    COMPUTE_PGM_RSRC_ONE_PRIORITY_SHIFT = 10
    COMPUTE_PGM_RSRC_ONE_PRIORITY_WIDTH = 2
    COMPUTE_PGM_RSRC_ONE_PRIORITY = 3072
    COMPUTE_PGM_RSRC_ONE_FLOAT_ROUND_MODE_32_SHIFT = 12
    COMPUTE_PGM_RSRC_ONE_FLOAT_ROUND_MODE_32_WIDTH = 2
    COMPUTE_PGM_RSRC_ONE_FLOAT_ROUND_MODE_32 = 12288
    COMPUTE_PGM_RSRC_ONE_FLOAT_ROUND_MODE_16_64_SHIFT = 14
    COMPUTE_PGM_RSRC_ONE_FLOAT_ROUND_MODE_16_64_WIDTH = 2
    COMPUTE_PGM_RSRC_ONE_FLOAT_ROUND_MODE_16_64 = 49152
    COMPUTE_PGM_RSRC_ONE_FLOAT_DENORM_MODE_32_SHIFT = 16
    COMPUTE_PGM_RSRC_ONE_FLOAT_DENORM_MODE_32_WIDTH = 2
    COMPUTE_PGM_RSRC_ONE_FLOAT_DENORM_MODE_32 = 196608
    COMPUTE_PGM_RSRC_ONE_FLOAT_DENORM_MODE_16_64_SHIFT = 18
    COMPUTE_PGM_RSRC_ONE_FLOAT_DENORM_MODE_16_64_WIDTH = 2
    COMPUTE_PGM_RSRC_ONE_FLOAT_DENORM_MODE_16_64 = 786432
    COMPUTE_PGM_RSRC_ONE_PRIV_SHIFT = 20
    COMPUTE_PGM_RSRC_ONE_PRIV_WIDTH = 1
    COMPUTE_PGM_RSRC_ONE_PRIV = 1048576
    COMPUTE_PGM_RSRC_ONE_ENABLE_DX10_CLAMP_SHIFT = 21
    COMPUTE_PGM_RSRC_ONE_ENABLE_DX10_CLAMP_WIDTH = 1
    COMPUTE_PGM_RSRC_ONE_ENABLE_DX10_CLAMP = 2097152
    COMPUTE_PGM_RSRC_ONE_DEBUG_MODE_SHIFT = 22
    COMPUTE_PGM_RSRC_ONE_DEBUG_MODE_WIDTH = 1
    COMPUTE_PGM_RSRC_ONE_DEBUG_MODE = 4194304
    COMPUTE_PGM_RSRC_ONE_ENABLE_IEEE_MODE_SHIFT = 23
    COMPUTE_PGM_RSRC_ONE_ENABLE_IEEE_MODE_WIDTH = 1
    COMPUTE_PGM_RSRC_ONE_ENABLE_IEEE_MODE = 8388608
    COMPUTE_PGM_RSRC_ONE_BULKY_SHIFT = 24
    COMPUTE_PGM_RSRC_ONE_BULKY_WIDTH = 1
    COMPUTE_PGM_RSRC_ONE_BULKY = 16777216
    COMPUTE_PGM_RSRC_ONE_CDBG_USER_SHIFT = 25
    COMPUTE_PGM_RSRC_ONE_CDBG_USER_WIDTH = 1
    COMPUTE_PGM_RSRC_ONE_CDBG_USER = 33554432
    COMPUTE_PGM_RSRC_ONE_RESERVED1_SHIFT = 26
    COMPUTE_PGM_RSRC_ONE_RESERVED1_WIDTH = 6
    #AMD_COMPUTE_PGM_RSRC_ONE_RESERVED1 = -67108864
end

@cenum SystemVGPRWorkitemID::UInt32 begin
    SYSTEM_VGPR_WORKITEM_ID_X = 0
    SYSTEM_VGPR_WORKITEM_ID_X_Y = 1
    SYSTEM_VGPR_WORKITEM_ID_X_Y_Z = 2
    SYSTEM_VGPR_WORKITEM_ID_UNDEFINED = 3
end

@cenum ComputePgmRsrcTwo::UInt32 begin
    COMPUTE_PGM_RSRC_TWO_ENABLE_SGPR_PRIVATE_SEGMENT_WAVE_BYTE_OFFSET_SHIFT = 0
    COMPUTE_PGM_RSRC_TWO_ENABLE_SGPR_PRIVATE_SEGMENT_WAVE_BYTE_OFFSET_WIDTH = 1
    COMPUTE_PGM_RSRC_TWO_ENABLE_SGPR_PRIVATE_SEGMENT_WAVE_BYTE_OFFSET = 1
    COMPUTE_PGM_RSRC_TWO_USER_SGPR_COUNT_SHIFT = 1
    COMPUTE_PGM_RSRC_TWO_USER_SGPR_COUNT_WIDTH = 5
    COMPUTE_PGM_RSRC_TWO_USER_SGPR_COUNT = 62
    COMPUTE_PGM_RSRC_TWO_ENABLE_TRAP_HANDLER_SHIFT = 6
    COMPUTE_PGM_RSRC_TWO_ENABLE_TRAP_HANDLER_WIDTH = 1
    COMPUTE_PGM_RSRC_TWO_ENABLE_TRAP_HANDLER = 64
    COMPUTE_PGM_RSRC_TWO_ENABLE_SGPR_WORKGROUP_ID_X_SHIFT = 7
    COMPUTE_PGM_RSRC_TWO_ENABLE_SGPR_WORKGROUP_ID_X_WIDTH = 1
    COMPUTE_PGM_RSRC_TWO_ENABLE_SGPR_WORKGROUP_ID_X = 128
    COMPUTE_PGM_RSRC_TWO_ENABLE_SGPR_WORKGROUP_ID_Y_SHIFT = 8
    COMPUTE_PGM_RSRC_TWO_ENABLE_SGPR_WORKGROUP_ID_Y_WIDTH = 1
    COMPUTE_PGM_RSRC_TWO_ENABLE_SGPR_WORKGROUP_ID_Y = 256
    COMPUTE_PGM_RSRC_TWO_ENABLE_SGPR_WORKGROUP_ID_Z_SHIFT = 9
    COMPUTE_PGM_RSRC_TWO_ENABLE_SGPR_WORKGROUP_ID_Z_WIDTH = 1
    COMPUTE_PGM_RSRC_TWO_ENABLE_SGPR_WORKGROUP_ID_Z = 512
    COMPUTE_PGM_RSRC_TWO_ENABLE_SGPR_WORKGROUP_INFO_SHIFT = 10
    COMPUTE_PGM_RSRC_TWO_ENABLE_SGPR_WORKGROUP_INFO_WIDTH = 1
    COMPUTE_PGM_RSRC_TWO_ENABLE_SGPR_WORKGROUP_INFO = 1024
    COMPUTE_PGM_RSRC_TWO_ENABLE_VGPR_WORKITEM_ID_SHIFT = 11
    COMPUTE_PGM_RSRC_TWO_ENABLE_VGPR_WORKITEM_ID_WIDTH = 2
    COMPUTE_PGM_RSRC_TWO_ENABLE_VGPR_WORKITEM_ID = 6144
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_ADDRESS_WATCH_SHIFT = 13
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_ADDRESS_WATCH_WIDTH = 1
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_ADDRESS_WATCH = 8192
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_MEMORY_VIOLATION_SHIFT = 14
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_MEMORY_VIOLATION_WIDTH = 1
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_MEMORY_VIOLATION = 16384
    COMPUTE_PGM_RSRC_TWO_GRANULATED_LDS_SIZE_SHIFT = 15
    COMPUTE_PGM_RSRC_TWO_GRANULATED_LDS_SIZE_WIDTH = 9
    COMPUTE_PGM_RSRC_TWO_GRANULATED_LDS_SIZE = 16744448
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_IEEE_754_FP_INVALID_OPERATION_SHIFT = 24
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_IEEE_754_FP_INVALID_OPERATION_WIDTH = 1
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_IEEE_754_FP_INVALID_OPERATION = 16777216
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_FP_DENORMAL_SOURCE_SHIFT = 25
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_FP_DENORMAL_SOURCE_WIDTH = 1
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_FP_DENORMAL_SOURCE = 33554432
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_IEEE_754_FP_DIVISION_BY_ZERO_SHIFT = 26
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_IEEE_754_FP_DIVISION_BY_ZERO_WIDTH = 1
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_IEEE_754_FP_DIVISION_BY_ZERO = 67108864
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_IEEE_754_FP_OVERFLOW_SHIFT = 27
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_IEEE_754_FP_OVERFLOW_WIDTH = 1
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_IEEE_754_FP_OVERFLOW = 134217728
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_IEEE_754_FP_UNDERFLOW_SHIFT = 28
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_IEEE_754_FP_UNDERFLOW_WIDTH = 1
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_IEEE_754_FP_UNDERFLOW = 268435456
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_IEEE_754_FP_INEXACT_SHIFT = 29
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_IEEE_754_FP_INEXACT_WIDTH = 1
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_IEEE_754_FP_INEXACT = 536870912
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_INT_DIVISION_BY_ZERO_SHIFT = 30
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_INT_DIVISION_BY_ZERO_WIDTH = 1
    COMPUTE_PGM_RSRC_TWO_ENABLE_EXCEPTION_INT_DIVISION_BY_ZERO = 1073741824
    COMPUTE_PGM_RSRC_TWO_RESERVED1_SHIFT = 31
    COMPUTE_PGM_RSRC_TWO_RESERVED1_WIDTH = 1
    #AMD_COMPUTE_PGM_RSRC_TWO_RESERVED1 = -2147483648
end

@cenum ElementByteSize::UInt32 begin
    ELEMENT_BYTE_SIZE_2 = 0
    ELEMENT_BYTE_SIZE_4 = 1
    ELEMENT_BYTE_SIZE_8 = 2
    ELEMENT_BYTE_SIZE_16 = 3
end

@cenum KernelCodeProperties::UInt32 begin
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_PRIVATE_SEGMENT_BUFFER_SHIFT = 0
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_PRIVATE_SEGMENT_BUFFER_WIDTH = 1
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_PRIVATE_SEGMENT_BUFFER = 1
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_DISPATCH_PTR_SHIFT = 1
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_DISPATCH_PTR_WIDTH = 1
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_DISPATCH_PTR = 2
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_QUEUE_PTR_SHIFT = 2
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_QUEUE_PTR_WIDTH = 1
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_QUEUE_PTR = 4
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_KERNARG_SEGMENT_PTR_SHIFT = 3
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_KERNARG_SEGMENT_PTR_WIDTH = 1
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_KERNARG_SEGMENT_PTR = 8
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_DISPATCH_ID_SHIFT = 4
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_DISPATCH_ID_WIDTH = 1
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_DISPATCH_ID = 16
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_FLAT_SCRATCH_INIT_SHIFT = 5
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_FLAT_SCRATCH_INIT_WIDTH = 1
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_FLAT_SCRATCH_INIT = 32
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_PRIVATE_SEGMENT_SIZE_SHIFT = 6
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_PRIVATE_SEGMENT_SIZE_WIDTH = 1
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_PRIVATE_SEGMENT_SIZE = 64
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_GRID_WORKGROUP_COUNT_X_SHIFT = 7
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_GRID_WORKGROUP_COUNT_X_WIDTH = 1
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_GRID_WORKGROUP_COUNT_X = 128
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_GRID_WORKGROUP_COUNT_Y_SHIFT = 8
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_GRID_WORKGROUP_COUNT_Y_WIDTH = 1
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_GRID_WORKGROUP_COUNT_Y = 256
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_GRID_WORKGROUP_COUNT_Z_SHIFT = 9
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_GRID_WORKGROUP_COUNT_Z_WIDTH = 1
    KERNEL_CODE_PROPERTIES_ENABLE_SGPR_GRID_WORKGROUP_COUNT_Z = 512
    KERNEL_CODE_PROPERTIES_RESERVED1_SHIFT = 10
    KERNEL_CODE_PROPERTIES_RESERVED1_WIDTH = 6
    KERNEL_CODE_PROPERTIES_RESERVED1 = 64512
    KERNEL_CODE_PROPERTIES_ENABLE_ORDERED_APPEND_GDS_SHIFT = 16
    KERNEL_CODE_PROPERTIES_ENABLE_ORDERED_APPEND_GDS_WIDTH = 1
    KERNEL_CODE_PROPERTIES_ENABLE_ORDERED_APPEND_GDS = 65536
    KERNEL_CODE_PROPERTIES_PRIVATE_ELEMENT_SIZE_SHIFT = 17
    KERNEL_CODE_PROPERTIES_PRIVATE_ELEMENT_SIZE_WIDTH = 2
    KERNEL_CODE_PROPERTIES_PRIVATE_ELEMENT_SIZE = 393216
    KERNEL_CODE_PROPERTIES_IS_PTR64_SHIFT = 19
    KERNEL_CODE_PROPERTIES_IS_PTR64_WIDTH = 1
    KERNEL_CODE_PROPERTIES_IS_PTR64 = 524288
    KERNEL_CODE_PROPERTIES_IS_DYNAMIC_CALLSTACK_SHIFT = 20
    KERNEL_CODE_PROPERTIES_IS_DYNAMIC_CALLSTACK_WIDTH = 1
    KERNEL_CODE_PROPERTIES_IS_DYNAMIC_CALLSTACK = 1048576
    KERNEL_CODE_PROPERTIES_IS_DEBUG_ENABLED_SHIFT = 21
    KERNEL_CODE_PROPERTIES_IS_DEBUG_ENABLED_WIDTH = 1
    KERNEL_CODE_PROPERTIES_IS_DEBUG_ENABLED = 2097152
    KERNEL_CODE_PROPERTIES_IS_XNACK_ENABLED_SHIFT = 22
    KERNEL_CODE_PROPERTIES_IS_XNACK_ENABLED_WIDTH = 1
    KERNEL_CODE_PROPERTIES_IS_XNACK_ENABLED = 4194304
    KERNEL_CODE_PROPERTIES_RESERVED2_SHIFT = 23
    KERNEL_CODE_PROPERTIES_RESERVED2_WIDTH = 9
    #KERNEL_CODE_PROPERTIES_RESERVED2 = -8388608
end

@cenum PowerTwo::UInt8 begin
    POWERTWO_1 = 0
    POWERTWO_2 = 1
    POWERTWO_4 = 2
    POWERTWO_8 = 3
    POWERTWO_16 = 4
    POWERTWO_32 = 5
    POWERTWO_64 = 6
    POWERTWO_128 = 7
    POWERTWO_256 = 8
end

@cenum EnabledControlDirective::UInt64 begin
    ENABLE_BREAK_EXCEPTIONS = 1
    ENABLE_DETECT_EXCEPTIONS = 2
    MAX_DYNAMIC_GROUP_SIZE = 4
    MAX_FLAT_GRID_SIZE = 8
    MAX_FLAT_WORKGROUP_SIZE = 16
    REQUIRED_DIM = 32
    REQUIRED_GRID_SIZE = 64
    REQUIRED_WORKGROUP_SIZE = 128
    REQUIRE_NO_PARTIAL_WORKGROUPS = 256
end

@cenum ExceptionKind::UInt16 begin
    EXCEPTION_KIND_INVALID_OPERATION = 1
    EXCEPTION_KIND_DIVISION_BY_ZERO = 2
    EXCEPTION_KIND_OVERFLOW = 4
    EXCEPTION_KIND_UNDERFLOW = 8
    EXCEPTION_KIND_INEXACT = 16
end

@cenum QueueProperties::UInt32 begin
    QUEUE_PROPERTIES_ENABLE_TRAP_HANDLER_SHIFT = 0
    QUEUE_PROPERTIES_ENABLE_TRAP_HANDLER_WIDTH = 1
    QUEUE_PROPERTIES_ENABLE_TRAP_HANDLER = 1
    QUEUE_PROPERTIES_IS_PTR64_SHIFT = 1
    QUEUE_PROPERTIES_IS_PTR64_WIDTH = 1
    QUEUE_PROPERTIES_IS_PTR64 = 2
    QUEUE_PROPERTIES_ENABLE_TRAP_HANDLER_DEBUG_SGPRS_SHIFT = 2
    QUEUE_PROPERTIES_ENABLE_TRAP_HANDLER_DEBUG_SGPRS_WIDTH = 1
    QUEUE_PROPERTIES_ENABLE_TRAP_HANDLER_DEBUG_SGPRS = 4
    QUEUE_PROPERTIES_ENABLE_PROFILING_SHIFT = 3
    QUEUE_PROPERTIES_ENABLE_PROFILING_WIDTH = 1
    QUEUE_PROPERTIES_ENABLE_PROFILING = 8
    QUEUE_PROPERTIES_USE_SCRATCH_ONCE_SHIFT = 4
    QUEUE_PROPERTIES_USE_SCRATCH_ONCE_WIDTH = 1
    QUEUE_PROPERTIES_USE_SCRATCH_ONCE = 16
    QUEUE_PROPERTIES__RESERVED1_SHIFT = 5
    QUEUE_PROPERTIES__RESERVED1_WIDTH = 27
    #_RESERVED1 = -32
end

@cenum SignalKind::Int64 begin
    SIGNAL_KIND_INVALID = 0
    SIGNAL_KIND_USER = 1
    SIGNAL_KIND_DOORBELL = -1
    SIGNAL_KIND_LEGACY_DOORBELL = -2
end

@cenum ImageChannelType::UInt8 begin
   IMAGE_CHANNEL_TYPE_SNORM_INT8 = 0
   IMAGE_CHANNEL_TYPE_SNORM_INT16 = 1
   IMAGE_CHANNEL_TYPE_UNORM_INT8 = 2
   IMAGE_CHANNEL_TYPE_UNORM_INT16 = 3
   IMAGE_CHANNEL_TYPE_UNORM_INT24 = 4
   IMAGE_CHANNEL_TYPE_SHORT_555 = 5
   IMAGE_CHANNEL_TYPE_SHORT_565 = 6
   IMAGE_CHANNEL_TYPE_INT_101010 = 7
   IMAGE_CHANNEL_TYPE_SIGNED_INT8 = 8
   IMAGE_CHANNEL_TYPE_SIGNED_INT16 = 9
   IMAGE_CHANNEL_TYPE_SIGNED_INT32 = 10
   IMAGE_CHANNEL_TYPE_UNSIGNED_INT8 = 11
   IMAGE_CHANNEL_TYPE_UNSIGNED_INT16 = 12
   IMAGE_CHANNEL_TYPE_UNSIGNED_INT32 = 13
   IMAGE_CHANNEL_TYPE_HALF_FLOAT = 14
   IMAGE_CHANNEL_TYPE_FLOAT = 15
end

@cenum ImageCapability::UInt32 begin
    IMAGE_CAPABILITY_NOT_SUPPORTED = 0
    IMAGE_CAPABILITY_READ_ONLY = 1
    IMAGE_CAPABILITY_WRITE_ONLY = 2
    IMAGE_CAPABILITY_READ_WRITE = 4
    IMAGE_CAPABILITY_READ_MODIFY_WRITE = 8
    IMAGE_CAPABILITY_ACCESS_INVARIANT_DATA_LAYOUT = 16
end

@cenum ImageDataLayout::UInt32 begin
    IMAGE_DATA_LAYOUT_OPAQUE = 0
    IMAGE_DATA_LAYOUT_LINEAR = 1
end

@cenum SamplerAddressingMode::UInt32 begin
    SAMPLER_ADDRESSING_MODE_UNDEFINED = 0
    SAMPLER_ADDRESSING_MODE_CLAMP_TO_EDGE = 1
    SAMPLER_ADDRESSING_MODE_CLAMP_TO_BORDER = 2
    SAMPLER_ADDRESSING_MODE_REPEAT = 3
    SAMPLER_ADDRESSING_MODE_MIRRORED_REPEAT = 4
end

@cenum SamplerCoordinateMode::UInt32 begin
    SAMPLER_COORDINATE_MODE_UNNORMALIZED = 0
    SAMPLER_COORDINATE_MODE_NORMALIZED = 1
end

@cenum SamplerFilterMode::UInt32 begin
    SAMPLER_FILTER_MODE_NEAREST = 0
    SAMPLER_FILTER_MODE_LINEAR = 1
end

@cenum AMDAgentInfo::UInt32 begin
    AMD_AGENT_INFO_CHIP_ID = 40960
    AMD_AGENT_INFO_CACHELINE_SIZE = 40961
    AMD_AGENT_INFO_COMPUTE_UNIT_COUNT = 40962
    AMD_AGENT_INFO_MAX_CLOCK_FREQUENCY = 40963
    AMD_AGENT_INFO_DRIVER_NODE_ID = 40964
    AMD_AGENT_INFO_MAX_ADDRESS_WATCH_POINTS = 40965
    AMD_AGENT_INFO_BDFID = 40966
    AMD_AGENT_INFO_MEMORY_WIDTH = 40967
    AMD_AGENT_INFO_MEMORY_MAX_FREQUENCY = 40968
    AMD_AGENT_INFO_PRODUCT_NAME = 40969
    AMD_AGENT_INFO_MAX_WAVES_PER_CU = 40970
    AMD_AGENT_INFO_NUM_SIMDS_PER_CU = 40971
    AMD_AGENT_INFO_NUM_SHADER_ENGINES = 40972
    AMD_AGENT_INFO_NUM_SHADER_ARRAYS_PER_SE = 40973
    AMD_AGENT_INFO_HDP_FLUSH = 40974
    AMD_AGENT_INFO_DOMAIN = 40975
end

@cenum AMDRegionInfo::UInt32 begin
    AMD_REGION_INFO_HOST_ACCESSIBLE = 40960
    AMD_REGION_INFO_BASE = 40961
    AMD_REGION_INFO_BUS_WIDTH = 40962
    AMD_REGION_INFO_MAX_CLOCK_FREQUENCY = 40963
end

@cenum CoherencyType::UInt32 begin
    COHERENCY_TYPE_COHERENT = 0
    COHERENCY_TYPE_NONCOHERENT = 1
end

@cenum AMDSignalAttribute::UInt32 begin
    AMD_SIGNAL_AMD_GPU_ONLY = 1
    AMD_SIGNAL_IPC = 2
end

@cenum Segment::UInt32 begin
    SEGMENT_GLOBAL = 0
    SEGMENT_READONLY = 1
    SEGMENT_PRIVATE = 2
    SEGMENT_GROUP = 3
end

@cenum MemoryPoolGlobalFlag::UInt32 begin
    MEMORY_POOL_GLOBAL_FLAG_KERNARG_INIT = 1
    MEMORY_POOL_GLOBAL_FLAG_FINE_GRAINED = 2
    MEMORY_POOL_GLOBAL_FLAG_COARSE_GRAINED = 4
end

@cenum MemoryPoolInfo::UInt32 begin
    MEMORY_POOL_INFO_SEGMENT = 0
    MEMORY_POOL_INFO_GLOBAL_FLAGS = 1
    MEMORY_POOL_INFO_SIZE = 2
    MEMORY_POOL_INFO_RUNTIME_ALLOC_ALLOWED = 5
    MEMORY_POOL_INFO_RUNTIME_ALLOC_GRANULE = 6
    MEMORY_POOL_INFO_RUNTIME_ALLOC_ALIGNMENT = 7
    MEMORY_POOL_INFO_ACCESSIBLE_BY_ALL = 15
    MEMORY_POOL_INFO_ALLOC_MAX_SIZE = 16
end

@cenum CopyDirection::UInt32 begin
    HostToHost = 0
    HostToDevice = 1
    DeviceToHost = 2
    DeviceToDevice = 3
end

@cenum MemoryPoolAccess::UInt32 begin
    MEMORY_POOL_ACCESS_NEVER_ALLOWED = 0
    MEMORY_POOL_ACCESS_ALLOWED_BY_DEFAULT = 1
    MEMORY_POOL_ACCESS_DISALLOWED_BY_DEFAULT = 2
end

@cenum LinkInfoType::UInt32 begin
    LINK_INFO_TYPE_HYPERTRANSPORT = 0
    LINK_INFO_TYPE_QPI = 1
    LINK_INFO_TYPE_PCIE = 2
    LINK_INFO_TYPE_INFINBAND = 3
    LINK_INFO_TYPE_XGMI = 4
end

@cenum AgentMemoryPoolInfo::UInt32 begin
    AGENT_MEMORY_POOL_INFO_ACCESS = 0
    AGENT_MEMORY_POOL_INFO_NUM_LINK_HOPS = 1
    AGENT_MEMORY_POOL_INFO_LINK_INFO = 2
end

@cenum PointerType::UInt32 begin
    POINTER_TYPE_UNKNOWN = 0
    POINTER_TYPE_HSA = 1
    POINTER_TYPE_LOCKED = 2
    POINTER_TYPE_GRAPHICS = 3
    POINTER_TYPE_IPC = 4
end

@cenum AMDEventType::UInt32 begin
    AMD_GPU_MEMORY_FAULT_EVENT = 0
end

@cenum MemoryFaultReason::Int32 begin
    MEMORY_FAULT_PAGE_NOT_PRESENT = 1
    MEMORY_FAULT_READ_ONLY = 2
    MEMORY_FAULT_NX = 4
    MEMORY_FAULT_HOST_ONLY = 8
    MEMORY_FAULT_DRAM_ECC = 16
    MEMORY_FAULT_IMPRECISE = 32
    MEMORY_FAULT_SRAM_ECC = 64
    MEMORY_FAULT_HANG = -2147483648
end

@cenum QueuePriority::UInt32 begin
    QUEUE_PRIORITY_LOW = 0
    QUEUE_PRIORITY_NORMAL = 1
    QUEUE_PRIORITY_HIGH = 2
end

@cenum WaveAction::UInt32 begin
    WAVE_ACTION_HALT = 1
    WAVE_ACTION_RESUME = 2
    WAVE_ACTION_KILL = 3
    WAVE_ACTION_DEBUG = 4
    WAVE_ACTION_TRAP = 5
    WAVE_ACTION_MAX = 6
end

@cenum HostAction::UInt32 begin
    HOST_ACTION_IGNORE = 1
    HOST_ACTION_EXIT = 2
    HOST_ACTION_NOTIFY = 4
end

@cenum WaveMode::UInt32 begin
    WAVE_MODE_SINGLE = 0
    WAVE_MODE_BROADCAST_PROCESS = 2
    WAVE_MODE_BROADCAST_PROCESS_CU = 3
    WAVE_MODE_MAX = 4
end

@cenum TrapType::UInt32 begin
    DEBUG_TRAP = 1
    MAX_TRAP = 2
end

@cenum AddressWatchMode::UInt32 begin
    ADDRESS_WATCH_READ = 0
    ADDRESS_WATCH_NON_READ = 1
    ADDRESS_WATCH_ATOMIC = 2
    ADDRESS_WATCH_ALL = 3
    ADDRESS_WATCH_NUM = 4
end

@cenum EventWaitStatus::UInt32 begin
    EVENT_WAIT_SUCCESS = 0
    EVENT_WAIT_ERROR = 4096
    EVENT_WAIT_TIMEOUT = 81923
end

@cenum ProgramInfo::UInt32 begin
    PROGRAM_INFO_MACHINE_MODEL = 0
    PROGRAM_INFO_PROFILE = 1
    PROGRAM_INFO_DEFAULT_FLOAT_ROUNDING_MODE = 2
end

@cenum FinalizerCallConvention::Int32 begin
    FINALIZER_CALL_CONVENTION_AUTO = -1
end

@cenum AQLProfileEventType::UInt32 begin
    AQLPROFILE_EVENT_TYPE_PMC = 0
    AQLPROFILE_EVENT_TYPE_SQTT = 1
end

@cenum AQLProfileBlockName::UInt32 begin
    AQLPROFILE_BLOCK_NAME_CPC = 0
    AQLPROFILE_BLOCK_NAME_CPF = 1
    AQLPROFILE_BLOCK_NAME_GDS = 2
    AQLPROFILE_BLOCK_NAME_GRBM = 3
    AQLPROFILE_BLOCK_NAME_GRBMSE = 4
    AQLPROFILE_BLOCK_NAME_SPI = 5
    AQLPROFILE_BLOCK_NAME_SQ = 6
    AQLPROFILE_BLOCK_NAME_SQCS = 7
    AQLPROFILE_BLOCK_NAME_SRBM = 8
    AQLPROFILE_BLOCK_NAME_SX = 9
    AQLPROFILE_BLOCK_NAME_TA = 10
    AQLPROFILE_BLOCK_NAME_TCA = 11
    AQLPROFILE_BLOCK_NAME_TCC = 12
    AQLPROFILE_BLOCK_NAME_TCP = 13
    AQLPROFILE_BLOCK_NAME_TD = 14
    AQLPROFILE_BLOCK_NAME_MCARB = 15
    AQLPROFILE_BLOCK_NAME_MCHUB = 16
    AQLPROFILE_BLOCK_NAME_MCMCBVM = 17
    AQLPROFILE_BLOCK_NAME_MCSEQ = 18
    AQLPROFILE_BLOCK_NAME_MCVML2 = 19
    AQLPROFILE_BLOCK_NAME_MCXBAR = 20
    AQLPROFILE_BLOCK_NAME_ATC = 21
    AQLPROFILE_BLOCK_NAME_ATCL2 = 22
    AQLPROFILE_BLOCK_NAME_GCEA = 23
    AQLPROFILE_BLOCK_NAME_RPB = 24
    AQLPROFILE_BLOCK_NAME_SDMA = 25
    AQLPROFILE_BLOCKS_NUMBER = 26
end

@cenum AQLProfileParameterName::UInt32 begin
    AQLPROFILE_PARAMETER_NAME_COMPUTE_UNIT_TARGET = 0
    AQLPROFILE_PARAMETER_NAME_VM_ID_MASK = 1
    AQLPROFILE_PARAMETER_NAME_MASK = 2
    AQLPROFILE_PARAMETER_NAME_TOKEN_MASK = 3
    AQLPROFILE_PARAMETER_NAME_TOKEN_MASK2 = 4
    AQLPROFILE_PARAMETER_NAME_SE_MASK = 5
    AQLPROFILE_PARAMETER_NAME_SAMPLE_RATE = 6
end

@cenum AQLProfileInfoType::UInt32 begin
    AQLPROFILE_INFO_COMMAND_BUFFER_SIZE = 0
    AQLPROFILE_INFO_PMC_DATA_SIZE = 1
    AQLPROFILE_INFO_PMC_DATA = 2
    AQLPROFILE_INFO_SQTT_DATA = 3
    AQLPROFILE_INFO_BLOCK_COUNTERS = 4
    AQLPROFILE_INFO_BLOCK_ID = 5
    AQLPROFILE_INFO_ENABLE_CMD = 6
end

@cenum CodeObjectStorageType::UInt32 begin
    CODE_OBJECT_STORAGE_TYPE_NONE = 0
    CODE_OBJECT_STORAGE_TYPE_FILE = 1
    CODE_OBJECT_STORAGE_TYPE_MEMORY = 2
end

@cenum LoadedCodeObjectKind::UInt32 begin
    LOADED_CODE_OBJECT_KIND_PROGRAM = 1
    LOADED_CODE_OBJECT_KIND_AGENT = 2
end

@cenum LoadedCodeObjectInfo::UInt32 begin
    LOADED_CODE_OBJECT_INFO_EXECUTABLE = 1
    LOADED_CODE_OBJECT_INFO_KIND = 2
    LOADED_CODE_OBJECT_INFO_AGENT = 3
    LOADED_CODE_OBJECT_INFO_CODE_OBJECT_STORAGE_TYPE = 4
    LOADED_CODE_OBJECT_INFO_CODE_OBJECT_STORAGE_MEMORY_BASE = 5
    LOADED_CODE_OBJECT_INFO_CODE_OBJECT_STORAGE_MEMORY_SIZE = 6
    LOADED_CODE_OBJECT_INFO_CODE_OBJECT_STORAGE_FILE = 7
    LOADED_CODE_OBJECT_INFO_LOAD_DELTA = 8
    LOADED_CODE_OBJECT_INFO_LOAD_BASE = 9
    LOADED_CODE_OBJECT_INFO_LOAD_SIZE = 10
end

ext: jl


text:
using JLD2
using Knet
using Distributions

const EPSILON = 1.0f-5

function normit(x)
    d = ndims(x) == 4 ? [1,2,3] : [1,]
    s = prod(size(x)[d])
    mu = sum(x, dims=d) ./ s
    x = x .- mu
    sigma = sqrt.(EPSILON .+ (sum(x .* x, dims=d)) ./ s)
    return x ./ sigma
end

function predict(w,x)
    x = conv4(w[1],x; stride=4) .+ w[2]
    x = normit(x)
    x = max.(0.3f0.*x,x)
    x = conv4(w[3],x; stride=2) .+ w[4]
    x = normit(x)
    x = max.(0.3f0.*x,x)
    x = mat(x)
    x = w[5]*x .+ w[6]
    x = normit(x)
    x = max.(0.3f0.*x,x)
    return abs.(w[7]*x .+ w[8])
end

function initialize_weights(number_of_actions)
    return [
        0.01f0*randn(Float32, 8, 8, 1, 16),
        0.01f0*randn(Float32, 1, 1, 16),
        0.01f0*randn(Float32, 4,4,16,32),
        0.01f0*randn(Float32, 1,1,32),
        0.01f0*randn(Float32, 256,2592),
        0.01f0*randn(Float32, 256,1),
        0.01f0*randn(Float32, number_of_actions,256),
        0.01f0*randn(Float32, number_of_actions,1),
    ]
end

loss(w,x,y) = sum(abs2, (y.-predict(w,x))) ./ length(y)
lossgradient = gradloss(loss)

function train(w, x,y; lr=.0001f0)
    dw, ll = lossgradient(w, x, y)
    for i in 1:length(w)
        w[i] -= lr .* dw[i]
    end
    return w, ll
end

function trainit(N,w,sensors,mqs, lr, batchsize, priorities)
    h = sort(priorities)[floor(Int64, 0.9*length(priorities))]
    _p = [max(h, p) for p in priorities]
    _at = Distributions.AliasTable(_p)
    tll = 0.0
    for j in 1:N
        if j % 100 == 0
            print('.'); flush(stdout)
        end
        x,y, bix = batchit(sensors, mqs, batchsize, _at);
        w,ll = train(w, x, y; lr=lr)
        priorities[bix] .= ll
        tll += ll
    end
    w, tll/N
end

function batchit(sensors, qs, batchsize, alias_table)
    @assert length(qs) > 0
    @assert length(sensors) == length(qs)
    NACT = length(qs[1])
    x = zeros(Float32, size(sensors[1])[1],size(sensors[1])[2],1,batchsize)
    y = zeros(Float32, NACT,batchsize)
    bix = rand(alias_table, batchsize)
    for (bnum, ix) in enumerate(bix)
        y[:,bnum] .= qs[ix]
        x[:,:,1,bnum] .= sensors[ix] ./ 255.0f0
    end
    KnetArray{Float32}(x), KnetArray{Float32}(y), bix
end

function savew(w)
    ww = [Array(wi) for wi in w]
    @save "w.jld2" ww
end

function loadw()
    @load "w.jld2" ww
    [KnetArray(wi) for wi in ww]
end

function network_predict(env::MyEnv, w, sensors)
    x = zeros(Float32, size(sensors)[1],size(sensors)[2],1,1)
    x[:,:,1,1] .= sensors./255.0f0
    v = Knet.Array(predict(w, Knet.KnetArray(x)))
    return reshape(v,length(v))
end

function moving_avg(y,w)
    ay = Float64[]
    for i in 1:length(y)
        push!(ay,mean(y[max(i-w,1):i]))
    end
    ay
end


text:
using Wavelets

wf = waveletfilter(WT.db6)
N = 1024*32;
x0 = rand(N);
L = int(log2(N)-2)
tn = 1000

println("iwt (N=",N,"), ", L, " levels")
f(x0, L, wf) = for i = 1:tn; idwt(x0, wf, L); end
f(x0, L, wf);
@time f(x0, L, wf);


text:
using FreqTables
using StatsBase
using Plots, StatsPlots
using Statistics
using DataFrames
Pkg.build("Conda")
ENV["PYTHON"]=""; Pkg.build("PyCall")
using PyPlot
using RCall
using CSV

foreach(readdir(Base.source_dir() * "\\methods"; join=true)) do f # load required structs
    if f[(end-2):end] == ".jl"
        println(f)
        include(f)
    end
end

function compare_approaches(k_range, N_range, percentile_cutoff_range, rho_range, niter)
    results = []
    ii = 1
    iters = length(k_range) * length(N_range) * length(percentile_cutoff_range) * niter * length(rho_range)
    for k in k_range
        for N in N_range
            for percentile_cutoff in percentile_cutoff_range
                for rho in rho_range
                    for iter in 1:niter
                        @info "iteration: $ii of $iters"
                        t = time()
                        ifp = method_IFP(k,k,N,percentile_cutoff, rho, verbose = false)
                        t_ipf = t - time()
                        t = time()
                        exc = method_exact(k,k,N,percentile_cutoff, rho, verbose = false)
                        t_exc = time() - t
                        t = time()
                        try
                            mcm = method_sampling_mcmc(k,k,N,percentile_cutoff, rho, verbose = false)
                            t_mcm = time() - t
                            push!(results,(k, N, percentile_cutoff, rho, ifp, exc, mcm, t_ipf, t_exc, t_mcm))
                        catch
                        end
                        ii += 1
                    end
                end
            end
        end
    end
    results
end

results = compare_approaches(3:6, 1000, 5:5:30, 0.00:0.05:0.50, 100)
results = results[[!isnothing(w) for w in results]]

kk = [w[1] for w in results]
nn = [w[2] for w in results]
pc = [w[3] for w in results]
rr = [w[4] for w in results]
ipf_mae = [w[5].mae for w in results]
ipf_vc = [w[5].vc for w in results]
ipf_bl = [w[5].blanks for w in results]
exa_mae = [w[6].mae for w in results]
exa_vc = [w[6].vc for w in results]
exa_bl = [w[6].blanks for w in results]
mcm_mae = [w[7].mae for w in results]
mcm_vc = [w[7].vc for w in results]
mcm_bl = [w[7].blanks for w in results]
ipf_t = [w[8] for w in results]
exa_t = [w[9] for w in results]
mcm_t = [w[10] for w in results]

results = DataFrame(K = kk, N = nn, PC = pc, RR = rr, IPF_blanks = ipf_bl, IPF_time = ipf_t, IPF_VC = ipf_vc, IPF_MAE = ipf_mae, EXA_blanks = exa_bl, EXA_time = exa_t, EXA_VC = exa_vc, EXA_MAE = exa_mae, MCM_blanks = mcm_bl, MCM_time = mcm_t, MCM_VC = mcm_vc, MCM_MAE = mcm_mae)

CSV.write("folder_path/CellSampling_results.csv", results, delim=";", decimal = ',')

max_stars_repo_name: Michal-Kot/FillMissingEntries

max_issues_repo_name: Michal-Kot/FillMissingEntries

max_forks_repo_name: Michal-Kot/FillMissingEntries


text:
function fact(n)
    r = 1
    while n > 1
        r = r * n
        n = n - 1
    end
    return r
end

for n = 0:10
    println(fact(n))
end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
include("../examples/rijke.jl")
include("../src/lss.jl")
using JLD
using OrdinaryDiffEq
function Rijke_tangent_sensitivity(n_spe)
	s = [6.9, 0.2]
	n = 2000
	d = N

	du_trj = zeros(d, d, n)
	dJ_trj = zeros(2,d,n)
	f_trj = zeros(d,n)
	X_trj = zeros(d,n)

	dJds = zeros(2)
	#vsh = zeros(d, n)
	nRunup = 1000000
	u = Rijke_ODE(rand(d),s,nRunup)

	for i = 1:n
		vel_p_i = view(u, 1:2*Ng)
		dJ_trj[1,1:2*Ng,i] = 0.5*vel_p_i
		dJ_trj[2,1:Ng,i] = cjpixf
		X_trj[:,i] = perturbation(u, s, 1.e-6)
		f!(view(f_trj,:,i), u, s, 1.)
		du_trj[:,:, i] = dRijke(u, s, 1.e-6)
		u = Rijke_ODE(u, s, 1)
	end
	y, xi = lss(du_trj, X_trj, zeros(d,n), s, 2)
	dJds = compute_sens(y, xi, dJ_trj, f_trj)
	#vsh[:,:] = y
	save(string("../data/rijke_tangent_sensitivity/",
		    "vsh_and_dJds_", string(n_spe),
		    ".jld"), "dJds", dJds)
end
function collect_sensitivities()
	pmap(Rijke_tangent_sensitivity, 6601:10000)
end

#=
	n_samples = 1000
	dJacds_ens, dJrayds_ens = zeros(2,n_samples),
							zeros(2,n_samples)
	up = zeros(d)
	um = zeros(d)
	eps = 1.e-3
	n = 100 # produces a growth of ≈ e
	for i = 1:n_samples
			@time dJacds_ens[:,i] = Zygote.gradient(s ->Jac(u, s, n), s)[1]
		dJrayds_ens[:,i] = Zygote.gradient(s ->Jray(u, s, n),
										   s)[1]
		u .= Rijke(u, [7.0, 0.2], n)
	end
=#

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
"""
	contact constraints
"""

struct ContactConstraints <: Constraints
	n
	ineq
end

function contact_constraints(model, T)
	n = model.nc * (T + 1) + model.nc * (T - 1) + 2 * model.nb * (T - 1) + 2 * (T - 1)
	ineq = con_ineq_contact(model, T)

	return ContactConstraints(n, ineq)
end

function constraints!(c, Z, con::ContactConstraints, model, idx, h, T)
	shift = 0

	# signed-distance function
	for t = 1:T
		x = view(Z, idx.x[t])

		if t == 1
			q = view(x, 1:model.nq)
			c[1:model.nc] = ϕ_func(model, q)
		end

		q⁺ = view(x, model.nq .+ (1:model.nq))
		c[model.nc + (t-1) * model.nc .+ (1:model.nc)] = ϕ_func(model, q⁺)
	end

	shift += model.nc * (T + 1)

	# friction cone
	for t = 1:T-1
		u = view(Z, idx.u[t])
		c[shift + (t - 1) * model.nc .+ (1:model.nc)] = friction_cone(model, u)
	end

	shift += model.nc * (T - 1)

	# maximum dissipation
	for t = 1:T-1
		x⁺ = view(Z, idx.x[t+1])
		u = view(Z, idx.u[t])
		c[shift + (t - 1) * 2 * model.nb .+ (1:2 * model.nb)] = maximum_dissipation(model,
			x⁺, u, h)
	end

	shift += 2 * model.nb * (T - 1)

	# impact complementarity
	for t = 1:T-1
		x⁺ = view(Z, idx.x[t + 1])
		u = view(Z, idx.u[t])

		q = view(x⁺, model.nq .+ (1:model.nq))
		s = view(u, model.idx_s)
		λ = view(u, model.idx_λ)

		c[shift + (t-1) * 1 + 1] = s[1] - (λ' * ϕ_func(model, q))[1]
	end

	shift += (T - 1)

	# friction cone complementarity
	for t = 1:T-1
		u = view(Z, idx.u[t])

		s = view(u, model.idx_s)
		ψ = view(u, model.idx_ψ)
		λ = view(u, model.idx_λ)
		b = view(u, model.idx_b)

		c[shift + (t-1) * 1 + 1] = s[1] - (ψ' * [model.μ * λ[1]; b])[1]
	end

	shift += (T - 1)

	# # friction parameter complementarity
	# for t = 1:T-1
	# 	u = view(Z, idx.u[t])
	#
	# 	s = view(u, model.idx_s)
	# 	b = view(u, model.idx_b)
	# 	η = view(u, model.idx_η)
	#
	# 	c[shift + (t-1) * 1 + 1] = s[1] - (η' * b)[1]
	# end
	#
	# shift += (T - 1)

    return nothing
end

function constraints_jacobian!(∇c, Z, con::ContactConstraints, model, idx, h, T)
	shift = 0
	c_shift = 0

	# signed-distance function
	ϕ(y) = ϕ_func(model, y)

	for t = 1:T
		x = view(Z, idx.x[t])

		if t == 1
			q = view(x, 1:model.nq)
			r_idx = 1:model.nc
			c_idx = idx.x[t][1:model.nq]
			len = length(r_idx) * length(c_idx)
			∇c[shift .+ (1:len)] = vec(ForwardDiff.jacobian(ϕ, q))
			shift += len
		end

		q⁺ = view(x, model.nq .+ (1:model.nq))
		r_idx = model.nc + (t-1) * model.nc .+ (1:model.nc)
		c_idx = idx.x[t][model.nq .+ (1:model.nq)]
		len = length(r_idx) * length(c_idx)
		∇c[shift .+ (1:len)] = vec(ForwardDiff.jacobian(ϕ, q⁺))
		shift += len
	end

	c_shift += model.nc * (T + 1)

	# friction cone
	fc(y) = friction_cone(model, y)
	for t = 1:T-1
		u = view(Z, idx.u[t])
		r_idx = c_shift + (t - 1) * model.nc .+ (1:model.nc)
		c_idx = idx.u[t]
		len = length(r_idx) * length(c_idx)
		∇c[shift .+ (1:len)] = vec(ForwardDiff.jacobian(fc, u))
		shift += len
	end

	c_shift += model.nc * (T - 1)

	# maximum dissipation
	for t = 1:T-1
		x⁺ = view(Z, idx.x[t+1])
		u = view(Z, idx.u[t])

		mdx(y) = maximum_dissipation(model, y, u, h)
		mdu(y) = maximum_dissipation(model, x⁺, y, h)

		r_idx = c_shift + (t - 1) * 2 * model.nb .+ (1:2 * model.nb)
		c_idx = idx.x[t + 1]
		len = length(r_idx) * length(c_idx)
		∇c[shift .+ (1:len)] = vec(ForwardDiff.jacobian(mdx, x⁺))
		shift += len

		c_idx = idx.u[t]
		len = length(r_idx) * length(c_idx)
		∇c[shift .+ (1:len)] = vec(ForwardDiff.jacobian(mdu, u))
		shift += len
	end

	c_shift += 2 * model.nb * (T - 1)

	# impact complementarity
	for t = 1:T-1
		x⁺ = view(Z, idx.x[t + 1])
		u = view(Z, idx.u[t])

		q = view(x⁺, model.nq .+ (1:model.nq))
		s = view(u, model.idx_s)
		λ = view(u, model.idx_λ)

		cq(y) = s - [λ' * ϕ_func(model, y)]
		cu(y) = y[model.idx_s] - [y[model.idx_λ]' * ϕ_func(model, q)]

		print()
		r_idx = c_shift + (t-1) * 1 + 1
		c_idx = idx.x[t + 1][model.nq .+ (1:model.nq)]
		len = length(r_idx) * length(c_idx)
		∇c[shift .+ (1:len)] = vec(ForwardDiff.jacobian(cq, q))
		shift += len

		c_idx = idx.u[t]
		len = length(r_idx) * length(c_idx)
		∇c[shift .+ (1:len)] = vec(ForwardDiff.jacobian(cu, u))
		shift += len
	end

	c_shift += (T - 1)

	# friction cone complementarity
	for t = 1:T-1
		u = view(Z, idx.u[t])

		s = view(u, model.idx_s)
		ψ = view(u, model.idx_ψ)

		cu(y) = y[model.idx_s] - [y[model.idx_ψ]' * [model.μ * y[model.idx_λ][1]; y[model.idx_b]]]

		r_idx = c_shift + (t-1) * 1 + 1
		c_idx = idx.u[t]
		len = length(r_idx) * length(c_idx)
		∇c[shift .+ (1:len)] = vec(ForwardDiff.jacobian(cu, u))
		shift += len
	end

	c_shift += (T - 1)

	return nothing
end

function constraints_sparsity(con::ContactConstraints, model, idx, T;
	shift_row = 0, shift_col = 0)

	row = []
    col = []

	con_shift = 0

	# signed-distance function
	for t = 1:T

		if t == 1
			r_idx = shift_row .+ (1:model.nc)
			c_idx = shift_col .+ idx.x[t][1:model.nq]
			row_col!(row, col, r_idx, c_idx)
		end

		r_idx = shift_row + model.nc + (t-1) * model.nc .+ (1:model.nc)
		c_idx = shift_col .+ idx.x[t][model.nq .+ (1:model.nq)]
		row_col!(row, col, r_idx, c_idx)
	end

	con_shift += model.nc * (T + 1)

	# friction cone
	for t = 1:T-1
		r_idx = shift_row + con_shift + (t - 1) * model.nc .+ (1:model.nc)
		c_idx = shift_col .+ idx.u[t]
		row_col!(row, col, r_idx, c_idx)
	end

	con_shift += model.nc * (T - 1)

	# maximum dissipation
	for t = 1:T-1
		r_idx = shift_row + con_shift + (t - 1) * 2 * model.nb .+ (1:2 * model.nb)
		c_idx = shift_col .+ idx.x[t + 1]
		row_col!(row, col, r_idx, c_idx)

		c_idx = shift_col .+ idx.u[t]
		row_col!(row, col, r_idx, c_idx)
	end

	con_shift += 2 * model.nb * (T - 1)

	# impact complementarity
	for t = 1:T-1
		r_idx = shift_row + con_shift + (t-1) * 1 + 1
		c_idx = shift_col .+ idx.x[t + 1][model.nq .+ (1:model.nq)]
		row_col!(row, col, r_idx, c_idx)

		c_idx = shift_col .+ idx.u[t]
		row_col!(row, col, r_idx, c_idx)
	end

	con_shift += (T - 1)

	# friction cone complementarity
	for t = 1:T-1
		r_idx = shift_row + con_shift + (t-1) * 1 + 1
		c_idx = shift_col .+ idx.u[t]
		row_col!(row, col, r_idx, c_idx)
	end

	con_shift += (T - 1)

    return collect(zip(row, col))
end

function con_ineq_contact(model, T)

	con_ineq = []
	shift = 0

	# signed-distance function
	push!(con_ineq, [i for i = 1:model.nc * (T + 1)])
	shift += model.nc * (T + 1)

	# friction cone
	push!(con_ineq, [i for i = shift .+ (1:model.nc * (T - 1))])
	shift += model.nc * (T - 1)

	# maximum dissipation
	# push!(con_ineq, [i for i = shift .+ (1:model.nb * (T - 1))])
	shift += 2 * model.nb * (T - 1)

	# impact complementarity
	push!(con_ineq, [i for i = shift .+ (1:(T - 1))])
	shift += (T - 1)

	# friction cone complementarity
	push!(con_ineq, [i for i = shift .+ (1:(T - 1))])
	shift += (T - 1)

    return vcat(con_ineq...)
end


text:
###############################################################################
#
#   padic.jl : flint padic numbers
#
###############################################################################

export FlintPadicField, padic, prime, teichmuller, log, sqrt

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

doc"""
    O(R::FlintPadicField, m::fmpz)
> Construct the value $0 + O(p^n)$ given $m = p^n$. An exception results if $m$
> is not found to be a power of `p = prime(R)`.
"""
function O(R::FlintPadicField, m::fmpz)
   if isone(m)
      N = 0
   else
      p = prime(R)
      if m == p
         N = 1
      else
         N = flog(m, p)
         p^(N) != m && error("Not a power of p in p-adic O()")
      end
   end
   d = padic(N)
   d.parent = R
   return d
end

doc"""
    O(R::FlintPadicField, m::fmpq)
> Construct the value $0 + O(p^n)$ given $m = p^n$. An exception results if $m$
> is not found to be a power of `p = prime(R)`.
"""
function O(R::FlintPadicField, m::fmpq)
   d = denominator(m)
   if isone(d)
      return O(R, numerator(m))
   end
   !isone(numerator(m)) && error("Not a power of p in p-adic O()")
   p = prime(R)
   if d == p
      N = -1
   else
     N = -flog(d, p)
     p^(-N) != d && error("Not a power of p in p-adic O()")
   end
   r = padic(N)
   r.parent = R
   return r
end

doc"""
    O(R::FlintPadicField, m::Integer)
> Construct the value $0 + O(p^n)$ given $m = p^n$. An exception results if $m$
> is not found to be a power of `p = prime(R)`.
"""
O(R::FlintPadicField, m::Integer) = O(R, fmpz(m))

elem_type(::Type{FlintPadicField}) = padic

doc"""
    base_ring(a::FlintPadicField)
> Returns `Union{}` as this field is not dependent on another field.
"""
base_ring(a::FlintPadicField) = Union{}

doc"""
    base_ring(a::padic)
> Returns `Union{}` as this field is not dependent on another field.
"""
base_ring(a::padic) = Union{}

doc"""
    parent(a::padic)
> Returns the parent of the given p-adic field element.
"""
parent(a::padic) = a.parent

isdomain_type(::Type{padic}) = true

isexact_type(R::Type{padic}) = false

function check_parent(a::padic, b::padic)
   parent(a) != parent(b) &&
      error("Incompatible padic rings in padic operation")
end

parent_type(::Type{padic}) = FlintPadicField

###############################################################################
#
#   Basic manipulation
#
###############################################################################

function Base.hash(a::padic, h::UInt)
   return xor(hash(lift(FlintQQ, a), h), xor(hash(prime(parent(a)), h), h))
end

doc"""
    prime(R::FlintPadicField)
> Return the prime $p$ for the given $p$-adic field.
"""
function prime(R::FlintPadicField)
   z = fmpz()
   ccall((:padic_ctx_pow_ui, :libflint), Void,
         (Ref{fmpz}, Int, Ref{FlintPadicField}), z, 1, R)
   return z
end

doc"""
    precision(a::padic)
> Return the precision of the given $p$-adic field element, i.e. if the element
> is known to $O(p^n)$ this function will return $n$.
"""
precision(a::padic) = a.N

doc"""
    valuation(a::padic)
> Return the valuation of the given $p$-adic field element, i.e. if the given
> element is divisible by $p^n$ but not a higher power of $p$ then the function
> will return $n$.
"""
valuation(a::padic) = a.v

doc"""
    lift(R::FlintRationalField, a::padic)
> Return a lift of the given $p$-adic field element to $\mathbb{Q}$.
"""
function lift(R::FlintRationalField, a::padic)
    ctx = parent(a)
    r = fmpq()
    ccall((:padic_get_fmpq, :libflint), Void,
          (Ref{fmpq}, Ref{padic}, Ref{FlintPadicField}), r, a, ctx)
    return r
end

doc"""
    lift(R::FlintIntegerRing, a::padic)
> Return a lift of the given $p$-adic field element to $\mathbb{Z}$.
"""
function lift(R::FlintIntegerRing, a::padic)
    ctx = parent(a)
    r = fmpz()
    ccall((:padic_get_fmpz, :libflint), Void,
          (Ref{fmpz}, Ref{padic}, Ref{FlintPadicField}), r, a, ctx)
    return r
end

doc"""
    zero(R::FlintPadicField)
> Return zero in the given $p$-adic field, to the default precision.
"""
function zero(R::FlintPadicField)
   z = padic(R.prec_max)
   ccall((:padic_zero, :libflint), Void, (Ref{padic},), z)
   z.parent = R
   return z
end

doc"""
    one(R::FlintPadicField)
> Return zero in the given $p$-adic field, to the default precision.
"""
function one(R::FlintPadicField)
   z = padic(R.prec_max)
   ccall((:padic_one, :libflint), Void, (Ref{padic},), z)
   z.parent = R
   return z
end

doc"""
    iszero(a::padic)
> Return `true` if the given p-adic field element is zero, otherwise return
> `false`.
"""
iszero(a::padic) = Bool(ccall((:padic_is_zero, :libflint), Cint,
                              (Ref{padic},), a))

doc"""
    isone(a::padic)
> Return `true` if the given p-adic field element is one, otherwise return
> `false`.
"""
isone(a::padic) = Bool(ccall((:padic_is_one, :libflint), Cint,
                             (Ref{padic},), a))

doc"""
    isunit(a::padic)
> Return `true` if the given p-adic field element is invertible, i.e. nonzero,
> otherwise return `false`.
"""
isunit(a::padic) = !Bool(ccall((:padic_is_zero, :libflint), Cint,
                              (Ref{padic},), a))

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function show(io::IO, x::padic)
   ctx = parent(x)
   cstr = ccall((:padic_get_str, :libflint), Ref{UInt8},
               (Ref{Void}, Ref{padic}, Ref{FlintPadicField}),
                   C_NULL, x, ctx)

   print(io, unsafe_string(cstr))

   ccall((:flint_free, :libflint), Void, (Ref{UInt8},), cstr)
   print(io, " + O(")
   print(io, prime(ctx))
   print(io, "^$(x.N))")
end

function show(io::IO, R::FlintPadicField)
   print(io, "Field of ", prime(R), "-adic numbers")
end

needs_parentheses(x::padic) = true

isnegative(x::padic) = false

show_minus_one(::Type{padic}) = true

###############################################################################
#
#   Canonicalisation
#
###############################################################################

canonical_unit(x::padic) = x

###############################################################################
#
#   Unary operators
#
###############################################################################

function -(x::padic)
   if iszero(x)
      return x
   end
   ctx = parent(x)
   z = padic(x.N)
   ccall((:padic_neg, :libflint), Void,
         (Ref{padic}, Ref{padic}, Ref{FlintPadicField}),
                     z, x, ctx)
   z.parent = ctx
   return z
end

###############################################################################
#
#   Binary operators
#
###############################################################################

function +(x::padic, y::padic)
   check_parent(x, y)
   ctx = parent(x)
   z = padic(min(x.N, y.N))
   z.parent = ctx
   ccall((:padic_add, :libflint), Void,
         (Ref{padic}, Ref{padic}, Ref{padic}, Ref{FlintPadicField}),
               z, x, y, ctx)
   return z
end

function -(x::padic, y::padic)
   check_parent(x, y)
   ctx = parent(x)
   z = padic(min(x.N, y.N))
   z.parent = ctx
   ccall((:padic_sub, :libflint), Void,
         (Ref{padic}, Ref{padic}, Ref{padic}, Ref{FlintPadicField}),
                  z, x, y, ctx)
   return z
end

function *(x::padic, y::padic)
   check_parent(x, y)
   ctx = parent(x)
   z = padic(min(x.N + y.v, y.N + x.v))
   z.parent = ctx
   ccall((:padic_mul, :libflint), Void,
         (Ref{padic}, Ref{padic}, Ref{padic}, Ref{FlintPadicField}),
               z, x, y, ctx)
   return z
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

+(a::padic, b::Integer) = a + parent(a)(b)

+(a::padic, b::fmpz) = a + parent(a)(b)

+(a::padic, b::fmpq) = a + parent(a)(b)

+(a::Integer, b::padic) = b + a

+(a::fmpz, b::padic) = b + a

+(a::fmpq, b::padic) = b + a

-(a::padic, b::Integer) = a - parent(a)(b)

-(a::padic, b::fmpz) = a - parent(a)(b)

-(a::padic, b::fmpq) = a - parent(a)(b)

-(a::Integer, b::padic) = parent(b)(a) - b

-(a::fmpz, b::padic) = parent(b)(a) - b

-(a::fmpq, b::padic) = parent(b)(a) - b

*(a::padic, b::Integer) = a*parent(a)(b)

*(a::padic, b::fmpz) = a*parent(a)(b)

*(a::padic, b::fmpq) = a*parent(a)(b)

*(a::Integer, b::padic) = b*a

*(a::fmpz, b::padic) = b*a

*(a::fmpq, b::padic) = b*a

###############################################################################
#
#   Comparison
#
###############################################################################

function ==(a::padic, b::padic)
   check_parent(a, b)
   ctx = parent(a)
   z = padic(min(a.N, b.N))
   ccall((:padic_sub, :libflint), Void,
         (Ref{padic}, Ref{padic}, Ref{padic}, Ref{FlintPadicField}),
               z, a, b, ctx)
   return Bool(ccall((:padic_is_zero, :libflint), Cint,
                (Ref{padic},), z))
end

function isequal(a::padic, b::padic)
   if parent(a) != parent(b)
      return false
   end
   return a.N == b.N && a == b
end

###############################################################################
#
#   Ad hoc comparison
#
###############################################################################

==(a::padic, b::Integer) = a == parent(a)(b)

==(a::padic, b::fmpz) = a == parent(a)(b)

==(a::padic, b::fmpq) = a == parent(a)(b)

==(a::Integer, b::padic) = parent(b)(a) == b

==(a::fmpz, b::padic) = parent(b)(a) == b

==(a::fmpq, b::padic) = parent(b)(a) == b

###############################################################################
#
#   Powering
#
###############################################################################

function ^(a::padic, n::Int)
   ctx = parent(a)
   z = padic(a.N + (n - 1)*a.v)
   z.parent = ctx
   ccall((:padic_pow_si, :libflint), Void,
                (Ref{padic}, Ref{padic}, Int, Ref{FlintPadicField}),
               z, a, n, ctx)
   return z
end

###############################################################################
#
#   Exact division
#
###############################################################################

function divexact(a::padic, b::padic)
   iszero(b) && throw(DivideError())
   check_parent(a, b)
   ctx = parent(a)
   z = padic(min(a.N - b.v, b.N - 2*b.v + a.v))
   z.parent = ctx
   ccall((:padic_div, :libflint), Cint,
         (Ref{padic}, Ref{padic}, Ref{padic}, Ref{FlintPadicField}),
               z, a, b, ctx)
   return z
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

divexact(a::padic, b::Integer) = a*(fmpz(1)//fmpz(b))

divexact(a::padic, b::fmpz) = a*(1//b)

divexact(a::padic, b::fmpq) = a*inv(b)

divexact(a::Integer, b::padic) = fmpz(a)*inv(b)

divexact(a::fmpz, b::padic) = inv((fmpz(1)//a)*b)

divexact(a::fmpq, b::padic) = inv(inv(a)*b)

###############################################################################
#
#   Inversion
#
###############################################################################

doc"""
    inv(a::padic)
> Returns $a^{-1}$. If $a = 0$ a `DivideError()` is thrown.
"""
function inv(a::padic)
   iszero(a) && throw(DivideError())
   ctx = parent(a)
   z = padic(a.N - 2*a.v)
   z.parent = ctx
   ccall((:padic_inv, :libflint), Cint,
         (Ref{padic}, Ref{padic}, Ref{FlintPadicField}), z, a, ctx)
   return z
end

###############################################################################
#
#   Divides
#
###############################################################################

doc"""
    divides(f::padic, g::padic)
> Returns a pair consisting of a flag which is set to `true` if $g$ divides
> $f$ and `false` otherwise, and a value $h$ such that $f = gh$ if
> such a value exists. If not, the value of $h$ is undetermined.
"""
function divides(a::padic, b::padic)
   iszero(b) && throw(DivideError())
   return true, divexact(a, b)
end

###############################################################################
#
#   GCD
#
###############################################################################

doc"""
    gcd(x::padic, y::padic)
> Returns the greatest common divisor of $x$ and $y$, i.e. the function returns
> $1$ unless both $a$ and $b$ are $0$, in which case it returns $0$.
"""
function gcd(x::padic, y::padic)
   check_parent(x, y)
   if iszero(x) && iszero(y)
      z = zero(parent(x))
   else
      z = one(parent(x))
   end
   return z
end

###############################################################################
#
#   Square root
#
###############################################################################

doc"""
    sqrt(a::padic)
> Return the $p$-adic square root of $a$. We define this only when the
> valuation of $a$ is even. The precision of the output will be
> precision$(a) -$ valuation$(a)/2$. If the square root does not exist, an
> exception is thrown.
"""
function sqrt(a::padic)
   (a.v % 2) != 0 && error("Unable to take padic square root")
   ctx = parent(a)
   z = padic(a.N - div(a.v, 2))
   z.parent = ctx
   res = Bool(ccall((:padic_sqrt, :libflint), Cint,
                    (Ref{padic}, Ref{padic}, Ref{FlintPadicField}), z, a, ctx))
   !res && error("Square root of p-adic does not exist")
   return z
end

###############################################################################
#
#   Special functions
#
###############################################################################

doc"""
    exp(a::padic)
> Return the $p$-adic exponential of $a$. We define this only when the
> valuation of $a$ is positive (unless $a = 0$). The precision of the output
> will be the same as the precision of the input. If the input is not valid an
> exception is thrown.
"""
function Base.exp(a::padic)
   !iszero(a) && a.v <= 0 && throw(DomainError())
   ctx = parent(a)
   z = padic(a.N)
   z.parent = ctx
   res = Bool(ccall((:padic_exp, :libflint), Cint,
                    (Ref{padic}, Ref{padic}, Ref{FlintPadicField}), z, a, ctx))
   !res && error("Unable to compute exponential")
   return z
end

doc"""
    log(a::padic)
> Return the $p$-adic logarithm of $a$. We define this only when the valuation
> of $a$ is zero (but not for $a == 0$). The precision of the output will be
> the same as the precision of the input. If the input is not valid an
> exception is thrown.
"""
function log(a::padic)
   (a.v > 0 || a.v < 0 || iszero(a)) && throw(DomainError())
   ctx = parent(a)
   z = padic(a.N)
   z.parent = ctx
   res = Bool(ccall((:padic_log, :libflint), Cint,
                    (Ref{padic}, Ref{padic}, Ref{FlintPadicField}), z, a, ctx))
   !res && error("Unable to compute logarithm")
   return z
end

doc"""
    teichmuller(a::padic)
> Return the Teichmuller lift of the $p$-adic value $a$. We require the
> valuation of $a$ to be nonnegative. The precision of the output will be the
> same as the precision of the input. For convenience, if $a$ is congruent to
> zero modulo $p$ we return zero. If the input is not valid an exception is
> thrown.
"""
function teichmuller(a::padic)
   a.v < 0 && throw(DomainError())
   ctx = parent(a)
   z = padic(a.N)
   z.parent = ctx
   ccall((:padic_teichmuller, :libflint), Void,
         (Ref{padic}, Ref{padic}, Ref{FlintPadicField}), z, a, ctx)
   return z
end

###############################################################################
#
#   Unsafe operators
#
###############################################################################

function zero!(z::padic)
   z.N = parent(z).prec_max
   ctx = parent(z)
   ccall((:padic_zero, :libflint), Void,
         (Ref{padic}, Ref{FlintPadicField}), z, ctx)
   return z
end

function mul!(z::padic, x::padic, y::padic)
   z.N = min(x.N + y.v, y.N + x.v)
   ctx = parent(x)
   ccall((:padic_mul, :libflint), Void,
         (Ref{padic}, Ref{padic}, Ref{padic}, Ref{FlintPadicField}),
               z, x, y, ctx)
   return z
end

function addeq!(x::padic, y::padic)
   x.N = min(x.N, y.N)
   ctx = parent(x)
   ccall((:padic_add, :libflint), Void,
         (Ref{padic}, Ref{padic}, Ref{padic}, Ref{FlintPadicField}),
               x, x, y, ctx)
   return x
end

function addeq!(z::padic, x::padic, y::padic)
   z.N = min(x.N, y.N)
   ctx = parent(x)
   ccall((:padic_add, :libflint), Void,
         (Ref{padic}, Ref{padic}, Ref{padic}, Ref{FlintPadicField}),
               z, x, y, ctx)
   return z
end

###############################################################################
#
#   Conversions and promotions
#
###############################################################################

promote_rule(::Type{padic}, ::Type{T}) where {T <: Integer} = padic

promote_rule(::Type{padic}, ::Type{fmpz}) = padic

promote_rule(::Type{padic}, ::Type{fmpq}) = padic

###############################################################################
#
#   Parent object overloads
#
###############################################################################

function (R::FlintPadicField)()
   z = padic(R.prec_max)
   z.parent = R
   return z
end

function (R::FlintPadicField)(n::fmpz)
   if isone(n)
      N = 0
   else
      p = prime(R)
      N, = remove(n, p)
   end
   z = padic(N + R.prec_max)
   ccall((:padic_set_fmpz, :libflint), Void,
         (Ref{padic}, Ref{fmpz}, Ref{FlintPadicField}), z, n, R)
   z.parent = R
   return z
end

function (R::FlintPadicField)(n::fmpq)
   m = denominator(n)
   if isone(m)
      return R(numerator(n))
   end
   p = prime(R)
   if m == p
      N = -1
   else
     N = -flog(m, p)
   end
   z = padic(N + R.prec_max)
   ccall((:padic_set_fmpq, :libflint), Void,
         (Ref{padic}, Ref{fmpq}, Ref{FlintPadicField}), z, n, R)
   z.parent = R
   return z
end

(R::FlintPadicField)(n::Integer) = R(fmpz(n))

function (R::FlintPadicField)(n::padic)
   parent(n) != R && error("Unable to coerce into p-adic field")
   return n
end

###############################################################################
#
#   FlintPadicField constructor
#
###############################################################################

# inner constructor is also used directly

doc"""
    FlintPadicField(p::Integer, prec::Int)
> Returns the parent object for the $p$-adic field for given prime $p$, where
> the default absolute precision of elements of the field is given by `prec`.
"""
function FlintPadicField(p::Integer, prec::Int)
   return FlintPadicField(fmpz(p), prec)
end


text:
function update_v!(state::State, data::Data, prior::Prior)
  for i in samplenames
    update_v!(i, state, data, prior)
  end
end

function update_v!(i::Char, state::State, data::Data, prior::Prior)
  lami = ref_lambda(state, i)
  zetai = ref_zeta(state, i)
  yi = ref_yfinite(data, i)
  vi = ref_v(state, i)

  nu = state.nu[lami]
  omega = state.omega[lami]
  psi = state.psi[lami]
  mu = state.mu[lami]

  shape = nu/2 .+ 1
  rate = (nu + zetai.^2 + ((yi - mu - psi .* zetai) .^ 2) ./ omega) / 2

  for n in eachindex(vi)
    vi[n] = rand(Gamma(shape[n], 1/rate[n]))
  end
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
using SafeTestsets

@safetestset "Atomic Solutions                                                                " begin include("atomic_solutions_tests.jl") end
@safetestset "Deterministic Solutions                                                         " begin include("deterministic_solutions_tests.jl") end


text:
using Arborist: wherevar, is_allowable_generated_body, insert_function_preamble
using MacroTools

@testset "wherevar" begin
    @test :T==wherevar(MacroTools.splitdef(:(foo(x::T) where {T} = x))[:whereparams][1])
    @test :T==wherevar(MacroTools.splitdef(:(foo(x::T) where {T<:Int} = x))[:whereparams][1])
end

@testset "insert_function_preamble" begin
    call_args_typetuple = Tuple{Int}
    ast = :(foo(x)=x)
    def = MacroTools.splitdef(ast)
    body = copy(def[:body])
    neobody = insert_function_preamble(
        body, def, call_args_typetuple
    ) |> MacroTools.flatten

    @test neobody == quote
        (x,) = call_args
        x
    end
end

@testset "is_allowable_generated_body" begin
    @testset "negatives" begin
        @test false == is_allowable_generated_body(
            :(x->2x)
        )

        @test false == is_allowable_generated_body(
            :((x for x in 1:20))
        )

        @test false == is_allowable_generated_body(
            quote
                y=0
                f = () -> y+=1  # NOT ALLOWED, closure
                f()
                return y
            end
        )
        @test false == is_allowable_generated_body(
            quote
                y=sum([1 for ii in 1:3]) ## comprehension
                return y
            end
        )

        @test false == is_allowable_generated_body(
            quote
                global x
                x=5
            end
        )
    end # @testset "negatives"

    @testset "positives" begin
        @test true == is_allowable_generated_body(
            :(2x)
        )

        @test true == is_allowable_generated_body(
            quote
                y=0
                f = Base.Fix1(+, y)
                f()
                return y
            end
        )
        @test true == is_allowable_generated_body(
            quote
                y=sum(identity, 1:3)
                return y
            end
        )
        @test true == is_allowable_generated_body(
            quote
                local x
                x=5
            end
        )
    end  # @testset positities
end


text:
"""
    nvertices(meshes::RefMesh)

Return the number of vertices of a reference mesh
"""
function Meshes.nvertices(mesh::RefMesh)
    Meshes.nvertices(mesh.mesh)
end

"""
    nelements(meshes::RefMeshes)

Return the number of elements of a reference mesh
"""
function Meshes.nelements(mesh::RefMesh)
    Meshes.nelements(mesh.mesh)
end

function normals(mesh::RefMesh{S,ME,M,N,T}) where {S,ME<:Meshes.SimpleMesh,M,N,T}
    if length(mesh.normals) == 0
        return SVector{length(Meshes.topology(mesh.mesh).connec)}(
            Meshes.Point3(Meshes.normal(Meshes.Triangle(mesh.mesh.points[[tri.indices...]]))) for tri in Meshes.topology(mesh.mesh).connec
        )
    else
        return mesh.normals
    end
    # TODO: Implement for RefMesh with GeometryBasics
end

"""
    normals_vertex(mesh::Meshes.SimpleMesh)

Compute per vertex normals and return them as a `StaticArrays.SVector`.

#! This is a naive approach because I have no time right know.
#! We just put the face mesh as a vertex mesh (and ovewritting values for common points)
# TODO: Use a real computation instead. See e.g.:
# https://stackoverflow.com/questions/45477806/general-method-for-calculating-smooth-vertex-normals-with-100-smoothness?noredirect=1&lq=1
"""
function normals_vertex(mesh::RefMesh)
    vertex_normals = fill(Meshes.Point3(0.0, 0.0, 0.0), Meshes.nvertices(mesh))
    for (i, tri) in enumerate(Meshes.topology(mesh.mesh).connec)
        vertex_normals[tri.indices[1]] = mesh.normals[i]
        vertex_normals[tri.indices[2]] = mesh.normals[i]
        vertex_normals[tri.indices[3]] = mesh.normals[i]
    end

    return SVector{length(vertex_normals)}(vertex_normals)
end

function normals_vertex(mesh::Meshes.SimpleMesh)
    vertex_normals = fill(Meshes.Point3(0.0, 0.0, 0.0), Meshes.nvertices(mesh))
    for (i, tri) in enumerate(Meshes.topology(mesh).connec)
        tri_norm = Meshes.normal(Meshes.Triangle(mesh.points[[tri.indices...]]))
        vertex_normals[tri.indices[1]] = tri_norm
        vertex_normals[tri.indices[2]] = tri_norm
        vertex_normals[tri.indices[3]] = tri_norm
    end

    return SVector{length(vertex_normals)}(vertex_normals)
end

"""
    nvertices(meshes::RefMeshes)

Return the number of vertices for each reference mesh as a vector of nvertices
"""
function Meshes.nvertices(meshes::RefMeshes)
    [Meshes.nvertices(i) for i in meshes.meshes]
    # TODO: Implement for RefMesh with GeometryBasics
end

"""
    nelements(meshes::RefMeshes)

Return the number of elements for each reference mesh as a vector of nelements
"""
function Meshes.nelements(meshes::RefMeshes)
    [Meshes.nelements(i) for i in meshes.meshes]
    # TODO: Implement for RefMesh with GeometryBasics
end


text:
using BitStreamLib

mutable struct Node
	value::UInt8
	leftptr::Union{Node, Nothing}
	rightptr::Union{Node, Nothing}
	Node() = new(0, nothing, nothing)
	Node(x, y, z) = new(x, y, z)
end

function min_heapify!(a::Array{Tuple{Int, Node}}, s::Int, i::Int)
	l = i * 2
	r = i * 2 + 1
	m = 0
	if l <= s && a[l][1] < a[i][1]
		m = l
	else
		m = i
	end

	if r <= s && a[r][1] < a[m][1]
		m = r
	end

	if m != i
		a[i], a[m] = a[m], a[i]
		min_heapify!(a, s, m)
	end
end

function build_queue(elements::Dict{UInt8, Int})::Array{Tuple{Int, Node}, 1}
	queue::Array{Tuple{Int, Node}, 1} = Array{Tuple{Int, Node}, 1}(undef, length(elements))
	s::Int = length(elements)

	el_vals = collect(values(elements))
	el_keys = collect(keys(elements))

	for i in 1 : s
		queue[i] = el_vals[i], Node(el_keys[i], nothing, nothing)
	end

	for i in cld(s, 2) : -1 : 1
		min_heapify!(queue, s, i)
	end

	return queue
end

function dequeue(queue::Array{Tuple{Int, Node}, 1}, s::Int)::Tuple{Int, Node}
	ret = queue[1]
	queue[1], queue[s] = queue[s], queue[1]
	min_heapify!(queue, s - 1, 1)
	return ret
end

function enqueue(queue::Array{Tuple{Int, Node}, 1}, s::Int, val::Tuple{Int, Node})
	if s == length(queue)
		push!(queue, val)
	else
		queue[s + 1] = val
	end
	i::Int = s + 1
	while i > 1 && queue[i][1] < queue[fld(i, 2)][1]
		p::Int = fld(i, 2)
		queue[i], queue[p] = queue[p], queue[i]
		i = p
	end
end

function make_huff_tree(data::Array{UInt8, 1})::Node
	count_list::Dict{UInt8, Int} = Dict()

	for i in data
		if !haskey(count_list, i)
			push!(count_list, i => 1)
		else
			count_list[i] += 1
		end
	end

	queue = build_queue(count_list)
	s = length(queue)

	while s > 1
		n1 = dequeue(queue, s)
		n2 = dequeue(queue, s - 1)
		n3 = (n1[1] + n2[1], Node(0, n1[2], n2[2]))
		enqueue(queue, s - 2, n3)
		s -= 1
	end
	return queue[1][2]
end

function make_huff_tree(data::String)::Node
	return make_huff_tree(translate(data))
end

function translate(data::String)::Array{UInt8, 1}
	int_data = Array{UInt8, 1}(undef, length(data))

	for i in 1 : length(data)
		int_data[i] = Int(data[i])
	end

	return int_data
end

function translate(data::Array{UInt8, 1})::String
	ret = ""

	for i in data
		ret *= Char(i)
	end

	return ret
end

function generate_mapping!(encode_table::Dict{UInt8, String}, cur_node::Node, cur_bin::Array{Bool, 1}, s::Int)
	if cur_node.leftptr == nothing
		bin = ""
		for i in 1 : s
			bin *= cur_bin[i] ? "1" : "0"
		end
		push!(encode_table, cur_node.value => bin)
		return
	end

	if s == length(cur_bin)
		push!(cur_bin, false)
	else
		cur_bin[s + 1] = false
	end
	generate_mapping!(encode_table, cur_node.leftptr, cur_bin, s + 1)

	cur_bin[s + 1] = true
	generate_mapping!(encode_table, cur_node.rightptr, cur_bin, s + 1)
end

function generate_mapping(huff_tree::Node)::Dict
	encode_table::Dict{UInt8, String} = Dict()

	generate_mapping!(encode_table, huff_tree, Array{Bool, 1}(undef, 4), 0)

	return encode_table
end

function get_decode_map(encode_table::Dict{UInt8, String})::Dict{String, UInt8}
	ret::Dict{String, UInt8} = Dict()

	for i in keys(encode_table)
		push!(ret, encode_table[i] => i)
	end

	return ret
end

function huff_encode(encode_table::Dict{UInt8, String}, data::String)::String
	ret = ""

	for i in data
		ret *= encode_table[UInt8(i)]
	end

	return ret
end

function huff_decode_str(decode_table::Dict{String, UInt8}, bin_data::String)::String
	ret = ""

	cur_bin = ""

	for i in bin_data
		cur_bin *= i
		if haskey(decode_table, cur_bin)
			ret *= Char(decode_table[cur_bin])
			cur_bin = ""
		end
	end

	return ret
end

function huff_encode(encode_table::Dict{UInt8, String}, data::Array{UInt8, 1})::BitStream
	ret = BitStream()
	for i in data
		temp_dat = encode_table[i]
		for j in temp_dat
			bit::Bool = j == '1'
			push!(ret, bit)
		end
	end
	return ret
end

function huff_encode(enc_tree::Node, data::Array{UInt8, 1})::BitStream
	return huff_encode(generate_mapping(enc_tree), data)
end

function huff_encode(enc_tree::Node, data::String)::BitStream
	return huff_encode(enc_tree, translate(data))
end

function huff_decode(decode_table::Dict{String, UInt8}, bin_data::BitStream)::Array{UInt8, 1}
	ret::Array{UInt8, 1} = Array{UInt8, 1}(undef, 0)

	cur_bin::String = ""

	for i in 1 : length(bin_data)
		cur_bin *= get(bin_data, i) ? '1' : '0'

		if haskey(decode_table, cur_bin)
			push!(ret, decode_table[cur_bin])
			cur_bin = ""
		end
	end

	return ret
end

function huff_decode(dec_tree::Node, bin_data::BitStream)::Array{UInt8, 1}
	return huff_decode(get_decode_map(generate_mapping(dec_tree)), bin_data)
end

function huff_decode_str(decode_table::Dict{String, UInt8}, bin_data::BitStream)::String
	return translate(huff_decode(decode_table, bin_data))
end

function read_file(file_name::String)::Array{UInt8, 1}
	strm = open(file_name)
	seekend(strm)
	eof = position(strm)
	seekstart(strm)

	ret::Array{UInt8, 1} = Array{UInt8, 1}(undef, eof)

	readbytes!(strm, ret)

	close(strm)

	return ret
end

function write_file(file_name::String, data::Array{UInt8, 1})
	io = open(file_name, "w")

	write(io, data)

	close(io)
end

function compress_tree(tree::Node, save::Array{UInt8, 1}, s::Int)
	if tree.leftptr == nothing
		save[s] = 0
		save[s + 1] = tree.value
	else
		save[s] = UInt8(fld(length(save), 2))
		push!(save, 0)
		push!(save, 0)
		compress_tree(tree.leftptr, save, length(save) - 1)
		save[s + 1] = UInt8(fld(length(save), 2))
		push!(save, 0)
		push!(save, 0)
		compress_tree(tree.rightptr, save, length(save) - 1)
	end
end

function compress_tree(tree::Node, save::Array{UInt16, 1}, s::Int)
	if tree.leftptr == nothing
		save[s] = 0
		save[s + 1] = tree.value
	else
		save[s] = UInt16(fld(length(save), 2))
		push!(save, 0)
		push!(save, 0)
		compress_tree(tree.leftptr, save, length(save) - 1)
		save[s + 1] = UInt16(fld(length(save), 2))
		push!(save, 0)
		push!(save, 0)
		compress_tree(tree.rightptr, save, length(save) - 1)
	end
end

function to_byte_array(data::Array{UInt16, 1})::Array{UInt8}
	ret::Array{UInt8, 1} = zeros(UInt8, length(data) * 2)

	for i in 1 : length(data)
		ret[(i - 1) * 2 + 1] = UInt8(data[i] & 255)
		ret[i * 2] = UInt8((data[i] & (255 << 8)) >> 8)
	end

	return ret
end

function compress_tree(tree::Node)::Union{Array{UInt8, 1}, Array{UInt16, 1}}
	if length(generate_mapping(tree)) > 128
		save2::Array{UInt16, 1} = UInt16[0, 0]
		compress_tree(tree, save2, 1)
		return save2
	else
		save::Array{UInt8, 1} = UInt8[0, 0]
		compress_tree(tree, save, 1)
		return save
	end
end

function decompress_tree(comp_tree::Union{Array{UInt8, 1}, Array{UInt16, 1}}, s::Int, cur_node::Node)
	if comp_tree[s] == 0
		cur_node.leftptr = nothing
		cur_node.rightptr = nothing
		cur_node.value = comp_tree[s + 1]
	else
		cur_node.leftptr = Node()
		decompress_tree(comp_tree, comp_tree[s] * 2 + 1, cur_node.leftptr)
		cur_node.rightptr = Node()
		decompress_tree(comp_tree, comp_tree[s + 1] * 2 + 1, cur_node.rightptr)
	end
end

function decompress_tree(comp_tree::Union{Array{UInt8, 1}, Array{UInt16, 1}})::Node
	root::Node = Node()
	decompress_tree(comp_tree, 1, root)
	return root
end

function write_to_file(filename::String, data::BitStream, tree::Array{UInt8, 1})
	io = open(filename, "w")

	write(io, Int64(length(data)))

	write(io, UInt16(length(tree)))

	write(io, tree)

	write(io, data.data)

	close(io)
end

function write_to_file(filename::String, data::BitStream, tree::Array{UInt16, 1})
	write_to_file(filename, data, to_byte_array(tree))
end

function write_to_file(filename::String, data::BitStream, tree::Node)
	write_to_file(filename, data, compress_tree(tree))
end

function compress_file(filein::String, fileout::String)
	data = read_file(filein)

	tree = make_huff_tree(data)

	comp_data = huff_encode(tree, data)

	write_to_file(fileout, comp_data, compress_tree(tree))
end

function decode_file(filename::String)::Array{UInt8, 1}
	io = open(filename, "r")

	bit_length::Int64 = read(io, Int64)

	tree_length::UInt16 = read(io, UInt16)

	comp_tree::Array{UInt8, 1} = Array{UInt8, 1}(undef, tree_length)

	readbytes!(io, comp_tree)

	comp_data::Array{UInt8, 1} = Array{UInt8, 1}(undef, cld(bit_length, 8))

	readbytes!(io, comp_data)

	close(io)

	comp_data_strm::BitStream = BitStream(comp_data, bit_length)

	println("hei1")

	tree::Node = decompress_tree(comp_tree)

	println("hei2")

	data = huff_decode(tree, comp_data_strm)

	return data
end

function decompress_file(filein::String, fileout::String)
	write_file(fileout, decode_file(filein))
end


text:
using PySerial
using Base.Test

@test typeof(list_ports()) == Array{Any,1}


text:
using RDKitMinimalLib
using GeometricFlux
using JSON
using GeometricFlux
using DataFrames
using LightGraphs
using Serialization

function process_atom_data(lines)
    matx = zeros(Float64, 4, length(lines))
    for (i,l) in enumerate(lines)
        coord_strs = strip(l[2:end]) |> x -> split(x, '\t')
        coord_strs = coord_strs[coord_strs .!= ""]
        coords = map(x -> parse(Float64, replace(x, "*^" => "e")
                               ), coord_strs)
        @assert length(coords) == 4
        matx[:,i] = coords
    end
    return matx
end

function get_df(filename)
    raw_lines = open(filename) do io
        readlines(io, keep=false)
    end
    n_atoms = parse(Int, raw_lines[1])

    props_meta = ["tag+i", "A", "B", "C", "mu", "alpha", "eps_homo",
                  "eps_lumo", "eps_gap", "R2", "zpve", "U0", "U",
                  "H", "G", "Cv"]
    raw_props = split(raw_lines[2], '\t')
    props_dict = Dict(zip(props_meta, raw_props))
    raw_atom_data = raw_lines[3:3+n_atoms-1]
    raw_numbers = raw_lines[3+n_atoms]
    smiles = split(raw_lines[4+n_atoms], '\t')

    mol_data = RDKitMinimalLib.get_mol(smiles[1]) |> RDKitMinimalLib.get_json |> JSON.parse
    moljson = mol_data["molecules"][1]

    mol_graph = SimpleGraph(n_atoms)

    num_bonds = length(moljson["bonds"])
    num_non_hydrogen = length(moljson["atoms"])

    for (i,bond) in enumerate(moljson["bonds"])
        add_edge!(mol_graph, bond["atoms"][1] + 1,  bond["atoms"][2] + 1)
    end

    Hcount = 0
    for (i,atom) in enumerate(moljson["atoms"])
        n_implicit_H = get(atom, "impHs", 0)
        for _ in 1:n_implicit_H
            add_edge!(mol_graph, i, num_non_hydrogen + 1 + Hcount)
            Hcount += 1
            num_bonds += 1
        end
    end
    @assert Hcount == n_atoms - num_non_hydrogen

    nf = process_atom_data(raw_atom_data)
    ef = zeros(1, num_bonds)
    bond_counter = 1
    for bond in moljson["bonds"]
        ef[bond_counter] = get(bond, "bo", 1)
        bond_counter += 1
    end

    for atom in moljson["atoms"]
        n_implicit_H = get(atom, "impHs", 0)
        for _ in 1:n_implicit_H
            ef[bond_counter] = 1
            bond_counter += 1
        end
    end
    mol_fg = FeaturedGraph(mol_graph, nf=nf, ef=ef)

    return Dict(
                props_dict...,
                "molecule" => mol_fg,
                "rawline" => raw_numbers,
           )
end

datarows = []

DATADIR = "/run/media/em/169747FC12AAEC72/datasets/qm9/"
DATADIR = "chemdata/"
for datafile in readdir(DATADIR)
    try
        push!(datarows, get_df(DATADIR * datafile))
    catch err
        println("Error in file $(datafile): $(err)")
    end
end

serialize("data.blob", datarows)


text:
struct MatsubaraPoleBasis <: AbstractBasis
    β::Float64
    poles::Vector{Float64}
end

function (basis::MatsubaraPoleBasis)(n::Vector{<:Integer})
    iv = (im * π / getbeta(basis)) .* n
    return 1 ./ (transpose(iv) .- basis.poles)
end

struct TauPoleBasis <: AbstractBasis
    β::Float64
    poles::Vector{Float64}
    statistics::Statistics
    wmax::Float64
end

getwmax(basis::TauPoleBasis) = basis.wmax

function TauPoleBasis(beta::Real, statistics::Statistics, poles::Vector{<:AbstractFloat})
    return TauPoleBasis(beta, poles, statistics, maximum(abs, poles))
end

function (basis::TauPoleBasis)(tau::Vector{<:AbstractFloat})
    all(τ -> 0 ≤ τ ≤ getbeta(basis), tau) ||
        throw(DomainError(tau, "tau must be in [0, beta]!"))

    x = (2 / getbeta(basis)) .* tau .- 1
    y = basis.poles ./ getwmax(basis)
    Λ = getbeta(basis) * getwmax(basis)
    if basis.statistics == fermion
        res = -LogisticKernel(Λ).(x, transpose(y))
    else
        K = RegularizedBoseKernel(Λ)
        res = -K.(x, transpose(y)) ./ transpose(y)
    end
    return transpose(res)
end

"""
Sparse pole representation
"""
struct SparsePoleRepresentation{T<:AbstractFloat} <: AbstractBasis
    basis::AbstractBasis
    poles::Vector{T}
    u::TauPoleBasis
    uhat::MatsubaraPoleBasis
    statistics::Statistics
    fitmat::Matrix{Float64}
    matrix::SVD
end

function Base.show(io::IO, obj::SparsePoleRepresentation{T}) where {T<:AbstractFloat}
    return print(io, "SparsePoleRepresentation for $(obj.basis) with poles at $(obj.poles)")
end

function SparsePoleRepresentation(
    basis::AbstractBasis, poles=default_omega_sampling_points(basis)
)
    y_sampling_points = poles ./ getwmax(basis)
    u = TauPoleBasis(getbeta(basis), basis.statistics, poles)
    uhat = MatsubaraPoleBasis(getbeta(basis), poles)
    weight = weight_func(basis.kernel, basis.statistics)(y_sampling_points)
    fitmat = -basis.s .* basis.v(poles) .* transpose(weight)
    matrix = svd(fitmat)
    return SparsePoleRepresentation(basis, poles, u, uhat, basis.statistics, fitmat, matrix)
end

getbeta(obj::SparsePoleRepresentation) = getbeta(obj.basis)

function Base.getproperty(obj::SparsePoleRepresentation, d::Symbol)
    if d === :v
        return nothing
    else
        return getfield(obj, d)
    end
end

function default_tau_sampling_points(obj::SparsePoleRepresentation)
    return default_tau_sampling_points(obj.basis)
end

function default_matsubara_sampling_points(obj::SparsePoleRepresentation)
    return default_matsubara_sampling_points(obj.basis)
end

iswellconditioned(::SparsePoleRepresentation) = false

"""
From IR to SPR

gl:
    Expansion coefficients in IR
"""
function from_IR(spr::SparsePoleRepresentation, gl::AbstractArray, dims=1)
    return mapslices(i -> spr.matrix \ i, gl; dims)
end

"""
From SPR to IR

g_spr:
    Expansion coefficients in SPR
"""
function to_IR(spr::SparsePoleRepresentation, g_spr::AbstractArray, dims=1)
    return mapslices(i -> spr.fitmat * i, g_spr; dims)
end


text:
using Test
using BenchmarkTools

function fink_filter_bank_AKS(c, L; nx=256, wd=1, pc=1, shift=false, Omega=false, safety_on=true)
    #c     - sets the scale sampling rate (1 is dyadic, 2 is half dyadic)
    #L     - number of angular bins (usually 8*pc or 16*pc)
    #wd    - width of the wavelets (default 1, wd=2 for a double covering)
    #pc    - plane coverage (default 1, full 2pi 2)
    #shift - shift in θ by 1/2 of the θ spacing
    #Omega - true= append Omega filter (all power beyond Nyquist) so the sum of filters is 1.0

    # -------- assertion errors to make sure arguments are reasonable
    @test wd <= L/2

    # -------- set parameters
    dθ   = pc*π/L
    θ_sh = shift ? dθ/2 : 0.0
    dx   = nx/2-1

    im_scale = convert(Int8,log2(nx))
    # -------- number of bins in radial direction (size scales)
    J = (im_scale-2)*c

    # -------- allocate output array of zeros
    filt      = zeros(nx, nx, J*L+(Omega ? 2 : 1))
    psi_index = zeros(Int32, J, L)
    theta     = zeros(Float64, L)
    j_value   = zeros(Float64, J)

    # -------- compute the required wd
    j_rad_exp = zeros(J)
    for j_ind = 1:J
        j = j_ind/c
        jrad  = im_scale-j-1
        j_rad_exp[j_ind] = 2^(jrad)
    end

    wd_j = max.(ceil.(L./(pc.*π.*j_rad_exp)),wd)

    if !safety_on
        wd_j.=wd
    end

    # loop over wd from small to large
    ## there is some uneeded redundancy added by doing this esp in l loop
    for wd in sort(unique(wd_j))
        # -------- allocate theta and logr arrays
        θ    = zeros(nx, nx)
        logr = zeros(nx, nx)

        wdθ  = wd*dθ
        norm = 1.0/sqrt(wd)
        # -------- loop over l
        for l = 0:L-1
            θ_l        = dθ*l+θ_sh
            theta[l+1] = θ_l

        # -------- allocate anggood BitArray
            anggood = falses(nx, nx)

        # -------- loop over pixels
            for x = 1:nx
                sx = mod(x+dx,nx)-dx -1    # define sx,sy so that no fftshift() needed
                for y = 1:nx
                    sy = mod(y+dx,nx)-dx -1
                    θ_pix  = mod(atan(sy, sx)+π -θ_l, 2*π)
                    θ_good = abs(θ_pix-π) <= wdθ

                    # If this is a pixel we might use, calculate log2(r)
                    if θ_good
                        anggood[y, x] = θ_good
                        θ[y, x]       = θ_pix
                        r2            = sx^2 + sy^2
                        logr[y, x]    = 0.5*log2(max(1,r2))
                    end
                end
            end
            angmask = findall(anggood)
        # -------- compute the wavelet in the Fourier domain
        # -------- the angular factor is the same for all j
            F_angular = norm .* cos.((θ[angmask].-π).*(L/(2*wd*pc)))

        # -------- loop over j for the radial part
        #    for (j_ind, j) in enumerate(1/c:1/c:im_scale-2)
            j_ind_w_wd = findall(wd_j.==wd)
            for j_ind = j_ind_w_wd
                j = j_ind/c
                j_value[j_ind] = j  # store for later
                jrad  = im_scale-j-1
                Δj    = abs.(logr[angmask].-jrad)
                rmask = (Δj .<= 1/c)

        # -------- radial part
                F_radial = cos.(Δj[rmask] .* (c*π/2))
                ind      = angmask[rmask]
        #      Let's have these be (J,L) if you reshape...
        #        f_ind    = (j_ind-1)*L+l+1
                f_ind    = j_ind + l*J
                filt[ind, f_ind] = F_radial .* F_angular[rmask]
                psi_index[j_ind,l+1] = f_ind
            end
        end
    end

    # -------- phi contains power near k=0 not yet accounted for
    filter_power = (sum(filt.*filt, dims=3))[:,:,1]

    # -------- for plane half-covered (pc=1), add other half-plane
    if pc == 1
        filter_power .+= circshift(filter_power[end:-1:1,end:-1:1],(1,1))
    end

    # -------- compute power required to sum to 1.0
    i0 = round(Int16,nx/2-2)
    i1 = round(Int16,nx/2+4)
    center_power = 1.0 .- fftshift(filter_power)[i0:i1,i0:i1]
    zind = findall(center_power .< 1E-15)
    center_power[zind] .= 0.0  # set small numbers to zero
    phi_cen = zeros(nx, nx)
    phi_cen[i0:i1,i0:i1] = sqrt.(center_power)

    # -------- before adding ϕ to filter bank, renormalize ψ if pc=1
    if pc==1 filt .*= sqrt(2.0) end  # double power for half coverage

    # -------- add result to filter array
    phi_index  = J*L+1
    filt[:,:,phi_index] .= fftshift(phi_cen)

    # -------- metadata dictionary
    info=Dict()
    info["npix"]         = nx
    info["j_value"]      = j_value
    info["theta_value"]  = theta
    info["psi_index"]    = psi_index
    info["phi_index"]    = phi_index
    info["pc"]           = pc
    info["wd"]           = wd_j
    info["fs_center_r"]  = j_rad_exp

    if Omega     # append a filter containing the rest (outside Nyquist)
        filter_power += filt[:,:,phi_index].^2
        edge_power    = 1.0 .- filter_power
        zind          = findall(edge_power .< 1E-15)
        edge_power[zind]     .= 0.0  # set small numbers to zero
        Omega_index           = J*L+2
        info["Omega_index"]   = Omega_index
        filt[:,:,Omega_index] = sqrt.(edge_power)
    end

    return filt, info
end

using Plots; gr()
using Measures
using FFTW
push!(LOAD_PATH, "/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/main")
using DHC_2DUtils

function plot_filter_bank_QA(filt, info; fname="filter_bank_QA.png")

    # -------- define plot1 to append plot to a list
    function plot1(ps, image; clim=nothing, bin=1.0, fsz=12, label=nothing)
        # ps    - list of plots to append to
        # image - image to heatmap
        # clim  - tuple of color limits (min, max)
        # bin   - bin down by this factor, centered on nx/2+1
        # fsz   - font size
        marg   = -1mm
        nx, ny = size(image)
        nxb    = nx/round(Integer, 2*bin)

        # -------- center on nx/2+1
        i0 = max(1,round(Integer, (nx/2+2)-nxb-1))
        i1 = min(nx,round(Integer, (nx/2)+nxb+1))
        lims = [i0,i1]
        subim = image[i0:i1,i0:i1]
        push!(ps, heatmap(image,aspect_ratio=:equal,clim=clim,
            xlims=lims, ylims=lims, size=(400,400),
            legend=false,xtickfontsize=fsz,ytickfontsize=fsz,#tick_direction=:out,
            rightmargin=marg,leftmargin=marg,topmargin=marg,bottommargin=marg))
        if label != nothing
            annotate!(lims'*[.96,.04],lims'*[.09,.91],text(label,:left,:white,16))
        end
        return
    end

    # -------- initialize array of plots
    ps   = []
    jval = info["j_value"]
    J    = length(jval)
    L    = length(info["theta_value"])
    # plot the first 3 and last 3 j values
    jplt = [1:3;J-2:J]

    ind = zeros(Int64,J,L)
    for j_ind=1:J
        for l=1:L
            ind[j_ind,l] = Int(j_ind-1 + (l-1)*J)
        end
    end
    ind .+= 1

    # -------- loop over j
    for j_ind in jplt
        # first 3 filters in Fourier space
        j = jval[j_ind]
        bfac = max(1,2.0^(j-1))
        for l = 1:3
            label=string("j=",j_ind," ℓ=",l-1)
            #label=latexstring("\\tilde{\\psi}_{$j_ind,$(l-1)}")
            plot1(ps, fftshift(filt[:,:,ind[j_ind,l]]), clim=(0,1),
                  bin=bfac, label=label)
        end
        fac = max(1,2.0^(5-j))
        # real part of config space of 3rd filter
        plot1(ps, fftshift(real(ifft(filt[:,:,ind[j_ind,3]]))),
              bin=fac, label=string("j=",j_ind," ℓ=",3-1))

        # sum all filters to form (partial) ring
        ring = dropdims(sum(filt[:,:,ind[j_ind,:]].^2, dims=3), dims=3)
        plot1(ps, fftshift(ring[:,:,1]), clim=(0,1),
              bin=bfac, label=string("j=",j_ind," ℓ=",0,":",L-1))
    end

    wavepow = fftshift(dropdims(sum(filt[:,:,ind].^2, dims=(3,4)), dims=(3,4)))
    plot1(ps, wavepow, clim=(-0.1,1), label=string("j=1:",J," ℓ=0:",L-1))

    if info["pc"]==1
        wavepow += circshift(wavepow[end:-1:1,end:-1:1],(1,1))
    end
    plot1(ps, wavepow, clim=(-0.1,1), bin=32, label="missing power")
    phi = filt[:,:,info["phi_index"]]
    phi_shift = fftshift(phi)
    plot1(ps, phi_shift.^2, clim=(-0.1,1), bin=32, label="ϕ")
    disc = wavepow+(phi_shift.^2)
    plot1(ps, fftshift(real(ifft(phi))), label="ϕ")
    plot1(ps, disc, clim=(-0.1,1), bin=32, label="all")

    myplot = plot(ps..., layout=(7,5), size=(1400,2000))
    savefig(myplot, fname)
end

filt, info  = fink_filter_bank_AKS(1, 8, pc=1, wd=1, safety_on=false)

@benchmark filt, info  = fink_filter_bank(1, 8, pc=1, wd=1)

@benchmark filt, info  = fink_filter_bank_AKS(1, 8, pc=1, wd=1)

@benchmark filt, info  = fink_filter_bank_AKS(1, 8, pc=1, wd=1, safety_on=false)

function fink_filter_hash_AKS(c, L; nx=256, wd=1, pc=1, shift=false, Omega=false)
    # -------- compute the filter bank
    filt, hash = fink_filter_bank_AKS(c, L; nx=nx, wd=wd, pc=pc, shift=shift, Omega=Omega)

    # -------- list of non-zero pixels
    flist = fink_filter_list(filt)

    # -------- pack everything you need into the info structure
    hash["filt_index"] = flist[1]
    hash["filt_value"] = flist[2]
    return hash
end

filt, info  = fink_filter_bank_AKS(1, 8, pc=2, wd=1)

hash = fink_filter_hash_AKS(1, 8, pc=1, wd=1)

filt, info  = fink_filter_bank(1, 8, pc=2, wd=1)

theme(:dark)

plot_filter_bank_QA(filt, info, fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/filt-8-pc1-wd1-aks.png")

filt, info  = fink_filter_bank_AKS(1, 8, pc=1, wd=1)

plot_filter_bank_QA(filt, info, fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/filt-8-pc1-wd1-aks.png")

filt, info  = fink_filter_bank(1, 8, pc=1, wd=1)

plot_filter_bank_QA(filt, info, fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/filt-8-pc1-wd1.png")

filt, info  = fink_filter_bank_AKS(1, 8, pc=1, wd=2)

plot_filter_bank_QA(filt, info, fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/filt-8-pc1-wd2-aks.png")

filt, info  = fink_filter_bank(1, 8, pc=1, wd=2)

plot_filter_bank_QA(filt, info, fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/filt-8-pc1-wd2.png")

hash = fink_filter_hash_AKS(1, 8, pc=1, wd=1)

DHC_compute(rand(256,256),hash,hash)

filt, info  = fink_filter_bank_AKS(1, 16, pc=1, wd=2)
info

plot_filter_bank_QA(filt, info, fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/filt-16-pc1-wd2-aks.png")

filt, info  = fink_filter_bank_AKS(1, 32, pc=1, wd=2)
info

plot_filter_bank_QA(filt, info, fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/filt-32-pc1-wd2-aks.png")

filt, info  = fink_filter_bank_AKS(1, 64, pc=1, wd=2)
info

plot_filter_bank_QA(filt, info, fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/filt-64-pc1-wd2-aks.png")

filt, info  = fink_filter_bank_AKS(2, 64, pc=1, wd=2)
info

plot_filter_bank_QA(filt, info, fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/filt-2-64-pc1-wd2-aks.png")

filt, info  = fink_filter_bank_AKS(2, 8, pc=1, wd=2)
info

plot_filter_bank_QA(filt, info, fname="/Users/saydjari/Dropbox/GradSchool_AKS/Doug/Projects/DHC/scratch_AKS/filt-2-8-pc1-wd2-aks.png")

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC