Skip to content

GSoC 2012 Application Stefan Krastanov: Vector Analysis

Jump to bottom Edit New page
Sean Vig edited this page Jul 18, 2013 · 2 revisions

GSoC 2012 Application Stefan Krastanov: Vector Analysis

Personal Background and Skills

My name is Stefan Krastanov and I am currently a MSc Physics student at ENS Lyon, France. I can be found on the mailing list and on the google code page ([email protected]), on GitHub (Krastanov) and occasionally on IRC.

I have a decent knowledge of programing and I am comfortable in a number of high-level languages. Lisp and Python fascinate me most. I believe that I have an in-depth knowledge of the pythonic idioms. I am also able to work idiomatically in C, Fortran and C++.

In the last few years the only serious coding that I have done was in python. With the exception of my past contributions to sympy, it was generally related to laboratory work (data acquisition and processing, preparation of simulations, testing, etc).

I work on a linux workstation, mainly in Vim, and I use Ipython for my python work. For version control of every personal project I use git. I have used this setup for a few years and I am proficient with it.

Past contributions to SymPy

The Project

My proposal is to work on a vector analysis framework for SymPy. It evolved from my wish to work on the tensor module (a project mentioned in the ideas page).

The discussion that led to this choice is on the mailing list:

Most of the project will be based on this publication http://groups.csail.mit.edu/mac/users/wisdom/AIM-2005-003.pdf. I am contacting the authors to request the lisp source code.

A detailed list of all the features that I want to implement can be found below, however, in a few words the project is about being able to represent scalar, vector and form fields in coordinate independent way in SymPy. The only way to work on vector analysis currently in sympy is to explicitly work with matrices and it is still quite clumsy.

Moreover, while not absolutely necessary for the creation of an abstract tensor module, this project will provide a unifying framework for future work in this direction.

I have the necessary mathematical background for the subject and I would have a great use for a tool that automates most of the task related to vector analysis. Spending the required 40h per week on the project will not be a problem.

Schedule (each bullet point is roughly a pull request)

  • Week 1 (13th May - 19th May)
    • manifolds and coordinate systems, scalar functions
  • Week 2 (20th May - 26nd May)
    • vector functions, partial derivatives, directional derivatives, nabla
  • Week 3-5 (27rd May - 16th June)
    • integral curves
    • one forms, higher order forms, tensor and wedge products
    • exterior derivatives
  • Week 6-10 (17th June - 21st July) (13 July is the midterm evaluation)
    • over-a-map operations
    • Lie derivatives
    • interior products
    • Midterm Evaluation
    • covariant derivatives, geodesic motion, connections
    • metrics
    • Hodge star
  • Week 11 (22nd July - 28th July)
    • writing trivial examples (tests already written during development) and notebooks with more involved examples showing real usage
    • add the already written docstrings to the sphinx documentation
  • Week 12-13 (29th July - 11th august) (it is not listed in Features bellow as it is mainly research)
    • Not all of these will get to a useful state:
      • Much of the code will depend on solve() in order to calculate inverse transformations. Where solve fails the user will be able to input the transformations explicitly, but it would be better if he does not need to. Hence, I will try to document all the cases that fail in order to have a guide to what future improvements of solve should be done. I may start to work on some of those improvements.
      • Create the basis for the tensor algebra framework (or just study it if it is already done by another GSoC project) and couple it with the vector analysis framework. In my project tensor products of forms or vectors would have been already implemented, but abstract index operation (i.e. like in xAct) is another thing. This is completely different formalism, however it is important to be able to switch between it and what would be implemented in my project.
      • Research a nice way to add operator fields (it would be already possible to build it from forms, however this is not the easiest way for beginners). It would be especially interesting to have assumption on those operators (unitary, orthogonal, symmetric, etc).
      • Check how plausible it is to use this module for infinite dimensional vector spaces.
      • Simplification of nested nablas?
      • The project is related to implementing abstract linear algebra (a project in which other students have interest). Check how well the two work together.

Features

An example with scalar fields over a trivial manifold:

# There is some implicit manifold behind all this.
# In the final code this will be explicit.

# define coordinate systems
R2 = RectangularCoordinateSystem(dim=2)
P2 = PolarCoordinateSystem(dim=2)

# from abstract point to coordinate representation and back
R2_chi = R2.to_coords # take a Point instance and return a column vector
R2_chi_inverse = R2.to_point # take a column vector and return a Point
P2_chi = P2.to_coords
P2_chi_inverse = P2.to_point

print P2_chi(R2_chi_inverse(Matrix([x, y])))
> [ sqrt(x**2 + y**2) ]
> [ atan(y,x)         ]

# define a point
point1 = R2_chi_inverse(Matrix[x, y])
point2 = P2_chi_inverse(Matrix[rho, theta])


# define a scalar function over the manifold
f_in_polar = a_SYMPY_function_that_takes_two_arguments
f = ManifoldFunction(f_in_polar, P2_chi)

print f(point1)
> x**2 + y**2
print f(point2)
> rho

# also some functions are already defined as methods of the coordinate system
print R2.x(point2)
> rho*sin(theta)

