text:
using Base.Meta

macro constants(array, stripprefix, expr)
    ret = Expr(:block)
    # Initialize the name lookup array
    push!(ret.args, :(const $array = Dict{UInt32,String}()))
    for e in expr.args
        if !isexpr(e, :const)
            continue
        end
        eq = e.args[1]
        @assert isexpr(eq, :(=))
        name = string(eq.args[1])
        name = replace(name, stripprefix => "")
        push!(ret.args, e)
        push!(ret.args, :($array[UInt32($(eq.args[1]))] = $name))
    end
    return esc(ret)
end

const EI_MAG0 = 0 #File identification byte 0 index
const ELFMAG0 = 0x7f #Magic number byte 0
const EI_MAG1 = 1 # File identification byte 1 index
const ELFMAG1 = UInt8('E') # Magic number byte 1
const EI_MAG2 = 2 # File identification byte 2 index
const ELFMAG2 = UInt8('L') # Magic number byte 2
const EI_MAG3 = 3 # File identification byte 3 index
const ELFMAG3 = UInt8('F') # Magic number byte 3

# e_indent[EIBCLASS]
@constants ELFCLASS "" begin
    const ELFCLASSNONE = 0 #Invalid class
    const ELFCLASS32 = 1 #32-bit objects
    const ELFCLASS64 = 2 #64-bit objects
end

# e_ident[EIDATA]
@constants ELFDATA "" begin
    const ELFDATANONE = 0
    const ELFDATA2LSB = 1
    const ELFDATA2MSB = 2
end

# e_version constants
@constants ELFVERSION "" begin
    const EV_NONE = 0       #Invalid Version
    const EV_CURRENT = 1    #Current Version
end

# e_type constants
@constants ET_TYPES "" begin
    const ET_NONE = 0   #No file type
    const ET_REL = 1   #Relocatable file
    const ET_EXEC = 2   #Executable file
    const ET_DYN = 3   #Shared object file
    const ET_CORE = 4   #Core file
end

@constants ET_RANGES "" begin
    const ET_LOOS = 0xfe00  #Operating system-specific
    const ET_HIOS = 0xfeff  #Operating system-specific
    const ET_LOPROC = 0xff00  #Processor-specific
    const ET_HIPROC = 0xffff  #Processor-specific
end

@constants EM_MACHINES "" begin
    # e_machine constants
    const EM_NONE = 0   #No machine
    const EM_M32 = 1   #AT&T WE 32100
    const EM_SPARC = 2   #SPARC
    const EM_386 = 3   #Intel 80386
    const EM_68K = 4   #Motorola 68000
    const EM_88K = 5   #Motorola 88000
    #6                  Reserved for future use (was EM_486)
    const EM_860 = 7   #Intel 80860
    const EM_MIPS = 8   #MIPS I Architecture
    const EM_S370 = 9   #IBM System/370 Processor
    const EM_MIPS_RS3_LE = 10  #MIPS RS3000 Little-endian
    #11-14              Reserved for future use
    const EM_PARISC = 15  #Hewlett-Packard PA-RISC
    #16  Reserved for future use
    const EM_VPP500 = 17  #Fujitsu VPP500
    const EM_SPARC32PLUS = 18  #Enhanced instruction set SPARC
    const EM_960 = 19  #Intel 80960
    const EM_PPC = 20  #PowerPC
    const EM_PPC64 = 21  #64-bit PowerPC
    const EM_S390 = 22  #IBM System/390 Processor
    const EM_SPU = 23  #IBM SPU/SPC
    #24-35   Reserved for future use
    const EM_V800 = 36  #NEC V800
    const EM_FR20 = 37  #Fujitsu FR20
    const EM_RH32 = 38  #TRW RH-32
    const EM_RCE = 39  #Motorola RCE
    const EM_ARM = 40  #ARM 32-bit architecture (AARCH32)
    const EM_ALPHA = 41  #Digital Alpha
    const EM_SH = 42  #Hitachi SH
    const EM_SPARCV9 = 43  #SPARC Version 9
    const EM_TRICORE = 44  #Siemens TriCore embedded processor
    const EM_ARC = 45  #Argonaut RISC Core, Argonaut Technologies Inc.
    const EM_H8_300 = 46  #Hitachi H8/300
    const EM_H8_300H = 47  #Hitachi H8/300H
    const EM_H8S = 48  #Hitachi H8S
    const EM_H8_500 = 49  #Hitachi H8/500
    const EM_IA_64 = 50  #Intel IA-64 processor architecture
    const EM_MIPS_X = 51  #Stanford MIPS-X
    const EM_COLDFIRE = 52  #Motorola ColdFire
    const EM_68HC12 = 53  #Motorola M68HC12
    const EM_MMA = 54  #Fujitsu MMA Multimedia Accelerator
    const EM_PCP = 55  #Siemens PCP
    const EM_NCPU = 56  #Sony nCPU embedded RISC processor
    const EM_NDR1 = 57  #Denso NDR1 microprocessor
    const EM_STARCORE = 58  #Motorola Star*Core processor
    const EM_ME16 = 59  #Toyota ME16 processor
    const EM_ST100 = 60  #STMicroelectronics ST100 processor
    const EM_TINYJ = 61  #Advanced Logic Corp. TinyJ embedded processor family
    const EM_X86_64 = 62  #AMD x86-64 architecture
    const EM_PDSP = 63  #Sony DSP Processor
    const EM_PDP10 = 64  #Digital Equipment Corp. PDP-10
    const EM_PDP11 = 65  #Digital Equipment Corp. PDP-11
    const EM_FX66 = 66  #Siemens FX66 microcontroller
    const EM_ST9PLUS = 67  #STMicroelectronics ST9+ 8/16 bit microcontroller
    const EM_ST7 = 68  #STMicroelectronics ST7 8-bit microcontroller
    const EM_68HC16 = 69  #Motorola MC68HC16 Microcontroller
    const EM_68HC11 = 70  #Motorola MC68HC11 Microcontroller
    const EM_68HC08 = 71  #Motorola MC68HC08 Microcontroller
    const EM_68HC05 = 72  #Motorola MC68HC05 Microcontroller
    const EM_SVX = 73  #Silicon Graphics SVx
    const EM_ST19 = 74  #STMicroelectronics ST19 8-bit microcontroller
    const EM_VAX = 75  #Digital VAX
    const EM_CRIS = 76  #Axis Communications 32-bit embedded processor
    const EM_JAVELIN = 77  #Infineon Technologies 32-bit embedded processor
    const EM_FIREPATH = 78  #Element 14 64-bit DSP Processor
    const EM_ZSP = 79  #LSI Logic 16-bit DSP Processor
    const EM_MMIX = 80  #Donald Knuth's educational 64-bit processor
    const EM_HUANY = 81  #Harvard University machine-independent object files
    const EM_PRISM = 82  #SiTera Prism
    const EM_AVR = 83  #Atmel AVR 8-bit microcontroller
    const EM_FR30 = 84  #Fujitsu FR30
    const EM_D10V = 85  #Mitsubishi D10V
    const EM_D30V = 86  #Mitsubishi D30V
    const EM_V850 = 87  #NEC v850
    const EM_M32R = 88  #Mitsubishi M32R
    const EM_MN10300 = 89  #Matsushita MN10300
    const EM_MN10200 = 90  #Matsushita MN10200
    const EM_PJ = 91  #picoJava
    const EM_OPENRISC = 92  #OpenRISC 32-bit embedded processor
    const EM_ARC_COMPACT = 93  #ARC International ARCompact processor (old spelling/synonym: EM_ARC_A5)
    const EM_XTENSA = 94  #Tensilica Xtensa Architecture
    const EM_VIDEOCORE = 95  #Alphamosaic VideoCore processor
    const EM_TMM_GPP = 96  #Thompson Multimedia General Purpose Processor
    const EM_NS32K = 97  #National Semiconductor 32000 series
    const EM_TPC = 98  #Tenor Network TPC processor
    const EM_SNP1K = 99  #Trebia SNP 1000 processor
    const EM_ST200 = 100 #STMicroelectronics (www.st.com) ST200 microcontroller
    const EM_IP2K = 101 #Ubicom IP2xxx microcontroller family
    const EM_MAX = 102 #MAX Processor
    const EM_CR = 103 #National Semiconductor CompactRISC microprocessor
    const EM_F2MC16 = 104 #Fujitsu F2MC16
    const EM_MSP430 = 105 #Texas Instruments embedded microcontroller msp430
    const EM_BLACKFIN = 106 #Analog Devices Blackfin (DSP) processor
    const EM_SE_C33 = 107 #S1C33 Family of Seiko Epson processors
    const EM_SEP = 108 #Sharp embedded microprocessor
    const EM_ARCA = 109 #Arca RISC Microprocessor
    const EM_UNICORE = 110 #Microprocessor series from PKU-Unity Ltd. and MPRC of Peking University
    const EM_EXCESS = 111 #eXcess: 16/32/64-bit configurable embedded CPU
    const EM_DXP = 112 #Icera Semiconductor Inc. Deep Execution Processor
    const EM_ALTERA_NIOS2 = 113 #Altera Nios II soft-core processor
    const EM_CRX = 114 #National Semiconductor CompactRISC CRX microprocessor
    const EM_XGATE = 115 #Motorola XGATE embedded processor
    const EM_C166 = 116 #Infineon C16x/XC16x processor
    const EM_M16C = 117 #Renesas M16C series microprocessors
    const EM_DSPIC30F = 118 #Microchip Technology dsPIC30F Digital Signal Controller
    const EM_CE = 119 #Freescale Communication Engine RISC core
    const EM_M32C = 120 #Renesas M32C series microprocessors
    #121-130 Reserved for future use
    const EM_TSK3000 = 131 #Altium TSK3000 core
    const EM_RS08 = 132 #Freescale RS08 embedded processor
    const EM_SHARC = 133 #Analog Devices SHARC family of 32-bit DSP processors
    const EM_ECOG2 = 134 #Cyan Technology eCOG2 microprocessor
    const EM_SCORE7 = 135 #Sunplus S+core7 RISC processor
    const EM_DSP24 = 136 #New Japan Radio (NJR) 24-bit DSP Processor
    const EM_VIDEOCORE3 = 137 #Broadcom VideoCore III processor
    const EM_LATTICEMICO32 = 138 #RISC processor for Lattice FPGA architecture
    const EM_SE_C17 = 139 #Seiko Epson C17 family
    const EM_TI_C6000 = 140 #The Texas Instruments TMS320C6000 DSP family
    const EM_TI_C2000 = 141 #The Texas Instruments TMS320C2000 DSP family
    const EM_TI_C5500 = 142 #The Texas Instruments TMS320C55x DSP family
    #143-159 Reserved for future use
    const EM_MMDSP_PLUS = 160 #STMicroelectronics 64bit VLIW Data Signal Processor
    const EM_CYPRESS_M8C = 161 #Cypress M8C microprocessor
    const EM_R32C = 162 #Renesas R32C series microprocessors
    const EM_TRIMEDIA = 163 #NXP Semiconductors TriMedia architecture family
    const EM_QDSP6 = 164 #QUALCOMM DSP6 Processor
    const EM_8051 = 165 #Intel 8051 and variants
    const EM_STXP7X = 166 #STMicroelectronics STxP7x family of configurable and extensible RISC processors
    const EM_NDS32 = 167 #Andes Technology compact code size embedded RISC processor family
    const EM_ECOG1 = 168 #Cyan Technology eCOG1X family
    const EM_ECOG1X = 168 #Cyan Technology eCOG1X family
    const EM_MAXQ30 = 169 #Dallas Semiconductor MAXQ30 Core Micro-controllers
    const EM_XIMO16 = 170 #New Japan Radio (NJR) 16-bit DSP Processor
    const EM_MANIK = 171 #M2000 Reconfigurable RISC Microprocessor
    const EM_CRAYNV2 = 172 #Cray Inc. NV2 vector architecture
    const EM_RX = 173 #Renesas RX family
    const EM_METAG = 174 #Imagination Technologies META processor architecture
    const EM_MCST_ELBRUS = 175 #MCST Elbrus general purpose hardware architecture
    const EM_ECOG16 = 176 #Cyan Technology eCOG16 family
    const EM_CR16 = 177 #National Semiconductor CompactRISC CR16 16-bit microprocessor
    const EM_ETPU = 178 #Freescale Extended Time Processing Unit
    const EM_SLE9X = 179 #Infineon Technologies SLE9X core
    const EM_L10M = 180 #Intel L10M
    const EM_K10M = 181 #Intel K10M
    #182 Reserved for future Intel use
    const EM_AARCH64 = 183 #ARM 64-bit architecture (AARCH64)
    #184 Reserved for future ARM use
    const EM_AVR32 = 185 #Atmel Corporation 32-bit microprocessor family
    const EM_STM8 = 186 #STMicroeletronics STM8 8-bit microcontroller
    const EM_TILE64 = 187 #Tilera TILE64 multicore architecture family
    const EM_TILEPRO = 188 #Tilera TILEPro multicore architecture family
    const EM_MICROBLAZE = 189 #Xilinx MicroBlaze 32-bit RISC soft processor core
    const EM_CUDA = 190 #NVIDIA CUDA architecture
    const EM_TILEGX = 191 #Tilera TILE-Gx multicore architecture family
    const EM_CLOUDSHIELD = 192 #CloudShield architecture family
    const EM_COREA_1ST = 193 #KIPO-KAIST Core-A 1st generation processor family
    const EM_COREA_2ND = 194 #KIPO-KAIST Core-A 2nd generation processor family
    const EM_ARC_COMPACT2 = 195 #Synopsys ARCompact V2
    const EM_OPEN8 = 196 #Open8 8-bit RISC soft processor core
    const EM_RL78 = 197 #Renesas RL78 family
    const EM_VIDEOCORE5 = 198 #Broadcom VideoCore V processor
    const EM_78KOR = 199 #Renesas 78KOR family
    const EM_56800EX = 200 #Freescale 56800EX Digital Signal Controller (DSC)
    const EM_BA1 = 201 #Beyond BA1 CPU architecture
    const EM_BA2 = 202 #Beyond BA2 CPU architecture
    const EM_XCORE = 203 #XMOS xCORE processor family
    const EM_MCHP_PIC = 204 #Microchip 8-bit PIC(r) family
