text:
export ft_eeprom_read, ft_eeprom_program
export Fteeprom232b, Fteeprom2232, Fteeprom232r
export Fteeprom2232h, Fteeprom4232h, Fteeprom232h
export Fteepromxseries

include("FT_DEVICE.jl")  # load constants

abstract Eeprom

function fteepromread_alleeprom{T<:Eeprom}(ft_handle::Culong, eepromdata::T)
  datasize = Cuint(sizeof(eepromdata))
  mfg = Array(UInt8,64); mfg[64] = 0
  mfgid = Array(UInt8,64); mfgid[64] = 0
  d = Array(UInt8,64); d[64] = 0
  sn = Array(UInt8,64); sn[64] = 0
  ee = Ref{T}(eepromdata)
  ft_status = ccall((:FT_EEPROM_Read, d2xx),
            Cuint,
            (Culong,Ref{T},Cuint,Ptr{UInt8},
              Ptr{UInt8},Ptr{UInt8},Ptr{UInt8}),
            ft_handle,ee,datasize,mfg,mfgid,d,sn)
  checkstatus(ft_status)
  mfg_string = bytestring(mfg[1:findfirst(mfg,0x00)-1])
  mfgid_string = bytestring(mfgid[1:findfirst(mfgid,0x00)-1])
  d_string = bytestring(d[1:findfirst(d,0x00)-1])
  sn_string = bytestring(sn[1:findfirst(sn,0x00)-1])
  return(eepromdata,mfg_string,mfgid_string,d_string,sn_string)
end

function fteepromprogram_alleeprom{T<:Eeprom}(ft_handle::Culong, eepromdata::T,
          mfg_string::ASCIIString, mfgid_string::ASCIIString,
          d_string::ASCIIString, sn_string::ASCIIString)
  datasize = sizeof(eepromdata)
  mfg = Array{UInt8,1}(mfg_string * "\0")
  mfgid = Array{UInt8,1}(mfgid_string * "\0")
  d = Array{UInt8,1}(d_string * "\0")
  sn = Array{UInt8,1}(sn_string * "\0")
  ee = Ref{T}(eepromdata)
  ft_status = ccall((:FT_EEPROM_Program, d2xx),
            Cuint,
            (Culong,Ref{T},Cuint,Ptr{UInt8},
              Ptr{UInt8},Ptr{UInt8},Ptr{UInt8}),
            ft_handle,ee,datasize,mfg,mfgid,d,sn)
  checkstatus(ft_status)
  return nothing
end

function showcommonelements(io::IO, eepromdata::Eeprom)
  f()=@printf(io,"deviceType = 0x%08x\n",eepromdata.deviceType);f()
  f()=@printf(io,"VendorId = 0x%04x\n",eepromdata.VendorId);f()
  f()=@printf(io,"ProductId = 0x%04x\n",eepromdata.ProductId);f()
  f()=@printf(io,"SerNumEnable = 0x%02x\n",eepromdata.SerNumEnable);f()
  f()=@printf(io,"MaxPower = 0x%04x\n",eepromdata.MaxPower);f()
end

# FT232B EEPROM structure for use with ft_eeprom_read and ft_eeprom_program
type Fteeprom232b <: Eeprom
  # Common header
  ### BEGIN common elements for all device EEPROMs ###
  deviceType :: Cuint # FT_DEVICE FTxxxx device type to be programmed
  # Device descriptor options
  VendorId :: Cshort # 0x0403
  ProductId :: Cshort # 0x6001
  SerNumEnable :: Cuchar # non-zero if serial number to be used
  # Config descriptor options
  MaxPower :: Cshort #  0 < MaxPower <= 500
  pad1::UInt8;pad2::UInt8;pad3::UInt8;pad4::UInt8  # needed to allign fields.  Why?
  ### END common elements for all device EEPROMs ###
  Fteeprom232b() = new(FT_DEVICE_232BM)
  Fteeprom232b(VendorId,ProductId,SerNumEnable,MaxPower) =
    new(FT_DEVICE_232BM,VendorId,ProductId,SerNumEnable,MaxPower)
end

function Base.show(io::IO, eepromdata::Fteeprom232b)
  showcommonelements(io,eepromdata)
end

function ft_eeprom_read(ft_handle::Culong, eepromdata::Fteeprom232b)
  @assert eepromdata.deviceType == FT_DEVICE_232BM
  fteepromread_alleeprom(ft_handle,eepromdata)
end

function ft_eeprom_program(ft_handle::Culong, eepromdata::Fteeprom232b,
                            mfg_string::ASCIIString,
                            mfgid_string::ASCIIString,
                            d_string::ASCIIString,
                            sn_string::ASCIIString)
  @assert eepromdata.deviceType == FT_DEVICE_232BM
  fteepromprogram_alleeprom(ft_handle,eepromdata,mfg_string,
                  mfgid_string,d_string,sn_string)
end

# FT2232 EEPROM structure for use with ft_eeprom_read and ft_eeprom_program
type Fteeprom2232 <: Eeprom
  ### BEGIN common elements for all device EEPROMs ###
  deviceType :: Cuint # FT_DEVICE FTxxxx device type to be programmed
  # Device descriptor options
  VendorId :: Cshort # 0x0403
  ProductId :: Cshort # 0x6001
  SerNumEnable :: Cuchar # non-zero if serial number to be used
  # Config descriptor options
  MaxPower :: Cshort #  0 < MaxPower <= 500
  pad1::UInt8;pad2::UInt8;pad3::UInt8;pad4::UInt8  # needed to allign fields.  Why?
  ### END common elements for all device EEPROMs ###
  # Drive options
  AIsHighCurrent :: Cuchar # non-zero if interface is high current
  BIsHighCurrent :: Cuchar # non-zero if interface is high current
  # Hardware options
  AIsFifo :: Cuchar # non-zero if interface is 245 FIFO
  AIsFifoTar :: Cuchar # non-zero if interface is 245 FIFO CPU target
  AIsFastSer :: Cuchar # non-zero if interface is Fast serial
  BIsFifo :: Cuchar # non-zero if interface is 245 FIFO
  BIsFifoTar :: Cuchar # non-zero if interface is 245 FIFO CPU target
  BIsFastSer :: Cuchar # non-zero if interface is Fast serial
  # Driver option
  ADriverType :: Cuchar
  BDriverType  :: Cuchar
  Fteeprom2232() = new(FT_DEVICE_2232C)
end

function Base.show(io::IO, eepromdata::Fteeprom2232)
  showcommonelements(io,eepromdata)
  f()=@printf(io,"AIsHighCurrent= 0x%02x\n",eepromdata.AIsHighCurrent);f()
  f()=@printf(io,"BIsHighCurrent= 0x%02x\n",eepromdata.BIsHighCurrent);f()
  f()=@printf(io,"AIsFifo= 0x%02x\n",eepromdata.AIsFifo);f()
  f()=@printf(io,"AIsFifoTar= 0x%02x\n",eepromdata.AIsFifoTar);f()
  f()=@printf(io,"AIsFastSer= 0x%02x\n",eepromdata.AIsFastSer);f()
  f()=@printf(io,"BIsFifo= 0x%02x\n",eepromdata.BIsFifo);f()
  f()=@printf(io,"BIsFifoTar= 0x%02x\n",eepromdata.BIsFifoTar);f()
  f()=@printf(io,"BIsFastSer= 0x%02x\n",eepromdata.BIsFastSer);f()
  f()=@printf(io,"ADriverType= 0x%02x\n",eepromdata.ADriverType);f()
  f()=@printf(io,"BDriverType= 0x%02x\n",eepromdata.BDriverType);f()
end

function ft_eeprom_read(ft_handle::Culong, eepromdata::Fteeprom2232)
  @assert eepromdata.deviceType == FT_DEVICE_2232C
  fteepromread_alleeprom(ft_handle,eepromdata)
end

function ft_eeprom_program(ft_handle::Culong, eepromdata::Fteeprom2232,
                            mfg_string::ASCIIString,
                            mfgid_string::ASCIIString,
                            d_string::ASCIIString,
                            sn_string::ASCIIString)
  @assert eepromdata.deviceType == FT_DEVICE_2232C
  fteepromprogram_alleeprom(ft_handle,eepromdata,mfg_string,
                  mfgid_string,d_string,sn_string)
end

# FT232R EEPROM structure for use with ft_eeprom_read and ft_eeprom_program
type Fteeprom232r <: Eeprom
  ### BEGIN common elements for all device EEPROMs ###
  deviceType :: Cuint # FT_DEVICE FTxxxx device type to be programmed
  # Device descriptor options
  VendorId :: Cshort # 0x0403
  ProductId :: Cshort # 0x6001
  SerNumEnable :: Cuchar # non-zero if serial number to be used
  # Config descriptor options
  MaxPower :: Cshort #  0 < MaxPower <= 500
  pad1::UInt8;pad2::UInt8;pad3::UInt8;pad4::UInt8  # needed to allign fields.  Why?
  ### END common elements for all device EEPROMs ###
  # Drive options
  IsHighCurrent :: Cuchar # non-zero if interface is high current
  # Hardware options
  UseExtOsc :: Cuchar # Use External Oscillator
  InvertTXD :: Cuchar # non-zero if invert TXD
  InvertRXD :: Cuchar # non-zero if invert RXD
  InvertRTS :: Cuchar # non-zero if invert RTS
  InvertCTS :: Cuchar # non-zero if invert CTS
  InvertDTR :: Cuchar # non-zero if invert DTR
  InvertDSR :: Cuchar # non-zero if invert DSR
  InvertDCD :: Cuchar # non-zero if invert DCD
  InvertRI :: Cuchar # non-zero if invert RI
  Cbus0 :: Cuchar # Cbus Mux control
  Cbus1 :: Cuchar # Cbus Mux control
  Cbus2 :: Cuchar # Cbus Mux control
  Cbus3 :: Cuchar # Cbus Mux control
  Cbus4 :: Cuchar # Cbus Mux control
# Driver option
  DriverType :: Cuchar #
  Fteeprom232r() = new(FT_DEVICE_232R)
end

function Base.show(io::IO, eepromdata::Fteeprom232r)
  showcommonelements(io,eepromdata)
  f()=@printf(io,"IsHighCurrent= 0x%02x\n",eepromdata.IsHighCurrent);f()
  f()=@printf(io,"UseExtOsc= 0x%02x\n",eepromdata.UseExtOsc);f()
  f()=@printf(io,"InvertTXD= 0x%02x\n",eepromdata.InvertTXD);f()
  f()=@printf(io,"InvertRXD= 0x%02x\n",eepromdata.InvertRXD);f()
  f()=@printf(io,"InvertRTS= 0x%02x\n",eepromdata.InvertRTS);f()
  f()=@printf(io,"InvertCTS= 0x%02x\n",eepromdata.InvertCTS);f()
  f()=@printf(io,"InvertDTR= 0x%02x\n",eepromdata.InvertDTR);f()
  f()=@printf(io,"InvertDSR= 0x%02x\n",eepromdata.InvertDSR);f()
  f()=@printf(io,"InvertDCD= 0x%02x\n",eepromdata.InvertDCD);f()
  f()=@printf(io,"InvertRI= 0x%02x\n",eepromdata.InvertRI);f()
  f()=@printf(io,"Cbus0= 0x%02x\n",eepromdata.Cbus0);f()
  f()=@printf(io,"Cbus1= 0x%02x\n",eepromdata.Cbus1);f()
  f()=@printf(io,"Cbus2= 0x%02x\n",eepromdata.Cbus2);f()
  f()=@printf(io,"Cbus3= 0x%02x\n",eepromdata.Cbus3);f()
  f()=@printf(io,"Cbus4= 0x%02x\n",eepromdata.Cbus4);f()
  f()=@printf(io,"DriverType= 0x%02x\n",eepromdata.DriverType);f()
end

function ft_eeprom_read(ft_handle::Culong, eepromdata::Fteeprom232r)
  @assert eepromdata.deviceType == FT_DEVICE_232R
  fteepromread_alleeprom(ft_handle,eepromdata)
end

function ft_eeprom_program(ft_handle::Culong, eepromdata::Fteeprom232r,
                            mfg_string::ASCIIString,
                            mfgid_string::ASCIIString,
                            d_string::ASCIIString,
                            sn_string::ASCIIString)
  @assert eepromdata.deviceType == FT_DEVICE_232R
  fteepromprogram_alleeprom(ft_handle,eepromdata,mfg_string,
                  mfgid_string,d_string,sn_string)
end

# FT2232H EEPROM structure for use with ft_eeprom_read and ft_eeprom_program
type Fteeprom2232h <: Eeprom
  ### BEGIN common elements for all device EEPROMs ###
  deviceType :: Cuint # FT_DEVICE FTxxxx device type to be programmed
  # Device descriptor options
  VendorId :: Cshort # 0x0403
  ProductId :: Cshort # 0x6001
  SerNumEnable :: Cuchar # non-zero if serial number to be used
  # Config descriptor options
  MaxPower :: Cshort #  0 < MaxPower <= 500
  pad1::UInt8;pad2::UInt8;pad3::UInt8;pad4::UInt8  # needed to allign fields.  Why?
  ### END common elements for all device EEPROMs ###
  # Drive options
  ALSlowSlew :: Cuchar # non-zero if AL pins have slow slew
  ALSchmittInput :: Cuchar # non-zero if AL pins are Schmitt input
  ALDriveCurrent :: Cuchar # valid values are 4mA, 8mA, 12mA, 16mA
  AHSlowSlew :: Cuchar # non-zero if AH pins have slow slew
  AHSchmittInput :: Cuchar # non-zero if AH pins are Schmitt input
  AHDriveCurrent :: Cuchar # valid values are 4mA, 8mA, 12mA, 16mA
  BLSlowSlew :: Cuchar # non-zero if BL pins have slow slew
  BLSchmittInput :: Cuchar # non-zero if BL pins are Schmitt input
  BLDriveCurrent :: Cuchar # valid values are 4mA, 8mA, 12mA, 16mA
  BHSlowSlew :: Cuchar # non-zero if BH pins have slow slew
  BHSchmittInput :: Cuchar # non-zero if BH pins are Schmitt input
  BHDriveCurrent :: Cuchar # valid values are 4mA, 8mA, 12mA, 16mA
  # Hardware options
  AIsFifo :: Cuchar # non-zero if interface is 245 FIFO
  AIsFifoTar :: Cuchar # non-zero if interface is 245 FIFO CPU target
  AIsFastSer :: Cuchar # non-zero if interface is Fast serial
  BIsFifo :: Cuchar # non-zero if interface is 245 FIFO
  BIsFifoTar :: Cuchar # non-zero if interface is 245 FIFO CPU target
  BIsFastSer :: Cuchar # non-zero if interface is Fast serial
  PowerSaveEnable :: Cuchar # non-zero if using BCBUS7 to save power for
  # self-powered designs
  # Driver option
  ADriverType :: Cuchar #
  BDriverType :: Cuchar #
  Fteeprom2232h() = new(FT_DEVICE_2232H)
