text:
# Definition of Boore & Thompson (2012) constant values to be used within subsequent functions
const m_ii_bt12 = [ 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0 ]
const r_jj_bt12 = [ 2.00, 3.17, 5.02, 7.96, 12.62, 20.00, 31.70, 50.24, 79.62, 126.19, 200.00, 317.00, 502.40, 796.20, 1262.00 ]
const num_m_ii_bt12 = 9
const num_r_jj_bt12 = 15

const coefs_wna_bt12 = [ 4.0    2.00  8.4312e-01 -2.8671e-02 2.0  1.7316e+00  1.1695e+00  2.1671e+00  9.6224e-01  9.9929e-01  9.4336e-01   4.0    2.00;
4.5    2.00  8.3064e-01 -5.3355e-08 2.0  1.6962e+00  1.3236e+00  2.0308e+00  9.5517e-01  1.0176e+00  9.6580e-01   4.5    2.00;
5.0    2.00  8.9701e-01 -1.0484e-02 2.0  1.7181e+00  1.2797e+00  1.9295e+00  1.0163e+00  1.0290e+00  9.6175e-01   5.0    2.00;
5.5    2.00  8.3574e-01 -1.2031e-06 2.0  1.3859e+00  1.2015e+00  1.9904e+00  8.1829e-01  1.0523e+00  9.8336e-01   5.5    2.00;
6.0    2.00  8.4485e-01 -1.2609e-02 2.0  1.2391e+00  9.9399e-01  2.0532e+00  7.4884e-01  1.0681e+00  9.7567e-01   6.0    2.00;
6.5    2.00  8.4348e-01 -1.4520e-06 2.0  1.1331e+00  9.3441e-01  2.0328e+00  7.1150e-01  1.0680e+00  9.7535e-01   6.5    2.00;
7.0    2.00  8.0078e-01 -1.0231e-03 2.0  9.9676e-01  8.0084e-01  2.1466e+00  5.9286e-01  1.0908e+00  9.7381e-01   7.0    2.00;
7.5    2.00  8.3352e-01 -2.0697e-03 2.0  9.4601e-01  8.2133e-01  2.1052e+00  6.1011e-01  1.0849e+00  9.5386e-01   7.5    2.00;
8.0    2.00  1.1818e+00 -3.7875e-01 2.0  7.6114e-01  2.4868e+00  3.4051e+00  5.1616e-01  1.1002e+00  9.3188e-01   8.0    2.00;
4.0    3.17  8.2110e-01 -4.7357e-03 2.0  1.8627e+00  1.6118e+00  2.0128e+00  9.6557e-01  1.0024e+00  9.6134e-01   4.0    3.17;
4.5    3.17  8.7601e-01 -5.9903e-07 2.0  1.8200e+00  1.4626e+00  1.9542e+00  1.0076e+00  1.0074e+00  9.5429e-01   4.5    3.17;
5.0    3.17  8.3971e-01 -3.0439e-03 2.0  1.5964e+00  1.3635e+00  1.9761e+00  8.9640e-01  1.0425e+00  9.7057e-01   5.0    3.17;
5.5    3.17  8.3455e-01 -2.8341e-06 2.0  1.3116e+00  1.1068e+00  2.0414e+00  7.7989e-01  1.0458e+00  9.8125e-01   5.5    3.17;
6.0    3.17  8.1645e-01 -8.3202e-08 2.0  1.1908e+00  9.8968e-01  2.0549e+00  6.9953e-01  1.0647e+00  9.7531e-01   6.0    3.17;
6.5    3.17  8.2440e-01 -2.2922e-19 2.0  1.1212e+00  9.2690e-01  2.0537e+00  6.7673e-01  1.0722e+00  9.7217e-01   6.5    3.17;
7.0    3.17  8.3242e-01 -4.0859e-03 2.0  1.0919e+00  1.0439e+00  2.0355e+00  6.6335e-01  1.0850e+00  9.6470e-01   7.0    3.17;
7.5    3.17  8.3241e-01 -8.3435e-03 2.0  1.1067e+00  1.1521e+00  1.9907e+00  6.4715e-01  1.0951e+00  9.4210e-01   7.5    3.17;
8.0    3.17  8.3994e-01 -2.6168e-02 2.0  9.6494e-01  9.0967e-01  2.2436e+00  5.8696e-01  1.1028e+00  9.2279e-01   8.0    3.17;
4.0    5.02  8.6329e-01 -4.0771e-02 2.0  1.8792e+00  1.7825e+00  2.0071e+00  9.5979e-01  1.0150e+00  9.6497e-01   4.0    5.02;
4.5    5.02  8.4825e-01 -5.8380e-05 2.0  1.8285e+00  1.6611e+00  1.9249e+00  9.6711e-01  1.0252e+00  9.7713e-01   4.5    5.02;
5.0    5.02  8.4446e-01 -3.4661e-22 2.0  1.5463e+00  1.4665e+00  1.9448e+00  8.6771e-01  1.0329e+00  9.7563e-01   5.0    5.02;
5.5    5.02  8.4978e-01 -6.2734e-09 2.0  1.5126e+00  1.4004e+00  1.9072e+00  8.4997e-01  1.0460e+00  9.7378e-01   5.5    5.02;
6.0    5.02  8.2129e-01 -4.5730e-06 2.0  1.2115e+00  1.0578e+00  2.0028e+00  7.0987e-01  1.0598e+00  9.7654e-01   6.0    5.02;
6.5    5.02  8.1544e-01 -1.9160e-08 2.0  1.0668e+00  8.4143e-01  2.0735e+00  6.5158e-01  1.0751e+00  9.7467e-01   6.5    5.02;
7.0    5.02  8.2414e-01 -4.0369e-03 2.0  1.0467e+00  7.1242e-01  2.1586e+00  6.2836e-01  1.0794e+00  9.5870e-01   7.0    5.02;
7.5    5.02  8.9473e-01 -5.3853e-02 2.0  1.0888e+00  1.4138e+00  2.0212e+00  6.7776e-01  1.0900e+00  9.5025e-01   7.5    5.02;
8.0    5.02  9.2147e-01 -8.9899e-02 2.0  9.5008e-01  1.1607e+00  2.4017e+00  6.2635e-01  1.0980e+00  9.2819e-01   8.0    5.02;
4.0    7.96  8.4779e-01 -2.8493e-02 2.0  2.0858e+00  2.6221e+00  1.8501e+00  9.5690e-01  1.0194e+00  9.7959e-01   4.0    7.96;
4.5    7.96  8.4677e-01 -2.1554e-14 2.0  1.9701e+00  2.2219e+00  1.7755e+00  9.8498e-01  1.0316e+00  9.9715e-01   4.5    7.96;
5.0    7.96  8.8283e-01 -4.1263e-07 2.0  1.9400e+00  1.8389e+00  1.8013e+00  1.0235e+00  1.0334e+00  9.7377e-01   5.0    7.96;
5.5    7.96  8.2785e-01 -8.1587e-17 2.0  1.4890e+00  1.4475e+00  1.9052e+00  8.1571e-01  1.0561e+00  9.8353e-01   5.5    7.96;
6.0    7.96  8.3393e-01 -4.2583e-08 2.0  1.4016e+00  1.4232e+00  1.8975e+00  7.8526e-01  1.0658e+00  9.7835e-01   6.0    7.96;
6.5    7.96  8.1429e-01 -2.8150e-05 2.0  1.0955e+00  9.3363e-01  2.0547e+00  6.6504e-01  1.0786e+00  9.8445e-01   6.5    7.96;
7.0    7.96  8.0350e-01 -3.1805e-06 2.0  1.0903e+00  1.0120e+00  2.0825e+00  6.3175e-01  1.0930e+00  9.6779e-01   7.0    7.96;
7.5    7.96  8.6455e-01 -4.9497e-02 2.0  1.0672e+00  1.3965e+00  2.0600e+00  6.2719e-01  1.0924e+00  9.5518e-01   7.5    7.96;
8.0    7.96  9.3068e-01 -1.0354e-01 2.0  1.0295e+00  1.7393e+00  2.1326e+00  6.3596e-01  1.0994e+00  9.2609e-01   8.0    7.96;
4.0   12.62  8.4343e-01 -9.4015e-03 2.0  1.7138e+00  2.7357e+00  1.8332e+00  8.7570e-01  1.0193e+00  9.9529e-01   4.0   12.62;
4.5   12.62  8.2603e-01 -1.5203e-02 2.0  1.7515e+00  2.3023e+00  1.8093e+00  8.6816e-01  1.0462e+00  9.9553e-01   4.5   12.62;
5.0   12.62  8.5078e-01 -6.8471e-17 2.0  1.6943e+00  1.8025e+00  1.7993e+00  9.1795e-01  1.0461e+00  1.0026e+00   5.0   12.62;
5.5   12.62  8.4406e-01 -4.2873e-07 2.0  1.6936e+00  1.6701e+00  1.7942e+00  8.9001e-01  1.0584e+00  9.8314e-01   5.5   12.62;
6.0   12.62  8.1873e-01 -9.7554e-08 2.0  1.3379e+00  1.3058e+00  1.9370e+00  7.4590e-01  1.0649e+00  9.8476e-01   6.0   12.62;
6.5   12.62  8.2130e-01 -1.1179e-06 2.0  1.0666e+00  9.5527e-01  2.0441e+00  6.5972e-01  1.0706e+00  9.9066e-01   6.5   12.62;
7.0   12.62  8.3901e-01 -3.8867e-03 2.0  1.0721e+00  9.7602e-01  2.0319e+00  6.7214e-01  1.0773e+00  9.7088e-01   7.0   12.62;
7.5   12.62  8.8154e-01 -5.3242e-02 2.0  1.0746e+00  1.4629e+00  2.0001e+00  6.4928e-01  1.0882e+00  9.4526e-01   7.5   12.62;
8.0   12.62  9.6254e-01 -1.5468e-01 2.0  8.4972e-01  1.1369e+00  2.6445e+00  5.5070e-01  1.0967e+00  9.2600e-01   8.0   12.62;
4.0   20.00  8.4020e-01 -4.0622e-03 2.0  1.7122e+00  3.7111e+00  1.7781e+00  8.3369e-01  1.0231e+00  1.0055e+00   4.0   20.00;
4.5   20.00  8.0074e-01 -1.1163e-03 2.0  1.6050e+00  2.5984e+00  1.7981e+00  8.0776e-01  1.0514e+00  1.0107e+00   4.5   20.00;
5.0   20.00  8.1957e-01 -2.6783e-08 2.0  1.5537e+00  1.9525e+00  1.8541e+00  8.1185e-01  1.0503e+00  9.9790e-01   5.0   20.00;
5.5   20.00  8.4597e-01 -2.1761e-28 2.0  1.6384e+00  1.8613e+00  1.7709e+00  8.6430e-01  1.0543e+00  9.8634e-01   5.5   20.00;
6.0   20.00  8.1936e-01 -2.8792e-08 2.0  1.2967e+00  1.3933e+00  1.9226e+00  7.3484e-01  1.0646e+00  9.9737e-01   6.0   20.00;
6.5   20.00  8.2149e-01 -2.1626e-07 2.0  1.2651e+00  1.1833e+00  1.8904e+00  7.2337e-01  1.0798e+00  9.7962e-01   6.5   20.00;
7.0   20.00  8.2885e-01 -1.4510e-02 2.0  1.0514e+00  8.6325e-01  2.0854e+00  6.3459e-01  1.0858e+00  9.7187e-01   7.0   20.00;
7.5   20.00  8.8599e-01 -6.7109e-02 2.0  1.0820e+00  1.4223e+00  2.1483e+00  6.5140e-01  1.0931e+00  9.5410e-01   7.5   20.00;
8.0   20.00  9.6561e-01 -1.5069e-01 2.0  9.2494e-01  1.6108e+00  2.3995e+00  5.7877e-01  1.0942e+00  9.3791e-01   8.0   20.00;
4.0   31.70  8.1172e-01 -4.8496e-04 2.0  1.5154e+00  4.2066e+00  1.7868e+00  7.5495e-01  1.0411e+00  1.0121e+00   4.0   31.70;
4.5   31.70  7.8971e-01 -1.1171e-07 2.0  1.3881e+00  2.8525e+00  1.9129e+00  7.1011e-01  1.0446e+00  1.0101e+00   4.5   31.70;
5.0   31.70  7.8907e-01 -3.0999e-07 2.0  1.3945e+00  1.9681e+00  1.8882e+00  7.2307e-01  1.0590e+00  1.0093e+00   5.0   31.70;
5.5   31.70  8.2863e-01 -4.3712e-29 2.0  1.5918e+00  1.8609e+00  1.7700e+00  8.3632e-01  1.0636e+00  1.0051e+00   5.5   31.70;
6.0   31.70  8.1063e-01 -2.5156e-07 2.0  1.3122e+00  1.4358e+00  1.8632e+00  7.3257e-01  1.0734e+00  1.0013e+00   6.0   31.70;
6.5   31.70  8.1517e-01 -5.7539e-06 2.0  1.1768e+00  1.1135e+00  1.8888e+00  6.9805e-01  1.0796e+00  9.8922e-01   6.5   31.70;
7.0   31.70  8.1577e-01 -1.2460e-04 2.0  1.1862e+00  1.0549e+00  1.9416e+00  6.8195e-01  1.0868e+00  9.6233e-01   7.0   31.70;
7.5   31.70  9.2526e-01 -1.1605e-01 2.0  1.0984e+00  1.6537e+00  2.1768e+00  6.4739e-01  1.0974e+00  9.5360e-01   7.5   31.70;
8.0   31.70  9.0070e-01 -1.0408e-01 2.0  9.3630e-01  1.4596e+00  2.2897e+00  5.6532e-01  1.0993e+00  9.3947e-01   8.0   31.70;
4.0   50.24  7.4376e-01 -9.0204e-03 2.0  1.1760e+00  5.3724e+00  1.9952e+00  5.7376e-01  1.0611e+00  1.0292e+00   4.0   50.24;
4.5   50.24  7.6692e-01 -2.8995e-07 2.0  1.2262e+00  3.1344e+00  1.9104e+00  6.3632e-01  1.0607e+00  1.0303e+00   4.5   50.24;
5.0   50.24  8.0410e-01 -1.7249e-07 2.0  1.2658e+00  2.0493e+00  1.8728e+00  7.1205e-01  1.0610e+00  1.0300e+00   5.0   50.24;
5.5   50.24  8.0289e-01 -7.3756e-10 2.0  1.3033e+00  1.9082e+00  1.8283e+00  7.1257e-01  1.0611e+00  1.0228e+00   5.5   50.24;
6.0   50.24  8.2058e-01 -3.3962e-07 2.0  1.4215e+00  1.7148e+00  1.7918e+00  7.7238e-01  1.0701e+00  9.9697e-01   6.0   50.24;
6.5   50.24  8.1998e-01 -3.7615e-08 2.0  1.3171e+00  1.5287e+00  1.8113e+00  7.3336e-01  1.0743e+00  9.8469e-01   6.5   50.24;
7.0   50.24  8.0401e-01 -1.8401e-03 2.0  1.2092e+00  1.1668e+00  1.9378e+00  6.7879e-01  1.0950e+00  9.7028e-01   7.0   50.24;
7.5   50.24  9.2124e-01 -1.1109e-01 2.0  1.0821e+00  1.7469e+00  2.2082e+00  6.4069e-01  1.0886e+00  9.5945e-01   7.5   50.24;
8.0   50.24  9.1213e-01 -1.0459e-01 2.0  1.0303e+00  1.7188e+00  2.1460e+00  6.1959e-01  1.0951e+00  9.3272e-01   8.0   50.24;
4.0   79.62  7.1877e-01 -8.4424e-07 2.0  9.8189e-01  5.7882e+00  1.9315e+00  4.8869e-01  1.0633e+00  1.0405e+00   4.0   79.62;
4.5   79.62  7.6806e-01 -7.7865e-10 2.0  1.1500e+00  4.2410e+00  1.8852e+00  5.8860e-01  1.0545e+00  1.0258e+00   4.5   79.62;
5.0   79.62  8.1221e-01 -1.9339e-06 2.0  1.3867e+00  3.2310e+00  1.7712e+00  7.2172e-01  1.0557e+00  1.0225e+00   5.0   79.62;
5.5   79.62  8.1179e-01 -3.6814e-15 2.0  1.3673e+00  2.3510e+00  1.7695e+00  7.3065e-01  1.0617e+00  1.0209e+00   5.5   79.62;
6.0   79.62  8.3805e-01 -3.9863e-07 2.0  1.4750e+00  1.7874e+00  1.7191e+00  8.0489e-01  1.0630e+00  1.0006e+00   6.0   79.62;
6.5   79.62  8.3461e-01 -3.2932e-04 2.0  1.5359e+00  1.6907e+00  1.7029e+00  8.1570e-01  1.0730e+00  9.7817e-01   6.5   79.62;
7.0   79.62  8.2516e-01 -1.0574e-03 2.0  1.4617e+00  1.6598e+00  1.7075e+00  7.8372e-01  1.0840e+00  9.7763e-01   7.0   79.62;
7.5   79.62  9.2678e-01 -1.2826e-01 2.0  1.0954e+00  1.8061e+00  2.2617e+00  6.3883e-01  1.0935e+00  9.5951e-01   7.5   79.62;
8.0   79.62  9.4781e-01 -1.3027e-01 2.0  1.1123e+00  1.9835e+00  2.2682e+00  6.5983e-01  1.0897e+00  9.3568e-01   8.0   79.62;
4.0  126.19  7.0065e-01 -3.0866e-03 2.0  7.7514e-01  4.3880e+00  2.1942e+00  3.9751e-01  1.0665e+00  1.0430e+00   4.0  126.19;
4.5  126.19  7.3944e-01 -6.1421e-03 2.0  8.4700e-01  3.1720e+00  2.1861e+00  4.6495e-01  1.0631e+00  1.0421e+00   4.5  126.19;
5.0  126.19  8.0513e-01 -3.1300e-07 2.0  1.4173e+00  3.7696e+00  1.7451e+00  7.2746e-01  1.0661e+00  1.0375e+00   5.0  126.19;
5.5  126.19  8.3644e-01 -1.2787e-02 2.0  1.6253e+00  3.0751e+00  1.6687e+00  8.0811e-01  1.0660e+00  1.0188e+00   5.5  126.19;
6.0  126.19  8.4568e-01 -1.0589e-02 2.0  1.6317e+00  2.5450e+00  1.6207e+00  8.3600e-01  1.0663e+00  1.0196e+00   6.0  126.19;
6.5  126.19  8.4142e-01 -9.2432e-03 2.0  1.6367e+00  2.2281e+00  1.6330e+00  8.4571e-01  1.0759e+00  1.0072e+00   6.5  126.19;
7.0  126.19  8.3694e-01 -1.8592e-12 2.0  1.4260e+00  1.8587e+00  1.6980e+00  7.8752e-01  1.0719e+00  9.9019e-01   7.0  126.19;
7.5  126.19  9.7787e-01 -1.5087e-01 2.0  1.3198e+00  2.7095e+00  2.0948e+00  7.4851e-01  1.0802e+00  9.6889e-01   7.5  126.19;
8.0  126.19  9.5551e-01 -1.3687e-01 2.0  1.0502e+00  1.6420e+00  2.2945e+00  6.6223e-01  1.0834e+00  9.5884e-01   8.0  126.19;
4.0  200.00  7.2408e-01 -3.7872e-06 2.0  7.9315e-01  3.3863e+00  2.1194e+00  4.3168e-01  1.0731e+00  1.0494e+00   4.0  200.00;
4.5  200.00  7.6883e-01 -1.0649e-10 2.0  8.4489e-01  3.0451e+00  2.0022e+00  5.1530e-01  1.0674e+00  1.0491e+00   4.5  200.00;
5.0  200.00  8.0446e-01 -1.1419e-07 2.0  1.2079e+00  3.5405e+00  1.7311e+00  6.7846e-01  1.0708e+00  1.0434e+00   5.0  200.00;
5.5  200.00  8.5125e-01 -2.0805e-02 2.0  1.6828e+00  3.7940e+00  1.5949e+00  8.1298e-01  1.0634e+00  1.0157e+00   5.5  200.00;
6.0  200.00  8.8936e-01 -3.4678e-02 2.0  1.9176e+00  3.6113e+00  1.5457e+00  9.0903e-01  1.0583e+00  1.0253e+00   6.0  200.00;
6.5  200.00  8.8360e-01 -4.3347e-02 2.0  2.1189e+00  3.0810e+00  1.5627e+00  9.3733e-01  1.0723e+00  1.0070e+00   6.5  200.00;
7.0  200.00  8.8401e-01 -4.3273e-02 2.0  2.2411e+00  3.1543e+00  1.5785e+00  9.5961e-01  1.0768e+00  9.8898e-01   7.0  200.00;
7.5  200.00  9.7709e-01 -1.3800e-01 2.0  1.3833e+00  2.6413e+00  2.0144e+00  7.8650e-01  1.0724e+00  9.7395e-01   7.5  200.00;
8.0  200.00  9.4766e-01 -1.2599e-01 2.0  1.0937e+00  1.7927e+00  2.1485e+00  6.8255e-01  1.0817e+00  9.6134e-01   8.0  200.00;
4.0  317.00  7.3621e-01 -2.6124e-07 2.0  6.7866e-01  1.5442e+00  2.3208e+00  4.0993e-01  1.0759e+00  1.0570e+00   4.0  317.00;
4.5  317.00  7.7046e-01 -2.0803e-06 2.0  8.4695e-01  2.4074e+00  1.9286e+00  5.2494e-01  1.0786e+00  1.0462e+00   4.5  317.00;
5.0  317.00  8.3255e-01 -8.9857e-03 2.0  1.2272e+00  3.5393e+00  1.6207e+00  6.9898e-01  1.0666e+00  1.0366e+00   5.0  317.00;
5.5  317.00  9.3328e-01 -9.4425e-02 2.0  1.7388e+00  5.6367e+00  1.5939e+00  8.2772e-01  1.0637e+00  1.0337e+00   5.5  317.00;
6.0  317.00  9.6786e-01 -1.0575e-01 2.0  2.0371e+00  4.4326e+00  1.5835e+00  9.3047e-01  1.0545e+00  1.0163e+00   6.0  317.00;
6.5  317.00  9.9191e-01 -1.2806e-01 2.0  2.2631e+00  4.7905e+00  1.5624e+00  9.7756e-01  1.0573e+00  1.0188e+00   6.5  317.00;
7.0  317.00  9.8463e-01 -1.2911e-01 2.0  2.4172e+00  4.1328e+00  1.5968e+00  1.0248e+00  1.0671e+00  1.0121e+00   7.0  317.00;
7.5  317.00  9.8211e-01 -1.4721e-01 2.0  1.4960e+00  2.3251e+00  2.0687e+00  8.3602e-01  1.0789e+00  9.7502e-01   7.5  317.00;
8.0  317.00  9.8036e-01 -1.3181e-01 2.0  1.4355e+00  2.1624e+00  2.0674e+00  8.5037e-01  1.0726e+00  9.6315e-01   8.0  317.00;
4.0  502.40  7.7145e-01 -2.1903e-08 2.0  7.5501e-01  1.2254e+00  1.9918e+00  4.8005e-01  1.0784e+00  1.0568e+00   4.0  502.40;
4.5  502.40  8.3924e-01 -1.7176e-02 2.0  8.4240e-01  1.7566e+00  1.7682e+00  5.8451e-01  1.0631e+00  1.0452e+00   4.5  502.40;
5.0  502.40  8.6317e-01 -3.0913e-02 2.0  1.2712e+00  3.3879e+00  1.4455e+00  7.1360e-01  1.0659e+00  1.0382e+00   5.0  502.40;
5.5  502.40  9.3729e-01 -1.1001e-01 2.0  1.9503e+00  6.6024e+00  1.4210e+00  8.4644e-01  1.0762e+00  1.0447e+00   5.5  502.40;
6.0  502.40  9.8756e-01 -1.2780e-01 2.0  2.3398e+00  5.8967e+00  1.4040e+00  9.5484e-01  1.0593e+00  1.0275e+00   6.0  502.40;
6.5  502.40  1.0255e+00 -1.6704e-01 2.0  2.7393e+00  5.7808e+00  1.5652e+00  1.0339e+00  1.0626e+00  1.0108e+00   6.5  502.40;
7.0  502.40  1.0233e+00 -1.6031e-01 2.0  2.9463e+00  4.5741e+00  1.5529e+00  1.0997e+00  1.0652e+00  1.0006e+00   7.0  502.40;
7.5  502.40  1.0200e+00 -1.2401e-01 2.0  3.0000e+00  3.9155e+00  1.6445e+00  1.1560e+00  1.0461e+00  9.6435e-01   7.5  502.40;
8.0  502.40  1.0002e+00 -1.1999e-01 2.0  2.5675e+00  3.3955e+00  1.6785e+00  1.1070e+00  1.0599e+00  9.6587e-01   8.0  502.40;
4.0  796.20  8.8270e-01 -8.3018e-02 2.0  1.2361e+00  3.2178e+00  1.2860e+00  6.5438e-01  1.0883e+00  1.0605e+00   4.0  796.20;
4.5  796.20  9.0131e-01 -9.0819e-02 2.0  1.4035e+00  3.7970e+00  1.3032e+00  7.1226e-01  1.0843e+00  1.0593e+00   4.5  796.20;
5.0  796.20  9.3558e-01 -8.3350e-02 2.0  1.5771e+00  4.0563e+00  1.2214e+00  7.8376e-01  1.0603e+00  1.0343e+00   5.0  796.20;
5.5  796.20  9.8071e-01 -1.2639e-01 2.0  2.1823e+00  5.4690e+00  1.2828e+00  8.9979e-01  1.0625e+00  1.0402e+00   5.5  796.20;
6.0  796.20  1.0292e+00 -1.6994e-01 2.0  2.9999e+00  6.5631e+00  1.2973e+00  1.0071e+00  1.0608e+00  1.0173e+00   6.0  796.20;
6.5  796.20  1.0716e+00 -2.0385e-01 2.0  3.0000e+00  6.0142e+00  1.4500e+00  1.0802e+00  1.0575e+00  1.0181e+00   6.5  796.20;
7.0  796.20  1.1170e+00 -2.2503e-01 2.0  3.0000e+00  5.5908e+00  1.6436e+00  1.1310e+00  1.0466e+00  1.0011e+00   7.0  796.20;
7.5  796.20  1.1243e+00 -2.3116e-01 2.0  3.0000e+00  4.7524e+00  1.7099e+00  1.1844e+00  1.0490e+00  9.9987e-01   7.5  796.20;
8.0  796.20  1.0006e+00 -1.1062e-01 2.0  2.8163e+00  3.0643e+00  1.6161e+00  1.2033e+00  1.0543e+00  9.8424e-01   8.0  796.20;
4.0 1262.00  1.1142e+00 -2.9561e-01 2.0  1.6692e+00  6.1510e+00  1.3141e+00  7.6115e-01  1.0829e+00  1.0502e+00   4.0 1262.00;
4.5 1262.00  1.1450e+00 -2.8700e-01 2.0  1.9553e+00  6.0141e+00  1.2838e+00  8.6221e-01  1.0634e+00  1.0392e+00   4.5 1262.00;
5.0 1262.00  1.1431e+00 -3.0083e-01 2.0  2.2941e+00  6.8268e+00  1.2609e+00  8.9057e-01  1.0726e+00  1.0475e+00   5.0 1262.00;
5.5 1262.00  1.1467e+00 -2.7310e-01 2.0  2.5178e+00  6.9746e+00  1.3125e+00  9.6257e-01  1.0554e+00  1.0372e+00   5.5 1262.00;
6.0 1262.00  1.1623e+00 -2.9893e-01 2.0  3.0000e+00  7.2278e+00  1.4351e+00  1.0457e+00  1.0629e+00  1.0252e+00   6.0 1262.00;
6.5 1262.00  1.2000e+00 -3.2059e-01 2.0  3.0000e+00  7.6028e+00  1.5529e+00  1.0983e+00  1.0535e+00  1.0295e+00   6.5 1262.00;
7.0 1262.00  1.2000e+00 -3.0723e-01 2.0  3.0000e+00  6.1282e+00  1.6668e+00  1.1508e+00  1.0466e+00  1.0024e+00   7.0 1262.00;
7.5 1262.00  1.2000e+00 -3.0499e-01 2.0  3.0000e+00  5.0982e+00  1.7273e+00  1.2181e+00  1.0481e+00  1.0067e+00   7.5 1262.00;
8.0 1262.00  1.2000e+00 -3.0348e-01 2.0  3.0000e+00  4.2696e+00  1.8291e+00  1.2610e+00  1.0492e+00  9.8681e-01   8.0 1262.00 ]

