Add benchmarks for pseudo remainder method

[1e9abhi1e10]

This PR adds the benchmarks for the pseudo-remainder method. Related Issue:

88

See sympy/sympy#25278

[1e9abhi1e10]

@oscarbenjamin if anything you want to modify in the code, then please tell me. Also, I little bit confused about the timings that how it should be calculated.

[oscarbenjamin]

I think this needs to be changed to 1.12: https://github.com/sympy/sympy_benchmarks/blob/3680cc194e5af27663c7e46215bf89789766291e/.github/workflows/run-benchmarks.yml#L34

[1e9abhi1e10]

I think this needs to be changed to 1.12:

https://github.com/sympy/sympy_benchmarks/blob/3680cc194e5af27663c7e46215bf89789766291e/.github/workflows/run-benchmarks.yml#L34

make a separate PR for this?

[oscarbenjamin]

make a separate PR for this?

It could be. The tests won't run on this one until that change is made.

[1e9abhi1e10]

there are almost 40 values of F and G from these benchmarks, should I include all ?

[oscarbenjamin]

there are almost 40 values of F and G from these benchmarks, should I include all ?

Maybe just pick some representative ones. It is good to time with smaller examples and also larger examples from each group. Also benchmarks should not take longer than around 1 second so don't use really large examples. We can still include the code to generate the examples of any size but we don't need to have the benchmark suite run all sizes every time it runs.

[oscarbenjamin]

Code should be added that can generate pairs of polynomials of different sizes as well as their gcd and what the prem should be. We can reuse this for other benchmarks like gcd later.

The Time class should have a params attribute that says which sizes to use for timings when running the benchmarks suite and the setup method should take a parameter that will be the values from params.

The purpose of storing the results in values is so that the tear down method can check that the result is correct (there is no point reporting timings if a wrong result is computed).

There should be time_ methods for different prem functions.

Many of these things are demonstrated by the example I linked before: https://github.com/sympy/sympy_benchmarks/blob/3680cc194e5af27663c7e46215bf89789766291e/benchmarks/integrate.py#L46

[oscarbenjamin]

If you look in the CI output under "Run benchmarks" you can find this output which shows the other polys benchmarks being run:

[ 42.24%] ··· ======== =============
               param1               
              -------- -------------
                 1       43.9±0.1μs 
                 10     74.9±0.09μs 
                100       699±2μs   
                500     12.5±0.05ms 
              ======== =============

[ 42.41%] ··· polys.TimePolyManyGens.time_is_linear                           ok
[ 42.41%] ··· ======== =============
               param1               
              -------- -------------
                 1      1.69±0.01μs 
                 10      19.9±0.1μs 
                100      569±0.3μs  
                500     12.0±0.01ms 
              ======== =============

I don't see any output for the prem benchmarks that you have added here so it looks like they are not being run by asv.

Do you see any output for the prem benchmarks when running asv locally?

[oscarbenjamin]

I'm not sure why the benchmark doesn't seem to run in CI.

[1e9abhi1e10]

Do you see any output for the prem benchmarks when running asv locally?

I am currently experiencing certain challenges while running asv locally, but I am determined to find solutions to overcome them. Do I need to use Python version 3.10 specifically for executing the benchmarks?

[oscarbenjamin]

I think that it should be possible with any version of Python.

The problem might be that the time_ methods do not have a parameter.

[oscarbenjamin]

Let's see if the added benchmarks run now.

[1e9abhi1e10]

Let's see if the added benchmarks run now.

Now the benchmarks for the prem method are running, but they are failing due to some reason. I will check it

[oscarbenjamin]

You will need to be able to run the benchmarks locally to get it working. Debugging via CI is impossible.

[1e9abhi1e10]

@oscarbenjamin Please have a look!

[oscarbenjamin]

This increases the total time taken to run the benchmarks from 15 minutes to 25 minutes which is too much. Some of the benchmarks are too slow and should be removed. Generally this benchmark suite cannot include anything that takes longer than a second and really we should limit to around 100ms.

For example time_prem_case1 with n=8 takes 4 seconds and with n=10 takes 17 seconds. Given that asv will run these multiple times in a loop this is too slow for what we can include in the benchmark suite. That means that time_prem_case1 should not be used with n=8 or n=10.

There are different ways to control this but one way is to make separate classes that have different values for params and try different sizes e.g.: https://github.com/sympy/sympy_benchmarks/blob/3680cc194e5af27663c7e46215bf89789766291e/benchmarks/solve.py#L498-L510

