text:
import Base.zero

export
    AV_TIME_BASE,
    AVMediaType,
    AVMEDIA_TYPE_UNKNOWN,
    AVMEDIA_TYPE_VIDEO,
    AVMEDIA_TYPE_AUDIO,
    AVMEDIA_TYPE_DATA,
    AVMEDIA_TYPE_SUBTITLE,
    AVMEDIA_TYPE_ATTACHMENT,
    AVMEDIA_TYPE_NB,
    AVPictureType,
    AV_PICTURE_TYPE_I,
    AV_PICTURE_TYPE_P,
    AV_PICTURE_TYPE_B,
    AV_PICTURE_TYPE_S,
    AV_PICTURE_TYPE_SI,
    AV_PICTURE_TYPE_SP,
    AV_PICTURE_TYPE_BI,
    AV_BUFFER_FLAG_READONLY,
    AVBuffer,
    AVBufferRef,
    AVBufferPool,
    AV_CH_FRONT_LEFT,
    AV_CH_FRONT_RIGHT,
    AV_CH_FRONT_CENTER,
    AV_CH_LOW_FREQUENCY,
    AV_CH_BACK_LEFT,
    AV_CH_BACK_RIGHT,
    AV_CH_FRONT_LEFT_OF_CENTER,
    AV_CH_FRONT_RIGHT_OF_CENTER,
    AV_CH_BACK_CENTER,
    AV_CH_SIDE_LEFT,
    AV_CH_SIDE_RIGHT,
    AV_CH_TOP_CENTER,
    AV_CH_TOP_FRONT_LEFT,
    AV_CH_TOP_FRONT_CENTER,
    AV_CH_TOP_FRONT_RIGHT,
    AV_CH_TOP_BACK_LEFT,
    AV_CH_TOP_BACK_CENTER,
    AV_CH_TOP_BACK_RIGHT,
    AV_CH_STEREO_LEFT,
    AV_CH_STEREO_RIGHT,
    AV_CH_WIDE_LEFT,
    AV_CH_WIDE_RIGHT,
    AV_CH_SURROUND_DIRECT_LEFT,
    AV_CH_SURROUND_DIRECT_RIGHT,
    AV_CH_LOW_FREQUENCY_2,
    AV_CH_LAYOUT_NATIVE,
    AV_CH_LAYOUT_MONO,
    AV_CH_LAYOUT_STEREO,
    AV_CH_LAYOUT_2POINT1,
    AV_CH_LAYOUT_2_1,
    AV_CH_LAYOUT_SURROUND,
    AV_CH_LAYOUT_3POINT1,
    AV_CH_LAYOUT_4POINT0,
    AV_CH_LAYOUT_4POINT1,
    AV_CH_LAYOUT_2_2,
    AV_CH_LAYOUT_QUAD,
    AV_CH_LAYOUT_5POINT0,
    AV_CH_LAYOUT_5POINT1,
    AV_CH_LAYOUT_5POINT0_BACK,
    AV_CH_LAYOUT_5POINT1_BACK,
    AV_CH_LAYOUT_6POINT0,
    AV_CH_LAYOUT_6POINT0_FRONT,
    AV_CH_LAYOUT_HEXAGONAL,
    AV_CH_LAYOUT_6POINT1,
    AV_CH_LAYOUT_6POINT1_BACK,
    AV_CH_LAYOUT_6POINT1_FRONT,
    AV_CH_LAYOUT_7POINT0,
    AV_CH_LAYOUT_7POINT0_FRONT,
    AV_CH_LAYOUT_7POINT1,
    AV_CH_LAYOUT_7POINT1_WIDE,
    AV_CH_LAYOUT_7POINT1_WIDE_BACK,
    AV_CH_LAYOUT_OCTAGONAL,
    AV_CH_LAYOUT_STEREO_DOWNMIX,
    AVMatrixEncoding,
    AV_MATRIX_ENCODING_NONE,
    AV_MATRIX_ENCODING_DOLBY,
    AV_MATRIX_ENCODING_DPLII,
    AV_MATRIX_ENCODING_DPLIIX,
    AV_MATRIX_ENCODING_DPLIIZ,
    AV_MATRIX_ENCODING_DOLBYEX,
    AV_MATRIX_ENCODING_DOLBYHEADPHONE,
    AV_MATRIX_ENCODING_NB,
    AV_DICT_MATCH_CASE,
    AV_DICT_IGNORE_SUFFIX,
    AV_DICT_DONT_STRDUP_KEY,
    AV_DICT_DONT_STRDUP_VAL,
    AV_DICT_DONT_OVERWRITE,
    AV_DICT_APPEND,
    AVDictionaryEntry,
    AVDictionary,
    AVFifoBuffer,
    AV_NUM_DATA_POINTERS,
    AV_FRAME_FLAG_CORRUPT,
    AVFrameSideDataType,
    AV_FRAME_DATA_PANSCAN,
    AV_FRAME_DATA_A53_CC,
    AV_FRAME_DATA_STEREO3D,
    AV_FRAME_DATA_MATRIXENCODING,
    AV_FRAME_DATA_DOWNMIX_INFO,
    AVFrameSideData,
    Array_8_Ptr,
    Array_8_Cint,
    Array_2_Ptr,
    Array_8_Uint64,
    AVRational,
    AVFrame,
    AV_LOG_QUIET,
    AV_LOG_PANIC,
    AV_LOG_FATAL,
    AV_LOG_ERROR,
    AV_LOG_WARNING,
    AV_LOG_INFO,
    AV_LOG_VERBOSE,
    AV_LOG_DEBUG,
    AV_LOG_SKIP_REPEATED,
    AVOptionType,
    AV_OPT_TYPE_FLAGS,
    AV_OPT_TYPE_INT,
    AV_OPT_TYPE_INT64,
    AV_OPT_TYPE_DOUBLE,
    AV_OPT_TYPE_FLOAT,
    AV_OPT_TYPE_STRING,
    AV_OPT_TYPE_RATIONAL,
    AV_OPT_TYPE_BINARY,
    AV_OPT_TYPE_CONST,
    AVOption,
    AVClass,
    AV_PIX_FMT_FLAG_BE,
    AV_PIX_FMT_FLAG_PAL,
    AV_PIX_FMT_FLAG_BITSTREAM,
    AV_PIX_FMT_FLAG_HWACCEL,
    AV_PIX_FMT_FLAG_PLANAR,
    AV_PIX_FMT_FLAG_RGB,
    AV_PIX_FMT_FLAG_PSEUDOPAL,
    AV_PIX_FMT_FLAG_ALPHA,
    PIX_FMT_BE,
    PIX_FMT_PAL,
    PIX_FMT_BITSTREAM,
    PIX_FMT_HWACCEL,
    PIX_FMT_PLANAR,
    PIX_FMT_RGB,
    PIX_FMT_PSEUDOPAL,
    PIX_FMT_ALPHA,
    AVComponentDescriptor,
    Array_4_AVComponentDescriptor,
    AVPixFmtDescriptor,
    AVPixelFormat,
    AV_PIX_FMT_NONE,
    AV_PIX_FMT_YUV420P,
    AV_PIX_FMT_YUYV422,
    AV_PIX_FMT_RGB24,
    AV_PIX_FMT_BGR24,
    AV_PIX_FMT_YUV422P,
    AV_PIX_FMT_YUV444P,
    AV_PIX_FMT_YUV410P,
    AV_PIX_FMT_YUV411P,
    AV_PIX_FMT_GRAY8,
    AV_PIX_FMT_MONOWHITE,
    AV_PIX_FMT_MONOBLACK,
    AV_PIX_FMT_PAL8,
    AV_PIX_FMT_YUVJ420P,
    AV_PIX_FMT_YUVJ422P,
    AV_PIX_FMT_YUVJ444P,
    AV_PIX_FMT_XVMC_MPEG2_MC,
    AV_PIX_FMT_XVMC_MPEG2_IDCT,
    AV_PIX_FMT_UYVY422,
    AV_PIX_FMT_UYYVYY411,
    AV_PIX_FMT_BGR8,
    AV_PIX_FMT_BGR4,
    AV_PIX_FMT_BGR4_BYTE,
    AV_PIX_FMT_RGB8,
    AV_PIX_FMT_RGB4,
    AV_PIX_FMT_RGB4_BYTE,
    AV_PIX_FMT_NV12,
    AV_PIX_FMT_NV21,
    AV_PIX_FMT_ARGB,
    AV_PIX_FMT_RGBA,
    AV_PIX_FMT_ABGR,
    AV_PIX_FMT_BGRA,
    AV_PIX_FMT_GRAY16BE,
    AV_PIX_FMT_GRAY16LE,
    AV_PIX_FMT_YUV440P,
    AV_PIX_FMT_YUVJ440P,
    AV_PIX_FMT_YUVA420P,
    AV_PIX_FMT_VDPAU_H264,
    AV_PIX_FMT_VDPAU_MPEG1,
    AV_PIX_FMT_VDPAU_MPEG2,
    AV_PIX_FMT_VDPAU_WMV3,
    AV_PIX_FMT_VDPAU_VC1,
    AV_PIX_FMT_RGB48BE,
    AV_PIX_FMT_RGB48LE,
    AV_PIX_FMT_RGB565BE,
    AV_PIX_FMT_RGB565LE,
    AV_PIX_FMT_RGB555BE,
    AV_PIX_FMT_RGB555LE,
    AV_PIX_FMT_BGR565BE,
    AV_PIX_FMT_BGR565LE,
    AV_PIX_FMT_BGR555BE,
    AV_PIX_FMT_BGR555LE,
    AV_PIX_FMT_VAAPI_MOCO,
    AV_PIX_FMT_VAAPI_IDCT,
    AV_PIX_FMT_VAAPI_VLD,
    AV_PIX_FMT_YUV420P16LE,
    AV_PIX_FMT_YUV420P16BE,
    AV_PIX_FMT_YUV422P16LE,
    AV_PIX_FMT_YUV422P16BE,
    AV_PIX_FMT_YUV444P16LE,
    AV_PIX_FMT_YUV444P16BE,
    AV_PIX_FMT_VDPAU_MPEG4,
    AV_PIX_FMT_DXVA2_VLD,
    AV_PIX_FMT_RGB444LE,
    AV_PIX_FMT_RGB444BE,
    AV_PIX_FMT_BGR444LE,
    AV_PIX_FMT_BGR444BE,
    AV_PIX_FMT_Y400A,
    AV_PIX_FMT_BGR48BE,
    AV_PIX_FMT_BGR48LE,
    AV_PIX_FMT_YUV420P9BE,
    AV_PIX_FMT_YUV420P9LE,
    AV_PIX_FMT_YUV420P10BE,
    AV_PIX_FMT_YUV420P10LE,
    AV_PIX_FMT_YUV422P10BE,
    AV_PIX_FMT_YUV422P10LE,
    AV_PIX_FMT_YUV444P9BE,
    AV_PIX_FMT_YUV444P9LE,
    AV_PIX_FMT_YUV444P10BE,
    AV_PIX_FMT_YUV444P10LE,
    AV_PIX_FMT_YUV422P9BE,
    AV_PIX_FMT_YUV422P9LE,
    AV_PIX_FMT_VDA_VLD,
    AV_PIX_FMT_GBRP,
    AV_PIX_FMT_GBRP9BE,
    AV_PIX_FMT_GBRP9LE,
    AV_PIX_FMT_GBRP10BE,
    AV_PIX_FMT_GBRP10LE,
    AV_PIX_FMT_GBRP16BE,
    AV_PIX_FMT_GBRP16LE,
    AV_PIX_FMT_YUVA422P,
    AV_PIX_FMT_YUVA444P,
    AV_PIX_FMT_YUVA420P9BE,
    AV_PIX_FMT_YUVA420P9LE,
    AV_PIX_FMT_YUVA422P9BE,
    AV_PIX_FMT_YUVA422P9LE,
    AV_PIX_FMT_YUVA444P9BE,
    AV_PIX_FMT_YUVA444P9LE,
    AV_PIX_FMT_YUVA420P10BE,
    AV_PIX_FMT_YUVA420P10LE,
    AV_PIX_FMT_YUVA422P10BE,
    AV_PIX_FMT_YUVA422P10LE,
    AV_PIX_FMT_YUVA444P10BE,
    AV_PIX_FMT_YUVA444P10LE,
    AV_PIX_FMT_YUVA420P16BE,
    AV_PIX_FMT_YUVA420P16LE,
    AV_PIX_FMT_YUVA422P16BE,
    AV_PIX_FMT_YUVA422P16LE,
    AV_PIX_FMT_YUVA444P16BE,
    AV_PIX_FMT_YUVA444P16LE,
    AV_PIX_FMT_VDPAU,
    AV_PIX_FMT_XYZ12LE,
    AV_PIX_FMT_XYZ12BE,
    AV_PIX_FMT_NV16,
    AV_PIX_FMT_NV20LE,
    AV_PIX_FMT_NV20BE,
    AV_PIX_FMT_NB,
    PIX_FMT_NONE,
    PIX_FMT_YUV420P,
    PIX_FMT_YUYV422,
    PIX_FMT_RGB24,
    PIX_FMT_BGR24,
    PIX_FMT_YUV422P,
    PIX_FMT_YUV444P,
    PIX_FMT_YUV410P,
    PIX_FMT_YUV411P,
    PIX_FMT_GRAY8,
    PIX_FMT_MONOWHITE,
    PIX_FMT_MONOBLACK,
    PIX_FMT_PAL8,
    PIX_FMT_YUVJ420P,
    PIX_FMT_YUVJ422P,
    PIX_FMT_YUVJ444P,
    PIX_FMT_XVMC_MPEG2_MC,
    PIX_FMT_XVMC_MPEG2_IDCT,
    PIX_FMT_UYVY422,
    PIX_FMT_UYYVYY411,
    PIX_FMT_BGR8,
    PIX_FMT_BGR4,
    PIX_FMT_BGR4_BYTE,
    PIX_FMT_RGB8,
    PIX_FMT_RGB4,
    PIX_FMT_RGB4_BYTE,
    PIX_FMT_NV12,
    PIX_FMT_NV21,
    PIX_FMT_ARGB,
    PIX_FMT_RGBA,
    PIX_FMT_ABGR,
    PIX_FMT_BGRA,
    PIX_FMT_GRAY16BE,
    PIX_FMT_GRAY16LE,
    PIX_FMT_YUV440P,
    PIX_FMT_YUVJ440P,
    PIX_FMT_YUVA420P,
    PIX_FMT_VDPAU_H264,
    PIX_FMT_VDPAU_MPEG1,
    PIX_FMT_VDPAU_MPEG2,
    PIX_FMT_VDPAU_WMV3,
    PIX_FMT_VDPAU_VC1,
    PIX_FMT_RGB48BE,
    PIX_FMT_RGB48LE,
    PIX_FMT_RGB565BE,
    PIX_FMT_RGB565LE,
    PIX_FMT_RGB555BE,
    PIX_FMT_RGB555LE,
    PIX_FMT_BGR565BE,
    PIX_FMT_BGR565LE,
    PIX_FMT_BGR555BE,
    PIX_FMT_BGR555LE,
    PIX_FMT_VAAPI_MOCO,
    PIX_FMT_VAAPI_IDCT,
    PIX_FMT_VAAPI_VLD,
    PIX_FMT_YUV420P16LE,
    PIX_FMT_YUV420P16BE,
    PIX_FMT_YUV422P16LE,
    PIX_FMT_YUV422P16BE,
    PIX_FMT_YUV444P16LE,
    PIX_FMT_YUV444P16BE,
    PIX_FMT_VDPAU_MPEG4,
    PIX_FMT_DXVA2_VLD,
    PIX_FMT_RGB444LE,
    PIX_FMT_RGB444BE,
    PIX_FMT_BGR444LE,
    PIX_FMT_BGR444BE,
    PIX_FMT_Y400A,
    PIX_FMT_BGR48BE,
    PIX_FMT_BGR48LE,
    PIX_FMT_YUV420P9BE,
    PIX_FMT_YUV420P9LE,
    PIX_FMT_YUV420P10BE,
    PIX_FMT_YUV420P10LE,
    PIX_FMT_YUV422P10BE,
    PIX_FMT_YUV422P10LE,
    PIX_FMT_YUV444P9BE,
    PIX_FMT_YUV444P9LE,
    PIX_FMT_YUV444P10BE,
    PIX_FMT_YUV444P10LE,
    PIX_FMT_YUV422P9BE,
    PIX_FMT_YUV422P9LE,
    PIX_FMT_VDA_VLD,
    PIX_FMT_GBRP,
    PIX_FMT_GBRP9BE,
    PIX_FMT_GBRP9LE,
    PIX_FMT_GBRP10BE,
    PIX_FMT_GBRP10LE,
    PIX_FMT_GBRP16BE,
    PIX_FMT_GBRP16LE,
    PIX_FMT_NB,
    PixelFormat,
    PIX_FMT_RGB32,
    PIX_FMT_RGB32_1,
    PIX_FMT_BGR32,
    PIX_FMT_BGR32_1,
    PIX_FMT_GRAY16,
    PIX_FMT_RGB48,
    PIX_FMT_RGB565,
    PIX_FMT_RGB555,
    PIX_FMT_RGB444,
    PIX_FMT_BGR48,
    PIX_FMT_BGR565,
    PIX_FMT_BGR555,
    PIX_FMT_BGR444,
    PIX_FMT_YUV420P9,
    PIX_FMT_YUV422P9,
    PIX_FMT_YUV444P9,
    PIX_FMT_YUV420P10,
    PIX_FMT_YUV422P10,
    PIX_FMT_YUV444P10,
    PIX_FMT_YUV420P16,
    PIX_FMT_YUV422P16,
    PIX_FMT_YUV444P16,
    PIX_FMT_GBRP9,
    PIX_FMT_GBRP10,
    PIX_FMT_GBRP16,
    AVSampleFormat,
    AV_SAMPLE_FMT_NONE,
    AV_SAMPLE_FMT_U8,
    AV_SAMPLE_FMT_S16,
    AV_SAMPLE_FMT_S32,
    AV_SAMPLE_FMT_FLT,
    AV_SAMPLE_FMT_DBL,
    AV_SAMPLE_FMT_U8P,
    AV_SAMPLE_FMT_S16P,
    AV_SAMPLE_FMT_S32P,
    AV_SAMPLE_FMT_FLTP,
    AV_SAMPLE_FMT_DBLP,
    AV_SAMPLE_FMT_NB,
    LIBAVUTIL_VERSION_MAJOR,
    LIBAVUTIL_VERSION_MINOR,
    LIBAVUTIL_VERSION_MICRO,
    LIBAVUTIL_BUILD,
    AVAudioFifo,
    AVDownmixType,
    AV_DOWNMIX_TYPE_UNKNOWN,
    AV_DOWNMIX_TYPE_LORO,
    AV_DOWNMIX_TYPE_LTRT,
    AV_DOWNMIX_TYPE_DPLII,
    AV_DOWNMIX_TYPE_NB,
    AVDownmixInfo,
    AV_OPT_FLAG_ENCODING_PARAM,
    AV_OPT_FLAG_DECODING_PARAM,
    AV_OPT_FLAG_METADATA,
    AV_OPT_FLAG_AUDIO_PARAM,
    AV_OPT_FLAG_VIDEO_PARAM,
    AV_OPT_FLAG_SUBTITLE_PARAM,
    AV_OPT_SEARCH_CHILDREN,
    AV_OPT_SEARCH_FAKE_OBJ,
    AV_STEREO3D_FLAG_INVERT,
    AVStereo3DType,
    AV_STEREO3D_2D,
    AV_STEREO3D_SIDEBYSIDE,
    AV_STEREO3D_TOPBOTTOM,
    AV_STEREO3D_FRAMESEQUENCE,
    AV_STEREO3D_CHECKERBOARD,
    AV_STEREO3D_SIDEBYSIDE_QUINCUNX,
    AV_STEREO3D_LINES,
    AV_STEREO3D_COLUMNS,
    AVStereo3D,
    Array_16_Uint32,
    AVXTEA