const coefs_ena_bt12 = [ 4.0    2.00  7.4499e-01 -3.0772e-02 2.0  1.3928e+00  8.3838e-01  2.4979e+00  6.7955e-01  1.0253e+00  9.6354e-01  4.0    2.00;
4.5    2.00  7.6659e-01 -3.0413e-02 2.0  1.3636e+00  9.7771e-01  2.3632e+00  6.9419e-01  1.0538e+00  9.6656e-01  4.5    2.00;
5.0    2.00  8.5088e-01 -2.7815e-06 2.0  1.2579e+00  9.8032e-01  2.0749e+00  7.8838e-01  1.0598e+00  9.5727e-01  5.0    2.00;
5.5    2.00  7.6469e-01 -1.0292e-07 2.0  1.1153e+00  8.5585e-01  2.1955e+00  6.0514e-01  1.0745e+00  9.7063e-01  5.5    2.00;
6.0    2.00  8.0539e-01 -4.1307e-17 2.0  1.0631e+00  8.3666e-01  2.1792e+00  6.3713e-01  1.0776e+00  9.5917e-01  6.0    2.00;
6.5    2.00  8.2454e-01 -1.1132e-16 2.0  1.0240e+00  7.5717e-01  2.1935e+00  6.4412e-01  1.0957e+00  9.4398e-01  6.5    2.00;
7.0    2.00  8.4196e-01 -8.3632e-03 2.0  9.8648e-01  7.3577e-01  2.2221e+00  6.4268e-01  1.1001e+00  9.1936e-01  7.0    2.00;
7.5    2.00  8.0918e-01 -1.9895e-06 2.0  9.1407e-01  7.9186e-01  2.2626e+00  5.7595e-01  1.1080e+00  9.0618e-01  7.5    2.00;
8.0    2.00  1.1574e+00 -4.1628e-01 2.0  7.1549e-01  3.3332e+00  4.1754e+00  4.0981e-01  1.1114e+00  8.7461e-01  8.0    2.00;
4.0    3.17  7.4249e-01 -3.4061e-02 2.0  1.3792e+00  8.1130e-01  2.5371e+00  6.6759e-01  1.0251e+00  9.6351e-01  4.0    3.17;
4.5    3.17  7.7349e-01 -2.2472e-02 2.0  1.3817e+00  9.9520e-01  2.3073e+00  7.2181e-01  1.0535e+00  9.6624e-01  4.5    3.17;
5.0    3.17  8.5089e-01 -4.4519e-07 2.0  1.2735e+00  9.9594e-01  2.0666e+00  7.9499e-01  1.0601e+00  9.5733e-01  5.0    3.17;
5.5    3.17  7.6277e-01 -3.8418e-08 2.0  1.1145e+00  8.4433e-01  2.2020e+00  6.0264e-01  1.0744e+00  9.7069e-01  5.5    3.17;
6.0    3.17  8.0431e-01 -1.3389e-17 2.0  1.0649e+00  8.3603e-01  2.1806e+00  6.3657e-01  1.0778e+00  9.5924e-01  6.0    3.17;
6.5    3.17  8.2327e-01 -9.8165e-21 2.0  1.0255e+00  7.5555e-01  2.1958e+00  6.4310e-01  1.0954e+00  9.4426e-01  6.5    3.17;
7.0    3.17  8.4445e-01 -1.2431e-02 2.0  9.8210e-01  7.3887e-01  2.2377e+00  6.3896e-01  1.1000e+00  9.1991e-01  7.0    3.17;
7.5    3.17  8.0631e-01 -4.7208e-07 2.0  9.0758e-01  7.6090e-01  2.2867e+00  5.6936e-01  1.1083e+00  9.0648e-01  7.5    3.17;
8.0    3.17  1.1238e+00 -3.3574e-01 2.0  7.3859e-01  2.2882e+00  3.4862e+00  4.8200e-01  1.1122e+00  8.7493e-01  8.0    3.17;
4.0    5.02  7.4005e-01 -3.5818e-02 2.0  1.3817e+00  8.0165e-01  2.5506e+00  6.6482e-01  1.0248e+00  9.6350e-01  4.0    5.02;
4.5    5.02  7.7668e-01 -1.8238e-02 2.0  1.3961e+00  1.0010e+00  2.2827e+00  7.3750e-01  1.0531e+00  9.6663e-01  4.5    5.02;
5.0    5.02  8.5014e-01 -3.3214e-07 2.0  1.2845e+00  1.0033e+00  2.0625e+00  7.9911e-01  1.0594e+00  9.5759e-01  5.0    5.02;
5.5    5.02  7.6180e-01 -2.1336e-17 2.0  1.1175e+00  8.4272e-01  2.2037e+00  6.0313e-01  1.0745e+00  9.7063e-01  5.5    5.02;
6.0    5.02  8.0269e-01 -2.0009e-17 2.0  1.0693e+00  8.3770e-01  2.1815e+00  6.3637e-01  1.0776e+00  9.5939e-01  6.0    5.02;
6.5    5.02  8.2171e-01 -3.0971e-31 2.0  1.0280e+00  7.5486e-01  2.1982e+00  6.4214e-01  1.0955e+00  9.4444e-01  6.5    5.02;
7.0    5.02  8.4613e-01 -1.5745e-02 2.0  9.8072e-01  7.4313e-01  2.2508e+00  6.3615e-01  1.0995e+00  9.1966e-01  7.0    5.02;
7.5    5.02  8.0493e-01 -6.6301e-08 2.0  9.0662e-01  7.5179e-01  2.2949e+00  5.6713e-01  1.1082e+00  9.0712e-01  7.5    5.02;
8.0    5.02  1.1064e+00 -3.1800e-01 2.0  7.4773e-01  2.1428e+00  3.4018e+00  4.8628e-01  1.1116e+00  8.7545e-01  8.0    5.02;
4.0    7.96  7.3601e-01 -3.8687e-02 2.0  1.3940e+00  7.9664e-01  2.5636e+00  6.6310e-01  1.0243e+00  9.6349e-01  4.0    7.96;
4.5    7.96  7.8586e-01 -7.4996e-03 2.0  1.4162e+00  1.0088e+00  2.2240e+00  7.7432e-01  1.0527e+00  9.6619e-01  4.5    7.96;
5.0    7.96  8.4919e-01 -2.4409e-07 2.0  1.2958e+00  1.0061e+00  2.0606e+00  8.0357e-01  1.0588e+00  9.5731e-01  5.0    7.96;
5.5    7.96  7.6152e-01 -1.7248e-17 2.0  1.1247e+00  8.4775e-01  2.2017e+00  6.0650e-01  1.0740e+00  9.7116e-01  5.5    7.96;
6.0    7.96  8.0101e-01 -5.6871e-34 2.0  1.0761e+00  8.4182e-01  2.1812e+00  6.3750e-01  1.0776e+00  9.5969e-01  6.0    7.96;
6.5    7.96  8.1921e-01 -6.4799e-35 2.0  1.0327e+00  7.5446e-01  2.2011e+00  6.4087e-01  1.0955e+00  9.4462e-01  6.5    7.96;
7.0    7.96  8.4706e-01 -2.0159e-02 2.0  9.7927e-01  7.4729e-01  2.2700e+00  6.3117e-01  1.0997e+00  9.1996e-01  7.0    7.96;
7.5    7.96  8.0254e-01 -1.0004e-19 2.0  9.0843e-01  7.4615e-01  2.3023e+00  5.6478e-01  1.1079e+00  9.0698e-01  7.5    7.96;
8.0    7.96  1.0991e+00 -3.1060e-01 2.0  7.5231e-01  2.0950e+00  3.3597e+00  4.8859e-01  1.1116e+00  8.7612e-01  8.0    7.96;
4.0   12.62  7.0000e-01 -2.4374e-01 2.0  1.7607e+00  6.9945e+00  2.9837e+00  4.0530e-01  1.0529e+00  1.0123e+00  4.0   12.62;
4.5   12.62  8.1100e-01 -5.4627e-02 2.0  1.3644e+00  2.1006e+00  1.9852e+00  6.6942e-01  1.0678e+00  1.0024e+00  4.5   12.62;
5.0   12.62  7.9093e-01 -9.1466e-03 2.0  1.4484e+00  1.7639e+00  1.8515e+00  7.2182e-01  1.0770e+00  9.8893e-01  5.0   12.62;
5.5   12.62  7.8662e-01 -6.4267e-07 2.0  1.1562e+00  1.0929e+00  1.9671e+00  6.4646e-01  1.0865e+00  9.8772e-01  5.5   12.62;
6.0   12.62  8.4063e-01 -3.7478e-07 2.0  1.1179e+00  8.6813e-01  1.9900e+00  7.1682e-01  1.0880e+00  9.6728e-01  6.0   12.62;
6.5   12.62  8.2699e-01 -9.6153e-08 2.0  1.0649e+00  9.0260e-01  2.0259e+00  6.5749e-01  1.0958e+00  9.5398e-01  6.5   12.62;
7.0   12.62  8.0730e-01 -1.0776e-06 2.0  9.7457e-01  7.1721e-01  2.2248e+00  5.9962e-01  1.1051e+00  9.3355e-01  7.0   12.62;
7.5   12.62  8.2952e-01 -1.1781e-02 2.0  9.4264e-01  9.2716e-01  2.2199e+00  6.0866e-01  1.1050e+00  9.2204e-01  7.5   12.62;
8.0   12.62  1.1075e+00 -3.2347e-01 2.0  7.9284e-01  2.4170e+00  3.4246e+00  4.9618e-01  1.1125e+00  8.8831e-01  8.0   12.62;
4.0   20.00  7.0972e-01 -1.5611e-01 2.0  1.0568e+00  1.5270e+01  3.4374e+00  3.3059e-01  1.0844e+00  1.0550e+00  4.0   20.00;
4.5   20.00  7.5373e-01 -3.5155e-02 2.0  1.0502e+00  3.2069e+00  2.1437e+00  5.0568e-01  1.0894e+00  1.0373e+00  4.5   20.00;
5.0   20.00  7.4294e-01 -4.4353e-03 2.0  9.7401e-01  1.3124e+00  2.0837e+00  5.1430e-01  1.0829e+00  1.0347e+00  5.0   20.00;
5.5   20.00  7.2999e-01 -2.9787e-07 2.0  1.0646e+00  8.7860e-01  2.0446e+00  5.2091e-01  1.0870e+00  1.0000e+00  5.5   20.00;
6.0   20.00  7.9094e-01 -2.4672e-17 2.0  1.0599e+00  9.0053e-01  1.9662e+00  5.9849e-01  1.0887e+00  9.8811e-01  6.0   20.00;
6.5   20.00  8.0349e-01 -4.4357e-07 2.0  1.1059e+00  1.1712e+00  1.9580e+00  6.2395e-01  1.0934e+00  9.6524e-01  6.5   20.00;
7.0   20.00  8.2080e-01 -6.2367e-05 2.0  1.1345e+00  1.2266e+00  1.9089e+00  6.6433e-01  1.1049e+00  9.4121e-01  7.0   20.00;
7.5   20.00  8.2713e-01 -1.6679e-02 2.0  9.7085e-01  9.6843e-01  2.1135e+00  5.9319e-01  1.1024e+00  9.1320e-01  7.5   20.00;
8.0   20.00  1.1435e+00 -3.3384e-01 2.0  8.6005e-01  2.9950e+00  3.0842e+00  5.6602e-01  1.1145e+00  8.9025e-01  8.0   20.00;
4.0   31.70  7.5523e-01 -6.1519e-02 2.0  7.0295e-01  1.6974e+01  2.7714e+00  3.5788e-01  1.0968e+00  1.0716e+00  4.0   31.70;
4.5   31.70  7.2579e-01 -7.4006e-02 2.0  8.4409e-01  7.7628e+00  2.6833e+00  3.5510e-01  1.0895e+00  1.0478e+00  4.5   31.70;
5.0   31.70  7.3171e-01 -1.8856e-02 2.0  9.7754e-01  2.6668e+00  2.1737e+00  4.6838e-01  1.0968e+00  1.0442e+00  5.0   31.70;
5.5   31.70  7.4865e-01 -4.8344e-07 2.0  9.7714e-01  1.1250e+00  2.1068e+00  5.1255e-01  1.0960e+00  1.0332e+00  5.5   31.70;
6.0   31.70  7.6646e-01 -3.0853e-07 2.0  1.0203e+00  9.3153e-01  2.0629e+00  5.4787e-01  1.0977e+00  1.0095e+00  6.0   31.70;
6.5   31.70  7.7282e-01 -7.6203e-07 2.0  9.9068e-01  9.3050e-01  2.0019e+00  5.4762e-01  1.1004e+00  9.9098e-01  6.5   31.70;
7.0   31.70  8.0833e-01 -2.0127e-03 2.0  9.9126e-01  1.1425e+00  1.9250e+00  6.0071e-01  1.1018e+00  9.7147e-01  7.0   31.70;
7.5   31.70  8.2974e-01 -1.7129e-02 2.0  9.6726e-01  1.0363e+00  2.0572e+00  5.9916e-01  1.1023e+00  9.2972e-01  7.5   31.70;
8.0   31.70  1.1337e+00 -3.4815e-01 2.0  8.2312e-01  2.8601e+00  3.2358e+00  5.1106e-01  1.1149e+00  8.9638e-01  8.0   31.70;
4.0   50.24  7.0000e-01 -3.3844e-03 2.0  6.2430e-01  2.7092e+01  2.1023e+00  3.3039e-01  1.1043e+00  1.0745e+00  4.0   50.24;
4.5   50.24  7.0000e-01 -1.0683e-03 2.0  7.0260e-01  9.4854e+00  2.2301e+00  3.6400e-01  1.0999e+00  1.0722e+00  4.5   50.24;
5.0   50.24  7.0000e-01 -1.8942e-02 2.0  7.8355e-01  4.0256e+00  2.4505e+00  3.7077e-01  1.0970e+00  1.0617e+00  5.0   50.24;
5.5   50.24  7.0000e-01 -1.6300e-02 2.0  8.8500e-01  1.9832e+00  2.3716e+00  4.0325e-01  1.0971e+00  1.0463e+00  5.5   50.24;
6.0   50.24  7.1272e-01 -1.4726e-07 2.0  8.9286e-01  8.1783e-01  2.3500e+00  4.3997e-01  1.0965e+00  1.0360e+00  6.0   50.24;
6.5   50.24  7.6839e-01 -6.4575e-08 2.0  1.0072e+00  9.1313e-01  1.9774e+00  5.3680e-01  1.1030e+00  9.9851e-01  6.5   50.24;
7.0   50.24  7.9821e-01 -2.1127e-05 2.0  9.9853e-01  1.0824e+00  1.8584e+00  5.8686e-01  1.0951e+00  9.8511e-01  7.0   50.24;
7.5   50.24  8.3695e-01 -3.1385e-02 2.0  9.8773e-01  1.0788e+00  2.0739e+00  5.9519e-01  1.0989e+00  9.3891e-01  7.5   50.24;
8.0   50.24  1.2000e+00 -4.0437e-01 2.0  8.5839e-01  4.1919e+00  3.1777e+00  5.4781e-01  1.1109e+00  9.0655e-01  8.0   50.24;
4.0   79.62  7.0000e-01 -8.0890e-18 2.0  5.5429e-01  3.4017e+01  2.0664e+00  3.1187e-01  1.1005e+00  1.0754e+00  4.0   79.62;
4.5   79.62  7.0000e-01 -1.3090e-07 2.0  6.3896e-01  1.3043e+01  2.1523e+00  3.4652e-01  1.1058e+00  1.0788e+00  4.5   79.62;
5.0   79.62  7.0000e-01 -1.2290e-02 2.0  7.3580e-01  6.5228e+00  2.3315e+00  3.6224e-01  1.1035e+00  1.0648e+00  5.0   79.62;
5.5   79.62  7.0000e-01 -5.1176e-03 2.0  8.5381e-01  2.3219e+00  2.3280e+00  4.0281e-01  1.0951e+00  1.0580e+00  5.5   79.62;
6.0   79.62  7.4394e-01 -7.1538e-03 2.0  8.5853e-01  1.0128e+00  2.3976e+00  4.6298e-01  1.0970e+00  1.0547e+00  6.0   79.62;
6.5   79.62  7.3383e-01 -7.5166e-06 2.0  9.1533e-01  7.6082e-01  2.2108e+00  4.7191e-01  1.1022e+00  1.0224e+00  6.5   79.62;
7.0   79.62  7.6594e-01 -1.1135e-05 2.0  1.0242e+00  1.0944e+00  1.8323e+00  5.4115e-01  1.1060e+00  9.8341e-01  7.0   79.62;
7.5   79.62  8.3184e-01 -3.7984e-02 2.0  1.0252e+00  1.1924e+00  1.9481e+00  5.9080e-01  1.1065e+00  9.5858e-01  7.5   79.62;
8.0   79.62  1.2000e+00 -4.2855e-01 2.0  8.2436e-01  4.5732e+00  3.5456e+00  5.0411e-01  1.1105e+00  9.2490e-01  8.0   79.62;
4.0  126.19  7.0000e-01 -9.3134e-03 2.0  6.5770e-01  2.4111e+01  2.0921e+00  3.4515e-01  1.0948e+00  1.0737e+00  4.0  126.19;
4.5  126.19  7.0000e-01 -6.4122e-03 2.0  7.0065e-01  1.1473e+01  2.2173e+00  3.6651e-01  1.0964e+00  1.0705e+00  4.5  126.19;
5.0  126.19  7.0000e-01 -3.8103e-02 2.0  8.1428e-01  5.6824e+00  2.5118e+00  3.6708e-01  1.0982e+00  1.0603e+00  5.0  126.19;
5.5  126.19  7.0000e-01 -2.0781e-02 2.0  8.8580e-01  2.2925e+00  2.3286e+00  4.0074e-01  1.0892e+00  1.0465e+00  5.5  126.19;
6.0  126.19  7.1685e-01 -2.3828e-03 2.0  9.4160e-01  1.1819e+00  2.1530e+00  4.6235e-01  1.0937e+00  1.0376e+00  6.0  126.19;
6.5  126.19  7.5639e-01 -6.9748e-06 2.0  1.0266e+00  1.0918e+00  1.9603e+00  5.4218e-01  1.0949e+00  1.0160e+00  6.5  126.19;
7.0  126.19  7.8490e-01 -1.8834e-05 2.0  1.1253e+00  1.1569e+00  1.8651e+00  6.0514e-01  1.0910e+00  9.7531e-01  7.0  126.19;
7.5  126.19  8.2386e-01 -3.8922e-02 2.0  1.0819e+00  1.3516e+00  1.9944e+00  5.9384e-01  1.0999e+00  9.5643e-01  7.5  126.19;
8.0  126.19  1.2000e+00 -4.3500e-01 2.0  8.3252e-01  3.9690e+00  3.5290e+00  4.9975e-01  1.1068e+00  9.2419e-01  8.0  126.19;
4.0  200.00  7.0000e-01 -9.9722e-09 2.0  5.7369e-01  2.2867e+01  2.0999e+00  3.2613e-01  1.0893e+00  1.0764e+00  4.0  200.00;
4.5  200.00  7.0000e-01 -2.0290e-15 2.0  6.9245e-01  1.1532e+01  2.0741e+00  3.6652e-01  1.0896e+00  1.0693e+00  4.5  200.00;
5.0  200.00  7.0000e-01 -1.7951e-02 2.0  7.5549e-01  5.9753e+00  2.4248e+00  3.6718e-01  1.0879e+00  1.0610e+00  5.0  200.00;
5.5  200.00  7.0000e-01 -1.2546e-02 2.0  8.5451e-01  2.9617e+00  2.2767e+00  4.0641e-01  1.0935e+00  1.0547e+00  5.5  200.00;
6.0  200.00  7.1629e-01 -4.5233e-07 2.0  9.3155e-01  1.7244e+00  2.1745e+00  4.5411e-01  1.0864e+00  1.0374e+00  6.0  200.00;
6.5  200.00  7.4775e-01 -1.1187e-05 2.0  1.0262e+00  1.3210e+00  1.9536e+00  5.2329e-01  1.0925e+00  1.0271e+00  6.5  200.00;
7.0  200.00  7.8534e-01 -5.3183e-05 2.0  1.0803e+00  1.3301e+00  1.8174e+00  6.0116e-01  1.0887e+00  1.0000e+00  7.0  200.00;
7.5  200.00  8.3010e-01 -4.3565e-02 2.0  1.0546e+00  1.3811e+00  1.9608e+00  5.9829e-01  1.1017e+00  9.7409e-01  7.5  200.00;
8.0  200.00  1.2000e+00 -4.2648e-01 2.0  8.5386e-01  4.6946e+00  3.5826e+00  5.2896e-01  1.0980e+00  9.4662e-01  8.0  200.00;
4.0  317.00  7.0000e-01 -3.0941e-07 2.0  5.4038e-01  2.0798e+01  2.1620e+00  3.1839e-01  1.0999e+00  1.0873e+00  4.0  317.00;
4.5  317.00  7.0000e-01 -2.0482e-02 2.0  6.1477e-01  1.5841e+01  2.4554e+00  3.1946e-01  1.0955e+00  1.0786e+00  4.5  317.00;
5.0  317.00  7.0000e-01 -3.9962e-02 2.0  7.0732e-01  8.4084e+00  2.7415e+00  3.3211e-01  1.0985e+00  1.0637e+00  5.0  317.00;
5.5  317.00  7.0390e-01 -7.4679e-03 2.0  7.6534e-01  3.2582e+00  2.3231e+00  3.8841e-01  1.0927e+00  1.0643e+00  5.5  317.00;
6.0  317.00  7.3097e-01 -3.7794e-04 2.0  8.8502e-01  2.1430e+00  2.1211e+00  4.6048e-01  1.0888e+00  1.0423e+00  6.0  317.00;
6.5  317.00  7.7222e-01 -4.3701e-03 2.0  1.0979e+00  2.0230e+00  1.8593e+00  5.6905e-01  1.0903e+00  1.0215e+00  6.5  317.00;
7.0  317.00  7.9944e-01 -2.1160e-03 2.0  1.1200e+00  1.6581e+00  1.8141e+00  6.2103e-01  1.0856e+00  1.0083e+00  7.0  317.00;
7.5  317.00  7.9780e-01 -2.1297e-06 2.0  1.2097e+00  1.5811e+00  1.7106e+00  6.5362e-01  1.0932e+00  9.8560e-01  7.5  317.00;
8.0  317.00  1.2000e+00 -4.1799e-01 2.0  9.6603e-01  5.5529e+00  3.0556e+00  5.6936e-01  1.1012e+00  9.5754e-01  8.0  317.00;
4.0  502.40  7.0000e-01 -1.7111e-02 2.0  4.6999e-01  1.8968e+01  2.8200e+00  2.6648e-01  1.0937e+00  1.0850e+00  4.0  502.40;
4.5  502.40  7.0000e-01 -2.6396e-02 2.0  5.1949e-01  1.6368e+01  2.8759e+00  2.7832e-01  1.0958e+00  1.0749e+00  4.5  502.40;
5.0  502.40  7.0000e-01 -2.1671e-02 2.0  6.1000e-01  8.6654e+00  2.6980e+00  3.1148e-01  1.0901e+00  1.0729e+00  5.0  502.40;
5.5  502.40  7.0002e-01 -7.1630e-03 2.0  6.9537e-01  4.5437e+00  2.4581e+00  3.5831e-01  1.0952e+00  1.0655e+00  5.5  502.40;
6.0  502.40  7.6134e-01 -6.3443e-06 2.0  8.9341e-01  3.7976e+00  1.9228e+00  5.0275e-01  1.0876e+00  1.0463e+00  6.0  502.40;
6.5  502.40  7.9007e-01 -2.5364e-03 2.0  1.0237e+00  2.5624e+00  1.7955e+00  5.8082e-01  1.0837e+00  1.0309e+00  6.5  502.40;
7.0  502.40  7.8695e-01 -2.0933e-03 2.0  1.1867e+00  2.1554e+00  1.6970e+00  6.3350e-01  1.0993e+00  1.0215e+00  7.0  502.40;
7.5  502.40  8.1809e-01 -4.3629e-03 2.0  1.4327e+00  2.5223e+00  1.5747e+00  7.2501e-01  1.0885e+00  9.9483e-01  7.5  502.40;
8.0  502.40  1.2000e+00 -3.9261e-01 2.0  1.1094e+00  5.2480e+00  2.8220e+00  6.5055e-01  1.0920e+00  9.6559e-01  8.0  502.40;
4.0  796.20  7.0000e-01 -2.3099e-02 2.0  4.5581e-01  1.7534e+01  3.1836e+00  2.5169e-01  1.0938e+00  1.0935e+00  4.0  796.20;
4.5  796.20  7.0000e-01 -3.7902e-02 2.0  4.8070e-01  1.8783e+01  3.4708e+00  2.4695e-01  1.0957e+00  1.0901e+00  4.5  796.20;
5.0  796.20  7.0000e-01 -8.6256e-03 2.0  5.4088e-01  8.8700e+00  2.5344e+00  3.0798e-01  1.0999e+00  1.0796e+00  5.0  796.20;
5.5  796.20  7.2574e-01 -6.1859e-07 2.0  6.3759e-01  6.1694e+00  2.2528e+00  3.8172e-01  1.0942e+00  1.0656e+00  5.5  796.20;
6.0  796.20  7.9007e-01 -1.1663e-07 2.0  8.9429e-01  5.2998e+00  1.8310e+00  5.4462e-01  1.0835e+00  1.0541e+00  6.0  796.20;
6.5  796.20  8.0321e-01 -5.4738e-06 2.0  1.0887e+00  4.4263e+00  1.6296e+00  6.2320e-01  1.0833e+00  1.0565e+00  6.5  796.20;
7.0  796.20  8.3257e-01 -1.9010e-02 2.0  1.3544e+00  3.6403e+00  1.5790e+00  7.0373e-01  1.0842e+00  1.0292e+00  7.0  796.20;
7.5  796.20  8.2613e-01 -6.7595e-03 2.0  1.5248e+00  2.8118e+00  1.5264e+00  7.6541e-01  1.0887e+00  1.0129e+00  7.5  796.20;
8.0  796.20  1.1507e+00 -3.2529e-01 2.0  1.4633e+00  5.9123e+00  2.1639e+00  7.6925e-01  1.0861e+00  9.8940e-01  8.0  796.20;
4.0 1262.00  7.0380e-01 -2.3859e-02 2.0  4.6750e-01  1.6009e+01  3.2471e+00  2.5732e-01  1.0983e+00  1.0854e+00  4.0 1262.00;
4.5 1262.00  7.2948e-01 -1.4944e-02 2.0  4.8377e-01  1.2697e+01  2.6412e+00  2.9720e-01  1.0910e+00  1.0799e+00  4.5 1262.00;
5.0 1262.00  7.1686e-01 -8.8205e-09 2.0  5.0825e-01  8.9544e+00  2.5879e+00  3.1741e-01  1.0940e+00  1.0730e+00  5.0 1262.00;
5.5 1262.00  7.5789e-01 -5.9107e-07 2.0  6.3038e-01  8.3660e+00  2.0527e+00  4.2616e-01  1.0955e+00  1.0773e+00  5.5 1262.00;
6.0 1262.00  7.8795e-01 -1.1888e-06 2.0  9.4374e-01  7.9155e+00  1.6171e+00  5.6637e-01  1.0922e+00  1.0692e+00  6.0 1262.00;
6.5 1262.00  8.3205e-01 -7.7356e-04 2.0  1.3982e+00  7.1069e+00  1.4151e+00  7.1961e-01  1.0743e+00  1.0470e+00  6.5 1262.00;
7.0 1262.00  8.2346e-01 -5.7461e-03 2.0  1.6840e+00  5.8126e+00  1.3549e+00  7.5587e-01  1.0847e+00  1.0338e+00  7.0 1262.00;
7.5 1262.00  8.4754e-01 -2.1169e-02 2.0  1.9081e+00  4.5097e+00  1.3691e+00  8.1498e-01  1.0837e+00  1.0215e+00  7.5 1262.00;
8.0 1262.00  8.5187e-01 -1.5704e-02 2.0  2.7561e+00  4.7263e+00  1.2542e+00  9.5985e-01  1.0839e+00  1.0185e+00  8.0 1262.00 ]


