Feeds:
Posts
Comments

Posts Tagged ‘Sympy’


This code is supposed to be (if some one does the work in the future) located in sage.tensor.differential_form_element.

The code presented below is a slight modification of Joris code for differential forms manipulation on SAGE.

Needed modules

from sage.symbolic.ring import SymbolicRing, SR
from sage.rings.ring_element import RingElement
from sage.algebras.algebra_element import AlgebraElement
from sage.rings.integer import Integer
from sage.combinat.permutation import Permutation

The advantage of using this is that, tensors defined here are an Algebra element, not just a python object as in the previous code.

The sage.combinat.permutation won’t be used (yet), but could be useful if tensor symmetries are defined.

TensorFormatter

class TensorFormatter:
    r"""
    This class contains all the functionality to print a tensor in a
    graphically pleasing way.  This class is called by the ``_latex_`` and
    ``_repr_`` methods of the Tensor class.
    """
    def __init__(self, space):
        r"""
        Construct a tensor formatter.  See
        ``TensorFormatter`` for more information.
        """
        self._space = space

    def repr(self, comp, fun):
        r"""
        String representation of a primitive tensor, i.e. a function
        times a tensor product of d's of the coordinate functions.

        INPUT:

        - ``comp`` -- a subscript of a differential form.

        - ``fun`` -- the component function of this form.

        EXAMPLES::

            sage: from sage.tensor.tensor_element import TensorFormatter
            sage: x, y, z = var('x, y, z')
            sage: U = CoordinatePatch((x, y, z))
            sage: D = TensorFormatter(U)
            sage: D.repr((0, 1), z^3)
            'z^3*dx@dy'

        """

        str = "@".join( \
            [('d%s' % self._space.coordinate(c).__repr__()) for c in comp])

        if fun == 1 and len(comp) > 0:
            # We have a non-trivial form whose component function is 1,
            # so we just return the formatted form part and ignore the 1.
            return str
        else:
            funstr = fun._repr_()

            if not self._is_atomic(funstr):
                funstr = '(' + funstr + ')'

            if len(str) > 0:
                return funstr + "*" + str
            else:
                return funstr

    def latex(self, comp, fun):
        r"""
        Latex representation of a primitive differential form, i.e. a function
        times a tensor product of d's of the coordinate functions.

        INPUT:

        - ``comp`` -- a subscript of a differential form.

        - ``fun`` -- the component function of this form.

        EXAMPLES::

            sage: from sage.tensor.tensor_element import TensorFormatter
            sage: x, y, z = var('x, y, z')
            sage: U = CoordinatePatch((x, y, z))
            sage: D = TensorFormatter(U)
            sage: D.latex((0, 1), z^3)
            'z^{3} d x \otimes d y'

        """

        from sage.misc.latex import latex

        str = " \otimes ".join( \
            [('d %s' % latex(self._space.coordinate(c))) for c in comp])

        if fun == 1 and len(comp) > 0:
            return str
        else:
            funstr = latex(fun)

            if not self._is_atomic(funstr):
                funstr = '(' + funstr + ')'

            return funstr + " " + str

    def _is_atomic(self, str):
        r"""
        Helper function to check whether a given string expression
        is atomic.

        EXAMPLES::

            sage: x, y, z = var('x, y, z')
            sage: U = CoordinatePatch((x, y, z))
            sage: from sage.tensor.tensor_element import TensorFormatter
            sage: D = TensorFormatter(U)
            sage: D._is_atomic('a + b')
            False
            sage: D._is_atomic('(a + b)')
            True
        """
        level = 0
        for n, c in enumerate(str):
            if c == '(':
                level += 1
            elif c == ')':
                level -= 1

            if c == '+' or c == '-':
                if level == 0 and n > 0:
                    return False
        return True

The only I’ve changed here is “DifferentialForm” by “Tensor” and “\wedge” by “\otimes”

The above code allows to write the tensor product in a basis, dx^1\otimes\cdots\otimes dx^n. The chosen symbol for denoting the tensor product was @.

Tensor Class