end

#EI_OSABI constants
@constants ELFOSABI "" begin
    const ELFOSABI_NONE = 0   #No extensions or unspecified
    const ELFOSABI_HPUX = 1   #Hewlett-Packard HP-UX
    const ELFOSABI_NETBSD = 2   #NetBSD
    const ELFOSABI_GNU = 3   #GNU
    const ELFOSABI_LINUX = 3   ##Linux historical - alias for ELFOSABI_GNU
    const ELFOSABI_SOLARIS = 6   #Sun Solaris
    const ELFOSABI_AIX = 7   #AIX
    const ELFOSABI_IRIX = 8   #IRIX
    const ELFOSABI_FREEBSD = 9   #FreeBSD
    const ELFOSABI_TRU64 = 10  #Compaq TRU64 UNIX
    const ELFOSABI_MODESTO = 11  #Novell Modesto
    const ELFOSABI_OPENBSD = 12  #Open BSD
    const ELFOSABI_OPENVMS = 13  #Open VMS
    const ELFOSABI_NSK = 14  #Hewlett-Packard Non-Stop Kernel
    const ELFOSABI_AROS = 15  #Amiga Research OS
end

# Special Section Indices
@constants SHN_INDICES "" begin
    const SHN_UNDEF = 0
    const SHN_LORESERVE = 0xff00
    const SHN_LOPROC = 0xff00
    const SHN_HIPROC = 0xff1f
    const SHN_LOOS = 0xff20
    const SHN_HIOS = 0xff3f
    const SHN_ABS = 0xfff1
    const SHN_COMMON = 0xfff2
    const SHN_XINDEX = 0xffff
    const SHN_HIRESERVE = 0xffff
end

# sh_type constants
@constants SHT_TYPES "" begin
    const SHT_NULL = 0
    const SHT_PROGBITS = 1
    const SHT_SYMTAB = 2
    const SHT_STRTAB = 3
    const SHT_RELA = 4
    const SHT_HASH = 5
    const SHT_DYNAMIC = 6
    const SHT_NOTE = 7
    const SHT_NOBITS = 8
    const SHT_REL = 9
    const SHT_SHLIB = 10
    const SHT_DYNSYM = 11
    const SHT_INIT_ARRAY = 14
    const SHT_FINI_ARRAY = 15
    const SHT_PREINIT_ARRAY = 16
    const SHT_GROUP = 17
    const SHT_SYMTAB_SHNDX = 18
    const SHT_GNU_HASH = 0x6ffffff6
    const SHT_GNU_LIBLIST = 0x6ffffff7
    const SHT_CHECKSUM = 0x6ffffff8
    const SHT_SUNW_move	= 0x6ffffffa
    const SHT_SUNW_COMDAT = 0x6ffffffb
    const SHT_SUNW_syminfo = 0x6ffffffc
    const SHT_GNU_verdef = 0x6ffffffd
    const SHT_GNU_verneed = 0x6ffffffe
    const SHT_GNU_versym = 0x6fffffff
    const SHT_X86_64_UNWIND = 0x70000001
end

@constants SHT_MASKS "" begin
    const SHT_LOOS = 0x60000000
    const SHT_HIOS = 0x6fffffff
    const SHT_LOPROC = 0x70000000
    const SHT_HIPROC = 0x7fffffff
    const SHT_LOUSER = 0x80000000
    const SHT_HIUSER = 0xffffffff
end

# Section Attribute Flags
@constants SHF_FLAGS "" begin
    const SHF_WRITE = 0x1
    const SHF_ALLOC = 0x2
    const SHF_EXECINSTR = 0x4
    const SHF_MERGE = 0x10
    const SHF_STRINGS = 0x20
    const SHF_INFO_LINK = 0x40
    const SHF_LINK_ORDER = 0x80
    const SHF_OS_NONCONFORMING = 0x100
    const SHF_GROUP = 0x200
    const SHF_TLS = 0x400
end

@constants SHF_MASKS "" begin
    const SHF_MASKOS = 0x0f000000
    const SHF_MASKPROC = 0xf0000000
end

# Section Group Flags
const GRP_COMDAT = 0x1
const GRP_MASKOS = 0x0ff00000
const GRP_MASKPROC = 0xf0000000

# Symbol Binding
@constants STB_BINDING "" begin
    const STB_LOCAL = 0
    const STB_GLOBAL = 1
    const STB_WEAK = 2
    const STB_LOOS = 10
    const STB_HIOS = 12
    const STB_LOPROC = 13
    const STB_HIPROC = 15
end

# Symbol Types
@constants STT_TYPES "" begin
    const STT_NOTYPE = 0
    const STT_OBJECT = 1
    const STT_FUNC = 2
    const STT_SECTION = 3
    const STT_FILE = 4
    const STT_COMMON = 5
    const STT_TLS = 6
    const STT_LOOS = 10
    const STT_HIOS = 12
    const STT_LOPROC = 13
    const STT_HIPROC = 15
end

# Symbol Visibility
@constants STV_VISIBILITY "" begin
    const STV_DEFAULT = 0
    const STV_INTERNAL = 1
    const STV_HIDDEN = 2
    const STV_PROTECTED = 3
end

# ptype constants
@constants P_TYPE "" begin
    const PT_NULL = 0
    const PT_LOAD = 1
    const PT_DYNAMIC = 2
    const PT_INTERP = 3
    const PT_NOTE = 4
    const PT_SHLIB = 5
    const PT_PHDR = 6
    const PT_TLS = 7
    const PT_LOOS = 0x60000000
    const PT_GNU_EH_FRAME = 0x6474e550
    const PT_GNU_STACK = 0x6474e551
    const PT_GNU_RELRO = 0x6474e552
    const PT_GNU_PROPERTY = 0x6474e553	# GNU property
    const PT_PAX_FLAGS = 0x65041580
    const PT_HIOS = 0x6fffffff
    const PT_LOPROC = 0x70000000
    const PT_HIPROC = 0x7fffffff
end

# p_flags constants
@constants PF "" begin
    const PF_X = 0x1 #Execute
    const PF_W = 0x2 #Write
    const PF_R = 0x4 #Read
    const PF_MASKOS = 0x0ff00000  #Unspecified
    const PF_MASKPROC = 0xf0000000  #Unspecified
end

# X86_64 relocations
@constants R_X86_64 "" begin
    const R_X86_64_NONE = 0
    const R_X86_64_64 = 1
    const R_X86_64_PC32 = 2
    const R_X86_64_GOT32 = 3
    const R_X86_64_PLT32 = 4
    const R_X86_64_COPY = 5
    const R_X86_64_GLOB_DAT = 6
    const R_X86_64_JUMP_SLOT = 7
    const R_X86_64_RELATIVE = 8
    const R_X86_64_GOTPCREL = 9
    const R_X86_64_32 = 10
    const R_X86_64_32S = 11
    const R_X86_64_16 = 12
    const R_X86_64_PC16 = 13
    const R_X86_64_8 = 14
    const R_X86_64_PC8 = 15
    const R_X86_64_DTPMOD64 = 16
    const R_X86_64_DTPOFF64 = 17
    const R_X86_64_TPOFF64 = 18
    const R_X86_64_TLSGD = 19
    const R_X86_64_TLSLD = 20
    const R_X86_64_DTPOFF32 = 21
    const R_X86_64_GOTTPOFF = 22
    const R_X86_64_TPOFF32 = 23
    const R_X86_64_PC64 = 24
    const R_X86_64_GOTOFF64 = 25
    const R_X86_64_GOTPC32 = 26
    const R_X86_64_GOT64 = 27
    const R_X86_64_GOTPCREL64 = 28
    const R_X86_64_GOTPC64 = 29
    const R_X86_64_GOTPLT64 = 30
    const R_X86_64_PLTOFF64 = 31
    const R_X86_64_SIZE32 = 32
    const R_X86_64_SIZE64 = 33
    const R_X86_64_GOTPC32_TLSDESC = 34
    const R_X86_64_TLSDESC_CALL = 35
    const R_X86_64_TLSDESC = 36
    const R_X86_64_IRELATIVE = 37
end

# PPC64 relocations
@constants R_PPC64 "" begin
    const R_PPC64_NONE = 0
    const R_PPC64_ADDR32 = 1
    const R_PPC64_ADDR24 = 2
    const R_PPC64_ADDR16 = 3
    const R_PPC64_ADDR16_LO = 4
    const R_PPC64_ADDR16_HI = 5
    const R_PPC64_ADDR16_HA = 6
    const R_PPC64_ADDR14 = 7
    const R_PPC64_ADDR14_BRTAKEN = 8
    const R_PPC64_ADDR14_BRNTAKEN = 9
    const R_PPC64_REL24 = 10
    const R_PPC64_REL14 = 11
    const R_PPC64_REL14_BRTAKEN = 12
    const R_PPC64_REL14_BRNTAKEN = 13
    const R_PPC64_GOT16 = 14
    const R_PPC64_GOT16_LO = 15
    const R_PPC64_GOT16_HI = 16
    const R_PPC64_GOT16_HA = 17
    const R_PPC64_GLOB_DAT = 20
    const R_PPC64_JMP_SLOT = 21
    const R_PPC64_RELATIVE = 22
    const R_PPC64_REL32 = 26
    const R_PPC64_ADDR64 = 38
    const R_PPC64_ADDR16_HIGHER = 39
    const R_PPC64_ADDR16_HIGHERA = 40
    const R_PPC64_ADDR16_HIGHEST = 41
    const R_PPC64_ADDR16_HIGHESTA = 42
    const R_PPC64_REL64 = 44
    const R_PPC64_TOC16 = 47
    const R_PPC64_TOC16_LO = 48
    const R_PPC64_TOC16_HI = 49
    const R_PPC64_TOC16_HA = 50
    const R_PPC64_TOC = 51
    const R_PPC64_ADDR16_DS = 56
    const R_PPC64_ADDR16_LO_DS = 57
    const R_PPC64_GOT16_DS = 58
    const R_PPC64_GOT16_LO_DS = 59
    const R_PPC64_TOC16_DS = 63
    const R_PPC64_TOC16_LO_DS = 64
    const R_PPC64_TLS = 67
    const R_PPC64_DTPMOD64 = 68
    const R_PPC64_TPREL16 = 69
    const R_PPC64_TPREL16_LO = 70
    const R_PPC64_TPREL16_HI = 71
    const R_PPC64_TPREL16_HA = 72
    const R_PPC64_TPREL64 = 73
    const R_PPC64_DTPREL16 = 74
    const R_PPC64_DTPREL16_LO = 75
    const R_PPC64_DTPREL16_HI = 76
    const R_PPC64_DTPREL16_HA = 77
    const R_PPC64_DTPREL64 = 78
    const R_PPC64_GOT_TLSGD16 = 79
    const R_PPC64_GOT_TLSGD16_LO = 80
    const R_PPC64_GOT_TLSGD16_HI = 81
    const R_PPC64_GOT_TLSGD16_HA = 82
    const R_PPC64_GOT_TLSLD16 = 83
    const R_PPC64_GOT_TLSLD16_LO = 84
    const R_PPC64_GOT_TLSLD16_HI = 85
    const R_PPC64_GOT_TLSLD16_HA = 86
    const R_PPC64_GOT_TPREL16_DS = 87
    const R_PPC64_GOT_TPREL16_LO_DS = 88
    const R_PPC64_GOT_TPREL16_HI = 89
    const R_PPC64_GOT_TPREL16_HA = 90
    const R_PPC64_GOT_DTPREL16_DS = 91
    const R_PPC64_GOT_DTPREL16_LO_DS = 92
    const R_PPC64_GOT_DTPREL16_HI = 93
    const R_PPC64_GOT_DTPREL16_HA = 94
    const R_PPC64_TPREL16_DS = 95
    const R_PPC64_TPREL16_LO_DS = 96
    const R_PPC64_TPREL16_HIGHER = 97
    const R_PPC64_TPREL16_HIGHERA = 98
    const R_PPC64_TPREL16_HIGHEST = 99
    const R_PPC64_TPREL16_HIGHESTA = 100
    const R_PPC64_DTPREL16_DS = 101
    const R_PPC64_DTPREL16_LO_DS = 102
    const R_PPC64_DTPREL16_HIGHER = 103
    const R_PPC64_DTPREL16_HIGHERA = 104
    const R_PPC64_DTPREL16_HIGHEST = 105
    const R_PPC64_DTPREL16_HIGHESTA = 106
    const R_PPC64_TLSGD = 107
    const R_PPC64_TLSLD = 108
    const R_PPC64_REL16 = 249
    const R_PPC64_REL16_LO = 250
    const R_PPC64_REL16_HI = 251
    const R_PPC64_REL16_HA = 252
end

@constants AUXV_TYPE "AT_" begin
    const AT_NULL = 0
    const AT_IGNORE = 1
    const AT_EXECFD = 2
    const AT_PHDR = 3
    const AT_PHENT = 4
    const AT_PHNUM = 5
    const AT_PAGESZ = 6
    const AT_BASE = 7
    const AT_FLAGS = 8
    const AT_ENTRY = 9
    const AT_NOTELF = 10
    const AT_UID = 11
    const AT_EUID = 12
    const AT_GID = 13
    const AT_EGID = 14
    const AT_PLATFORM = 15
    const AT_HWCAP = 16
    const AT_CLKTCK = 17
    const AT_SYSINFO = 32
    const AT_SYSINFO_EHDR = 33
end

