text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

using Tensors, Materials, Test

# TODO: Civilized way to place this wall of numbers at the end of the source file?
# TODO: Or maybe read the expected results from a text file, like the other tests do?
let stresses_expected = [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [10.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [20.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [30.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [40.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [50.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [60.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [70.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [80.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [90.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [100.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [110.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [119.995, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [129.743, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [137.812, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [143.250, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [146.767, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [149.313, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [151.463, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [153.465, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [155.402, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [157.301, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [159.164, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [160.991, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [162.780, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [164.529, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [166.236, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [167.902, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [169.527, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [171.110, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [172.653, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [174.155, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [175.618, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [177.042, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [178.429, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [179.779, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [181.094, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [182.373, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [183.619, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [184.833, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [186.014, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [10.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [20.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [30.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [40.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [50.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [60.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [70.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [80.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [90.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [100.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [110.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [120.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [130.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [140.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [150.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [159.999, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [169.920, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [179.177, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [187.533, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [195.357, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [202.848, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [209.935, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [216.120, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [219.606, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [216.262, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [209.090, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [203.945, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [201.241, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [200.038, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [199.668, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [199.760, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [200.115, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [200.622, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [201.219, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [201.867, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [202.542, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [203.231, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [203.922, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [204.611, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [205.293, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [205.967, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [206.630, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [207.282, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [207.922, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [208.551, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [209.168, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [209.774, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [210.368, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [210.952, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [211.525, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [212.087, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [212.640, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [213.183, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [213.717, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [214.242, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [214.758, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [215.266, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [215.766, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [216.259, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [216.743, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [217.221, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [217.691, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [218.155, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [218.612, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [219.063, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [219.507, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [219.946, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [220.379, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [220.807, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [221.229, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [221.646, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [222.059, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [222.466, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [222.868, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [223.266, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [223.660, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [224.049, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [224.435, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [224.816, 0.0, 0.0, 0.0, 0.0, 0.0]],
    strains_expected = [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                        [0.00005, -0.00002, -0.00002, 0.0, 0.0, 0.0],
                        [0.00010, -0.00003, -0.00003, 0.0, 0.0, 0.0],
                        [0.00015, -0.00004, -0.00004, 0.0, 0.0, 0.0],
                        [0.00020, -0.00006, -0.00006, 0.0, 0.0, 0.0],
                        [0.00025, -0.00008, -0.00008, 0.0, 0.0, 0.0],
                        [0.00030, -0.00009, -0.00009, 0.0, 0.0, 0.0],
                        [0.00035, -0.00010, -0.00010, 0.0, 0.0, 0.0],
                        [0.00040, -0.00012, -0.00012, 0.0, 0.0, 0.0],
                        [0.00045, -0.00014, -0.00014, 0.0, 0.0, 0.0],
                        [0.00050, -0.00015, -0.00015, 0.0, 0.0, 0.0],
                        [0.00055, -0.00017, -0.00017, 0.0, 0.0, 0.0],
                        [0.00060, -0.00018, -0.00018, 0.0, 0.0, 0.0],
                        [0.00065, -0.00020, -0.00020, 0.0, 0.0, 0.0],
                        [0.00070, -0.00021, -0.00021, 0.0, 0.0, 0.0],
                        [0.00075, -0.00023, -0.00023, 0.0, 0.0, 0.0],
                        [0.00080, -0.00025, -0.00025, 0.0, 0.0, 0.0],
                        [0.00085, -0.00028, -0.00028, 0.0, 0.0, 0.0],
                        [0.00090, -0.00030, -0.00030, 0.0, 0.0, 0.0],
                        [0.00095, -0.00032, -0.00032, 0.0, 0.0, 0.0],
                        [0.00100, -0.00034, -0.00034, 0.0, 0.0, 0.0],
                        [0.00105, -0.00037, -0.00037, 0.0, 0.0, 0.0],
                        [0.00110, -0.00039, -0.00039, 0.0, 0.0, 0.0],
                        [0.00115, -0.00041, -0.00041, 0.0, 0.0, 0.0],
                        [0.00120, -0.00044, -0.00044, 0.0, 0.0, 0.0],
                        [0.00125, -0.00046, -0.00046, 0.0, 0.0, 0.0],
                        [0.00130, -0.00048, -0.00048, 0.0, 0.0, 0.0],
                        [0.00135, -0.00051, -0.00051, 0.0, 0.0, 0.0],
                        [0.00140, -0.00053, -0.00053, 0.0, 0.0, 0.0],
                        [0.00145, -0.00055, -0.00055, 0.0, 0.0, 0.0],
                        [0.00150, -0.00058, -0.00058, 0.0, 0.0, 0.0],
                        [0.00155, -0.00060, -0.00060, 0.0, 0.0, 0.0],
                        [0.00160, -0.00062, -0.00062, 0.0, 0.0, 0.0],
                        [0.00165, -0.00065, -0.00065, 0.0, 0.0, 0.0],
                        [0.00170, -0.00067, -0.00067, 0.0, 0.0, 0.0],
                        [0.00175, -0.00070, -0.00070, 0.0, 0.0, 0.0],
                        [0.00180, -0.00072, -0.00072, 0.0, 0.0, 0.0],
                        [0.00185, -0.00074, -0.00074, 0.0, 0.0, 0.0],
                        [0.00190, -0.00077, -0.00077, 0.0, 0.0, 0.0],
                        [0.00195, -0.00079, -0.00079, 0.0, 0.0, 0.0],
                        [0.00200, -0.00081, -0.00081, 0.0, 0.0, 0.0],
                        [0.00107, -0.00053, -0.00053, 0.0, 0.0, 0.0],
                        [0.00107, -0.00053, -0.00053, 0.0, 0.0, 0.0],
                        [0.00112, -0.00055, -0.00055, 0.0, 0.0, 0.0],
                        [0.00117, -0.00056, -0.00056, 0.0, 0.0, 0.0],
                        [0.00122, -0.00058, -0.00058, 0.0, 0.0, 0.0],
                        [0.00127, -0.00059, -0.00059, 0.0, 0.0, 0.0],
                        [0.00132, -0.00061, -0.00061, 0.0, 0.0, 0.0],
                        [0.00137, -0.00062, -0.00062, 0.0, 0.0, 0.0],
                        [0.00142, -0.00064, -0.00064, 0.0, 0.0, 0.0],
                        [0.00147, -0.00065, -0.00065, 0.0, 0.0, 0.0],
                        [0.00152, -0.00067, -0.00067, 0.0, 0.0, 0.0],
                        [0.00157, -0.00068, -0.00068, 0.0, 0.0, 0.0],
                        [0.00162, -0.00070, -0.00070, 0.0, 0.0, 0.0],
                        [0.00167, -0.00071, -0.00071, 0.0, 0.0, 0.0],
                        [0.00172, -0.00073, -0.00073, 0.0, 0.0, 0.0],
                        [0.00177, -0.00074, -0.00074, 0.0, 0.0, 0.0],
                        [0.00182, -0.00076, -0.00076, 0.0, 0.0, 0.0],
                        [0.00187, -0.00077, -0.00077, 0.0, 0.0, 0.0],
                        [0.00192, -0.00079, -0.00079, 0.0, 0.0, 0.0],
                        [0.00197, -0.00081, -0.00081, 0.0, 0.0, 0.0],
                        [0.00202, -0.00082, -0.00082, 0.0, 0.0, 0.0],
                        [0.00207, -0.00084, -0.00084, 0.0, 0.0, 0.0],
                        [0.00212, -0.00086, -0.00086, 0.0, 0.0, 0.0],
                        [0.00217, -0.00088, -0.00088, 0.0, 0.0, 0.0],
                        [0.00222, -0.00089, -0.00089, 0.0, 0.0, 0.0],
                        [0.00227, -0.00092, -0.00092, 0.0, 0.0, 0.0],
                        [0.00232, -0.00094, -0.00094, 0.0, 0.0, 0.0],
                        [0.00237, -0.00098, -0.00098, 0.0, 0.0, 0.0],
                        [0.00242, -0.00101, -0.00101, 0.0, 0.0, 0.0],
                        [0.00247, -0.00103, -0.00103, 0.0, 0.0, 0.0],
                        [0.00252, -0.00106, -0.00106, 0.0, 0.0, 0.0],
                        [0.00257, -0.00109, -0.00109, 0.0, 0.0, 0.0],
                        [0.00262, -0.00111, -0.00111, 0.0, 0.0, 0.0],
                        [0.00267, -0.00113, -0.00113, 0.0, 0.0, 0.0],
                        [0.00272, -0.00116, -0.00116, 0.0, 0.0, 0.0],
                        [0.00277, -0.00118, -0.00118, 0.0, 0.0, 0.0],
                        [0.00282, -0.00121, -0.00121, 0.0, 0.0, 0.0],
                        [0.00287, -0.00123, -0.00123, 0.0, 0.0, 0.0],
                        [0.00292, -0.00126, -0.00126, 0.0, 0.0, 0.0],
                        [0.00297, -0.00128, -0.00128, 0.0, 0.0, 0.0],
                        [0.00302, -0.00131, -0.00131, 0.0, 0.0, 0.0],
                        [0.00307, -0.00133, -0.00133, 0.0, 0.0, 0.0],
                        [0.00312, -0.00135, -0.00135, 0.0, 0.0, 0.0],
                        [0.00317, -0.00138, -0.00138, 0.0, 0.0, 0.0],
                        [0.00322, -0.00140, -0.00140, 0.0, 0.0, 0.0],
                        [0.00327, -0.00143, -0.00143, 0.0, 0.0, 0.0],
                        [0.00332, -0.00145, -0.00145, 0.0, 0.0, 0.0],
                        [0.00337, -0.00148, -0.00148, 0.0, 0.0, 0.0],
                        [0.00342, -0.00150, -0.00150, 0.0, 0.0, 0.0],
                        [0.00347, -0.00152, -0.00152, 0.0, 0.0, 0.0],
                        [0.00352, -0.00155, -0.00155, 0.0, 0.0, 0.0],
                        [0.00357, -0.00157, -0.00157, 0.0, 0.0, 0.0],
                        [0.00362, -0.00160, -0.00160, 0.0, 0.0, 0.0],
                        [0.00367, -0.00162, -0.00162, 0.0, 0.0, 0.0],
                        [0.00372, -0.00165, -0.00165, 0.0, 0.0, 0.0],
                        [0.00377, -0.00167, -0.00167, 0.0, 0.0, 0.0],
                        [0.00382, -0.00170, -0.00170, 0.0, 0.0, 0.0],
                        [0.00387, -0.00172, -0.00172, 0.0, 0.0, 0.0],
                        [0.00392, -0.00174, -0.00174, 0.0, 0.0, 0.0],
                        [0.00397, -0.00177, -0.00177, 0.0, 0.0, 0.0],
                        [0.00402, -0.00179, -0.00179, 0.0, 0.0, 0.0],
                        [0.00407, -0.00182, -0.00182, 0.0, 0.0, 0.0],
                        [0.00412, -0.00184, -0.00184, 0.0, 0.0, 0.0],
                        [0.00417, -0.00187, -0.00187, 0.0, 0.0, 0.0],
                        [0.00422, -0.00189, -0.00189, 0.0, 0.0, 0.0],
                        [0.00427, -0.00192, -0.00192, 0.0, 0.0, 0.0],
                        [0.00432, -0.00194, -0.00194, 0.0, 0.0, 0.0],
                        [0.00437, -0.00197, -0.00197, 0.0, 0.0, 0.0],
                        [0.00442, -0.00199, -0.00199, 0.0, 0.0, 0.0],
                        [0.00447, -0.00201, -0.00201, 0.0, 0.0, 0.0],
                        [0.00452, -0.00204, -0.00204, 0.0, 0.0, 0.0],
                        [0.00457, -0.00206, -0.00206, 0.0, 0.0, 0.0],
                        [0.00462, -0.00209, -0.00209, 0.0, 0.0, 0.0],
                        [0.00467, -0.00211, -0.00211, 0.0, 0.0, 0.0],
                        [0.00472, -0.00214, -0.00214, 0.0, 0.0, 0.0],
                        [0.00477, -0.00216, -0.00216, 0.0, 0.0, 0.0],
                        [0.00482, -0.00219, -0.00219, 0.0, 0.0, 0.0],
                        [0.00487, -0.00221, -0.00221, 0.0, 0.0, 0.0],
                        [0.00492, -0.00224, -0.00224, 0.0, 0.0, 0.0],
                        [0.00497, -0.00226, -0.00226, 0.0, 0.0, 0.0],
                        [0.00502, -0.00229, -0.00229, 0.0, 0.0, 0.0]],

    parameters = DSAParameterState(E =    200.0e3,
                                   nu =      0.3,
                                   R0 =   100.0,
                                   Kn =   100.0,
                                   nn =    10.0,
                                   C1 = 10000.0,
                                   D1 =   100.0,
                                   C2 = 50000.0,
                                   D2 =  1000.0,
                                   Q =     50.0,
                                   b =      0.1,
                                   w =      1e-5,
                                   P1 =   200.0,
                                   P2 =     1e-1,
                                   m =      0.66,
                                   m1 =     6.0,
                                   m2 =     6.0,
                                   M1 =  6000.0,
                                   M2 =  6000.0,
                                   ba =     1e4,
                                   xi =      1.0),
    material    = DSA(parameters=parameters),
    dtime       = 0.25,
    n_steps     = 100,
    n_interrupt = 40,  # for interrupt-and-hold test

    dstrain11   = 2e-4 * dtime, # Corresponds to 10 MPa elastic stress response

    tostrain(tens::Symm2) = copy(tovoigt(tens; offdiagscale=2.0)),
    tostress(tens::Symm2) = copy(tovoigt(tens)),

    times     = [material.drivers.time],
    stresses  = [tostress(material.variables.stress)],
    strains   = [tostrain(material.drivers.strain)],
    Ras       = [copy(material.variables.Ra)],
    tas       = [copy(material.variables.ta)],
    cumeqs    = [copy(material.variables.cumeq)]

    function snapshot!()
        push!(times, material.drivers.time)
        push!(stresses, tostress(material.variables.stress))
        push!(strains, tostrain(material.drivers.strain))
        push!(Ras, material.variables.Ra)
        push!(tas, material.variables.ta)
        push!(cumeqs, copy(material.variables.cumeq))
    end

    # TODO: This doesn't actually test anything.
    # # Uninterrupted test
    # material2 = DSA(parameters = parameters)
    # times2    = [material2.drivers.time]
    # stresses2 = [tostress(material2.variables.stress)]
    # strains2  = [tostrain(material2.drivers.strain)]
    # for i in 1:n_steps
    #     uniaxial_increment!(material2, dstrain11, dtime)
    #     update_material!(material2)
    #     push!(times2, material2.drivers.time)
    #     push!(stresses2, copy(tovoigt(material2.variables.stress)))
    #     push!(strains2, copy(tovoigt(material2.drivers.strain; offdiagscale = 2.0)))
    # end

    # Interrupted test
    for i in 1:n_interrupt
        uniaxial_increment!(material, dstrain11, dtime)
        update_material!(material)
        snapshot!()
    end

    # Interrupt and hold
    # Drive to zero stress
    strain_at_stop = material.drivers.strain[1,1]
    let dstress11 = -material.variables.stress[1,1]
        stress_driven_uniaxial_increment!(material, dstress11, dtime)
    end
    update_material!(material)
    snapshot!()
    # Hold for 3600 seconds
    stress_driven_uniaxial_increment!(material, 0.0, 3600)
    update_material!(material)
    snapshot!()

    # Continue test
    dstrain_extra = strain_at_stop - material.drivers.strain[1,1]
    n_extra_steps = Int(ceil(dstrain_extra / dstrain11))
    for i in (n_interrupt + 1):(n_steps + n_extra_steps)
        uniaxial_increment!(material, dstrain11, dtime)
        update_material!(material)
        snapshot!()
    end

    for i in 1:length(times)
        @test isapprox(stresses[i], stresses_expected[i]; atol = 1e-3)
        @test isapprox(strains[i],  strains_expected[i]; atol = 1e-5)
    end
    dcumeq = cumeqs[end] - cumeqs[end - 1]
    @test dcumeq > 0
end

# Plotting
# using PyPlot
# x11     = [a[1] for a in strains]
# y11     = [a[1] for a in stresses]
# x112    = [a[1] for a in strains2]
# y112    = [a[1] for a in stresses2]
# RasNorm = [Ra / parameters.P1 for Ra in Ras]
# tasNorm = [ta / maximum(tas) for ta in tas]
# fig = figure("test_DSA.jl", figsize = (5, 12)) # Create a new blank figure
# subplot(211)
# plot(x11,y11, label = "interrupted")
# plot(x112,y112,linestyle = "--", label = "uninterrupted")
# title("test_DSA.jl")
# xlabel("Strain, \$\\varepsilon_{11}\$")
# ylabel("Stress, \$\\sigma_{11}\$")
# legend()
# subplot(212)
# plot(times, RasNorm, label = "\$R_a\$")
# plot(times, tasNorm, linestyle = "--", label = "\$t_a\$")
# xlim([3600.0, maximum(times)])
# title("Normalized Evolution of \$R_a\$ & \$t_a\$")
# xlabel("Time")
# ylabel("Ra, ta")
# legend()
# fig.canvas.draw() # Update the figure
# show()
# gcf()


text:
const IV_TERMINALS = [Symbol(".=="), Symbol(".<="), Symbol(".>=")] # Calls that should terminate the tree search
const IV_PARAMETERIZED = Dict(:any => 0, :all => 0, :all_before => 1, :all_after => 1, :all_between => 2, :any_between=>2)

# Negation operator applied to bool or :anybool
flex_not(b::Union{Symbol, Bool}) =  (b == :anybool ? :anybool : !b)

all_before(τ, i) = all(τ[1:i])
all_after(τ, i) = all(τ[i:end])

all_between(τ, i, j) = all(τ[min(i,j):min(length(τ), max(i,j))])
any_between(τ, i, j) = any(τ[min(i,j):min(length(τ), max(i,j))])

function all_between_inv(out, i, j, N, rng::AbstractRNG)
    l = min(i,j)
    l > N && throw(InfeasibleConstraint("Couldn't find feasible constraint due to out of bounds range"))
    h = min(max(i,j), N)
    arr = Array{Any}(undef, N)
    fill!(arr, :anybool)
    if out == true
         arr[l:h] .= true
    elseif out == false
        arr[rand(rng, l:h)] = false
    end
    (arr,)
end

function any_between_inv(out, i, j, N, rng::AbstractRNG)
    l = min(i,j)
    l > N && throw(InfeasibleConstraint("Couldn't find feasible constraint due to out of bounds range"))
    h = min(max(i,j), N)
    arr = Array{Any}(undef, N)
    fill!(arr, :anybool)
    if out == true
        arr[rand(rng, l:h)] = true
    elseif out == false
        arr[l:h] .= false
    end
    (arr,)
end

function all_before_inv(out, i, N, rng::AbstractRNG)
    arr = Array{Any}(undef, N)
    fill!(arr, :anybool)
    if out == true
         arr[1:i] .= true
    elseif out == false
        arr[rand(rng, 1:i)] = false
    end
    (arr,)
end

function all_after_inv(out, i, N, rng::AbstractRNG)
    arr = Array{Any}(undef, N)
    fill!(arr, :anybool)
    if out == true
         arr[i:N] .= true
    elseif out == false
        arr[rand(rng, i:N)] = false
    end
    (arr,)
end

function not_inv(out, rng::AbstractRNG)
    out == :anybool && return :anybool
    (length(out) > 1 ) ? ([flex_not(v) for v in out],) : flex_not(out)
end

# Inverse of the "and" operator which takes two boolean inputs
# If output is true then both inputs are true
# If the output is false then at least one of the inputs is false
function and_inv(out, rng::AbstractRNG)
    if out == :anybool
        return (:anybool, :anybool)
    elseif out
        return (true, true)
    else
        return rand(rng, [(false, :anybool), (:anybool, false)])
    end
end

# Inverse of the "or" operator which takes two boolean inputs
# If the output is true then at least one of the inputs is true
# If the output is false then both if the inputs are false
function or_inv(out, rng::AbstractRNG)
    if out == :anybool
        return (:anybool, :anybool)
    elseif out
        return rand(rng, [(true, :anybool), (:anybool, true)])
    else
        return (false, false)
    end
end

# Inverse of the "all" or "globally" operator which takes a vector of bools
# If the output is true then all inputs in the time series are true
# If the output is false then there is at least on false in the time series
function all_inv(out, N, rng::AbstractRNG)
    arr = Array{Any}(undef, N)
    fill!(arr, :anybool)
    if out == true
         fill!(arr, true)
    elseif out == false
        arr[rand(rng, 1:N)] = false
    end
    (arr,)
end

# Inverse of the "any" or "eventually" operator which takes a vector of bools
# If the output is true then there is at least on true in the time series
# If the output is false then al
function any_inv(out, N, rng::AbstractRNG)
    arr = Array{Any}(undef, N)
    fill!(arr, :anybool)
    if out == true
        arr[rand(rng, 1:N)] = true
    elseif out == false
        fill!(arr, false)
    end
    (arr,)
end

# This function applies the inverse of a scalar boolean operation for each time
# `op_inv` is a function that specifies the scalar boolean operator inverse
# This function assumes that `op` takes in 2 arguments and outputs 1
function bitwise_op_inv(out, op_inv, rng::AbstractRNG)
    N = length(out)
    arr1, arr2 = Array{Any}(undef, N), Array{Any}(undef, N)
    for i in 1:N
        arr1[i], arr2[i] = op_inv(out[i], rng)
    end
    return (arr1, arr2)
end

# Defines the invers of the bitwise "and" operator
bitwise_and_inv(out, rng::AbstractRNG) = (out isa Array) ? bitwise_op_inv(out, and_inv, rng) : and_inv(out, rng)

# Defines the inverse of the bitwise "or" operator
bitwise_or_inv(out, rng::AbstractRNG) = (out isa Array) ? bitwise_op_inv(out, or_inv, rng) : or_inv(out, rng)

# Mapping an operation to its inverse
bool_inverses = Dict(
        :&& => bitwise_and_inv,
        :|| => bitwise_or_inv,
        :any => any_inv,
        :all => all_inv,
        :! => not_inv,
        :all_before => all_before_inv,
        :all_after => all_after_inv,
        :all_between => all_between_inv,
        :any_between => any_between_inv,
        :.& => bitwise_and_inv,
        :.| => bitwise_or_inv
    )

# Get a list of expressions and outcomes that need to be satisfied for the top level expression to be true
function sample_constraints(expr, N, rng::AbstractRNG)
    constraints = []
    sample_constraints!(expr, true, constraints, N, rng)
    constraints
end

# This version passes around a list of constraints (expression, values) pairs that need to be satisfied
# `expr` is an expression to sample a trajectory from
# `truthval` - The desired output of the expression. Could be a scalar bool or a time series of bools
# `constraints` - The list of constraints and their corresponding truth values
# `N` - The length of the time series
function sample_constraints!(expr, truthval, constraints, N, rng::AbstractRNG)
    all(truthval .== :anybool) && return
    if expr.head == :call && expr.args[1] in IV_TERMINALS
        # Add constrain expression to the list of constraints and return
        push!(constraints, [expr, truthval])
        return
    end
    if expr.head == :call
        # This could be :any, :all, :&, :|
        # Get the inverse of these operations and recurse
        # Special handing for expanding operators that need to be passed `N`
        op = expr.args[1]
        inv = bool_inverses[op]
        if op in keys(IV_PARAMETERIZED)
            ps = IV_PARAMETERIZED[op]
            if length(truthval) == 1
                inv_res = inv(truthval, expr.args[3:3+(ps - 1)]..., N, rng)
                sample_constraints!(expr.args[2], inv_res[1], constraints, N, rng)
            else
                for i=1:N
                    inv_res = [fill(:anybool, i-1)..., inv(truthval[i], expr.args[3:3+(ps - 1)]..., N-i+1, rng)[1]...]
                    sample_constraints!(expr.args[2], inv_res, constraints, N, rng)
                end
            end
        else
            inv_res = inv(truthval, rng)
            for i in 1:length(inv_res)
                sample_constraints!(expr.args[i+1], inv_res[i], constraints, N, rng)
            end
        end

    else
        # Here "head" contains the operator and args contains the expressions
        # Get the inverse and recurse directly
        inv_res = bool_inverses[expr.head](truthval, rng)
        for i in 1:length(inv_res)
            sample_constraints!(expr.args[i], inv_res[i], constraints, N, rng)
        end
    end
end

## The stuff below here is for implies clauses
function parse_implications(expr)
    implications = Dict{Expr, Expr}()
    parse_implications!(expr, implications)
    implications
end

function parse_implications!(expr, implications)
    if expr.head == :vect
        implications[expr.args[1]] = expr.args[2]
        return
    elseif expr.head == :&&
        parse_implications!(expr.args[1], implications)
        parse_implications!(expr.args[2], implications)
    else
        throw(error("Unregonized head ", expr.head))
    end
end

function and_expressions(exprs)
    new_expr = exprs[1]
    for i=2:length(exprs)
        new_expr = Expr(:.&, new_expr, exprs[i])
    end
    Expr(:call, :all, new_expr)
end

ext: jl


text:
using BigO
using Documenter

DocMeta.setdocmeta!(BigO, :DocTestSetup, :(using BigO); recursive=true)

makedocs(;
    modules=[BigO],
    authors="Endre Márk Borza",
    repo="https://github.com/endremborza/BigO.jl/blob/{commit}{path}#{line}",
    sitename="BigO.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://endremborza.github.io/BigO.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/endremborza/BigO.jl",
    devbranch="main"
)


text:
__precompile__(false)
module Shearlab

using Images
using Wavelets
using CairoMakie
#using ArrayFire
using DSP
using FFTW

export
			WT, filt_gen, resize_image, load_image, imageplot,
			rescale, plot_wavelet,subsampling, upsampling,
			cconvol, reverse, clamp, perform_wavorth_transform, snr,
			padarray, fliplr, upsample, fix, dshear, checkfiltersizes, Filterconfigs,
			Filterswedgebandlow, getwedgebandpasslowpassfilters2D, Preparedfilters,
			preparefilters2D, getshearletidxs2D, Shearletsystems2D, getshearletsystem2D,
			Shearlets2D, getshearlets2D,
			sheardec2D, shearrec2D,
			Serialpreparedfilters, prepareserial2D, Sheardecserial2D, sheardecserial2D,
			shearrecserial2D, finishserial2D#,
			# AFArray

include("util/operations.jl")
include("2D/getshearlets2D.jl")
include("2D/shearletdecrec2D.jl")
include("2D/shearletserialdecrec2D.jl")
include("fast_wavelet/filters.jl")
include("fast_wavelet/operations_wavelet.jl")
include("fast_wavelet/imaging.jl")

# Alias for functions to have same names as in ShearLab3D
const SLupsample = upsample
const SLdshare = dshear
const SLcheckFilterSizes = checkfiltersizes
const SLWedgeBandpassAndLowpassFilters2D = getwedgebandpasslowpassfilters2D
const SLprepareFilters2D = preparefilters2D
const SLgetShearletIdxs2D = getshearletidxs2D
const SLgetShearletSystem2D = getshearletsystem2D
const SLsheardec2D = sheardec2D
const SLshearrec2D = shearrec2D
const SLprepareSerial2D = prepareserial2D
const SLgetShearlets2D = getshearlets2D
const SLsheardecSerial2D = sheardecserial2D
const SLshearrecSerial2D = shearrecserial2D
const SLfinishSerial2D = finishserial2D

end #module


text:
using Test
using ThreeBodyDecay

@testset "jp_str macro" begin
    @test jp"3+" == str2jp("3+")
    @test jp"3^+" == jp"3+"
    @test jp"1/2+" == jp(1//2, '+')
    @test jp"3/2-" == jp(3//2, '-')
    @test jp"3+" == jp(3, '+')
end

@testset "jp ⊗ jp" begin
    Swave = jp(1//2,'+') ⊗ jp(1,'-')
    Pwave = [sw ⊗ jp(1,'-') for sw in Swave]
    Dwave = [sw ⊗ jp(2,'+') for sw in Swave]
    #
    @test length(Swave) == 2
    @test length(vcat(Pwave...)) == 5
    @test length(Set(vcat(Pwave...))) == 3
end

let
    @test length(possible_ls(jp"3/2-", jp"3-"; jp=jp"1/2+")) == 4
    #
    lsLSv = possible_lsLS(1, jp"1/2+",
        [jp"1+", jp"1/2+", jp"0-", jp"1/2-"])
    @test size(lsLSv) == (1,2)
    #
    lsLSv = possible_lsLS(1, 1, '+',
        ThreeBodySpins(2,1,0; two_h0=1),
        ThreeBodyParities('+', '+', '-'; P0='-'))
    @test size(lsLSv) == (1,2)
end

@testset "x2 over2" begin
    @test x2(1//2) == 1
    @test x2(3) == 6
    @test x2((2,3//2)) == (4,3)
    #
    @test over2(3) == 3//2
    @test over2(2) == 1
    @test over2((2,3)) == (1,3//2)
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
using Plots, Printf

"""
	shallow_water_2D(; nx, ny, dam_x=true, do_visu=true)

2D shallow water equations solver for an instantaneous dam break.
The Lax-Friedrichs Method was applied to the continuity equation.
Geometry (length of 40 meters, width of 20 meters) and initial conditions
(half of domain has initial water level of 5 meters, other half is dry)
matches BASEMENT version 2.8 Test Case H_1 "Dam break in a closed channel."
Solution of momentum equations requires division by H. Therefore a minimum
H must be set throughout the domain, to avoid numerical instability.

# Arguments
    - `nx`: number of discretised cells for x dimension.
    - `ny`: number of discretised cells for y dimension.
    - `dam_x`: if true, 1D dam break in x-direction, else 1D dam break in y-direction.
    - `do_visu`: if true, each physical time step will be ploted.

# Return values
    - `H`: The solution arrray (Water surface height in m).
"""
@views function shallow_water_2D(;
    # Numerics
    nx = 512,
    ny = 256,
    # Modeling
    dam_x = true,
    # Visualisation
    do_visu = true)
    # Physics
    Lx, Ly = 40.0, 20.0
    g = 9.81
    u_max = 5 #from review of results with very small timesteps
    H_init = 5.0
    time = 0
    # Numerics
    nout = 100
    H_min = 1.0
    # Derived numerics
    dx, dy = Lx / nx, Ly / ny
    nt = 10000#cld(ttot, dt)
    xc, yc = LinRange(dx / 2, Lx - dx / 2, nx), LinRange(dy / 2, Ly - dy / 2, ny)
    # Initial Conditions
    H = zeros(Float64, nx, ny) .+ H_min
    if (dam_x)
        H[1:round(Int64(nx / 2)), :] = zeros(round(Int64(nx / 2)), ny) .+ H_init #1D dam break in x-direction
    else
        H[:, 1:round(Int64(ny / 2))] = zeros(Float64, nx, round(Int64(ny / 2))) .+ H_init #1D dam break in y direction
    end
    #H[1:round(Int64(nx/2)), 1:round(Int64(ny/2))] = zeros(Float64, round(Int64(nx/2)), round(Int64(ny/2))) .+ H_init #2D dam break
    #H[round(Int64(3*nx/8)):round(Int64(5*nx/8)), round(Int64(3*ny/8)):round(Int64(5*ny/8))] = zeros(Float64, round(Int64(nx/4))+1, round(Int64(ny/4))+1) .+ H_init #2D water column
    # Array initialisation
    H_new = copy(H)
    u = zeros(Float64, nx, ny)
    u_new = copy(u)
    v = zeros(Float64, nx, ny)
    v_new = copy(v)
    F = zeros(Float64, nx + 1, ny)
    G = zeros(Float64, nx, ny + 1)
    t_tic = 0.0
    niter = 0
    # Time loop
    for it = 0:nt
        if (it == 11)
            t_tic = Base.time()
            niter = 0
        end
        dt = dx / (max(maximum(abs.(u)), maximum(abs.(v)), u_max)) / 4.1
        time = time + dt
        λx = dx / 2 / dt #x-direction signal speed for Lax-Friedrichs scheme
        λy = dy / 2 / dt #y-direction signal speed for Lax-Friedrichs scheme
        #Continuity
        F[2:end-1, :] .= 1 / 2 .* (H[1:end-1, :] .* u[1:end-1, :] .+ H[2:end, :] .* u[2:end, :]) .- λx / 2 .* (H[2:end, :] .- H[1:end-1, :])
        G[:, 2:end-1] .= 1 / 2 .* (H[:, 1:end-1] .* v[:, 1:end-1] .+ H[:, 2:end] .* v[:, 2:end]) .- λy / 2 .* (H[:, 2:end] .- H[:, 1:end-1])
        H_new[2:end-1, :] .= H[2:end-1, :] .- dt / dx .* (F[3:end-1, :] .- F[2:end-2, :]) .- dt / dy .* (G[2:end-1, 2:end] .- G[2:end-1, 1:end-1])
        #Momentum in x-direction
        F[2:end-1, :] .= 1 / 2 .* ((H[1:end-1, :] .* ((u[1:end-1, :]) .^ 2) .+ 0.5 * g .* ((H[1:end-1, :]) .^ 2)) .+ (H[2:end, :] .* ((u[2:end, :]) .^ 2) .+ 0.5 * g .* ((H[2:end, :]) .^ 2)))
        G[:, 2:end-1] .= 1 / 2 .* ((H[:, 1:end-1] .* u[:, 1:end-1] .* v[:, 1:end-1]) .+ (H[:, 2:end] .* u[:, 2:end] .* v[:, 2:end]))
        u_new[2:end-1, :] .= (1 / 2 .* (H[1:end-2, :] .* u[1:end-2, :] .+ H[3:end, :] .* u[3:end, :]) .- dt / dx .* (F[3:end-1, :] .- F[2:end-2, :]) .- dt / dy .* (G[2:end-1, 2:end] .- G[2:end-1, 1:end-1])) ./ H_new[2:end-1, :]
        #Momentum in y-direction
        F[2:end-1, :] .= 1 / 2 .* ((H[1:end-1, :] .* u[1:end-1, :] .* v[1:end-1, :]) .+ (H[2:end, :] .* u[2:end, :] .* v[2:end, :]))
        G[:, 2:end-1] .= 1 / 2 .* ((H[:, 1:end-1] .* ((v[:, 1:end-1]) .^ 2) .+ 0.5 * g .* ((H[:, 1:end-1]) .^ 2)) .+ (H[:, 2:end] .* ((v[:, 2:end]) .^ 2) .+ 0.5 * g .* ((H[:, 2:end]) .^ 2)))
        v_new[:, 2:end-1] .= (1 / 2 .* (H[:, 1:end-2] .* v[:, 1:end-2] .+ H[:, 3:end] .* v[:, 3:end]) .- dt / dx .* (F[2:end, 2:end-1] .- F[1:end-1, 2:end-1]) .- dt / dy .* (G[:, 3:end-1] .- G[:, 2:end-2])) ./ H_new[:, 2:end-1]
        #No-flux boundary conditions
        H_new[1, :] .= H_new[2, :]
        H_new[end, :] .= H_new[end-1, :]
        u_new[1, :] .= 0
        u_new[end, :] .= 0
        v_new[:, 1] .= 0
        v_new[:, end] .= 0

        H .= H_new
        u .= u_new
        v .= v_new

        niter += 1
        if do_visu && (it % nout == 0)
            p1 = plot(xc, H[:, round(Int64(ny / ny))], xlims = (xc[1], xc[end]), ylims = (-3, 8),
                xlabel = "Lx at y=10 (m)", ylabel = "water surface elevation (m)", label = "h",
                title = "time = $(round(time, sigdigits=3)) s", linewidth = :1.0, framestyle = :box)
            plot!(xc, u[:, round(Int64(ny / ny))], label = "u", linewidth = :1.0)
            plot!(xc, v[:, round(Int64(ny / ny))], label = "v", linewidth = :1.0)
            display(p1)
        end
    end
    t_toc = Base.time() - t_tic
    A_eff = (3 * 2) / 1e9 * nx * sizeof(Float64)  # Effective main memory access per iteration [GB]
    t_it = t_toc / niter                      # Execution time per iteration [s]
    T_eff = A_eff / t_it                       # Effective memory throughput [GB/s]
    @printf("Time = %1.3f sec, T_eff = %1.4f GB/s (niter = %d)\n", t_toc, round(T_eff, sigdigits = 3), niter)
    return H
end

# shallow_water_2D(; nx = 256, ny = 128, do_visu = true)


text:
function FundamentalNutationArguments(arg0::IERSConventions, arg1::TimeScale, arg2::InputStream, arg3::JString)
    return FundamentalNutationArguments((IERSConventions, TimeScale, InputStream, JString), arg0, arg1, arg2, arg3)
end

function FundamentalNutationArguments(arg0::IERSConventions, arg1::TimeScale, arg2::InputStream, arg3::JString, arg4::TimeScales)
    return FundamentalNutationArguments((IERSConventions, TimeScale, InputStream, JString, TimeScales), arg0, arg1, arg2, arg3, arg4)
end

function FundamentalNutationArguments(arg0::IERSConventions, arg1::TimeScale, arg2::List)
    return FundamentalNutationArguments((IERSConventions, TimeScale, List), arg0, arg1, arg2)
end

function FundamentalNutationArguments(arg0::IERSConventions, arg1::TimeScale, arg2::List, arg3::TimeScales)
    return FundamentalNutationArguments((IERSConventions, TimeScale, List, TimeScales), arg0, arg1, arg2, arg3)
end

function evaluate_all(obj::FundamentalNutationArguments, arg0::AbsoluteDate)
    return jcall(obj, "evaluateAll", BodiesElements, (AbsoluteDate,), arg0)
end

function evaluate_all(obj::FundamentalNutationArguments, arg0::FieldAbsoluteDate)
    return jcall(obj, "evaluateAll", FieldBodiesElements, (FieldAbsoluteDate,), arg0)
end


text:
@testset "242.valid-anagram.jl" begin
    @test is_anagram("nagaram", "anagram") == true
    @test is_anagram("rat", "car") == false
end


text:
function BFS(G::DictGraph, start::Int)::SortedDict
    bsfDict = SortedDict(i => [-1, typemax(Int), -1] for i in keys(G.relation))
    bsfDict[start] = [0, 0, -1]

    𝑄 = Array{Int, 1}()

    push!(𝑄, start)

    while !isempty(𝑄)
        𝑢 = pop!(𝑄)
        for 𝑣 in G[𝑢]
            if bsfDict[𝑣[1]][1] == -1
                bsfDict[𝑣[1]][1] = 0
                bsfDict[𝑣[1]][2] = bsfDict[𝑢][2] + 1
                bsfDict[𝑣[1]][3] = 𝑢
                push!(𝑄, 𝑣.u)
            end
        end
        bsfDict[𝑢][1] = 1
    end

    return bsfDict
end

function BFSv2(G::DictGraph, start::Int)::SortedDict
    bsfDict = SortedDict(i => [-1, typemax(Int), -1] for i in keys(G.relation))
    bsfDict[start] = [0, 0, -1]

    𝑄 = Array{Int, 1}()

    push!(𝑄, start)

    while !isempty(𝑄)                                       #while the stack is not empty
        𝑢 = pop!(𝑄)

        for 𝑣 in G[𝑢]                              #for each vertex v that belongs to G.Adj[u]
            if bsfDict[𝑣[1]][2] == typemax(Int)
                bsfDict[𝑣[1]][1] = 0
                bsfDict[𝑣[1]][2] = bsfDict[𝑢][2] + 1
                bsfDict[𝑣[1]][3] = 𝑢
                push!(𝑄, 𝑣[1])
            end
        end

        bsfDict[𝑢][1] = 1
    end
end


text:
using ClearSky
using Test

using ClearSky: MOLPARAM

@testset "Molparam" begin include("test_molparam.jl") end
#@testset "Gray" begin include("test_gray.jl") end


text:
tupvec(vec) = Tuple([entry] for entry in vec)
coords_to_axes(src, dst) = [src[1], dst[1]], [src[2], dst[2]]

function calculate_guide_pos(axis, arrow_pos, dist_from_graph)
    if axis == "x"
        base_src_tern = (arrow_pos[1], arrow_pos[2], 0)
        base_dst_tern = (arrow_pos[2], arrow_pos[1], 0)
        oriented_dfg = [0, -1] .* dist_from_graph
    elseif axis == "y"
        base_src_tern = (0, arrow_pos[1], arrow_pos[2])
        base_dst_tern = (0, arrow_pos[2], arrow_pos[1])
        oriented_dfg = [1, 1] .* dist_from_graph
    elseif axis == "z"
        base_src_tern = (arrow_pos[2], 0, arrow_pos[1])
        base_dst_tern = (arrow_pos[1], 0, arrow_pos[2])
        oriented_dfg = [-1, 1] .* dist_from_graph
    else
        throw(ArgumentError("Invalid axis (\"x, y, z\""))
    end
    arrow_src = collect(tern2cart(base_src_tern)) .+ 2 .* oriented_dfg
    arrow_dst = collect(tern2cart(base_dst_tern)) .+ 2 .* oriented_dfg
    textpos = ((arrow_src .+ arrow_dst) ./ 2) .+ oriented_dfg
    return (arrow_src, arrow_dst, textpos)
end

function calculate_ticks(axis, ticks, tick_length, tickmult)
    tick_endpoints = Vector{Vector{Float64}}(undef, length(ticks))
    tick_startpoints = Vector{Tuple{Float64, Float64}}(undef, length(ticks))

    function start_structure(axis, tick_pos)
        if axis == "x"
            return 1 - tick_pos, tick_pos, 0
        elseif axis == "y"
            return 1 - tick_pos, 0, tick_pos
        elseif axis == "z"
            return 0, tick_pos, 1 - tick_pos
        else
            throw(ArgumentError("Invalid axis (\"x, y, z\""))
        end
    end

    function end_structure(axis, tick_pos)
        if axis == "x"
            return 1 - tick_pos, 0, tick_pos
        elseif axis == "y"
            return 0, 1 - tick_pos, tick_pos
        elseif axis == "z"
            return 1 - tick_pos, tick_pos, 0
        else
            throw(ArgumentError("Invalid axis (\"x, y, z\""))
        end
    end

    for (i, tick_pos) in enumerate(ticks)
        start_point = tern2cart(start_structure(axis, tick_pos)...)
        end_point = tern2cart(end_structure(axis, tick_pos)...)
        tick_offset = normalize(collect(start_point .- end_point)) ./ (1 / tick_length)
        tick_startpoints[i] = start_point
        tick_endpoints[i] = start_point .+ tickmult .* tick_offset
    end

    return tick_startpoints, tick_endpoints
end


text:
include("Potts.jl")
using IterTools:product

function pth_order_tensor(r)
    a =  Tuple([2 for i in 1:r])
    M =  zeros(a)
    for i in CartesianIndices(M)
        s = prod(2*[Tuple(i)...] .- 3)
        M[i] = s
    end
    return M
end

n =15
q = 2

Terms = Dict{}()

M =  pth_order_tensor(2)

for i in 1:n-1
    Terms[(i,i+1)] = (0.3 + rand())*M

end

Terms[(1,n)] = (0.3 + rand())*M

Terms[(1,3,5,7,9)] = 0.5*pth_order_tensor(5)

F = FactorGraph([2,5],n,q,Terms);

n_samples =10^6

s = raw_sampler_potts(F, n_samples, false)
samples=[2*a .- 3  for a in s]
sample_path ="saved_model/samples.csv"
writedlm( sample_path,hcat(samples...))


text:
using Plots
pyplot()

include("../heavy_chain.jl")
include("../heavy_chain_plot.jl")

#some parameters
xmin = -1
xmax = 1
dt = 1/1000
all_t = 0:dt:40
all_u = [NaN for i = 1:length(all_t)]
x_coords = collect(linspace(0, 1.2, 200))

#create chain
chain = Chain(x_coords)

#define input u(t)
#u_in(t) = mod(t, 3)>1.0 ? 1.0:-1.0

#simulate and animate
anim = Animation()
iter=1
u_in = 1.0
for t = all_t
  #take frame
  if mod(t, 0.04)<dt*0.5 || mod(t, 0.04)>(0.04-dt*0.5)
    println("Frame at t = $(t)")
    l = Plots.@layout [a{0.7h}; b]
    fig = plot(layout = l)
    heavy_chain_plot!(fig[1], chain, xmin, xmax, t)
    plot!(fig[2], all_t, all_u, lab = "u")
    frame(anim)
  end

  #simulate one step
  if chain.w[1]>0.2
    u_in = -(chain.w[1]+0.7)^4
  elseif chain.w[1]<-0.2
    u_in = (chain.w[1]-0.7)^4
  end
  simulate!(chain, u_in, dt)
  all_u[iter] = u_in
  iter+=1
end
gif(anim, "animation.mp4", fps = round(Int, 1/0.04))

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
using BCTRNN
using DiffEqSensitivity
using OrdinaryDiffEq
import DiffEqFlux: FastChain, FastDense
import Flux: ClipValue, ADAM

# Not in Project.toml
using Plots
gr()

include("sine_wave_dataloader.jl")

function train_sine_fc(epochs, solver=Tsit5();
  sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)),
  T=Float32, model_size=5,
  mtkize=false, gen_jac=false, lr=0.02, kwargs...)

  train_dl = generate_2d_data(T)

  f_in = 2
  f_out = 1

  im = BCTRNN.InputAllToAll()
  sm = BCTRNN.SynsFullyConnected()
  om = BCTRNN.OutputAll()
  wiring = BCTRNN.WiringConfig(f_in, model_size, im,sm,om)

  model = FastChain(BCTRNN.Mapper(f_in),
                    BCTRNN.LTC(wiring, solver, sensealg; T, mtkize, gen_jac, kwargs...),
                    FastDense(wiring.n_out, f_out))

  hs = []
  for (k,v) in wiring.matrices
    push!(hs, heatmap(v, title=k))
  end
  display(plot(hs..., layout=length(hs)))

  cb = BCTRNN.MyCallback(T; ecb=mycb, nepochs=epochs, nsamples=length(train_dl))
  #opt = GalacticOptim.Flux.Optimiser(ClipValue(0.5), ADAM(0.02))
  opt = BCTRNN.ClampBoundOptim(BCTRNN.get_bounds(model,T)..., ClipValue(T(1.0)), ADAM(T(lr)))
  BCTRNN.optimize(model, BCTRNN.loss_seq, cb, opt, train_dl, epochs, T), model
end

@time train_sine_fc(100, Tsit5(); model_size=10, mtkize=false) # 29.983295 seconds (134.07 M allocations: 6.758 GiB, 5.35% gc time)

@time train_sine_fc(100, Tsit5(); model_size=10, mtkize=false, sensealg=InterpolatingAdjoint(autojacvec=ZygoteVJP()))
@time train_sine_fc(100, Tsit5(); model_size=10, mtkize=false, sensealg=InterpolatingAdjoint(checkpointing=true, autojacvec=ZygoteVJP()))

@time train_sine_fc(100, Tsit5(); model_size=10, mtkize=true, gen_jac=true, sensealg=InterpolatingAdjoint())
@time train_sine_fc(100, Tsit5(); model_size=10, mtkize=true, gen_jac=true, sensealg=InterpolatingAdjoint(autojacvec=ZygoteVJP()))
@time train_sine_fc(100, Tsit5(); model_size=10, mtkize=true, gen_jac=true, sensealg=InterpolatingAdjoint(autojacvec=false))

@time train_sine_fc(100, Tsit5(); model_size=10, mtkize=false, sensealg=BacksolveAdjoint())
@time train_sine_fc(100, Tsit5(); model_size=10, mtkize=false, sensealg=BacksolveAdjoint(autojacvec=ZygoteVJP())) # slow
@time train_sine_fc(100, Tsit5(); model_size=10, mtkize=false, sensealg=BacksolveAdjoint(checkpointing=true, autojacvec=ZygoteVJP())) # slow
@time train_sine_fc(100, Tsit5(); model_size=10, mtkize=false, sensealg=BacksolveAdjoint(checkpointing=true, autojacvec=false)) # slow
@time train_sine_fc(100, Tsit5(); model_size=10, mtkize=true, gen_jac=true, sensealg=BacksolveAdjoint()) # slow
#@time train_sine_fc(100, Tsit5(); model_size=10, mtkize=false, sensealg=BacksolveAdjoint(autojacvec=EnzymeVJP())) # fails crashes

@time train_sine_fc(100, Tsit5(); model_size=10, mtkize=false, sensealg=QuadratureAdjoint())
@time train_sine_fc(100, Tsit5(); model_size=10, mtkize=false, sensealg=QuadratureAdjoint(autojacvec=ZygoteVJP()))
@time train_sine_fc(100, Tsit5(); model_size=10, mtkize=false, sensealg=QuadratureAdjoint(compile=true, autojacvec=ZygoteVJP()))
@time train_sine_fc(100, Tsit5(); model_size=10, mtkize=true, gen_jac=true, sensealg=QuadratureAdjoint(compile=true, autojacvec=false))
#@time train_sine_fc(100, Tsit5(); model_size=10, mtkize=false, sensealg=QuadratureAdjoint(autojacvec=EnzymeVJP()))

function train_sine(epochs, solver=Tsit5();
  sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)),
  T=Float32,
  n_sensory=2, n_inter=5, n_command=5, n_motor=1,
  sensory_inter=2, inter_command=3, command_command=2, command_motor=2,
  mtkize=false, lr=0.02, kwargs...)

  train_dl = generate_2d_data(T)

  f_in = 2
  f_out = 1

  im = BCTRNN.InputAllToFirstN(n_sensory)
  sm = BCTRNN.SynsNCP(; n_sensory, n_inter, n_command, n_motor,
    sensory_inter, inter_command, command_command, command_motor
  )
  om = BCTRNN.OutputAll()
  wiring = BCTRNN.WiringConfig(f_in, im,sm,om)

  model = FastChain(BCTRNN.Mapper(f_in),
                    BCTRNN.LTC(wiring, solver, sensealg; T, mtkize, kwargs...),
                    FastDense(wiring.n_out, f_out))

  hs = []
  for (k,v) in wiring.matrices
    push!(hs, heatmap(v, title=k))
  end
  display(plot(hs..., layout=length(hs)))

  cb = BCTRNN.MyCallback(T; ecb=mycb, nepochs=epochs, nsamples=length(train_dl))
  #opt = GalacticOptim.Flux.Optimiser(ClipValue(0.5), ADAM(0.02))
  opt = BCTRNN.ClampBoundOptim(BCTRNN.get_bounds(model,T)..., ClipValue(T(1.0)), ADAM(T(lr)))
  BCTRNN.optimize(model, BCTRNN.loss_seq, cb, opt, train_dl, epochs, T), model
end

#  60.991400 seconds (168.50 M allocations: 19.765 GiB, 5.07% gc time)
@time train_sine(100, Tsit5(), mtkize=false) #  36.071273 seconds (190.71 M allocations: 9.256 GiB, 5.49% gc time)
@time train_sine(100, Tsit5(), n_sensory=3, n_motor=4, mtkize=false)
@time train_sine(200, Tsit5(); model_size=3, abstol=1e-4, reltol=1e-3)

@time train_sine(100, AutoTsit5(Rosenbrock23(autodiff=false)); model_size=10)
@time train_sine(20, AutoTsit5(Rosenbrock23(autodiff=false)); model_size=30)

#@time train_sine(20, SOSRA(); model_size=20)
#@time train_sine(200, SRA3(); model_size=7)


text:
const expected_errs =
    ["SBML files with rules are not supported",
    "Model contains no reactions.",
    "are not yet implemented.",
    "Please make reaction irreversible or rearrange kineticLaw to the form `term1 - term2`."]

function lower_one(fn; verbose=false)
    k = 0
    n_dvs = 0
    n_ps = 0
    time = 0.0
    diffeq_retcode = :nothing
    expected_err = false
    err = ""
    try
        SBMLToolkit.checksupport(fn)
        ml = SBML.readSBML(fn, doc -> begin
            set_level_and_version(3, 2)(doc)
            convert_simplify_math(doc)
        end)
        k = 1
        rs = ReactionSystem(ml)
        k = 2
        sys = ODESystem(ml)
        n_dvs = length(states(sys))
        n_ps = length(parameters(sys))
        k = 3
        prob  = ODEProblem(sys, Pair[], (0, 1.))
        k = 4
        sol = solve(prob, TRBDF2(), dtmax=0.5; force_dtmin=false, unstable_check=unstable_check = (dt, u, p, t) -> any(isnan, u))
        diffeq_retcode = sol.retcode
        if diffeq_retcode == :Success
            k = 5
            time = @belapsed solve($prob, Rosenbrock23())
        end
    catch e
        err = string(e)
        if sum([occursin(e, err) for e in expected_errs]) > 0
            expected_err = true
        end
        if length(err) > 1000 # cutoff since I got ArgumentError: row size (9088174) too large
            err = err[1:1000]
        end
    finally
        verbose && @info("$(basename(fn)) done with a code $k and error msg: $err")
        return (fn, k, n_dvs, n_ps, time, err, expected_err, diffeq_retcode)
    end
end

function lower_fns(fns; verbose=false, write_fn=nothing)
    df = DataFrame(file=String[], retcode=Int[], n_dvs=Int[], n_ps=Int[], time=Float64[], error=String[], expected_error=Bool[], diffeq_retcode=Symbol[])
    for fn in fns
        row = lower_one(fn; verbose=verbose)
        push!(df, row)
    end
    write_fn !== nothing && CSV.write(joinpath(logdir, write_fn), df)
    df
end


text:
#Wave equation Differential Equation
function square_grid_KG!(df,f,param,t)

    #Get Operators and variables
    Opps        = param[1]
    vars        = param[2]

    av  = Opps.aux_vex

    RHS     = Opps.DiffO_1
    BC      = Opps.BC_Mat
    Alpha   = Opps.Abs_Mat
    Lz      = Opps.Lz

    #Fields of our problem (Julia changes the rhs if you change the lhs)
    ψ   = f.x[1]
    ϕ   = f.x[2]

    dψ  = df.x[1]
    dϕ  = df.x[2]

    #Differential Equation
    if (vars.dissipation)

        #If Rotation is on,
        if( vars.Ω != 0.0 )

            #Update absorption matrix
            alpha!(Alpha , vars , Opps , t)

            #First Derivative
            dψ .=  ϕ

            #Second Derivative
            mul!(dϕ , RHS , ψ )
            av =  Alpha * (ϕ  + Lz * ψ)
            dϕ .=  dϕ .- av

        else

            # No need to get absorption matrix
            #alpha!( Alpha , vars , Opps , t)

            #First Derivative
            dψ .=  ϕ

            #Second derivative
            mul!(dϕ , RHS , ψ )
            av =  Alpha * (ϕ  + Lz * ψ)
            dϕ .=  dϕ .- av

        end

    else

        #First Derivative
        dψ .=  ϕ

        #Second Derivative
        mul!(dϕ , RHS , ψ )

    end

    #Boundary Conditions
    if(vars.Boundaries != :periodic )

        mul!(av , BC , dϕ )
       dϕ .= av

    end

end

function solve_wave_equation_2D(p_in::Param)

    #
    # Setup Routine
    #

    #Welcome message
    welcome_message();

    #Start time
    init_time = time()
    max_runtime = time_string_to_seconds(p_in.max_runtime)

    #Setup routine setup check
    if( !setup_folder(p_in) )
        print_exit_message(  time() - init_time )
        return false
    end

    #Check if parameters are OK
    println("Checking Parameters...")
    p = p_in

    if( p.resume_simulation == :true)
        p = get_parameters_from_file(p.input_file_name)

        if( p_in.t_sim_final <= p.t_sim_final )
            println("ERROR: Simulation already at final time")
            print_exit_message(  time() - init_time )
            return false
        end

        p = Param(  p,

                    resume_simulation = p_in.resume_simulation,
                    input_file_name = p_in.input_file_name,

                    t_sim_init  = p.t_sim_final,
                    t_sim_final = p_in.t_sim_final,

                    folder = p_in.folder,
                    fname = p_in.fname
            )
    end

    print_parameters_to_screen(p)

    println("OK")

    #Storage file
    println("Creating HDF5 file...")

    fdata   = "./data"
    fname   = fdata * p.folder *  "$(p.fname).h5"
    fid = h5open(fname, "w")

    save_parameters_to_file( p , fid )

    println("OK")

    #Create operators structures
    println("Setting up Operators...")
    Operators = set_up_operator(p)
    println("OK")

    #Time variables
    tmin = p.t_sim_init
    tmax = p.t_sim_final
    tspan       = ( tmin , tmax )

    iter        = 0
    iter_save   = 0
    every       = floor(Int , p.out_every_t/p.deltat + 1.0)
    save_time   = 0

    #Initial Conditions
    println("Setting up inicial configuration...")
    ψ , dψ = get_gaussian_pulse( p , Operators)

    #If resume simulation
    if( p.resume_simulation == :true)

        #Get max iteration
        _ , _ , itermax = get_time_variables( p.input_file_name )

        #Get field at previous simulation
        _,_,_, ψmat , dψmat = get_fields( p.input_file_name , itermax)

        #Reshape to one column vector
        ψ   = reshape( ψmat , p.xnodes*p.ynodes)
        dψ  = reshape( dψmat , p.xnodes*p.ynodes)
    end

    U = ArrayPartition( Operators.BC_Mat * ψ , Operators.BC_Mat *dψ )
    println("OK")

    #Parameters
    my_params = ( Operators , p)

    #Problem ODE and integrator
    prob = ODEProblem(  square_grid_KG! , U , tspan , my_params )

    integrator = init( prob , RK4() , save_everystep=false , dt=p.deltat , adaptive=false )

    #Create a group for the functions
    g1 =  create_group(fid , "ψ");
    g2 =  create_group(fid , "dψ");

    #Save first configuration of field

    dset_name = string("iteration_",iter_save)

    data1 = integrator.u.x[1]
    dset1, dtype1  = create_dataset(g1, dset_name , data1 )
    write_dataset(dset1, dtype1 , data1)

    data2 = integrator.u.x[2]
    dset2, dtype2  = create_dataset(g2, dset_name , data2 )
    write_dataset(dset2, dtype2 , data2)

    #Add HDF5.attributes time to both datasets
    HDF5.attributes(dset1)["time"] =  integrator.t
    HDF5.attributes(dset2)["time"] =  integrator.t

    #Progress bar and time evolution
    println("Starting up simulation...")
    prog1 = Progress( Int(floor(p.t_sim_final-p.t_sim_init))*100 )

    for (u,t) in tuples(integrator)

        #Increase iteration counter
        iter += 1

        #Step integrator
        DifferentialEquations.step!(integrator)

        #Update bar
        ProgressMeter.update!(prog1, Int(floor(integrator.t*100)))

        #If save iteration
        if iter % every == 0

            iter_save += 1
            save_time = integrator.t

            #data set names
            dset_name = string("iteration_",iter_save)

            #data set for both displacement and derivative
            data1 = integrator.u.x[1]
            dset1, dtype1  = create_dataset(g1, dset_name , data1 )
            write_dataset(dset1, dtype1 , data1)

            data2 = integrator.u.x[2]
            dset2, dtype2  = create_dataset(g2, dset_name , data2 )
            write_dataset(dset2, dtype2 , data2)

            #Add HDF5.attributes time to both datasets
            HDF5.attributes(dset1)["time"] =  integrator.t
            HDF5.attributes(dset2)["time"] =  integrator.t

        end

        #Check if time has exceeded
        if( time_since(init_time) > max_runtime )
            println("\nMaximum allowed runtime reached , Ending Simulation!")
            break
        end

    end

    #Add time information to file
    HDF5.attributes(fid)["max_iter"]    = iter_save
    HDF5.attributes(fid)["t_min"]       = p.t_sim_init
    HDF5.attributes(fid)["t_max"]       = save_time
    HDF5.attributes(fid)["out_every_t"] = p.out_every_t
    HDF5.attributes(fid)["dt"]          = p.deltat

    close(fid)
    println("Saving Files...")

    print_exit_message(  time() - init_time )

    return true
end


text:
function test_contexts()
    lims = [0.0 1.0; 99.0 100.0]
    cdist = UniformContext(lims)
    X = sample(1, cdist)
    @test lims[1, 1] < X[1, 1] < lims[1, 2]
    @test lims[2, 1] < X[1, 2] < lims[2, 2]
end

test_contexts()


text:
include("trace.jl")

"""
    DynamicDSLFunction{T} <: GenerativeFunction{T,DynamicDSLTrace}

A generative function based on a shallowly embedding modeling language based on Julia functions.

Constructed using the `@gen` keyword.
Most methods in the generative function interface involve a end-to-end execution of the function.
"""
struct DynamicDSLFunction{T} <: GenerativeFunction{T,DynamicDSLTrace}
    params_grad::Dict{Symbol,Any}
    params::Dict{Symbol,Any}
    arg_types::Vector{Type}
    has_defaults::Bool
    arg_defaults::Vector{Union{Some{Any},Nothing}}
    julia_function::Function
    has_argument_grads::Vector{Bool}
    accepts_output_grad::Bool
end

function DynamicDSLFunction(arg_types::Vector{Type},
                     arg_defaults::Vector{Union{Some{Any},Nothing}},
                     julia_function::Function,
                     has_argument_grads, ::Type{T},
                     accepts_output_grad::Bool) where {T}
    params_grad = Dict{Symbol,Any}()
    params = Dict{Symbol,Any}()
    has_defaults = any(arg -> arg != nothing, arg_defaults)
    DynamicDSLFunction{T}(params_grad, params, arg_types,
                has_defaults, arg_defaults,
                julia_function,
                has_argument_grads, accepts_output_grad)
end

function DynamicDSLTrace(gen_fn::T, args) where {T<:DynamicDSLFunction}
    # pad args with default values, if available
    if gen_fn.has_defaults && length(args) < length(gen_fn.arg_defaults)
        defaults = gen_fn.arg_defaults[length(args)+1:end]
        defaults = map(x -> something(x), defaults)
        args = Tuple(vcat(collect(args), defaults))
    end
    DynamicDSLTrace{T}(gen_fn, args)
end

accepts_output_grad(gen_fn::DynamicDSLFunction) = gen_fn.accepts_output_grad

mutable struct GFUntracedState
    params::Dict{Symbol,Any}
end

function (gen_fn::DynamicDSLFunction)(args...)
    state = GFUntracedState(gen_fn.params)
    gen_fn.julia_function(state, args...)
end

function exec(gen_fn::DynamicDSLFunction, state, args::Tuple)
    gen_fn.julia_function(state, args...)
end

# whether there is a gradient of score with respect to each argument
# it returns 'nothing' for those arguemnts that don't have a derivatice
has_argument_grads(gen::DynamicDSLFunction) = gen.has_argument_grads

"Global reference to the GFI state for the dynamic modeling language."
const state = gensym("state")

"Implementation of @trace for the dynamic modeling language."
function dynamic_trace_impl(expr::Expr)
    @assert expr.head == :gentrace "Not a Gen trace expression."
    call, addr = expr.args[1], expr.args[2]
    if (call.head != :call) error("syntax error in @trace at $(call)") end
    fn = call.args[1]
    args = Expr(:tuple, call.args[2:end]...)
    if addr != nothing
        addr = something(addr)
        return Expr(:call, GlobalRef(@__MODULE__, :traceat), state, fn, args, addr)
    else
        return Expr(:call, GlobalRef(@__MODULE__, :splice), state, fn, args)
    end
end

# Defaults for untraced execution
@inline traceat(state::GFUntracedState, gen_fn::GenerativeFunction, args, key) =
    gen_fn(args...)

@inline traceat(state::GFUntracedState, dist::Distribution, args, key) =
    random(dist, args...)

@inline splice(state::GFUntracedState, gen_fn::DynamicDSLFunction, args::Tuple) =
    gen_fn(args...)

########################
# trainable parameters #
########################

"Implementation of @param for the dynamic modeling language."
function dynamic_param_impl(expr::Expr)
    @assert expr.head == :genparam "Not a Gen param expression."
    name = expr.args[1]
    Expr(:(=), name, Expr(:call, GlobalRef(@__MODULE__, :read_param), state, QuoteNode(name)))
end

function read_param(state, name::Symbol)
    if haskey(state.params, name)
        state.params[name]
    else
        throw(UndefVarError(name))
    end
end

##################
# AddressVisitor #
##################

struct AddressVisitor
    visited::DynamicSelection
end

AddressVisitor() = AddressVisitor(DynamicSelection())

function visit!(visitor::AddressVisitor, addr)
    if addr in visitor.visited
        error("Attempted to visit address $addr, but it was already visited")
    end
    push!(visitor.visited, addr)
end

function all_visited(visited::Selection, choices::ChoiceMap)
    allvisited = true
    for (key, _) in get_values_shallow(choices)
        allvisited = allvisited && (key in visited)
    end
    for (key, submap) in get_submaps_shallow(choices)
        if !(key in visited)
            subvisited = visited[key]
            allvisited = allvisited && all_visited(subvisited, submap)
        end
    end
    allvisited
end

function get_unvisited(visited::Selection, choices::ChoiceMap)
    unvisited = choicemap()
    for (key, _) in get_values_shallow(choices)
        if !(key in visited)
            set_value!(unvisited, key, get_value(choices, key))
        end
    end
    for (key, submap) in get_submaps_shallow(choices)
        if !(key in visited)
            subvisited = visited[key]
            sub_unvisited = get_unvisited(subvisited, submap)
            set_submap!(unvisited, key, sub_unvisited)
        end
    end
    unvisited
end

get_visited(visitor) = visitor.visited

function check_no_submap(constraints::ChoiceMap, addr)
    if !isempty(get_submap(constraints, addr))
        error("Expected a value at address $addr but found a sub-assignment")
    end
end

function check_no_value(constraints::ChoiceMap, addr)
    if has_value(constraints, addr)
        error("Expected a sub-assignment at address $addr but found a value")
    end
end

function gen_fn_changed_error(addr)
    error("Generative function changed at address: $addr")
end

include("simulate.jl")
include("generate.jl")
include("propose.jl")
include("assess.jl")
include("project.jl")
include("update.jl")
include("regenerate.jl")
include("backprop.jl")

export DynamicDSLFunction

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
"""
Test multi-processing
"""
# Load library to main worker
@time using Distributed, RetroSignalModel, CSV
import RetroSignalModel as rs
# Add workers
addprocs(1, exeflags="--project=$(Base.active_project())");
@show nprocs();

# Load library to workers
@time @everywhere using RetroSignalModel

# Parameter searching
@time rs.paramSearching(rs.rtgM4(); num_sim=3, distributed=true, saveall=true)


text:
using MolecularIntegrals
using BenchmarkTools

function time_ethane()
    bfs = build_basis(ethane,"6-31G")
    @time MolecularIntegrals.all_twoe_ints(bfs);
    #@time MolecularIntegrals.all_twoe_ints(bfs,MolecularIntegrals.coulomb_hgp);
end

function vrr_timings()
    x=y=z=0.0
    xyz = [x,y,z]
    xyza = xyz + [0.1,0.05,0.025]
    xa,ya,za = xyza
    ex = 1.0
    A = B = xyz
    ax,ay,az = bx,by,bz = xyz
    C = D = xyza
    cx,cy,cz = dx,dy,dz = xyza
    amax = cmax = 4

    #print("# vrr_widearray($amax,$cmax) ")
    #@btime MolecularIntegrals.vrr_widearray($amax,$cmax, $ex,$ex,$ex,$ex, $xyz,$xyz,$xyza,$xyza);
    vrrs = zeros(Float64,amax+cmax+1,nao[amax],nao[cmax])
    print("vrr!($amax,$cmax)      ")
    @btime MolecularIntegrals.vrr!($vrrs, $amax,$cmax, $ex,$ex,$ex,$ex, $xyz,$xyz,$xyza,$xyza);
    return nothing
end

# vrr2(4,4)            17.560 ms (166510 allocations: 6.46 MiB)
# vrr_widearray(4,4)   3.572 ms (104351 allocations: 3.85 MiB)
# vrr(4,4)        3.864 ms (81789 allocations: 1.70 MiB)

function hand_vrr_timings()
    x=y=z=0.0
    xyz = [x,y,z]
    xyza = xyz + [0.1,0.05,0.025]
    xa,ya,za = xyza
    ex = 1.0
    A = B = xyz
    ax,ay,az = bx,by,bz = xyz
    C = D = xyza
    cx,cy,cz = dx,dy,dz = xyza
    amax = cmax = 1
    #print("# vrr_widearray($amax,$cmax) ")
    #@btime MolecularIntegrals.vrr_widearray($amax,$cmax, $ex,$ex,$ex,$ex, $xyz,$xyz,$xyza,$xyza);
    print("# vrr($amax,$cmax)      ")
    @btime MolecularIntegrals.vrr($amax,$cmax, $ex,$ex,$ex,$ex, $xyz,$xyz,$xyza,$xyza);
    print("# vrr_pp               ")
    @btime MolecularIntegrals.vrr_pp($ex,$ex,$ex,$ex, $xyz,$xyz,$xyza,$xyza);
end

# vrr_widearray(1,1)   23.402 μs (519 allocations: 20.45 KiB)
# vrr5(1,1)             33.262 μs (384 allocations: 15.94 KiB)
# vrr_pp                 1.586 μs (16 allocations: 2.66 KiB)

function hrr_timings()
    x=y=z=0.0
    xyz = [x,y,z]
    xyza = xyz + [0.1,0.05,0.025]
    xa,ya,za = xyza
    ex = 1
    A = B = xyz
    ax,ay,az = bx,by,bz = xyz
    C = D = xyza
    cx,cy,cz = dx,dy,dz = xyza
    aexpn=bexpn=cexpn=dexpn = ex
    ashell,bshell,cshell,dshell = (2,2,2,2)
    print("# hrr($ashell,$bshell,$cshell,$dshell) ")
    @btime MolecularIntegrals.hrr($ashell,$bshell,$cshell,$dshell, $aexpn,$bexpn,$cexpn,$dexpn, $A,$B,$C,$D);
    #print("# hrr_dict($ashell,$bshell,$cshell,$dshell)  ")
    #@btime MolecularIntegrals.hrr_dict($ashell,$bshell,$cshell,$dshell, $aexpn,$bexpn,$cexpn,$dexpn, $A,$B,$C,$D);
    nothing
end

# hrr2(2,2,2,2)   28.719 ms (282814 allocations: 12.60 MiB)
# hrr1(2,2,2,2)   20.398 ms (289908 allocations: 8.86 MiB)
# hrr5(2,2,2,2)   19.829 ms (257436 allocations: 6.44 MiB)

function ethane_timing()
    for bname in ["sto3g","6-31G","cc-pvdz"]
        bfs = build_basis(ethane,bname)
        println("Ethane/$bname: nbf=$(length(bfs))")
        fetcher = MolecularIntegrals.eri_fetcher(bfs)
        #print("Huzinaga ")
        #@btime MolecularIntegrals.all_twoe_ints($bfs);
        print("HGP      ")
        @btime MolecularIntegrals.all_twoe_ints_chrr($bfs);
        print("Rys      ")
        @btime MolecularIntegrals.all_twoe_ints_rys($bfs);
    end
    return nothing
end

#vrr_timings()
#hand_vrr_timings()
#hrr_timings()
ethane_timing()


text:
module MyLearn2ClassifyMany

using Random
# Activation functions (not needed if MyLearn2Classify is already imported)
linear(z) = z
dlinear(z) = 1.0
dtanh(z) = 1 - tanh(z)^2
sigmoid(z::Real) = 1.0 / (1.0 + exp(-z))
function dsigmoid(z::Real)
    sigmoid_z = 1.0 / (1.0 + exp(-z))
    return sigmoid_z * (1 - sigmoid_z)
end

"""
    Y = g(X, W, b, f_a)

Inputs:
* `X` is a d x N Array{Float64,2} -- even when n or d equal 1
* `W` is an n x d Array{Float64,2}
* `b` is n element Array{Float64,1}

Output: `Y`, the 1 x N matrix of outputs
"""
function g(X::AbstractArray, W::AbstractMatrix, b::AbstractVector, f_a::Function)
    return sum(f_a.(W * X .+ b); dims=1) # type Array{Float64,2}
end

"""
    dW, db, loss = grad_loss_1layer(f_a, df_a, x, y, W, b)

Inputs:
* `x` is a d x N Array{Float64,2} -- even when d = 1
* `y` is a n x N  Array{Float64,2}
* `W` is an n x d Array{Float64,2}
* `b` is n element Array{Float64,1}

Outputs:
* `dW`: gradients with respect to weights
* `db`: gradients with respect to biases
* `loss`: loss function value
"""
function grad_loss_1layer(
        f_a::Function,
        df_a::Function,
        x::AbstractMatrix,
        y::AbstractMatrix,
        W::AbstractMatrix,
        b::AbstractVector
    )

    n, d = size(W)
    N = size(y, 2) # assume y is matrix of size n x N

    dW = zeros(n, d)
    db = zeros(n)
    loss = 0.0

    for k in 1:N
        for p in 1:n
            error = y[p, k] - f_a(W[p, :]' * x[:, k] + b[p])
            common_term = error * df_a(W[p, :]' * x[:, k] + b[p])
            for q in 1:d
                dW[p, q] = dW[p, q] - 2 / N * common_term * x[q, k]
            end
            db[p] = db[p] - 2 / N * common_term
            loss += 1 / N * error^2
        end
    end
    return dW, db, loss
end

"""
    dW, db, loss = grad_loss_1layer_1output(f_a, df_a, X, y, W, b)

Inputs:
* `X` is a d x N Array{Float64,2} -- even when d= 1
* `y` is a 1 x N element Array{Float64,2}
* `W` is an n x d Array{Float64,2}
* `b` is n element Array{Float64,1}

Outputs:
* `dW`: vector of gradients with respect to weights
* `db`: gradient with respect to bias
* `loss`: loss function value
"""
function grad_loss_1layer_1output(
        f_a::Function,
        df_a::Function,
        X::AbstractMatrix,
        y::AbstractMatrix,
        W::AbstractMatrix,
        b::AbstractVector
    )

    n, d = size(W)
    N = size(X, 2)

    dW = zeros(n, d)
    db = zeros(n)
    loss = 0.0

    for k in 1:N
        error = y[k] - sum(f_a.(W * X[:, k] + b))
        for p in 1:n
            for q in 1:d
                dW[p, q] -= 2 / N * error * df_a(W[p, :]' * X[:, k] + b[p]) * X[q, k]
            end
            db[p] -= 2 / N * error * df_a(W[p, :]' * X[:, k] + b[p])
        end

        loss += 1/N * error^2
    end
    return dW, db, loss
end

"""
    W, b, loss = learn2classify_asgd_1layer(f_a, df_a, grad_loss, X, y, W0, b0[,
        mu=1e-3, iters=500, batch_size=10, seed=0])

Perform accelerated stochastic gradient descent and return optimized weights
    `W`, bias `b`, and vector of loss function values `loss`.

**Inputs**
* `f_a`, `df_a`: activation function and its derivative, respectively
* `grad_loss`: gradient of the loss function
* `X`: d x N matrix representing N samples having dimension d
* `y`: 1 x N row matrix containing scalar function values for the N input samples
* `W0`: n x d initial weight matrix, where n is the hidden layer dimension
* `b0`: length-n initial bias vector
* `mu`: gradient descent step size
* `iters`: number of gradient descent iterations to perform
* `batch_size`: number of samples to use for each iteration
* `seed`: specify a random seed for batch selection

**Outputs**
* `W`: optimized matrix of weights
* `b`: optimized bias vector
* `loss`: vector of loss function values for all iterations
"""
function learn2classify_asgd_1layer(
        f_a::Function,
        df_a::Function,
        grad_loss::Function,
        X::AbstractMatrix,
        y::AbstractMatrix,
        W0::AbstractMatrix,
        b0::AbstractVector,
        mu::Number=1e-3,
        iters::Integer=500,
        batch_size::Integer=10,
        seed::Integer=0
    )

    (seed != 0) && seed!(seed) # use seed if provided

    d = size(W0, 2) #number of inputs
    n = size(W0, 1) # number of neurons
    N = size(X, 2) # number of training samples

    W = W0
    b = b0

    loss = zeros(iters)

    lambda_k = 0
    q_k = W
    p_k = b
    for i in 1:iters
        batch_idx = randperm(N)
        batch_idx = batch_idx[1:min(batch_size, N)]

        dW, db, loss_i = grad_loss(f_a, df_a, X[:, batch_idx], y[:, batch_idx], W, b)

        q_kp1 = W - mu * dW
        p_kp1 = b - mu * db

        lambda_kp1 = (1 + sqrt(1 + 4 * lambda_k^2)) / 2
        gamma_k = (1 - lambda_k) / lambda_kp1

        W = (1 - gamma_k) * q_kp1 + gamma_k * q_k
        b = (1 - gamma_k) * p_kp1 + gamma_k * p_k

        q_k = q_kp1
        p_k = p_kp1
        lambda_k = lambda_kp1

        loss[i] = loss_i
    end
    return W, b, loss
end

export g, grad_loss_1layer, grad_loss_1layer_1output, learn2classify_asgd_1layer, linear, dlinear, dtanh, sigmoid, dsigmoid

end # module

max_stars_repo_name: jevogel/ml

max_issues_repo_name: jevogel/ml

max_forks_repo_name: jevogel/ml


text:
include("common.jl")
jllversion=v"1.7.0"
build_julia(ARGS, v"1.6.0"; jllversion)
build_julia(ARGS, v"1.7.0-beta3"; jllversion)

max_stars_repo_name: singularitti/Yggdrasil

max_issues_repo_name: singularitti/Yggdrasil

max_forks_repo_name: singularitti/Yggdrasil


text:
__precompile__()

module CUDArt

import CUDAdrv
using Compat

export
    # pointer symbols
    CudaPtr, rawpointer, CUDA_NULL,
    # array symbols
    AbstractCudaArray, AbstractCudaVector, AbstractCudaMatrix,
    CudaArray, CudaVector, CudaMatrix, CudaVecOrMat,
    CudaPitchedArray, HostArray, to_host,
    # other symbols
    device, devices, devcount, device_reset, attribute, capability,
    driver_version, runtime_version,
    pitchel, pitchbytes,
    device_synchronize, synchronize,
    Stream, null_stream, cudasleep,
    destroy, free, cudafinalizer, wait_free

import Base: ==, -, +, getindex, setindex!,
    length, size, ndims, eltype, similar, pointer, stride,
    copy, convert, reinterpret, show, summary,
    copy!, get!, fill!, wait, unsafe_convert, vec

# Prepare the CUDA runtime API bindings
include("libcudart-6.5.jl")
import .CUDArt_gen
const rt = CUDArt_gen

const ext = joinpath(dirname(@__DIR__), "deps", "ext.jl")
isfile(ext) || error("CUDArt.jl has not been built, please run Pkg.build(\"CUDArt\").")
include(ext)
if !configured
    # default (non-functional) values for critical variables,
    # making it possible to _load_ the package at all times.
    const libcuda_path = nothing
    const libnvml_path = nothing
    const nvidiasmi_path = nothing
end
const libcudart = libcudart_path
const libnvml = libnvml_path
const nvidia_smi = nvidiasmi_path

include("version.jl")
include("types.jl")
include("device.jl")
include("stream.jl")
include("pointer.jl")
include("arrays.jl")
include("compile.jl")
include("execute.jl")

include("precompile.jl")

function __init__()
    if !configured
        warn("CUDArt.jl has not been successfully built, and will not work properly.")
        warn("Please run Pkg.build(\"CUDArt\") and restart Julia.")
        return
    end

    c_async_send_cudastream[] = cfunction(async_send_cudastream, Void, (rt.cudaStream_t, rt.cudaError_t, Ptr{Void}))
end

end


text:
ProjDir = @__DIR__

include("./gen6-6-8.jl")

fig = plot_models([m6_6s, m6_7s, m6_8s], :waic)
plot(fig)


text:
##################################################
# Polyacetylene
##################################################

# Find the actual crossing point
x0, y0 = TheoryOfCISS.find_crossing(
    Polyacetylene(N = 38),
    xsymbol = :θ,
    x0 = 1.5,
    y0 = 4.9,
    threshold=1e-6,
    attenuation=1e-2,
    correct_overlaps = false
);

# Calculate the D and C vectors
vs = TheoryOfCISS.atpoint(
    polarization,
    x0,
    y0,
    Polyacetylene(N = 38),
    xsymbol = :θ,
    ysymbol = :E,
    calc = TheoryOfCISS.calc_data1,
    correct_overlaps = false,
);

##################################################
# Helicene
##################################################

let
    f = polarization
    xsymbol = :δz
    x0 = 0.8
    ysymbol = :E
    y0 = -3.5
    mol = Helicene(N=7)

    x0, y0 = TheoryOfCISS.find_crossing(
        x0, y0, mol,
        xsymbol = xsymbol,
        threshold=1e-6,
        attenuation=1e-2
    )

    TheoryOfCISS.atpoint(
        f, x0, y0, mol,
        xsymbol = xsymbol,
        ysymbol = ysymbol,
        calc = TheoryOfCISS.calc_data1,
        threshold=1e-6,
    )
end

###

let
    f = x -> 100real(x.s) / real(x.t)
    xsymbol = :δz
    x0 = 0.8
    ysymbol = :E
    y0 = -3.5
    mol = Helicene(N=7)

    x0, y0 = TheoryOfCISS.find_crossing(
        x0, y0, mol,
        xsymbol = xsymbol,
        threshold=1e-6,
        attenuation=1e-2
    )

    TheoryOfCISS.atpoint(
        f, x0, y0, mol,
        xsymbol = xsymbol,
        ysymbol = ysymbol,
        calc = TheoryOfCISS.calc_data6,
        threshold=1e-6,
    )
end


text:
# This is a simple example showing how forward and backward Euler
# could be wrapped

module InternalEuler

using DiffEqBase, LinearAlgebra

# make a algorithm type
abstract type EulerAlgs <: DiffEqBase.AbstractODEAlgorithm end

struct FwdEulerAlg <: EulerAlgs end

struct BwdEulerAlg <: EulerAlgs end

function DiffEqBase.solve(prob::DiffEqBase.AbstractODEProblem{uType,tType,isinplace},
                          Alg::FwdEulerAlg;
                          dt=(prob.tspan[2]-prob.tspan[1])/100,
                          tstops=tType[],
                          kwargs... # ignored kwargs
                          ) where {uType,tType,isinplace}
    u0 = prob.u0
    f = prob.f
    tspan = prob.tspan
    p = prob.p

    if isempty(tstops)
        tstops = tspan[1]:dt:tspan[2]
    end
    @assert tstops[1]==tspan[1]

    nt = length(tstops)
    out = Vector{uType}(undef, nt)
    out[1] = copy(u0)
    tmp = copy(u0)
    for i=2:nt
        t = tstops[i]
        dt = t-tstops[i-1]
        if isinplace
            f(tmp,out[i-1],p,t)
            out[i] = out[i-1] + dt*tmp
        else
            out[i] = out[i-1] + dt*f(out[i-1],p,t)
        end
    end
    # make solution type
    DiffEqBase.build_solution(prob, Alg, tstops, out)
end

function DiffEqBase.solve(prob::DiffEqBase.AbstractODEProblem{uType,tType,isinplace},
                          Alg::BwdEulerAlg;
                          dt=(prob.tspan[2]-prob.tspan[1])/100,
                          tstops=tType[],
                          tol=1e-5,
                          maxiter=100,
                          kwargs... # ignored kwargs
                          ) where {uType,tType,isinplace}
    u0 = prob.u0
    f = prob.f
    tspan = prob.tspan
    p = prob.p
    # TODO: fix numparameters as it picks up the Jacobian
#    @assert !isinplace "Only out of place functions supported"
    @assert DiffEqBase.has_jac(f) "Provide Jacobian as f(::Val{:jac}, ...)"
    jac = (u,p,t) -> f(Val{:jac}(),u,p,t)

    if isempty(tstops)
        tstops = tspan[1]:dt:tspan[2]
    end
    @assert tstops[1]==tspan[1]

    nt = length(tstops)
    out = Vector{uType}(nt)
    out[1] = copy(u0)
    for i=2:nt
        t = tstops[i]
        dt = t-tstops[i-1]
        out[i] = newton(t, dt, out[i-1], p, f, jac, tol, maxiter)
    end
    # make solution type
    DiffEqBase.build_solution(prob, Alg, tstops, out)
end

function newton(t, dt, u_last, p, f, jac, tol, maxiter)
    res = (u) -> u - u_last - dt*f(u,p,t)
    u = u_last + dt*f(u_last,p,t) # forward Euler step as first guess
    for i=1:maxiter
        du = -(I - dt*jac(u,p,t))\res(u)
        u += du
        norm(du, Inf)<tol && return u
    end
    error("Newton not converged")
end

end


text:
using DataFrames
using NamedArrays
using FunctionEnvelope
using SchumakerSpline

bounds = (0.0, 1.0)
x_name = :alpha
evaluations = DataFrame()

function func_inc1(x)
    feval = NamedArray([x, -sqrt(x)])
    setnames!(feval, ["utility", "x_E"], 1)
    return feval
end

function func_inc2(x)
    feval = NamedArray([sqrt(x), sqrt(x)])
    setnames!(feval, ["utility", "x_E"], 1)
    return feval
end

function func_inc3(x)
    feval = NamedArray([3*(x^2), sqrt(x)])
    setnames!(feval, ["utility", "x_E"], 1)
    return feval
end

function dec1(x)
    feval = NamedArray([-sqrt(x), sqrt(x)])
    setnames!(feval, ["utility", "x_E"], 1)
    return feval
end
function dec2(x)
    feval = NamedArray([-3*x^2 + 1, sqrt(x)])
    setnames!(feval, ["utility", "x_E"], 1)
    return feval
end

const_Array = NamedArray([0.1, sqrt(0.2)])
setnames!(const_Array, ["utility", "x_E"], 1)

# Testing basic evaluations
dd = evaluate_function(func_inc1,0.3,:func_inc1, :example_x)
size(dd) == (1, 4)
output, dd2 = evaluate_function_and_add_to_dd(dd, func_inc2,0.4,:func_inc2, :example_x, :utility)
size(dd2) == (2, 4)
output, dd3 = evaluate_function_and_add_to_dd(dd2, dec1,0.4,:dec1, :example_x, :utility)
size(dd3) == (3, 4)
output, dd4 = evaluate_function_and_add_to_dd(dd3, const_Array,0.4,:const_Array, :example_x, :utility)
size(dd4) == (4, 4)

# Testing the construction of a dataframe of functions
dd_function = make_function_dataframe([:func_inc1, :func_inc2, :func_inc3, :dec1, :dec2, :constant],
                                      [func_inc1, func_inc2, func_inc3, dec1, dec2, const_Array],
                                      [:Increasing, :Increasing, :Increasing, :Decreasing, :Decreasing, :Constant])
dd_function[2,:func] == func_inc2

# Testing the crawler.
step = 0.11
limits = (0.0001, 0.9999)
step_right_to_left = true
x_name = :alpha
output_for_envelope = :utility

# Getting Initial evaluations
dd_eval = initial_evaluations(dd_function, limits, step_right_to_left, x_name, output_for_envelope)

# Testing
target = 0.8
max_vals =  maximum_possible_value(target, dd_function, dd_eval, output_for_envelope, x_name)

# Initial evaluations
grid = reverse(0.0001:step:0.9999)
dd_eval = crawler(dd_function, grid, x_name, output_for_envelope; do_all_at_ends = true)

x_array = sort(unique(dd_eval[!,x_name]))
bracketing_parameter = 0.8
max_interval_width = 0.01
func_name = :func_inc2
interval = Tuple{Float64,Float64}([0.01, 0.8])
type = :Increasing
constant = 0.2
dd_eval_rows = size(dd_eval)[1]
# This should have one intercept at x = 0.2^2 = 0.04
interc, dd_eval = get_intercept_with_constant(dd_eval, func_name, type, constant, dd_function, interval, x_name, output_for_envelope, bracketing_parameter, max_interval_width)
size(dd_eval)[1] > dd_eval_rows
abs(interc - 0.04) < max_interval_width
# This should have one intercept at x = 0.2
func_name = :func_inc1
interc, dd_eval = get_intercept_with_constant(dd_eval, func_name, type, constant, dd_function, interval, x_name, output_for_envelope, bracketing_parameter, max_interval_width)
abs(interc - 0.2) < max_interval_width

# These should intercept at x = 1/3
increasing_function = :func_inc3
decreasing_function = :func_inc1
interval = Tuple{Float64,Float64}([0.1,0.6])
interp, dd_eval = get_intercept(dd_eval, increasing_function, decreasing_function, dd_function, interval,
                                x_name, output_for_envelope, bracketing_parameter, max_interval_width)
abs(interp - (1/3)) < max_interval_width

# These should intercept at x = (-1 + sqrt(13))/(6)
increasing_function = :func_inc1
decreasing_function = :dec2
interval = Tuple{Float64,Float64}([0.1,0.6])
interp, dd_eval = get_intercept(dd_eval, increasing_function, decreasing_function, dd_function, interval,
                                x_name, output_for_envelope, bracketing_parameter, max_interval_width)
analytical_intercept = (-1 + sqrt(13))/(6) # Just do quadratic formula to get this.
abs(interp - analytical_intercept) < max_interval_width

# Based on the above functions
dd_interval, dd_eval = get_upper_envelope(dd_eval, dd_function, x_name, output_for_envelope, bracketing_parameter, max_interval_width; grid = reverse(collect(grid)))
env = Envelope(dd_interval, dd_eval, dd_function, x_name, output_for_envelope)

plt = plot(env; points = true, legend = :none)
plt
plot(plt; legend=:best)

# Testing nice syntax for evaluation.
abs(env(1.0) - 2.998) < 1e-2


text:
using ReservoirComputing, MLJLinearModels

#model parameters
const approx_res_size = 30
const radius = 1.2
const activation = tanh
const degree = 6
const sigma = 0.1
const beta = 0.0
const alpha = 1.0
const nla_type = NLADefault()
const in_size = 6
const out_size = 3
const γ = 0.5
const prior_model_size = 3
const u0 = [1.0,0.0,0.0]
const tspan = (0.0,1000.0)
const datasize = 50
W_in = ReservoirComputing.physics_informed_input(approx_res_size, in_size, sigma, γ, prior_model_size)
W = ReservoirComputing.init_reservoir_givendeg(approx_res_size, radius, degree)
const extended_states = false
const h_steps = 2

const train_len = 50
const predict_len = 12
data = ones(Float64, in_size-out_size, 100).+γ

train = data[:, 1:train_len]

train = data[:, 1:1+train_len-1]
test = data[:, train_len:train_len+predict_len-1]

#user physics function
function lorenz()
end

#physics data generator for training and prediction
trange = collect(range(tspan[1], tspan[2], length = train_len))
dt = trange[2]-trange[1]
tsteps = push!(trange, dt + trange[end])
tspan_new = (tspan[1], dt+tspan[2])

function prior_model(u0, tspan_new, tsteps, model = lorenz)
    sol = ones(length(u0), length(tsteps)).*γ
    return sol
end
physics_model_data = prior_model(u0, tspan_new, tsteps)

#constructor 1
hesn = HESN(W,
    train,
    prior_model,
    u0,
    tspan,
    datasize,
    W_in,
    activation = activation,
    alpha = alpha,
    nla_type = nla_type,
    extended_states = extended_states)

#test constructor
@test isequal(Integer(floor(approx_res_size/in_size)*in_size), hesn.res_size)
@test isequal(train, hesn.train_data)
@test isequal(prior_model(u0, tspan_new, tsteps), hesn.physics_model_data)
@test isequal(vcat(train, physics_model_data[:, 1:end-1]), vcat(hesn.train_data, hesn.physics_model_data[:,1:end-1]))
@test isequal(alpha, hesn.alpha)
@test isequal(activation, hesn.activation)
@test isequal(nla_type, hesn.nla_type)
@test isequal(prior_model, hesn.prior_model)
@test isequal(datasize, hesn.datasize)
@test isequal(u0, hesn.u0)
@test isequal(tspan, hesn.tspan)
@test isequal(dt, hesn.dt)
@test size(hesn.W) == (hesn.res_size, hesn.res_size)
@test size(hesn.W_in) == (hesn.res_size, hesn.in_size)
@test size(hesn.states) == (hesn.res_size, train_len)

#test dimension mismatch of hesn constructor
bad_in_size = 4
W_in = ReservoirComputing.physics_informed_input(approx_res_size, bad_in_size, sigma, γ, prior_model_size)
@test_throws DimensionMismatch hesn2 = HESN(W, train, prior_model, u0, tspan, datasize, W_in, activation = activation, alpha = alpha, nla_type = nla_type, extended_states = extended_states)

bad_res_size = approx_res_size-1
W_in = ReservoirComputing.physics_informed_input(bad_res_size, in_size, sigma, γ, prior_model_size)
@test_throws DimensionMismatch hesn2 = HESN(W, train, prior_model, u0, tspan, datasize, W_in, activation = activation, alpha = alpha, nla_type = nla_type, extended_states = extended_states)

#test train
linear_model = Ridge(beta, Analytical())
W_out = HESNtrain(linear_model, hesn)
@test size(W_out) == (out_size, hesn.res_size+size(hesn.physics_model_data,1))

#test predict
output = HESNpredict(hesn, predict_len, W_out)
@test size(output) == (out_size, predict_len)


text:
# for use with static, axis-symmetric metrics
# create new abstract type for easy re-definition
abstract type AbstractAutoDiffMetricParams{T} <: AbstractMetricParams{T} end
abstract type AbstractAutoDiffStaticAxisSymmetricParams{T} <:
              AbstractAutoDiffMetricParams{T} end

# interface
metric_components(m::AbstractAutoDiffStaticAxisSymmetricParams{T}, rθ) where {T} =
    error("Not implemented for $(typeof(m)).")

"""
    inverse_metric_components(g_comp)

Calculates ``g^{tt}``, ``g^{rr}``, ``g^{\\theta\\theta}``, ``g^{\\phi\\phi}``, ``g^{t\\phi}`` of a static,
axis-symmetric metric from ``g_{tt}``, ``g_{rr}``, ``g_{\\theta\\theta}``, ``g_{\\phi\\phi}``, ``g_{t\\phi}``
using a symbolically computed inverse matrix method.

## Developer notes

To recreate:

```julia
using Symbolics
@variables g[1:5] # non zero metric components
metric = [
    g[1] 0 0 g[5]
    0 g[2] 0 0
    0 0 g[3] 0
    g[5] 0 0 g[4]
]
inv(metric)
```
"""
@fastmath function inverse_metric_components(g_comp)
    @inbounds let g1 = g_comp[1],
        g2 = g_comp[2],
        g3 = g_comp[3],
        g4 = g_comp[4],
        g5 = g_comp[5]

        (
            (g2 * g3 * g4) / (g1 * g2 * g3 * g4 - (g5^2) * g2 * g3),
            (g1 * g3 * g4 - (g5^2) * g3) / (g1 * g2 * g3 * g4 - (g5^2) * g2 * g3),
            (g1 * g2 * g4 - (g5^2) * g2) / (g1 * g2 * g3 * g4 - (g5^2) * g2 * g3),
            (g1 * g2 * g3) / (g1 * g2 * g3 * g4 - (g5^2) * g2 * g3),
            (-g2 * g3 * g5) / (g1 * g2 * g3 * g4 - (g5^2) * g2 * g3)
        )
    end
end

"""
    calc_symb_geodesic_equation(g1inv, j1, j2, v)

Symbolically pre-computed geodesic equations for a static, axis-symmetric metric.

## Developer notes

To recreate:

```julia
using Symbolics, Tullio
@variables ginv[1:5], j1[1:5], j2[1:5], v[1:4] # non zero metric components
inverse_metric = [
    ginv[1] 0 0 ginv[5]
    0 ginv[2] 0 0
    0 0 ginv[3] 0
    ginv[5] 0 0 ginv[4]
]
j1_mat = [
    j1[1] 0 0 j1[5]
    0 j1[2] 0 0
    0 0 j1[3] 0
    j1[5] 0 0 j1[4]
]
j2_mat = [
    j2[1] 0 0 j2[5]
    0 j2[2] 0 0
    0 0 j2[3] 0
    j2[5] 0 0 j2[4]
]
j0 = zeros(Float64, (4, 4))
jacobian = (j0, j1_mat, j2_mat, j0)
# christoffel symbols
@tullio Γ[i, k, l] :=
    1 / 2 *
    inverse_metric[i, m] *
    (jacobian[l][m, k] + jacobian[k][m, l] - jacobian[m][k, l])
# compute geodesic equation
@tullio δxδλ[i] := -v[j] * Γ[i, j, k] * v[k]
```
"""
@fastmath function compute_geodesic_equation(ginv, j1, j2, v::SVector{4,T}) where {T}
    @inbounds @SVector [
        -2(0.5ginv[5] * j1[4] + 0.5ginv[1] * j1[5]) * v[2] * v[4] -
        2(0.5ginv[1] * j1[1] + 0.5ginv[5] * j1[5]) * v[1] * v[2] -
        2(0.5ginv[1] * j2[1] + 0.5ginv[5] * j2[5]) * v[1] * v[3] -
        2(0.5ginv[5] * j2[4] + 0.5ginv[1] * j2[5]) * v[3] * v[4],
        0.5(v[1]^2) * ginv[2] * j1[1] +
        0.5(v[3]^2) * ginv[2] * j1[3] +
        0.5(v[4]^2) * ginv[2] * j1[4] +
        ginv[2] * j1[5] * v[1] * v[4] - 0.5(v[2]^2) * ginv[2] * j1[2] -
        ginv[2] * j2[2] * v[2] * v[3],
        0.5(v[1]^2) * ginv[3] * j2[1] +
        0.5(v[2]^2) * ginv[3] * j2[2] +
        0.5(v[4]^2) * ginv[3] * j2[4] +
        ginv[3] * j2[5] * v[1] * v[4] - 0.5(v[3]^2) * ginv[3] * j2[3] -
        ginv[3] * j1[3] * v[2] * v[3],
        -2(0.5ginv[4] * j1[4] + 0.5ginv[5] * j1[5]) * v[2] * v[4] -
        2(0.5ginv[4] * j2[4] + 0.5ginv[5] * j2[5]) * v[3] * v[4] -
        2(0.5ginv[5] * j1[1] + 0.5ginv[4] * j1[5]) * v[1] * v[2] -
        2(0.5ginv[4] * j2[5] + 0.5ginv[5] * j2[1]) * v[1] * v[3]
    ]
end

@fastmath function constrain_time(g_comp, v, μ=0.0, positive::Bool = true)
    @inbounds begin
        discriminant = (
            -g_comp[1] * g_comp[2] * v[2]^2 - g_comp[1] * g_comp[3] * v[3]^2 -
            g_comp[1] * μ^2 - (g_comp[1] * g_comp[4] - g_comp[5]^2) * v[4]^2
        )
        if positive
            -(g_comp[5] * v[4] + √discriminant) / g_comp[1]
        else
            -(g_comp[5] * v[4] - √discriminant) / g_comp[1]
        end
    end
end

function constrain(m::AbstractAutoDiffStaticAxisSymmetricParams{T}, u, v; μ::T = 0.0) where {T}
    rθ = @SVector [u[2], u[3]]
    g_comps = metric_components(m, rθ)
    constrain_time(g_comps, v, μ)
end

@inbounds function geodesic_eq(
    m::AbstractAutoDiffStaticAxisSymmetricParams{T},
    u,
    v
) where {T}
    # get the only position components we need for this metric type
    rθ = @SVector [u[2], u[3]]
    # calculate all non-zero components
    g_comps = metric_components(m, rθ)
    # use AD to get derivatives
    jacs = ForwardDiff.vector_mode_jacobian(x -> metric_components(m, x), rθ)
    # calculate all non-zero inverse matric components
    ginv_comps = inverse_metric_components(g_comps)
    # calculate acceleration
    compute_geodesic_equation(ginv_comps, jacs[:, 1], jacs[:, 2], v)
end

export AbstractAutoDiffStaticAxisSymmetricParams


text:
using MPI
using Oceananigans
using Oceananigans.Distributed

MPI.Init()

comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)
Nranks = MPI.Comm_size(comm)

topology = (Periodic, Periodic, Flat)
arch = MultiArch(CPU(); topology, ranks=(1, Nranks, 1))
grid = RectilinearGrid(arch; topology, size=(16, 16), halo=(3, 3), extent=(2π, 2π))

model = NonhydrostaticModel(; grid)

uᵢ = rand(size(grid)...)
vᵢ = rand(size(grid)...)
set!(model, u=uᵢ, v=vᵢ)

u, v, w = model.velocities
ζ = ∂x(v) - ∂y(u)

simulation = Simulation(model, Δt=0.01, stop_iteration=3)

simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.velocities, (; ζ)),
                                                      schedule = IterationInterval(1),
                                                      with_halos = true,
                                                      filename = "test_output_writing_rank$rank",
                                                      overwrite_existing = true)

run!(simulation)


text:
using PortAudio, SampledSignals, LibSndFile
using HTTP
using WAV
using Statistics
using Plots
using Images
using PyCall
using DSP

function spectronew(val)
    fs = 16000
    audio_test = SampleBuf(val, 16000)
    n = length(audio_test.data)
    nw = n÷25
    spec = DSP.spectrogram(audio_test.data[:], nw, nw÷10; fs=fs)
    png(
        Plots.heatmap(spec.time, spec.freq, pow2db.(spec.power), xaxis = false, yaxis=false, legend = :none, ticks=false, margins = -1.2Plots.cm, size=(28.5,29)),
        "testcase",
    )
    img = load("testcase.png")
    img = imresize(img, (28,28))
    img = reshape(convert(Array{Float32}, Gray.(img)), (28, 28, 1));
    return img
end

function sendsingledataup(data)
   r  = HTTP.post("http://localhost:8003/adddata", [], string(spectronew(data)))
   # println(String(r.body))
end

labeldict = Dict("left" => 0, "right" => 1, "up" => 2,  "down" => 3,""=>4)

function makeamove(mode)
    py"""
    import pyautogui as autogui
    def makemode(mode):
        #autogui.moveTo(200, 200, duration=1)
        if mode == 0:
            print('Mode: Left')
            #autogui.click(clicks=2, button='left')
            autogui.press('left')
        elif mode == 1:
            print('Mode: Right')
            autogui.press('right')
        elif mode == 2:
            print('Mode: Up')
            #autogui.click(clicks=2, button='up')
            autogui.press('up')
        elif mode == 3:
            print('Mode: Down')
            autogui.press('down')
        else:
            print('Mode: NA')
    """
    py"makemode"(mode)
end

py"""
import collections
import contextlib
import sys
import wave

import webrtcvad

def read_wave(path):
    with contextlib.closing(wave.open(path, 'rb')) as wf:
        num_channels = wf.getnchannels()
        assert num_channels == 1
        sample_width = wf.getsampwidth()
        assert sample_width == 2
        sample_rate = wf.getframerate()
        assert sample_rate in (8000, 16000, 32000, 48000)
        pcm_data = wf.readframes(wf.getnframes())
        return pcm_data, sample_rate

def write_wave(path, audio, sample_rate):
    with contextlib.closing(wave.open(path, 'wb')) as wf:
        wf.setnchannels(1)
        wf.setsampwidth(2)
        wf.setframerate(sample_rate)
        wf.writeframes(audio)

class Frame(object):
    def __init__(self, bytes, timestamp, duration):
        self.bytes = bytes
        self.timestamp = timestamp
        self.duration = duration

def frame_generator(frame_duration_ms, audio, sample_rate):
    n = int(sample_rate * (frame_duration_ms / 1000.0) * 2)
    offset = 0
    timestamp = 0.0
    duration = (float(n) / sample_rate) / 2.0
    while offset + n < len(audio):
        yield Frame(audio[offset:offset + n], timestamp, duration)
        timestamp += duration
        offset += n

def vad_collector(sample_rate, frame_duration_ms,
                  padding_duration_ms, vad, frames):
    num_padding_frames = int(padding_duration_ms / frame_duration_ms)
    ring_buffer = collections.deque(maxlen=num_padding_frames)
    triggered = False

    voiced_frames = []
    for frame in frames:
        is_speech = vad.is_speech(frame.bytes, sample_rate)

        #sys.stdout.write('1' if is_speech else '0')
        if not triggered:
            ring_buffer.append((frame, is_speech))
            num_voiced = len([f for f, speech in ring_buffer if speech])
            if num_voiced > 0.9 * ring_buffer.maxlen:
                triggered = True
                # sys.stdout.write('+(%s)' % (ring_buffer[0][0].timestamp,))
                for f, s in ring_buffer:
                    voiced_frames.append(f)
                ring_buffer.clear()
        else:
            voiced_frames.append(frame)
            ring_buffer.append((frame, is_speech))
            num_unvoiced = len([f for f, speech in ring_buffer if not speech])
            if num_unvoiced > 0.9 * ring_buffer.maxlen:
                # sys.stdout.write('-(%s)' % (frame.timestamp + frame.duration))
                triggered = False
                yield b''.join([f.bytes for f in voiced_frames])
                ring_buffer.clear()
                voiced_frames = []
    # if triggered:
        # print("Here I arrive")
        # sys.stdout.write('-(%s)' % (frame.timestamp + frame.duration))
    #sys.stdout.write('\n')
    if voiced_frames:
        yield b''.join([f.bytes for f in voiced_frames])

def doyourvaejobplease(path):
    audio, sample_rate = read_wave(path)
    vad = webrtcvad.Vad(int(1))
    frames = frame_generator(30, audio, sample_rate)
    frames = list(frames)
    segments = vad_collector(sample_rate, 30, 300, vad, frames)
    for i, segment in enumerate(segments):
        path = './audios/chunk-%002d.wav' % (i,)
        write_wave(path, segment, sample_rate)
"""

using Glob
while true
    stream = PortAudioStream("default", 1, 0)
    buf = read(stream, 2s).data[:]
    wavwrite(buf,  "test8final.wav", Fs=48000, nbits=16,  compression=WAVE_FORMAT_PCM)
    py"doyourvaejobplease"("test8final.wav")
    close(stream)
    listofwavs = glob("*.wav","./audios")
    interestingdata = map(x->wavread(x)[1][:], listofwavs)
    for i in interestingdata
        datara = ceil(Int, length(i)/16000)
        zerosadd =  datara*16000 - length(i)
        zeroestobeadded = zeros(Float64, zerosadd)
        i = vcat(i, zeroestobeadded)
        batch = map(x->i[16000*(x-1)+1:16000*x] ,1:datara)
    	for j in batch
            sendsingledataup(j)
        end
    end

    # # rm("audios")
    map(x->rm(x), listofwavs)
    interestingdata = nothing
    listofwavs = nothing
    buf = nothing
    r = HTTP.get("http://localhost:8003/digitreg", [])
    mode = labeldict[String(r.body)]
    makeamove(mode)
    #print(varinfo())
    sleep(0.1)
end

max_stars_repo_name: ashwani-rathee/Sardor

max_issues_repo_name: ashwani-rathee/Sardor

max_forks_repo_name: ashwani-rathee/Sardor