text:
# RLWE LUTs (look-up tables)
# Copyright (C): Jaka Smrekar (vinctux) <vinctux@outlook.com>, 2016
#
# Reference:       https://www.microsoft.com/en-us/research/wp-content/uploads/2016/05/RLWE-1.pdf
# x86 C Reference: https://github.com/dconnolly/msr-latticecrypto/blob/master/ntt_constants.c

# Precomputed table Ω, storing primitive
# roots. Used in plain NTT.
global const Ω = [12288, 1479, 8246, 4134, 6429, 1351, 7678, 7935, 5559]

# Precomputed table Ψ_{rev, k^{-1}}, storing
# scaled powers of ψ in bit-reversed order.
# Support:  N = {512, 1024}, at Q = 12289
const Ψ = Dict{Integer, Array{Int16, 1}}(
        512  => [
             8193,   493,  6845,  9908,  1378, 10377,  7952,   435, 10146,  1065,   404,  7644,
             1207,  3248, 11121,  5277,  2437,  3646,  2987,  6022,  9867,  6250, 10102,  9723,
             1002,  7278,  4284,  7201,   875,  3780,  1607,  4976,  8146,  4714,   242,  1537,
             3704,  9611,  5019,   545,  5084, 10657,  4885, 11272,  3066, 12262,  3763, 10849,
             2912,  5698, 11935,  4861,  7277,  9808, 11244,  2859,  7188,  1067,  2401, 11847,
              390, 11516,  8511,  3833,  2780,  7094,  4895,  1484,  2305,  5042,  8236,  2645,
             7875,  9442,  2174,  7917,  1689,  3364,  4057,  3271, 10863,  4654,  1777, 10626,
             3636,  7351,  9585,  6998,   160,  3149,  4437, 12286, 10123,  3915,  7370, 12176,
             4048,  2249,  2884,  1153,  9103,  6882,  2126, 10659,  3510,  5332,  2865,  9919,
             9320,  8311,  9603,  9042,  3016, 12046,  9289, 11618,  7098,  3136,  9890,  3400,
             2178,  1544,  5559,   420,  8304,  4905,   476,  3531,  9326,  4896,  9923,  3051,
             3091,    81,  1000,  4320,  1177,  8034,  9521, 10654, 11563,  7678, 10436, 12149,
             3014,  9088,  5086,  1326, 11119,  2319, 11334,   790,  2747,  7443,  3135,  3712,
             1062,  9995,  7484,  8736,  9283,  2744, 11726,  2975,  9664,   949,  7468,  9650,
             7266,  5828,  6561,  7698,  3328,  6512,  1351,  7311,  8155,  5736,   722, 10984,
             4043,  7143, 10810,     1,  8668,  2545,  3504,  8747, 11077,  1646,  9094,  5860,
             1759,  8582,  3694,  7110,  8907, 11934,  8058,  9741,  9558,  3932,  5911,  4890,
             3637,  8830,  5542, 12144,  5755,  7657,  7901, 11029, 11955,  9863, 10861,  1696,
             3284,  2881,  7197,  2089,  9000,  2013,   729,  9048, 11809,  2842, 11267,     9,
             6498,   544,  2468,   339,  1381,  2525,  8112,  3584,  6958,  4989, 10616,  8011,
             5374,  9452, 12159,  4354,  9893,  7837,  3296,  8340,  7222,  2197,   118,  2476,
             5767,   827,  8541, 11336,  8855,  8760,  9381,   218,  9928, 10446,  9259,  4115,
             6142,  2447,  3963, 11713,  1954,  2051,  1805,  2882,   453,  6381, 11871,  8517,
             4774,  6860,  4737,  1293,   156,  9522,  8320,  3991,  5876,  2281, 10258,  6956,
             1489,  2500,  1583,  6347, 11026, 12240,  6374,  1483,   350,  1512, 10474,  6906,
             9087,  7796,  5369,  2057, 10314,  3757,  9364, 11942,  7535, 10431,   426,  3315,
             2738,  6421,  2655,  6554,   723,   174,  1693,  9280,  6099,   295,  5766, 11637,
             8527,  2919,  8273,  8212,  3728,  8240,  6299,  1159,  1146, 11341, 11964, 10885,
             5297,  6170,  3956,  1360, 11089,  7105,  9734,  6167, 10695,  1962,  5106,  6328,
             9597,   168,  7991,  8960,  6370,  7856,  3834,  5257, 10542,  9166,  9235,  5486,
             6507,  1566,  2948,  9786, 11606,  9830,  8633, 12225,  8049,  8719, 11454,  6224,
             8243,   709,  1319,  9139,  1958,  7967, 10211, 11177,  8210,  1058, 11848, 11367,
            11239,  7753,  5445,  3860,  9606,  1190,  8471,  6118,  3789,   147,  5456,  7840,
             7540,  5537,  4789,  4467,  4075,  5315,  4324,  4916, 10120, 11767,  7210,  9027,
             1973,  5574, 11011,  2344,  8775,  1041,  1018,  6364, 11821,  8301, 11907,   316,
             6950,  5446,  6093,  3710, 10256,  3998, 10367,  8410,  1254, 11316,  5435,  1359,
             7083,  5529,  9090, 12233,  8724, 11635, 10587,  1987,  6427,  6136,  6874,  3643,
              400,  1728,  4948,  6137,  5057,  7591,  3445,  7509,  2049,  7377, 10968,   192,
             5241,  9369,  9162,  8120,   787,  8807,  1010,  6821,  6415,   677,  6234,  3336,
            12237,  9115,  1323,  2766, 12138, 10162,  8332,  9450,  2505,  5906, 10710, 11858,
             4782,  6403,  9260,  5594,  8076, 11785,   605,  9987,  3600,  3263,  7665,  6077,
              421,  8209,  6068,  3602, 11286,  3532, 12048, 12231,  7280,  1956, 11404,  6008,
              8851,  2844,   975,  4212,  5681, 8812, 12147, 11184
        ],
        1024 => [
             8193,   493,  6845,  9908,  1378, 10377,  7952,   435, 10146,  1065,   404,  7644,
             1207,  3248, 11121,  5277,  2437,  3646,  2987,  6022,  9867,  6250, 10102,  9723,
             1002,  7278,  4284,  7201,   875,  3780,  1607,  4976,  8146,  4714,   242,  1537,
             3704,  9611,  5019,   545,  5084, 10657,  4885, 11272,  3066, 12262,  3763, 10849,
             2912,  5698, 11935,  4861,  7277,  9808, 11244,  2859,  7188,  1067,  2401, 11847,
              390, 11516,  8511,  3833,  2780,  7094,  4895,  1484,  2305,  5042,  8236,  2645,
             7875,  9442,  2174,  7917,  1689,  3364,  4057,  3271, 10863,  4654,  1777, 10626,
             3636,  7351,  9585,  6998,   160,  3149,  4437, 12286, 10123,  3915,  7370, 12176,
             4048,  2249,  2884,  1153,  9103,  6882,  2126, 10659,  3510,  5332,  2865,  9919,
             9320,  8311,  9603,  9042,  3016, 12046,  9289, 11618,  7098,  3136,  9890,  3400,
             2178,  1544,  5559,   420,  8304,  4905,   476,  3531,  9326,  4896,  9923,  3051,
             3091,    81,  1000,  4320,  1177,  8034,  9521, 10654, 11563,  7678, 10436, 12149,
             3014,  9088,  5086,  1326, 11119,  2319, 11334,   790,  2747,  7443,  3135,  3712,
             1062,  9995,  7484,  8736,  9283,  2744, 11726,  2975,  9664,   949,  7468,  9650,
             7266,  5828,  6561,  7698,  3328,  6512,  1351,  7311,  8155,  5736,   722, 10984,
             4043,  7143, 10810,     1,  8668,  2545,  3504,  8747, 11077,  1646,  9094,  5860,
             1759,  8582,  3694,  7110,  8907, 11934,  8058,  9741,  9558,  3932,  5911,  4890,
             3637,  8830,  5542, 12144,  5755,  7657,  7901, 11029, 11955,  9863, 10861,  1696,
             3284,  2881,  7197,  2089,  9000,  2013,   729,  9048, 11809,  2842, 11267,     9,
             6498,   544,  2468,   339,  1381,  2525,  8112,  3584,  6958,  4989, 10616,  8011,
             5374,  9452, 12159,  4354,  9893,  7837,  3296,  8340,  7222,  2197,   118,  2476,
             5767,   827,  8541, 11336,  3434,  3529,  2908, 12071,  2361,  1843,  3030,  8174,
             6147,  9842,  8326,   576, 10335, 10238, 10484,  9407, 11836,  5908,   418,  3772,
             7515,  5429,  7552, 10996, 12133,  2767,  3969,  8298,  6413, 10008,  2031,  5333,
            10800,  9789, 10706,  5942,  1263,    49,  5915, 10806, 11939, 10777,  1815,  5383,
             3202,  4493,  6920, 10232,  1975,  8532,  2925,   347,  4754,  1858, 11863,  8974,
             9551,  5868,  9634,  5735, 11566, 12115, 10596,  3009,  6190, 11994,  6523,   652,
             3762,  9370,  4016,  4077,  8561,  4049,  5990, 11130, 11143,   948,   325,  1404,
             6992,  6119,  8333, 10929,  1200,  5184,  2555,  6122,  1594, 10327,  7183,  5961,
             2692, 12121,  4298,  3329,  5919,  4433,  8455,  7032,  1747,  3123,  3054,  6803,
             5782, 10723,  9341,  2503,   683,  2459,  3656,    64,  4240,  3570,  835,   6065,
             4046, 11580, 10970,  3150, 10331,  4322,  2078,  1112,  4079, 11231,   441,   922,
             1050,  4536,  6844,  8429,  2683, 11099,  3818,  6171,  8500, 12142,  6833,  4449,
             4749,  6752,  7500,  7822,  8214,  6974,  7965,  7373,  2169,   522,  5079,  3262,
            10316,  6715,  1278,  9945,  3514, 11248, 11271,  5925,   468,  3988,   382, 11973,
             5339,  6843,  6196,  8579,  2033,  8291,  1922,  3879, 11035,   973,  6854, 10930,
             5206,  6760,  3199,    56,  3565,   654,  1702, 10302,  5862,  6153,  5415,  8646,
            11889, 10561,  7341,  6152,  7232,  4698,  8844,  4780, 10240,  4912,  1321, 12097,
             7048,  2920,  3127,  4169, 11502,  3482, 11279,  5468,  5874, 11612,  6055,  8953,
               52,  3174, 10966,  9523,   151,  2127,  3957,  2839,  9784,  6383,  1579,   431,
             7507,  5886,  3029,  6695,  4213,   504, 11684,  2302,  8689,  9026,  4624,  6212,
            11868,  4080,  6221,  8687,  1003,  8757,   241,    58,  5009, 10333,   885,  6281,
             3438,  9445, 11314,  8077,  6608,  3477,   142,  1105,  8841,   343,  4538,  1908,
             1208,  4727,  7078, 10423, 10125,  6873, 11573, 10179,   416,   814,  1705,  2450,
             8700,   717,  9307,  1373,  8186,  2429, 10568, 10753,  7228, 11071,   438,  8774,
             5993,  3278,  4209,  6877,  3449,  1136,  3708,  3238,  2926,  1826,  4489,  3171,
             8024,  8611,  1928,   464,  3205,  8930,  7080,  1092, 10900, 10221, 11943,  4404,
             9126,  4032,  7449,  6127,  8067, 10763,   125,   540,  8921,  8062,   612,  8051,
            12229,  9572,  9089, 10754, 10029,    68,  6453,  7723,  4781,  4924,  1014,   448,
             3942,  5232,  1327,  8682,  3744,  7326,  3056,  9761,  5845,  5588,   412,  7187,
             3975,  4883,  3087,  6454,  2257,  7784,  5676,  1417,  8400, 11710,  5596,  5987,
             9175,  2769,  5966,   212,  6555, 11113,  5508, 11014,  1125,  4860, 10844,  1131,
             4267,  6636,  2275,  9828,  5063,  4176,  3765,  1518,  8794,  4564, 10224,  5826,
             3534,  3961,  4145, 10533,   506, 11034,  6505, 10897,  2674, 10077,  3338,  9013,
             3511,  6811, 11111,  2776,  1165,  2575,  8881, 10347,   377,  4578, 11914, 10669,
            10104,   392, 10453,   425,  9489,   193,  2231,  6197,  1038, 11366,  6204,  8122,
             2894,  3654, 10975, 10545,  6599,  2455, 11951,  3947,    20,  5002,  5163,  4608,
             8946,  8170, 10138,  1522,  8665, 10397,  3344,  5598, 10964,  6565, 11260,  1945,
            11041,  9847,  7174,  4939,  2148,  6330,  3959,  5797,  4913,  3528,  8054,  3825,
             8914,  9998,  4335,  8896,  9342,  3982,  6680, 11653,  7790,  6617,  1737,   622,
            10485, 10886,  6195,  7100,  1687,   406, 12143,  5268,  9389, 12050,   994,  7735,
             5464,  7383,  4670,   512,   364,  9929,  3028,  5216,  5518,  1226,  7550,  8038,
             7043,  7814, 11053,  3017,  3121,  7584,  2600, 11232,  6780, 12085,  5219,  1409,
             9600,  4605,  8151, 12109,   463,  8882,  8308, 10821,  9247, 10945,  9806,  2054,
             6203,  6643,  3120,  6105,  8348,  8536,  6919,  8753, 11007,  8717,  9457,  2021,
             9060,  4730,  3929, 10583,  3723,   845,  1936,     7,  5054,  3154,  3285,  4360,
             3805, 11522,  2213,  4153, 12239, 12073,  5526,   769,  4099,  3944,  5604,  5530,
            11024,  9282,  2171,  3480,  7434,  8520,  3232, 11996,  9656,  1406,  2945,  5349,
             7207,  4590, 11607, 11309,  5202,   844,  7082,  4050,  8016,  9068,  9694,  8452,
             7000,  5662,   567,  2941,  8619,  3808,  4987,  2373,  5135,    63,  7605,  3360,
            11839, 10345,   578,  6921,  7628,   510,  5386,  2622,  7806,  5703, 10783,  9224,
            11379,  5900,  4719, 11538,  3502,  5789, 10631,  5618,   826,  5043,  3090, 10891,
             9951,  7596,  2293, 11872,  6151,  3469,  4443,  8871,  1555,  1802,  5103,  1891,
             1223,  2334,  7878,  1590,   881,   365,  1927, 11274,  4510,  9652,  2946,  6828,
             1280,   614, 10918, 12265,  7250,  6742,  9804, 11385,  2276, 11307,  2593,   879,
             7899,  8071,  3454,  8531,  3795,  9021,  5776,  1849,  7766,  7988,   457,     8,
              530,  9663,  7785, 11511,  3578,  7592, 10588,  3466,  8972,  9757,  3332,   139,
             2046,  2940, 10808,  9332,   874,  2301,  5650, 12119,   150,   648,  8000,  9982,
             9416,  2827,  2434, 11498,  6481, 12268,  9754, 11169, 11823, 11259,  3821, 10608,
             2929,  6263,  4649,  6320,  9687, 10388,   502,  5118,  8496,  6226, 10716,  8443,
             7624,  6883,  9269,  6616,  8620,  5287,   944,  7519,  6125,  1882, 11249, 10254,
             5410,  1251,  1790,  5275,  8449, 10447,  4113,    72,  2828,  4352,  7455,  2712,
            11048,  7911,  3451,  4094,  6508,  3045, 11194,  2643,  1783,  7211,  4974,  7724,
             9811,  9449,  3019,  4194,  2730,  6878, 10421,  2253,  4518,  9195,  7469, 11129,
             9173, 12100,  1763,  2209,  9617,  5170,   865,  1279,  1694, 10759,  8420,  4423,
            10555,  3815,  5832, 10939
            ]
)