@constants DYNAMIC_TYPE "" begin
    const DT_NULL = 0
    const DT_NEEDED = 1
    const DT_PLTRELSZ = 2
    const DT_PLTGOT = 3
    const DT_HASH = 4
    const DT_STRTAB = 5
    const DT_SYMTAB = 6
    const DT_RELA = 7
    const DT_RELASZ = 8
    const DT_RELAENT = 9
    const DT_STRSZ = 10
    const DT_SYMENT = 11
    const DT_INIT = 12
    const DT_FINI = 13
    const DT_SONAME = 14
    const DT_RPATH = 15
    const DT_SYMBOLIC = 16
    const DT_REL = 17
    const DT_RELSZ = 18
    const DT_RELENT = 19
    const DT_PLTREL = 20
    const DT_DEBUG = 21
    const DT_TEXTREL = 22
    const DT_JMPREL = 23
    const DT_BIND_NOW = 24
    const DT_INIT_ARRAY = 25
    const DT_FINI_ARRAY = 26
    const DT_INIT_ARRAYSZ = 27
    const DT_FINI_ARRAYSZ = 28
    const DT_RUNPATH = 29
    const DT_FLAGS = 30
    const DT_ENCODING = 32
    const OLD_DT_LOOS = 0x60000000
    const DT_LOOS = 0x6000000d
    const DT_HIOS = 0x6ffff000
    const DT_VALRNGLO = 0x6ffffd00
    const DT_VALRNGHI = 0x6ffffdff
    const DT_ADDRRNGLO = 0x6ffffe00
    const DT_ADDRRNGHI = 0x6ffffeff
    const DT_VERSYM = 0x6ffffff0
    const DT_RELACOUNT = 0x6ffffff9
    const DT_RELCOUNT = 0x6ffffffa
    const DT_FLAGS_1 = 0x6ffffffb
    const DT_VERDEF = 0x6ffffffc
    const DT_VERDEFNUM = 0x6ffffffd
    const DT_VERNEED = 0x6ffffffe
    const DT_VERNEEDNUM = 0x6fffffff
    const OLD_DT_HIOS = 0x6fffffff
    const DT_LOPROC = 0x70000000
    const DT_HIPROC = 0x7fffffff
    const DT_VALRNGLO = 0x6ffffd00
    const DT_GNU_PRELINKED = 0x6ffffdf5# Prelinking timestamp
    const DT_GNU_CONFLICTSZ = 0x6ffffdf6# Size of conflict section */
    const DT_GNU_LIBLISTSZ = 0x6ffffdf7# Size of library list */
    const DT_CHECKSUM = 0x6ffffdf8
    const DT_PLTPADSZ = 0x6ffffdf9
    const DT_MOVEENT = 0x6ffffdfa
    const DT_MOVESZ = 0x6ffffdfb
    const DT_FEATURE_1 = 0x6ffffdfc# Feature selection (DTF_*).  */
    const DT_POSFLAG_1 = 0x6ffffdfd# Flags for DT_* entries, effecting
    const DT_SYMINSZ = 0x6ffffdfe# Size of syminfo table (in bytes) */
    const DT_SYMINENT = 0x6ffffdff# Entry size of syminfo */
    const DT_VALRNGHI = 0x6ffffdff
    const DT_VALNUM = 12
    const DT_ADDRRNGLO = 0x6ffffe00
    const DT_GNU_HASH = 0x6ffffef5# GNU-style hash table.  */
    const DT_TLSDESC_PLT = 0x6ffffef6
    const DT_TLSDESC_GOT = 0x6ffffef7
    const DT_GNU_CONFLICT = 0x6ffffef8# Start of conflict section */
    const DT_GNU_LIBLIST = 0x6ffffef9# Library list */
    const DT_CONFIG = 0x6ffffefa# Configuration information.  */
    const DT_DEPAUDIT = 0x6ffffefb# Dependency auditing.  */
    const DT_AUDIT = 0x6ffffefc# Object auditing.  */
    const DT_PLTPAD = 0x6ffffefd# PLT padding.  */
    const DT_MOVETAB = 0x6ffffefe# Move table.  */
    const DT_SYMINFO = 0x6ffffeff# Syminfo table.  */
    const DT_ADDRRNGHI = 0x6ffffeff
    const DT_VERSYM = 0x6ffffff0
    const DT_RELACOUNT = 0x6ffffff9
    const DT_RELCOUNT = 0x6ffffffa
    const DT_FLAGS_1 = 0x6ffffffb# State flags, see DF_1_* below.  */
    const DT_VERDEF = 0x6ffffffc# Address of version definition
    const DT_VERDEFNUM = 0x6ffffffd# Number of version definitions */
    const DT_VERNEED = 0x6ffffffe# Address of table with needed
    const DT_VERNEEDNUM = 0x6fffffff# Number of needed versions */
end

@constants NOTE_TYPES "" begin
    const NT_GNU_ABI_TAG = 1
    const NT_GNU_HWCAP = 2
    const NT_GNU_BUILD_ID = 3
    const NT_GNU_GOLD_VERSION = 4
    const NT_GNU_PROPERTY_TYPE_0 = 5
end

@constants GNU_PROPERTY "" begin
    const GNU_PROPERTY_STACK_SIZE = UInt32(1)
    const GNU_PROPERTY_NO_COPY_ON_PROTECTED	= UInt32(2)
    const GNU_PROPERTY_LOPROC = 0xc0000000
    const GNU_PROPERTY_HIPROC = 0xdfffffff
    const GNU_PROPERTY_LOUSER = 0xe0000000
    const GNU_PROPERTY_HIUSER = 0xffffffff
    const GNU_PROPERTY_X86_ISA_1_USED = 0xc0010002
    const GNU_PROPERTY_X86_ISA_1_NEEDED = 0xc0008002
    const GNU_PROPERTY_X86_FEATURE_1_AND = 0xc0000002
end

@constants GNU_PROPERTY_X86_ISA_1 "GNU_PROPERTY_X86_ISA_1_" begin
    const GNU_PROPERTY_X86_ISA_1_BASELINE = UInt32(1)
    const GNU_PROPERTY_X86_ISA_1_V2 = UInt32(2)
    const GNU_PROPERTY_X86_ISA_1_V3 = UInt32(4)
    const GNU_PROPERTY_X86_ISA_1_V4	= UInt32(8)
end

const GNU_PROPERTY_X86_ISA_1_to_str = Dict(GNU_PROPERTY_X86_ISA_1_BASELINE=>"CMOV, CX8 (cmpxchg8b), FPU (fld), MMX, OSFXSR (fxsave), SCE (syscall), SSE, SSE2", GNU_PROPERTY_X86_ISA_1_V2 => "CMPXCHG16B (cmpxchg16b), LAHF-SAHF (lahf), POPCNT (popcnt), SSE3, SSSE3, SSE4.1, SSE4.2", GNU_PROPERTY_X86_ISA_1_V3 => "AVX, AVX2, BMI1, BMI2, F16C, FMA, LZCNT, MOVBE, XSAVE", GNU_PROPERTY_X86_ISA_1_V4 => "AVX512F, AVX512BW, AVX512CD, AVX512DQ and AVX512VL")

@constants GNU_PROPERTY_X86_FEATURE_1 "GNU_PROPERTY_X86_FEATURE_1_" begin
    const GNU_PROPERTY_X86_FEATURE_1_IBT = UInt32(1)
    const GNU_PROPERTY_X86_FEATURE_1_SHSTK = UInt32(2)
end

@constants ELF_NOTE_OS "ELF_NOTE_OS_" begin
    const ELF_NOTE_OS_LINUX =	0
    const ELF_NOTE_OS_GNU = 1
    const ELF_NOTE_OS_SOLARIS2 = 2
    const ELF_NOTE_OS_FREEBSD = 3
end

@constants VER_NDX "VER_NDX_" begin
	const VER_NDX_LOCAL = 0x0000
	const VER_NDX_GLOBAL = 0x0001
	const VER_NDX_ELIMINATE= 0xff01
end

@constants VER_DEF "VER_DEF_" begin
    const VER_DEF_NONE = 0
    const VER_DEF_CURRENT = 1
    const VER_DEF_NUM = 2
end

@constants VER_FLG "VER_FLG_" begin
    const VER_FLG_BASE = 0x1
    const VER_FLG_WEAK = 0x2
end

@constants VER_NEED "VER_NEED_" begin
    const VER_NEED_NONE = 0
    const VER_NEED_CURRENT = 1
    const VER_NEED_NUM = 2
end


text:
module DestructuringTests
using  ExpressionPatterns.Destructuring
using  ExpressionPatterns.Dispatch
using  Base.Test

macro anyvec(x) anyvec(x) end

@metafunction anyvec(x) = x
@metafunction anyvec([*{elts}]) = :(Any[$(map(anyvec, elts)...)])

## letds

# No destructuring.
@letds x = 10 y = 20 @test (x,y) == (10,20)

# Basic destructuring.
@letds (x,y) = :(x,y) @test (x,y) == (:x,:y)

# Destructuring with nested expressions.
@letds f((a+b),c{d}) = :(func(10+20,A{T})) begin
  @test (f,a,b,c,d) == (:func,10,20,:A,:T)
end

# Destructuring with slurps.
@letds (for x in [*{elts}]; *{body}  end) =
      :(for i in [1,2,3,4]; print(i) end) begin

  @test (x,elts,body) == (:i, Any[1,2,3,4], Any[:(print(i))])
end

# Destructuring with alterning slurps.
@letds [*{a,b}] = :[1,2,3,4] @test (a,b) == ([1,3],[2,4])

# Destructuring with multiple slurps
@letds [*{before3}, 3, *{after3}] = :[1,2,3,4,5] begin
  @test (before3,after3) == ([1,2],[4,5])
end

# Destructuring with nested slurps.
@letds (begin *{f(*{elts})}  end) =
      :(begin f(1,2); g(2,3,4) end) begin

  @test f    == [:f, :g]
  @test elts == @anyvec [[1,2],[2,3,4]]
end

# Triple nested slurps! this is fun!
@letds (begin *{[*{f(*{(*{elts},)})}]} end) =
      :(begin [f((x,y),(a,b)), g((1,2),(4,5))]
      	      [h((1,2,3,4))] end) begin

  @test f    == @anyvec [[:f, :g], [:h]]
  @test elts == @anyvec [[[[:x,:y],[:a,:b]],[[1,2],[4,5]]],[[[1,2,3,4]]]]
end

## macrods

@macrods nodestruct(x,y,z) (x,y,z)
@test   @nodestruct(a,10+20,f(x,y)) == (:a,:(10+20), :(f(x,y)))

@macrods basicdestruct(x+y) (x,y)
@test   @basicdestruct(1+2) == (1,2)

@macrods nestedexpr(f(x::T1, y::T2) = body) (f,x,y,T1,T2,body)
@test   @nestedexpr(add(a::Int, b::Int) = a+b) == (:add, :a,:b,:Int,:Int,:(a+b))

@macrods withslurp(for x in  coll;     *{body}   end)    :(map($coll)   do $x; $(body...) end)
@test   @withslurp(for i in [1,2,3]; i*=2; i+10; end) ==  [12, 14, 16]

#alternating slurps
@macrods switch(x, begin *{case(c),expr} end) (x,c,expr)
ex=
@switch i begin
  case(1); print(1);
  case(2); println(2);
end
@test ex == (:i, [1,2], [:(print(1)), :(println(2))])

@macrods nestedslurp([*{(*{tuple_elts},), [*{vector_elts}]}]) (tuple_elts, vector_elts)
@test   @nestedslurp([(1,2,3), [2,3,4], (3,2,5), [7,4,3]]) ==
        (Any[[1,2,3],[3,2,5]], Any[[2,3,4], [7,4,3]])

## anonds | funds

@test @anonds((x,y) -> x+y)(1,2) == 3

@test @anonds(f(x,y) -> :($f($x)($y)))(:(g(1,2))) == :(g(1)(2))

@test @anonds(*{a,b,c} -> (a,b,c))(1,2,3,4,5,6) == ([1,4],[2,5],[3,6])

@funds olddict((K=>V)(*{keys=>values},)) = :(Dict{$K,$V}($([:($key=>$value) for (key,value) in zip(keys,values)]...)))
@test  olddict(:((Symbol=>Int)(:a=>1, :b=>2, :c=>3))) == :(Dict{Symbol,Int}(:a=>1, :b=>2, :c=>3))

end

ext: jl


text:
module aoc_01
    using Test
    export read_input, part1, part2
    export test_part1, test_part2, test
    export input, demo

    const input_r = r""

    function read_input(filename = "input.txt")
        lines = readlines(filename)
        # for line in eachline(filename)
        for line in lines
        end
        # [parse_line(line) for line in eachline(filename)]
        lines
    end

    function parse_line(line)
    end

    const input = read_input("input.txt")
    # const demo = read_input("demo.txt")

    function part1(input = input)
    end

    function test_part1()
        @testset "Part 1" begin
            @test true
        end
        part1( demo )
    end

    function part2(input = input)
    end

    function test_part2()
        @testset "Part 2" begin
            @test true
        end
        part2( demo )
    end

    function test()
        test_part1()
        test_part2()
    end
end


text:
# All implementations of data loaders/creaters are in this file or loaded from here

export SyntheticGlacier, ITMIXGlacier, RGIGlacier

# generic loader functions
include("generic.jl")

"""
Synthetic glaciers:
- Bench
"""
struct SyntheticGlacier <: GlacierID
    id::Symbol
end
getname(gid::SyntheticGlacier) = string(gid.id)
getregion(::SyntheticGlacier) = Region.Synthetic

const SyntheticBenchLoader = :SynthBench
include("synthetic-bench.jl")

"""
If there is nothing to load, return `SomeData(nothing)`.

Probably only makes sense for iv and thickness, thus only those are implemented.
"""
const NothingLoader = :Nothing
load!(::DataSet{IVData, NothingLoader}, args...) = SomeData(nothing), Matrix{Bool}(0,0)
load!(::DataSet{ThicknessData, NothingLoader}, args...) = SomeData(nothing)

"All ITMIX glaciers ever defined.  Check to code for which versions"
const ITMIXnames = [:Academy, :Aqqutikitsoq, :Austfonna, :AustreGroenfjordbreen, # v2
                    :Brewster, :ChhotaShigri, # v2
                    :Columbia, :Devon, :Elbrus, :Freya, :Hellstugubreen, :Kesselwandferner,
                    :Mocho, :NorthGlacier, :SouthGlacier, :Starbuck, :Synthetic1, :Synthetic2,
                    :Synthetic3, :Tasman, :Unteraar, :Urumqi, :Washmawapta]
const itmix2only = [:AustreGroenfjordbreen, :ChhotaShigri]

