Implement @precompute decorator

[saidctb]

Previously, the @python decorator was utilized to signal to the compiler that a function should be executed at compile time, and its result should be directly incorporated into the code. However, this mechanism was removed in the issue #668 because it was not handled properly.
I want to reintroduce this mechanism and make it more general than before as I think it can be useful in many instances, where we need to precompute some expressions and hardcode them directly in the code to make the code faster.
Additionally, it can be used in unit testing and generate an executable that does everything instead of separating the implementation from the unit testing.
Key Features:

• The function that contains the @precompute decorator should be able to call other functions in the same file or access the imported files. previously the function with the @python didn't have access to other functions in the same file.
• The results from calling a function with the @precompute decorator should be declared as Final, to ensure immutability.
• The translated code will not contain those functions unless we can handle callbacks. In that case, Fortran/C can call the Python implementation of those functions.
• The functions with a precompute decorator can only be called using literals and constant arrays, we can also handle arguments that are declared Final because we can extract their value.
  Here I constructed many examples to illustrate its usefulness.
• Compute matrix-matrix product and use precompute decorator to generate test cases:

import numpy as np

# Matrix Matrix product
def matmat(A:'float[:,:]', B:'float[:,:]'):
    """
    Compute the matrix matrix product.

    """
    ...

    return C
#======================================================
#UNIT TESTING
# Use the precompute decorator to create the test cases in the compilation time
@precompute
def generate_matrices(n:'int'):
    """ Generate matrices and compute their product using the numpy library"""
    A = np.random.rand(n,n)
    B = np.random.rand(n,n)
    C = A@B
    return A,B,C

if __name__ == '__main__':

    A1,B1,C1 = generate_matrices(5)
    A2,B2,C2 = generate_matrices(8)
    A3,B3,C3 = generate_matrices(12)

    err1 = abs(C1-matmat(A1,B1)).sum()
    err2 = abs(C2-matmat(A2,B2)).sum()
    err3 = abs(C3-matmat(A3,B3)).sum()
    for r in (err1, err2, err3):
        if abs(r)<tol:
            print("PASSED!")
        else:
            print('FAILED!')

• Compute numerical integrals and check the results using the Sympy library:

import numpy as np

# Define the functions to be integrated
def f1(x:'float'):
    return x**2 + np.sin(x)

def f2(x:'float'):
    return np.cos(x)
    
def f3(x:'float'):
    return np.tan(x)**2
    
def f4(x:'float'):
    return x**2 + x**3

# numerical integration
def numerical_integral(func:'(float)(float)', a:'float', b:'float'):
    """
    Compute the definite integral of a function numerically using the trapezoidal rule.
    """
    ...
    return result

#======================================================
#UNIT TESTING
def numerical_integrals(a:'float',b:'float'):
    integral1 = numerical_integral(f1, a, b)
    integral2 = numerical_integral(f2, a, b)
    integral3 = numerical_integral(f3, a, b)
    integral4 = numerical_integral(f4, a, b)
    return integral1,integral2, integral3, integral4

# Use the precompute decorator to compute exact integrals
@precompute
def exact_integrals(a:'float',b:'float'):
    from sympy.abc import x
    integral1 = sp.integrate(f1(x), (x, a, b))
    integral2 = sp.integrate(f2(x), (x, a, b))
    integral3 = sp.integrate(f3(x), (x, a, b))
    integral4 = sp.integrate(f4(x), (x, a, b))
    return integral1,integral2, integral3, integral4

if __name__ == '__main__':

    num = numerical_integrals(0.,1.)
    ext = exact_integrals(0.,1.)
    for r1,r2 in zip(num, ext):
        if abs(r1-r2)<tol:
            print("PASSED!")
        else:
            print('FAILED!')

• Compute the minimum of a function using gradient descent algorithm and compare with an exact solution:

import numpy as np

# Define the function to be optimized
def f(x:'float'):
    return (x - 2) ** 2 + np.sin(x)

