text:
letterA = (Thebes.Point3D[Thebes.Point3D(0.0, 0.0, 0.0), Thebes.Point3D(1.0, 0.0, 0.0), Thebes.Point3D(0.4, -2.0, 0.0), Thebes.Point3D(1.25, -2.0, 0.0), Thebes.Point3D(0.5, -2.5, 0.0), Thebes.Point3D(1.3125, -2.5, 0.0), Thebes.Point3D(1.0, -5.0, 0.0), Thebes.Point3D(1.5, -4.0, 0.0), Thebes.Point3D(3.0, 0.0, 0.0), Thebes.Point3D(2.0, 0.0, 0.0), Thebes.Point3D(2.6, -2.0, 0.0), Thebes.Point3D(1.75, -2.0, 0.0), Thebes.Point3D(2.5, -2.5, 0.0), Thebes.Point3D(1.6875, -2.5, 0.0), Thebes.Point3D(2.0, -5.0, 0.0), Thebes.Point3D(0.0, 0.0, 1.0), Thebes.Point3D(1.0, 0.0, 1.0), Thebes.Point3D(0.4, -2.0, 1.0), Thebes.Point3D(1.25, -2.0, 1.0), Thebes.Point3D(0.5, -2.5, 1.0), Thebes.Point3D(1.3125, -2.5, 1.0), Thebes.Point3D(1.0, -5.0, 1.0), Thebes.Point3D(1.5, -4.0, 1.0), Thebes.Point3D(3.0, 0.0, 1.0), Thebes.Point3D(2.0, 0.0, 1.0), Thebes.Point3D(2.6, -2.0, 1.0), Thebes.Point3D(1.75, -2.0, 1.0), Thebes.Point3D(2.5, -2.5, 1.0), Thebes.Point3D(1.6875, -2.5, 1.0), Thebes.Point3D(2.0, -5.0, 1.0)], Array{Int64,1}[[3, 4, 2, 1], [5, 6, 4, 3], [7, 8, 6, 5], [12, 11, 9, 10], [14, 13, 11, 12], [8, 15, 13, 14], [6, 14, 12, 4], [7, 15, 8], [16, 17, 19, 18], [18, 19, 21, 20], [20, 21, 23, 22], [25, 24, 26, 27], [27, 26, 28, 29], [29, 28, 30, 23], [19, 27, 29, 21], [23, 30, 22], [1, 16, 18, 3], [3, 18, 20, 5], [5, 20, 22, 7], [7, 22, 30, 15], [15, 30, 28, 13], [13, 28, 26, 11], [11, 26, 24, 9], [9, 24, 25, 10], [10, 25, 27, 12], [12, 27, 19, 4], [4, 19, 17, 2], [2, 17, 16, 1], [8, 23, 21, 6], [6, 21, 29, 14], [14, 29, 23, 8]])
letterB = (Thebes.Point3D[Thebes.Point3D(0.0, 0.0, 0.0), Thebes.Point3D(2.0, 0.0, 0.0), Thebes.Point3D(3.0, -0.5, 0.0), Thebes.Point3D(3.0, -1.5, 0.0), Thebes.Point3D(2.0, -2.5, 0.0), Thebes.Point3D(3.0, -3.5, 0.0), Thebes.Point3D(3.0, -4.5, 0.0), Thebes.Point3D(2.0, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(1.0, -1.0, 0.0), Thebes.Point3D(2.0, -1.0, 0.0), Thebes.Point3D(2.0, -2.0, 0.0), Thebes.Point3D(1.0, -2.0, 0.0), Thebes.Point3D(1.0, -3.0, 0.0), Thebes.Point3D(2.0, -3.0, 0.0), Thebes.Point3D(2.0, -4.0, 0.0), Thebes.Point3D(1.0, -4.0, 0.0), Thebes.Point3D(0.0, 0.0, 1.0), Thebes.Point3D(2.0, 0.0, 1.0), Thebes.Point3D(3.0, -0.5, 1.0), Thebes.Point3D(3.0, -1.5, 1.0), Thebes.Point3D(2.0, -2.5, 1.0), Thebes.Point3D(3.0, -3.5, 1.0), Thebes.Point3D(3.0, -4.5, 1.0), Thebes.Point3D(2.0, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0), Thebes.Point3D(1.0, -1.0, 1.0), Thebes.Point3D(2.0, -1.0, 1.0), Thebes.Point3D(2.0, -2.0, 1.0), Thebes.Point3D(1.0, -2.0, 1.0), Thebes.Point3D(1.0, -3.0, 1.0), Thebes.Point3D(2.0, -3.0, 1.0), Thebes.Point3D(2.0, -4.0, 1.0), Thebes.Point3D(1.0, -4.0, 1.0)], Array{Int64,1}[[10, 11, 2, 1], [11, 12, 5, 4, 3, 2], [5, 12, 13, 14, 15], [5, 15, 16, 8, 7, 6], [8, 16, 17, 9], [9, 17, 14, 13, 10, 1], [18, 19, 28, 27], [19, 20, 21, 22, 29, 28], [32, 31, 30, 29, 22], [23, 24, 25, 33, 32, 22], [26, 34, 33, 25], [18, 27, 30, 31, 34, 26], [1, 2, 19, 18], [2, 3, 20, 19], [3, 4, 21, 20], [4, 5, 22, 21], [5, 6, 23, 22], [6, 7, 24, 23], [7, 8, 25, 24], [8, 9, 26, 25], [9, 1, 18, 26], [11, 10, 27, 28], [10, 13, 30, 27], [12, 11, 28, 29], [13, 12, 29, 30], [17, 16, 33, 34], [16, 15, 32, 33], [15, 14, 31, 32], [14, 17, 34, 31]])
letterC = (Thebes.Point3D[Thebes.Point3D(1.0, 0.0, 0.0), Thebes.Point3D(2.0, 0.0, 0.0), Thebes.Point3D(3.0, -1.0, 0.0), Thebes.Point3D(3.0, -1.5, 0.0), Thebes.Point3D(2.0, -1.5, 0.0), Thebes.Point3D(2.0, -1.0, 0.0), Thebes.Point3D(1.5, -0.5, 0.0), Thebes.Point3D(1.0, -1.0, 0.0), Thebes.Point3D(1.0, -4.0, 0.0), Thebes.Point3D(1.5, -4.5, 0.0), Thebes.Point3D(2.0, -4.0, 0.0), Thebes.Point3D(2.0, -3.5, 0.0), Thebes.Point3D(3.0, -3.5, 0.0), Thebes.Point3D(3.0, -4.0, 0.0), Thebes.Point3D(2.0, -5.0, 0.0), Thebes.Point3D(1.0, -5.0, 0.0), Thebes.Point3D(0.0, -4.0, 0.0), Thebes.Point3D(0.0, -1.0, 0.0), Thebes.Point3D(1.0, 0.0, 1.0), Thebes.Point3D(2.0, 0.0, 1.0), Thebes.Point3D(3.0, -1.0, 1.0), Thebes.Point3D(3.0, -1.5, 1.0), Thebes.Point3D(2.0, -1.5, 1.0), Thebes.Point3D(2.0, -1.0, 1.0), Thebes.Point3D(1.5, -0.5, 1.0), Thebes.Point3D(1.0, -1.0, 1.0), Thebes.Point3D(1.0, -4.0, 1.0), Thebes.Point3D(1.5, -4.5, 1.0), Thebes.Point3D(2.0, -4.0, 1.0), Thebes.Point3D(2.0, -3.5, 1.0), Thebes.Point3D(3.0, -3.5, 1.0), Thebes.Point3D(3.0, -4.0, 1.0), Thebes.Point3D(2.0, -5.0, 1.0), Thebes.Point3D(1.0, -5.0, 1.0), Thebes.Point3D(0.0, -4.0, 1.0), Thebes.Point3D(0.0, -1.0, 1.0)], Array{Int64,1}[[6, 5, 4, 3], [7, 6, 3, 2], [7, 2, 1], [7, 1, 18, 8], [8, 18, 17, 9], [9, 17, 16, 10], [10, 16, 15], [15, 14, 11, 10], [14, 13, 12, 11], [21, 22, 23, 24], [20, 21, 24, 25], [19, 20, 25], [26, 36, 19, 25], [27, 35, 36, 26], [28, 34, 35, 27], [33, 34, 28], [28, 29, 32, 33], [29, 30, 31, 32], [1, 2, 20, 19], [2, 3, 21, 20], [3, 4, 22, 21], [4, 5, 23, 22], [5, 6, 24, 23], [6, 7, 25, 24], [7, 8, 26, 25], [8, 9, 27, 26], [9, 10, 28, 27], [10, 11, 29, 28], [11, 12, 30, 29], [12, 13, 31, 30], [13, 14, 32, 31], [14, 15, 33, 32], [15, 16, 34, 33], [16, 17, 35, 34], [17, 18, 36, 35], [18, 1, 19, 36]])
letterD = (Thebes.Point3D[Thebes.Point3D(0.0, 0.0, 0.0), Thebes.Point3D(2.0, 0.0, 0.0), Thebes.Point3D(3.0, -1.0, 0.0), Thebes.Point3D(3.0, -4.0, 0.0), Thebes.Point3D(2.0, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(1.0, -1.0, 0.0), Thebes.Point3D(1.5, -1.0, 0.0), Thebes.Point3D(2.0, -1.5, 0.0), Thebes.Point3D(2.0, -3.5, 0.0), Thebes.Point3D(1.5, -4.0, 0.0), Thebes.Point3D(1.0, -4.0, 0.0), Thebes.Point3D(0.0, 0.0, 1.0), Thebes.Point3D(2.0, 0.0, 1.0), Thebes.Point3D(3.0, -1.0, 1.0), Thebes.Point3D(3.0, -4.0, 1.0), Thebes.Point3D(2.0, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0), Thebes.Point3D(1.0, -1.0, 1.0), Thebes.Point3D(1.5, -1.0, 1.0), Thebes.Point3D(2.0, -1.5, 1.0), Thebes.Point3D(2.0, -3.5, 1.0), Thebes.Point3D(1.5, -4.0, 1.0), Thebes.Point3D(1.0, -4.0, 1.0)], Array{Int64,1}[[7, 8, 2, 1], [8, 9, 3, 2], [9, 10, 4, 3], [10, 11, 5, 4], [11, 12, 6, 5], [12, 7, 1, 6], [13, 14, 20, 19], [14, 15, 21, 20], [15, 16, 22, 21], [16, 17, 23, 22], [17, 18, 24, 23], [18, 13, 19, 24], [1, 2, 14, 13], [2, 3, 15, 14], [3, 4, 16, 15], [4, 5, 17, 16], [5, 6, 18, 17], [6, 1, 13, 18], [8, 7, 19, 20], [9, 8, 20, 21], [10, 9, 21, 22], [11, 10, 22, 23], [12, 11, 23, 24], [7, 12, 24, 19]])
letterE = (Thebes.Point3D[Thebes.Point3D(0.0, 0.0, 0.0), Thebes.Point3D(3.0, 0.0, 0.0), Thebes.Point3D(3.0, -1.0, 0.0), Thebes.Point3D(1.0, -1.0, 0.0), Thebes.Point3D(1.0, -2.0, 0.0), Thebes.Point3D(2.0, -2.0, 0.0), Thebes.Point3D(2.0, -3.0, 0.0), Thebes.Point3D(1.0, -3.0, 0.0), Thebes.Point3D(1.0, -4.0, 0.0), Thebes.Point3D(3.0, -4.0, 0.0), Thebes.Point3D(3.0, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(0.0, 0.0, 1.0), Thebes.Point3D(3.0, 0.0, 1.0), Thebes.Point3D(3.0, -1.0, 1.0), Thebes.Point3D(1.0, -1.0, 1.0), Thebes.Point3D(1.0, -2.0, 1.0), Thebes.Point3D(2.0, -2.0, 1.0), Thebes.Point3D(2.0, -3.0, 1.0), Thebes.Point3D(1.0, -3.0, 1.0), Thebes.Point3D(1.0, -4.0, 1.0), Thebes.Point3D(3.0, -4.0, 1.0), Thebes.Point3D(3.0, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0)], Array{Int64,1}[[4, 3, 2, 1], [12, 9, 8, 5, 4, 1], [8, 7, 6, 5], [12, 11, 10, 9], [13, 14, 15, 16], [13, 16, 17, 20, 21, 24], [17, 18, 19, 20], [21, 22, 23, 24], [1, 2, 14, 13], [2, 3, 15, 14], [3, 4, 16, 15], [4, 5, 17, 16], [5, 6, 18, 17], [6, 7, 19, 18], [7, 8, 20, 19], [8, 9, 21, 20], [9, 10, 22, 21], [10, 11, 23, 22], [11, 12, 24, 23], [12, 1, 13, 24]])
letterF = (Thebes.Point3D[Thebes.Point3D(0.0, 0.0, 0.0), Thebes.Point3D(1.0, 0.0, 0.0), Thebes.Point3D(1.0, -2.0, 0.0), Thebes.Point3D(2.0, -2.0, 0.0), Thebes.Point3D(2.0, -3.0, 0.0), Thebes.Point3D(1.0, -3.0, 0.0), Thebes.Point3D(1.0, -4.0, 0.0), Thebes.Point3D(3.0, -4.0, 0.0), Thebes.Point3D(3.0, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(0.0, 0.0, 1.0), Thebes.Point3D(1.0, 0.0, 1.0), Thebes.Point3D(1.0, -2.0, 1.0), Thebes.Point3D(2.0, -2.0, 1.0), Thebes.Point3D(2.0, -3.0, 1.0), Thebes.Point3D(1.0, -3.0, 1.0), Thebes.Point3D(1.0, -4.0, 1.0), Thebes.Point3D(3.0, -4.0, 1.0), Thebes.Point3D(3.0, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0)], Array{Int64,1}[[1, 10, 7, 6, 3, 2], [6, 5, 4, 3], [10, 9, 8, 7], [12, 13, 16, 17, 20, 11], [13, 14, 15, 16], [17, 18, 19, 20], [1, 2, 12, 11], [2, 3, 13, 12], [3, 4, 14, 13], [4, 5, 15, 14], [5, 6, 16, 15], [6, 7, 17, 16], [7, 8, 18, 17], [8, 9, 19, 18], [9, 10, 20, 19], [10, 1, 11, 20]])
letterG = (Thebes.Point3D[Thebes.Point3D(1.0, 0.0, 0.0), Thebes.Point3D(2.0, 0.0, 0.0), Thebes.Point3D(3.0, -1.0, 0.0), Thebes.Point3D(3.0, -2.0, 0.0), Thebes.Point3D(1.5, -2.0, 0.0), Thebes.Point3D(1.5, -1.5, 0.0), Thebes.Point3D(2.0, -1.5, 0.0), Thebes.Point3D(2.0, -1.0, 0.0), Thebes.Point3D(1.5, -0.5, 0.0), Thebes.Point3D(1.0, -1.0, 0.0), Thebes.Point3D(1.0, -4.0, 0.0), Thebes.Point3D(1.5, -4.5, 0.0), Thebes.Point3D(2.0, -4.0, 0.0), Thebes.Point3D(2.0, -3.5, 0.0), Thebes.Point3D(3.0, -3.5, 0.0), Thebes.Point3D(3.0, -4.0, 0.0), Thebes.Point3D(2.0, -5.0, 0.0), Thebes.Point3D(1.0, -5.0, 0.0), Thebes.Point3D(0.0, -4.0, 0.0), Thebes.Point3D(0.0, -1.0, 0.0), Thebes.Point3D(1.0, 0.0, 1.0), Thebes.Point3D(2.0, 0.0, 1.0), Thebes.Point3D(3.0, -1.0, 1.0), Thebes.Point3D(3.0, -2.0, 1.0), Thebes.Point3D(1.5, -2.0, 1.0), Thebes.Point3D(1.5, -1.5, 1.0), Thebes.Point3D(2.0, -1.5, 1.0), Thebes.Point3D(2.0, -1.0, 1.0), Thebes.Point3D(1.5, -0.5, 1.0), Thebes.Point3D(1.0, -1.0, 1.0), Thebes.Point3D(1.0, -4.0, 1.0), Thebes.Point3D(1.5, -4.5, 1.0), Thebes.Point3D(2.0, -4.0, 1.0), Thebes.Point3D(2.0, -3.5, 1.0), Thebes.Point3D(3.0, -3.5, 1.0), Thebes.Point3D(3.0, -4.0, 1.0), Thebes.Point3D(2.0, -5.0, 1.0), Thebes.Point3D(1.0, -5.0, 1.0), Thebes.Point3D(0.0, -4.0, 1.0), Thebes.Point3D(0.0, -1.0, 1.0)], Array{Int64,1}[[7, 6, 5, 4], [8, 7, 4, 3], [9, 8, 3, 2], [1, 9, 2], [10, 9, 1, 20], [11, 10, 20, 19], [12, 11, 19, 18], [12, 18, 17], [13, 12, 17, 16], [16, 15, 14, 13], [24, 25, 26, 27], [23, 24, 27, 28], [22, 23, 28, 29], [22, 29, 21], [40, 21, 29, 30], [39, 40, 30, 31], [38, 39, 31, 32], [37, 38, 32], [36, 37, 32, 33], [33, 34, 35, 36], [1, 2, 22, 21], [2, 3, 23, 22], [3, 4, 24, 23], [4, 5, 25, 24], [5, 6, 26, 25], [6, 7, 27, 26], [7, 8, 28, 27], [8, 9, 29, 28], [9, 10, 30, 29], [10, 11, 31, 30], [11, 12, 32, 31], [12, 13, 33, 32], [13, 14, 34, 33], [14, 15, 35, 34], [15, 16, 36, 35], [16, 17, 37, 36], [17, 18, 38, 37], [18, 19, 39, 38], [19, 20, 40, 39], [20, 1, 21, 40]])
letterH = (Thebes.Point3D[Thebes.Point3D(0.0, 0.0, 0.0), Thebes.Point3D(1.0, 0.0, 0.0), Thebes.Point3D(1.0, -2.0, 0.0), Thebes.Point3D(2.0, -2.0, 0.0), Thebes.Point3D(2.0, 0.0, 0.0), Thebes.Point3D(3.0, 0.0, 0.0), Thebes.Point3D(3.0, -5.0, 0.0), Thebes.Point3D(2.0, -5.0, 0.0), Thebes.Point3D(2.0, -3.0, 0.0), Thebes.Point3D(1.0, -3.0, 0.0), Thebes.Point3D(1.0, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(0.0, 0.0, 1.0), Thebes.Point3D(1.0, 0.0, 1.0), Thebes.Point3D(1.0, -2.0, 1.0), Thebes.Point3D(2.0, -2.0, 1.0), Thebes.Point3D(2.0, 0.0, 1.0), Thebes.Point3D(3.0, 0.0, 1.0), Thebes.Point3D(3.0, -5.0, 1.0), Thebes.Point3D(2.0, -5.0, 1.0), Thebes.Point3D(2.0, -3.0, 1.0), Thebes.Point3D(1.0, -3.0, 1.0), Thebes.Point3D(1.0, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0)], Array{Int64,1}[[12, 11, 10, 3, 2, 1], [10, 9, 4, 3], [9, 8, 7, 6, 5, 4], [13, 14, 15, 22, 23, 24], [15, 16, 21, 22], [16, 17, 18, 19, 20, 21], [1, 2, 14, 13], [2, 3, 15, 14], [3, 4, 16, 15], [4, 5, 17, 16], [5, 6, 18, 17], [6, 7, 19, 18], [7, 8, 20, 19], [8, 9, 21, 20], [9, 10, 22, 21], [10, 11, 23, 22], [11, 12, 24, 23], [12, 1, 13, 24]])
letterI = (Thebes.Point3D[Thebes.Point3D(0.0, 0.0, 0.0), Thebes.Point3D(3.0, 0.0, 0.0), Thebes.Point3D(3.0, -1.0, 0.0), Thebes.Point3D(2.0, -1.0, 0.0), Thebes.Point3D(2.0, -4.0, 0.0), Thebes.Point3D(3.0, -4.0, 0.0), Thebes.Point3D(3.0, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(0.0, -4.0, 0.0), Thebes.Point3D(1.0, -4.0, 0.0), Thebes.Point3D(1.0, -1.0, 0.0), Thebes.Point3D(0.0, -1.0, 0.0), Thebes.Point3D(0.0, 0.0, 1.0), Thebes.Point3D(3.0, 0.0, 1.0), Thebes.Point3D(3.0, -1.0, 1.0), Thebes.Point3D(2.0, -1.0, 1.0), Thebes.Point3D(2.0, -4.0, 1.0), Thebes.Point3D(3.0, -4.0, 1.0), Thebes.Point3D(3.0, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0), Thebes.Point3D(0.0, -4.0, 1.0), Thebes.Point3D(1.0, -4.0, 1.0), Thebes.Point3D(1.0, -1.0, 1.0), Thebes.Point3D(0.0, -1.0, 1.0)], Array{Int64,1}[[12, 11, 4, 3, 2, 1], [11, 10, 5, 4], [10, 9, 8, 7, 6, 5], [13, 14, 15, 16, 23, 24], [16, 17, 22, 23], [17, 18, 19, 20, 21, 22], [1, 2, 14, 13], [2, 3, 15, 14], [3, 4, 16, 15], [4, 5, 17, 16], [5, 6, 18, 17], [6, 7, 19, 18], [7, 8, 20, 19], [8, 9, 21, 20], [9, 10, 22, 21], [10, 11, 23, 22], [11, 12, 24, 23], [12, 1, 13, 24]])
letterJ = (Thebes.Point3D[Thebes.Point3D(1.0, 0.0, 0.0), Thebes.Point3D(1.5, 0.0, 0.0), Thebes.Point3D(2.5, -1.0, 0.0), Thebes.Point3D(2.5, -4.0, 0.0), Thebes.Point3D(3.0, -4.0, 0.0), Thebes.Point3D(3.0, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(0.0, -4.0, 0.0), Thebes.Point3D(1.75, -4.0, 0.0), Thebes.Point3D(1.75, -1.0, 0.0), Thebes.Point3D(1.25, -0.5, 0.0), Thebes.Point3D(0.75, -1.0, 0.0), Thebes.Point3D(0.75, -2.0, 0.0), Thebes.Point3D(0.0, -2.0, 0.0), Thebes.Point3D(0.0, -1.0, 0.0), Thebes.Point3D(1.0, 0.0, 1.0), Thebes.Point3D(1.5, 0.0, 1.0), Thebes.Point3D(2.5, -1.0, 1.0), Thebes.Point3D(2.5, -4.0, 1.0), Thebes.Point3D(3.0, -4.0, 1.0), Thebes.Point3D(3.0, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0), Thebes.Point3D(0.0, -4.0, 1.0), Thebes.Point3D(1.75, -4.0, 1.0), Thebes.Point3D(1.75, -1.0, 1.0), Thebes.Point3D(1.25, -0.5, 1.0), Thebes.Point3D(0.75, -1.0, 1.0), Thebes.Point3D(0.75, -2.0, 1.0), Thebes.Point3D(0.0, -2.0, 1.0), Thebes.Point3D(0.0, -1.0, 1.0)], Array{Int64,1}[[14, 13, 12, 15], [15, 12, 11, 2, 1], [2, 11, 10, 3], [3, 10, 9, 4], [9, 8, 7, 6, 5, 4], [30, 27, 28, 29], [16, 17, 26, 27, 30], [18, 25, 26, 17], [19, 24, 25, 18], [19, 20, 21, 22, 23, 24], [1, 2, 17, 16], [2, 3, 18, 17], [3, 4, 19, 18], [4, 5, 20, 19], [5, 6, 21, 20], [6, 7, 22, 21], [7, 8, 23, 22], [8, 9, 24, 23], [9, 10, 25, 24], [10, 11, 26, 25], [11, 12, 27, 26], [12, 13, 28, 27], [13, 14, 29, 28], [14, 15, 30, 29], [15, 1, 16, 30]])
letterK = (Thebes.Point3D[Thebes.Point3D(0.0, 0.0, 0.0), Thebes.Point3D(1.0, 0.0, 0.0), Thebes.Point3D(1.0, -2.0, 0.0), Thebes.Point3D(2.0, 0.0, 0.0), Thebes.Point3D(3.0, 0.0, 0.0), Thebes.Point3D(1.5, -2.5, 0.0), Thebes.Point3D(3.0, -5.0, 0.0), Thebes.Point3D(2.0, -5.0, 0.0), Thebes.Point3D(1.0, -3.0, 0.0), Thebes.Point3D(1.0, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(0.0, 0.0, 1.0), Thebes.Point3D(1.0, 0.0, 1.0), Thebes.Point3D(1.0, -2.0, 1.0), Thebes.Point3D(2.0, 0.0, 1.0), Thebes.Point3D(3.0, 0.0, 1.0), Thebes.Point3D(1.5, -2.5, 1.0), Thebes.Point3D(3.0, -5.0, 1.0), Thebes.Point3D(2.0, -5.0, 1.0), Thebes.Point3D(1.0, -3.0, 1.0), Thebes.Point3D(1.0, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0)], Array{Int64,1}[[11, 10, 9, 3, 2, 1], [9, 6, 3], [6, 5, 4, 3], [9, 8, 7, 6], [12, 13, 14, 20, 21, 22], [14, 17, 20], [14, 15, 16, 17], [17, 18, 19, 20], [1, 2, 13, 12], [2, 3, 14, 13], [3, 4, 15, 14], [4, 5, 16, 15], [5, 6, 17, 16], [6, 7, 18, 17], [7, 8, 19, 18], [8, 9, 20, 19], [9, 10, 21, 20], [10, 11, 22, 21], [11, 1, 12, 22]])
letterL = (Thebes.Point3D[Thebes.Point3D(0.0, 0.0, 0.0), Thebes.Point3D(3.0, 0.0, 0.0), Thebes.Point3D(3.0, -1.0, 0.0), Thebes.Point3D(1.0, -1.0, 0.0), Thebes.Point3D(1.0, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(0.0, 0.0, 1.0), Thebes.Point3D(3.0, 0.0, 1.0), Thebes.Point3D(3.0, -1.0, 1.0), Thebes.Point3D(1.0, -1.0, 1.0), Thebes.Point3D(1.0, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0)], Array{Int64,1}[[4, 3, 2, 1], [6, 5, 4, 1], [7, 8, 9, 10], [7, 10, 11, 12], [1, 2, 8, 7], [2, 3, 9, 8], [3, 4, 10, 9], [4, 5, 11, 10], [5, 6, 12, 11], [6, 1, 7, 12]])
letterM = (Thebes.Point3D[Thebes.Point3D(0.0, 0.0, 0.0), Thebes.Point3D(0.5, 0.0, 0.0), Thebes.Point3D(1.0, -3.0, 0.0), Thebes.Point3D(1.5, -2.5, 0.0), Thebes.Point3D(2.0, -3.0, 0.0), Thebes.Point3D(2.5, 0.0, 0.0), Thebes.Point3D(3.0, 0.0, 0.0), Thebes.Point3D(2.5, -5.0, 0.0), Thebes.Point3D(1.5, -3.5, 0.0), Thebes.Point3D(0.5, -5.0, 0.0), Thebes.Point3D(0.0, 0.0, 1.0), Thebes.Point3D(0.5, 0.0, 1.0), Thebes.Point3D(1.0, -3.0, 1.0), Thebes.Point3D(1.5, -2.5, 1.0), Thebes.Point3D(2.0, -3.0, 1.0), Thebes.Point3D(2.5, 0.0, 1.0), Thebes.Point3D(3.0, 0.0, 1.0), Thebes.Point3D(2.5, -5.0, 1.0), Thebes.Point3D(1.5, -3.5, 1.0), Thebes.Point3D(0.5, -5.0, 1.0)], Array{Int64,1}[[10, 3, 2, 1], [10, 9, 4, 3], [9, 8, 5, 4], [8, 7, 6, 5], [11, 12, 13, 20], [13, 14, 19, 20], [14, 15, 18, 19], [15, 16, 17, 18], [1, 2, 12, 11], [2, 3, 13, 12], [3, 4, 14, 13], [4, 5, 15, 14], [5, 6, 16, 15], [6, 7, 17, 16], [7, 8, 18, 17], [8, 9, 19, 18], [9, 10, 20, 19], [10, 1, 11, 20]])
letterN = (Thebes.Point3D[Thebes.Point3D(0.0, 0.0, 0.0), Thebes.Point3D(1.0, 0.0, 0.0), Thebes.Point3D(1.0, -2.0, 0.0), Thebes.Point3D(2.0, 0.0, 0.0), Thebes.Point3D(3.0, 0.0, 0.0), Thebes.Point3D(3.0, -5.0, 0.0), Thebes.Point3D(2.0, -5.0, 0.0), Thebes.Point3D(2.0, -3.0, 0.0), Thebes.Point3D(1.0, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(0.0, 0.0, 1.0), Thebes.Point3D(1.0, 0.0, 1.0), Thebes.Point3D(1.0, -2.0, 1.0), Thebes.Point3D(2.0, 0.0, 1.0), Thebes.Point3D(3.0, 0.0, 1.0), Thebes.Point3D(3.0, -5.0, 1.0), Thebes.Point3D(2.0, -5.0, 1.0), Thebes.Point3D(2.0, -3.0, 1.0), Thebes.Point3D(1.0, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0)], Array{Int64,1}[[10, 9, 3, 2, 1], [9, 8, 4, 3], [8, 7, 6, 5, 4], [11, 12, 13, 19, 20], [13, 14, 18, 19], [14, 15, 16, 17, 18], [1, 2, 12, 11], [2, 3, 13, 12], [3, 4, 14, 13], [4, 5, 15, 14], [5, 6, 16, 15], [6, 7, 17, 16], [7, 8, 18, 17], [8, 9, 19, 18], [9, 10, 20, 19], [10, 1, 11, 20]])
letterO = (Thebes.Point3D[Thebes.Point3D(1.0, 0.0, 0.0), Thebes.Point3D(2.0, 0.0, 0.0), Thebes.Point3D(3.0, -1.0, 0.0), Thebes.Point3D(3.0, -4.0, 0.0), Thebes.Point3D(2.0, -5.0, 0.0), Thebes.Point3D(1.0, -5.0, 0.0), Thebes.Point3D(0.0, -4.0, 0.0), Thebes.Point3D(0.0, -1.0, 0.0), Thebes.Point3D(1.0, -1.5, 0.0), Thebes.Point3D(1.5, -1.0, 0.0), Thebes.Point3D(2.0, -1.5, 0.0), Thebes.Point3D(2.0, -3.5, 0.0), Thebes.Point3D(1.5, -4.0, 0.0), Thebes.Point3D(1.0, -3.5, 0.0), Thebes.Point3D(1.0, 0.0, 1.0), Thebes.Point3D(2.0, 0.0, 1.0), Thebes.Point3D(3.0, -1.0, 1.0), Thebes.Point3D(3.0, -4.0, 1.0), Thebes.Point3D(2.0, -5.0, 1.0), Thebes.Point3D(1.0, -5.0, 1.0), Thebes.Point3D(0.0, -4.0, 1.0), Thebes.Point3D(0.0, -1.0, 1.0), Thebes.Point3D(1.0, -1.5, 1.0), Thebes.Point3D(1.5, -1.0, 1.0), Thebes.Point3D(2.0, -1.5, 1.0), Thebes.Point3D(2.0, -3.5, 1.0), Thebes.Point3D(1.5, -4.0, 1.0), Thebes.Point3D(1.0, -3.5, 1.0)], Array{Int64,1}[[10, 2, 1], [10, 11, 3, 2], [11, 12, 4, 3], [12, 13, 5, 4], [13, 6, 5], [13, 14, 7, 6], [14, 9, 8, 7], [9, 10, 1, 8], [15, 16, 24], [16, 17, 25, 24], [17, 18, 26, 25], [18, 19, 27, 26], [19, 20, 27], [20, 21, 28, 27], [21, 22, 23, 28], [22, 15, 24, 23], [1, 2, 16, 15], [2, 3, 17, 16], [3, 4, 18, 17], [4, 5, 19, 18], [5, 6, 20, 19], [6, 7, 21, 20], [7, 8, 22, 21], [8, 1, 15, 22], [10, 9, 23, 24], [11, 10, 24, 25], [12, 11, 25, 26], [13, 12, 26, 27], [14, 13, 27, 28], [9, 14, 28, 23]])
letterP = (Thebes.Point3D[Thebes.Point3D(0.0, 0.0, 0.0), Thebes.Point3D(1.0, 0.0, 0.0), Thebes.Point3D(1.0, -2.0, 0.0), Thebes.Point3D(2.0, -2.0, 0.0), Thebes.Point3D(3.0, -3.0, 0.0), Thebes.Point3D(3.0, -4.0, 0.0), Thebes.Point3D(2.0, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(1.0, -3.0, 0.0), Thebes.Point3D(1.5, -3.0, 0.0), Thebes.Point3D(2.0, -3.5, 0.0), Thebes.Point3D(1.5, -4.0, 0.0), Thebes.Point3D(1.0, -4.0, 0.0), Thebes.Point3D(0.0, 0.0, 1.0), Thebes.Point3D(1.0, 0.0, 1.0), Thebes.Point3D(1.0, -2.0, 1.0), Thebes.Point3D(2.0, -2.0, 1.0), Thebes.Point3D(3.0, -3.0, 1.0), Thebes.Point3D(3.0, -4.0, 1.0), Thebes.Point3D(2.0, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0), Thebes.Point3D(1.0, -3.0, 1.0), Thebes.Point3D(1.5, -3.0, 1.0), Thebes.Point3D(2.0, -3.5, 1.0), Thebes.Point3D(1.5, -4.0, 1.0), Thebes.Point3D(1.0, -4.0, 1.0)], Array{Int64,1}[[8, 13, 9, 3, 2, 1], [9, 10, 4, 3], [10, 11, 5, 4], [11, 6, 5], [11, 12, 7, 6], [8, 7, 12, 13], [14, 15, 16, 22, 26, 21], [16, 17, 23, 22], [17, 18, 24, 23], [18, 19, 24], [19, 20, 25, 24], [26, 25, 20, 21], [1, 2, 15, 14], [2, 3, 16, 15], [3, 4, 17, 16], [4, 5, 18, 17], [5, 6, 19, 18], [6, 7, 20, 19], [7, 8, 21, 20], [8, 1, 14, 21], [9, 13, 26, 22], [13, 12, 25, 26], [12, 11, 24, 25], [11, 10, 23, 24], [10, 9, 22, 23]])
letterQ = (Thebes.Point3D[Thebes.Point3D(1.0, 0.0, 0.0), Thebes.Point3D(1.5, 0.0, 0.0), Thebes.Point3D(2.0, -0.5, 0.0), Thebes.Point3D(2.5, 0.0, 0.0), Thebes.Point3D(3.0, 0.0, 0.0), Thebes.Point3D(3.0, -0.5, 0.0), Thebes.Point3D(2.5, -1.0, 0.0), Thebes.Point3D(3.0, -1.5, 0.0), Thebes.Point3D(3.0, -4.0, 0.0), Thebes.Point3D(2.0, -5.0, 0.0), Thebes.Point3D(1.0, -5.0, 0.0), Thebes.Point3D(0.0, -4.0, 0.0), Thebes.Point3D(0.0, -1.0, 0.0), Thebes.Point3D(1.5, -1.0, 0.0), Thebes.Point3D(2.0, -1.5, 0.0), Thebes.Point3D(2.0, -3.5, 0.0), Thebes.Point3D(1.5, -4.0, 0.0), Thebes.Point3D(1.0, -3.5, 0.0), Thebes.Point3D(1.0, -1.5, 0.0), Thebes.Point3D(1.0, 0.0, 1.0), Thebes.Point3D(1.5, 0.0, 1.0), Thebes.Point3D(2.0, -0.5, 1.0), Thebes.Point3D(2.5, 0.0, 1.0), Thebes.Point3D(3.0, 0.0, 1.0), Thebes.Point3D(3.0, -0.5, 1.0), Thebes.Point3D(2.5, -1.0, 1.0), Thebes.Point3D(3.0, -1.5, 1.0), Thebes.Point3D(3.0, -4.0, 1.0), Thebes.Point3D(2.0, -5.0, 1.0), Thebes.Point3D(1.0, -5.0, 1.0), Thebes.Point3D(0.0, -4.0, 1.0), Thebes.Point3D(0.0, -1.0, 1.0), Thebes.Point3D(1.5, -1.0, 1.0), Thebes.Point3D(2.0, -1.5, 1.0), Thebes.Point3D(2.0, -3.5, 1.0), Thebes.Point3D(1.5, -4.0, 1.0), Thebes.Point3D(1.0, -3.5, 1.0), Thebes.Point3D(1.0, -1.5, 1.0)], Array{Int64,1}[[14, 3, 2, 1], [14, 7, 6, 5, 4, 3], [14, 15, 8, 7], [15, 16, 9, 8], [16, 17, 10, 9], [17, 11, 10], [17, 18, 12, 11], [18, 19, 13, 12], [19, 14, 1, 13], [20, 21, 22, 33], [22, 23, 24, 25, 26, 33], [26, 27, 34, 33], [27, 28, 35, 34], [28, 29, 36, 35], [29, 30, 36], [30, 31, 37, 36], [31, 32, 38, 37], [32, 20, 33, 38], [1, 2, 21, 20], [2, 3, 22, 21], [3, 4, 23, 22], [4, 5, 24, 23], [5, 6, 25, 24], [6, 7, 26, 25], [7, 8, 27, 26], [8, 9, 28, 27], [9, 10, 29, 28], [10, 11, 30, 29], [11, 12, 31, 30], [12, 13, 32, 31], [13, 1, 20, 32], [19, 18, 37, 38], [18, 17, 36, 37], [17, 16, 35, 36], [16, 15, 34, 35], [15, 14, 33, 34], [14, 19, 38, 33]])
letterR = (Thebes.Point3D[Thebes.Point3D(0.0, 0.0, 0.0), Thebes.Point3D(1.0, 0.0, 0.0), Thebes.Point3D(1.0, -1.5, 0.0), Thebes.Point3D(2.0, 0.0, 0.0), Thebes.Point3D(3.0, 0.0, 0.0), Thebes.Point3D(1.5, -2.0, 0.0), Thebes.Point3D(2.0, -2.0, 0.0), Thebes.Point3D(3.0, -3.0, 0.0), Thebes.Point3D(3.0, -4.0, 0.0), Thebes.Point3D(2.0, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(1.0, -3.0, 0.0), Thebes.Point3D(1.5, -3.0, 0.0), Thebes.Point3D(2.0, -3.5, 0.0), Thebes.Point3D(1.5, -4.0, 0.0), Thebes.Point3D(1.0, -4.0, 0.0), Thebes.Point3D(0.0, 0.0, 1.0), Thebes.Point3D(1.0, 0.0, 1.0), Thebes.Point3D(1.0, -1.5, 1.0), Thebes.Point3D(2.0, 0.0, 1.0), Thebes.Point3D(3.0, 0.0, 1.0), Thebes.Point3D(1.5, -2.0, 1.0), Thebes.Point3D(2.0, -2.0, 1.0), Thebes.Point3D(3.0, -3.0, 1.0), Thebes.Point3D(3.0, -4.0, 1.0), Thebes.Point3D(2.0, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0), Thebes.Point3D(1.0, -3.0, 1.0), Thebes.Point3D(1.5, -3.0, 1.0), Thebes.Point3D(2.0, -3.5, 1.0), Thebes.Point3D(1.5, -4.0, 1.0), Thebes.Point3D(1.0, -4.0, 1.0)], Array{Int64,1}[[11, 16, 12, 3, 2, 1], [6, 5, 4, 3], [6, 3, 12], [12, 13, 7, 6], [13, 14, 8, 7], [9, 8, 14], [10, 9, 14, 15], [11, 10, 15, 16], [17, 18, 19, 28, 32, 27], [19, 20, 21, 22], [28, 19, 22], [22, 23, 29, 28], [23, 24, 30, 29], [30, 24, 25], [31, 30, 25, 26], [32, 31, 26, 27], [1, 2, 18, 17], [2, 3, 19, 18], [3, 4, 20, 19], [4, 5, 21, 20], [5, 6, 22, 21], [6, 7, 23, 22], [7, 8, 24, 23], [8, 9, 25, 24], [9, 10, 26, 25], [10, 11, 27, 26], [11, 1, 17, 27], [16, 15, 31, 32], [15, 14, 30, 31], [14, 13, 29, 30], [13, 12, 28, 29], [12, 16, 32, 28]])
letterS = (Thebes.Point3D[Thebes.Point3D(0.0, 0.0, 0.0), Thebes.Point3D(3.0, 0.0, 0.0), Thebes.Point3D(3.0, -3.0, 0.0), Thebes.Point3D(1.0, -3.0, 0.0), Thebes.Point3D(1.0, -4.0, 0.0), Thebes.Point3D(3.0, -4.0, 0.0), Thebes.Point3D(3.0, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(0.0, -2.0, 0.0), Thebes.Point3D(2.0, -2.0, 0.0), Thebes.Point3D(2.0, -1.0, 0.0), Thebes.Point3D(0.0, -1.0, 0.0), Thebes.Point3D(0.0, 0.0, 1.0), Thebes.Point3D(3.0, 0.0, 1.0), Thebes.Point3D(3.0, -3.0, 1.0), Thebes.Point3D(1.0, -3.0, 1.0), Thebes.Point3D(1.0, -4.0, 1.0), Thebes.Point3D(3.0, -4.0, 1.0), Thebes.Point3D(3.0, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0), Thebes.Point3D(0.0, -2.0, 1.0), Thebes.Point3D(2.0, -2.0, 1.0), Thebes.Point3D(2.0, -1.0, 1.0), Thebes.Point3D(0.0, -1.0, 1.0)], Array{Int64,1}[[12, 11, 2, 1], [11, 10, 3, 2], [10, 9, 4, 3], [9, 8, 5, 4], [8, 7, 6, 5], [13, 14, 23, 24], [14, 15, 22, 23], [15, 16, 21, 22], [16, 17, 20, 21], [17, 18, 19, 20], [1, 2, 14, 13], [2, 3, 15, 14], [3, 4, 16, 15], [4, 5, 17, 16], [5, 6, 18, 17], [6, 7, 19, 18], [7, 8, 20, 19], [8, 9, 21, 20], [9, 10, 22, 21], [10, 11, 23, 22], [11, 12, 24, 23], [12, 1, 13, 24]])
letterT = (Thebes.Point3D[Thebes.Point3D(1.0, 0.0, 0.0), Thebes.Point3D(2.0, 0.0, 0.0), Thebes.Point3D(2.0, -4.0, 0.0), Thebes.Point3D(3.0, -4.0, 0.0), Thebes.Point3D(3.0, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(0.0, -4.0, 0.0), Thebes.Point3D(1.0, -4.0, 0.0), Thebes.Point3D(1.0, 0.0, 1.0), Thebes.Point3D(2.0, 0.0, 1.0), Thebes.Point3D(2.0, -4.0, 1.0), Thebes.Point3D(3.0, -4.0, 1.0), Thebes.Point3D(3.0, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0), Thebes.Point3D(0.0, -4.0, 1.0), Thebes.Point3D(1.0, -4.0, 1.0)], Array{Int64,1}[[8, 3, 2, 1], [8, 7, 6, 5, 4, 3], [9, 10, 11, 16], [11, 12, 13, 14, 15, 16], [1, 2, 10, 9], [2, 3, 11, 10], [3, 4, 12, 11], [4, 5, 13, 12], [5, 6, 14, 13], [6, 7, 15, 14], [7, 8, 16, 15], [8, 1, 9, 16]])
letterU = (Thebes.Point3D[Thebes.Point3D(1.0, 0.0, 0.0), Thebes.Point3D(2.0, 0.0, 0.0), Thebes.Point3D(3.0, -1.0, 0.0), Thebes.Point3D(3.0, -5.0, 0.0), Thebes.Point3D(2.0, -5.0, 0.0), Thebes.Point3D(2.0, -1.0, 0.0), Thebes.Point3D(1.0, -1.0, 0.0), Thebes.Point3D(1.0, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(0.0, -1.0, 0.0), Thebes.Point3D(1.0, 0.0, 1.0), Thebes.Point3D(2.0, 0.0, 1.0), Thebes.Point3D(3.0, -1.0, 1.0), Thebes.Point3D(3.0, -5.0, 1.0), Thebes.Point3D(2.0, -5.0, 1.0), Thebes.Point3D(2.0, -1.0, 1.0), Thebes.Point3D(1.0, -1.0, 1.0), Thebes.Point3D(1.0, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0), Thebes.Point3D(0.0, -1.0, 1.0)], Array{Int64,1}[[10, 7, 6, 3, 2, 1], [5, 4, 3, 6], [9, 8, 7, 10], [11, 12, 13, 16, 17, 20], [13, 14, 15, 16], [17, 18, 19, 20], [1, 2, 12, 11], [2, 3, 13, 12], [3, 4, 14, 13], [4, 5, 15, 14], [5, 6, 16, 15], [6, 7, 17, 16], [7, 8, 18, 17], [8, 9, 19, 18], [9, 10, 20, 19], [10, 1, 11, 20]])
letterV = (Thebes.Point3D[Thebes.Point3D(1.5, 0.0, 0.0), Thebes.Point3D(3.0, -2.0, 0.0), Thebes.Point3D(3.0, -5.0, 0.0), Thebes.Point3D(2.0, -5.0, 0.0), Thebes.Point3D(2.0, -2.0, 0.0), Thebes.Point3D(1.5, -1.5, 0.0), Thebes.Point3D(1.0, -2.0, 0.0), Thebes.Point3D(1.0, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(0.0, -2.0, 0.0), Thebes.Point3D(1.5, 0.0, 1.0), Thebes.Point3D(3.0, -2.0, 1.0), Thebes.Point3D(3.0, -5.0, 1.0), Thebes.Point3D(2.0, -5.0, 1.0), Thebes.Point3D(2.0, -2.0, 1.0), Thebes.Point3D(1.5, -1.5, 1.0), Thebes.Point3D(1.0, -2.0, 1.0), Thebes.Point3D(1.0, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0), Thebes.Point3D(0.0, -2.0, 1.0)], Array{Int64,1}[[6, 5, 2, 1], [5, 4, 3, 2], [7, 6, 1, 10], [10, 9, 8, 7], [11, 12, 15, 16], [12, 13, 14, 15], [20, 11, 16, 17], [17, 18, 19, 20], [1, 2, 12, 11], [2, 3, 13, 12], [3, 4, 14, 13], [4, 5, 15, 14], [5, 6, 16, 15], [6, 7, 17, 16], [7, 8, 18, 17], [8, 9, 19, 18], [9, 10, 20, 19], [10, 1, 11, 20]])
letterW = (Thebes.Point3D[Thebes.Point3D(0.5, 0.0, 0.0), Thebes.Point3D(1.5, -1.5, 0.0), Thebes.Point3D(2.5, 0.0, 0.0), Thebes.Point3D(3.0, -5.0, 0.0), Thebes.Point3D(2.5, -5.0, 0.0), Thebes.Point3D(2.0, -2.0, 0.0), Thebes.Point3D(1.5, -2.5, 0.0), Thebes.Point3D(1.0, -2.0, 0.0), Thebes.Point3D(0.5, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(0.5, 0.0, 1.0), Thebes.Point3D(1.5, -1.5, 1.0), Thebes.Point3D(2.5, 0.0, 1.0), Thebes.Point3D(3.0, -5.0, 1.0), Thebes.Point3D(2.5, -5.0, 1.0), Thebes.Point3D(2.0, -2.0, 1.0), Thebes.Point3D(1.5, -2.5, 1.0), Thebes.Point3D(1.0, -2.0, 1.0), Thebes.Point3D(0.5, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0)], Array{Int64,1}[[1, 8, 7, 2], [7, 6, 3, 2], [6, 5, 4, 3], [10, 9, 8, 1], [12, 17, 18, 11], [12, 13, 16, 17], [13, 14, 15, 16], [11, 18, 19, 20], [1, 2, 12, 11], [2, 3, 13, 12], [3, 4, 14, 13], [4, 5, 15, 14], [5, 6, 16, 15], [6, 7, 17, 16], [7, 8, 18, 17], [8, 9, 19, 18], [9, 10, 20, 19], [10, 1, 11, 20]])
letterX = (Thebes.Point3D[Thebes.Point3D(0.0, 0.0, 0.0), Thebes.Point3D(0.5, 0.0, 0.0), Thebes.Point3D(1.5, -2.0, 0.0), Thebes.Point3D(2.5, 0.0, 0.0), Thebes.Point3D(3.0, 0.0, 0.0), Thebes.Point3D(2.0, -2.5, 0.0), Thebes.Point3D(3.0, -5.0, 0.0), Thebes.Point3D(2.5, -5.0, 0.0), Thebes.Point3D(1.5, -3.0, 0.0), Thebes.Point3D(0.5, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(1.0, -2.5, 0.0), Thebes.Point3D(0.0, 0.0, 1.0), Thebes.Point3D(0.5, 0.0, 1.0), Thebes.Point3D(1.5, -2.0, 1.0), Thebes.Point3D(2.5, 0.0, 1.0), Thebes.Point3D(3.0, 0.0, 1.0), Thebes.Point3D(2.0, -2.5, 1.0), Thebes.Point3D(3.0, -5.0, 1.0), Thebes.Point3D(2.5, -5.0, 1.0), Thebes.Point3D(1.5, -3.0, 1.0), Thebes.Point3D(0.5, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0), Thebes.Point3D(1.0, -2.5, 1.0)], Array{Int64,1}[[3, 2, 1, 12], [6, 5, 4, 3], [6, 3, 12, 9], [9, 8, 7, 6], [12, 11, 10, 9], [24, 13, 14, 15], [15, 16, 17, 18], [21, 24, 15, 18], [18, 19, 20, 21], [21, 22, 23, 24], [1, 2, 14, 13], [2, 3, 15, 14], [3, 4, 16, 15], [4, 5, 17, 16], [5, 6, 18, 17], [6, 7, 19, 18], [7, 8, 20, 19], [8, 9, 21, 20], [9, 10, 22, 21], [10, 11, 23, 22], [11, 12, 24, 23], [12, 1, 13, 24]])
letterY = (Thebes.Point3D[Thebes.Point3D(1.25, 0.0, 0.0), Thebes.Point3D(1.75, 0.0, 0.0), Thebes.Point3D(1.75, -2.75, 0.0), Thebes.Point3D(3.0, -5.0, 0.0), Thebes.Point3D(2.5, -5.0, 0.0), Thebes.Point3D(1.5, -3.5, 0.0), Thebes.Point3D(0.5, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(1.25, -2.75, 0.0), Thebes.Point3D(1.25, 0.0, 1.0), Thebes.Point3D(1.75, 0.0, 1.0), Thebes.Point3D(1.75, -2.75, 1.0), Thebes.Point3D(3.0, -5.0, 1.0), Thebes.Point3D(2.5, -5.0, 1.0), Thebes.Point3D(1.5, -3.5, 1.0), Thebes.Point3D(0.5, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0), Thebes.Point3D(1.25, -2.75, 1.0)], Array{Int64,1}[[9, 6, 3, 2, 1], [6, 5, 4, 3], [9, 8, 7, 6], [10, 11, 12, 15, 18], [12, 13, 14, 15], [15, 16, 17, 18], [1, 2, 11, 10], [2, 3, 12, 11], [3, 4, 13, 12], [4, 5, 14, 13], [5, 6, 15, 14], [6, 7, 16, 15], [7, 8, 17, 16], [8, 9, 18, 17], [9, 1, 10, 18]])
letterZ = (Thebes.Point3D[Thebes.Point3D(0.0, 0.0, 0.0), Thebes.Point3D(3.0, 0.0, 0.0), Thebes.Point3D(3.0, -1.0, 0.0), Thebes.Point3D(1.0, -1.0, 0.0), Thebes.Point3D(3.0, -4.0, 0.0), Thebes.Point3D(3.0, -5.0, 0.0), Thebes.Point3D(0.0, -5.0, 0.0), Thebes.Point3D(0.0, -4.0, 0.0), Thebes.Point3D(2.0, -4.0, 0.0), Thebes.Point3D(0.0, -1.0, 0.0), Thebes.Point3D(0.0, 0.0, 1.0), Thebes.Point3D(3.0, 0.0, 1.0), Thebes.Point3D(3.0, -1.0, 1.0), Thebes.Point3D(1.0, -1.0, 1.0), Thebes.Point3D(3.0, -4.0, 1.0), Thebes.Point3D(3.0, -5.0, 1.0), Thebes.Point3D(0.0, -5.0, 1.0), Thebes.Point3D(0.0, -4.0, 1.0), Thebes.Point3D(2.0, -4.0, 1.0), Thebes.Point3D(0.0, -1.0, 1.0)], Array{Int64,1}[[1, 4, 3, 2], [1, 10, 9, 6, 5, 4], [9, 8, 7, 6], [12, 13, 14, 11], [14, 15, 16, 19, 20, 11], [16, 17, 18, 19], [1, 2, 12, 11], [2, 3, 13, 12], [3, 4, 14, 13], [4, 5, 15, 14], [5, 6, 16, 15], [6, 7, 17, 16], [7, 8, 18, 17], [8, 9, 19, 18], [9, 10, 20, 19], [10, 1, 11, 20]])


text:
#############################
# Resampling uncertain values
#############################
include("uncertain_values/test_resampling_certain_value.jl")
include("uncertain_values/test_resampling_uncertainvalues.jl")
include("uncertain_values/test_resampling_uncertainvalues_kde.jl")
include("uncertain_values/test_resampling_UncertainScalarPopulation.jl")

ext: jl


text:
using Thebes, Luxor

include(dirname(pathof(Thebes)) * "/../data/moreobjects.jl")

platonics = [:boxtorus, :concave, :cone, :crossshape, :cube, :cuboctahedron, :dodecahedron , :geodesic, :helix2,
:icosahedron, :icosidodecahedron, :octahedron, :octtorus, :rhombicosidodecahedron,
:rhombicuboctahedron, :rhombitruncated_cubeoctahedron, :rhombitruncated_icosidodecahedron,
:snub_cube, :snub_dodecahedron, :sphere2, :tet3d, :tetrahedron, :triangle, :truncated_cube,
:truncated_dodecahedron, :truncated_icosahedron, :truncated_octahedron, :truncated_tetrahedron]

function main()
    Drawing(800, 800, "cubes.svg")
    origin()
    background("orange")
    #setopacity(0.75)
    o = platonics[rand(1:end)]

    @layer begin
        sethue("magenta")
        carpet(300)
    end
    axes3D()
    setline(10)
    for x in -100:50:100
        for y in -100:50:100
            object = make(moreobjects[2])
            scaleby!(object, 10, 10, 10)
            moveby!(object, x, y, 50rand())
            rotateby!(object, object.vertices[1], rand(), rand(), rand())
            pin(object)

            # dark version
            moveby!(object, 0, 0, -object.vertices[1].z)
            scaleby!(object, 1, 1, 0)
            pin(object)
        end
    end
    finish()
end

main()


text:
using VCFTools
using MendelImpute
using GeneticVariation

function run(data::String, width::Int)
    tgtfile = "./tgt_masked_unphased_" * data
    reffile = "./ref_" * data
    outfile = "./mendel_imputed_" * data
    phase(tgtfile, reffile, outfile = outfile, width = width, fast_method=true)
end

data  = ARGS[1]
width = ARGS[2]
run(data, width)

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
function torn_system_jacobian_sparsity(sys)
    s = structure(sys)
    @unpack fullvars, graph, partitions = s

    # The sparsity pattern of `nlsolve(f, u, p)` w.r.t `p` is difficult to
    # determine in general. Consider the "simplest" case, a linear system. We
    # have
    #                   A u = p.
    # Clearly, the sparsity of `u` depends on the sparsity of both `p` and `A`
    # in a non-trivial way. However, in the generic case, `u` is dense even when
    # `A` and `p` are sparse. For instance
    #
    # ```julia
    # julia> using Random, SparseArrays
    #
    # julia> A = sprand(MersenneTwister(1234), 100, 100, 0.1);
    #
    # julia> p = sprand(MersenneTwister(12345), 100, 0.05);
    #
    # julia> count(x->abs(x) < 1e-5, A \ Vector(p))
    # 0
    # ```
    #
    # Let 𝑇 be the set of tearing variables and 𝑉 be the set of all *states* in
    # the residual equations. In the following code, we are going to assume the
    # connection between 𝑇 (the `u` in from above) and 𝑉 ∖ 𝑇 (the `p` in from
    # above) has full incidence.
    #
    # Note that as we are reducing algebraic equations numerically, it could be
    # the case that a later partition (a BLT block) contains tearing variables
    # from other partitions.
    #
    # We know that partitions are BLT ordered. Hence, the tearing variables in
    # each partition is unique, and all states in a partition must be
    # either differential variables or algebraic tearing variables that are
    # from previous partitions. Hence, we can build the dependency chain as we
    # traverse the partitions.

    # `avars2dvars` maps a algebraic variable to its differential variable
    # dependencies.
    avars2dvars = Dict{Int,Set{Int}}()
    c = 0
    for partition in partitions
        @unpack e_residual, v_residual = partition
        # initialization
        for tvar in v_residual
            avars2dvars[tvar] = Set{Int}()
        end
        for teq in e_residual
            c += 1
            for var in 𝑠neighbors(graph, teq)
                # Skip the tearing variables in the current partition, because
                # we are computing them from all the other states.
                Graphs.insorted(var, v_residual) && continue
                deps = get(avars2dvars, var, nothing)
                if deps === nothing # differential variable
                    @assert !isalgvar(s, var)
                    for tvar in v_residual
                        push!(avars2dvars[tvar], var)
                    end
                else # tearing variable from previous partitions
                    @assert isalgvar(s, var)
                    for tvar in v_residual
                        union!(avars2dvars[tvar], avars2dvars[var])
                    end
                end
            end
        end
    end

    dvrange = diffvars_range(s)
    dvar2idx = Dict(v=>i for (i, v) in enumerate(dvrange))
    I = Int[]; J = Int[]
    eqidx = 0
    for ieq in 𝑠vertices(graph)
        isalgeq(s, ieq) && continue
        eqidx += 1
        for ivar in 𝑠neighbors(graph, ieq)
            if isdiffvar(s, ivar)
                push!(I, eqidx)
                push!(J, dvar2idx[ivar])
            elseif isalgvar(s, ivar)
                for dvar in avars2dvars[ivar]
                    push!(I, eqidx)
                    push!(J, dvar2idx[dvar])
                end
            end
        end
    end
    sparse(I, J, true)
end

"""
    partitions_dag(s::SystemStructure)

Return a DAG (sparse matrix) of partitions to use for parallelism.
"""
function partitions_dag(s::SystemStructure)
    @unpack partitions, graph = s

    # `partvars[i]` contains all the states that appear in `partitions[i]`
    partvars = map(partitions) do partition
        ipartvars = Set{Int}()
        for req in partition.e_residual
            union!(ipartvars, 𝑠neighbors(graph, req))
        end
        ipartvars
    end

    I, J = Int[], Int[]
    n = length(partitions)
    for (i, partition) in enumerate(partitions)
        for j in i+1:n
            # The only way for a later partition `j` to depend on an earlier
            # partition `i` is when `partvars[j]` contains one of tearing
            # variables of partition `i`.
            if !isdisjoint(partvars[j], partition.v_residual)
                # j depends on i
                push!(I, i)
                push!(J, j)
            end
        end
    end

    sparse(I, J, true, n, n)
end

function gen_nlsolve(sys, eqs, vars; checkbounds=true)
    @assert !isempty(vars)
    @assert length(eqs) == length(vars)
    rhss = map(x->x.rhs, eqs)
    # We use `vars` instead of `graph` to capture parameters, too.
    allvars = unique(collect(Iterators.flatten(map(ModelingToolkit.vars, rhss))))
    params = setdiff(allvars, vars)

    u0map = defaults(sys)
    # splatting to tighten the type
    u0 = [map(var->get(u0map, var, 1e-3), vars)...]
    # specialize on the scalar case
    isscalar = length(u0) == 1
    u0 = isscalar ? u0[1] : SVector(u0...)

    fname = gensym("fun")
    f = Func(
        [
         DestructuredArgs(vars, inbounds=!checkbounds)
         DestructuredArgs(params, inbounds=!checkbounds)
        ],
        [],
        isscalar ? rhss[1] : MakeArray(rhss, SVector)
    ) |> SymbolicUtils.Code.toexpr

    solver_call = LiteralExpr(quote
            $numerical_nlsolve(
                               $fname,
                               # initial guess
                               $u0,
                               # "captured variables"
                               ($(params...),)
                              )
        end)

    [
     fname ← @RuntimeGeneratedFunction(f)
     DestructuredArgs(vars, inbounds=!checkbounds) ← solver_call
    ]
end

function get_torn_eqs_vars(sys; checkbounds=true)
    s = structure(sys)
    partitions = s.partitions
    vars = s.fullvars
    eqs = equations(sys)

    torn_eqs  = map(idxs-> eqs[idxs], map(x->x.e_residual, partitions))
    torn_vars = map(idxs->vars[idxs], map(x->x.v_residual, partitions))

    gen_nlsolve.((sys,), torn_eqs, torn_vars, checkbounds=checkbounds)
end

function build_torn_function(
        sys;
        expression=false,
        jacobian_sparsity=true,
        checkbounds=false,
        kw...
    )

    rhss = []
    for eq in equations(sys)
        isdiffeq(eq) && push!(rhss, eq.rhs)
    end

    out = Sym{Any}(gensym("out"))
    odefunbody = SetArray(
        !checkbounds,
        out,
        rhss
    )

    s = structure(sys)
    states = map(i->s.fullvars[i], diffvars_range(s))
    syms = map(Symbol, states)
    pre = get_postprocess_fbody(sys)

    expr = SymbolicUtils.Code.toexpr(
        Func(
             [
              out
              DestructuredArgs(states, inbounds=!checkbounds)
              DestructuredArgs(parameters(sys), inbounds=!checkbounds)
              independent_variables(sys)
             ],
             [],
             pre(Let(
                 collect(Iterators.flatten(get_torn_eqs_vars(sys, checkbounds=checkbounds))),
                 odefunbody
                ))
            )
    )
    if expression
        expr
    else
        observedfun = let sys = sys, dict = Dict()
            function generated_observed(obsvar, u, p, t)
                obs = get!(dict, value(obsvar)) do
                    build_observed_function(sys, obsvar, checkbounds=checkbounds)
                end
                obs(u, p, t)
            end
        end

        ODEFunction{true}(
                          @RuntimeGeneratedFunction(expr),
                          sparsity = torn_system_jacobian_sparsity(sys),
                          syms = syms,
                          observed = observedfun,
                         )
    end
end

"""
    find_solve_sequence(partitions, vars)

given a set of `vars`, find the groups of equations we need to solve for
to obtain the solution to `vars`
"""
function find_solve_sequence(partitions, vars)
    subset = filter(x -> !isdisjoint(x.v_residual, vars), partitions)
    isempty(subset) && return []
    vars′ = mapreduce(x->x.v_residual, union, subset)
    if vars′ == vars
        return subset
    else
        return find_solve_sequence(partitions, vars′)
    end
end

function build_observed_function(
        sys, ts;
        expression=false,
        output_type=Array,
        checkbounds=true
    )

    if (isscalar = !(ts isa AbstractVector))
        ts = [ts]
    end
    ts = Symbolics.scalarize.(value.(ts))

    vars = Set()
    foreach(Base.Fix1(vars!, vars), ts)
    ivs = independent_variables(sys)
    dep_vars = collect(setdiff(vars, ivs))

    s = structure(sys)
    @unpack partitions, fullvars, graph = s
    diffvars = map(i->fullvars[i], diffvars_range(s))
    algvars = map(i->fullvars[i], algvars_range(s))

    required_algvars = Set(intersect(algvars, vars))
    obs = observed(sys)
    observed_idx = Dict(map(x->x.lhs, obs) .=> 1:length(obs))
    # FIXME: this is a rather rough estimate of dependencies.
    maxidx = 0
    sts = Set(states(sys))
    for (i, s) in enumerate(dep_vars)
        idx = get(observed_idx, s, nothing)
        if idx === nothing
            if !(s in sts)
                throw(ArgumentError("$s is either an observed nor a state variable."))
            end
            continue
        end
        idx > maxidx && (maxidx = idx)
    end
    vs = Set()
    for idx in 1:maxidx
        vars!(vs, obs[idx].rhs)
        union!(required_algvars, intersect(algvars, vs))
        empty!(vs)
    end

    varidxs = findall(x->x in required_algvars, fullvars)
    subset = find_solve_sequence(partitions, varidxs)
    if !isempty(subset)
        eqs = equations(sys)

        torn_eqs  = map(idxs-> eqs[idxs.e_residual], subset)
        torn_vars = map(idxs->fullvars[idxs.v_residual], subset)

        solves = gen_nlsolve.((sys,), torn_eqs, torn_vars; checkbounds=checkbounds)
    else
        solves = []
    end

    subs = []
    for sym in vars
        eqidx = get(observed_idx, sym, nothing)
        eqidx === nothing && continue
        push!(subs, sym ← obs[eqidx].rhs)
    end
    pre = get_postprocess_fbody(sys)

    ex = Func(
        [
         DestructuredArgs(diffvars, inbounds=!checkbounds)
         DestructuredArgs(parameters(sys), inbounds=!checkbounds)
         independent_variables(sys)
        ],
        [],
        pre(Let(
            [
             collect(Iterators.flatten(solves))
             map(eq -> eq.lhs←eq.rhs, obs[1:maxidx])
             subs
            ],
            isscalar ? ts[1] : MakeArray(ts, output_type)
           ))
    ) |> Code.toexpr

    expression ? ex : @RuntimeGeneratedFunction(ex)
end

struct ODAEProblem{iip}
end

ODAEProblem(args...; kw...) = ODAEProblem{true}(args...; kw...)

"""
    ODAEProblem{iip}(sys, u0map, tspan, parammap = DiffEqBase.NullParameters(); kw...)

This constructor acts similar to the one for [`ODEProblem`](@ref) with the following changes:
`ODESystem`s can sometimes be further reduced if `structural_simplify` has
already been applied to them. This is done this constructor.
In these cases, the constructor uses the knowledge of the strongly connected
components calculated during the process of simplification as the basis for
building pre-simplified nonlinear systems in the implicit solving.
In summary: these problems are structurally modified, but could be
more efficient and more stable. Note, the returned object is still of type [`ODEProblem`](@ref).
"""
function ODAEProblem{iip}(
                          sys,
                          u0map,
                          tspan,
                          parammap=DiffEqBase.NullParameters();
                          kw...
                         ) where {iip}
    s = structure(sys)
    @unpack fullvars = s
    dvs = map(i->fullvars[i], diffvars_range(s))
    ps = parameters(sys)
    defs = defaults(sys)

    u0 = ModelingToolkit.varmap_to_vars(u0map, dvs; defaults=defs)
    p = ModelingToolkit.varmap_to_vars(parammap, ps; defaults=defs)

    ODEProblem{iip}(build_torn_function(sys; kw...), u0, tspan, p; kw...)
end


text:
"""
    struct Convex{T1<:AbstractArray,T2<:Union{Flux.Zeros, AbstractVector},T3}
        W::T1
        U::T1
        b::T2
        σ::T3
    end

Input Convex Neural Network (ICNN) layer by Amos et al.
"""
struct Convex{T1<:AbstractArray,T2<:Union{Flux.Zeros,AbstractVector},T3} <: AbstractLayer
    W::T1
    U::T1
    b::T2
    σ::T3
end

function Convex(
    z_in::T,
    x_in::T,
    out::T,
    activation = identity::Function;
    fw = randn::Function,
    fb = zeros::Function,
    precision = Float32,
) where {T<:Integer}
    return Convex(
        fw(precision, out, z_in),
        fw(precision, out, x_in),
        fb(precision, out),
        activation,
    )
end

function (m::Convex)(x::AbstractArray)
    W, b, σ = m.W, m.b, m.σ
    sz = size(x)
    x = reshape(x, sz[1], :) # reshape to handle dims > 1 as batch dimensions
    x = σ.(W * x .+ b)

    return reshape(x, :, sz[2:end]...)
end

function (m::Convex)(z::AbstractArray, x::AbstractArray)
    W, U, b, σ = m.W, m.U, m.b, m.σ
    sz = size(z)
    sx = size(x)
    z = reshape(z, sz[1], :)
    x = reshape(x, sx[1], :)
    z = σ.(softplus.(W) * z + U * x .+ b)

    return reshape(z, :, sz[2:end]...)
end

Flux.@functor Convex

function Base.show(io::IO, model::Convex{T1,T2,T3}) where {T1,T2,T3}
    print(
        io,
        "Input convex layer{$T1,$T2,$T3}\n",
        "nonnegative weights for: $(model.W |> size)\n",
        "input weights: $(model.U |> size)\n",
        "bias: $(model.b |> size)\n",
        "activation: $(model.σ)\n",
    )
end

struct ICNN{T} <: AbstractChain
    layers::T
end

function ICNN(
    din::TI,
    dout::TI,
    layer_sizes::TT,
    activation = identity::Function;
    fw = randn::Function,
    fb = zeros::Function,
    precision = Float32,
) where {TI<:Integer,TT<:Union{Tuple,AbstractVector}}

    layers = (Dense(din, layer_sizes[1], activation),)

    if length(layer_sizes) > 1
        i = 1
        for out in layer_sizes[2:end]
            layers = (
                layers...,
                Convex(
                    layer_sizes[i],
                    din,
                    out,
                    activation;
                    fw = fw,
                    fb = fb,
                    precision = precision,
                ),
            )
            i += 1
        end
    end
    layers = (
        layers...,
        Convex(
            layer_sizes[end],
            din,
            dout,
            identity;
            fw = fw,
            fb = fb,
            precision = precision,
        ),
    )

    return ICNN(layers)

end

(m::ICNN)(x::AbstractArray) = begin
    z = m.layers[1](x)
    for i = 2:length(m.layers)
        z = m.layers[i](z, x)
    end
    return z
end

Flux.@functor ICNN

function Base.show(io::IO, model::ICNN{T}) where {T}
    print(
        io,
        "Input convex neural network{$T}\n",
        "Layers: 1 Dense layer + $(length(model.layers)-1) convex layers\n",
    )
end

"""
    struct FastConvex{I<:Integer,F1,F2} <: DiffEqFlux.FastLayer
        zin::I
        xin::I
        out::I
        σ::F1
        initial_params::F2
    end

Fast ICNN layer
"""
struct FastConvex{I<:Integer,F1,F2} <: DiffEqFlux.FastLayer
    zin::I
    xin::I
    out::I
    σ::F1
    initial_params::F2
end

function FastConvex(
    zin::Integer,
    xin::Integer,
    out::Integer,
    σ = identity;
    fw = randn,
    fb = zeros,
    precision = Float32,
)
    initial_params() =
        vcat(vec(fw(precision, out, zin)), vec(fw(precision, out, xin)), fb(precision, out))
    return FastConvex{typeof(out),typeof(σ),typeof(initial_params)}(
        zin,
        xin,
        out,
        σ,
        initial_params,
    )
end

function (f::FastConvex)(z, x, p)
    f.σ.(
        softplus.(reshape(p[1:(f.out*f.zin)], f.out, f.zin)) * z .+
        reshape(p[f.out*f.zin+1:(f.out*f.zin+f.out*f.xin)], f.out, f.xin) * x .+
        p[(f.out*f.zin+f.out*f.xin+1):end],
    )
end

DiffEqFlux.paramlength(f::FastConvex) = f.out * (f.zin + f.xin + 1)
DiffEqFlux.initial_params(f::FastConvex) = f.initial_params()

struct FastICNN{T} <: DiffEqFlux.FastLayer
    layers::T
end

function FastICNN(
    din::TI,
    dout::TI,
    layer_sizes::TT,
    activation = identity::Function;
    fw = randn::Function,
    fb = zeros::Function,
    precision = Float32,
) where {TI<:Integer,TT<:Union{Tuple,AbstractVector}}

    layers = (FastDense(din, layer_sizes[1], activation),)

    if length(layer_sizes) > 1
        i = 1
        for out in layer_sizes[2:end]
            layers = (
                layers...,
                FastConvex(
                    layer_sizes[i],
                    din,
                    out,
                    activation;
                    fw = fw,
                    fb = fb,
                    precision = precision,
                ),
            )
            i += 1
        end
    end
    layers = (
        layers...,
        FastConvex(
            layer_sizes[end],
            din,
            dout,
            identity;
            fw = fw,
            fb = fb,
            precision = precision,
        ),
    )

    return FastICNN(layers)

end

DiffEqFlux.initial_params(c::FastICNN) = vcat(initial_params.(c.layers)...)

function (m::FastICNN)(x::AbstractArray, p)
    z = m.layers[1](x, p[1:DiffEqFlux.paramlength(m.layers[1])])
    counter = DiffEqFlux.paramlength(m.layers[1])
    for i = 2:length(m.layers)
        z = m.layers[i](z, x, p[counter+1:counter+DiffEqFlux.paramlength(m.layers[i])])
        counter += DiffEqFlux.paramlength(m.layers[i])
    end

    return z
end


text:
@testset "policies" begin

    @testset "RandomPolicy" begin
        p = RandomPolicy(DiscreteSpace(3))
        obs = (reward = 0.0, terminal = false, state = 1)

        Random.seed!(p, 321)
        actions = [p(obs) for _ in 1:100]
        Random.seed!(p, 321)
        new_actions = [p(obs) for _ in 1:100]
        @test actions == new_actions
    end

    @testset "WeightedRandomPolicy" begin
        @testset "1D" begin
            weight = [1, 2, 3]
            ratio = weight ./ sum(weight)
            N = 1000
            actions = [:a, :b, :c]
            p = WeightedRandomPolicy(weight, actions = actions, seed = 123)

            samples = [p(nothing, MINIMAL_ACTION_SET) for _ in 1:N]
            stats = countmap(samples)
            for (a, r) in zip(actions, ratio)
                @test isapprox(r, stats[a] / N, atol = 0.05)
            end

            legal_actions = [:aa, :cc]
            legal_actions_mask = [true, false, true]
            obs = (
                reward = 0.0,
                terminal = false,
                state = nothing,
                legal_actions = legal_actions,
                legal_actions_mask = legal_actions_mask,
            )
            samples = [p(obs) for _ in 1:N]
            stats = countmap(samples)

            weighted_ratio = ratio[legal_actions_mask] ./ sum(ratio[legal_actions_mask])
            for i in 1:length(legal_actions)
                @test isapprox(stats[legal_actions[i]] / N, weighted_ratio[i], atol = 0.05)
            end
        end

        @testset "2D" begin
            n_state = 2
            weight = reshape(1:6, 3, n_state)
            ratio = weight ./ sum(weight, dims = 1)
            N = 1000
            actions = [:a, :b, :c]

            p = WeightedRandomPolicy(weight, actions = actions, seed = 123)

            for state in 1:n_state
                samples = [p((state = state,), MINIMAL_ACTION_SET) for _ in 1:N]
                stats = countmap(samples)
                for (a, r) in zip(actions, ratio[:, state])
                    @test isapprox(r, stats[a] / N, atol = 0.05)
                end
            end

            for state in 1:n_state
                legal_actions = [:aa, :cc]
                legal_actions_mask = [true, false, true]
                obs = (
                    reward = 0.0,
                    terminal = false,
                    state = state,
                    legal_actions = legal_actions,
                    legal_actions_mask = legal_actions_mask,
                )
                samples = [p(obs) for _ in 1:N]
                stats = countmap(samples)

                weighted_ratio =
                    ratio[legal_actions_mask, state] ./
                    sum(ratio[legal_actions_mask, state])
                for i in 1:length(legal_actions)
                    @test isapprox(
                        stats[legal_actions[i]] / N,
                        weighted_ratio[i],
                        atol = 0.05,
                    )
                end
            end
        end
    end
end


text:
export ParsingSystem, Move

abstract type ParsingSystem end;
abstract type Move end;

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
module AdditiveRungeKuttaMethod
export AdditiveRungeKutta
export ARK2GiraldoKellyConstantinescu
export ARK548L2SA2KennedyCarpenter, ARK437L2SA1KennedyCarpenter

# Naive formulation that uses equation 3.8 from Giraldo, Kelly, and Constantinescu (2013) directly.
# Seems to cut the number of solver iterations by half but requires Nstages - 1 additional storage.
struct NaiveVariant end
additional_storage(::NaiveVariant, Q, Nstages) = (Lstages = ntuple(i -> similar(Q), Nstages),)

# Formulation that does things exactly as in Giraldo, Kelly, and Constantinescu (2013).
# Uses only one additional vector of storage regardless of the number of stages.
struct LowStorageVariant end
additional_storage(::LowStorageVariant, Q, Nstages) = (Qtt = similar(Q),)

using GPUifyLoops
include("AdditiveRungeKuttaMethod_kernels.jl")

using StaticArrays

using ..ODESolvers
ODEs = ODESolvers
using ..SpaceMethods
using ..LinearSolvers
using ..MPIStateArrays: device, realview

"""
    op! = EulerOperator(f!, ϵ)

Construct a linear operator which performs an explicit Euler step ``Q + α f(Q)``,
where `f!` and `op!` both operate inplace, with extra arguments passed through, i.e.
```
op!(LQ, Q, args...)
```
is equivalent to
```
f!(dQ, Q, args...)
LQ .= Q .+ ϵ .* dQ
```
"""
mutable struct EulerOperator{F,FT}
  f!::F
  ϵ::FT
end

function (op::EulerOperator)(LQ, Q, args...)
  op.f!(LQ, Q, args..., increment=false)
  @. LQ = Q + op.ϵ * LQ
end

"""
    AdditiveRungeKutta(f, l, linearsolver, RKAe, RKAi, RKB, RKC, Q;
                       split_nonlinear_linear, variant, dt, t0 = 0)

This is a time stepping object for implicit-explicit time stepping of a
decomposed differential equation. When `split_nonlinear_linear == false`
the equation is assumed to be decomposed as

```math
  \\dot{Q} = [l(Q, t)] + [f(Q, t) - l(Q, t)]
```

where `Q` is the state, `f` is the full tendency and
`l` is the chosen linear operator. When `split_nonlinear_linear == true`
the assumed decomposition is

```math
  \\dot{Q} = [l(Q, t)] + [f(Q, t)]
```

where `f` is now only the nonlinear tendency. For both decompositions the
linear operator `l` is integrated implicitly whereas the remaining part
is integrated explicitly. Other arguments are the required time step size `dt`
and the optional initial time `t0`. The resulting linear systems are solved
using the provided `linearsolver` solver. This time stepping object is intended
to be passed to the `solve!` command.

The constructor builds an additive Runge--Kutta scheme
based on the provided `RKAe`, `RKAi`, `RKB` and `RKC` coefficient arrays.
Additionally `variant` specifies which of the analytically equivalent but numerically
different formulations of the scheme is used.

The available concrete implementations are:

  - [`ARK2GiraldoKellyConstantinescu`](@ref)
  - [`ARK548L2SA2KennedyCarpenter`](@ref)
  - [`ARK437L2SA1KennedyCarpenter`](@ref)
"""
mutable struct AdditiveRungeKutta{T, RT, AT, LT, V, VS, Nstages, Nstages_sq} <: ODEs.AbstractODESolver
  "time step"
  dt::RT
  "time"
  t::RT
  "rhs function"
  rhs!
  "rhs linear operator"
  rhs_linear!
  "implicit operator, pre-factorized"
  implicitoperator!
  "linear solver"
  linearsolver::LT
  "Storage for solution during the AdditiveRungeKutta update"
  Qstages::NTuple{Nstages, AT}
  "Storage for RHS during the AdditiveRungeKutta update"
  Rstages::NTuple{Nstages, AT}
  "Storage for the linear solver rhs vector"
  Qhat::AT
  "RK coefficient matrix A for the explicit scheme"
  RKA_explicit::SArray{NTuple{2, Nstages}, RT, 2, Nstages_sq}
  "RK coefficient matrix A for the implicit scheme"
  RKA_implicit::SArray{NTuple{2, Nstages}, RT, 2, Nstages_sq}
  "RK coefficient vector B (rhs add in scaling)"
  RKB::SArray{Tuple{Nstages}, RT, 1, Nstages}
  "RK coefficient vector C (time scaling)"
  RKC::SArray{Tuple{Nstages}, RT, 1, Nstages}
  split_nonlinear_linear::Bool
  "Variant of the ARK scheme"
  variant::V
  "Storage dependent on the variant of the ARK scheme"
  variant_storage::VS

  function AdditiveRungeKutta(rhs!,
                              rhs_linear!,
                              linearsolver::AbstractLinearSolver,
                              RKA_explicit, RKA_implicit, RKB, RKC,
                              split_nonlinear_linear, variant,
                              Q::AT; dt=nothing, t0=0) where {AT<:AbstractArray}

    @assert dt != nothing

    T = eltype(Q)
    LT = typeof(linearsolver)
    RT = real(T)

    Nstages = length(RKB)

    Qstages = (Q, ntuple(i -> similar(Q), Nstages - 1)...)
    Rstages = ntuple(i -> similar(Q), Nstages)
    Qhat = similar(Q)

    V = typeof(variant)
    variant_storage = additional_storage(variant, Q, Nstages)
    VS = typeof(variant_storage)

    # The code throughout assumes SDIRK implicit tableau so we assert that
    # here.
    for is = 2:Nstages
      @assert RKA_implicit[is, is] ≈ RKA_implicit[2, 2]
    end

    α = dt * RKA_implicit[2, 2]
    # Here we are passing NaN for the time since prefactorization assumes the
    # operator is time independent.  If that is not the case the NaN will
    # surface.
    implicitoperator! = prefactorize(EulerOperator(rhs_linear!, -α),
                                     linearsolver, Q, nothing, T(NaN))

    new{T, RT, AT, LT, V, VS, Nstages, Nstages ^ 2}(RT(dt), RT(t0),
                                                    rhs!, rhs_linear!, implicitoperator!, linearsolver,
                                                    Qstages, Rstages, Qhat,
                                                    RKA_explicit, RKA_implicit, RKB, RKC,
                                                    split_nonlinear_linear,
                                                    variant, variant_storage)
  end
end

function AdditiveRungeKutta(spacedisc::AbstractSpaceMethod,
                            spacedisc_linear::AbstractSpaceMethod,
                            linearsolver::AbstractLinearSolver,
                            RKA_explicit, RKA_implicit, RKB, RKC,
                            split_nonlinear_linear, variant,
                            Q::AT; dt=nothing, t0=0) where {AT<:AbstractArray}
  rhs! = (x...; increment) -> SpaceMethods.odefun!(spacedisc, x..., increment = increment)
  rhs_linear! = (x...; increment) -> SpaceMethods.odefun!(spacedisc_linear, x..., increment = increment)
  AdditiveRungeKutta(rhs!, rhs_linear!, linearsolver,
                     RKA_explicit, RKA_implicit, RKB, RKC,
                     split_nonlinear_linear, variant,
                     Q; dt=dt, t0=t0)
end

# this will only work for iterative solves
# direct solvers use prefactorization
ODEs.isadjustable(ark::AdditiveRungeKutta) = ark.implicitoperator! isa EulerOperator
function ODEs.updatedt!(ark::AdditiveRungeKutta, dt)
  @assert ODEs.isadjustable(ark)
  ark.dt = dt
  α = dt * ark.RKA_implicit[2, 2]
  ark.implicitoperator! = EulerOperator(ark.rhs_linear!, -α)
end
ODEs.updatetime!(ark::AdditiveRungeKutta, time) = (ark.t = time)

function ODEs.dostep!(Q, ark::AdditiveRungeKutta, p, timeend::Real,
                      adjustfinalstep::Bool)
  time, dt = ark.t, ark.dt
  if adjustfinalstep && time + dt > timeend
    dt = timeend - time
  end
  @assert dt > 0

  ODEs.dostep!(Q, ark, p, time, dt)

  if dt == ark.dt
    ark.t += dt
  else
    ark.t = timeend
  end

end

function ODEs.dostep!(Q, ark::AdditiveRungeKutta, p, time::Real, dt::Real,
                      slow_δ = nothing, slow_rv_dQ = nothing,
                      slow_scaling = nothing)
  ODEs.dostep!(Q, ark, ark.variant, p, time, dt, slow_δ, slow_rv_dQ, slow_scaling)
end

function ODEs.dostep!(Q, ark::AdditiveRungeKutta, variant::NaiveVariant,
                      p, time::Real, dt::Real,
                      slow_δ = nothing, slow_rv_dQ = nothing,
                      slow_scaling = nothing)
  implicitoperator!, linearsolver = ark.implicitoperator!, ark.linearsolver
  RKA_explicit, RKA_implicit = ark.RKA_explicit, ark.RKA_implicit
  RKB, RKC = ark.RKB, ark.RKC
  rhs!, rhs_linear! = ark.rhs!, ark.rhs_linear!
  Qstages, Rstages = ark.Qstages, ark.Rstages
  Qhat = ark.Qhat
  split_nonlinear_linear = ark.split_nonlinear_linear
  Lstages = ark.variant_storage.Lstages

  rv_Q = realview(Q)
  rv_Qstages = realview.(Qstages)
  rv_Lstages = realview.(Lstages)
  rv_Rstages = realview.(Rstages)
  rv_Qhat = realview(Qhat)

  Nstages = length(RKB)

  threads = 256
  blocks = div(length(rv_Q) + threads - 1, threads)

  # calculate the rhs at first stage to initialize the stage loop
  rhs!(Rstages[1], Qstages[1], p, time + RKC[1] * dt, increment = false)

  if dt != ark.dt
    α = dt * RKA_implicit[2, 2]
    implicitoperator! = EulerOperator(rhs_linear!, -α)
  end

  rhs_linear!(Lstages[1], Qstages[1], p, time + RKC[1] * dt, increment = false)

  # note that it is important that this loop does not modify Q!
  for istage = 2:Nstages
    stagetime = time + RKC[istage] * dt

    # this kernel also initializes Qstages[istage] with an initial guess
    # for the linear solver
    @launch(device(Q), threads = threads, blocks = blocks,
            stage_update!(variant, rv_Q, rv_Qstages, rv_Lstages, rv_Rstages, rv_Qhat,
                          RKA_explicit, RKA_implicit, dt, Val(istage),
                          Val(split_nonlinear_linear), slow_δ, slow_rv_dQ))

    #solves Q_tt = Qhat + dt * RKA_implicit[istage, istage] * rhs_linear!(Q_tt)
    α = dt * RKA_implicit[istage, istage]
    linearoperator! = function(LQ, Q)
      rhs_linear!(LQ, Q, p, stagetime; increment = false)
      @. LQ = Q - α * LQ
    end
    linearsolve!(implicitoperator!, linearsolver, Qstages[istage], Qhat, p, stagetime)

    rhs!(Rstages[istage], Qstages[istage], p, stagetime, increment = false)
    rhs_linear!(Lstages[istage], Qstages[istage], p, stagetime, increment = false)
  end

  # compose the final solution
  @launch(device(Q), threads = threads, blocks = blocks,
          solution_update!(variant, rv_Q, rv_Lstages, rv_Rstages, RKB, dt,
                           Val(Nstages), Val(split_nonlinear_linear),
                           slow_δ, slow_rv_dQ, slow_scaling))
end

function ODEs.dostep!(Q, ark::AdditiveRungeKutta, variant::LowStorageVariant,
                      p, time::Real, dt::Real,
                      slow_δ = nothing, slow_rv_dQ = nothing,
                      slow_scaling = nothing)
  implicitoperator!, linearsolver = ark.implicitoperator!, ark.linearsolver
  RKA_explicit, RKA_implicit = ark.RKA_explicit, ark.RKA_implicit
  RKB, RKC = ark.RKB, ark.RKC
  rhs!, rhs_linear! = ark.rhs!, ark.rhs_linear!
  Qstages, Rstages = ark.Qstages, ark.Rstages
  Qhat = ark.Qhat
  split_nonlinear_linear = ark.split_nonlinear_linear
  Qtt = ark.variant_storage.Qtt

  rv_Q = realview(Q)
  rv_Qstages = realview.(Qstages)
  rv_Rstages = realview.(Rstages)
  rv_Qhat = realview(Qhat)
  rv_Qtt = realview(Qtt)

  Nstages = length(RKB)

  threads = 256
  blocks = div(length(rv_Q) + threads - 1, threads)

  # calculate the rhs at first stage to initialize the stage loop
  rhs!(Rstages[1], Qstages[1], p, time + RKC[1] * dt, increment = false)

  if dt != ark.dt
    α = dt * RKA_implicit[2, 2]
    implicitoperator! = EulerOperator(rhs_linear!, -α)
  end

  # note that it is important that this loop does not modify Q!
  for istage = 2:Nstages
    stagetime = time + RKC[istage] * dt

    # this kernel also initializes Qtt for the linear solver
    @launch(device(Q), threads = threads, blocks = blocks,
            stage_update!(variant, rv_Q, rv_Qstages, rv_Rstages, rv_Qhat, rv_Qtt,
                          RKA_explicit, RKA_implicit, dt, Val(istage),
                          Val(split_nonlinear_linear), slow_δ, slow_rv_dQ))

    #solves Q_tt = Qhat + dt * RKA_implicit[istage, istage] * rhs_linear!(Q_tt)
    linearsolve!(implicitoperator!, linearsolver, Qtt, Qhat, p, stagetime)

    #update Qstages
    Qstages[istage] .+= Qtt

    rhs!(Rstages[istage], Qstages[istage], p, stagetime, increment = false)
  end

  if split_nonlinear_linear
    for istage = 1:Nstages
      stagetime = time + RKC[istage] * dt
      rhs_linear!(Rstages[istage], Qstages[istage], p, stagetime, increment = true)
    end
  end

  # compose the final solution
  @launch(device(Q), threads = threads, blocks = blocks,
          solution_update!(variant, rv_Q, rv_Rstages, RKB, dt, Val(Nstages),
                           slow_δ, slow_rv_dQ, slow_scaling))
end

"""
    ARK2GiraldoKellyConstantinescu(f, l, linearsolver, Q; dt, t0,
                                   split_nonlinear_linear, variant, paperversion)

This function returns an [`AdditiveRungeKutta`](@ref) time stepping object,
see the documentation of [`AdditiveRungeKutta`](@ref) for arguments definitions.
This time stepping object is intended to be passed to the `solve!` command.

`paperversion=true` uses the coefficients from the paper, `paperversion=false`
uses coefficients that make the scheme (much) more stable but less accurate

This uses the second-order-accurate 3-stage additive Runge--Kutta scheme of
Giraldo, Kelly and Constantinescu (2013).

### References
    @article{giraldo2013implicit,
      title={Implicit-explicit formulations of a three-dimensional nonhydrostatic unified model of the atmosphere ({NUMA})},
      author={Giraldo, Francis X and Kelly, James F and Constantinescu, Emil M},
      journal={SIAM Journal on Scientific Computing},
      volume={35},
      number={5},
      pages={B1162--B1194},
      year={2013},
      publisher={SIAM}
    }
"""
function ARK2GiraldoKellyConstantinescu(F, L,
                                        linearsolver::AbstractLinearSolver,
                                        Q::AT; dt=nothing, t0=0,
                                        split_nonlinear_linear=false,
                                        variant=LowStorageVariant(),
                                        paperversion=false) where {AT<:AbstractArray}

  @assert dt != nothing

  T = eltype(Q)
  RT = real(T)

  a32 = RT(paperversion ? (3 + 2sqrt(2)) / 6 : 1 // 2)
  RKA_explicit = [RT(0)           RT(0)   RT(0);
                  RT(2 - sqrt(2)) RT(0)   RT(0);
                  RT(1 - a32)     RT(a32) RT(0)]

  RKA_implicit = [RT(0)               RT(0)               RT(0);
                  RT(1 - 1 / sqrt(2)) RT(1 - 1 / sqrt(2)) RT(0);
                  RT(1 / (2sqrt(2)))  RT(1 / (2sqrt(2)))  RT(1 - 1 / sqrt(2))]

  RKB = [RT(1 / (2sqrt(2))), RT(1 / (2sqrt(2))), RT(1 - 1 / sqrt(2))]
  RKC = [RT(0), RT(2 - sqrt(2)), RT(1)]

  Nstages = length(RKB)

  AdditiveRungeKutta(F, L, linearsolver,
                     RKA_explicit, RKA_implicit, RKB, RKC,
                     split_nonlinear_linear,
                     variant,
                     Q; dt=dt, t0=t0)
end

"""
    ARK548L2SA2KennedyCarpenter(f, l, linearsolver, Q; dt, t0,
                                split_nonlinear_linear, variant)

This function returns an [`AdditiveRungeKutta`](@ref) time stepping object,
see the documentation of [`AdditiveRungeKutta`](@ref) for arguments definitions.
This time stepping object is intended to be passed to the `solve!` command.

This uses the fifth-order-accurate 8-stage additive Runge--Kutta scheme of
Kennedy and Carpenter (2013).

### References

    @article{kennedy2019higher,
      title={Higher-order additive Runge--Kutta schemes for ordinary differential equations},
      author={Kennedy, Christopher A and Carpenter, Mark H},
      journal={Applied Numerical Mathematics},
      volume={136},
      pages={183--205},
      year={2019},
      publisher={Elsevier}
    }
"""
function ARK548L2SA2KennedyCarpenter(F, L,
                                     linearsolver::AbstractLinearSolver,
                                     Q::AT; dt=nothing, t0=0,
                                     split_nonlinear_linear=false,
                                     variant=LowStorageVariant()) where {AT<:AbstractArray}

  @assert dt != nothing

  T = eltype(Q)
  RT = real(T)

  Nstages = 8
  gamma = RT(2 // 9)

  # declared as Arrays for mutability, later these will be converted to static arrays
  RKA_explicit = zeros(RT, Nstages, Nstages)
  RKA_implicit = zeros(RT, Nstages, Nstages)
  RKB = zeros(RT, Nstages)
  RKC = zeros(RT, Nstages)

  # the main diagonal
  for is = 2:Nstages
    RKA_implicit[is, is] = gamma
  end

  RKA_implicit[3, 2] = RT(2366667076620 //  8822750406821)
  RKA_implicit[4, 2] = RT(-257962897183 // 4451812247028)
  RKA_implicit[4, 3] = RT(128530224461 // 14379561246022)
  RKA_implicit[5, 2] = RT(-486229321650 // 11227943450093)
  RKA_implicit[5, 3] = RT(-225633144460 // 6633558740617)
  RKA_implicit[5, 4] = RT(1741320951451 // 6824444397158)
  RKA_implicit[6, 2] = RT(621307788657 // 4714163060173)
  RKA_implicit[6, 3] = RT(-125196015625 // 3866852212004)
  RKA_implicit[6, 4] = RT(940440206406 // 7593089888465)
  RKA_implicit[6, 5] = RT(961109811699 // 6734810228204)
  RKA_implicit[7, 2] = RT(2036305566805 // 6583108094622)
  RKA_implicit[7, 3] = RT(-3039402635899 // 4450598839912)
  RKA_implicit[7, 4] = RT(-1829510709469 // 31102090912115)
  RKA_implicit[7, 5] = RT(-286320471013 // 6931253422520)
  RKA_implicit[7, 6] = RT(8651533662697 // 9642993110008)

  RKA_explicit[3, 1] = RT(1 // 9)
  RKA_explicit[3, 2] = RT(1183333538310 // 1827251437969)
  RKA_explicit[4, 1] = RT(895379019517 // 9750411845327)
  RKA_explicit[4, 2] = RT(477606656805 // 13473228687314)
  RKA_explicit[4, 3] = RT(-112564739183 // 9373365219272)
  RKA_explicit[5, 1] = RT(-4458043123994 // 13015289567637)
  RKA_explicit[5, 2] = RT(-2500665203865 // 9342069639922)
  RKA_explicit[5, 3] = RT(983347055801 // 8893519644487)
  RKA_explicit[5, 4] = RT(2185051477207 // 2551468980502)
  RKA_explicit[6, 1] = RT(-167316361917 // 17121522574472)
  RKA_explicit[6, 2] = RT(1605541814917 // 7619724128744)
  RKA_explicit[6, 3] = RT(991021770328 // 13052792161721)
  RKA_explicit[6, 4] = RT(2342280609577 // 11279663441611)
  RKA_explicit[6, 5] = RT(3012424348531 // 12792462456678)
  RKA_explicit[7, 1] = RT(6680998715867 // 14310383562358)
  RKA_explicit[7, 2] = RT(5029118570809 // 3897454228471)
  RKA_explicit[7, 3] = RT(2415062538259 // 6382199904604)
  RKA_explicit[7, 4] = RT(-3924368632305 // 6964820224454)
  RKA_explicit[7, 5] = RT(-4331110370267 // 15021686902756)
  RKA_explicit[7, 6] = RT(-3944303808049 // 11994238218192)
  RKA_explicit[8, 1] = RT(2193717860234 // 3570523412979)
  RKA_explicit[8, 2] = RKA_explicit[8, 1]
  RKA_explicit[8, 3] = RT(5952760925747 // 18750164281544)
  RKA_explicit[8, 4] = RT(-4412967128996 // 6196664114337)
  RKA_explicit[8, 5] = RT(4151782504231 // 36106512998704)
  RKA_explicit[8, 6] = RT(572599549169 // 6265429158920)
  RKA_explicit[8, 7] = RT(-457874356192 // 11306498036315)

  RKB[2] = 0
  RKB[3] = RT(3517720773327 // 20256071687669)
  RKB[4] = RT(4569610470461 // 17934693873752)
  RKB[5] = RT(2819471173109 // 11655438449929)
  RKB[6] = RT(3296210113763 // 10722700128969)
  RKB[7] = RT(-1142099968913 // 5710983926999)

  RKC[2] = RT(4 // 9)
  RKC[3] = RT(6456083330201 // 8509243623797)
  RKC[4] = RT(1632083962415 // 14158861528103)
  RKC[5] = RT(6365430648612 // 17842476412687)
  RKC[6] = RT(18 // 25)
  RKC[7] = RT(191 // 200)

  for is = 2:Nstages
    RKA_implicit[is, 1] = RKA_implicit[is, 2]
  end

  for is = 1:Nstages-1
    RKA_implicit[Nstages, is] = RKB[is]
  end

  RKB[1] = RKB[2]
  RKB[8] = gamma

  RKA_explicit[2, 1] = RKC[2]
  RKA_explicit[Nstages, 1] = RKA_explicit[Nstages, 2]

  RKC[1] = 0
  RKC[Nstages] = 1

  # conversion to static arrays
  RKA_explicit = SMatrix{Nstages, Nstages}(RKA_explicit)
  RKA_implicit = SMatrix{Nstages, Nstages}(RKA_implicit)
  RKB = SVector{Nstages}(RKB)
  RKC = SVector{Nstages}(RKC)

  ark = AdditiveRungeKutta(F, L, linearsolver,
                           RKA_explicit, RKA_implicit, RKB, RKC,
                           split_nonlinear_linear,
                           variant,
                           Q; dt=dt, t0=t0)
end

"""
    ARK437L2SA1KennedyCarpenter(f, l, linearsolver, Q; dt, t0,
                                split_nonlinear_linear, variant)

This function returns an [`AdditiveRungeKutta`](@ref) time stepping object,
see the documentation of [`AdditiveRungeKutta`](@ref) for arguments definitions.
This time stepping object is intended to be passed to the `solve!` command.

This uses the fourth-order-accurate 7-stage additive Runge--Kutta scheme of
Kennedy and Carpenter (2013).

### References
    @article{kennedy2019higher,
      title={Higher-order additive Runge--Kutta schemes for ordinary differential equations},
      author={Kennedy, Christopher A and Carpenter, Mark H},
      journal={Applied Numerical Mathematics},
      volume={136},
      pages={183--205},
      year={2019},
      publisher={Elsevier}
    }
"""
function ARK437L2SA1KennedyCarpenter(F, L,
                                     linearsolver::AbstractLinearSolver,
                                     Q::AT; dt=nothing, t0=0,
                                     split_nonlinear_linear=false,
                                     variant=LowStorageVariant()) where {AT<:AbstractArray}

  @assert dt != nothing

  T = eltype(Q)
  RT = real(T)

  Nstages = 7
  gamma = RT(1235 // 10000)

  # declared as Arrays for mutability, later these will be converted to static arrays
  RKA_explicit = zeros(RT, Nstages, Nstages)
  RKA_implicit = zeros(RT, Nstages, Nstages)
  RKB = zeros(RT, Nstages)
  RKC = zeros(RT, Nstages)

  # the main diagonal
  for is = 2:Nstages
    RKA_implicit[is, is] = gamma
  end

  RKA_implicit[3, 2] = RT(624185399699 // 4186980696204)
  RKA_implicit[4, 2] = RT(1258591069120 // 10082082980243)
  RKA_implicit[4, 3] = RT(-322722984531 // 8455138723562)
  RKA_implicit[5, 2] = RT(-436103496990 // 5971407786587)
  RKA_implicit[5, 3] = RT(-2689175662187 // 11046760208243)
  RKA_implicit[5, 4] = RT(4431412449334 // 12995360898505)
  RKA_implicit[6, 2] = RT(-2207373168298 // 14430576638973)
  RKA_implicit[6, 3] = RT(242511121179 // 3358618340039)
  RKA_implicit[6, 4] = RT(3145666661981 // 7780404714551)
  RKA_implicit[6, 5] = RT(5882073923981 // 14490790706663)
  RKA_implicit[7, 2] = 0
  RKA_implicit[7, 3] = RT(9164257142617 // 17756377923965)
  RKA_implicit[7, 4] = RT(-10812980402763 // 74029279521829)
  RKA_implicit[7, 5] = RT(1335994250573 // 5691609445217)
  RKA_implicit[7, 6] = RT(2273837961795 // 8368240463276)

  RKA_explicit[3, 1] = RT(247 // 4000)
  RKA_explicit[3, 2] = RT(2694949928731 // 7487940209513)
  RKA_explicit[4, 1] = RT(464650059369 // 8764239774964)
  RKA_explicit[4, 2] = RT(878889893998 // 2444806327765)
  RKA_explicit[4, 3] = RT(-952945855348 // 12294611323341)
  RKA_explicit[5, 1] = RT(476636172619 // 8159180917465)
  RKA_explicit[5, 2] = RT(-1271469283451 // 7793814740893)
  RKA_explicit[5, 3] = RT(-859560642026 // 4356155882851)
  RKA_explicit[5, 4] = RT(1723805262919 // 4571918432560)
  RKA_explicit[6, 1] = RT(6338158500785 // 11769362343261)
  RKA_explicit[6, 2] = RT(-4970555480458 // 10924838743837)
  RKA_explicit[6, 3] = RT(3326578051521 // 2647936831840)
  RKA_explicit[6, 4] = RT(-880713585975 // 1841400956686)
  RKA_explicit[6, 5] = RT(-1428733748635 // 8843423958496)
  RKA_explicit[7, 2] = RT(760814592956 // 3276306540349)
  RKA_explicit[7, 3] = RT(-47223648122716 // 6934462133451)
  RKA_explicit[7, 4] = RT(71187472546993 // 9669769126921)
  RKA_explicit[7, 5] = RT(-13330509492149 // 9695768672337)
  RKA_explicit[7, 6] = RT(11565764226357 // 8513123442827)

  RKB[2] = 0
  RKB[3] = RT(9164257142617 // 17756377923965)
  RKB[4] = RT(-10812980402763 // 74029279521829)
  RKB[5] = RT(1335994250573 // 5691609445217)
  RKB[6] = RT(2273837961795 // 8368240463276)
  RKB[7] = RT(247 // 2000)

  RKC[2] = RT(247 // 2000)
  RKC[3] = RT(4276536705230 // 10142255878289)
  RKC[4] = RT(67 // 200)
  RKC[5] = RT(3 // 40)
  RKC[6] = RT(7 // 10)

  for is = 2:Nstages
    RKA_implicit[is, 1] = RKA_implicit[is, 2]
  end

  for is = 1:Nstages-1
    RKA_implicit[Nstages, is] = RKB[is]
  end

  RKB[1] = RKB[2]

  RKA_explicit[2, 1] = RKC[2]
  RKA_explicit[Nstages, 1] = RKA_explicit[Nstages, 2]

  RKC[1] = 0
  RKC[Nstages] = 1

  # conversion to static arrays
  RKA_explicit = SMatrix{Nstages, Nstages}(RKA_explicit)
  RKA_implicit = SMatrix{Nstages, Nstages}(RKA_implicit)
  RKB = SVector{Nstages}(RKB)
  RKC = SVector{Nstages}(RKC)

  ark = AdditiveRungeKutta(F, L, linearsolver,
                           RKA_explicit, RKA_implicit, RKB, RKC,
                           split_nonlinear_linear,
                           variant,
                           Q; dt=dt, t0=t0)
end

end

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
# Copyright (c) 2012: John Myles White and other contributors.
#
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation
# files (the "Software"), to deal in the Software without
# restriction, including without limitation the rights to use,
# copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
# OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
# BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

immutable LevenbergMarquardt <: Optim.Optimizer end

"""
    `levenberg_marquardt(f!, g!, initial_x, m; <keyword arguments>`

Returns the argmin over x of `sum(f(x).^2)` using the Levenberg-Marquardt
algorithm, and an estimate of the Jacobian of `f` at x.

The function `f!(x, storage)` should take an input `x` vector of length n and return an output
vector of length m. The function `g!(x, storage)` is the Jacobian of f, and should return an
m x n Jacobian matrix. `initial_x` is an initial guess for the solution.

Implements box constraints as described in Kanzow, Yamashita, Fukushima (2004; J
Comp & Applied Math).

# Keyword arguments
* `tolX::Real=1e-8`: search tolerance in x
* `tolG::Real=1e-12`: search tolerance in gradient
* `maxIter::Integer=100`: maximum number of iterations
* `min_step_quality=1e-3`: for steps below this quality, the trust region is shrinked
* `good_step_quality=0.75`: for steps above this quality, the trust region is expanded
* `lambda::Real=10.0`: (inverse of) initial trust region radius
* `lambda_increase=10.`: `lambda` is multiplied by this factor after step below min quality
* `lambda_decrease=0.1`: `lambda` is multiplied by this factor after good quality steps
* `show_trace::Bool=false`: print a status summary on each iteration if true
* `lower,upper=[]`: bound solution to these limits
"""
function levenberg_marquardt{T}(f!::Function, g!::Function, initial_x::AbstractVector{T},
                                m::Int; tolX::Real = 1e-8, tolG::Real = 1e-12,
                                maxIter::Integer = 100, lambda::Real = 10.0,
                                lambda_increase::Real = 10., lambda_decrease::Real = 0.1,
                                min_step_quality::Real = 1e-3, good_step_quality::Real = 0.75,
                                show_trace::Bool = false, lower::Vector{T} = Array{T}(0),
                                upper::Vector{T} = Array{T}(0))

    # check parameters
    ((isempty(lower) || length(lower)==length(initial_x)) && (isempty(upper) || length(upper)==length(initial_x))) ||
            throw(ArgumentError("Bounds must either be empty or of the same length as the number of parameters."))
    ((isempty(lower) || all(initial_x .>= lower)) && (isempty(upper) || all(initial_x .<= upper))) ||
            throw(ArgumentError("Initial guess must be within bounds."))
    (0 <= min_step_quality < 1) || throw(ArgumentError(" 0 <= min_step_quality < 1 must hold."))
    (0 < good_step_quality <= 1) || throw(ArgumentError(" 0 < good_step_quality <= 1 must hold."))
    (min_step_quality < good_step_quality) || throw(ArgumentError("min_step_quality < good_step_quality must hold."))

    # other constants
    const MAX_LAMBDA = 1e16 # minimum trust region radius
    const MIN_LAMBDA = 1e-16 # maximum trust region radius
    const MIN_DIAGONAL = 1e-6 # lower bound on values of diagonal matrix used to regularize the trust region step

    converged = false
    x_converged = false
    g_converged = false
    need_jacobian = true
    iterCt = 0
    x = copy(initial_x)
    delta_x = copy(initial_x)
    f_calls = 0
    g_calls = 0

    # Create buffers
    n = length(x)
    fcur = Vector{T}(m)
    trial_f = Vector{T}(m)
    JJ = Matrix{T}(n, n)
    J = Matrix{T}(m, n)
    n_buffer = Vector{T}(n)

    f!(x, fcur)
    f_calls += 1
    residual = sum(abs2, fcur)

    # Maintain a trace of the system.
    tr = Optim.OptimizationTrace{LevenbergMarquardt}()
    if show_trace
        d = Dict("lambda" => lambda)
        os = Optim.OptimizationState{LevenbergMarquardt}(iterCt, sum(abs2, fcur), NaN, d)
        push!(tr, os)
        println(os)
    end

    while (~converged && iterCt < maxIter)
        if need_jacobian
            g!(x, J)
            g_calls += 1
            need_jacobian = false
        end
        # we want to solve:
        #    argmin 0.5*||J(x)*delta_x + f(x)||^2 + lambda*||diagm(J'*J)*delta_x||^2
        # Solving for the minimum gives:
        #    (J'*J + lambda*diagm(DtD)) * delta_x == -J' * f(x), where DtD = sum(abs2, J,1)
        # Where we have used the equivalence: diagm(J'*J) = diagm(sum(abs2, J,1))
        # It is additionally useful to bound the elements of DtD below to help
        # prevent "parameter evaporation".
        DtD = vec(sum(abs2, J, 1))
        for i in 1:length(DtD)
            if DtD[i] <= MIN_DIAGONAL
                DtD[i] = MIN_DIAGONAL
            end
        end
        # delta_x = ( J'*J + lambda * Diagonal(DtD) ) \ ( -J'*fcur )
        At_mul_B!(JJ, J, J)
        @simd for i in 1:n
            @inbounds JJ[i, i] += lambda * DtD[i]
        end
        At_mul_B!(n_buffer, J, fcur)
        scale!(n_buffer, -1)
        delta_x = JJ \ n_buffer

        # apply box constraints
        if !isempty(lower)
            @simd for i in 1:n
               @inbounds delta_x[i] = max(x[i] + delta_x[i], lower[i]) - x[i]
            end
        end
        if !isempty(upper)
            @simd for i in 1:n
               @inbounds delta_x[i] = min(x[i] + delta_x[i], upper[i]) - x[i]
            end
        end

        # if the linear assumption is valid, our new residual should be:
        predicted_residual = sum(abs2, J*delta_x + fcur)

        # try the step and compute its quality
        f!(x + delta_x, trial_f)
        f_calls += 1
        trial_residual = sum(abs2, trial_f)
        # step quality = residual change / predicted residual change
        rho = (trial_residual - residual) / (predicted_residual - residual)
        if rho > min_step_quality
            x += delta_x
            fcur .= trial_f
            residual = trial_residual
            if rho > good_step_quality
                # increase trust region radius
                lambda = max(lambda_decrease*lambda, MIN_LAMBDA)
            end
            need_jacobian = true
        else
            # decrease trust region radius
            lambda = min(lambda_increase*lambda, MAX_LAMBDA)
        end
        iterCt += 1

        # show state
        if show_trace
            g_norm = norm(J' * fcur, Inf)
            d = Dict("g(x)" => g_norm, "dx" => delta_x, "lambda" => lambda)
            os = Optim.OptimizationState{LevenbergMarquardt}(iterCt, sum(abs2, fcur), g_norm, d)
            push!(tr, os)
            println(os)
        end

        # check convergence criteria:
        # 1. Small gradient: norm(J^T * fcur, Inf) < tolG
        # 2. Small step size: norm(delta_x) < tolX
        if norm(J' * fcur, Inf) < tolG
            g_converged = true
        elseif norm(delta_x) < tolX*(tolX + norm(x))
            x_converged = true
        end
        converged = g_converged | x_converged
    end

    Optim.MultivariateOptimizationResults(
        LevenbergMarquardt(),  # method
        initial_x,             # initial_x
        x,                     # minimizer
        sum(abs2, fcur),       # minimum
        iterCt,                # iterations
        !converged,            # iteration_converged
        x_converged,           # x_converged
        0.0,                   # x_tol
        0.0,
        false,                 # f_converged
        0.0,                   # f_tol
        0.0,
        g_converged,           # g_converged
        tolG,                  # g_tol
        0.0,
        false,                 # f_increased
        tr,                    # trace
        f_calls,               # f_calls
        g_calls,               # g_calls
        0                      # h_calls
    )
end


text:
"""
    beta_rng(α, β, shape=1; seed=nothing)

Generate a `shape` element array of random variables from a Beta(`α`, `β`) distribution. Optionally you can set a specific seed.

# Notes

The Beta distribution is given:

``f(x,α,β) = \\frac{x^{α-1 (1-x)^{β-1}}}{Β(α,β)} \\quad x \\in [0,1]``

where ``Β(α,β) = \\frac{Γ(α)Γ(β)}{Γ(α+β)}``

# Examples

```julia-repl
julia> beta_rng(1,2)
1-element Vector{Float64}:
 0.44456674672905633

julia> beta_rng(1, 2, (2,2))
2×2 Matrix{Float64}:
 0.0792132  0.595657
 0.737615   0.649721

```

# References

D.P. Kroese, T. Taimre, Z.I. Botev. Handbook of Monte Carlo Methods.
  Wiley Series in Probability and Statistics, John Wiley & Sons, New York, 2011.

Law, A. Simulation modeling and analysis, 5th Ed. McGraw Hill Education, Tuscon, 2013.
"""
function beta_rng(α::Real, β::Real, shape::Union{Int, Tuple{Vararg{Int}}}=1; seed::Union{Int, Nothing}=nothing)
    seed_setter(seed)
    Y₁ = gamma_rng(α, 1, shape, seed=seed)
    Y₂ = gamma_rng(β, 1, shape, seed=seed)
    X = Y₁ ./ (Y₁ .+ Y₂)
    return X
end


text:
# Bi-objective linear assignment problem (bilap)
#
# Example 9.38 (from Ulungu and Teghem, 1994), page 255 of
# Multicriteria Optimization (2nd edt), M. Ehrgott, Springer 2005.

# ---- Packages to use
using vOptGeneric, JuMP, GLPK

# ---- Values of the instance to solve
C1 = [  5  1  4  7 ;  # coefficients's vector of the objective 1
        6  2  2  6 ;
        2  8  4  4 ;
        3  5  7  1   ]

C2 = [  3  6  4  2 ;  # coefficients's vector of the objective 2
        1  3  8  3 ;
        5  2  2  3 ;
        4  2  3  5   ]

n  = size(C2,1)       # number of lines/columns

# ---- setting the model
bilap = vModel( GLPK.Optimizer )

@variable( bilap, x[1:n,1:n], Bin )

@addobjective( bilap , Min, sum( C1[i,j]*x[i,j] for i=1:n,j=1:n ) )
@addobjective( bilap , Min, sum( C2[i,j]*x[i,j] for i=1:n,j=1:n ) )

@constraint( bilap , cols[i=1:n], sum(x[i,j] for j=1:n) == 1 )
@constraint( bilap , rows[j=1:n], sum(x[i,j] for i=1:n) == 1 )

# ---- Invoking the solver (Chalmet method)
vSolve( bilap, method=:chalmet, step=0.5 )

# ---- Querying the results
Y_N = getY_N( bilap )

# ---- Displaying the results
printX_E( bilap )

max_stars_repo_name: vOptSolver/vOptGeneric

max_issues_repo_name: vOptSolver/vOptGeneric

max_forks_repo_name: vOptSolver/vOptGeneric


text:
# NOTE: this is copied from Pluto.run with launch_browser url changed
using Pluto: Sockets, HTTP, ServerSession, http_router_for, WorkspaceManager, MSG_DELIM, unpack, process_ws_message

function run_with_notebook_open(session::ServerSession, nb::Notebook)
    pluto_router = http_router_for(session)
    host = session.options.server.host
    port = session.options.server.port

    hostIP = parse(Sockets.IPAddr, host)
    if port === nothing
        port, serversocket = Sockets.listenany(hostIP, UInt16(1234))
    else
        try
            serversocket = Sockets.listen(hostIP, UInt16(port))
        catch e
            @error "Port with number $port is already in use. Use Pluto.run() to automatically select an available port."
            return
        end
    end

    kill_server = Ref{Function}(identity)

    servertask = @async HTTP.serve(hostIP, UInt16(port), stream=true, server=serversocket) do http::HTTP.Stream
        # messy messy code so that we can use the websocket on the same port as the HTTP server

        if HTTP.WebSockets.is_upgrade(http.message)
            try
                requestURI = http.message.target |> HTTP.URIs.unescapeuri |> HTTP.URI
                @assert endswith(requestURI.path, string(session.secret))

                HTTP.WebSockets.upgrade(http) do clientstream
                    if !isopen(clientstream)
                        return
                    end
                    try
                    while !eof(clientstream)
                        # This stream contains data received over the WebSocket.
                        # It is formatted and MsgPack-encoded by send(...) in PlutoConnection.js
                        try
                            parentbody = let
                                # For some reason, long (>256*512 bytes) WS messages get split up - `readavailable` only gives the first 256*512
                                data = UInt8[]
                                while !endswith(data, MSG_DELIM)
                                    if eof(clientstream)
                                        if isempty(data)
                                            return
                                        end
                                        @warn "Unexpected eof after" data
                                        append!(data, MSG_DELIM)
                                        break
                                    end
                                    append!(data, readavailable(clientstream))
                                end
                                # TODO: view to avoid memory allocation
                                unpack(data[1:end - length(MSG_DELIM)])
                            end
                            process_ws_message(session, parentbody, clientstream)
                        catch ex
                            if ex isa InterruptException
                                kill_server[]()
                            elseif ex isa HTTP.WebSockets.WebSocketError || ex isa EOFError
                                # that's fine!
                            elseif ex isa InexactError
                                # that's fine! this is a (fixed) HTTP.jl bug: https://github.com/JuliaWeb/HTTP.jl/issues/471
                                # TODO: remove this switch
                            else
                                bt = stacktrace(catch_backtrace())
                                @warn "Reading WebSocket client stream failed for unknown reason:" exception = (ex, bt)
                            end
                        end
                    end
                    catch ex
                        if ex isa InterruptException
                            kill_server[]()
                        else
                            bt = stacktrace(catch_backtrace())
                            @warn "Reading WebSocket client stream failed for unknown reason:" exception = (ex, bt)
                        end
                    end
                end
            catch ex
                if ex isa InterruptException
                    kill_server[]()
                elseif ex isa Base.IOError
                    # that's fine!
                elseif ex isa ArgumentError && occursin("stream is closed", ex.msg)
                    # that's fine!
                else
                    bt = stacktrace(catch_backtrace())
                    @warn "HTTP upgrade failed for unknown reason" exception = (ex, bt)
                end
            end
        else
            request::HTTP.Request = http.message
            request.body = read(http)
            HTTP.closeread(http)

            # If a "token" url parameter is passed in from binder, then we store it to add to every URL (so that you can share the URL to collaborate).
            params = HTTP.queryparams(HTTP.URI(request.target))
            if haskey(params, "token") && session.binder_token === nothing
                session.binder_token = params["token"]
            end

            request_body = IOBuffer(HTTP.payload(request))
            if eof(request_body)
                # no request body
                response_body = HTTP.handle(pluto_router, request)
            else
                @warn "HTTP request contains a body, huh?" request_body
            end

            request.response::HTTP.Response = response_body
            request.response.request = request
            try
                HTTP.startwrite(http)
                write(http, request.response.body)
                HTTP.closewrite(http)
            catch e
                if isa(e, Base.IOError) || isa(e, ArgumentError)
                    # @warn "Attempted to write to a closed stream at $(request.target)"
                else
                    rethrow(e)
                end
            end
        end
    end

    address = if session.options.server.root_url === nothing
        hostPretty = (hostStr = string(hostIP)) == "127.0.0.1" ? "localhost" : hostStr
        portPretty = Int(port)
        "http://$(hostPretty):$(portPretty)/"
    else
        session.options.server.root_url
    end
    Sys.set_process_title("Pluto server - $address")

    notebook_url = address * "edit?id=$(nb.notebook_id)"
    if session.options.server.launch_browser && open_in_default_browser(notebook_url)
        println("Opening $notebook_url in your default browser... ~ have fun!")
    else
        println("Go to $notebook_url in your browser to start your notebook ~ have fun!")
    end
    println()
    println("Press Ctrl+C in this terminal to stop Pluto")
    println()

    kill_server[] = () -> @sync begin
        println("\n\nClosing Pluto... Restart Julia for a fresh session. \n\nHave a nice day! 🎈")
        @async close(serversocket)
        # TODO: HTTP has a kill signal?
        # TODO: put do_work tokens back
        for client in values(session.connected_clients)
            @async close(client.stream)
        end
        empty!(session.connected_clients)
        for (notebook_id, ws) in WorkspaceManager.workspaces
            @async WorkspaceManager.unmake_workspace(ws)
        end
    end

    try
        # create blocking call and switch the scheduler back to the server task, so that interrupts land there
        wait(servertask)
    catch e
        if e isa InterruptException
            kill_server[]()
        elseif e isa TaskFailedException
            # nice!
        else
            rethrow(e)
        end
    end
end


text:
export @minimize

"""
`@minimize cost [st ctr] [with slv_opt]`

Minimize a given problem with cost function `cost`, constraints `ctr` and solver options `slv_opt`.

# Example

```julia
julia> A, b = randn(10,4), randn(10);

julia> @minimize ls(A*x-b) + 0.5*norm(x);
    it |      gamma |        fpr |        tau |        FBE |
 ------|------------|------------|------------|------------|
     1 | 2.9152e-02 | 2.7656e+00 | 1.0000e+00 | 5.5181e+00 |
     9 | 2.9152e-02 | 9.9682e-05 | 1.0000e+00 | 4.4086e+00 |

julia> @minimize ls(A*x-b) st x >= 0.;
    it |      gamma |        fpr |        tau |        FBE |
 ------|------------|------------|------------|------------|
     1 | 5.8304e-02 | 1.0068e+00 | 1.0000e+00 | 6.6282e+00 |
     3 | 5.8304e-02 | 9.5210e-16 | 1.0000e+00 | 6.5654e+00 |

julia> it, slv = @minimize ls(A*x-b) st norm(x) == 2.0 with PG(maxit = 5);
    it |      gamma |        fpr |
 ------|------------|------------|
     1 | 6.1373e-02 | 2.2090e+01 |
     5 | 3.0686e-02 | 5.5190e-01 |

```

Returns as output a tuple containing the number of iterations and the constructed solver.

"""
macro minimize(cf::Union{Expr, Symbol})
	cost = esc(cf)
    return :(solve(problem($(cost)), default_solver()))
end

macro minimize(cf::Union{Expr, Symbol}, s::Symbol, cstr::Union{Expr, Symbol})
	cost = esc(cf)
	if s == :(st)
		constraints = esc(cstr)
		return :(solve(problem($(cost), $(constraints)), default_solver()))
	elseif s == :(with)
		solver = esc(cstr)
        return :(solve(problem($(cost)), $(solver)))
	else
		error("wrong symbol after cost function! use `st` or `with`")
	end
end

macro minimize(cf::Union{Expr, Symbol}, s::Symbol, cstr::Union{Expr, Symbol}, w::Symbol, slv::Union{Expr, Symbol})
	cost = esc(cf)
	s != :(st) && error("wrong symbol after cost function! use `st`")
	constraints = esc(cstr)
	w != :(with) && error("wrong symbol after constraints! use `with`")
	solver = esc(slv)
    return :(solve(problem($(cost), $(constraints)), $(solver)))
end