text:
# 4D Vertex structure

struct Vertex4D
    dx::Float32
    dy::Float32
    dz::Float32
    dw::Float32
    xsvp::Int64
    ysvp::Int64
    zsvp::Int64
    wsvp::Int64

    function Vertex4D(xsv::Int, ysv::Int, zsv::Int, wsv::Int)
        ssv = Float32((xsv + ysv + zsv + wsv) * UNSKEW_4D)
        new(-xsv - ssv, -ysv - ssv, -zsv - ssv, -wsv - ssv,
            (xsv * PRIME_X)%Int64,
            (ysv * PRIME_Y)%Int64,
            (zsv * PRIME_Z)%Int64,
            (wsv * PRIME_W)%Int64)
    end
end
Vertex4D(v) = Vertex4D(((v >> 0) & 3) - 1, ((v >> 2) & 3) - 1,
                       ((v >> 4) & 3) - 1, ((v >> 6) & 3) - 1)

# Lookup Tables & Gradients

const vertex_ind = [
     7,17,22,24,25,26,28,29,35,36,39,40,53,54,57,58,64,65,67,68,
     7,17,22,25,26,28,29,35,36,40,53,54,58,65,68,
     1, 3, 5, 7,15,17,22,25,26,29,36,40,54,58,65,68,
     1, 7, 8,17,18,22,23,25,26,29,36,40,54,58,65,68,69,
     7,17,24,25,26,28,29,35,39,40,53,57,58,67,68,
     3, 7,17,25,26,28,29,35,36,39,40,53,54,57,58,68,
     3, 7,17,25,26,28,29,36,40,54,58,68,
     3, 7, 8,17,18,25,26,28,29,36,40,54,58,68,69,
     2, 3, 6, 7,16,17,24,25,28,29,39,40,57,58,67,68,
     3, 7,17,25,26,28,29,39,40,57,58,68,
     3, 7,17,25,26,28,29,40,58,68,
     3, 7, 8,17,18,25,26,28,29,30,40,58,68,69,
     2, 7, 9,17,19,24,25,27,28,29,39,40,57,58,67,68,71,
     3, 7, 9,17,19,25,26,28,29,39,40,57,58,68,71,
     3, 7, 9,17,19,25,26,28,29,31,40,58,68,71,
     3, 7,10,17,20,25,26,28,29,30,31,40,58,68,69,71,72,
     7,22,24,25,26,28,35,36,39,40,53,64,65,67,68,
     5, 7,22,25,26,28,29,35,36,39,40,53,54,64,65,68,
     5, 7,22,25,26,29,35,36,40,54,65,68,
     5, 7, 8,22,23,25,26,29,35,36,40,54,65,68,69,
     6, 7,24,25,26,28,29,35,36,39,40,53,57,64,67,68,
     7,25,26,28,29,35,36,39,40,53,58,65,67,68,
     7,25,26,28,29,35,36,39,40,54,58,65,68,69,
     7, 8,25,26,29,36,40,41,54,58,65,68,69,
     6, 7,24,25,28,29,35,39,40,57,67,68,
     7,25,26,28,29,35,36,39,40,57,58,67,68,71,
     7,25,26,28,29,35,36,39,40,58,68,
     7, 8,25,26,28,29,36,40,41,58,68,69,
     6, 7, 9,24,25,27,28,29,35,39,40,57,67,68,71,
     7, 9,25,28,29,39,40,43,57,58,67,68,71,
     7, 9,25,26,28,29,39,40,43,58,68,71,
     7,10,25,26,28,29,40,41,43,58,68,69,71,72,
     4, 5, 6, 7,21,22,24,25,35,36,39,40,64,65,67,68,
     5, 7,22,25,26,35,36,39,40,64,65,68,
     5, 7,22,25,26,35,36,40,65,68,
     5, 7, 8,22,23,25,26,35,36,37,40,65,68,69,
     6, 7,24,25,28,35,36,39,40,64,67,68,
     7,25,26,28,29,35,36,39,40,64,65,67,68,76,
     7,25,26,28,29,35,36,39,40,65,68,
     7, 8,25,26,29,35,36,40,41,65,68,69,
     6, 7,24,25,28,35,39,40,67,68,
     7,25,26,28,29,35,36,39,40,67,68,
     7,25,26,28,29,35,36,39,40,68,
     7, 8,25,26,28,29,35,36,39,40,41,68,69,
     6, 7, 9,24,25,27,28,35,39,40,42,67,68,71,
     7, 9,25,28,29,35,39,40,43,67,68,71,
     7, 9,25,26,28,29,35,36,39,40,43,68,71,
     7,25,26,28,29,36,39,40,41,43,58,68,69,71,72,
     4, 7,11,22,24,25,32,34,35,36,39,40,64,65,67,68,76,
     5, 7,11,22,25,26,32,35,36,39,40,64,65,68,76,
     5, 7,11,22,25,26,32,35,36,40,44,65,68,76,
     5, 7,12,22,25,26,33,35,36,37,40,44,65,68,69,76,77,
     6, 7,11,24,25,28,34,35,36,39,40,64,67,68,76,
     7,11,25,35,36,39,40,46,64,65,67,68,76,
     7,11,25,26,35,36,39,40,46,65,68,76,
     7,12,25,26,35,36,40,41,46,65,68,69,76,77,
     6, 7,11,24,25,28,34,35,39,40,45,67,68,76,
     7,11,25,28,35,36,39,40,46,67,68,76,
     7,11,25,26,28,29,35,36,39,40,46,68,76,
     7,25,26,29,35,36,39,40,41,46,65,68,69,76,77,
     6, 7,13,24,25,28,35,38,39,40,42,45,67,68,71,76,78,
     7,13,25,28,35,39,40,43,46,67,68,71,76,78,
     7,25,28,29,35,36,39,40,43,46,67,68,71,76,78,
     7,25,26,28,29,35,36,39,40,41,43,46,68,69,71,76,79,
    17,22,24,25,26,28,35,53,54,57,58,64,65,67,68,
    15,17,22,25,26,28,29,35,36,53,54,57,58,64,65,68,
    15,17,22,25,26,29,36,53,54,58,65,68,
    15,17,18,22,23,25,26,29,36,53,54,58,65,68,69,
    16,17,24,25,26,28,29,35,39,53,54,57,58,64,67,68,
    17,25,26,28,29,35,40,53,54,57,58,65,67,68,
    17,25,26,28,29,36,40,53,54,57,58,65,68,69,
    17,18,25,26,29,36,40,54,58,59,65,68,69,
    16,17,24,25,28,29,39,53,57,58,67,68,
    17,25,26,28,29,39,40,53,54,57,58,67,68,71,
    17,25,26,28,29,40,53,54,57,58,68,
    17,18,25,26,28,29,40,54,58,59,68,69,
    16,17,19,24,25,27,28,29,39,53,57,58,67,68,71,
    17,19,25,28,29,39,40,57,58,61,67,68,71,
    17,19,25,26,28,29,40,57,58,61,68,71,
    17,20,25,26,28,29,40,58,59,61,68,69,71,72,
    21,22,24,25,26,28,35,36,39,53,54,57,64,65,67,68,
    22,25,26,28,35,36,40,53,54,58,64,65,67,68,
    22,25,26,29,35,36,40,53,54,58,64,65,68,69,
    22,23,25,26,29,36,40,54,58,65,66,68,69,
    24,25,26,28,35,39,40,53,57,58,64,65,67,68,
    25,26,28,29,35,36,39,40,53,54,57,58,64,65,67,68,
     7,17,22,25,26,28,29,35,36,40,53,54,58,65,68,69,
    25,26,29,36,40,54,58,65,68,69,
    24,25,28,29,35,39,40,53,57,58,64,67,68,71,
     7,17,24,25,26,28,29,35,39,40,53,57,58,67,68,71,
     7,17,25,26,28,29,35,36,39,40,53,54,57,58,65,67,68,69,71,
    25,26,28,29,36,40,54,58,65,68,69,
    24,25,27,28,29,39,40,57,58,67,68,70,71,
    25,28,29,39,40,57,58,67,68,71,
    25,26,28,29,39,40,57,58,67,68,71,
    25,26,28,29,40,58,68,69,71,72,
    21,22,24,25,35,36,39,53,64,65,67,68,
    22,25,26,35,36,39,40,53,54,64,65,67,68,76,
    22,25,26,35,36,40,53,54,64,65,68,
    22,23,25,26,35,36,40,54,65,66,68,69,
    24,25,28,35,36,39,40,53,57,64,65,67,68,76,
     7,22,24,25,26,28,35,36,39,40,53,64,65,67,68,76,
     7,22,25,26,28,29,35,36,39,40,53,54,58,64,65,67,68,69,76,
    25,26,29,35,36,40,54,58,65,68,69,
    24,25,28,35,39,40,53,57,64,67,68,
     7,24,25,26,28,29,35,36,39,40,53,57,58,64,65,67,68,71,76,
     7,25,26,28,29,35,36,39,40,58,65,67,68,
     7,25,26,28,29,35,36,39,40,54,58,65,68,69,
    24,25,27,28,35,39,40,57,67,68,70,71,
    25,28,29,35,39,40,57,58,67,68,71,
     7,25,26,28,29,35,36,39,40,57,58,67,68,71,
     7,25,26,28,29,36,39,40,58,68,69,71,72,
    21,22,24,25,32,34,35,36,39,53,64,65,67,68,76,
    22,25,32,35,36,39,40,64,65,67,68,74,76,
    22,25,26,32,35,36,40,64,65,68,74,76,
    22,25,26,33,35,36,40,65,66,68,69,74,76,77,
    24,25,34,35,36,39,40,64,65,67,68,75,76,
    25,35,36,39,40,64,65,67,68,76,
    25,26,35,36,39,40,64,65,67,68,76,
    25,26,35,36,40,65,68,69,76,77,
    24,25,28,34,35,39,40,64,67,68,75,76,
    25,28,35,36,39,40,64,65,67,68,76,
     7,25,26,28,29,35,36,39,40,64,65,67,68,76,
     7,25,26,29,35,36,39,40,65,68,69,76,77,
    24,25,28,35,38,39,40,67,68,70,71,75,76,78,
    25,28,35,39,40,67,68,71,76,78,
     7,25,28,29,35,36,39,40,67,68,71,76,78,
    25,26,28,29,35,36,39,40,68,69,71,76,79,
    14,15,16,17,21,22,24,25,53,54,57,58,64,65,67,68,
    15,17,22,25,26,53,54,57,58,64,65,68,
    15,17,22,25,26,53,54,58,65,68,
    15,17,18,22,23,25,26,53,54,55,58,65,68,69,
    16,17,24,25,28,53,54,57,58,64,67,68,
    17,25,26,28,29,53,54,57,58,64,65,67,68,86,
    17,25,26,28,29,53,54,57,58,65,68,
    17,18,25,26,29,53,54,58,59,65,68,69,
    16,17,24,25,28,53,57,58,67,68,
    17,25,26,28,29,53,54,57,58,67,68,
    17,25,26,28,29,53,54,57,58,68,
    17,18,25,26,28,29,53,54,57,58,59,68,69,
    16,17,19,24,25,27,28,53,57,58,60,67,68,71,
    17,19,25,28,29,53,57,58,61,67,68,71,
    17,19,25,26,28,29,53,54,57,58,61,68,71,
    17,25,26,28,29,40,54,57,58,59,61,68,69,71,72,
    21,22,24,25,35,53,54,57,64,65,67,68,
    22,25,26,35,36,53,54,57,58,64,65,67,68,86,
    22,25,26,35,36,53,54,58,64,65,68,
    22,23,25,26,36,53,54,58,65,66,68,69,
    24,25,28,35,39,53,54,57,58,64,65,67,68,86,
    17,22,24,25,26,28,35,53,54,57,58,64,65,67,68,86,
    17,22,25,26,28,29,35,36,40,53,54,57,58,64,65,67,68,69,86,
    25,26,29,36,40,53,54,58,65,68,69,
    24,25,28,35,39,53,57,58,64,67,68,
    17,24,25,26,28,29,35,39,40,53,54,57,58,64,65,67,68,71,86,
    17,25,26,28,29,40,53,54,57,58,65,67,68,
    17,25,26,28,29,36,40,53,54,57,58,65,68,69,
    24,25,27,28,39,53,57,58,67,68,70,71,
    25,28,29,39,40,53,57,58,67,68,71,
    17,25,26,28,29,39,40,53,54,57,58,67,68,71,
    17,25,26,28,29,40,54,57,58,68,69,71,72,
    21,22,24,25,35,53,64,65,67,68,
    22,25,26,35,36,53,54,64,65,67,68,
    22,25,26,35,36,53,54,64,65,68,
    22,23,25,26,35,36,53,54,64,65,66,68,69,
    24,25,28,35,39,53,57,64,65,67,68,
    22,24,25,26,28,35,36,39,40,53,54,57,58,64,65,67,68,76,86,
    22,25,26,35,36,40,53,54,58,64,65,67,68,
    22,25,26,29,35,36,40,53,54,58,64,65,68,69,
    24,25,28,35,39,53,57,64,67,68,
    24,25,28,35,39,40,53,57,58,64,65,67,68,
    25,26,28,29,35,36,39,40,53,54,57,58,64,65,67,68,
    25,26,28,29,35,36,40,53,54,58,65,67,68,69,
    24,25,27,28,35,39,53,57,64,67,68,70,71,
    24,25,28,29,35,39,40,53,57,58,64,67,68,71,
    25,26,28,29,35,39,40,53,57,58,65,67,68,71,
    25,26,28,29,36,39,40,54,57,58,65,67,68,69,71,72,
    21,22,24,25,32,34,35,53,64,65,67,68,73,76,
    22,25,32,35,36,53,64,65,67,68,74,76,
    22,25,26,32,35,36,53,54,64,65,68,74,76,
    22,25,26,35,36,40,54,64,65,66,68,69,74,76,77,
    24,25,34,35,39,53,64,65,67,68,75,76,
    25,35,36,39,40,53,64,65,67,68,76,
    22,25,26,35,36,39,40,53,54,64,65,67,68,76,
    22,25,26,35,36,40,54,64,65,68,69,76,77,
    24,25,28,34,35,39,53,57,64,67,68,75,76,
    24,25,28,35,36,39,40,53,57,64,65,67,68,76,
    25,26,28,35,36,39,40,53,58,64,65,67,68,76,
    25,26,29,35,36,39,40,54,58,64,65,67,68,69,76,77,
    24,25,28,35,39,40,57,64,67,68,70,71,75,76,78,
    24,25,28,35,39,40,57,64,67,68,71,76,78,
    25,28,29,35,36,39,40,57,58,64,65,67,68,71,76,78,
    25,26,28,29,35,36,39,40,58,65,67,68,69,71,76,
    14,17,22,24,25,47,50,52,53,54,57,58,64,65,67,68,86,
    15,17,22,25,26,47,50,53,54,57,58,64,65,68,86,
    15,17,22,25,26,47,50,53,54,58,65,68,80,86,
    15,17,22,25,26,48,51,53,54,55,58,65,68,69,80,86,87,
    16,17,24,25,28,47,52,53,54,57,58,64,67,68,86,
    17,25,47,53,54,57,58,64,65,67,68,82,86,
    17,25,26,47,53,54,57,58,65,68,82,86,
    17,25,26,48,53,54,58,59,65,68,69,82,86,87,
    16,17,24,25,28,47,52,53,57,58,67,68,81,86,
    17,25,28,47,53,54,57,58,67,68,82,86,
    17,25,26,28,29,47,53,54,57,58,68,82,86,
    17,25,26,29,53,54,57,58,59,65,68,69,82,86,87,
    16,17,24,25,28,49,53,56,57,58,60,67,68,71,81,86,88,
    17,25,28,49,53,57,58,61,67,68,71,82,86,88,
    17,25,28,29,53,54,57,58,61,67,68,71,82,86,88,
    17,25,26,28,29,53,54,57,58,59,61,68,69,71,82,86,89,
    21,22,24,25,35,50,52,53,54,57,64,65,67,68,86,
    22,25,50,53,54,57,58,64,65,67,68,84,86,
    22,25,26,50,53,54,58,64,65,68,84,86,
    22,25,26,51,53,54,58,65,66,68,69,84,86,87,
    24,25,52,53,54,57,58,64,65,67,68,85,86,
    25,53,54,57,58,64,65,67,68,86,
    25,26,53,54,57,58,64,65,67,68,86,
    25,26,53,54,58,65,68,69,86,87,
    24,25,28,52,53,57,58,64,67,68,85,86,
    25,28,53,54,57,58,64,65,67,68,86,
    17,25,26,28,29,53,54,57,58,64,65,67,68,86,
    17,25,26,29,53,54,57,58,65,68,69,86,87,
    24,25,28,53,56,57,58,67,68,70,71,85,86,88,
    25,28,53,57,58,67,68,71,86,88,
    17,25,28,29,53,54,57,58,67,68,71,86,88,
    25,26,28,29,53,54,57,58,68,69,71,86,89,
    21,22,24,25,35,50,52,53,64,65,67,68,83,86,
    22,25,35,50,53,54,64,65,67,68,84,86,
    22,25,26,35,36,50,53,54,64,65,68,84,86,
    22,25,26,36,53,54,58,64,65,66,68,69,84,86,87,
    24,25,35,52,53,57,64,65,67,68,85,86,
    25,35,53,54,57,58,64,65,67,68,86,
    22,25,26,35,36,53,54,57,58,64,65,67,68,86,
    22,25,26,36,53,54,58,64,65,68,69,86,87,
    24,25,28,35,39,52,53,57,64,67,68,85,86,
    24,25,28,35,39,53,54,57,58,64,65,67,68,86,
    25,26,28,35,40,53,54,57,58,64,65,67,68,86,
    25,26,29,36,40,53,54,57,58,64,65,67,68,69,86,87,
    24,25,28,39,53,57,58,64,67,68,70,71,85,86,88,
    24,25,28,39,53,57,58,64,67,68,71,86,88,
    25,28,29,39,40,53,54,57,58,64,65,67,68,71,86,88,
    25,26,28,29,40,53,54,57,58,65,67,68,69,71,86,
    21,22,24,25,35,53,62,63,64,65,67,68,73,76,83,86,90,
    22,25,35,53,62,64,65,67,68,74,76,84,86,90,
    22,25,35,36,53,54,64,65,67,68,74,76,84,86,90,
    22,25,26,35,36,53,54,64,65,66,68,69,74,76,84,86,91,
    24,25,35,53,63,64,65,67,68,75,76,85,86,90,
    25,35,53,64,65,67,68,76,86,90,
    22,25,35,36,53,54,64,65,67,68,76,86,90,
    25,26,35,36,53,54,64,65,68,69,76,86,91,
    24,25,35,39,53,57,64,65,67,68,75,76,85,86,90,
    24,25,35,39,53,57,64,65,67,68,76,86,90,
    25,35,36,39,40,53,54,57,58,64,65,67,68,76,86,90,
    25,26,35,36,40,53,54,58,64,65,67,68,69,76,86,
    24,25,28,35,39,53,57,64,67,68,70,71,75,76,85,86,92,
    25,28,35,39,53,57,64,67,68,71,76,86,92,
    25,28,35,39,40,53,57,58,64,65,67,68,71,76,86,
    25,26,28,29,35,36,39,40,53,54,57,58,64,65,67,68,69,71,76,86,
]

