Can't figure out how polynomial conversion works in the M2/Sage interface

[saad khalid]

Hey everyone:

I was hoping some of you could provide some insight where my knowledge is lacking. I'm trying to add to the M2/Sage interface by adding a conversion for the M2 Divide class. Sage can already convert polynomials, so my hope was to just have it treat the Divide class as two polynomials, one in the numerator and the other in the denominator. Here is an example of an object in the divide class:

macaulay2.eval("""
K = toField(QQ[zet]/(zet^8 - zet^7 +zet^5 - zet^4 +zet^3 -zet + 1))
A=matrix{{zet^1,0},{0,zet^13}}
needsPackage "InvariantRing"
G=generateGroup({A},K)
P = molienSeries G
 """)

P is a "Divide" object. Here is an example of Sage converting a polynomial from M2 to Sage:

macaulay2.eval("""
needsPackage "Points";
M = matrix{{1,2,3},{4,5,6}}
R = QQ[x,y,MonomialOrder=>Lex];
(Q,inG,G) = points(M,R)
G#0
ring G#0
""")
G0 = macaulay2('G#0').to_sage(); G0

G in M2 is a list of polynomials, and G#0 is the first one. Here's what it looks like:

 3      2
y  - 15y  + 74y - 120

I convert this to a sage object, G0, and it appears as:

y^3 - 15*y^2 + 74*y - 120

The part that I really need is where Sage changes M2's two-line display for exponents on variables to just the standard Sage notation. I'm looking at the to_sage()
method in the source code but I can't seem to figure out how it's doing it. Am I missing something obvious? G#0 is a polynomial ring I believe, and here 
is the code in the to_sage function for PolynomialRings:
            elif cls_str == "PolynomialRing":
                from sage.rings.all import PolynomialRing
                from sage.rings.polynomial.term_order import inv_macaulay2_name_mapping

                #Get the base ring
                base_ring = self.coefficientRing().to_sage()

                #Get a string list of generators
                gens = str(self.gens())[1:-1]

                # Check that we are dealing with default degrees, i.e. 1's.
                if self.degrees().any("x -> x != {1}").to_sage():
                    raise ValueError("cannot convert Macaulay2 polynomial ring with non-default degrees to Sage")
                #Handle the term order
                external_string = self.external_string()
                order = None
                if "MonomialOrder" not in external_string:
                    order = "degrevlex"
                else:
                    for order_name in inv_macaulay2_name_mapping:
                        if order_name in external_string:
                            order = inv_macaulay2_name_mapping[order_name]
                if len(gens) > 1 and order is None:
                    raise ValueError("cannot convert Macaulay2's term order to a Sage term order")

                return PolynomialRing(base_ring, order=order, names=gens)

Where in this code does that conversion happen? It's like parsing ascii art...

Thanks!

[Dima Pasechnik]

On Saturday, July 2, 2016 at 7:13:41 PM UTC+1, saad khalid wrote:

Hey everyone:

I was hoping some of you could provide some insight where my knowledge is lacking. I'm trying to add to the M2/Sage interface by adding a conversion for the M2 Divide class. Sage can already convert polynomials, so my hope was to just have it treat the Divide class as two polynomials, one in the numerator and the other in the denominator.

Why? Sage can work with multivariate rational functions just fine:

sage: R.<x,y>=QQ[]

sage: f=x*y

sage: g=x^2+y^2

sage: f/g

x*y/(x^2 + y^2)

sage: type(f/g)

<type 'sage.rings.fraction_field_element.FractionFieldElement'>

sage: 

 

Here is an example of an object in the divide class:

macaulay2.eval("""
K = toField(QQ[zet]/(zet^8 - zet^7 +zet^5 - zet^4 +zet^3 -zet + 1))
A=matrix{{zet^1,0},{0,zet^13}}
needsPackage "InvariantRing"
G=generateGroup({A},K)
P = molienSeries G
 """)

M2 can provide you numerator and denominator of a Divide object, and then

you can use the existing Sage code to convert these two polynomials, and then form the fraction.

HTH,

Dmitrii

[saad khalid]