"""
ITMIX glacier data
"""
struct ITMIXGlacier <: GlacierID
    id::Symbol # the name as symbol
    version::Int # ITMIX version
    function ITMIXGlacier(name, version)
        @assert 0<version<3
        @assert name in ITMIXnames "$name not in BITEModel.ITMIXnames"
        new(name, version)
    end
end

getname(gid::ITMIXGlacier) = string(gid.id)
getregion(gid::ITMIXGlacier) = Region._Region(gid.version+200)
function getrgi(gid::ITMIXGlacier)
    error("Function not up to date")
    r = iscalving_isicecap[gid.id][3]
    r==nothing && error("Glacier $gid has no RGI number.")
    return r
end
const ITMIXLoader = :ITMIX

ITMIXv1_glaciers = [ITMIXGlacier(n,1) for n in ITMIXnames if !(n in itmix2only)]
ITMIXv2_glaciers = [ITMIXGlacier(n,2) for n in ITMIXnames]
include("itmix.jl")

"""
As defined in the RGI 6.0 technical document.
"""
@with_kw struct RGIattrs
    # bgndate::Date
    # enddate::Date
    lon::F
    lat::F
    subregion::Int
    area::F
    zmin::F
    zmax::F
    # zmed::F
    # slope::F
    aspect::F
    # lmax::F
    status::Int # 0,1 ok, 2 not proper outline
    form::Int
    termtype::Int
    surging::Int
    name::String
end
function RGIattrs(id::Union{String,Symbol})
    reg = getrgi_region(string(id))
    nr = getrgi_nr(string(id))
    RGIattrs(read_attrib_file(reg, nr)[2:end]...)
end

"""
Randolph glacier inventory (RGI) glaciers

http://www.glims.org/RGI/index.html

Use RGI 6.0 numbers, e.g. RGI60-07.01455

Note that RGI itself only provides the outlines.  However Matthias
also prepared DEM files.
"""
struct RGIGlacier <: GlacierID
    id::Symbol
    attrs::RGIattrs
    function RGIGlacier(id::Union{String,Symbol}, attrs::RGIattrs)
        @assert getrgi_version(string(id))==v"6.0" "Need RGI version 6.0 IDs"
        new(uppercase(string(id)), attrs)
    end
    RGIGlacier(id::Union{String,Symbol}) = RGIGlacier(id, RGIattrs(id))
end
function RGIGlacier(nr::Int, region::Int)
    r = VAWTools.int2str2(region)
    id = VAWTools.int2str5(nr)
    RGIGlacier("RGI60-$r.$id")
end
# getid uses default
getrgi(gid::RGIGlacier) = string(getid(gid)) # same as ID in this case, but string
getregion(gid::RGIGlacier) = Region._Region(getrgi_region(gid))
"RGILoader"
const RGILoader = :RGI
getname(gid::RGIGlacier)::String = gid.attrs.name
Base.show(stream::IO, gid::RGIGlacier) = print(stream, "RGIGlacier($(gid.id),...)")

# RGI specifics:
function parse_rgi(rgi::AbstractString)
    vers = VersionNumber(parse(Int,split(rgi,'-')[1][4:end])÷10)
    reg = parse(Int,split(split(rgi,'-')[2],'.')[1])
    nr = parse(Int,split(split(rgi,'-')[2],'.')[2])
    @assert 0<reg<20 "Not a valid RGI region!"
    return vers, reg, nr
end
parse_rgi(gid::RGIGlacier) = parse_rgi(string(gid.id))
getrgi_version(rgi) = parse_rgi(rgi)[1]
getrgi_region(rgi) = parse_rgi(rgi)[2]
getrgi_nr(rgi) = parse_rgi(rgi)[3]

include("rgi.jl")

"""
Antarctic Peninsula glaciers with ID as Huss&Farinotti uses.
"""
struct PeninsulaGlacier <: GlacierID
    id::Int
end
getname(gid::PeninsulaGlacier) = string(gid.id)
getregion(::PeninsulaGlacier) = Region.AntarcticPeninsula
"Data loader from H&F Peninsula"
const PeninsulaLoader = :Peninsula
include("peninsula.jl")

"""
GloGEM bdot data-loader

Glacier ID according to RGI 6.0
"""
const GloGEMLoader = :GloGEM
include("glogem.jl")

"""
GlaThiDa thickness data-loader

Glacier ID according to RGI 6.0
"""
const GlaThiDaLoader = :GlaThiDa
include("glathida.jl")

"""
Misc CCI loaders (IV and dh/dt)

"""
const CCILoader = :CCI
include("CCI.jl")

"""
Glamos data-loader

(No glaciers associated with it)
"""
const GlamosLoader = :GLAMOS

"""
No data loader

Always returns SomeData(nothing)
"""
const VoidLoader = :VOID
load!(ds::DataSet{<:Any, VoidLoader}, varargs...) = SomeData(nothing)


text:
using Oceananigans
using Plots
using JLD2

include("Analysis.jl")

using .Analysis

filenames = [
             "oceananigans_unstable_bickley_jet_Nh32_WENO5_fields.jld2",
             "oceananigans_unstable_bickley_jet_Nh64_WENO5_fields.jld2",
             "oceananigans_unstable_bickley_jet_Nh128_WENO5_fields.jld2",
             "oceananigans_unstable_bickley_jet_Nh256_WENO5_fields.jld2",
             "oceananigans_unstable_bickley_jet_Nh512_WENO5_fields.jld2",
             "oceananigans_unstable_bickley_jet_Nh1024_WENO5_fields.jld2",
            # "oceananigans_unstable_bickley_jet_Nh2048_WENO5_fields.jld2"
            ]

resolutions = [32, 64, 128, 256, 512, 1024, 2048]

results = [Analysis.oceananigans_statistics(filename) for filename in filenames]

 t = map(r -> r[1], results)
Δe = map(r -> r[2], results)
Δχ = map(r -> r[3], results)
c★ = map(r -> r[4], results)

kwargs = Dict(
              :xlabel => "Time",
             )

energy_plot   = Analysis.plot_resolutions(t, Δe, resolutions;
                                          ylabel = "Kinetic energy",
                                          legend = :bottomleft,
                                          kwargs...)

max_c_plot    = Analysis.plot_resolutions(t, c★, resolutions;
                                          ylabel = "\$ \\max |c| \$",
                                          legend = :bottomleft,
                                          kwargs...)

variance_plot = Analysis.plot_resolutions(t, Δχ, resolutions;
                                          ylabel = "Tracer variance",
                                          legend = :bottomleft,
                                          kwargs...)

display(plot(energy_plot, variance_plot, max_c_plot, layout=(3, 1), size=(800, 600)))

filename = filenames[end]

file = jldopen(filename)
iterations = parse.(Int, keys(file["timeseries/t"]))

c0 = Analysis.CellFileField(file, :c, 0)
c1 = Analysis.CellFileField(file, :c, iterations[end])

close(file)

h_plot = histogram(interior(c0)[:], bins=20, label="\$ t = 0 \$")
histogram!(h_plot, interior(c1)[:], bins=20, label="\$ t = 200 \$")

display(h_plot)