const offsets_4D =
    UInt16[  20,  35,  51,  68,  83,  99, 111, 126, 142, 154, 164, 178, 195, 210, 224, 241,
            256, 272, 284, 299, 315, 329, 343, 356, 368, 382, 393, 405, 420, 433, 445, 459,
            475, 487, 497, 511, 523, 537, 548, 560, 570, 581, 591, 604, 618, 630, 643, 658,
            675, 690, 704, 721, 736, 749, 761, 775, 789, 801, 814, 829, 846, 860, 875, 892,
            907, 923, 935, 950, 966, 980, 994,1007,1019,1033,1044,1056,1071,1084,1096,1110,
           1126,1140,1154,1167,1181,1197,1213,1223,1237,1253,1272,1283,1296,1306,1317,1327,
           1339,1353,1364,1376,1390,1406,1425,1436,1447,1466,1479,1493,1505,1516,1530,1543,
           1558,1571,1583,1597,1610,1620,1631,1641,1653,1664,1678,1691,1705,1715,1728,1741,
           1757,1769,1779,1793,1805,1819,1830,1842,1852,1863,1873,1886,1900,1912,1925,1940,
           1952,1966,1977,1989,2003,2019,2038,2049,2060,2079,2092,2106,2118,2129,2143,2156,
           2166,2177,2187,2200,2211,2230,2243,2257,2267,2280,2296,2310,2323,2337,2351,2367,
           2381,2393,2406,2421,2433,2444,2458,2471,2484,2498,2512,2528,2543,2556,2572,2587,
           2604,2619,2633,2650,2665,2678,2690,2704,2718,2730,2743,2758,2775,2789,2804,2821,
           2836,2849,2861,2875,2888,2898,2909,2919,2931,2942,2956,2969,2983,2993,3006,3019,
           3033,3045,3058,3073,3085,3096,3110,3123,3136,3150,3164,3180,3195,3208,3224,3239,
           3256,3270,3285,3302,3316,3326,3339,3352,3367,3380,3396,3411,3428,3441,3456,3476]

const lookup_4D =
    Vertex4D.([0x01,0x04,0x05,0x10,0x11,0x14,0x15,0x16,0x19,0x1a,0x25,0x26,0x29,0x40,0x41,0x44,
               0x45,0x46,0x49,0x4a,0x50,0x51,0x52,0x54,0x55,0x56,0x58,0x59,0x5a,0x5b,0x5e,0x61,
               0x62,0x64,0x65,0x66,0x67,0x68,0x69,0x6a,0x6b,0x6d,0x6e,0x76,0x79,0x7a,0x85,0x86,
               0x89,0x91,0x92,0x94,0x95,0x96,0x97,0x98,0x99,0x9a,0x9b,0x9d,0x9e,0xa1,0xa4,0xa5,
               0xa6,0xa7,0xa9,0xaa,0xab,0xad,0xae,0xaf,0xb5,0xb6,0xb9,0xba,0xbb,0xbe,0xbf,0xd6,
               0xd9,0xda,0xe5,0xe6,0xe9,0xea,0xeb,0xee,0xef,0xfa,0xfb,0xfe])


