text:
using Test
using LabelNumerals
using RomanNumerals

@testset "RomanNumeral test" begin
    # Constructor tests
    @test LabelNumeral(RomanNumeral, 46) == LabelNumeral(RomanNumeral, "XLVI")
    @test LabelNumeral(RomanNumeral, 1999) == LabelNumeral(RomanNumeral, "MCMXCIX")
    @test LabelNumeral(rn"XXXX") == LabelNumeral(rn"XL")
    @test LabelNumeral(RomanNumeral,1) == LabelNumeral(rn"I")
    @test LabelNumeral(rn"xx") == LabelNumeral(rn"XX")

    @test_throws Meta.ParseError LabelNumeral(RomanNumeral, "nope")
    @test_throws Meta.ParseError LabelNumeral(RomanNumeral, "XLX")

    # arithmetic tests
    @test LabelNumeral(rn"I") + LabelNumeral(rn"IX") == LabelNumeral(rn"X")
    @test LabelNumeral(rn"L") - LabelNumeral(RomanNumeral, 1) == LabelNumeral(rn"XLIX")
    @test LabelNumeral(rn"X") > LabelNumeral(rn"IX")
    @test LabelNumeral(rn"IX") <= LabelNumeral(rn"X")
    @test LabelNumeral(rn"IX") < LabelNumeral(rn"X")
    @test isless(LabelNumeral(rn"IX"), LabelNumeral(rn"X"))
    @test begin
        max(LabelNumeral(rn"C"), LabelNumeral(rn"X"), LabelNumeral(rn"M")) ==
            LabelNumeral(rn"M")
    end
    @test begin
        min(LabelNumeral(rn"M"), LabelNumeral(rn"X"), LabelNumeral(rn"C")) ==
            LabelNumeral(rn"X")
    end
    @test string(LabelNumeral(rn"X"; prefix="A-",caselower=true)) == "A-x"
end

@testset "Integer test" begin
    # Constructor tests
    @test LabelNumeral(Int, 46) == LabelNumeral(Int, "46")
    @test LabelNumeral(Int, 1999) == LabelNumeral(Int, "1999")
    @test LabelNumeral(1000) == LabelNumeral("1000")
    @test LabelNumeral(Int,1) == LabelNumeral("1")

    # arithmetic tests
    @test LabelNumeral(1) + LabelNumeral(9) == LabelNumeral(10)
    @test LabelNumeral(10) - LabelNumeral(1) == LabelNumeral(9)
    @test LabelNumeral(10) > LabelNumeral(9)
    @test LabelNumeral(9) <= LabelNumeral(10)
    @test LabelNumeral(9) < LabelNumeral(10)
    @test isless(LabelNumeral(9), LabelNumeral(10))
    @test begin
        max(LabelNumeral(100), LabelNumeral(10), LabelNumeral(1000)) ==
            LabelNumeral(1000)
    end
    @test begin
        min(LabelNumeral(1000), LabelNumeral(10), LabelNumeral(100)) ==
            LabelNumeral(10)
    end
    @test string(LabelNumeral(10; prefix="A-",caselower=true)) == "A-10"
end

@testset "AlphaNumeral test" begin
    # Constructor tests
    @test LabelNumeral(AlphaNumeral, 46) == LabelNumeral(AlphaNumeral, "TT")
    @test LabelNumeral(AlphaNumeral, 156) == LabelNumeral(AlphaNumeral, "ZZZZZZ")
    @test LabelNumeral(AlphaNumeral, 1) == LabelNumeral(AlphaNumeral, "A")
    @test hash(LabelNumeral(AlphaNumeral, "BBB")) ==
        hash(LabelNumeral(AlphaNumeral, "CCC") - LabelNumeral(AlphaNumeral,"A"))

    println(LabelNumeral(AlphaNumeral, "BBB"))

    @test convert(Bool, LabelNumeral(AlphaNumeral,100)) == true
    @test convert(BigInt, LabelNumeral(AlphaNumeral,100)) == BigInt(100)
    @test convert(Int, LabelNumeral(AlphaNumeral,100)) == 100
    @test convert(LabelNumeral{AlphaNumeral}, 100) == LabelNumeral(AlphaNumeral, 100)
    @test LabelNumeral(AlphaNumeral, 100) + 10 == 110

    @test_throws DomainError LabelNumeral(AlphaNumeral, "ABC")
    @test_throws DomainError LabelNumeral(AlphaNumeral, "157")
    @test_throws DomainError LabelNumeral(AlphaNumeral, "0")
    @test_throws DomainError LabelNumeral(AlphaNumeral, 157)
    @test_throws DomainError LabelNumeral(AlphaNumeral, 0)

    # arithmetic tests
    @test LabelNumeral(an"II") + LabelNumeral(an"A") == LabelNumeral(an"JJ")
    @test LabelNumeral(an"JJ") - LabelNumeral(an"A") == LabelNumeral(an"II")
    @test LabelNumeral(an"JJ") > LabelNumeral(an"II")
    @test LabelNumeral(an"Z") <= LabelNumeral(an"BB")
    @test LabelNumeral(an"Z") < LabelNumeral(an"BB")
    @test isless(LabelNumeral(an"Z"), LabelNumeral(an"BB"))
    @test begin
        max(LabelNumeral(an"I"), LabelNumeral(an"AA"), LabelNumeral(an"A")) ==
            LabelNumeral(an"AA")
    end
    @test begin
        min(LabelNumeral(an"I"), LabelNumeral(an"AA"), LabelNumeral(an"A")) ==
            LabelNumeral(an"A")
    end
    @test string(LabelNumeral(an"J"; prefix="A-",caselower=true)) == "A-j"
end

A2N = Dict(
    "One" => 1,
    "Two" => 2,
    "Three" => 3,
    "Four" => 4,
    "Five" => 5,
    "Six" => 6,
    "Seven" => 7,
    "Eight" => 8,
    "Nine" => 9,
    "Ten" => 10,
    "Eleven"=> 11,
    #"Twelve" => 12,
    "Thirteen" => 13,
    "Fourteen" => 14,
    "Fifteen" => 15,
    "Sixteen" => 16,
    "Seventeen" => 17,
    "Eighteen" => 18,
    "Nineteen" => 19,
    "Twenty" => 20
)

@testset "LookupNumeral test" begin
    registerLookupNumerals(A2N, 1, 20)
    # Constructor tests
    @test LabelNumeral(LookupNumeral, 10) == LabelNumeral(LookupNumeral, "Ten")
    @test LabelNumeral(LookupNumeral, 11) == LabelNumeral(LookupNumeral, "Eleven")
    @test LabelNumeral(LookupNumeral, 1) == LabelNumeral(LookupNumeral, "One")
    @test hash(LabelNumeral(LookupNumeral, "Ten")) ==
        hash(LabelNumeral(LookupNumeral, "Eleven") - LabelNumeral(LookupNumeral,"One"))

    println(LabelNumeral(LookupNumeral, "Ten"))

    @test convert(Bool, LabelNumeral(LookupNumeral,10)) == true
    @test convert(BigInt, LabelNumeral(LookupNumeral,11)) == BigInt(11)
    @test convert(Int, LabelNumeral(LookupNumeral,10)) == 10
    @test convert(LabelNumeral{LookupNumeral}, 10) == LabelNumeral(LookupNumeral, 10)
    @test LabelNumeral(LookupNumeral, 9) + 10 == 19

    @test_throws DomainError LabelNumeral(LookupNumeral, "ABC")
    @test_throws DomainError LabelNumeral(LookupNumeral, "157")
    @test_throws DomainError LabelNumeral(LookupNumeral, "0")
    @test_throws DomainError LabelNumeral(LookupNumeral, 157)
    @test_throws DomainError LabelNumeral(LookupNumeral, 0)
    @test_throws DomainError begin
        LabelNumeral(ln"Ten") + LabelNumeral(ln"Two")
    end

    # arithmetic tests
    @test LabelNumeral(ln"One") + LabelNumeral(ln"Two") == LabelNumeral(ln"Three")
    @test LabelNumeral(ln"Three") - LabelNumeral(ln"Two") == LabelNumeral(ln"One")
    @test LabelNumeral(ln"Four") > LabelNumeral(ln"One")
    @test LabelNumeral(ln"Ten") <= LabelNumeral(ln"Eleven")
    @test LabelNumeral(ln"Ten") < LabelNumeral(ln"Eleven")
    @test isless(LabelNumeral(ln"Ten"), LabelNumeral(ln"Eleven"))
    @test begin
        max(LabelNumeral(ln"One"), LabelNumeral(ln"Five"), LabelNumeral(ln"Ten")) ==
            LabelNumeral(ln"Ten")
    end
    @test begin
        min(LabelNumeral(ln"Ten"), LabelNumeral(ln"Five"), LabelNumeral(ln"One")) ==
            LabelNumeral(ln"One")
    end
    @test string(LabelNumeral(ln"Ten"; prefix="A-",caselower=true)) == "A-ten"
end

@testset "AlphaNumNumeral test" begin
    # Constructor tests
    @test LabelNumeral(AlphaNumNumeral, 46) == LabelNumeral(AlphaNumNumeral, "BU")
    @test LabelNumeral(AlphaNumNumeral, 156) == LabelNumeral(AlphaNumNumeral, "GA")
    @test LabelNumeral(AlphaNumNumeral, 1) == LabelNumeral(AlphaNumNumeral, "B")
    @test hash(LabelNumeral(AlphaNumNumeral, "BBB")) ==
        hash(LabelNumeral(AlphaNumNumeral, "BBB") - LabelNumeral(AlphaNumNumeral,"A"))

    println(LabelNumeral(AlphaNumNumeral, "BBB"))

    @test convert(Bool, LabelNumeral(AlphaNumNumeral,100)) == true
    @test convert(BigInt, LabelNumeral(AlphaNumNumeral,100)) == BigInt(100)
    @test convert(Int, LabelNumeral(AlphaNumNumeral,100)) == 100
    @test convert(LabelNumeral{AlphaNumNumeral}, 100) == LabelNumeral(AlphaNumNumeral, 100)
    @test LabelNumeral(AlphaNumNumeral, 100) + 10 == 110

    @test_throws DomainError LabelNumeral(AlphaNumNumeral, "A23BC")
    @test_throws DomainError LabelNumeral(AlphaNumNumeral, "157")
    @test_throws DomainError LabelNumeral(AlphaNumNumeral, "0")

    # arithmetic tests
    @test LabelNumeral(ann"II") + LabelNumeral(ann"A") == LabelNumeral(ann"II")
    @test LabelNumeral(ann"JJ") - LabelNumeral(ann"B") == LabelNumeral(ann"JI")
    @test LabelNumeral(ann"JJ") > LabelNumeral(ann"II")
    @test LabelNumeral(ann"Z") <= LabelNumeral(ann"BB")
    @test LabelNumeral(ann"Z") < LabelNumeral(ann"BB")
    @test isless(LabelNumeral(ann"Z"), LabelNumeral(ann"BB"))
    @test begin
        max(LabelNumeral(ann"I"), LabelNumeral(ann"Z"), LabelNumeral(ann"A")) ==
            LabelNumeral(ann"Z")
    end
    @test begin
        min(LabelNumeral(ann"I"), LabelNumeral(ann"Z"), LabelNumeral(ann"A")) ==
            LabelNumeral(ann"A")
    end
    @test string(LabelNumeral(ann"J"; prefix="A-",caselower=true)) == "A-j"
end

@testset "Search Label" begin
    @test length(findLabels("XXX"; pfxList=["X",""])) == 6
    @test length(findLabels("X-XX"; pfxList=["X-",""])) == 3
    @test length(findLabels("A10"; pfxList=["A",""])) == 1
end


text:
using LinearAlgebra
using Zygote: hessian_dual, hessian_reverse

@testset "hessian: $hess" for hess in [hessian_dual, hessian_reverse]

  @test hess(x -> x[1]*x[2], randn(2)) ≈ [0 1; 1 0]
  @test hess(((x,y),) -> x*y, randn(2)) ≈ [0 1; 1 0]  # original docstring version
  @test hess(x -> sum(x.^3), [1 2; 3 4]) ≈ Diagonal([6, 18, 12, 24])
  @test hess(sin, pi/2) ≈ -1

  @test_throws Exception hess(sin, im*pi)
  @test_throws Exception hess(x -> x+im, pi)
  @test_throws Exception hess(identity, randn(2))
end

@testset "jacobian(f, args...)" begin
  @test jacobian(identity, [1,2])[1] == [1 0; 0 1]

  j1 = jacobian((a,x) -> a.^2 .* x, [1,2,3], 1)
  @test j1[1] ≈ Diagonal([2,4,6])
  @test j1[2] ≈ [1, 4, 9]
  @test j1[2] isa Vector

  j2 = jacobian((a,t) -> sum(a .* t[1]) + t[2], [1,2,3], (4,5))  # scalar output is OK
  @test j2[1] == [4 4 4]
  @test j2[1] isa Matrix
  @test j2[2] === nothing  # input other than Number, Array is ignored

  j3 = jacobian((a,d) -> prod(a, dims=d), [1 2; 3 4], 1)
  @test j3[1] ≈ [3 1 0 0; 0 0 4 2]
  @test j3[2] ≈ [0, 0]  # pullback is always Nothing, but array already allocated

  j4 = jacobian([1,2,-3,4,-5]) do xs
    map(x -> x>0 ? x^3 : 0, xs)  # pullback gives Nothing for some elements x
  end
  @test j4[1] ≈ Diagonal([3,12,0,48,0])

  j5 = jacobian((x,y) -> hcat(x[1], y), fill(pi), exp(1))  # zero-array
  @test j5[1] isa Matrix
  @test vec(j5[1]) == [1, 0]

  @test_throws ArgumentError jacobian(identity, [1,2,3+im])
  @test_throws ArgumentError jacobian(sum, [1,2,3+im])  # scalar, complex

  f6(x,y) = abs2.(x .* y)
  g6 = gradient(first∘f6, [1+im, 2], 3+4im)
  j6 = jacobian((x,y) -> abs2.(x .* y), [1+im, 2], 3+4im)
  @test j6[1][1,:] ≈ g6[1]
  @test j6[2][1] ≈ g6[2]
end

@testset "jacobian(loss, ::Params)" begin
  xs = [1 2; 3 4]
  ys = [5,7,9];
  Jxy = jacobian(() -> ys[1:2] .+ sum(xs.^2), Params([xs, ys]))
  @test Jxy[ys] ≈ [1 0 0; 0 1 0]
  @test Jxy[xs] ≈ [2 6 4 8; 2 6 4 8]
end

ext: jl


text:
################################################################################################
#
# OpenCV_hpp.jl
#
# Julia constants from OpenCV headers:
#
# cvdef.h
# core.h
# core.hpp
# types_c.h
# imgproc.hpp
# imgcodec.hpp
# videoio.hpp
# highgui.hpp
# features2d.hpp
# objectdetect.hpp
# calib3d.hpp (not in use)
#
################################################################################################

# C++ constant from <cfloat>
const DBL_MAX = 1.79769e+308
# constant declarations (not all)

# core.hpp
# Operations on arrays
const SORT_EVERY_ROW    = 0
const SORT_EVERY_COLUMN = 1
const SORT_ASCENDING    = 0
const SORT_DESCENDING   = 16

const COVAR_SCRAMBLED = 0
const COVAR_NORMAL    = 1
const COVAR_USE_AVG   = 2
const COVAR_SCALE     = 4
const COVAR_ROWS      = 8
const COVAR_COLS      = 16

# matrix decomposition types
const DECOMP_LU       = 0
const DECOMP_SVD      = 1
const DECOMP_EIG      = 2
const DECOMP_CHOLESKY = 3
const DECOMP_QR       = 4
const DECOMP_NORMAL   = 16