[1e9abhi1e10]

Should I add examples of f and g for specific values of n? If yes then will I add around 16 different values for different cases and specific values of n.

Like these: https://github.com/sympy/sympy_benchmarks/blob/3680cc194e5af27663c7e46215bf89789766291e/benchmarks/solve.py#L423-L433

[oscarbenjamin]

Should I add examples of f and g for specific values of n?

No, it is better to keep the code that generates them.

Like these:

Those are only done like that because the code that generates them needs PyDy which is a separate library. The examples are hard-coded to avoid having the benchmarks depend on other libraries.

[oscarbenjamin]

This is still too slow. It takes the benchmark suite from 15 minutes to 20 minutes. We don't need 25% of the test run to just be spent on prem and we will want to benchmark more methods later.

[1e9abhi1e10]

@oscarbenjamin In the next step, should I add benchmarks for the pdiv method? Also, should I include those benchmarks in this pull request (PR), or should I create a new one?

[1e9abhi1e10]

I think the time_ functions for prem method is falling due to teardown function.

[oscarbenjamin]

The asv output indicates that this is not working:

[ 41.78%] ··· ...REMLinearDenseQuadraticGCD.time_prem_methods             failed
[ 41.78%] ··· ================== ======== ======== ========
              --                           param2          
              ------------------ --------------------------
                    param1          1        3        5    
              ================== ======== ======== ========
                     prem         failed   failed   failed 
                  Poly_prem       failed   failed   failed 
               PolyElement_prem   failed   failed   failed 
              ================== ======== ======== ========

[ 41.95%] ··· ...mePREMQuadraticNonMonicGCD.time_prem_methods             failed
[ 41.95%] ··· ================== ======== ======== ========
              --                           param2          
              ------------------ --------------------------
                    param1          1        3        5    
              ================== ======== ======== ========
                     prem         failed   failed   failed 
                  Poly_prem       failed   failed   failed 
               PolyElement_prem   failed   failed   failed 
              ================== ======== ======== ========

The table looks nice but we need numbers rather than "failed" :)

[1e9abhi1e10]

Sorry, I am trying to squash the commits but there is a problem with the extension. Soon I will fix it!

[1e9abhi1e10]

Why the asv output is structured in this way?

[ 42.11%] ··· ...imePREMSparseGCDHighDegree.time_prem_methods                 ok
[ 42.11%] ··· ================== ======== =============
                    param1        param2               
              ------------------ -------- -------------
                     prem           1        965±20μs  
                     prem           3        1.86±0ms  
                     prem           5      3.38±0.05ms 
                     prem           8      6.31±0.03ms 
                  Poly_prem         1        68.4±2μs  
                  Poly_prem         3      2.63±0.05ms 
                  Poly_prem         5       11.6±0.4ms 
                  Poly_prem         8       45.9±0.6ms 
               PolyElement_prem     1      3.99±0.01ms 
               PolyElement_prem     3      12.1±0.01ms 
               PolyElement_prem     5      26.3±0.05ms 
               PolyElement_prem     8       46.0±0.1ms 
              ================== ======== =============

It should be like this

[100.00%] ··· polys.TimePREMSparseGCDHighDegree.time_prem_methods                                                                                                                                          ok
[100.00%] ··· ================== ============ ============= ============ =============
              --                                         param2
              ------------------ -----------------------------------------------------
                    param1            1             3            5             8
              ================== ============ ============= ============ =============
                     prem          545±20μs    1.11±0.03ms   2.04±0.2ms   3.76±0.04ms
                  Poly_prem        39.8±3μs    1.66±0.03ms   7.90±0.4ms    29.6±0.8ms
               PolyElement_prem   2.71±0.1ms    7.85±0.3ms    18.2±1ms     31.3±0.9ms
              ================== ============ ============= ============ =============

Should I use tuple of a tuple for params?

[oscarbenjamin]

Why the asv output is structured in this way?

Possibly there are too many cases for a horizontal layout. It make a difference to make the names shorter e.g. PolyElement_prem could just be sparse or something. We don't need prem to be part of the name because it is already shown as part of the prem benchmark.

It might be that asv tries to detect the "terminal width" or something. Maybe there is an option to asv to control that.

[1e9abhi1e10]

Should this problem be resolved?