text:
abstract type WordVectorsEmbedder{T,S} <: AbstractEmbedder{T,S} end

"""
Constant that represents embeddings libraries used in text embedding.
"""
const EmbeddingsLibrary{T,S} = Union{
    ConceptNet{<:Languages.Language,S,T},
    Word2Vec.WordVectors{S,T,<:Integer},
    Glowe.WordVectors{S,T,<:Integer},
    EmbeddingsAnalysis.CompressedWordVectors{<:QuantizedArrays.AbstractQuantization,
                                             <:Unsigned,<:Distances.PreMetric,T,S,<:Integer}
}

"""
    document2vec(embedder, document)

Word-embeddings approach to document embedding. It embeds documents
using word embeddings libraries and some algorithm for combining
these (depending on the type of `embedder`).

# Arguments
  * `embedder::WordVectorsEmbedder` is the embedder
  * `document::Vector{String}` the document to be embedded,
     where each vector element corresponds to a sentence
"""
function __document2vec(embedder::WordVectorsEmbedder{T,S},
                        document::Vector{String};
                        kwargs...  # for the unused arguments
                       ) where {T,S}
    # Initializations
    n = length(document)
    m = size(embedder.embeddings)[1]  # number of vector components
    embedded_words = Vector{Vector{String}}(undef, n)
    doc_word_embeddings = Vector{Matrix{T}}(undef, n)

    # Get word embeddings for each sentence
    @inbounds for i in 1:n
        words = tokenize(document[i], method=:stringanalysis)
        _embs, _mtoks = word_embeddings(embedder.embeddings,
                                        words;
                                        keep_size=false,
                                        max_compound_word_length=1,
                                        wildcard_matching=true,
                                        print_matched_words=false)
        doc_word_embeddings[i] = _embs
        embedded_words[i] = words[setdiff(1:length(words), _mtoks)]
    end

    # Remove empty embeddings
    filter!(!isempty, doc_word_embeddings)
    filter!(!isempty, embedded_words)

    # Create sentence embeddings
    sntembs = sentences2vec(embedder,
                            doc_word_embeddings;
                            embedded_words=embedded_words)
    # If no words were embedded, sntembs is a `D`×0 Matrix
    # which, when squashed, becomes a `D`-element zero Vector
    # (`D` is the dimensionality of the document embedding)
    return squash(sntembs)
