An experimental Python-to-C transpiler and domain specific language for embedded high-performance computing

Hi!

Julia 1.7 will feature a mechanism called "BLAS trampoline", which allows for loading a BLAS library at runtime (see here). This PR makes use of that and adds MKL as a dependency (plus project management/instantiation stuff). Someone running the benchmark doesn't need to have MKL installed already, it's installed during the activation of the julia project before the benchmarks are run. This will only work with julia 1.7+ though (which should be released relatively soon), so if that's a concern I can add a note somewhere saying that the benchmarks require it.

I've also fixed the benchmarking script, since the existing version in the repo tried to run a test_riccati.jl.in file, which didn't exist. I've also made it so that the NM and NREP parameter is passed on the command line instead of hardcoding via a file.

With this and the removal of some unnecessary (I think) copies, I get these results locally:

Note that I've only run python3 run_benchmark_julia.py as well as python3 run_benchmark_numpy.py, as I don't know how to run your BLASFEO benchmarks. I suspect the remaining difference in benchmarking times is due to this, as well as my laptop not being as fast as the machine you've run your benchmarks on. Nevertheless, since Julia+MKL can match performance of prometeo/BLASFEO once the matrix size becomes large enough i.e. we're memory bound, I suspect the speedup is real.

Cheers!

Hi! I was looking forward to installing the module but came across this error while trying to install the module on Python v3.9.10:

C:\Users\HP>py -3.9 -m pip install prometeo-dsl
ERROR: Ignored the following versions that require a different python version: 0.0.1 Requires-Python >=3.6, <=3.9; 0.0.10 Requires-Python >=3.6, <=3.9; 0.0.2 Requires-Python >=3.6, <=3.9; 0.0.3 Requires-Python >=3.6, <=3.9; 0.0.4 Requires-Python >=3.6, <=3.9; 0.0.5 Requires-Python >=3.6, <=3.9; 0.0.6 Requires-Python >=3.6, <=3.9; 0.0.7 Requires-Python >=3.6, <=3.9; 0.0.9 Requires-Python >=3.6, <=3.9
ERROR: Could not find a version that satisfies the requirement prometeo-dsl (from versions: none)
ERROR: No matching distribution found for prometeo-dsl

The screenshot for the same:

Hey, after reading a little bit about this project, I love it but I'm a little confused one a few parts; does it transpile to exe? Or what does it transpile the python code to so users can run? Also; would I have to write my python programs in such a way so prometeo is able to access and execute it?

TY.