# norm types
const NORM_INF       = 1
const NORM_L1        = 2
const NORM_L2        = 4
const NORM_L2SQR     = 5
const NORM_HAMMING   = 6
const NORM_HAMMING2  = 7
const NORM_TYPE_MASK = 7
const NORM_RELATIVE  = 8
const NORM_MINMAX    = 32

# comparison types
const CMP_EQ = 0
const CMP_GT = 1
const CMP_GE = 2
const CMP_LT = 3
const CMP_LE = 4
const CMP_NE = 5

const GEMM_1_T = 1
const GEMM_2_T = 2
const GEMM_3_T = 4

const DFT_INVERSE        = 1
const DFT_SCALE          = 2
const DFT_ROWS           = 4
const DFT_COMPLEX_OUTPUT = 16
const DFT_REAL_OUTPUT    = 32
const DCT_INVERSE        = DFT_INVERSE
const DCT_ROWS           = DFT_ROWS

# Various border types, image boundaries are denoted with '|'
const BORDER_CONSTANT    = 0 # iiiiii|abcdefgh|iiiiiii  with some specified 'i'
const BORDER_REPLICATE   = 1 # aaaaaa|abcdefgh|hhhhhhh
const BORDER_REFLECT     = 2 # fedcba|abcdefgh|hgfedcb
const BORDER_WRAP        = 3 # cdefgh|abcdefgh|abcdefg
const BORDER_REFLECT_101 = 4 # gfedcb|abcdefgh|gfedcba
const BORDER_TRANSPARENT = 5 # uvwxyz|absdefgh|ijklmno
const BORDER_REFLECT101  = BORDER_REFLECT_101
const BORDER_DEFAULT     = BORDER_REFLECT_101
const BORDER_ISOLATED    = 16 # do not look outside of ROI

# Mat rtypes (core.hpp)
const CV_8U  = 0       # unsigned int8    (0:255)
const CV_8S  = 1       # signed int8      (-128:127)
const CV_16U = 2       # 16bit unsigned
const CV_16S = 3       # 16bit signed
const CV_32S = 4       # 32bit signed int
const CV_32F = 5       # 32bit float
const CV_64F = 6       # 64bit double

const CV_MatType = Int
const CV_8UC1 = 0
const CV_8UC2 = 8
const CV_8UC3 = 16
const CV_8UC4 = 24
const CV_8SC1 = 1
const CV_8SC2 = 9
const CV_8SC3 = 17
const CV_8SC4 = 25
const CV_16UC1 = 2
const CV_16UC2 = 10
const CV_16UC3 = 18
const CV_16UC4 = 26
const CV_16SC1 = 3
const CV_16SC2 = 11
const CV_16SC3 = 19
const CV_16SC4 = 27
const CV_32SC1 = 4
const CV_32SC2 = 12
const CV_32SC3 = 20
const CV_32SC4 = 28
const CV_32FC1 = 5
const CV_32FC2 = 13
const CV_32FC3 = 21
const CV_32FC4 = 29
const CV_64FC1 = 6
const CV_64FC2 = 14
const CV_64FC3 = 22
const CV_64FC4 = 30
#end of CV_MatType

# Make a Dict() for lookup of image formats (CV_Mat types)
CV_MAT_TYPE= Dict( 0 => "CV_8UC1",
                   8 => "CV_8UC2",
                  16 => "CV_8UC3",
                  24 => "CV_8UC4",
                   1 => "CV_8SC1",
                   9 => "CV_8SC2",
                  17 => "CV_8SC3",
                  25 => "CV_8SC4",
                   2 => "CV_16UC1",
                  10 => "CV_16UC2",
                  18 => "CV_16UC3",
                  26 => "CV_16UC4",
                   3 => "CV_16SC1",
                  11 => "CV_16SC2",
                  19 => "CV_16SC3",
                  27 => "CV_16SC4",
                   4 => "CV_32SC1",
                  12 => "CV_32SC2",
                  20 => "CV_32SC3",
                  28 => "CV_32SC4",
                   5 => "CV_32FC1",
                  13 => "CV_32FC2",
                  21 => "CV_32FC3",
                  29 => "CV_32FC4",
                   6 => "CV_64FC1",
                  14 => "CV_64FC2",
                  22 => "CV_64FC3",
                  30 => "CV_64FC4")

# K-means
const KMEANS_RANDOM_CENTERS     = 0
const KMEANS_PP_CENTERS         = 2
const KMEANS_USE_INITIAL_LABELS = 1

# Line thickness
const FILLED  = -1
const LINE_4  = 4
const LINE_8  = 8
const LINE_AA = 16

# Font type and style
const FONT_HERSHEY_SIMPLEX        = 0
const FONT_HERSHEY_PLAIN          = 1
const FONT_HERSHEY_DUPLEX         = 2
const FONT_HERSHEY_COMPLEX        = 3
const FONT_HERSHEY_TRIPLEX        = 4
const FONT_HERSHEY_COMPLEX_SMALL  = 5
const FONT_HERSHEY_SCRIPT_SIMPLEX = 6
const FONT_HERSHEY_SCRIPT_COMPLEX = 7
const FONT_ITALIC                 = 16

# cvdef.h
const CV_CPU_NONE   = 0
const CV_CPU_MMX    = 1
const CV_CPU_SSE    = 2
const CV_CPU_SSE2   = 3
const CV_CPU_SSE3   = 4
const CV_CPU_SSSE3  = 5
const CV_CPU_SSE4_1 = 6
const CV_CPU_SSE4_2 = 7
const CV_CPU_POPCNT = 8
const CV_CPU_AVX    = 10
const CV_CPU_NEON   = 11

# core_c.h
const CV_AUTOSTEP        = 0x7fffffff
const CV_MAX_ARR         = 10
const CV_NO_DEPTH_CHECK  = 1
const CV_NO_CN_CHECK     = 2
const CV_NO_SIZE_CHECK   = 4
const CV_CMP_EQ   =  0
const CV_CMP_GT   =  1
const CV_CMP_GE   =  2
const CV_CMP_LT   =  3
const CV_CMP_LE   =  4
const CV_CMP_NE   =  5
const CV_CHECK_RANGE  = 1
const CV_CHECK_QUIET  = 2
const CV_RAND_UNI     = 0
const CV_RAND_NORMAL  = 1

const CV_REDUCE_SUM = 0
const CV_REDUCE_AVG = 1
const CV_REDUCE_MAX = 2
const CV_REDUCE_MIN = 3

const CV_SORT_EVERY_ROW    = 0
const CV_SORT_EVERY_COLUMN = 1
const CV_SORT_ASCENDING    = 0
const CV_SORT_DESCENDING   = 16

const CV_GEMM_A_T         = 1
const CV_GEMM_B_T         = 2
const CV_GEMM_C_T         = 4
const CV_SVD_MODIFY_A     = 1
const CV_SVD_U_T          = 2
const CV_SVD_V_T          = 4
const CV_LU               = 0
const CV_SVD              = 1
const CV_SVD_SYM          = 2
const CV_CHOLESKY         = 3
const CV_QR               = 4
const CV_NORMAL           = 16
const CV_COVAR_SCRAMBLED  = 0
const CV_COVAR_NORMAL     = 1
const CV_COVAR_USE_AVG    = 2
const CV_COVAR_SCALE      = 4
const CV_COVAR_ROWS       = 8
const CV_COVAR_COLS       = 16

const CV_PCA_DATA_AS_ROW  = 0
const CV_PCA_DATA_AS_COL  = 1
const CV_PCA_USE_AVG      = 2
const CV_C                = 1
const CV_L1               = 2
const CV_L2               = 4
const CV_NORM_MASK        = 7
const CV_RELATIVE         = 8
const CV_DIFF             = 16
const CV_MINMAX           = 32

const CV_DIFF_C        = CV_DIFF
const CV_DIFF_L1       = CV_DIFF
const CV_DIFF_L2       = CV_DIFF
const CV_RELATIVE_C    = CV_RELATIVE
const CV_RELATIVE_L1   = CV_RELATIVE
const CV_RELATIVE_L2   = CV_RELATIVE

const CV_REDUCE_SUM    = 0
const CV_REDUCE_AVG    = 1
const CV_REDUCE_MAX    = 2
const CV_REDUCE_MIN    = 3

const CV_DXT_FORWARD  = 0
const CV_DXT_INVERSE  = 1
const CV_DXT_SCALE    = 2
const CV_DXT_INV_SCALE = (CV_DXT_INVERSE + CV_DXT_SCALE)
const CV_DXT_INVERSE_SCALE = CV_DXT_INV_SCALE
const CV_DXT_ROWS     = 4
const CV_DXT_MUL_CONJ = 8

const CV_FRONT = 1
const CV_BACK = 0

const CV_GRAPH_VERTEX      =   1
const CV_GRAPH_TREE_EDGE   =   2
const CV_GRAPH_BACK_EDGE   =   4
const CV_GRAPH_FORWARD_EDGE =  8
const CV_GRAPH_CROSS_EDGE   =  16
const CV_GRAPH_ANY_EDGE     =  30
const CV_GRAPH_NEW_TREE     =  32
const CV_GRAPH_BACKTRACKING =  64
const CV_GRAPH_OVER         =  -1
const CV_GRAPH_ALL_ITEMS    = -1

# flags for graph vertices and edges */
const CV_GRAPH_ITEM_VISITED_FLAG = (1 << 30)
# const CV_IS_GRAPH_VERTEX_VISITED(vtx) = (((CvGraphVtx*)(vtx))->flags & CV_GRAPH_ITEM_VISITED_FLAG)
# const CV_IS_GRAPH_EDGE_VISITED(edge) \ (((CvGraphEdge*)(edge))->flags & CV_GRAPH_ITEM_VISITED_FLAG)
const CV_GRAPH_SEARCH_TREE_NODE_FLAG =  (1 << 29)
const CV_GRAPH_FORWARD_EDGE_FLAG     =  (1 << 28)

# ocl.hpp
const TYPE_DEFAULT     = (1 << 0)
const TYPE_CPU         = (1 << 1)
const TYPE_GPU         = (1 << 2)
const TYPE_ACCELERATOR = (1 << 3)
const TYPE_DGPU        = TYPE_GPU + (1 << 16)
const TYPE_IGPU        = TYPE_GPU + (1 << 17)
const TYPE_ALL         = 0xFFFFFFFF

# types_c.h
const CV_RNG_COEFF = 4164903690       #U = unsigned
const CV_PI = 3.1415926535897932384626433832795
const CV_LOG2 = 0.69314718055994530941723212145818
const CV_AUTO_STEP = 0x7fffffff
const CV_CN_MAX = 512
const CV_CN_SHIFT = 3
const CV_HIST_ARRAY = 0
const CV_HIST_SPARSE =  1
const CV_HIST_TREE  = CV_HIST_SPARSE
const CV_HIST_UNIFORM =  1
const CV_MAGIC_MASK =  0xFFFF0000
const CV_MAT_MAGIC_VAL =  0x42420000
const CV_MATND_MAGIC_VAL = 0x42430000
const CV_MAX_DIM = 32
const CV_NODE_EMPTY = 32
const CV_NODE_FLOW = 8
const CV_NODE_INT = 1
const CV_NODE_INTEGER = CV_NODE_INT
const CV_NODE_MAP =  6
const CV_NODE_NAMED = 64
const CV_NODE_NONE =  0
const CV_NODE_REAL =  2
const CV_NODE_REF =  4
const CV_NODE_SEQ = 5
const CV_NODE_SEQ_SIMPLE =   256
const CV_NODE_STR =   3
const CV_NODE_TYPE_MASK  =  7
const CV_NODE_USER =   16
const CV_SEQ_ELTYPE_BITS =  12
const CV_SEQ_ELTYPE_CODE = CV_8UC1
const CV_SEQ_ELTYPE_CONNECTED_COMP = 0
const CV_SEQ_ELTYPE_GENERIC =  0
const CV_SEQ_ELTYPE_GRAPH_EDGE = 0
const CV_SEQ_ELTYPE_GRAPH_VERTEX = 0
const CV_SEQ_ELTYPE_INDEX = CV_32SC1
const CV_SEQ_ELTYPE_POINT = CV_32SC2
const CV_SEQ_ELTYPE_POINT3D = CV_32FC3
const CV_SEQ_ELTYPE_TRIAN_ATR = 0
const CV_SEQ_KIND_BITS = 2
const CV_SEQ_MAGIC_VAL =  0x42990000
const CV_SET_MAGIC_VAL = 0x42980000
const CV_SPARSE_MAT_MAGIC_VAL  = 0x42440000
const CV_STORAGE_APPEND =  2
const CV_STORAGE_MAGIC_VAL  = 0x42890000
const CV_STORAGE_READ =  0
const CV_STORAGE_WRITE  = 1
const CV_STORAGE_WRITE_BINARY =  CV_STORAGE_WRITE
const CV_STORAGE_WRITE_TEXT  = CV_STORAGE_WRITE
const CV_SUBMAT_FLAG_SHIFT  = 15
const CV_TERMCRIT_EPS  = 2
const CV_TERMCRIT_ITER  = 1
const CV_USRTYPE1 =  7
const CV_WHOLE_SEQ_END_INDEX =  0x3fffffff

# Errors
const CV_StsOk = 0
const CV_StsBackTrace = -1
const CV_StsError = -2
const CV_StsInternal = -3
const CV_StsNoMem = -4
const CV_StsBadArg = -5
const CV_StsBadFunc = -6
const CV_StsNoConv = -7
const CV_StsAutoTrace = -8
const CV_HeaderIsNull = -9
const CV_BadImageSize = -10
const CV_BadOffset = -11
const CV_BadDataPtr = -12
const CV_BadStep = -13
const CV_BadModelOrChSeq = -14
const CV_BadNumChannels = -15
const CV_BadNumChannel1U = -16
const CV_BadDepth = -17
const CV_BadAlphaChannel = -18
const CV_BadOrder = -19
const CV_BadOrigin = -20
const CV_BadAlign = -21
const CV_BadCallBack = -22
const CV_BadTileSize = -23
const CV_BadCOI = -24
const CV_BadROISize = -25
const CV_MaskIsTiled = -26
const CV_StsNullPtr = -27,
const CV_StsVecLengthErr = -28
const CV_StsFilterStructContentErr = -29
const CV_StsKernelStructContentErr = -30
const CV_StsFilterOffsetErr = -31
const CV_StsBadSize = -201
const CV_StsDivByZero = -202
const CV_StsInplaceNotSupported = -203
const CV_StsObjectNotFound = -204
const CV_StsUnmatchedFormats = -205
const CV_StsBadFlag = -206
const CV_StsBadPoint = -207
const CV_StsBadMask = -208
const CV_StsUnmatchedSizes = -209
const CV_StsUnsupportedFormat = -210
const CV_StsOutOfRange = -211
const CV_StsParseError = -212
const CV_StsNotImplemented = -213
const CV_StsBadMemBlock = -214
const CV_StsAssert = -215
const CV_GpuNotSupported = -216
const CV_GpuApiCallError = -217
const CV_GpuNppCallError = -218
const CV_GpuCufftCallError = -219

# imgproc.hpp

# morphological operations
const MORPH_ERODE    = 0
const MORPH_DILATE   = 1
const MORPH_OPEN     = 2
const MORPH_CLOSE    = 3
const MORPH_GRADIENT = 4
const MORPH_TOPHAT   = 5
const MORPH_BLACKHAT = 6

# shape of the structuring element
const MORPH_RECT    = 0
const MORPH_CROSS   = 1
const MORPH_ELLIPSE = 2

# interpolation algorithm
const INTER_NEAREST        = 0 # nearest neighbor interpolation
const INTER_LINEAR         = 1 # bilinear interpolation
const INTER_CUBIC          = 2 # bicubic interpolation
const INTER_AREA           = 3 # area-based (or super) interpolation
const INTER_LANCZOS4       = 4 #  Lanczos interpolation over 8x8 neighborhood