end

function document2vec(embedder::WordVectorsEmbedder{T,S},
                      document::Vector{String};
                      kwargs...  # for the unused arguments
                     ) where {T,S}
    # An embedded document is a zero vector only if it was not properly embedded
    embedded_document = __document2vec(embedder, document; kwargs...)
    is_embedded = !iszero(embedded_document)
    # Check for OOV (out-of-vocabulary) policy
    if embedder.config.oov_policy == :large_vector && !is_embedded
        embedded_document .= T(DEFAULT_OOV_VAL)
    end
    return embedded_document, is_embedded
end

"""
    word_embeddings(word_vectors, document_tokens [;kwargs])

Returns a matrix corresponding to the word embeddings of `document_tokens`
as well as the indices of missing i.e. not-embedded tokens.

# Arguments
  * `word_vectors::EmbeddingsLibrary` wordvectors object; can be
     a `Word2Vec.WordVectors`, `Glowe.WordVectors` or `ConceptnetNumberbatch.ConceptNet`
  * `document_tokens::Vector{String}` the words to be embedded,
     where each vector element corresponds to a word

# Keyword arguments
  * `keep_size::Bool` a `false` value discards vectors for words not found
     while a `true` value (default) places a zero vector in the embeddings
     matrix
  * `print_matched_words::Bool` if `true`, the words that were and that were
     not embedded are printed (default `false`)
  * `kwargs...` the rest of the keyword arguments are `ConceptNet` specific
     and can be found by inspecting the help of `ConceptnetNumberbatch.embed_document`
"""
function word_embeddings(word_vectors::Union{Word2Vec.WordVectors{S1,T,H},
                                             Glowe.WordVectors{S1,T,H},
                                             EmbeddingsAnalysis.CompressedWordVectors{Q,U,D,T,S1,H}},
                         document_tokens::Vector{S2};
                         keep_size::Bool=true,
                         print_matched_words::Bool=false,
                         kwargs...
                        ) where {S1<:AbstractString, T<:Real,
                                 H<:Integer, S2<:AbstractString,
                                 Q<:QuantizedArrays.AbstractQuantization,
                                 U<:Unsigned, D<:Distances.PreMetric}
    # Initializations
    n = size(word_vectors)[1]
    p = length(document_tokens)
    found_positions = Int[] # column positions in word embedding matrix of found tokens
    found_tokens = Int[]    # positions in the token vector on found tokens

    # Search for matching words (exact match)
    for (i, token) in enumerate(document_tokens)
        if haskey(word_vectors.vocab_hash, token)
            push!(found_positions, word_vectors.vocab_hash[token])
            push!(found_tokens, i)
        end
    end
    missing_tokens = setdiff(1:p, found_tokens)  # positions of tokens that were not found

    # Construct document structure
    m = ifelse(keep_size, p, length(found_tokens))
    embedding = zeros(T, n, m)
    for (pos, tok, j) in zip(found_positions, found_tokens, 1:m)
        embedding[:, ifelse(keep_size, tok, j)] = word_vectors.vectors[:, pos]
    end

    # Return document embedding and missing tokens
    if print_matched_words
        println("Embedded words: $(tokens[found_tokens])")
        println("Mismatched words: $(tokens[missing_tokens])")
    end
    return embedding, missing_tokens
end

# Replicate ConceptnetNumberbatch embed_document function
function word_embeddings(conceptnet::ConceptNet{L,K,E},
                         document_tokens::Vector{S};
                         language=Languages.English(),
                         keep_size::Bool=true,
                         compound_word_separator::String="_",
                         max_compound_word_length::Int=1,
                         wildcard_matching::Bool=false,
                         print_matched_words::Bool=false
                        ) where {L<:Language, K, E<:Real, S<:AbstractString}
    return ConceptnetNumberbatch.embed_document(conceptnet, document_tokens,
                language=language, keep_size=keep_size,
                compound_word_separator=compound_word_separator,
                max_compound_word_length=max_compound_word_length,
                wildcard_matching=wildcard_matching,
                print_matched_words=print_matched_words)
end

"""
    squash(m)

Function that creates a single mean vector from a matrix `m`
and performs some normalization operations as well.
"""
squash(m::Matrix{T}) where T<:AbstractFloat = begin
    # Calculate document embedding
    ##############################
    v = vec(sum(m, dims=2))
    return v./(norm(v,2)+eps(T))
end

"""
    squash(vv, m)

Function that creates a single mean vector from a vector of vectors `vv`
where each vector has a length `m` and performs some normalization
operations as well.
"""
squash(vv::Vector{Vector{T}}, m::Int) where T<:AbstractFloat = begin
    v = zeros(T,m)
    for w in vv
        v+= w
    end
    return v./(norm(v,2)+eps(T))
end

ext: jl


text:
using PastaQ
using ITensors
using Random

Random.seed!(1234)

N = 4     # Number of qubits
depth = 4 # Depth of the quantum circuit

# Generate random quantum circuit built out of
# layers of single-qubit random rotations + `CX`
# gates, alternating between even and of odd layers.
println("Random circuit of depth $depth on $N qubits:")
circuit = randomcircuit(N, depth; twoqubitgates="CX", onequbitgates="Rn")
display(circuit)
println()

# 1. Unitary quantum circuit
# Returns the MPS at the output of the quantum circuit:
# `|ψ⟩ = Û|0,0,…,0⟩`
# where `Û` is the unitary circuit.
println("Applying random circuit to compute |ψ⟩ = U|0,0,…,0⟩...")
ψ = runcircuit(circuit)
@show maxlinkdim(ψ)
println()

# A representation of the unitary operation as a MPO
# is obtained using the flag `process=true`:
println("Approximating random circuit as an MPO U...")
U = runcircuit(circuit; process=true)
@show maxlinkdim(U)
println()