# Precomputed table Ψ^{-1}_{rev, k^{-1}}, storing
# scaled powers of ψ^{-1} in bit-reversed order.
# Support:  N = {512, 1024}, at Q = 12289
const Ψ_ = Dict{Integer, Array{Int16, 1}}(
        512  => [
             8193, 11796,  2381,  5444, 11854,  4337,  1912, 10911,  7012,  1168,  9041, 11082,
             4645, 11885, 11224,  2143,  7313, 10682,  8509, 11414,  5088,  8005,  5011, 11287,
             2566,  2187,  6039,  2422,  6267,  9302,  8643,  9852,  8456,  3778,   773, 11899,
              442,  9888, 11222,  5101,  9430,  1045,  2481,  5012,  7428,   354,  6591,  9377,
             1440,  8526,    27,  9223,  1017,  7404,  1632,  7205, 11744,  7270,  2678,  8585,
            10752, 12047,  7575,  4143,  8758, 11813,  7384,  3985, 11869,  6730, 10745, 10111,
             8889,  2399,  9153,  5191,   671,  3000,   243,  9273,  3247,  2686,  3978,  2969,
             2370,  9424,  6957,  8779,  1630, 10163,  5407,  3186, 11136,  9405, 10040,  8241,
              113,  4919,  8374,  2166,     3,  7852,  9140, 12129,  5291,  2704,  4938,  8653,
             1663, 10512,  7635,  1426,  9018,  8232,  8925, 10600,  4372, 10115,  2847,  4414,
             9644,  4053,  7247,  9984, 10805,  7394,  5195,  9509,   953,  3748, 11462,  6522,
             9813, 12171, 10092,  5067,  3949,  8993,  4452,  2396,  7935,   130,  2837,  6915,
             4278,  1673,  7300,  5331,  8705,  4177,  9764, 10908, 11950,  9821, 11745,  5791,
            12280,  1022,  9447,   480,  3241, 11560, 10276,  3289, 10200,  5092,  9408,  9005,
            10593,  1428,  2426,   334,  1260,  4388,  4632,  6534,   145,  6747,  3459,  8652,
             7399,  6378,  8357,  2731,  2548,  4231,   355,  3382,  5179,  8595,  3707, 10530,
             6429,  3195, 10643,  1212,  3542,  8785,  9744,  3621, 12288,  1479,  5146,  8246,
             1305, 11567,  6553,  4134,  4978, 10938,  5777,  8961,  4591,  5728,  6461,  5023,
             2639,  4821, 11340,  2625,  9314,   563,  9545,  3006,  3553,  4805,  2294, 11227,
             8577,  9154,  4846,  9542, 11499,   955,  9970,  1170, 10963,  7203,  3201,  9275,
              140,  1853,  4611,   726,  1635,  2768,  4255, 11112,  7969, 11289, 12208,  9198,
             9238,  2366,  7393,  2963,  1105,   142,  3477,  6608,  8077, 11314,  9445,  3438,
             6281,   885, 10333,  5009,    58,   241,  8757,  1003,  8687,  6221,  4080, 11868,
             6212,  4624,  9026,  8689,  2302, 11684,   504,  4213,  6695,  3029,  5886,  7507,
              431,  1579,  6383,  9784,  2839,  3957,  2127,   151,  9523, 10966,  3174,    52,
             8953,  6055, 11612,  5874,  5468, 11279,  3482, 11502,  4169,  3127,  2920,  7048,
            12097,  1321,  4912, 10240,  4780,  8844,  4698,  7232,  6152,  7341, 10561, 11889,
             8646,  5415,  6153,  5862, 10302,  1702,   654,  3565,    56,  3199,  6760,  5206,
            10930,  6854,   973, 11035,  3879,  1922,  8291,  2033,  8579,  6196,  6843,  5339,
            11973,   382,  3988,   468,  5925, 11271, 11248,  3514,  9945,  1278,  6715, 10316,
             3262,  5079,   522,  2169,  7373,  7965,  6974,  8214,  7822,  7500,  6752,  4749,
             4449,  6833, 12142,  8500,  6171,  3818, 11099,  2683,  8429,  6844,  4536,  1050,
              922,   441, 11231,  4079,  1112,  2078,  4322, 10331,  3150, 10970, 11580,  4046,
             6065,   835,  3570,  4240,    64,  3656,  2459,   683,  2503,  9341, 10723,  5782,
             6803,  3054,  3123,  1747,  7032,  8455,  4433,  5919,  3329,  4298, 12121,  2692,
             5961,  7183, 10327,  1594,  6122,  2555,  5184,  1200, 10929,  8333,  6119,  6992,
             1404,   325,   948, 11143, 11130,  5990,  4049,  8561,  4077,  4016,  9370,  3762,
              652,  6523, 11994,  6190,  3009, 10596, 12115, 11566,  5735,  9634,  5868,  9551,
             8974, 11863,  1858,  4754,   347,  2925,  8532,  1975, 10232,  6920,  4493,  3202,
             5383,  1815, 10777, 11939, 10806,  5915,    49,  1263,  5942, 10706,  9789, 10800,
             5333,  2031, 10008,  6413,  8298,  3969,  2767, 12133, 10996,  7552,  5429,  7515,
             3772,   418,  5908, 11836,  9407, 10484, 10238, 10335,   576,  8326,  9842,  6147,
             8174,  3030,  1843,  2361, 12071,  2908,  3529,  3434
            ],

        1024 => [
             8193, 11796,  2381,  5444, 11854,  4337,  1912, 10911,  7012,  1168,  9041, 11082,
             4645, 11885, 11224,  2143,  7313, 10682,  8509, 11414,  5088,  8005,  5011, 11287,
             2566,  2187,  6039,  2422,  6267,  9302,  8643,  9852,  8456,  3778,   773, 11899,
              442,  9888, 11222,  5101,  9430,  1045,  2481,  5012,  7428,   354,  6591,  9377,
             1440,  8526,    27,  9223,  1017,  7404,  1632,  7205, 11744,  7270,  2678,  8585,
            10752, 12047,  7575,  4143,  8758, 11813,  7384,  3985, 11869,  6730, 10745, 10111,
             8889,  2399,  9153,  5191,   671,  3000,   243,  9273,  3247,  2686,  3978,  2969,
             2370,  9424,  6957,  8779,  1630, 10163,  5407,  3186, 11136,  9405, 10040,  8241,
              113,  4919,  8374,  2166,     3,  7852,  9140, 12129,  5291,  2704,  4938,  8653,
             1663, 10512,  7635,  1426,  9018,  8232,  8925, 10600,  4372, 10115,  2847,  4414,
             9644,  4053,  7247,  9984, 10805,  7394,  5195,  9509,   953,  3748, 11462,  6522,
             9813, 12171, 10092,  5067,  3949,  8993,  4452,  2396,  7935,   130,  2837,  6915,
             4278,  1673,  7300,  5331,  8705,  4177,  9764, 10908, 11950,  9821, 11745,  5791,
            12280,  1022,  9447,   480,  3241, 11560, 10276,  3289, 10200,  5092,  9408,  9005,
            10593,  1428,  2426,   334,  1260,  4388,  4632,  6534,   145,  6747,  3459,  8652,
             7399,  6378,  8357,  2731,  2548,  4231,   355,  3382,  5179,  8595,  3707, 10530,
             6429,  3195, 10643,  1212,  3542,  8785,  9744,  3621, 12288,  1479,  5146,  8246,
             1305, 11567,  6553,  4134,  4978, 10938,  5777,  8961,  4591,  5728,  6461,  5023,
             2639,  4821, 11340,  2625,  9314,   563,  9545,  3006,  3553,  4805,  2294, 11227,
             8577,  9154,  4846,  9542, 11499,   955,  9970,  1170, 10963,  7203,  3201,  9275,
              140,  1853,  4611,   726,  1635,  2768,  4255, 11112,  7969, 11289, 12208,  9198,
             9238,  2366,  7393,  2963, 11184, 12147,  8812,  5681,  4212,   975,  2844,  8851,
             6008, 11404,  1956,  7280, 12231, 12048,  3532, 11286,  3602,  6068,  8209,   421,
             6077,  7665,  3263,  3600,  9987,   605, 11785,  8076,  5594,  9260,  6403,  4782,
            11858, 10710,  5906,  2505,  9450,  8332, 10162, 12138,  2766,  1323,  9115, 12237,
             3336,  6234,   677,  6415,  6821,  1010,  8807,   787,  8120,  9162,  9369,  5241,
              192, 10968,  7377,  2049,  7509,  3445,  7591,  5057,  6137,  4948,  1728,   400,
             3643,  6874,  6136,  6427,  1987, 10587, 11635,  8724, 12233,  9090,  5529,  7083,
             1359,  5435, 11316,  1254,  8410, 10367,  3998, 10256,  3710,  6093,  5446,  6950,
              316, 11907,  8301, 11821,  6364,  1018,  1041,  8775,  2344, 11011,  5574,  1973,
             9027,  7210, 11767, 10120,  4916,  4324,  5315,  4075,  4467,  4789,  5537,  7540,
             7840,  5456,   147,  3789,  6118,  8471,  1190,  9606,  3860,  5445,  7753, 11239,
            11367, 11848,  1058,  8210, 11177, 10211,  7967,  1958,  9139,  1319,   709,  8243,
             6224, 11454,  8719,  8049, 12225,  8633,  9830, 11606,  9786,  2948,  1566,  6507,
             5486,  9235,  9166, 10542,  5257,  3834,  7856,  6370,  8960,  7991,   168,  9597,
             6328,  5106,  1962, 10695,  6167,  9734,  7105, 11089,  1360,  3956,  6170,  5297,
            10885, 11964, 11341,  1146,  1159,  6299,  8240,  3728,  8212,  8273,  2919,  8527,
            11637,  5766,   295,  6099,  9280,  1693,   174,   723,  6554,  2655,  6421,  2738,
             3315,   426, 10431,  7535, 11942,  9364,  3757, 10314,  2057,  5369,  7796,  9087,
             6906, 10474,  1512,   350,  1483,  6374, 12240, 11026,  6347,  1583,  2500,  1489,
             6956, 10258,  2281,  5876,  3991,  8320,  9522,   156,  1293,  4737,  6860,  4774,
             8517, 11871,  6381,   453,  2882,  1805,  2051,  1954, 11713,  3963,  2447,  6142,
             4115,  9259, 10446,  9928,   218,  9381,  8760,  8855,  1350,  6457,  8474,  1734,
             7866,  3869,  1530, 10595, 11010, 11424,  7119,  2672, 10080, 10526,   189,  3116,
             1160,  4820,  3094,  7771, 10036,  1868,  5411,  9559,  8095,  9270,  2840,  2478,
             4565,  7315,  5078, 10506,  9646,  1095,  9244,  5781,  8195,  8838,  4378,  1241,
             9577,  4834,  7937,  9461, 12217,  8176,  1842,  3840,  7014, 10499, 11038,  6879,
             2035,  1040, 10407,  6164,  4770, 11345,  7002,  3669,  5673,  3020,  5406,  4665,
             3846,  1573,  6063,  3793,  7171, 11787,  1901,  2602,  5969,  7640,  6026,  9360,
             1681,  8468,  1030,   466,  1120,  2535,    21,  5808,   791,  9855,  9462,  2873,
             2307,  4289, 11641, 12139,   170,  6639,  9988, 11415,  2957,  1481,  9349, 10243,
            12150,  8957,  2532,  3317,  8823,  1701,  4697,  8711,   778,  4504,  2626, 11759,
            12281, 11832,  4301,  4523, 10440,  6513,  3268,  8494,  3758,  8835,  4218,  4390,
            11410,  9696,   982, 10013,   904,  2485,  5547,  5039,    24,  1371, 11675, 11009,
             5461,  9343,  2637,  7779,  1015, 10362, 11924, 11408, 10699,  4411,  9955, 11066,
            10398,  7186, 10487, 10734,  3418,  7846,  8820,  6138,   417,  9996,  4693,  2338,
             1398,  9199,  7246, 11463,  6671,  1658,  6500,  8787,   751,  7570,  6389,   910,
             3065,  1506,  6586,  4483,  9667,  6903, 11779,  4661,  5368, 11711,  1944,   450,
             8929,  4684, 12226,  7154,  9916,  7302,  8481,  3670,  9348, 11722,  6627,  5289,
             3837,  2595,  3221,  4273,  8239,  5207, 11445,  7087,   980,   682,  7699,  5082,
             6940,  9344, 10883,  2633,   293,  9057,  3769,  4855,  8809, 10118,  3007,  1265,
             6759,  6685,  8345,  8190, 11520,  6763,   216,    50,  8136, 10076,   767,  8484,
             7929,  9004,  9135,  7235, 12282, 10353, 11444,  8566,  1706,  8360,  7559,  3229,
            10268,  2832,  3572,  1282,  3536,  5370,  3753,  3941,  6184,  9169,  5646,  6086,
            10235,  2483,  1344,  3042,  1468,  3981,  3407, 11826,   180,  4138,  7684,  2689,
            10880,  7070,   204,  5509,  1057,  9689,  4705,  9168,  9272,  1236,  4475,  5246,
             4251,  4739, 11063,  6771,  7073,  9261,  2360, 11925, 11777,  7619,  4906,  6825,
             4554, 11295,   239,  2900,  7021,   146, 11883, 10602,  5189,  6094,  1403,  1804,
            11667, 10552,  5672,  4499,   636,  5609,  8307,  2947,  3393,  7954,  2291,  3375,
             8464,  4235,  8761,  7376,  6492,  8330,  5959, 10141,  7350,  5115,  2442,  1248,
            10344,  1029,  5724,  1325,  6691,  8945,  1892,  3624, 10767,  2151,  4119,  3343,
             7681,  7126,  7287, 12269,  8342,   338,  9834,  5690,  1744,  1314,  8635,  9395,
             4167,  6085,   923, 11251,  6092, 10058, 12096,  2800, 11864,  1836, 11897,  2185,
             1620,   375,  7711, 11912,  1942,  3408,  9714, 11124,  9513,  1178,  5478,  8778,
             3276,  8951,  2212,  9615,  1392,  5784,  1255, 11783,  1756,  8144,  8328,  8755,
             6463,  2065,  7725,  3495, 10771,  8524,  8113,  7226,  2461, 10014,  5653,  8022,
            11158,  1445,  7429, 11164,  1275,  6781,  1176,  5734, 12077,  6323,  9520,  3114,
             6302,  6693,   579,  3889, 10872,  6613,  4505, 10032,  5835,  9202,  7406,  8314,
             5102, 11877,  6701,  6444,  2528,  9233,  4963,  8545,  3607, 10962,  7057,  8347,
            11841, 11275,  7365,  7508,  4566,  5836, 12221,  2260,  1535,  3200,  2717,    60,
             4238, 11677,  4227,  3368, 11749, 12164,  1526,  4222,  6162,  4840,  8257,  3163,
             7885,   346,  2068,  1389, 11197,  5209,  3359,  9084, 11825, 10361,  3678,  4265,
             9118,  7800, 10463,  9363,  9051,  8581, 11153,  8840,  5412,  8080,  9011,  6296,
             3515, 11851,  1218,  5061,  1536,  1721,  9860,  4103, 10916,  2982, 11572,  3589,
             9839, 10584, 11475, 11873,  2110,   716,  5416,  2164,  1866,  5211,  7562, 11081,
            10381,  7751, 11946,  3448
            ]
)