const INTER_MAX            = 7 # mask for interpolation codes
const WARP_FILL_OUTLIERS   = 8
const WARP_INVERSE_MAP     = 16

const INTER_BITS      = 5
const INTER_BITS2     = INTER_BITS * 2
const INTER_TAB_SIZE  = 1 << INTER_BITS
const INTER_TAB_SIZE2 = INTER_TAB_SIZE * INTER_TAB_SIZE

# Distance types for Distance Transform and M-estimators
const DIST_USER    = -1  # User defined distance
const DIST_L1      = 1   # distance = |x1-x2| + |y1-y2|
const DIST_L2      = 2   # the simple euclidean distance
const DIST_C       = 3   # distance = max(|x1-x2|,|y1-y2|)
const DIST_L12     = 4   # L1-L2 metric: distance = 2(sqrt(1+x*x/2) - 1))
const DIST_FAIR    = 5   # distance = c^2(|x|/c-log(1+|x|/c)), c = 1.3998
const DIST_WELSCH  = 6   # distance = c^2/2(1-exp(-(x/c)^2)), c = 2.9846
const DIST_HUBER   = 7   # distance = |x|<c ? x^2/2 : c(|x|-c/2), c=1.345

# Mask size for distance transform
const DIST_MASK_3       = 3
const DIST_MASK_5       = 5
const DIST_MASK_PRECISE = 0

# type of the threshold operation
const THRESH_BINARY     = 0 # value = value > threshold ? max_value : 0
const THRESH_BINARY_INV = 1 # value = value > threshold ? 0 : max_value
const THRESH_TRUNC      = 2 # value = value > threshold ? threshold : value
const THRESH_TOZERO     = 3 #  value = value > threshold ? value : 0
const THRESH_TOZERO_INV = 4 #  value = value > threshold ? 0 : value
const THRESH_MASK       = 7
const THRESH_OTSU       = 8  #  use Otsu algorithm to choose the optimal threshold value

# adaptive threshold algorithm
const ADAPTIVE_THRESH_MEAN_C     = 0
const ADAPTIVE_THRESH_GAUSSIAN_C = 1

const PROJ_SPHERICAL_ORTHO  = 0
const PROJ_SPHERICAL_EQRECT = 1

# class of the pixel in GrabCut algorithm
const GC_BGD    = 0 #  background
const GC_FGD    = 1 # foreground
const GC_PR_BGD = 2 # most probably background
const GC_PR_FGD = 3  # most probably foreground

# GrabCut algorithm flags
const GC_INIT_WITH_RECT  = 0
const GC_INIT_WITH_MASK  = 1
const GC_EVAL            = 2

# distanceTransform algorithm flags
const DIST_LABEL_CCOMP = 0
const DIST_LABEL_PIXEL = 1

# floodfill algorithm flags
const FLOODFILL_FIXED_RANGE = 1 << 16
const FLOODFILL_MASK_ONLY   = 1 << 17

# type of the template matching operation
const TM_SQDIFF        = 0
const TM_SQDIFF_NORMED = 1
const TM_CCORR         = 2
const TM_CCORR_NORMED  = 3
const TM_CCOEFF        = 4
const TM_CCOEFF_NORMED = 5

# connected components algorithm output formats
const CC_STAT_LEFT   = 0
const CC_STAT_TOP    = 1
const CC_STAT_WIDTH  = 2
const CC_STAT_HEIGHT = 3
const CC_STAT_AREA   = 4
const CC_STAT_MAX    = 5

# mode of the contour retrieval algorithm
const RETR_EXTERNAL  = 0 # retrieve only the most external (top-level) contours
const RETR_LIST      = 1 # retrieve all the contours without any hierarchical information
const RETR_CCOMP     = 2 #  retrieve the connected components (that can possibly be nested)
const RETR_TREE      = 3 #  retrieve all the contours and the whole hierarchy
const RETR_FLOODFILL = 4

# the contour approximation algorithm
const CHAIN_APPROX_NONE      = 1
const CHAIN_APPROX_SIMPLE    = 2
const CHAIN_APPROX_TC89_L1   = 3
const CHAIN_APPROX_TC89_KCOS = 4

# Variants of a Hough transform
const HOUGH_STANDARD      = 0
const HOUGH_PROBABILISTIC = 1
const HOUGH_MULTI_SCALE   = 2
const HOUGH_GRADIENT      = 3

# Variants of Line Segment Detector
const LSD_REFINE_NONE = 0
const LSD_REFINE_STD  = 1
const LSD_REFINE_ADV  = 2

# Histogram comparison methods
const HISTCMP_CORREL        = 0
const HISTCMP_CHISQR        = 1
const HISTCMP_INTERSECT     = 2
const HISTCMP_BHATTACHARYYA = 3
const HISTCMP_HELLINGER     = HISTCMP_BHATTACHARYYA
const HISTCMP_CHISQR_ALT    = 4
const HISTCMP_KL_DIV        = 5

# the color conversion code
const COLOR_BGR2BGRA     = 0
const COLOR_RGB2RGBA     = COLOR_BGR2BGRA

const COLOR_BGRA2BGR     = 1
const COLOR_RGBA2RGB     = COLOR_BGRA2BGR

const COLOR_BGR2RGBA     = 2
const COLOR_RGB2BGRA     = COLOR_BGR2RGBA

const COLOR_RGBA2BGR     = 3
const COLOR_BGRA2RGB     = COLOR_RGBA2BGR

const COLOR_BGR2RGB      = 4
const COLOR_RGB2BGR      = COLOR_BGR2RGB

const COLOR_BGRA2RGBA    = 5
const COLOR_RGBA2BGRA    = COLOR_BGRA2RGBA

const COLOR_BGR2GRAY     = 6
const COLOR_RGB2GRAY     = 7
const COLOR_GRAY2BGR     = 8
const COLOR_GRAY2RGB     = COLOR_GRAY2BGR
const COLOR_GRAY2BGRA    = 9
const COLOR_GRAY2RGBA    = COLOR_GRAY2BGRA
const COLOR_BGRA2GRAY    = 10
const COLOR_RGBA2GRAY    = 11
const COLOR_BGR2BGR565   = 12
const COLOR_RGB2BGR565   = 13
const COLOR_BGR5652BGR   = 14
const COLOR_BGR5652RGB   = 15
const COLOR_BGRA2BGR565  = 16
const COLOR_RGBA2BGR565  = 17
const COLOR_BGR5652BGRA  = 18
const COLOR_BGR5652RGBA  = 19

const COLOR_GRAY2BGR565  = 20
const COLOR_BGR5652GRAY  = 21

const COLOR_BGR2BGR555   = 22
const COLOR_RGB2BGR555   = 23
const COLOR_BGR5552BGR   = 24
const COLOR_BGR5552RGB   = 25
const COLOR_BGRA2BGR555  = 26
const COLOR_RGBA2BGR555  = 27
const COLOR_BGR5552BGRA  = 28
const COLOR_BGR5552RGBA  = 29

const COLOR_GRAY2BGR555  = 30
const COLOR_BGR5552GRAY  = 31

const COLOR_BGR2XYZ      = 32
const COLOR_RGB2XYZ      = 33
const COLOR_XYZ2BGR      = 34
const COLOR_XYZ2RGB      = 35

const COLOR_BGR2YCrCb    = 36
const COLOR_RGB2YCrCb    = 37
const COLOR_YCrCb2BGR    = 38
const COLOR_YCrCb2RGB    = 39

const COLOR_BGR2HSV      = 40
const COLOR_RGB2HSV      = 41

const COLOR_BGR2Lab      = 44
const COLOR_RGB2Lab      = 45

const COLOR_BGR2Luv      = 50
const COLOR_RGB2Luv      = 51
const COLOR_BGR2HLS      = 52
const COLOR_RGB2HLS      = 53

const COLOR_HSV2BGR      = 54,
const COLOR_HSV2RGB      = 55,

const COLOR_Lab2BGR      = 56
const COLOR_Lab2RGB      = 57
const COLOR_Luv2BGR      = 58
const COLOR_Luv2RGB      = 59
const COLOR_HLS2BGR      = 60
const COLOR_HLS2RGB      = 61

const COLOR_BGR2HSV_FULL = 66
const COLOR_RGB2HSV_FULL = 67
const COLOR_BGR2HLS_FULL = 68
const COLOR_RGB2HLS_FULL = 69

const COLOR_HSV2BGR_FULL = 70
const COLOR_HSV2RGB_FULL = 71
const COLOR_HLS2BGR_FULL = 72
const COLOR_HLS2RGB_FULL = 73

const COLOR_LBGR2Lab     = 74
const COLOR_LRGB2Lab     = 75
const COLOR_LBGR2Luv     = 76
const COLOR_LRGB2Luv     = 77

const COLOR_Lab2LBGR     = 78
const COLOR_Lab2LRGB     = 79
const COLOR_Luv2LBGR     = 80
const COLOR_Luv2LRGB     = 81

const COLOR_BGR2YUV      = 82
const COLOR_RGB2YUV      = 83
const COLOR_YUV2BGR      = 84
const COLOR_YUV2RGB      = 85

# YUV 4:2:0 family to RGB
const COLOR_YUV2RGB_NV12  = 90
const COLOR_YUV2BGR_NV12  = 91
const COLOR_YUV2RGB_NV21  = 92
const COLOR_YUV2BGR_NV21  = 93
const COLOR_YUV420sp2RGB  = COLOR_YUV2RGB_NV21
const COLOR_YUV420sp2BGR  = COLOR_YUV2BGR_NV21

const COLOR_YUV2RGBA_NV12 = 94
const COLOR_YUV2BGRA_NV12 = 95
const COLOR_YUV2RGBA_NV21 = 96
const COLOR_YUV2BGRA_NV21 = 97
const COLOR_YUV420sp2RGBA = COLOR_YUV2RGBA_NV21
const COLOR_YUV420sp2BGRA = COLOR_YUV2BGRA_NV21

const COLOR_YUV2RGB_YV12  = 98
const COLOR_YUV2BGR_YV12  = 99
const COLOR_YUV2RGB_IYUV  = 100
const COLOR_YUV2BGR_IYUV  = 101
const COLOR_YUV2RGB_I420  = COLOR_YUV2RGB_IYUV
const COLOR_YUV2BGR_I420  = COLOR_YUV2BGR_IYUV
const COLOR_YUV420p2RGB   = COLOR_YUV2RGB_YV12
const COLOR_YUV420p2BGR   = COLOR_YUV2BGR_YV12

const COLOR_YUV2RGBA_YV12 = 102
const COLOR_YUV2BGRA_YV12 = 103
const COLOR_YUV2RGBA_IYUV = 104
const COLOR_YUV2BGRA_IYUV = 105
const COLOR_YUV2RGBA_I420 = COLOR_YUV2RGBA_IYUV
const COLOR_YUV2BGRA_I420 = COLOR_YUV2BGRA_IYUV
const COLOR_YUV420p2RGBA  = COLOR_YUV2RGBA_YV12
const COLOR_YUV420p2BGRA  = COLOR_YUV2BGRA_YV12

const COLOR_YUV2GRAY_420  = 106
const COLOR_YUV2GRAY_NV21 = COLOR_YUV2GRAY_420
const COLOR_YUV2GRAY_NV12 = COLOR_YUV2GRAY_420
const COLOR_YUV2GRAY_YV12 = COLOR_YUV2GRAY_420
const COLOR_YUV2GRAY_IYUV = COLOR_YUV2GRAY_420
const COLOR_YUV2GRAY_I420 = COLOR_YUV2GRAY_420
const COLOR_YUV420sp2GRAY = COLOR_YUV2GRAY_420
const COLOR_YUV420p2GRAY  = COLOR_YUV2GRAY_420

# YUV 4:2:2 family to RGB
const COLOR_YUV2RGB_UYVY = 107
const COLOR_YUV2BGR_UYVY = 108
#COLOR_YUV2RGB_VYUY = 109
#COLOR_YUV2BGR_VYUY = 110
const COLOR_YUV2RGB_Y422 = COLOR_YUV2RGB_UYVY
const COLOR_YUV2BGR_Y422 = COLOR_YUV2BGR_UYVY
const COLOR_YUV2RGB_UYNV = COLOR_YUV2RGB_UYVY
const COLOR_YUV2BGR_UYNV = COLOR_YUV2BGR_UYVY

const COLOR_YUV2RGBA_UYVY = 111
const COLOR_YUV2BGRA_UYVY = 112
#COLOR_YUV2BGRA_VYUY = 114
const COLOR_YUV2RGBA_Y422 = COLOR_YUV2RGBA_UYVY
const COLOR_YUV2BGRA_Y422 = COLOR_YUV2BGRA_UYVY
const COLOR_YUV2RGBA_UYNV = COLOR_YUV2RGBA_UYVY
const COLOR_YUV2BGRA_UYNV = COLOR_YUV2BGRA_UYVY

const COLOR_YUV2RGB_YUY2 = 115
const COLOR_YUV2BGR_YUY2 = 116
const COLOR_YUV2RGB_YVYU = 117
const COLOR_YUV2BGR_YVYU = 118
const COLOR_YUV2RGB_YUYV = COLOR_YUV2RGB_YUY2
const COLOR_YUV2BGR_YUYV = COLOR_YUV2BGR_YUY2
const COLOR_YUV2RGB_YUNV = COLOR_YUV2RGB_YUY2
const COLOR_YUV2BGR_YUNV = COLOR_YUV2BGR_YUY2

const COLOR_YUV2RGBA_YUY2 = 119
const COLOR_YUV2BGRA_YUY2 = 120
const COLOR_YUV2RGBA_YVYU = 121
const COLOR_YUV2BGRA_YVYU = 122
const COLOR_YUV2RGBA_YUYV = COLOR_YUV2RGBA_YUY2
const COLOR_YUV2BGRA_YUYV = COLOR_YUV2BGRA_YUY2
const COLOR_YUV2RGBA_YUNV = COLOR_YUV2RGBA_YUY2
const COLOR_YUV2BGRA_YUNV = COLOR_YUV2BGRA_YUY2

const COLOR_YUV2GRAY_UYVY = 123
const COLOR_YUV2GRAY_YUY2 = 124
#CV_YUV2GRAY_VYUY    = CV_YUV2GRAY_UYVY
const COLOR_YUV2GRAY_Y422 = COLOR_YUV2GRAY_UYVY
const COLOR_YUV2GRAY_UYNV = COLOR_YUV2GRAY_UYVY
const COLOR_YUV2GRAY_YVYU = COLOR_YUV2GRAY_YUY2
const COLOR_YUV2GRAY_YUYV = COLOR_YUV2GRAY_YUY2
const COLOR_YUV2GRAY_YUNV = COLOR_YUV2GRAY_YUY2

# alpha premultiplication
const COLOR_RGBA2mRGBA    = 125
const COLOR_mRGBA2RGBA    = 126

# RGB to YUV 4:2:0 family
const COLOR_RGB2YUV_I420  = 127
const COLOR_BGR2YUV_I420  = 128
const COLOR_RGB2YUV_IYUV  = COLOR_RGB2YUV_I420
const COLOR_BGR2YUV_IYUV  = COLOR_BGR2YUV_I420

const COLOR_RGBA2YUV_I420 = 129
const COLOR_BGRA2YUV_I420 = 130
const COLOR_RGBA2YUV_IYUV = COLOR_RGBA2YUV_I420
const COLOR_BGRA2YUV_IYUV = COLOR_BGRA2YUV_I420
const COLOR_RGB2YUV_YV12  = 131
const COLOR_BGR2YUV_YV12  = 132
const COLOR_RGBA2YUV_YV12 = 133
const COLOR_BGRA2YUV_YV12 = 134