[ 42.11%] ··· ...imePREMSparseGCDHighDegree.time_prem_methods                 ok
[ 42.11%] ··· ======== ======== =============
               param1   param2               
              -------- -------- -------------
                prem      1        931±5μs   
                prem      3      1.85±0.01ms 
                prem      5      3.20±0.02ms 
                prem      8      6.14±0.05ms 
                poly      1       64.3±0.5μs 
                poly      3      2.47±0.02ms 
                poly      5      11.3±0.03ms 
                poly      8       43.6±0.2ms 
               sparse     1      4.00±0.02ms 
               sparse     3      12.0±0.04ms 
               sparse     5       26.2±0.1ms 
               sparse     8       45.9±0.2ms 
              ======== ======== =============

[oscarbenjamin]

Should this problem be resolved?

I don't know if it is an easy thing to change but if not then that is fine.

It looks like param_names can be used to control the names that are displayed: https://asv.readthedocs.io/en/stable/benchmarks.html#benchmark-attributes

I don't know if there is anything in the docs about changing the layout.

[oscarbenjamin]

This is looking a lot better but the structure is still not quite right for what we want to do when adding more cases later. I reorganised it a bit like this:

from sympy import *

class _GCDExample:
    """A benchmark example with two polynomials and their gcd."""

    def __init__(self, n):
        f, g, d, syms = self.make_poly(n)
        self.f = f
        self.g = g
        self.d = d
        self.x = syms[0]
        self.y = syms[1:]
        self.syms = syms
        self.ring = QQ[syms]

    def make_poly(self, n):
        raise NotImplementedError("Subclass should implement this.")

    def to_poly(self, expr):
        return Poly(expr, self.syms)

    def to_ring(self, expr):
        return self.ring(expr)

    def as_expr(self):
        return (self.f, self.g, self.d, self.syms)

    def as_poly(self):
        return (self.to_poly(self.f), self.to_poly(self.g), self.to_poly(self.d))

    def as_ring(self):
        return (self.to_ring(self.f), self.to_ring(self.g), self.to_ring(self.d), self.ring)

class _SparseGCDHighDegree(_GCDExample):
    """A pair of polynomials in n symbols with a high degree sparse GCD.

    >>> for n in [1, 2, 3, 4, 5]:
    ...     f, g, d, syms = _SparseGCDHighDegree(n).as_expr()
    ...     print(d)
    x**2 + y1**2 + 1
    x**3 + y1**3 + y2**3 + 1
    x**4 + y1**4 + y2**4 + y3**4 + 1
    x**5 + y1**5 + y2**5 + y3**5 + y4**5 + 1
    x**6 + y1**6 + y2**6 + y3**6 + y4**6 + y5**6 + 1
    """

    def make_poly(self, n):
        x, *y = syms = symbols("x, y1:{}".format(n+1))
        d = 1 + x ** (n + 1) + sum([y[i] ** (n + 1) for i in range(n)])
        f = d * (-2 + x ** (n + 1) + sum([y[i] ** (n + 1) for i in range(n)]))
        g = d * (2 + x ** (n + 1) + sum([y[i] ** (n + 1) for i in range(n)]))
        return f, g, d, syms

class _QuadraticNonMonicGCD(_GCDExample):
    """A pair of quadratic polynomials with a non-monic GCD.

    >>> for n in [1, 2, 3, 4, 5]:
    ...     f, g, d, syms = _QuadraticNonMonicGCD(n).as_expr()
    ...     print(d)
    x**2*y1**2 + 1
    x**2*y1**2 + y2**2 + 1
    x**2*y1**2 + y2**2 + y3**2 + 1
    x**2*y1**2 + y2**2 + y3**2 + y4**2 + 1
    x**2*y1**2 + y2**2 + y3**2 + y4**2 + y5**2 + 1
    """
    def make_poly(self, n):
        x, *y = syms = symbols("x, y1:{}".format(n+1))
        d = 1 + x ** 2 * y[0] ** 2 + sum([y[i] ** 2 for i in range(1, n)])
        f = d * (-1 + x ** 2 - y[0] ** 2 + sum([y[i] ** 2 for i in range(1, n)]))
        g = d * (2 + x * y[0] + sum(y[1:n])) ** 2
        return f, g, d, syms