const FF_LAMBDA_SHIFT = 7
const FF_LAMBDA_SCALE = 1 << FF_LAMBDA_SHIFT
const FF_QP2LAMBDA = 118
const FF_LAMBDA_MAX = 256 * 128 - 1
const FF_QUALITY_SCALE = FF_LAMBDA_SCALE

const AV_NOPTS_VALUE = reinterpret(Int64, 0x8000000000000000)

const AV_TIME_BASE = 1000000

immutable AVRational
    num::Cint
    den::Cint
end

const AV_TIME_BASE_Q = AVRational( 1 , AV_TIME_BASE )

# begin enum AVMediaType
const AVMediaType=Cint
const AVMEDIA_TYPE_UNKNOWN = Int32(-1)
const AVMEDIA_TYPE_VIDEO = Int32(0)
const AVMEDIA_TYPE_AUDIO = Int32(1)
const AVMEDIA_TYPE_DATA = Int32(2)
const AVMEDIA_TYPE_SUBTITLE = Int32(3)
const AVMEDIA_TYPE_ATTACHMENT = Int32(4)
const AVMEDIA_TYPE_NB = Int32(5)
# end enum AVMediaType

# begin enum AVPictureType
const AVPictureType=UInt32
const AV_PICTURE_TYPE_I = UInt32(1)
const AV_PICTURE_TYPE_P = UInt32(2)
const AV_PICTURE_TYPE_B = UInt32(3)
const AV_PICTURE_TYPE_S = UInt32(4)
const AV_PICTURE_TYPE_SI = UInt32(5)
const AV_PICTURE_TYPE_SP = UInt32(6)
const AV_PICTURE_TYPE_BI = UInt32(7)
# end enum AVPictureType

const AV_BUFFER_FLAG_READONLY = 1 << 0

const AVBuffer=Void

immutable AVBufferRef
    buffer::Ptr{AVBuffer}
    data::Ptr{UInt8}
    size::Cint
end

const AVBufferPool=Void

const AV_CH_FRONT_LEFT = 0x00000001
const AV_CH_FRONT_RIGHT = 0x00000002
const AV_CH_FRONT_CENTER = 0x00000004
const AV_CH_LOW_FREQUENCY = 0x00000008
const AV_CH_BACK_LEFT = 0x00000010
const AV_CH_BACK_RIGHT = 0x00000020
const AV_CH_FRONT_LEFT_OF_CENTER = 0x00000040
const AV_CH_FRONT_RIGHT_OF_CENTER = 0x00000080
const AV_CH_BACK_CENTER = 0x00000100
const AV_CH_SIDE_LEFT = 0x00000200
const AV_CH_SIDE_RIGHT = 0x00000400
const AV_CH_TOP_CENTER = 0x00000800
const AV_CH_TOP_FRONT_LEFT = 0x00001000
const AV_CH_TOP_FRONT_CENTER = 0x00002000
const AV_CH_TOP_FRONT_RIGHT = 0x00004000
const AV_CH_TOP_BACK_LEFT = 0x00008000
const AV_CH_TOP_BACK_CENTER = 0x00010000
const AV_CH_TOP_BACK_RIGHT = 0x00020000
const AV_CH_STEREO_LEFT = 0x20000000
const AV_CH_STEREO_RIGHT = 0x40000000
const AV_CH_WIDE_LEFT = 0x0000000080000000
const AV_CH_WIDE_RIGHT = 0x0000000100000000
const AV_CH_SURROUND_DIRECT_LEFT = 0x0000000200000000
const AV_CH_SURROUND_DIRECT_RIGHT = 0x0000000400000000
const AV_CH_LOW_FREQUENCY_2 = 0x0000000800000000
const AV_CH_LAYOUT_NATIVE = 0x8000000000000000
const AV_CH_LAYOUT_MONO = AV_CH_FRONT_CENTER
const AV_CH_LAYOUT_STEREO = AV_CH_FRONT_LEFT | AV_CH_FRONT_RIGHT
const AV_CH_LAYOUT_2POINT1 = AV_CH_LAYOUT_STEREO | AV_CH_LOW_FREQUENCY
const AV_CH_LAYOUT_2_1 = AV_CH_LAYOUT_STEREO | AV_CH_BACK_CENTER
const AV_CH_LAYOUT_SURROUND = AV_CH_LAYOUT_STEREO | AV_CH_FRONT_CENTER
const AV_CH_LAYOUT_3POINT1 = AV_CH_LAYOUT_SURROUND | AV_CH_LOW_FREQUENCY
const AV_CH_LAYOUT_4POINT0 = AV_CH_LAYOUT_SURROUND | AV_CH_BACK_CENTER
const AV_CH_LAYOUT_4POINT1 = AV_CH_LAYOUT_4POINT0 | AV_CH_LOW_FREQUENCY
const AV_CH_LAYOUT_2_2 = (AV_CH_LAYOUT_STEREO | AV_CH_SIDE_LEFT) | AV_CH_SIDE_RIGHT
const AV_CH_LAYOUT_QUAD = (AV_CH_LAYOUT_STEREO | AV_CH_BACK_LEFT) | AV_CH_BACK_RIGHT
const AV_CH_LAYOUT_5POINT0 = (AV_CH_LAYOUT_SURROUND | AV_CH_SIDE_LEFT) | AV_CH_SIDE_RIGHT
const AV_CH_LAYOUT_5POINT1 = AV_CH_LAYOUT_5POINT0 | AV_CH_LOW_FREQUENCY
const AV_CH_LAYOUT_5POINT0_BACK = (AV_CH_LAYOUT_SURROUND | AV_CH_BACK_LEFT) | AV_CH_BACK_RIGHT
const AV_CH_LAYOUT_5POINT1_BACK = AV_CH_LAYOUT_5POINT0_BACK | AV_CH_LOW_FREQUENCY
const AV_CH_LAYOUT_6POINT0 = AV_CH_LAYOUT_5POINT0 | AV_CH_BACK_CENTER
const AV_CH_LAYOUT_6POINT0_FRONT = (AV_CH_LAYOUT_2_2 | AV_CH_FRONT_LEFT_OF_CENTER) | AV_CH_FRONT_RIGHT_OF_CENTER
const AV_CH_LAYOUT_HEXAGONAL = AV_CH_LAYOUT_5POINT0_BACK | AV_CH_BACK_CENTER
const AV_CH_LAYOUT_6POINT1 = AV_CH_LAYOUT_5POINT1 | AV_CH_BACK_CENTER
const AV_CH_LAYOUT_6POINT1_BACK = AV_CH_LAYOUT_5POINT1_BACK | AV_CH_BACK_CENTER
const AV_CH_LAYOUT_6POINT1_FRONT = AV_CH_LAYOUT_6POINT0_FRONT | AV_CH_LOW_FREQUENCY
const AV_CH_LAYOUT_7POINT0 = (AV_CH_LAYOUT_5POINT0 | AV_CH_BACK_LEFT) | AV_CH_BACK_RIGHT
const AV_CH_LAYOUT_7POINT0_FRONT = (AV_CH_LAYOUT_5POINT0 | AV_CH_FRONT_LEFT_OF_CENTER) | AV_CH_FRONT_RIGHT_OF_CENTER
const AV_CH_LAYOUT_7POINT1 = (AV_CH_LAYOUT_5POINT1 | AV_CH_BACK_LEFT) | AV_CH_BACK_RIGHT
const AV_CH_LAYOUT_7POINT1_WIDE = (AV_CH_LAYOUT_5POINT1 | AV_CH_FRONT_LEFT_OF_CENTER) | AV_CH_FRONT_RIGHT_OF_CENTER
const AV_CH_LAYOUT_7POINT1_WIDE_BACK = (AV_CH_LAYOUT_5POINT1_BACK | AV_CH_FRONT_LEFT_OF_CENTER) | AV_CH_FRONT_RIGHT_OF_CENTER
const AV_CH_LAYOUT_OCTAGONAL = ((AV_CH_LAYOUT_5POINT0 | AV_CH_BACK_LEFT) | AV_CH_BACK_CENTER) | AV_CH_BACK_RIGHT
const AV_CH_LAYOUT_STEREO_DOWNMIX = AV_CH_STEREO_LEFT | AV_CH_STEREO_RIGHT

# begin enum AVMatrixEncoding
const AVMatrixEncoding=UInt32
const AV_MATRIX_ENCODING_NONE = UInt32(0)
const AV_MATRIX_ENCODING_DOLBY = UInt32(1)
const AV_MATRIX_ENCODING_DPLII = UInt32(2)
const AV_MATRIX_ENCODING_DPLIIX = UInt32(3)
const AV_MATRIX_ENCODING_DPLIIZ = UInt32(4)
const AV_MATRIX_ENCODING_DOLBYEX = UInt32(5)
const AV_MATRIX_ENCODING_DOLBYHEADPHONE = UInt32(6)
const AV_MATRIX_ENCODING_NB = UInt32(7)
# end enum AVMatrixEncoding

const AV_DICT_MATCH_CASE = 1
const AV_DICT_IGNORE_SUFFIX = 2
const AV_DICT_DONT_STRDUP_KEY = 4
const AV_DICT_DONT_STRDUP_VAL = 8
const AV_DICT_DONT_OVERWRITE = 16
const AV_DICT_APPEND = 32

immutable AVDictionaryEntry
    key::Ptr{UInt8}
    value::Ptr{UInt8}
end

const AVDictionary=Void

immutable AVFifoBuffer
    buffer::Ptr{UInt8}
    rptr::Ptr{UInt8}
    wptr::Ptr{UInt8}
    _end::Ptr{UInt8}
    rndx::UInt32
    wndx::UInt32
end

const AV_NUM_DATA_POINTERS = 8
const AV_FRAME_FLAG_CORRUPT = 1 << 0

# begin enum AVFrameSideDataType
const AVFrameSideDataType=UInt32
const AV_FRAME_DATA_PANSCAN = UInt32(0)
const AV_FRAME_DATA_A53_CC = UInt32(1)
const AV_FRAME_DATA_STEREO3D = UInt32(2)
const AV_FRAME_DATA_MATRIXENCODING = UInt32(3)
const AV_FRAME_DATA_DOWNMIX_INFO = UInt32(4)
# end enum AVFrameSideDataType

immutable AVFrameSideData
    _type::AVFrameSideDataType
    data::Ptr{UInt8}
    size::Cint
    metadata::Ptr{AVDictionary}
end

immutable Array_8_Ptr
    d1::Ptr{UInt8}
    d2::Ptr{UInt8}
    d3::Ptr{UInt8}
    d4::Ptr{UInt8}
    d5::Ptr{UInt8}
    d6::Ptr{UInt8}
    d7::Ptr{UInt8}
    d8::Ptr{UInt8}
end

zero(::Type{Array_8_Ptr}) = Array_8_Ptr(fill(C_NULL,8)...)

immutable Array_8_Cint
    d1::Cint
    d2::Cint
    d3::Cint
    d4::Cint
    d5::Cint
    d6::Cint
    d7::Cint
    d8::Cint
end

zero(::Type{Array_8_Cint}) = Array_8_Cint(fill(zero(Cint),8)...)

immutable Array_2_Ptr
    d1::Ptr{Void}
    d2::Ptr{Void}
end

zero(::Type{Array_2_Ptr}) = Array_2_Ptr(fill(C_NULL,2)...)

immutable Array_8_Uint64
    d1::UInt64
    d2::UInt64
    d3::UInt64
    d4::UInt64
    d5::UInt64
    d6::UInt64
    d7::UInt64
    d8::UInt64
end

zero(::Type{Array_8_Uint64}) = Array_8_Uint64(fill(zero(UInt64),8)...)

immutable AVFrame
    data::Array_8_Ptr
    linesize::Array_8_Cint
    extended_data::Ptr{Ptr{UInt8}}
    width::Cint
    height::Cint
    nb_samples::Cint
    format::Cint
    key_frame::Cint
    pict_type::AVPictureType
    base::Array_8_Ptr
    sample_aspect_ratio::AVRational
    pts::Int64
    pkt_pts::Int64
    pkt_dts::Int64
    coded_picture_number::Cint
    display_picture_number::Cint
    quality::Cint
    reference::Cint
    qscale_table::Ptr{Int8}
    qstride::Cint
    qscale_type::Cint
    mbskip_table::Ptr{UInt8}
    motion_val::Array_2_Ptr
    mb_type::Ptr{UInt32}
    dct_coeff::Ptr{Int16}
    ref_index::Array_2_Ptr
    opaque::Ptr{Void}
    error::Array_8_Uint64
    _type::Cint
    repeat_pict::Cint
    interlaced_frame::Cint
    top_field_first::Cint
    palette_has_changed::Cint
    buffer_hints::Cint
    pan_scan::Ptr{Void}  #Ptr{AVPanScan}
    reordered_opaque::Int64
    hwaccel_picture_private::Ptr{Void}
    owner::Ptr{Void} #Ptr{AVCodecContext}
    thread_opaque::Ptr{Void}
    motion_subsample_log2::UInt8
    sample_rate::Cint
    channel_layout::UInt64
    buf::Array_8_Ptr
    extended_buf::Ptr{Ptr{AVBufferRef}}
    nb_extended_buf::Cint
    side_data::Ptr{Ptr{AVFrameSideData}}
    nb_side_data::Cint
    flags::Cint
end

const AV_LOG_QUIET = -8
const AV_LOG_PANIC = 0
const AV_LOG_FATAL = 8
const AV_LOG_ERROR = 16
const AV_LOG_WARNING = 24
const AV_LOG_INFO = 32
const AV_LOG_VERBOSE = 40
const AV_LOG_DEBUG = 48

# Skipping MacroDefinition: av_dlog ( pctx , ... )

const AV_LOG_SKIP_REPEATED = 1

# begin enum AVOptionType
const AVOptionType=UInt32
const AV_OPT_TYPE_FLAGS = UInt32(0)
const AV_OPT_TYPE_INT = UInt32(1)
const AV_OPT_TYPE_INT64 = UInt32(2)
const AV_OPT_TYPE_DOUBLE = UInt32(3)
const AV_OPT_TYPE_FLOAT = UInt32(4)
const AV_OPT_TYPE_STRING = UInt32(5)
const AV_OPT_TYPE_RATIONAL = UInt32(6)
const AV_OPT_TYPE_BINARY = UInt32(7)
const AV_OPT_TYPE_CONST = UInt32(128)
# end enum AVOptionType

immutable AVOption
    name::Ptr{UInt8}
    help::Ptr{UInt8}
    offset::Cint
    _type::AVOptionType
    default_val::Void
    min::Cdouble
    max::Cdouble
    flags::Cint
    unit::Ptr{UInt8}
end

immutable AVClass
    class_name::Ptr{UInt8}
    item_name::Ptr{Void}
    option::Ptr{AVOption}
    version::Cint
    log_level_offset_offset::Cint
    parent_log_context_offset::Cint
    child_next::Ptr{Void}
    child_class_next::Ptr{Void}
end

const AV_PIX_FMT_FLAG_BE = 1 << 0
const AV_PIX_FMT_FLAG_PAL = 1 << 1
const AV_PIX_FMT_FLAG_BITSTREAM = 1 << 2
const AV_PIX_FMT_FLAG_HWACCEL = 1 << 3
const AV_PIX_FMT_FLAG_PLANAR = 1 << 4
const AV_PIX_FMT_FLAG_RGB = 1 << 5
const AV_PIX_FMT_FLAG_PSEUDOPAL = 1 << 6
const AV_PIX_FMT_FLAG_ALPHA = 1 << 7
const PIX_FMT_BE = AV_PIX_FMT_FLAG_BE
const PIX_FMT_PAL = AV_PIX_FMT_FLAG_PAL
const PIX_FMT_BITSTREAM = AV_PIX_FMT_FLAG_BITSTREAM
const PIX_FMT_HWACCEL = AV_PIX_FMT_FLAG_HWACCEL
const PIX_FMT_PLANAR = AV_PIX_FMT_FLAG_PLANAR
const PIX_FMT_RGB = AV_PIX_FMT_FLAG_RGB
const PIX_FMT_PSEUDOPAL = AV_PIX_FMT_FLAG_PSEUDOPAL
const PIX_FMT_ALPHA = AV_PIX_FMT_FLAG_ALPHA