text:
"""
    MTLDS_variational(mt_enc, mt_post, hphi, emission, d, flag_mt_emission,
    flag_mt_x0)

A deterministic MT-LDS model (subtype of `MTDSModel`) following the equations:

`` x_t = A_z x_{t-1} + B_z u_{t} ``

`` y_t = C_z x_y + d_z + s .* ϵ_t``

for ``ϵ_t \\sim 𝑵(0, I)`` and ``x_0 = 0``. Note the following restrictions vs a
full LDS:

1. Latent system is deterministic.
2. The inputs cannot affect the emission process (i.e. no ``D_z u_t`` term).
3. The latent state ``x_0`` is constant and equal to zero.
4. The emission covariance is diagonal.

The limitations exist because this was all I needed for my experiments, but it
is not difficult to extend the model to relax each of these, especially (2-4).

## Struct description
The model is variational in ``z``, and contains an encoder (`mt_enc`) which is
an RNN which encodes the sequences `Y` and `U` (which should be of size ``d × T ×
n_b``, where ``d`` is the size of each vector, ``T`` is the sequence length, and
``n_b`` is the batch size) into a single vector. One can instead perform standard
(non-amortized VI by using a LookupTable here. The `mt_post` network is a
network which returns a tuple of (mean, (diagonal) sqrt variance) defining a
variational Gaussian posterior. The method `encode` will return two tuples, the
first of which is a sample from this posterior (using reparameterization),
elements 2 and 3 are the variational posterior mean and stdev components. (The
second tuple is empty -- this is reserved for ``x_0`` if it has a separate
encoder.) Note that these variational networks need not be used, e.g. if one is
using an alternative MCO approach.

The function `hphi` is a network ``h_ϕ: \\text{dom}(z) → \\text{dom}(θ)`` which
takes a ``k × n_b`` matrix ``z`` and returns the LDS parameter vector θ. The
`emission` is either a pre-specified network which is optimized without any MT
modulation, or an Int which specifies the emission dimension, for which a MT
affine emission is constructed (``C_z + d_z``). Finally ``flag_mt_emission`` is
a boolean which is true in the latter case (MT affine emission) and false
otherwise, and `flag_mt_x0` is a placeholder for whether the ``x_0`` is
modulated via ``z``.

## Methods
As per all `MTDSModel` objects, the following interface is available:

- `encode(m::MTLDS_variational, Y::AbstractArray{T}, U::AbstractArray{T})`,
  which encodes the sequences `Y` and `U` into a sample and variational
  posterior parameters of ``z``.
- `forward(m::MTLDS_variational, z::AbstractVecOrMat{T}, U::AbstractArray{T})`,
  the forward model, which constructs the LDS for the latent variable `z` and
  calculates the output mean.
- `reconstruct(m::MTDSModel, Y, U)` which performs both `encode` and `forward`
  operations, effectively reconstructing `Y` through the bottleneck of the
  inference of ``z``.

### Constructor
See `create_mtlds` for a more helpful constructor.

"""
struct MTLDS_variational{U,V,W,S,T} <: MTDSModel
    mt_enc::U
    mt_post::V
    hphi::W
    emission::S
    logstd::T
    d::Int
    flag_mt_emission::Bool
    flag_mt_x0::Bool
end

Flux.@treelike MTLDS_variational
# common superclass properties
is_amortized(m::MTLDS_variational) = true
is_amortized(m::MTLDS_variational{U,V,W,S,T}) where {U <: LookupTable, V,W,S,T} = false
has_x0_encoding(m::MTLDS_variational) = false

function Base.show(io::IO, l::MTLDS_variational)
    d_x0, d_mt = l.d, size(l.mt_post.Dense1.W, 1)
    d_y = get_final_layer_dim(l.emission)
    if (d_y === nothing) || l.mt_enc isa LookupTable
        d_uy = d_y === nothing ? "" : format(", in=?, out={:d}", d_y)
    else
        d_u = size(l.mt_enc.cell.Wi, 2) - something(d_y, 0)
        d_uy = format(", in={:d}, out={:d}", d_u, d_y)
    end
    emission_type = l.flag_mt_emission ? "" : ", *non-MT emission function*"
    h0_type = l.flag_mt_x0 ? "'MT'" : "'non-MT'"
    enc_type = l.mt_enc isa LookupTable ? "'Standard VI'" : "'Amortized VI'"
    print(io, "MTLDS_variational(state=", d_x0, d_uy, ", d_mt=", d_mt, ", ", emission_type,
    "x0_type=", h0_type,  ", enc_type=", enc_type,")")
end

is_amortized(m::MTLDS_variational{U,V,W,S,T}) where {U <: LookupTable, V, W, S, T} = false

forward(m::MTLDS_variational, z, U; T_steps=size(U, 2)) = forward(m, nothing, z, U; T_steps=T_steps)
function forward(m::MTLDS_variational, x0, z::AbstractVecOrMat{T}, U::AbstractArray{T}; T_steps=size(Y,2)) where T <: AbstractFloat
    # note that T_steps is currently ignored (this is no great issue).
    !(x0 === nothing) && @warn "MTLDS_variational currently doesn't accept x0 -- fixed at 0. Ignoring..."
    @argcheck size(z,2) == size(U,3)
    θ = m.hphi(z)
    # currently don't modulate the initial state.
    return _mtlds_model_forward_parvector(m, θ, U)
end

    function _mtlds_model_forward_parvector(m::MTLDS_variational, θ::AbstractVector{T},
            U::AbstractMatrix{T}, h0::AbstractVector{T}=zeros(T, m.d)) where T <: AbstractFloat
        n_u, tT = size(U)
        _A, _B = A(θ, m.d), B(θ, m.d, n_u)
        lds = Flux.Recur(LDSCell_simple_u(_A, _B, h0))
        lds_state = [lds(U[:,t]) for t in 1:tT]
        X = hcat(lds_state...)
        if m.flag_mt_emission
            _C, d = C(θ, m.d, n_u, size(Y, 1)), dmt(θ, m.d, n_u, size(Y, 1))
            return _C * X .+ d
        end
        return m.emission(X)
    end

    function _mtlds_model_forward_parvector(m::MTLDS_variational, θ::AbstractMatrix{T},
            U::AbstractArray{T}, h0::AbstractVecOrMat{T}=zeros(T, m.d)) where T <: AbstractFloat
        n_u, tT, n_b = size(U)
        h0 = ensure_matrix(h0)
        _Bs = B(θ, m.d, n_u)
        _As = cat([A(θ[:,b], m.d) for b in 1:n_b]..., dims=3)
        lds = Flux.Recur(LDSCell_batch_u(_As, _Bs, h0))

        if m.flag_mt_emission
            _C, d = C(θ, m.d, n_u, m.emission), dmt(θ, m.d, n_u, m.emission)
            lds_state = [unsqueeze(lds(U[:,t,:]), 2) for t in 1:tT]
            X = cat(lds_state..., dims=2)
            return batch_matmul(_C, X) .+ unsqueeze(d, 2)
        else
            Ys = [unsqueeze(m.emission(lds(U[:,t,:])), 2) for t in 1:tT]
            return cat(Ys..., dims=2)
        end
    end

encode(m::MTLDS_variational, Y::AbstractArray{T}, U::AbstractArray{T}; T_enc=size(Y, 2),
    stochastic=true) where T = encode(m, vcat(Y, U); T_enc=T_enc, stochastic=stochastic)

encode(m::MTLDS_variational, Y::MaskedArray{T}, U::AbstractArray{T}; T_enc=size(Y, 2),
    stochastic=true) where T =
    encode(m, vcat(stacked(Y), U); T_enc=T_enc, stochastic=stochastic)

function encode(m::MTLDS_variational, data_enc::AbstractArray; T_enc=size(data_enc, 2), stochastic=true)
    m.flag_mt_x0 && error("Unable to modulate x0 currently. Change encoder code, and `_mtlds_model_forward`.")
    smpz, μ_z, σ_z = posterior_sample(m.mt_enc, m.mt_post, data_enc, T_enc, stochastic)
    return (smpz, μ_z, σ_z), (nothing, nothing, nothing)   # second tuple is for if x0 has separate encoder.
end

# Cannot define function on an abstract type, so we have to def these inside each subtype.
(m::MTLDS_variational)(z::AbstractArray{T,2}, U::AbstractArray{T,3}) where T = forward(m, z, U)
(m::MTLDS_variational)(Y::AbstractArray{T,3}, U::AbstractArray{T,3}) where T = reconstruct(m, Y, U)

"""
    create_mtlds(d_x, d_in, d_y, d_enc_state, d_mt; encoder=:LSTM,
    emission=:linear, d_hidden=64)

Create `MTLDS_variational` with the dimensions `d_x` of the latent system,
`d_in` of the input vectors, `d_y` of the output vectors, `d_enc_state` of
the RNN encoder, `d_mt` of the latent hierarchical variable ``z``.

### Optional arguments:

- `encoder` is the encoding RNN for the variational posterior. Can be in {`:RNN`
, `:GRU`, `:LSTM`, `:LookupTable`}, with the latter providing standard VI.
- `emission` specifies the system emission equation, which by
  default is `:linear`, which constructs a MT affine output. Otherwise, any
  callable function that maps a state vector (matrix) → outputs can be given, but
  it will not use any MT adaptation.
- `d_hidden`: number of hidden units in the one-hidden-layer MLP used for ``h_ϕ``.
- `spherical_var`: whether the emission covariance is spherical (and hence has a
  single parameter) (default: true), or is an axis-aligned ellipsoid with `n_y`
  parameters.
- `hphi`: default `:MLP`, can be `:Linear`. This determines whether the prior over
  θ is linear or nonlinear. If == `:MLP`, `d_hidden` specifies layers, o.w. it's ignored.
"""
function create_mtlds(d_x, d_in, d_y, d_enc_state, d_mt; encoder=:LSTM, emission=:linear,
    d_hidden=64, spherical_var=false, hphi=:MLP)

    if encoder == :LSTM
        init_enc = LSTM(d_y+d_in, d_enc_state)
    elseif encoder == :GRU
        init_enc = GRU(d_y+d_in, d_enc_state)
    elseif encoder == :RNN
        init_enc = RNN(d_y+d_in, d_enc_state)
    elseif encoder == :LookupTable
        init_enc = LookupTable(d_mt, 0.01f0, -1f0)  # 2/3: initial posterior noise for mu/logstd
    else
        error("encoder must be specified as :LSTM, :GRU, :RNN, :LookupTable")
    end

    enc_post = MultiDense(Dense(d_enc_state, d_mt, identity), Dense(d_enc_state, d_mt, σ))
    enc_post.Dense2.b.data .= -2   # initialize posteriors to be low variance

    n_dec = d_x*(d_x+1)/2 + d_x*d_in*2
    if emission == :linear
        emission = d_y  ### C(θ), d(θ)
        n_dec += d_y * d_x + d_y
        flag_mt_emission = true
    else
        (supertype(typeof(emission)) === Any) || error("unexpected emission object: expecting fn or :linear")
        flag_mt_emission = false
    end

    if hphi == :MLP
        hphi = mlp(d_mt, d_hidden..., Int(n_dec))
    elseif hphi == :Linear
        hphi = Dense(d_mt, Int(n_dec))
    else
        error("hphi must be specified as :MLP, :Linear")
    end

    emission_std = spherical_var ? Tracker.param([0.0f0]) : Tracker.param(zeros(Float32, d_y))
    MTLDS_variational(init_enc, enc_post, hphi, emission, emission_std, d_x, flag_mt_emission, false)
end

ext: jl


text:
function get_a(obj::FieldButcherArrayProvider)
    return jcall(obj, "getA", Vector{Vector{RealFieldElement}}, ())
end

function get_b(obj::FieldButcherArrayProvider)
    return jcall(obj, "getB", Vector{RealFieldElement}, ())
end

function get_c(obj::FieldButcherArrayProvider)
    return jcall(obj, "getC", Vector{RealFieldElement}, ())