Well, here's what i tried and it doesn't seem to work after I build sage again. I added this code to the to_sage() function,

            elif cls_str == "Divide":
                self_Div = self
                div_Numerator = macaulay2('numerator self_Div')
                div_Denominator = macaulay2('denominator self_Div')
                div_Numerator = div_Numerator.sage_polystring()
                div_Denominator = div_Denominator.sage_polystring()
                sage_Div = div_Numerator/div_Denominator
                return sage_Div

 
I added this after the main else statement. The code looks like rubbish to me but I didn't know what else to do, so i thought I'd show something and ask for improvements, if that's okay... When running the following code now in Sage, here is my error message:

sage: macaulay2.eval("""
....: K = toField(QQ[zet]/(zet^8 - zet^7 +zet^5 - zet^4 +zet^3 -zet + 1))
....: A=matrix{{zet^1,0},{0,zet^14}}
....: needsPackage "InvariantRing"
....: G=generateGroup({A},K)
....: P = molienSeries G
....:  """)
K

PolynomialRing

| zet 0                                |
| 0   -zet^7+zet^6-zet^4+zet^3-zet^2+1 |

        2       2
Matrix K  <--- K

InvariantRing

Package

{| 1 0 |, | -zet^7-zet^2 0     |, | zet^7 0                             |, | zet^7-zet^6-zet^3+zet^2-1 0     |, | -zet^7+zet^6-zet^4+zet^3-zet^2+1 0   |, | zet 0                                |, | -zet^6-zet 0     |, | zet^5 0        |, | zet^2 0                        |, | zet^6 0                         |, | zet^7-zet^5+zet^4-zet^3+zet-1 0     |, | zet^4 0          |, | -zet^7+zet^5-zet^4-zet+1 0     |, | zet^3 0            |, | -zet^5-1 0     |}
 | 0 1 |  | 0            zet^3 |  | 0     zet^7-zet^5+zet^4-zet^3+zet-1 |  | 0                         zet^6 |  | 0                                zet |  | 0   -zet^7+zet^6-zet^4+zet^3-zet^2+1 |  | 0          zet^4 |  | 0     -zet^5-1 |  | 0     -zet^7+zet^5-zet^4-zet+1 |  | 0     zet^7-zet^6-zet^3+zet^2-1 |  | 0                             zet^7 |  | 0     -zet^6-zet |  | 0                        zet^2 |  | 0     -zet^7-zet^2 |  | 0        zet^5 |

List

         2    3    4    5    6    7    8    9    10    11    12    13    14
1 - T + T  - T  + T  - T  + T  - T  + T  - T  + T   - T   + T   - T   + T
---------------------------------------------------------------------------
        2          3    4    5    7    8           2           2    3    4
 (1 - T) (1 - T + T  - T  + T  - T  + T )(1 + T + T )(1 + T + T  + T  + T )

Expression of class Divide

sage: G = macaulay2('P').to_sage()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/home/saad/sage/local/lib/python2.7/site-packages/sage/all_cmdline.pyc in <module>()
----> 1 G = macaulay2('P').to_sage()

/home/saad/sage/local/lib/python2.7/site-packages/sage/interfaces/macaulay2.pyc in to_sage(self)
   1151             elif cls_str == "Divide":
   1152                 self_Div = self
-> 1153                 div_Numerator = macaulay2('numerator self_Div')
   1154                 div_Denominator = macaulay2('denominator self_Div')
   1155                 div_Numerator = div_Numerator.sage_polystring()

/home/saad/sage/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in __call__(self, x, name)
    241
    242         if isinstance(x, six.string_types):
--> 243             return cls(self, x, name=name)
    244         try:
    245             return self._coerce_from_special_method(x)

/home/saad/sage/local/lib/python2.7/site-packages/sage/interfaces/expect.pyc in __init__(self, parent, value, is_name, name)
   1379             except (RuntimeError, ValueError) as x:
   1380                 self._session_number = -1
-> 1381                 raise_(TypeError, x, sys.exc_info()[2])
   1382             except BaseException:
   1383                 self._session_number = -1