text:
# Problem with ML ML estimator
# The initial point is valid;
# L-BFGS fails during line search for an SDE(2) model

using EstimateSDE
using EstimateSDE:AddSdeRoot,ResampleAr,InitialModel,estimate_trace,
  has_statstate, ml_estimate_init
using Random
using Distributions
using Printf
using DataStructures
using TypedTables
using StatsPlots
using JuliaDB
using .Threads

## Settings
proc0 = testprocess(:C)
proc = SdeProcess(proc0.eq) # Use stationary initial conditions
N = 200
Tav = 1.0
Tkld = 0.5
p = 2
Trs = [Tkld;Tav]
nspec = 1000

rng = MersenneTwister(539872437) # fails!
t,y = rand(rng,proc,N,Tav)
problem = SdeProblemLike(:stat,:elimQ,t,y,p,true,nothing)

## Maximum Likelihood estimate
tor = MlEstimator()
init_tors = []
push!(init_tors,AddSdeRoot(-10/Tav))
[push!(init_tors,ResampleAr(Tr,m,-10/Tav)) for Tr in Trs, m in [:nearest,:linear]]
est_ml,trace = estimate_trace(tor,problem,init_tors)

## Failed initiation
trace_failed = trace[3]
proc_ini = trace_failed[:proc_ini]
est,opt = ml_estimate_init(:dynamichmc,tor,problem,proc_ini)

ext: jl


text:
##==================================================
##
## Script name: Causal Matrix Completion
##
##
## Purpose of script: Given a Matrix of observed Y(0), indicator for missing
## entries and parameters for the algorithm, run the causal matrix Completion
## method as described in Athey et al. (JASA 2021)
## Author: Eray Turkel
##
## Date Last Edited: 09-22-2021
##
## Email: eturkel@stanford.edu
##
## Input:  TBC

using LinearAlgebra: pinv, svd, norm, I, dot, diagm,Diagonal, QRIteration
using Arpack: svds
using Statistics: mean
using Random

function shrink_matrix(A, lambda)
    # Implements the matrix shrinkage operator.
    # Denoting the SVD of A= U S V',
    # Replaces the matrix of eigenvalues S with S_shrink
    # where S_shrink replaces the i'th singular value
    # in the matrix S with max(sigma_i(A)- lambda, 0).
    N = size(A)[1];
    T = size(A)[2];
    rank_max=min(N,T);
    F = svd(A);
    S = F.S;
    num_singular=size(S)[1];
    for i in 1:num_singular
        S[i]=max(S[i]-lambda,0);
    end
    return  F.U* Diagonal(S) * F.Vt;
end

function FEestimator(Y, L, O)
    #Extract the time and unit fixed effects
    #Given Y, current estimate of L, and Observed Matrix O
    est_matrix=( (Y-L) .* O) ;
    N=size(Y)[1];
    T=size(Y)[2];
    time_columns=[ceil(i/N)==j for i in 1:(N*T), j in 1:T];
    unit_columns=[(i%N)==(j%N) for i in 1:(N*T), j in 1:N];
    X=hcat(time_columns,unit_columns);
    X=X[(vec(O).!=0),:];
    y=vec(est_matrix)[(vec(O).!=0),];
    fixed_effects=inv(X' * X) * (X' * vec(y));
    time_FE= fixed_effects[1:T]
    unit_FE= fixed_effects[(T+1):(T+N)]
    return (T=time_FE, U=unit_FE);
end

function softImpute_noFE(Y,O,lambda, n_iter,verbose= true)
    # Y is the Y(0) matrix of size N by T
    # O is the matrix of size N by T marking observed entries,
    # with O_(i,t)=1 if Y(0)_(i,t) is observed, O_(i,t)=0 otherwise.
    # lambda is the regularization parameter
    N=size(Y)[1];
    T=size(Y)[2];
    N_Obs=(sum(sum(O,dims=2),dims=1))[1];
    O_orth=((O .- 1) .* (-1));
    P_o_Y=(Y .* O);
    L_prev= P_o_Y;
    L_next= shrink_matrix(((Y) .* O)+ (L_prev .* O_orth), lambda*(N_Obs/2));
    for i in 1:n_iter
        L_prev=deepcopy(L_next)
        L_next= shrink_matrix(((Y) .* O)+ (L_prev .* O_orth), lambda*(N_Obs/2));
        distance= norm((L_next-L_prev));
        if (verbose)
            println("Current Dist:");
            println(distance);
        end
    end
    return L_next
end

function softImpute(Y,O,lambda, n_iter, verbose= true)
    # Y is the Y(0) matrix of size N by T
    # O is the matrix of size N by T marking observed entries,
    # with O_(i,t)=1 if Y(0)_(i,t) is observed, O_(i,t)=0 otherwise.
    # lambda is the regularization parameter
    N=size(Y)[1];
    T=size(Y)[2];
    N_Obs=(sum(sum(O,dims=2),dims=1))[1];
    O_orth=((O .- 1) .* (-1));
    P_o_Y=(Y .* O);
    L_prev= P_o_Y;
    fixed_effects= FEestimator(Y,zeros(N,T),O);
    Time= ones(N,1) * (fixed_effects.T)';
    Unit= (fixed_effects.U) * ones(1,T);
    L_next= shrink_matrix(((Y-Time-Unit) .* O)+ (L_prev .* O_orth), lambda*(N_Obs/2));
    for i in 1:n_iter
        L_prev= deepcopy(L_next);
        fixed_effects= FEestimator(Y,L_prev,O);
        Time= ones(N,1) * (fixed_effects.T)';
        Unit= (fixed_effects.U) * ones(1,T);
        L_new= shrink_matrix(((Y-Time-Unit) .* O)+ ((L_prev) .* O_orth), lambda*(N_Obs/2));
        L_next= deepcopy(L_new);
        distance= norm((L_next-L_prev));
        if (verbose)
            println("Current Dist:");
            println(distance);
        end
    end
    return (L=L_next, T=Time, U=Unit, Y_est= L_next+ Time+ Unit)
end

function MC_estimator(Y,O,reg_lambda, n_iter, FE_estimate= false, verbose=true)
    if (FE_estimate)
        return softImpute(Y,O,reg_lambda, n_iter, verbose);
    else
        return softImpute_noFE(Y,O,reg_lambda, n_iter, verbose);
    end
end

unit_fixed_rand= 8   .+ randn(1000)*0.5;
time_fixed_rand= 0.5 .+ randn(50)*0.1;
#Create random low rank matrix for the true dgp
L_rand = randn(1000,50)*0.5;
L_low = svd(L_rand);
Sigmas= L_low.S;
num_sigmas=size(Sigmas)[1];
Sigmas[(num_sigmas):num_sigmas]=zeros(num_sigmas-(num_sigmas)+1);
L_rand=  L_low.U* Diagonal(Sigmas) * L_low.V;
eps = randn(1000,50)*0.05;

Y= unit_fixed_rand * ones(1,50) + ones(1000,1) * time_fixed_rand' + L_rand +eps;

N=size(Y)[1];
T=size(Y)[2];

O_1=[ j<15 for i in 1:(150), j in 1:T];
O_2=[ j<25 for i in 1:(150), j in 1:T];
O_3=[ j<35 for i in 1:(150), j in 1:T];
O_4=[ j<(T+1) for i in 1:(550), j in 1:T];

O=vcat(O_1,O_2,O_3,O_4);

(MC_estimator(Y,O,0.001,1000))

(MC_estimator(Y,O,0.001,50,true))


text:
import Base: isempty

export Translation,
       an_element,
       constraints_list,
       linear_map

"""
    Translation{N, S<:LazySet{N}, VN<:AbstractVector{N}} <: AbstractAffineMap{N, S}

Type that represents a lazy translation.

The translation of set `X` along vector `v` is the map:

```math
x ↦ x + v,\\qquad x ∈ X
```

A translation is a special case of an affine map ``A x + b, x ∈ X`` where the
linear map ``A`` is the identity matrix and the translation vector ``b`` is
``v``.

### Fields

- `X` -- set
- `v` -- vector that defines the translation

### Notes

The translation preserves convexity: if `X` is convex, then any translation of
`X` is convex as well.

### Example

```jldoctest translation
julia> X = BallInf([2.0, 2.0, 2.0], 1.0);

julia> v = [1.0, 0.0, 0.0]; # translation along dimension 1

julia> tr = Translation(X, v);

julia> typeof(tr)
Translation{Float64, BallInf{Float64, Vector{Float64}}, Vector{Float64}}

julia> tr.X
BallInf{Float64, Vector{Float64}}([2.0, 2.0, 2.0], 1.0)

julia> tr.v
3-element Vector{Float64}:
 1.0
 0.0
 0.0
```
The sum operator `+` is overloaded to create translations:

```jldoctest translation
julia> X + v == Translation(X, v)
true
```
And so does the Minkowski sum operator, `⊕`:

```jldoctest translation
julia> X ⊕ v == Translation(X, v)
true
```

The translation of a translation is performed immediately:

```jldoctest translation
julia> tr = (X+v)+v
Translation{Float64, BallInf{Float64, Vector{Float64}}, Vector{Float64}}(BallInf{Float64, Vector{Float64}}([2.0, 2.0, 2.0], 1.0), [2.0, 0.0, 0.0])

julia> tr.v
3-element Vector{Float64}:
 2.0
 0.0
 0.0
```

The dimension of a translation is obtained with the `dim` function:

```jldoctest translation
julia> dim(tr)
3
```

For the support vector (resp. support function) along vector `d`, use `σ` and
`ρ` respectively:

```jldoctest translation
julia> σ([1.0, 0.0, 0.0], tr)
3-element Vector{Float64}:
 5.0
 2.0
 2.0

julia> ρ([1.0, 0.0, 0.0], tr)
5.0
```
See the docstring of each of these functions for details.

The `an_element` function is useful to obtain an element of a translation:

```jldoctest translation
julia> e = an_element(tr)
3-element Vector{Float64}:
 4.0
 2.0
 2.0
```

The lazy linear map of a translation is an affine map, since the following
simplification rule applies: ``M * (X ⊕ v) = (M * X) ⊕ (M * v)``:

```jldoctest translation
julia> using LinearAlgebra: I

julia> M = Matrix(2.0I, 3, 3);

julia> Q = M * tr;

julia> Q isa AffineMap && Q.M == M && Q.X == tr.X && Q.v == 2 * tr.v
true
```

Use the `isempty` method to query if the translation is empty; it falls back
to the `isempty` method of the wrapped set:

```jldoctest translation
julia> isempty(tr)
false
```

The list of constraints of the translation of a polyhedron (in general, a set
whose `constraints_list` is available) can be computed from a lazy translation:

```jldoctest translation
julia> constraints_list(tr)
6-element Vector{HalfSpace{Float64, LazySets.Arrays.SingleEntryVector{Float64}}}:
 HalfSpace{Float64, LazySets.Arrays.SingleEntryVector{Float64}}([1.0, 0.0, 0.0], 5.0)
 HalfSpace{Float64, LazySets.Arrays.SingleEntryVector{Float64}}([0.0, 1.0, 0.0], 3.0)
 HalfSpace{Float64, LazySets.Arrays.SingleEntryVector{Float64}}([0.0, 0.0, 1.0], 3.0)
 HalfSpace{Float64, LazySets.Arrays.SingleEntryVector{Float64}}([-1.0, 0.0, 0.0], -3.0)
 HalfSpace{Float64, LazySets.Arrays.SingleEntryVector{Float64}}([0.0, -1.0, 0.0], -1.0)
 HalfSpace{Float64, LazySets.Arrays.SingleEntryVector{Float64}}([0.0, 0.0, -1.0], -1.0)
```
"""
struct Translation{N, S<:LazySet{N}, VN<:AbstractVector{N}} <: AbstractAffineMap{N, S}
    X::S
    v::VN

    # default constructor with dimension check
    function Translation(X::S, v::VN) where {N, VN<:AbstractVector{N}, S<:LazySet{N}}
        @assert dim(X) == length(v) "cannot create a translation of a set of " *
            "dimension $(dim(X)) along a vector of length $(length(v))"

        return new{N, S, VN}(X, v)
    end
end