#Demosaicing
const COLOR_BayerBG2BGR = 46
const COLOR_BayerGB2BGR = 47
const COLOR_BayerRG2BGR = 48
const COLOR_BayerGR2BGR = 49

const COLOR_BayerBG2RGB = COLOR_BayerRG2BGR
const COLOR_BayerGB2RGB = COLOR_BayerGR2BGR
const COLOR_BayerRG2RGB = COLOR_BayerBG2BGR
const COLOR_BayerGR2RGB = COLOR_BayerGB2BGR

const COLOR_BayerBG2GRAY = 86
const COLOR_BayerGB2GRAY = 87
const COLOR_BayerRG2GRAY = 88
const COLOR_BayerGR2GRAY = 89

# Demosaicing using Variable Number of Gradients
const COLOR_BayerBG2BGR_VNG = 62
const COLOR_BayerGB2BGR_VNG = 63
const COLOR_BayerRG2BGR_VNG = 64
const COLOR_BayerGR2BGR_VNG = 65

const COLOR_BayerBG2RGB_VNG = COLOR_BayerRG2BGR_VNG
const COLOR_BayerGB2RGB_VNG = COLOR_BayerGR2BGR_VNG
const COLOR_BayerRG2RGB_VNG = COLOR_BayerBG2BGR_VNG
const COLOR_BayerGR2RGB_VNG = COLOR_BayerGB2BGR_VNG

# Edge-Aware Demosaicing
const COLOR_BayerBG2BGR_EA  = 135
const COLOR_BayerGB2BGR_EA  = 136
const COLOR_BayerRG2BGR_EA  = 137
const COLOR_BayerGR2BGR_EA  = 138

const COLOR_BayerBG2RGB_EA  = COLOR_BayerRG2BGR_EA
const COLOR_BayerGB2RGB_EA  = COLOR_BayerGR2BGR_EA
const COLOR_BayerRG2RGB_EA  = COLOR_BayerBG2BGR_EA
const COLOR_BayerGR2RGB_EA  = COLOR_BayerGB2BGR_EA

const COLOR_COLORCVT_MAX  = 139

# COLORMAP
const COLORMAP_AUTUMN  = 0
const COLORMAP_BONE    = 1
const COLORMAP_JET     = 2
const COLORMAP_WINTER  = 3
const COLORMAP_RAINBOW = 4
const COLORMAP_OCEAN   = 5
const COLORMAP_SUMMER  = 6
const COLORMAP_SPRING  = 7
const COLORMAP_COOL    = 8
const COLORMAP_HSV     = 9
const COLORMAP_PINK    = 10
const COLORMAP_HOT     = 11

# types of intersection between rectangles
const INTERSECT_NONE = 0
const INTERSECT_PARTIAL  = 1
const INTERSECT_FULL  = 2

# highgui.hpp
const WindowProperty = UInt32
const WINDOW_NORMAL     = 0x00000000 # the user can resize the window (no constraint)
const WINDOW_AUTOSIZE   = 0x00000001 # the user cannot resize the window, the size is constrainted by the image displayed
const WINDOW_OPENGL     = 0x00001000 # window with opengl support
const WINDOW_FULLSCREEN = 1          # change the window to fullscreen
const WINDOW_FREERATIO  = 0x00000100 # the image expends as much as it can (no ratio constraint)
const WINDOW_KEEPRATIO  = 0x00000000 # the ratio of the image is respected

const getWindowPropertyFlag = Int64
const WND_PROP_FULLSCREEN   = 0  # fullscreen property    (can be WINDOW_NORMAL or WINDOW_FULLSCREEN)
const WND_PROP_AUTOSIZE     = 1  # autosize property      (can be WINDOW_NORMAL or WINDOW_AUTOSIZE)
const WND_PROP_ASPECT_RATIO = 2  # window's aspect ration (can be set to WINDOW_FREERATIO or WINDOW_KEEPRATIO);
const WND_PROP_OPENGL       = 3  # opengl support

const mouseEventFlag = Int64
const EVENT_MOUSEMOVE      = 0
const EVENT_LBUTTONDOWN    = 1
const EVENT_RBUTTONDOWN    = 2
const EVENT_MBUTTONDOWN    = 3
const EVENT_LBUTTONUP      = 4
const EVENT_RBUTTONUP      = 5
const EVENT_MBUTTONUP      = 6
const EVENT_LBUTTONDBLCLK  = 7
const EVENT_RBUTTONDBLCLK  = 8
const EVENT_MBUTTONDBLCLK  = 9
const EVENT_MOUSEHWHEEL    = 11
const EVENT_FLAG_SHIFTKEY  = 16
const EVENT_FLAG_ALTKEY    = 32

# videoio.hpp
const captureFlag = Int64
const CAP_ANY          = 0     # autodetect
const CAP_VFW          = 200   # platform native
const CAP_V4L          = 200
const CAP_V4L2         = CAP_V4L
const CAP_FIREWARE     = 300   # IEEE 1394 drivers
const CAP_FIREWIRE     = CAP_FIREWARE
const CAP_IEEE1394     = CAP_FIREWARE
const CAP_DC1394       = CAP_FIREWARE
const CAP_CMU1394      = CAP_FIREWARE
const CAP_QT           = 500   # QuickTime
const CAP_UNICAP       = 600   # Unicap drivers
const CAP_DSHOW        = 700   # DirectShow (via videoInput)
const CAP_PVAPI        = 800   # PvAPI, Prosilica GigE SDK
const CAP_OPENNI       = 900   # OpenNI (for Kinect)
const CAP_OPENNI_ASUS  = 910   # OpenNI (for Asus Xtion)
const CAP_ANDROID      = 1000  # Android
const CAP_XIAPI        = 1100  # XIMEA Camera API
const CAP_AVFOUNDATION = 1200  # AVFoundation framework for iOS (OS X Lion will have the same API)
const CAP_GIGANETIX    = 1300  # Smartek Giganetix GigEVisionSDK
const CAP_MSMF         = 1400  # Microsoft Media Foundation (via videoInput)
const CAP_INTELPERC    = 1500  # Intel Perceptual Computing SDK
const CAP_OPENNI2      = 1600  # OpenNI2 (for Kinect)

# generic properties (based on DC1394 properties)
const capturePROPflag = Int64
const CAP_PROP_POS_MSEC       = 0
const CAP_PROP_POS_FRAMES     = 1
const CAP_PROP_POS_AVI_RATIO  = 2
const CAP_PROP_FRAME_WIDTH    = 3
const CAP_PROP_FRAME_HEIGHT   = 4
const CAP_PROP_FPS            = 5
const CAP_PROP_FOURCC         = 6
const CAP_PROP_FRAME_COUNT    = 7
const CAP_PROP_FORMAT         = 8
const CAP_PROP_MODE           = 9
const CAP_PROP_BRIGHTNESS    = 10
const CAP_PROP_CONTRAST      = 11
const CAP_PROP_SATURATION    = 12
const CAP_PROP_HUE           = 13
const CAP_PROP_GAIN          = 14
const CAP_PROP_EXPOSURE      = 15
const CAP_PROP_CONVERT_RGB   = 16
const CAP_PROP_WHITE_BALANCE_BLUE_U = 17
const CAP_PROP_RECTIFICATION = 18
const CAP_PROP_MONOCROME     = 19
const CAP_PROP_SHARPNESS     = 20
const CAP_PROP_AUTO_EXPOSURE = 21        # DC1394: exposure control
const CAP_PROP_GAMMA         = 22
const CAP_PROP_TEMPERATURE   = 23
const CAP_PROP_TRIGGER       = 24
const CAP_PROP_TRIGGER_DELAY = 25
const CAP_PROP_WHITE_BALANCE_RED_V = 26
const CAP_PROP_ZOOM          = 27
const CAP_PROP_FOCUS         = 28
const CAP_PROP_GUID          = 29
const CAP_PROP_ISO_SPEED     = 30
const CAP_PROP_BACKLIGHT     = 32
const CAP_PROP_PAN           = 33
const CAP_PROP_TILT          = 34
const CAP_PROP_ROLL          = 35
const CAP_PROP_IRIS          = 36
const CAP_PROP_SETTINGS      = 37

# FourCC codecs for VideoWriter
const CV_FOURCC_IYUV = [pointer("I"),pointer("Y"),pointer("U"),pointer("V")]   #for yuv420p into an uncompressed AVI
const CV_FOURCC_DIV3 = [pointer("D"),pointer("I"),pointer("V"),pointer("3")]   #for DivX MPEG-4 codec
const CV_FOURCC_MP42 = [pointer("M"),pointer("P"),pointer("4"),pointer("2")]   #for MPEG-4 codec
const CV_FOURCC_DIVX = [pointer("D"),pointer("I"),pointer("V"),pointer("X")]   #for DivX codec
const CV_FOURCC_PIM1 = [pointer("P"),pointer("I"),pointer("M"),pointer("1")]   #for MPEG-1 codec
const CV_FOURCC_I263 = [pointer("I"),pointer("2"),pointer("6"),pointer("3")]   #for ITU H.263 codec
const CV_FOURCC_MPEG = [pointer("M"),pointer("P"),pointer("E"),pointer("G")]   #for MPEG-1 codec

# imgcodecs.hpp
const imreadFlag = Int64
const IMREAD_UNCHANGED  = -1         # 8bit, color or not
const IMREAD_GRAYSCALE  = 0          # 8bit, gray
const IMREAD_COLOR      = 1          # ?, color
const IMREAD_ANYDEPTH   = 2          # any depth, ?
const IMREAD_ANYCOLOR   = 4          # ?, any color
const IMREAD_LOAD_GDAL  = 8          # Use gdal driver

const imwriteFlag = Int64
const IMWRITE_JPEG_QUALITY        = 1
const IMWRITE_JPEG_PROGRESSIVE    = 2
const IMWRITE_JPEG_OPTIMIZE       = 3
const IMWRITE_JPEG_RST_INTERVAL   = 4
const IMWRITE_JPEG_LUMA_QUALITY   = 5
const IMWRITE_JPEG_CHROMA_QUALITY = 6
const IMWRITE_PNG_COMPRESSION     = 16
const IMWRITE_PNG_STRATEGY        = 17
const IMWRITE_PNG_BILEVEL         = 18
const IMWRITE_PXM_BINARY          = 32
const IMWRITE_WEBP_QUALITY        = 64

const IMWRITE_PNG_STRATEGY_DEFAULT      = 0
const IMWRITE_PNG_STRATEGY_FILTERED     = 1
const IMWRITE_PNG_STRATEGY_HUFFMAN_ONLY = 2
const IMWRITE_PNG_STRATEGY_RLE          = 3
const IMWRITE_PNG_STRATEGY_FIXED        = 4

# objectdetect.hpp
# CvHaarClassifierCascade
const CASCADE_DO_CANNY_PRUNING    = 1
const CASCADE_SCALE_IMAGE         = 2
const CASCADE_FIND_BIGGEST_OBJECT = 4
const CASCADE_DO_ROUGH_SEARCH     = 8

# HOGDescriptor
const L2Hys = 0
const DEFAULT_NLEVELS = 64

# features2d.hpp
const DEFAULT = 0
const DRAW_OVER_OUTIMG = 1
const NOT_DRAW_SINGLE_POINTS = 2
const DRAW_RICH_KEYPOINTS = 4

# calib3d.hpp
const LMEDS  = 4 #  least-median algorithm
const RANSAC = 8 #  RANSAC algorithm

const SOLVEPNP_ITERATIVE = 0
const SOLVEPNP_EPNP      = 1  #  F.Moreno-Noguer et al, "EPnP: Efficient Perspective-n-Point Camera Pose Estimation"
const SOLVEPNP_P3P       = 2  #  X.S. Gao et al, "Complete Solution Classification for the Perspective-Three-Point Problem"
const SOLVEPNP_DLS       = 3  #  Joel A. Hesch et al, "A Direct Least-Squares (DLS) Method for PnP"
const SOLVEPNP_UPNP      = 4  #  A.Penate-Sanchez et al. "Exhaustive Linearization for Robust Camera Pose and Focal Length Estimation"

const CALIB_CB_ADAPTIVE_THRESH = 1
const CALIB_CB_NORMALIZE_IMAGE = 2
const CALIB_CB_FILTER_QUADS    = 4
const CALIB_CB_FAST_CHECK      = 8

const CALIB_CB_SYMMETRIC_GRID  = 1
const CALIB_CB_ASYMMETRIC_GRID = 2
const CALIB_CB_CLUSTERING      = 4

const CALIB_USE_INTRINSIC_GUESS = 0x00001
const CALIB_FIX_ASPECT_RATIO    = 0x00002
const CALIB_FIX_PRINCIPAL_POINT = 0x00004
const CALIB_ZERO_TANGENT_DIST   = 0x00008
const CALIB_FIX_FOCAL_LENGTH    = 0x00010
const CALIB_FIX_K1              = 0x00020
const CALIB_FIX_K2              = 0x00040
const CALIB_FIX_K3              = 0x00080
const CALIB_FIX_K4              = 0x00800
const CALIB_FIX_K5              = 0x01000
const CALIB_FIX_K6              = 0x02000
const CALIB_RATIONAL_MODEL      = 0x04000
const CALIB_THIN_PRISM_MODEL    = 0x08000
const CALIB_FIX_S1_S2_S3_S4     = 0x10000

# only for stereo
const CALIB_FIX_INTRINSIC       = 0x00100
const CALIB_SAME_FOCAL_LENGTH   = 0x00200

# for stereo rectification
const CALIB_ZERO_DISPARITY      = 0x00400

# the algorithm for finding fundamental matrix
const FM_7POINT = 1 #  7-point algorithm
const FM_8POINT = 2 #  8-point algorithm
const FM_LMEDS  = 4 # least-median algorithm
const FM_RANSAC = 8  # RANSAC algorithm

# video. Object Tracking
const OPTFLOW_USE_INITIAL_FLOW     = 4
const OPTFLOW_LK_GET_MIN_EIGENVALS = 8
const OPTFLOW_FARNEBACK_GAUSSIAN   = 256

const MOTION_TRANSLATION = 0
const MOTION_EUCLIDEAN   = 1
const MOTION_AFFINE      = 2
const MOTION_HOMOGRAPHY  = 3

# OpenGL interop constants
# cv::ogl::RenderModes
# include "opengl.hpp"
const POINTS          = 0x0000
const LINES           = 0x0001
const LINE_LOOP       = 0x0002
const LINE_STRIP      = 0x0003
const TRIANGLES       = 0x0004
const TRIANGLE_STRIP  = 0x0005
const TRIANGLE_FAN    = 0x0006
const QUADS           = 0x0007
const QUAD_STRIP      = 0x0008
const POLYGON         = 0x0009

# Access
const READ_ONLY             = 0x88B8
const WRITE_ONLY            = 0x88B9
const READ_WRITE            = 0x88BA
# Target
const ARRAY_BUFFER          = 0x8892
const ELEMENT_ARRAY_BUFFER  = 0x8893
const PIXEL_PACK_BUFFER     = 0x88EB
const PIXEL_UNPACK_BUFFER   = 0x88EC


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Exodus"
version = v"0.1.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/gsjaardema/seacas.git", "1452c325ab5d507d000397c713a5b0c4dc57bf81")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/seacas

mkdir build
cd build

### The SEACAS code will install in ${INSTALL_PATH}/bin, ${INSTALL_PATH}/lib, and ${INSTALL_PATH}/include.
INSTALL_PATH=${prefix}

FORTRAN=NO
NETCDF_PATH=${prefix}
PNETCDF_PATH=${prefix}
HDF5_PATH=${prefix}

### Set to ON for parallel compile; otherwise OFF for serial (default)
if [ "${MPI}" = "" ]; then
  netcdf_parallel=$($NETCDF_PATH/bin/nc-config --has-parallel)
  if [ "${netcdf_parallel}" == "yes" ]; then
    MPI=YES
  else
    MPI=NO
  fi