List of Features (i.e. rules to implement) and tentative class structure:

  • defining a manifold with arbitrary (in form and number) coordinates systems
    • Manifold(Basic) (it is Basic so it can be used in expressions like cartesian products of manifolds)
    • CoordinateSystem(Basic)
      • a number of common ones already implemented (eg Cartesian2D or Polar)
      • having PointToCoordinates and CoordinatesToPoint methods that are sympy functions
    • Point(Symbol) having a Manifold and a CoordinateSystem among its args
  • scalar fields
    • ScalarField a sympy function(al) that has a point2coordinate function and a coordinate2scalar function among its args
  • the rest will be implemented in classes akin to ScalarField or containers akin to Mul and Add or operators like Derivative
  • vector fields
    • basis fields
    • directional derivatives of scalar fields
      • including the special cases of div, rot, laplacian, ...
    • integral curves
      • differential equation for specified coordinates
      • Taylor expansion independent of coordinates
      • evolution operator
    • commutators (of vectors/directional derivatives)
  • one-form fields
    • basis fields and dual fields
    • differentials of scalar functions
  • higher rank forms
    • zero order form
    • tensor product of higher rank forms
    • wedge product
      • one-forms
      • higher rank forms (not very trivial if they are not products of one forms)
  • exterior derivative
    • closed, exact
  • over-a-map operations
    • pullback of a function, vector or form field
    • pushforward of a vector field
    • differential of a map
  • Lie derivatives
  • interior products
  • covariant derivatives, Christoffel's symbol
  • parallel transport, geodesic motion, curvature
  • connections, torsion
  • metrics
  • Hodge star operator
  • Other objects not directly related to vector calculus like Vector, Operator, VectorSpace, BilinearForm, Tensor, etc (those are not fields, no calculus, just algebra). Those may be implemented by other gsoc projects (there have been at least two propositions on the mailing list). If not, only skeleton implementation of those will be done, as discussed on the mailing list. This will permit work on base-independent linear algebra (not the same as coordinate-independent vector calculus). For example it will permit solving http://scicomp.stackexchange.com/questions/74/symbolic-software-packages-for-matrix-expressions. I add these objects in the notes here so it is clearer what will not be part of the project (for example Gram–Schmidt_process on abstract vectors will not be in the project as it involves linear algebra and not vector calculus). Nonetheless, the basic objects and operations on them will be added to ensure that there will not be future conflicts.

More examples

The examples are not written in OOP style in order to have the sequence of operations clearly shown.

Transformations of coordinates:

# define coordinate systems
R2 = RectangularCoordinateSystem(dim=2)
P2 = PolarCoordinateSystem(dim=2)

# or calculate it directly
jacobian_matrix(R2, P2) # returns a function that takes two symbols (eg. rho and theta) and returns the matrix

# get an explicit expression of x and y in terms of rho and theta
x_and_y_in_terms_of_rho_and_theta = R2.to_coords(P2.to_point(Matrix([rho, theta])))

Get the differential equation in rectangular coordinates for the integral curve of a vector field defined in polar coordinates:

R2 = RectangularCoordinateSystem(dim=2)
P2 = PolarCoordinateSystem(dim=2)
A = VectorField(name='A', coord_system=P2, components=[sympy_function_of_rho_and_theta_returning_component_1,
                                                       sympy_function_of_rho_and_theta_returning_component_2
                                                       ])
eq = integral_curve_equation(field=A, coord=R2, vars=[x,y])

Common manifolds will be already implemented in the code, but the user will be able to define new ones. Here is a very rough example of how this will be done:

my_manifold = Manifold('R3')
my_manifold.add_coord_system('cart')
my_manifold.add_coord_system('spherical')

my_manifold.status()
> cart and spherical are unrepresentable

my_manifold.connect('cart', 'spherical', [r(x,y,z), theta(x,y,z), phi(x,y,z)])
my_manifold.status()
> `solve()` was unable to create spherical->cart transformation, specify one explicitly if you need it

It should also be possible to define cartesian products of manifolds (eg a torus):

my_torus = circle_manifold*circle_manifold
my_torus.coord_systems
> [ list of all possible products of coordinate charts ]

Other examples without code snippets:

  • calculate the exterior derivative of div(A)dVolume and show that it is AndSurface
  • given a path and a vector, calculate the vector obtained by parallel transport of the original vector over the path
  • given a velocity vector and a starting point get the geodesic path
  • from a vector basis get a dual form basis and vice-versa
  • given a metric construct metric-dual vector/form basis to a given form/vector basis
  • calculate hodge duals of forms

Implementation Style:

Most of the objects will derive from Symbol or Function. The manifolds and spaces will derive from Basic or even object. Something should be done for vector valued functions. The derivatives, integrals, wedge products, etc. will have unevaluated representation, subclassing the appropriate operators in sympy (usually Expr).

Contraction of a form and a vector will probably subclass Expr (otherwise Mul will be both scalar product and contraction).

Mul and Add will be used (no subclassing of Expr for them) in some way akin to evaluate=False. .flatten will be used for the scalar part of the expressions because it is the closest thing to a canonicalizer that we have, but the rest will be canonicalized by some global function (vector_canonicalizer for example). At the end of the project this will be moved inside __mul__ and __add__ (ugly solution but to me it seems that it is the best before having some kind of AST crawlers, this was discussed on the mailing list).

The rest is just implementing straightforward mathematical rules. Almost exclusively the "naive" algorithms will be used.

Add a custom footer
Add a custom sidebar
Clone this wiki locally