# 2. Quantum circuit with noise
# Apply a noise model `noise` (specified by appropriate
# parameters) to each quantum gate in `gates`.
# Returns a mixed density operator as MPO:
# `ρ = ε(|0,0,…⟩⟨0,0,…|)`
# where `ε` is the quantum channel.
# Here, the noise is a single-qubit amplitude damping
# channel with decay rate `γ=0.01`..
println(
  "Running the circuit with amplitude damping to compute the state ρ = ε(|0,0,…⟩⟨0,0,…|)..."
)
ρ = runcircuit(circuit; noise=("amplitude_damping", (γ=0.01,)))
@show maxlinkdim(ρ)
println()

# A representation of the quantum channel as a MPO
# is obtained using the flag `process=true`, which
# returns the Choi matrix `Λ` of the channel:`:
println(
  "Running the circuit with amplitude damping to compute the Choi matrix Λ of the quantum channel...",
)
Λ = runcircuit(circuit; process=true, noise=("amplitude_damping", (γ=0.01,)))
@show maxlinkdim(Λ)


text:
# Author: Lukasz Bratos

module Methods

export mbisekcji, mstycznych, msiecznych

"""
Funkcja rozwiazuje rownanie postaci f(x) = 0 metoda bisekcji.
function mbisekcji(f, a::Float64, b::Float64, delta::Float64, epsilon::Float64)
    Dane:
    f – funkcja f(x) zadana jako anonimowa funkcja,
    a, b – końce przedziału początkowego,
    delta, epsilon - dokładności obliczeń,
    Wyniki:
    (r, v, it, err) – czwórka, gdzie
    r – przybliżenie pierwiastka równania f(x) = 0,
    v – wartość f(r),
    it – liczba wykonanych iteracji,
    err – sygnalizacja błędu
        0 - brak błędu
        1 - funkcja nie zmienia znaku w przedziale [a,b]
"""
function mbisekcji(f, a::Float64, b::Float64, delta::Float64, epsilon::Float64)
    fa = f(a)
    fb = f(b)
    e = b - a
    err = 0
    it = 0
    r = Float64(0.0)
    v = Float64(0.0)

    if sign(fa) == sign(fb)
        err = 1
        return (r, v, it, err)
    end

    while e > epsilon
        it += 1

        e = e / 2
        r = a + e
        v = f(r)

        if abs(e) < delta || abs(v) < epsilon
            return (r, v, it, err)
        end

        if sign(v) != sign(fa)
            b = r
            fb = v
        else
            a = r
            fa = v
        end
    end
end

"""
Funkcja rozwiazuje rownanie postaci f(x) = 0 metoda Newtona (stycznych).
function mstycznych(f,pf,x0::Float64, delta::Float64, epsilon::Float64, maxit::Int)
    Dane:
    f, pf– funkcją f(x) oraz pochodną f′(x) zadane jako anonimowe funkcje,
    x0 – przybliżenie początkowe,
    delta, epsilon – dokładności obliczeń,
    maxit – maksymalna dopuszczalna liczba iteracji,
    Wyniki:
    (r, v, it, err) – czwórka, gdzie
    r – przybliżenie pierwiastka równania f(x) = 0,
    v – wartość f(r),
    it – liczba wykonanych iteracji,
    err – sygnalizacja błędu
        0 - metoda zbieżna
        1 - nie osiągnięto wymaganej dokładności w maxit iteracji,
        2 - pochodna bliska zeru
"""
function mstycznych(f, pf, x0::Float64, delta::Float64, epsilon::Float64, maxit::Int)
    v = f(x0)
    err = 0

    if abs(v) < epsilon
        return (x0, v, 0, err)
    end

    for it in 1:maxit
        if abs(pf(x0)) < epsilon
            # pochodna bliska zeru
            err = 2
            return (x0, v, it, err)
        end

        x1 = x0 - v / pf(x0)
        v = f(x1)

        if abs(x1 - x0) < delta || abs(v) < epsilon
            return (x1, v, it, err)
        end

        x0 = x1
    end

    # nie osiagnieto wymaganej dokladnosci w maxit iteracji
    err = 1
    return (x0, v, maxit, err)
end

"""
Funkcja rozwiazuje rownanie postaci f(x) = 0 metoda siecznych.
function msiecznych(f, x0::Float64, x1::Float64, delta::Float64, epsilon::Float64,maxit::Int)
    Dane:
    f – funkcja f(x)zadana jako anonimowa funkcja,
    x0, x1 – przybliżenia początkowe,
    delta, epsilon – dokładności obliczeń,
    maxit – maksymalna dopuszczalna liczba iteracji,
    Wyniki:
    (r, v, it, err)– czwórka, gdzie
    r – przybliżenie pierwiastka równania f(x) = 0,
    v – wartość f(r),
    it – liczba wykonanych iteracji,
    err – sygnalizacja błędu
    0 - metoda zbieżna
    1 - nie osiągnięto wymaganej dokładności w maxit iteracji
"""
function msiecznych(f, x0::Float64, x1::Float64, delta::Float64, epsilon::Float64, maxit::Int)
    fx0 = f(x0)
    fx1 = f(x1)
    err = 0

    for it in 1:maxit
        if abs(fx0) > abs(fx1)
            x0, x1 = x1, x0
            fx0, fx1 = fx1, fx0
        end

        s = (x1 - x0) / (fx1 - fx0)
        x1 = x0
        fx1 = fx0
        x0 = x0 - fx0 * s
        fx0 = f(x0)

        if abs(x1 - x0) < delta || abs(fx0) < epsilon
            return (x0, fx0, it, err)
        end
    end

    # nie osiagnieto wymaganej dokladnosci w maxit iteracji
    err = 1
    return (x0, fx0, maxit, err)
end

end


text:
"""
In Block
"""

export InBlock

mutable struct InBlock <: AbstractBlock
    inport::Vector{InPort}
    outport::Vector{OutPort}

    function InBlock()
        @createblock new(Vector{InPort}(), Vector{OutPort}()) 1 1
    end
end

"""
IO
"""

function Base.show(io::IO, blk::InBlock)
    println(io, "InBlock()")
end

"""
to expr
"""

function _toexpr(blk::InBlock)
    _toexpr(blk.inport[1])
end


text:
"""
    CoordinateBoundaryConditions(left, right)

A set of two `BoundaryCondition`s to be applied along a coordinate x, y, or z.

The `left` boundary condition is applied on the negative or lower side of the coordinate
while the `right` boundary condition is applied on the positive or higher side.
"""
mutable struct CoordinateBoundaryConditions{L, R}
     left :: L
    right :: R
end


text:
using Distributions

