text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -32.35079088597315 <= q <= 33.04707617725232)

                     add_NL_constraint(m, :(log(1 + exp(0.7516814527569986 + 0.07211062597860263*log(1 + exp(0.800956367833014 + 0.1575422515944389*log(1 + exp(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4]))) + 0.30181017969165413*log(1 + exp(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4]))) + 0.5024447348308962*log(1 + exp(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4]))) + -0.7463877606567912*log(1 + exp(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))) + -0.28099609060437203*log(1 + exp(-0.8241959750660866 + 0.48752067922981546*log(1 + exp(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4]))) + -0.8942226763869914*log(1 + exp(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4]))) + -0.1882956347370066*log(1 + exp(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4]))) + 0.0613423792438379*log(1 + exp(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))) + 0.19401527858092393*log(1 + exp(-0.014178711859067938 + -0.6496362145686132*log(1 + exp(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4]))) + 0.5744456236018198*log(1 + exp(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4]))) + 0.06990002879133872*log(1 + exp(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4]))) + 0.11738131892224501*log(1 + exp(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))) + -0.43991506890783594*log(1 + exp(-0.38559964579408224 + -0.6882459841770032*log(1 + exp(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4]))) + 0.5869261088566509*log(1 + exp(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4]))) + 0.7286086900917508*log(1 + exp(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4]))) + 0.9385542676650069*log(1 + exp(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))))) + log(1 + exp(0.5442086655639207 + -0.859033717055464*log(1 + exp(0.800956367833014 + 0.1575422515944389*log(1 + exp(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4]))) + 0.30181017969165413*log(1 + exp(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4]))) + 0.5024447348308962*log(1 + exp(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4]))) + -0.7463877606567912*log(1 + exp(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))) + -0.6994107397450926*log(1 + exp(-0.8241959750660866 + 0.48752067922981546*log(1 + exp(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4]))) + -0.8942226763869914*log(1 + exp(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4]))) + -0.1882956347370066*log(1 + exp(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4]))) + 0.0613423792438379*log(1 + exp(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))) + -0.21447765318835765*log(1 + exp(-0.014178711859067938 + -0.6496362145686132*log(1 + exp(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4]))) + 0.5744456236018198*log(1 + exp(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4]))) + 0.06990002879133872*log(1 + exp(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4]))) + 0.11738131892224501*log(1 + exp(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))) + 0.9247182813449877*log(1 + exp(-0.38559964579408224 + -0.6882459841770032*log(1 + exp(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4]))) + 0.5869261088566509*log(1 + exp(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4]))) + 0.7286086900917508*log(1 + exp(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4]))) + 0.9385542676650069*log(1 + exp(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))))) + log(1 + exp(-0.33457802196328057 + 0.8511359109817849*log(1 + exp(0.800956367833014 + 0.1575422515944389*log(1 + exp(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4]))) + 0.30181017969165413*log(1 + exp(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4]))) + 0.5024447348308962*log(1 + exp(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4]))) + -0.7463877606567912*log(1 + exp(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))) + 0.9056512987027685*log(1 + exp(-0.8241959750660866 + 0.48752067922981546*log(1 + exp(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4]))) + -0.8942226763869914*log(1 + exp(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4]))) + -0.1882956347370066*log(1 + exp(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4]))) + 0.0613423792438379*log(1 + exp(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))) + -0.901983268277478*log(1 + exp(-0.014178711859067938 + -0.6496362145686132*log(1 + exp(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4]))) + 0.5744456236018198*log(1 + exp(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4]))) + 0.06990002879133872*log(1 + exp(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4]))) + 0.11738131892224501*log(1 + exp(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))) + -0.02854306019224717*log(1 + exp(-0.38559964579408224 + -0.6882459841770032*log(1 + exp(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4]))) + 0.5869261088566509*log(1 + exp(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4]))) + 0.7286086900917508*log(1 + exp(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4]))) + 0.9385542676650069*log(1 + exp(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))))) + log(1 + exp(0.85223750045835 + -0.7468516096711739*log(1 + exp(0.800956367833014 + 0.1575422515944389*log(1 + exp(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4]))) + 0.30181017969165413*log(1 + exp(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4]))) + 0.5024447348308962*log(1 + exp(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4]))) + -0.7463877606567912*log(1 + exp(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))) + -0.1269236474608202*log(1 + exp(-0.8241959750660866 + 0.48752067922981546*log(1 + exp(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4]))) + -0.8942226763869914*log(1 + exp(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4]))) + -0.1882956347370066*log(1 + exp(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4]))) + 0.0613423792438379*log(1 + exp(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))) + -0.5407492362268407*log(1 + exp(-0.014178711859067938 + -0.6496362145686132*log(1 + exp(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4]))) + 0.5744456236018198*log(1 + exp(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4]))) + 0.06990002879133872*log(1 + exp(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4]))) + 0.11738131892224501*log(1 + exp(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))) + -0.8670343803826994*log(1 + exp(-0.38559964579408224 + -0.6882459841770032*log(1 + exp(-0.7991039641167363 + 0.28486813994855353*$(x[1]) + 0.682420711072186*$(x[2]) + 0.05952972537514656*$(x[3]) + 0.998655851604032*$(x[4]))) + 0.5869261088566509*log(1 + exp(-0.9827986828893067 + 0.8570099551669048*$(x[1]) + -0.029287184975956393*$(x[2]) + 0.28818320654817864*$(x[3]) + -0.9376552616355038*$(x[4]))) + 0.7286086900917508*log(1 + exp(0.33306588949676286 + 0.667638331016386*$(x[1]) + 0.7691056200843742*$(x[2]) + 0.5908191290581821*$(x[3]) + 0.06769592108634237*$(x[4]))) + 0.9385542676650069*log(1 + exp(-0.7350489056562997 + 0.545921427040915*$(x[1]) + -0.3116746162014481*$(x[2]) + 0.844944868364411*$(x[3]) + -0.9086564944285445*$(x[4]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# Networks found by Bose-Nelson algorithm, except for N=16
# http://pages.ripco.net/~jgamble/nw.html
# http://www.cs.brandeis.edu/~hugues/sorting_networks.html
# https://github.com/JeffreySarnoff/SortingNetworks.jl
sorting_network_parameters = Dict(
    2 => (((1,2),),),
    4 => (((1,2), (3,4)), ((1,3), (2,4)), ((2,3),)),
    8 => (
        ((1, 2), (3, 4), (5, 6), (7, 8)), ((1, 3), (2, 4), (5, 7), (6, 8)),
        ((2, 3), (6, 7), (1, 5), (4, 8)), ((2, 6), (3, 7)),
        ((2, 5), (4, 7)), ((3, 5), (4, 6)), ((4, 5),)
    ),
    16 => (
        ((4, 11), (12, 15), (5, 14), (3, 13), (1, 7), (9, 10), (2,8), (6, 16)),
        ((1, 2), (3, 5), (7, 8), (13, 14), (4, 6), (9, 12), (11, 16), (10, 15)),
        ((1, 4), (7, 11), (14, 15), (2, 6), (8, 16), (3, 9), (10, 13), (5, 12)),
        ((1, 3), (8, 14), (15, 16), (2, 5), (6, 12), (4, 9), (11, 13), (7, 10)),
        ((2, 3), (4, 7), (8, 9), (12, 14), (6, 10), (13, 15), (5, 11)),
        ((3, 7), (12, 13), (2, 4), (6, 8), (9, 10), (14, 15)),
        ((3, 4), (5, 7), (11, 12), (13, 14)),
        ((9, 11), (5, 6), (7, 8), (10, 12)),
        ((4, 5), (6, 7), (8, 9), (10, 11), (12, 13)),
        ((7, 8), (9, 10))
    ),
    32 => (
        ((1, 2), (3, 4), (5, 6), (7, 8), (9, 10), (11, 12), (13, 14), (15, 16), (17, 18), (19, 20), (21, 22), (23, 24),
         (25, 26), (27, 28), (29, 30), (31, 32)),
        ((1, 3), (2, 4), (5, 7), (6, 8), (9, 11), (10, 12), (13, 15), (14, 16), (17, 19), (18, 20), (21, 23), (22, 24),
         (25, 27), (26, 28), (29, 31), (30, 32)),
        ((2, 3), (6, 7), (1, 5), (4, 8), (10, 11), (14, 15), (9, 13), (12, 16), (18, 19), (22, 23), (17, 21), (20, 24),
         (26, 27), (30, 31), (25, 29), (28, 32)),
        ((2, 6), (3, 7), (10, 14), (11, 15), (1, 9), (8, 16), (18, 22), (19, 23), (26, 30), (27, 31), (17, 25),
         (24, 32)),
        ((2, 5), (4, 7), (10, 13), (12, 15), (18, 21), (20, 23), (26, 29), (28, 31), (1, 17), (16, 32)),
        ((3, 5), (4, 6), (11, 13), (12, 14), (2, 10), (7, 15), (19, 21), (20, 22), (27, 29), (28, 30), (18, 26), (23, 31)),
        ((4, 5), (12, 13), (2, 9), (3, 11), (6, 14), (8, 15), (20, 21), (28, 29), (18, 25), (19, 27), (22, 30), (24, 31)),
        ((4, 12), (3, 9), (5, 13), (8, 14), (20, 28), (19, 25), (21, 29), (24, 30), (2, 18), (15, 31)),
        ((4, 11), (6, 13), (20, 27), (22, 29), (2, 17), (3, 19), (14, 30), (16, 31)),
        ((4, 10), (7, 13), (20, 26), (23, 29), (3, 17), (16, 30)),
        ((4, 9), (8, 13), (6, 10), (7, 11), (20, 25), (24, 29), (22, 26), (23, 27)),
        ((5, 9), (8, 12), (21, 25), (24, 28), (4, 20), (13, 29)),
        ((6, 9), (8, 11), (22, 25), (24, 27), (4, 19), (5, 21), (12, 28), (14, 29)),
        ((7, 9), (8, 10), (23, 25), (24, 26), (4, 18), (6, 22), (11, 27), (15, 29)),
        ((8, 9), (24, 25), (4, 17), (6, 21), (7, 23), (10, 26), (12, 27), (16, 29)),
        ((8, 24), (7, 21), (5, 17), (6, 18), (9, 25), (12, 26), (15, 27), (16, 28)),
        ((8, 23), (6, 17), (7, 19), (10, 25), (14, 26), (16, 27)), ((8, 22), (7, 17), (11, 25), (16, 26)),
        ((8, 21), (12, 25)), ((8, 20), (13, 25)), ((8, 19), (14, 25), (12, 20), (13, 21)),
        ((8, 18), (15, 25), (11, 19), (14, 22)), ((8, 17), (16, 25), (10, 18), (12, 19), (14, 21), (15, 23)),
        ((9, 17), (12, 18), (16, 24), (15, 21)), ((10, 17), (16, 23), (14, 18), (15, 19)), ((11, 17), (16, 22)),
        ((12, 17), (16, 21)), ((13, 17), (16, 20)), ((14, 17), (16, 19)), ((15, 17), (16, 18)), ((16, 17),)
    )
)

ext: jl


text:
struct Framework
    plugins
    Framework(plugins, hooks=[]; opts...) = new(PluginStack(plugins, hooks; opts...))
end

struct EmptyPlugin{Id} <: Plugin
    EmptyPlugin{Id}(;opts...) where Id = new{Id}()
end
Plugins.symbol(::EmptyPlugin) = :empty
for i = 1:1000
    Plugins.register(EmptyPlugin{i})
end

mutable struct CounterPlugin{Id} <: Plugin
    hook1count::Int
    hook2count::Int
    hook3count::Int
    CounterPlugin{Id}() where Id = new{Id}(0, 0, 0)
end
Plugins.symbol(::CounterPlugin) = :counter
for i = 1:100
    Plugins.register(CounterPlugin{i})
end

@inline hook1(plugin::CounterPlugin, framework) = begin
    plugin.hook1count += 1
    return false
end

hook2_handler(plugin::CounterPlugin, framework) = begin
    plugin.hook2count += 1
    return false
end

@inline hook3(plugin::CounterPlugin, framework, p1, p2) = begin
    plugin.hook3count += p2
    return false
end

chain_of_empties(length=20, startat=0) = [EmptyPlugin{i + startat} for i = 1:length]

callmanytimes(framework, hook, times=1e5) = for i=1:times hook(framework) end

mutable struct FrameworkTestPlugin <: Plugin
    calledwithframework
    FrameworkTestPlugin() = new("Never called")
end
Plugins.register(FrameworkTestPlugin)
hook1(plugin::FrameworkTestPlugin, framework) = plugin.calledwithframework = framework

mutable struct EventTestPlugin <: Plugin
    calledwithframework
    calledwithevent
    EventTestPlugin() = new("Never called", "Never called")
end
Plugins.register(EventTestPlugin)
event_handler(plugin::EventTestPlugin, framework, event) = begin
    plugin.calledwithframework = framework
    plugin.calledwithevent = event
end

struct ConfigurablePlugin <: Plugin
    config::String
    ConfigurablePlugin(;config = "default") = new(config)
end
Plugins.register(ConfigurablePlugin)
checkconfig_handler(plugin::ConfigurablePlugin, framework, event) = begin
    if event.config !== plugin.config
        throw("Not the same!")
    end
    return plugin.config
end

struct PropagationStopperPlugin <: Plugin
end
Plugins.register(PropagationStopperPlugin)
propagationtest(plugin::PropagationStopperPlugin, framework, data) = data == 42
propagationtest_nodata(plugin::PropagationStopperPlugin, framework) = true

struct PropagationCheckerPlugin <: Plugin
end
Plugins.register(PropagationCheckerPlugin)
propagationtest(plugin::PropagationCheckerPlugin, framework, data) = data === 32 || throw("Not 32!")
propagationtest_nodata(plugin::PropagationCheckerPlugin, framework) = throw("Not stopped!")

mutable struct DynamicPlugin <: Plugin
    lastdata
end
Plugins.register(DynamicPlugin)
dynamismtest(plugin::DynamicPlugin, framework, data) = plugin.lastdata = data

# Hook cache test:
mutable struct SharedState
    plugins::PluginStack
    shared_counter::Int
end

struct App{TCache}
    state::SharedState
    hooks::TCache
    function App(plugins, hooklist)
        state = SharedState(PluginStack(plugins, hooklist), 0)
        cache = hook_cache(state.plugins)
        return new{typeof(cache)}(state, cache)
    end
end

const OP_CYCLES = 1e7

function op(app::App)
    counters = [counter for counter in app.state.plugins if counter isa CounterPlugin]
    @info "op: A sample operation on the app, involving hook1() calls in a semi-realistic setting."
    @info "op: $(length(counters)) CounterPlugins found, $(length(app.state.plugins)) plugins in total, each CounterPlugin incrementing a private counter."

    start_ts = time_ns()
    for i in 1:OP_CYCLES
        app.hooks.hook3(app, i, 1)
    end
    end_ts = time_ns()

    for i = 1:length(counters)
        @test counters[i].hook3count == OP_CYCLES
    end

    time_diff = end_ts - start_ts
    avg_calltime = time_diff / OP_CYCLES
    @info "op: $OP_CYCLES hook1() calls took $(time_diff / 1e9) secs. That is $avg_calltime nanosecs per call on average, or $(avg_calltime / length(counters)) ns per in-plugin counter increment."
end

mutable struct LifeCycleTestPlugin <: Plugin
    setupcalledwith
    shutdowncalledwith
    deferredinitcalledwith
    LifeCycleTestPlugin() = new()
end
Plugins.register(LifeCycleTestPlugin)
Plugins.setup!(plugin::LifeCycleTestPlugin, framework) = plugin.setupcalledwith = framework
Plugins.shutdown!(plugin::LifeCycleTestPlugin, framework) = begin
    plugin.shutdowncalledwith = framework
    if framework === 42
        error("shutdown called with 42")
    end
end
deferred_init(plugin::Plugin, ::Any) = true
deferred_init(plugin::LifeCycleTestPlugin, data) = plugin.deferredinitcalledwith = data

@testset "Plugins.jl basics" begin
    @testset "Plugin chain" begin
        a1 = Framework([CounterPlugin{1}, EmptyPlugin])
        @show innerplugin = a1.plugins[:empty]
        @show counter = a1.plugins[:counter]
        a1_hook1s = hooklist(a1.plugins, hook1)
        @test length(a1_hook1s) === 1
        callmanytimes(a1, a1_hook1s)
        @info "$(length(a1.plugins))-length chain, $(length(a1_hook1s)) counter (1e5 cycles):"
        @time callmanytimes(a1, a1_hook1s)
        hooklist(a1.plugins, hook2_handler)(a1)
        @test counter.hook1count == 2e5
        @test counter.hook2count == 1
    end

    @testset "Same plugin twice" begin
        a2 = Framework([CounterPlugin{1}, CounterPlugin{2}])
        innercounter = a2.plugins[2]
        outercounter = a2.plugins[1]
        a2_hook1s = hooklist(a2.plugins, hook1)
        @test length(a2_hook1s) === 2
        callmanytimes(a2, a2_hook1s)
        @info "$(length(a2.plugins))-length chain, $(length(a2_hook1s)) counters (1e5 cycles):"
        @time callmanytimes(a2, a2_hook1s)
        hooklist(a2.plugins, hook2_handler)(a2)
        @test innercounter.hook1count == 2e5
        @test innercounter.hook2count == 1
        @test outercounter.hook1count == 2e5
        @test outercounter.hook2count == 1
    end

    @testset "Chain of empty Plugins to skip" begin
        chainedapp = Framework(vcat([CounterPlugin{1}], chain_of_empties(20), [CounterPlugin{2}], chain_of_empties(20, 21)))
        innerplugin = chainedapp.plugins[22]
        outerplugin = chainedapp.plugins[1]
        chainedapp_hook1s = hooklist(chainedapp.plugins, hook1)
        callmanytimes(chainedapp, chainedapp_hook1s)
        @info "$(length(chainedapp.plugins))-length chain,  $(length(chainedapp_hook1s))  counters (1e5 cycles):"
        @time callmanytimes(chainedapp, chainedapp_hook1s)
        @test outerplugin.hook1count == 2e5
        @test outerplugin.hook2count == 0
        @test innerplugin.hook1count == 2e5
        @test innerplugin.hook2count == 0
    end

    @testset "Unhandled hook returns false" begin
        app = Framework([EmptyPlugin{1}])
        @test hooklist(app.plugins, hook1)() == false
        @test call_optional(hooklist(app.plugins, hook1)) == nothing
    end

    @testset "Framework goes through" begin
        frameworktestapp = Framework([EmptyPlugin{1}, FrameworkTestPlugin])
        hooklist(frameworktestapp.plugins, hook1)(frameworktestapp)
        @test frameworktestapp.plugins[2].calledwithframework === frameworktestapp
    end

    @testset "Event object" begin
        eventtestapp = Framework([EmptyPlugin{1}, EventTestPlugin])
        event = (name="test event", data=42)
        hooklist(eventtestapp.plugins, event_handler)(eventtestapp, event)
        @test eventtestapp.plugins[2].calledwithframework === eventtestapp
        @test eventtestapp.plugins[2].calledwithevent === event
    end

    @testset "Multiple apps with same chain, differently configured" begin
        app2config = "app2config"
        app1 = Framework([EmptyPlugin{1}, ConfigurablePlugin])
        app2 = Framework([EmptyPlugin{2}, ConfigurablePlugin]; config=app2config)
        event1 = (config ="default",)
        event2 = (config = app2config,)
        hooklist(app1.plugins, checkconfig_handler)(app1, event1)
        @test_throws String hooklist(app1.plugins, checkconfig_handler)(app1, event2)
        @test call_optional(hooklist(app1.plugins, checkconfig_handler), app1, event1) == event1.config
        hooklist(app2.plugins, checkconfig_handler)(app2, event2)
        @test_throws String hooklist(app2.plugins, checkconfig_handler)(app2, event1)
        @test call_optional(hooklist(app2.plugins, checkconfig_handler), app2, event2) == event2.config
    end

    @testset "Stopping Propagation" begin
        spapp = Framework([EmptyPlugin{1}, PropagationStopperPlugin, EmptyPlugin{2}, PropagationCheckerPlugin])
        hooklist(spapp.plugins, propagationtest)(spapp, 42) === true # It is stopped so the checker does not throw
        hooklist(spapp.plugins, propagationtest)(spapp, 32) === false # Not stopped but accepted by the checker
        @test_throws String hooklist(spapp.plugins, propagationtest)(spapp, 41)

        @test hooklist(spapp.plugins, propagationtest_nodata)(spapp) === true
    end

    @testset "HookList iteration" begin
        iapp = Framework([EmptyPlugin{1}, CounterPlugin{2}, EmptyPlugin{2}, CounterPlugin{1}, EmptyPlugin{3}])
        c1 = iapp.plugins[4]
        c2 = iapp.plugins[2]
        hookers = [c2, c1]
        @test length(hooklist(iapp.plugins, hook1)) === 2
        i = 1
        for hook in hooklist(iapp.plugins, hook1)
            @test hookers[i] === hook.plugin
            i += 1
        end
    end

    @testset "Accessing plugins directly" begin
        app = Framework([EmptyPlugin{1}, CounterPlugin{1}])
        empty = app.plugins[1]
        counter = app.plugins[2]
        @test get(app.plugins, :empty) === empty
        @test get(app.plugins, :counter) === counter
        @test app.plugins[:empty] === empty
        @test length(app.plugins) == 2
    end

    @testset "Hook cache" begin
        counters = [CounterPlugin{i} for i=2:40]
        empties = [EmptyPlugin{i} for i=1:100]
        pluginarr = [CounterPlugin{1}, empties..., counters...]
        @info "Measuring time to first hook call with $(length(pluginarr)) uniquely typed plugins, $(length(counters) + 1) implementig the hook."
        @time begin
            simpleapp = SharedState(PluginStack(pluginarr, [hook1]), 0)
            simpleapp_hooks = hooks(simpleapp)
            simpleapp_hooks.hook1(simpleapp)
            @test simpleapp.plugins[1].hook1count == 1
        end
    end

    @testset "Hook cache as type parameter" begin
        counters = [CounterPlugin{i} for i=2:2]
        empties = [EmptyPlugin{i} for i=1:100]
        app = App([CounterPlugin{1}, empties..., counters...], [hook3])
        op(app)
    end

    @testset "Lifecycle Hooks" begin
        app = Framework([EmptyPlugin{1}, LifeCycleTestPlugin])
        plugin = app.plugins[2]
        @test setup!(app.plugins, app).allok == true
        @test plugin.setupcalledwith === app

        # Create a non-standard lifecycle hook
        lifecycle_hook = Plugins.create_lifecyclehook(deferred_init)
        @test string(lifecycle_hook) == "deferred_init"
        @test lifecycle_hook(app.plugins, "42").allok === true
        @test plugin.deferredinitcalledwith === "42"

        @test Plugins.shutdown!(app.plugins, app).allok === true
        @test plugin.shutdowncalledwith === app
        notallok = Plugins.shutdown!(app.plugins, 42)
        @test notallok.allok === false
        @test (notallok.results[2] isa Tuple) === true
        @test (notallok.results[2][1] isa Exception) === true
        @test (stacktrace(notallok.results[2][2]) isa AbstractVector{StackTraces.StackFrame}) === true
        @test plugin.shutdowncalledwith === 42
    end

    @testset "Modifying plugins" begin
        c2 = CounterPlugin{2}()
        app = Framework([EmptyPlugin, CounterPlugin{1}], [hook1])
        c1 = app.plugins[2]
        cache = hooks(app)
        cache.hook1(app)
        push!(app.plugins, c2)
        cache = hooks(app)
        cache.hook1(app)
        @test c2.hook1count == 1
        @test c1.hook1count == 2

        @test pop!(app.plugins) === c2
        cache = hooks(app)
        cache.hook1(app)
        @test c2.hook1count == 1
        @test c1.hook1count == 3

        @test popfirst!(app.plugins) isa EmptyPlugin
        cache = hooks(app)
        cache.hook1(app)
        @test c2.hook1count == 1
        @test c1.hook1count == 4

        pushfirst!(app.plugins, c2)
        cache = hooks(app)
        cache.hook1(app)
        @test c2.hook1count == 2
        @test c1.hook1count == 5
        @test app.plugins[1] === c2

        pop!(app.plugins)
        @test length(app.plugins) == 1
        @test isempty(app.plugins) == false
        pop!(app.plugins)
        @test length(app.plugins) == 0
        @test isempty(app.plugins) == true
        @test_throws ArgumentError pop!(app.plugins)

        app = Framework([EmptyPlugin{1}, CounterPlugin{3}], [hook1])
        c3 = app.plugins[:counter]
        @test isempty(app.plugins) == false
        cache = hooks(app)
        cache.hook1(app)
        @test c3.hook1count == 1
        empty!(app.plugins)
        @test isempty(app.plugins) == true
        cache = hooks(app)
        cache.hook1(app)
        @test c3.hook1count == 1

        c5 = CounterPlugin{5}()
        app = Framework([EmptyPlugin{1}, CounterPlugin{4}], [hook1])
        c4 = app.plugins[:counter]
        @test isempty(app.plugins) == false
        cache = hooks(app)
        cache.hook1(app)
        @test c4.hook1count == 1
        @test c5.hook1count == 0
        app.plugins[2] = c5
        cache = hooks(app)
        cache.hook1(app)
        @test c4.hook1count == 1
        @test c5.hook1count == 1
    end
end


text:
"""
function function_value(f::Function, fs::ScalarFunctionSpace{dim,T}, cell::Int,q_point::Int) where {dim,T}
Evaluate function f(x): x ∈ ℝⁿ ↦ ℝ with x given by cell number `cell` and quadrature point with index `qpoint`
    using fs FunctionSpace data
"""
function function_value(f::Function, fs::ScalarFunctionSpace{dim,T}, cell::Int,q_point::Int) where {dim,T}
    coords = get_cell_coordinates(cell, fs.mesh)
    return f(spatial_coordinate(fs, q_point, coords))
end

"""
function spatial_coordinate(fs::ScalarFunctionSpace{dim}, q_point::Int, x::AbstractVector{Vec{dim,T}})
Map coordinates of quadrature point `q_point` of Scalar Function Space `fs`
into domain with vertices `x`
"""
function spatial_coordinate(fs::ScalarFunctionSpace{dim}, q_point::Int, x::AbstractVector{Vec{dim,T}}) where {dim,T}
    n_base_funcs = getngeobasefunctions(fs)
    @assert length(x) == n_base_funcs
    vec = zero(Vec{dim,T})
    @inbounds for i in 1:n_base_funcs
        vec += geometric_value(fs, q_point, i) * x[i]
    end
    return vec
end

# Trial Functions
struct TrialFunction{dim,T,shape,N1}
    fs::DiscreteFunctionSpace
    m_values::Array{T,N1}
    components::Int
end

@inline getnbasefunctions(u::TrialFunction) = getnbasefunctions(u.fs)
@inline getfunctionspace(u::TrialFunction) = u.fs
@inline getnlocaldofs(u::TrialFunction) = getnlocaldofs(getfunctionspace(u))
@inline getncomponents(u::TrialFunction) = u.components

function TrialFunction(fs::ScalarTraceFunctionSpace{dim,T}) where {dim,T}
    mesh = getmesh(fs)
    m_values = fill(zero(T) * T(NaN), getncells(mesh), getnbasefunctions(fs), n_faces_per_cell(mesh))
    return TrialFunction{dim,T,getshape(getfiniteelement(fs)),3}(fs, m_values, 1)
end

function TrialFunction(fs::ScalarFunctionSpace{dim,T}) where {dim,T}
    mesh = getmesh(fs)
    m_values = fill(zero(T) * T(NaN), getncells(mesh), getnbasefunctions(fs))
    return TrialFunction{dim,T,getshape(getfiniteelement(fs)),2}(fs, m_values, 1)
end

function TrialFunction(fs::VectorFunctionSpace{dim,T}) where {dim,T}
    mesh = getmesh(fs)
    m_values = fill(zero(T) * T(NaN), getncells(mesh), getnbasefunctions(fs))
    return TrialFunction{dim,T,getshape(getfiniteelement(fs)),2}(fs, m_values, dim)
end

function TrialFunction(fs::DiscreteFunctionSpace{dim}, components::Int, m_values::Array{T,N}) where {dim,T,N}
    mesh = getmesh(fs)
    @assert size(m_values,1) == getncells(mesh)
    @assert size(m_values,2) == getnbasefunctions(fs)
    return TrialFunction{dim,T,getshape(getfiniteelement(fs)),N}(fs, m_values, components)
end

reference_coordinate(fs::ScalarFunctionSpace{dim,T},
x_ref::Vec{dim,T}, x::Vec{dim,T}) where {dim,T} = fs.Jinv[]⋅(x-x_ref)

"""
function value(u_h::TrialFunction{dim,T}, cell::Int, x::Vec{dim,T})
    get trial function value on cell `cell` at point `x`
"""
function value(u_h::TrialFunction{dim,T}, node::Int, cell::Int) where {dim,T}
    u = zero(T)
    mesh = getmesh(u_h.fs)
    x = mesh.nodes[node].x
    ξ = reference_coordinate(u_h.fs, mesh.nodes[mesh.cells[cell].nodes[1]].x, x)
    for i in 1:getnbasefunctions(u_h.fs)
        u  += u_h.m_values[cell, i]*value(getfiniteelement(u_h.fs), i, ξ)
    end
    return u
end

function nodal_avg(u_h::TrialFunction{dim,T}) where {dim,T}
    mesh = getmesh(u_h.fs)
    nodalu_h = Vector{Float64}(undef, length(mesh.nodes))
    share_count = zeros(Int,length(mesh.nodes))
    fill!(nodalu_h,0)
    @inbounds for (cell_idx, cell) in enumerate(CellIterator(mesh))
        reinit!(u_h, cell)
        for node in getnodes(cell)
            u = value(u_h, node, cell)
            nodalu_h[node] += u
            share_count[node] += 1
        end
    end
    return nodalu_h./share_count
end

function errornorm(u_h::TrialFunction{dim,T}, u_ex::Function, norm_type::String="L2") where {dim,T}
    mesh = getmesh(u_h.fs)
    Etu_h = zero(T)
    if norm_type == "L2"
        n_basefuncs_s = getnbasefunctions(u_h)
        @inbounds for (cell_idx, cell) in enumerate(CellIterator(mesh))
            reinit!(u_h, cell)
            Elu_h = zero(T)
            for q_point in 1:getnquadpoints(u_h.fs)
                dΩ = getdetJdV(u_h.fs, q_point)
                u = zero(T)
                for i in 1:getnbasefunctions(u_h.fs)
                    u  += u_h.m_values[cell.current_cellid[], i]*shape_value(u_h.fs, q_point, i)
                end
                # Integral (u_h - u_ex) dΩ
                Elu_h += (u-u_ex(spatial_coordinate(u_h.fs, q_point, cell.coords)))^2*dΩ
            end
            Etu_h += Elu_h
        end
    else
        throw("Norm $norm_type not available")
    end
    return Etu_h
end


text:
###############################################################################
#
#   fmpz.jl : BigInts
#
###############################################################################

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

export fmpz, FlintZZ, FlintIntegerRing, parent, show, convert, hash, fac, bell,
       binom, isprime, fdiv, cdiv, tdiv, div, rem, mod, gcd, xgcd, lcm, invmod,
       powmod, abs, divrem, isqrt, popcount, prevpow2, nextpow2, ndigits, dec,
       bin, oct, hex, base, one, zero, divexact, fits, sign, nbits, deepcopy,
       tdivpow2, fdivpow2, cdivpow2, flog, clog, cmpabs, clrbit!, setbit!,
       combit!, crt, divisible, divisor_lenstra, fdivrem, tdivrem, fmodpow2,
       gcdinv, isprobabprime, issquare, jacobi, remove, root, size, isqrtrem,
       sqrtmod, trailing_zeros, sigma, eulerphi, fib, moebiusmu, primorial,
       risingfac, numpart, canonical_unit, needs_parentheses, isnegative,
       show_minus_one, parseint, addeq!, mul!, isunit, isequal, num, den,
       iszero, rand

###############################################################################
#
#   Data type and parent methods
#
###############################################################################

parent_type(::Type{fmpz}) = FlintIntegerRing

doc"""
    parent(a::fmpz)
> Returns the unique Flint integer parent object `FlintZZ`.
"""
parent(a::fmpz) = FlintZZ

elem_type(::Type{FlintIntegerRing}) = fmpz

doc"""
    base_ring(a::FlintIntegerRing)
> Returns `Union{}` as this ring is not dependent on another ring.
"""
base_ring(a::FlintIntegerRing) = Union{}

doc"""
    base_ring(a::fmpz)
> Returns `Union{}` as the parent ring is not dependent on another ring.
"""
base_ring(a::fmpz) = Union{}

isdomain_type(::Type{fmpz}) = true

################################################################################
#
#   Hashing
#
################################################################################

# Similar to hash for BigInt found in julia/base

function _fmpz_is_small(a::fmpz)
   return __fmpz_is_small(a.d)
end

function _fmpz_limbs(a::fmpz)
   return __fmpz_limbs(a.d)
end

function hash_integer(a::fmpz, h::UInt)
   return _hash_integer(a.d, h)
end

function hash(a::fmpz, h::UInt)
   return hash_integer(a, h)
end

function __fmpz_is_small(a::Int)
   return (unsigned(a) >> (Sys.WORD_SIZE - 2) != 1)
end

function __fmpz_limbs(a::Int)
   if __fmpz_is_small(a)
      return 0
   end
   b = unsafe_load(convert(Ref{Cint}, unsigned(a)<<2), 2)
   return b
end

function _hash_integer(a::Int, h::UInt)
   s = __fmpz_limbs(a)
   s == 0 && return Base.hash_integer(a, h)
   # get the pointer after the first two Cint
   d = convert(Ptr{Ref{UInt}}, unsigned(a) << 2) + 2*sizeof(Cint)
   p = unsafe_load(d)
   b = unsafe_load(p)
   h = xor(Base.hash_uint(xor(ifelse(s < 0, -b, b), h)), h)
   for k = 2:abs(s)
      h = xor(Base.hash_uint(xor(unsafe_load(p, k), h)), h)
   end
   return h
end

###############################################################################
#
#   Basic manipulation
#
###############################################################################

function deepcopy_internal(a::fmpz, dict::ObjectIdDict)
   z = fmpz()
   ccall((:fmpz_set, :libflint), Void, (Ref{fmpz}, Ref{fmpz}), z, a)
   return z
end

doc"""
    one(R::FlintIntegerRing)
> Return the integer $1$.
"""
one(R::FlintIntegerRing) = fmpz(1)

doc"""
    zero(R::FlintIntegerRing)
> Return the integer $1$.
"""
zero(R::FlintIntegerRing) = fmpz(0)

doc"""
    sign(a::fmpz)
> Returns the sign of $a$, i.e. $+1$, $0$ or $-1$.
"""
sign(a::fmpz) = Int(ccall((:fmpz_sgn, :libflint), Cint, (Ref{fmpz},), a))

doc"""
    fits(::Type{Int}, a::fmpz)
> Returns `true` if the given integer fits into an `Int`, otherwise returns
> `false`.
"""
fits(::Type{Int}, a::fmpz) = ccall((:fmpz_fits_si, :libflint), Bool,
                                   (Ref{fmpz},), a)

doc"""
    fits(::Type{UInt}, a::fmpz)
> Returns `true` if the given integer fits into a `UInt`, otherwise returns
> `false`.
"""
fits(::Type{UInt}, a::fmpz) = sign(a) < 0 ? false :
              ccall((:fmpz_abs_fits_ui, :libflint), Bool, (Ref{fmpz},), a)

doc"""
    size(a::fmpz)
> Returns the number of limbs required to store the absolute value of $a$.
"""
size(a::fmpz) = Int(ccall((:fmpz_size, :libflint), Cint, (Ref{fmpz},), a))

doc"""
    isunit(a::fmpz)
> Return `true` if the given integer is a unit, i.e. $\pm 1$, otherwise return
> `false`.
"""
isunit(a::fmpz) = ccall((:fmpz_is_pm1, :libflint), Bool, (Ref{fmpz},), a)

doc"""
    iszero(a::fmpz)
> Return `true` if the given integer is zero, otherwise return `false`.
"""
iszero(a::fmpz) = ccall((:fmpz_is_zero, :libflint), Bool, (Ref{fmpz},), a)

doc"""
    isone(a::fmpz)
> Return `true` if the given integer is one, otherwise return `false`.
"""
isone(a::fmpz) = ccall((:fmpz_is_one, :libflint), Bool, (Ref{fmpz},), a)

doc"""
    denominator(a::fmpz)
> Returns the denominator of $a$ thought of as a rational. Always returns $1$.
"""
function denominator(a::fmpz)
   return fmpz(1)
end

doc"""
    numerator(a::fmpz)
> Returns the numerator of $a$ thought of as a rational. Always returns $a$.
"""
function numerator(a::fmpz)
   return a
end

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

string(x::fmpz) = dec(x)

show(io::IO, x::fmpz) = print(io, string(x))

show(io::IO, a::FlintIntegerRing) = print(io, "Integer Ring")

needs_parentheses(x::fmpz) = false

isnegative(x::fmpz) = x < 0

show_minus_one(::Type{fmpz}) = false

###############################################################################
#
#   Canonicalisation
#
###############################################################################

canonical_unit(x::fmpz) = x < 0 ? fmpz(-1) : fmpz(1)

###############################################################################
#
#   Unary operators and functions, e.g. -fmpz(12), ~fmpz(12)
#
###############################################################################

function -(x::fmpz)
    z = fmpz()
    ccall((:__fmpz_neg, :libflint), Void, (Ref{fmpz}, Ref{fmpz}), z, x)
    return z
end

function ~(x::fmpz)
    z = fmpz()
    ccall((:fmpz_complement, :libflint), Void, (Ref{fmpz}, Ref{fmpz}), z, x)
    return z
end

function abs(x::fmpz)
    z = fmpz()
    ccall((:fmpz_abs, :libflint), Void, (Ref{fmpz}, Ref{fmpz}), z, x)
    return z
end

###############################################################################
#
#   Binary operators and functions
#
###############################################################################

# Metaprogram to define functions +, -, *, gcd, lcm,
#                                 &, |, $

for (fJ, fC) in ((:+, :add), (:-,:sub), (:*, :mul),
                 (:&, :and), (:|, :or), (:$, :xor))
    @eval begin
        function ($fJ)(x::fmpz, y::fmpz)
            z = fmpz()
            ccall(($(string(:fmpz_, fC)), :libflint), Void,
                  (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, y)
            return z
        end
    end
end

# Metaprogram to define functions fdiv, cdiv, tdiv, div, mod

for (fJ, fC) in ((:fdiv, :fdiv_q), (:cdiv, :cdiv_q), (:tdiv, :tdiv_q),
                 (:div, :tdiv_q))
    @eval begin
        function ($fJ)(x::fmpz, y::fmpz)
            iszero(y) && throw(DivideError())
            z = fmpz()
            ccall(($(string(:fmpz_, fC)), :libflint), Void,
                  (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, y)
            return z
        end
    end
end

function divexact(x::fmpz, y::fmpz)
    iszero(y) && throw(DivideError())
    z = fmpz()
    ccall((:fmpz_divexact, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, y)
    z
end

function rem(x::fmpz, c::fmpz)
    iszero(c) && throw(DivideError())
    q = fmpz()
    r = fmpz()
    ccall((:fmpz_tdiv_qr, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), q, r, x, abs(c))
    return r
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

function +(x::fmpz, c::Int)
    z = fmpz()
    if c >= 0
       ccall((:fmpz_add_ui, :libflint), Void,
             (Ref{fmpz}, Ref{fmpz}, Int), z, x, c)
    else
       ccall((:fmpz_sub_ui, :libflint), Void,
             (Ref{fmpz}, Ref{fmpz}, Int), z, x, -c)
    end
    return z
end

+(c::Int, x::fmpz) = x + c

function -(x::fmpz, c::Int)
    z = fmpz()
    if c >= 0
       ccall((:fmpz_sub_ui, :libflint), Void,
             (Ref{fmpz}, Ref{fmpz}, Int), z, x, c)
    else
       ccall((:fmpz_add_ui, :libflint), Void,
             (Ref{fmpz}, Ref{fmpz}, Int), z, x, -c)
    end
    return z
end

function -(c::Int, x::fmpz)
    z = fmpz()
    if c >= 0
       ccall((:fmpz_sub_ui, :libflint), Void,
             (Ref{fmpz}, Ref{fmpz}, Int), z, x, c)
    else
       ccall((:fmpz_add_ui, :libflint), Void,
             (Ref{fmpz}, Ref{fmpz}, Int), z, x, -c)
    end
    ccall((:__fmpz_neg, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}), z, z)
    return z
end

function *(x::fmpz, c::Int)
    z = fmpz()
    ccall((:fmpz_mul_si, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Int), z, x, c)
    return z
end

*(c::Int, x::fmpz) = x * c

+(a::fmpz, b::Integer) = a + fmpz(b)

+(a::Integer, b::fmpz) = fmpz(a) + b

-(a::fmpz, b::Integer) = a - fmpz(b)

-(a::Integer, b::fmpz) = fmpz(a) - b

*(a::fmpz, b::Integer) = a*fmpz(b)

*(a::Integer, b::fmpz) = fmpz(a)*b

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

function divexact(x::fmpz, y::Int)
    y == 0 && throw(DivideError())
    z = fmpz()
    ccall((:fmpz_divexact_si, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Int), z, x, y)
    z
end

divexact(x::fmpz, y::Integer) = divexact(x, fmpz(y))

divexact(x::Integer, y::fmpz) = divexact(fmpz(x), y)

###############################################################################
#
#   Ad hoc division
#
###############################################################################

function rem(x::fmpz, c::Int)
    c == 0 && throw(DivideError())
    r = ccall((:fmpz_tdiv_ui, :libflint), Int, (Ref{fmpz}, Int), x, abs(c))
    return sign(x) > 0 ? r : -r
end

function tdivpow2(x::fmpz, c::Int)
    c < 0 && throw(DomainError())
    z = fmpz()
    ccall((:fmpz_tdiv_q_2exp, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Int), z, x, c)
    return z
end

function fdivpow2(x::fmpz, c::Int)
    c < 0 && throw(DomainError())
    z = fmpz()
    ccall((:fmpz_fdiv_q_2exp, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Int), z, x, c)
    return z
end

function fmodpow2(x::fmpz, c::Int)
    c < 0 && throw(DomainError())
    z = fmpz()
    ccall((:fmpz_fdiv_r_2exp, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Int), z, x, c)
    return z
end

function cdivpow2(x::fmpz, c::Int)
    c < 0 && throw(DomainError())
    z = fmpz()
    ccall((:fmpz_cdiv_q_2exp, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Int), z, x, c)
    return z
end

function div(x::fmpz, c::Int)
    c == 0 && throw(DivideError())
    z = fmpz()
    ccall((:fmpz_tdiv_q_si, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Int), z, x, c)
    return z
end

function tdiv(x::fmpz, c::Int)
    c == 0 && throw(DivideError())
    z = fmpz()
    ccall((:fmpz_tdiv_q_si, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Int), z, x, c)
    return z
end

function fdiv(x::fmpz, c::Int)
    c == 0 && throw(DivideError())
    z = fmpz()
    ccall((:fmpz_fdiv_q_si, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Int), z, x, c)
    return z
end

function cdiv(x::fmpz, c::Int)
    c == 0 && throw(DivideError())
    z = fmpz()
    ccall((:fmpz_cdiv_q_si, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Int), z, x, c)
    return z
end

###############################################################################
#
#   Division with remainder
#
###############################################################################

function divrem(x::fmpz, y::fmpz)
    iszero(y) && throw(DivideError())
    z1 = fmpz()
    z2 = fmpz()
    ccall((:fmpz_tdiv_qr, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z1, z2, x, y)
    z1, z2
end

function tdivrem(x::fmpz, y::fmpz)
    iszero(y) && throw(DivideError())
    z1 = fmpz()
    z2 = fmpz()
    ccall((:fmpz_tdiv_qr, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z1, z2, x, y)
    z1, z2
end

function fdivrem(x::fmpz, y::fmpz)
    iszero(y) && throw(DivideError())
    z1 = fmpz()
    z2 = fmpz()
    ccall((:fmpz_fdiv_qr, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z1, z2, x, y)
    z1, z2
end

###############################################################################
#
#   Powering
#
###############################################################################

function ^(x::fmpz, y::Int)
    if y < 0; throw(DomainError()); end
    if isone(x); return x; end
    if x == -1; return isodd(y) ? x : -x; end
    if y > typemax(UInt); throw(DomainError()); end
    if y == 0; return one(FlintZZ); end
    if y == 1; return x; end
    z = fmpz()
    ccall((:fmpz_pow_ui, :libflint), Void, (Ref{fmpz}, Ref{fmpz}, UInt), z, x, UInt(y))
    return z
end

###############################################################################
#
#   Comparison
#
###############################################################################

function cmp(x::fmpz, y::fmpz)
    Int(ccall((:fmpz_cmp, :libflint), Cint,
              (Ref{fmpz}, Ref{fmpz}), x, y))
end

==(x::fmpz, y::fmpz) = cmp(x,y) == 0

<=(x::fmpz, y::fmpz) = cmp(x,y) <= 0

>=(x::fmpz, y::fmpz) = cmp(x,y) >= 0

<(x::fmpz, y::fmpz) = cmp(x,y) < 0

>(x::fmpz, y::fmpz) = cmp(x,y) > 0

function cmpabs(x::fmpz, y::fmpz)
    Int(ccall((:fmpz_cmpabs, :libflint), Cint,
              (Ref{fmpz}, Ref{fmpz}), x, y))
end

isless(x::fmpz, y::fmpz) = x < y

###############################################################################
#
#   Ad hoc comparison
#
###############################################################################

function cmp(x::fmpz, y::Int)
    Int(ccall((:fmpz_cmp_si, :libflint), Cint, (Ref{fmpz}, Int), x, y))
end

==(x::fmpz, y::Int) = cmp(x,y) == 0

<=(x::fmpz, y::Int) = cmp(x,y) <= 0

>=(x::fmpz, y::Int) = cmp(x,y) >= 0

<(x::fmpz, y::Int) = cmp(x,y) < 0

>(x::fmpz, y::Int) = cmp(x,y) > 0

==(x::Int, y::fmpz) = cmp(y,x) == 0

<=(x::Int, y::fmpz) = cmp(y,x) >= 0

>=(x::Int, y::fmpz) = cmp(y,x) <= 0

<(x::Int, y::fmpz) = cmp(y,x) > 0

>(x::Int, y::fmpz) = cmp(y,x) < 0

function cmp(x::fmpz, y::UInt)
    Int(ccall((:fmpz_cmp_ui, :libflint), Cint, (Ref{fmpz}, UInt), x, y))
end

==(x::fmpz, y::UInt) = cmp(x,y) == 0

<=(x::fmpz, y::UInt) = cmp(x,y) <= 0

>=(x::fmpz, y::UInt) = cmp(x,y) >= 0

<(x::fmpz, y::UInt) = cmp(x,y) < 0

>(x::fmpz, y::UInt) = cmp(x,y) > 0

==(x::UInt, y::fmpz) = cmp(y,x) == 0

<=(x::UInt, y::fmpz) = cmp(y,x) >= 0

>=(x::UInt, y::fmpz) = cmp(y,x) <= 0

<(x::UInt, y::fmpz) = cmp(y,x) > 0

>(x::UInt, y::fmpz) = cmp(y,x) < 0

###############################################################################
#
#   Shifting
#
###############################################################################

doc"""
    <<(x::fmpz, c::Int)
> Return $2^cx$ where $c \geq 0$.
"""
function <<(x::fmpz, c::Int)
    c < 0 && throw(DomainError())
    c == 0 && return x
    z = fmpz()
    ccall((:fmpz_mul_2exp, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Int), z, x, c)
    return z
end

doc"""
    >>(x::fmpz, c::Int)
> Return $x/2^c$, discarding any remainder, where $c \geq 0$.
"""
function >>(x::fmpz, c::Int)
    c < 0 && throw(DomainError())
    c == 0 && return x
    z = fmpz()
    ccall((:fmpz_fdiv_q_2exp, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Int), z, x, c)
    return z
end

###############################################################################
#
#   Modular arithmetic
#
###############################################################################

doc"""
    mod(x::fmpz, y::fmpz)
> Return the remainder after division of $x$ by $y$. The remainder will be the
> least nonnegative remainder.
"""
function mod(x::fmpz, y::fmpz)
   z = fmpz()
   ccall((:fmpz_mod, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, y)
   return z
end

doc"""
    mod(x::fmpz, y::Int)
> Return the remainder after division of $x$ by $y$. The remainder will be the
> least nonnegative remainder.
"""
function mod(x::fmpz, c::Int)
    c == 0 && throw(DivideError())
    if c > 0
        return ccall((:fmpz_fdiv_ui, :libflint), Int, (Ref{fmpz}, Int), x, c)
    else
        r = ccall((:fmpz_fdiv_ui, :libflint), Int, (Ref{fmpz}, Int), x, -c)
        return r == 0 ? 0 : r + c
    end
end

doc"""
    powmod(x::fmpz, p::fmpz, m::fmpz)
> Return $x^p (\mod m)$. The remainder will be in the range $[0, m)$
"""
function powmod(x::fmpz, p::fmpz, m::fmpz)
    m <= 0 && throw(DomainError())
    if p < 0
       x = invmod(x, m)
       p = -p
    end
    r = fmpz()
    ccall((:fmpz_powm, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Ref{fmpz}),
          r, x, p, m)
    return r
end

doc"""
    powmod(x::fmpz, p::Int, m::fmpz)
> Return $x^p (\mod m)$. The remainder will be in the range $[0, m)$
"""
function powmod(x::fmpz, p::Int, m::fmpz)
    m <= 0 && throw(DomainError())
    if p < 0
       x = invmod(x, m)
       p = -p
    end
    r = fmpz()
    ccall((:fmpz_powm_ui, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Int, Ref{fmpz}),
          r, x, p, m)
    return r
end

doc"""
    invmod(x::fmpz, m::fmpz)
> Return $x^{-1} (\mod m)$. The remainder will be in the range $[0, m)$
"""
function invmod(x::fmpz, m::fmpz)
    m <= 0 && throw(DomainError())
    z = fmpz()
    if isone(m)
        return fmpz(0)
    end
    if ccall((:fmpz_invmod, :libflint), Cint,
             (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, m) == 0
       error("Impossible inverse in invmod")
    end
    return z
end

doc"""
    sqrtmod(x::fmpz, m::fmpz)
> Return a square root of $x (\mod m)$ if one exists. The remainder will be in
> the range $[0, m)$. We require that $m$ is prime, otherwise the algorithm may
> not terminate.
"""
function sqrtmod(x::fmpz, m::fmpz)
    m <= 0 && throw(DomainError())
    z = fmpz()
    if (ccall((:fmpz_sqrtmod, :libflint), Cint,
              (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, m) == 0)
        error("no square root exists")
    end
    return z
end

doc"""
    crt(r1::fmpz, m1::fmpz, r2::fmpz, m2::fmpz, signed=false)
> Find $r$ such that $r \equiv r_1 (\mod m_1)$ and $r \equiv r_2 (\mod m_2)$.
> If `signed = true`, $r$ will be in the range $-m_1m_2/2 < r \leq m_1m_2/2$.
> If `signed = false` the value will be in the range $0 \leq r < m_1m_2$.
"""
function crt(r1::fmpz, m1::fmpz, r2::fmpz, m2::fmpz, signed=false)
   z = fmpz()
   ccall((:fmpz_CRT, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Cint),
          z, r1, m1, r2, m2, signed)
   return z
end

doc"""
    crt(r1::fmpz, m1::fmpz, r2::Int, m2::Int, signed=false)
> Find $r$ such that $r \equiv r_1 (\mod m_1)$ and $r \equiv r_2 (\mod m_2)$.
> If `signed = true`, $r$ will be in the range $-m_1m_2/2 < r \leq m_1m_2/2$.
> If `signed = false` the value will be in the range $0 \leq r < m_1m_2$.
"""
function crt(r1::fmpz, m1::fmpz, r2::Int, m2::Int, signed=false)
   z = fmpz()
   r2 < 0 && throw(DomainError())
   m2 < 0 && throw(DomainError())
   ccall((:fmpz_CRT_ui, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Int, Int, Cint),
          z, r1, m1, r2, m2, signed)
   return z
end

###############################################################################
#
#   Integer logarithm
#
###############################################################################

doc"""
    flog(x::fmpz, c::fmpz)
> Return the floor of the logarithm of $x$ to base $c$.
"""
function flog(x::fmpz, c::fmpz)
    c <= 0 && throw(DomainError())
    x <= 0 && throw(DomainError())
    return ccall((:fmpz_flog, :libflint), Int,
                 (Ref{fmpz}, Ref{fmpz}), x, c)
end

doc"""
    clog(x::fmpz, c::fmpz)
> Return the ceiling of the logarithm of $x$ to base $c$.
"""
function clog(x::fmpz, c::fmpz)
    c <= 0 && throw(DomainError())
    x <= 0 && throw(DomainError())
    return ccall((:fmpz_clog, :libflint), Int,
                 (Ref{fmpz}, Ref{fmpz}), x, c)
end

doc"""
    flog(x::fmpz, c::Int)
> Return the floor of the logarithm of $x$ to base $c$.
"""
function flog(x::fmpz, c::Int)
    c <= 0 && throw(DomainError())
    return ccall((:fmpz_flog_ui, :libflint), Int,
                 (Ref{fmpz}, Int), x, c)
end

doc"""
    clog(x::fmpz, c::Int)
> Return the ceiling of the logarithm of $x$ to base $c$.
"""
function clog(x::fmpz, c::Int)
    c <= 0 && throw(DomainError())
    return ccall((:fmpz_clog_ui, :libflint), Int,
                 (Ref{fmpz}, Int), x, c)
end

###############################################################################
#
#   GCD and LCM
#
###############################################################################

doc"""
    gcd(x::fmpz, y::fmpz)
> Return the greatest common divisor of $x$ and $y$. The returned result will
> always be nonnegative and will be zero iff $x$ and $y$ are zero.
"""
function gcd(x::fmpz, y::fmpz)
   z = fmpz()
   ccall((:fmpz_gcd, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, y)
   return z
end

doc"""
    gcd(x::Array{fmpz, 1})
> Return the greatest common divisor of the elements of $x$. The returned
> result will always be nonnegative and will be zero iff all elements of $x$
> are zero.
"""
function gcd(x::Array{fmpz, 1})
   if length(x) == 0
      error("Array must not be empty")
   elseif length(x) == 1
      return x[1]
   end

   z = fmpz()
   ccall((:fmpz_gcd, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x[1], x[2])

   for i in 3:length(x)
      ccall((:fmpz_gcd, :libflint), Void,
            (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, z, x[i])
      if isone(z)
         return z
      end
   end

   return z
end

doc"""
    lcm(x::fmpz, y::fmpz)
> Return the least common multiple of $x$ and $y$. The returned result will
> always be nonnegative and will be zero iff $x$ and $y$ are zero.
"""
function lcm(x::fmpz, y::fmpz)
   z = fmpz()
   ccall((:fmpz_lcm, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, y)
   return z
end

doc"""
    lcm(x::Array{fmpz, 1})
> Return the least common multiple of the elements of $x$. The returned result
> will always be nonnegative and will be zero iff the elements of $x$ are zero.
"""
function lcm(x::Array{fmpz, 1})
   if length(x) == 0
      error("Array must not be empty")
   elseif length(x) == 1
      return x[1]
   end

   z = fmpz()
   ccall((:fmpz_lcm, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x[1], x[2])

   for i in 3:length(x)
      ccall((:fmpz_lcm, :libflint), Void,
            (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, z, x[i])
   end

   return z
end

###############################################################################
#
#   Extended GCD
#
###############################################################################

doc"""
    gcdx(a::fmpz, b::fmpz)
> Return a tuple $g, s, t$ such that $g$ is the greatest common divisor of $a$
> and $b$ and integers $s$ and $t$ such that $g = as + bt$.
"""
function gcdx(a::fmpz, b::fmpz)
    if b == 0 # shortcut this to ensure consistent results with gcdx(a,b)
        return a < 0 ? (-a, -one(FlintZZ), zero(FlintZZ)) : (a, one(FlintZZ), zero(FlintZZ))
    end
    g = fmpz()
    s = fmpz()
    t = fmpz()
    ccall((:fmpz_xgcd, :libflint), Void,
        (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Ref{fmpz}),
        g, s, t, a, b)
    g, s, t
end

doc"""
    gcdinv(a::fmpz, b::fmpz)
> Return a tuple $g, s$ where $g$ is the greatest common divisor of $a$ and
> $b$ and where $s$ is the inverse of $a$ modulo $b$ if $g = 1$. This function
> can be used to detect impossible inverses, i.e. where $a$ and $b$ are not
> coprime, and to yield the common factor of $a$ and $b$ if they are not
> coprime. We require $b \geq a \geq 0$.
"""
function gcdinv(a::fmpz, b::fmpz)
   a < 0 && throw(DomainError())
   b < a && throw(DomainError())
   g = fmpz()
   s = fmpz()
   ccall((:fmpz_gcdinv, :libflint), Void,
        (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Ref{fmpz}),
        g, s, a, b)
   return g, s
end

###############################################################################
#
#   Roots
#
###############################################################################

doc"""
    isqrt(x::fmpz)
> Return the floor of the square root of $x$.
"""
function isqrt(x::fmpz)
    x < 0 && throw(DomainError())
    z = fmpz()
    ccall((:fmpz_sqrt, :libflint), Void, (Ref{fmpz}, Ref{fmpz}), z, x)
    return z
end

doc"""
    isqrtrem(x::fmpz)
> Return a tuple $s, r$ consisting of the floor $s$ of the square root of $x$
> and the remainder $r$, i.e. such that $x = s^2 + r$. We require $x \geq 0$.
"""
function isqrtrem(x::fmpz)
    x < 0 && throw(DomainError())
    s = fmpz()
    r = fmpz()
    ccall((:fmpz_sqrtrem, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), s, r, x)
    return s, r
end

doc"""
    root(x::fmpz, n::Int)
> Return the floor of the $n$-the root of $x$. We require $n > 0$ and that
> $x \geq 0$ if $n$ is even.
"""
function root(x::fmpz, n::Int)
   x < 0 && iseven(n) && throw(DomainError())
   n <= 0 && throw(DomainError())
   z = fmpz()
   ccall((:fmpz_root, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz}, Int), z, x, n)
   return z
end

###############################################################################
#
#   Factorization
#
###############################################################################

function _factor(a::fmpz)
   # This is a hack around https://github.com/JuliaLang/julia/issues/19963
   # Remove this once julia 6.0 is required
   if a == 1 || a == -1
     return Dict{fmpz, Int}(), a
   end

   F = fmpz_factor()
   ccall((:fmpz_factor, :libflint), Void, (Ref{fmpz_factor}, Ref{fmpz}), F, a)
   res = Dict{fmpz, Int}()
   for i in 1:F.num
     z = fmpz()
     ccall((:fmpz_factor_get_fmpz, :libflint), Void,
           (Ref{fmpz}, Ref{fmpz_factor}, Int), z, F, i - 1)
     res[z] = unsafe_load(F.exp, i)
   end
   return res, canonical_unit(a)
end

function factor(a::fmpz)
   fac, z = _factor(a)
   return Fac(z, fac)
end

###############################################################################
#
#   Number theoretic/combinatorial
#
###############################################################################

doc"""
    divisible(x::fmpz, y::fmpz)
> Return `true` if $x$ is divisible by $y$, otherwise return `false`. We
> require $x \neq 0$.
"""
function divisible(x::fmpz, y::fmpz)
   iszero(y) && throw(DivideError())
   Bool(ccall((:fmpz_divisible, :libflint), Cint,
              (Ref{fmpz}, Ref{fmpz}), x, y))
end

doc"""
    divisible(x::fmpz, y::Int)
> Return `true` if $x$ is divisible by $y$, otherwise return `false`. We
> require $x \neq 0$.
"""
function divisible(x::fmpz, y::Int)
   y == 0 && throw(DivideError())
   Bool(ccall((:fmpz_divisible_si, :libflint), Cint,
              (Ref{fmpz}, Int), x, y))
end

doc"""
    issquare(x::fmpz)
> Return `true` if $x$ is a square, otherwise return `false`.
"""
issquare(x::fmpz) = Bool(ccall((:fmpz_is_square, :libflint), Cint,
                               (Ref{fmpz},), x))

doc"""
    is_prime(x::UInt)
> Return `true` if $x$ is a prime number, otherwise return `false`.
"""
is_prime(x::UInt) = Bool(ccall((:n_is_prime, :libflint), Cint, (UInt,), x))

doc"""
    isprime(x::fmpz)
> Return `true` if $x$ is a prime number, otherwise return `false`.
"""
# flint's fmpz_is_prime doesn't work yet
isprime(x::fmpz) = Bool(ccall((:fmpz_is_probabprime, :libflint), Cint,
                              (Ref{fmpz},), x))

doc"""
    isprobabprime(x::fmpz)
> Return `true` if $x$ is a very probably a prime number, otherwise return
> `false`. No counterexamples are known to this test, but it is conjectured
> that infinitely many exist.
"""
isprobabprime(x::fmpz) = Bool(ccall((:fmpz_is_probabprime, :libflint), Cint,
                                    (Ref{fmpz},), x))

doc"""
    remove(x::fmpz, y::fmpz)
> Return the tuple $n, z$ such that $x = y^nz$ where $y$ and $z$ are coprime.
"""
function remove(x::fmpz, y::fmpz)
   iszero(y) && throw(DivideError())
   z = fmpz()
   num = ccall((:fmpz_remove, :libflint), Int,
               (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, y)
   return num, z
end

remove(x::fmpz, y::Integer) = remove(x, fmpz(y))

remove(x::Integer, y::fmpz) = remove(fmpz(x), y)

remove(x::Integer, y::Integer) = remove(fmpz(x), fmpz(y))

doc"""
    valuation(x::fmpz, y::fmpz)
> Return the largest $n$ such that $y^n$ divides $x$.
"""
function valuation(x::fmpz, y::fmpz)
   n, _ = remove(x, y)
   return n
end

valuation(x::fmpz, y::Integer) = valuation(x, fmpz(y))

valuation(x::Integer, y::fmpz) = valuation(fmpz(x), y)

valuation(x::Integer, y::Integer) = valuation(fmpz(x), fmpz(y))

doc"""
    divisor_lenstra(n::fmpz, r::fmpz, m::fmpz)
> If $n$ has a factor which lies in the residue class $r (\mod m)$ for
> $0 < r < m < n$, this function returns such a factor. Otherwise it returns
> $0$. This is only efficient if $m$ is at least the cube root of $n$. We
> require gcd$(r, m) = 1$ and this condition is not checked.
"""
function divisor_lenstra(n::fmpz, r::fmpz, m::fmpz)
   r <= 0 && throw(DomainError())
   m <= r && throw(DomainError())
   n <= m && throw(DomainError())
   z = fmpz()
   if !Bool(ccall((:fmpz_divisor_in_residue_class_lenstra, :libflint),
       Cint, (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, n, r, m))
      z = 0
   end
   return z
end

doc"""
    fac(x::Int)
> Return the factorial of $x$, i.e. $x! = 1.2.3\ldots x$. We require
> $x \geq 0$.
"""
function fac(x::Int)
    x < 0 && throw(DomainError())
    z = fmpz()
    ccall((:fmpz_fac_ui, :libflint), Void, (Ref{fmpz}, UInt), z, x)
    return z
end

doc"""
    risingfac(x::fmpz, y::Int)
> Return the rising factorial of $x$, i.e. $x(x + 1)(x + 2)\ldots (x + n - 1)$.
> If $n < 0$ we throw a `DomainError()`.
"""
function risingfac(x::fmpz, y::Int)
    y < 0 && throw(DomainError())
    z = fmpz()
    ccall((:fmpz_rfac_ui, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, UInt), z, x, y)
    return z
end

doc"""
    risingfac(x::Int, y::Int)
> Return the rising factorial of $x$, i.e. $x(x + 1)(x + 2)\ldots (x + n - 1)$.
> If $n < 0$ we throw a `DomainError()`.
"""
function risingfac(x::Int, y::Int)
    y < 0 && throw(DomainError())
    z = fmpz()
    if x < 0
       if y <= -x # we don't pass zero
          z = isodd(y) ? -risingfac(-x - y + 1, y) :
                          risingfac(-x - y + 1, y)
       end
    else
       ccall((:fmpz_rfac_uiui, :libflint), Void,
             (Ref{fmpz}, UInt, UInt), z, x, y)
    end
    return z
end

doc"""
    primorial(x::Int)
>  Return the primorial of $n$, i.e. the product of all primes less than or
> equal to $n$. If $n < 0$ we throw a `DomainError()`.
"""
function primorial(x::Int)
    x < 0 && throw(DomainError())
    z = fmpz()
    ccall((:fmpz_primorial, :libflint), Void,
          (Ref{fmpz}, UInt), z, x)
    return z
end

doc"""
    fib(x::Int)
>  Return the $n$-th Fibonacci number $F_n$. We define $F_1 = 1$, $F_2 = 1$ and
> $F_{i + 1} = F_i + F_{i - 1}$ for all $i > 2$. We require $n \geq 0$. For
> convenience, we define $F_0 = 0$.
"""
function fib(x::Int)
    x < 0 && throw(DomainError())
    z = fmpz()
    ccall((:fmpz_fib_ui, :libflint), Void,
          (Ref{fmpz}, UInt), z, x)
    return z
end

doc"""
    bell(x::Int)
> Return the Bell number $B_n$.
"""
function bell(x::Int)
    x < 0 && throw(DomainError())
    z = fmpz()
    ccall((:arith_bell_number, :libflint), Void,
          (Ref{fmpz}, UInt), z, x)
    return z
end

doc"""
    binom(n::Int, k::Int)
> Return the binomial coefficient $\frac{n!}{(n - k)!k!}$. If $n, k < 0$ or
> $k > n$ we return $0$.
"""
function binom(n::Int, k::Int)
    n < 0 && return fmpz(0)
    k < 0 && return fmpz(0)
    z = fmpz()
    ccall((:fmpz_bin_uiui, :libflint), Void,
          (Ref{fmpz}, UInt, UInt), z, n, k)
    return z
end

doc"""
    moebiusmu(x::fmpz)
> Returns the Moebius mu function of $x$ as an \code{Int}. The value
> returned is either $-1$, $0$ or $1$. If $x < 0$ we throw a `DomainError()`.
"""
function moebiusmu(x::fmpz)
   x < 0 && throw(DomainError())
   return Int(ccall((:fmpz_moebius_mu, :libflint), Cint,
                    (Ref{fmpz},), x))
end

doc"""
    jacobi(x::fmpz, y::fmpz)
> Return the value of the Jacobi symbol $\left(\frac{x}{y}\right)$. If
> $y \leq x$ or $x < 0$, we throw a `DomainError()`.
"""
function jacobi(x::fmpz, y::fmpz)
   y <= x && throw(DomainError())
   x < 0 && throw(DomainError())
   return Int(ccall((:fmpz_jacobi, :libflint), Cint,
                    (Ref{fmpz}, Ref{fmpz}), x, y))
end

doc"""
    sigma(x::fmpz, y::Int)
> Return the value of the sigma function, i.e. $\sum_{0 < d \;| x} d^y$. If
> $y < 0$ we throw a `DomainError()`.
"""
function sigma(x::fmpz, y::Int)
   y < 0 && throw(DomainError())
   z = fmpz()
   ccall((:fmpz_divisor_sigma, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz}, Int), z, x, y)
   return z
end

doc"""
    eulerphi(x::fmpz)
> Return the value of the Euler phi function at $x$, i.e. the number of
> positive integers less than $x$ that are coprime with $x$.
"""
function eulerphi(x::fmpz)
   x < 0 && throw(DomainError())
   z = fmpz()
   ccall((:fmpz_euler_phi, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz}), z, x)
   return z
end

doc"""
    numpart(x::Int)
> Return the number of partitions of $x$. This function is not available on
> Windows 64.
"""
function numpart(x::Int)
   if (is_windows() ? true : false) && Int == Int64
      error("not yet supported on win64")
   end
   x < 0 && throw(DomainError())
   z = fmpz()
   ccall((:partitions_fmpz_ui, :libarb), Void,
         (Ref{fmpz}, UInt), z, x)
   return z
end

doc"""
    numpart(x::fmpz)
> Return the number of partitions of $x$. This function is not available on
> Windows 64.
"""
function numpart(x::fmpz)
   if (is_windows() ? true : false) && Int == Int64
      error("not yet supported on win64")
   end
   x < 0 && throw(DomainError())
   z = fmpz()
   ccall((:partitions_fmpz_fmpz, :libarb), Void,
         (Ref{fmpz}, Ref{fmpz}, Int), z, x, 0)
   return z
end

###############################################################################
#
#   Number bases/digits
#
###############################################################################

doc"""
    bin(n::fmpz)
> Return $n$ as a binary string.
"""
bin(n::fmpz) = base(n, 2)

doc"""
    oct(n::fmpz)
> Return $n$ as a octal string.
"""
oct(n::fmpz) = base(n, 8)

doc"""
    dec(n::fmpz)
> Return $n$ as a decimal string.
"""
dec(n::fmpz) = base(n, 10)

doc"""
    hex(n::fmpz) = base(n, 16)
> Return $n$ as a hexadecimal string.
"""
hex(n::fmpz) = base(n, 16)

doc"""
    base(n::fmpz, b::Integer)
> Return $n$ as a string in base $b$. We require $2 \leq b \leq 62$.
"""
function base(n::fmpz, b::Integer)
    2 <= b <= 62 || error("invalid base: $b")
    p = ccall((:fmpz_get_str,:libflint), Ref{UInt8},
              (Ref{UInt8}, Cint, Ref{fmpz}), C_NULL, b, n)
    s = unsafe_string(p)
    ccall((:flint_free, :libflint), Void, (Ref{UInt8},), p)
    return s
end

function ndigits_internal(x::fmpz, b::Integer = 10)
    # fmpz_sizeinbase might return an answer 1 too big
    n = Int(ccall((:fmpz_sizeinbase, :libflint), UInt,
                  (Ref{fmpz}, Int32), x, b))
    abs(x) < fmpz(b)^(n - 1) ? n - 1 : n
end

doc"""
    ndigits(x::fmpz, b::Integer = 10)
> Return the number of digits of $x$ in the base $b$ (default is $b = 10$).
"""
ndigits(x::fmpz, b::Integer = 10) = iszero(x) ? 1 : ndigits_internal(x, b)

doc"""
    nbits(x::fmpz)
> Return the number of binary bits of $x$. We return zero if $x = 0$.
"""
nbits(x::fmpz) = iszero(x) ? 0 : Int(ccall((:fmpz_sizeinbase, :libflint), UInt,
                  (Ref{fmpz}, Int32), x, 2))  # docu states: always correct
                                #if base is power of 2

###############################################################################
#
#   Bit fiddling
#
###############################################################################

doc"""
    popcount(x::fmpz)
> Return the number of ones in the binary representation of $x$.
"""
popcount(x::fmpz) = Int(ccall((:fmpz_popcnt, :libflint), UInt,
                              (Ref{fmpz},), x))

doc"""
    prevpow2(x::fmpz)
> Return the previous power of $2$ up to including $x$.
"""
prevpow2(x::fmpz) = x < 0 ? -prevpow2(-x) :
                            (x <= 2 ? x : one(FlintZZ) << (ndigits(x, 2) - 1))

doc"""
    nextpow2(x::fmpz)
> Return the next power of $2$ that is at least $x$.
"""
nextpow2(x::fmpz) = x < 0 ? -nextpow2(-x) :
                            (x <= 2 ? x : one(FlintZZ) << ndigits(x - 1, 2))

doc"""
    trailing_zeros(x::fmpz)
> Count the trailing zeros in the binary representation of $x$.
"""
trailing_zeros(x::fmpz) = ccall((:fmpz_val2, :libflint), Int,
                                (Ref{fmpz},), x)

###############################################################################
#
#   Bitwise operations (unsafe)
#
###############################################################################

doc"""
    clrbit!(x::fmpz, c::Int)
> Clear bit $c$ of $x$, where the least significant bit is the $0$-th bit. Note
> that this function modifies its input in-place.
"""
function clrbit!(x::fmpz, c::Int)
    c < 0 && throw(DomainError())
    ccall((:fmpz_clrbit, :libflint), Void, (Ref{fmpz}, Int), x, c)
end

doc"""
    setbit!(x::fmpz, c::Int)
> Set bit $c$ of $x$, where the least significant bit is the $0$-th bit. Note
> that this function modifies its input in-place.
"""
function setbit!(x::fmpz, c::Int)
    c < 0 && throw(DomainError())
    ccall((:fmpz_setbit, :libflint), Void, (Ref{fmpz}, Int), x, c)
end

doc"""
    combit!(x::fmpz, c::Int)
> Complement bit $c$ of $x$, where the least significant bit is the $0$-th bit.
> Note that this function modifies its input in-place.
"""
function combit!(x::fmpz, c::Int)
    c < 0 && throw(DomainError())
    ccall((:fmpz_combit, :libflint), Void, (Ref{fmpz}, Int), x, c)
end

###############################################################################
#
#   Unsafe operators
#
###############################################################################

function mul!(z::fmpz, x::fmpz, y::fmpz)
   ccall((:fmpz_mul, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, y)
   return z
end

function addmul!(z::fmpz, x::fmpz, y::fmpz, c::fmpz)
   ccall((:fmpz_addmul, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, y)
   return z
end

function addeq!(z::fmpz, x::fmpz)
   ccall((:fmpz_add, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, z, x)
   return z
end

function add!(z::fmpz, x::fmpz, y::fmpz)
   ccall((:fmpz_add, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, y)
   return z
end

function zero!(z::fmpz)
   ccall((:fmpz_zero, :libflint), Void,
         (Ref{fmpz},), z)
   return z
end

###############################################################################
#
#   Parent object overloads
#
###############################################################################

(::FlintIntegerRing)() = fmpz()

(::FlintIntegerRing)(a::Integer) = fmpz(a)

(::FlintIntegerRing)(a::AbstractString) = fmpz(a)

(::FlintIntegerRing)(a::fmpz) = a

(::FlintIntegerRing)(a::Float64) = fmpz(a)

(::FlintIntegerRing)(a::Float32) = fmpz(Float64(a))

(::FlintIntegerRing)(a::Float16) = fmpz(Float64(a))

(::FlintIntegerRing)(a::BigFloat) = fmpz(BigInt(a))

###############################################################################
#
#   String parser
#
###############################################################################

function parse(::Type{fmpz}, s::String, base::Int = 10)
    s = string(s)
    sgn = s[1] == '-' ? -1 : 1
    i = 1 + (sgn == -1)
    z = fmpz()
    err = ccall((:fmpz_set_str, :libflint),
               Int32, (Ref{fmpz}, Ref{UInt8}, Int32),
               z, string(SubString(s, i)), base)
    err == 0 || error("Invalid big integer: $(repr(s))")
    return sgn < 0 ? -z : z
end

###############################################################################
#
#   Random generation
#
###############################################################################

function rand(R::FlintIntegerRing, n::UnitRange{Int})
   return R(rand(n))
end

###############################################################################
#
#   Constructors
#
###############################################################################

fmpz(s::AbstractString) = parse(fmpz, s)

fmpz(z::Integer) = fmpz(BigInt(z))

fmpz(z::Float16) = fmpz(Float64(z))

fmpz(z::Float32) = fmpz(Float64(z))

fmpz(z::BigFloat) = fmpz(BigInt(z))

###############################################################################
#
#   Conversions and promotions
#
###############################################################################

convert(::Type{fmpz}, a::Integer) = fmpz(a)

function convert(::Type{BigInt}, a::fmpz)
   r = BigInt()
   ccall((:fmpz_get_mpz, :libflint), Void, (Ref{BigInt}, Ref{fmpz}), r, a)
   return r
end

function convert(::Type{Int}, a::fmpz)
   (a > typemax(Int) || a < typemin(Int)) && throw(InexactError())
   return ccall((:fmpz_get_si, :libflint), Int, (Ref{fmpz},), a)
end

function convert(::Type{UInt}, a::fmpz)
   (a > typemax(UInt) || a < 0) && throw(InexactError())
   return ccall((:fmpz_get_ui, :libflint), UInt, (Ref{fmpz}, ), a)
end

function convert(::Type{Float64}, n::fmpz)
    # rounds to zero
    ccall((:fmpz_get_d, :libflint), Float64, (Ref{fmpz},), n)
end

convert(::Type{Float32}, n::fmpz) = Float32(Float64(n))

convert(::Type{Float16}, n::fmpz) = Float16(Float64(n))

convert(::Type{BigFloat}, n::fmpz) = BigFloat(BigInt(n))

Base.promote_rule(::Type{fmpz}, ::Type{T}) where {T <: Integer} = fmpz

promote_rule(::Type{fmpz}, ::Type{T}) where {T <: Integer} = fmpz


text:
module PrecessionNutation

using AstronomicalTime
using ERFA
using StaticArrays
using MuladdMacro

import EarthOrientation: precession_nutation00

include("helper.jl")
include("tables.jl")

export precession_nutation_erfa, precession_nutation_06

function precession_nutation_erfa(ep::TTEpoch)
    x, y = ERFA.xy06(julian1(ep), julian2(ep))
    s = ERFA.s06(julian1(ep), julian2(ep), x, y)
    x, y, s
end

### IAU 2006 precession and IAU 2000A_R06 nutation
function precession_nutation_IAU2006_IAU2000A_R06(ep_tt::TTEpoch; revision = 0)
    tt_jc = jc(ep_tt)

    if revision === 0
        # P03
        X = @evalpoly(tt_jc, -16_617.0, +2_004_191_898.00,    -429_782.90, -198_618.34,     +7.578,  +5.928_5) # μas
        Y = @evalpoly(tt_jc,  -6_951.0,        -25_896.00, -22_407_274.70,   +1_900.59, +1_112.526,  +0.135_8) # μas
        s = @evalpoly(tt_jc,     +94.0,         +3_808.65,        -122.68,  -72_574.11,    +27.980, +15.620_0) # μas

    elseif revision === 1
        # P03_rev1
        X = @evalpoly(tt_jc, -16_617.0, +2_004_191_804.00,    -429_755.80, -198_618.29,     +7.575,  +5.928_5) # μas
        Y = @evalpoly(tt_jc,  -6_951.0,        -24_867.00, -22_407_272.70,   +1_900.26, +1_112.525,  +0.135_8) # μas

        # TODO: confirm values for s
        s = @evalpoly(tt_jc,     +94.0,         +3_808.65,        -122.68,  -72_574.11,    +27.980, +15.620_0) # μas

    elseif revision === 2
        # P03_rev2
        X = @evalpoly(tt_jc, -16_617.0, +2_004_192_130.00,    -429_775.20, -198_618.39,     +7.576,  +5.928_5) # μas
        Y = @evalpoly(tt_jc,  -6_951.0,        -25_817.00, -22_407_280.10,   +1_900.46, +1_112.526,  +0.135_8) # μas

        # TODO: confirm values for s
        s = @evalpoly(tt_jc,     +94.0,         +3_808.65,        -122.68,  -72_574.11,    +27.980, +15.620_0) # μas

    else
        error("Revision $revision not implemented")
    end

    fund_args = fundamental_arguments(tt_jc)

    X_harm_coeff = zeros(MVector{5, Float64})
    Y_harm_coeff = zeros(MVector{5, Float64})
    s_harm_coeff = zeros(MVector{5, Float64})

    @inbounds for frequency in frequencies

        arg = 0.0
        @inbounds @simd for k in frequency.coefficients_idx
            @muladd arg += frequency.coefficients[k] * fund_args[k]
        end
        si = sin(arg)
        co = cos(arg)

        @inbounds for component in frequency.components
            @muladd ampl = component.s * si + component.c * co

            if component.variable == 1
                X_harm_coeff[component.poweridx] += ampl
            elseif component.variable == 2
                Y_harm_coeff[component.poweridx] += ampl
            else
                s_harm_coeff[component.poweridx] += ampl
            end
        end
    end

    X += @evalpoly(tt_jc, X_harm_coeff[1], X_harm_coeff[2], X_harm_coeff[3], X_harm_coeff[4], X_harm_coeff[5])
    Y += @evalpoly(tt_jc, Y_harm_coeff[1], Y_harm_coeff[2], Y_harm_coeff[3], Y_harm_coeff[4], Y_harm_coeff[5])
    s += @evalpoly(tt_jc, s_harm_coeff[1], s_harm_coeff[2], s_harm_coeff[3], s_harm_coeff[4], s_harm_coeff[5])

    X = μas_to_rad(X)
    Y = μas_to_rad(Y)
    s = μas_to_rad(s)

    X, Y, s - X * Y / 2
end

const precession_nutation_06 = precession_nutation_IAU2006_IAU2000A_R06

end # module


text:
struct DeepTreeEnv
    depth::Int64
    adj_list::Dict{Int64, Vector{Int64}}
end

function DeepTreeEnv(depth::Int64, model::String, rng::MersenneTwister)
    @assert depth > 3
    num_nodes = compute_num_nodes(depth)
    list_of_nodes = collect(1:num_nodes)
    if model == "true"

    elseif model == "optimistic"

    else

    end
end

function compute_num_nodes(depth::Int64)
    sum([2^i for i in 0:depth])
end


text:
function PlotStairs(label_id, values::AbstractArray{T}, count::Integer, xscale::Real = 1.0, x0::Real = 0.0, offset::Integer = 0, stride::Integer = sizeof(T)) where {T<:ImPlotData}
    LibCImPlot.PlotStairs(label_id, values, count, xscale, x0, offset, stride)
end

function PlotStairs(label_id, values::AbstractArray{T}, count::Integer, xscale::Real = 1.0, x0::Real = 0.0, offset::Integer = 0, stride::Integer = sizeof(Float64)) where {T<:Real}
    LibCImPlot.PlotStairs(label_id, Float64.(values), count, xscale, x0, offset, stride)
end

function PlotStairs(label_id, x::AbstractArray{T}, y::AbstractArray{T}, count::Integer, offset::Integer = 0, stride::Integer = sizeof(T)) where {T<:ImPlotData}
    LibCImPlot.PlotStairs(label_id, x, y, count, offset, stride)
end

function PlotStairs(label_id, x::AbstractArray{T}, y::AbstractArray{T}, count::Integer, offset::Integer = 0, stride::Integer = sizeof(Float64)) where {T<:Real}
    LibCImPlot.PlotStairs(label_id, Float64.(x), Float64.(y), count, offset, stride)
end

function PlotStairs(y::AbstractArray{T};
    label_id::String="",
    count::Integer=length(y),
    xscale::Real = 1.0, x0::Real = 0.0,
    offset::Integer=0,
    stride::Integer=1
) where {T <: ImPlotData}

    LibCImPlot.PlotStairs(label_id, y, count, xscale, x0, offset, stride * sizeof(T))
end

function PlotStairs(x::AbstractArray{T}, y::AbstractArray{T};
    count::Integer = min(length(x), length(y)),
    offset::Integer = 0,
    stride::Integer = 1,
    label_id::String = ""
) where {T <: ImPlotData}

    LibCImPlot.PlotStairs(label_id, x, y, count, offset, stride * sizeof(T))
end


text:
using LinearAlgebra
using Distances

# unittest performance
# versus python hsic_varioustests_npy.py unittest:
#    20sec elapsed / 52sec user (i.e. multithreaded)
# using distmat1: 10.3sec
# using Distances.jl pairwise 13.2sec
# using distmatvec: 10.7sec

# performance discussion see https://discourse.julialang.org/t/sum-of-hadamard-products/3531/5

# for comparison, do the non-vectorized calculation
# points are in columns, dimensions(variables) across rows
# note transposed from previous convention!
function distmatslow(X,Y)
  m = size(X,2)
  @assert m==size(Y,2) #"size mismatch"
  D = Array{Float32}(undef,m,m)
  for ix = 1:m
    @inbounds for iy = 1:m
      x = X[:,ix]
      y = Y[:,iy]
      d = norm(x-y)
      D[ix,iy] = d*d
    end
  end
  return D
end

function distmat(X,Y)
  m = size(X,2)
  @assert m==size(Y,2) #"size mismatch"
  D = Array{Float32}(undef,m,m)
  for ix = 1:m
    @inbounds for iy = ix:m
      x = X[:,ix]
      y = Y[:,iy]
      d = x-y
      d2 = d'*d
      D[ix,iy] = d2
      D[iy,ix] = d2
    end
  end
  return D
end

#
function distmatvec(X,Y)
  """
  points are in columns, dimensions(variables) down rows
  """
  m = size(X,2)
  @assert m==size(Y,2) #"size mismatch"

  XY = sum(X .* Y,dims=1)
  #XY = XY.reshape(m,1)
  #R1 = n_.tile(XY,(1,Y.shape[0]))
  R1 = repeat(XY',1,m)  # outer=(1,m) not suported by autograd
  R2 = repeat(XY,m,1)
  #xy_ = X * Y'
  D = R1 + R2 - 2.f0 * X' * Y
  D
end

# python-ish version requires broadcasting different shapes - "outer sum"
function distmat1(X,Y)
  A = sum(X .* X,dims=1)'
  B = sum(Y .* Y,dims=1)
  C = X' * Y
  A .+ B .- 2.f0 * C
end

function eye(n)
  Matrix{Float32}(I,n,n)
end

function hsic(X,Y,sigma)
  #=
  # 1/m^2 Tr Kx H Ky H
  X,Y have data in COLUMNS, each row is a dimension
  # NB is transposed from python
  =#
  #println("X=\n", X'[1:2,:])
  #Yt = Y;; println("Y=\n", Yt[1:2:])
  m = size(X,2)
  println("hsic between ",m,"points")
  H = eye(m) - (1.f0 / m) * ones(Float32,m,m)
  #Dxx = distmatvec(X,X)
  #Dyy = distmatvec(Y,Y)
  Dxx = distmat1(X,X)
  Dyy = distmat1(Y,Y)
  #Dxx = pairwise(SqEuclidean(),X,X,dims=2)
  #Dyy = pairwise(SqEuclidean(),Y,Y,dims=2)
  sigma2 = 2.f0 * sigma*sigma
  Kx = exp.( -Dxx / sigma2 )
  Ky = exp.( -Dyy / sigma2 )
  Kxc = Kx * H
  Kyc = Ky * H
  thehsic = (1.f0 / (m*m)) * sum(Kxc' .* Kyc)
  return thehsic	# type float32
end

# Pkg.add("NPZ")
using NPZ
# aug19: this matches the output of the unittest in hsic_varioustests_npy
function unittest()
  X = Float32[0.1 0.2 0.3;
                5 4 3]
  Y = Float32[1 2 3;
               2 2 2]
  X = transpose(X)
  Y = transpose(Y)
  println(distmatslow(X,X))
  println(distmat(X,X))
  println(distmatslow(Y,Y))
  println(distmat(Y,Y))
  println("hsic(X,Y,0.5)=",hsic(X,Y,0.5f0))

  # larger test
  data = npzread("/tmp/_data.npz")
  # todo convert arrays to float32
  X = data["arr_0"]
  Y = data["arr_1"]
  X = convert(Array{Float32,2},X')
  Y = convert(Array{Float32,2},Y')
  println("X=\n", X'[1:2,:])
  println("Y=\n", Y'[1:2,:])

  println("\nindependent hsic(X,Y,1)=",hsic(X,Y,1.f0))

  println("\nidentical hsic(X,X,1)=",hsic(X,copy(X),1.f0))

  Y2 = X .* X
  println("Y2=",Y2'[1:2,:])
  println("\nnonlinear hsic(X,Y*Y,1)=",hsic(X,Y2,1.f0))
end


text:
using KernelAbstractions.Extras.LoopInfo: @unroll

#####
##### Periodic boundary conditions
#####

@kernel function fill_periodic_west_and_east_halo!(c, H::Int, N)
  j, k = @index(Global, NTuple)
  @unroll for i = 1:H
      @inbounds begin
          c[i, j, k] = c[N+i, j, k]     # west
          c[N+H+i, j, k] = c[H+i, j, k] # east
      end
  end
end

@kernel function fill_periodic_south_and_north_halo!(c, H::Int, N)
  i, k = @index(Global, NTuple)
  @unroll for j = 1:H
      @inbounds begin
          c[i, j, k] = c[i, N+j, k]     # south
          c[i, N+H+j, k] = c[i, H+j, k] # north
      end
  end
end

@kernel function fill_periodic_bottom_and_top_halo!(c, H::Int, N)
  i, j = @index(Global, NTuple)
  @unroll for k = 1:H
      @inbounds begin
          c[i, j, k] = c[i, j, N+k]        # top
          c[i, j, N+H+k] = c[i, j, H+k] # bottom
      end
  end
end

function fill_west_and_east_halo!(c, ::PBC, ::PBC, arch, dep, grid, args...; kw...)
  c_parent = parent(c)
  yz_size = size(c_parent)[[2, 3]]
  event = launch!(arch, grid, yz_size, fill_periodic_west_and_east_halo!, c_parent, grid.Hx, grid.Nx; dependencies=dep, kw...)
  return event
end

function fill_south_and_north_halo!(c, ::PBC, ::PBC, arch, dep, grid, args...; kw...)
  c_parent = parent(c)
  xz_size = size(c_parent)[[1, 3]]
  event = launch!(arch, grid, xz_size, fill_periodic_south_and_north_halo!, c_parent, grid.Hy, grid.Ny; dependencies=dep, kw...)
  return event
end

function fill_bottom_and_top_halo!(c, ::PBC, ::PBC, arch, dep, grid, args...; kw...)
  c_parent = parent(c)
  xy_size = size(c_parent)[[1, 2]]
  event = launch!(arch, grid, xy_size, fill_periodic_bottom_and_top_halo!, c_parent, grid.Hz, grid.Nz; dependencies=dep, kw...)
  return event
end

####
#### Implement single periodic directions (for Distributed memory architectures)
####

@kernel function fill_periodic_west_halo!(c, H::Int, N)
  j, k = @index(Global, NTuple)
  @unroll for i = 1:H
      @inbounds begin
          c[i, j, k] = c[N+i, j, k]     # west
      end
  end
end

@kernel function fill_periodic_east_halo!(c, H::Int, N)
  j, k = @index(Global, NTuple)
  @unroll for i = 1:H
      @inbounds begin
          c[N+H+i, j, k] = c[H+i, j, k] # east
      end
  end
end

@kernel function fill_periodic_south_halo!(c, H::Int, N)
  i, k = @index(Global, NTuple)
  @unroll for j = 1:H
      @inbounds begin
          c[i, j, k] = c[i, N+j, k]     # south
      end
  end
end

@kernel function fill_periodic_north_halo!(c, H::Int, N)
  i, k = @index(Global, NTuple)
  @unroll for j = 1:H
      @inbounds begin
          c[i, N+H+j, k] = c[i, H+j, k] # north
      end
  end
end

@kernel function fill_periodic_bottom_halo!(c, H::Int, N)
  i, j = @index(Global, NTuple)
  @unroll for k = 1:H
      @inbounds begin
          c[i, j, k] = c[i, j, N+k]        # bottom
      end
  end
end

@kernel function fill_periodic_top_halo!(c, H::Int, N)
  i, j = @index(Global, NTuple)
  @unroll for k = 1:H
      @inbounds begin
          c[i, j, N+H+k] = c[i, j, H+k] # top
      end
  end
end

function fill_west_halo!(c, ::PBC, arch, dep, grid, args...; kw...)
  c_parent = parent(c)
  yz_size = size(c_parent)[[2, 3]]
  event = launch!(arch, grid, yz_size, fill_periodic_west_halo!, c_parent, grid.Hx, grid.Nx; dependencies=dep, kw...)
  return event
end

function fill_east_halo!(c, ::PBC, arch, dep, grid, args...; kw...)
  c_parent = parent(c)
  yz_size = size(c_parent)[[2, 3]]
  event = launch!(arch, grid, yz_size, fill_periodic_east_halo!, c_parent, grid.Hx, grid.Nx; dependencies=dep, kw...)
  return event
end

function fill_south_halo!(c, ::PBC, arch, dep, grid, args...; kw...)
  c_parent = parent(c)
  xz_size = size(c_parent)[[1, 3]]
  event = launch!(arch, grid, xz_size, fill_periodic_south_halo!, c_parent, grid.Hy, grid.Ny; dependencies=dep, kw...)
  return event
end

function fill_north_halo!(c, ::PBC, arch, dep, grid, args...; kw...)
  c_parent = parent(c)
  xz_size = size(c_parent)[[1, 3]]
  event = launch!(arch, grid, xz_size, fill_periodic_north_halo!, c_parent, grid.Hy, grid.Ny; dependencies=dep, kw...)
  return event
end

function fill_bottom_halo!(c, ::PBC, arch, dep, grid, args...; kw...)
  c_parent = parent(c)
  xy_size = size(c_parent)[[1, 2]]
  event = launch!(arch, grid, xy_size, fill_periodic_bottom_halo!, c_parent, grid.Hz, grid.Nz; dependencies=dep, kw...)
  return event
end

function fill_top_halo!(c, ::PBC, arch, dep, grid, args...; kw...)
  c_parent = parent(c)
  xy_size = size(c_parent)[[1, 2]]
  event = launch!(arch, grid, xy_size, fill_periodic_top_halo!, c_parent, grid.Hz, grid.Nz; dependencies=dep, kw...)
  return event
end


text:
using Printf
import Flux: σ
using ModelingToolkit
using GalacticOptim
using Optim
using DiffEqFlux
using NeuralPDE
using Quadrature, Cubature, Cuba
using Plots

@parameters t,x
@variables c(..)
@derivatives Dt'~t
@derivatives Dxx''~x
@derivatives Dx'~x

# Parameters

v = 1
R = 0
D = 0.1 # diffusion
t_max = 2.0
x_min = -1.0
x_max = 1.0

# Equations, initial and boundary conditions
eqs = [ Dt(c(t, x)) ~ D * Dxx(c(t,x)) - Dx(c(t,x)) ]

bcs = [
        c(0, x) ~ cos(π*x) + 1.0,
        c(t, x_min) ~ c(t, x_max)

]

# Space and time domains
domains = [t ∈ IntervalDomain(0.0,t_max),
        x ∈ IntervalDomain(x_min,x_max)
]

# Discretization
nx = 32
dx = (x_max-x_min) / (nx - 1)
dt = 0.01

# Neural network
dim = length(domains)
output = length(eqs)
hidden = 8

chain = FastChain( FastDense(dim, hidden, σ),
                    FastDense(hidden, hidden, σ),
                    FastDense(hidden, 1))

strategy = GridTraining(dx=[dt,dx])

discretization = PhysicsInformedNN(chain, strategy=strategy)

pde_system = PDESystem(eqs, bcs, domains, [t,x], [c])
prob = discretize(pde_system,discretization)

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

res = GalacticOptim.solve(prob,Optim.BFGS();cb=cb,maxiters=1200)

# Plots

phi = discretization.phi

initθ = discretization.initθ

acum =  [0;accumulate(+, length.(initθ))]
sep = [acum[i]+1 : acum[i+1] for i in 1:length(acum)-1]
minimizers = [res.minimizer[s] for s in sep]
ts,xs = [domain.domain.lower:dx:domain.domain.upper for domain in domains]

anim = @animate for (i, t) in enumerate(0:dt:t_max)
    @info "Animating frame $i..."
    c_predict = reshape([phi([t, x], res.minimizer)[1] for x in xs], length(xs))
    title = @sprintf("Advection-diffusion t = %.3f", t)
    plot(xs, c_predict, label="", title=title , ylims=(0., 2))
end

gif(anim, "advection_diffusion_pinn.gif", fps=15)

c_predict = reshape([ phi([0, x], res.minimizer)[1] for x in xs], length(xs))
plot(xs,c_predict)

# Plot correct solution
using JLD2
file = jldopen("advection_diffusion/simulation/cosine_advection_diffusion.jld2")
iterations = parse.(Int, keys(file["timeseries/t"]))

anim = @animate for (i, iter) in enumerate(iterations)
    @info "Animating frame $i..."
    Hx = file["grid/Hx"]
    x = file["grid/xC"][1+Hx:end-Hx]
    t = file["timeseries/t/$iter"]
    c = file["timeseries/c/$iter"][:]

    title = @sprintf("Advection-diffusion t = %.3f", t)
    p = plot(x, c .+ 1, linewidth=2, title=title, label="Oceananigans",
             xlabel="x", ylabel="Tracer", xlims=(-1, 1), ylims=(0, 2))

    c_predict = reshape([phi([t, x], res.minimizer)[1] for x in xs], length(xs))

    plot!(p, x, c_predict, linewidth=2, label="Neural PDE")
end

gif(anim, "advection_diffusion_comparison.gif", fps=15)

max_stars_repo_path: Old Stuff/1_convection_diffusion.jl

max_stars_repo_name: drozzy/NeuralOceananigans

max_issues_repo_path: Old Stuff/1_convection_diffusion.jl

max_issues_repo_name: drozzy/NeuralOceananigans

max_forks_repo_path: Old Stuff/1_convection_diffusion.jl

max_forks_repo_name: drozzy/NeuralOceananigans


text:
module AutoARIMA

using LinearAlgebra
using Optim
using Polynomials
using RecipesBase
using StaticArrays
using Statistics
using HypothesisTests

export seriesA,seriesB,seriesB2,seriesC,seriesD,seriesE,seriesF,seriesG,dowj,wine,lake
export autocovariance,autocovariance_matrix,autocorrelation,autocorrelation_matrix,partial_autocorrelation
export isinversible, isstationary
export innovations, levinson_durbin, least_squares, yule_walker, hannan_rissanen
export boxcox, inv_boxcox, guerrero, difference, integrate
export simulate, forecast, fit, toarma, k, residuals
export aic, aicc, bic, mse, rmse, mae, mape
export correlogram,partial_correlogram
export AbstractModel, AbstractParams
export ARParams, MAParams, ARMAParams, ARIMAParams, MSARIMAParams
export ARModel, MAModel, ARMAModel
export MA∞, AR∞
export boxjenkins

include("datasets.jl")
include("stats.jl")
include("transforms.jl")
include("abstract.jl")
include("criteria.jl")
include("simulate.jl")
include("recipes.jl")
include("ls.jl")
include("arma.jl")
include("ar.jl")
include("ma.jl")
include("arima.jl")
include("sarima.jl")
include("auto.jl")

end


text:
module GLM_

import MLJBase

export OLSRegressor #, OLS, LinearRegression

import ..GLM # strange syntax for lazy-loading

const OLSFitResult = GLM.LinearModel

OLSFitResult(coefs::Vector, b=nothing) = OLSFitResult(coefs, b)

####
#### OLSRegressor
####

mutable struct OLSRegressor <: MLJBase.Deterministic{OLSFitResult}
    fit_intercept::Bool
#    allowrankdeficient::Bool
end

function OLSRegressor(;fit_intercept=true)
    return OLSRegressor(fit_intercept)
end

# synonyms
# const OLS = OLSRegressor
# const LinearRegression = OLSRegressor

####
#### fit/predict OLSRegressor
####

function MLJBase.fit(model::OLSRegressor, verbosity::Int, X, y::Vector)

    Xmatrix = MLJBase.matrix(X)
    features = MLJBase.schema(X).names

    if model.fit_intercept
        fitresult = GLM.lm(hcat(Xmatrix, ones(eltype(Xmatrix), size(Xmatrix, 1), 1)), y)
    else
        fitresult = GLM.lm(Xmatrix, y)
    end

    coefs = GLM.coef(fitresult)

    ## TODO: add feature importance curve to report using `features`
    report = Dict(:coef => coefs[1:end-Int(model.fit_intercept)]
                , :intercept => ifelse(model.fit_intercept, coefs[end], nothing)
                , :deviance => GLM.deviance(fitresult)
                , :dof_residual => GLM.dof_residual(fitresult)
                , :stderror => GLM.stderror(fitresult)
                , :vcov => GLM.vcov(fitresult))
    cache = nothing

    return fitresult, cache, report
end

function MLJBase.predict(model::OLSRegressor, fitresult::OLSFitResult, Xnew)
    Xmatrix = MLJBase.matrix(Xnew)
    model.fit_intercept && (Xmatrix = hcat(Xmatrix, ones(eltype(Xmatrix), size(Xmatrix, 1), 1)))
    return GLM.predict(fitresult, Xmatrix)
end

# metadata:
MLJBase.load_path(::Type{<:OLSRegressor}) = "MLJModels.GLM_.OLSRegressor"
MLJBase.package_name(::Type{<:OLSRegressor}) = "GLM"
MLJBase.package_uuid(::Type{<:OLSRegressor}) = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
MLJBase.package_url(::Type{<:OLSRegressor}) = "https://github.com/JuliaStats/GLM.jl"
MLJBase.is_pure_julia(::Type{<:OLSRegressor}) = :yes
MLJBase.input_kinds(::Type{<:OLSRegressor}) = [:continuous, ]
MLJBase.output_kind(::Type{<:OLSRegressor}) = :continuous
MLJBase.output_quantity(::Type{<:OLSRegressor}) = :univariate

end # module

max_stars_repo_name: ablaom/MLJModels

max_issues_repo_name: ablaom/MLJModels

max_forks_repo_name: ablaom/MLJModels


text:
function run_simulation(set_sizes; fun = conjunctive_set, kwargs...)
    results = DataFrame[]
    for n in set_sizes
        experiment = Experiment(set_size = n, populate_visicon = fun,
            n_trials = 10^4)
        run_condition!(experiment; kwargs...)
        df = DataFrame(experiment.data)
        hit_rate = mean(df.target_present .== df.response)
        g = groupby(df, [:target_present,:response])
        temp = combine(g, :rt => mean)
        temp[!,:distractors] .= n
        temp[!,:hit_rate] .= hit_rate
        push!(results, temp)
    end
    return vcat(results...)
end