immutable AVComponentDescriptor
    plane::UInt16
    step_minus1::UInt16
    offset_plus1::UInt16
    shift::UInt16
    depth_minus1::UInt16
end

immutable Array_4_AVComponentDescriptor
    d1::AVComponentDescriptor
    d2::AVComponentDescriptor
    d3::AVComponentDescriptor
    d4::AVComponentDescriptor
end

zero(::Type{Array_4_AVComponentDescriptor}) = Array_4_AVComponentDescriptor(fill(zero(AVComponentDescriptor),4)...)

immutable AVPixFmtDescriptor
    name::Ptr{UInt8}
    nb_components::UInt8
    log2_chroma_w::UInt8
    log2_chroma_h::UInt8
    flags::UInt8
    comp::Array_4_AVComponentDescriptor
end

# begin enum AVPixelFormat
const AVPixelFormat=Cint
const AV_PIX_FMT_NONE = Int32(-1)
const AV_PIX_FMT_YUV420P = Int32(0)
const AV_PIX_FMT_YUYV422 = Int32(1)
const AV_PIX_FMT_RGB24 = Int32(2)
const AV_PIX_FMT_BGR24 = Int32(3)
const AV_PIX_FMT_YUV422P = Int32(4)
const AV_PIX_FMT_YUV444P = Int32(5)
const AV_PIX_FMT_YUV410P = Int32(6)
const AV_PIX_FMT_YUV411P = Int32(7)
const AV_PIX_FMT_GRAY8 = Int32(8)
const AV_PIX_FMT_MONOWHITE = Int32(9)
const AV_PIX_FMT_MONOBLACK = Int32(10)
const AV_PIX_FMT_PAL8 = Int32(11)
const AV_PIX_FMT_YUVJ420P = Int32(12)
const AV_PIX_FMT_YUVJ422P = Int32(13)
const AV_PIX_FMT_YUVJ444P = Int32(14)
const AV_PIX_FMT_XVMC_MPEG2_MC = Int32(15)
const AV_PIX_FMT_XVMC_MPEG2_IDCT = Int32(16)
const AV_PIX_FMT_UYVY422 = Int32(17)
const AV_PIX_FMT_UYYVYY411 = Int32(18)
const AV_PIX_FMT_BGR8 = Int32(19)
const AV_PIX_FMT_BGR4 = Int32(20)
const AV_PIX_FMT_BGR4_BYTE = Int32(21)
const AV_PIX_FMT_RGB8 = Int32(22)
const AV_PIX_FMT_RGB4 = Int32(23)
const AV_PIX_FMT_RGB4_BYTE = Int32(24)
const AV_PIX_FMT_NV12 = Int32(25)
const AV_PIX_FMT_NV21 = Int32(26)
const AV_PIX_FMT_ARGB = Int32(27)
const AV_PIX_FMT_RGBA = Int32(28)
const AV_PIX_FMT_ABGR = Int32(29)
const AV_PIX_FMT_BGRA = Int32(30)
const AV_PIX_FMT_GRAY16BE = Int32(31)
const AV_PIX_FMT_GRAY16LE = Int32(32)
const AV_PIX_FMT_YUV440P = Int32(33)
const AV_PIX_FMT_YUVJ440P = Int32(34)
const AV_PIX_FMT_YUVA420P = Int32(35)
const AV_PIX_FMT_VDPAU_H264 = Int32(36)
const AV_PIX_FMT_VDPAU_MPEG1 = Int32(37)
const AV_PIX_FMT_VDPAU_MPEG2 = Int32(38)
const AV_PIX_FMT_VDPAU_WMV3 = Int32(39)
const AV_PIX_FMT_VDPAU_VC1 = Int32(40)
const AV_PIX_FMT_RGB48BE = Int32(41)
const AV_PIX_FMT_RGB48LE = Int32(42)
const AV_PIX_FMT_RGB565BE = Int32(43)
const AV_PIX_FMT_RGB565LE = Int32(44)
const AV_PIX_FMT_RGB555BE = Int32(45)
const AV_PIX_FMT_RGB555LE = Int32(46)
const AV_PIX_FMT_BGR565BE = Int32(47)
const AV_PIX_FMT_BGR565LE = Int32(48)
const AV_PIX_FMT_BGR555BE = Int32(49)
const AV_PIX_FMT_BGR555LE = Int32(50)
const AV_PIX_FMT_VAAPI_MOCO = Int32(51)
const AV_PIX_FMT_VAAPI_IDCT = Int32(52)
const AV_PIX_FMT_VAAPI_VLD = Int32(53)
const AV_PIX_FMT_YUV420P16LE = Int32(54)
const AV_PIX_FMT_YUV420P16BE = Int32(55)
const AV_PIX_FMT_YUV422P16LE = Int32(56)
const AV_PIX_FMT_YUV422P16BE = Int32(57)
const AV_PIX_FMT_YUV444P16LE = Int32(58)
const AV_PIX_FMT_YUV444P16BE = Int32(59)
const AV_PIX_FMT_VDPAU_MPEG4 = Int32(60)
const AV_PIX_FMT_DXVA2_VLD = Int32(61)
const AV_PIX_FMT_RGB444LE = Int32(62)
const AV_PIX_FMT_RGB444BE = Int32(63)
const AV_PIX_FMT_BGR444LE = Int32(64)
const AV_PIX_FMT_BGR444BE = Int32(65)
const AV_PIX_FMT_Y400A = Int32(66)
const AV_PIX_FMT_BGR48BE = Int32(67)
const AV_PIX_FMT_BGR48LE = Int32(68)
const AV_PIX_FMT_YUV420P9BE = Int32(69)
const AV_PIX_FMT_YUV420P9LE = Int32(70)
const AV_PIX_FMT_YUV420P10BE = Int32(71)
const AV_PIX_FMT_YUV420P10LE = Int32(72)
const AV_PIX_FMT_YUV422P10BE = Int32(73)
const AV_PIX_FMT_YUV422P10LE = Int32(74)
const AV_PIX_FMT_YUV444P9BE = Int32(75)
const AV_PIX_FMT_YUV444P9LE = Int32(76)
const AV_PIX_FMT_YUV444P10BE = Int32(77)
const AV_PIX_FMT_YUV444P10LE = Int32(78)
const AV_PIX_FMT_YUV422P9BE = Int32(79)
const AV_PIX_FMT_YUV422P9LE = Int32(80)
const AV_PIX_FMT_VDA_VLD = Int32(81)
const AV_PIX_FMT_GBRP = Int32(82)
const AV_PIX_FMT_GBRP9BE = Int32(83)
const AV_PIX_FMT_GBRP9LE = Int32(84)
const AV_PIX_FMT_GBRP10BE = Int32(85)
const AV_PIX_FMT_GBRP10LE = Int32(86)
const AV_PIX_FMT_GBRP16BE = Int32(87)
const AV_PIX_FMT_GBRP16LE = Int32(88)
const AV_PIX_FMT_YUVA422P = Int32(89)
const AV_PIX_FMT_YUVA444P = Int32(90)
const AV_PIX_FMT_YUVA420P9BE = Int32(91)
const AV_PIX_FMT_YUVA420P9LE = Int32(92)
const AV_PIX_FMT_YUVA422P9BE = Int32(93)
const AV_PIX_FMT_YUVA422P9LE = Int32(94)
const AV_PIX_FMT_YUVA444P9BE = Int32(95)
const AV_PIX_FMT_YUVA444P9LE = Int32(96)
const AV_PIX_FMT_YUVA420P10BE = Int32(97)
const AV_PIX_FMT_YUVA420P10LE = Int32(98)
const AV_PIX_FMT_YUVA422P10BE = Int32(99)
const AV_PIX_FMT_YUVA422P10LE = Int32(100)
const AV_PIX_FMT_YUVA444P10BE = Int32(101)
const AV_PIX_FMT_YUVA444P10LE = Int32(102)
const AV_PIX_FMT_YUVA420P16BE = Int32(103)
const AV_PIX_FMT_YUVA420P16LE = Int32(104)
const AV_PIX_FMT_YUVA422P16BE = Int32(105)
const AV_PIX_FMT_YUVA422P16LE = Int32(106)
const AV_PIX_FMT_YUVA444P16BE = Int32(107)
const AV_PIX_FMT_YUVA444P16LE = Int32(108)
const AV_PIX_FMT_VDPAU = Int32(109)
const AV_PIX_FMT_XYZ12LE = Int32(110)
const AV_PIX_FMT_XYZ12BE = Int32(111)
const AV_PIX_FMT_NV16 = Int32(112)
const AV_PIX_FMT_NV20LE = Int32(113)
const AV_PIX_FMT_NV20BE = Int32(114)
const AV_PIX_FMT_NB = Int32(115)
const PIX_FMT_NONE = Int32(-1)
const PIX_FMT_YUV420P = Int32(0)
const PIX_FMT_YUYV422 = Int32(1)
const PIX_FMT_RGB24 = Int32(2)
const PIX_FMT_BGR24 = Int32(3)
const PIX_FMT_YUV422P = Int32(4)
const PIX_FMT_YUV444P = Int32(5)
const PIX_FMT_YUV410P = Int32(6)
const PIX_FMT_YUV411P = Int32(7)
const PIX_FMT_GRAY8 = Int32(8)
const PIX_FMT_MONOWHITE = Int32(9)
const PIX_FMT_MONOBLACK = Int32(10)
const PIX_FMT_PAL8 = Int32(11)
const PIX_FMT_YUVJ420P = Int32(12)
const PIX_FMT_YUVJ422P = Int32(13)
const PIX_FMT_YUVJ444P = Int32(14)
const PIX_FMT_XVMC_MPEG2_MC = Int32(15)
const PIX_FMT_XVMC_MPEG2_IDCT = Int32(16)
const PIX_FMT_UYVY422 = Int32(17)
const PIX_FMT_UYYVYY411 = Int32(18)
const PIX_FMT_BGR8 = Int32(19)
const PIX_FMT_BGR4 = Int32(20)
const PIX_FMT_BGR4_BYTE = Int32(21)
const PIX_FMT_RGB8 = Int32(22)
const PIX_FMT_RGB4 = Int32(23)
const PIX_FMT_RGB4_BYTE = Int32(24)
const PIX_FMT_NV12 = Int32(25)
const PIX_FMT_NV21 = Int32(26)
const PIX_FMT_ARGB = Int32(27)
const PIX_FMT_RGBA = Int32(28)
const PIX_FMT_ABGR = Int32(29)
const PIX_FMT_BGRA = Int32(30)
const PIX_FMT_GRAY16BE = Int32(31)
const PIX_FMT_GRAY16LE = Int32(32)
const PIX_FMT_YUV440P = Int32(33)
const PIX_FMT_YUVJ440P = Int32(34)
const PIX_FMT_YUVA420P = Int32(35)
const PIX_FMT_VDPAU_H264 = Int32(36)
const PIX_FMT_VDPAU_MPEG1 = Int32(37)
const PIX_FMT_VDPAU_MPEG2 = Int32(38)
const PIX_FMT_VDPAU_WMV3 = Int32(39)
const PIX_FMT_VDPAU_VC1 = Int32(40)
const PIX_FMT_RGB48BE = Int32(41)
const PIX_FMT_RGB48LE = Int32(42)
const PIX_FMT_RGB565BE = Int32(43)
const PIX_FMT_RGB565LE = Int32(44)
const PIX_FMT_RGB555BE = Int32(45)
const PIX_FMT_RGB555LE = Int32(46)
const PIX_FMT_BGR565BE = Int32(47)
const PIX_FMT_BGR565LE = Int32(48)
const PIX_FMT_BGR555BE = Int32(49)
const PIX_FMT_BGR555LE = Int32(50)
const PIX_FMT_VAAPI_MOCO = Int32(51)
const PIX_FMT_VAAPI_IDCT = Int32(52)
const PIX_FMT_VAAPI_VLD = Int32(53)
const PIX_FMT_YUV420P16LE = Int32(54)
const PIX_FMT_YUV420P16BE = Int32(55)
const PIX_FMT_YUV422P16LE = Int32(56)
const PIX_FMT_YUV422P16BE = Int32(57)
const PIX_FMT_YUV444P16LE = Int32(58)
const PIX_FMT_YUV444P16BE = Int32(59)
const PIX_FMT_VDPAU_MPEG4 = Int32(60)
const PIX_FMT_DXVA2_VLD = Int32(61)
const PIX_FMT_RGB444LE = Int32(62)
const PIX_FMT_RGB444BE = Int32(63)
const PIX_FMT_BGR444LE = Int32(64)
const PIX_FMT_BGR444BE = Int32(65)
const PIX_FMT_Y400A = Int32(66)
const PIX_FMT_BGR48BE = Int32(67)
const PIX_FMT_BGR48LE = Int32(68)
const PIX_FMT_YUV420P9BE = Int32(69)
const PIX_FMT_YUV420P9LE = Int32(70)
const PIX_FMT_YUV420P10BE = Int32(71)
const PIX_FMT_YUV420P10LE = Int32(72)
const PIX_FMT_YUV422P10BE = Int32(73)
const PIX_FMT_YUV422P10LE = Int32(74)
const PIX_FMT_YUV444P9BE = Int32(75)
const PIX_FMT_YUV444P9LE = Int32(76)
const PIX_FMT_YUV444P10BE = Int32(77)
const PIX_FMT_YUV444P10LE = Int32(78)
const PIX_FMT_YUV422P9BE = Int32(79)
const PIX_FMT_YUV422P9LE = Int32(80)
const PIX_FMT_VDA_VLD = Int32(81)
const PIX_FMT_GBRP = Int32(82)
const PIX_FMT_GBRP9BE = Int32(83)
const PIX_FMT_GBRP9LE = Int32(84)
const PIX_FMT_GBRP10BE = Int32(85)
const PIX_FMT_GBRP10LE = Int32(86)
const PIX_FMT_GBRP16BE = Int32(87)
const PIX_FMT_GBRP16LE = Int32(88)
const PIX_FMT_NB = Int32(89)
# end enum AVPixelFormat

const PixelFormat = AVPixelFormat

macro AV_PIX_FMT_NE(be,le)
    Symbol("AV_PIX_FMT_"*string(le))
end

const AV_PIX_FMT_RGB32 = @AV_PIX_FMT_NE(ARGB,BGRA)
const AV_PIX_FMT_RGB32_1 = @AV_PIX_FMT_NE(RGBA,ABGR)
const AV_PIX_FMT_BGR32 = @AV_PIX_FMT_NE(ABGR,RGBA)
const AV_PIX_FMT_BGR32_1 = @AV_PIX_FMT_NE(BGRA,ARGB)
const AV_PIX_FMT_GRAY16 = @AV_PIX_FMT_NE(GRAY16BE,GRAY16LE)
const AV_PIX_FMT_RGB48 = @AV_PIX_FMT_NE(RGB48BE,RGB48LE)
const AV_PIX_FMT_RGB565 = @AV_PIX_FMT_NE(RGB565BE,RGB565LE)
const AV_PIX_FMT_RGB555 = @AV_PIX_FMT_NE(RGB555BE,RGB555LE)
const AV_PIX_FMT_RGB444 = @AV_PIX_FMT_NE(RGB444BE,RGB444LE)
const AV_PIX_FMT_BGR48 = @AV_PIX_FMT_NE(BGR48BE,BGR48LE)
const AV_PIX_FMT_BGR565 = @AV_PIX_FMT_NE(BGR565BE,BGR565LE)
const AV_PIX_FMT_BGR555 = @AV_PIX_FMT_NE(BGR555BE,BGR555LE)
const AV_PIX_FMT_BGR444 = @AV_PIX_FMT_NE(BGR444BE,BGR444LE)
const AV_PIX_FMT_YUV420P9 = @AV_PIX_FMT_NE(YUV420P9BE,YUV420P9LE)
const AV_PIX_FMT_YUV422P9 = @AV_PIX_FMT_NE(YUV422P9BE,YUV422P9LE)
const AV_PIX_FMT_YUV444P9 = @AV_PIX_FMT_NE(YUV444P9BE,YUV444P9LE)
const AV_PIX_FMT_YUV420P10 = @AV_PIX_FMT_NE(YUV420P10BE,YUV420P10LE)
const AV_PIX_FMT_YUV422P10 = @AV_PIX_FMT_NE(YUV422P10BE,YUV422P10LE)
const AV_PIX_FMT_YUV444P10 = @AV_PIX_FMT_NE(YUV444P10BE,YUV444P10LE)
const AV_PIX_FMT_YUV420P16 = @AV_PIX_FMT_NE(YUV420P16BE,YUV420P16LE)
const AV_PIX_FMT_YUV422P16 = @AV_PIX_FMT_NE(YUV422P16BE,YUV422P16LE)
const AV_PIX_FMT_YUV444P16 = @AV_PIX_FMT_NE(YUV444P16BE,YUV444P16LE)
const AV_PIX_FMT_GBRP9 = @AV_PIX_FMT_NE(GBRP9BE,GBRP9LE)
const AV_PIX_FMT_GBRP10 = @AV_PIX_FMT_NE(GBRP10BE,GBRP10LE)
const AV_PIX_FMT_GBRP16 = @AV_PIX_FMT_NE(GBRP16BE,GBRP16LE)

