Introducing Symbolic Binary operators and i32 to S Casting function

[anutosh491]

Towards #1907 and #1846 (comment)

What all has been addressed

1. SymbolicAdd namespace was removed and macro for supporting all symbolic binary operators was introduced.
2. The S() casting function was introduced .
   A random example would be

(lf) anutosh491@spbhat68:~/lpython/lpython$ cat examples/expr2.py 
from sympy import Symbol

def main0():
    x: S = Symbol('x')
    y: S = Symbol('y')
    z: S = x + y
    a: S = S(20)
    b: S = z ** a
    c: S = a / z
    print(b)
    print(c)

main0()

# Not implemented yet in LPython:
#if __name__ == "__main__":
#    main()
(lf) anutosh491@spbhat68:~/lpython/lpython$ lpython --backend=c --enable-symengine examples/expr2.py 
(x + y)**20
20/(x + y)

3. Intial Implementation for SymbolicExpression ttype  #1846 (comment)

[anutosh491]

Once I have an initial review and I get an idea of what sort of tests would be most appropriate here , I'll add few integration tests in this PR . I guess that would be enough for the PR , I will then work towards chaining of operators and introducing a few elementary functions like sin , cos , abs in the following PR.

[certik]

I would say go ahead and add tests for all the new features that you implemented.

I also recommend to split PRs into two kinds: "new features but no refactorings" and "refactorings but no new features". That way it's easier to review and see if anything broke. Right now we need a lot more tests to check every corner case. Then we can have more confidence in refactoring that nothing broke.

[anutosh491]

I also recommend to split PRs into two kinds: "new features but no refactorings" and "refactorings but no new features".

I will surely keep this in mind for future work . I will add tests soon .

[anutosh491]

I am actually concerned about two things (which could be addressed in subsequent prs or this one ) @certik

As we went through Symengine's C wrapper file, most of the functions defined take in basic variable/variables and assign it to another basic variable like basic_add etc.

1. So while dealing with integers there are few main functions

//! Assign to s, a long.
CWRAPPER_OUTPUT_TYPE integer_set_si(basic s, long i);
//! Assign to s, a ulong.
CWRAPPER_OUTPUT_TYPE integer_set_ui(basic s, unsigned long i);
//! Assign to s, an integer that has base 10 representation c.
CWRAPPER_OUTPUT_TYPE integer_set_str(basic s, const char *c);

So a program like the following would work because we would have integer_set_si(y, 20); in the C generated code

def test_symbolic_operations():
    x: S = Symbol('x')
    y: S = S(20)
    z: S = x + y
    print(z)

But something like the following might not because there is no variable to assign the integer to

def test_symbolic_operations():
    x: S = Symbol('x')
    z: S = x + S(20)
    print(z)

Well we could have a workaround for this using the following function void basic_const_set(basic s, const char *c);

void test_symbolic_operations() {
    basic x, z;
    basic_new_stack(x);
    basic_new_stack(z);
    basic_symbol_set(x, "x");

    basic_const_set(z, "20");
    basic_add(z, x, z);

    char* result = basic_str(z);
    printf("%s\n", result);

    basic_str_free(result);
    basic_del_stack(x);
    basic_del_stack(z);
}

where we assign 20 as a constant to z and then assign x + z to z

[anutosh491]

2. But I'm not sure if we have a workaround for

def test_symbolic_operations():
    z: S = S(30) + S(20)
    print(z)

Here we might need a dummy variable I guess !?

[certik]

This PR looks good, can you please rebase it on top of the latest main? (resolving conflicts) Then I'll do a final review and merge it.

Yes, I think you need to create a temporary variable when traversing the expression, to be able to handle things like z = (x+y+z+x*y)**2 etc.

I would also get assert working, so that you can start testing it like this:

assert z == S(0)

To ensure things work (currently we only print it, so if it prints incorrectly, the test would not reveal it).

[certik]

Here is how to handle expression evaluation robustly:

diff --git a/src/libasr/codegen/asr_to_c_cpp.h b/src/libasr/codegen/asr_to_c_cpp.h
index 903ed5eee..cb9cc0d41 100644
--- a/src/libasr/codegen/asr_to_c_cpp.h
+++ b/src/libasr/codegen/asr_to_c_cpp.h
@@ -83,6 +83,34 @@ struct CPPDeclarationOptions: public DeclarationOptions {
     }
 };
 