isoperationtype(::Type{<:Translation}) = true
isconvextype(::Type{Translation{N, S, VN}}) where {N, S, VN} = isconvextype(S)

# constructor from a Translation: perform the translation immediately
Translation(tr::Translation{N}, v::AbstractVector{N}) where {N} = Translation(tr.X, tr.v + v)

"""
    +(X::LazySet, v::AbstractVector)

Convenience constructor for a translation.

### Input

- `X` -- set
- `v` -- vector

### Output

The symbolic translation of ``X`` along vector ``v``.
"""
+(X::LazySet, v::AbstractVector) = Translation(X, v)

# translation from the left
+(v::AbstractVector, X::LazySet) = Translation(X, v)

"""
    ⊕(X::LazySet, v::AbstractVector)

Unicode alias constructor ⊕ (`oplus`) for the lazy translation operator.
"""
⊕(X::LazySet, v::AbstractVector) = Translation(X, v)

# translation from the left
⊕(v::AbstractVector, X::LazySet) = Translation(X, v)

# the translation of a lazy linear map is a (lazy) affine map
Translation(lm::LinearMap, v::AbstractVector) = AffineMap(lm.M, lm.X, v)

# the linear map of a translation is a (lazy) affine map:
# M * (X ⊕ v) = (M * X) ⊕ (M * v)
LinearMap(M::AbstractMatrix, tr::Translation) = AffineMap(M, tr.X, M * tr.v)

# --- AbstractAffineMap interface functions ---

function matrix(tr::Translation{N}) where {N}
    return Diagonal(fill(one(N), dim(tr)))
end

function vector(tr::Translation)
    return tr.v
end

function set(tr::Translation)
    return tr.X
end

# ============================
# LazySet interface functions
# ============================

"""
    σ(d::AbstractVector, tr::Translation)

Return the support vector of a translation.

### Input

- `d`  -- direction
- `tr` -- translation

### Output

The support vector in the given direction.
If the direction has norm zero, the result depends on the wrapped set.
"""
function σ(d::AbstractVector, tr::Translation)
    return tr.v + σ(d, tr.X)
end

"""
    ρ(d::AbstractVector, tr::Translation)

Return the support function of a translation.

### Input

- `d`  -- direction
- `tr` -- translation

### Output

The support function in the given direction.
"""
function ρ(d::AbstractVector, tr::Translation)
    return dot(d, tr.v) + ρ(d, tr.X)
end

"""
    an_element(tr::Translation)

Return some element of a translation.

### Input

- `tr` -- translation

### Output

An element in the translation.

### Notes

This function first asks for `an_element` function of the wrapped set, then
translates this element according to the given translation vector.
"""
function an_element(tr::Translation)
    return an_element(tr.X) + tr.v
end

function isboundedtype(::Type{<:Translation{N, S}}) where {N, S}
    return isboundedtype(S)
end

"""
    constraints_list(tr::Translation)

Return the list of constraints of the translation of a set.

### Input

- `tr` -- lazy translation of a polyhedron

### Output

The list of constraints of the translation.

### Notes

We assume that the set wrapped by the lazy translation `X` offers a method
`constraints_list(⋅)`.

### Algorithm

Let the translation be defined by the set of points `y` such that `y = x + v` for
all `x ∈ X`. Then, each defining halfspace `a⋅x ≤ b` is transformed to
`a⋅y ≤ b + a⋅v`.
"""
function constraints_list(tr::Translation)
    return _constraints_list_translation(tr.X, tr.v)
end

function _constraints_list_translation(X::LazySet, v::AbstractVector)
    @assert applicable(constraints_list, X) "this function requires that " *
        "the `constraints_list` method is applicable"

    constraints_X = constraints_list(X)
    constraints_TX = similar(constraints_X)
    @inbounds for (i, ci) in enumerate(constraints_X)
        constraints_TX[i] = HalfSpace(ci.a, ci.b + dot(ci.a, v))
    end
    return constraints_TX
end

"""
    ∈(x::AbstractVector, tr::Translation)

Check whether a given point is contained in the translation of a set.

### Input

- `x`  -- point/vector
- `tr` -- translation of a set

### Output

`true` iff ``x ∈ tr``.

### Algorithm

This implementation relies on the set membership function for the wrapped set
`tr.X`, since ``x ∈ X ⊕ v`` iff ``x - v ∈ X``.
"""
function ∈(x::AbstractVector, tr::Translation)
    return x - tr.v ∈ tr.X
end

"""
    linear_map(M::AbstractMatrix, tr::Translation)

Concrete linear map of a polyhedron in constraint representation.

### Input

- `M`  -- matrix
- `tr` -- translation of a set

### Output

A concrete set corresponding to the linear map.
The type of the result depends on the type of the set wrapped by `tr`.

### Algorithm

We compute `translate(linear_map(M, tr.X), M * tr.v)`.
"""
function linear_map(M::AbstractMatrix, tr::Translation)
    @assert dim(tr) == size(M, 2) "a linear map of size $(size(M)) cannot be " *
        "applied to a set of dimension $(dim(tr))"

    return translate(linear_map(M, tr.X), M * tr.v)
end

function concretize(tr::Translation)
    return translate(concretize(tr.X), tr.v)
end


text:
# ROW REDUCTION
function rowreduce_M{T}(M::Array{T},C,inds,fetched,qready)
  N = size(M,1)
  for i = 1:N
    i_ind = findfirst(inds,i)
    if i_ind !=0 # M ones
      v = M[:,i_ind]
      put!(C,v)
      M[i,i_ind] = 1
      for j = i+1:N
        M[j,i_ind] = 0
      end
      i_ind_p = i_ind + 1
    else
      v = fetch(C)::Vector{T}
      i_ind_p = fld(i-first(inds),step(inds))+2
    end
    put!(fetched,i)

    for k_ind = i_ind_p:length(inds) # M loop
      scale = M[i,k_ind]/v[i]
      M[i,k_ind] = scale
      for j = i+1:N
        M[j,k_ind] -= scale*v[j]
      end
    end

    take!(qready)
  end
  M
end

function rowreduce_S{T}(S::Array{T},C,inds,fetched,qready)
  N = size(S,1)
  for i = 1:N
    v = fetch(C)::Vector{T}
    put!(fetched,i)
    for k_ind in 1:length(inds) # S loop
      scale = S[i,k_ind]/v[i]
      S[i,k_ind] = scale
      for j = i+1:N
        S[j,k_ind] -= scale*v[j]
      end
    end
    take!(qready)
  end
  S
end

# BACK SUBSTITUTION
function backsub_S{T}(S::Array{T},C,inds,fetched,qready)
  # S = fetch(S_r)
  N = size(S,1)

  for i = N:-1:1
    v = fetch(C)::Vector{T}
    put!(fetched,i)
    for k_ind in 1:length(inds) # S loop
      scale = S[i,k_ind]
      for j = 1:i-1
        S[j,k_ind] -= scale*v[j]
      end
    end
    take!(qready)
  end
  S
end

function generateremoterefs(workerprocs)
  ar = Array(typeof(RemoteRef()),size(workerprocs))
  for (i,w) in enumerate(workerprocs)
    ar[i] = RemoteRef(w)
  end
  ar
end

function parallel_invert(M,workerprocs=workers();verbose=true)
  @assert size(M,1)==size(M,2)
  N = size(M,1)
  S = eye(eltype(M),N)
  nws = length(workerprocs)
  enum_workers = enumerate(workerprocs)

  verbose && println("Channels...")

  @compat C = RemoteRef(myid())
  @compat fetched_M = generateremoterefs(workerprocs)
  @compat fetched_S = generateremoterefs(workerprocs)
  @compat qready_M = generateremoterefs(workerprocs)
  @compat qready_S = generateremoterefs(workerprocs)

  # ROW REDUCTION
  verbose && println("Reducing rows...")
  @sync begin
    @async begin
      for i = 1:N
        @sync begin
          for pn = 1:nws
            @async begin
              take!(fetched_M[pn])
            end
            @async begin
              take!(fetched_S[pn])
            end
          end
        end
        take!(C)
        for pn = 1:nws
          put!(qready_M[pn],i)
          put!(qready_S[pn],i)
        end
      end
    end
    for (pn,pid) = enum_workers
      @async M[:,pn:nws:N] = remotecall_fetch(rowreduce_M,pid,M[:,pn:nws:N],C,pn:nws:N,fetched_M[pn],qready_M[pn])
      @async S[:,pn:nws:N] = remotecall_fetch(rowreduce_S,pid,S[:,pn:nws:N],C,pn:nws:N,fetched_S[pn],qready_S[pn])
    end
  end

  for pn = 1:nws
    close(fetched_M[pn])
    close(qready_M[pn])
  end

  verbose && println("Onto back substitution...")
  # BACK SUBSTITUTION
  @sync begin
    @async begin
      for i = N:-1:1
        put!(C,M[:,i])
        @sync begin
          for pn = 1:nws
            @async begin
              take!(fetched_S[pn])
            end
          end
        end
        take!(C)
        for pn = 1:nws
          put!(qready_S[pn],i)
        end
      end
    end
    for (pn,pid) = enum_workers
      @async S[:,pn:nws:N] = remotecall_fetch(backsub_S,pid,S[:,pn:nws:N],C,pn:nws:N,fetched_S[pn],qready_S[pn])
    end

  end
  close(C)
  for pn = 1:nws
    close(fetched_S[pn])
    close(qready_S[pn])
  end

  S
end


text:
mutable struct SymmetricTensor{T} <: AbstractArray{T,4}
    w::Int
    h::Int
    f::Int
    contents::Array{T,3}

    function SymmetricTensor{T}(w::Int, h::Int, f::Int) where T
        contents = Array{T}(undef, w, h, div(f * (f + 1), 2))
        new{T}(w, h, f, contents)
    end
end

function Base.size(o::SymmetricTensor)
    return (o.w, o.h, o.f, o.f)
end

function Base.getindex(o::SymmetricTensor, i::Int, j::Int, k::Int, l::Int)
    o.contents[i, j, convert_indices(o.f, order_indices(k, l)...)]
end

function Base.setindex!(o::SymmetricTensor, x, i::Int, j::Int, k::Int, l::Int)
    o.contents[i, j, convert_indices(o.f, order_indices(k, l)...)] = x
end

function order_indices(i::Int, j::Int)
    return j < i ? (j, i) : (i, j)
end

function convert_indices(n::Int, i::Int, j::Int)
    return n * (i - 1) + j - div(i * (i - 1), 2)
end

"""
```
descriptor = RegionCovariance(feature_img)
```
"""
struct RegionCovariance
    P::Array{Float64,3}
    Q::SymmetricTensor{Float64}
    nrows::Int
    ncols::Int
    nfeatures::Int

    function RegionCovariance(F::AbstractArray{Float64,3})
        nrows, ncols, nfeatures = length.(axes(F))
        P = featurewise_integral_images(F, nrows, ncols, nfeatures)
        Q = feature_product_integral_images(F, nrows, ncols, nfeatures)
        new(P, Q, nrows, ncols, nfeatures)
    end

    # Calculate the integral image of each feature layer of `F`. This is equivalent to
    # calling `integral_image(Fᵢ)` in a loop, where `Fᵢ` is a view of ith feature layer of
    # `F`, but reproducing the functionality here results in significantly fewer
    # allocations.
    function featurewise_integral_images(F::AbstractArray{Float64,3},
                                         nrows::Int, ncols::Int, nfeatures::Int)
        P = Array{Float64,3}(undef, nrows, ncols, nfeatures)

        @inbounds for f in 1:nfeatures
            for r in 1:nrows, c in 1:ncols
                P[r,c,f] = F[r,c,f]
                P[r,c,f] += r > 1 ? P[r-1,c,f] : 0.0
                P[r,c,f] += c > 1 ? P[r,c-1,f] : 0.0
                P[r,c,f] -= r > 1 && c > 1 ? P[r-1,c-1,f] : 0.0
            end
        end

        return P
    end

    # Given a feature image with 𝑛 features, where each pixel is a corresponds to a feature
    # vector 𝐱 ∈ ℝⁿ, calculate the 𝑛 × 𝑛 matrix 𝐴 for each 𝐱 such that 𝐴[𝑖,𝑗] = 𝐱ᵢ ⋅ 𝐱ⱼ, and
    # then the integral image of that matrix. Note that 𝐴 in this case is symmetric: a
    # custom data structure is used to leverage this symmetry and save space.
    function feature_product_integral_images(F::AbstractArray{Float64,3},
                                             nrows::Int, ncols::Int, nfeatures::Int)
        Q = SymmetricTensor{Float64}(nrows, ncols, nfeatures)

        @inbounds for fᵢ in 1:nfeatures
            for fⱼ in fᵢ:nfeatures
                for r in 1:nrows, c in 1:ncols
                    Q[r,c,fᵢ,fⱼ] = F[r,c,fᵢ] * F[r,c,fⱼ]
                    Q[r,c,fᵢ,fⱼ] += r > 1 ? Q[r-1,c,fᵢ,fⱼ] : 0.0
                    Q[r,c,fᵢ,fⱼ] += c > 1 ? Q[r,c-1,fᵢ,fⱼ] : 0.0
                    Q[r,c,fᵢ,fⱼ] -= r > 1 && c > 1 ? Q[r-1,c-1,fᵢ,fⱼ] : 0.0
                end
            end
        end

        return Q
    end
end

function covariance_matrix(d::RegionCovariance, rows::AbstractRange, cols::AbstractRange)
    𝐂 = Array{Float64}(undef, d.nfeatures, d.nfeatures)
    return covariance_matrix!(𝐂, d, rows, cols)
end

function covariance_matrix!(𝐂::AbstractArray{Float64,2}, d::RegionCovariance,
                            rows::AbstractRange, cols::AbstractRange)
    row₀, col₀ = first.((rows, cols))
    row₁, col₁ = last.((rows, cols))
    return covariance_matrix!(𝐂, d, row₀, col₀, row₁, col₁)
end

function covariance_matrix(d::RegionCovariance,
                           row₀::Integer, col₀::Integer,
                           row₁::Integer, col₁::Integer)
    𝐂 = Array{Float64}(undef, d.nfeatures, d.nfeatures)
    return covariance_matrix!(𝐂, d, row₀, col₀, row₁, col₁)
end