fi

CFLAGS="-Wall -Wunused -pedantic -std=c11"
CXXFLAGS="-Wall -Wunused -pedantic"

GENERATOR="Unix Makefiles"
SHARED="YES"
BUILD_TYPE="RELEASE"

### Set to YES to enable the building of a thread-safe version of the Exodus and IOSS libraries.
THREADSAFE=${THREADSAFE:-NO}

if [ "$THREADSAFE" == "YES" ] ; then
  THREAD_SAFE_OPT="-DSEACASProj_EXTRA_LINK_FLAGS=-lpthread"
fi

function check_enable()
{
  local path=$1
  if [ -e "${path}" ]
  then
    echo "YES"
  else
    echo "NO"
  fi
}

HAVE_NETCDF=$(check_enable "${NETCDF_PATH}/include/netcdf.h")

### Define to NO to *enable* exodus deprecated functions
OMIT_DEPRECATED=${OMIT_DEPRECATED:-NO}

NUMPROCS=${NUMPROCS:-4}

# BUG needs to work with cray too.
if [ "${MPI}" == "YES" ]; then
  if [ "${USE_SRUN}" == "YES" ]
  then
    MPI_EXEC=$(which srun)
    MPI_SYMBOLS="-D MPI_EXEC=${MPI_EXEC} -D MPI_EXEC_NUMPROCS_FLAG=-N  -DMPI_EXEC_DEFAULT_NUMPROCS:STRING=${NUMPROCS} -DMPI_EXEC_MAX_NUMPROCS:STRING=${NUMPROCS}"
    MPI_BIN=$(dirname "${MPI_EXEC}")
  else
    MPI_EXEC=$(which mpiexec)
    MPI_SYMBOLS="-D MPI_EXEC=${MPI_EXEC}  -DMPI_EXEC_DEFAULT_NUMPROCS:STRING=${NUMPROCS} -DMPI_EXEC_MAX_NUMPROCS:STRING=${NUMPROCS}"
    MPI_BIN=$(dirname "${MPI_EXEC}")
  fi
  CXX=mpicxx
  CC=mpicc
  FC=mpif77
fi

### You can add these below if you want more verbosity...
#-D CMAKE_VERBOSE_MAKEFILE:BOOL=ON \
#-D SEACASProj_VERBOSE_CONFIGURE=ON \

### You can add these below to regenerate the flex and bison files for
### aprepro and aprepro_lib May have to touch aprepro.l aprepro.y
### aprepro.ll and aprepro.yy to have them regenerate
#-D GENERATE_FLEX_FILES=ON \
#-D GENERATE_BISON_FILES=ON \

###------------------------------------------------------------------------
cmake -G "${GENERATOR}" \
    -D CMAKE_INSTALL_PREFIX:PATH=${prefix} \
    -D CMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -D CMAKE_CXX_COMPILER:FILEPATH=${CXX} \
    -D CMAKE_C_COMPILER:FILEPATH=${CC} \
    -D CMAKE_Fortran_COMPILER:FILEPATH=${FC} \
    -D CMAKE_CXX_FLAGS="${CXXFLAGS} ${CXX_WARNING_FLAGS}" \
    -D CMAKE_C_FLAGS="${CFLAGS} ${C_WARNING_FLAGS}" \
    -D CMAKE_Fortran_FLAGS="${FFLAGS} ${F77_WARNING_FLAGS}" \
    -D CMAKE_INSTALL_RPATH:PATH=${INSTALL_PATH}/lib \
    -D BUILD_SHARED_LIBS:BOOL=${SHARED} \
    -D CMAKE_BUILD_TYPE=${BUILD_TYPE} \
    -D SEACASProj_ENABLE_SEACASExodus=YES \
    -D SEACASProj_ENABLE_SEACASExodus_for=${FORTRAN} \
    -D SEACASProj_ENABLE_SEACASExoIIv2for32=${FORTRAN} \
    -D SEACASProj_ENABLE_TESTS=YES \
    -D SEACASExodus_ENABLE_STATIC:BOOL=${STATIC} \
    -D SEACASProj_SKIP_FORTRANCINTERFACE_VERIFY_TEST:BOOL=YES \
    -D SEACASProj_HIDE_DEPRECATED_CODE:BOOL=${OMIT_DEPRECATED_CODE} \
    -D SEACASProj_ENABLE_Fortran=${FORTRAN} \
    \
    -D TPL_ENABLE_Netcdf:BOOL=${HAVE_NETCDF} \
    -D TPL_ENABLE_MPI:BOOL=${MPI} \
    -D TPL_ENABLE_Pthread:BOOL=${THREADSAFE} \
    ${THREAD_SAFE_OPT} \
    -D SEACASExodus_ENABLE_THREADSAFE:BOOL=${THREADSAFE} \
    \
    ${MPI_SYMBOLS} \
    \
    -D MPI_BIN_DIR:PATH=${MPI_BIN} \
    -D NetCDF_ROOT:PATH=${NETCDF_PATH} \
    -D HDF5_ROOT:PATH=${HDF5_PATH} \
    -D HDF5_NO_SYSTEM_PATHS=YES \
    -D PNetCDF_ROOT:PATH=${PNETCDF_PATH} \
    \
    ..

echo ""
echo "INSTALL_PATH: ${INSTALL_PATH}"
echo "  "
echo "          CC: ${CC}"
echo "         CXX: ${CXX}"
echo "          FC: ${FC}"
echo "         MPI: ${MPI}"
echo "      SHARED: ${SHARED}"
echo "  BUILD_TYPE: ${BUILD_TYPE}"
echo "  THREADSAFE: ${THREADSAFE}"
echo "  "
echo " HAVE_NETCDF: ${HAVE_NETCDF}"
echo ""

make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc = "glibc"),
    Platform("aarch64", "linux"; libc = "glibc"),
    Platform("x86_64", "macos"),
    Platform("aarch64", "macos"),
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libexodus", :libexodus)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a"))
    Dependency(PackageSpec(name="NetCDF_jll", uuid="7243133f-43d8-5620-bbf4-c2c921802cf3"))
    Dependency(PackageSpec(name="HDF5_jll", uuid="0234f1f7-429e-5d53-9886-15a909be8d59"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
using BlockBandedMatrices, LinearAlgebra, MatrixFactorizations

@testset "BlockBandedMatrix QR/QL" begin
    @testset "Square QR" begin
        N = 5
        A = BlockBandedMatrix{Float64}(undef, (1:N,1:N), (2,1))
        A.data .= randn.()

        F = qr(A)
        @test F.factors ≈ qr(Matrix(A)).factors
        @test F.τ ≈ LinearAlgebra.qrfactUnblocked!(Matrix(A)).τ

        Q,R = F
        Q̃,R̃ = qr(Matrix(A))
        @test R ≈ R̃
        @test Q ≈ Q̃

        b = randn(size(A,1))
        @test Q'b ≈ Q̃'b
        @test Q*b ≈ Q̃*b

        @test F\b ≈ ldiv!(F, copy(b)) ≈ Matrix(A)\b ≈ A\b
    end

    @testset "Thin QR" begin
        N = 5
        A = BlockBandedMatrix{Float64}(undef, (1:N+1,1:N), (2,1))
        A.data .= randn.()

        F = qr(A)
        @test F.factors ≈ qr(Matrix(A)).factors
        @test F.τ ≈ LinearAlgebra.qrfactUnblocked!(Matrix(A)).τ

        Q,R = F
        Q̃,R̃ = qr(Matrix(A))
        @test R ≈ R̃
        @test Q ≈ Q̃

        b = randn(size(A,1))
        @test Q'b ≈ Q̃'b
        @test Q*b ≈ Q̃*b

        @test Q'b ≈ lmul!(Q', copy(b)) ≈ Q̃'*b

        @test F\b ≈ ldiv!(F, copy(b))[1:15] ≈ Matrix(A)\b ≈ A\b
    end

    @testset "Wide QR" begin
        N = 5
        A = BlockBandedMatrix{Float64}(undef, (1:N,1:N+1), (2,1))
        A.data .= randn.()

        F = qr(A)
        @test F.factors ≈ qr(Matrix(A)).factors
        @test F.τ ≈ LinearAlgebra.qrfactUnblocked!(Matrix(A)).τ

        Q,R = F
        Q̃,R̃ = qr(Matrix(A))
        @test R ≈ R̃
        @test Q ≈ Q̃

        b = randn(size(A,1))
        @test Q'b ≈ Q̃'b
        @test Q*b ≈ Q̃*b
        Q'b
        # @test_broken DimensionMismatch ldiv!(F, copy(b))

        @test_broken F\b ≈ ldiv!(F, copy(b)) ≈ Matrix(A)\b ≈ A\b
    end

    @testset "Square QL" begin
        N = 5
        A = BlockBandedMatrix{Float64}(undef, (1:N,1:N), (2,1))
        A.data .= randn.()

        F = ql(A)
        @test F.factors ≈ ql(Matrix(A)).factors
        @test F.τ ≈ MatrixFactorizations.qlfactUnblocked!(Matrix(A)).τ

        Q,L = F
        Q̃,L̃ = ql(Matrix(A))
        @test L ≈ L̃
        @test Q ≈ Q̃

        b = randn(size(A,1))
        @test Q'b ≈ Q̃'b
        @test Q*b ≈ Q̃*b

        b = randn(size(A,1))
        @test F\b ≈ ldiv!(F, copy(b)) ≈ Matrix(A)\b ≈ A\b
    end

    @testset "Thin QL" begin
        N = 5
        A = BlockBandedMatrix{Float64}(undef, (1:N+1,1:N), (2,1))
        A.data .= randn.()

        # F = ql(A)
        # @test F.factors ≈ ql(Matrix(A)).factors
        # @test F.τ ≈ MatrixFactorizations.qlfactUnblocked!(Matrix(A)).τ

        # Q,L = F
        # Q̃,L̃ = ql(Matrix(A))
        # @test L ≈ L̃
        # @test Q ≈ Q̃

        # b = randn(size(A,1))
        # @test Q'b ≈ Q̃'b
        # @test Q*b ≈ Q̃*b

        # b = randn(size(A,1))
        # @test F\b ≈ ldiv!(F, copy(b)) ≈ Matrix(A)\b ≈ A\b
    end

    @testset "Wide QL" begin
        N = 5
        A = BlockBandedMatrix{Float64}(undef, (1:N,1:N+1), (2,1))
        A.data .= randn.()

        @test_throws ArgumentError ql(A)
    end

    @testset "Complex QR/QL" begin
        @testset "Simple" begin
            A=BlockBandedMatrix{Float64}(I, ([1,1],[1,1]), (0,0))
            Af=qr(A)
            B=BlockBandedMatrix{ComplexF64}(I, ([1,1],[1,1]), (0,0))
            Bf=qr(B)
            @test Af.factors == Bf.factors
            @test Af.τ == Bf.τ
        end

        @testset "Off-diagonals" begin
            for qrl in [qr,ql]
                for T in [Float64, ComplexF64]
                    A = BlockBandedMatrix{T}(undef, ([3,2,1,3],[3,2,1,3]), (1,1))
                    A.data .= rand(T, size(A.data))
                    Af = qrl(A)
                    @test Af.factors isa BlockBandedMatrix

                    Am = Matrix(A)
                    Amf = qrl(Am)

                    @test Af.factors ≈ Amf.factors
                end
            end
        end
    end
end

# N = 500;
# A = BlockBandedMatrix{Float64}(undef, (1:N,1:N), (2,1));
# A.data .= randn.();
# @time F = qr(A); # 11s
# b = randn(size(A,1));
# @time F\b; # 0.6s


text: raw(t::Unitful.Time) = Float64(ustrip(s(t)))


text:
module TurbulenceClosures

export
    AbstractEddyViscosityClosure,
    IsotropicDiffusivity,
    AnisotropicDiffusivity,
    AnisotropicBiharmonicDiffusivity,
    TwoDimensionalLeith,
    ConstantSmagorinsky,
    SmagorinskyLilly,
    AnisotropicMinimumDissipation,
    HorizontallyCurvilinearAnisotropicDiffusivity,

    DiffusivityFields,
    calculate_diffusivities!,

    ∇_κ_∇c,
    ∇_κ_∇T,
    ∇_κ_∇S,
    ∂ⱼ_2ν_Σ₁ⱼ,
    ∂ⱼ_2ν_Σ₂ⱼ,
    ∂ⱼ_2ν_Σ₃ⱼ,

    cell_diffusion_timescale

using CUDA
using KernelAbstractions

import Oceananigans.Utils: with_tracers

using Oceananigans
using Oceananigans.Grids
using Oceananigans.Operators
using Oceananigans.BoundaryConditions
using Oceananigans.Fields
using Oceananigans.Buoyancy
using Oceananigans.Utils

using Oceananigans.Architectures: AbstractArchitecture, device, @hascuda

#####
##### Molecular viscosity and thermal diffusivity definitions
##### Approximate viscosities and thermal diffusivities for seawater at 20ᵒC and 35 psu,
##### according to Sharqawy et al., "Thermophysical properties of seawater: A review of
##### existing correlations and data" (2010).
#####

const ν₀ = 1.05e-6
const κ₀ = 1.46e-7

#####
##### Abstract types
#####

"""
    AbstractTurbulenceClosure

Abstract supertype for turbulence closures.
"""
abstract type AbstractTurbulenceClosure end

"""
    AbstractTensorDiffusivity <: AbstractTurbulenceClosure

Abstract supertype for turbulence closures that are defined by a tensor viscosity and
tensor diffusivities.
"""
abstract type AbstractTensorDiffusivity <: AbstractTurbulenceClosure end

#####
##### 'Tupled closure' implementation
#####

for stress_div in (:∂ⱼ_2ν_Σ₁ⱼ, :∂ⱼ_2ν_Σ₂ⱼ, :∂ⱼ_2ν_Σ₃ⱼ)
    @eval begin
        @inline function $stress_div(i, j, k, grid::AbstractGrid{FT}, clock, closure_tuple::Tuple, U,
                                     K_tuple, args...) where FT

            stress_div_ijk = zero(FT)

            ntuple(Val(length(closure_tuple))) do α
                @inbounds closure = closure_tuple[α]
                @inbounds K = K_tuple[α]
                stress_div_ijk += $stress_div(i, j, k, grid, clock, closure, U, K, args...)
            end

            return stress_div_ijk
        end
    end
end

@inline function ∇_κ_∇c(i, j, k, grid::AbstractGrid{FT}, clock, closure_tuple::Tuple,
                        c, tracer_index, K_tuple, args...) where FT

    flux_div_ijk = zero(FT)

    ntuple(Val(length(closure_tuple))) do α
        @inbounds closure = closure_tuple[α]
        @inbounds K = K_tuple[α]
        flux_div_ijk +=  ∇_κ_∇c(i, j, k, grid, clock, closure, c, tracer_index, K, args...)
    end

    return flux_div_ijk
end

function calculate_diffusivities!(K_tuple::Tuple, arch, grid, closure_tuple::Tuple, args...)
    ntuple(Val(length(closure_tuple))) do α
        @inbounds closure = closure_tuple[α]
        @inbounds K = K_tuple[α]
        calculate_diffusivities!(K, arch, grid, closure, args...)
    end

    return nothing
end

with_tracers(tracers, closure_tuple::Tuple) =
    Tuple(with_tracers(tracers, closure) for closure in closure_tuple)

#####
##### Include module code
#####

include("turbulence_closure_utils.jl")
include("diffusion_operators.jl")
include("viscous_dissipation_operators.jl")
include("velocity_tracer_gradients.jl")
include("abstract_eddy_viscosity_closures.jl")

include("closure_tuples.jl")

include("turbulence_closure_implementations/nothing_closure.jl")
include("turbulence_closure_implementations/isotropic_diffusivity.jl")
include("turbulence_closure_implementations/anisotropic_diffusivity.jl")
include("turbulence_closure_implementations/horizontally_curvilinear_anistropic_diffusivity.jl")
include("turbulence_closure_implementations/anisotropic_biharmonic_diffusivity.jl")

# Eddy viscosity closures
include("turbulence_closure_implementations/leith_enstrophy_diffusivity.jl")
include("turbulence_closure_implementations/smagorinsky_lilly.jl")
include("turbulence_closure_implementations/anisotropic_minimum_dissipation.jl")

include("diffusivity_fields.jl")
include("turbulence_closure_diagnostics.jl")

#####
##### Some value judgements here
#####

"""
    ConstantSmagorinsky

An alias for `SmagorinskyLilly`.
"""
const ConstantSmagorinsky = SmagorinskyLilly

end # module


text:
using Statistics,Serialization,StatsPlots, DataFrames

# parameter values
paras = ["beta_l5","alpha_o050","beta_l5_alpha_l025", "baseline_GER", "beta_l5_alpha_l025_alpha_o025"]

path_to_data_prefix = "../../data/working_paper_2/"

filename_prefix = "../../figures/working_paper_2/"

mkpath(filename_prefix)

# define period of plotting for economic plots
start_week = 3 # cut of first week
end_week = 102
# set periods for covid plots
inf_start_day = 1
inf_end_day = 600 # cut of after period 400
inf_start_week = 3
inf_end_week = 86 # 58*7 = 406 and 42*7 = 294

include("$(path_to_data_prefix)baseline_GER/covid_par_ini.jl") # necessary for some basic parameter values (nhh)

datapoint = fld(T,datat)

# scatterplot
gdploss_means = []
totcas_means = []
pubacc_means = []

# boxplots
number_batches = 20
gdploss_all_means = DataFrame()
totcas_all_means = DataFrame()
pubacc_all_means = DataFrame()
badpoltime_all_means = DataFrame()
polsw_all_means = DataFrame()
totinfmean_all_means = DataFrame()
totinfstd_all_means_plus = DataFrame()
totinfstd_all_means_minus = DataFrame()
infmean_all_means = DataFrame()
infstd_all_means_plus = DataFrame()
infstd_all_means_minus = DataFrame()
totcasmean_all_means = DataFrame()
totcasstd_all_means_plus = DataFrame()
totcasstd_all_means_minus = DataFrame()
gdpmean_all_means = DataFrame()
gdpstd_all_means_plus = DataFrame()
gdpstd_all_means_minus = DataFrame()
unempmean_all_means = DataFrame()
unempstd_all_means_plus = DataFrame()
unempstd_all_means_minus = DataFrame()
pubaccmean_all_means = DataFrame()
pubaccmean_all_means_plus = DataFrame()
pubaccmean_all_means_minus = DataFrame()
totfirmsmean_all_means = DataFrame()

for (i,p) in enumerate(paras)
    local path_to_data = string(path_to_data_prefix, p , "/")
    # store data scatterplot
    local worker_results = deserialize(open("$(path_to_data)batchdata.dat"))
    gdploss = mean((results -> results[:togdploss]).(worker_results))
    totcas = mean(100*(results -> results[:totcas]).(worker_results))
    pubacch = (results -> results[:pubacctraj]).(worker_results)
    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    global nn = size(gdpmean)[1]
    local pubacc = mean([pubacch[i][datapoint+1]*datat/(365*gdpmean[2]) for  i = 1:size(pubacch)[1]])
    push!(gdploss_means,gdploss)
    push!(totcas_means,totcas)
    push!(pubacc_means,pubacc)
    # store data boxplot
    gdploss_all = (results -> results[:togdploss]).(worker_results)
    totcas_all = 100*(results -> results[:totcas]).(worker_results)
    pubacc_all = [pubacch[i][datapoint+1]*datat/(365*gdpmean[2]) for  i = 1:size(pubacch)[1]]
    gdploss_all_means[!, Symbol.(paras[i])] = gdploss_all
    totcas_all_means[!, Symbol.(paras[i])] = totcas_all
    pubacc_all_means[!, Symbol.(paras[i])] = pubacc_all
    badpoltime_all = (results -> results[:badpoltime]).(worker_results)
    polsw_all = (results -> results[:polswitchcount]).(worker_results)
    badpoltime_all_means[!, Symbol.(paras[i])] = badpoltime_all
    polsw_all_means[!, Symbol.(paras[i])] = polsw_all

    # store data dynamics
    totinfmean_all = vec(mean(hcat((results -> results[:totinftraj]).(worker_results)...), dims=2))
    totinfstd_all = vec(std(hcat((results -> results[:totinftraj]).(worker_results)...), dims=2))
    totinfmean_all_means[!, Symbol.(paras[i])] = totinfmean_all
    totinfstd_all_means_plus[!, Symbol.(paras[i])] = totinfmean_all + totinfstd_all
    totinfstd_all_means_minus[!, Symbol.(paras[i])] = totinfmean_all - totinfstd_all
    infmean = vec(mean(100*cat((results -> results[:inftraj]).(worker_results)...,dims=3), dims=3)[2:nn,3])
    infstd = vec(std(100*cat((results -> results[:inftraj]).(worker_results)...,dims=3), dims=3)[2:nn,3])
    infmean_all_means[!, Symbol.(paras[i])] = infmean
    infstd_all_means_plus[!, Symbol.(paras[i])] = infmean + infstd
    infstd_all_means_minus[!, Symbol.(paras[i])] = infmean - infstd
    global totcasmean = vec(mean(cat(100*(results -> results[:totcastraj]).(worker_results)...,dims=3), dims=3)[2:nn,3])
    global totcasstd = vec(std(cat(100*(results -> results[:totcastraj]).(worker_results)...,dims=3), dims=3)[2:nn,3])
    totcasmean_all_means[!, Symbol.(paras[i])] = totcasmean
    totcasstd_all_means_plus[!, Symbol.(paras[i])] = totcasmean + totcasstd
    totcasstd_all_means_minus[!, Symbol.(paras[i])] = totcasmean - totcasstd
    # gdp
    gdpmean_all = vec(mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2))
    gdpstd_all = vec(std(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2))
    gdpmean_all_means[!, Symbol.(paras[i])] = gdpmean_all
    gdpstd_all_means_plus[!, Symbol.(paras[i])] = gdpmean_all + gdpstd_all
    gdpstd_all_means_minus[!, Symbol.(paras[i])] = gdpmean_all - gdpstd_all
    # unemployment
    unempmean_all = vec(mean(hcat((results -> results[:unempltraj]).(worker_results)...), dims=2))
    unempstd_all = vec(std(hcat((results -> results[:unempltraj]).(worker_results)...), dims=2))
    unempmean_all_means[!, Symbol.(paras[i])] = unempmean_all
    unempstd_all_means_plus[!, Symbol.(paras[i])] = unempmean_all + unempstd_all
    unempstd_all_means_minus[!, Symbol.(paras[i])] = unempmean_all - unempstd_all
    # public account, adjusted
    pubh = (datat/(365*gdpmean[2]))*hcat((results -> results[:pubacctraj]).(worker_results)...)
    pubmean = vec(mean(pubh, dims=2))
    pubstd = vec(std(pubh, dims=2))
    pubaccmean_all_means[!, Symbol.(paras[i])] = pubmean
    pubaccmean_all_means_plus[!, Symbol.(paras[i])] = pubmean + pubstd
    pubaccmean_all_means_minus[!, Symbol.(paras[i])] = pubmean - pubstd
    #totfirms
    totfirmsmean_all = vec(mean(hcat((results -> results[:totfirmtraj]).(worker_results)...), dims=2))
    totfirmsstd_all = vec(std(hcat((results -> results[:totfirmtraj]).(worker_results)...), dims=2))
    totfirmsmean_all_means[!, Symbol.(paras[i])] = totfirmsmean_all