class _TimeOP:
    """
    Benchmarks comparing Poly implementations of a given operation.
    """
    param_names = [
        'size',
        'impl',
    ]

    def setup(self, n, impl):

        examples = self.GCDExampleCLS(n)

        expected = self.expected(*examples.as_expr())

        if impl == 'expr':
            func = self.get_func_expr(*examples.as_expr())
            expected = examples.to_expr(expected)
        elif impl == 'dense':
            func = self.get_func_poly(*examples.as_poly())
            expected = examples.to_poly(expected)
        elif impl == 'sparse':
            func = self.get_func_sparse(*examples.as_ring())
            expected = examples.to_ring(expected)

        self.func = func
        self.expected = expected

    def time_op(self, n, impl):
        self.returned = self.func()

    def teardown(self, n, impl):
        assert self.expected == self.returned

class _TimePREM(_TimeOP):

    def expected(self, f, g, d, syms):
        prem_f_g_x = rem(f * LC(g, x) ** (degree(f, x) - degree(g, x) + 1), g, x)
        return prem(f, g, x)

    def get_func_expr(self, f, g, d, syms):
        x = syms[0]
        return lambda: prem(f, g, x)

    def get_func_poly(self, f, g, d):
        return lambda: f.prem(g)

    def get_func_sparse(self, f, g, d, ring):
        return lambda: f.prem(g)

class TimePREM_SparseGCDHighDegree(_TimePREM):
    GCDExampleCLS = _SparseGCDHighDegree
    params = [(1, 3, 5), ('expr', 'dense', 'sparse')]

The key point about this is that if we want to test something else like gcd or pquo etc then it is very easy to add that. We just add a _TimeGCD class and everything else carries over for free without duplicating lots of repetitive code:

class _TimeGCD(_TimeOP):

    def expected(self, f, g, d, syms):
        return d

    def get_func_expr(self, f, g, d, syms):
        return lambda: gcd(f, g)

    def get_func_poly(self, f, g, d):
        return lambda: f.gcd(g)

    def get_func_sparse(self, f, g, d, ring):
        return lambda: f.gcd(g)

Note how this _TimeGCD class contains the minimum possible code that might not be reused for anything other than the specific case of timing GCD across the three implementations. Everything else is factored out.

[oscarbenjamin]

So in the end the output looks like:

[ 41.78%] ··· polys.TimePREM_LinearDenseQuadraticGCD.time_op                  ok
[ 41.78%] ··· ====== ============= ============= =============
              --                        impl                  
              ------ -----------------------------------------
               size       expr         dense         sparse   
              ====== ============= ============= =============
                1       1.54±0ms    74.6±0.08μs    222±0.5μs  
                3     6.52±0.01ms     1.53±0ms      2.43±0ms  
                5     20.2±0.04ms   8.79±0.04ms   12.4±0.01ms 
              ====== ============= ============= =============

[ 41.95%] ··· polys.TimePREM_QuadraticNonMonicGCD.time_op                     ok
[ 41.95%] ··· ====== ============= ============= =============
              --                        impl                  
              ------ -----------------------------------------
               size       expr         dense         sparse   
              ====== ============= ============= =============
                1       773±2μs      105±0.4μs     411±0.9μs  
                3       3.00±0ms    5.24±0.01ms   5.99±0.02ms 
                5     8.33±0.01ms   10.6±0.03ms   12.9±0.05ms 
              ====== ============= ============= =============

[ 42.11%] ··· polys.TimePREM_SparseGCDHighDegree.time_op                      ok
[ 42.11%] ··· ====== ============= ============= =============
              --                        impl                  
              ------ -----------------------------------------
               size       expr         dense         sparse   
              ====== ============= ============= =============
                1       723±2μs      53.9±0.1μs    151±0.8μs  
                3       1.49±0ms    2.16±0.01ms   2.57±0.01ms 
                5     2.65±0.01ms   9.62±0.02ms   10.9±0.06ms 
                8     5.17±0.01ms   36.6±0.05ms   40.0±0.05ms 
              ====== ============= ============= =============

[ 42.28%] ··· polys.TimePREM_SparseNonMonicQuadratic.time_op                  ok
[ 42.28%] ··· ====== ========= ============= =============
              --                      impl                
              ------ -------------------------------------
               size     expr       dense         sparse   
              ====== ========= ============= =============
                1     467±1μs    62.4±0.1μs    191±0.2μs  
                3     583±1μs     1.79±0ms      2.08±0ms  
                5     693±2μs   4.88±0.01ms   5.30±0.01ms 
                8     858±3μs   9.64±0.03ms   10.2±0.01ms 
              ====== ========= ============= =============

[oscarbenjamin]

I think this looks good.

Thanks, and let's get it merged.

[1e9abhi1e10]

Thanks a lot!