end

function get_order(obj::EmbeddedRungeKuttaFieldIntegrator)
    return jcall(obj, "getOrder", jint, ())
end


text:
export AbstractSampler
abstract type AbstractSampler end


text:
# License is MIT: LICENSE.md

using ModuleInterfaceTools

using Test

@api use StrLiterals, StrEntities

@testset "LaTeX Entities" begin
    @test f"\<dagger>" == "†"
    #@test f"\<mscrl>" == "𝓁" # \U1f4c1
    @test f"\<c_l>" == "𝓁" # \U1f4c1
    @test f"\<nleqslant>" == "⩽̸" # \u2a7d\u338
end
@testset "Emoji Entities" begin
    @test f"\:sailboat:" == "\u26f5"
    @test f"\:ring:"     == "\U1f48d"
    @test f"\:flag-us:"  == "\U1f1fa\U1f1f8"
end
@testset "Unicode Entities" begin
    @test f"\N{end of text}" == "\x03" # \3
    @test f"\N{TIBETAN LETTER -A}" == "\u0f60"
    @test f"\N{LESS-THAN OR SLANTED EQUAL TO}" == "\u2a7d"
    @test f"\N{REVERSED HAND WITH MIDDLE FINGER EXTENDED}" == "\U1f595"
end
@testset "HTML Entities" begin
    @test f"\≰"    == "\u2270"
    @test f"\𝒫"   == "\U1d4ab"
    @test f"\≨"   == "\u2268\ufe00"
end
@testset "Unicode constants" begin
    @test f"\u{3}"     == "\x03"
    @test f"\u{f60}"   == "\u0f60"
    @test f"\u{2a7d}"  == "\u2a7d"
    @test f"\u{1f595}" == "\U0001f595"
end


text:
"""
    test_Poly

"""
ROOT = pwd()
include(joinpath(ROOT, "src_PolynomialOptim", "PolynomialOptim.jl"))

x = Variable("x", Complex)
y = Variable("y", Complex)
z = Variable("z", Real)
b = Variable("b", Bool)

## Polynomial algebra
pt = Point(Dict(x=>2, y=>1+im, z=>0+im, b=>2.1))
print(pt)

p1 = x^2 + 3*y + conj(x) + conj(z) + b
println(p1)
evaluate(p1, pt) == 10 + 3im

p2 = y^6 - y^6 + (-y*x*b) / 4 + π*conj(b)
println(p2)

## Poly operators
@assert evaluate(x, pt) == 2
@assert evaluate(y, pt) == 1+im
@assert evaluate(conj(y), pt) == 1-im
@assert evaluate(z, pt) == 0
@assert evaluate(b, pt) == 1

@assert evaluate(x*y, pt) == 2+2im
@assert evaluate(p2, pt) == (-1-im)/2 + π

@assert evaluate(real(y), pt) == 1
@assert evaluate(imag(y), pt) == 1

@assert evaluate(x*y*conj(y), pt) == 2*(1+im)*(1-im)
@assert evaluate(x*y + conj(y), pt) == 3+im
@assert evaluate(x*y + 1 + 1im, pt) == 3+3im

## Point algebra
pt = Point(Dict(x=>2, y=>1+im, z=>0+im, b=>2.1))
pt1 = Point(Dict(z=>0+im))
pt2 = Point(Dict(x=>3, y=>-1+2im, b=>-1))
@assert merge(pt, pt1, pt2) == Point(Dict(x=>5, y=>3im, z=>0+im, b=>1))

@assert pt + pt1 - 3*pt2 == Point(Dict(x=>-7, y=>4-5im, z=>0+im, b=>1))

@assert norm(pt2, Inf) == 3
@assert norm(pt2, 1) == 6
@assert norm(pt2, 2) == √(14)

## Poly cplx2real
p = x + y
println("p: $p")

p1 = copy(p)
println("p1: $p1")

add!(p1, y*z)
println("p: $p")
println("p1: $p1")

real_p, imag_p = cplx2real(p1)

expo = Exponent(Dict(x=>Degree(1,0), y=>Degree(1,0), z=>Degree(2,0)))
println(expo)
real_p, imag_p = cplx2real(expo)
println(expo)

println(real_p)
println(imag_p)

println("p1: $p1")
println("real part: $real_p")
println("imag part: $imag_p")

### Exponent generation...
x = Variable("x", Complex)
y = Variable("y", Complex)
z = Variable("z", Complex)
μ = Variable("μ", Complex)

input_vars = [x, y, z, μ]
d = 3

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
# policy file better sanity fewer contacts
pinft = 0.5*pinf_o

pshopt = [[1,1,1],[1,1,1]] - pspar*[[0.15,0.5,0],[0.15,0.5,0]]

shopmeett = [8,28,10]

tauadj = tauadj_o

socialmaxyyt = socialmaxyy_o
socialmaxoyt = socialmaxoy_o
socialmaxoot = socialmaxoo_o
socialmaxyot = socialmaxyo_o


text:
function bootstrap(x, t, n_trees)
    n_features = size(x)[2]
    n_features_forest = Int(floor(sqrt(n_features)))
    bootstrapped_x = Array{Float64}(undef, length(t), n_features_forest, n_trees)
    bootstrapped_t = Array{Any}(undef, n_trees, length(t))
    using_feature = Array{Int}(undef, n_trees, n_features_forest)
    newaxis = [CartesianIndex()]
    for i in 1 : n_trees
        ind = rand(1:length(t), length(t))
        col = sample(1 : n_features, n_features_forest)
        k = x[ind, col]
        bootstrapped_x[:, :, i] = k[:, :, newaxis]
        bootstrapped_t[i, :] = t[ind]
        using_feature[i, :] = col
    end
    return bootstrapped_x, bootstrapped_t, using_feature
end

"""
    RandomForest(nt; alpha = 0.01)
RandomForest Model. `nt` is the number of trees, and `alpha` is the same as `alpha` in DecisionTree.

# Example
```jldoctest classification
julia> model = RandomForest(10)
RandomForest(0.01, 10, DecisionTree[], Vector{Any}[], #undef)

julia> fit!(model, x, t)
10×1 Matrix{Int64}:
 1
 2
 2
 2
 2
 1
 1
 1
 1
 1

julia> println(predict(model, x))
Any[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 1, 2, 2, 1, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2]
```
"""
mutable struct RandomForest{N}
    α::Float64
    n_trees::Int64
    forest::Array{DecisionTree, 1}
    classes::Array{Array{Any, 1}, 1}
    using_feature::Array{Int64, 2}
    RandomForest(nt::Int64; alpha::Float64 = 0.01) = new{nt}(alpha, nt, Array{DecisionTree}(undef, 0), Array{Array}(undef, 0))
end

function fit!(model::RandomForest, x, t)
    forest = Array{DecisionTree}(undef, model.n_trees)
    classes = Array{Array}(undef, model.n_trees)
    bootstrapped_x, bootstrapped_t, using_feature = bootstrap(x, t, model.n_trees)
    for i in 1:model.n_trees
        tree = DecisionTree(alpha = model.α)
        fit!(tree, bootstrapped_x[:, :, i], bootstrapped_t[i, :])
        forest[i], classes[i] = tree, tree.classes
    end
    model.forest = forest
    model.classes = classes
    model.using_feature = using_feature
end

function predict(model::RandomForest, x)
    solution = Array{Any}(undef, model.n_trees, size(x, 1))
    for (tree, feature, i) in zip(model.forest, 1:size(model.using_feature, 1), 1:model.n_trees)
        solution[i, :] = predict(tree, x[:, model.using_feature[feature, :]])
    end
    predicts = Array{Any}(undef, size(x, 1))
    for i in 1:size(x, 1)
        class, counts = count_unique(solution[:, i])
        predicts[i] = class[argmax(counts)]
    end
    return predicts
end


text:
"""
	resample(x, resolution::Number; crs, method)
	resample(x; to, method)
    resample(xs...; to=first(xs), method)

`resample` uses `ArchGDAL.gdalwarp` to resample a [`Raster`](@ref) or
[`RasterStack`](@ref) to a new `resolution` and optionally new `crs`,
or to snap to the bounds, resolution and crs of the object `to`.

# Arguments

- `x`: the object to resample.
- `resolution`: a `Number` specifying the resolution for the output.
    If the keyword argument `crs` (described below) is specified, `resolution` must be in units of the `crs`.

# Keywords

- `to`: an `AbstractRaster` whos resolution, crs and bounds will be snapped to.
    For best results it should roughly cover the same extent, or a subset of `A`.
- `crs`: A `GeoFormatTypes.GeoFormat` such as `EPSG(x)` or `WellKnownText(string)` specifying an
    output crs (`A` will be reprojected to `crs` in addition to being resampled). Defaults to `crs(A)`
- `method`: A `Symbol` or `String` specifying the method to use for resampling.
    From the docs for [`gdalwarp`](https://gdal.org/programs/gdalwarp.html#cmdoption-gdalwarp-r):
    * `:near`: nearest neighbour resampling (default, fastest algorithm, worst interpolation quality).
    * `:bilinear`: bilinear resampling.
    * `:cubic`: cubic resampling.
    * `:cubicspline`: cubic spline resampling.
    * `:lanczos`: Lanczos windowed sinc resampling.
    * `:average`: average resampling, computes the weighted average of all non-NODATA contributing pixels.
        rms root mean square / quadratic mean of all non-NODATA contributing pixels (GDAL >= 3.3)
    * `:mode`: mode resampling, selects the value which appears most often of all the sampled points.
    * `:max`: maximum resampling, selects the maximum value from all non-NODATA contributing pixels.
    * `:min`: minimum resampling, selects the minimum value from all non-NODATA contributing pixels.
    * `:med`: median resampling, selects the median value of all non-NODATA contributing pixels.
    * `:q1`: first quartile resampling, selects the first quartile value of all non-NODATA contributing pixels.
    * `:q3`: third quartile resampling, selects the third quartile value of all non-NODATA contributing pixels.
    * `:sum`: compute the weighted sum of all non-NODATA contributing pixels (since GDAL 3.1)

    Where NODATA values are set to `missingval`.

Note: `missingval` of `missing` does not work with GDAL. Use `replace_missing(A, newmissingval)` to
assign a missing value before using `resample` if the current value is `missing`. This will be
automated in future versions.

# Example

Resample a WorldClim layer to match an EarthEnv layer:

```jldoctest
using Rasters, Plots
A = Raster(WorldClim{Climate}, :prec; month=1)
B = Raster(EarthEnv{HabitatHeterogeneity}, :evenness)

a = plot(A)
b = plot(resample(A; to=B))

savefig(a, "build/resample_example_before.png")
savefig(b, "build/resample_example_after.png")
# output
```

### Before `resample`:

![before resample](resample_example_before.png)

### After `resample`:

![after resample](resample_example_after.png)

$EXPERIMENTAL
"""
function resample end
resample(xs::RasterStackOrArray...; kw...) = resample(xs; kw...)
function resample(ser::AbstractRasterSeries, args...; kw...)
    map(x -> resample(x, args...; kw...), ser)
end
function resample(xs::Union{Tuple,NamedTuple}; to=first(xs), kw...)
    map(x -> resample(x; to, kw...), xs)
end
function resample(A::RasterStackOrArray, resolution::Number;
    crs::GeoFormat=crs(A), method=:near, kw...
)
    wkt = convert(String, convert(WellKnownText, crs))
    flags = Dict(
        :t_srs => wkt,
        :tr => [resolution, resolution],
        :r => method,
    )
    return warp(A, flags; kw...)