/home/saad/sage/local/lib/python2.7/site-packages/sage/interfaces/expect.pyc in __init__(self, parent, value, is_name, name)
   1374         else:
   1375             try:
-> 1376                 self._name = parent._create(value, name=name)
   1377             # Convert ValueError and RuntimeError to TypeError for
   1378             # coercion to work properly.

/home/saad/sage/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in _create(self, value, name)
    431     def _create(self, value, name=None):
    432         name = self._next_var_name() if name is None else name
--> 433         self.set(name, value)
    434         return name
    435

/home/saad/sage/local/lib/python2.7/site-packages/sage/interfaces/macaulay2.pyc in set(self, var, value)
    318         ans = Expect.eval(self, cmd)
    319         if ans.find("stdio:") != -1:
--> 320             raise RuntimeError("Error evaluating Macaulay2 code.\nIN:%s\nOUT:%s"%(cmd, ans))
    321
    322     def _object_class(self):

TypeError: Error evaluating Macaulay2 code.
IN:sage1=numerator self_Div;
OUT:stdio:13:7:(3): error: no method found for applying numerator to:
     argument   :  self    (of class IndexedVariable)
                       Div

It's certainly a different error message from what I was getting before when I tried this, but i don't know what it means.

[Dima Pasechnik]

On Sunday, July 3, 2016 at 1:21:18 AM UTC+1, saad khalid wrote:

Well, here's what i tried and it doesn't seem to work after I build sage again. I added this code to the to_sage() function,

            elif cls_str == "Divide":
                self_Div = self
                div_Numerator = macaulay2('numerator self_Div')
                div_Denominator = macaulay2('denominator self_Div')
                div_Numerator = div_Numerator.sage_polystring()
                div_Denominator = div_Denominator.sage_polystring()
                sage_Div = div_Numerator/div_Denominator
                return sage_Div

 

self_Div is not known to M2 in any way, yet you try calling in M2 'numerator self_Div'.

You need a way to get the counterpart of 'self' in M2 somehow.

[Dima Pasechnik]

On Sunday, July 3, 2016 at 11:05:06 AM UTC+1, Dima Pasechnik wrote:

On Sunday, July 3, 2016 at 1:21:18 AM UTC+1, saad khalid wrote:

Well, here's what i tried and it doesn't seem to work after I build sage again. I added this code to the to_sage() function,

            elif cls_str == "Divide":
                self_Div = self
                div_Numerator = macaulay2('numerator self_Div')
                div_Denominator = macaulay2('denominator self_Div')
                div_Numerator = div_Numerator.sage_polystring()
                div_Denominator = div_Denominator.sage_polystring()
                sage_Div = div_Numerator/div_Denominator
                return sage_Div

 

self_Div is not known to M2 in any way, yet you try calling in M2 'numerator self_Div'.

You need a way to get the counterpart of 'self' in M2 somehow.

Instead of

               self_Div = self
                div_Numerator = macaulay2('numerator self_Div')
                div_Denominator = macaulay2('denominator self_Div')
 

use

                div_Numerator = self.numerator()
                div_Denominator = self.denominator()

and it will be almost it (div_Denominator will be Product, and will need extra processing

because of this)

[saad khalid]

Hmm, so I tried this instead but it also didn't work and gave a similar error:

elif cls_str == "Divide":
        div_Str = repr_str
                div_Numerator = macaulay2('numerator div_Str')
                div_Denominator = macaulay2('denominator div_Str')

                div_Numerator = div_Numerator.sage_polystring()
                div_Denominator = div_Denominator.sage_polystring()
                sage_Div = div_Numerator/div_Denominator
                return sage_Div

I guess I don't understand what each of those variables contains... What is self in this situation? I thought that self would be the Divide object, since that is what's being passed to the to_sage function. But when that didn't work, i assumed that repr_str would contain the string of the Divide object, but that's also not working.

[Dima Pasechnik]

as I just posted, you can apply M2 functions to self.

self is of M2 type Divide, and so you can do

self.numerator()

etc.

[saad khalid]