text:
export iauFk52h
"""
Transform FK5 (J2000.0) star data into the Hipparcos system.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given (all FK5, equinox J2000.0, epoch J2000.0):
   r5      double    RA (radians)
   d5      double    Dec (radians)
   dr5     double    proper motion in RA (dRA/dt, rad/Jyear)
   dd5     double    proper motion in Dec (dDec/dt, rad/Jyear)
   px5     double    parallax (arcsec)
   rv5     double    radial velocity (km/s, positive = receding)

Returned (all Hipparcos, epoch J2000.0):
   rh      double    RA (radians)
   dh      double    Dec (radians)
   drh     double    proper motion in RA (dRA/dt, rad/Jyear)
   ddh     double    proper motion in Dec (dDec/dt, rad/Jyear)
   pxh     double    parallax (arcsec)
   rvh     double    radial velocity (km/s, positive = receding)

Notes:

   1> This function transforms FK5 star positions and proper motions
   into the system of the Hipparcos catalog.

   2. The proper motions in RA are dRA/dt rather than
   cos(Dec)*dRA/dt, and are per year rather than per century.

   3. The FK5 to Hipparcos transformation is modeled as a pure
   rotation and spin;  zonal errors in the FK5 catalog are not
   taken into account.

   4. See also iauH2fk5, iauFk5hz, iauHfk5z.

Called:
   iauStarpv    star catalog data to space motion pv-vector
   iauFk5hip    FK5 to Hipparcos rotation and spin
   iauRxp       product of r-matrix and p-vector
   iauPxp       vector product of two p-vectors
   iauPpp       p-vector plus p-vector
   iauPvstar    space motion pv-vector to star catalog data

Reference:

   F.Mignard & M.Froeschle, Astron.Astrophys., 354, 732-739 (2000).

This revision:  2017 October 12

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauFk52h(double r5, double d5,
#               double dr5, double dd5, double px5, double rv5,
#               double *rh, double *dh,
#               double *drh, double *ddh, double *pxh, double *rvh)

function iauFk52h(r5::Real, d5::Real, dr5::Real, dd5::Real, px5::Real, rv5::Real)
   ref_rh  = Ref{Float64}(0.0)
   ref_dh  = Ref{Float64}(0.0)
   ref_drh = Ref{Float64}(0.0)
   ref_ddh = Ref{Float64}(0.0)
   ref_pxh = Ref{Float64}(0.0)
   ref_rvh = Ref{Float64}(0.0)

   ccall((:iauFk52h, libsofa_c), Cdouble,
         (Cdouble, Cdouble, Cdouble, Cdouble,
         Cdouble, Cdouble,
         Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},
         Ref{Cdouble}, Ref{Cdouble}, ),
         convert(Float64, r5), convert(Float64, d5),
         convert(Float64, dr5), convert(Float64, dd5),
         convert(Float64, px5), convert(Float64, rv5),
         ref_rh, ref_dh, ref_drh, ref_ddh, ref_pxh, ref_rvh)

   return ref_rh[], ref_dh[], ref_drh[], ref_ddh[], ref_pxh[], ref_rvh[]
end


text:
function rand!(low::Real, high::Real, out::NDArray)
  _random_uniform(low, high, out)
end
function rand{N}(low::Real, high::Real, shape::NTuple{N, Int})
  rand(low, high, shape, cpu())
end
function rand{N}(low::Real, high::Real, shape::NTuple{N, Int}, ctx::Context)
  out = empty(shape, ctx)
  rand!(low, high, out)
end

function randn!(mean::Real, stdvar::Real, out::NDArray)
  _random_gaussian(mean, stdvar, out)
end
function randn{N}(mean::Real, stdvar::Real, shape::NTuple{N,Int})
  randn(mean, stdvar, shape, cpu())
end
function randn{N}(mean::Real, stdvar::Real, shape::NTuple{N,Int}, ctx::Context)
  out = empty(shape, ctx)
  randn!(mean, stdvar, out)
end

function srand!(seed_state::Int)
  @mxcall(:MXRandomSeed, (Cint,), seed_state)
end


text:
function finde(A,n)
    Antwort = zeros(size(A,1)*size(A,2),2);
    counter = 1;
    for i = 1:size(A,1)
        for j = 1:size(A,2)
            if A[i,j] == n
                Antwort[counter,1] = i;
                Antwort[counter,2] = j;
                counter = counter+1;
            end
        end
    end
    X = Antwort[:,1];
    Y = Antwort[:,2];
    #display(X)
    #display(Y)
    #display(counter)
    Antwort1 = resize(Antwort[:,1],counter-1);
    Antwort2 = resize(Antwort[:,2],counter-1);
    Antwortneu = hcat(Antwort1,Antwort2);
return Antwortneu;

end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

# Implementation(s) provided in separate file, loaded via require:
function _h5io_keys end
function _h5io_objtype end
function _h5io_open end

function _h5io_add_path_to_dest(dest_with_subpath::Tuple{Any,AbstractString}, path::AbstractString)
    dest, old_path = dest_with_subpath
    @assert !isempty(old_path)
    @assert !isempty(path)
    (dest, "$old_path/$path")
end

_h5io_add_path_to_dest(dest, path::AbstractString) = (dest, path)

function _h5io_write!(dest_with_subpath::Tuple{Any,AbstractString}, data::AbstractArray{<:Real})
    dest, path = dest_with_subpath
    dest[path] = data
    nothing
end

function _h5io_write!(dest_with_subpath::Tuple{Any,AbstractString}, data::AbstractArray{<:Nothing})
    nothing
end

function _h5io_write!(dest_with_subpath::Tuple{Any,AbstractString}, data::VectorOfSimilarVectors)
    _h5io_write!(dest_with_subpath, Array(flatview(data)))
end

function _h5io_write!(dest, data::NamedTuple{names}) where {names}
    for name in names
        _h5io_write!(_h5io_add_path_to_dest(dest, String(name)), Base.getproperty(data, name))
    end
    nothing
end

function _h5io_write!(dest, data)
    Tables.istable(data) || throw(ArgumentError("data is not a table."))
    cols = Tables.columns(data)::NamedTuple
    @assert cols isa NamedTuple
    _h5io_write!(dest, cols)
end

function _h5io_read end

function _h5io_read(src_with_subpath::Tuple{Any,AbstractString})
    src, subpath = src_with_subpath
    _h5io_read(src[subpath])
end

_h5io_read(src) = _h5io_read(src, _h5io_objtype(src))

_h5io_read(src, ::Val{:dataset}) = _h5io_read_postprocess(read(src))

function _h5io_read(src, ::Val{:datafile})
    names = (sort(_h5io_keys(src))...,)
    values = map(k -> _h5io_read((src, k)), names)
    name_symbols = map(Symbol, names)
    x = NamedTuple{name_symbols}(values)
    _h5io_read_postprocess(x)
end

_h5io_read_postprocess(V::AbstractVector{<:Real}) = V

_h5io_read_postprocess(VV::AbstractMatrix{<:Real}) = VectorOfSimilarVectors(VV)

_h5io_read_postprocess(nt::NamedTuple) = TypedTables.Table(nt)

_h5io_read_postprocess(nt::NamedTuple{(:info, :logdensity, :params, :weight)}) =
    DensitySampleVector((nt.params, nt.logdensity, nt.weight, nt.info, Array{Nothing}(undef, size(nt.info)...)))

_h5io_read_postprocess(nt::NamedTuple{(:chaincycle, :chainid, :sampletype, :stepno)}) =
    MCMCSampleIDVector((nt.chainid, nt.chaincycle, nt.stepno, nt.sampletype))


text:
import InfrastructureSystems

function main(args)
    if length(args) != 2
        println("Usage: julia generate_structs.jl INPUT_FILE OUTPUT_DIRECTORY")
        exit(1)
    end

    InfrastructureSystems.generate_structs(args[1], args[2])
end

main(ARGS)


text:
function gate(::GateName"pauli_channel", dims::Tuple = (2,);
              pauli_ops = ["Id","X","Y","Z"],
              error_probabilities = prepend!(zeros(length(pauli_ops)^length(dims)-1),1))
  @assert sum(error_probabilities) ≈ 1
  @assert all(dims .== 2)
  N = length(dims)
  length(error_probabilities) > (1 << 10) && error("Hilbert space too large")
  error_probabilities ./= sum(error_probabilities)
  kraus = zeros(Complex{Float64},1<<N,1<<N,length(pauli_ops)^N)
  basis = vec(reverse.(Iterators.product(fill(pauli_ops,N)...)|>collect))
  for (k,ops) in enumerate(basis)
    kraus[:,:,k] = √error_probabilities[k] * reduce(kron,gate.(ops))
  end
  return kraus
end

gate(::GateName"bit_flip", dims::Tuple = (2,); p::Number) =
  gate("pauli_channel", dims; error_probabilities = prepend!(p/(2^length(dims)-1) * ones(2^length(dims)-1), 1-p), pauli_ops = ["Id","X"])

gate(::GateName"phase_flip", dims::Tuple = (2,); p::Number) =
  gate("pauli_channel", dims; error_probabilities = prepend!(p/(2^length(dims)-1) * ones(2^length(dims)-1), 1-p), pauli_ops = ["Id","Z"])

function gate(::GateName"AD", dims::Tuple = (2,); γ::Real)
  N = length(dims)
  @assert all(dims .== 2)
  kraus = zeros(2,2,2)
  kraus[:,:,1] = [1 0
                  0 sqrt(1-γ)]
  kraus[:,:,2] = [0 sqrt(γ)
                  0 0]
  N == 1 && return kraus
  k1 = kraus[:,:,1]
  k2 = kraus[:,:,2]
  T = vec(Iterators.product(fill([k1,k2],N)...)|>collect)
  K = zeros(Float64,1<<N,1<<N,1<<N)
  for x in 1:1<<N
    K[:,:,x] = reduce(kron,T[x])
  end
  return K
end
#
# To accept the gate name "amplitude_damping"
gate(::GateName"amplitude_damping", dims::Tuple = (2,); kwargs...) = gate("AD", dims; kwargs...)

function gate(::GateName"PD", dims::Tuple = (2,); γ::Real)
  N = length(dims)
  @assert all(dims .== 2)
  kraus = zeros(2,2,2)
  kraus[:,:,1] = [1 0
                  0 sqrt(1-γ)]
  kraus[:,:,2] = [0 0
                  0 sqrt(γ)]

  N == 1 && return kraus
  k1 = kraus[:,:,1]
  k2 = kraus[:,:,2]
  T = vec(Iterators.product(fill([k1,k2],N)...)|>collect)
  K = zeros(Float64,1<<N,1<<N,1<<N)
  for x in 1:1<<N
    K[:,:,x] = reduce(kron,T[x])
  end
  return K
end

# To accept the gate name "phase_damping"
gate(::GateName"phase_damping", dims::Tuple = (2,); kwargs...) = gate("PD", dims; kwargs...)
#
# To accept the gate name "dephasing"
gate(::GateName"dephasing", dims::Tuple = (2,); kwargs...) = gate("PD", dims; kwargs...)

# make general n-qubit
gate(::GateName"DEP", dims::Tuple = (2,); p::Number) =
  gate("pauli_channel", dims; error_probabilities = prepend!(p/(4^length(dims)-1) * ones(4^length(dims)-1), 1-p), pauli_ops = ["Id","X","Y","Z"])

# To accept the gate name "depolarizing"
gate(::GateName"depolarizing", dims::Tuple = (2,); kwargs...) =
  gate("DEP", dims; kwargs...)

function insertnoise(circuit::Vector{<:Vector{<:Any}}, noisemodel::Tuple; gate = nothing)#idlenoise::Bool = false)
  max_g_size = maxgatesize(circuit)

  # single noise model for all
  if noisemodel[1] isa String
    tmp = []
    for k in 1:max_g_size
      push!(tmp, k => noisemodel)
    end
    #if idlenoise
    #  push!(tmp, "idle" => noisemodel)
    #end
    noisemodel = Tuple(tmp)
  end

  noisycircuit = []
  for layer in circuit
    noisylayer = []
    for g in layer
      push!(noisylayer, g)
      applynoise = (isnothing(gate) ? true :
                    gate isa String ? g[1] == gate : g[1] in gate)
      #@applytogate = gate isa String ? g[1] == gate : g[1] in gate
      #if isnothing(gate) || applytogate
      if applynoise
        nq = g[2]
        # n-qubit gate
        if nq isa Tuple
          gatenoiseindex = findfirst(x -> x == length(nq), first.(noisemodel))
          isnothing(gatenoiseindex) && error("Noise model not defined for $(length(nq))-qubit gates!")
          gatenoise = last(noisemodel[gatenoiseindex])
          # Check whether the n-qubit Kraus channel has been defined
          # n -qubit Kraus operator
          noisecheck = PastaQ.gate(gatenoise[1], Tuple(repeat([2], length(nq))); gatenoise[2]...)
          # if the single-qubit copy is return, use productnoise, and throw a warning
          if size(noisecheck,1) < 1<<length(nq)
            @warn "$(length(nq))-qubit Kraus operators for the $(gatenoise[1]) noise not defined. Applying tensor-product noise instead.\n"
            # tensor-product noise
            for q in nq
              push!(noisylayer,(gatenoise[1], q, gatenoise[2]))
            end
          else
            # correlated n-qubit noise
            push!(noisylayer,(gatenoise[1], nq, gatenoise[2]))
          end
        # 1-qubit gate
        else
          gatenoiseindex = findfirst(x -> x == 1, first.(noisemodel))
          isnothing(gatenoiseindex) && error("Noise model not defined for 1-qubit gates!")
          gatenoise = last(noisemodel[gatenoiseindex])
          push!(noisylayer, (gatenoise[1], nq, gatenoise[2]))
        end
      end
    end
    #if idlenoise
    #  gatenoiseindex = findfirst(x -> x == "idle", first.(noisemodel))
    #  isnothing(gatenoiseindex) && error("Noise model not defined for idling qubits!")
    #  gatenoise = last(noisemodel[gatenoiseindex])
    #  busy_qubits = vcat([collect(g[2]) for g in layer]...)
    #  idle_qubits = filter(y -> y ∉ busy_qubits, 1:nqubits(circuit))
    #  for n in idle_qubits
    #    push!(noisylayer, (gatenoise[1], n, gatenoise[2]))
    #  end
    #end
    push!(noisycircuit, noisylayer)
  end
  return noisycircuit
end

insertnoise(circuit::Vector{<:Any}, noisemodel::Tuple; kwargs...) =
  insertnoise([circuit], noisemodel; kwargs...)[1]

function maxgatesize(circuit::Vector{<:Vector{<:Any}})
  maxsize = 0
  for layer in circuit
    for g in layer
      maxsize = length(g[2]) > maxsize ? length(g[2]) : maxsize
    end
  end
  return maxsize
end

maxgatesize(circuit::Vector{<:Any}) =
  maxgatesize([circuit])


text:
# Useful overlays
"""
    my_overlay
Overlaying hallucinated trajectory on the ground truth
# Fields
- `color::Colorant`
- `scene::Scene`
"""
struct my_overlay
    scene::Scene
    color::Colorant # Needs to be of form colorant"Colorname"
end

function AutomotiveVisualization.add_renderable!(rendermodel::RenderModel,
     overlay::my_overlay)
    add_renderable!(rendermodel, overlay.scene,
    car_color=overlay.color)
    return rendermodel
end

"""
    curve_pts_overlay
Displays circles at the curve points that constitute the lanes of a road:
- color: the color of the point
- size: the size of the point

# Examples
```julia
render(scene,road,[curvepts_overlay(roadway_ext,colorant"yellow",0.05)])
```
"""
struct curvepts_overlay
    roadway
    color::Colorant # eg: colorant"yellow"
    size::Float64
end

function AutomotiveVisualization.add_renderable!(rendermodel::RenderModel,
     overlay::curvepts_overlay, scene::Scene, roadway::R) where {S,D,I,R}

    num_segments = length(roadway.segments)

    # Loop over segments
    for seg_i in 1:num_segments
        seg = roadway.segments[seg_i]
        num_lanes = length(seg.lanes)

        # Within a segment, loop over lanes
        for lane_i in 1:num_lanes
            lane = seg.lanes[lane_i]

            # Within a lane, loop over the curve points
            num_curvepts = length(lane.curve)
            for cpt_i in 1:num_curvepts
                cpt = lane.curve[cpt_i]

                add_instruction!(rendermodel,render_circle,
                    (cpt.pos.x,cpt.pos.y,overlay.size,overlay.color))
            end
        end
    end

end

"""
    LaneOverlay

- Intended to debug lanes and how vehicle 4 is teleporting from lane a to lane b2

# Examples
```julia
lane_overlay = LaneOverlay(roadway[LaneTag(1,1)], RGBA(0.0,0.0,1.0,0.5))
render(scene, roadway, [lane_overlay], cam=FitToContentCamera(0.))
```
"""
struct LaneOverlay
    lane::Lane
    color::Colorant
end

function AutomotiveVisualization.add_renderable!(rendermodel::RenderModel,
     overlay::LaneOverlay, scene::Scene, roadway::Roadway)
    add_instruction!(rendermodel, overlay.lane, roadway, color_asphalt=overlay.color) # this display a lane with the specified color
    return rendermodel
end

"""
    MergeOverlay