const PIX_FMT_RGB32 = AV_PIX_FMT_RGB32
const PIX_FMT_RGB32_1 = AV_PIX_FMT_RGB32_1
const PIX_FMT_BGR32 = AV_PIX_FMT_BGR32
const PIX_FMT_BGR32_1 = AV_PIX_FMT_BGR32_1
const PIX_FMT_GRAY16 = AV_PIX_FMT_GRAY16
const PIX_FMT_RGB48 = AV_PIX_FMT_RGB48
const PIX_FMT_RGB565 = AV_PIX_FMT_RGB565
const PIX_FMT_RGB555 = AV_PIX_FMT_RGB555
const PIX_FMT_RGB444 = AV_PIX_FMT_RGB444
const PIX_FMT_BGR48 = AV_PIX_FMT_BGR48
const PIX_FMT_BGR565 = AV_PIX_FMT_BGR565
const PIX_FMT_BGR555 = AV_PIX_FMT_BGR555
const PIX_FMT_BGR444 = AV_PIX_FMT_BGR444
const PIX_FMT_YUV420P9 = AV_PIX_FMT_YUV420P9
const PIX_FMT_YUV422P9 = AV_PIX_FMT_YUV422P9
const PIX_FMT_YUV444P9 = AV_PIX_FMT_YUV444P9
const PIX_FMT_YUV420P10 = AV_PIX_FMT_YUV420P10
const PIX_FMT_YUV422P10 = AV_PIX_FMT_YUV422P10
const PIX_FMT_YUV444P10 = AV_PIX_FMT_YUV444P10
const PIX_FMT_YUV420P16 = AV_PIX_FMT_YUV420P16
const PIX_FMT_YUV422P16 = AV_PIX_FMT_YUV422P16
const PIX_FMT_YUV444P16 = AV_PIX_FMT_YUV444P16
const PIX_FMT_GBRP9 = AV_PIX_FMT_GBRP9
const PIX_FMT_GBRP10 = AV_PIX_FMT_GBRP10
const PIX_FMT_GBRP16 = AV_PIX_FMT_GBRP16

# begin enum AVSampleFormat
const AVSampleFormat=Cint
const AV_SAMPLE_FMT_NONE = Int32(-1)
const AV_SAMPLE_FMT_U8 = Int32(0)
const AV_SAMPLE_FMT_S16 = Int32(1)
const AV_SAMPLE_FMT_S32 = Int32(2)
const AV_SAMPLE_FMT_FLT = Int32(3)
const AV_SAMPLE_FMT_DBL = Int32(4)
const AV_SAMPLE_FMT_U8P = Int32(5)
const AV_SAMPLE_FMT_S16P = Int32(6)
const AV_SAMPLE_FMT_S32P = Int32(7)
const AV_SAMPLE_FMT_FLTP = Int32(8)
const AV_SAMPLE_FMT_DBLP = Int32(9)
const AV_SAMPLE_FMT_NB = Int32(10)
# end enum AVSampleFormat

# Skipping MacroDefinition: AV_VERSION_INT ( a , b , c ) ( a << 16 | b << 8 | c )
# Skipping MacroDefinition: AV_VERSION_DOT ( a , b , c ) a ## . ## b ## . ## c
# Skipping MacroDefinition: AV_VERSION ( a , b , c ) AV_VERSION_DOT ( a , b , c )

const LIBAVUTIL_VERSION_MAJOR = 53
const LIBAVUTIL_VERSION_MINOR = 3
const LIBAVUTIL_VERSION_MICRO = 0

# Skipping MacroDefinition: LIBAVUTIL_VERSION_INT AV_VERSION_INT ( LIBAVUTIL_VERSION_MAJOR , LIBAVUTIL_VERSION_MINOR , LIBAVUTIL_VERSION_MICRO )
# Skipping MacroDefinition: LIBAVUTIL_VERSION AV_VERSION ( LIBAVUTIL_VERSION_MAJOR , LIBAVUTIL_VERSION_MINOR , LIBAVUTIL_VERSION_MICRO )

# const LIBAVUTIL_BUILD = LIBAVUTIL_VERSION_INT

# Skipping MacroDefinition: LIBAVUTIL_IDENT "Lavu" AV_STRINGIFY ( LIBAVUTIL_VERSION )
# Skipping MacroDefinition: FF_API_PIX_FMT ( LIBAVUTIL_VERSION_MAJOR < 54 )
# Skipping MacroDefinition: FF_API_CONTEXT_SIZE ( LIBAVUTIL_VERSION_MAJOR < 54 )
# Skipping MacroDefinition: FF_API_PIX_FMT_DESC ( LIBAVUTIL_VERSION_MAJOR < 54 )
# Skipping MacroDefinition: FF_API_AV_REVERSE ( LIBAVUTIL_VERSION_MAJOR < 54 )
# Skipping MacroDefinition: FF_API_AUDIOCONVERT ( LIBAVUTIL_VERSION_MAJOR < 54 )
# Skipping MacroDefinition: FF_API_CPU_FLAG_MMX2 ( LIBAVUTIL_VERSION_MAJOR < 54 )
# Skipping MacroDefinition: FF_API_LLS_PRIVATE ( LIBAVUTIL_VERSION_MAJOR < 54 )
# Skipping MacroDefinition: FF_API_AVFRAME_LAVC ( LIBAVUTIL_VERSION_MAJOR < 54 )
# Skipping MacroDefinition: FF_API_VDPAU ( LIBAVUTIL_VERSION_MAJOR < 54 )
# Skipping MacroDefinition: FF_API_XVMC ( LIBAVUTIL_VERSION_MAJOR < 54 )

const AVAudioFifo=Void

# begin enum AVDownmixType
const AVDownmixType=UInt32
const AV_DOWNMIX_TYPE_UNKNOWN = UInt32(0)
const AV_DOWNMIX_TYPE_LORO = UInt32(1)
const AV_DOWNMIX_TYPE_LTRT = UInt32(2)
const AV_DOWNMIX_TYPE_DPLII = UInt32(3)
const AV_DOWNMIX_TYPE_NB = UInt32(4)
# end enum AVDownmixType

immutable AVDownmixInfo
    preferred_downmix_type::AVDownmixType
    center_mix_level::Cdouble
    center_mix_level_ltrt::Cdouble
    surround_mix_level::Cdouble
    surround_mix_level_ltrt::Cdouble
    lfe_mix_level::Cdouble
end

# Skipping MacroDefinition: DECLARE_ALIGNED ( n , t , v ) t __attribute__ ( ( aligned ( n ) ) ) v
# Skipping MacroDefinition: DECLARE_ASM_CONST ( n , t , v ) static const t av_used __attribute__ ( ( aligned ( n ) ) ) v
# Skipping MacroDefinition: av_malloc_attrib __attribute__ ( ( __malloc__ ) )
# Skipping MacroDefinition: av_alloc_size ( ... )

const AV_OPT_FLAG_ENCODING_PARAM = 1
const AV_OPT_FLAG_DECODING_PARAM = 2
const AV_OPT_FLAG_METADATA = 4
const AV_OPT_FLAG_AUDIO_PARAM = 8
const AV_OPT_FLAG_VIDEO_PARAM = 16
const AV_OPT_FLAG_SUBTITLE_PARAM = 32
const AV_OPT_SEARCH_CHILDREN = 0x0001
const AV_OPT_SEARCH_FAKE_OBJ = 0x0002
const AV_STEREO3D_FLAG_INVERT = 1 << 0

# begin enum AVStereo3DType
const AVStereo3DType=UInt32
const AV_STEREO3D_2D = UInt32(0)
const AV_STEREO3D_SIDEBYSIDE = UInt32(1)
const AV_STEREO3D_TOPBOTTOM = UInt32(2)
const AV_STEREO3D_FRAMESEQUENCE = UInt32(3)
const AV_STEREO3D_CHECKERBOARD = UInt32(4)
const AV_STEREO3D_SIDEBYSIDE_QUINCUNX = UInt32(5)
const AV_STEREO3D_LINES = UInt32(6)
const AV_STEREO3D_COLUMNS = UInt32(7)
# end enum AVStereo3DType

immutable AVStereo3D
    _type::AVStereo3DType
    flags::Cint
end

immutable Array_16_Uint32
    d1::UInt32
    d2::UInt32
    d3::UInt32
    d4::UInt32
    d5::UInt32
    d6::UInt32
    d7::UInt32
    d8::UInt32
    d9::UInt32
    d10::UInt32
    d11::UInt32
    d12::UInt32
    d13::UInt32
    d14::UInt32
    d15::UInt32
    d16::UInt32
end

zero(::Type{Array_16_Uint32}) = Array_16_Uint32(fill(zero(UInt32),16)...)

immutable AVXTEA
    key::Array_16_Uint32
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
function truncation_projectors(H::Vector{ITensor}, ψ::Vector{ITensor}; kwargs...)
  return truncation_projectors(MPO(H), MPS(ψ); kwargs...)
end

maxdim_arg(T::Type) = round(Int, Sys.total_memory() / sizeof(T), RoundDown)
_maxdim_arg(a) = maxdim_arg(mapreduce(eltype, promote_type, a))
maxdim_arg(a::AbstractArray) = _maxdim_arg(a)
maxdim_arg(a::Union{MPS,MPO}) = _maxdim_arg(a)

function truncation_projectors(H::MPO, ψ::MPS; split_tags=("" => ""), split_plevs=(0 => 1), cutoff=1e-8, maxdim=maxdim_arg(H), tol=1e-6, maxiter=10)
  left_tags, right_tags = split_tags
  left_plev, right_plev = split_plevs
  N = length(ψ)
  @assert length(H) == N
  ψd = dag(ψ)
  Hd = dag(H)
  ψd = sim(linkinds, ψd)
  Hd = sim(linkinds, Hd)
  Hd, ψd = sim(siteinds, commoninds, Hd, ψd)
  linkindsH = linkinds(H)
  new_tags = tags.(linkindsH)
  new_plevs = plev.(linkindsH)
  L = Vector{ITensor}(undef, N)
  L[1] = ψ[1] * H[1] * Hd[1] * ψd[1]
  normalize!(L[1])
  for n in 2:N
    L[n] = L[n - 1] * ψ[n] * H[n] * Hd[n] * ψd[n]
    normalize!(L[n])
  end
  R = Vector{ITensor}(undef, N)
  R[N - 1] = ψ[N] * H[N] * Hd[N] * ψd[N]
  normalize!(R[N - 1])
  for n in reverse(2:(N - 1))
    R[n - 1] = R[n] * ψ[n] * H[n] * Hd[n] * ψd[n]
    normalize!(R[n - 1])
  end
  Us = Vector{ITensor}(undef, N - 1)
  Uds = Vector{ITensor}(undef, N - 1)
  for n in 1:(N - 1)
    link_ψn = commonind(ψ[n], ψ[n + 1])
    link_Hn = commonind(H[n], H[n + 1])
    links_n = (link_ψn, link_Hn)
    left_links_n = setprime(addtags(links_n, left_tags), left_plev)
    right_links_n = setprime(addtags(links_n, right_tags), right_plev)
    Ln = replaceinds(L[n], links_n => left_links_n)
    Rn = replaceinds(R[n], links_n => right_links_n)
    ρn = Ln * Rn
    U, S, V, spec = svd(ρn, left_links_n; maxdim=maxdim, cutoff=cutoff, lefttags=new_tags[n], leftplev=new_plevs[n])
    Ud = replaceinds(dag(U), left_links_n => right_links_n)
    overlap_old = (ρn * dag(U) * dag(Ud))[]
    err = 1.0
    iter = 0
    while err > tol && iter ≤ maxiter
      ΛU = dag(ρn * dag(Ud))
      U, P = polar(ΛU, left_links_n)
      u = commonind(U, P)
      setprime!(U, new_plevs[n]; inds=u)
      Ud = replaceinds(dag(U), left_links_n => right_links_n)
      overlap_new = (ρn * U * Ud)[]
      err = abs(overlap_new - overlap_old) / overlap_new
      iter += 1
    end
    if err > tol
      @warn "In $iter iterations, computing projector between sites $n and $(n + 1) gave a truncation error of $err, which is larger than the tolerance $tol."
    end
    Us[n] = U
    Uds[n] = Ud
  end
  return Us .=> Uds
end

# Assume prime level convention for site indices of MPO of dag(s) -> s'
function truncation_projectors_assume_primes(H::MPO, ψ::MPS; cutoff=1e-8)
  N = length(ψ)
  @assert length(H) == N
  ψd = dag(ψ)''
  # TODO: add swapprime(siteinds, H, 0 => 1)
  Hd = swapprime(dag(H), 0 => 1)'
  L = Vector{ITensor}(undef, N)
  L[1] = ψ[1] * H[1] * Hd[1] * ψd[1]
  for n in 2:N
    L[n] = L[n - 1] * ψ[n] * H[n] * Hd[n] * ψd[n]
  end
  R = Vector{ITensor}(undef, N)
  R[N - 1] = ψ[N] * H[N] * Hd[N] * ψd[N]
  for n in reverse(2:(N - 1))
    R[n - 1] = R[n] * ψ[n] * H[n] * Hd[n] * ψd[n]
  end

  for n in 1:(N - 1)
    left_tags = "left"
    right_tags = "right"
    link_ψn = commonind(ψ[n], ψ[n + 1])
    link_Hn = commonind(H[n], H[n + 1])
    links_n = (link_ψn, link_Hn)
    left_links_n = addtags(links_n, left_tags)
    right_links_n = addtags(links_n, right_tags)
    Ln = addtags(L[n], left_tags; inds=links_n)
    Rn = addtags(R[n], right_tags; inds=links_n)
    ρn = Ln * Rn
    U, S, V, spec = svd(ρn, left_links_n; cutoff=cutoff)
    Ud = replacetags(dag(U), left_tags => right_tags)
    for _ in 1:10
      ΛU = dag(ρn * Ud)
      U, _ = polar(ΛU, left_links_n)
      Ud = replacetags(dag(U), left_tags => right_tags)
    end
  end
end

ext: jl


text:
function generate_inputdata1d(u, data, polydeg, datatyp::String)
    f(x) = u(x)
    h = polydeg +1
    a = data[1]
    b = data[end]
    if datatyp == "NNPP" || datatyp == "NNPPh"
        modal = zeros(1,h)
        a1 = a[1]
        b1 = b[1]
        n,_ = gauss_lobatto_nodes_weights(h)
        _, inverse_vandermonde_legendre = vandermonde_legendre(n)
        nab = zeros(h)
        for j=1:h
            nab[j]=(b1-a1)/2*n[j]+(a1+b1)/2
        end

        indicator=zeros(1,h)
        for i in 1:h
            indicator[1,i] = u(nab[i])
        end
        multiply_dimensionwise!(modal, inverse_vandermonde_legendre, indicator)

        total_energy = 0.0
        for i in 1:h
            total_energy += modal[1, i]^2
        end
        total_energy_clip1 = 0.0
        for i in 1:(h-1)
            total_energy_clip1 += modal[1, i]^2
        end
        total_energy_clip2 = 0.0
        for i in 1:(h-2)
            total_energy_clip2 += modal[1, i]^2
        end

        x1 = (total_energy - total_energy_clip1)/total_energy
        x2 = (total_energy_clip1 - total_energy_clip2)/total_energy_clip1

        if datatyp == "NNPP"
            return x1, x2
        elseif datatyp == "NNPPh"
            return x1, x2, h
        end
    elseif datatyp == "NNRH"
        coef = zeros(h)
        help = zeros(h)
        n,w = gauss_lobatto_nodes_weights(h)
        _, inverse_vandermonde_legendre = vandermonde_legendre(n)

        #Transformation: [-1,1] to [a,b]
        nab=zeros(h)
        wab=zeros(h)
        for j=1:h
            nab[j]=(b-a)/2*n[j]+(a+b)/2
            wab[j]=(b-a)/2*w[j]
        end

        for i=1:h
            for l=1:h
                poly,_= legendre_polynomial_and_derivative(i-1,n[l])./sqrt(i-0.5)
                coef[i] += f(nab[l])*poly*wab[l]
                help[i] += poly*wab[l]*poly
            end
            coef[i] = coef[i]/help[i]
        end

        indicator = zeros(1,h)
        for i in 1:h
            indicator[1,i] = u(nab[i])
        end

        function uh(x)
            uh_approx=0
            for i = 1:h
                t = x*2/(b-a)-(a+b)/(b-a)
                phi,_= legendre_polynomial_and_derivative(i-1,t)./sqrt(i-0.5)
                uh_approx += indicator[i]*phi
            end
            uh_approx
        end
        meanuh = mean(uh,nab)

        return uh(a), uh(b), meanuh
    end

end

function generate_inputXi1d(element_id, input_data, c2e, mesh, datatyp)
    if datatyp == "NNPP" || datatyp == "NNPPh"

        return input_data[:, element_id]

    elseif datatyp == "NNRH"
        leaf_cell_ids = leaf_cells(mesh.tree)
        # get neighbors id
        cell_id = leaf_cell_ids[element_id]
        neighbor_ids = Array{Int64}(undef, 2)
        for direction in 1:n_directions(mesh.tree)

            if has_neighbor(mesh.tree, cell_id, direction)
                neighbor_cell_id = mesh.tree.neighbor_ids[direction, cell_id]
                if has_children(mesh.tree, neighbor_cell_id) # Cell has small neighbor
                    if direction == 1
                        neighbor_ids[direction] = c2e[mesh.tree.child_ids[2, neighbor_cell_id]]
                    else
                        neighbor_ids[direction] = c2e[mesh.tree.child_ids[1, neighbor_cell_id]]
                    end
                else # Cell has same refinement level neighbor
                    neighbor_ids[direction] = c2e[neighbor_cell_id]
                end
            else # Cell is small and has large neighbor
                parent_id = mesh.tree.parent_ids[cell_id]
                neighbor_cell_id = mesh.tree.neighbor_ids[direction, parent_id]
                neighbor_ids[direction] = c2e[neighbor_cell_id]
            end
        end

        # input data der neighbors und Zelle
        return input_data[3, neighbor_ids[1]], input_data[3, element_id], input_data[3, neighbor_ids[2]], input_data[1, element_id], input_data[2, element_id]
    end
end


text:
module CMSSPModel

# Stdlib Requirements
using Logging
using Printf
using Random
using Statistics
using LinearAlgebra

# Non-stdlib essentials
using StaticArrays

# POMDP Stuff
using POMDPs
using POMDPModelTools
using POMDPModels
using POMDPPolicies
using POMDPSimulators
using LocalApproximationValueIteration

# Requirements
# POMDPs.isterminal
# POMDPs.generate_sr # For continuous regions
# POMDPs.reward # For cost - decomposed