# Gradient descent optimization
def gradient_descent(func, initial_guess, learning_rate=0.01, max_iter=1000, tol=1e-6):
    x = initial_guess
    for _ in range(max_iter):
        grad = (func(x + tol) - func(x)) / tol  # Approximate gradient
        x_new = x - learning_rate * grad
        if abs(x_new - x) < tol:
            break
        x = x_new
    return x

#======================================================
# UNIT TESTING
@precompute
def symbolic_optimization():
    import sympy as sp
    x = sp.Symbol('x')
    gradient = sp.diff(f(x), x)
    critical_points = sp.solve(gradient, x)
    return float(critical_points)

@precompute
def exact_critical_points():
    return symbolic_optimization()

if __name__ == '__main__':
    initial_guess = 0  # Initial guess for numerical optimization

    # Numerical optimization using gradient descent
    numerical_result = gradient_descent(f, initial_guess)

    # Exact optimization
    exact_result = exact_critical_points()

    # Unit test
    if np.isclose(numerical_result, exact_result):
        print("PASSED!")
    else:
        print("FAILED!")

[EmilyBourne]

What is the difference between this @precompute decorator and the @sympy decorator?

[saidctb]

I don't remember having @lambdify I was thinking about lambdas

[EmilyBourne]

It is defined here:

pyccel/pyccel/decorators.py

Lines 25 to 34 in b39fe9a

	def lambdify(f):
	args = f.__code__.co_varnames
	from sympy import symbols
	args = symbols(args)
	expr = f(*args)
	def wrapper(*vals):
	return expr.subs(zip(args,vals)).doit()
	return wrapper

Like I say it looks like your description of what @sympy should do. Sympy's lambdify converts a sympy expression into a function which also seems to be what you are describing above. So there is definitely some redundancy between these two

[EmilyBourne]

But there is an example here: https://github.com/pyccel/pyccel/blob/devel/tests/symbolic/scripts/lambdas.py#L6-L11
Which I do not understand at all as g1 should already be a function

[saidctb]

yes, I don't think we need it, we should stick with the more general @sympy

[saidctb]

But there is an example here: https://github.com/pyccel/pyccel/blob/devel/tests/symbolic/scripts/lambdas.py#L6-L11 Which I do not understand at all as g1 should already be a function

Yes, its behavior does not make sense in Python, we do not handle lambdas correctly,
we instead assume that a lambda expression is a symbolic expression, which is False in Python

[EmilyBourne]

However, this mechanism was removed in the issue #668 because it was not handled properly.

This is not why it was removed. It was also removed because we couldn't understand the difference between @python and @sympy. We could think of use cases for @sympy (most of which looked like the examples you provide above) but we couldn't work out what cases additionally required a @python decorator. Of the two we then picked the one with the name that was closest to being something understandable. I prefer @precompute as a name for this decorator

[saidctb]

their behaviour is different because one returns a symbolic expressions as explained above

[EmilyBourne]

All your examples seem to be about testing. Why do you need the tests to be translated?

[saidctb]

There are other cases that I can provide, We can use it for example to prepare some data by reading them from a file and manipulating them using different libraries and then use their result to do some computations, something like

@precompute
def prepare_data():
    ...
    return A
def manipulate_data(A:'float[:]'):
    ...
A = prepare_data()
manipulate_data(A)

This will allow us to generate an executable that is more portable contrary to having a shared library.
This is needed generally in Embedded programming where you have limited memory and you need a small executable

[EmilyBourne]

At the moment @sympy does not work in Python because it behaves differently we should wrap the function that has a sympy decorator and evaluate its results by doing

def sympy(f):
    def wrap(x):
        expr = f(x)
        return expr.subs(Symbol('x'), x)
     return wrap
``

This implementation requires:

1. That f take the same arguments as the sympy expression
2. That the names of the arguments to f and the arguments in the sympy expression match
3. That the arguments passed to f not be used.

Why not implement the decorator like:

import sympy as sp
import inspect

def sympy(f):
    sig = inspect.signature(f)
    args = [sp.Symbol(a) for a in sig.parameters]
    return sp.lambdify(args, f(*args))