end

function Base.show(io::IO, eepromdata::Fteeprom2232h)
  showcommonelements(io,eepromdata)
  f()=@printf(io,"ALSlowSlew= 0x%02x\n",eepromdata.ALSlowSlew);f()
  f()=@printf(io,"ALSchmittInput= 0x%02x\n",eepromdata.ALSchmittInput);f()
  f()=@printf(io,"ALDriveCurrent= 0x%02x\n",eepromdata.ALDriveCurrent);f()
  f()=@printf(io,"AHSlowSlew= 0x%02x\n",eepromdata.AHSlowSlew);f()
  f()=@printf(io,"AHSchmittInput= 0x%02x\n",eepromdata.AHSchmittInput);f()
  f()=@printf(io,"AHDriveCurrent= 0x%02x\n",eepromdata.AHDriveCurrent);f()
  f()=@printf(io,"BLSlowSlew= 0x%02x\n",eepromdata.BLSlowSlew);f()
  f()=@printf(io,"BLSchmittInput= 0x%02x\n",eepromdata.BLSchmittInput);f()
  f()=@printf(io,"BLDriveCurrent= 0x%02x\n",eepromdata.BLDriveCurrent);f()
  f()=@printf(io,"BHSlowSlew= 0x%02x\n",eepromdata.BHSlowSlew);f()
  f()=@printf(io,"BHSchmittInput= 0x%02x\n",eepromdata.BHSchmittInput);f()
  f()=@printf(io,"BHDriveCurrent= 0x%02x\n",eepromdata.BHDriveCurrent);f()
  f()=@printf(io,"AIsFifo= 0x%02x\n",eepromdata.AIsFifo);f()
  f()=@printf(io,"AIsFifoTar= 0x%02x\n",eepromdata.AIsFifoTar);f()
  f()=@printf(io,"AIsFastSer= 0x%02x\n",eepromdata.AIsFastSer);f()
  f()=@printf(io,"BIsFifo= 0x%02x\n",eepromdata.BIsFifo);f()
  f()=@printf(io,"BIsFifoTar= 0x%02x\n",eepromdata.BIsFifoTar);f()
  f()=@printf(io,"BIsFastSer= 0x%02x\n",eepromdata.BIsFastSer);f()
  f()=@printf(io,"BIsFastSer= 0x%02x\n",eepromdata.BIsFastSer);f()
  f()=@printf(io,"PowerSaveEnable= 0x%02x\n",eepromdata.PowerSaveEnable);f()
  f()=@printf(io,"ADriverType= 0x%02x\n",eepromdata.ADriverType);f()
  f()=@printf(io,"BDriverType= 0x%02x\n",eepromdata.BDriverType);f()
end

function ft_eeprom_read(ft_handle::Culong, eepromdata::Fteeprom2232h)
  @assert eepromdata.deviceType == FT_DEVICE_2232H
  fteepromread_alleeprom(ft_handle,eepromdata)
end

function ft_eeprom_program(ft_handle::Culong, eepromdata::Fteeprom2232h,
                            mfg_string::ASCIIString,
                            mfgid_string::ASCIIString,
                            d_string::ASCIIString,
                            sn_string::ASCIIString)
  @assert eepromdata.deviceType == FT_DEVICE_2232H
  fteepromprogram_alleeprom(ft_handle,eepromdata,mfg_string,
                  mfgid_string,d_string,sn_string)
end

# FT4232H EEPROM structure for use with ft_eeprom_read and ft_eeprom_program
type Fteeprom4232h <: Eeprom
  ### BEGIN common elements for all device EEPROMs ###
  deviceType :: Cuint # FT_DEVICE FTxxxx device type to be programmed
  # Device descriptor options
  VendorId :: Cshort # 0x0403
  ProductId :: Cshort # 0x6001
  SerNumEnable :: Cuchar # non-zero if serial number to be used
  # Config descriptor options
  MaxPower :: Cshort #  0 < MaxPower <= 500
  pad1::UInt8;pad2::UInt8;pad3::UInt8;pad4::UInt8  # needed to allign fields.  Why?
  ### END common elements for all device EEPROMs ###
  # Drive options
  ASlowSlew :: Cuchar # non-zero if A pins have slow slew
  ASchmittInput :: Cuchar # non-zero if A pins are Schmitt input
  ADriveCurrent :: Cuchar # valid values are 4mA, 8mA, 12mA, 16mA
  BSlowSlew :: Cuchar # non-zero if B pins have slow slew
  BSchmittInput :: Cuchar # non-zero if B pins are Schmitt input
  BDriveCurrent :: Cuchar # valid values are 4mA, 8mA, 12mA, 16mA
  CSlowSlew :: Cuchar # non-zero if C pins have slow slew
  CSchmittInput :: Cuchar # non-zero if C pins are Schmitt input
  CDriveCurrent :: Cuchar # valid values are 4mA, 8mA, 12mA, 16mA
  DSlowSlew :: Cuchar # non-zero if D pins have slow slew
  DSchmittInput :: Cuchar # non-zero if D pins are Schmitt input
  DDriveCurrent :: Cuchar # valid values are 4mA, 8mA, 12mA, 16mA
  # Hardware options
  ARIIsTXDEN :: Cuchar # non-zero if port A uses RI as RS485 TXDEN
  BRIIsTXDEN :: Cuchar # non-zero if port B uses RI as RS485 TXDEN
  CRIIsTXDEN :: Cuchar # non-zero if port C uses RI as RS485 TXDEN
  DRIIsTXDEN :: Cuchar # non-zero if port D uses RI as RS485 TXDEN
  # Driver option
  ADriverType :: Cuchar #
  BDriverType :: Cuchar #
  CDriverType :: Cuchar #
  DDriverType :: Cuchar #
  Fteeprom4232h() = new(FT_DEVICE_4232H)
end

function Base.show(io::IO, eepromdata::Fteeprom4232h)
  showcommonelements(io,eepromdata)
  f()=@printf(io,"ASlowSlew= 0x%02x\n",eepromdata.ASlowSlew);f()
  f()=@printf(io,"ASchmittInput= 0x%02x\n",eepromdata.ASchmittInput);f()
  f()=@printf(io,"ADriveCurrent= 0x%02x\n",eepromdata.ADriveCurrent);f()
  f()=@printf(io,"BSlowSlew= 0x%02x\n",eepromdata.BSlowSlew);f()
  f()=@printf(io,"BSchmittInput= 0x%02x\n",eepromdata.BSchmittInput);f()
  f()=@printf(io,"BDriveCurrent= 0x%02x\n",eepromdata.BDriveCurrent);f()
  f()=@printf(io,"CSlowSlew= 0x%02x\n",eepromdata.CSlowSlew);f()
  f()=@printf(io,"CSchmittInput= 0x%02x\n",eepromdata.CSchmittInput);f()
  f()=@printf(io,"CDriveCurrent= 0x%02x\n",eepromdata.CDriveCurrent);f()
  f()=@printf(io,"DSlowSlew= 0x%02x\n",eepromdata.DSlowSlew);f()
  f()=@printf(io,"DSchmittInput= 0x%02x\n",eepromdata.DSchmittInput);f()
  f()=@printf(io,"DDriveCurrent= 0x%02x\n",eepromdata.DDriveCurrent);f()
  f()=@printf(io,"ARIIsTXDEN= 0x%02x\n",eepromdata.ARIIsTXDEN);f()
  f()=@printf(io,"BRIIsTXDEN= 0x%02x\n",eepromdata.BRIIsTXDEN);f()
  f()=@printf(io,"CRIIsTXDEN= 0x%02x\n",eepromdata.CRIIsTXDEN);f()
  f()=@printf(io,"DRIIsTXDEN= 0x%02x\n",eepromdata.DRIIsTXDEN);f()
  f()=@printf(io,"ADriverType= 0x%02x\n",eepromdata.ADriverType);f()
  f()=@printf(io,"BDriverType= 0x%02x\n",eepromdata.BDriverType);f()
  f()=@printf(io,"CDriverType= 0x%02x\n",eepromdata.CDriverType);f()
  f()=@printf(io,"DDriverType= 0x%02x\n",eepromdata.DDriverType);f()
end

function ft_eeprom_read(ft_handle::Culong, eepromdata::Fteeprom4232h)
  @assert eepromdata.deviceType == FT_DEVICE_4232H
  fteepromread_alleeprom(ft_handle,eepromdata)
end

function ft_eeprom_program(ft_handle::Culong, eepromdata::Fteeprom4232h,
                            mfg_string::ASCIIString,
                            mfgid_string::ASCIIString,
                            d_string::ASCIIString,
                            sn_string::ASCIIString)
  @assert eepromdata.deviceType == FT_DEVICE_4232H
  fteepromprogram_alleeprom(ft_handle,eepromdata,mfg_string,
                  mfgid_string,d_string,sn_string)
end

# FT232H EEPROM structure for use with ft_eeprom_read and ft_eeprom_program
type Fteeprom232h <: Eeprom
  ### BEGIN common elements for all device EEPROMs ###
  deviceType :: Cuint # FT_DEVICE FTxxxx device type to be programmed
  # Device descriptor options
  VendorId :: Cshort # 0x0403
  ProductId :: Cshort # 0x6001
  SerNumEnable :: Cuchar # non-zero if serial number to be used
  # Config descriptor options
  MaxPower :: Cshort #  0 < MaxPower <= 500
  pad1::UInt8;pad2::UInt8;pad3::UInt8;pad4::UInt8  # needed to allign fields.  Why?
  ### END common elements for all device EEPROMs ###
  # Drive options
  ACSlowSlew :: Cuchar # non-zero if AC bus pins have slow slew
  ACSchmittInput :: Cuchar # non-zero if AC bus pins are Schmitt input
  ACDriveCurrent :: Cuchar # valid values are 4mA, 8mA, 12mA, 16mA
  ADSlowSlew :: Cuchar # non-zero if AD bus pins have slow slew
  ADSchmittInput :: Cuchar # non-zero if AD bus pins are Schmitt input
  ADDriveCurrent :: Cuchar # valid values are 4mA, 8mA, 12mA, 16mA
  # CBUS options
  Cbus0 :: Cuchar # Cbus Mux control
  Cbus1 :: Cuchar # Cbus Mux control
  Cbus2 :: Cuchar # Cbus Mux control
  Cbus3 :: Cuchar # Cbus Mux control
  Cbus4 :: Cuchar # Cbus Mux control
  Cbus5 :: Cuchar # Cbus Mux control
  Cbus6 :: Cuchar # Cbus Mux control
  Cbus7 :: Cuchar # Cbus Mux control
  Cbus8 :: Cuchar # Cbus Mux control
  Cbus9 :: Cuchar # Cbus Mux control
  # FT1248 options
  FT1248Cpol :: Cuchar # FT1248 clock polarity - clock idle high (1) or clock idle  :: Cuchar #(0)
  FT1248Lsb :: Cuchar # FT1248 data is LSB (1) or MSB (0)
  FT1248FlowControl :: Cuchar # FT1248 flow control enable
  # Hardware options
  IsFifo :: Cuchar # non-zero if interface is 245 FIFO
  IsFifoTar :: Cuchar # non-zero if interface is 245 FIFO CPU target
  IsFastSer :: Cuchar # non-zero if interface is Fast serial
  IsFT1248 :: Cuchar # non-zero if interface is FT1248
  PowerSaveEnable :: Cuchar ## Driver option
  DriverType :: Cuchar #  :: Cuchar #FT X Series EEPROM structure for use with ft_eeprom_read and ft_eeprom_program
  Fteeprom232h() = new(FT_DEVICE_232H)
end

function Base.show(io::IO, eepromdata::Fteeprom232h)
  showcommonelements(io,eepromdata)
  f()=@printf(io,"ACSlowSlew = 0x%02x\n",eepromdata.ACSlowSlew);f()
  f()=@printf(io,"ACSchmittInput = 0x%02x\n",eepromdata.ACSchmittInput);f()
  f()=@printf(io,"ACDriveCurrent = 0x%02x\n",eepromdata.ACDriveCurrent);f()
  f()=@printf(io,"ADSlowSlew = 0x%02x\n",eepromdata.ADSlowSlew);f()
  f()=@printf(io,"ADSchmittInput = 0x%02x\n",eepromdata.ADSchmittInput);f()
  f()=@printf(io,"ADDriveCurrent = 0x%02x\n",eepromdata.ADDriveCurrent);f()
  f()=@printf(io,"Cbus0 = 0x%02x\n",eepromdata.Cbus0);f()
  f()=@printf(io,"Cbus1 = 0x%02x\n",eepromdata.Cbus1);f()
  f()=@printf(io,"Cbus2 = 0x%02x\n",eepromdata.Cbus2);f()
  f()=@printf(io,"Cbus3 = 0x%02x\n",eepromdata.Cbus3);f()
  f()=@printf(io,"Cbus4 = 0x%02x\n",eepromdata.Cbus4);f()
  f()=@printf(io,"Cbus5 = 0x%02x\n",eepromdata.Cbus5);f()
  f()=@printf(io,"Cbus6 = 0x%02x\n",eepromdata.Cbus6);f()
  f()=@printf(io,"Cbus7 = 0x%02x\n",eepromdata.Cbus7);f()
  f()=@printf(io,"Cbus8 = 0x%02x\n",eepromdata.Cbus8);f()
  f()=@printf(io,"Cbus9 = 0x%02x\n",eepromdata.Cbus9);f()
  f()=@printf(io,"FT1248Cpol = 0x%02x\n",eepromdata.FT1248Cpol);f()
  f()=@printf(io,"FT1248Lsb = 0x%02x\n",eepromdata.FT1248Lsb);f()
  f()=@printf(io,"FT1248FlowControl = 0x%02x\n",eepromdata.FT1248FlowControl);f()
  f()=@printf(io,"IsFifo = 0x%02x\n",eepromdata.IsFifo);f()
  f()=@printf(io,"IsFifoTar = 0x%02x\n",eepromdata.IsFifoTar);f()
  f()=@printf(io,"IsFastSer = 0x%02x\n",eepromdata.IsFastSer);f()
  f()=@printf(io,"IsFT1248 = 0x%02x\n",eepromdata.IsFT1248);f()
  f()=@printf(io,"PowerSaveEnable = 0x%02x\n",eepromdata.PowerSaveEnable);f()
  f()=@printf(io,"DriverType = 0x%02x\n",eepromdata.DriverType);f()
end