end

hcapfrac = zeros(datapoint+1)
for th = 1:datapoint+1
    hcapfrac[th] = hcap / (nhh * icufrac)
end

# Infections
pl8 = @df infmean_all_means[inf_start_week:inf_end_week, :] plot(datat*collect(inf_start_week:inf_end_week), cols(1), linewidth = [2],   legend = false, xlabel = "Days", ylabel = "infected [%]", label = false, linecolor = 1)
@df infstd_all_means_plus[inf_start_week:inf_end_week, :] plot!(datat*collect(inf_start_week:inf_end_week), cols(1), linewidth = [1], linestyle= [:dot], legend = false,linecolor = 1 )
@df infstd_all_means_minus[inf_start_week:inf_end_week, :] plot!(datat*collect(inf_start_week:inf_end_week), cols(1), linestyle= [:dot], linewidth = [1], legend = false,linecolor = 1)

global i=0
for p in paras
    global i = i+1
    @df infmean_all_means[inf_start_week:inf_end_week, :] plot!(datat*collect(inf_start_week:inf_end_week), cols(i), linewidth = [2],  legend = false, label = false, linecolor = i)
    @df infstd_all_means_plus[inf_start_week:inf_end_week, :] plot!(datat*collect(inf_start_week:inf_end_week), cols(i), linewidth = [1], linestyle= [:dot], linecolor = i)
    @df infstd_all_means_minus[inf_start_week:inf_end_week, :] plot!(datat*collect(inf_start_week:inf_end_week), cols(i), linestyle= [:dot], linewidth = [1], linecolor = i)
end
savefig(pl8,"$(filename_prefix)policy_comparison_infection_dynamics.pdf")

# Casulties
pl9 = @df totcasmean_all_means[inf_start_week:inf_end_week, :] plot(datat*collect(inf_start_week:inf_end_week),cols(1), linewidth = [2],   legend = false, xlabel = "Days", ylabel = "casualties [%]", label = false, linecolor =1)
@df totcasstd_all_means_plus[inf_start_week:inf_end_week, :] plot!(datat*collect(inf_start_week:inf_end_week),cols(1), linewidth = [1], linestyle= [:dot], legend = false,linecolor =1)
@df totcasstd_all_means_minus[inf_start_week:inf_end_week, :] plot!(datat*collect(inf_start_week:inf_end_week),cols(1), linestyle= [:dot], linewidth = [1], legend = false,linecolor = 1)
global i=0
for p in paras
    global i = i+1
    @df totcasmean_all_means[inf_start_week:inf_end_week, :] plot!(datat*collect(inf_start_week:inf_end_week), cols(i), linewidth = [2],  legend = false, label = false, linecolor = i)
    @df totcasstd_all_means_plus[inf_start_week:inf_end_week, :] plot!(datat*collect(inf_start_week:inf_end_week), cols(i), linewidth = [1], linestyle= [:dot], linecolor = i)
    @df totcasstd_all_means_minus[inf_start_week:inf_end_week, :] plot!(datat*collect(inf_start_week:inf_end_week), cols(i), linestyle= [:dot], linewidth = [1], linecolor = i)
end
savefig(pl9,"$(filename_prefix)policy_comparison_casulties_dynamics.pdf")

# GDP
pl10 = @df gdpmean_all_means[start_week:end_week, :] plot(datat*collect(start_week:end_week), cols(1),  linewidth = [2],  legend = false, xlabel = "Days", ylabel = "GDP",label = false, linecolor =1)
@df gdpstd_all_means_plus[start_week:end_week, :] plot!(datat*collect(start_week:end_week), cols(1), linewidth = [1], linestyle= [:dot], legend = false,linecolor = 1, label = ["" "" ""])
@df gdpstd_all_means_minus[start_week:end_week, :] plot!(datat*collect(start_week:end_week), cols(1), linestyle= [:dot], linewidth = [1], legend = false,linecolor = 1, label = ["" "" ""])
global i=0
for p in paras
    global i = i+1
    @df gdpmean_all_means[start_week:end_week, :] plot!(datat*collect(start_week:end_week), cols(i), linewidth = [2],  legend = false, label = false, linecolor = i)
    @df gdpstd_all_means_plus[start_week:end_week, :] plot!(datat*collect(start_week:end_week), cols(i), linewidth = [1], linestyle= [:dot], linecolor = i)
    @df gdpstd_all_means_minus[start_week:end_week, :] plot!(datat*collect(start_week:end_week), cols(i), linestyle= [:dot], linewidth = [1], linecolor = i)
end
savefig(pl10,"$(filename_prefix)policy_comparison_gdp_dynamics.pdf")

# Unemployment
pl11 = @df unempmean_all_means[start_week:end_week, :] plot(datat*collect(start_week:end_week), cols(1),   linewidth = [2], xlabel = "Days", ylabel = "unempl. [%]", legend = false, linecolor =1)
@df unempstd_all_means_plus[start_week:end_week, :] plot!(datat*collect(start_week:end_week), cols(1), linewidth = [1], linestyle= [:dot],  linecolor = 1,legend = false, label = ["" "" ""])
@df unempstd_all_means_minus[start_week:end_week, :] plot!(datat*collect(start_week:end_week), cols(1), linestyle= [:dot], linewidth = [1], linecolor = 1,legend = false, label = ["" "" ""])
global i=0
for p in paras
    global i = i+1
    @df unempmean_all_means[start_week:end_week, :] plot!(datat*collect(start_week:end_week), cols(i), linewidth = [2],  legend = false, label = false, linecolor = i)
    @df unempstd_all_means_plus[start_week:end_week, :] plot!(datat*collect(start_week:end_week), cols(i), linewidth = [1], linestyle= [:dot], linecolor = i)
    @df unempstd_all_means_minus[start_week:end_week, :] plot!(datat*collect(start_week:end_week), cols(i), linestyle= [:dot], linewidth = [1], linecolor = i)
end
savefig(pl11,"$(filename_prefix)policy_comparison_unemployment_dynamics.pdf")

# public account
pl12 = @df pubaccmean_all_means[start_week:end_week, :] plot(datat*collect(start_week:end_week), cols(1), title = "Public Account",  linewidth = [2],  legend = false, linecolor =1)
@df pubaccmean_all_means_plus[start_week:end_week, :] plot!(datat*collect(start_week:end_week), cols(1), linewidth = [1], linestyle= [:dot], linecolor = 1, legend = false,label = ["" "" ""])
@df pubaccmean_all_means_minus[start_week:end_week, :] plot!(datat*collect(start_week:end_week), cols(1), linestyle= [:dot], linewidth = [1], linecolor = 1,legend = false, label = ["" "" ""])
global i=0
for p in paras
    global i = i+1
    @df pubaccmean_all_means[start_week:end_week, :] plot!(datat*collect(start_week:end_week), cols(i), linewidth = [2],  legend = false, label = false, linecolor = i)
    @df pubaccmean_all_means_plus[start_week:end_week, :] plot!(datat*collect(start_week:end_week), cols(i), linewidth = [1], linestyle= [:dot], linecolor = i)
    @df pubaccmean_all_means_minus[start_week:end_week, :] plot!(datat*collect(start_week:end_week), cols(i), linestyle= [:dot], linewidth = [1], linecolor = i)
end
savefig(pl12,"$(filename_prefix)policy_comparison_public_account_dynamics.pdf")


text:
import Celeste.Stripe82Score

function test_infer_field_and_score_object()
    rcf = RunCamcolField(4263, 5, 119)
    objid = "1237663784734490644"
    ParallelRun.infer_field(rcf, datadir, datadir; objid=objid)
    truthfile = joinpath(datadir, "coadd_for_4263_5_119.fit")
    Stripe82Score.score_object_disk(rcf, objid, datadir, truthfile, datadir)
end

function test_score_field()
    results_filename = "celeste-004263-5-0119.jld"

    if !isfile(joinpath(datadir, results_filename))
        results_url = "https://www.dropbox.com/s/8jcdmluahsf38tm/celeste-004263-5-0119.jld"
        run(`wget --quiet -P $datadir $results_url`)
    end

    rcf = RunCamcolField(4263, 5, 119)
    truthfile = joinpath(datadir, "coadd_for_4263_5_119.fit")
    Stripe82Score.score_field_disk(rcf, datadir, truthfile, datadir)