Sorry if I've missed something, but is there a method you have to allow tests to run on the transpiled C to ensure it's working? (and presumably that it gives the same results as the python code, unless you've got some way of being 100% sure that is guaranteed due to your language restrictions)

This PR implements prometeo's heap computation based on the meta-information stored by the parser and Dijkstra's algorithm on the call graph that computes the worst-case heap usage. In order to do this, first, the meta-information is used to resolve the values of dims and dimv variables.

This is part of a brainstorming with the development team on how to define meta-information structures for classes and methods within prometeo's parser. Rough draft:

Python example

class ClassA:
    def __init__(self):
        self.attr1 : int = 0

    def method1(self, arg1 : int) -> int:
        return self.attr1 + arg1

class ClassB:
    def __init__(self):
        self.attr1 : int = 0
        self.attr2 : ClassA = ClassA()

    def method1(self, arg1 : int) -> ClassA:
        return self.attr2

Meta-information:

[
    'ClassA' : [
        'constructor_args' : 
            [('self', 'self')],
        'attributes' : 
            [('attr1', 'int', 0)],
        'methods' : [
            'method1' : 
                'args' : 
                [('self','self'), ('arg1', 'int')],
                'return_type' : 'int'
        ]
    ]
]

[
    'ClassB' : [
        'constructor_args' : 
            [('self', 'self')],
        'attributes' : 
            [('attr1', 'int', 0), ('attr2', 'ClassA', 'None')],
        'methods' : [
            'method1' : 
                'args' : 
                [('self','self'), ('arg1', 'int')],
                'return_type' : 'ClassA'
        ]
    ]
]

This PR is a first step towards improved type checking, which is currently a bit messy especially for nested attributes, BinOps and UnaryOps, Lists and so on. The PR also contains preliminary work on CasADi code-gen.

replacing https://github.com/zanellia/prometeo/pull/9 - the implementation of instance attributes required quite some fixes to the heap computation machinery.

This PR implements prometeo's heap computation based on the meta-information stored by the parser and Dijkstra's algorithm on the call graph that computes the worst-case heap usage. In order to do this, first, the meta-information is used to resolve the values of dims and dimv variables. In summary, the heap computation works (roughly) as follows:

1. prometeo's parser goes through the Python code and builds a record of the constructors being defined and called: these are the only functions that give rise to memory escape in the DSL. pmat() and pvec() are the only built-in constructors, so the heap usage (64- and 8-byte aligned) is directly computed. When a constructor associated to a user class MyClass is defined the memory footprint is added to the heap usage of its __init__ method and MyClass() is registered as a constructor.

2. after the Python code has been transpiled to C, another parser goes thorough the Python AST once more and builds a call graph, and, from the call graph, a reachability map, that describes which function scopes can be reached from every function scope - cycles are allowed, as long as there is no memory allocation in them.

3. The worst case heap usage is computed by calculating the shortest path on a properly constructed graph (slightly modified version of the call graph) whose vertexes represent function scopes and whose edges are given by (minus) the memory allocated within each function. This graph is built such that, whenever a method that is registered as constructor at point 1 is called, the amount of memory associated with that scope is added to the caller's memory usage (i.e. accounting for memory escape)

Addressing https://github.com/zanellia/prometeo/issues/6. The AST parser now builds an additional data structure called meta_info. For example

from prometeo import *

n : dims = 10

class ClassA():
    A : pmat = pmat(n,n)

    def method1(self, arg1 : int) -> int:
        return 0

class ClassB():
    attr1 : ClassA = ClassA()

class ClassC():
    attr2 : ClassB = ClassB()

def main() -> int:
    A : pmat = pmat(n,n)
    D : ClassC = ClassC()
    return 0

generates

{
   "global": {},
   "global@ClassA": {
      "attr": {
         "A": "pmat"
      },
      "methods": {
         "method1": {
            "args": {
               "arg1": "int"
            },
            "return_type": "int"
         }
      }
   },
   "global@ClassB": {
      "attr": {
         "attr1": "ClassA"
      },
      "methods": {}
   },
   "global@ClassC": {
      "attr": {
         "attr2": "ClassB"
      },
      "methods": {}
   }
}

Moving to numpy-like syntax: A[n][m] -> A[n,m].

Additional changes:

• dump static analysis info to JSON files
• implement factorizations and solves (LU and Cholesky)
• fix assignments across pmats, pvecs and scalars

macOs BigSur 11.4 MacBookPro 13-inch, M1, 2020

I cloned the repo, installed with pip install -e ., ran ipython, and tried import prometeo. I get the error:

Python 3.8.8 (default, Apr 13 2021, 12:59:45)
Type 'copyright', 'credits' or 'license' for more information
IPython 7.22.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: import prometeo
---------------------------------------------------------------------------
OSError                                   Traceback (most recent call last)
<ipython-input-1-371e560288bc> in <module>
----> 1 import prometeo

~/repos/lib/prometeo/prometeo/__init__.py in <module>
      1 from . import cgen
----> 2 from . import linalg
      3 from . import mem
      4 from . import auxl
      5 from . import cmdline

~/repos/lib/prometeo/prometeo/linalg/__init__.py in <module>
----> 1 from .pmat_blasfeo_wrapper import *
      2 from .pvec_blasfeo_wrapper import *
      3 from .pmat import *
      4 from .pvec import *
      5 from .blasfeo_wrapper import *

~/repos/lib/prometeo/prometeo/linalg/pmat_blasfeo_wrapper.py in <module>
----> 1 from .blasfeo_wrapper import *
      2 from ctypes import *
      3
      4
      5 bw.blasfeo_dgeex1.restype = c_double

~/repos/lib/prometeo/prometeo/linalg/blasfeo_wrapper.py in <module>
      2 import os
      3
----> 4 bw = CDLL(os.path.dirname(__file__) + '/../lib/blasfeo/libblasfeo.so')
      5
      6 class blasfeo_dmat(Structure):

~/opt/anaconda3/lib/python3.8/ctypes/__init__.py in __init__(self, name, mode, handle, use_errno, use_last_error, winmode)
    379
    380         if handle is None:
--> 381             self._handle = _dlopen(self._name, mode)
    382         else:
    383             self._handle = handle

OSError: dlopen(/Users/vas/repos/lib/prometeo/prometeo/linalg/../lib/blasfeo/libblasfeo.so, 6): image not found

The current interface to BLAS and LAPACK routines should be made more flexible, in the sense that different levels of granularity should be accessible though the same method without adding unnecessary complexity to simple calls. Below is a draft signature for gemm suggested by @giaf:

gemm(A[1:2,1:2], B.T, C[3:4,1:2], D=None, alpha=1.0, beta=1.0) ,

where A, B, C and D can also be transposed through e.g. A.T and views are possible too, for example:

gemm(A[1:2,1:2], B.T, C[3:4,1:2], D, beta=2.0).

This is easily implementable in Python, but will require some work in prometeo's parser. In particular in the argument inspection here

https://github.com/zanellia/prometeo/blob/f45a6347ea524164cc812fad80c9db1f4f94e17f/prometeo/cgen/code_gen_c.py#L1826)

we'd need to do something like this:

nargs = len(args)

# extract keyword arguments
...

# check expression in views (if any) using https://github.com/zanellia/prometeo/blob/f45a6347ea524164cc812fad80c9db1f4f94e17f/prometeo/cgen/code_gen_c.py#L244

# extract and unparse views of matrix arguments (if any)
first_index = astu.unparse(node.targets[0].slice.value.elts[0]).strip('\n')
second_index = astu.unparse(node.targets[0].slice.value.elts[1]).strip('\n')

# check for transposition (look if arg has an attribute an if it is 'T')

# pass this info to a generic `process_arg()` method

As of now sized type checking is carried out in a rather case-specific fashion. At some point, the typed_record (https://github.com/zanellia/prometeo/blob/master/prometeo/cgen/code_gen_c.py#L324) and the var_dim_record (https://github.com/zanellia/prometeo/blob/master/prometeo/cgen/code_gen_c.py#L328) generated during the static analysis should be merged into a single dictionary and they should support any valid sized type e.g

n : dims = 10
A : pmat = pmat(n,n)

class ClassA():
    A : pmat = pmat(n,n)

class ClassB():
    A : pmat = pmat(n,n)
    B : ClassA = ClassA()

As of now 1) and 2) work fine, but to support 3) and more sophisticated user-defined sized types the code needs some restructuring.