function ft_eeprom_read(ft_handle::Culong, eepromdata::Fteeprom232h)
  @assert eepromdata.deviceType == FT_DEVICE_232H
  fteepromread_alleeprom(ft_handle,eepromdata)
end

function ft_eeprom_program(ft_handle::Culong, eepromdata::Fteeprom232h,
                            mfg_string::ASCIIString,
                            mfgid_string::ASCIIString,
                            d_string::ASCIIString,
                            sn_string::ASCIIString)
  @assert eepromdata.deviceType == FT_DEVICE_232H
  fteepromprogram_alleeprom(ft_handle,eepromdata,mfg_string,
                  mfgid_string,d_string,sn_string)
end

type Fteepromxseries <: Eeprom
  ### BEGIN common elements for all device EEPROMs ###
  deviceType :: Cuint # FT_DEVICE FTxxxx device type to be programmed
  # Device descriptor options
  VendorId :: Cshort # 0x0403
  ProductId :: Cshort # 0x6001
  SerNumEnable :: Cuchar # non-zero if serial number to be used
  # Config descriptor options
  MaxPower :: Cshort #  0 < MaxPower <= 500
  pad1::UInt8;pad2::UInt8;pad3::UInt8;pad4::UInt8  # needed to allign fields.  Why?
  ### END common elements for all device EEPROMs ###
  # Drive options
  ACSlowSlew :: Cuchar # non-zero if AC bus pins have slow slew
  ACSchmittInput :: Cuchar # non-zero if AC bus pins are Schmitt input
  ACDriveCurrent :: Cuchar # valid values are 4mA, 8mA, 12mA, 16mA
  ADSlowSlew :: Cuchar # non-zero if AD bus pins have slow slew
  ADSchmittInput :: Cuchar # non-zero if AD bus pins are Schmitt input
  ADDriveCurrent :: Cuchar # valid values are 4mA, 8mA, 12mA, 16mA
  # CBUS options
  Cbus0 :: Cuchar # Cbus Mux control
  Cbus1 :: Cuchar # Cbus Mux control
  Cbus2 :: Cuchar # Cbus Mux control
  Cbus3 :: Cuchar # Cbus Mux control
  Cbus4 :: Cuchar # Cbus Mux control
  Cbus5 :: Cuchar # Cbus Mux control
  Cbus6 :: Cuchar # Cbus Mux control
  # UART signal options
  InvertTXD :: Cuchar # non-zero if invert TXD
  InvertRXD :: Cuchar # non-zero if invert RXD
  InvertRTS :: Cuchar # non-zero if invert RTS
  InvertCTS :: Cuchar # non-zero if invert CTS
  InvertDTR :: Cuchar # non-zero if invert DTR
  InvertDSR :: Cuchar # non-zero if invert DSR
  InvertDCD :: Cuchar # non-zero if invert DCD
  InvertRI :: Cuchar # non-zero if invert RI
  # Battery Charge Detect options
  BCDEnable :: Cuchar # Enable Battery Charger Detection
  BCDForceCbusPWREN :: Cuchar # asserts the power enable signal on CBUS when charging port detected
  BCDDisableSleep :: Cuchar # forces the device never to go into sleep mode
  # I2C options
  I2CSlaveAddress :: Cushort # I2C slave device address
  I2CDeviceId :: Cuint # I2C device ID
  I2CDisableSchmitt :: Cuchar # Disable I2C Schmitt trigger
  # FT1248 options
  FT1248Cpol :: Cuchar # FT1248 clock polarity - clock idle high (1) orclock idle low (0)
  FT1248Lsb :: Cuchar # FT1248 data is LSB (1) or MSB (0)
  FT1248FlowControl :: Cuchar # FT1248 flow control enable
  # Hardware options
  RS485EchoSuppress :: Cuchar
  PowerSaveEnable :: Cuchar ## Driver option
  DriverType :: Cuchar
  Fteepromxseries() = new(FT_DEVICE_X_SERIES)
end

function Base.show(io::IO, eepromdata::Fteepromxseries)
  showcommonelements(io,eepromdata)
  f()=@printf(io,"ACSlowSlew= 0x%02x\n",eepromdata.ACSlowSlew);f()
  f()=@printf(io,"ACSchmittInput= 0x%02x\n",eepromdata.ACSchmittInput);f()
  f()=@printf(io,"ACDriveCurrent= 0x%02x\n",eepromdata.ACDriveCurrent);f()
  f()=@printf(io,"ADSlowSlew= 0x%02x\n",eepromdata.ADSlowSlew);f()
  f()=@printf(io,"ADSchmittInput= 0x%02x\n",eepromdata.ADSchmittInput);f()
  f()=@printf(io,"ADDriveCurrent= 0x%02x\n",eepromdata.ADDriveCurrent);f()
  f()=@printf(io,"Cbus0= 0x%02x\n",eepromdata.Cbus0);f()
  f()=@printf(io,"Cbus1= 0x%02x\n",eepromdata.Cbus1);f()
  f()=@printf(io,"Cbus2= 0x%02x\n",eepromdata.Cbus2);f()
  f()=@printf(io,"Cbus3= 0x%02x\n",eepromdata.Cbus3);f()
  f()=@printf(io,"Cbus4= 0x%02x\n",eepromdata.Cbus4);f()
  f()=@printf(io,"Cbus5= 0x%02x\n",eepromdata.Cbus5);f()
  f()=@printf(io,"Cbus6= 0x%02x\n",eepromdata.Cbus6);f()
  f()=@printf(io,"InvertTXD= 0x%02x\n",eepromdata.InvertTXD);f()
  f()=@printf(io,"InvertRXD= 0x%02x\n",eepromdata.InvertRXD);f()
  f()=@printf(io,"InvertRTS= 0x%02x\n",eepromdata.InvertRTS);f()
  f()=@printf(io,"InvertCTS= 0x%02x\n",eepromdata.InvertCTS);f()
  f()=@printf(io,"InvertDTR= 0x%02x\n",eepromdata.InvertDTR);f()
  f()=@printf(io,"InvertDSR= 0x%02x\n",eepromdata.InvertDSR);f()
  f()=@printf(io,"InvertDCD= 0x%02x\n",eepromdata.InvertDCD);f()
  f()=@printf(io,"InvertRI= 0x%02x\n",eepromdata.InvertRI);f()
  f()=@printf(io,"BCDEnable= 0x%02x\n",eepromdata.BCDEnable);f()
  f()=@printf(io,"BCDForceCbusPWREN= 0x%02x\n",eepromdata.BCDForceCbusPWREN);f()
  f()=@printf(io,"BCDDisableSleep= 0x%02x\n",eepromdata.BCDDisableSleep);f()
  f()=@printf(io,"I2CSlaveAddress= 0x%02x\n",eepromdata.I2CSlaveAddress);f()
  f()=@printf(io,"I2CDeviceId= 0x%02x\n",eepromdata.I2CDeviceId);f()
  f()=@printf(io,"I2CDisableSchmitt= 0x%02x\n",eepromdata.I2CDisableSchmitt);f()
  f()=@printf(io,"FT1248Cpol= 0x%02x\n",eepromdata.FT1248Cpol);f()
  f()=@printf(io,"FT1248Lsb= 0x%02x\n",eepromdata.FT1248Lsb);f()
  f()=@printf(io,"FT1248FlowControl= 0x%02x\n",eepromdata.FT1248FlowControl);f()
  f()=@printf(io,"RS485EchoSuppress= 0x%02x\n",eepromdata.RS485EchoSuppress);f()
  f()=@printf(io,"PowerSaveEnable= 0x%02x\n",eepromdata.PowerSaveEnable);f()
  f()=@printf(io,"DriverType= 0x%02x\n",eepromdata.DriverType);f()
end

function ft_eeprom_read(ft_handle::Culong, eepromdata::Fteepromxseries)
  @assert eepromdata.deviceType == FT_DEVICE_X_SERIES
  fteepromread_alleeprom(ft_handle,eepromdata)
end

function ft_eeprom_program(ft_handle::Culong, eepromdata::Fteepromxseries,
                            mfg_string::ASCIIString,
                            mfgid_string::ASCIIString,
                            d_string::ASCIIString,
                            sn_string::ASCIIString)
  @assert eepromdata.deviceType == FT_DEVICE_X_SERIES
  fteepromprogram_alleeprom(ft_handle,eepromdata,mfg_string,
                  mfgid_string,d_string,sn_string)
end

function ft_eeprom_read(ft_handle::Culong)
  (devicetype,id,sn,d) = ft_getdeviceinfo(ft_handle)
  if devicetype == FT_DEVICE_232BM
    return ft_eeprom_read(ft_handle,Fteeprom232b())
  elseif devicetype == FT_DEVICE_2232C
    return ft_eeprom_read(ft_handle,Fteeprom2232())
  elseif devicetype == FT_DEVICE_232R
    return ft_eeprom_read(ft_handle,Fteeprom232r())
  elseif devicetype == FT_DEVICE_2232H
    return ft_eeprom_read(ft_handle,Fteeprom2232h())
  elseif devicetype == FT_DEVICE_4232H
    return ft_eeprom_read(ft_handle,Fteeprom4232h())
  elseif devicetype == FT_DEVICE_232H
    return ft_eeprom_read(ft_handle,Fteeprom232h())
  elseif devicetype == FT_DEVICE_X_SERIES
    return ft_eeprom_read(ft_handle,Fteepromxseries())
  else
    error("Device type ",devicetype," unknown")
  end
end


text:
@testset "NEQ Kirschmer" begin
  R,x = FlintZZ["x"]
  K,a = number_field(x^2+3x+1,"a")
  kt,t = K["t"]
  E, b = number_field( t^2 + (a+5)* (-2*a+2) )
  @test Hecke.isnorm( E, -a*(-2*a+2) )[1]
end

ext: jl


text:
#! /usr/bin/julia

# Rosetta Code, Combinations and permutations

function Base.binomial{T<:FloatingPoint}(n::T, k::T)
    exp(lfact(n) - lfact(n - k) - lfact(k))
end

function Base.factorial{T<:FloatingPoint}(n::T, k::T)
    exp(lfact(n) - lfact(k))
end

⊞{T<:Real}(n::T, k::T) = binomial(n, k)
⊠{T<:Real}(n::T, k::T) = factorial(n, n-k)

function picknk{T<:Integer}(lo::T, hi::T)
    n = rand(lo:hi)
    k = rand(1:n)
    return (n, k)
end

nsamp = 10

print("Tests of the combinations (⊞) and permutations (⊠) operators for ")
println("integer values.")
println()
lo, hi = 1, 12
print(nsamp, " samples of n ⊠ k with n in [", lo, ", ", hi, "] ")
println("and k in [1, n].")
for i in 1:nsamp
    (n, k) = picknk(lo, hi)
    println(@sprintf "    %2d ⊠ %2d = %18d" n k n ⊠ k)
end

lo, hi = 10, 60
println()
print(nsamp, " samples of n ⊞ k with n in [", lo, ", ", hi, "] ")
println("and k in [1, n].")
for i in 1:nsamp
    (n, k) = picknk(lo, hi)
    println(@sprintf "    %2d ⊞ %2d = %18d" n k n ⊞ k)
end

println()
print("Tests of the combinations (⊞) and permutations (⊠) operators for ")
println("(big) float values.")
println()
lo, hi = 5, 15000
print(nsamp, " samples of n ⊠ k with n in [", lo, ", ", hi, "] ")
println("and k in [1, n].")
for i in 1:nsamp
    (n, k) = picknk(lo, hi)
    n = BigFloat(n)
    k = BigFloat(k)
    println(@sprintf "    %7.1f ⊠ %7.1f = %10.6e" n k n ⊠ k)
end

lo, hi = 100, 1000
println()
print(nsamp, " samples of n ⊞ k with n in [", lo, ", ", hi, "] ")
println("and k in [1, n].")
for i in 1:nsamp
    (n, k) = picknk(lo, hi)
    n = BigFloat(n)
    k = BigFloat(k)
    println(@sprintf "    %7.1f ⊞ %7.1f = %10.6e" n k n ⊞ k)
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
module GPARatScale
    using Stheno, TemporalGPs, Distributions, Optim, Zygote, Random
    # Data generation
    include(joinpath("data", "toy_data.jl"))

    # GPAR
    include(joinpath("gp", "optimized.jl"))
    # Temporal GP
    include(joinpath("gp", "temporal_gp_inference.jl"))
    # Util
    include("util.jl")
end

max_stars_repo_name: willtebbutt/GPAR-at-scale

max_issues_repo_name: willtebbutt/GPAR-at-scale

max_forks_repo_name: willtebbutt/GPAR-at-scale


text:
######################################################################################################
# Merged the implementation of Wouter Koolen and Xuedong Shang                                       #
# * https://bitbucket.org/wmkoolen/tidnabbil/src/master/purex_games_paper/                           #
# * https://github.com/xuedong/LinBAI.jl/blob/master/linear/pep_linear.jl                            #
# For the function `glrt` for linear and lipschitz structures, we use functions in                   #
# `utilities/envelope.jl` to implement, which is based on the Envelope Theorem.                      #
######################################################################################################
# A Pure Exploration problem (pep) is parameterised by
# - a domain, \mathcal M, representing the prior knowledge about the structure
# - a query, as embodied by a correct-answer function istar
# To specify a PE problem here, we need to compute the following things:
# - nanswers: number of possible answers
# - istar: correct answer for feasible μ
# - glrt: value and best response (λ and ξ) to (N, ̂μ) or (w, μ)
# - oracle: characteristic time and oracle weights at μ
using IterTools;
include("envelope.jl"); # an elegant way to compute alt_min for Lipschitz structure

"""
Classical.
All functions are implemented by Wouter Koolen.
Note that we can also use Envelope Theorem to simplify `glrt` and to implement `oracle`.
"""
struct BestArm
    expfam; # common exponential family
end
nanswers(pep::BestArm, μ) = length(μ);
istar(pep::BestArm, μ) = argmax(μ);
getexpfam(pep::BestArm, k) = pep.expfam;
function glrt(pep::BestArm, w, μ)
    @assert length(size(μ)) == 1
    ⋆ = argmax(μ); # index of best arm among μ
    val, k, θ = minimum(
        begin
        # transport ⋆ and k to weighted midpoint
        θ = (w[⋆]*μ[⋆]+w[k]*μ[k])/(w[⋆]+w[k]);
        w[⋆]*d(pep.expfam, μ[⋆], θ) + w[k]*d(pep.expfam, μ[k], θ), k, θ
        end
        for k in eachindex(μ)
        if k != ⋆
    );
    λ = copy(μ);
    λ[⋆] = θ;
    λ[k] = θ;
    val, (k, λ), (⋆, μ);