This code is incomplete due to:

  • I’ve not defined the TensorsAlgebra, which should be done in parallel.
  • There are a lot of attributes not presented in this class.
  • class Tensor(AlgebraElement):
        r"""
        Tensor class.
        """
    
        def __init__(self, parent, degree, fun = None):
            r"""
            Construct a tensor.
    
            INPUT:
    
            - ``parent`` -- Parent algebra of tensors.
    
            - ``degree`` -- Degree of the tensor.
    
            - ``fun`` (default: None) -- Initialize this differential form with the given function.  If the degree is not zero, this argument is silently ignored.
    
            EXAMPLES::
    
                sage: x, y, z = var('x, y, z')
                sage: F = Tensors(); F
                Algebra of tensors in the variables x, y, z
                sage: f = Tensor(F, 0, sin(z)); f
                sin(z)
    
            """
    
            from sage.tensor.tensorss import Tensors
            if not isinstance(parent, Tensors):
                raise TypeError, "Parent not an algebra of tensors."
    
            RingElement.__init__(self, parent)
    
            self._degree = degree
            self._components = {}
    
            if degree == 0 and fun is not None:
                self.__setitem__([], fun)
    
        def __getitem__(self, subscript):
            r"""
            Return a given component of the tensor.
    
            INPUT:
    
            - ``subscript``: subscript of the component.  Must be an integer
            or a list of integers.
    
            EXAMPLES::
    
                sage: x, y, z = var('x, y, z')
                sage: F = Tensors(); F
                Algebra of tensors in the variables x, y, z
                sage: f = Tensor(F, 0, sin(x*y)); f
                sin(x*y)
                sage: f[()]
                sin(x*y)
            """
    
            if isinstance(subscript, (Integer, int)):
                subscript = (subscript, )
            else:
                subscript = tuple(subscript)
    
            dim = self.parent().base_space().dim()
            if any([s >= dim for s in subscript]):
                raise ValueError, "Index out of bounds."
    
            if len(subscript) != self._degree:
                raise TypeError, "%s is not a subscript of degree %s" %\
                    (subscript, self._degree)
    
            """sign, subscript = sort_subscript(subscript)"""
    
            if subscript in self._components:
                return sign*self._components[subscript]
            else:
                return 0
    
        def __setitem__(self, subscript, fun):
            r"""
            Modify a given component of the tensor.
    
            INPUT:
    
            - ``subscript``: subscript of the component.  Must be an integer or a list of integers.
    
            EXAMPLES::
    
                sage: F = Tensors(); F
                Algebra of tensors in the variables x, y, z
                sage: f = Tensor(F, 2)
                sage: f[1, 2] = x; f
                x*dy@dz
            """
    
            if isinstance(subscript, (Integer, int)):
                subscript = (subscript, )
            else:
                subscript = tuple(subscript)
    
            dim = self.parent().base_space().dim()
            if any([s >= dim for s in subscript]):
                raise ValueError, "Index out of bounds."
    
            if len(subscript) != self._degree:
                raise TypeError, "%s is not a subscript of degree %s" %\
                    (subscript, self._degree)
    
            """sign, subscript = sort_subscript(subscript)"""
            self._components[subscript] = SR(fun)

    Ok, so again I’ve changed “DifferentialForm(s)” by “Tensor(s)”, drop the permutation of indices (’cause tensors do not need to be neither symmetric nor anti-symmetric.

    I’ll keep working with this code… It’s all by now. Oh! I’ll post next week the rest of the code based in Sergey’s GRPy.

    Enjoy.

    Dox

    Read Full Post »


    I’ve just updated the SAGE worksheet which uses the definitions described in the previous posts.

  • There are some explanations in text format
  • The code has been hidden… because is long.
  • Moreover… I’ve discover something really amazing! Joris Vankerschaver‘s code of the differential form package in SAGE. Thus, I could use some ideas from Joris’ code to improve GRmodule. Nice, Isn’t it?

    Let’s hope I could so something nice this weekend!

    Don’t forget check the worksheet, and post some comment for feedback! 😉

    Enjoy!

    Dox

    Read Full Post »


    Hi everyone!

    This time the Christoffel connection will be defined.

    The code

    As usual, here is the code:

    class Christoffel(Tensor):
        '''The class to represent Christoffel Symbols of the second kind. Please
            note that while it inherits from Tensor, Christoffel symbols are
            NOT tensors'''
    
        def __init__(self,metr,symbol='C',rank=(1,2),sh=(1,-1,-1)):
    
            # The metric
            self.g_down = metr
    
            # Since we have a metric we do indeed have a coordinate system
            self.rep  = self.g_down.rep
    
            self.g_up = metr.inverse
    
            # Please note that this call will trigger a call to allocate in
            # the Tensor class, but the allocate will actually be the allocate
            # defined below
            super(Christoffel,self).__init__(symbol,rank,sh,coords=metr.coord)
    
        def allocate(self,rank):
            Tensor.allocate(self,rank)
            # Now that we have allocated things, time to actually compute things
            for i in xrange(self.dim):
                for k in xrange(self.dim):
                    for l in xrange(self.dim):
                        sum = 0
                        for m in xrange(self.dim):
                            term1 = diff(self.g_down[m,k],self.g_down.coords[l])
                            term2 = diff(self.g_down[m,l],self.g_down.coords[k])
                            term3 = diff(self.g_down[k,l],self.g_down.coords[m])
    
                            tr = self.g_up[i,m] * (term1+term2-term3)
    
                            sum += tr
                        res = sum/2
                        self.components[i,k,l] = res
            self.getNonZero()

    This code is almost a copy of Sergey’s one, except for the use of xrange instead of np.arange, and the fact that I’ve dropped the minus signs denoting the shape of the tensors.

    Sage implementation

    This time I won’t present a Python file, but a SAGE file, GRmodule.

    Enjoy!

    Dox

    Read Full Post »


    Hello everyone. I’ll do it quick, ’cause I’m really tired, and by the way, today… Feb. 8th, I’m turning 30 years old 🙂

    First of all, I create a page in my google site for the GRmodule files. Note that the name is GRmodule and not GR-module, because the dash is not admissible in a module name (beginner mistake!)

    In the last post…

    The Metric class was defined, and the first steps of the implementation were walked.

    Inverting the metric

    Once the components of the metric are given, one can call the invert method for assigning the inverse of the metric. Don’t know why, but I’d like to keep track of what’s going on! so, I change a little bit the code from Sergei, and my invert method returns the inverse metric tensor,

      def invert(self):
          '''Find the inverse of the metric and store the result in a
          Metric object self.inverse'''
    
          '''Create a unit matrix of dimension dim'''
          temp = sympy.eye(self.dim)
    
          '''Assign the values of the metric to temp'''
          for key in self.components.keys():
              id = tuple(np.abs(key))
              temp[id] = self.components[key]
    
          '''invert the matrix with inv() from sympy'''
          inv = temp.inv()
          '''convert the matrix in a dictionary'''
          inverse = self._dictkeycopy(self.components)
          for i in range(self.dim):
              for j in range(self.dim):
                  inverse[i,j] = inv[i,j]
          self.inverse = Metric(self.coord,rank=(2,0),sh=(1,1),symbol='g_inv')
          self.inverse.components = inverse
          return self.inverse

    First a temporal matrix is created, and just for assuring it’s invertible, one creates a unit matrix of dimension dim. That’s what the sympy.eye does!

    Secondly, the values of the metric are assigned to temp.

    The temp matrix is inverted. The sympy command to inver a matrix is inv.

    Next, the inverse matrix is converted in a dictionary (just like in previous cases).

    Finally, the characteristic of tensor is given to the inverse matrix object… and it’s returned!, i.e., it can be assigned.

    In the implementation file the components of the metric are assigned, the metric is invert… and the result of the inversion is printed.

    How to run it?!

  • Download the latest files and store them in the same folder, e.g., GR.
  • Open a terminal and move to the GR folder cd path/to/GR
  • In the terminal type $ python Proof-GR-module.py > Result.txt
  • Open the Result.txt file, i.e., $ emacs Result.txt
  • Enjoy life!

    DOX

    Happy B-day to me!!! 🙂

     

    Read Full Post »


    In the last two post I introduce the classes formalTensor (which has been not that useful) and Tensor, the latter include a series of attributes.

    Now is time to get started with GR 😛

    Metric Class

    I left this class just as Sergei defined, it looked weird to me at the beginning, but again I decipher it 😉 with the help of Reinteract, b.t.w., those who want to play around with Reinteract, it’s in Ubuntu repositories, so just type in your console,

    $ sudo apt-get install reinteract

    Ok. Continuing… here is the code,

    class Metric(Tensor):
      '''Represents a metric. Note that coordinates now MUST be provided'''
    
      def __init__(self,coord,rank=(0,2),sh=(-1,-1),symbol='g'):
          self.coord=coord
          super(Metric,self).__init__(symbol,rank,sh,coords=coord)
    
      def invert(self):
          '''Find the inverse of the metric and store the result in a
          Metric object self.inverse'''
    
          '''Create a unit matrix of dimension dim'''
          temp = sympy.eye(self.dim)
    
          '''Assign the values of the metric to temp'''
          for key in self.components.keys():
              id = tuple(np.abs(key))
              temp[id] = self.components[key]
    
          '''invert the matrix with inv() from sympy'''
          inv = temp.inv()
          '''convert the matrix in a dictionary'''
          inverse = self._dictkeycopy(self.components)
          for i in range(self.dim):
              for j in range(self.dim):
                  inverse[i,j] = inv[i,j]
          self.inverse = Metric(self.coord,rank=(2,0),sh=(1,1),symbol='g_inv')
          self.inverse.components = inverse

    Inherit…

    The first I notice, as a non-expert programmer (or non-programmer at all 😛 ) was the word Tensor inside the brackets of the Metric class. This means that Metric is a Tensor… just as a Tensor is an Object. Thus, all methods defined in the Tensor class are applicable to the Metric class.

    Coordinate system needed

    The __init__ method define the rank, shape and symbol of the metric, which are always the same. However, the set of coordinates must be given.

    Example of implementation

    Coordinated should be entered like a tuple, and the name of coordinates MUST be “declared” as sympy.Symbol

     t = Symbol('t')
    r = Symbol('r')
    th = Symbol('theta')
    ph = Symbol('phi')
    g = Metric((t,r,th,ph))

    Given the this data, the Metric class calls the Tensor class and create a covariant rank 2 tensor. And we should give the components. For example, the easiest Schwarzschild metric (the metric is easier than the name!!!),

    g[0,0] = -(1 - 2*M/r)
    g[1,1] = 1/(1 - 2*M/r)
    g[2,2] = r**2
    g[3,3] = r**2*(sin(th))**2

    WAIT!!!! Don’t forget to declare the mass parameter, so, again…

    M = Symbol('M')
    g[0,0] = -(1 - 2*M/r)
    g[1,1] = 1/(1 - 2*M/r)
    g[2,2] = r**2
    g[3,3] = r**2*(sin(th))**2

    … Sorry guys, it’s again too late! and is already Monday! :-/

    Ok, tomorrow we’ll continue with the inverse metric 😉

    Get GRmodule.py
    and the proof-file.py.

    Cya tomorrow!

    Dox

    Read Full Post »


    A few days ago I post in sage-devel google group my desire of contribute in the development of a general relativity (perhaps beyond) module.

    A search in Google leads me to a class, named GRPy, created by Sergei Ossokine. My experience with it was not as pleasant as I’d like, but it’s a beginning for starting my own Module. 🙂 Yeah!!!

    However, since I’m not a programmer, there are lot of things to learn in the way.

    The requirements

    We are suppose to do symbolic calculations, so a requirement is the python package for doing so, say sympy

    import sympy

    For managing arrays, one of the most useful packages is numpy, so we could call it,

    import numpy as np

    this line does not import the module… but if we want to use a command of the module, i.e. array, we should call it as, np.array.

    I just learnt that there exist a module for iterate, and in particular allows to define Cartesian products. It is called itertools. Let’s import it,

    import itertools

    Finally, from the copy module, let call the deepcopy function,

    from copy import deepcopy

    which copy and object and all structures inside it!

    Define a Tensor Class

    I’ll follow to Sergei at the beginning… so let’s start as he did.

    A tensor is an object, so is should inherit the characteristics of one. Moreover, a tensor has an identifier symbol, and a rank

    class formalTensor(object):
    ... def __init__(self,symbol,rank):
    ...    self.symbol = sympy.Symbol(symbol)
    ...    self.rank = rank

    I still don’t now how useful is to define this formalTensor class. But the next one is the really important,

    class Tensor(object):
      '''A class to represent a tensor in a particular basis'''
      def __init__(self,symbol,rank,shape,sym=None,coords=None,**args):
        self.symbol = sympy.Symbol(symbol) # The symbol to represent this Tensor
        self.coords = coords # The coordinate system we are using for the representation
        self.shape = shape # The shape of the tensor, for example (-1,1) for k_{a}^{b
        self.rank = rank # Rank
        self.contr = rank[0]
        self.cov = rank[1]
    
        # We need to know the dimensionality of our spacetime. For the moment
        # we will deduce it from the coordinates provide, or if none are given, then
        # the assumption will be made that it's stored in the optional arguments
        if coords is not None:
            self.dim = len(self.coords)
        else:
            self.dim = args['dim']
    
        if coords is not None:
          self.allocate(rank)
          self.rep = True
        self.symbolic = formalTensor(symbol,rank)
    
      def __setitem__(self,idx,val):
        self.components[idx]=val
    
      def __getitem__(self,idx):
        return self.components[idx]
    
      def allocate(self,rank):
        '''Allocate the dictionary(hash table) necessary to store the components
          Note that covariant indices are negative! (except for 0 of course)'''
        n = rank[0] + rank[1]
        indc = list(itertools.product(range(self.dim),repeat=n))
        mastr=[]
        for i in range(len(indc)):
            temp=[]
            for k in range(len(indc[i])):
    
                if self.shape[k] == -1:
    
                    temp.append(-indc[i][k])
                else:
                    temp.append(indc[i][k])
            mastr.append(tuple(temp))
        self.components = dict(zip(mastr,[0 for i in range(len(indc))]))
    
      def _dictkeycopy(self, hay):
        keys = hay.keys()
        return dict(zip(keys,[0]*len(keys)))
    
      def getNonZero(self):
        '''Returns only non-zero components of the tensor, if the coordinate
        system is provided'''
        if self.rep:
          nonzerok =[]
          nonzerov = []
          for key in self.components.keys():
            if self.components[key] !=0:
              nonzerok.append(key)
              nonzerov.append(self.components[key])
          d = dict(zip(nonzerok,nonzerov))
          keys = d.keys()
          keys.sort()
          self.nonzero = [(key,d[key]) for key in keys]
          return self.nonzero
        else:
          print "Attempted to get components that have not been initialized!"
    
      def __str__(self):
        '''Print a "nice" human - readable representation of the tensor'''
        self.getNonZero()
    
        # We will print only non-zero components unless all the components are zero
        ttl=""
        if self.nonzero:
    
          for i in range(len(self.nonzero)):
              ttl += str(self.nonzero[i][0]) + " "+str(self.nonzero[i][1])+"\n"
        else:
          ttl ="All the components of the tensor are 0!"
        return ttl

    WOW… long way to go explain!

    Tensors have,

  • Symbol: for identification.
  • Rank: or number of indices… even if it is a connection
  • Shape: a list in the form (a,b) for denoting a tensor a times contravariant and b times covariant.
  • Symmetry: not implemented!!! 😦
  • Coordinate system: A list with the name of coordinates, example: (t, r, theta, phi).
  • Any other arguents: that’s what **arg stands for. Example dim=3 would set the dimension to 3.
  • The constructor

    Any class has a constructor named __init__, where the variables of the object are defined and initialized. Additionally, one should call the variables using the prefix self. For a nice introduction to classes see for example A primer on scientific programming with Python.

    Said so, the first few line can be followed, let’s jump to the first non-so-trivial thing.

    Special Methods

    Methods are the possible actions defined for objects in our class. Those methods whose names start with underscore are known as special methods.

    Above, two of these are used.

  • __setitem__: which set a given component of the tensor,
  • __getitem__: which return the value of a given component.
  • Ok, it’s midnight in UK and I’m really tired.

    Tomorrow I’ll explain the allocate method.

    I’m uploading the python file to my page (as usual… at the end of the page, in the attachment section).

    So far… GR-module.py (version 1)

    Enjoy!

    DOX.

    Read Full Post »