- Draw line joining vehicle and its associated merge vehicle
- Inspired by MergeNeighborsOverlay from AutonomousMerging.jl
"""
@with_kw mutable struct MergeOverlay
    env::MergingEnvironment = MergingEnvironment()
end
function AutomotiveVisualization.add_renderable!(rendermodel::RenderModel,
     overlay::MergeOverlay,scene::Scene, roadway::Roadway)
    for veh in scene
        ego_veh = veh
        merge_veh = find_merge_vehicle(overlay.env,scene,veh)
        if merge_veh != nothing
            A = get_front(ego_veh)
            B = get_rear(merge_veh)
            add_instruction!(rendermodel, render_line_segment, (A.x, A.y, B.x, B.y, colorant"blue", 0.5))
        end
    end
end

max_stars_repo_name: sisl/AutomotiveInteraction

max_issues_repo_name: sisl/AutomotiveInteraction

max_forks_repo_name: sisl/AutomotiveInteraction


text:
module TestMetaInfo_VariableTimestep

using Test
using Mimi

import Mimi:
    compdef, reset_compdefs, first_period, last_period

reset_compdefs()

@defcomp ch4forcing1 begin
    c_N2Oconcentration = Parameter(index=[time],unit="ppbv")
    c_CH4concentration = Parameter(index=[time],unit="ppbv")
    f0_CH4baseforcing = Parameter(unit="W/m2")
    fslope_CH4forcingslope = Parameter(unit="W/m2")
    c0_baseN2Oconc = Parameter(unit="ppbv")
    c0_baseCH4conc = Parameter(unit="ppbv")
    f_CH4forcing = Variable(index=[time],unit="W/m2")
    over_baseoverlap = Variable(unit="W/m2")
    over = Variable(index=[time],unit="W/m2")

    function run_timestep(p, v, d, t)
        #from p.16 in Hope 2009
        if is_first(t)
            #calculate baseline forcing overlap in first time period
            v.over_baseoverlap = -0.47 * log(1 + 2.0e-5 * (p.c0_baseN2Oconc * p.c0_baseCH4conc)^0.75 + 5.3e-15 * p.c0_baseCH4conc * (p.c0_baseCH4conc * p.c0_baseN2Oconc)^1.52)
        end

        v.over[t] = -0.47 * log(1 + 2.0e-5 * (p.c_CH4concentration[t] * p.c0_baseN2Oconc)^0.75 + 5.3e-15 * p.c_CH4concentration[t] * (p.c0_baseN2Oconc * p.c_CH4concentration[t])^1.52)
        v.f_CH4forcing[t] = p.f0_CH4baseforcing + p.fslope_CH4forcingslope * (sqrt(p.c_CH4concentration[t]) - sqrt(p.c0_baseCH4conc)) + v.over[t] - v.over_baseoverlap
    end
end

@Mimi.defmodel test_model begin
    index[time] = [2010:2100; 2105:5:2200]
    component(ch4forcing1)
    component(ch4forcing1, ch4forcing2) # add another one with a different name
end

c1 = compdef(test_model, :ch4forcing1)
c2 = compdef(test_model, :ch4forcing2)

@test c1 == compdef(:ch4forcing1)
@test_throws ErrorException compdef(:missingcomp)

@test c2.comp_id.module_name == :TestMetaInfo_VariableTimestep
@test c2.comp_id.comp_name == :ch4forcing1
@test c2.name == :ch4forcing2

vars = Mimi.variable_names(c2)
@test length(vars) == 3

pars = Mimi.parameter_names(c2)
@test length(pars) == 6

@test first_period(test_model.md, c1) == 2010
@test last_period(test_model.md, c1) == 2200

end # module


text:
@testset "$(rpad("Timeseries",80))" begin
    ntime = 10
    Δt = 0.1

    ts = TimeSeries{eltype(Δt)}(ntime, Δt, 1)
    @test typeof(ts) <: AbstractArray
    @test firstindex(ts)   == 0
    @test firstindex(ts,1) == firstindex(ts.t,1) - 1
    @test firstindex(ts,2) == 1
    @test lastindex(ts)    == ntime
    @test lastindex(ts,1)  == lastindex(ts.t,1) - 1
    @test lastindex(ts,2)  == 1
    @test axes(ts)   == (0:ntime,)
    @test axes(ts,1) == 0:ntime
    @test axes(ts,2) == 1:1
    @test size(ts.t) == (ntime+1,)

    compute_timeseries!(ts, 0.)
    t = collect(0:Δt:ntime*Δt)
    @test ts.t ≈ t atol=eps()

    ts1 = TimeSeries(ntime, Δt)
    compute_timeseries!(ts1, 0.)
    @test ts1 == ts
end


text:
fmt1 = "%.2f"
fmt2 = "%s"

tbl_dbl_width = begin
  r1 = Range(Box{Int}(Point(0,0), Point(2,1)), 1,1)
  r2 = Range(Box{Int}(Point(0,1), Point(1,2)), 2,1)
  r3 = Range(Box{Int}(Point(1,1), Point(2,2)), 3,2)

  layout = push(push(Layout([r1]), r2), r3)
  styles = [
    expandstyle(fmt2),
    uniformstyle(fmt1)
  ]
  data = Any[ExpandDataElement("Head",1,2), reshape(["< 15"],1,1), reshape([.5],1,1)]

  Table(data, layout, styles)
end

tbl_dbl_height = begin
  r1 = Range(Box{Int}(Point(0,0), Point(1,2)), 1,1)
  r2 = Range(Box{Int}(Point(1,0), Point(2,1)), 2,1)
  r3 = Range(Box{Int}(Point(1,1), Point(2,2)), 3,2)

  layout = push(push(Layout([r1]), r2), r3)
  styles = [
    expandstyle(fmt2),
    uniformstyle(fmt1)
  ]
  data = Any[ExpandDataElement("Head",2,1), reshape(["< 15"],1,1), reshape([.5],1,1)]

  Table(data, layout, styles)
end

tbl_hier_width = begin
  r1 = Range(Box{Int}(Point(0,0),Point(4,1)),1,1)
  r2 = Range(Box{Int}(Point(0,1),Point(4,3)),2,1)

  layout = push(Layout([r1]),r2)
  styles = [
    expandstyle(fmt2)
    cartesianstyle(fmt2)
  ]

  data = Any[
    ExpandDataElement(["Age / Sex"], 1, 4),
    CartesianDataElement(
      Vector[["< 15", ">= 15"],
             ["M", "F"]],true)]

  Table(data, layout, styles)
end

tbl_hier_height = begin
  r1 = Range(Box{Int}(Point(0,0),Point(1,4)),1,1)
  r2 = Range(Box{Int}(Point(1,0),Point(3,4)),2,1)

  layout = push(Layout([r1]),r2)
  styles = [
    expandstyle(fmt2)
    cartesianstyle(fmt2)
  ]

  data = Any[
    ExpandDataElement(["Age / Sex"], 4, 1),
    CartesianDataElement(
      Vector[["< 15", ">= 15"],
             ["M", "F"]],false)]

  Table(data, layout, styles)
end

tbl_pivottable = begin
  srand(100)
  r1 = Range(Box{Int}(Point(0,0),Point(2,2)),1,1) # Empty
  r2 = Range(Box{Int}(Point(2,0),Point(6,1)),2,1) # Col Header
  r3 = Range(Box{Int}(Point(2,1),Point(6,3)),3,2) # Col Levels
  r4 = Range(Box{Int}(Point(0,2),Point(2,3)),4,3) # Row Header
  r5 = Range(Box{Int}(Point(0,3),Point(2,5)),5,2)
  r6 = Range(Box{Int}(Point(2,3),Point(6,5)),6,4)

  layout = push(push(push(push(push(Layout([r1]), r2), r3), r4), r5), r6)
  styles = [
    expandstyle("%s"),
    cartesianstyle("%s"),
    uniformstyle("%s"),
    uniformstyle("%.2f"),
  ]

  data=Any[
    ExpandDataElement("",2,2),
    ExpandDataElement("Age / Sex",1,4),
    CartesianDataElement(
    Vector[["M", "F"],["\$\\leq\$ 15", "> 15"]], true),
    reshape(["Sex", "stats"],1,2),
    CartesianDataElement(
        Vector[["LT 6", "GT 6"],["mean"]], false),
    rand(2,4)
  ]

  Table(data, layout, styles)
end

tblnames = [:tbl_dbl_width;
        :tbl_dbl_height;
        :tbl_hier_width;
        :tbl_hier_height;
        :tbl_pivottable]


text:
using Test

for testscen in 1:2
    if testscen == 2
        continue # this validation data is missing
    end

    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = page_model()
    include("../src/components/RCPSSPScenario.jl")
    include("../src/components/EquityWeighting.jl")

    rcpsspscenario = addrcpsspscenario(m, scenario)
    equityweighting = add_comp!(m, EquityWeighting)

    equityweighting[:tct_percap_totalcosts_total] = readpagedata(m, "test/validationdata/$valdir/tct_per_cap_totalcostspercap.csv")
    equityweighting[:act_adaptationcosts_total] = readpagedata(m, "test/validationdata/$valdir/act_adaptationcosts_tot.csv")
    equityweighting[:act_percap_adaptationcosts] = readpagedata(m, "test/validationdata/$valdir/act_percap_adaptationcosts.csv")
    equityweighting[:cons_percap_aftercosts] = readpagedata(m, "test/validationdata/$valdir/cons_percap_aftercosts.csv")
    equityweighting[:cons_percap_consumption] = readpagedata(m, "test/validationdata/$valdir/cons_percap_consumption.csv")
    equityweighting[:cons_percap_consumption_0] = readpagedata(m, "test/validationdata/cons_percap_consumption_0.csv")
    equityweighting[:rcons_percap_dis] = readpagedata(m, "test/validationdata/$valdir/rcons_per_cap_DiscRemainConsumption.csv")
    equityweighting[:yagg_periodspan] = readpagedata(m, "test/validationdata/yagg_periodspan.csv")
    equityweighting[:pop_population] = readpagedata(m, "test/validationdata/$valdir/pop_population.csv")
    equityweighting[:y_year_0] = 2015.
    equityweighting[:y_year] = Mimi.dim_keys(m.md, :time)
    equityweighting[:grw_gdpgrowthrate] = rcpsspscenario[:grw_gdpgrowthrate]
    equityweighting[:popgrw_populationgrowth] = rcpsspscenario[:popgrw_populationgrowth]

    p = load_parameters(m)
    set_leftover_params!(m, p)

    run(m)

    # Generated data
    df = m[:EquityWeighting, :df_utilitydiscountfactor]
    wtct_percap = m[:EquityWeighting, :wtct_percap_weightedcosts]
    pct_percap = m[:EquityWeighting, :pct_percap_partiallyweighted]
    dr = m[:EquityWeighting, :dr_discountrate]
    dfc = m[:EquityWeighting, :dfc_consumptiondiscountrate]
    pct = m[:EquityWeighting, :pct_partiallyweighted]
    pcdt = m[:EquityWeighting, :pcdt_partiallyweighted_discounted]
    wacdt = m[:EquityWeighting, :wacdt_partiallyweighted_discounted]
    aact = m[:EquityWeighting, :aact_equityweightedadaptation_discountedaggregated]

    wit = m[:EquityWeighting, :wit_equityweightedimpact]
    addt = m[:EquityWeighting, :addt_equityweightedimpact_discountedaggregated]
    addt_gt = m[:EquityWeighting, :addt_gt_equityweightedimpact_discountedglobal]
    te = m[:EquityWeighting, :te_totaleffect]

    # Recorded data
    df_compare = readpagedata(m, "test/validationdata/df_utilitydiscountfactor.csv")
    wtct_percap_compare = readpagedata(m, "test/validationdata/$valdir/wtct_percap_weightedcosts.csv")
    pct_percap_compare = readpagedata(m, "test/validationdata/$valdir/pct_percap_partiallyweighted.csv")
    dr_compare = readpagedata(m, "test/validationdata/$valdir/dr_discountrate.csv")
    dfc_compare = readpagedata(m, "test/validationdata/$valdir/dfc_consumptiondiscountrate.csv")
    pct_compare = readpagedata(m, "test/validationdata/$valdir/pct_partiallyweighted.csv")
    pcdt_compare = readpagedata(m, "test/validationdata/$valdir/pcdt_partiallyweighted_discounted.csv")
    wacdt_compare = readpagedata(m, "test/validationdata/$valdir/wacdt_partiallyweighted_discounted.csv")
    aact_compare = readpagedata(m, "test/validationdata/$valdir/aact_equityweightedadaptation_discountedaggregated.csv")

    wit_compare = readpagedata(m, "test/validationdata/$valdir/wit_equityweightedimpact.csv")
    addt_compare = readpagedata(m, "test/validationdata/$valdir/addt_equityweightedimpact_discountedaggregated.csv")
    addt_gt_compare = (scenario == "NDCs" ? 9.960706559386551e8 : 3.445096140740376e8)
    te_compare = (scenario == "NDCs" ? 1.0320923880568126e9 : 5.0496576104580307e8)

    @test df ≈ df_compare rtol=1e-8
    @test wtct_percap ≈ wtct_percap_compare rtol=1e-7
    @test pct_percap ≈ pct_percap_compare rtol=1e-7
    @test dr ≈ dr_compare rtol=1e-5
    @test dfc ≈ dfc_compare rtol=1e-7
    @test pct ≈ pct_compare rtol=1e-3
    @test pcdt ≈ pcdt_compare rtol=1e-3
    @test wacdt ≈ wacdt_compare rtol=1e-4

    @test aact ≈ aact_compare rtol=1e-3

    @test wit ≈ wit_compare rtol=1e-3
    @test addt ≈ addt_compare rtol=1e-2
    @test addt_gt ≈ addt_gt_compare rtol=1e-2

    @test te ≈ te_compare rtol=1e-2
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
function testInverse(p::Float64, C::Float64, m::Int64)
  ϵ::Float64 = (p - C)/2.0
  testFunction(BernoulliFactory.inverse, m, C/p, makef(p), C, ϵ, GLOBAL_RNG)
end

Random.seed!(192837465)

testInverse(0.8, 0.4, ntrials)
testInverse(0.5, 0.2, ntrials)
testInverse(0.9, 0.7, ntrials)
testInverse(0.15, 0.1, ntrials)


text:
"A simple contruction helper that runs some sanity checks on the types"
@inline function game_sanity_checks(uids, TD, TC)
    @assert(isempty(intersect(uids...)),
            "Invalid uids: Two players can not control the same input")
    @assert(sum(length(uis) for uis in uids) == n_controls(TD),
            "Not all inputs have been assigned to players.")
    @assert(all(isbits(uir) for uir in uids),
            "Invalid uids: all ranges should be isbits to make things fast.")
    @assert(all(eltype(uir) == Int for uir in uids),
            "Invalid uids: the elements of the u_idx_range should be integers.")
end

"""
$(TYPEDEF)

A representation of a general game with potentially non-linear dynamics and
non-quadratic costs.
"""
struct GeneralGame{TH<:Integer,TU<:Tuple,TD<:ControlSystem, TC} <: AbstractGame
    h::TH
    uids::TU
    dyn::TD
    cost::TC

    function GeneralGame(h::TH, uids::TU, dyn::TD, cost::TC) where {TH,TU,TD,TC}
        game_sanity_checks(uids, TD, TC)
        return new{TH,TU,TD,TC}(h, uids, dyn, cost)
    end
end

dyntype(::Type{<:GeneralGame{TH,TU,TD}}) where {TH,TU,TD} = TD
dynamics(g::GeneralGame) = g.dyn
player_costs(g::GeneralGame) = g.cost
horizon(g::GeneralGame) = g.h
uindex(g::GeneralGame) = g.uids

"""
    $(TYPEDSIGNATURES)

Transforms a `GeneralGame` `g` with feedback linearizable dynamics (i.e.
`LinearizationStyle(dynamics(g)) isa FeedbackLinearization`) to it's linearized
form.
"""
function transform_to_feedbacklin(g::GeneralGame, x0)
    @assert LinearizationStyle(dynamics(g)) isa FeedbackLinearization
    # transform the dynamics
    lin_dyn = feedbacklin(dynamics(g))
    # approximate the cost by a cost in ξ coordinates
    ξ_cost = map(enumerate(player_costs(g))) do (i, c)
        transformed_cost(dynamics(g), c, n_players(g))
    end |> SVector{n_players(g)}

    # transformed game
    gξ = GeneralGame(horizon(g), uindex(g), lin_dyn, ξ_cost)
    # transform initial conditions
    ξ0 = ξ_from(dynamics(g), x0)
    if λ_issingular(dynamics(g), ξ0)
        @warn "State conversion map is singular at provided initial conditions."
    end

    return gξ, ξ0
end

"""
$(TYPEDEF)

A struct to represent a multi-player linear system where each player controls a
range of inputs.

# Parameters

- `TD`:   the type of the dynamics
- `TC`:   the type of the player costs

# Fields