function covariance_matrix!(𝐂::AbstractArray{Float64,2}, d::RegionCovariance,
                            row₀::Integer, col₀::Integer,
                            row₁::Integer, col₁::Integer)
    region_size = (row₁-row₀+1)*(col₁-col₀+1)

    𝐐 = sum_feature_products_over_region(d.Q, Val{d.nfeatures}(), row₀, col₀, row₁, col₁)
    𝐩 = sum_over_region_by_feature(d.P, Val{d.nfeatures}(), row₀, col₀, row₁, col₁)
    result = (𝐐 - ((𝐩 * 𝐩') / region_size)) / (region_size-1)

    𝐂 .= result
end

function sum_feature_products_over_region(Q::SymmetricTensor{Float64}, ::Val{nfeatures},
                                          row₀::Integer, col₀::Integer, row₁::Integer,
                                          col₁::Integer) where nfeatures
    𝐐 = zeros(MMatrix{nfeatures,nfeatures,Float64})

    @inbounds for fᵢ in 1:nfeatures
        for fⱼ in fᵢ:nfeatures
            𝐐[fᵢ,fⱼ] = Q[row₁,col₁,fᵢ,fⱼ]
            𝐐[fᵢ,fⱼ] -= row₀ > 1 ? Q[row₀-1,col₁,fᵢ,fⱼ] : 0.0
            𝐐[fᵢ,fⱼ] -= col₀ > 1 ? Q[row₁,col₀-1,fᵢ,fⱼ] : 0.0
            𝐐[fᵢ,fⱼ] += row₀ > 1 && col₀ > 1 ? Q[row₀-1,col₀-1,fᵢ,fⱼ] : 0.0
            𝐐[fⱼ,fᵢ] = 𝐐[fᵢ,fⱼ]
        end
    end

    return 𝐐
end

function sum_over_region_by_feature(P::AbstractArray{Float64,3}, ::Val{nfeatures},
                                    row₀::Integer, col₀::Integer, row₁::Integer,
                                    col₁::Integer) where nfeatures
    𝐩 = zeros(MVector{nfeatures,Float64})

    @inbounds for f in 1:nfeatures
        𝐩[f] = P[row₁,col₁,f]
        𝐩[f] -= row₀ > 1 ? P[row₀-1,col₁,f] : 0.0
        𝐩[f] -= col₀ > 1 ? P[row₁,col₀-1,f] : 0.0
        𝐩[f] += row₀ > 1 && col₀ > 1 ? P[row₀-1,col₀-1,f] : 0.0
    end

    return 𝐩
end


text:
"A box in R^n."
struct BoxSpace{T <: Number,N} <: AbstractSpace
    low::Array{T,N}
    high::Array{T,N}
end

size(s::BoxSpace) = size(s.low)

"""
    BoxSpace(low::Number, high::Number, size::Tuple{Vararg{Int}}=(1,))

```julia
BoxSpace(-1, 1)
BoxSpace(-1, 1, (2,3))
```
"""
BoxSpace(low::Number, high::Number, size::Tuple{Vararg{Int}}=(1,)) = BoxSpace(fill(low, size), fill(high, size))

"""
    BoxSpace(low::Array{<:Number}, high::Array{<:Number})

```julia
BoxSpace([0, 0, 0], [1, 2, 3])
```
"""
BoxSpace(low::Array{<:Number}, high::Array{<:Number}) = size(low) == size(high) && BoxSpace(low, high)

==(x::BoxSpace, y::BoxSpace) = x.low == y.low && x.high == y.high

sample(s::BoxSpace) = map((l, h) -> l + rand() * (h - l), s.low, s.high)

occursin(xs::Array{<:Number}, s::BoxSpace) = size(s) == size(xs) &&
    all(map((a, b, c) -> a ≤ b ≤ c, s.low, xs, s.high))
occursin(x::Number, s::BoxSpace) = all(map((l, h) -> l ≤ x ≤ h, s.low, s.high))


text:
### The function eps() outputs the machine epsilon of the argument.
### Argument currently set to 1e-14.
### Resulting EPSILON is approx. 1.6e-30.
const EPSILON = eps(1e-14)

### EPSILON32 is 1f-30.
const EPSILON32 = "0.000000000000000000000000000001f"

### Euler–Mascheroni constant.
### γ = 0.5772156649015...
import Base.MathConstants.eulergamma

### Numerical approximation to the digamma function.
### Based on eq. (12), without looking at the accompanying source
### code, of: K. S. Kölbig, "Programs for computing the logarithm of
### the gamma function, and the digamma function, for complex
### argument," Computer Phys. Commun. vol. 4, pp. 221–226 (1972).
const DIGAMMA_c =
"""
inline float
digamma(float x)

		{
		float p = 0.0f;

		if (x < 7)
		{
			int n = 7 - floor(x);
			for (int v=1; v < n; v++)
				p -= 1 / (x + v);

			p -= 1 / x;
			x += n;
		}

		float t = 1 / x;
		p += log(x) - 0.5f * t;
		t *= t;
		p -= t * 0.08333333333333333f
				- 0.008333333333333333f * t
				+ 0.003968253968253968f * t*t
				- 0.004166666666666667f * t*t*t
				+ 0.007575757575757576f * t*t*t*t
				- 0.021092796092796094f * t*t*t*t*t
				+ 0.08333333333333333f * t*t*t*t*t*t
				- 0.4432598039215686f * t*t*t*t*t*t*t;

		return p;
		}
		"""

### Parallel Gauss-Jordan elimination (reduced row echelon form) algorithm for solving the linear system Ax=b.
### Partial pivoting unnecessary for symmetric positive-definite matrices. Ref,
### "Accuracy and Stability of Numerical Algorithms" Higham, 2002.
const LINSOLVE_c =
"""
inline void
linsolve(long K, long z, local float *A, local float *b)

		{
		for (long j=0; j<K-1; j++)
		{
			if ((j <= z) && (z < K-1))
			{
				float c = -A[K * j + (z + 1)] / A[K * j + j];

			    for (long l=j+1; l<K; l++)
			        A[K * l + (z + 1)] += c * A[K * l + j];

			    b[z + 1] += c * b[j];
			}
			barrier(CLK_LOCAL_MEM_FENCE);
		}

		for (long j=K-1; j>0; j--)
		{
			if (z < j)
				b[z] -= b[j] * A[K * j + z] / A[K * j + j];

			barrier(CLK_LOCAL_MEM_FENCE);
		}

		b[z] *= 1 / A[K * z + z];
		}
		"""

"Type alias for a vector of vectors."
VectorList{T} = Vector{Vector{T}}

"Type alias for a vector of matrices."
MatrixList{T} = Vector{Matrix{T}}

"Prevent overflow to ±Inf."
finite(x::Union{AbstractFloat, Array{<:AbstractFloat}}) = sign.(x) .* min.(abs.(x), floatmax.(x))

### LogSumExp function.
### Overflow safe.
import Distributions.logsumexp

function additive_logistic(x::Matrix{<:Real}; dims::Integer)
	"Additive logistic function of a real-valued matrix over the given dimension."
	"Overflow safe."

	if dims in [1,2]
		x = exp.(x .- maximum(x, dims=dims))
		x = x ./ sum(x, dims=dims)
	end

	return x
end

function additive_logistic(x::Vector{<:Real})
	"Additive logistic function of a real-valued vector."
	"Overflow safe."

	x = exp.(x .- maximum(x))
	x = x / sum(x)

	return x
end

function additive_logistic(x::Matrix{<:Real})
	"Additive logistic function of a real-valued matrix."
	"Overflow safe."

	x = exp.(x .- maximum(x))
	x = x / sum(x)

	return x
end

function isstochastic(P::Matrix{<:Real}; dims::Integer)
	"Check to see if X is a stochastic matrix."
	"if dims = 1, check for left stochastic matrix."
	"if dims = 2, check for right stochastic matrix."

	if dims == 1
		x = all([isprobvec(P[:,j]) for j in 1:size(P, 2)])
	elseif dims == 2
		x = all([isprobvec(P[i,:]) for i in 1:size(P, 1)])
	else
		x = P
	end

	return x
end

### Keep until pull request is merged.
import Distributions.xlogy
import Distributions.binomlogpdf
Distributions.xlogy(x::T, y::T) where {T<:Real} = x != zero(T) ? x * log(y) : zero(log(y))
Distributions.binomlogpdf(n::Real, p::Real, k::Real) = (isinteger(k) & (zero(k) <= k <= n)) ? convert(typeof(float(p)), loggamma(n + 1) - loggamma(k + 1) - loggamma(n - k + 1) + xlogy(k, p) + xlogy(n - k, 1 - p)) : convert(typeof(float(p)), -Inf)

### Keep until pull request is merged.
function Distributions.entropy(d::Dirichlet)
    α = d.alpha
    α0 = d.alpha0
    k = length(α)

    if k == 1
    	en = 0.0

    else
	    en = d.lmnB + (α0 - k) * digamma(α0)
	    for j in 1:k
	        @inbounds αj = α[j]
	        en -= (αj - 1.0) * digamma(αj)
	    end
	end

    return en
end


text:
function test()
  x = normal(0, 1)
  y = y + 5
  mean(y)
end

function test()
  x(ω) = normal(ω, 0, 1)
  mean(x)
end

struct SymbolicΩ
  x::Dunno
end

normal(ω::SymbolicΩ, args...) = (f = normal, args = args)
mean()

# Idea 1 is to basically do a symbolic exection
# Create a symbolic omega type
# redefine all ops on it
# +(x::Symbolic, y::Symbolic) =
# Then if I want to compute an expectaiton I analyze the graph


text:
begin
	using Plots
	using Printf
	using Unitful
	using DataFrames
	using CSV
	using Images
	using Interpolations
	using QuadGK
	using UnitfulEquivalences
	using Markdown
	using InteractiveUtils
	using PlutoUI
	using LsqFit
	using Distributions
	using Statistics
	using StatsBase
end
begin
	using Peaks
	using Glob
	using FFTW
	using DSP
end
import Unitful:
    nm, μm, mm, cm, m, km,
    mg, g, kg,
    ps, ns, μs, ms, s, minute, hr, d, yr, Hz, kHz, MHz, GHz,
    eV,
    μJ, mJ, J,
	μW, mW, W,
    A, N, mol, mmol, V, L, M

import PhysicalConstants.CODATA2018: N_A

function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end
lfi = ingredients("../src/LaserLab.jl")
PlutoUI.TableOfContents(title="Table of Contents", indent=true)

max_stars_repo_name: proscrite/LaserLab

max_issues_repo_name: proscrite/LaserLab

max_forks_repo_name: proscrite/LaserLab


text:
module UnionTyped
using Try
using TryExperimental
using TryExperimental: Maybe
f(xs) = Ok(xs) |> Try.and_then(xs -> trygetindex(xs, 1)) |> Maybe.ok
end  # module UnionTyped

module ConcretelyTyped
using Try
using TryExperimental
using TryExperimental: ConcreteResult, Maybe
function trygetfirst(xs)::ConcreteResult{eltype(xs),BoundsError}
    trygetindex(xs, 1)
end
f(xs) = Ok(xs) |> Try.and_then(trygetfirst) |> Maybe.ok
end  # module


text:
using JuliaInterpreter
using Test

@testset "core" begin
    @test JuliaInterpreter.is_quoted_type(QuoteNode(Int32), :Int32)
    @test !JuliaInterpreter.is_quoted_type(QuoteNode(Int32), :Int64)
    @test !JuliaInterpreter.is_quoted_type(QuoteNode(Int32(0)), :Int32)
    @test !JuliaInterpreter.is_quoted_type(Int32, :Int32)

    function buildexpr()
        items = [7, 3]
        ex = quote
            X = $items
            for x in X
                println(x)
            end
        end
        return ex
    end
    frame = JuliaInterpreter.enter_call(buildexpr)
    lines = JuliaInterpreter.framecode_lines(frame.framecode.src)
    # Test that the :copyast ends up on the same line as the println
    if isdefined(Base.IRShow, :show_ir_stmt)   # only works on Julia 1.6 and higher
        @test any(str->occursin(":copyast", str) && occursin("println", str), lines)
    end
end


text:
#module ResNetModule
export ResNet

import Knet
using Knet.Layers21: Conv, BatchNorm, Linear, Block, Add
using Knet.Ops21: relu, pool, mean
using Artifacts

"""
    ResNet(; nblocks=(2,2,2,2), block=ResNetBottleneck, groups=1, bottleneck=1, classes=1000)
    ResNet(name::String; pretrained=true)

Return a ResNet model based on keyword arguments or a name specifying a predefined
structure. Load pretrained weights if `pretrained=true`, randomly initialize otherwise. The
models satisfy the indexing interface, e.g. model[1] is the preprocessing layer and
model[1:end-1] omits the imagenet classification top.

A ResNet model `r` applied to an input of size `(width,height,channels=3,images)` returns a
`(classes,1)` vector of class scores. The model is a sequence of 6 functions: `r[1]` is the
input block, `r[2:5]` are four stages of residual blocks, and `r[6]` is the output
block. One can manipulate a ResNet model using the array interface, e.g. `r[1:5](img)` would
give the output of the last residual block, `pushfirst!(r, preprocess)` would add the
`preprocess` function to the beginning, etc.

The input and output blocks are the same for every ResNet model:

    Input:
        Conv(7×7, 3=>64, padding=3, stride=2, BatchNorm(), relu)
        Op(pool; window=3, stride=2, padding=1)

    Output:
        Op(pool; op=mean, window=typemax(Int))
        reshape2d
        Linear(classes, bias)

Stage i in 1:4 consists of `nblocks[i]` blocks of type `block` which can be
`ResNetBottleneck` (default) or `ResNetBasic`. `ResNetBasic` is a residual block with two
3×3 convolutions. `ResNetBottleneck` is a residual block with a 1×1 convolution that brings
the number of channels down, a 3×3 convolution (possibly grouped if `groups > 1`), and
another 1×1 convolution that brings the number of channels up. The keyword argument
`bottleneck` controls the ratio of the output channels to the intermediate channels. The
number of channels is doubled at each stage: blocks in the i'th stage have `2^(i+5)` output
channels for `ResNetBasic` and `2^(i+7)` output channels for `ResNetBottleneck`.  The
spatial resolution is halved at each stage: the first 3×3 convolution at each stage except
the first has `stride=2`.

The pretrained models come from
[torchvision](https://pytorch.org/vision/stable/models.html). A generic ResNet model can
handle any size input and any number of classes, but for pretrained models the training
settings should be replicated for best results: The number of classes is 1000 and the
expected size of the input is (224,224,3,N). The input images are resized to 256 pixels at
the short edge and center cropped to size (224,224). The pixel values are transformed from
the (0,1) range to:

    (rgb_pixel_0_1 .- [0.485, 0.456, 0.406]) ./ [0.229, 0.224, 0.225]

The class labels, as well as training, validation and test sets can be found at
[Kaggle](https://www.kaggle.com/c/imagenet-object-localization-challenge). Here are all the
pretrained models with name, settings, size in bytes, inference time compared to resnet50,
and top-1 validation accuracy.

    name              size   time  top1   ref1    settings
    ----              ----   ----  ----   ----    --------
    resnet18           45M   0.41  68.94  69.758  (nblocks=(2,2,2,2), block=ResNetBasic)
    resnet34           84M   0.57  72.76  73.314  (nblocks=(3,4,6,3), block=ResNetBasic)
    resnet50           98M   1.00  75.28  76.130  (nblocks=(3,4,6,3), bottleneck=4)
    resnet101         171M   1.58  76.83  77.374  (nblocks=(3,4,23,3), bottleneck=4)
    resnet152         231M   2.19  77.80  78.312  (nblocks=(3,8,36,3), bottleneck=4)
    wide_resnet50_2   264M   2.11  77.75  78.468  (nblocks=(3,4,6,3), bottleneck=2)
    wide_resnet101_2  485M   3.00  78.19  78.848  (nblocks=(3,4,23,3), bottleneck=2)
    resnext50_32x4d    96M   1.34  77.06  77.618  (nblocks=(3,4,6,3), groups=32, bottleneck=2)
    resnext101_32x8d  340M   3.45  78.69  79.312  (nblocks=(3,4,23,3), groups=32)

Note: The top1 accuracy figures are slightly different from the ones given by torchvision
(ref1), which is probably due to the differences in preprocessing, in particular `imresize`
in Julia gives different results compared to `Resize` in `torchvision.transforms`.

References:
* He, Kaiming et al. "Deep Residual Learning for Image Recognition." 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (2016): 770-778.
* Zagoruyko, Sergey and Nikos Komodakis. "Wide Residual Networks." ArXiv abs/1605.07146 (2016)
* Xie, Saining et al. "Aggregated Residual Transformations for Deep Neural Networks." 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (2017): 5987-5995.
* https://pytorch.org/vision/master/models.html

"""
function ResNet(; nblocks = (2,2,2,2), block = ResNetBottleneck, groups = 1, bottleneck = 1, classes = 1000)
    s = Block(; name="$block$nblocks")
    push!(s, Op(imagenet_preprocess; normalize=resnet_pt_normalize))
    push!(s, ResNetInput())
    x, y = 64, (block === ResNetBasic ? 64 : 256)
    for (stage, nblock) in enumerate(nblocks)
        if stage > 1; y *= 2; end
        b = y ÷ bottleneck
        blocks = Block(; name="Stage$stage")
        for iblock in 1:nblock
            stride = (stage > 1 && iblock == 1) ? 2 : 1
            push!(blocks, block(x, b, y; stride, groups))
            x = y
        end
        push!(s, blocks)
    end
    push!(s, ResNetOutput(y, classes))
    resnet_init(s)
end

function ResNetBottleneck(x, b, y; activation=relu, groups=1, padding=1, stride=1, o...)
    Add(
        Block(
            ConvBN(1, 1, x, b; activation),
            ConvBN(3, 3, b÷groups, b; activation, groups, padding, stride),
            ConvBN(1, 1, b, y),
        ),
        (x != y ? ConvBN(1, 1, x, y; stride) : identity);
        activation)
end

function ResNetBasic(x, b, y; stride=1, padding=1, activation=relu, o...)
    Add(
        Block(
            ConvBN(3, 3, x, b; activation, padding, stride),
            ConvBN(3, 3, b, y; padding),
        ),
        (x != y ? ConvBN(1, 1, x, y; stride) : identity);
        activation)
end

function ResNetInput()
    Block(
        ConvBN(7, 7, 3, 64; stride=2, padding=3, activation=relu),
        Op(pool; window=3, stride=2, padding=1);
        name = "Input"
    )
end

function ResNetOutput(xchannels, classes)
    Block(
        Op(pool; op=mean, window=typemax(Int)),
        reshape2d,
        Linear(xchannels, classes; binit=zeros);
        name = "Output"
    )
end

# Frequently used combo
ConvBN(x...; o...) = Conv(x...; o..., normalization=BatchNorm())

# Torchvision normalization
resnet_pt_normalize(x) = (x .- [0.485, 0.456, 0.406]) ./ [0.229, 0.224, 0.225]

# Run a single image so weights get initialized
resnet_init(m) = (m(Knet.atype(zeros(Float32,224,224,3,1))); m)

# Pretrained models from torchvision
resnetmodels = Dict{String,NamedTuple}(
    "resnet18" => (nblocks=(2,2,2,2), block=ResNetBasic),
    "resnet34" => (nblocks=(3,4,6,3), block=ResNetBasic),
    "resnet50" => (nblocks=(3,4,6,3), bottleneck=4),
    "resnet101" => (nblocks=(3,4,23,3), bottleneck=4),
    "resnet152" => (nblocks=(3,8,36,3), bottleneck=4),
    "wide_resnet50_2" => (nblocks=(3,4,6,3), bottleneck=2),
    "wide_resnet101_2" => (nblocks=(3,4,23,3), bottleneck=2),
    "resnext50_32x4d" => (nblocks=(3,4,6,3), groups=32, bottleneck=2),
    "resnext101_32x8d" => (nblocks=(3,4,23,3), groups=32),
)

function ResNet(s::String; pretrained=true)
    @assert haskey(resnetmodels, s)  "Please choose from known ResNet models:\n$(collect(keys(resnetmodels)))"
    kwargs = resnetmodels[s]
    model = ResNet(; kwargs...)
    pretrained && setweights!(model, joinpath(@artifact_str(s), "$s.jld2"))
    return model
end

#end # ResNetModule


text:
include("euler/euler.jl")
using .Primes: smallest_prime_factor
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int)::Int
    while true
        prime = smallest_prime_factor(n)
        if prime < n
            n ÷= prime
        else
            return n
        end
    end
end

compute(13_195)

compute(600_851_475_143)

@benchmark compute(600_851_475_143)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using Test, Logging, Random
using DecisionProgramming

rng = MersenneTwister(4)

@info "Testing random_diagram"
@test_throws DomainError random_diagram(rng, -1, 1, 1, 1, 1)
@test_throws DomainError random_diagram(rng, 1, -1, 1, 1, 1)
@test_throws DomainError random_diagram(rng, 0, 0, 1, 1, 1)
@test_throws DomainError random_diagram(rng, 1, 1, 0, 1, 1)
@test_throws DomainError random_diagram(rng, 1, 1, 1, 0, 1)
@test_throws DomainError random_diagram(rng, 1, 1, 1, 1, 0)

for (n_C, n_D) in [(1, 0), (0, 1)]
    rng = RandomDevice()
    C, D, V = random_diagram(rng, n_C, n_D, 1, 1, 1)
    @test isa(C, Vector{ChanceNode})
    @test isa(D, Vector{DecisionNode})
    @test isa(V, Vector{ValueNode})
    @test length(C) == n_C
    @test length(D) == n_D
    @test length(V) == 1
    @test all(!isempty(v.I_j) for v in V)
end

@info "Testing random States"
@test_throws DomainError States(rng, [0], 10)
@test isa(States(rng, [2, 3], 10), States)

@info "Testing random Probabilities"
S = States([2, 3, 2])
c = ChanceNode(3, [1, 2])
@test isa(Probabilities(rng, c, S; n_inactive=0), Probabilities)
@test isa(Probabilities(rng, c, S; n_inactive=1), Probabilities)
@test isa(Probabilities(rng, c, S; n_inactive=2*3*(2-1)), Probabilities)
@test_throws DomainError Probabilities(rng, c, S; n_inactive=2*3*(2-1)+1)

@info "Testing random Consequences"
S = States([2, 3])
v = ValueNode(3, [1, 2])
@test isa(Consequences(rng, v, S; low=-1.0, high=1.0), Consequences)
@test_throws DomainError Consequences(rng, v, S; low=1.1, high=1.0)

@info "Testing random LocalDecisionStrategy"
S = States([2, 3, 2])
d = DecisionNode(3, [1, 2])
@test isa(LocalDecisionStrategy(rng, d, S), LocalDecisionStrategy)


text:
# This file is helpful for creating jobfiles for a cluster
#   supercomputing environment.

# --- TASK VARIABLES ---
# how many job files to create (targets will be split in files)
# each file will correspond to one cluster or the ones specified
nnodes = 1
procs = 9
# for cluster (supercomputing) environments you need to specify queue
# each environment has different queues
# queue = "normal"
queue = "skx-normal"
# slurm command
cmd = "sbatch"
# julia path, usually in cluster must me specified, if julia is
# on a personal computer and its in the path just write julia
juliapath = "julia"
# juliapath = "/home1/05863/mgarciat/julia-1./julia"
# jobname used by slurm launcher
jobname = "rideaustin"
# max running time
maxjobtime = "24:00:00"

# -- MODEL FIT VARIABLES ---
# variables that will be passed to model_script.jl
workspace = pwd()  # for home computer
modelscript = joinpath(workspace, "4_hyperparameter_search.jl")
ngens = 48
gensize = 4
num_threads = 10
walltime = 300.0
# walltime = 300.0
cvsplits = 5
datadir = joinpath(workspace, "productivity_splits/")
# name for the resulting bash file
launcherfile = joinpath(workspace, "slurm_launcher.sh")
jobdir = joinpath(workspace, "job_scripts")
allocation = "Measuring-exposure-t"

# this is the header of the script
header = """#!/bin/bash
#
#SBATCH -N $nnodes
#SBATCH -n $procs
#SBATCH -p $queue
#SBATCH -J $jobname
#SBATCH -t $maxjobtime
#SBATCH -A $allocation
#
export WALLTIME=$(float(walltime))
export JULIAPATH="$juliapath"
export SCRIPT="$modelscript"
export NGENS=$ngens
export GENSIZE=$gensize
export CVSPLITS=$cvsplits
export JULIA_NUM_THREADS=$num_threads
#
srun \$JULIAPATH \$SCRIPT \$WALLTIME \$NGENS \$GENSIZE \$CVSPLITS"""

if isfile(launcherfile)
    rm(launcherfile)
end
if isdir(jobdir)
    rm(jobdir, recursive=true)
end
mkdir(jobdir)

datafiles = readdir(datadir)
N = length(datafiles)  # must be 2^tree_depth - 1 splits
open(launcherfile, "a") do launcherio
    num_files = if N % procs == 0
        N ÷ procs
    else
        N ÷ procs + 1
    end
    for j in 1:num_files
        # cd into folder and sbatch then go back
        jobfile = joinpath(jobdir, "job-$(j).sh")
        open(jobfile, "w") do jobio
            joblist = []
            for i in 1:procs
                idx = procs * (j - 1) + i
                if idx > N
                    break
                end
                file = datafiles[idx]
                push!(joblist, file)
            end
            joblist = join(joblist, " ")
            write(jobio, "$header $joblist\n")
        end
        write(launcherio, "$cmd \"$jobfile\"\n")
    end
end

max_stars_repo_name: mauriciogtec/gfen-reproduce

max_issues_repo_name: mauriciogtec/gfen-reproduce

max_forks_repo_name: mauriciogtec/gfen-reproduce


text:
# FT_STATUS (DWORD)
export FT_OK, FT_INVALID_HANDLE, FT_DEVICE_NOT_FOUND
export FT_DEVICE_NOT_OPENED, FT_IO_ERROR, FT_INSUFFICIENT_RESOURCES
export FT_INVALID_PARAMETER, FT_INVALID_BAUD_RATE, FT_DEVICE_NOT_OPENED_FOR_ERASE
export FT_DEVICE_NOT_OPENED_FOR_WRITE, FT_FAILED_TO_WRITE_DEVICE, FT_EEPROM_READ_FAILED
export FT_EEPROM_WRITE_FAILED, FT_EEPROM_ERASE_FAILED, FT_EEPROM_NOT_PRESENT
export FT_EEPROM_NOT_PROGRAMMED, FT_INVALID_ARGS, FT_NOT_SUPPORTED
export FT_OTHER_ERROR
export ftstatusdict

const FT_OK = 0
const FT_INVALID_HANDLE = 1
const FT_DEVICE_NOT_FOUND = 2
const FT_DEVICE_NOT_OPENED = 3
const FT_IO_ERROR = 4
const FT_INSUFFICIENT_RESOURCES = 5
const FT_INVALID_PARAMETER = 6
const FT_INVALID_BAUD_RATE = 7
const FT_DEVICE_NOT_OPENED_FOR_ERASE = 8
const FT_DEVICE_NOT_OPENED_FOR_WRITE = 9
const FT_FAILED_TO_WRITE_DEVICE = 10
const FT_EEPROM_READ_FAILED = 11
const FT_EEPROM_WRITE_FAILED = 12
const FT_EEPROM_ERASE_FAILED = 13
const FT_EEPROM_NOT_PRESENT = 14
const FT_EEPROM_NOT_PROGRAMMED = 15
const FT_INVALID_ARGS = 16
const FT_NOT_SUPPORTED = 17
const FT_OTHER_ERROR = 18

ftstatusdict = Dict(
  FT_OK                           => "FT_OK",
  FT_INVALID_HANDLE               => "FT_INVALID_HANDLE",
  FT_DEVICE_NOT_FOUND             => "FT_DEVICE_NOT_FOUND",
  FT_DEVICE_NOT_OPENED            => "FT_DEVICE_NOT_OPENED",
  FT_IO_ERROR                     => "FT_IO_ERROR",
  FT_INSUFFICIENT_RESOURCES       => "FT_INSUFFICIENT_RESOURCES",
  FT_INVALID_PARAMETER            => "FT_INVALID_PARAMETER",
  FT_INVALID_BAUD_RATE            => "FT_INVALID_BAUD_RATE",
  FT_DEVICE_NOT_OPENED_FOR_ERASE  => "FT_DEVICE_NOT_OPENED_FOR_ERASE",
  FT_DEVICE_NOT_OPENED_FOR_WRITE  => "FT_DEVICE_NOT_OPENED_FOR_WRITE",
  FT_FAILED_TO_WRITE_DEVICE       => "FT_FAILED_TO_WRITE_DEVICE",
  FT_EEPROM_READ_FAILED           => "FT_EEPROM_READ_FAILED",
  FT_EEPROM_WRITE_FAILED          => "FT_EEPROM_WRITE_FAILED",
  FT_EEPROM_ERASE_FAILED          => "FT_EEPROM_ERASE_FAILED",
  FT_EEPROM_NOT_PRESENT           => "FT_EEPROM_NOT_PRESENT",
  FT_EEPROM_NOT_PROGRAMMED        => "FT_EEPROM_NOT_PROGRAMMED",
  FT_INVALID_ARGS                 => "FT_INVALID_ARGS",
  FT_NOT_SUPPORTED                => "FT_NOT_SUPPORTED",
  FT_OTHER_ERROR                  => "FT_OTHER_ERROR")

type Ftd2xxError <: Exception
  ft_status :: UInt64
end
function Base.showerror(io::IO, e::Ftd2xxError)
  fts = e.ft_status
  errortext = get(ftstatusdict,fts,"Code Not Found")
  print(io, "Ftd2xxError: FT_STATUS=$(fts), $errortext")
end

function checkstatus(ft_status::Cuint)
  if ft_status != 0
    throw(Ftd2xxError(ft_status))
  end
  return 0
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAppsV1DaemonSetCondition <: SwaggerModel
    lastTransitionTime::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: lastTransitionTime
    message::Any # spec type: Union{ Nothing, String } # spec name: message
    reason::Any # spec type: Union{ Nothing, String } # spec name: reason
    status::Any # spec type: Union{ Nothing, String } # spec name: status
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiAppsV1DaemonSetCondition(;lastTransitionTime=nothing, message=nothing, reason=nothing, status=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1DaemonSetCondition, Symbol("lastTransitionTime"), lastTransitionTime)
        setfield!(o, Symbol("lastTransitionTime"), lastTransitionTime)
        validate_property(IoK8sApiAppsV1DaemonSetCondition, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        validate_property(IoK8sApiAppsV1DaemonSetCondition, Symbol("reason"), reason)
        setfield!(o, Symbol("reason"), reason)
        validate_property(IoK8sApiAppsV1DaemonSetCondition, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        validate_property(IoK8sApiAppsV1DaemonSetCondition, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiAppsV1DaemonSetCondition

const _property_map_IoK8sApiAppsV1DaemonSetCondition = Dict{Symbol,Symbol}(Symbol("lastTransitionTime")=>Symbol("lastTransitionTime"), Symbol("message")=>Symbol("message"), Symbol("reason")=>Symbol("reason"), Symbol("status")=>Symbol("status"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiAppsV1DaemonSetCondition = Dict{Symbol,String}(Symbol("lastTransitionTime")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("message")=>"String", Symbol("reason")=>"String", Symbol("status")=>"String", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiAppsV1DaemonSetCondition }) = collect(keys(_property_map_IoK8sApiAppsV1DaemonSetCondition))
Swagger.property_type(::Type{ IoK8sApiAppsV1DaemonSetCondition }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAppsV1DaemonSetCondition[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1DaemonSetCondition }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1DaemonSetCondition[property_name]

function check_required(o::IoK8sApiAppsV1DaemonSetCondition)
    (getproperty(o, Symbol("status")) === nothing) && (return false)
    (getproperty(o, Symbol("type")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAppsV1DaemonSetCondition }, name::Symbol, val)
end


text:
twice(x)=2x
half(x)=x/2
addone(x)=x+1

x=1
y=x |> twice |> addone |> half
@show y

z=begin
   x |>
   twice |>
   addone |>
   half
end

@show z

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Cassette
using InteractiveUtils
using MixedModeDebugger
using Test

using MixedModeDebugger: make_interpretted_mode_overdub, remove_interpretted_mode_overdub
@testset "make_interpretted_mode_overdub" begin
    @testset "Happy path" begin
        # these should not error.
        @testset "+(1)" begin
            overdubs_before = methods(Cassette.overdub)
            @test make_interpretted_mode_overdub(@which +(1)) isa Any
            @test methods(Cassette.overdub) != overdubs_before
            @test remove_interpretted_mode_overdub(@which +(1)) isa Any
            @test_broken methods(Cassette.overdub) == overdubs_before
        end

        @testset "+(1, 2)" begin
            overdubs_before = methods(Cassette.overdub)
            @test make_interpretted_mode_overdub(@which +(1, 2)) isa Any
            @test methods(Cassette.overdub) != overdubs_before
            @test remove_interpretted_mode_overdub(@which +(1, 2)) isa Any
            @test_broken methods(Cassette.overdub) == overdubs_before
        end

        @testset "rand(1)" begin
            overdubs_before = methods(Cassette.overdub)
            @test make_interpretted_mode_overdub(@which rand(1)) isa Any
            @test methods(Cassette.overdub) != overdubs_before
            @test remove_interpretted_mode_overdub(@which rand(1)) isa Any
            @test_broken methods(Cassette.overdub) == overdubs_before
        end

        @testset "rand(1:10, (10,20))" begin
            overdubs_before = methods(Cassette.overdub)
            @test make_interpretted_mode_overdub(@which rand(1:10, (10,20))) isa Any
            @test methods(Cassette.overdub) != overdubs_before
            @test remove_interpretted_mode_overdub(@which rand(1:10, (10,20))) isa Any
            @test_broken methods(Cassette.overdub) == overdubs_before
        end
    end
end


text:
"""
Partial function application.
"""
function partial(f, args...; kwargs...)
    function new_f(fargs...; fkwargs...)
        return f(args..., fargs...; kwargs..., fkwargs...)
    end
    return new_f
end

"""
Doesn't fail iff the two blocks are equivalent.
Blocks are equivalent if they represent equal time-series, even if the representations of
those time series are different (e.g. a range vs a vector)
"""
function test_equivalent(x::Block, y::Block)
    # Convert the blocks to a standard representation.
    x = Block(collect(x))
    y = Block(collect(y))
    @test isequal(x, y)
end

"""
Evaluate in several different ways, and ensure that they are all equivalent.
"""
function _evaluate(node::Node, t0::DateTime, t1::DateTime)
    # Standard evaluation in one go.
    block = evaluate(node, t0, t1)

    # Evaluation with a batch interval.
    for divisions in [1, 2, 3, 4, 5, 13]
        # We don't care exactly what batch interval we get, just pick something that roughly
        # works.
        batch_interval, _ = divrem(t1 - t0, divisions)
        test_block = evaluate(node, t0, t1; batch_interval)
        test_equivalent(test_block, block)
    end

    # Re-evaluation, and copying state.
    state = start_at([node], t0)
    copied_state = duplicate(state)
    evaluate_until!(state, t1)
    test_equivalent(only(state.evaluated_node_to_blocks[node]), block)
    evaluate_until!(copied_state, t1)
    test_equivalent(only(copied_state.evaluated_node_to_blocks[node]), block)

    return block
end

"""
Map the given function over each of the block's times & values.
"""
_map(f, block::Block) = Block([time => f(time, value) for (time, value) in block])

# Common functionality for testing binary operators that also perform alignment.
# Disable formatting, so as to permit more consistent layout of timeseries.
#! format: off
const b1 = Block([
    DateTime(2000, 1, 1) => 1,
    DateTime(2000, 1, 2) => 2,
    DateTime(2000, 1, 3) => 3,
    DateTime(2000, 1, 4) => 4,
])

const b2 = Block([
    DateTime(2000, 1, 2) => 5,
    DateTime(2000, 1, 3) => 6,
    DateTime(2000, 1, 5) => 8,
])

const b3 = Block([
    DateTime(2000, 1, 1) => 15,
])

const b4 = Block([
    DateTime(2000, 1, 1) => 1,
    DateTime(2000, 1, 2) => 2,
    DateTime(2000, 1, 3) => 3,
    DateTime(2000, 1, 4) => 4,
    DateTime(2000, 1, 5) => 5,
    DateTime(2000, 1, 6) => 6,
    DateTime(2000, 1, 7) => 7,
])

const b_boolean = Block([
    DateTime(2000, 1, 1) => true,
    DateTime(2000, 1, 2) => false,
    DateTime(2000, 1, 3) => true,
    DateTime(2000, 1, 4) => true,
])

const n1 = block_node(b1)
const n2 = block_node(b2)
const n3 = block_node(b3)
const n4 = block_node(b4)
const n_boolean = TimeDag.block_node(b_boolean)

const _T_START = DateTime(2000, 1, 1)
const _T_END = DateTime(2001, 1, 1)

_eval(n) = _evaluate(n, _T_START, _T_END)
# Evaluation when not wanting to perform tests.
_eval_fast(n) = evaluate(n, _T_START, _T_END)

"""
    _get_rand_svec_block(rng::AbstractRNG, dim::Int, n_obs::Int) -> Block

Get a block of value type `SVector{dim,Float64}`, of length `n_obs`, with random values.
"""
function _get_rand_svec_block(rng::AbstractRNG, dim::Int, n_obs::Int)
    times = _T_START:Day(1):(_T_START + Day(n_obs - 1))
    values = [SVector{dim}(rand(rng, dim)) for _ in 1:n_obs]
    return Block(times, values)
end

"""
    _get_rand_vec_block(rng::AbstractRNG, dim::Int, n_obs::Int) -> Block

Get a block of value type `Vector{Float64}`, of length `n_obs`, with random values.
"""
function _get_rand_vec_block(rng::AbstractRNG, dim::Int, n_obs::Int)
    times = _T_START:Day(1):(_T_START + Day(n_obs - 1))
    values = [rand(rng, dim) for _ in 1:n_obs]
    return Block(times, values)
end

"""
    _get_rand_block(rng::AbstractRNG, n_obs::Int) -> Block

Get a block of value type `Float64`, of length `n_obs`, with random values.
"""
function _get_rand_block(rng::AbstractRNG, n_obs::Int)
    times = _T_START:Day(1):(_T_START + Day(n_obs - 1))
    values = rand(rng, n_obs)
    return Block(times, values)
end

function _test_unary_op(f_timedag, block::Block, f=f_timedag; time_agnostic::Bool=true)
    ft = time_agnostic ? (_, x) -> f(x) : f

    # Check basic evaluation.
    value = first(block.values)
    node = block_node(block)
    @test _eval(f_timedag(node)) == _map(ft, block)

    if time_agnostic
        # Check constant propagation (only relevant for time agnostic nodes).
        @test f_timedag(constant(value)) == constant(f(value))
        @test f_timedag(constant(value)) === constant(f(value))
    end

    # Two instances of the NodeOp instance should compare equal for equal
    # type parameters.
    T = typeof(value)
    @test T != Float32
    for ta in (true, false)
        @test TimeDag.SimpleUnary{f,ta,T}() == TimeDag.SimpleUnary{f,ta,T}()
        @test TimeDag.SimpleUnary{f,ta,Float32}() != TimeDag.SimpleUnary{f,ta,T}()
    end
end

function _test_binary_op(f_timedag, f=f_timedag; time_agnostic::Bool=true)
    ft = time_agnostic ? (_, x, y) -> f(x, y) : f

    @testset "Common (fast) alignment" begin
        # Should apply for *any* user choice of alignment.
        for alignment in (UNION, LEFT, INTERSECT)
            n = f_timedag(n1, n1, alignment)
            block = _eval(n)
            @test block == Block([
                DateTime(2000, 1, 1) => ft(DateTime(2000, 1, 1), 1, 1),
                DateTime(2000, 1, 2) => ft(DateTime(2000, 1, 2), 2, 2),
                DateTime(2000, 1, 3) => ft(DateTime(2000, 1, 3), 3, 3),
                DateTime(2000, 1, 4) => ft(DateTime(2000, 1, 4), 4, 4),
            ])
            # For fast alignment, we expect *identical* timestamps as on the input.
            @test block.times === b1.times
        end
    end

    @testset "union alignment" begin
        n = f_timedag(n1, n2)
        @test n === f_timedag(n1, n2, UNION)
        block = _eval(n)
        @test block == Block([
            DateTime(2000, 1, 2) => ft(DateTime(2000, 1, 2), 2, 5),
            DateTime(2000, 1, 3) => ft(DateTime(2000, 1, 3), 3, 6),
            DateTime(2000, 1, 4) => ft(DateTime(2000, 1, 4), 4, 6),
            DateTime(2000, 1, 5) => ft(DateTime(2000, 1, 5), 4, 8),
        ])
    end

    @testset "intersect alignment" begin
        n = f_timedag(n1, n2, INTERSECT)

        @test _eval(n) == Block([
            DateTime(2000, 1, 2) => ft(DateTime(2000, 1, 2), 2, 5),
            DateTime(2000, 1, 3) => ft(DateTime(2000, 1, 3), 3, 6),
        ])
    end

    @testset "left alignment" begin
        n = f_timedag(n1, n2, LEFT)
        @test _eval(n) == Block([
            DateTime(2000, 1, 2) => ft(DateTime(2000, 1, 2), 2, 5),
            DateTime(2000, 1, 3) => ft(DateTime(2000, 1, 3), 3, 6),
            DateTime(2000, 1, 4) => ft(DateTime(2000, 1, 4), 4, 6),
        ])

        # Catch edge-case in which there was a bug.
        @test _eval(f_timedag(n2, n3, LEFT)) == Block([
            DateTime(2000, 1, 2) => ft(DateTime(2000, 1, 2), 5, 15),
            DateTime(2000, 1, 3) => ft(DateTime(2000, 1, 3), 6, 15),
            DateTime(2000, 1, 5) => ft(DateTime(2000, 1, 5), 8, 15),
        ])
        @test _eval(f_timedag(n3, n2, LEFT)) == Block{Int64}()
    end

    if time_agnostic
        @testset "constant propagation" begin
            value = 2.0  # valid for all ops we're currently testing.
            @test f_timedag(constant(value), constant(value)) === constant(f(value, value))
            @test f_timedag(value, constant(value)) === constant(f(value, value))
            @test f_timedag(constant(value), value) === constant(f(value, value))
        end
    end

    # Two instances of the NodeOp instance should compare equal for equal
    # type parameters.
    T = Float64
    for A in (UnionAlignment, IntersectAlignment, LeftAlignment)
        for ta in (true, false)
            @test TimeDag.SimpleBinary{f,ta,T,A}() == TimeDag.SimpleBinary{f,ta,T,A}()
            @test TimeDag.SimpleBinary{f,ta,Float32,A}() != TimeDag.SimpleBinary{f,ta,T,A}()
        end
    end

    @testset "initial values" begin
        # Initial value of wrong type.
        @test_throws ArgumentError f_timedag(n1, n2; initial_values=(-1, -2.0))

        @testset "union alignment" begin
            n = f_timedag(n1, n2; initial_values=(-1, -2))
            @test n === f_timedag(n1, n2, UNION; initial_values=(-1, -2))

            @test _eval(n) == Block([
                DateTime(2000, 1, 1) => ft(DateTime(2000, 1, 1), 1, -2),
                DateTime(2000, 1, 2) => ft(DateTime(2000, 1, 2), 2, 5),
                DateTime(2000, 1, 3) => ft(DateTime(2000, 1, 3), 3, 6),
                DateTime(2000, 1, 4) => ft(DateTime(2000, 1, 4), 4, 6),
                DateTime(2000, 1, 5) => ft(DateTime(2000, 1, 5), 4, 8),
            ])

            n = f_timedag(n2, n1; initial_values=(-1, -2))
            @test _eval(n) == Block([
                DateTime(2000, 1, 1) => ft(DateTime(2000, 1, 1), -1, 1),
                DateTime(2000, 1, 2) => ft(DateTime(2000, 1, 2), 5, 2),
                DateTime(2000, 1, 3) => ft(DateTime(2000, 1, 3), 6, 3),
                DateTime(2000, 1, 4) => ft(DateTime(2000, 1, 4), 6, 4),
                DateTime(2000, 1, 5) => ft(DateTime(2000, 1, 5), 8, 4),
            ])
        end

        @testset "left alignment" begin
            n = f_timedag(n1, n2, LEFT; initial_values=(-1, -2))
            # Initial left value should be ignored.
            @test n === f_timedag(n1, n2, LEFT; initial_values=(-42, -2))
            @test _eval(n) == Block([
                DateTime(2000, 1, 1) => ft(DateTime(2000, 1, 1), 1, -2),
                DateTime(2000, 1, 2) => ft(DateTime(2000, 1, 2), 2, 5),
                DateTime(2000, 1, 3) => ft(DateTime(2000, 1, 3), 3, 6),
                DateTime(2000, 1, 4) => ft(DateTime(2000, 1, 4), 4, 6),
            ])
        end

        # Intersect alignment
        @testset "intersect alignment" begin
            n = f_timedag(n1, n2, INTERSECT; initial_values=(-1, -2))
            # Initial values should always be ignored.
            @test n === f_timedag(n1, n2, INTERSECT)
        end
    end
end

function _test_ternary_op(f_timedag, f=f_timedag; time_agnostic::Bool=true)
    ft = time_agnostic ? (_, x, y, z) -> f(x, y, z) : f

    @testset "Common (fast) alignment" begin
        # Should apply for *any* user choice of alignment.
        for alignment in (UNION, LEFT, INTERSECT)
            n = f_timedag(n1, n1, n1, alignment)
            block = _eval(n)
            @test block == Block([
                DateTime(2000, 1, 1) => ft(DateTime(2000, 1, 1), 1, 1, 1),
                DateTime(2000, 1, 2) => ft(DateTime(2000, 1, 2), 2, 2, 2),
                DateTime(2000, 1, 3) => ft(DateTime(2000, 1, 3), 3, 3, 3),
                DateTime(2000, 1, 4) => ft(DateTime(2000, 1, 4), 4, 4, 4),
            ])
            # For fast alignment, we expect *identical* timestamps as on the input.
            @test block.times === b1.times
        end
    end

    @testset "union alignment" begin
        n = f_timedag(n1, n2, n3)
        @test n === f_timedag(n1, n2, n3, UNION)
        block = _eval(n)
        @test block == Block([
            DateTime(2000, 1, 2) => ft(DateTime(2000, 1, 2), 2, 5, 15),
            DateTime(2000, 1, 3) => ft(DateTime(2000, 1, 3), 3, 6, 15),
            DateTime(2000, 1, 4) => ft(DateTime(2000, 1, 4), 4, 6, 15),
            DateTime(2000, 1, 5) => ft(DateTime(2000, 1, 5), 4, 8, 15),
        ])
    end

    @testset "intersect alignment" begin
        n = f_timedag(n1, n2, n4, INTERSECT)

        @test _eval(n) == Block([
            DateTime(2000, 1, 2) => ft(DateTime(2000, 1, 2), 2, 5, 2),
            DateTime(2000, 1, 3) => ft(DateTime(2000, 1, 3), 3, 6, 3),
        ])
    end

    @testset "left alignment" begin
        n = f_timedag(n1, n2, n3, LEFT)
        @test _eval(n) == Block([
            DateTime(2000, 1, 2) => ft(DateTime(2000, 1, 2), 2, 5, 15),
            DateTime(2000, 1, 3) => ft(DateTime(2000, 1, 3), 3, 6, 15),
            DateTime(2000, 1, 4) => ft(DateTime(2000, 1, 4), 4, 6, 15),
        ])

        # Catch edge-case in which there was a bug.
        @test _eval(f_timedag(n2, n3, n1, LEFT)) == Block([
            DateTime(2000, 1, 2) => ft(DateTime(2000, 1, 2), 5, 15, 2),
            DateTime(2000, 1, 3) => ft(DateTime(2000, 1, 3), 6, 15, 3),
            DateTime(2000, 1, 5) => ft(DateTime(2000, 1, 5), 8, 15, 4),
        ])
        @test _eval(f_timedag(n3, n2, n1, LEFT)) == Block{Int64}()
    end

    if time_agnostic
        @testset "constant propagation" begin
            value = 2.0  # valid for all ops we're currently testing.
            cv = constant(value)
            cf = constant(f(value, value, value))
            @test f_timedag(cv, cv, cv) === cf
            @test f_timedag(cv, value, value) === cf
            @test f_timedag(value, cv, value) === cf
            @test f_timedag(value, value, cv) === cf
            @test f_timedag(cv, cv, value) === cf
            @test f_timedag(cv, value, cv) === cf
            @test f_timedag(value, cv, cv) === cf
        end
    end

    # Two instances of the NodeOp instance should compare equal for equal
    # type parameters.
    T = Float64
    for A in (UnionAlignment, IntersectAlignment, LeftAlignment)
        for ta in (true, false)
            @test TimeDag.SimpleNary{f,ta,3,T,A}() == TimeDag.SimpleNary{f,ta,3,T,A}()
            @test TimeDag.SimpleNary{f,ta,3,Float32,A}() != TimeDag.SimpleNary{f,ta,3,T,A}()
        end
    end

    @testset "initial values" begin
        # Initial value of wrong type.
        @test_throws ArgumentError f_timedag(n1, n2, n3; initial_values=(-1, -2.0, -3))

        @testset "union alignment" begin
            n = f_timedag(n1, n2, n3; initial_values=(-1, -2, -3))
            @test n === f_timedag(n1, n2, n3, UNION; initial_values=(-1, -2, -3))

            @test _eval(n) == Block([
                DateTime(2000, 1, 1) => ft(DateTime(2000, 1, 1), 1, -2, 15),
                DateTime(2000, 1, 2) => ft(DateTime(2000, 1, 2), 2, 5, 15),
                DateTime(2000, 1, 3) => ft(DateTime(2000, 1, 3), 3, 6, 15),
                DateTime(2000, 1, 4) => ft(DateTime(2000, 1, 4), 4, 6, 15),
                DateTime(2000, 1, 5) => ft(DateTime(2000, 1, 5), 4, 8, 15),
            ])

            n = f_timedag(n2, n1, n3; initial_values=(-1, -2, -3))
            @test _eval(n) == Block([
                DateTime(2000, 1, 1) => ft(DateTime(2000, 1, 1), -1, 1, 15),
                DateTime(2000, 1, 2) => ft(DateTime(2000, 1, 2), 5, 2, 15),
                DateTime(2000, 1, 3) => ft(DateTime(2000, 1, 3), 6, 3, 15),
                DateTime(2000, 1, 4) => ft(DateTime(2000, 1, 4), 6, 4, 15),
                DateTime(2000, 1, 5) => ft(DateTime(2000, 1, 5), 8, 4, 15),
            ])
        end

        @testset "left alignment" begin
            n = f_timedag(n1, n2, n3, LEFT; initial_values=(-1, -2, -3))
            # TODO Support this.
            # # Initial left value should be ignored.
            # @test n === f_timedag(n1, n2, n3, LEFT; initial_values=(-42, -2, -3))
            @test _eval(n) == Block([
                DateTime(2000, 1, 1) => ft(DateTime(2000, 1, 1), 1, -2, 15),
                DateTime(2000, 1, 2) => ft(DateTime(2000, 1, 2), 2, 5, 15),
                DateTime(2000, 1, 3) => ft(DateTime(2000, 1, 3), 3, 6, 15),
                DateTime(2000, 1, 4) => ft(DateTime(2000, 1, 4), 4, 6, 15),
            ])
        end

        # Intersect alignment
        @testset "intersect alignment" begin
            n = f_timedag(n1, n2, n3, INTERSECT; initial_values=(-1, -2, -3))
            # Initial values should always be ignored.
            @test n === f_timedag(n1, n2, n3, INTERSECT)
        end
    end
end

#! format: on


text:
baremodule Historic

macro define end
# macro define_static end

function clear end
function events end
function flattable end
function recorders end

# Tools
function taskprefix end
function taskdata end
function taskid end
function objid end
function uniqueid end

const defaultprefix = taskprefix
const defaultdata = taskdata

module Internal

import ..Historic: @define
using ..Historic: Historic

using Base.Meta: isexpr
using Logging: Logging, AbstractLogger, current_logger, handle_message
using UUIDs: UUID, uuid4

include("BlockLinkedLists.jl")
using .BlockLinkedLists: BlockLinkedList, blocks

include("utils.jl")
include("define.jl")
include("recording.jl")
include("logging.jl")
include("defaultlogger.jl")
include("id.jl")
include("tools.jl")
# include("blockinglogger.jl")
include("accessors.jl")
include("tables.jl")

function __init__()
    init_unique_id()
end

end  # module Internal

@define Scratch
using .Scratch: @record

Internal.define_docstrings()

end  # baremodule Historic


text:
module LossTests

using Photon, Test

function test_distance(loss)
    p = randn(10,8)
    y = randn(10,8)
    r = loss(p,y)
    @test r > 0
    @test loss(p,y) == loss(y,p)

    r = loss(p,p)
    @test r == 0

    p = randn(10,1,8)
    y = randn(1,10,8)
    r = loss(p,y)
    @test r > 0

    p = KorA(randn(10,8))
    y = KorA(randn(10,8))
    r = loss(p,y)
    @test loss(p,y) == loss(y,p)
    @test r > 0

    r = loss(p,p)
    @test r == 0

end

function test_nondistance(loss)
    p = randn(10,8)
    y = randn(10,8)
    r = loss(p,y)
    @test r isa Number

    p = randn(10,1,8)
    y = randn(1,10,8)
    r = loss(p,y)
    @test r isa Number

    p = KorA(randn(10,8))
    y = KorA(randn(10,8))
    r = loss(p,y)
    @test r isa Number

    @test loss(p,p) <= loss(p,y)

end

function test_classification(loss)
    p = Float32.(rand(0:0.0001:1.0,10,16))
    y = Float32.(rand(0:1, 10,16))
    r = loss(p, y)
    @test r > 0

    r = loss(KorA(p), KorA((y)))
    @test r > 0

end

@testset "Loss" begin
    test_distance(L1Loss())
    test_distance(L1Loss(reduction=sum))
    test_distance(L2Loss())
    test_distance(LNLoss(3))

    test_nondistance(PseudoHuberLoss())

    test_classification(HingeLoss())
    test_classification(CrossEntropyLoss())
    test_classification(BCELoss())
    test_classification(FocalLoss())
end

end


text:
export Array_, similar_, copy_

struct CacheKey{T}
    name::Symbol
    l::Int
    c::Int

    CacheKey(n::Symbol, A::AT, sz::NTuple{M,Int} = size(A)) where AT <: AbstractArray{T,M} where {T,M} =
        new{AT}(n, prod(sz), checksum(sz))
    CacheKey(n::Symbol, ::Val{T}, sz::NTuple{M,Int}) where {T,M} = new{Array{T,M}}(n, prod(sz), checksum(sz))
end

checksum(size::NTuple{N,Int}) where N = sum(ntuple(i -> size[i] * i^2, Val(N)))

using LRUCache

const cache = LRU{CacheKey, AbstractArray}(;maxsize = 100) # very crude: fixed size, for now

# function Base.getindex(lru::LRU, key::CacheKey{AT}) where AT
#     node = lru.ht[key]
#     LRUCache.move_to_front!(lru.q, node)
#     return node.v::AT # to improve type stability?
# end

"""
    similar_(name, A)     ≈ similar(A)
    Array_{T}(name, size) ≈ Array{T}(undef, size)
New arrays for intermediate results, drawn from an LRU cache when `length(A) >= 2000`.
The cache's key uses `name::Symbol` as well as type & size to ensure different uses don't collide.

    copy_(name, A) = copyto!(similar_(name, A), A)
Just like that.
"""
similar_(A::AbstractArray) = similar_(:array_, A)

function similar_(name::Symbol, A::TA)::TA where TA<:AbstractArray
    if length(A) < 2000
        similar(A)
    else
        get(cache, CacheKey(name, A)) do
            similar(A)
        end
    end
end

struct Array_{T} end

@doc @doc(similar_)
Array_(any...) = Array_{Float64}(any...)

Array_{T}(sz::Vararg{Int}) where {T} =  Array_{T}(:Array_, sz)

Array_{T}(name::Symbol, sz::Vararg{Int}) where {T} =  Array_{T}(name, sz)

Array_{T}(sz::Tuple) where {T} = Array{T}(:Array_, sz)

function Array_{T}(name::Symbol, sz::NTuple{N,Int}) where {T,N}
    key = CacheKey(name, Val(T), sz)
    if prod(sz) < 2000
        Array{T}(undef, sz);
    else
        get(cache, key) do
            Array{T, N}(undef, sz)
        end
    end
end

@doc @doc(similar_)
copy_(A::AbstractArray) = copy_(:copy_, A)

copy_(name::Symbol, A::AbstractArray) = copyto!(similar_(name, A), A)


text:
module SMLMFrameConnection

using SMLMData
using Hungarian
using NearestNeighbors
using Optim
using StatsBase
using Statistics

include("structdefinitions.jl")
include("precluster.jl")
include("defineidealFC.jl")
include("organizeclusters.jl")
include("computeclusterinfo.jl")
include("estimateparams.jl")
include("estimatedensities.jl")
include("create_costmatrix.jl")
include("solveLAP.jl")
include("linkclusters.jl")
include("compress_connectID.jl")
include("combinelocalizations.jl")
include("connectlocalizations.jl")
include("frameconnect.jl")

end