end

test_score_field()
test_infer_field_and_score_object()


text:
VERSION >= v"0.7.0-" && using InteractiveUtils
versioninfo()

if VERSION < v"0.7.0-"
    Pkg.clone(pwd(), "LibSndFile")
    Pkg.build("LibSndFile")
    # for now we need SampledSignals master
    Pkg.checkout("SampledSignals")
else
    using Pkg
    # for now we need to `clone` because there's no way to specify the
    # package name for `add`
    Pkg.clone(pwd(), "LibSndFile")
    Pkg.build("LibSndFile")
    Pkg.add(PackageSpec(name="SampledSignals", rev="master"))
end
# manually install test dependencies so we can run the test script directly, which avoids
# clobberling our environment
Pkg.add("Compat")
Pkg.add("FileIO")


text:
struct OrgParseError
    content::AbstractString
    point::Integer
    allowedtypes::Vector{<:Type}
end

function Base.showerror(io::IO, ex::OrgParseError)
    beforepoint = if ex.point > 10 "…" else "" end *
        replace(ex.content[max(1, prevind(ex.content, ex.point, 10)):ex.point-1], '\n' => "\\n")
    afterpoint = replace(ex.content[ex.point:min(end, nextind(ex.content, ex.point, 45))], '\n' => "\\n") *
        if ex.point + 45 < length(ex.content) "…" else "" end
    print(io, "Org parse failed at index $(ex.point) (char $(length(ex.content[1:ex.point]))) of string:\n")
    printstyled(io, " ", beforepoint, color=:light_black)
    printstyled(io, "", afterpoint, "\n", color=:yellow)
    printstyled(io, " "^length(beforepoint), " ^", bold=true, color=:red)
    print(io, "\nThe immediate remaining content did not match any of the allowed components:\n")
    print(io, " • ", join(string.(ex.allowedtypes), "\n • "), "\n")
    Base.Experimental.show_error_hints(io, ex)
end

function parseorg(content::AbstractString, typematchers::Dict{Char, <:AbstractVector{<:Type}},
                  typefallbacks::AbstractVector{<:Type};
                  debug::Bool=false, partial::Bool=false, maxobj::Integer=0)
    point, objects = 1, OrgComponent[]
    points = [point]
    clen = lastindex(content) # this does not change, help the compiler
    while point <= clen
        if debug print("\n\e[36m$(lpad(point, 4))\e[37m") end
        obj::Union{OrgComponent, Nothing} = nothing
        types = DataType[]
        # use the next non-whitespace char to guide type matchers
        char = if content[point] == ' ' || content[point] == '\t'
            nextnonspace = findfirst(c -> c != ' ' && c != '\t',
                                     @inbounds @view content[point:clen])
            if !isnothing(nextnonspace)
                content[point + nextnonspace - 1]
            else
                content[point]
            end
        else
            content[point]
        end
        if char in keys(typematchers)
            types = typematchers[char]
            for type in types
                # profiling indicates that @inbounds @view content[point:clen] is about a
                # third faster than @inbounds @view content[point:end] for large strings
                res = consume(type, @inbounds @view content[point:clen])
                if !isnothing(res)
                    width, obj = res
                    point += width
                    break
                end
            end
        end
        if isnothing(obj)
            for type in typefallbacks
                res = consume(type, @inbounds @view content[point:clen])
                if !isnothing(res)
                    width, obj = res
                    point += width
                    break
                end
            end
        end
        if isnothing(obj)
            if partial
                return point-1, objects
            elseif content isa SubString
                throw(OrgParseError(content.string, content.offset + point,
                                    vcat(types, typefallbacks)))
            else
                throw(OrgParseError(content, point, vcat(types, typefallbacks)))
            end
        end
        if obj isa TextPlain &&
            length(objects) > 0 &&
            objects[end] isa TextPlain
            objects[end] *= obj
        else
            push!(objects, obj)
        end
        push!(points, point)
        if debug
            print(rpad(" ─[\e[33m$(typeof(objects[end]))\e[37m" *
                "(\e[35m$(length(objects))\e[37m)]─", 42, '─'),
                  "> \e[36m$(lpad(point, 4))\e[37m\n",
                  "     \e[32m'", rpad(content[points[end-1]:prevind(content, points[end])] * "'", 28),
                  "\e[37m  ", objects[end])
        end
        if partial && 0 < length(objects) == maxobj
            return point-1, objects
        end
    end
    if debug print("\n\n") end
    if partial
        point-1, objects
    else
        objects
    end
end

function parseorg(content::AbstractString, typefallbacks::AbstractVector{<:Type};
                  debug=false, partial=false)
    parseorg(content, Dict{Char, Vector{DataType}}(),
             typefallbacks; debug, partial)
end

# Parsing objects, with restrictions

filtermatcher(allowed::Vector{<:Type}, matchers=org_object_matchers) =
    filter(p -> !isempty(p.second),
           Dict{Char, Vector{<:Type}}(key => filter(v -> v ∈ allowed, value)
                                      for (key, value) in matchers))

const org_object_resticted_matchers =
    Dict(o => filtermatcher(org_object_restrictions[o])
         for o in keys(org_object_restrictions))

const org_object_resticted_fallbacks =
    Dict(o => [org_object_fallbacks ∩ org_object_restrictions[o]; TextPlainForced]
         for o in keys(org_object_restrictions))

function parseobjects(context::Type, content::AbstractString; debug=false, partial=false)
    @assert context in keys(org_object_restrictions)
    Vector{Object}(
        parseorg(content, org_object_resticted_matchers[context],
                 org_object_resticted_fallbacks[context]; debug, partial))
end

# parsing utilities

function forwardsbalenced(content::AbstractString, point::Integer=1, limit::Integer=lastindex(content);
                          bracketpairs::Dict{Char, Char}=Dict{Char, Char}(), escapechars::Vector{Char}=Char[],
                          quotes::Vector{Char}=Char[], spacedquotes::Vector{Char}=Char[])
    open = content[point]
    if !(open in keys(bracketpairs))
        throw(ErrorException("index $point does not lie on a recognised open bracket: $(join(keys(bracketpairs), ", "))"))
    end
    close = bracketpairs[open]
    currentquote = nothing
    depth = 1
    while point < limit
        point = nextind(content, point)
        # escaped chars
        if content[point] in escapechars &&
            point < limit
            point += 1
            continue
        end
        # unquote
        if !isnothing(currentquote)
            if content[point] == currentquote
                if currentquote in quotes
                    currentquote = nothing
                    continue
                elseif currentquote in spacedquotes &&
                    (point == limit || content[point+1] == ' ')
                    currentquote = nothing
                    continue
                end
            else
                continue
            end
        end
        # quotes
        if content[point] in quotes
            currentquote = content[point]
            continue
        end
        # quotes with space
        if content[point] in spacedquotes &&
            (point == 1 || content[point-1] == ' ')
            currentquote = content[point]
            continue
        end
        # brackets
        if content[point] == open
            depth += 1
        elseif content[point] == close
            depth -= 1
        end
        if depth == 0
            return point
        end
    end
end


text:
abstract type AbstractUIElement end
abstract type AbstractUIHTMLElement <: AbstractUIElement end

struct UIElement <: AbstractUIHTMLElement
    html::String
end

struct UIHeader <: AbstractUIHTMLElement
    html::String
end

struct UIFooter <: AbstractUIHTMLElement
    html::String
end

struct UIModel <: AbstractUIElement
    id::String
    default::String
end

struct UIWatch <: AbstractUIElement
    id::String
    code::String
end

UIHeader(el::UIElement) = el
UIFooter(el::UIElement) = el
UIModel(el::UIModel) = el

UIHeader(els::Tuple) = unroll(els)
UIFooter(els::Tuple) = unroll(els)
UIModel(els::Tuple) = unroll(els)

tag_wrap(tag, contents...) = (UIElement("<$tag>"), contents, UIElement("</$tag>"))

function extract_uitype(t::Type, el::AbstractUIElement)
    if typeof(el) <: t
        (el)
    else
        ()
    end
end

function extract_uitype(t::Type, content::Tuple)
    out = ()
    for c in content
        if typeof(c) <: t
            out = (out..., c)
        end
    end
    out
end

convert_uielement(t::Type, el::AbstractUIHTMLElement) = t(el.html)
convert_uielement(t::Type, el::UIModel) = el # no op

function convert_uielement(t::Type, content::Tuple)
    map(c -> convert_uielement(t, c), content)
end

function unroll(content::Tuple)
    out = ()
    for c in content
        if typeof(c) <: AbstractUIElement
            out = (out..., c)
        elseif typeof(c) <: AbstractString
            out = (out..., UIElement(c))
        else
            out = (out..., unroll(c)...)
        end
    end
    out
end

function unroll(content::String)
    # Hit if UI is just a string
    UIElement(content)
end


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using GLMakie, NetCDF
using YAXArrays.Datasets: open_mfdataset
GLMakie.activate!() # HIDE
let
    dataset = open_mfdataset("./data/ETOPO1_halfdegree.nc")
    lon = dataset.lon.values
    lat = dataset.lat.values
    data = dataset.ETOPO1avg.data[:,:]
    function toCartesian(lon, lat; r = 1, cxyz = (0,0,0) )
        lat, lon = lat*π/180, lon*π/180
        x = cxyz[1] + r * cos(lat) * cos(lon)
        y = cxyz[2] + r * cos(lat) * sin(lon)
        z = cxyz[3] + r *sin(lat)
        x,y,z
    end
    function lonlat3D(lon, lat, data; cxyz = (0,0,0))
        xyzw = zeros(size(data)..., 3)
        for (i,lon) in enumerate(lon), (j,lat) in enumerate(lat)
            x, y, z = toCartesian(lon, lat; cxyz = cxyz)
            xyzw[i,j,1] = x
            xyzw[i,j,2] = y
            xyzw[i,j,3] = z
        end
        xyzw[:,:,1], xyzw[:,:,2], xyzw[:,:,3]
    end
    # this is needed in order to have a closed surface
    lonext  = cat(collect(lon), lon[1], dims = 1)
    dataext = begin
        tmpdata = zeros(size(lon)[1]+1,  size(lat)[1])
        tmpdata[1:size(lon)[1],:] = data
        tmpdata[size(lon)[1]+1,:] = data[1,:]
        tmpdata
    end
    #now let's plot that surface!
    xetopo, yetopo, zetopo = lonlat3D(lonext, lat, dataext)

    fig = Figure(resolution =(950,900),fontsize = 25,backgroundcolor = :black)
    ax = Axis3(fig, aspect = :data, azimuth = -0.2π, elevation = -0.05π)
    pltobj=surface!(ax,xetopo, yetopo, zetopo, color = dataext,
        colormap = :hot, backlight = 1.0f0, colorrange =(-6000,5000))
    cbar  =Colorbar(fig,pltobj,label = "ETOPO1 [m]",ticklabelsize = 20,
        tickalign =1, vertical = true, ticksize=15, labelcolor = :white,
        ticklabelcolor = :white, tickcolor = :grey90, height = Relative(0.5))
    hidespines!(ax)
    hidedecorations!(ax)
    fig[1,1] = ax
    fig[1,2] = cbar
    fig
    save(joinpath(@__DIR__, "output", "topographySphere.png"), fig, px_per_unit = 2.0) # HIDE
end

using Pkg # HIDE
Pkg.status(["GLMakie","NetCDF", "YAXArrays"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using DataFrames, ROOT, ROOTDataFrames
using Base.Test

df = TreeDataFrame([string(Pkg.dir(), "/ROOTDataFrames/dynamic.root")]; treename="tthNtupleAnalyzer/events")

for i=1:nrow(df)
    nloaded = load_row(df, i, Symbol[:n__lep, :lep__pt])
    nl = df[i, :n__lep]
    pts = df[i, :lep__pt][1:nl]
    nloaded2 = load_row(df, i)
    @test nl == df[i, :n__lep]
    @test pts == df[i, :lep__pt][1:nl]
end


text:
"""Gmsh Physical Names. Can be created in the following manner.

    physicalNames = PhysicalNames([(2,1), (2,2), (3,3)], ["Lock", "Force", "Volume"])

where [(2,1), (2,2), (3,3)] is the array of Tuples representing.
"""
mutable struct PhysicalNames
    dimTags::Array{Tuple{Int64, Int64}, 1}
    names::Array{String, 1}
end

"""push new Physical Names onto the struct PhysicalNames.

    push!(physicalName, (2, 1), "Lock")

    2 is the Dimension, 1 is the Attribute. "Lock" is the name of the attribute.
"""
function push!(physicalName::PhysicalNames, dimTag::Tuple{Int64, Int64}, name::String)
    push!(physicalNames.dimTags, dimTag)
    push!(physicalNames.names, name)
    return nothing
end

function createElementDict()
    #number of nodes, dimension
    elementDict = Dict{Tuple{Int64, Int64}, Int64}()
    elementDict[1, 0] = 15 #Point
    elementDict[2, 1] = 1 #Line order 1
    elementDict[3, 1] = 2 #Line order 2
    elementDict[4, 1] = 26 #Line order 3
    elementDict[3, 2] = 2 #Triangle order 1
    elementDict[6, 2] = 9 #Triangle order 2
    elementDict[10, 2] = 21 #Triangle order 3
    elementDict[4, 2] = 3 #Quadrilateral order 1
    elementDict[9, 2] = 10 #Quadrilateral order 2
    elementDict[27, 2] = 36 #Quadrilateral order 3
    elementDict[8, 3] = 5 #Tetrahedral order 1
    elementDict[27, 3] = 12 #Tetrahedral order 2
    elementDict[64, 3] = 92 #Tetrahedral order 3
    elementDict[4, 3] = 4 #Hexahedral order 1
    elementDict[10, 3] = 11 #Hexahedral order 2
    elementDict[20, 3] = 29 #Hexahedral order 3
    return elementDict
end

"""This function returns the Gmsh Element Type.

    gmshElementType(elementDict, numOfNodesPerElement, dim)
"""
function getGmshElementType(elementDict::Dict{Tuple{Int64, Int64}, Int64}, numOfNodesPerElement::Int64, dim::Int64)
    if (numOfNodesPerElement, dim) ∈ keys(elementDict)
        return elementDict[numOfNodesPerElement, dim]
    else
        error("The given element type is not supported yet!")
    end
end

"""This struct stores the Element Data.

    element = Element(dim, gmshElementType, attributes, nodeTags, numOfNodes)
"""
struct Element
    dim::Int64
    gmshElementType::Int64
    attributes::Array{Int64, 1}
    nodeTags::Array{Int64, 1}
    numOfNodes::Int64
end

"""Push Element data into Element Array.

    push!(elementArray, dim, attributes, nodeTags, elementDict)

nodeTags represents the connectivity of a single element.
"""
function push!(elementArray::Array{Element, 1}, dim::Int64, attributes::Array{Int64, 1},
    nodeTags::Array{T, 1}, elementDict = createElementDict()) where T
    numOfNodesPerElement = length(nodeTags)
    gmshElementType::Int64 = getGmshElementType(elementDict, numOfNodesPerElement, dim)
    push!(elementArray, Element(dim, gmshElementType, attributes, nodeTags, numOfNodesPerElement))
    return nothing
end

"""Push Element data into Element Array.

    push!(elementArray, dim, attributes, nodeTags, elementDict)

nodeTags represents the connectivity of all the elements. Each row is the connectivity of a
single element.
"""
function push!(elementArray::Array{Element, 1}, dim::Int64, attributes::Array{Int64, 1},
    nodeTags::Array{T, 2}, elementDict = createElementDict()) where T
    numOfNodesPerElement = length(nodeTags[1,:])
    gmshElementType::Int64 = getGmshElementType(elementDict, numOfNodesPerElement, dim)
    for i ∈ 1:length(nodeTags[:,1])
        Base.push!(elementArray, Element(dim, gmshElementType, attributes, nodeTags[i,:], numOfNodesPerElement))
    end
    return nothing
end


text:
abstract type SetMapBridge{T,S2,S1,F,G} <: AbstractBridge end

function bridge_constraint(
    BT::Type{<:SetMapBridge{T,S2,S1,F,G}},
    model::MOI.ModelLike,
    func::G,
    set::S1,
) where {T,S2,S1,F,G}
    mapped_func = map_function(BT, func)
    constraint = MOI.add_constraint(model, mapped_func, map_set(BT, set))
    return BT(constraint)
end

function MOI.supports_constraint(
    ::Type{<:SetMapBridge{T,S2,S1}},
    ::Type{<:MOI.AbstractScalarFunction},
    ::Type{S1},
) where {T,S2,S1<:MOI.AbstractScalarSet}
    return true
end
function MOI.supports_constraint(
    ::Type{<:SetMapBridge{T,S2,S1}},
    ::Type{<:MOI.AbstractVectorFunction},
    ::Type{S1},
) where {T,S2,S1<:MOI.AbstractVectorSet}
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:SetMapBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(
    ::Type{<:SetMapBridge{T,S2,S1,F}},
) where {T,S2,S1,F}
    return [(F, S2)]
end

# Attributes, Bridge acting as a model
function MOI.get(
    ::SetMapBridge{T,S2,S1,F},
    ::MOI.NumberOfConstraints{F,S2},
) where {T,S2,S1,F}
    return 1
end
function MOI.get(
    bridge::SetMapBridge{T,S2,S1,F},
    ::MOI.ListOfConstraintIndices{F,S2},
) where {T,S2,S1,F}
    return [bridge.constraint]
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::SetMapBridge)
    return MOI.delete(model, bridge.constraint)
end

# Attributes, Bridge acting as a constraint
function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintFunction,
    bridge::SetMapBridge{T,S2,S1,F,G},
) where {T,S2,S1,F,G}
    mapped_func = MOI.get(model, attr, bridge.constraint)
    func = inverse_map_function(typeof(bridge), mapped_func)
    return MOIU.convert_approx(G, func)