In particular, if we analyze

from prometeo import *

n : dims = 10

class ClassA():
    A : pmat = pmat(n,n)

def main() -> int:
    A : pmat = pmat(n,n)
    return 0

it will result in the following typed_record:

{
    "global": {
        "n": "dims"
    },
    "global@ClassA": {
        "A": "pmat"
    },
    "global@main": {
        "A": "pmat"
    }
}

and the following var_dim_record:

{
    "global": {},
    "global@ClassA": {
        "A": [
            "n",
            "n"
        ]
    },
    "global@main": {
        "A": [
            "n",
            "n"
        ]
    }
}

High-priority

• [ ] Support for basic types: int, double, pmat and pvec
  • [x] declaration of variables of supported types
  • [x] specification of function arguments of supported types
  • [x] simple operations ('+', '-', '*', '/') on int and double types
• [ ] Support for high-level basic linear algebra operations on pmat variables
  • [x] C = A + B
  • [x] C = A - B
  • [x] C = A * B
  • [x] C = A\B
• [x] Support for subscripting of pmat variables
  • [x] get value (val = A[i,j])
  • [x] set value (A[i,j] = val)
• [x] Support for functions calls
• [x] Support for structures
  • [x] structure declaration
  • [x] structure access
• [ ] Support of arrays of all supported types
  • [x] array of pmat declaration through Python Lists
  • [x] array of other supported types declaration through Python Lists
  • [x] array of structures declaration through Python Lists
  • [x] accessing and setting arrays of pmats
• [x] Support for function definitions (no mem escape)
• [x] Extend AST parser to keep track of heap usage
• [ ] AST-based verification of user code
• [x] Unit testing

Low-priority/long-term

• [x] Support for simple classes
• [ ] Support for CasADi functions

Improvements

• [x] Remove dims from type hints in declarations (keep only in function signature)
• [x] Add linear algebra routines with explicit return
• [x] Add support for slicing