+class SymEngineStack {
+    std::vector<std::string> stack;
+    size_t stack_top = -1;
+    std::string &symengine_src;
+
+    SymEngineStack(std::string symengine_src) : symengine_src{symengine_src} {}
+
+    std::string pop() {
+        LCOMPILERS_ASSERT(stack_pop > 0 && stack_pop < stack.size())
+        std::string top = stack[stack_top];
+        stack_top--;
+        return top;
+    }
+
+    std::string push() {
+        if (stack_top < stack.size()) {
+            // Create a unique name for a symengine variable:
+            std::string var = "stack" + std::to_string(stack.size());
+            stack.push_back(var);
+            symengine_src += "basic " + var + ";";
+            symengine_src += "symengine_basic_new_stack(" + var + ");";
+        }
+        stack_top++;
+        return stack[stack_top];
+    }
+
+}
+
 template <class Struct>
 class BaseCCPPVisitor : public ASR::BaseVisitor<Struct>
 {
@@ -105,6 +133,9 @@ public:
     std::map<uint64_t, std::string> const_var_names;
     std::map<int32_t, std::string> gotoid2name;
     std::map<std::string, std::string> emit_headers;
+    std::string symengine_src;
+    SymEngineStack symengine_stack(symengine_src);
+
 
     // Output configuration:
     // Use std::string or char*
@@ -2384,7 +2415,6 @@ R"(#include <stdio.h>
             SET_INTRINSIC_NAME(Exp, "exp");
             SET_INTRINSIC_NAME(Exp2, "exp2");
             SET_INTRINSIC_NAME(Expm1, "expm1");
-            SET_INTRINSIC_NAME(SymbolicSymbol, "Symbol");
             SET_INTRINSIC_NAME(SymbolicInteger, "Integer");
             SET_INTRINSIC_NAME(SymbolicPi, "pi");
             case (static_cast<int64_t>(ASRUtils::IntrinsicFunctions::SymbolicAdd)):
@@ -2394,11 +2424,21 @@ R"(#include <stdio.h>
             case (static_cast<int64_t>(ASRUtils::IntrinsicFunctions::SymbolicPow)): {
                 LCOMPILERS_ASSERT(x.n_args == 2);
                 this->visit_expr(*x.m_args[0]);
-                std::string arg1 = src;
                 this->visit_expr(*x.m_args[1]);
-                std::string arg2 = src;
-                out = arg1 + "," + arg2;
-                src = out;
+                std::string arg2 = symengine_stack.pop()
+                std::string arg1 = symengine_stack.pop()
+                std::string target = symengine_stack.push()
+                symengine_src += "symengine_add(" + target + "," + arg1 + ","
+                    + arg2 + ");\n";
+                src = target;
+                break;
+            }
+            SET_INTRINSIC_NAME(SymbolicSymbol, "Symbol"); {
+                LCOMPILERS_ASSERT(x.n_args == 1);
+                this->visit_expr(*x.m_args[0]);
+                std::string target = symengine_stack.push()
+                symengine_src += "symengine_symbol(" + target + "," + src + ");\n";
+                src = target;
                 break;
             }
             default : {

There are two issues to fix:

• symengine_src must be eventually emitted to src
• the stack must be emptied when the function exits, so that we start creating new variables for the next function (we can't reuse variables from the old function)

The first issue is fixed by:

--- a/src/libasr/codegen/asr_to_c.cpp
+++ b/src/libasr/codegen/asr_to_c.cpp
@@ -1233,6 +1233,10 @@ R"(    // Initialise Numpy
         for (size_t i=0; i<x.n_values; i++) {
             this->visit_expr(*x.m_values[i]);
             ASR::ttype_t* value_type = ASRUtils::expr_type(x.m_values[i]);
+            if (value_type->type == ASR::ttypeType::SymbolicExpression) {
+                out += symengine_src;
+                symengine_src = "";
+            }
             if( ASRUtils::is_array(value_type) ) {
                 src += "->data";
             }

[certik]

I think that this looks good as is. Thanks!

[anutosh491]

I'm seeing a rather queue based implementation , I'll try my best to break my thinking down for you . Let's try to solve this program

from sympy import Symbol, pi
def main0():
    x: S = Symbol('x')
    y: S = Symbol('y')
    z: S = pi + x + y
main0()

The C generated code for this should be

void main0()
{
    basic x;
    basic_new_stack(x);
    basic y;
    basic_new_stack(y);
    basic z;
    basic_new_stack(z);
    symbol_set(x, "x");
    symbol_set(y, "y");
    basic stack0;
    basic_new_stack(stack0);
    basic stack1;
    basic_new_stack(stack1);
    basic_const_pi(stack1);
    basic_add(stack0, stack1, x);   
    basic_add(z, stack0, y);
}

I'll discuss the logic in the comment below .

[anutosh491]

Let's think for a while that we are using a container and not a stack .

1. As we traverse the ASR , the first thing we deal with are the annotation statements (x: S , y: S, z: S) . This generates the following

    basic x;
    basic_new_stack(x);
    basic y;
    basic_new_stack(y);
    basic z;
    basic_new_stack(z);

By the end of this our container would have [x, y, z].

2. Once we are done with handling annotations , we essentially deal with assignments & other operations and this would be in the same order of how annotations work. Hence we would first deal with assignments related to x , then y, then z. So for the following

    symbol_set(x, "x");
    symbol_set(y, "y");

We would need the first two elements of the container [x, y, z], so we pop them out . Now the container has [z]

3. Coming to the last part which is the addition operation . As per hat I notice pi + x + y would be dealt as ((pi+ x) + y) i.e left-associative
   i) The first element would end up being the target . Container has [z], this element would be target , we pop it out
   ii) We visit pi + x . We need a target but container is empty , which means we need a temporary variable stack0 .We push this into the container ([stack0]) and use it as target by popping this out.
   iii) Then we visit pi . Container is empty , we use temporary variable stack1 . This would generate

    basic stack0;
    basic_new_stack(stack0);
    basic stack1;
    basic_new_stack(stack1);
    basic_const_pi(stack1);

iV) Finally we push x and y in our container. And we can use that to get

    basic_add(stack0, stack1, x);   
    basic_add(z, stack0, y);