with the function defined as:

@sympy
def f(x : float):
    from sympy import tan
    return tan(x).diff()

?

[saidctb]

yes, we can do it like that

[saidctb]

Another use case is to store the results of the factorial functions in an array and use the result in another function to improve performance.

import numpy as np

# Function to compute factorial
def factorial(n):
    return 1 if n <= 1 else n * factorial(n - 1)

# Function to compute binomial coefficient (n choose k)
def binomial_coefficient(n, k):
    if n<20:
        fact_n = fact_array[n]
    else:
        fact_n = factorial(n)

    if n<20:
        fact_k = fact_array[k]
    else:
        fact_k = factorial(k)

    if n-k<20:
        fact_nk = fact_array[n-k]
    else:
        fact_nk = factorial(n-k)
    # Compute binomial coefficient using the computed factorials
    return fact_n // (fact_k * fact_nk)

@precompute
def compute_first_n_factorials(n):
    from functools import lru_cache
    fact = lru_cache(factorial)
    return np.array([fact(i) for i in range(n)])

n:'Final'  = 20
fact_array = compute_first_n_factorials(n)

if __name__ == "__main__":
    n = 10
    k = 3
    binomial_coeff = binomial_coefficient(n, k)
    print(f"Binomial coefficient C({n}, {k}) is: {binomial_coeff}")
``

[saidctb]

@EmilyBourne , you also have this example

[saidctb]

This case cannot be replicated by the use of @sympy or lambdify because both work on symbolic expressions

[EmilyBourne]

Summary of related decorators:

@python

Equivalent to the proposed @precompute or the C++ keyword constexpr. May help reduce the size of the executable. According to the discussion on #668 we decided that this was not a priority at the moment so we need a strong use case to change this opinion as there are lots of high priority issues in progress.

Possible use cases:

1. Embedded programming
2. binomial_coefficients

@sympy

Allows the execution of Python code to obtain a symbolic expression. Creates a function which is described by the symbolic result of the original function.
E.g:

@sympy 
def f(x:'float'):
    from sympy import Symbol, tan
    x_symbol = Symbol('x')
    return tan(x_symbol).diff()

should generate:

@sympy
def f(x:float):
    import sympy as sp
    return 1/(sp.cos(x)**2)

This is useful but, as written, there are several restrictions which must be tested during the translation. Namely:

1. That f take the same arguments as the sympy expression
2. That the names of the arguments to f and the arguments in the sympy expression match
3. That the arguments passed to f not be used.

With these restrictions we can write the decorator as:

def sympy(f):
    def wrap(*args):
        expr = f(*args)
        sig = inspect.signature(f)
        subs = {sp.Symbol(a_name): a for a_name, a in zip(sig.parameters, args)}
        return expr.subs(subs)
     return wrap

I am not sure if this version allows arrays to be passed. It should be reasonably costly to evaluate expr in this way.

Alternatively we could modify the implementation slightly. E.g:

@sympy 
def f(x:'float'):
    from sympy import tan
    return tan(x).diff()

should generate:

def f(x:float):
    import numpy as np
    return 1/(np.cos(x)**2)

This removes the restrictions mentioned above. The decorator can then be written as:

import sympy as sp
import inspect

def sympy(f):
    sig = inspect.signature(f)
    args = [sp.Symbol(a) for a in sig.parameters]
    return sp.lambdify(args, f(*args))

@lambdify

This decorator should be removed. Its name is highly misleading and it relies on the code not being executable. An example can be found in the tests:

pyccel/tests/symbolic/scripts/lambdas.py

Lines 4 to 11 in b39fe9a

	f1 = lambda x: x**2 + 1
	f2 = lambda x,y: x**2 + f1(y)*f1(x)
	g1 = lambda x: f1(x)**2 + 1
	# lambda expressions can be printed
	#$ header m1(double)
	m1 = lambdify(g1)