export
    TPDistribution,
    is_inf_hor,
    CMSSPState,
    CMSSPAction,
    CMSSP,
    modetype,
    continuoustype,
    modeactiontype,
    controlactiontype,
    get_modeswitch_actions,
    get_control_actions,
    set_modeswitch_mdp!,
    mode_index,
    mode_actionindex,
    BridgeSample,
    OpenLoopVertex

include("cmssp_base.jl")

end # End of submodule


text:
using ZMQ
using LibSerialPort

# Message accessors for ZMQ and SerialPort
get_message(socket::ZMQ.Socket) = String(unsafe_string(ZMQ.recv(socket)))
get_message(sp::SerialPort) = readuntil(sp, '\n', 100)

"""
Return a ZMQ subscriber socket for the provided address.
All messages except those beginning with filter_prefix will be ignored.
If filter_prefix is an empty string (the default), no messages are ignored.
"""
function get_data_source(context::ZMQ.Context, address::String, filter_prefix::String="")

    # Create a subscriber socket and attempt connection to publisher
    socket = Socket(context, SUB)
    try
        ZMQ.connect(socket, address)
    catch
        println("Could not open ZMQ socket at $address")
        quit()
    end

    # Subscribe to published messages, filtering on provided filter_prefix string
    ZMQ.set_subscribe(socket, filter_prefix)

    println("Connected to data publisher at $address")
    return socket
end

"""
Return a LibSerialPort.SerialPort
"""
function get_data_source(address::String, baud::Integer)
    sp = try
        open(address, baud)
    catch
        println("Could not open serial port at $address")
        quit()
    end
    return sp
end

max_stars_repo_name: andrewadare/realtime

max_issues_repo_name: andrewadare/realtime

max_forks_repo_name: andrewadare/realtime


text:
export diageye

using StatsFuns: logistic
using StatsFuns: softmax, softmax!
using LoopVectorization

import LinearAlgebra

diageye(::Type{T}, n::Int) where { T <: Real } = Matrix{T}(I, n, n)
diageye(n::Int)                                = diageye(Float64, n)

function normalize_sum(x::Array{Float64,1})
    x ./ sum(x)
end

const sigmoid = logistic

dtanh(x) = 1 - tanh(x)^2

"""
    mirrorlog(x)

Returns `log(1 - x)`.
"""
mirrorlog(x) = log(1 - x)

"""
    xtlog(x)

Returns `x * log(x)`.
"""
xtlog(x) = x * log(x)

"""
    negate_inplace!(A)

Returns `-A`, modifying and reusing `A` storage if possible.

See also: [`mul_inplace!`](@ref)
"""
function negate_inplace! end

negate_inplace!(A::AbstractArray) = -A
negate_inplace!(A::Real)          = -A
negate_inplace!(A::Array)         = map!(-, A, A)

"""
    mul_inplace!(alpha, A)

Returns `alpha * A`, modifying and reusing `A` storage if possible.

See also: [`negate_inplace!`](@ref)
"""
function mul_inplace! end

mul_inplace!(alpha, A::AbstractArray) = alpha * A
mul_inplace!(alpha, A::Real)          = alpha * A
mul_inplace!(alpha::T, A::Array{T}) where { T <: Real } = lmul!(alpha, A)

"""
    rank1update(A, x)
    rank1update(A, x, y)

Helper function for A + x * y'. Uses optimised BLAS version for AbstractFloats and fallbacks to a generic implementation in case of differentiation
"""
function rank1update end

rank1update(A::AbstractMatrix, x::AbstractVector)                    = rank1update(eltype(A), eltype(x), eltype(x), A, x, x)
rank1update(A::AbstractMatrix, x::AbstractVector, y::AbstractVector) = rank1update(eltype(A), eltype(x), eltype(y), A, x, y)

rank1update(A::Real, x::Real)          = rank1update(A, x, x)
rank1update(A::Real, x::Real, y::Real) = A + x * y

function rank1update(::Type{ T }, ::Type{ T }, ::Type{T}, A::Matrix, x::Vector, y::Vector) where { T <: LinearAlgebra.BlasFloat }
    return LinearAlgebra.BLAS.ger!(one(T), x, y, copy(A))
end

function rank1update(::Type{ T1 }, ::Type{ T2 }, ::Type{ T3 }, A::AbstractMatrix, x::AbstractVector, y::AbstractVector) where { T1 <: Real, T2 <: Real, T3 <: Real }
    T = promote_type(T1, T2, T3)
    B = Matrix{T}(undef, size(A))
    return rank1update!(B, A, x, y)
end

function rank1update!(B::AbstractMatrix, A::AbstractMatrix, x::AbstractVector, y::AbstractVector)
    sz = size(A)
    @inbounds for k2 in 1:sz[2]
        yk2 = y[k2]
        @inbounds for k1 in 1:sz[1]
            B[k1,k2] = A[k1,k2] + x[k1] * yk2
        end
    end
    return B
end

"""
    mul_trace(A, B)

Computes tr(A * B) without allocating A * B.
"""
function mul_trace end

mul_trace(A::Real, B::Real) = A * B

function mul_trace(A::AbstractMatrix, B::AbstractMatrix)
    sA, sB = size(A), size(B)
    @assert (sA === sB) && (length(sA) === 2) && (first(sA) === last(sA))
    result = zero(promote_type(eltype(A), eltype(B)))
    n = first(sA)
    @turbo for i in 1:n
        for j in 1:n
            result += A[i, j] * B[j, i]
        end
    end
    return result
end