end
# solve for x such that d(μ1, μx) + x*d(μi, μx) == v
# where μx = (μ1+x*μi)/(1+x)
function X(expfam, μ1, μi, v)
    kl1i = d(expfam, μ1, μi); # range of V(x) is [0, kl1i]
    @assert 0 ≤ v ≤ kl1i "0 ≤ $v ≤ $kl1i";
    α = binary_search(
        z -> let μz = (1-z)*μ1+z*μi
        (1-z)*d(expfam, μ1, μz) + z*d(expfam, μi, μz) - (1-z)*v
        end,
        0, 1, ϵ = kl1i*1e-10);
    α/(1-α), (1-α)*μ1+α*μi
end
# oracle problem
function oracle(pep, μs)
    μstar = maximum(μs);
    if all(μs .== μstar) # yes, this happens
        return Inf, ones(length(μs))/length(μs);
    end
    # determine upper range for subsequent binary search
    hi = minimum(
        d(pep.expfam, μstar, μ)
        for μ in μs
        if μ != μstar
    );
    val = binary_search(
        z -> sum(
            let μx = X(pep.expfam, μstar, μ, z)[2];
            d(pep.expfam, μstar, μx) / d(pep.expfam, μ, μx)
            end
            for μ in μs
            if μ != μstar
            ) - 1,
        0, hi);
    ws = [(μ == μstar) ? 1. : X(pep.expfam, μstar, μ, val)[1] for μ in μs];
    Σ = sum(ws);
    Σ/val, ws ./ Σ;
end
# oracle problem for best μ in confidence interval
function optimistic_oracle(pep::BestArm, hμ, N)
    t = sum(N);
    μdn = [ddn(pep.expfam, hμ[k], log(t)/N[k]) for k in eachindex(hμ)];
    μup = [dup(pep.expfam, hμ[k], log(t)/N[k]) for k in eachindex(hμ)];
    # try to make each arm the best-looking arm so far, then move everybody else down as far as possible
    minimum(oracle(pep, ((k == j) ? μup[k] : μdn[k] for k in eachindex(hμ))) for j in eachindex(hμ) if μup[j] > maximum(μdn));
end

"""
Linear Structure.
Most functions are implemented by Xuedong Shang.
We use Envelope Theorem to simplify `glrt` and to implement `oracle_linear_bai` and `oracle_linear_threshold`.
"""
function sherman_morrison(Vinv, u, v)
    num = (Vinv*u)*transpose(transpose(Vinv)*v)
    denum = 1 + transpose(v)*Vinv*u
    return Vinv .- num / denum
end
function sherman_morrison(Vinv, u)
    Vinv_u = Vinv*u
    num = Vinv_u*transpose(Vinv_u)
    denum = 1 + transpose(u)*Vinv_u
    return Vinv .- num / denum
end
function build_gaps(arms)
    gaps = Vector{Float64}[]
    for pair in subsets(arms, 2)
        gap1 = pair[1] - pair[2]
        push!(gaps, gap1)
        gap2 = pair[2] - pair[1]
        push!(gaps, gap2)
        #@show pair
    end
    return gaps
end

struct LinearBestArm
    expfam; # common exponential family
    arms;  # array of size K of arms in R^d