Here we assume that a lambda function does not work like a lambda function in Python. Rather it works like a sympy expression. The lambdify function (despite being defined in decorators.py it seems this is a function not a decorator) is then used like sympy's lambdify function but on these pseudo-lambda functions.

Currently the implementation of the decorator matches the suggested implementation for @sympy.

[saidctb]

Even if it was not removed, the sympy decorator returns a function that can be evaluated, the precompute does the evaluation

[saidctb]

This is the difference

@precompute
def f1(x):
    return x**3
@sympy 
def f2(x):
    return Symbol('x')**3

a1 = f1(2.)
a2 = f2(2.)

This should be translated into something like

a1 = 8.
a2 = f2(2.)
x2 = 3.14
a3 = f2(x2)
#f1(x2) will raise an error, because we don't know the value of x2

one evaluates the function with arguments passed as literals, the other transforms a symbolic expression that is the result of the function into a callable function

[EmilyBourne]

I understood, you explained this above. But as discussed here, as long as lambdify still exists it gives the impression that we should have:

@precompute
def f1(x):
    return x**3
@sympy
def f2(y):
     from sympy import x, lambdify
     f = x**3
     return f.subs(x,y)
@lambdify
def f3(x):
     return x**3
a1 = f1(2.)
a2 = f2(2.)
a3 = f3(2.)

which should be translated into something like

a1 = 8.
a2 = 8.
a2 = f3(2.)

which is why I say "if the lambdify decorator is removed". Not because it actually changes anything but because its presence and misleading name makes @sympy ambiguous

[saidctb]

lambdify does not execute the code contrary to sympy and precompute , with sympy and precompute we call the Python interpreter to execute the function, with lambdify we do the analysis ourselves and generate the lambdified function

[saidctb]

Since the beginning lambdify has been different from the others, the issue is that lambdify is restricted compared to @sympy so there is no need to keep it

[saidctb]

Another Example in the domain of finite elements.
We want to perform the analysis on a specific mesh instead of writing something general that works on any mesh
which is slower in general.

# Function to read node and element information from a mesh file
@precompute
def read_mesh(filename):
    ...
    return nodes, elements

# Function that performs Finite Element Analysis using node and element data
def perform_fea(parameters):
   # Read node and element information from the mesh file
    filename: 'Final' = 'mesh_data.txt'  # Provide the filename of the mesh data file
    nodes, elements = read_mesh(filename)
    # Perform FEA calculations using the provided node and element data
    print("Performing Finite Element Analysis...")
    # Placeholder for actual FEA calculations
    # Access nodes and elements to perform calculations

if __name__ == '__main__':
    # Perform Finite Element Analysis
    parameters = ...
    perform_fea(parameters)
``

[saidctb]

I think it is fine because this is a special case, if the user does not like it he would not need to use it
The whole reason why we want to precompute some stuff at the compilation stage is to optimize the code
and generate a less general code, which will allow the Fortran compiler to further optimize it.
Also this gives the user more freedom to use other libraries that we do not support

[EmilyBourne]

Also this gives the user more freedom to use other libraries that we do not support

It is exactly this case that I am thinking of where you may not want the result to be const

[saidctb]

can you give me an example where the use of const result will be a problem?

[EmilyBourne]

@precompute
def complex_function_using_library_not_supported_by_pyccel():
     ...
     return generated_result

starting_point = complex_function_using_library_not_supported_by_pyccel()
starting_point += dt

[saidctb]

okay then the user will need to do this to force the use of Final

degree:'Final',nodes, elements = read_mesh(...)
``

[EmilyBourne]

I am checking the old tests to make sure I have fully understood. I see this:

pyccel/tests/symbolic/scripts/decorator.py

Lines 21 to 36 in b39fe9a

	@sympy
	def f1_sympy():
	#This function must return a sympy expression
	#that depends on the arguments of the function
	from sympy.abc import x,y
	expr = x*y + 2
	return expr
	@sympy
	def f2_sympy(x):
	#This function must return a sympy expression
	#that depends on the arguments of the function
	from sympy import diff
	u = x**2 + 2*x
	u_dx = diff(u, x)
	return u_dx