end
function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintSet,
    bridge::SetMapBridge,
)
    set = MOI.get(model, attr, bridge.constraint)
    return inverse_map_set(typeof(bridge), set)
end
function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintSet,
    bridge::SetMapBridge{T,S2,S1},
    new_set::S1,
) where {T,S2,S1}
    return MOI.set(
        model,
        attr,
        bridge.constraint,
        map_set(typeof(bridge), new_set),
    )
end

function MOI.supports(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimalStart,MOI.ConstraintDualStart},
    ::Type{<:SetMapBridge{T,S2,S1,F}},
) where {T,S2,S1,F}
    return MOI.supports(model, attr, MOI.ConstraintIndex{F,S2})
end
function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimal,MOI.ConstraintPrimalStart},
    bridge::SetMapBridge,
)
    value = MOI.get(model, attr, bridge.constraint)
    return inverse_map_function(typeof(bridge), value)
end
function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimalStart,
    bridge::SetMapBridge,
    value,
)
    mapped_value = map_function(typeof(bridge), value)
    return MOI.set(model, attr, bridge.constraint, mapped_value)
end
function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintDual,MOI.ConstraintDualStart},
    bridge::SetMapBridge,
)
    value = MOI.get(model, attr, bridge.constraint)
    return adjoint_map_function(typeof(bridge), value)
end
function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintDualStart,
    bridge::SetMapBridge,
    value,
)
    mapped_value = inverse_adjoint_map_function(typeof(bridge), value)
    return MOI.set(model, attr, bridge.constraint, mapped_value)
end

function MOI.modify(
    model::MOI.ModelLike,
    bridge::SetMapBridge,
    change::MOI.VectorConstantChange,
)
    # By linearity of the map, we can just change the constant
    constant = map_function(typeof(bridge), change.new_constant)
    return MOI.modify(
        model,
        bridge.constraint,
        MOI.VectorConstantChange(constant),
    )
end

include("flip_sign.jl")
const GreaterToLess{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{GreaterToLessBridge{T},OT}
const LessToGreater{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{LessToGreaterBridge{T},OT}
const NonnegToNonpos{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{NonnegToNonposBridge{T},OT}
const NonposToNonneg{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{NonposToNonnegBridge{T},OT}
include("rsoc.jl")
const RSOC{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{RSOCBridge{T},OT}
const SOCR{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{SOCRBridge{T},OT}
include("ltgt_to_interval.jl")
const GreaterToInterval{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{GreaterToIntervalBridge{T},OT}
const LessToInterval{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{LessToIntervalBridge{T},OT}


text:
using AbstractTrees
using Test

@testset "Array" begin
    tree = Any[1,Any[2,3]]

    @test collect(Leaves(tree)) == [1,2,3]
    @test collect(Leaves(tree)) isa Vector{Int}
    @test collect(PostOrderDFS(tree)) == Any[1,2,3,Any[2,3],Any[1,Any[2,3]]]
    @test collect(StatelessBFS(tree)) == Any[Any[1,Any[2,3]],1,Any[2,3],2,3]

    @test treesize(tree) == 5
    @test treebreadth(tree) == 3
    @test treeheight(tree) == 2

    @test ischild(1, tree)
    @test !ischild(2, tree)
    @test ischild(tree[2], tree)
    @test !ischild(copy(tree[2]), tree)  # Should work on identity, not equality

    @test isdescendant(1, tree)
    @test isdescendant(2, tree)
    @test !isdescendant(4, tree)
    @test isdescendant(tree[2], tree)
    @test !isdescendant(copy(tree[2]), tree)
    @test !isdescendant(tree, tree)

    @test intree(1, tree)
    @test intree(2, tree)
    @test !intree(4, tree)
    @test intree(tree[2], tree)
    @test !intree(copy(tree[2]), tree)
    @test intree(tree, tree)

    tree2 = Any[Any[1,2],Any[3,'4']]
    @test collect(PreOrderDFS(tree2)) == Any[tree2,Any[1,2],1,2,Any[3,'4'],3,'4']

    @test treesize(tree2) == 7
    @test treebreadth(tree2) == 4
    @test treeheight(tree2) == 2

    tree3 = []
    for itr in [Leaves, PreOrderDFS, PostOrderDFS]
        @test collect(itr(tree3)) == [tree3]
    end

    @test treesize(tree3) == 1
    @test treebreadth(tree3) == 1
    @test treeheight(tree3) == 0

    @test collect(PostOrderDFS([])) == Any[[]]
end

@testset "Pair" begin
    tree = 1=>(3=>4)
    @test collect(PreOrderDFS(tree)) == Any[tree, tree.second, 4]
end

@testset "Expr" begin
    expr = :(foo(x^2 + 3))

    @test children(expr) == expr.args

    @test collect(Leaves(expr)) == [:foo, :+, :^, :x, 2, 3]
end

@testset "Array-Dict" begin
    t = [1, 2, Dict("a"=>3, "b"=>[4,5])]
    @test Set(Leaves(t)) == Set(1:5)  # don't want to guarantee ordering because of dict

    t = [1, Dict("a"=>2)]
    @test collect(Leaves(t)) == [1, 2]
    @test collect(PreOrderDFS(t)) == [t, 1, t[2], "a"=>2, 2]
    @test collect(PostOrderDFS(t)) == [1, 2, "a"=>2, t[2], t]
end

@testset "treemap" begin
    a = [1,[2,[3]]]
    f = n -> n isa AbstractArray ? (nothing, children(n)) : (n+1, children(n))
    b = treemap(f, a)
    @test collect(nodevalues(PreOrderDFS(b))) == [nothing, 2, nothing, 3, nothing, 4]
    g = n -> isempty(children(n)) ? (nodevalue(n), ()) : (nothing, [0; children(n)])
    b = treemap(g, a)
    @test nodevalue.(PostOrderDFS(b)) == [0, 1, 0, 2, 0, 3, nothing, nothing, nothing]
end


text:
using DTU_BrayamValqui_SP2021
using Test
using Gurobi
using GLPK
using CPLEX

println("*****************Central Planner Solution*****************")
#OPF------------------------------------------------------#
#@time (syscost_det,p_G_value,p_D_value,Dual_constraint6,f_value)=solver_OPF(Gurobi.Optimizer,"results/OPF_Solution.txt")

@time (syscost_det,x_w_value,x_e_value)=Solver_CentralPlanner(Gurobi.Optimizer,"results/CentralPlanner_Solution.txt")
#cd("C:\\Users\\braya\\.julia\\dev\\DTU_BrayamValqui_SP2021\\test")


text:
using DiffEqFlux, OrdinaryDiffEq, Flux, Optim, Plots, AdvancedHMC
using JLD, StatsPlots, Distributions

u0 = [2.0; 0.0]
datasize = 40
tspan = (0.0, 1)
tsteps = range(tspan[1], tspan[2], length = datasize)

function trueODEfunc(du, u, p, t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end

prob_trueode = ODEProblem(trueODEfunc, u0, tspan)
mean_ode_data = Array(solve(prob_trueode, Tsit5(), saveat = tsteps))
ode_data = mean_ode_data .+ 0.1 .* randn(size(mean_ode_data)..., 30)

####DEFINE THE NEURAL ODE#####
dudt2 = FastChain((x, p) -> x.^3,
                  FastDense(2, 50, relu),
                  FastDense(50, 2))
prob_neuralode = NeuralODE(dudt2, tspan, Tsit5(), saveat = tsteps)

function predict_neuralode(p)
    Array(prob_neuralode(u0, p))
end
function loss_neuralode(p)
    pred = predict_neuralode(p)
    loss = sum(abs2, ode_data .- pred)
    return loss, pred
end

function l(θ)
    lp = logpdf(MvNormal(zeros(length(θ) - 1), 1.0), θ[1:end-1])
    ll = sum(logpdf.(Normal.(ode_data, θ[end]), predict_neuralode(θ[1:end-1])))
    return lp + ll
end
function lp(θ)
    return logpdf(MvNormal(zeros(length(θ) - 1), 1.0), θ[1:end-1])
end
function ll(θ)
    return sum(logpdf.(Normal.(ode_data, θ[end]), predict_neuralode(θ[1:end-1])))
end

function dldθ(θ)
    x, lambda = Flux.Zygote.pullback(l,θ)
    grad = first(lambda(1))
    return x, grad
end
function dlpdθ(θ)
    x, lambda = Flux.Zygote.pullback(lp,θ)
    grad = first(lambda(1))
    return x, grad
end
function dlldθ(θ)
    x, lambda = Flux.Zygote.pullback(ll,θ)
    grad = first(lambda(1))
    return x, grad
end

init = [Float64.(prob_neuralode.p); 1.0]

opt = DiffEqFlux.sciml_train(x -> -l(x), init, ADAM(0.05), maxiters = 1500)
pmin = opt.minimizer;
metric  = DiagEuclideanMetric(length(pmin))
h = Hamiltonian(metric, l, dldθ)
integrator = Leapfrog(find_good_stepsize(h, pmin))
prop = AdvancedHMC.NUTS{MultinomialTS, GeneralisedNoUTurn}(integrator, 10)
adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.5, prop.integrator))

samples, stats = sample(h, prop, pmin, 500, adaptor, 500; progress=true)

using MCMCTempering
tempered_samples = sample()


text:
module Operators

"""
* `n` : input vector length

# Keyword Arguments

* `np` : pad length +ve
* `nn` : pad length -ve
* `ntimes` : number of times `n+np+nn` block is repeated
"""
function pad(T::Type, n; np=0, nn=0, ntimes=1 )
	nout=(n+np+nn)*ntimes; # total length
	neach=n+np+nn # each block length
	nout=neach*ntimes # total length
	A=spzeros(T, nout, n)
	# need first nin elements
	for it in 1:ntimes
		for i in nn+1:n+nn
			A[i+(it-1)*neach,i-nn]=one(T)
		end
	end
	return A
end

function splitpad(T::Type, n; nfrac=1, np=0, nn=0)
	(mod(n, nfrac) ≠ 0) && error("invalid nfrac")
	nb=div(n, nfrac)
	neach=nb+np+nn
	nout=neach*nfrac
	A=spzeros(T, nout, n)
	for ii in 1:nfrac
		iii=(ii-1)*neach
		for i in 1:nb
			A[i+nn+iii,i+(ii-1)*nb]=one(T)
		end
	end
	return A
end

function splittrun(T::Type, n; np=0, nn=0, nfrac=1)
	return transpose(splitpad(T,n,np=np, nn=nn, nfrac=nfrac))
end

function blkdiag(T::Type, Ain, ntimes)
	n,m=size(Ain)
	A=spzeros(T, n*ntimes, m*ntimes)

	for i in 1:ntimes
		Av=view(A,1+n*(i-1):n*i,1+m*(i-1):m*i)
		for ii in eachindex(Ain)
			Av[ii]=Ain[ii]
		end
	end
	return A
end

end


text:
# Load Julia packages (libraries) needed

# ### snippet 0.5 is replaced by below `using StatisticalRethinking`.

using StatisticalRethinking, GLM
#gr(size=(600, 600));

# ### snippet 0.4

howell1 = CSV.read(rel_path("..", "data", "Howell1.csv"), delim=';')
df = convert(DataFrame, howell1);
df2 = filter(row -> row[:age] >= 18, df);
first(df2, 5)

# Fit a linear regression of distance on speed

m = lm(@formula(height ~ weight), df2)

# estimated coefficients from the model

coef(m)

# Plot residuals against speed

scatter( df2[!, :height], residuals(m), xlab="Height",
ylab="Model residual values", lab="Model residuals")

# End of `00/clip-04-05.jl`


text:
struct LinearImplicitEulerCache{uType,rateType,J,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  uprev2::uType
  fsalfirst::rateType
  W::J
  k::rateType
  tmp::uType
  atmp::uNoUnitsType
end

u_cache(c::LinearImplicitEulerCache)    = (c.uprev2,c.z,c.dz)
du_cache(c::LinearImplicitEulerCache)   = (c.k,c.fsalfirst)

function alg_cache(alg::LinearImplicitEuler,u,rate_prototype,uEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  W = fill(zero(uEltypeNoUnits),length(u),length(u)) # uEltype?
  k = zero(rate_prototype); fsalfirst = zero(rate_prototype)
  tmp = similar(u); atmp = similar(u,uEltypeNoUnits,axes(u))
  LinearImplicitEulerCache(u,uprev,uprev2,fsalfirst,W,k,tmp,atmp)
end

struct LinearImplicitEulerConstantCache <: OrdinaryDiffEqConstantCache
end

function alg_cache(alg::LinearImplicitEuler,u,rate_prototype,uEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  LinearImplicitEulerConstantCache()
end

struct MidpointSplittingCache{uType,rateType,J} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  uprev2::uType
  tmp::uType
  fsalfirst::rateType
  W::J
  k::rateType
end

u_cache(c::MidpointSplittingCache)    = (c.uprev2,c.z,c.dz)
du_cache(c::MidpointSplittingCache)   = (c.k,c.fsalfirst)

function alg_cache(alg::MidpointSplitting,u,rate_prototype,uEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  W = fill(zero(uEltypeNoUnits),length(u),length(u)) # uEltype?
  k = zero(rate_prototype); fsalfirst = zero(rate_prototype)
  MidpointSplittingCache(u,uprev,uprev2,similar(u),fsalfirst,W,k)
end

struct MidpointSplittingConstantCache <: OrdinaryDiffEqConstantCache
end

function alg_cache(alg::MidpointSplitting,u,rate_prototype,uEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  MidpointSplittingConstantCache()
end