end
nanswers(pep::LinearBestArm, μ) = length(pep.arms);
narms(pep::LinearBestArm, μ) = length(pep.arms);
istar(pep::LinearBestArm, μ) = argmax([arm'μ for arm in pep.arms]);
armstar(pep::LinearBestArm, μ) = pep.arms[argmax([arm'μ for arm in pep.arms])];
getexpfam(pep::LinearBestArm, k) = pep.expfam;
############################################ Required by CG-C and Lk-C implemented by Xuedong Shang ############################################
function alt_min_cgc(pep::LinearBestArm, w, µ, k) # by Xuedong Shang
    sum_w = sum(w)
    w = w/sum(w)  # avoid dividing by small quantities
    arm = pep.arms[k]
    arm_star = armstar(pep, µ)
    @assert arm != arm_star
    K = length(pep.arms)
    direction = arm .- arm_star
    d = length(direction)
    # Construction of the matrix V = (sum_k w_k a_k a_k^T)^{-1}
    sum_arms_matrix = zeros(d,d)
    for j in 1:K
        sum_arms_matrix .+= w[j].* (pep.arms[j]*transpose(pep.arms[j]))
    end
    Vinv = inv(sum_arms_matrix)
    # Closest point
    η = sum_w * (direction'µ) / ((direction')*Vinv*direction)
    λ = µ .- η/sum_w * Vinv * direction
    # Divergence to that point
    val = .5 * sum_w * (direction'µ)^2 / ((direction')*Vinv*direction)
    return val, λ, k
end
function alt_min_cgc(pep::LinearBestArm, w, µ) # by Xuedong Shang
    minimum(alt_min_cgc(pep, w, µ, i) for i in 1:nanswers(pep, µ) if i != istar(pep, µ))
end
function glrt_cgc(pep::LinearBestArm, w::Vector, μ) # by Xuedong Shang
    @assert length(size(μ)) == 1
    star = istar(pep, µ)
    val, λ, k = alt_min_cgc(pep, w, µ)
    #Return: distance to closest alternative, best arm and closest λ for that alternative, best arm and vector for closest point in model
    val, (k, λ), (star, μ);
end
function glrt_cgc(pep::LinearBestArm, P::Matrix, μ) # by Xuedong Shang
    # P is a matrix of size nb_answers*nb_arms (both of which may be different from length(µ)).
    # In the BAI case, nb_answers = nb_arms != length(µ)
    nb_answers = size(P)[1]
    @assert length(size(μ)) == 1
    val = 0
    λs = zeros(nb_answers, length(µ))  # nb_answers * length(µ0)
    star = istar(pep, µ)
    for i in 1:nb_answers
        if i != star
            λs[i, :] = copy(µ)  # µ belongs to ¬i
        else
            val_i, λ_i, _ = alt_min_cgc(pep, P[i,:], µ)
            #println("i $i ; val_i $val_i ; λ_i $λ_i")
            λs[i, :] = λ_i
            val += val_i
        end
    end
    val, λs, (star, μ);
end
################################################################################################################################################
############################################## Our elegant implementation by Envelope Theorem ##################################################
function glrt(pep::LinearBestArm, N, μ, Vxinv)
    @assert length(size(μ)) == 1
    ⋆ = argmax(μ); # index of best arm among μ
    val, _, _, _ = alt_min_linear_bai(N, μ, pep.arms, Vxinv);
    return val, ⋆, μ;
end
function oracle_linear_bai(pep::LinearBestArm, μ)
    K = length(pep.arms); dim = length(μ); reward = [μ'pep.arms[k] for k=1:K]; μi = argmax(reward);
    @assert all(hμ_in_lambda(reward, μi, K)) "μ'arm has >=2 best arms.";
    x = [1.0/K for k=1:K]; Vxinv = zeros(Float64, dim,dim);
    for k in 1:K
        Vxinv += pep.arms[k]*(pep.arms[k]')/K;
    end
    for t=1:10000
        if is_complete_square(floor(Int, t/K))
            z = [1.0/K for k=1:K];
        else
            f, ∇f, fidx = compute_f_∇f_linear_bai(x, μ, t^(-0.9)/K, pep.arms, Vxinv);
            Σ = [[(i==j) ? 1 : 0 for j=1:K]-x for i=1:K];
            A = [[Σ[i]'∇f[j] for i=1:K] for j in fidx]; # construct payoff matrix
            z = solveZeroSumGame(A, K, length(fidx));
        end
        x = t/(t+1.0) * x + z/(t+1.0);
        for k=1:K
            Vxinv = sherman_morrison(Vxinv*t/(t+1.0), z[k]*pep.arms[k]/(t+1.0));
        end
    end
    value, _, _, _ = alt_min_linear_bai(x, μ, pep.arms, Vxinv);
    return 1/value, x;
end
#########################################################################################################################################
struct LinearThreshold
    expfam; # common exponential family
    arms;  # array of size K of arms in R^d
    threshold;  # threshold
end
nanswers(pep::LinearThreshold, μ) = 2^length(µ);
narms(pep::LinearThreshold, µ) = length(pep.arms);
istar(pep::LinearThreshold, μ) = 1 + sum([(µ[k] > pep.threshold) * 2^(k-1) for k in 1:length(µ)]);
getexpfam(pep::LinearThreshold, k) = pep.expfam;
############################################ Required by CG-C and Lk-C implemented by Xuedong Shang ############################################
function alt_min_cgc(pep::LinearThreshold, w, µ, k)
    sum_w = sum(w)
    w = w/sum(w)  # avoid dividing by small quantities
    K = narms(pep, µ)
    # Construction of the matrix V = (sum_k w_k a_k a_k^T)^{-1}
    sum_arms_matrix = zeros(length(µ),length(µ))
    for j in 1:K
        sum_arms_matrix .+= w[j].* (pep.arms[j]*transpose(pep.arms[j]))
    end
    Vinv = inv(sum_arms_matrix)
    # Closest point
    η = (µ[k] - pep.threshold) / Vinv[k,k]
    λ = µ .- η * Vinv[:,k]
    # Divergence to that point
    val = .5 * sum_w * (µ[k] - pep.threshold)^2 / Vinv[k,k]
    return val, λ, k
end
function alt_min_cgc(pep::LinearThreshold, w, µ)
    minimum(alt_min_cgc(pep, w, µ, i) for i in 1:length(µ))
end
function glrt_cgc(pep::LinearThreshold, w::Vector, μ)
    @assert length(size(μ)) == 1
    star = istar(pep, µ)
    val, λ, k = alt_min_cgc(pep, w, µ)
    answer = star
    if µ[k]>pep.threshold
        answer -= 2^(k-1)
    else
        answer += 2^(k-1)
    end
    #Return: distance to closest alternative, answer and closest λ for that alternative, answer and vector for closest point in model
    val, (answer, λ), (star, μ);
end
function glrt_cgc(pep::LinearThreshold, P::Matrix, μ)
    # P is a matrix of size nb_answers*nb_arms (both of which may be different from length(µ)).
    # In the BAI case, nb_answers = nb_arms != length(µ)
    nb_answers = nanswers(pep, µ)
    @assert length(size(μ)) == 1
    val = 0
    λs = zeros(nb_answers, length(µ))
    star = istar(pep, µ)
    for i in 1:nb_answers
        if i != star
            λs[i, :] = copy(µ)  # µ belongs to ¬i
        else
            val_i, λ_i, _ = alt_min_cgc(pep, P[i,:], µ)
            λs[i, :] = λ_i
            val += val_i
        end
    end
    val, λs, (star, μ);
end
################################################################################################################################################
############################################## Our elegant implementation by Envelope Theorem ##################################################
function glrt(pep::LinearThreshold, N, μ, Vxinv)
    @assert length(size(μ)) == 1
    answer = [((pep.arms[k]'μ>pep.threshold) ? 1 : 0) for k=1:length(pep.arms)];
    val, _, _ = alt_min_linear_threshold(N, μ, pep.arms, Vxinv);
    return val, answer, μ;
end
function oracle_linear_threshold(pep::LinearThreshold, μ)
    K = length(pep.arms); dim = length(μ); reward = [μ'pep.arms[k] for k=1:K]; μi = argmax(reward);
    @assert all(hμ_in_lambda_threshold(reward, K, pep.threshold)) "the difference between μ'arm and threshold should be greater than $(eps())";
    x = [1.0/K for k=1:K]; Vxinv = zeros(Float64, dim,dim);
    for k in 1:K
        Vxinv += pep.arms[k]*(pep.arms[k]')/K;
    end
    for t=1:10000
        if is_complete_square(floor(Int, t/K))
            z = [1.0/K for k=1:K];
        else
            f, ∇f, fidx = compute_f_∇f_linear_threshold(pep, x, μ, t^(-0.9)/K, Vxinv);
            Σ = [[(i==j) ? 1 : 0 for j=1:K]-x for i=1:K];
            A = [[Σ[i]'∇f[j] for i=1:K] for j in fidx]; # construct payoff matrix
            z = solveZeroSumGame(A, K, length(fidx));
        end
        x = t/(t+1.0) * x + z/(t+1.0);
        for k=1:K
            Vxinv = sherman_morrison(Vxinv*t/(t+1.0), z[k]*pep.arms[k]/(t+1.0));
        end
    end
    value, _, _ = alt_min_linear_threshold(x, μ, pep.arms, Vxinv);
    return 1/value, x;
end
################################################################################################################################################

"""
Lipschitz Structure.
We use Envelope Theorem to simplify `glrt` function and to implement `oracle_lipschitz_bai`.
"""
struct LipschitzBestArm
    expfam; # common exponential family
    arms;  # array of size K of arms
    μ; # reward of the arms
    L; # lipschitz constant
end
nanswers(pep::LipschitzBestArm, μ) = length(μ);
istar(pep::LipschitzBestArm, μ) = argmax(μ);
getexpfam(pep::LipschitzBestArm, k) = pep.expfam;
getlipschconst(pep::LipschitzBestArm) = pep.L;
############################################## Our elegant implementation by Envelope Theorem ##################################################
function glrt(pep::LipschitzBestArm, N, μ, L)
    @assert length(size(μ)) == 1
    ⋆ = argmax(μ); # index of best arm among μ
    val, _, _, _ = alt_min_lipschitz(N, μ, pep.arms, L);
    val, ⋆, μ;
end
function oracle_lipschitz_bai(pep::LipschitzBestArm)
    K = length(pep.arms); μi = argmax(pep.μ);
    @assert all(hμ_in_lambda(pep.μ, μi, K) && is_lipschitz(pep.μ, pep.arms, pep.L)) "μ has >=2 best arms or μ violates $(pep.L)-Lipschitz assumption.";
    x = [1.0/K for k=1:K];
    for t=1:100000 # depending on the number of arms
        if is_complete_square(floor(Int, t/K))
            z = [1.0/K for k=1:K];
        else
            f, ∇f, fidx = compute_f_∇f_lipschitz_bai(x, pep.μ, t^(-0.9)/K, pep.arms, pep.L);
            Σ = [[(i==j) ? 1 : 0 for j=1:K]-x for i=1:K];
            A = [[Σ[i]'∇f[j] for i=1:K] for j in fidx]; # construct payoff matrix
            z = solveZeroSumGame(A, K, length(fidx));
        end
        x = t/(t+1.0) * x + z/(t+1.0);
    end
    value, _, _, _ = alt_min_lipschitz(x, pep.μ, pep.arms, pep.L);
    return 1/value, x;
end
#########################################################################################################################################


text:
"""
`Basics.excludeDoubles()  ... exclude 'double entries' from various list; see below:

  + `(confList::Array{Configuration,1})`
    ... to exlude from the (non-relativistic) confList all 'doubles', ie. configurations that occured before in the list;
        a new confList::Array{Configuration,1} is returned.

  + `(csfList::Array{CsfR,1})`
    ... to exlude from the (relativistic) csfList all 'doubles', ie. CSF that occured already before in the list;
        a new csfList::Array{CsfR,1} is returned.
"""
function Basics.excludeDoubles(confList::Array{Configuration,1})
    confListNew = [ deepcopy(confList[1]) ]

    for  ic  in 2:length(confList)
        add = true
        for  id in 1:ic-1    if   confList[ic] == confList[id]    add = false;   break   end    end
        if  add   push!( confListNew, confList[ic])   end
    end

    return( confListNew )
end

function Basics.excludeDoubles(confList::Array{CsfR,1})
    error("Not yet implemented !")

    csfListNew = [ deepcopy(csfList[1]) ]

    for  ic  in 2:length(csfList)
        add = true
    	  for  id in 1:ic-1
    		if   csfList[ic] == csfList[id]    add = false;   break   end
    	  end
    	  if  add   push!( csfListNew, csfList[ic])	end
    end
    return( csfListNew )
end


text:
models = (
    MultiTaskLassoRegressor,
    MultiTaskLassoCVRegressor,
    MultiTaskElasticNetRegressor,
    MultiTaskElasticNetCVRegressor
    )

fparams = (
    MultiTaskLassoRegressor=(:coef, :intercept),
    MultiTaskLassoCVRegressor=(:coef, :intercept, :alpha, :mse_path, :alphas),
    MultiTaskElasticNetRegressor=(:coef, :intercept),
    MultiTaskElasticNetCVRegressor=(:coef, :intercept, :alpha, :mse_path, :l1_ratio)
    )

@testset "Fit/Predict" begin
    # check roughly that models return sensible results
    X, y, ls = simple_regression_2()
    for mod in models
        m = mod()
        # ratio of rmse to LS
        r, f = test_regression_2(m, X, y, ls)
        @test r < 200
        @test keys(MB.fitted_params(m, f)) == getproperty(fparams, Symbol(mod))
    end
end

@testset "Docstrings" begin
    for mod in models
        m = mod()
        @test !isempty(MB.docstring(m))
    end
end


text:
# ---
# title: 869. Reordered Power of 2
# id: problem869
# author: Indigo
# date: 2021-07-11
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/reordered-power-of-2/description/>
# hidden: true
# ---
#
# Starting with a positive integer `N`, we reorder the digits in any order
# (including the original order) such that the leading digit is not zero.
#
# Return `true` if and only if we can do this in a way such that the resulting
# number is a power of 2.
#
#
#
# **Example 1:**
#
#
#
#     Input: 1
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: 10
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: 16
#     Output: true
#
#
# **Example 4:**
#
#
#
#     Input: 24
#     Output: false
#
#
# **Example 5:**
#
#
#
#     Input: 46
#     Output: true
#
#
#
#
# **Note:**
#
#   1. `1 <= N <= 10^9`
#
#
## @lc code=start
using LeetCode

function reordered_power_of2(n::Integer)
    cntr1 = OffsetArray(fill(0, 10), -1)
    cntr2 = OffsetArray(fill(0, 10), -1)
    function count_num!(cntr, n::Integer)
        fill!(cntr, 0)
        for d in digits(n)
            cntr[d] += 1
        end
        cntr
    end
    count_num!(cntr1, n)
    for i in 0:31
        (cntr1 == count_num!(cntr2, 1 << i)) && return true
    end
    return false
end
## @lc code=end


text:
#=
## Integration with `NetworkDynamics.jl`
In this example we model a Network based on the [Kuroamoto model](https://en.wikipedia.org/wiki/Kuramoto_model) using the [`NetworkDynamics.jl`](https://github.com/FHell/NetworkDynamics.jl) package.
=#

using NetworkDynamics
using Graphs
using BlockSystems
using ModelingToolkit
using OrdinaryDiffEq
using Plots

#=
In the Kuramoto model each vertex $i$ has it's own intrinsic angular frequency $\omega_i$.
The angle of frequency is given by
```math
\dot\phi = \omega_i + \sum e
```
where $\sum e$ sums over all the connected edges.
=#
@parameters t ω edgesum(t)
@variables ϕ(t)
D = Differential(t)

kmvert = IOBlock([D(ϕ) ~ ω + edgesum],
                 [edgesum],
                 [ϕ],
                 name=:kmvert)

function generate_ode_vertex(iob::IOBlock, aggregator; kwargs...)
    gen = generate_io_function(iob; kwargs...)
    f = (du, u, edges, p, t) -> gen.f_ip(du, u, aggregator(edges), p, t)
    vars = ModelingToolkit.tosymbol.(gen.states)
    ODEVertex(f = f, dim = length(vars), sym = vars)
end

## allocation free oop aggregator. might be more difficult for more-dimensional systems
## unfortunately there are no vector variables in MDK and we can't model the aggregators
## as an IOSystem.
function sumedges(edges)
    r = 0.0
    for e in edges
        r += e[1]
    end
    return r
end

vertex = generate_ode_vertex(kmvert, sumedges)
nothing #hide

#=
The edges are defined as
```math
e_{i,j} = K \cdot \mathrm{sin}(\phi_i - \phi_j)
```
For the (static) edges we can't used IOSystem function building in the moment.
Well, we could but than we'd get an `ODEEdge` with zero-massmatrix.
=#
@parameters v₁(t) v₂(t) K
edge_ip = build_function([K*sin(v₁[1] - v₂[1])], [v₁], [v₂], K, t, expression=Val{false})[2]
edge = StaticEdge(f = edge_ip, dim = 1, coupling=:antisymmetric)

nothing #hide
#=
Lets runt the [same example](https://fhell.github.io/NetworkDynamics.jl/dev/heterogeneous_system/#Heterogenous-parameters) as in `NetworkDynamics.jl`
=#
N = 8
g = watts_strogatz(N,2,0) # ring network
nd = network_dynamics(vertex, edge, g)

ω = collect(1:N)./N
ω .-= sum(ω)/N
K = 3.
p = (ω, K); # p[1] vertex parameters, p[2] edge parameters

x0 = collect(1:N)./N
x0 .-= sum(x0)./N

#src using BenchmarkTools
#src @btime $nd($x0, $x0, $p, 0.0)

tspan = (0., 10.)
prob = ODEProblem(nd, x0, tspan, p)
sol = solve(prob, Tsit5())
plot(sol, ylabel="ϕ")


text:
module BrokemiaHelperCaveWall

using ..Ahorn, Maple

@mapdef Entity "BrokemiaHelper/caveWall" CaveWall(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight,
    tiletype::String="3", disableTransitionFading::Bool=false)

const placements = Ahorn.PlacementDict(
    "Cave Wall (BrokemiaHelper)" => Ahorn.EntityPlacement(
        CaveWall,
        "rectangle",
        Dict{String, Any}(),
        Ahorn.tileEntityFinalizer
    ),
)

Ahorn.editingOptions(entity::CaveWall) = Dict{String, Any}(
    "tiletype" => Ahorn.tiletypeEditingOptions()
)

Ahorn.minimumSize(entity::CaveWall) = 8, 8
Ahorn.resizable(entity::CaveWall) = true, true

Ahorn.selection(entity::CaveWall) = Ahorn.getEntityRectangle(entity)

Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::CaveWall, room::Maple.Room) = Ahorn.drawTileEntity(ctx, room, entity, alpha=0.7)

end

max_stars_repo_name: Brokemia/BrokemiaHelper

max_issues_repo_name: Brokemia/BrokemiaHelper

max_forks_repo_name: Brokemia/BrokemiaHelper


text:
mutable struct FASTQRecord
   sequence::BioSequence{DNAAlphabet{2}}
   quality::Vector{UInt8}
   raw::FASTQ.Record

   FASTQRecord() = new( BioSequence{DNAAlphabet{2}}(),
                        Vector{UInt8}(), FASTQ.Record() )
end

Base.read!( parser::FASTQ.Reader, rec::FASTQRecord ) = read!( parser, rec.raw )

function Base.fill!( rec::FASTQRecord, offset=33 )
   rec.sequence = BioSequences.BioSequence{BioSequences.DNAAlphabet{2}}( rec.raw.data,
                                                                         first(rec.raw.sequence),
                                                                         last(rec.raw.sequence) )
   rec.quality  = rec.raw.data[rec.raw.quality] - convert(UInt8, offset)
   rec
end

Base.string( rec::FASTQRecord ) = String(rec.raw.data[rec.raw.identifier])

function reverse_complement!( rec::FASTQRecord )
   reverse_complement!( rec.sequence )
   reverse!( rec.quality )
end


text:
### Before parsing

unsafeattribute(::Nothing, ::Core.AbstractString) = nothing
unsafeattribute(x::LightXML.XMLElement, name::Core.AbstractString) = attribute(x, name)

function parse_scalar(xel, name::String, T; default::Union{String,Nothing}=nothing)
    scalarstr = unsafeattribute(xel, name)
    if scalarstr === nothing
        if default === nothing
            @error "no parsable scalar found"
        else
            scalarstr = default
        end
    end

    return parse(T, scalarstr)
end

function parse_vector(xel, name::String, T; default::Union{String,Nothing}=nothing)
    vectorstr = unsafeattribute(xel, name)
    if vectorstr === nothing
        if default === nothing
            @error "no parsable vector found"
        else
            vectorstr = default
        end
    end

    return parse.(T,split(vectorstr))
end

function parse_inertiamatrix(xinertia, T)
    if xinertia === nothing
        J = zeros(T, 3, 3)
    else
        ixx = parse_scalar(xinertia, "ixx", T)
        ixy = parse_scalar(xinertia, "ixy", T, default = "0")
        ixz = parse_scalar(xinertia, "ixz", T, default = "0")
        iyy = parse_scalar(xinertia, "iyy", T)
        iyz = parse_scalar(xinertia, "iyz", T, default = "0")
        izz = parse_scalar(xinertia, "izz", T)
        J = [ixx ixy ixz; ixy iyy iyz; ixz iyz izz]
    end

    return J
end

function parse_pose(xpose, T)
    if xpose === nothing
        x, q = zeros(T, 3), one(UnitQuaternion{T})
    else
        x = parse_vector(xpose, "xyz", T, default = "0 0 0")
        rpy = parse_vector(xpose, "rpy", T, default = "0 0 0")
        q = UnitQuaternion(RotZYX(rpy[3], rpy[2], rpy[1]))
    end

    return x, q
end

function parse_inertia(xinertial, T)
    if xinertial === nothing
        x = zeros(T, 3)
        q = one(UnitQuaternion{T})
        m = zero(T)
        J = zeros(T, 3, 3)
    else
        x, q = parse_pose(find_element(xinertial, "origin"), T)
        J = parse_inertiamatrix(find_element(xinertial, "inertia"), T)
        m = parse_scalar(find_element(xinertial, "mass"), "value", T, default = "0")
    end

    return x, q, m, J
end

function parse_robotmaterials(xroot, T)
    xmaterials = get_elements_by_tagname(xroot, "material")

    mdict = Dict{String,Vector{T}}()

    for xmaterial in xmaterials
        name = attribute(xmaterial, "name")
        colornode = find_element(xmaterial, "color")
        colorvec = parse_vector(colornode, "rgba", T)
        mdict[name] = colorvec
    end

    return mdict
end

function parse_xmaterial(xmaterial, materialdict, T)
    if xmaterial === nothing
        color = RGBA(0.75, 0.75, 0.75)
    else
        name = unsafeattribute(xmaterial, "name")
        if haskey(materialdict, name)
            colorvec = materialdict[name]
        else
            colorvec = [0.75; 0.75; 0.75; 1.0]
        end
        colornode = find_element(xmaterial, "color")
        colorvec = parse_vector(colornode, "rgba", T, default = string(colorvec[1]," ",colorvec[2]," ",colorvec[3]," ",colorvec[4]))

        color = RGBA(colorvec...)
    end

    return color
end

function parse_shape(xvisual, materialdict, T)
    if xvisual === nothing
        shape = nothing
    else
        xgeometry = find_element(xvisual, "geometry")
        @assert xgeometry !== nothing

        color = parse_xmaterial(find_element(xvisual, "material"), materialdict, T)
        x, q = parse_pose(find_element(xvisual, "origin"), T)

        shapenodes = LightXML.XMLElement[]
        for node in child_nodes(xgeometry)  # node is an instance of XMLNode
            if is_elementnode(node)
                push!(shapenodes, XMLElement(node))
            end
        end

        if length(shapenodes) == 0
            shape = nothing
        else
            if length(shapenodes) > 1
                @info "Multiple geometries."
            end

            shapenode = shapenodes[1]
            if name(shapenode) == "box"
                xyz = parse_vector(shapenode, "size", T, default = "1 1 1")
                shape = Box(xyz..., zero(T), color = color, xoffset = x, qoffset = q)
            elseif name(shapenode) == "cylinder"
                r = parse_scalar(shapenode, "radius", T, default = "0.5")
                l = parse_scalar(shapenode, "length", T, default = "1")
                shape = Cylinder(r, l, zero(T), color = color, xoffset = x, qoffset = q)
            elseif name(shapenode) == "sphere"
                r = parse_scalar(shapenode, "radius", T, default = "0.5")
                shape = Sphere(r, zero(T), color = color, xoffset = x, qoffset = q)
            elseif name(shapenode) == "mesh"
                path = attribute(shapenode, "filename")
                scale = parse_vector(shapenode, "scale", T, default = "1 1 1")
                shape = Mesh(path, zero(T), zeros(T, 3, 3), scale=scale, color = color, xoffset = x, qoffset = q)
            else
                @info "Unkown geometry."
                shape = nothing
            end
        end

    end

    return shape
end

function parse_link(xlink, materialdict, T)
    xvisual = find_element(xlink, "visual")

    x, q, m, J = parse_inertia(find_element(xlink, "inertial"), T)
    shape = parse_shape(find_element(xlink, "visual"), materialdict, T)
    name = attribute(xlink, "name")

    if shape === nothing
        link = Body(m, J, name=name)
    else
        link = shape
        link.m = m
        link.J = J
        link.name = name
    end

    link.state.xc = x
    link.state.qc = q

    return link
end

function parse_links(xlinks, materialdict, T)
    ldict = Dict{String,AbstractBody{T}}()

    for xlink in xlinks
        link = parse_link(xlink, materialdict, T)
        ldict[link.name] = link
    end

    return ldict
end

function parse_joint(xjoint, plink, clink, T)
    jointtype = attribute(xjoint, "type")
    x, q = parse_pose(find_element(xjoint, "origin"), T)
    axis = parse_vector(find_element(xjoint, "axis"), "xyz", T, default = "1 0 0")
    p1 = x
    name = attribute(xjoint, "name")

    # TODO limits for revolute joint?
    if jointtype == "revolute" || jointtype == "continuous"
        joint = EqualityConstraint(Revolute(plink, clink, axis; p1=p1, qoffset = q), name=name)
    elseif jointtype == "prismatic"
        joint = EqualityConstraint(Prismatic(plink, clink, axis; p1=p1, qoffset = q), name=name)
    elseif jointtype == "planar"
        joint = EqualityConstraint(Planar(plink, clink, axis; p1=p1, qoffset = q), name=name)
    elseif jointtype == "fixed"
        joint = EqualityConstraint(Fixed(plink, clink; p1=p1, qoffset = q), name=name)
    elseif jointtype == "floating"
        joint = EqualityConstraint(Floating(plink, clink), name=name)
    elseif jointtype == "ball"
        joint = EqualityConstraint(Spherical(plink, clink; p1=p1), name=name)
    else
        @error "Unknown joint type"
    end

    return joint
end

function parse_loop_joint(xjoint, link1, link2, T)
    find_element(xjoint, "link1")
    find_element(xjoint, "link2")

    jointtype = attribute(xjoint, "type")
    axis = parse_vector(find_element(xjoint, "axis"), "xyz", T, default = "1 0 0")
    x1, q1 = parse_pose(find_element(xjoint, "link1"), T)
    x2, q2 = parse_pose(find_element(xjoint, "link2"), T) # The orientation q2 of the second body is ignored because it is determined by the mechanism's structure
    p1 = x1
    p2 = x2
    name = attribute(xjoint, "name")

    # TODO limits for revolute joint?
    if jointtype == "revolute" || jointtype == "continuous"
        joint = EqualityConstraint(Revolute(link1, link2, axis; p1=p1, p2=p2, qoffset = q1), name=name)
    elseif jointtype == "prismatic"
        joint = EqualityConstraint(Prismatic(link1, link2, axis; p1=p1, p2=p2, qoffset = q1), name=name)
    elseif jointtype == "planar"
        joint = EqualityConstraint(Planar(link1, link2, axis; p1=p1, p2=p2, qoffset = q1), name=name)
    elseif jointtype == "fixed"
        joint = EqualityConstraint(Fixed(link1, link2; p1=p1, p2=p2, qoffset = q1), name=name)
    elseif jointtype == "floating"
        joint = EqualityConstraint(Floating(link1, link2), name=name)
    elseif jointtype == "ball"
        joint = EqualityConstraint(Spherical(plink, clink; p1=p1, p2=p2), name=name)
    else
        @error "Unknown joint type"
    end

    return joint
end

function parse_joints(xjoints, ldict, floating, T)
    origins = Origin{T}[]
    links = Body{T}[]
    joints = EqualityConstraint{T}[]
    floatingname = ""

    for name in keys(ldict)
        childflag = false
        for xjoint in xjoints
            xchild = find_element(xjoint, "child")
            childname = attribute(xchild, "link")
            if childname == name
                childflag = true
                break
            end
        end
        if childflag
            push!(links, ldict[name])
        else
            origin = Origin{T}()
            if floating # keep current link and create new origin
                push!(links, ldict[name])
                floatingname = name
            else # make current link origin
                origin = Origin(ldict[name])
            end
            push!(origins, origin)
        end
    end

    @assert length(origins) == 1 "Multiple origins"
    origin = origins[1]

    for xjoint in xjoints
        xplink = find_element(xjoint, "parent")
        plink = ldict[attribute(xplink, "link")]
        if plink.id == origin.id
            plink = origin
        end
        xclink = find_element(xjoint, "child")
        clink = ldict[attribute(xclink, "link")]

        joint = parse_joint(xjoint, plink, clink, T)
        push!(joints, joint)
    end

    if floating
        originjoint = EqualityConstraint(Floating(origin, ldict[floatingname]), name="auto_generated_floating_joint")
        push!(joints, originjoint)
    end

    return origin, links, joints
end

# TODO This might be missing the detection of a direct loop, i.e. only two links connected by two joints
# TODO Also only works for a single loop closure in a cycle (so no ladders)
function parse_loop_joints(xloopjoints, origin, joints, ldict, T)
    loopjoints = EqualityConstraint{T}[]

    for xloopjoint in xloopjoints
        xlink1 = find_element(xloopjoint, "link1")
        xlink2 = find_element(xloopjoint, "link2")
        link1 = ldict[attribute(xlink1, "link")]
        link2 = ldict[attribute(xlink2, "link")]

        predlist = Tuple{Int64,Int64}[]
        jointlist = [(joints[i].id,joints[i].parentid,joints[i].childids) for i=1:length(joints)]
        linkid = link1.id

        while true # create list of predecessor joints and parent links for link1
            for (i,jointdata) in enumerate(jointlist)
                if linkid ∈ jointdata[3]
                    push!(predlist,(jointdata[1],jointdata[2]))
                    linkid = jointdata[2]
                    deleteat!(jointlist,i)
                    break
                end
            end
            if linkid == origin.id
                break
            end
        end

        jointlist = [(joints[i].id,joints[i].parentid,joints[i].childids) for i=1:length(joints)]
        linkid = link2.id
        joint1id = 0
        joint2id = 0
        foundflag = false

        while true # check which predecessor link of link2 is also a predecessor link of link1
            for (i,jointdata) in enumerate(jointlist)
                if linkid ∈ jointdata[3]
                    joint2id = jointdata[1]
                    linkid = jointdata[2]
                    deleteat!(jointlist,i)
                    break
                end
            end
            for el in predlist
                if linkid == el[2]
                    joint1id = el[1]
                    foundflag = true
                    break
                end
            end
            foundflag && break
        end

        # Find and remove joints to combine them
        joint1 = 0
        joint2 = 0
        for (i,joint) in enumerate(joints)
            if joint.id == joint1id
                joint1 = joint
                deleteat!(joints,i)
                break
            end
        end
        # if joint1id == joint2id # already combined joint
        #     joint2 = joint
        for (i,joint) in enumerate(joints)
            if joint.id == joint2id
                joint2 = joint
                deleteat!(joints,i)
                break
            end
        end

        joint = cat(joint1,joint2)
        push!(joints,joint)
        loopjoint = parse_loop_joint(xloopjoint, link1, link2, T)
        push!(loopjoints, loopjoint)
    end

    return joints, loopjoints
end

function parse_urdf(filename, floating, ::Type{T}) where T
    xdoc = LightXML.parse_file(filename)
    xroot = LightXML.root(xdoc)
    @assert LightXML.name(xroot) == "robot"

    xlinks = get_elements_by_tagname(xroot, "link")
    xjoints = get_elements_by_tagname(xroot, "joint")
    xloopjoints = get_elements_by_tagname(xroot, "loop_joint")

    materialdict = parse_robotmaterials(xroot, T)
    ldict = parse_links(xlinks, materialdict, T)
    origin, links, joints = parse_joints(xjoints, ldict, floating, T)

    joints, loopjoints = parse_loop_joints(xloopjoints, origin, joints, ldict, T)

    free(xdoc)

    return origin, links, joints, loopjoints
end

### After parsing

function set_parsed_values!(mechanism::Mechanism{T}, loopjoints) where T
    graph = mechanism.graph
    xjointlist = Dict{Int64,SVector{3,T}}() # stores id, x in world frame
    qjointlist = Dict{Int64,UnitQuaternion{T}}() # stores id, q in world frame

    for id in graph.rdfslist # from root to leaves
        component = getcomponent(mechanism, id)
        !(component isa Body) && continue # only for bodies

        body = component
        xbodylocal = body.state.xc
        qbodylocal = body.state.qc
        shape = body.shape

        parentid = get_parentid(mechanism, id, loopjoints)
        constraint = geteqconstraint(mechanism, parentid)

        grandparentid = constraint.parentid
        if grandparentid === nothing # predecessor is origin
            parentbody = mechanism.origin

            xparentbody = SA{T}[0; 0; 0]
            qparentbody = one(UnitQuaternion{T})

            xparentjoint = SA{T}[0; 0; 0]
            qparentjoint = one(UnitQuaternion{T})
        else
            parentbody = getbody(mechanism, grandparentid)

            grandgrandparentid = get_parentid(mechanism, grandparentid, loopjoints)
            parentconstraint = geteqconstraint(mechanism, grandgrandparentid)

            xparentbody = parentbody.state.xc # in world frame
            qparentbody = parentbody.state.qc # in world frame

            xparentjoint = xjointlist[parentconstraint.id] # in world frame
            qparentjoint = qjointlist[parentconstraint.id] # in world frame
        end

        ind1 = findfirst(x->x==id,constraint.childids)
        ind2 = ind1+1

        # urdf joint's x and q in parent's (parentbody) frame
        xjointlocal = vrotate(xparentjoint + vrotate(constraint.constraints[ind1].vertices[1], qparentjoint) - xparentbody, inv(qparentbody))
        qjointlocal = qparentbody \ qparentjoint * constraint.constraints[ind2].qoffset

        # store joint's x and q in world frame
        xjoint = xparentbody + vrotate(xjointlocal, qparentbody)
        qjoint = qparentbody * qjointlocal
        xjointlist[constraint.id] = xjoint
        qjointlist[constraint.id] = qjoint

        # difference to parent body (parentbody)
        qoffset = qjointlocal * qbodylocal

        # actual joint properties
        p1 = xjointlocal # in parent's (parentbody) frame
        p2 = vrotate(-xbodylocal, inv(qbodylocal)) # in body frame (xbodylocal and qbodylocal are both relative to the same (joint) frame -> rotationg by inv(body.q) gives body frame)
        constraint.constraints[ind1].vertices = (p1, p2)

        V3 = vrotate(constraint.constraints[ind2].V3', qjointlocal) # in parent's (parentbody) frame
        V12 = (svd(skew(V3)).Vt)[1:2,:]
        constraint.constraints[ind2].V3 = V3'
        constraint.constraints[ind2].V12 = V12
        constraint.constraints[ind2].qoffset = qoffset # in parent's (parentbody) frame

        # actual body properties
        setPosition!(body) # set everything to zero
        setPosition!(parentbody, body, p1 = p1, p2 = p2, Δq = qoffset)
        xbody = body.state.xc
        qbody = body.state.qc

        # shape relative
        if !(typeof(shape) <: EmptyShape)
            shape.xoffset = vrotate(xjoint + vrotate(shape.xoffset, qjoint) - xbody, inv(qbody))
            shape.qoffset = qoffset \ qjointlocal * shape.qoffset
        end
    end
    for (i,constraint) in enumerate(loopjoints)
        @assert length(constraint.childids) == 2 # Loop joint connects only two bodies

        parentid1 = constraint.parentid
        parentid2 = constraint.childids[1]
        if parentid1 === nothing # predecessor is origin
            parentbody1 = mechanism.origin

            xparentbody1 = SA{T}[0; 0; 0]
            qparentbody1 = one(UnitQuaternion{T})

            xparentjoint1 = SA{T}[0; 0; 0]
            qparentjoint1 = one(UnitQuaternion{T})
        else
            parentbody1 = getbody(mechanism, parentid1)

            grandparentid1 = get_parentid(mechanism, parentid1, loopjoints)
            parentconstraint1 = geteqconstraint(mechanism, grandparentid1)

            xparentbody1 = parentbody1.state.xc # in world frame
            qparentbody1 = parentbody1.state.qc # in world frame

            xparentjoint1 = xjointlist[parentconstraint1.id] # in world frame
            qparentjoint1 = qjointlist[parentconstraint1.id] # in world frame
        end
        parentbody2 = getbody(mechanism, parentid2)

        grandparentid2 = get_parentid(mechanism, parentid2, loopjoints)
        parentconstraint2 = geteqconstraint(mechanism, grandparentid2)

        xparentbody2 = parentbody2.state.xc # in world frame
        qparentbody2 = parentbody2.state.qc # in world frame

        xparentjoint2 = xjointlist[parentconstraint2.id] # in world frame
        qparentjoint2 = qjointlist[parentconstraint2.id] # in world frame

        ind1 = 1
        ind2 = ind1+1

        # urdf joint's x and q in parent's (parentbody) frame
        xjointlocal1 = vrotate(xparentjoint1 + vrotate(constraint.constraints[ind1].vertices[1], qparentjoint1) - xparentbody1, inv(qparentbody1))
        xjointlocal2 = vrotate(xparentjoint2 + vrotate(constraint.constraints[ind1].vertices[2], qparentjoint2) - xparentbody2, inv(qparentbody2))
        qjointlocal1 = qparentbody1 \ qparentjoint1 * constraint.constraints[ind2].qoffset

        # difference to parent body (parentbody)
        qoffset1 = qjointlocal1 * qparentbody2 #  qparentbody2 = body in for loop above

        # actual joint properties
        p1 = xjointlocal1 # in parent's (parentbody1) frame
        p2 = xjointlocal2 # in parent's (parentbody2) frame
        constraint.constraints[ind1].vertices = (p1, p2)

        V3 = vrotate(constraint.constraints[ind2].V3', qjointlocal1) # in parent's (parentbody1) frame
        V12 = (svd(skew(V3)).Vt)[1:2,:]
        constraint.constraints[ind2].V3 = V3'
        constraint.constraints[ind2].V12 = V12
        constraint.constraints[ind2].qoffset = qoffset1 # in parent's (parentbody1) frame
    end
end

function get_parentid(mechanism, id, loopjoints)
    graph = mechanism.graph
    conns = connections(graph, id)
    for connsid in conns
        constraint = geteqconstraint(mechanism, connsid)
        if constraint ∉ loopjoints && id ∈ constraint.childids
            return connsid
        end
    end

    return nothing
end


text:
using Regression
using Distributions
using PyPlot

function tranform(x::Array{Float64}, dim::Int64=10)
    len = length(x)
    d = Array(Float64, len, dim)
    for i = 1:len
        for j = 1:dim
            d[i,j] = (x[i])^(j-1)
        end
    end
    d
end

function f_e(x,thetas)
    r = 0
    for i = 1:length(thetas)
        r += thetas[i]*x^(i-1)
    end
    r
end

# prepare simulation data
f(x) = 0.22*x^2-1.3*x^3+0.1*x^5
noise = Normal(0,0.9)
xs_o = rand(100)*4
ys = map(f, xs_o) + rand(noise, 100)

dim = 6                      # feature vector with n dimension
model = LinearPred(dim)      # the equatin (2)
sqrloss = SqrLoss()          # the loss function (3)
reg = SqrL2Reg(2.0)          # squared L2 regularizer with coefficient 1
xs = tranform(xs_o,dim)

# this is (4)
rmodel = riskmodel(model,sqrloss)
# assume xs is the 10 by m matrix(list of feature vectors)
# and the ys is a vector of labels either -1 or 1
objective_function = Regression.RegRiskFun(rmodel, reg, xs', ys)

thetas = Array(Float64, dim)
# some adjustable parameters for optimization
options = Regression.Options()
# use gradient descent
solver = Regression.GD()
# run the optimization
Regression.solve!(solver, objective_function, thetas, options, Regression.no_op)
# clean up weak features
thetas = map(x-> abs(x) < 0.0001? 0 : x, thetas)

# visualize the data and the line
clf()
scatter(xs_o,ys,c="g")
min_x,max_x = foldr(min, Inf, xs_o),foldr(max, -Inf, xs_o)
plot_x = [min_x:0.01:max_x]
plot_y = map(x->f_e(x,thetas), plot_x)
plot(plot_x, plot_y, c="b")


text:
# http://homepages.inf.ed.ac.uk/imurray2/pub/10ess/elliptical_slice.m
# https://github.com/jobovy/bovy_mcmc/blob/master/bovy_mcmc/elliptical_slice.py

function elliptical_slice(xx, prior::Array{Float64,}, log_like_fn, cur_log_like)
    D = length(xx);
    cur_log_like = log_like_fn(xx);
    angle_range = 0;
    hh = log(rand()) + cur_log_like;

    # Set up a bracket of angles and pick a first proposal.
    # "phi = (theta'-theta)" is a change in angle.
    if angle_range <= 0
        # Bracket whole ellipse with both edges at first proposed point
        phi = rand() * 2*pi;
        phi_min = phi - 2*pi;
        phi_max = phi;
    else
        # Randomly center bracket on current point
        phi_min = -angle_range*rand();
        phi_max = phi_min + angle_range;
        phi = rand()*(phi_max - phi_min) + phi_min;
    end

    xx_prop = xx  # TODO: Need this due to scoping rules?

    # Slice sampling loop
    while true
        # Compute xx for proposed angle difference and check if it's on the slice
        xx_prop = xx*cos(phi) + nu*sin(phi);
        cur_log_like = log_like_fn(xx_prop);
        if cur_log_like > hh
            # New point is on slice, ** EXIT LOOP **
            break;
        end
        # Shrink slice to rejected point
        if phi > 0
            phi_max = phi;
        elseif phi < 0
            phi_min = phi;
        else
            error("BUG DETECTED: Shrunk to current position and still not acceptable.");
        end
        # Propose new angle difference
        phi = rand() * (phi_max - phi_min) + phi_min;
    end
    return xx_prop, cur_log_like
end

# User specified Cholesky of prior covariance:
function elliptical_slice(xx, Sigma::Array{Float64,2}, llk, cur_llk)
    nu = reshape(Sigma'*randn(D, 1), size(xx));
    elliptical_slice(xx, nu, llk, cur_llk)
end


text:
module UnitfulHartree
using Unitful
import Unitful

macro lintpragma(s) end
@lintpragma("Ignore use of undeclared variable ElectronMass")
@lintpragma("Ignore use of undeclared variable mₑ")
@lintpragma("Ignore use of undeclared variable ElementaryCharge")
@lintpragma("Ignore use of undeclared variable e₀")
@lintpragma("Ignore use of undeclared variable CoulombForceConstant")
@lintpragma("Ignore use of undeclared variable kₑ")
@lintpragma("Ignore use of undeclared variable ReducedPlanckConstant")
@lintpragma("Ignore use of undeclared variable ħ")
@lintpragma("Ignore use of undeclared variable BohrRadius")
@lintpragma("Ignore use of undeclared variable a₀")
@lintpragma("Ignore use of undeclared variable HartreeEnergy")
@lintpragma("Ignore use of undeclared variable Eₕ")
@lintpragma("Ignore use of undeclared variable RydbergEnergy")
@lintpragma("Ignore use of undeclared variable Ry")
@lintpragma("Ignore use of undeclared variable ClassicalElectronRadius")
@lintpragma("Ignore use of undeclared variable rₑ")
@lintpragma("Ignore use of undeclared variable Density")
@lintpragma("Ignore use of undeclared variable ρ")
@lintpragma("Ignore use of undeclared variable ContractedDensityGradient")
@lintpragma("Ignore use of undeclared variable σₑ")
@lintpragma("Ignore use of undeclared variable FirstDensityDerivative")
@lintpragma("Ignore use of undeclared variable ∂ϵ_∂ρ")
@lintpragma("Ignore use of undeclared variable FirstGradientDerivative")
@lintpragma("Ignore use of undeclared variable ∂ϵ_∂σ")
@lintpragma("Ignore use of undeclared variable SecondDensityDerivative")
@lintpragma("Ignore use of undeclared variable ∂²ϵ_∂ρ²")
@lintpragma("Ignore use of undeclared variable SecondGradientDerivative")
@lintpragma("Ignore use of undeclared variable ∂²ϵ_∂σ²")
@lintpragma("Ignore use of undeclared variable SecondDensityGradientDerivative")
@lintpragma("Ignore use of undeclared variable ∂²ϵ_∂ρ∂σ")
@lintpragma("Ignore use of undeclared variable ThirdDensityDerivative")
@lintpragma("Ignore use of undeclared variable ∂³ϵ_∂ρ³")
@lintpragma("Ignore use of undeclared variable ThirdGradientDerivative")
@lintpragma("Ignore use of undeclared variable ∂³ϵ_∂σ³")
@lintpragma("Ignore use of undeclared variable ThirdDensity2GradientDerivative")
@lintpragma("Ignore use of undeclared variable ∂³ϵ_∂ρ²∂σ")
@lintpragma("Ignore use of undeclared variable ThirdDensityGradient2Derivative")
@lintpragma("Ignore use of undeclared variable ∂³ϵ_∂ρ∂σ²")

@unit mₑ "mₑ" ElectronMass 9.1093835611e-31*u"kg" false
@unit e₀ "eₒ" ElementaryCharge 1.602176620898e-19*u"C" false
@unit kₑ "kₑ" CoulombForceConstant 1/(4π)u"ϵ0" false
@unit ħ "ħ" ReducedPlanckConstant Unitful.ħ false
@unit a₀ "a₀" BohrRadius 1ħ^2/(1kₑ*mₑ*e₀^2) false
@unit Eₕ "Eₕ" HartreeEnergy mₑ*e₀^4*kₑ^2/(1ħ^2) true
@unit Ry "Ry" RydbergEnergy 0.5Eₕ true
@unit rₑ "rₑ" ClassicalElectronRadius 1e₀^2*kₑ/(1mₑ*Unitful.c^2) false
const α = 1e₀^2*1kₑ/(1Unitful.c*ħ)
const mₚ = 1836.15mₑ
const μ_b = e₀*ħ/(2mₑ)
const ϵ₀ = 1/(4π*kₑ)

@unit ρ         "ρ"         Density                         (1a₀)^-3       false
@unit σₑ        "σₑ"        ContractedDensityGradient       (1a₀)^-8       false
@unit ∂ϵ_∂ρ     "∂ϵ_∂ρ"     FirstDensityDerivative          (1Eₕ)*(1a₀)^3  false
@unit ∂ϵ_∂σ     "∂ϵ_∂σ"     FirstGradientDerivative         (1Eₕ)*(1a₀)^8  false
@unit ∂²ϵ_∂ρ²   "∂²ϵ_∂ρ²"   SecondDensityDerivative         (1Eₕ)*(1a₀)^6  false
@unit ∂²ϵ_∂σ²   "∂²ϵ_∂σ²"   SecondGradientDerivative        (1Eₕ)*(1a₀)^16 false
@unit ∂²ϵ_∂ρ∂σ  "∂²ϵ_∂ρ∂σ"  SecondDensityGradientDerivative (1Eₕ)*(1a₀)^11 false
@unit ∂³ϵ_∂ρ³   "∂³ϵ_∂ρ³"   ThirdDensityDerivative          (1Eₕ)*(1a₀)^9  false
@unit ∂³ϵ_∂σ³   "∂³ϵ_∂σ³"   ThirdGradientDerivative         (1Eₕ)*(1a₀)^24 false
@unit ∂³ϵ_∂ρ²∂σ "∂³ϵ_∂ρ²∂σ" ThirdDensity2GradientDerivative (1Eₕ)*(1a₀)^14 false
@unit ∂³ϵ_∂ρ∂σ² "∂³ϵ_∂ρ∂σ²" ThirdDensityGradient2Derivative (1Eₕ)*(1a₀)^19 false
const ϵ = Eₕ
const rho = ρ
const sig = σₑ
const Eh = Eₕ

# Some gymnastics required here because if we precompile, we cannot add to
# Unitful.basefactors at compile time and expect the changes to persist to runtime.
const localunits = Unitful.basefactors
function __init__()
    merge!(Unitful.basefactors, localunits)
    Unitful.register(UnitfulHartree)
end

end # module

max_stars_repo_name: mdavezac/DFTShims

max_issues_repo_name: mdavezac/DFTShims

max_forks_repo_name: mdavezac/DFTShims


text:
mutable struct LidToOHId <: AbstractVector{Int32}
  nlids::Int
  noids::Int32
end

function Base.copy(a::LidToOHId)
  LidToOHId(
    copy(a.nlids),
    copy(a.noids))
end

Base.size(a::LidToOHId) = (a.nlids,)
Base.IndexStyle(::Type{<:LidToOHId}) = IndexLinear()

@inline function Base.getindex(a::LidToOHId,lid::Integer)
  @boundscheck begin
    if !( 1<=lid && lid<=length(a) )
      throw(BoundsError(a,lid))
    end
  end
  if lid<=a.noids
    Int32(lid)
  else
    a.noids-Int32(lid)
  end
end

@inline function Base.push!(a::LidToOHId,ohid::Integer)
  @boundscheck begin
    if (a.noids - ohid) != (length(a) + 1)
      throw(DomainError(ohid))
    end
  end
  a.nlids += 1
  a
end

struct LidToGid <: AbstractVector{Int}
  oid_to_gid::UnitRange{Int}
  hid_to_gid::Vector{Int}
end

function Base.copy(a::LidToGid)
  LidToGid(
    copy(a.oid_to_gid),
    copy(a.hid_to_gid))
end

Base.size(a::LidToGid) = (length(a.oid_to_gid)+length(a.hid_to_gid),)
Base.IndexStyle(::Type{<:LidToGid}) = IndexLinear()

@inline function Base.getindex(a::LidToGid,lid::Integer)
  @boundscheck begin
    if !( 1<=lid && lid<=length(a) )
      throw(BoundsError(a,lid))
    end
  end
  noids = length(a.oid_to_gid)
  if lid <= noids
    oid = lid
    a.oid_to_gid[oid]
  else
    hid = lid - noids
    a.hid_to_gid[hid]
  end
end

@inline function Base.push!(a::LidToGid,gid::Integer)
  push!(a.hid_to_gid,gid)
  a
end

struct LidToPart <: AbstractVector{Int32}
  noids::Int
  part::Int32
  hid_to_part::Vector{Int32}
end

function Base.copy(a::LidToPart)
  LidToPart(
    copy(a.noids),
    copy(a.part),
    copy(a.hid_to_part))
end

Base.size(a::LidToPart) = (a.noids+length(a.hid_to_part),)
Base.IndexStyle(::Type{<:LidToPart}) = IndexLinear()

@inline function Base.getindex(a::LidToPart,lid::Integer)
  @boundscheck begin
    if !( 1<=lid && lid<=length(a) )
      throw(BoundsError(a,gid))
    end
  end
  if lid <= a.noids
    a.part
  else
    hid = lid - a.noids
    a.hid_to_part[hid]
  end
end

@inline function Base.push!(a::LidToPart,part::Integer)
  push!(a.hid_to_part,part)
  a
end

struct GidToLid{T} <: AbstractDict{Int,Int32}
  oid_to_gid::UnitRange{Int}
  oid_to_lid::T
  gid_to_lid::Dict{Int,Int32} # Only for ghost
end

function Base.copy(a::GidToLid)
  GidToLid(
    copy(a.oid_to_gid),
    copy(a.oid_to_lid),
    copy(a.gid_to_lid))
end

Base.length(a::GidToLid) = length(a.oid_to_gid) + length(a.gid_to_lid)
Base.haskey(a::GidToLid,gid::Integer) = (gid in a.oid_to_gid) || haskey(a.gid_to_lid,gid)

function Base.keys(a::GidToLid)
  Base.Iterators.flatten((a.oid_to_gid,keys(a.gid_to_lid)))
end

function Base.values(a::GidToLid)
  Base.Iterators.flatten(
    (Int32(1):Int32(length(a.oid_to_gid)),values(a.gid_to_lid)))
end

@inline function Base.getindex(a::GidToLid,gid::Integer)
  @boundscheck begin
    if ! haskey(a,gid)
      throw(KeyError(gid))
    end
  end
  if gid in a.oid_to_gid
    oid = 1 + gid - first(a.oid_to_gid)
    a.oid_to_lid[oid]
  else
    a.gid_to_lid[gid]
  end
end

@inline function Base.setindex!(a::GidToLid,lid::Integer,gid::Integer)
  if ! haskey(a,gid)
    a.gid_to_lid[gid] = lid
  end
  gid
end

function Base.iterate(a::GidToLid)
  itr = zip(keys(a),values(a))
  next = iterate(itr)
  if next == nothing
    return nothing
  end
  item, state = next
  item[1]=>item[2], (itr,state)
end

function Base.iterate(a::GidToLid,(itr,state))
  next = iterate(itr,state)
  if next == nothing
    return nothing
  end
  item, state = next
  item[1]=>item[2], (itr,state)
end

struct LinearGidToPart <:AbstractVector{Int32}
  ngids::Int
  part_to_firstgid::Vector{Int}
end

function Base.copy(a::LinearGidToPart)
  LinearGidToPart(a.ngids,copy(a.part_to_firstgid))
end

Base.size(a::LinearGidToPart) = (a.ngids,)
Base.IndexStyle(::Type{<:LinearGidToPart}) = IndexLinear()

function Base.getindex(a::LinearGidToPart,gid::Integer)
  @boundscheck begin
    if !( 1<=gid && gid<=length(a) )
      throw(BoundsError(a,gid))
    end
  end
  Int32(searchsortedlast(a.part_to_firstgid,gid))
end

struct CartesianGidToPart{N} <: AbstractVector{Int32}
  ngids::NTuple{N,Int}
  part_to_firstgid::NTuple{N,Vector{Int}}
end

function Base.copy(a::CartesianGidToPart)
  CartesianGidToPart(a.ngids,map(copy,a.part_to_firstgid))
end

Base.size(a::CartesianGidToPart) = (prod(a.ngids),)
Base.IndexStyle(::Type{<:CartesianGidToPart}) = IndexLinear()

function Base.getindex(a::CartesianGidToPart,gid::Integer)
  cgid = Tuple(CartesianIndices(a.ngids)[gid])
  cpart = map(searchsortedlast,a.part_to_firstgid,cgid)
  nparts = map(length,a.part_to_firstgid)
  part = LinearIndices(nparts)[CartesianIndex(cpart)]
  Int32(part)
end

struct IndexSet <: AbstractIndexSet
  part::Int
  lid_to_gid::Vector{Int}
  lid_to_part::Vector{Int32}
  oid_to_lid::Vector{Int32}
  hid_to_lid::Vector{Int32}
  lid_to_ohid::Vector{Int32}
  gid_to_lid::Dict{Int,Int32}
end

function IndexSet(
  part::Int,
  lid_to_gid::Vector{Int},
  lid_to_part::Vector{Int32},
  oid_to_lid::Vector{Int32},
  hid_to_lid::Vector{Int32},
  lid_to_ohid::Vector{Int32})

  gid_to_lid = Dict{Int,Int32}()
  for (lid,gid) in enumerate(lid_to_gid)
    gid_to_lid[gid] = lid
  end
  IndexSet(
    part,
    lid_to_gid,
    lid_to_part,
    oid_to_lid,
    hid_to_lid,
    lid_to_ohid,
    gid_to_lid)
end

function IndexSet(
  part::Int,
  lid_to_gid::Vector{Int},
  lid_to_part::Vector{Int32},
  oid_to_lid::Vector{Int32},
  hid_to_lid::Vector{Int32})

  lid_to_ohid = zeros(Int32,length(lid_to_gid))
  lid_to_ohid[oid_to_lid] = 1:length(oid_to_lid)
  lid_to_ohid[hid_to_lid] = -(1:length(hid_to_lid))

  IndexSet(
    part,
    lid_to_gid,
    lid_to_part,
    oid_to_lid,
    hid_to_lid,
    lid_to_ohid)
end

function IndexSet(
  part::Int,
  lid_to_gid::Vector{Int},
  lid_to_part::Vector{Int32})

  oid_to_lid = collect(Int32,findall(owner->owner==part,lid_to_part))
  hid_to_lid = collect(Int32,findall(owner->owner!=part,lid_to_part))
  IndexSet(
    part,
    lid_to_gid,
    lid_to_part,
    oid_to_lid,
    hid_to_lid)
end

function Base.copy(a::IndexSet)
  IndexSet(
    copy(a.part),
    copy(a.lid_to_gid),
    copy(a.lid_to_part),
    copy(a.oid_to_lid),
    copy(a.hid_to_lid),
    copy(a.lid_to_ohid),
    copy(a.gid_to_lid))
end

struct ExtendedIndexRange <: AbstractIndexSet
  part::Int
  lid_to_gid::Vector{Int}
  lid_to_part::Vector{Int32}
  oid_to_lid::Vector{Int32}
  hid_to_lid::Vector{Int32}
  lid_to_ohid::Vector{Int32}
  gid_to_lid::GidToLid{Vector{Int32}}
end

function ExtendedIndexRange(
  part::Integer,
  lid_to_gid::Vector{Int},
  lid_to_part::Vector{Int32},
  firstgid::Integer)

  oid_to_lid = collect(Int32,findall(owner->owner==part,lid_to_part))
  hid_to_lid = collect(Int32,findall(owner->owner!=part,lid_to_part))
  noids = length(oid_to_lid)
  lid_to_ohid = zeros(Int32,length(lid_to_gid))
  lid_to_ohid[oid_to_lid] = 1:length(oid_to_lid)
  lid_to_ohid[hid_to_lid] = -(1:length(hid_to_lid))
  offset = firstgid-1
  oid_to_gid = Int(1+offset):Int(noids+offset)
  gid_to_lid_for_ghost = Dict{Int,Int32}()
  for lid in hid_to_lid
    gid_to_lid_for_ghost[lid_to_gid[lid]] = lid
  end
  gid_to_lid = GidToLid(oid_to_gid,oid_to_lid,gid_to_lid_for_ghost)
  ExtendedIndexRange(
    part,
    lid_to_gid,
    lid_to_part,
    oid_to_lid,
    hid_to_lid,
    lid_to_ohid,
    gid_to_lid)
end

function Base.copy(a::ExtendedIndexRange)
  ExtendedIndexRange(
    copy(a.part),
    copy(a.lid_to_gid),
    copy(a.lid_to_part),
    copy(a.oid_to_lid),
    copy(a.hid_to_lid),
    copy(a.lid_to_ohid),
    copy(a.gid_to_lid))
end

struct IndexRange <: AbstractIndexSet
  part::Int
  lid_to_gid::LidToGid
  lid_to_part::LidToPart
  oid_to_lid::UnitRange{Int32}
  hid_to_lid::Vector{Int32}
  lid_to_ohid::LidToOHId
  gid_to_lid::GidToLid{UnitRange{Int32}}
end

function Base.copy(a::IndexRange)
  IndexRange(
    copy(a.part),
    copy(a.lid_to_gid),
    copy(a.lid_to_part),
    copy(a.oid_to_lid),
    copy(a.hid_to_lid),
    copy(a.lid_to_ohid),
    copy(a.gid_to_lid))
end

function IndexRange(
  part::Integer,
  noids::Integer,
  firstgid::Integer)

  offset = firstgid-1
  oid_to_gid = Int(1+offset):Int(noids+offset)
  hid_to_gid = Int[]
  lid_to_gid = LidToGid(oid_to_gid,hid_to_gid)
  hid_to_part = Int32[]
  lid_to_part = LidToPart(noids,part,hid_to_part)
  oid_to_lid = Int32(1):Int32(noids)
  hid_to_lid = Int32[]
  nlids = Int(noids)
  lid_to_ohid = LidToOHId(nlids,Int32(noids))
  gid_to_lid_for_ghost = Dict{Int,Int32}()
  gid_to_lid = GidToLid(oid_to_gid,oid_to_lid,gid_to_lid_for_ghost)
  IndexRange(
    part,
    lid_to_gid,
    lid_to_part,
    oid_to_lid,
    hid_to_lid,
    lid_to_ohid,
    gid_to_lid)
end

function IndexRange(
  part::Integer,
  noids::Integer,
  firstgid::Integer,
  hid_to_gid::Vector{Int},
  hid_to_part::Vector{Int32})

  offset = firstgid-1
  oid_to_gid = Int(1+offset):Int(noids+offset)
  lid_to_gid = LidToGid(oid_to_gid,hid_to_gid)
  lid_to_part = LidToPart(noids,part,hid_to_part)
  oid_to_lid = Int32(1):Int32(noids)
  nhids = length(hid_to_part)
  hid_to_lid = collect(Int32, (1:nhids) .+ noids )
  gid_to_lid_for_ghost = Dict{Int,Int32}()
  for (hid,gid) in enumerate(hid_to_gid)
    lid = hid + noids
    gid_to_lid_for_ghost[gid] = lid
  end
  nlids = Int(noids+nhids)
  lid_to_ohid = LidToOHId(nlids,Int32(noids))
  gid_to_lid = GidToLid(oid_to_gid,oid_to_lid,gid_to_lid_for_ghost)
  IndexRange(
    part,
    lid_to_gid,
    lid_to_part,
    oid_to_lid,
    hid_to_lid,
    lid_to_ohid,
    gid_to_lid)
end


text:
#################### ChainSummary ####################

#################### Types and Constructors ####################

struct ChainSummary
  value::Array{Float64, 3}
  rownames::Vector{AbstractString}
  colnames::Vector{AbstractString}
  header::AbstractString

  function ChainSummary(value::Array{Float64, 3},
                        rownames::Vector{AbstractString},
                        colnames::Vector{AbstractString},
                        header::AbstractString)
    dim = size(value)
    length(rownames) == dim[1] ||
      throw(DimensionMismatch("numbers of rownames and rows differ"))
    length(colnames) == dim[2] ||
      throw(DimensionMismatch("numbers of colnames and columns differ"))
    new(value, rownames, colnames, header)
  end
end

function ChainSummary(value::Array{Float64, 3}, rownames::Vector{T},
                     colnames::Vector{U}, header::AbstractString) where {T<:AbstractString, U<:AbstractString}
  ChainSummary(copy(value), AbstractString[rownames...],
               AbstractString[colnames...], header)
end

function ChainSummary(value::Matrix{Float64}, rownames::Vector{T},
                     colnames::Vector{U}, header::AbstractString) where {T<:AbstractString, U<:AbstractString}
  dim = size(value)
  ChainSummary(reshape(value, dim[1], dim[2], 1), AbstractString[rownames...],
               AbstractString[colnames...], header)
end

#################### Base Methods ####################

## write n ' ' characters to io
wrtsp(io::IO, n) = while (n -= 1) >= 0 write(io, ' ') end

function Base.show(io::IO, s::ChainSummary)
  ## rowname widths
  rnwid = map(length, s.rownames)
  mxrnwid = maximum(rnwid)
  ## column name widths
  cnwid = map(length, s.colnames)
  ## s.value as right alignable strings
  charv = mapslices(showoff, s.value, dims = [1])
  colwid = 1 .+ max.(cnwid, vec(maximum(map(length, charv), dims = [1, 3] )))
  m, n, f = size(charv)
  for k in 1:f
    ## write the column headers centered on the column widths
    wrtsp(io, mxrnwid)
    for j in 1:n
      ## do not count the leading space
      nspace = colwid[j] - cnwid[j] - 1
      ## divide by 2 rounding down
      nright = nspace >> 1
      wrtsp(io, 1 + nspace - nright)
      print(io, s.colnames[j])
      wrtsp(io, nright)
    end
    println(io)
    for i in 1:m
      wrtsp(io, mxrnwid - rnwid[i])
      print(io, s.rownames[i])
      for j in 1:n
        wrtsp(io, colwid[j] - length(charv[i, j, k]))
        print(io, charv[i, j, k])
      end
      println(io)
    end
    println(io)
  end
end


text:
include("RotationAngle.jl")
include("ProjectedRotationAngle.jl")


text:
# io.jl: functions to do IO more effiently

@doc """
### Utils.BufferedIO

  This type provides a means to buffer IO in memory before writing it to a file.
  Data written to the object is stored in the IOBuffer until flush() is called,
  when the buffer is (efficiently) dumped into the file
"""->
mutable struct BufferedIO{T <: IO}  <: IO
  fstream::T
  fbuf::IOBuffer
end

@doc """
### Utils.BufferedIO

  Constructor for BufferedIO type.  Takes an IOStream and creates an IOBuffer.
  If no IOStream is given, a dummy stream is created.  If a dummy stream is
  used, flush must never be called on the returned BufferedIO object

  Inputs:

    f: an IOStream object, defaults to a dummy stream

  Outputs:

    a BufferedIO object
"""->
function BufferedIO(f::IO=DevNull)
  buf = IOBuffer()
  fbuf = BufferedIO{typeof(f)}(f, buf)

  # register atexit hook to make sure any buffered data is flushed before
  # julia exits
  # this causes errors to be printed when julia exists, possible #10431
  atexit( () -> if isopen(fbuf.fstream)
                  close(fbuf)
                end
        )
  return fbuf
end

"""
  Alternative constructor for BufferedIO, emulating the open() function.
  This function creates the underlying file using open() and then creates
  a BufferedIO around it.

  Inputs:

    fname: AbstractString, name of file to open
    mode: file open mode, see documentation of open(), defaults to append

  Outputs:

    a BufferedIO object
"""
function BufferedIO(fname::AbstractString, mode::AbstractString="a")

  f = open(fname, mode)
  return BufferedIO(f)
end
# only provide write functionality, for now
import Base.write, Base.flush, Base.close, Base.isopen

# I think this is all that needs to be implemented, because it is the job of
# functions like println and print to convert things to arrays of UInt8s
write(io::BufferedIO, x::UInt8) = write(io.fbuf, x)

"""
  `Base` function `flush` extended for BufferedIO
"""
function flush(io::BufferedIO)
  write(io.fstream, take!(io.fbuf))
  flush(io.fstream)
end

"""
  `Base` function `close` extended for BufferedIO
"""
function close(io::BufferedIO)
  if isopen(io.fstream)
    flush(io)
    close(io.fstream)
  end
end

"""
  `Base` function `isopen` extended for BufferedIO
"""
function isopen(io::BufferedIO)
  return isopen(io.fstream)
end

#------------------------------------------------------------------------------
# create some useful streams

"""
  Buffered version of STDOUT.  This should *always* be used instead of STDOUT
"""
global const BSTDOUT = BufferedIO(STDOUT)

"""
  Buffered version of STDERR.  This should *always* be used instead of STDERR
"""
global const BSTDERR = BufferedIO(STDERR)

#=
#------------------------------------------------------------------------------
# create file io to numbered file with time stamp, unbuffered

"""
  This IO type writes the output to a numbered file (see the constructor),
  without buffering, prefixing the output with a time stamp using MPI_Wtime
  Note that only println() includes the time stamp prefix, print() does not.
"""
type DebugFileIO  <: IO
  fstream::IO  # TODO: make this an IOStream (concrete)
end

@doc """
### Utils.DebugFileIO

  Constructor for DebugFileIO type.  Creates a numbered file by appending
  the supplied rank to the file name (see get_parallel_fname in ODLCommonTools)

  Inputs:
    fname: file name, including extension
    mode: mode flags for opening the file (see docs for open() )
    rank: the rank of the process

  Outputs:
    a DebugFileIO object
"""->
function DebugFileIO(fname::AbstractString, mode::AbstractString, rank::Integer )

  f = open( get_parallel_fname(fname, rank), mode)
  fbuf = DebugFileIO(f)

  # register atexit hook to make sure any buffered data is flushed before
  # julia exits
  # this causes errors to be printed when julia exists, possible #10431
  atexit( () -> if isopen(fbuf.fstream)
                  close(fbuf)
                end
        )
  return fbuf
end

write(io::DebugFileIO, x::UInt8) = write(io.fstream, x)
flush(io::DebugFileIO) = flush(io.fstream)
close(io::DebugFileIO) = close(io.fstream)

import Base.println

function println(io::DebugFileIO, args...)
  tval = MPI.Wtime()
  println(io.fstream, @printf("%16e ", tval), args...)
end
=#

#------------------------------------------------------------------------------
# other random stuff

"""
  This function prints 2D arrays so they look like matrices (m x n format)
"""
function printArray(io::IO, A::AbstractMatrix)

  m, n = size(A)
  for i=1:m
    if i == 1
      print(io, "[")
    else
      print(io, " ")
    end
    for j=1:n
      print(io, A[i, j])

      if j == n
        if i == m
          print(io, "]")
        else
          print(io, ";")
        end
        print(io, "\n")
      else
        print(io, ", ")
      end  # end if
    end  # end j
  end  # end i

  return io
end

function printArray(A::AbstractMatrix)
  return printArray(STDOUT, A)
end