$(TYPEDFIELDS)
"""
struct LQGame{TU,TD<:Union{LTVSystem, LTISystem}, TC} <: AbstractGame
    uids::TU
    "The full linear system dynamics. A vector (time) of `LinearSystem`s."
    dyn::TD
    "The cost representation. A vector (time) of vector (player) of
    `QuadraticPlayerCost`"
    pcost::TC

    function LQGame(uids::TU, dyn::TD, pcost::TC) where {TU,TD,TC}
        game_sanity_checks(uids, TD, eltype(TC))
        @assert eltype(eltype(TC)) <: QuadraticPlayerCost "LQGames require quadratic costs."
        return new{TU,TD,TC}(uids, dyn, pcost)
    end
end

dynamics(g::LQGame) = g.dyn
player_costs(g::LQGame) = g.pcost
horizon(g::LQGame) = length(g.pcost)
uindex(g::LQGame) = g.uids

function lq_approximation(solver, g::GeneralGame, op::SystemTrajectory)
    lqg = lqgame_preprocess_alloc(g)
    lq_approximation!(lqg, solver, g, op)
    return lqg
end

function lq_approximation!(lqg::LQGame, solver, g::GeneralGame,
                           op::SystemTrajectory)
    for (k, (xₖ, uₖ)) in enumerate(zip(op.x, op.u))
        # discrete linearization along the operating point
        t = time_disc2cont(op, k)
        linearize!(lqg, dynamics(g), xₖ, uₖ, t, k)
        # quadratiation of the cost along the operating point
        player_costs(lqg)[k] = map(player_costs(g)) do pcₖⁱ
            c = quadraticize!(qcache(solver), pcₖⁱ, g, xₖ, uₖ, t)
            return regularize(solver, c)
        end
    end
end

@inline function linearize!(lqg::LQGame, dyn::ControlSystem, xₖ, uₖ, t, k)
    return linearize!(LinearizationStyle(dyn), lqg, dyn, xₖ, uₖ, t, k)
end

@inline function linearize!(::LinearizationStyle, lqg, dyn, xₖ, uₖ, t, k)
    return dynamics(lqg)[k] = linearize_discrete(dyn, xₖ, uₖ, t)
end
# TODO: Maybe we can do this less subtle. Currenty we simply do nothing because the
# LQGame already holds the right linear dynamics (created in
# lqgame_preprocess_alloc)
@inline linearize!(::TrivialLinearization, lqg, args...) = dynamics(lqg)


text:
using BayesianNonparametrics
using Test
using Distributions

# write your own tests here
include("distributionTests.jl")
include("dpmTests.jl")


text:
#
# Functions to create DiffEq problems
#
using OrdinaryDiffEq

function _problem_setup(sim, tend, inputs)
    pybamm = pyimport("pybamm")

    input_parameter_order = isnothing(inputs) ? nothing : collect(keys(inputs))
    p = isnothing(inputs) ? nothing : collect(values(inputs))

    sim.build()
    fn_str, u0_str = sim.built_model.generate_julia_diffeq(
        input_parameter_order=input_parameter_order,
        dae_type="implicit",
        get_consistent_ics_solver=pybamm.CasadiSolver()
    )

    # PyBaMM-generated functions
    sim_fn! = runtime_eval(Meta.parse(fn_str))
    sim_u0! = runtime_eval(Meta.parse(u0_str))

    # Evaluate initial conditions
    len_y = convert(Int, sim.built_model.len_rhs_and_alg)
    u0 = Array{Float64}(undef, len_y)
    sim_u0!(u0, p)

    # Scale the time
    tau = sim.built_model.timescale.evaluate()
    tspan = (0, tend/tau)
    sim_fn!, u0, tspan, p
end

function get_ode_problem(sim, tend, inputs)
    sim_fn!, u0, tspan, p = _problem_setup(sim, tend, inputs)

    # Create problem, isinplace is explicitly true as cannot be inferred from
    # runtime_eval function
    ODEProblem{true}(sim_fn!, u0, tspan, p)
end

# Defaults
get_ode_problem(sim, tend::Real) = get_ode_problem(sim, tend, nothing)
get_ode_problem(sim, inputs::AbstractDict) = get_ode_problem(sim, 3600, inputs)
get_ode_problem(sim) = get_ode_problem(sim, 3600, nothing)

function get_dae_problem(sim, tend, inputs)
    sim_fn!, u0, tspan, p = _problem_setup(sim, tend, inputs)

    # Create vector of 1s and 0s to indicate differential and algebraic variables
    len_rhs = convert(Int, sim.built_model.len_rhs)
    len_alg = convert(Int, sim.built_model.len_alg)
    differential_vars = vcat(ones(len_rhs), zeros(len_alg))

    # Create problem, isinplace is explicitly true as cannot be inferred from
    # runtime_eval function
    du0 = zeros(size(u0))
    DAEProblem{true}(sim_fn!, du0, u0, tspan, p, differential_vars=differential_vars)
end

# Defaults
get_dae_problem(sim, tend::Real) = get_dae_problem(sim, tend, nothing)
get_dae_problem(sim, inputs::AbstractDict) = get_dae_problem(sim, 3600, inputs)
get_dae_problem(sim) = get_dae_problem(sim, 3600, nothing)


text:
mutable struct SwapRateTrigger <: ExerciseStrategy
  rateTimes::Vector{Float64}
  swapTriggers::Vector{Float64}
  exerciseTimes::Vector{Float64}
  currentIndex::Int
  rateIndex::Vector{Int}
end

function SwapRateTrigger(rateTimes::Vector{Float64}, swapTriggers::Vector{Float64}, exerciseTimes::Vector{Float64})
  check_increasing_times(rateTimes)
  length(rateTimes) > 1 || error("Rate times must contain at least two values")

  check_increasing_times(exerciseTimes)
  length(swapTriggers) == length(exerciseTimes) || error("swapTriggers / exerciseTimes mismatch")

  rateIndex = Vector{Int}(undef, length(exerciseTimes))
  j = 1
  @inbounds @simd for i in eachindex(exerciseTimes)
    while j <= length(rateTimes) && rateTimes[j] < exerciseTimes[i]
      j += 1
    end

    rateIndex[i] = j
  end

  return SwapRateTrigger(rateTimes, swapTriggers, exerciseTimes, 1, rateIndex)
end

relevant_times(srt::SwapRateTrigger) = srt.exerciseTimes

clone(srt::SwapRateTrigger) = SwapRateTrigger(copy(srt.rateTimes), copy(srt.swapTriggers), copy(srt.exerciseTimes), srt.currentIndex, copy(srt.rateIndex))


text:
export atomium
###
# The atomium: a famous building in Belgium
###
cylinder(::Type{T} = Float64) where {T} = UnitDisk{T}() × UnitInterval{T}()

cylinder(radius::T, length::T) where {T} = radius .* UnitDisk{T}() × (0 .. length)

function atomium()
    sphere1 = ball(0.25)
    spheres = union(sphere1,sphere1 .+ SVector( 0.6, 0.6, 0.6))
    spheres= union(spheres, sphere1 .+ SVector( 0.6, 0.6,-0.6))
    spheres= union(spheres, sphere1 .+ SVector( 0.6,-0.6, 0.6))
    spheres= union(spheres, sphere1 .+ SVector( 0.6,-0.6,-0.6))
    spheres= union(spheres, sphere1 .+ SVector(-0.6, 0.6, 0.6))
    spheres= union(spheres, sphere1 .+ SVector(-0.6, 0.6,-0.6))
    spheres= union(spheres, sphere1 .+ SVector(-0.6,-0.6, 0.6))
    spheres= union(spheres, sphere1 .+ SVector(-0.6,-0.6,-0.6))
    cyl1 = cylinder(0.10, 1.2)
    spheres= union(spheres, cyl1 .+ SVector(0.6, 0.6, -0.6));
    spheres= union(spheres, cyl1 .+ SVector(0.6,-0.6, -0.6));
    spheres= union(spheres, cyl1 .+ SVector(-0.6, 0.6,-0.6));
    spheres= union(spheres, cyl1 .+ SVector(-0.6,-0.6,-0.6));
    cyl2 = DomainSets.rotate(cyl1, 0.0, pi/2.0, 0.0)
    spheres= union(spheres, cyl2 .+ SVector(-0.6,  0.6, 0.6))
    spheres= union(spheres, cyl2 .+ SVector(-0.6, -0.6, 0.6))
    spheres= union(spheres, cyl2 .+ SVector(-0.6,  0.6,-0.6))
    spheres= union(spheres, cyl2 .+ SVector(-0.6, -0.6,-0.6))
    cyl2b = DomainSets.rotate(cyl1, pi/2.0, 0.0, 0.0)
    spheres= union(spheres, cyl2b .+ SVector( 0.6, 0.6,-0.6))
    spheres= union(spheres, cyl2b .+ SVector(-0.6, 0.6, 0.6))
    spheres= union(spheres, cyl2b .+ SVector( 0.6, 0.6, 0.6))
    spheres= union(spheres, cyl2b .+ SVector(-0.6, 0.6,-0.6))
    cyl3 = cylinder(0.10, 1.2*sqrt(3))
    cyl3 = DomainSets.rotate(cyl3, 0.0, acos(1/sqrt(3)), 0.0)
    cyl3 = DomainSets.rotate(cyl3, 0.0, 0.0, pi/4)
    spheres= union(spheres, cyl3 .+ SVector( -0.6, -0.6, -0.6))
    cyl4 = cylinder(0.10, 1.2*sqrt(3))
    cyl4 = DomainSets.rotate(cyl4, 0.0, -acos(1/sqrt(3)), 0.0)
    cyl4 = DomainSets.rotate(cyl4, 0.0, 0.0, pi/4)
    spheres= union(spheres, cyl4 .+ SVector(  0.6, 0.6,-0.6))
    cyl5 = cylinder(0.10, 1.2*sqrt(3))
    cyl5 = DomainSets.rotate(cyl5, 0.0, acos(1/sqrt(3)), 0.0)
    cyl5 = DomainSets.rotate(cyl5, 0.0, 0.0, -pi/4)
    spheres= union(spheres, cyl5 .+ SVector( -0.6, +0.6, -0.6))
    cyl6 = cylinder(0.10, 1.2*sqrt(3))
    cyl6 = DomainSets.rotate(cyl6, 0.0, -acos(1/sqrt(3)), 0.0)
    cyl6 = DomainSets.rotate(cyl6, 0.0, 0.0, -pi/4)
    spheres= union(spheres, cyl6 .+ SVector( 0.6, -0.6, -0.6))
    atomium = spheres
end


text:
ALL_TASK_OFFSETS[("linux-x86", v"1.8.0")] = Dict(
    :END => 256,
    :_isexception => 30,
    :_state => 28,
    :bufsz => 248,
    :copy_stack => 252,
    :ctx => 88,
    :donenotify => 12,
    :eh => 84,
    :excstack => 80,
    :gcstack => 64,
    :logstate => 20,
    :next => 0,
    :prio => 78,
    :ptls => 72,
    :queue => 4,
    :result => 16,
    :rngState0 => 32,
    :rngState1 => 40,
    :rngState2 => 48,
    :rngState3 => 56,
    :sizeof_ctx => 156,
    :start => 24,
    :sticky => 29,
    :stkbuf => 244,
    :tid => 76,
    :tls => 8,
    :tls_base_context => 3416,
    :world_age => 68,
)


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCell
include("useful_routines.jl")

vectosymmconverter = CutCell.vector_to_symmetric_matrix_converter()
E1 = [
    1.0 0.0
    0.0 0.0
    0.0 1.0
]

E2 = [
    0.0 0.0
    0.0 1.0
    1.0 0.0
]

@test allapprox(vectosymmconverter[1], E1)
@test allapprox(vectosymmconverter[2], E2)

matrix = CutCell.make_row_matrix(E1, [1, 2, 3])
expE1 = hcat(E1, 2 * E1, 3 * E1)
@test allapprox(matrix, expE1)

matrix = CutCell.interpolation_matrix([1.0, 2.0, 3.0], 2)
testm = [
    1.0 0.0 2.0 0.0 3.0 0.0
    0.0 1.0 0.0 2.0 0.0 3.0
]
@test allapprox(matrix, testm)

basis = TensorProductBasis(2, 1)
cellmap = CutCell.CellMap([0.0, 0.0], [3.0, 1.0])
grad = CutCell.transform_gradient(gradient(basis, 0.0, 0.0), CutCell.jacobian(cellmap))
@test allapprox(grad[1, :], [-1 / 6, -1 / 2])
@test allapprox(grad[2, :], [-1 / 6, 1 / 2])
@test allapprox(grad[3, :], [1 / 6, -1 / 2])
@test allapprox(grad[4, :], [1 / 6, 1 / 2])

l, m = 1.0, 2.0
stiffness = plane_strain_voigt_hooke_matrix(1.0, 2.0)
cellmap = CutCell.CellMap([0.0, 0.0], [1.0, 1.0])
quad = tensor_product_quadrature(2, 2)
bf = CutCell.bilinear_form(basis, quad, stiffness, CutCell.jacobian(cellmap))
@test size(bf) == (8, 8)

rhsfunc(x) = [1.0, 1.0]
rhs = CutCell.linear_form(rhsfunc, basis, quad, cellmap)
@test allapprox(rhs, 0.25 * ones(8))

basis = TensorProductBasis(2, 2)
quad = tensor_product_quadrature(2, 4)
cellmap = CutCell.CellMap([-1, -1.0], [1.0, 1.0])
M = CutCell.mass_matrix(basis, quad, cellmap, 2)
rhsfunc(x) = [x[1]^2 + x[2]^2, 2x[1] * x[2]]
R = CutCell.linear_form(rhsfunc, basis, quad, cellmap)

sol = reshape(M \ R, 2, :)
coords = hcat([cellmap(basis.points[:, i]) for i = 1:9]...)
exactsol = hcat([rhsfunc(coords[:, i]) for i = 1:9]...)
@test allapprox(sol, exactsol, 1e3eps())

K = CutCell.bilinear_form(basis, quad, stiffness, CutCell.jacobian(cellmap))
bodyforce(x) = -4 * (l + 2m) * [1.0, 0.0]
R = CutCell.linear_form(bodyforce, basis, quad, cellmap)

boundarynodeids = [1, 4, 7, 8, 9, 6, 3, 2]
bcvals = exactsol[:, boundarynodeids]
CutCell.apply_dirichlet_bc!(K, R, boundarynodeids, bcvals)

sol = reshape(K \ R, 2, :)
@test allapprox(sol, exactsol, 1e3eps())

K = CutCell.bilinear_form(basis, quad, stiffness, cellmap)
bodyforce(x) = -4 * (l + 2m) * [1.0, 0.0]
R = CutCell.linear_form(bodyforce, basis, quad, cellmap)

boundarynodeids = [1, 4, 7, 8, 9, 6, 3, 2]
bcvals = exactsol[:, boundarynodeids]
CutCell.apply_dirichlet_bc!(K, R, boundarynodeids, bcvals)

sol = reshape(K \ R, 2, :)
@test allapprox(sol, exactsol, 1e3eps())

polyorder = 1
numqp = 2
basis = TensorProductBasis(2, polyorder)
quad = tensor_product_quadrature(2, numqp)
lambda, mu = 1.0, 2.0
penalty = 1e0
dx = 0.1
e11 = dx / 1.0
e22 = -lambda / (lambda + 2mu) * e11
dy = e22

cellmap = CutCell.CellMap([0.0, 0.0], [1.0, 1.0])
stiffness = plane_strain_voigt_hooke_matrix(lambda, mu)
facequads = CutCell.face_quadratures(numqp)
facedetjac = CutCell.face_determinant_jacobian(cellmap)
bf = CutCell.bilinear_form(basis, quad, stiffness, cellmap)

ftop1 = CutCell.face_traction_component_operator(
    basis,
    facequads[1],
    [0.,1.],
    [0.,-1.],
    stiffness,
    facedetjac[1],
    cellmap,
)
ftop2 = CutCell.face_traction_component_operator(
    basis,
    facequads[2],
    [1.,0.],
    [1.,0.],
    stiffness,
    facedetjac[2],
    cellmap,
)
ftop4 = CutCell.face_traction_component_operator(
    basis,
    facequads[4],
    [1.,0.],
    [-1.0, 0],
    stiffness,
    facedetjac[4],
    cellmap,
)

bcop1 = penalty*CutCell.component_mass_matrix(basis,facequads[1],[0.,1.],facedetjac[1])
bcop2 = penalty*CutCell.component_mass_matrix(basis,facequads[2],[1.,0.],facedetjac[2])
bcop4 = penalty*CutCell.component_mass_matrix(basis,facequads[4],[1.,0.],facedetjac[4])
bcrhs2 = penalty*CutCell.component_linear_form(dx,basis,facequads[2],[1.,0.],facedetjac[2])

K = bf - ftop1 - ftop2 - ftop4 + bcop1 + bcop2 + bcop4
rhs = bcrhs2

sol = K\rhs
disp = reshape(sol,2,:)
testdisp = [0. 0. dx dx
            0. dy 0. dy]
@test allapprox(disp,testdisp,1e3eps())


text:
# Get approximated hessian (posterior precision matrix)
function vi_get_apprhess!(H::MatIO, β::Real, J::MatI, Λ0::MatI)
    unsafe_copy!(H, Λ0) # H ← Λ0
    BLAS.syrk!('U', 'T', β, J, 1.0, H)
    @inbounds for j in axes(H, 1), i in 1:j-1
        H[j,i] = H[i,j]
    end
    return nothing
end

function vi_get_apprhess!(H::MatIO, λy::VecI, J::MatI, Λ0::MatI, ΛJ::MatB)
    for j in axes(J, 2)
        @simd for i in eachindex(λy)
            @inbounds ΛJ[i,j] = λy[i] * J[i,j]
        end
    end
    unsafe_copy!(H, Λ0)                      # H ← Λ0
    BLAS.gemm!('T', 'N', 1.0, ΛJ, J, 1.0, H) # H ← JΛ * J + Λ0
    return nothing
end

function vi_get_apprhess!(H::MatIO, Λy::MatI, J::MatI, Λ0::MatI, ΛJ::MatB)
    unsafe_copy!(H, Λ0)                       # H  ← Λ0
    BLAS.symm!('L', 'U', 1.0, Λy, J, 0.0, ΛJ) # ΛJ ← Λy * J
    BLAS.gemm!('T', 'N', 1.0, ΛJ, J, 1.0,  H) # H  ← JΛ * J + Λ0
    return nothing
end

# Get 1st-order variational inference gradient
function vi_get_gradient!(d::VecIO, β::Real, J::MatI, r::VecI, Λ0::MatI, Δμ::VecI)
    BLAS.symv!('U', 1.0, Λ0, Δμ,  0.0, d) # d ← Λ0  * Δμ
    BLAS.gemv!('T',   β,  J,  r, -1.0, d) # d ← βJ' * rs - Λ0  * Δμ
    return nothing
end

function vi_get_gradient!(d::VecIO, ΛJ::MatI, r::VecI, Λ0::MatI, Δμ::VecI)
    BLAS.symv!('U', 1.0, Λ0, Δμ,  0.0, d) # d ← Λ0  * Δμ
    BLAS.gemv!('T', 1.0, ΛJ,  r, -1.0, d) # d ← ΛJ' * rs - Λ0  * Δμ
    return nothing
end

# Perform BLAS.trsm!
function vi_get_blastrsm!(m::MatO, A::MatI)
    BLAS.trsm!('L', 'L', 'N', 'N', 1.0, A, m)
    BLAS.trsm!('L', 'L', 'T', 'N', 1.0, A, m)
    return nothing
end

# Get evidence (free energy) of variational inference
function vi_get_evidence!(Λy::Union{Real, VecI, MatI}, rs::VecI, ny::Int, Λ0::Union{Real, VecI, MatI}, Δμ::VecI, nd::Int, Fs::MatB, Λb::MatB)
    lm_get_cholesky!(Fs)
    vi_get_blastrsm!(Λb, Fs)
    return lm_get_squaresum(Λy, rs, ny) + lm_get_squaresum(Λ0, Δμ, nd) + 0.5 * (logdet(Fs, nd) + tr(Λb, nd))
end

struct VariationalInferenceOptimizer
    @def prop Vector{Float64} rs rt μs μt dμ δμ Δμ
    @def prop Matrix{Float64} Js Jt ΛJ Λs Λt Λb Ft Fb
    @def prop Int nd ny

    function VariationalInferenceOptimizer(nd::Int, ny::Int)
        @def vars Vector{Float64}(undef, ny) rs rt
        @def vars Vector{Float64}(undef, nd) μs μt dμ δμ Δμ
        @def vars Matrix{Float64}(undef, ny, nd) Js Jt ΛJ
        @def vars Matrix{Float64}(undef, nd, nd) Λs Λt Λb Ft Fb
        return new(rs, rt, μs, μt, dμ, δμ, Δμ, Js, Jt, ΛJ, Λs, Λt, Λb, Ft, Fb, nd, ny)
    end
end

function vi_get_optimize!(xs::VecIO, fn::AbstractObjective, μ0::VecI, Λ0::MatI, βy::Real, τ::Real, h::Real, itmax::Int, vi::VariationalInferenceOptimizer)
    #### Allocations ####
    @get vi nd ny rs rt μs μt dμ δμ Δμ
    @get vi Js Jt Λs Λt Λb Ft Fb
    #### Initialization
    ih1 = inv(h)
    ih2 = ih1 * ih1
    toN = eachindex(1:nd)
    unsafe_copy!(μs, μ0)
    unsafe_copy!(Λs, Λ0)
    #### Compute precision matrix
    lm_get_residual!(rs, μs, fn)
    lm_get_jacobian!(Js, μs, fn)
    vi_get_apprhess!(Λt, βy, Js, Λ0)
    # Controlling parameters of trust region
    a = τ * diagmax(Λt, nd); b = 2.0
    #### Compute free energy
    fill!(Δμ, 0.0)
    unsafe_copy!(Ft, Λs)
    unsafe_copy!(Fb, Λt)
    Fnow = vi_get_evidence!(βy, rs, ny, Λ0, Δμ, nd, Ft, Fb)
    #### Iteration
    it = 0
    while it < itmax
        it += 1
        # Damped Gauss-Newton approximation
        unsafe_copy!(Ft, Λt)
        λ = a * diagmax(Ft, nd)
        @simd for i in toN
            @inbounds Ft[i,i] += λ
        end
        cholesky_state = lm_get_cholesky!(Ft)
        # Levenberg-Marquardt step by 1st order linearization
        vi_get_gradient!(dμ, βy, Js, rs, Λ0, Δμ)
        lm_get_blastrsv!(dμ, Ft)
        # Finite-difference of 2nd order directional derivative
        @simd for i in toN
            @inbounds μt[i] = μs[i] + h * dμ[i]
        end
        lm_get_residual!(rt, μt, fn)
        lm_get_gradient!(rt, rs, Js, dμ, ih1, ih2)
        # Geodesic acceleration step by 2nd order linearization
        lm_get_gradient!(δμ, βy, Js, rt)
        lm_get_blastrsv!(δμ, Ft)
        gratio = 2.0 * nrm2(δμ) / nrm2(dμ)
        gratio > 0.75 && (a *= b; b *= 2.0; continue)
        # (Levenberg-Marquardt velocity) + (Geodesic acceleration)
        BLAS.axpy!(1.0, δμ, dμ) # dμ ← dμ + δμ
        # Compute predicted gain of the least square cost by linear approximation
        LinApprox = lm_get_squaresum(Λt, dμ, nd) + λ * dot(dμ, dμ, nd)
        # Compute gain ratio
        @simd for i in toN
            @inbounds μt[i] = μs[i] + dμ[i]
        end
        lm_get_residual!(rt, μt, fn)
        lm_get_jacobian!(Jt, μt, fn)
        vi_get_apprhess!(Λb, βy, Jt, Λ0)
        @simd for i in toN
            @inbounds δμ[i] = μt[i] - μ0[i]
        end
        unsafe_copy!(Ft, Λt)
        unsafe_copy!(Fb, Λb)
        Fnew = vi_get_evidence!(βy, rt, ny, Λ0, δμ, nd, Ft, Fb)

        Δ = Fnow - Fnew; ρ = Δ / LinApprox
        if ρ > 0.0
            Fnow = Fnew
            unsafe_copy!(μs, μt)
            unsafe_copy!(Λs, Λt)
            (maximum(δμ) < 1e-15 || Δ < 1e-10) && break # Check localized convergence
            unsafe_copy!(rs, rt)
            unsafe_copy!(Λt, Λb)
            unsafe_copy!(Js, Jt)
            unsafe_copy!(Δμ, δμ)
            a = ρ < 0.9367902323681495 ? a * (1.0 - cubic(2.0 * ρ - 1.0)) : a / 3.0
            b = 2.0
        else
            a > 1e10 && break # Check localized divergence
            a *= b; b *= 2.0
        end
    end
    unsafe_copy!(xs, μs); return xs
end

function vi_get_optimize!(xs::VecIO, fn::AbstractObjective, μ0::VecI, Λ0::MatI, Λy::Union{VecI, MatI}, τ::Real, h::Real, itmax::Int, vi::VariationalInferenceOptimizer)
    #### Allocations ####
    @get vi nd ny rs rt μs μt dμ δμ Δμ
    @get vi Js Jt ΛJ Λs Λt Λb Ft Fb
    #### Initialization
    ih1 = inv(h)
    ih2 = ih1 * ih1
    toN = eachindex(1:nd)
    unsafe_copy!(μs, μ0)
    unsafe_copy!(Λs, Λ0)
    #### Compute precision matrix
    lm_get_residual!(rs, μs, fn)
    lm_get_jacobian!(Js, μs, fn)
    vi_get_apprhess!(Λt, Λy, Js, Λ0, ΛJ) # ΛJ = Λy * Js + Λ0
    # Controlling parameters of trust region
    a = τ * diagmax(Λt, nd); b = 2.0
    #### Compute free energy
    fill!(Δμ, 0.0)
    unsafe_copy!(Ft, Λs)
    unsafe_copy!(Fb, Λt)
    Fnow = vi_get_evidence!(Λy, rs, ny, Λ0, Δμ, nd, Ft, Fb)
    #### Iteration
    it = 0
    while it < itmax
        it += 1
        # Damped Gauss-Newton approximation
        unsafe_copy!(Ft, Λt)
        λ = a * diagmax(Ft, nd)
        @simd for i in toN
            @inbounds Ft[i,i] += λ
        end
        cholesky_state = lm_get_cholesky!(Ft)
        # Levenberg-Marquardt step by 1st order linearization
        vi_get_gradient!(dμ, ΛJ, rs, Λ0, Δμ)
        lm_get_blastrsv!(dμ, Ft)
        # Finite-difference of 2nd order directional derivative
        @simd for i in toN
            @inbounds μt[i] = μs[i] + h * dμ[i]
        end
        lm_get_residual!(rt, μt, fn)
        lm_get_gradient!(rt, rs, Js, dμ, ih1, ih2)
        # Geodesic acceleration step by 2nd order linearization
        lm_get_gradient!(δμ, ΛJ, rt)
        lm_get_blastrsv!(δμ, Ft)
        gratio = 2.0 * nrm2(δμ) / nrm2(dμ)
        gratio > 0.75 && (a *= b; b *= 2.0; continue)
        # (Levenberg-Marquardt velocity) + (Geodesic acceleration)
        BLAS.axpy!(1.0, δμ, dμ) # dμ ← dμ + δμ
        # Compute predicted gain of the least square cost by linear approximation
        LinApprox = lm_get_squaresum(Λt, dμ, nd) + λ * dot(dμ, dμ, nd)
        # Compute gain ratio
        @simd for i in toN
            @inbounds μt[i] = μs[i] + dμ[i]
        end
        lm_get_residual!(rt, μt, fn)
        lm_get_jacobian!(Jt, μt, fn)
        vi_get_apprhess!(Λb, Λy, Jt, Λ0, ΛJ)
        @simd for i in toN
            @inbounds δμ[i] = μt[i] - μ0[i]
        end
        unsafe_copy!(Ft, Λt)
        unsafe_copy!(Fb, Λb)
        Fnew = vi_get_evidence!(Λy, rt, ny, Λ0, δμ, nd, Ft, Fb)

        Δ = Fnow - Fnew; ρ = Δ / LinApprox
        if ρ > 0.0
            Fnow = Fnew
            unsafe_copy!(μs, μt)
            unsafe_copy!(Λs, Λt)
            (maximum(δμ) < 1e-15 || Δ < 1e-10) && break # Check localized convergence
            unsafe_copy!(rs, rt)
            unsafe_copy!(Λt, Λb)
            unsafe_copy!(Js, Jt)
            unsafe_copy!(Δμ, δμ)
            a = ρ < 0.9367902323681495 ? a * (1.0 - cubic(2.0 * ρ - 1.0)) : a / 3.0
            b = 2.0
        else
            a > 1e10 && break # Check localized divergence
            a *= b; b *= 2.0
        end
    end
    unsafe_copy!(xs, μs); return xs
end

optimize!(xs::VecIO, fn::AbstractObjective, μ0::VecI, Λ0::MatI, Λy::Union{Real, VecI, MatI}, vi::VariationalInferenceOptimizer; τ::Real=1e-3, h::Real=0.1, itmax::Int=100) = vi_get_optimize!(xs, fn, μ0, Λ0, Λy, τ, h, itmax, vi)