end
function resample(A::RasterStackOrArray; to, method=:near, kw...)
    all(hasdim(to, (XDim, YDim))) || throw(ArgumentError("`to` mush have both XDim and YDim dimensions to resize with GDAL"))
    if sampling(to, XDim) isa Points
        to = set(to, dims(to, XDim) => Intervals(Start()))
    end
    if sampling(to, YDim) isa Points
        to = set(to, dims(to, YDim) => Intervals(Start()))
    end

    wkt = convert(String, convert(WellKnownText, crs(to)))
    xres, yres = map(abs ∘ step, span(to, (XDim, YDim)))
    (xmin, xmax), (ymin, ymax) = bounds(to, (XDim, YDim))
    flags = Dict(
        :t_srs => wkt,
        :tr => [yres, xres],
        :te => [xmin, ymin, xmax, ymax],
        :r => method,
    )
    return warp(A, flags; kw...)
end


text:
module PrairieIO
using LightXML
using Images
using FileIO
using FixedPointNumbers
using Compat
using SimpleTraits,Unitful

import Compat: UTF8String,String,ASCIIString

include("importPrairie.jl")

export
   importPrairie,
   getPrairieFrames

end # module


text:
module MAGMA
using CUDAdrv, CUDAapi, CUDAnative, CuArrays
using CEnum

using LinearAlgebra: triu, tril, dot, checksquare, chkstride1

# export wrappers in svds
export magma_gels!, magma_gesvd!, magma_gesdd!, magmaInit, magmaFinalize, magma_gebrd!

# export wrappers in linearsystemsolver
export magma_gesv!, magma_getri!, magma_getrs!, magma_getrf!

# export some wrappers in clang auto-generation
export magma_init, magma_finalize

include("clang/libmagma_common.jl")
include("clang/libmagma_api.jl")

"""
the path to magma binary
"""
const libmagma = "/usr/local/magma/lib/libmagma.so"

macro magmafunc(function_name)
	return Expr(:quote, Symbol("magma_", function_name))
end

macro magmafunc_gpu(function_name)
	return Expr(:quote, Symbol("magma_", function_name, "_gpu"))
end

# ! for parts of subroutines which can not do the query
macro magmafunc_nb(function_name)
	return Expr(:quote, Symbol("magma_get_", function_name, "_nb"))
end

const magmaTypeTuple= (Float32, Float64, ComplexF32, ComplexF64)
const magmaTypeList = ["Float32", "Float64", "ComplexF32", "ComplexF64"]
const magmaTypeDict = Dict(Float32=>"s", Float64=>"d", ComplexF32=>"c", ComplexF64=>"z",
"Float32"=>"s", "Float64"=>"d", "ComplexF32"=>"c", "ComplexF64"=>"z",)

# include the files of subroutines
include("dense/dense.jl")

end  # modul MAGMA


text:
module ShiftedArrays

import Base: checkbounds, getindex, setindex!, parent, size, axes
export ShiftedArray, ShiftedVector, shifts, default
export CircShiftedArray, CircShiftedVector
export lag, lead

include("shiftedarray.jl")
include("circshiftedarray.jl")
include("lag.jl")
include("circshift.jl")
include("fftshift.jl")

end


text:
const Min_Max = NamedTuple{(:min, :max), Tuple{Float64, Float64}}
const Up_Down = NamedTuple{(:up, :down), Tuple{Float64, Float64}}
const From_To_Float = NamedTuple{(:from, :to), Tuple{Float64, Float64}}
const FromTo_ToFrom_Float = NamedTuple{(:from_to, :to_from), Tuple{Float64, Float64}}

"From http://www.pserc.cornell.edu/matpower/MATPOWER-manual.pdf Table B-4"

IS.@scoped_enum(GeneratorCostModels, PIECEWISE_LINEAR = 1, POLYNOMIAL = 2,)

IS.@scoped_enum(AngleUnits, DEGREES = 1, RADIANS = 2,)

IS.@scoped_enum(BusTypes, ISOLATED = 1, PQ = 2, PV = 3, REF = 4, SLACK = 5,)

IS.@scoped_enum(
    LoadModels,
    ConstantImpedance = 1, # Z
    ConstantCurrent = 2,   # I
    ConstantPower = 3,     # P
)

"From https://www.eia.gov/survey/form/eia_923/instructions.pdf"

IS.@scoped_enum(
    PrimeMovers,
    BA = 1,  # Energy Storage, Battery
    BT = 2,  # Turbines Used in a Binary Cycle (including those used for geothermal applications)
    CA = 3,  # Combined-Cycle – Steam Part
    CC = 4,  # Combined-Cycle - Aggregated Plant *augmentation of EIA
    CE = 5,  # Energy Storage, Compressed Air
    CP = 6,  # Energy Storage, Concentrated Solar Power
    CS = 7,  # Combined-Cycle Single-Shaft Combustion turbine and steam turbine share a single generator
    CT = 8,  # Combined-Cycle Combustion Turbine Part
    ES = 9,  # Energy Storage, Other (Specify on Schedule 9, Comments)
    FC = 10,  # Fuel Cell
    FW = 11,  # Energy Storage, Flywheel
    GT = 12,  # Combustion (Gas) Turbine (including jet engine design)
    HA = 13,  # Hydrokinetic, Axial Flow Turbine
    HB = 14,  # Hydrokinetic, Wave Buoy
    HK = 15,  # Hydrokinetic, Other
    HY = 16,  # Hydraulic Turbine (including turbines associated with delivery of water by pipeline)
    IC = 17,  # Internal Combustion (diesel, piston, reciprocating) Engine
    PS = 18,  # Energy Storage, Reversible Hydraulic Turbine (Pumped Storage)
    OT = 19,  # Other – Specify on SCHEDULE 9.
    ST = 20,  # Steam Turbine (including nuclear, geothermal and solar steam; does not include combined-cycle turbine)
    PVe = 21,  # Photovoltaic *renaming from EIA PV to PVe to avoid conflict with BusType.PV
    WT = 22,  # Wind Turbine, Onshore
    WS = 23,  # Wind Turbine, Offshore
)

"AER Aggregated Fuel Code From https://www.eia.gov/survey/form/eia_923/instructions.pdf"

IS.@scoped_enum(
    ThermalFuels,
    COAL = 1,  # COL    # Anthracite Coal and Bituminous Coal
    WASTE_COAL = 2,  # WOC    # Waste/Other Coal (includes anthracite culm, gob, fine coal, lignite waste, waste coal)
    DISTILLATE_FUEL_OIL = 3,  # DFO # Distillate Fuel Oil (Diesel, No. 1, No. 2, and No. 4
    WASTE_OIL = 4,  # WOO    # Waste Oil Kerosene and JetFuel Butane, Propane,
    PETROLEUM_COKE = 5,  # PC  # Petroleum Coke
    RESIDUAL_FUEL_OIL = 6,     # RFO # Residual Fuel Oil (No. 5, No. 6 Fuel Oils, and Bunker Oil)
    NATURAL_GAS = 7,  # NG    # Natural Gas
    OTHER_GAS = 8,  # OOG    # Other Gas and blast furnace gas
    NUCLEAR = 9,  # NUC # Nuclear Fission (Uranium, Plutonium, Thorium)
    AG_BIPRODUCT = 10,  # ORW    # Agricultural Crop Byproducts/Straw/Energy Crops
    MUNICIPAL_WASTE = 11,  # MLG    # Municipal Solid Waste – Biogenic component
    WOOD_WASTE = 12,  # WWW     # Wood Waste Liquids excluding Black Liquor (BLQ) (Includes red liquor, sludge wood, spent sulfite liquor, and other wood-based liquids)
    GEOTHERMAL = 13,  # GEO     # Geothermal
    OTHER = 14,  # OTH     # Other
)

IS.@scoped_enum(StateTypes, Differential = 1, Algebraic = 2, Hybrid = 3,)

PS_MAX_LOG = parse(Int, get(ENV, "PS_MAX_LOG", "50"))
DEFAULT_BASE_MVA = 100.0

const POWER_SYSTEM_STRUCT_DESCRIPTOR_FILE =
    joinpath(dirname(pathof(PowerSystems)), "descriptors", "power_system_structs.json")

const DEFAULT_SYSTEM_FREQUENCY = 60.0

const INFINITE_TIME = 1e6
const START_COST = 1e8


text:
import DataFrames

"""
"""
function fix_column_types!(
        df::DataFrames.AbstractDataFrame;
        categorical_feature_names::AbstractVector{Symbol} = Symbol[],
        continuous_feature_names::AbstractVector{Symbol} = Symbol[],
        categorical_label_names::AbstractVector{Symbol} = Symbol[],
        continuous_label_names::AbstractVector{Symbol} = Symbol[],
        float_type::Type{<:AbstractFloat} = Float64,
        )::Nothing
    function make_categorical_column(
            old_vector::AbstractVector,
            )::Vector{Any}
        new_vector = Vector{Any}(undef, length(old_vector))
        for i = 1:length(old_vector)
            old_vector_ith_element = old_vector[i]
            if DataFrames.ismissing(old_vector_ith_element)
                new_vector[i] = DataFrames.missing
            else
                new_vector[i] = string(old_vector_ith_element)
            end
        end
        return new_vector
    end
    function make_continuous_column(
            old_vector::AbstractVector,
            )::Vector{Any}
        new_vector = Vector{Any}(undef, length(old_vector))
        for i = 1:length(old_vector)
            old_vector_ith_element = old_vector[i]
            if DataFrames.ismissing(old_vector_ith_element)
                new_vector[i] = DataFrames.missing
            else
                new_vector[i] = float_type(old_vector_ith_element)
            end
        end
        return new_vector
    end
    for x in categorical_feature_names
        transform_columns!(df, make_categorical_column, x)
    end
    for x in continuous_feature_names
        transform_columns!(df, make_continuous_column, x)
    end
    for x in categorical_label_names
        transform_columns!(df, make_categorical_column, x)
    end
    for x in continuous_label_names
        transform_columns!(df, make_continuous_column, x)
    end
    for x in DataFrames.names(df)
        transform_columns!(df, fix_type, x)
    end
    return nothing
end

"""
"""
function check_column_types(
        df::DataFrames.AbstractDataFrame;
        categorical_feature_names::AbstractVector{Symbol} = Symbol[],
        continuous_feature_names::AbstractVector{Symbol} = Symbol[],
        categorical_label_names::AbstractVector{Symbol} = Symbol[],
        continuous_label_names::AbstractVector{Symbol} = Symbol[],
        )::Nothing
    for column_name in DataFrames.names(df)
        column_eltype=eltype(
            collect(
                DataFrames.skipmissing(
                    df[column_name]
                    )
                )
            )
        if column_name in categorical_feature_names
            if column_eltype<:AbstractString
            else
                error(
                    string(
                        "Column \"",
                        column_name,
                        "\" has eltype() \"",
                        column_eltype,
                        "\". However, this column is categorical,",
                        "and therefore its eltype() must be <: ",
                        "AbstractString.",
                        )
                    )
            end
        elseif column_name in categorical_label_names
            if column_eltype<:AbstractString
            else
                error(
                    string(
                        "Column \"",
                        column_name,
                        "\" has eltype() \"",
                        column_eltype,
                        "\". However, this column is categorical,",
                        "and therefore its eltype() must be <: ",
                        "AbstractString.",
                        )
                    )
            end
        elseif column_name in continuous_feature_names
            if column_eltype<:AbstractFloat
            else
                error(
                    string(
                        "Column \"",
                        column_name,
                        "\" has eltype() \"",
                        column_eltype,
                        "\". However, this column is continuous,",
                        "and therefore its eltype() must be <:",
                        "AbstractFloat.",
                        )
                    )
            end
        elseif column_name in continuous_label_names
            if column_eltype<:AbstractFloat
            else
                error(
                    string(
                        "Column \"",
                        column_name,
                        "\" has eltype() \"",
                        column_eltype,
                        "\". However, this column is continuous,",
                        "and therefore its eltype() must be <:",
                        "AbstractFloat.",
                        )
                    )
            end
        else
            if column_eltype<:AbstractString
            elseif column_eltype<:AbstractFloat
                error(
                    string(
                        "Column \"",
                        column_name,
                        "\" has eltype() \"",
                        column_eltype,
                        "\". However, we only allow AbstractStrings and ",
                        "AbstractFloats. Use T <: AbstractString if",
                        "it is a categorical column. Use T <: ",
                        "AbstractFloat if it is a continuous column."
                        )
                    )
            else
            end
        end
    end
    return nothing