Sorry! I hadn't seen that, that's great, thank you! The denominator is proving to be a bit tricky to deal with, mostly because I can't decide what I want to do with it. I did it though, but I'm not happy with it. I'm also an amatuer at regex(this is the first time I've used it, I figured it was about time I started learning now). But, I Think it should work? So, I made another method named sage_prodstring that does a similar job to sage_polystring. Here is the code for it:

    def sage_prodstring(self):
    """
        If this Macaulay2 element is a Product, return a string
        representation of this Product that is suitable for
        evaluation in Python. Needed internally for using to_sage on objects of class Divide.

        EXAMPLES::


        """
    external_string = self.external_String()
    prod_String = re.findall("new Power from \{(.+?),(.+?)\}", external_string)
    final_Prod = ""
    for i in range(0, len(prod_String)):
        final_Prod += "(" + prod_String[i][0] + ")" + '^' + prod_String[i][1] + "*"
    final_Prod = final_Prod[:-1]
    return final_Prod

I don't really know what to put in the Examples part. If you have any recommendations, I'd love to know. Also, I feel like this implementation is inneficient, since it first runs external_string, which is a Sage method that parses Macaulay's ascii art into a str, and then my method parses that str into something that you could actually use as python code. I feel like, ideally, I wouldn't need to run external_string, and could directly convert the ascii output from macaulay directly into the python code output I wanted, but i couldn't figure out exactly how to do that... I hope this is okay. I don't think i can test it yet until I fill in the spot for the Example though.

I then used this method in the to_sage() function for the denominator, as seen here:

            elif cls_str == "Divide":

                div_Numerator = self.numerator()
                div_Denominator = self.denominator()

                div_Numerator = div_Numerator.sage_polystring()
        div_Denominator = div_Denominator.sage_prodstring()
                sage_Div = "(" +  div_Numerator ")" + "/" + "(" + div_Denominator + ")"
                return sage_Div

How does it look? Thanks for all the help thus far.

[saad khalid]

So, after fixing a small typo, here is the output from an example:

sage: macaulay2.eval("""
....: K = toField(QQ[zet]/(zet^8 - zet^7 +zet^5 - zet^4 +zet^3 -zet + 1))
....: A=matrix{{zet^1,0},{0,zet^14}}
....: needsPackage "InvariantRing"
....: G=generateGroup({A},K)
....: P = molienSeries G
....:  """)

sage: macaulay2('P').to_sage()

'(1-T+T**2-T**3+T**4-T**5+T**6-T**7+T**8-T**9+T**10-T**11+T**12-T**13+T**14)/((1-T)^2*(1-T+T^3-T^4+T^5-T^7+T^8)^1*(1+T+T^2)^1*(1+T+T^2+T^3+T^4)^1)'

That's the output from the last command, anyways. Now, what I don't understand is how to make it into like... a function, or a Sage object or something like that. Like, I'd want to be able to compute its Taylor series just by doing something like "macaulay2('P').to_sage().taylor()" or something, but I'm not sure how to convert the string.

[leif]