You have 1 function with 1 argument which matches my description above, but also a second function f1_sympy which takes no arguments. Do you know how this second case was supposed to work?

[saidctb]

The behavior here was not clear, because we didn't care about the code working in Python, but the result was always the same, In the case where we don't have arguments we deduced them from the expression

[saidctb]

For me, it seems that the behavior of the @precompute is clear and distinct from @sympy, and it is needed in many cases to optimize the code. unless there is a clear reason not to support it, we should implement it.

[EmilyBourne]

For me, it seems that the behavior of the @precompute is clear and distinct from @sympy

I agree.

and it is needed in many cases to optimize the code

I disagree. It is helpful but not indispensable.

unless there is a clear reason not to support it, we should implement it.

I agree but I would like to see Yaman's opinion on this (he is on holiday this week).

That said I don't see this as a higher priority issue than the currently ongoing work. If you want to take care of this and Yaman agrees then I think we can create a dedicated issue but the PR that adds this will also need to add clear and detailed documentation in this file:

https://github.com/pyccel/pyccel/blob/devel/docs/decorators.md

[saidctb]

Sure, this is a new feature, it will improve the performance of the code, I don't see it differently than the @sympy in terms of importance.
but yes, I can take care of it when I have time

[EmilyBourne]

Sure, this is a new feature, it will improve the performance of the code, I don't see it differently than the @sympy in terms of importance.

This I agree with completely. Fixing @sympy is also a long way down my priority list 😅

[saidctb]

Another use case where we can extract the length of the SIMD vector using a Python function that calculates it depending on the architecture.
This will allow the user to mimic the notion of macros in C/C++ and create a code specific to the architecture.

@precompute
def get_simd_length():
    # compute the simd length specific to the architecture 
    import platform

    def calculate_simd_vector_total_size(architecture):
        if 'x86' in architecture:
            # Assuming AVX architecture
            vector_register_size_bits = 256
        elif 'arm' in architecture:
            # Assuming ARM NEON architecture
            vector_register_size_bits = 128
        else:
            raise ValueError("Unsupported architecture")

        return vector_register_size_bits

    # Detect architecture
    architecture = platform.machine()

    # Calculate total size of SIMD vector
    total_size_bits = calculate_simd_vector_total_size(architecture)
    return total_size_bits//64

def performe_computation(y:'float[:]')
    simd_lenght:'Final' = get_simd_lenght()
    x:'Final' = np.empty(simd_length)
    for i in range(y.shape[0]//simd_length):
        x[:] = y[i*simd_lenght:(i+1*simd_length)]
        #$ omp simd
        for j in range(simd_length):
            ...

[saidctb]

@EmilyBourne , you have this case where we mimic the notion of macros in C/C++, and generate code specific to the architecture

[saidctb]

Another example is to check if a library exists and use it otherwise refer to manual implementation

@precompute
def check_blas():
    """
    Checks if BLAS is available.
    
    Returns:
        bool: True if BLAS is available, False otherwise.
    """
    try:
        import scipy.linalg.blas
        return True
    except ImportError:
        return False

def _blas_matrix_vector_product(A, x):
    from scipy.linalg.blas import dtrmv
    return dtrmv(A, x)

def _manual_matrix_vector_product(A, x):
    """
    Performs matrix-vector product Ax using manual implementation.
    
    Parameters:
        A (ndarray): Matrix of shape (m, n).
        x (ndarray): Vector of shape (n,).
        
    Returns:
        ndarray: Resulting vector of shape (m,).
    """
    ...

def matrix_vector_product(A, x):
    """
    Performs matrix-vector product Ax using BLAS if available; otherwise, uses manual implementation.
    
    Parameters:
        A (ndarray): Matrix of shape (m, n).
        x (ndarray): Vector of shape (n,).
        
    Returns:
        ndarray: Resulting vector of shape (m,).
    """
    if check_blas():
        return _blas_matrix_vector_product(A, x)
    else:
        return _manual_matrix_vector_product(A, x)
    ```