So basically we are splitting pi + x + y to stack0 + y where stack0 is stack1 + x where stack1 is pi .

[anutosh491]

I do think this would work , I'll give it a shot and get back to you . As of now , I am seeing positive results .

[certik]

@anutosh491 go ahead and see if this works. Make sure you also test a case like (x+y+z) + (x+z) and x+(y+(z+x)) and similar. Once you implement mul, do (x+y+z)*(x+z)*(y+z)*(x+y) etc.

[anutosh491]

Yes, I tried these complex cases out (which are nested with brackets like the ones you've mentioned above).
I can see that the approach is working for both

1. right recursive blocks - x+(y+(z+x))
2. left recursive blocks - (((x + y) + z) + x)
3. Misc cases - (x+y+z) + (x+z)

I just have a minute doubt that I'll point out on the pr but other I feel the approach is good . I'll raise a pr soon !

[namannimmo10]

Ahh! Adding symbolics to LPython already! 🚀

[certik]

@namannimmo10 absolutely. We are tying up all loose ends, this addition makes LPython complete in terms of scope and I think it will be a big hit for SymPy.

[namannimmo10]

That's awesome! I would love to see LPython supporting both numerics and symbolics with fast compilation. Would be incredible to use over the years!
My team is working on GPU methods and using CuPy. Since Python is slow, we write raw CUDA kernels to accelerate the methods. The speed of LPython (with --fast) is much faster than C, so it would be great if LPython could compile CuPy code at some point. So basically you won't need to deal with low level code directly at all to accelerate the methods - everything would happen in the backend, simply write code in Python, compile it with LPython fast and get the same speed.

[certik]

Yes, good idea. I opened up an issue for this: #2075, if you have some CuPy code that you could share, that would be great.

[anutosh491]

Adding support for compiling CuPy code would be awesome @namannimmo10