saad khalid wrote:
> Sorry! I hadn't seen that, that's great, thank you! The denominator is
> proving to be a bit tricky to deal with, mostly because I can't decide
> what I want to do with it. I did it though, but I'm not happy with it.
> I'm also an amatuer at regex(this is the first time I've used it, I
> figured it was about time I started learning now). But, I Think it
> should work? So, I made another method named sage_prodstring that does a
> similar job to sage_polystring. Here is the code for it:
>
> |
> defsage_prodstring(self):
> """
> If this Macaulay2 element is a Product, return a string
> representation of this Product that is suitable for
> evaluation in Python. Needed internally for using to_sage on
> objects of class Divide.
>
> EXAMPLES::
>
>
> """
> external_string =self.external_String()
> prod_String =re.findall("new Power from
> \{(.+?),(.+?)\}",external_string)
> final_Prod =""
> fori inrange(0,len(prod_String)):

> final_Prod +="("+prod_String[i][0]+")"+'^'+prod_String[i][1]+"*"
> final_Prod =final_Prod[:-1]
> returnfinal_Prod
> |

Side note: Never use += on a string in a loop. Here it's even worse,
since you are concatenating strings on the right hand side as well.

Hint: "*".join(["alpha","beta","gamma"]) gives 'alpha*beta*gamma' for
example.

In your case the list to join is something like

[ "(%s)^%s" % (f[0], f[1]) for f in prod_String ]

so you could replace the five last lines by

return "*".join( [ "(%s)^%s" % (f[0], f[1]) for f in prod_String ] )


(Haven't tested that, your EXAMPLES or TESTS should do though...)


-leif

[Dima Pasechnik]

On Monday, July 4, 2016 at 6:09:59 AM UTC+1, saad khalid wrote:

So, after fixing a small typo, here is the output from an example:

sage: macaulay2.eval("""
....: K = toField(QQ[zet]/(zet^8 - zet^7 +zet^5 - zet^4 +zet^3 -zet + 1))
....: A=matrix{{zet^1,0},{0,zet^14}}
....: needsPackage "InvariantRing"
....: G=generateGroup({A},K)
....: P = molienSeries G
....:  """)

sage: macaulay2('P').to_sage()

'(1-T+T**2-T**3+T**4-T**5+T**6-T**7+T**8-T**9+T**10-T**11+T**12-T**13+T**14)/((1-T)^2*(1-T+T^3-T^4+T^5-T^7+T^8)^1*(1+T+T^2)^1*(1+T+T^2+T^3+T^4)^1)'

sage: R.<T>=QQ[]

sage: Rf=FractionField(R)

sage: f=Rf('(1+T)/(1-T)^2'); f # your string comes here 

(T + 1)/(T^2 - 2*T + 1)

sage: [f.derivative(T,j)(0)/factorial(j) for j in [0..10]]

[1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21]

# you can instead create a symbolic function:

sage: var('T')

T

sage: ff(T)=SR('(1+T)/(1-T)^2')

sage: ff.taylor(T,0,10)

T |--> 21*T^10 + 19*T^9 + 17*T^8 + 15*T^7 + 13*T^6 + 11*T^5 + 9*T^4 + 7*T^3 + 5*T^2 + 3*T + 1

[saad khalid]

Thanks! That makes much more sense and seems simpler, I'm glad I know that now!

Side note:  Never use += on a string in a loop.  Here it's even worse,
since you are concatenating strings on the right hand side as well.

Is it because it's slow?

Also, do you think I should have the to_sage function automatically convert it to a symbolic function or fracfield, or should I just leave it as outputting the string version?

[Vincent Delecroix]

Could this kind of discussion can be done on sage-support? As it name
suggests sage-devel is dedicated to development discussions. If we want
to keep developers reading it, it is better to not flood it with users
usage questions.

[saad khalid]

Sorry about that, I thought sage-devel was for talking/asking questions related to doing development on Sage. Since I was editting the source, I assumed it would be okay to post it here, but if it's not, I'll definitely remember that in the future.
I created a trac for review:
https://trac.sagemath.org/ticket/20936#comment:1

[Dima Pasechnik]

IMHO this is borderline between devel and support - it is about extending Sage interface to Macaulay2 with a new type.

[leif]

saad khalid wrote:
> Sorry about that, I thought sage-devel was for talking/asking questions
> related to doing development on Sage. Since I was editting the source, I
> assumed it would be okay to post it here, but if it's not, I'll
> definitely remember that in the future.

Well it is, but neither from the subject nor from your first post was
clear what you intended, or more precisely, what your motivation was.

I also first thought you had some specific problem to solve (unrelated
to Sage, but /using/ Sage and/or M2), and asked for help because Sage or
its interface to M2 simply didn't offer what you needed. But I
hesitated to kick you to sage-support... ;-)

Only after reading some further posts, I got the impression that you're
actually intending to /contribute to Sage/, and just gave an example of
application.

> I created a trac for review:
> https://trac.sagemath.org/ticket/20936#comment:1

Yep, that's where such a discussion should end up (or lead to). And
trac is by far better suited for reviewing code than a mailing list is.

Occasionally, people announce new tickets here, return for very specific
questions, for design decisions, or to attract reviewers, but the main
discussion and development happens on trac.


-leif