end


text:
const solver = CSDP.CSDPSolver(printlevel=0)

import MathProgBase
const MPB_test_path = joinpath(dirname(pathof(MathProgBase)), "..", "test")

@testset "Linear tests" begin
    include(joinpath(MPB_test_path, "linproginterface.jl"))
    linprogsolvertest(solver, 1e-6)
end

@testset "Conic tests" begin
    include(joinpath(MPB_test_path, "conicinterface.jl"))
    # FIXME fails on Windows 32 bits... Maybe I should put linear vars/cons
    # in a diagonal matrix in SemidefiniteModels.jl instead of many 1x1 blocks
    @static if !Sys.iswindows() || Sys.WORD_SIZE != 32
        @testset "Conic linear tests" begin
            coniclineartest(solver, duals=true, tol=1e-6)
        end

# :Error for SOC1 on Mac OS: https://travis-ci.org/JuliaOpt/CSDP.jl/jobs/601461403#L389-L391
#@testset "Conic SOC tests" begin
#    conicSOCtest(CSDP.CSDPSolver(printlevel=0, write_prob="soc.prob"), duals=true, tol=1e-6)
#end

# CSDP returns :Suboptimal for SOCRotated1
#        @testset "Conic SOC rotated tests" begin
#            conicSOCRotatedtest(solver, duals=true, tol=1e-6)
#        end
    end

    @testset "Conic SDP tests" begin
        conicSDPtest(solver, duals=false, tol=1e-6)
    end
end


text:
# bbox is a "bounding box" representing outer limits of latitude & longitude: [bottomlat leftlon; toplat rightlon]
# rasterdensity is number of pixels per degree, so rasterdensity = 1/res with res = raster resolution measured in degrees
# TO DO: make these consistent with lon, lat convention in rest of code, use resolution instead of rasterdensity
function bbox2ranges(bbox, rasterdensity)
    latindexes, lonindexes = bbox2indexes(bbox, rasterdensity)
    latindex = latindexes[1] : latindexes[2]
    lonindex = lonindexes[1] : lonindexes[2]
    return latindex, lonindex
end