function residuals_subperiods(fm::FactorModel, break_period, variable_index)  # gets residuals from OLS regression of x on factors for both subsamples for variable_index (i.e. column of the x matrix)
    #fm1 = FactorModel(fm.x[1:break_period, :], "ICp1")  # select number of factors according to ICp1
    #fm2 = FactorModel(fm.x[break_period+1:end, :], "ICp1")  # select number of factors according to ICp1
    x1, y1 = fm.factors[1:break_period, 1:fm.number_of_factors], fm.x[1:break_period, variable_index]
    x2, y2 = fm.factors[break_period+1:end, 1:fm.number_of_factors], fm.x[break_period+1:end, variable_index]
    #factors1, loadings1, number_of_factors1 = calculate_factors(fm.x[1:break_period, :])
    #factors2, loadings2, number_of_factors2 = calculate_factors(fm.x[break_period+1:end, :])
    #return( (fm.x[1:break_period, variable_index] - (factors1[:, 1:fm.number_of_factors]*loadings1[:, fm.number_of_factors]')[:, variable_index]), (fm.x[break_period+1:end, variable_index] - (factors2[:, 1:fm.number_of_factors]*loadings2[:, fm.number_of_factors]')[:, variable_index]) )  # TODO: number of factors is likely to be higher if a structural break is present. --> number of factors should be reestimated for subperiods
    return( ((eye(break_period) - x1*inv(x1'x1)x1')y1), ((eye(length(y2)) - x2*inv(x2'x2)x2')y2) )   # use OLS for subperiods as in Breitung Eickmeier 2011 or Factor Model as in the 2009 version of the same paper?
end

# static factor Chow-tests Breitung, Eickmeier (2011)

# LR test
function LR_test(fm::FactorModel, break_period::Int64, variable_index::Int64)  # TODO: in 2009 version of their paper Breitung and Eickmeier note that we can also use residuals from OLS regressions because they are asymptotically the same as PCA
    T = size(fm.x, 1)
    likelihood_ratio = T * (log(sum(fm.factor_residuals[:, variable_index].^2)) - log(apply(+, map(res->sum(res.^2), residuals_subperiods(fm, break_period, variable_index)))))
    return likelihood_ratio
end

function Wald_test(fm::FactorModel, break_period::Int64, variable_index::Int64)
    design_matrix = hcat(fm.factors[:, 1:fm.number_of_factors], fm.factors[:, 1:fm.number_of_factors].*[zeros(break_period), ones(size(fm.x, 1)-break_period)])
    ols_estimates = inv(design_matrix'design_matrix)design_matrix'fm.x[:, variable_index]
    residuals = fm.x[:, variable_index] - design_matrix * ols_estimates
    # calculate variance-covariance matrix according to White(1980) TODO: maybe this should be maximum likelihood estimation of var-cov matrix (i.e. negative inverse of information matrix)
    coefficient_covariance = inv(design_matrix'design_matrix)*(design_matrix'*diagm(residuals.^2)*design_matrix)*inv(design_matrix'design_matrix)
    wald_stat = (ols_estimates[fm.number_of_factors+1:end]'inv(coefficient_covariance[fm.number_of_factors+1:end, fm.number_of_factors+1:end])ols_estimates[fm.number_of_factors+1:end])[1]
    return wald_stat
end

function LM_test(fm::FactorModel, break_period::Int64, variable_index::Int64)
    design_matrix = hcat(fm.factors[:, 1:fm.number_of_factors], fm.factors[:, 1:fm.number_of_factors].*[zeros(break_period), ones(size(fm.x, 1)-break_period)])
    ols_estimates = inv(design_matrix'design_matrix)design_matrix'fm.factor_residuals[:, variable_index]
    residuals = fm.factor_residuals[:, variable_index] - design_matrix * ols_estimates
    R_squared = 1 - sum(residuals.^2) / sum(fm.factor_residuals[:, variable_index].^2)
    lagrange_multiplier = size(fm.x, 1) * R_squared
end

# GLS version of LM-statistic: can be used for approximate factor models where the error term results from an AR(p) model (see Breitung, Eickmeier 2011)
function LM_test_gls(fm::FactorModel, break_period::Int64, variable_index::Int64, max_lags::Int64=10)
    # get gls transformed series for period before and after the break
    lag_polynomial = estimate_lag_polynomial(fm, variable_index, max_lags)[1]
    white_variance1, white_variance2 = estimate_white_variance_subperiods(fm, variable_index, length(lag_polynomial), break_period)
    gls_x = vcat(  # transform both subperiods separately
        gls_transform(fm.x[1:break_period, variable_index], lag_polynomial, white_variance1),
        gls_transform(fm.x[break_period+1:end, variable_index], lag_polynomial, white_variance2)
    )
    gls_factors = vcat(  # same goes for factors
        apply(hcat, [gls_transform(fm.factors[1:break_period, r], lag_polynomial, white_variance1) for r in 1:fm.number_of_factors]),
        apply(hcat, [gls_transform(fm.factors[break_period+1:end, r], lag_polynomial, white_variance2) for r in 1:fm.number_of_factors])
    )
    gls_factors_star = vcat(zeros(break_period, fm.number_of_factors), gls_factors[break_period+1:end, :])
    design_matrix = hcat(gls_factors, gls_factors_star)
    ols_estimates = inv(design_matrix'design_matrix)design_matrix'gls_x
    residuals = gls_x - design_matrix * ols_estimates
    R_squared = 1 - sum(residuals.^2) / sum(gls_x.^2)
    lagrange_multiplier = length(gls_x) * R_squared
    # TODO: this gives way too high a rejection rate
end

function gls_transform(vec::Array{Float64, 1}, polynomial::Array{Float64, 1}, variance_estimate::Float64=1)  # gls_transform a vector given a lag polynomial and an estiamtion of the variance (e.g. \hat \rho y_{it} on page 75 of Breitung, Eickmeier (2011))
    lagged_vecs = lagged_matrix(vec, [0:length(polynomial)])
    (lagged_vecs[:, 1] - sum([lagged_vecs[:, i+1]*polynomial[i] for i in 1:length(polynomial)])) / variance_estimate
end

function estimate_lag_polynomial(fm::FactorModel, variable_index::Int64, max_lags::Int64=10)  # used for GLS transformation of the model
    # Breitung, Eickmeier 2011 page 75
    # we could also use a HAC estimator (see page 75 of Breitung, Eickmeier 2011 but this can be less efficient in small samples)
    lag_num = indmin([bic_criterion(regress_on_lags(fm.factor_residuals[:, variable_index], p)[2], p) for p in 1:max_lags])
    polynomial, residuals = regress_on_lags(fm.factor_residuals[:, variable_index], lag_num)

end

function estimate_white_variance_subperiods(fm::FactorModel, variable_index::Int64, lag_number::Int64, break_period::Int64)
    coefs, residuals = regress_on_lags(fm.factor_residuals[:, variable_index], lag_number)
    x1, x2 = reshape(fm.factor_residuals[1+lag_number:break_period, variable_index], (break_period-lag_number, 1)), reshape(fm.factor_residuals[break_period+1:end, variable_index], (size(fm.x, 1)-break_period, 1))
    resid1, resid2 = residuals[1:break_period-lag_number], residuals[break_period-lag_number+1:end]
    (inv(x1'x1)*(x1'*diagm(resid1.^2)*x1)*inv(x1'x1))[1], (inv(x2'x2)*(x2'*diagm(resid2.^2)*x2)*inv(x2'x2))[1]
end

function regress_on_lags(vec::Array{Float64, 1}, p::Int64)  # returns coefficients and residuals from a regression of vec on p of its lags
    lags = lagged_matrix(vec, [0:p])
    y, x = lags[:, 1], lags[:, 2:end]
    coefs = inv(x'x)*x'y
    resid = y - x * coefs
    coefs, resid
end

function bic_criterion(residuals, k)
    T = length(residuals)
    log(var(residuals)) + k*log(T)/T # from Bai Ng (2008)
end

function quandt_andrews(fm::FactorModel, test_statistic::Function, variable_index::Int64; obs_margin_ratio=0.15)
    first_obs, last_obs = int(ceil(obs_margin_ratio*size(fm.x, 1)))+1, int(ceil((1-obs_margin_ratio)*size(fm.x, 1)))-1  # first and last obs to be tested such that there are some obs at the margins
    test_stats = [test_statistic(fm, t, variable_index) for t in first_obs:last_obs]  # note that LM_test_gls takes default value for max_lags
    stats_order = sortperm(test_stats, rev=true)
    return (first_obs+stats_order-1, test_stats[stats_order])
end

function quandt_andrews_critical_value(r::Int64, p_value::Float64=0.05)  # returns critical value of sup LM statistci for a given number of r = factors / degrees of freedom  taken from Andrews (2003)
    quandt_andrews_critical_values = [[0.01 => 12.16, 0.05 => 8.68, 0.1 => 7.12],
        [0.01 => 15.56, 0.05 => 11.72, 0.1 => 10.0],
        [0.01 => 18.07, 0.05 => 14.14, 0.1 => 12.28],
        [0.01 => 20.47, 0.05 => 16.36, 0.1 => 14.34],
        [0.01 => 22.66, 0.05 => 18.32, 0.1 => 16.30],
        [0.01 => 24.74, 0.05 => 20.24, 0.1 => 18.11],
        [0.01 => 26.72, 0.05 => 22.06, 0.1 => 19.87]
    ]  # TODO: add rest from table or write a simulation for critical values
    if r > length(quandt_andrews_critical_values)
        warn(string("r>", length(quandt_andrews_critical_values), " not implemented yet. Pull requests are welcome."))
    end
    return quandt_andrews_critical_values[r][p_value]
end

function chow_test(fm::FactorModel, test_statistic::Function, break_period::Int64, variable_index::Int64, significance_level::Float64=0.05)  # test one of the chow tests
    critical_value = quantile(Distributions.Chisq(fm.number_of_factors), 1-significance_level)
    test_statistic(fm, break_period, variable_index) > critical_value
end

function bootstrap_test(fm::FactorModel, test_statistic::Function, break_period::Int64, variable_index::Int64, bootstrap::Function=residual_bootstrap, significance_level::Float64=0.05; B=1000)  # same as chow_test but bootstraps rejections
    stat = test_statistic(fm, break_period, variable_index)
    bootstraps = bootstrap(fm, B, boot_fm->test_statistic(boot_fm, break_period, variable_index))
    p_value = mean(bootstraps .> stat)
    p_value < significance_level
end


text:
"""
    safe_div(x, y)

Safely divide `x` by `y`. If `y` is zero, return `x` directly.
"""
safe_div(x, y) = ifelse(iszero(y), x, x/y)

"""
    maximum_dims(dims)

Return the maximum value for each dimension. An array of dimensions `dims` is accepted.
The maximum of each dimension in the element is computed.
"""
maximum_dims(dims::AbstractArray{<:Integer}) = (maximum(dims), )
maximum_dims(dims::AbstractArray{NTuple{N, T}}) where {N,T} = ntuple(i -> maximum(x->x[i], dims), N)
maximum_dims(dims::AbstractArray{CartesianIndex{N}}) where {N} = ntuple(i -> maximum(x->x[i], dims), N)

function reverse_indices(idx::Array{T}) where T
    rev = Dict{T,Vector{CartesianIndex}}()
    for (ind, val) = pairs(idx)
        rev[val] = get(rev, val, CartesianIndex[])
        push!(rev[val], ind)
    end
    rev
end

function count_indices(idx::AbstractArray)
    counts = zero.(idx)
    for i in unique(idx)
        counts += sum(idx .== i) * (idx .== i)
    end
    return counts
end

function divide_by_counts!(xs, idx::AbstractArray, dims)
    colons = Base.ntuple(_->Colon(), dims)
    counts = count_indices(idx)
    for i in CartesianIndices(counts)
        view(xs, colons..., i) ./= counts[i]
    end
    return xs
end


text:
function helper_fn(z, attempt)
	method, N, M, d = z
	proc_id = myid();

	helper_start_time = time()

	time_limit_overall = 6 * 60 * 60 # 6 hours
	time_limit_sub = 30 * 60 # 30 minutes
	depth_limit = d

	grand_choice_probs = readcsv(string("../tcm_data/synthetic_data_ECG_models/MNL_ECG2_Neq", N, "_grand_model_predict.csv"))
	grand_assortments = convert(Array{Int64}, readcsv(string("../tcm_data/synthetic_data_ECG_models/MNL_ECG2_Neq", N, "_grand_assortments.csv")))

	assortments = copy(grand_assortments)
	choice_probs = copy(grand_choice_probs)

	test_assortments = copy(grand_assortments)
	test_choice_probs = copy(grand_choice_probs)

	println("Before if")

	if (M < 2^(N-1))
		ass_ord_inds = convert(Array{Int64}, readcsv(string("../tcm_data/synthetic_data_ECG_models/Neq", N, "_randperm.csv")) )
		@show typeof(ass_ord_inds)
		assortments = grand_assortments[ ass_ord_inds[1:M], :]
		choice_probs = grand_choice_probs[ ass_ord_inds[1:M], :]

		test_assortments = grand_assortments[ ass_ord_inds[(M+1):end], :]
		test_choice_probs = grand_choice_probs[ ass_ord_inds[(M+1):end], :]
	end

	println("After if")
	num_train_assortments = size(assortments,1)
	num_test_assortments = size(test_assortments, 1)

	# FOREST, L1, COLD-START, SUBPROBLEM VIA TOP DOWN INDUCTION
	if (method == "HCG")
		initial_left = Array{Int64,1}[]
		initial_right = Array{Int64,1}[]
		initial_product = Array{Int64,1}[]
		initial_isLeaf = Array{Bool,1}[]

		println("Starting forest L1 HCG")
		lambda, forest_left, forest_right, forest_product, forest_isLeaf, primal_objval, elapsed_time = tcm_estimateForest_L1(N, assortments, choice_probs, time_limit_overall, depth_limit, initial_left, initial_right, initial_product, initial_isLeaf)
		L1_error = primal_objval / size(assortments,1)

		K = length(lambda)
		forest_predict, forest_predict_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, assortments)
		forest_train_error = sum( abs.( forest_predict - choice_probs  )  ) / size(assortments,1)

		forest_all_predict, forest_all_predict_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, grand_assortments)
		forest_all_error = sum( abs.( forest_all_predict - grand_choice_probs  )  ) / size(grand_assortments,1)

		forest_test_predict, forest_test_predict_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, test_assortments)
		forest_test_error = sum( abs.( forest_test_predict - test_choice_probs  )  ) / size(test_assortments,1)

		row_string = string(N,  ",",
							num_train_assortments, ",",
							d, ",",
							"forest-cold-HCG", ",",
							proc_id, ",",
							attempt, ",",
							num_test_assortments, ",",
							primal_objval, ",",
							L1_error, ",",
							forest_train_error, ",",
							forest_test_error, ",",
							forest_all_error, ",",
							K, ",",
							elapsed_time)
		row_string = string(row_string, "\n")

		@spawnat 1 begin
			print(forest_outcsvhandle, row_string)
			flush(forest_outcsvhandle);
		end
	end

	# FOREST, L1, COLD-START, SUBPROBLEM VIA MIO
	if (method == "ECG")
		initial_left = Array{Int64,1}[]
		initial_right = Array{Int64,1}[]
		initial_product = Array{Int64,1}[]
		initial_isLeaf = Array{Bool,1}[]

		println("Starting forest L1 ECG")
		lambda_MIO, forest_left_MIO, forest_right_MIO, forest_product_MIO, forest_isLeaf_MIO, primal_objval_MIO, elapsed_time_MIO = tcm_estimateForest_L1_MIO(N, assortments, choice_probs, time_limit_overall, time_limit_sub, depth_limit, initial_left, initial_right, initial_product, initial_isLeaf)

		L1_error_MIO = primal_objval_MIO / size(assortments,1)

		K_MIO = length(lambda_MIO)
		forest_MIO_predict, forest_MIO_predict_by_class = tcm_predictForest(N, K_MIO, lambda_MIO, forest_left_MIO, forest_right_MIO, forest_product_MIO, forest_isLeaf_MIO, assortments)
		forest_MIO_train_error = sum( abs.( forest_MIO_predict - choice_probs  )  ) / size(assortments,1)

		forest_MIO_all_predict, forest_MIO_all_predict_by_class = tcm_predictForest(N, K_MIO, lambda_MIO, forest_left_MIO, forest_right_MIO, forest_product_MIO, forest_isLeaf_MIO, grand_assortments)
		forest_MIO_all_error = sum( abs.( forest_MIO_all_predict - grand_choice_probs  )  ) / size(grand_assortments,1)

		forest_MIO_test_predict, forest_MIO_test_predict_by_class = tcm_predictForest(N, K_MIO, lambda_MIO, forest_left_MIO, forest_right_MIO, forest_product_MIO, forest_isLeaf_MIO, test_assortments)
		forest_MIO_test_error = sum( abs.( forest_MIO_test_predict - test_choice_probs  )  ) / size(test_assortments,1)

		row_string = string(N,  ",",
							num_train_assortments, ",",
							d, ",",
							"forest-cold-ECG", ",",
							proc_id, ",",
							attempt, ",",
							num_test_assortments, ",",
							primal_objval_MIO, ",",
							L1_error_MIO, ",",
							forest_MIO_train_error, ",",
							forest_MIO_test_error, ",",
							forest_MIO_all_error, ",",
							K_MIO, ",",
							elapsed_time_MIO)
		row_string = string(row_string, "\n")

		@spawnat 1 begin
			print(forest_outcsvhandle, row_string)
			flush(forest_outcsvhandle);
		end
	end

	if (method == "RTS")
		initial_left = Array{Int64,1}[]
		initial_right = Array{Int64,1}[]
		initial_product = Array{Int64,1}[]
		initial_isLeaf = Array{Bool,1}[]

		K_vec = [100, 200, 500, 1000, 2000, 5000, 10000, 20000, 100000 ]

		for K in K_vec

			println("Starting forest L1 RTS")
			lambda_RTS, forest_left_RTS, forest_right_RTS, forest_product_RTS, forest_isLeaf_RTS, primal_objval_RTS, elapsed_time_RTS = tcm_estimateForest_L1_RTS(N, assortments, choice_probs, depth_limit, K)

			L1_error_RTS = primal_objval_RTS / size(assortments,1)

			forest_RTS_predict, forest_RTS_predict_by_class = tcm_predictForest(N, K, lambda_RTS, forest_left_RTS, forest_right_RTS, forest_product_RTS, forest_isLeaf_RTS, assortments)
			forest_RTS_train_error = sum( abs.( forest_RTS_predict - choice_probs  )  ) / size(assortments,1)

			forest_RTS_all_predict, forest_RTS_all_predict_by_class = tcm_predictForest(N, K, lambda_RTS, forest_left_RTS, forest_right_RTS, forest_product_RTS, forest_isLeaf_RTS, grand_assortments)
			forest_RTS_all_error = sum( abs.( forest_RTS_all_predict - grand_choice_probs  )  ) / size(grand_assortments,1)

			forest_RTS_test_predict, forest_RTS_test_predict_by_class = tcm_predictForest(N, K, lambda_RTS, forest_left_RTS, forest_right_RTS, forest_product_RTS, forest_isLeaf_RTS, test_assortments)
			forest_RTS_test_error = sum( abs.( forest_RTS_test_predict - test_choice_probs  )  ) / size(test_assortments,1)

			row_string = string(N,  ",",
								num_train_assortments, ",",
								d, ",",
								"forest-cold-RTS", ",",
								proc_id, ",",
								attempt, ",",
								num_test_assortments, ",",
								primal_objval_RTS, ",",
								L1_error_RTS, ",",
								forest_RTS_train_error, ",",
								forest_RTS_test_error, ",",
								forest_RTS_all_error, ",",
								K, ",",
								elapsed_time_RTS)
			row_string = string(row_string, "\n")

			# @show row_string

			@spawnat 1 begin
				print(forest_outcsvhandle, row_string)
				flush(forest_outcsvhandle);
			end
		end
	end

	helper_end_time = time()

	return helper_end_time - helper_start_time;
end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
using JLD
using Plots
using LaTeXStrings
using TexTables
import Base.Iterators.product
import Distributions.mean
import Distributions.std
#reset past font scaling
Plots.scalefontsizes()
#scale fonts
Plots.scalefontsizes(1.5)
#Extracting stored varibles
for (maternOrder, dset) in product(["12", "32"], ["Far", "Close", "Circle"])
  ld = load("./out/jld/IncludePred$(dset)$(maternOrder).jld")
  for key in keys(ld)
    symb = Symbol(key * dset * maternOrder)
    data = ld[key]
    @eval $(symb) = $(data)
  end
end

# --------------------
# Start out by creating plots of the data
# --------------------
NTrain = 10000
NTest = 10
d = 2
N = NTrain + NTest

# Prediction points close
xTrain = rand(d, NTrain)
xTest = rand(d, NTest)

plotClose = plot(xlims=(0.0,1.5),
                 ylims=(0.0,1.5),
                 legend=:topleft)
scatter!(plotClose,
         xTrain[1, :],
         xTrain[2, :],
         markershape=:circle,
         color=:blue,
         label=L"\mathrm{Training Points}")

scatter!(plotClose,
         xTest[1, :],
         xTest[2, :],
         markershape=:diamond,
         color=:orange,
         markersize=7.0,
         label=L"\mathrm{Prediction Points}")
savefig(plotClose, "./out/plots/predClose.pdf")

# Prediction points far
xTrain = rand(d, NTrain)
xTest = 1.1 * ones(d,NTest) .+ 0.1 * rand(d, NTest)

plotFar = plot(xlims=(0.0,1.5),
               ylims=(0.0,1.5),
               legend=:topleft)
scatter!(plotFar,
         xTrain[1, :],
         xTrain[2, :],
         markershape=:circle,
         color=:blue,
         label=L"\mathrm{Training Points}")

scatter!(plotFar,
         xTest[1, :],
         xTest[2, :],
         markershape=:diamond,
         color=:orange,
         markersize=7.0,
         label=L"\mathrm{Prediction Points}")
savefig(plotFar, "./out/plots/predFar.pdf")

h = 1 / NTrain
xTrain = h : h : 1
xTrain = vcat(sin.(2 * π * xTrain)',
              2 * cos.(2 * π * xTrain)')
xTest = 0.1 * randn(2,NTest)

# Reducing the size of training set to improve visualization
NTrain = NTrain / 20
# Prediction points circle
h = 1 / NTrain
xTrain = h : h : 1
xTrain = vcat(sin.(2 * π * xTrain)',
              2 * cos.(2 * π * xTrain)')
xTest = 0.1 * randn(2,NTest)

plotCircle = plot(legend=:topleft)
scatter!(plotCircle,
         xTrain[1, :],
         xTrain[2, :],
         markershape=:circle,
         color=:blue,
         label=L"\mathrm{Training\ Points}")

scatter!(plotCircle,
         xTest[1, :],
         xTest[2, :],
         markershape=:diamond,
         color=:orange,
         markersize=7.0,
         label=L"\mathrm{Prediction\ Points}")
savefig(plotCircle, "./out/plots/predCircle.pdf")

# --------------------
# Create tables
# --------------------
outputTables = Array{IndexedTable{1,1},1}(undef, 6)

for (k, (dset, maternOrder)) in enumerate(
                                product(["Close", "Far", "Circle"], ["12", "32"]))

  @eval meanErrorNoPred = $(Symbol("meanErrorNoPred" * dset * maternOrder))
  @eval stdErrorNoPred = $(Symbol("stdErrorNoPred" * dset * maternOrder))
  @eval KLTrueApproxNoPred = $(Symbol("KLTrueApproxNoPred" * dset * maternOrder))
  @eval KLApproxTrueNoPred = $(Symbol("KLApproxTrueNoPred" * dset * maternOrder))
  @eval meanErrorPredFirst = $(Symbol("meanErrorPredFirst" * dset * maternOrder))
  @eval stdErrorPredFirst = $(Symbol("stdErrorPredFirst" * dset * maternOrder))
  @eval KLTrueApproxPredFirst = $(Symbol("KLTrueApproxPredFirst" * dset * maternOrder))
  @eval KLApproxTruePredFirst = $(Symbol("KLApproxTruePredFirst" * dset * maternOrder))
  @eval meanErrorPredLast = $(Symbol("meanErrorPredLast" * dset * maternOrder))
  @eval stdErrorPredLast = $(Symbol("stdErrorPredLast" * dset * maternOrder))
  @eval KLTrueApproxPredLast = $(Symbol("KLTrueApproxPredLast" * dset * maternOrder))
  @eval KLApproxTruePredLast = $(Symbol("KLApproxTruePredLast" * dset * maternOrder))

  meanTabNoPred = TableCol("\$\\mathrm{$(dset)}, \\nu = $(maternOrder[1]) / $(maternOrder[2])\$",
                     [Symbol("\$\\left|\\mu-\\mu_{\\NP}\\right|\$")],
                     [mean(meanErrorNoPred)],
                     [std(meanErrorNoPred)],
                     )

  meanTabPredFirst = TableCol("\$\\mathrm{$(dset)}, \\nu = $(maternOrder[1]) / $(maternOrder[2])\$",
                     [Symbol("\$\\left|\\mu-\\mu_{\\PF}\\right|\$")],
                     [mean(meanErrorPredFirst)],
                     [std(meanErrorPredFirst)],
                     )

  meanTabPredLast = TableCol("\$\\mathrm{$(dset)}, \\nu = $(maternOrder[1]) / $(maternOrder[2])\$",
                     [Symbol("\$\\left|\\mu-\\mu_{\\PL}\\right|\$")],
                     [mean(meanErrorPredLast)],
                     [std(meanErrorPredLast)],
                     )

  stdTabNoPred = TableCol("\$\\mathrm{$(dset)}, \\nu = $(maternOrder[1]) / $(maternOrder[2])\$",
                       [Symbol(:"\$\\left|\\sigma^2-\\sigma_{\\NP}^2\\right|\$")],
                       [mean(stdErrorNoPred)],
                       # [std(stdErrorNoPred)],
                       )

  stdTabPredFirst = TableCol("\$\\mathrm{$(dset)}, \\nu = $(maternOrder[1]) / $(maternOrder[2])\$",
                       [Symbol("\$\\left|\\sigma^2-\\sigma_{\\PF}^2\\right|\$")],
                       [mean(stdErrorPredFirst)],
                       # [std(stdErrorPredFirst)],
                       )

  stdTabPredLast = TableCol("\$\\mathrm{$(dset)}, \\nu = $(maternOrder[1]) / $(maternOrder[2])\$",
                       [Symbol("\$\\left|\\sigma^2-\\sigma_{\\PL}^2\\right|\$")],
                       [mean(stdErrorPredLast)],
                       # [std(stdErrorPredLast)],
                       )

  KLTrueAppTabNoPred = TableCol("\$\\mathrm{$(dset)}, \\nu = $(maternOrder[1]) / $(maternOrder[2])\$",
                     [Symbol("\$\\KL*{\\KM}{\\KM_{\\NP}}\$")],
                     [mean(KLTrueApproxNoPred)],
                     # [std(KLTrueApproxNoPred)],
                     )

  KLTrueAppTabPredFirst = TableCol("\$\\mathrm{$(dset)}, \\nu = $(maternOrder[1]) / $(maternOrder[2])\$",
                     [Symbol("\$\\KL*{\\KM}{\\KM_{\\PF}}\$")],
                     [mean(KLTrueApproxPredFirst)],
                     # [std(KLTrueApproxPredFirst)],
                     )

  KLTrueAppTabPredLast = TableCol("\$\\mathrm{$(dset)}, \\nu = $(maternOrder[1]) / $(maternOrder[2])\$",
                     [Symbol("\$\\KL*{\\KM}{\\KM_{\\PL}}\$")],
                     [mean(KLTrueApproxPredLast)],
                     # [std(KLTrueApproxPredLast)],
                     )

  KLAppTrueTabNoPred = TableCol("\$\\mathrm{$(dset)}, \\nu = $(maternOrder[1]) / $(maternOrder[2])\$",
                     [Symbol("\$\\KL*{\\KM_{\\NP}}{\\KM}\$")],
                     [mean(KLApproxTrueNoPred)],
                     # [std(KLApproxTrueNoPred)],
                     )

  KLAppTrueTabPredFirst = TableCol("\$\\mathrm{$(dset)}, \\nu = $(maternOrder[1]) / $(maternOrder[2])\$",
                     [Symbol("\$\\KL*{\\KM_{\\PF}}{\\KM}\$")],
                     [mean(KLApproxTruePredFirst)],
                     # [std(KLApproxTruePredFirst)],
                     )

  KLAppTrueTabPredLast = TableCol("\$\\mathrm{$(dset)}, \\nu = $(maternOrder[1]) / $(maternOrder[2])\$",
                     [Symbol("\$\\KL*{\\KM_{\\PL}}{\\KM}\$")],
                     [mean(KLApproxTruePredLast)],
                     # [std(KLApproxTruePredLast)],
                     )

  # Setting the format

  setfield!.(meanTabNoPred.data.vals, :format, "{:,2e}")
  setfield!.(meanTabNoPred.data.vals, :format_se, "{:,2e}")

  setfield!.(meanTabPredFirst.data.vals, :format, "{:,2e}")
  setfield!.(meanTabPredFirst.data.vals, :format_se, "{:,2e}")

  setfield!.(meanTabPredLast.data.vals, :format, "{:,2e}")
  setfield!.(meanTabPredLast.data.vals, :format_se, "{:,2e}")

  setfield!.(stdTabNoPred.data.vals, :format, "{:,2e}")
  # setfield!.(stdTabNoPred.data.vals, :format_se, "{:,2e}")

  setfield!.(stdTabPredFirst.data.vals, :format, "{:,2e}")
  # setfield!.(stdTabPredFirst.data.vals, :format_se, "{:,2e}")

  setfield!.(stdTabPredLast.data.vals, :format, "{:,2e}")
  # setfield!.(stdTabPredLast.data.vals, :format_se, "{:,2e}")

  setfield!.(KLTrueAppTabNoPred.data.vals, :format, "{:,2e}")
  # setfield!.(KLTrueAppTabNoPred.data.vals, :format_se, "{:,2e}")

  setfield!.(KLTrueAppTabPredFirst.data.vals, :format, "{:,2e}")
  # setfield!.(KLTrueAppTabPredFirst.data.vals, :format_se, "{:,2e}")

  setfield!.(KLTrueAppTabPredLast.data.vals, :format, "{:,2e}")
  # setfield!.(KLTrueAppTabPredLast.data.vals, :format_se, "{:,2e}")

  setfield!.(KLAppTrueTabNoPred.data.vals, :format, "{:,2e}")
  # setfield!.(KLAppTrueTabNoPred.data.vals, :format_se, "{:,2e}")

  setfield!.(KLAppTrueTabPredFirst.data.vals, :format, "{:,2e}")
  # setfield!.(KLAppTrueTabPredFirst.data.vals, :format_se, "{:,2e}")

  setfield!.(KLAppTrueTabPredLast.data.vals, :format, "{:,2e}")
  # setfield!.(KLAppTrueTabPredLast.data.vals, :format_se, "{:,2e}")

  outputTables[k] = vcat(meanTabNoPred,
                         meanTabPredFirst,
                         meanTabPredLast,
                         stdTabNoPred,
                         stdTabPredFirst,
                         stdTabPredLast,
                         KLTrueAppTabNoPred,
                         KLTrueAppTabPredFirst,
                         KLTrueAppTabPredLast,
                         KLAppTrueTabNoPred,
                         KLAppTrueTabPredFirst,
                         KLAppTrueTabPredLast)
end
tab = hcat(outputTables[:]...)
f = open("./out/tables/IncludePredTables.tex", "w")
write(f, to_tex(tab))
close(f)


text:
using DelimitedFiles

cola_init() = register(DataDep(
    "GLUE-CoLA",
    """
    The Corpus of Linguistic Acceptability (CoLA) task (GLUE version)
    """,
    "https://dl.fbaipublicfiles.com/glue/data/CoLA.zip",
    "f212fcd832b8f7b435fb991f101abf89f96b933ab400603bf198960dfc32cbff";
    post_fetch_method = fn -> begin
      mv(fn, "CoLA.zip")
      DataDeps.unpack("CoLA.zip")
      innerdir = "CoLA"
      innerfiles = readdir(innerdir)
      mv.(joinpath.(innerdir, innerfiles), innerfiles)
      rm(innerdir)
    end
))

struct CoLA <: Dataset end

function trainfile(::CoLA)
    sets = readdlm(datadep"GLUE-CoLA/train.tsv", '\t', String)
    [selectdim(sets, 2, i) for i = (4,2)]
end

function devfile(::CoLA)
    sets = readdlm(datadep"GLUE-CoLA/dev.tsv", '\t', String)
    [selectdim(sets, 2, i) for i = (4,2)]
end

function testfile(::CoLA)
    sets, headers = readdlm(datadep"GLUE-CoLA/test.tsv", '\t', String; header=true)
    [selectdim(sets, 2, 2)]
end

get_labels(::CoLA) = ("0", "1")


text:
# ------------ Branch trace ------------ #

struct BranchTrace{T <: RecordSite, B <: RecordSite} <: Trace
    cond::T
    branch::B
    BranchTrace(cond::T, branch::B) where {T, B} = new{T, B}(cond, branch)
end

# If-else branch site
struct ConditionalBranchCallSite{C, A, J, L, R} <: CallSite
    trace::BranchTrace
    score::Float64
    cond_kernel::C
    cond_args::J
    cond::Bool
    branch::A
    branch_args::L
    ret::R
end
get_score(cbs::ConditionalBranchCallSite) = cbs.score


text:
# ---
# title: 1416. Restore The Array
# id: problem1416
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/restore-the-array/description/>
# hidden: true
# ---
#
# A program was supposed to print an array of integers. The program forgot to
# print whitespaces and the array is printed as a string of digits and all we
# know is that all integers in the array were in the range `[1, k]` and there
# are no leading zeros in the array.
#
# Given the string `s` and the integer `k`. There can be multiple ways to
# restore the array.
#
# Return _the number of possible array_ that can be printed as a string `s`
# using the mentioned program.
#
# The number of ways could be very large so return it **modulo** `10^9 + 7`
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "1000", k = 10000
#     Output: 1
#     Explanation: The only possible array is [1000]
#
#
# **Example 2:**
#
#
#
#     Input: s = "1000", k = 10
#     Output: 0
#     Explanation: There cannot be an array that was printed this way and has all integer >= 1 and <= 10.
#
#
# **Example 3:**
#
#
#
#     Input: s = "1317", k = 2000
#     Output: 8
#     Explanation: Possible arrays are [1317],[131,7],[13,17],[1,317],[13,1,7],[1,31,7],[1,3,17],[1,3,1,7]
#
#
# **Example 4:**
#
#
#
#     Input: s = "2020", k = 30
#     Output: 1
#     Explanation: The only possible array is [20,20]. [2020] is invalid because 2020 > 30. [2,020] is ivalid because 020 contains leading zeros.
#
#
# **Example 5:**
#
#
#
#     Input: s = "1234567890", k = 90
#     Output: 34
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 10^5`.
#   * `s` consists of only digits and doesn't contain leading zeros.
#   * `1 <= k <= 10^9`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
$(SIGNATURES)

Create voronoi diagram of point set.

Returns a mesh of triangles.
"""
function voronoi(points::Vector{Point{3, T}}) where T <: AbstractFloat
    result = tetrahedralize(JLTetGenIO(points), "Qw")
    Mesh{Triangle}(result)
end

"""
$(SIGNATURES)

Tetrahedralize a mesh of polygons with optional facet markers.
Returns a mesh of tetrahdra.
"""
function TetGen.tetrahedralize(
        mesh::Mesh{3, Float64, <: TetGen.Ngon}, command = "Qp";
        marker = :markers, holes = Point{3, Float64}[]
    )
    f = faces(mesh)
    kw_args = Any[:facets => metafree(f),:holes => holes]
    if hasproperty(f, marker)
        push!(kw_args, :facetmarkers => getproperty(f, marker))
    end
    tio = JLTetGenIO(coordinates(mesh); kw_args...)
    result = tetrahedralize(tio, command)
    return Mesh{Tetrahedron}(result)
end

"""
$(SIGNATURES)

Tetrahedralize a domain described by a mesh of triangles.
Returns a mesh of tetrahdra.
"""
function TetGen.tetrahedralize(mesh::Mesh{3, Float64, <: TetGen.Triangle}, command = "Qp")
    tio = JLTetGenIO(coordinates(mesh); facets = faces(mesh))
    result = tetrahedralize(tio, command)
    Mesh{Tetrahedron}(result)
end


text:
# """
#     Object
#
# Overarching type for objects displayable on `Axes3D`.
# """
# abstract type Object3D <: Object end


text:
#=
day14:
- Julia version: 1.7.0
- Author: Paul.Mealor
- Date: 2021-12-14
=#

using AoC, Test
using StructEquality # https://docs.juliahub.com/StructEquality/TwsrV/1.0.0/

@def_structequal struct Rule{T}
    target::T
    results::Vector{T}
end

function parseRule(line)
    captures = match(r"^((.)(.)) -> (.)", line).captures
    return Rule(String(captures[1]),[captures[2]*captures[4],captures[4]*captures[3]])
end
@test parseRule("AB -> C") == Rule("AB",["AC","CB"])

function toPairs(line::T) where T
    result = Dict{T,BigInt}()
    for i in 1:length(line)-1
        pair = line[i:i+1]
        result[pair] = get(result,pair,0) + 1
    end
    return result
end
@test toPairs("ABCD") == Dict("AB"=>1,"BC"=>1,"CD"=>1)

parseInput(lines) = (toPairs(lines[1]), parseRule.(lines[3:end]))
@test parseInput(["ABCD","","AB -> C", "AC -> D"]) == (Dict("AB"=>1,"BC"=>1,"CD"=>1), [Rule("AB",["AC","CB"]),Rule("AC",["AD","DC"])])

doStep(input::Tuple) = doStep(input...)
function doStep(state::T, rules) where T
    newState = T()
    for rule in rules
        elementCount = get(state,rule.target,0)
        if elementCount != 0
            for sequence in rule.results
                newState[sequence] = get(newState,sequence,0) + elementCount
            end
        end
    end
    return newState
end
@test doStep(parseInput(exampleLines(14,1))) == Dict("NC"=>1,"CN"=>1,"NB"=>1,"BC"=>1,"CH"=>1,"HB"=>1)

doSteps(input::Tuple, n) = doSteps(input..., n)
function doSteps(state, rules, n)
    for i in 1:n
        state = doStep(state, rules)
    end
    return state
end
@test doSteps(parseInput(exampleLines(14,1)),1) == toPairs("NCNBCHB")
@test doSteps(parseInput(exampleLines(14,1)),2) == toPairs("NBCCNBBBCBHCB")
@test doSteps(parseInput(exampleLines(14,1)),3) == toPairs("NBBBCNCCNBBNBNBBCHBHHBCHB")
@test doSteps(parseInput(exampleLines(14,1)),4) == toPairs("NBBNBNBBCCNBCNCCNBBNBBNBBBNBBNBBCBHCBHHNHCBBCBHCB")

function answer(lines, steps)
    lastElement = lines[1][end]
    state = doSteps(parseInput(lines), steps)
    counts = Dict{Char, BigInt}(lastElement=>1)
    for (pair,count) in state
        counts[pair[1]] = get(counts,pair[1],0)+count
    end
    mostElementCount = argmax(pair->pair[2], counts)[2]
    leastElementCount = argmin(pair->pair[2], counts)[2]
    return mostElementCount - leastElementCount
end
part1(lines) = answer(lines, 10)
part2(lines) = answer(lines, 40)
@test part1(exampleLines(14,1)) == 1588
@test part2(exampleLines(14,1)) == 2188189693529
@show "testing complete"

@show lines(14) |> ll -> @time part1(ll)
@show lines(14) |> ll -> @time part2(ll)
@show lines(14) |> ll -> @time answer(ll,400)
@show lines(14) |> ll -> @time answer(ll,4000)


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

struct SearchApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_searchMessages(_api::SearchApi, query::String; count=nothing, highlight=nothing, page=nothing, sort=nothing, sort_dir=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/search.messages", ["slackAuth"])
    Swagger.set_param(_ctx.query, "count", count)  # type Int32
    Swagger.set_param(_ctx.query, "highlight", highlight)  # type Bool
    Swagger.set_param(_ctx.query, "page", page)  # type Int32
    Swagger.set_param(_ctx.query, "query", query)  # type String
    Swagger.set_param(_ctx.query, "sort", sort)  # type String
    Swagger.set_param(_ctx.query, "sort_dir", sort_dir)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Searches for messages matching a query.
Param: query::String (required)
Param: count::Int32
Param: highlight::Bool
Param: page::Int32
Param: sort::String
Param: sort_dir::String
Return: Dict{String, Any}
"""
function searchMessages(_api::SearchApi, query::String; count=nothing, highlight=nothing, page=nothing, sort=nothing, sort_dir=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_searchMessages(_api, query; count=count, highlight=highlight, page=page, sort=sort, sort_dir=sort_dir, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function searchMessages(_api::SearchApi, response_stream::Channel, query::String; count=nothing, highlight=nothing, page=nothing, sort=nothing, sort_dir=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_searchMessages(_api, query; count=count, highlight=highlight, page=page, sort=sort, sort_dir=sort_dir, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export searchMessages


text:
using ProtoBuf

print_hdr(tname) = println("testing $tname...")

type V1
    f1::Int32
    f2::Bool
    V1() = (a=new(); clear(a); a)
    V1(f1,f2) = new(f1,f2)
end

type V2
    f1::Int64
    f2::Bool
    f3::Int64
    V2() = (a=new(); clear(a); a)
    V2(f1,f2,f3) = new(f1,f2,f3)
end

function check_samestruct()
    iob = PipeBuffer();
    writeval = V1(1, true);
    writeproto(iob, writeval);
    readval = readproto(iob, V1());
    @assert ProtoBuf.protoeq(readval, writeval)

    iob = PipeBuffer();
    writeval = V2(1, true, 20);
    writeproto(iob, writeval);
    readval = readproto(iob, V2());
    @assert ProtoBuf.protoeq(readval, writeval)
end

# write V1, read V2
# write V2, read V1
function check_V1_v2()
    iob = PipeBuffer();
    writeval = V1(1, true);
    writeproto(iob, writeval);
    readval = readproto(iob, V2());
    checkval = V2(1,true,0)
    @assert ProtoBuf.protoeq(readval, checkval)

    iob = PipeBuffer();
    writeval = V2(1, true, 20);
    writeproto(iob, writeval);
    readval = readproto(iob, V1());
    checkval = V1(1,true)
    @assert ProtoBuf.protoeq(readval, checkval)

    iob = PipeBuffer();
    writeval = V2(typemax(Int64), true, 20);
    writeproto(iob, writeval);
    readval = readproto(iob, V1());
    checkval = V1(0,true)
    @assert ProtoBuf.protoeq(readval, checkval)
end

print_hdr("serialization across type versions")
check_samestruct()
check_V1_v2()


text:
using Test
include("makeff.jl")

@testset "baseff" begin
  @test :(f(x)) == baseff(:x, :f, 1)
  @test :(f(x,x)) == baseff(:x, :f, 2)
  @test :(f(x,x,x)) == baseff(:x, :f, 3)
  @test :(P(x,x,x)) == baseff(:x, :P, 3)
end

@testset "replace1" begin
  @test :(P(f(x,x,x,x), f(x,x,x,x))) == replace1(:x, :(P(x,x)), :(f(x,x,x,x)))
  @test :(P(f(x,x,x,x), g(x), f(x,x,x,x))) == replace1(:x, :(P(x,g(x),x)), :(f(x,x,x,x)))

end

@testset "makedd" begin
  @test :(f(x)) == makedd(:x, :(f(x)), :x, 0, 1)
  @test :(f(f(x))) == makedd(:x, :(f(x)), :x, 1, 1)
  @test :(f(f(f(x)))) == makedd(:x, :(f(x)), :x, 2, 1)
  @test :(f(f(f(f(x))))) == makedd(:x, :(f(x)), :x, 3, 1)
end

@testset "makeff(1,*)" begin
  @test :(P()) == makeff(:x, :P, :f, 1, 0)
  @test :(P(f(f(x)))) == makeff(:x, :P, :f, 1, 1)
  @test :(P(f(f(x,x),f(x,x)),f(f(x,x),f(x,x)))) == makeff(:x, :P, :f, 1, 2)
  @test :(P(f(f(x,x,x),f(x,x,x),f(x,x,x)),f(f(x,x,x),f(x,x,x),f(x,x,x)),f(f(x,x,x),f(x,x,x),f(x,x,x)))) == makeff(:x, :P, :f, 1, 3)
end

@testset "makeff(*,1)" begin
  @test :(P(x)) == makeff(:x, :P, :f, 0, 1)
  @test :(P(f(f(x)))) == makeff(:x, :P, :f, 1, 1)
  @test :(P(f(f(x)))) == makeff(:x, :P, :f, 2, 1)
  @test :(P(f(f(f(x))))) == makeff(:x, :P, :f, 3, 1)
  @test :(P(f(f(f(f(x)))))) == makeff(:x, :P, :f, 4, 1)
end

@testset "makeff(*,*)" begin
  @test :(P(f(f(x,x),f(x,x)),f(f(x,x),f(x,x)))) == makeff(:x, :P, :f, 2, 2)
  @test :(P(f(f(x,x,x),f(x,x,x),f(x,x,x)),f(f(x,x,x),f(x,x,x),f(x,x,x)),f(f(x,x,x),f(x,x,x),f(x,x,x)))) == makeff(:x, :P, :f, 2, 3)

  @test :(P(f(f(f(x,x),f(x,x)),f(f(x,x),f(x,x))),f(f(f(x,x),f(x,x)),f(f(x,x),f(x,x))))) == makeff(:x, :P, :f, 3, 2)
  @test :(P(f(f(f(x, x, x), f(x, x, x), f(x, x, x)), f(f(x, x, x), f(x, x, x), f(x, x, x)), f(f(x, x, x), f(x, x, x), f(x, x, x))), f(f(f(x, x, x), f(x, x, x), f(x, x, x)), f(f(x, x, x), f(x, x, x), f(x, x, x)), f(f(x, x, x), f(x, x, x), f(x, x, x))), f(f(f(x, x, x), f(x, x, x), f(x, x, x)), f(f(x, x, x), f(x, x, x), f(x, x, x)), f(f(x, x, x), f(x, x, x), f(x, x, x))))) == makeff(:x, :P, :f, 3, 3)

end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
#=
jim.jl
jiffy image display
2019-02-23 Jeff Fessler, University of Michigan
=#

export jim

using UnitfulRecipes
using Plots: heatmap, plot, plot!, Plot
import Plots # gui
using MosaicViews: mosaicview
using FFTViews: FFTView
using OffsetArrays: OffsetMatrix
import OffsetArrays # no_offset_view
using AxisArrays: AxisArray, axisnames, axisvalues
#using MIRTjim: prompt

# global default key/values
const jim_table = Dict([
 :aspect_ratio => :equal,
 :clim => nothing,
 :color => :grays,
 :colorbar => :legend,
 :gui => false,
 :prompt => false,
 :line3plot => true, # lines around sub image for 3d mosaic?
 :line3type => (:yellow),
 :ncol => 0,
 :nrow => 0,
 :padval => nothing,
 :mosaic_npad => 1,
 :tickdigit => 1,
 :title => "",
 :xlabel => nothing,
 :ylabel => nothing,
 :fft0 => false,
 :yflip => nothing, # defer to minimum value of y - see default below
 :yreverse => nothing, # defer to whether y is non-ascending
 :abswarn => isinteractive(), # warn when taking abs of complex images?
 :infwarn => isinteractive(), # warn when image has Inf?
 :nanwarn => isinteractive(), # warn when image has NaN?
])

jim_def = deepcopy(jim_table)

jim_stack = Any[] # for push! and pop!

# min,max of an iterable, excluding Inf, NaN
# caution: complex Unitful arrays are <:Number not <:Complex
function _dogood(z, fun::Function, init::Real)
    T = eltype(first(z))
    z = Iterators.filter(isfinite, z)
    z = Iterators.map(x -> x / oneunit(T), z)
    if first(z) isa Complex
        z = Iterators.map(abs, z)
    end
    return fun(z; init) * oneunit(real(T))
end
_maxgood(z::AbstractArray{<:Number}) = _dogood(z, maximum, -Inf)
_maxgood(z::AbstractArray{<:AbstractArray}) = maximum(_maxgood, z)
_mingood(z::AbstractArray{<:Number}) = _dogood(z, minimum, Inf)
_mingood(z::AbstractArray{<:AbstractArray}) = minimum(_mingood, z)

"""
    nothing_else(x, y)
Return `y` if `x` is nothing, else return `x`.
"""
function nothing_else(x, y)
    return x == nothing ? y : x
end

# parsimonious axis ticks by default
function _ticks(x::AbstractVector{<:Number})
    if x[1] isa Real
        minfloor = x -> floor(minimum(x), digits = jim_def[:tickdigit])
        maxceil = x -> ceil(maximum(x), digits = jim_def[:tickdigit])
    else # Unitful
        minfloor = x -> floor(eltype(x), minimum(x), digits = jim_def[:tickdigit])
        maxceil = x -> ceil(eltype(x), maximum(x), digits = jim_def[:tickdigit])
    end
    z0 = zero(x[1]) # units
    ticks = (minimum(x) < z0 < maximum(x)) ?
        [minfloor(x), z0, maxceil(x)] :
        [minfloor(x), maxceil(x)]
    return ticks
end

# subtle issues with default labels depending on Unitful or not
_label(s::Symbol, x::AbstractVector{<:Real}) =
    isnothing(jim_def[s]) ? "" : jim_def[s]
_label(s::Symbol, x::AbstractVector{<:Number}) = jim_def[s]

_fft0_axis(n::Int) = (-n÷2):(n÷2 - iseven(n))

"""
    jim(z, ...)

A jiffy image display of `z` using `heatmap`.

in
- `z` image, can be 2D or higher, if higher then it uses `mosaicviews`

option
- `aspect_ratio`; default `:equal`
- `clim`; default `(minimum(z),maximum(z))`
- `color` (colormap, e.g. `:hsv`); default `:grays`
- `colorbar` (e.g. `:none`); default `:legend`
- `gui` call `Plots.gui()` immediately?; default `false`
- `prompt` call `prompt()` immediately?; default `false`
- `ncol` for mosaicview for 3D and higher arrays; default `0` does auto select
   "Number of tiles in column direction."
- `nrow` for mosaicview for 3D and higher arrays; default `0` does auto select
- `padval` padding value for mosaic view; default `minimum(z)`
- `line3plot` lines around sub image for 3d mosaic; default `true`
- `line3type` line type around sub image for 3d mosaic; default `(:yellow)`
- `mosaic_npad` # of pixel padding for mosaic view; default `1`
- `fft0` if true use FFTView to display (2D only); default `false`
- `title`; default `""`
- `xlabel`; default `nothing` (or units if applicable)
- `ylabel`; default `nothing`
- `yflip`; default `true` if `minimum(y) >= 0`
- `yreverse`; default `true` if `y[1] > y[end]`
- `x` values for x axis; default `collect(axes(z)[1])`
- `y` values for y axis; default `collect(axes(z)[2])`
- `xticks`; default `[minimum(x),maximum(x)]` (usually)
- `yticks`; default `[minimum(y),maximum(y)]`

out
- returns plot handle, type `Plots.Plot`
"""
function jim(
    z::AbstractMatrix{<:Number} ;
    abswarn::Bool = jim_def[:abswarn],
    kwargs...,
)

    !any(isfinite, z) && throw("no finite values")

    if all(z -> imag(z) == zero(real(z[1])), z)
        z = real(z) # remove any zero imaginary part
    else
        abswarn && (@warn "magnitude at $(caller_name())")
        z = abs.(z) # due to Unitful complex types not being <: Complex
    end
    return _jim(z ; kwargs...)
end

# 2D RealU matrix
function _jim(z::AbstractMatrix{<:RealU} ;
    aspect_ratio = jim_def[:aspect_ratio],
    clim = nothing_else(jim_def[:clim], (_mingood(z), _maxgood(z))),
    color = jim_def[:color],
    colorbar = jim_def[:colorbar],
    title::AbstractString = jim_def[:title],
    fft0::Bool = jim_def[:fft0],
    gui::Bool = jim_def[:gui],
    prompt::Bool = jim_def[:prompt],
    x::AbstractVector{<:Number} = fft0 ? _fft0_axis(size(z,1)) : axes(z,1),
    y::AbstractVector{<:Number} = fft0 ? _fft0_axis(size(z,2)) : axes(z,2),
    xy::Tuple = (x,y),
    xticks = _ticks(x),
    yticks = _ticks(y),
    xlabel::Union{Nothing,AbstractString} = _label(:xlabel, x),
    ylabel::Union{Nothing,AbstractString} = _label(:ylabel, y),
    yflip::Bool = nothing_else(jim_def[:yflip], minimum(y) >= zero(y[1])),
    yreverse::Bool = nothing_else(jim_def[:yreverse], y[1] > y[end]),
#   abswarn::Bool = jim_def[:abswarn], # not used here
    infwarn::Bool = jim_def[:infwarn],
    nanwarn::Bool = jim_def[:nanwarn],
    kwargs...
)

    # because some backends require y to be in ascending order
    if yreverse
        y = reverse(y)
        z = reverse(z, dims=2)
    end

    if fft0
        z = FFTView(z)[x,y]
    end

    # warnings for non-number values
    infwarn && any(isinf, z) && @warn("$(sum(isinf, z)) ±Inf")
    nanwarn && any(isnan, z) && @warn("$(sum(isnan, z)) NaN")

    if _mingood(z) ≈ _maxgood(z) # uniform or nearly uniform image
        tmp = (_mingood(z) == _maxgood(z)) ? "Uniform $(z[1])" :
            "Nearly uniform $((_mingood(z),_maxgood(z)))"

        p = plot( ; aspect_ratio,
            xlim = (x[1], x[end]),
            ylim = (y[1], y[end]),
            title,
            yflip,
            xlabel,
            ylabel,
            xticks,
            yticks,
            annotate = (x[(end+1)÷2], y[(end+1)÷2], tmp, :red),
            kwargs...
        )

    else

        p = heatmap(xy..., z' ;
            transpose = false,
            aspect_ratio,
            clim,
            color,
            colorbar,
            title,
            yflip,
            xlabel,
            ylabel,
            xticks,
            yticks,
            kwargs...
        )
    end

    gui && Plots.gui()
    prompt && MIRTjim.prompt()
    return p
end # jim

# OffsetArray / OffsetMatrix
# https://github.com/JuliaPlots/Plots.jl/issues/2410
_axes(z,j) = axes(z,j).parent .+ axes(z,j).offset
function jim(z::OffsetMatrix{<:Number} ;
    x = _axes(z,1),
    y = _axes(z,2),
    kwargs...
)
    jim(OffsetArrays.no_offset_view(z) ; x, y, kwargs...)
end
#=
This approach fails because z' is no longer an OffsetMatrix
function Plots.heatmap(z::OffsetMatrix{<:Number}; kwargs...)
    x = axes(z,1); x = x.parent .+ x.offset
    y = axes(z,2); y = y.parent .+ y.offset
    heatmap(x, y, z; kwargs...)
end
function Plots.heatmap(x, y, z::OffsetMatrix{<:Number}; kwargs...)
    heatmap(x, y, OffsetArrays.no_offset_view(z); kwargs...)
end
=#

# AxisArray / AxisMatrix
const AxisMatrix{T} = AxisArray{T,2,D,Ax} where {T,D,Ax}
function jim(z::AxisMatrix{<:Number} ;
    x = axisvalues(z)[1],
    y = axisvalues(z)[2],
    xlabel = String(axisnames(z)[1]),
    ylabel = String(axisnames(z)[2]),
    kwargs...
)
    jim(parent(z) ; x, y, xlabel, ylabel, kwargs...)
end

"""
    jim(z, title::String ; kwargs...)
"""
jim(z::AbstractArray, title::AbstractString ; kwargs...) =
    jim(z ; title, kwargs...)

"""
    jim(x, y, z ; kwargs...)

The `x` and `y` axes can be Unitful thanks to UnitfulRecipes.
"""
jim(x::AbstractVector{<:Number}, y::AbstractVector{<:Number}, z ; kwargs...) =
    jim(z ; x, y, kwargs...)

"""
    jim(x, y, z, title::String ; kwargs...)

Allow `title` as positional argument for convenience.
"""
jim(x::AbstractVector, y, z, title::AbstractString ; kwargs...) =
    jim(x, y, z ; title, kwargs...)

"""
    jim(key::Symbol, value::Any)
Set default value for one of the keys.
"""
function jim(key::Symbol, value)
    global jim_def
    !haskey(jim_def, key) && throw(ArgumentError("no key $key"))
    jim_def[key] = value
end

"""
    jim(plot1, plot2, ... ; gui=?, prompt=?, kwargs...)
Subplot-type layout, where `kwargs` are passed to `plot`.
- `gui` call `Plots.gui()` immediately?; default `false`
- `prompt` call `prompt()` immediately?; default `false`
"""
function jim(p::Plot, args... ;
    gui::Bool = jim_def[:gui],
    prompt::Bool = jim_def[:prompt],
    kwargs...,
)
    out = plot(p, args... ; kwargs...)
    gui && Plots.gui()
    prompt && MIRTjim.prompt()
    return out
end

"""
    jim(:test::Symbol)

`jim(:keys)` return default keys.

`jim(:defs)` return `Dict` of default keys / vals.

`jim(:key)` return `Dict[key]` if possible.

`jim(:reset)` reset to defaults.

`jim(:push!)` `push!` current defaults to the stack.

`jim(:pop!)` `pop!` defaults from the stack.

`jim(:blank)` return blank plot.
"""
function jim(test::Symbol)
    global jim_def
    if test === :keys
        return sort(collect(keys(jim_def)))
    end
    if test === :blank
        return plot(legend=false, grid=false, foreground_color_subplot=:white)
    end
    if test === :defs
        return jim_def
    end
    if test === :push!
        return push!(jim_stack, deepcopy(jim_def))
    end
    if test === :pop!
        return jim_def = deepcopy(pop!(jim_stack))
    end
    if test === :reset
        return jim_def = deepcopy(jim_table)
    end
    if haskey(jim_def, test)
        return jim_def[test]
    end

    throw("symbol $test")
end


max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/TreeKnit

max_forks_repo_name: PierreBarrat/TreeKnit


text:
module pretrained
import ..pretrained_dir
include(joinpath(pretrained_dir, "v2_tiny", "voc2007", "v2_tiny_voc.jl"))
end #module