"""
    v_a_vT(v, a)

Computes v*a*v^T with a single allocation.
"""
function v_a_vT(v::AbstractVector, a::Real)
    T      = promote_type(eltype(v), typeof(a))
    result = zeros(T, length(v), length(v))
    mul!(result, v, v', a, one(T))
    return result
end

function v_a_vT(v, a)
    result = v*v'
    result *= a
    return result
end

"""
    v_a_vT(v1, a, v2)

Computes v1*a*v2^T with a single allocation.
"""
function v_a_vT(v1::AbstractVector, a::Real, v2::AbstractVector)
    T      = promote_type(eltype(v1), typeof(a), eltype(v2))
    result = zeros(T, length(v1), length(v2))
    mul!(result, v1, v2', a, one(T))
    return result
end

function v_a_vT(v1, a, v2)
    result = v1*v2'
    result *= a
    return result
end


text:
function assemble_residuals!(duh, dΩ, dω, Y, qAPVM, ϕ, F, n)
  bᵤ(v) = ∫(-1.0*qAPVM*(v⋅⟂(F,n)))dΩ + ∫(DIV(v)*ϕ)dω
  bₕ(q) = ∫(-q*DIV(F))dω
  bₕᵤ((v,q)) = bᵤ(v) + bₕ(q)
  Gridap.FESpaces.assemble_vector!(bₕᵤ, duh, Y)
end

function shallow_water_rosenbrock_time_step!(
  y₂, ϕ, F, q₁, q₂, duh₁, duh₂, H1h, H1hchol, y_wrk,  # in/out args
  model, dΩ, dω, Y, V, Q, R, S, f, g, y₁, y₀,         # in args
  RTMMchol, L2MMchol, Amat, Bchol, Blfchol,           # more in args
  dt, τ, leap_frog,
  assem=SparseMatrixAssembler(SparseMatrixCSC{Float64,Int},Vector{Float64},R,S))                                   # ...yet more in args
  # energetically balanced second order rosenbrock shallow water solver
  # reference: eqns (24) and (39) of
  # https://github.com/BOM-Monash-Collaborations/articles/blob/main/energetically_balanced_time_integration/EnergeticallyBalancedTimeIntegration_SW.tex

  n = get_normal_vector(Triangulation(model))
  dt₁ = dt
  if leap_frog
    dt₁ = 2.0*dt
  end

  y₀v, y₁v, y₂v = get_free_dof_values(y₀,y₁,y₂)

  # multifield terms
  u₁, h₁ = y₁

  # 1.1: the mass flux
  compute_mass_flux!(F,dΩ,V,RTMMchol,u₁*h₁)
  # 1.2: the bernoulli function
  compute_bernoulli_potential!(ϕ,dΩ,Q,L2MMchol,u₁⋅u₁,h₁,g)
  # 1.3: the potential vorticity
  compute_potential_vorticity!(q₁,H1h,H1hchol,dΩ,R,S,h₁,u₁,f,n,assem)
  # 1.4: assemble the momentum and continuity equation residuals
  assemble_residuals!(duh₁, dΩ, dω, Y, q₁ - τ*u₁⋅∇(q₁), ϕ, F, n)

  # Solve for du₁, dh₁ over a MultiFieldFESpace
  solve!(duh₁, Blfchol, duh₁)

  # update
  y₂v .=  y₀v .+ dt₁ .* duh₁

  u₂, h₂ = y₂
  # 2.1: the mass flux
  compute_mass_flux!(F,dΩ,V,RTMMchol,u₁*(2.0*h₁ + h₂)/6.0+u₂*(h₁ + 2.0*h₂)/6.0)
  # 2.2: the bernoulli function
  compute_bernoulli_potential!(ϕ,dΩ,Q,L2MMchol,(u₁⋅u₁ + u₁⋅u₂ + u₂⋅u₂)/3.0,0.5*(h₁ + h₂),g)
  # 2.3: the potential vorticity
  compute_potential_vorticity!(q₂,H1h,H1hchol,dΩ,R,S,h₂,u₂,f,n,assem)
  # 2.4: assemble the momentum and continuity equation residuals
  assemble_residuals!(duh₂, dΩ, dω, Y, 0.5*(q₁ - τ*u₁⋅∇(q₁) + q₂ - τ*u₂⋅∇(q₂)), ϕ, F, n)

  # subtract A*[du₁,dh₁] from [du₂,dh₂] vector
  mul!(y_wrk, Amat, duh₁)
  duh₂ .= duh₂ .- y_wrk

  # solve for du₂, dh₂
  solve!(duh₂, Bchol, duh₂)

  # update yⁿ⁺¹
  y₂v .= y₁v .+ dt .* duh₂
end

function compute_mean_depth!(wrk, L2MM, h)
  # compute the mean depth over the sphere, for use in the approximate Jacobian
  mul!(wrk, L2MM, get_free_dof_values(h))
  h_int = sum(wrk)
  a_int = sum(L2MM,dims=[1,2])
  h_avg = h_int/a_int
  h_avg
end

function shallow_water_rosenbrock_time_stepper(
  model, order, degree,
  h₀, u₀, f₀, g, H₀,
  λ, dt, τ, N;
  mass_matrix_solver::Gridap.Algebra.LinearSolver=Gridap.Algebra.BackslashSolver(),
  jacobian_matrix_solver::Gridap.Algebra.LinearSolver=Gridap.Algebra.BackslashSolver(),
  leap_frog=false,
  write_diagnostics=true,
  write_diagnostics_freq=1,
  dump_diagnostics_on_screen=true,
  write_solution=false,
  write_solution_freq=N/10,
  output_dir="nswe_eq_ncells_$(num_cells(model))_order_$(order)_rosenbrock")

  # Forward integration of the shallow water equations
  Ω = Triangulation(model)
  dΩ = Measure(Ω, degree)
  dω = Measure(Ω, degree, ReferenceDomain())

  # Setup the trial and test spaces
  R, S, U, V, P, Q = setup_mixed_spaces(model, order)

  Y = MultiFieldFESpace([V, Q])
  X = MultiFieldFESpace([U, P])

  # assemble the mass matrices (RTMM mass matrix not needed)
  H1MM, _, L2MM, H1MMchol, RTMMchol, L2MMchol =
    setup_and_factorize_mass_matrices(dΩ, R, S, U, V, P, Q;
                                      mass_matrix_solver=mass_matrix_solver)

  # Project the initial conditions onto the trial spaces
  b₁(q)   = ∫(q*h₀)dΩ
  b₂(v)   = ∫(v⋅u₀)dΩ
  b₃(s)   = ∫(s*f₀)*dΩ
  rhs1    = assemble_vector(b₃, S)
  f       = FEFunction(S, copy(rhs1))
  solve!(get_free_dof_values(f),H1MMchol,get_free_dof_values(f))

  # assemble the approximate MultiFieldFESpace Jacobian
  n = get_normal_vector(Ω)

  Amat((u,p),(v,q)) = ∫(-dt*λ*f*(v⋅⟂(u,n)))dΩ + ∫(dt*λ*g*(DIV(v)*p))dω - ∫(dt*λ*H₀*(q*DIV(u)))dω
  Mmat((u,p),(v,q)) = ∫(u⋅v)dΩ + ∫(p*q)dΩ # block mass matrix
  A = assemble_matrix(Amat, X,Y)
  M = assemble_matrix(Mmat, X,Y)
  B = M - A

  Bchol = numerical_setup(symbolic_setup(jacobian_matrix_solver,B),B)
  Mchol = numerical_setup(symbolic_setup(mass_matrix_solver,M),M)
  # leap frog matrices, using 2×dt
  lf = 1.0
  if leap_frog
    lf = 2.0
  end
  Blf = M - lf*A
  Blfchol = numerical_setup(symbolic_setup(jacobian_matrix_solver,Blf),Blf)

  # multifield initial condtions
  b₄((v,q)) = b₁(q) + b₂(v)
  rhs2  = assemble_vector(b₄, Y)
  solve!(rhs2, Mchol, rhs2)
  yn  = FEFunction(Y, rhs2)

  un, hn = yn

  hnv, fv, ynv = get_free_dof_values(hn,f,yn)

  # work arrays
  h_tmp = copy(hnv)
  w_tmp = copy(fv)

  # build the potential vorticity lhs operator once just to initialise
  bmm(a,b) = ∫(a*hn*b)dΩ
  H1h      = assemble_matrix(bmm, R, S)
  H1hchol  = numerical_setup(symbolic_setup(mass_matrix_solver,H1h),H1h)

  function run_simulation(pvd=nothing)
    diagnostics_file = joinpath(output_dir,"nswe__rosenbrock_diagnostics.csv")

    ϕ      = clone_fe_function(Q,hn)
    F      = clone_fe_function(V,un)
    wn     = clone_fe_function(S,f)
    q1     = clone_fe_function(S,f)
    q2     = clone_fe_function(S,f)

    # mulifield fe functions
    ym1     = clone_fe_function(Y,yn)
    ym2     = clone_fe_function(Y,yn)
    duh1    = copy(ynv)
    duh2    = copy(ynv)
    y_wrk   = copy(ynv)

    if (write_diagnostics)
      initialize_csv(diagnostics_file,"time", "mass", "vorticity", "kinetic", "potential", "power")
    end

    # first step, no leap frog
    istep = 1
    shallow_water_rosenbrock_time_step!(yn, ϕ, F, q1, q2, duh1, duh2, H1h, H1hchol, y_wrk,
                                        model, dΩ, dω, Y, V, Q, R, S, f, g, ym1, ym2,
                                        RTMMchol, L2MMchol, A, Bchol, Blfchol, dt, τ, false)

    if (write_diagnostics && write_diagnostics_freq>0 && mod(istep, write_diagnostics_freq) == 0)
      compute_diagnostic_vorticity!(wn, dΩ, S, H1MMchol, un, get_normal_vector(Ω))
      dump_diagnostics_shallow_water!(h_tmp, w_tmp,
                                      model, dΩ, dω, S, L2MM, H1MM,
                                      hn, un, wn, ϕ, F, g, istep, dt,
                                      diagnostics_file,
                                      dump_diagnostics_on_screen)
    end
    # time step iteration loop
    for istep in 2:N
      aux=ym2
      ym2=ym1
      ym1=yn
      yn=aux
      shallow_water_rosenbrock_time_step!(yn, ϕ, F, q1, q2, duh1, duh2, H1h, H1hchol, y_wrk,
                                          model, dΩ, dω, Y, V, Q, R, S, f, g, ym1, ym2,
                                          RTMMchol, L2MMchol, A, Bchol, Blfchol, dt, τ, leap_frog)

      # IMPORTANT NOTE: We need to extract un, hn out of yn at each iteration because
      #                 the association of yn with its object instance changes at the beginning of
      #                 each iteration
      un, hn = yn
      if (write_diagnostics && write_diagnostics_freq>0 && mod(istep, write_diagnostics_freq) == 0)
        compute_diagnostic_vorticity!(wn, dΩ, S, H1MMchol, un, get_normal_vector(Ω))
        dump_diagnostics_shallow_water!(h_tmp, w_tmp,
                                        model, dΩ, dω, S, L2MM, H1MM,
                                        hn, un, wn, ϕ, F, g, istep, dt,
                                        diagnostics_file,
                                        dump_diagnostics_on_screen)
      end
      if (write_solution && write_solution_freq>0 && mod(istep, write_solution_freq) == 0)
        compute_diagnostic_vorticity!(wn, dΩ, S, H1MMchol, un, get_normal_vector(Ω))
        pvd[dt*Float64(istep)] = new_vtk_step(Ω,joinpath(output_dir,"n=$(istep)"),["hn"=>hn,"un"=>un,"wn"=>wn])
      end
    end
    hn, un
  end
  if (write_diagnostics || write_solution)
    rm(output_dir,force=true,recursive=true)
    mkdir(output_dir)
  end
  if (write_solution)
    pvdfile=joinpath(output_dir,"nswe_eq_ncells_$(num_cells(model))_order_$(order)_rosenbrock")
    paraview_collection(run_simulation,pvdfile)
  else
    run_simulation()
  end
end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
mutable struct MultiCoefplot
    name::Union{Null,String}
    note::Union{Null,AbstractCaption}
    xtitle::Union{Null,String}
    ytitle::Union{Null,String}
    dict::OrderedDict{Symbol,Coefplot}
    legends::OrderedDict{Symbol,Union{Legend,Null}}
    legends_options::LegendsOption
    other_components::Vector{Any}
    vertical::Bool

    function MultiCoefplot(dict::OrderedDict{Symbol,Coefplot},
                      xtitle::Union{Null,String}=Null(),ytitle::Union{Null,String}=Null(),
                      name::Union{Null,String}=Null(), note::Union{Null,AbstractCaption}=Null(),
                      legends_options::LegendsOption=LegendsOption(),other_components::Vector{Any}=Any[])
        @assert length(dict) > 1 "Can't make a MultiCoefplot out of a singleton"
        new_dict = deepcopy(dict)
        autolegends = OrderedDict()

        base_locs = Dict(label => i for (i,label) in enumerate(get_unique_syms(new_dict)))

        for (c_symbol, c, color) in zip(keys(new_dict),values(new_dict),COLOR_PALATTE)
            check_if_all_singlecoefplot_options_conform(c)
            setcolor!(c, color)
            c.vertical = false # reset all coefplots to non vertical
            for (k,v) in c.dict
                @assert haskey(base_locs,k)
                v.thiscoef_loc = base_locs[k]
            end
            representitive_scoptions = first(values(c.dict)).options
            this_legend = Legend(string(c_symbol),representitive_scoptions)
            push!(autolegends, c_symbol => this_legend)
        end
        new(name, note, xtitle, ytitle, new_dict, autolegends, legends_options, other_components, false)
    end
end

function get_unique_syms(m::MultiCoefplot)
    syms = Symbol[]
    for c in values(m.dict)
        append!(syms, keys(c.dict))
    end
    unique!(syms)
    return syms
end

function get_unique_syms(d::OrderedDict{Symbol,Coefplot})
    syms = Symbol[]
    for c in values(d)
        append!(syms, keys(c.dict))
    end
    unique!(syms)
    return syms
end

function setxtitle!(m::MultiCoefplot, x::String)
    m.xtitle = x
    return m
end

function setytitle!(m::MultiCoefplot, y::String)
    m.ytitle = y
    return m
end

function setname!(m::MultiCoefplot, name::String)
    m.name = name
    return m
end

function includenote!(m::MultiCoefplot, note::String)
    m.note = AbstractCaption(note,  default_note_options())
    return m
end

function setlegends!(m::MultiCoefplot, ps::Pair{Symbol,T} ...) where T <: Union{String, Missing}
    for p in ps
        if ismissing(p.second)
            m.legends[p.first] = Null()
        else
            c = m.dict[p.first]
            check_if_all_singlecoefplot_options_conform(c)
            representitive_scoptions = first(values(c.dict)).options
            m.legends[p.first] = Legend(p.second, representitive_scoptions)
        end
    end
end

function addcomponent!(m::MultiCoefplot, v::T) where T <: Union{VLine,VBand,HLine,HBand}
    push!(m.other_components,v)
    return m
end

function clearcomponents!(m::MultiCoefplot)
    m.other_components = Vector{Any}()
    return m
end

function transpose!(m::MultiCoefplot)
    m.vertical = ~m.vertical
    for c in values(m.dict)
        c.vertical = m.vertical
    end
    xtitle = deepcopy(m.xtitle)
    m.xtitle = deepcopy(m.ytitle)
    m.ytitle = xtitle
    return m
end

function MultiCoefplot(pair::Pair{Symbol,Coefplot}...)
    MultiCoefplot(OrderedDict(pair))
end

function to_plotable(m::MultiCoefplot, interval::Union{Missing,Real}=missing)
    options = gen_other_option_from_mcoefplot(m)
    legend_nonmissing = collect(values(filter(x -> ~isnull(x.second), m.legends)))
    legend_labels = PGFPlotsX.Options(Symbol("legend entries") => join([lgd.l for lgd in values(legend_nonmissing)],","))
    merge!(options,gen_legend_options(m.legends_options),legend_labels)
    scs = gen_scoefplots_from_mcoefplot(m, interval)
    Plotable(options, m.note, scs, m.other_components, values(legend_nonmissing))
end

function gen_other_option_from_mcoefplot(m::MultiCoefplot)
    specified_options_list = get_coefplot_options.(values(m.dict))
    merged_option_1 = merge(specified_options_list...)

    options_list = gen_other_option_from_coefplot.(values(m.dict);seperate_line=true)
    merged_option_2 = merge_merge([[:xtick, :xticklabels], [:ytick, :yticklabels]], options_list...)
    merged_option = merge(merged_option_1,merged_option_2)

    has_x_ticks = haskey(merged_option.dict,Symbol("extra x ticks"))
    has_y_ticks = haskey(merged_option.dict,Symbol("extra y ticks"))
    if has_x_ticks
        xtick = deepcopy(merged_option[:xtick])
        sorted_xtick = sort(xtick)
        merged_option[Symbol("extra x ticks")] = [(sorted_xtick[i]+sorted_xtick[i+1])/2 for i in 1:(length(sorted_xtick)-1)]
    end
    if has_y_ticks
        ytick = deepcopy(merged_option[:ytick])
        sorted_ytick = sort(ytick)
        merged_option[Symbol("extra x ticks")] = [(sorted_ytick[i]+sorted_ytick[i+1])/2 for i in 1:(length(sorted_ytick)-1)]
    end

    merged_option[:xmin] = minimum([option[:xmin] for option in options_list])
    merged_option[:ymin] = minimum([option[:ymin] for option in options_list])
    merged_option[:xmax] = maximum([option[:xmax] for option in options_list])
    merged_option[:ymax] = maximum([option[:ymax] for option in options_list])

    labels_and_titles = PGFPlotsX.Options()
    if ~isnull(m.xtitle)
        labels_and_titles[:xlabel] = m.xtitle
        labels_and_titles[Symbol("label style")] = "{font=\\footnotesize}"
    end
    if ~isnull(m.ytitle)
        labels_and_titles[:ylabel] = m.ytitle
        labels_and_titles[Symbol("label style")] = "{font=\\footnotesize}"
    end
    if ~isnull(m.name)
        labels_and_titles[:title] = m.name
        labels_and_titles[Symbol("title style")] = "{font=\\large}"
    end

    return merge!(merged_option,labels_and_titles) # give labels_and_titles better priority
end

function gen_scoefplots_from_mcoefplot(m::MultiCoefplot,interval::Union{Missing,Real}=missing) # need to rework. identify which singlecoefplot should be together. This implementation is very fragile
    v = PGFPlotsX.TikzElement[]
    l = length(m.dict)
    if ismissing(interval)
        interval = 1/(l+1)
    end
    d = - interval * (l-1)/2
    for coefplot in values(m.dict)
        shifted_c = shift(coefplot,d)
        d = d + interval
        append!(v,collect(shifted_c))
    end
    return v
end


text:
using ITensors
using Test
using Combinatorics: permutations
import Random: seed!

# Enable debug checking for these tests
ITensors.enable_debug_checks()

seed!(12345)

digits(::Type{T},x...) where {T} = T(sum([x[length(x)-k+1]*10^(k-1) for k=1:length(x)]))

@testset "Dense ITensor basic functionality" begin

@testset "ITensor constructors" begin
  i = Index(2, "i")
  j = Index(2, "j")
  k = Index(2, "k")
  l = Index(2, "l")

  @testset "Default" begin
    A = ITensor()
    @test storage(A) isa NDTensors.Dense{Float64}
  end

  @testset "Undef with index" begin
    A = ITensor(undef, i)
    @test storage(A) isa NDTensors.Dense{Float64}
  end

  @testset "Default with indices" begin
    A = ITensor(i,j)
    @test storage(A) isa NDTensors.Dense{Float64}
  end

  @testset "Index set operations" begin
    A = randomITensor(i, j)
    B = randomITensor(j, k)
    C = randomITensor(k, l)
    @test hascommoninds(A, B)
    @test hascommoninds(B, C)
    @test !hascommoninds(A, C)
  end

  @testset "Get element with end" begin
    a = Index(2)
    b = Index(3)
    A = randomITensor(a, b)
    @test A[end, end] == A[a => 2, b => 3]
    @test A[2, end] == A[a => 2, b => 3]
    @test A[1, end] == A[a => 1, b => 3]
    @test A[end, 2] == A[a => 2, b => 2]
    @test A[a => end, b => end] == A[a => 2, b => 3]
    @test A[a => 2, b => end] == A[a => 2, b => 3]
    @test A[a => 1, b => end] == A[a => 1, b => 3]
    @test A[a => end, b => 3] == A[a => 2, b => 3]
    @test A[a => end, b => 2] == A[a => 2, b => 2]
    @test A[b => end, a => end] == A[a => 2, b => 3]
    @test A[b => 2, a => end] == A[a => 2, b => 2]
    @test A[b => 1, a => end] == A[a => 2, b => 1]
    @test A[b => end, a => 2] == A[a => 2, b => 3]
    @test A[b => end, a => 1] == A[a => 1, b => 3]
  end

  @testset "Set element with end" begin
    _i = Index(2, "i")
    _j = Index(3, "j")

    A = ITensor(_i, _j)
    A[_i => end, _j => end] = 2.5
    @test A[_i => dim(_i), _j => dim(_j)] == 2.5

    A = ITensor(_i, _j)
    A[_j => end, _i => end] = 3.5
    @test A[_i => dim(_i), _j => dim(_j)] == 3.5

    A = ITensor(_i, _j)
    A[_j => end, _i => 1] = 4.5
    @test A[_i => 1, _j => dim(_j)] == 4.5
  end

  @testset "Random" begin
    A = randomITensor(i, j)

    # Test hasind, hasinds
    @test hasind(A, i)
    @test hasind(i)(A)

    @test hasinds(A, i)
    @test hasinds(A, j)
    @test hasinds(A, [i, j])
    @test hasinds([i, j])(A)
    @test hasinds(A, IndexSet(j))
    @test hasinds(A, j, i)
    @test hasinds(A, (i, j))
    @test hasinds(A, IndexSet(i, j))
    @test hasinds(j, i)(A)
    @test hasinds(i)(A)
    @test hasinds(IndexSet(j))(A)
    @test hasinds((i, j))(A)
    @test hasinds(IndexSet(i, j))(A)

    @test storage(A) isa NDTensors.Dense{Float64}

    @test ndims(A) == order(A) == 2 == length(inds(A))
    @test size(A) == dims(A) == (2,2)
    @test dim(A) == 4

    At = randomITensor(Index(2), Index(3))
    @test maxdim(At) == 3
    @test mindim(At) == 2
    @test dim(At, 1) == 2
    @test dim(At, 2) == 3

    B = randomITensor(IndexSet(i,j))
    @test storage(B) isa NDTensors.Dense{Float64}
    @test ndims(B) == order(B) == 2 == length(inds(B))
    @test size(B) == dims(B) == (2,2)

    A = randomITensor()
    @test eltype(A) == Float64
    @test ndims(A) == 0
end

@testset "ITensor iteration" begin
  A = randomITensor(i, j)
  Is = eachindex(A)
  @test length(Is) == dim(A)
  sumA = 0.0
  for I in Is
    sumA += A[I]
  end
  @test sumA ≈ sum(ITensors.data(A))
  sumA = 0.0
  for a in A
    sumA += a
  end
  @test sumA ≈ sum(A)
  @test sumA ≈ sum(A)
end

  @testset "From matrix" begin
    M = [1 2; 3 4]
    A = itensor(M,i,j)
    @test storage(A) isa NDTensors.Dense{Float64}

    @test M ≈ Matrix(A,i,j)
    @test M' ≈ Matrix(A,j,i)
    @test_throws DimensionMismatch vector(A)

    @test size(A,1) == size(M,1) == 2
    @test_throws BoundsError size(A,3)
    @test_throws BoundsError size(A,0)
    @test_throws ErrorException size(M,0)
    # setstorage changes the internal data but not indices
    N = [5 6; 7 8]
    A = itensor(M, i, j)
    B = ITensors.setstorage(A, NDTensors.Dense(vec(N)))
    @test N == Matrix(B, i, j)
    @test storage(A) isa NDTensors.Dense{Float64}
    @test storage(B) isa NDTensors.Dense{Int}

    M = [1 2 3; 4 5 6]
    @test_throws DimensionMismatch itensor(M,i,j)
  end

  @testset "To Matrix" begin
    TM = randomITensor(i,j)

    M1 = matrix(TM)
    for ni in i, nj in j
      @test M1[val(ni),val(nj)] ≈ TM[ni,nj]
    end

    M2 = Matrix(TM,j,i)
    for ni in i, nj in j
      @test M2[val(nj),val(ni)] ≈ TM[ni,nj]
    end

    T3 = randomITensor(i,j,k)
    @test_throws DimensionMismatch Matrix(T3,i,j)
  end

  @testset "To Vector" begin
    TV = randomITensor(i)

    V = vector(TV)
    for ni in i
      @test V[val(ni)] ≈ TV[ni]
    end
    V = Vector(TV)
    for ni in i
      @test V[val(ni)] ≈ TV[ni]
    end
    V = Vector(TV, i)
    for ni in i
      @test V[val(ni)] ≈ TV[ni]
    end
    V = Vector{ComplexF64}(TV)
    for ni in i
      @test V[val(ni)] ≈ complex(TV[ni])
    end

    T2 = randomITensor(i,j)
    @test_throws DimensionMismatch vector(T2)
  end

  @testset "Complex" begin
    A = ITensor(Complex,i,j)
    @test storage(A) isa NDTensors.Dense{Complex}
  end

  @testset "Random complex" begin
    A = randomITensor(ComplexF64,i,j)
    @test storage(A) isa NDTensors.Dense{ComplexF64}
  end

  @testset "From complex matrix" begin
    M = [1+2im 2; 3 4]
    A = itensor(M,i,j)
    @test storage(A) isa NDTensors.Dense{ComplexF64}
  end

end

@testset "Convert to complex" begin
  i = Index(2,"i")
  j = Index(2,"j")
  A = randomITensor(i,j)
  B = complex(A)
  for ii ∈ dim(i), jj ∈ dim(j)
    @test complex(A[i=>ii,j=>jj]) == B[i=>ii,j=>jj]
  end
end

@testset "Complex Number Operations" begin
  i = Index(3,"i")
  j = Index(4,"j")

  A = randomITensor(ComplexF64,i,j)

  rA = real(A)
  iA = imag(A)
  @test norm(rA+1im*iA - A) < 1E-8
  @test eltype(rA) == Float64
  @test eltype(iA) == Float64

  cA = conj(A)
  @test eltype(cA) == ComplexF64
  @test norm(cA) ≈ norm(A)

  B = randomITensor(Float64,i,j)

  cB = conj(B)
  @test eltype(cB) == Float64
  @test norm(cB) ≈ norm(B)
end

@testset "similar" begin
  i = Index(2,"i")
  j = Index(2,"j")
  A = randomITensor(i,j)
  B = similar(A)
  @test inds(B) == inds(A)
  Ac = similar(A, ComplexF32)
  @test storage(Ac) isa NDTensors.Dense{ComplexF32}
end

@testset "fill!" begin
  i = Index(2,"i")
  j = Index(2,"j")
  A = randomITensor(i,j)
  fill!(A, 1.0)
  @test all(ITensors.data(A) .== 1.0)
end

@testset "fill! using broadcast" begin
  i = Index(2,"i")
  j = Index(2,"j")
  A = randomITensor(i,j)
  A .= 1.0
  @test all(ITensors.data(A) .== 1.0)
end

@testset "copyto!" begin
  i = Index(2,"i")
  j = Index(2,"j")
  M = [1 2; 3 4]
  A = itensor(M,i,j)
  N = 2*M
  B = itensor(N,i,j)
  copyto!(A, B)
  @test A == B
  @test ITensors.data(A) == vec(N)
  A = itensor(M,i,j)
  B = itensor(N,j,i)
  copyto!(A, B)
  @test A == B
  @test ITensors.data(A) == vec(transpose(N))
end

@testset "Unary -" begin
  i = Index(2,"i")
  j = Index(2,"j")
  M = [1 2; 3 4]
  A = itensor(M,i,j)
  @test -A == itensor(-M, i, j)
end

@testset "dot" begin
  i = Index(2,"i")
  a = [1.0; 2.0]
  b = [3.0; 4.0]
  A = itensor(a,i)
  B = itensor(b,i)
  @test dot(A, B) == 11.0
end

@testset "mul!" begin
  i = Index(2; tags="i")
  j = Index(2; tags="j")
  k = Index(2; tags="k")

  A = randomITensor(i, j)
  B = randomITensor(j, k)
  C = randomITensor(i, k)
  mul!(C, A, B)
  @test C ≈ A*B

  A = randomITensor(i, j)
  B = randomITensor(j, k)
  C = randomITensor(k, i)
  mul!(C, A, B)
  @test C ≈ A*B

  A = randomITensor(i, j)
  B = randomITensor(k, j)
  C = randomITensor(i, k)
  mul!(C, A, B)
  @test C ≈ A*B

  A = randomITensor(i, j)
  B = randomITensor(k, j)
  C = randomITensor(k, i)
  mul!(C, A, B)
  @test C ≈ A*B

  A = randomITensor(j, i)
  B = randomITensor(j, k)
  C = randomITensor(i, k)
  mul!(C, A, B)
  @test C ≈ A*B

  A = randomITensor(j, i)
  B = randomITensor(j, k)
  C = randomITensor(k, i)
  mul!(C, A, B)
  @test C ≈ A*B

  A = randomITensor(j, i)
  B = randomITensor(k, j)
  C = randomITensor(i, k)
  mul!(C, A, B)
  @test C ≈ A*B

  A = randomITensor(j, i)
  B = randomITensor(k, j)
  C = randomITensor(k, i)
  mul!(C, A, B)
  @test C ≈ A*B

  A = randomITensor(i, j)
  B = randomITensor(k, j)
  C = randomITensor(k, i)
  α = 2
  β = 3
  R = mul!(copy(C), A, B, α, β)
  @test α*A*B+β*C ≈ R

  @testset "In-place bugs" begin

    @testset "Bug 1" begin
      l1 = Index(3, "l=1")
      l2 = Index(3, "l=2")
      s = Index(2, "s")

      A = randomITensor(s', s)
      B = randomITensor(l1, s, l2)

      C = randomITensor(l1, s', l2)

      C .= A .* B

      @test C ≈ A * B
    end

    @testset "Bug 2" begin
      is = [Index(n+1, "i$n") for n in 1:6]

      for ais in permutations((1, 2, 3)),
          bis in permutations((2, 3, 4)),
          cis in permutations((1, 4))
        A = randomITensor(ntuple(i->is[ais[i]], Val(length(ais))))
        B = randomITensor(ntuple(i->is[bis[i]], Val(length(bis))))
        C = randomITensor(ntuple(i->is[cis[i]], Val(length(cis))))

        C .= A .* B

        @test C ≈ A * B
      end

      for ais in permutations((1, 2, 3)),
          bis in permutations((2, 3, 4, 5)),
          cis in permutations((1, 4, 5))
        A = randomITensor(ntuple(i->is[ais[i]], Val(length(ais))))
        B = randomITensor(ntuple(i->is[bis[i]], Val(length(bis))))
        C = randomITensor(ntuple(i->is[cis[i]], Val(length(cis))))

        C .= A .* B

        @test C ≈ A * B
      end
    end

  end

  @testset "In-place outer bug" begin
    l1 = Index(3, "l=1")
    s = Index(2, "s")

    A = randomITensor(l1)
    B = randomITensor(s)
    C = randomITensor(s, l1)

    C .= A .* B

    @test C ≈ A * B
  end

  @testset "In-place contractions" begin
    i1 = Index(2, "i1")
    i2 = Index(2, "i2")
    i3 = Index(2, "i3")
    i4 = Index(2, "i4")
    i5 = Index(2, "i5")
    i6 = Index(2, "i6")
    j1 = Index(2, "j1")
    j2 = Index(2, "j2")
    j3 = Index(2, "j3")

    #A = randomITensor(s', s)
    #B = randomITensor(l1, s, l2)

    #C = randomITensor(l1, s', l2)

    C .= A .* B
    @test C ≈ A * B
  end

end

@testset "exponentiate" begin
  s1 = Index(2,"s1")
  s2 = Index(2,"s2")
  i1 = Index(2,"i1")
  i2 = Index(2,"i2")
  Amat = rand(2,2,2,2)
  A = itensor(Amat,i1,i2,s1,s2)

  Aexp = exp(A,(i1,i2),(s1,s2))
  Amatexp = reshape(exp(reshape(Amat,4,4)),2,2,2,2)
  Aexp_from_mat = itensor(Amatexp,i1,i2,s1,s2)
  @test Aexp ≈ Aexp_from_mat

  #test that exponentiation works when indices need to be permuted
  Aexp = exp(A,(s1,s2),(i1,i2))
  Amatexp = Matrix(exp(reshape(Amat,4,4))')
  Aexp_from_mat = itensor(reshape(Amatexp,2,2,2,2),s1,s2,i1,i2)
  @test Aexp ≈ Aexp_from_mat

  #test exponentiation when hermitian=true is used
  Amat = reshape(Amat, 4,4)
  Amat = reshape(Amat+Amat'+randn(4,4)*1e-10,2,2,2,2)
  A = itensor(Amat,i1,i2,s1,s2)
  Aexp = exp(A,(i1,i2),(s1,s2),ishermitian=true)
  Amatexp = reshape(parent(exp(LinearAlgebra.Hermitian(reshape(Amat,4,4)))),
                    2,2,2,2)
  Aexp_from_mat = itensor(Amatexp,i1,i2,s1,s2)
  @test Aexp ≈ Aexp_from_mat
  Aexp = exp(A,(i1,i2),(s1,s2); ishermitian=true)
  Amatexp = reshape(parent(exp(LinearAlgebra.Hermitian(reshape(Amat,4,4)))),
                    2,2,2,2)
  Aexp_from_mat = itensor(Amatexp,i1,i2,s1,s2)
  @test Aexp ≈ Aexp_from_mat
end

@testset "onehot (setelt)" begin
  i = Index(2,"i")

  T = onehot(i(1))
  @test T[i(1)] ≈ 1.0
  @test T[i(2)] ≈ 0.0

  T = setelt(i(2))
  @test T[i(1)] ≈ 0.0
  @test T[i(2)] ≈ 1.0

  # Test onehot taking Pair{Index,Int}
  T = onehot(i=>2)
  @test T[i(1)] ≈ 0.0
  @test T[i(2)] ≈ 1.0

  j = Index(2,"j")

  T = onehot(j=>2,i=>1)
  @test T[j=>1,i=>1] ≈ 0.0
  @test T[j=>2,i=>1] ≈ 1.0
  @test T[j=>1,i=>2] ≈ 0.0
  @test T[j=>2,i=>2] ≈ 0.0
end

@testset "add, subtract, and axpy" begin
  i = Index(2,"i")
  a = [1.0; 2.0]
  b = [3.0; 4.0]
  A = itensor(a,i)
  B = itensor(b,i)
  c = [5.0; 8.0]
  @test A + B == itensor([4.0; 6.0], i)
  @test axpy!(2.0, A, B) == itensor(c, i)
  a = [1.0; 2.0]
  b = [3.0; 4.0]
  A = itensor(a,i)
  B = itensor(b,i)
  c = [5.0; 8.0]
  @test (B .+= 2.0 .* A) == itensor(c, i)
  a = [1.0; 2.0]
  b = [3.0; 4.0]
  A = itensor(a,i)
  B = itensor(b,i)
  c = [8.0; 12.0]
  @test (A .= 2.0 .* A .+ 2.0 .* B) == itensor(c, i)
  a = [1.0; 2.0]
  A = itensor(a,i)
  B = ITensor(2.0)
  @test_throws DimensionMismatch A + B
  a = [1.0; 2.0]
  A = itensor(a,i)
  B = ITensor()
  @test_throws DimensionMismatch A + B
  a = [1.0; 2.0]
  a = [1.0; 2.0]
  A = itensor(a,i)
  B = ITensor(2.0)
  @test_throws DimensionMismatch A - B
  a = [1.0; 2.0]
  A = itensor(a,i)
  B = ITensor()
  @test_throws DimensionMismatch A - B
  a = [1.0; 2.0]
  A = itensor(a,i)
  B = ITensor(2.0)
  @test_throws DimensionMismatch B - A
  a = [1.0; 2.0]
  A = itensor(a,i)
  B = ITensor()
  @test_throws DimensionMismatch B - A
  a = [1.0; 2.0]
  b = [3.0; 4.0]
  A = itensor(a,i)
  B = itensor(b,i)
  c = [2.0; 2.0]
  @test B - A == itensor(c, i)
  @test A - B == -itensor(c, i)
end

@testset "mul! and rmul!" begin
  i = Index(2,"i")
  a = [1.0; 2.0]
  b = [2.0; 4.0]
  A = itensor(a,i)
  A2, A3 = copy(A), copy(A)
  B = itensor(b,i)
  @test mul!(A2, A, 2.0) == B == (A2 .= 0 .* A2 .+ 2 .* A)
  @test rmul!(A, 2.0) == B == ITensors.scale!(A3, 2)
  #make sure mul! works also when A2 has NaNs in it
  A = itensor([1.0; 2.0],i)
  A2 = itensor([NaN; 1.],i)
  @test mul!(A2, A, 2.0) == B

  i = Index(2,"i")
  j = Index(2,"j")
  M = [1 2; 3 4]
  A = itensor(M,i,j)
  N = 2*M
  B = itensor(N,j,i)
  @test ITensors.data(mul!(B, A, 2.0)) == 2.0*vec(transpose(M))
end

@testset "Construct from Array" begin
  i = Index(2,"index_i")
  j = Index(2,"index_j")

  M = [1. 2;
       3 4]
  T = itensor(M, i, j)
  T[i => 1, j => 1] = 3.3
  @test M[1, 1] == 3.3
  @test T[i => 1, j => 1] == 3.3
  @test storage(T) isa Dense{Float64}

  M = [1 2;
       3 4]
  T = itensor(M, i, j)
  T[i => 1, j => 1] = 3.3
  @test M[1, 1] == 1
  @test T[i => 1, j => 1] == 3.3
  @test storage(T) isa Dense{Float64}

  M = [1 2;
       3 4]
  T = itensor(Int, M, i, j)
  T[i => 1, j => 1] = 6
  @test M[1, 1] == 6
  @test T[i => 1, j => 1] == 6
  @test storage(T) isa Dense{Int}

  M = [1. 2;
       3 4]
  T = ITensor(M, i, j)
  T[i => 1, j => 1] = 3.3
  @test M[1, 1] == 1
  @test T[i => 1, j => 1] == 3.3
end

@testset "ITensor Array constructor view behavior" begin
  d = 2
  i = Index(d)

  # view
  A = randn(Float64, d, d)
  T = itensor(A, i', dag(i))
  @test storage(T) isa NDTensors.Dense{Float64}
  A[1, 1] = 2.0
  T[1, 1] == 2.0

  # view
  A = rand(Int, d, d)
  T = itensor(Int, A, i', dag(i))
  @test storage(T) isa NDTensors.Dense{Int}
  A[1, 1] = 2
  T[1, 1] == 2

  # no view
  A = rand(Int, d, d)
  T = itensor(A, i', dag(i))
  @test storage(T) isa NDTensors.Dense{Float64}
  A[1, 1] = 2
  T[1, 1] ≠ 2

  # no view
  A = randn(Float64, d, d)
  T = ITensor(A, i', dag(i))
  @test storage(T) isa NDTensors.Dense{Float64}
  A[1, 1] = 2
  T[1, 1] ≠ 2

  # no view
  A = rand(Int, d, d)
  T = ITensor(Int, A, i', dag(i))
  @test storage(T) isa NDTensors.Dense{Int}
  A[1, 1] = 2
  T[1, 1] ≠ 2

  # no view
  A = rand(Int, d, d)
  T = ITensor(A, i', dag(i))
  @test storage(T) isa NDTensors.Dense{Float64}
  A[1, 1] = 2
  T[1, 1] ≠ 2
end

@testset "Convert to Array" begin
  i = Index(2,"i")
  j = Index(3,"j")
  T = randomITensor(i,j)

  A = Array{Float64}(T,i,j)
  for I in CartesianIndices(T)
    @test A[I] == T[I]
  end

  T11 = T[1,1]
  T[1,1] = 1
  @test T[1,1] == 1
  @test T11 != 1
  @test A[1,1] == T11

  A = Matrix{Float64}(T,i,j)
  for I in CartesianIndices(T)
    @test A[I] == T[I]
  end

  A = Matrix(T,i,j)
  for I in CartesianIndices(T)
    @test A[I] == T[I]
  end

  A = Array(T,i,j)
  for I in CartesianIndices(T)
    @test A[I] == T[I]
  end

  T = randomITensor(i)
  A = Vector(T)
  for I in CartesianIndices(T)
    @test A[I] == T[I]
  end
end

@testset "Test isapprox for ITensors" begin
  m,n = rand(0:20,2)
  i = Index(m)
  j = Index(n)
  realData = rand(m,n)
  complexData = complex(realData)
  A = itensor(realData, i,j)
  B = itensor(complexData, i,j)
  @test A≈B
  @test B≈A
  A = permute(A,j,i)
  @test A≈B
  @test B≈A
end

@testset "permute, always_copy = false" begin
  i = Index(2)
  A = ITensor(i, i')
  Ap = permute(A, i, i')
  A[i => 1, i' => 1] = 1
  @test A[i => 1, i' => 1] == 1
  @test Ap[i => 1, i' => 1] == 1
end

@testset "permute, always_copy = true" begin
  i = Index(2)
  A = ITensor(i, i')
  Ap = permute(A, i, i'; always_copy = true)
  A[i => 1, i' => 1] = 1
  @test A[i => 1, i' => 1] == 1
  @test Ap[i => 1, i' => 1] == 0
end

@testset "ITensor tagging and priming" begin
  s1 = Index(2,"Site,s=1")
  s2 = Index(2,"Site,s=2")
  l = Index(3,"Link")
  ltmp = settags(l,"Temp")
  A1 = randomITensor(s1,l,l')
  A2 = randomITensor(s2,l',l'')
  @testset "firstind(::ITensor,::String)" begin
    @test s1==firstind(A1, "Site")
    @test s1==firstind(A1, "s=1")
    @test s1==firstind(A1, "s=1,Site")
    @test l==firstind(A1; tags="Link", plev=0)
    @test l'==firstind(A1; plev=1)
    @test l'==firstind(A1; tags="Link", plev=1)
    @test s2==firstind(A2, "Site")
    @test s2==firstind(A2, "s=2")
    @test s2==firstind(A2, "Site")
    @test s2==firstind(A2, plev=0)
    @test s2==firstind(A2; tags="s=2", plev=0)
    @test s2==firstind(A2; tags="Site", plev=0)
    @test s2==firstind(A2; tags="s=2,Site", plev=0)
    @test l'==firstind(A2; plev=1)
    @test l'==firstind(A2; tags="Link", plev=1)
    @test l''==firstind(A2; plev=2)
    @test l''==firstind(A2; tags="Link", plev=2)
  end
  @testset "addtags(::ITensor,::String,::String)" begin
    s1u = addtags(s1, "u")
    lu = addtags(l, "u")

    A1u = addtags(A1, "u")
    @test hasinds(A1u,s1u,lu,lu')

    A1u = addtags(A1, "u", "Link")
    @test hasinds(A1u,s1,lu,lu')

    A1u = addtags(A1, "u"; tags="Link")
    @test hasinds(A1u,s1,lu,lu')

    A1u = addtags(A1, "u"; plev=0)
    @test hasinds(A1u,s1u,lu,l')

    A1u = addtags(A1, "u"; tags="Link", plev=0)
    @test hasinds(A1u,s1,lu,l')

    A1u = addtags(A1, "u"; tags="Link", plev=1)
    @test hasinds(A1u,s1,l,lu')
  end
  @testset "removetags(::ITensor,::String,::String)" begin
    A2r = removetags(A2,"Site")
    @test hasinds(A2r,removetags(s2,"Site"),l',l'')

    A2r = removetags(A2,"Link";plev=1)
    @test hasinds(A2r,s2,removetags(l,"Link")',l'')

    A2r = replacetags(A2,"Link","Temp";plev=1)
    @test hasinds(A2r,s2,ltmp',l'')
  end
  @testset "replacetags(::ITensor,::String,::String)" begin
    s2tmp = replacetags(s2, "Site", "Temp")

    @test s2tmp == replacetags(s2, "Site" => "Temp")

    ltmp = replacetags(l, "Link", "Temp")

    A2r = replacetags(A2, "Site", "Temp")
    @test hasinds(A2r,s2tmp,l',l'')

    A2r = replacetags(A2, "Site" => "Temp")
    @test hasinds(A2r, s2tmp, l', l'')

    A2r = replacetags(A2,"Link","Temp")
    @test hasinds(A2r,s2,ltmp',ltmp'')

    A2r = replacetags(A2, "Site" => "Link", "Link" => "Site")
    @test hasinds(A2r, replacetags(s2, "Site" => "Link"),
                       replacetags(l', "Link" => "Site"),
                       replacetags(l'', "Link" => "Site"))
  end
  @testset "prime(::ITensor,::String)" begin
    A2p = prime(A2)
    @test A2p==A2'
    @test hasinds(A2p,s2',l'',l''')

    A2p = prime(A2,2)
    A2p = A2''
    @test hasinds(A2p,s2'',l''',l'''')

    A2p = prime(A2,"s=2")
    @test hasinds(A2p,s2',l',l'')
  end

  @testset "mapprime" begin
    @test hasinds(mapprime(A2,1,7),s2,l^7,l'')
    @test hasinds(mapprime(A2,0,1),s2',l',l'')
  end

  @testset "replaceprime" begin
    @test hasinds(mapprime(A2, 1 => 7), s2, l^7, l'')
    @test hasinds(mapprime(A2, 0 => 1), s2', l', l'')
    @test hasinds(mapprime(A2, 1 => 7, 0 => 1), s2', l^7, l'')
    @test hasinds(mapprime(A2, 1 => 2, 2 => 1), s2, l'', l')
    @test hasinds(mapprime(A2, 1 => 0, 0 => 1), s2', l, l'')
  end

  @testset "setprime" begin
    @test hasinds(setprime(A2,2,s2), s2'', l', l'')
    @test hasinds(setprime(A2,0,l''), s2, l', l)
  end

  @testset "swapprime" begin
    @test hasinds(swapprime(A2,1,3),l''',s2,l'')
  end
end

@testset "ITensor other index operations" begin

  s1 = Index(2, "Site,s=1")
  s2 = Index(2, "Site,s=2")
  l = Index(3, "Link")
  A1 = randomITensor(s1, l, l')
  A2 = randomITensor(s2, l', l'')

  @testset "ind(::ITensor)" begin
    @test ind(A1, 1) == s1
    @test ind(A1, 2) == l
  end

  @testset "replaceind and replaceinds" begin
    rA1 = replaceind(A1, s1, s2)
    @test hasinds(rA1, s2, l, l')
    @test hasinds(A1, s1, l, l')

    # Pair notation (like Julia's replace function)
    rA1 = replaceind(A1, s1 => s2)
    @test hasinds(rA1, s2, l, l')
    @test hasinds(A1, s1, l, l')

    replaceind!(A1, s1, s2)
    @test hasinds(A1, s2, l, l')

    rA2 = replaceinds(A2, (s2, l'), (s1, l))
    @test hasinds(rA2, s1, l, l'')
    @test hasinds(A2, s2, l', l'')

    # Pair notation (like Julia's replace function)
    rA2 = replaceinds(A2, s2 => s1, l' => l)
    @test hassameinds(rA2, (s1, l, l''))
    @test hassameinds(A2, (s2, l', l''))

    # Test ignoring indices that don't exist
    rA2 = replaceinds(A2, s1 => l, l' => l)
    @test hassameinds(rA2, (s2, l, l''))
    @test hassameinds(A2, (s2, l', l''))

    replaceinds!(A2, (s2, l'), (s1, l))
    @test hasinds(A2, s1, l, l'')
  end

  @testset "replaceinds fixed errors" begin
    l = Index(3; tags="l")
    s = Index(2; tags="s")
    l̃, s̃ = sim(l), sim(s)
    A = randomITensor(s, l)
    Ã = replaceinds(A, (l, s), (l̃, s̃))
    @test ind(A, 1) == s
    @test ind(A, 2) == l
    @test ind(Ã, 1) == s̃
    @test ind(Ã, 2) == l̃
    @test_throws ErrorException replaceinds(A, (l, s), (s̃, l̃))
  end

  @testset "swapinds and swapinds!" begin
    s = Index(2)
    t = Index(2)
    Ast = randomITensor(s, s', t, t')
    Ats = swapinds(Ast, (s, s'), (t, t'))
    @test Ast != Ats
    @test Ast == swapinds(Ats, (s, s'), (t, t'))

    swapinds!(Ats, (s, s'), (t, t'))
    @test Ast == Ats
  end

end #End "ITensor other index operations"

@testset "Converting Real and Complex Storage" begin

  @testset "Add Real and Complex" begin
    i = Index(2,"i")
    j = Index(2,"j")
    TC = randomITensor(ComplexF64,i,j)
    TR = randomITensor(Float64,i,j)

    S1 = TC+TR
    S2 = TR+TC
    @test typeof(S1.storage) == NDTensors.Dense{ComplexF64,Vector{ComplexF64}}
    @test typeof(S2.storage) == NDTensors.Dense{ComplexF64,Vector{ComplexF64}}
    for ii=1:dim(i),jj=1:dim(j)
      @test S1[i=>ii,j=>jj] ≈ TC[i=>ii,j=>jj]+TR[i=>ii,j=>jj]
      @test S2[i=>ii,j=>jj] ≈ TC[i=>ii,j=>jj]+TR[i=>ii,j=>jj]
    end
  end

end

@testset "ITensor, NDTensors.Dense{$SType} storage" for SType ∈ (Float64,
                                                                 ComplexF64)
  mi,mj,mk,ml,mα = 2,3,4,5,6,7
  i = Index(mi,"i")
  j = Index(mj,"j")
  k = Index(mk,"k")
  l = Index(ml,"l")
  α = Index(mα,"alpha")
  @testset "Set and get values with IndexVals" begin
    A = ITensor(SType,i,j,k)
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      A[k=>kk,j=>jj,i=>ii] = digits(SType,ii,jj,kk)
    end
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      @test A[j=>jj,k=>kk,i=>ii]==digits(SType,ii,jj,kk)
    end
    @test_throws DimensionMismatch A[1]
  end
  @testset "Test permute(ITensor,Index...)" begin
    A = randomITensor(SType,i,k,j)
    permA = permute(A,k,j,i)
    @test k==inds(permA)[1]
    @test j==inds(permA)[2]
    @test i==inds(permA)[3]
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      @test A[k=>kk,i=>ii,j=>jj]==permA[i=>ii,j=>jj,k=>kk]
    end
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      @test A[k=>kk,i=>ii,j=>jj]==permA[i=>ii,j=>jj,k=>kk]
    end
    # TODO: I think this was doing slicing, but what is the output
    # of slicing an ITensor?
    #@testset "getindex and setindex with vector of IndexVals" begin
    #    k_inds = [k=>kk for kk ∈ 1:dim(k)]
    #    for ii ∈ 1:dim(i), jj ∈ 1:dim(j)
    #      @test A[k_inds,i=>ii,j=>jj]==permA[i=>ii,j=>jj,k_inds...]
    #    end
    #    for ii ∈ 1:dim(i), jj ∈ 1:dim(j)
    #        A[k_inds,i=>ii,j=>jj]=collect(1:length(k_inds))
    #    end
    #    permA = permute(A,k,j,i)
    #    for ii ∈ 1:dim(i), jj ∈ 1:dim(j)
    #      @test A[k_inds,i=>ii,j=>jj]==permA[i=>ii,j=>jj,k_inds...]
    #    end
    #end
  end
  @testset "Set and get values with Ints" begin
    A = ITensor(SType,i,j,k)
    A = permute(A,k,i,j)
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      A[kk,ii,jj] = digits(SType,ii,jj,kk)
    end
    A = permute(A,i,j,k)
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      @test A[ii,jj,kk]==digits(SType,ii,jj,kk)
    end
  end
  @testset "Test scalar(ITensor)" begin
    x = SType(34)
    A = ITensor(x)
    @test x==scalar(A)
    A = ITensor(SType,i,j,k)
    @test_throws DimensionMismatch scalar(A)
  end
  @testset "Test norm(ITensor)" begin
    A = randomITensor(SType,i,j,k)
    @test norm(A)≈sqrt(scalar(dag(A)*A))
  end
  @testset "Test add ITensors" begin
    A = randomITensor(SType,i,j,k)
    B = randomITensor(SType,k,i,j)
    C = A+B
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      @test C[i=>ii,j=>jj,k=>kk]==A[j=>jj,i=>ii,k=>kk]+B[i=>ii,k=>kk,j=>jj]
    end
    @test array(permute(C,i,j,k))==array(permute(A,i,j,k))+array(permute(B,i,j,k))
  end

  @testset "Test factorizations of an ITensor" begin

    A = randomITensor(SType,i,j,k,l)

    @testset "Test SVD of an ITensor" begin
      U,S,V,spec,u,v = svd(A,(j,l))
      @test storage(S) isa NDTensors.Diag{Float64,Vector{Float64}}
      @test A≈U*S*V
      @test U*dag(prime(U,u))≈δ(SType,u,u') atol=1e-13
      @test V*dag(prime(V,v))≈δ(SType,v,v') atol=1e-13
    end

    @testset "Test SVD of an ITensor with different algorithms" begin
      U, S, V, spec, u, v = svd(A, j, l; alg = "recursive")
      @test storage(S) isa NDTensors.Diag{Float64,Vector{Float64}}
      @test A ≈ U * S * V
      @test U * dag(prime(U, u)) ≈ δ(SType, u, u') atol = 1e-13
      @test V * dag(prime(V, v)) ≈ δ(SType, v, v') atol = 1e-13

      U, S, V, spec, u, v = svd(A, j,l; alg = "divide_and_conquer")
      @test storage(S) isa NDTensors.Diag{Float64,Vector{Float64}}
      @test A ≈ U * S * V
      @test U * dag(prime(U, u)) ≈ δ(SType, u, u') atol = 1e-13
      @test V * dag(prime(V, v)) ≈ δ(SType, v, v') atol = 1e-13

      U, S, V, spec, u, v = svd(A, j,l; alg = "qr_iteration")
      @test storage(S) isa NDTensors.Diag{Float64,Vector{Float64}}
      @test A ≈ U * S * V
      @test U * dag(prime(U, u)) ≈ δ(SType, u, u') atol = 1e-13
      @test V * dag(prime(V, v)) ≈ δ(SType, v, v') atol = 1e-13

      @test_throws ErrorException svd(A, j,l; alg = "bad_alg")
    end

    #@testset "Test SVD of a DenseTensor internally" begin
    #  Lis = commoninds(A,IndexSet(j,l))
    #  Ris = uniqueinds(A,Lis)
    #  Lpos,Rpos = NDTensors.getperms(inds(A),Lis,Ris)
    #  # XXX this function isn't used anywhere in ITensors
    #  # (it is no longer needed because of the combiner)
    #  Ut,St,Vt,spec = svd(NDTensors.tensor(A), Lpos, Rpos)
    #  U = itensor(Ut)
    #  S = itensor(St)
    #  V = itensor(Vt)
    #  u = commonind(U, S)
    #  v = commonind(V, S)
    #  @test storage(S) isa NDTensors.Diag{Float64,Vector{Float64}}
    #  @test A≈U*S*V
    #  @test U*dag(prime(U,u))≈δ(SType,u,u') atol=1e-13
    #  @test V*dag(prime(V,v))≈δ(SType,v,v') atol=1e-13
    #end

    @testset "Test SVD truncation" begin
        ii = Index(4)
        jj = Index(4)
        T = randomITensor(ComplexF64,ii,jj)
        U,S,V = svd(T,ii;maxdim=2)
        u,s,v = svd(matrix(T))
        @test norm(U*S*V-T) ≈ sqrt(s[3]^2+s[4]^2)
    end

    @testset "Test QR decomposition of an ITensor" begin
      Q,R,q = qr(A,(i,l))
      q = commonind(Q,R)
      @test A ≈ Q*R atol=1e-13
      @test Q*dag(prime(Q,q)) ≈ δ(SType,q,q') atol=1e-13
    end

    @testset "Regression test for QR decomposition of an ITensor with all indices on one side" begin
      a = Index(2, "a")
      b = Index(2, "b")
      Vab = randomITensor(a, b)
      Q, R = qr(Vab, (a, b))
      @test hasinds(Q, (a, b))
      @test Vab ≈ Q * R atol = 1e-13
    end

    @testset "Test polar decomposition of an ITensor" begin
      U,P,u = polar(A,(k,l))
      @test A ≈ U*P atol=1e-13
      #Note: this is only satisfied when left dimensions
      #are greater than right dimensions
      UUᵀ =  U*dag(prime(U,u))

      # TODO: use a combiner to combine the u indices to make
      # this test simpler
      for ii ∈ 1:dim(u[1]),
          jj ∈ 1:dim(u[2]),
          iip ∈ 1:dim(u[1]),
          jjp ∈ 1:dim(u[2])
        val = UUᵀ[u[1](ii),u[2](jj),u[1]'(iip),u[2]'(jjp)]
        if ii==iip && jj==jjp
          @test val ≈ one(SType) atol=1e-13
        else
          @test val ≈ zero(SType) atol=1e-13
        end
      end
    end

    @testset "Test Hermitian eigendecomposition of an ITensor" begin
      is = IndexSet(i,j)
      T = randomITensor(SType, is..., prime(is)...)
      T = T + swapprime(dag(T), 0, 1)
      D, U, spec, l, r = eigen(T; ishermitian=true)
      @test T ≈ prime(U) * D * dag(U) atol=1e-13
      UUᴴ =  U * prime(dag(U), r)
      @test UUᴴ ≈ δ(r, r')
    end

    @testset "Test factorize of an ITensor" begin

      @testset "factorize default" begin
        L,R = factorize(A, (j,l))
        l = commonind(L, R)
        @test A ≈ L*R
        @test L*dag(prime(L, l)) ≈ δ(SType, l, l')
        @test R*dag(prime(R, l)) ≉ δ(SType, l, l')
      end

      @testset "factorize ortho left" begin
        L,R = factorize(A, (j,l); ortho="left")
        l = commonind(L, R)
        @test A ≈ L*R
        @test L*dag(prime(L, l)) ≈ δ(SType, l, l')
        @test R*dag(prime(R, l)) ≉ δ(SType, l, l')
      end

      @testset "factorize ortho right" begin
        L,R = factorize(A, (j,l); ortho="right")
        l = commonind(L, R)
        @test A ≈ L*R
        @test L*dag(prime(L, l)) ≉ δ(SType, l, l')
        @test R*dag(prime(R, l)) ≈ δ(SType, l, l')
      end

      @testset "factorize ortho none" begin
        L,R = factorize(A, (j,l); ortho="none")
        l = commonind(L, R)
        @test A ≈ L*R
        @test L*dag(prime(L, l)) ≉ δ(SType, l, l')
        @test R*dag(prime(R, l)) ≉ δ(SType, l, l')
      end

      @testset "factorize when ITensor has primed indices" begin
        A = randomITensor(i, i')
        L, R = factorize(A, i)
        l = commonind(L, R)
        @test A ≈ L * R
        @test L * dag(prime(L, l)) ≈ δ(SType, l, l')
        @test R * dag(prime(R, l)) ≉ δ(SType, l, l')

        @test_throws ErrorException factorize(A, i; which_decomp="svd", svd_alg = "bad_alg")
      end

    end # End factorize tests

    @testset "Test error for bad decomposition inputs" begin
      @test_throws ErrorException svd(A)
      @test_throws ErrorException svd(A, inds(A))
      @test_throws ErrorException eigen(A, inds(A), inds(A))
      #@test_throws ErrorException factorize(A)
    end

  end # End ITensor factorization testset

end # End Dense storage test

@testset "dag copy behavior" begin
  i = Index(4,"i")

  v1 = randomITensor(i)
  cv1 = dag(v1)
  cv1[1] = -1
  @test v1[1] ≈ cv1[1]

  v2 = randomITensor(i)
  cv2 = dag(v2;always_copy=true)
  orig_elt = v2[1]
  cv2[1] = -1
  @test v2[1] ≈ orig_elt

  v3 = randomITensor(ComplexF64,i)
  orig_elt = v3[1]
  cv3 = dag(v3)
  cv3[1] = -1
  @test v3[1] ≈ orig_elt

  v4 = randomITensor(ComplexF64,i)
  cv4 = dag(v4;always_copy=true)
  orig_elt = v4[1]
  cv4[1] = -1
  @test v4[1] ≈ orig_elt
end

@testset "filter ITensor indices" begin
  i = Index(2, "i")
  A = randomITensor(i, i')
  @test hassameinds(filterinds(A; plev = 0), (i,))
  @test hassameinds(inds(A; plev = 0), (i,))
  is = inds(A)
  @test hassameinds(filterinds(is; plev = 0), (i,))
  @test hassameinds(inds(is; plev = 0), (i,))
end

@testset "product" begin
  s1 = Index(2, "s1")
  s2 = Index(2, "s2")
  s3 = Index(2, "s3")

  rA = Index(3, "rA")
  lA = Index(3, "lA")

  rB = Index(3, "rB")
  lB = Index(3, "lB")

  # operator * operator
  A = randomITensor(s1', s2', dag(s1), dag(s2), lA, rA)
  B = randomITensor(s1', s2', dag(s1), dag(s2), lB, rB)
  AB = product(A, B)
  @test hassameinds(AB, (s1', s2', s1, s2, lA, rA, lB, rB))
  @test AB ≈ mapprime(prime(A; inds = (s1', s2', s1, s2)) * B, 2 => 1)

  # operator * operator, common dangling indices
  A = randomITensor(s1', s2', dag(s1), dag(s2), lA, rA)
  B = randomITensor(s1', s2', dag(s1), dag(s2), dag(lA), dag(rA))
  AB = product(A, B)
  @test hassameinds(AB, (s1', s2', s1, s2))
  @test AB ≈ mapprime(prime(A; inds = (s1', s2', s1, s2)) * B, 2 => 1)

  # operator * operator, apply_dag, common dangling indices
  A = randomITensor(s1', s2', dag(s1), dag(s2), lA, rA)
  B = randomITensor(s1', s2', dag(s1), dag(s2), lB, rB)
  ABAdag = product(A, B; apply_dag = true)
  AB = mapprime(prime(A; inds = (s1', s2', s1, s2)) * B, 2 => 1)
  Adag = swapprime(dag(A), 0 => 1; inds = (s1', s2', s1, s2))
  @test hassameinds(ABAdag, (s1', s2', s1, s2, lB, rB))
  @test ABAdag ≈ mapprime(prime(AB; inds = (s1', s2', s1, s2)) * Adag, 2 => 1)

  # operator * operator, more complicated
  A = randomITensor(s1', s2', dag(s1), dag(s2), lA, rA)
  B = randomITensor(s1', s2', s3', dag(s1), dag(s2), dag(s3), lB, rB, dag(rA))
  AB = product(A, B)
  @test hassameinds(AB, (s1', s2', s3', s1, s2, s3, lA, lB, rB))
  @test AB ≈ mapprime(prime(A; inds = (s1', s2', s1, s2)) * B, 2 => 1)

  # state * operator (1)
  A = randomITensor(dag(s1), dag(s2), lA, rA)
  B = randomITensor(s1', s2', dag(s1), dag(s2), lB, rB)
  AB = product(A, B)
  @test hassameinds(AB, (s1, s2, lA, rA, lB, rB))
  @test AB ≈ mapprime(prime(A; inds = (s1, s2)) * B)

  # state * operator (2)
  A = randomITensor(dag(s1'), dag(s2'), lA, rA)
  B = randomITensor(s1', s2', dag(s1), dag(s2), lB, rB)
  @test_throws ErrorException product(A, B)

  # operator * state (1)
  A = randomITensor(s1', s2', dag(s1), dag(s2), lA, rA)
  B = randomITensor(s1', s2', lB, rB)
  @test_throws ErrorException product(A, B)

  # operator * state (2)
  A = randomITensor(s1', s2', dag(s1), dag(s2), lA, rA)
  B = randomITensor(s1, s2, lB, rB, dag(lA))
  AB = product(A, B)
  @test hassameinds(AB, (s1, s2, rA, lB, rB))
  @test AB ≈ noprime(A * B)

  # state * state (1)
  A = randomITensor(dag(s1), dag(s2), lA, rA)
  B = randomITensor(s1, s2, lB, rB)
  AB = product(A, B)
  @test hassameinds(AB, (lA, rA, lB, rB))
  @test AB ≈ A * B

  # state * state (2)
  A = randomITensor(dag(s1'), dag(s2'), lA, rA)
  B = randomITensor(s1, s2, lB, dag(rA))
  AB = product(A, B)
  @test hassameinds(AB, (s1', s2', s1, s2, lA, lB))
  @test AB ≈ A * B

  # state * state (3)
  A = randomITensor(dag(s1'), dag(s2'), lA, rA)
  B = randomITensor(s1, s2, lB, rB)
  @test_throws ErrorException product(A, B)

  # state * state (4)
  A = randomITensor(dag(s1), dag(s2), lA, rA)
  B = randomITensor(s1', s2', lB, rB)
  @test_throws ErrorException product(A, B)

  # state * state (5)
  A = randomITensor(dag(s1'), dag(s2'), lA, rA)
  B = randomITensor(s1', s2', lB, rB)
  @test_throws ErrorException product(A, B)

end

@testset "hastags" begin
  i = Index(2, "i, x")
  j = Index(2, "j, x")
  A = randomITensor(i, j)
  @test hastags(A, "i")
  @test anyhastags(A, "i")
  @test !allhastags(A, "i")
  @test allhastags(A, "x")
end

end # End Dense ITensor basic functionality

# Disable debug checking once tests are completed
ITensors.disable_debug_checks()

nothing