function bbox2indexes(bbox, rasterdensity)
    latindexes = round.(Int, reverse(rasterdensity*(90 .- bbox[:,1]), dims=1) + [1 0]')
    lonindexes = round.(Int, rasterdensity*(bbox[:,2] .- (-180)) + [1 0]')
    return latindexes, lonindexes
end

function roundbbox(bbox, rasterdensity)
    newbbox = zeros(2,2)
    newbbox[1,:] = floor.(bbox[1,:]*rasterdensity)/rasterdensity
    newbbox[2,:] = ceil.(bbox[2,:]*rasterdensity)/rasterdensity
    return newbbox
end

filterrange(range, lowhigh) = range[(range .>= lowhigh[1]) .& (range .<= lowhigh[2])]

function getlats(bbox, rasterdensity, halfshift)
	bottomlat, toplat = bbox[:,1]
	res = 1/rasterdensity
	shift = halfshift ? res/2 : 0.0
	return bottomlat + shift : res : toplat - shift
	# return toplat - shift : -res : bottomlat + shift
end

function getlons(bbox, rasterdensity, halfshift)
	leftlon, rightlon = bbox[:,2]
	res = 1/rasterdensity
	shift = halfshift ? res/2 : 0.0
	return leftlon + shift : res : rightlon - shift
end

# Boolean function, true for a leap year (for Gregorian calendar, assumes year >= 1583).
# The Gregorian calendar has 97 leap years every 400 years:
#     Every year divisible by 4 is a leap year.
#     However, every year divisible by 100 is not a leap year.
#     However, every year divisible by 400 is a leap year after all.
# So, 1700, 1800, 1900, 2100, and 2200 are not leap years, but 1600, 2000, and 2400 are leap years.
leapyear(year) = divisible(year,4) && (!divisible(year,100) || divisible(year,400))
divisible(x,y) = x % y == 0

function meshgrid2d(x,y)
	xx = repeat(x', length(y), 1)
	yy = repeat(y, 1, length(x))
	return xx, yy
end

# Borrowed from Matlab, see fspecial('disk', rad)
function diskfilterkernel(rad)
	crad  = ceil(Int, rad-0.5)
	x,y = meshgrid2d(-crad:crad,-crad:crad)
	maxxy = max.(abs.(x),abs.(y))
	minxy = min.(abs.(x),abs.(y))
	m1 = (rad^2 .< (maxxy .+ 0.5).^2 + (minxy .- 0.5).^2) .* (minxy .- 0.5) +
	  (rad^2 .>= (maxxy .+ 0.5).^2 + (minxy .- 0.5).^2) .* sqrt.(max.(0, rad^2 .- (maxxy .+ 0.5).^2))
	m2 = (rad^2 .> (maxxy .- 0.5).^2 + (minxy .+ 0.5).^2) .* (minxy .+ 0.5) +
	  (rad^2 .<= (maxxy .- 0.5).^2 + (minxy .+ 0.5).^2) .* sqrt.(max.(rad^2 .- (maxxy .- 0.5).^2))
	sgrid = (rad^2 * (0.5*(asin.(m2/rad) .- asin.(m1/rad)) + 0.25*(sin.(2*asin.(m2/rad)) - sin.(2*asin.(m1/rad)))) +
	     		- (maxxy .- 0.5).*(m2-m1) + (m1 .- minxy .+ 0.5)) .*
	      ((((rad^2 .< (maxxy .+ 0.5).^2 + (minxy .+ 0.5).^2) .&
	      	(rad^2 .> (maxxy .- 0.5).^2 + (minxy .- 0.5).^2)) .|
	      	((minxy .== 0) .& (maxxy .- 0.5 .< rad) .& (maxxy .+ 0.5 .>= rad))))
	sgrid = sgrid .+ ((maxxy .+ 0.5).^2 + (minxy .+ 0.5).^2 .< rad^2)
	sgrid[crad+1,crad+1] = min.(pi*rad^2, pi/2)
	if ((crad>0) && (rad > crad-0.5) && (rad^2 < (crad-0.5)^2+0.25))
		m1a  = sqrt(rad^2 - (crad - 0.5).^2)
		m1n = m1a/rad
		sg0 = 2*(rad^2*(0.5*asin(m1n) + 0.25*sin(2*asin(m1n)))-m1a*(crad-0.5))
		sgrid[2*crad+1,crad+1] = sg0
		sgrid[crad+1,2*crad+1] = sg0
		sgrid[crad+1,1]        = sg0
		sgrid[1,crad+1]        = sg0
		sgrid[2*crad,crad+1]   = sgrid[2*crad,crad+1] - sg0
		sgrid[crad+1,2*crad]   = sgrid[crad+1,2*crad] - sg0
		sgrid[crad+1,2]        = sgrid[crad+1,2]      - sg0
		sgrid[2,crad+1]        = sgrid[2,crad+1]      - sg0
	end
	sgrid[crad+1,crad+1] = min(sgrid[crad+1,crad+1], 1)
	h = centered(sgrid/sum(sgrid))
end

meandrop(x; dims=dims) = dropdims(mean(x, dims=dims), dims=dims)
sumdrop(x; dims=dims) = dropdims(sum(x, dims=dims), dims=dims)

shifthalfcell(data, kernel=centered([.25 .25; .25 .25])) = imfilter(data, kernel)

function shiftallcells!(windCF)
    yearlength = size(windCF, 1)
    updateprogress = Progress(yearlength, 1)
    kernel = centered([.25 .25; .25 .25])
    for i=1:yearlength
        windCF[i,:,:] = shifthalfcell(windCF[i,:,:], kernel)
        next!(updateprogress)
    end
end

function max!(x::AbstractArray, val)
    @inbounds for i in eachindex(x)
        x[i] = max(x[i], val)
    end
	x
end

selfmap!(f,x) = map!(f,x,x)

# Apply a function fn to the data matrix by chunks to avoid memory issues
function gridsplit(data::AbstractArray, fn::Function, resulttype::DataType; nmax=9000, overlap=250)
    rows, cols = size(data)
    nparts = length(1:nmax:rows) * length(1:nmax:cols)
    nparts > 1 && println("\nSplitting data matrix into $nparts parts to avoid memory issues...")
    result = zeros(resulttype, size(data))
    part = 0
    for r = 1:nmax:rows, c = 1:nmax:cols
        part += 1
        nparts > 1 && println("\nPart $part/$nparts:")
        rowrange_in = max(r-overlap, 1):min(r-1+nmax+overlap, rows)
        colrange_in = max(c-overlap, 1):min(c-1+nmax+overlap, cols)
        rowrange_out = r:min(r-1+nmax, rows)
        colrange_out = c:min(c-1+nmax, cols)
        rowdelta = rowrange_out[1] - rowrange_in[1] + 1
        coldelta = colrange_out[1] - colrange_in[1] + 1
        rowrange_delta = rowdelta:(rowdelta + length(rowrange_out) - 1)
        colrange_delta = coldelta:(coldelta + length(colrange_out) - 1)
        result[rowrange_out,colrange_out] = fn(data[rowrange_in,colrange_in])[rowrange_delta,colrange_delta]
    end
    return result
end

# function gridsurface(x, y; color=color, nmax=9000, args...)
#     rows, cols = size(color)
#     rowchunks, colchunks = ceil(Int, rows/nmax), ceil(Int, cols/nmax)
#     rowsplit, colsplit = ceil(Int, rows/rowchunks), ceil(Int, cols/colchunks)
#     nparts = rowchunks * colchunks
#     nparts > 1 && println("\nSplitting surface plat into $nparts parts to avoid Makie GPU limits...")
#     part = 0
#     scene = Scene()
#     for r = 1:rowsplit:rows, c = 1:colsplit:cols
#         part += 1
#         nparts > 1 && println("\nPart $part/$nparts:")
#         rowrange = r:min(r-1+rowsplit, rows)
#         colrange = c:min(c-1+colsplit, cols)
#         surface!(x[rowrange,colrange], y[rowrange,colrange]; color=color[rowrange,colrange], args...)
#     end
#     return scene
# end

row2lon(row::Int, res) = (row - 0.5) * res - 180
col2lat(col::Int, res) = 90 - (col - 0.5) * res
lon2row(lon::Real, res) = floor(Int, mod(180+lon, 360) / res) + 1
lat2col(lat::Real, res) = floor(Int, min(180-res/2, 90-lat) / res) + 1

function rowcol2lonlat(rowcol::Tuple{Int,Int}, res)
    row, col = rowcol
    @assert 1 <= row <= 360/res
    @assert 1 <= col <= 180/res
    return row2lon(row,res), col2lat(col,res)
end

function lonlat2rowcol(lonlat::Tuple{<:Real,<:Real}, res)
    lon, lat = lonlat
    @assert -90 <= lat <= 90
    return lon2row(lon,res), lat2col(lat, res)
end

# convert indexes of datasets with 0.01 degree resolution to indexes of ERA5 resolution (0.28125 degrees)
function eraranges(lonrange, latrange, res, erares)
    topleft = rowcol2lonlat((lonrange[1],latrange[1]), res)
    bottomright = rowcol2lonlat((lonrange[end],latrange[end]), res)
    row1, col1 = lonlat2rowcol(topleft, erares)
    row2, col2 = lonlat2rowcol(bottomright, erares)
    nlons = round(Int, 360/erares)
    eralonranges = row2 >= row1 ? [row1:row2] : [row1:nlons, 1:row2]
    eralatrange = col1:col2
    return eralonranges, eralatrange
end

function eraindexlookup(lons, lats, eralonranges, eralatrange)
    eralonrange = length(eralonranges) == 1 ? eralonranges[1] : [eralonranges[1]; eralonranges[2]]
    lonindexlookup = zeros(Int,1280)
    latindexlookup = zeros(Int,640)
    for (i,r) in enumerate(eralonrange)
        lonindexlookup[r] = i
    end
    for (i,r) in enumerate(eralatrange)
        latindexlookup[r] = i
    end
    eralonindexes = lonindexlookup[lon2row.(lons,0.28125)]
    eralatindexes = latindexlookup[lat2col.(lats,0.28125)]
    return eralonindexes, eralatindexes
end

# Automatically detect the "bounding box" of nonzero data in a matrix.
# Returns a tuple of indexes of the box containing data, (lon,lat) or (row,col).
function getbboxranges(regions::AbstractMatrix, padding::Int=0)
    data = (regions .> 0) .& (regions .!= NOREGION)
    lonrange = dataindexes_lon(vec(any(data, dims=2)), padding)     # all longitudes with region data
    latrange = dataindexes_lat(vec(any(data, dims=1)), padding)     # all latitudes with region data
    return lonrange, latrange
end

# Given a vector of Bool indicating nonzero data along latitudes, return a range
# of indexes indicating the area where there is data.
# Optionally add some padding elements to the data area (to ensure that there is
# some offshore territory around the regions).
function dataindexes_lat(latdata::AbstractVector, padding::Int=0)
    len = length(latdata)
    first, last = extrema(findall(latdata))         # first & last indexes of the data region
    return max(1, first-padding):min(len, last+padding)   # lat indexes of region elements
end

# Given a vector of Bool indicating nonzero data along longitudes, return a vector
# of indexes indicating the area where there is data. Do this by eliminating the
# longest contiguous sequence of zero data, considering wraparound (at lon=180).
# Optionally add some padding elements to the data area (to ensure that there is
# some offshore territory around the regions).
function dataindexes_lon(londata::AbstractVector, padding::Int=0)
    len = length(londata)
    seq = longest_circular_sequence(londata, false) # first & last indexes of the longest contiguous sequence of NONregion elements
    first, last = seq[2]+1, seq[1]-1                # the rest are region elements (including "holes" in the sequence)
    last = last >= first ? last : last + len        # make sure last > first so the loop below works
    return [mod1(i,len) for i = first-padding:last+padding]   # lon indexes of region elements
end

# Return the start and end indexes of the longest contiguous sequence of v such that v[i] == x,
# allowing the sequence to wrap around the end.
function longest_circular_sequence(v::AbstractVector, x)
    longest = Int[]
    current = Int[]
    haswrapped = false
    i = 1
    len = length(v)
    sequencelength(seq) = isempty(seq) ? 0 : mod(seq[2]-seq[1], len) + 1
    while true
        if v[i] == x
            if isempty(current)
                current = [i,i]
            else
                current[2] = i
            end
        else
            longest = sequencelength(current) > sequencelength(longest) ? current : longest
            haswrapped && return longest
            current = Int[]
        end
        if !haswrapped && i==len
            haswrapped = true
            i = 1
            sequencelength(current) == len && return current
        else
            i += 1
        end
    end
end

function coordmap(len::Int, croprange)
    coords = zeros(Int, len)
    indexes = (1:len)[croprange]
    for (i,n) = enumerate(indexes)
        coords[n] = i
    end
    return coords
end

# area of a grid cell in km2 for a given latitude (in degrees) and raster resolution
rastercellarea(lat, res) = cosd(lat) * (2*6371*π/(360/res))^2

function eralonlat(options, lonrange, latrange)
    @unpack res, erares = options

    lons = (-180+res/2:res:180-res/2)[lonrange]         # longitude values (pixel center)
    lats = (90-res/2:-res:-90+res/2)[latrange]          # latitude values (pixel center)
    cellarea = rastercellarea.(lats, res)

    eralonranges, eralatrange = eraranges(lonrange, latrange, res, erares)
    eralonrange = length(eralonranges) == 1 ? eralonranges[1] : [eralonranges[1]; eralonranges[2]]
    eralons = (-180+erares/2:erares:180-erares/2)[eralonrange]     # longitude values (pixel center)
    eralats = (90-erares/2:-erares:-90+erares/2)[eralatrange]      # latitude values (pixel center)

    lonmap = coordmap(round(Int, 360/res), lonrange)      # map full longitude indexes to cropped indexes
    latmap = coordmap(round(Int, 180/res), latrange)      # ditto latitude

    return eralons, eralats, lonmap, latmap, cellarea
end

function uncrop(croppedarray, lonrange, latrange, res)
    nlon, nlat = round(Int, 360/res), round(Int, 180/res)
    full = zeros(eltype(croppedarray), nlon, nlat)
    full[lonrange, latrange] = croppedarray
    return full
end

function drawmap(mapdata)
    skip = ceil(Int, maximum(size(mapdata))/3000)
    mirrormap = reverse(mapdata[1:skip:end,1:skip:end], dims=2)
    # display(heatmap(mirrormap, size=(1200, 900), c=[cgrad(:viridis)[x] for x in 0.0:0.2:1.0]))
    # display(heatmap(mirrormap, size=(1200, 900), c=cgrad(:viridis, [0, 0.5, 1])))
    heatmap(mirrormap)
    # display(countmap(mirrormap[:]))
end

# I can't believe Makie can't do this more conveniently.
function plotlines(x,y)
    scene = Scene()
    for i = 1:size(y,2)
        lines!(scene, x, y[:,i], color=ColorSchemes.Set1_9[i])
    end
    scene
end

function drawregionmap(regionname)
    plotly()
    regions, offshoreregions, regionlist, lonrange, latrange = loadregions(regionname)
    skip = ceil(Int, maximum(size(mapdata))/3000)
    reg = reverse(regions[1:skip:end,1:skip:end]', dims=1)
    display(heatmap(reg .+ (reg.>0).*20, size=(1200, 900)))
    reg = reverse(offshoreregions[1:4:end,1:4:end]', dims=1)
    display(heatmap(reg .+ (reg.>0).*20, size=(1200, 900)))
end

function resize_categorical(regions, regionlist, lonrange, latrange; skipNOREGION=false)
    res = 0.01          # resolution of auxiliary datasets [degrees per pixel]
    erares = 0.28125    # resolution of ERA5 datasets [degrees per pixel]

    eralons, eralats, lonmap, latmap, cellarea = eralonlat(Dict(:res=>res, :erares=>erares), lonrange, latrange)
    numreg = length(regionlist)
    nlons, nlats = length(eralons), length(eralats)

    RegType = eltype(regions)
    smallregions = zeros(RegType, nlons, nlats)
    count0 = Dict(i => 0 for i = 0:numreg)
    count0[NOREGION] = 0

    for (j, eralat) in enumerate(eralats)
        colrange = latmap[lat2col(eralat+erares/2, res):lat2col(eralat-erares/2, res)-1]
        for (i, eralon) in enumerate(eralons)
            # get all high resolution row and column indexes within this ERA5 cell
            rowrange = lonmap[lon2row(eralon-erares/2, res):lon2row(eralon+erares/2, res)-1]
            count = copy(count0)
            foundnonzero = false
            for c in colrange, r in rowrange
                (c == 0 || r == 0) && continue
                reg = regions[r,c]
                skipNOREGION && (reg == 0 || reg == NOREGION) && continue
                count[reg] += 1
                foundnonzero = true
            end
            smallregions[i,j] = foundnonzero ? findmax(count)[2] : 0
        end
    end

    return smallregions
end

function matlab2elin(; gisregion="Europe8", year=2018)
    filenamesuffix = ""
    _, _, regionlist, _, _ = loadregions(gisregion)

    # CF_pvrooftop, capacity_pvrooftop
    region = string.(regionlist)
    tech = ["PVplantA", "PVplantB", "PVroof", "CSPA", "CSPB", "WindonshoreA", "WindonshoreB", "Windoffshore"]
    classname = ["PVPA", "PVPB", "PVR", "CSPA", "CSPB", "WONA", "WONB", "WOFF"]
    capvar = ["capacity_pvplantA", "capacity_pvplantB", "capacity_pvrooftop", "capacity_cspplantA", "capacity_cspplantB",
                "capacity_onshoreA", "capacity_onshoreB", "capacity_offshore"]
    cfvar = ["CFtime_pvplantA", "CFtime_pvplantB", "CFtime_pvrooftop", "CFtime_cspplantA", "CFtime_cspplantB",
                "CFtime_windonshoreA", "CFtime_windonshoreB", "CFtime_windoffshore"]

    winddata = matread(in_datafolder("output", "GISdata_wind$(year)_$gisregion$filenamesuffix.mat"))
    solardata = matread(in_datafolder("output", "GISdata_solar$(year)_$gisregion$filenamesuffix.mat"))
    data = merge(winddata, solardata)

    # read number of classes from wind & solar GIS output
    nwindclasses = [size(winddata[varname], 2) for varname in capvar[6:8]]
    nsolarclasses = [size(solardata[varname], 2) for varname in capvar[1:5]]
    nclasses = [nsolarclasses; nwindclasses]

    for t = 1:length(tech)
        open(in_datafolder("output", "capacity_$(tech[t]).inc"), "w") do f
            for (r,reg) in enumerate(region)
                for c = 1:nclasses[t]
                    val = data[capvar[t]][r,c]
                    !isnan(val) && val > 0 && @printf(f, "%-3s . %s%-2d %12.6f\n", reg, classname[t], c, val)
                end
            end
        end
        open(in_datafolder("output", "cf_$(tech[t]).inc"), "w") do f
            for (r,reg) in enumerate(region)
                for c = 1:nclasses[t]
                    for h = 1:8760
                        val = data[cfvar[t]][h,r,c]
                        !isnan(val) && val > 0 && @printf(f, "%-3s . %s%-2d . h%04d %10.6f\n", reg, classname[t], c, h, val)
                    end
                end
            end
        end
    end
end

function dms2deg(dms)
    sg = sign(dms)
    dms = abs(dms)
    d = sg*floor(dms/10000)
    dms = dms-10000*d
    m = floor(dms/100)
    return d + m/60 + (dms-100*m)/3600
end

# Capital recovery factor for investments
CRF(r, T) = r / (1 - 1/(1+r)^T)

in_datafolder(names...) = joinpath(getconfig("datafolder"), names...)

max_stars_repo_name: henrikhodel/GlobalEnergyGIS

max_issues_repo_name: henrikhodel/GlobalEnergyGIS

max_forks_repo_name: henrikhodel/GlobalEnergyGIS


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module SubFoo1
thissub = "Foo1.SubFoo1"
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
function spectrumsweep(p, μlist, selfy = false)
    pn = reconstruct(p, μN = μlist[1])
    h = rectangle_randombounds_sc(pn, 0, 0.0, selfy = selfy)
    sp = spectrum(h, method = ArpackPackage(nev = 16, sigma = 0.0im))
    energieslist = zeros(Float64, length(sp.energies), length(μlist))
    energieslist[:,1] = sp.energies
    for i in 2:length(μlist)
        println(i/length(μlist))
        pn = reconstruct(p, μN = μlist[i])
        h = rectangle_randombounds_sc(pn, 0, 0.0, selfy = selfy)
        sp = spectrum(h,
            method = ArpackPackage(nev = 16, sigma = 0.0im))
        energieslist[:,i] = sp.energies
    end
    return energieslist, μlist
end

function spectrumsweepλ(p, λlist, selfy = false)
    pn = reconstruct(p, λ = λlist[1])
    h = rectangle_randombounds_sc(pn, 0, 0.0, selfy = selfy)
    sp = spectrum(h, method = ArpackPackage(nev = 16, sigma = 0.0im))
    energieslist = zeros(Float64, length(sp.energies), length(λlist))
    energieslist[:,1] = sp.energies
    for i in 2:length(λlist)
        println(i/length(λlist))
        pn = reconstruct(p, λ = λlist[i])
        h = rectangle_randombounds_sc(pn, 0, 0.0)
        sp = spectrum(h,
            method = ArpackPackage(nev = 16, sigma = 0.0im))
        energieslist[:,i] = sp.energies
    end
    return energieslist, λlist
end

function spectrumsweepb(p, blist, selfy = false)
    pn = reconstruct(p, EZ = SA[0, blist[1], 0])
    h = rectangle_randombounds_sc(pn, 0., 0.0, selfy = selfy)
    sp = spectrum(h, method = ArpackPackage(nev = 32, sigma = 0.0im))
    energieslist = zeros(Float64, length(sp.energies), length(blist))
    energieslist[:,1] = sp.energies
    for i in 2:length(blist)
        println(i/length(blist))
        pn = reconstruct(p, EZ = SA[0, blist[i], 0])
        h = rectangle_randombounds_sc(pn, 0, 0.0, selfy = selfy)
        sp = spectrum(h,
            method = ArpackPackage(nev = 32, sigma = 0.0im))
        energieslist[:,i] = sp.energies
    end
    return energieslist, blist
end

function spectrumsweepα(p, list, selfy = false)
    pn = reconstruct(p, α = list[1])
    h = rectangle_randombounds_sc(pn, 0, 0.0, selfy = selfy)
    @time sp = spectrum(h, method = ArpackPackage(nev = 32, sigma = 0.0im))
    energieslist = SharedArray(zeros(Float64, length(sp.energies), length(list)))
    energieslist[:,1] = sp.energies
    for i in 2:length(list) #@sync @distributed
        println(i/length(list))
        pn = reconstruct(p, α = list[i])
        h = rectangle_randombounds_sc(pn, 0, 0.0, selfy = selfy)
        sp = spectrum(h,
            method = ArpackPackage(nev = 32, sigma = 0.0im))
        energieslist[:,i] = sp.energies
    end
    return energieslist, list
end

"given somre presets and a list of angles (in degrees) computes the spectrum of the
lowest energy states as a funcion of θ"
function splittingvsrotation(p, list, selfy = true; kw...)
    list = list .* π/180
    println("first h...")
    h = rectangle_randombounds_sc(p, list[1], 0.0, selfy = selfy; kw...)
    println("first sp...")
    sp = spectrum(h, method = ArpackPackage(nev = 16, sigma = 0.0im))
    energieslist = SharedArray(zeros(Float64, length(sp.energies), length(list)))
    energieslist[:,1] = sp.energies
    @sync @distributed for i in 2:length(list)
        println(i/length(list))
        println("h ...")
        hn = rectangle_randombounds_sc(p, list[i], 0.0, selfy = selfy; kw...)
        println("sp...")
        spn = spectrum(hn,
            method = ArpackPackage(nev = 16, sigma = 0.0im))
        energieslist[:,i] = spn.energies
    end
    return energieslist, list
end

function splittingvsdisorder(p, list, selfy = true; kw...)
    println("first h...")
    h = rectangle_randombounds_sc(p, 0.0, list[1], selfy = selfy; kw...)
    println("first sp...")
    sp = spectrum(h, method = ArpackPackage(nev = 16, sigma = 0.0im))
    energieslist = SharedArray(zeros(Float64, length(sp.energies), length(list)))
    energieslist[:,1] = sp.energies
    @sync @distributed for i in 2:length(list)
        println(i/length(list))
        println("$(list[i])...")
        hn = rectangle_randombounds_sc(p, 0.0, list[i], selfy = selfy; kw...)
        println("sp...")
        spn = spectrum(hn,
            method = ArpackPackage(nev = 16, sigma = 0.0im))
        energieslist[:,i] = spn.energies
    end
    return energieslist, list
end

# function merge_max2(a, b)
#     x, y = Baselet.sort((a..., b...); by = last, rev = true)
#     return (x, y)
# end


text:
# This script solved the TD-EFIE and TD-MFIE 'at the same time'. It verifies
# the code dealing with systems of time domain equations in the special case
# of a block diagonal system.
using CompScienceMeshes, BEAST
Γ = readmesh(joinpath(@__DIR__,"sphere2.in"))

X = raviartthomas(Γ)
Y = buffachristiansen(Γ)

Δt, Nt = 0.6, 200
δ = timebasisdelta(Δt, Nt)
T0 = timebasiscxd0(Δt, Nt)
T1 = timebasisshiftedlagrange(Δt,Nt,1)
T2 = timebasisshiftedlagrange(Δt, Nt, 2)
T3 = timebasisshiftedlagrange(Δt, Nt, 3)

duration = 20 * Δt
delay = 1.5 * duration
amplitude = 1.0
gaussian = creategaussian(duration, delay, amplitude)

direction, polarisation = ẑ, x̂
E = BEAST.planewave(polarisation, direction, gaussian, 1.0)
H = direction × E

SL = TDMaxwell3D.singlelayer(speedoflight=1.0, numdiffs=0)
DL = TDMaxwell3D.doublelayer(speedoflight=1.0)
I = Identity()
N = NCross()

@hilbertspace k l
@hilbertspace j m
emfie = @discretise(
    (0.5(N⊗I) + 1.0DL)[k,j] + SL[l,m] == -1.0H[k] - E[l],
    k∈Y⊗δ, l∈X⊗δ, j∈X⊗T1, m∈X⊗T1)
xemfie = solve(emfie)


text:
#=
Functions for working with periodograms of scalar data.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date : 2014-08-21

References
----------

http://quant-econ.net/jl/estspec.html

=#
import DSP

"""
Smooth the data in x using convolution with a window of requested size and type.

##### Arguments

- `x::Array`: An array containing the data to smooth
- `window_len::Int(7)`: An odd integer giving the length of the window
- `window::AbstractString("hanning")`: A string giving the window type. Possible values
are `flat`, `hanning`, `hamming`, `bartlett`, or `blackman`

##### Returns

- `out::Array`: The array of smoothed data
"""
function smooth(x::Array, window_len::Int, window::AbstractString="hanning")
    if length(x) < window_len
        throw(ArgumentError("Input vector length must be >= window length"))
    end

    if window_len < 3
       throw(ArgumentError("Window length must be at least 3."))
    end

    if iseven(window_len)
        window_len += 1
        println("Window length must be odd, reset to $window_len")
    end

    windows = Dict("hanning" => DSP.hanning,
                   "hamming" => DSP.hamming,
                   "bartlett" => DSP.bartlett,
                   "blackman" => DSP.blackman,
                   "flat" => DSP.rect  # moving average
                   )

    # Reflect x around x[0] and x[-1] prior to convolution
    k = round(Int, window_len / 2)
    xb = x[1:k]   # First k elements
    xt = x[end-k+1:end]  # Last k elements
    s = [reverse(xb); x; reverse(xt)]

    # === Select window values === #
    if !haskey(windows, window)
        msg = "Unrecognized window type '$window'"
        print(msg * " Defaulting to hanning")
        window = "hanning"
    end

    w = windows[window](window_len)

    return conv(w ./ sum(w), s)[window_len+1:end-window_len]
end

"Version of `smooth` where `window_len` and `window` are keyword arguments"
function smooth(x::Array; window_len::Int=7, window::AbstractString="hanning")
    smooth(x, window_len, window)
end

function periodogram(x::Vector)
    n = length(x)
    I_w = abs(fft(x)).^2 ./ n
    w = 2pi * (0:n-1) ./ n  # Fourier frequencies

    # int rounds to nearest integer. We want to round up or take 1/2 + 1 to
    # make sure we get the whole interval from [0, pi]
    ind = iseven(n) ? round(Int, n / 2  + 1) : ceil(Int, n / 2)
    w, I_w = w[1:ind], I_w[1:ind]
    return w, I_w
end

function periodogram(x::Vector, window::AbstractString, window_len::Int=7)
    w, I_w = periodogram(x)
    I_w = smooth(I_w, window_len=window_len, window=window)
    return w, I_w
end

"""
Computes the periodogram

    I(w) = (1 / n) | sum_{t=0}^{n-1} x_t e^{itw} |^2

at the Fourier frequences w_j := 2 pi j / n, j = 0, ..., n - 1, using the fast
Fourier transform.  Only the frequences w_j in [0, pi] and corresponding values
I(w_j) are returned.  If a window type is given then smoothing is performed.

##### Arguments

- `x::Array`: An array containing the data to smooth
- `window_len::Int(7)`: An odd integer giving the length of the window
- `window::AbstractString("hanning")`: A string giving the window type. Possible values
are `flat`, `hanning`, `hamming`, `bartlett`, or `blackman`

##### Returns

- `w::Array{Float64}`: Fourier frequencies at which the periodogram is evaluated
- `I_w::Array{Float64}`: The periodogram at frequences `w`

"""
periodogram

"""
Compute periodogram from data `x`, using prewhitening, smoothing and recoloring.
The data is fitted to an AR(1) model for prewhitening, and the residuals are
used to compute a first-pass periodogram with smoothing.  The fitted
coefficients are then used for recoloring.

##### Arguments

- `x::Array`: An array containing the data to smooth
- `window_len::Int(7)`: An odd integer giving the length of the window
- `window::AbstractString("hanning")`: A string giving the window type. Possible values
are `flat`, `hanning`, `hamming`, `bartlett`, or `blackman`

##### Returns

- `w::Array{Float64}`: Fourier frequencies at which the periodogram is evaluated
- `I_w::Array{Float64}`: The periodogram at frequences `w`

"""
function ar_periodogram(x::Array, window::AbstractString="hanning", window_len::Int=7)
    # run regression
    x_current, x_lagged = x[2:end], x[1:end-1]  # x_t and x_{t-1}
    coefs = collect(linreg(x_lagged, x_current))

    # get estimated values and compute residual
    est = [ones(x_lagged) x_lagged] * coefs
    e_hat = x_current - est

    phi = coefs[2]

    # compute periodogram on residuals
    w, I_w = periodogram(e_hat, window, window_len)

    # recolor and return
    I_w = I_w ./ abs(1 - phi .* exp(im.*w)).^2

    return w, I_w
end