The Vectors and Electromagnetism GSoC projects

[Stefan Krastanov]

Today we had the first discussion with Prasoon and Sachin about their projects.

We did not progress much but at least we outlined the two general approaches that we can use for these modules (specifically for creating vector fields). I will give them somewhat arbitrary names here:

- the `mechanics` way - having a Vector class that keeps all the information about the field and it is not part of expression trees in the way Basic and Expr are. For instance Vector(something along cartesian.x)+Vector(something along spherical.r) will result in Vector(complex internal structure).

- the `diffgeom` way - having base/unit vectors and building all the rest in terms of their linear combinations (all this expressed as sympy expressions).

Prasoon and Sachin did not have the time to look at the example problem that was given in the previous email yet (no harm done there, there is still some time before the official starting date). Probably this will be the subject of our next discussion.

The next discussion was scheduled for tomorrow. After that I suggest that we keep most of the discussions to the mailing list and the gihub wiki and meet on irc / realtime wikis / google docs / etc  once a week.

Stefan

[Aaron Meurer]

The discussion was at http://piratepad.net/KBviCWUlA3.

I'm curious what you think of this kind of discussion, as opposed to
IRC. Google docs is also an option (it has a chat). I think the
downside is that unlike IRC, which is logged at
http://colabti.org/irclogger/irclogger_logs/sympy, it's a little
harder to search through these discussions afterwords.

Aaron Meurer

> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy?hl=en-US.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Stefan Krastanov]

Here is a quick summary from today:

- probably scalar fields will be represented simply by SymPy expressions where some of the symbols will have special meaning (the coordinates)

- probably vectors will be represented like in mechanics (one object, not necessarily a sympy expression)

- using reference systems translated and rotated with respect to each other is rather unclear at the moment: before continuing with the irc meetings I would suggest that the students provide a _nice_ wiki presenting the answers to the questions in the previous thread and also the following questions:

Bellow is the question in "mathematical" terms. Transform it in whatever way you find appropriate to fit your suggested APIs:

- in 3D

- I have three points A, B and C.
- I use each of them as the zero of three different coordinate systems

- The A and B systems are both Cartesian but rotated by theta_AB around axis_AB

- The C system is spherical (r, phi, theta). The theta=0 axis is rotated wrt the z axis of A by the Euler angles alpha, beta, gamma

- I define a scalar field in A, another scalar field in B and a vector field in C

- I want the sum of the scalar fields

- I want the gradient of that sum

- I want the convective derivative of the vector field from C wrt the gradient from the question above

- I want to express the entities from the above 4 question in each of the three coordinate systems.

- For all this please explicitly choose some fields for the examples and calculate the expected results by hand (and add them to the example session as mock results).

I think that this will really stress test the suggested API. The only thing missing is the time dependence needed in mechanics. I strongly suggest that we first finish the considerations above before continuing.

@Prasoon and Sachin, when will you be able to provide a detailed wiki page with an example session for what is asked here? There is really no need to hurry (officially GSoC has not started yet) so please take your time (a week?).

Stefan

[Stefan Krastanov]

And to be mean I will also write this down using the diffgeom module :)

Seriously though, while the diffgeom module uses a completely different approach that is not compatible with the needs of the `mechanics` module (on which we are focusing), having such a comparison will be useful to show bad choices both in the diffgeom module and in the present project.

[Jason Moore]

Thanks for sharing this. Just one comment. I wouldn't let the mechanics way prevent you all from doing this the "right" way. With the right way being whatever a more general approach is. The mechanics module was designed from a very narrow point of view in terms of the mathematics. Without thinking about this very hard, it seems that whatever we do in mechanics has got to be some small subset of the generalities of vector calculus but we make sure to really carefully deal with the details that most affect our systems (rotations, vectors, reference frames, time differentiation).

I'm looking forward to seeing some of the problems done out by hand that are outside the scope of the way we think about this from the mechanics perspective. Maybe then I can comment better on the how mechanics would work with a generalization. The stress test you wrote out is a nice start.

Jason
moorepants.info
+01 530-601-9791

[Sachin Joglekar]

@Stefan : I would recommend having a separate class for ScalarFields. Even I wasn't sure of the need for this till yesterday, when we came across the problem of how the user would define a scalar field in any coordinate system he wants(which is not the global frame). In such cases, I propose something like the following-

rho = ScalarField(6*x**2*y, c_rect1)

rho.express(c_sph)

If we find a better way to this this I am all for it, but for now the above seems elegant and convenient.

In any case, I tried my hand at expressing a few of the steps put forward by Stefan in a mock SymPy session. I would request all of you to have a look at it and express your views (please elaborate on the reasons for any editing if done in the current code). It is a WIP obviously.

The wiki link - https://github.com/sympy/sympy/wiki/Vectors-EM-framework

[Gilbert Gede]

Sachin,

I like where you are going with this. If I'm interpreting it correctly, each CoordinateSystem has to be attached to a ReferenceFrame, and is fixed (although possibly rotated and/or translated upon coordinate system definition) with respect to that ReferenceFrame? Prasoon, Stefan, others - what are your thoughts on this?

I think there might be some better/easier ways to access basis vectors and some other issues, but that discussion can come after a consensus is reached on the CoordinateSystem/ReferenceFrame distinctions.

-Gilbert

[Sachin Joglekar]

@Gilbert, we could also let a CoordSystem have a motion, with something like 

system_a.set_vel(translational = ..., angular = ....) <- this would be with respect to some system defined in that frame only.

Then coordinates of this system, when expressed in some other system, would be functions of coordinate variables AND the 'time' variable of the global reference frame.

[Stefan Krastanov]

@Sachin, it would be best to finish what we have here before we starting discussion of motion.

I will add a TODO to the wiki page with the issues that I see (in about an hour).

[Stefan Krastanov]

Updated with a todo list: https://github.com/sympy/sympy/wiki/Vectors-EM-framework

Try to make the wikipage easy to read. Think of this as if it was documentation for the module.

@Prasoon, if what Sachin has written is incompatible with what you have in mind, just start a section below what he has written.

[Sachin Joglekar]

Stefan, could you have a look at the page again?

[Prasoon Shukla]

Okay then. I have read everyone's comments and also the stress test. At the moment, I believe that using the things we have discussed in the last two days, we can have a mock session for this stress test without much difficulty.

I still think that some of the things I mentioned in the proposal aren't very clear. What I think I will do is that I'll add some explanatory notes about what's what on the wiki page first. Hopefully, that will clear out some of the things that are unclear at this point.

Once that part is done, I shall try to provide the examples that would further clarify any remaining doubts.

Also, I will make a separate implementation than what Sachin has made (in a separate section at the bottom of the page) just so that ideas don't end up getting mixed. Once we are clear on all things, then we can make a unified wiki page that details the API.

[Stefan Krastanov]

Gilbert, could you check what Sachin has written on the wiki about `ScalarField`, namely the part that it is necessary when using different origins for the time variable.

IMO this is unnecessary complication, but maybe I am missing something.

--

[Gilbert Gede]

I believe that the part about different origin for the time variable has been taken out? I agree though, it seemed to add an unnecessary level of complication.

But I left some comments and questions throughout the wiki.

[Prasoon Shukla]

@Everyone:

I have added a separate heading towards the end of the wiki page. I would have done so yesterday but my internet had given up on me yesterday.

Anyway, a few things I'd like to mention before you read it.

First, I did not read the existing wiki page before adding content to it. This I primarily did so that I do not get influenced by the content already mentioned on the wiki page. I will read the existing page now and then we can work on merging the solutions that arise from the discussion.

Second, I believe that there were quite a few confusions related to the working of the CoordSys class and the RefFrame class. I have tried to discuss that issue a bit as well.

Then, I have tried providing a mock session as well. The point about Euler angles confused me so I have not incorporated that in the additions to the wiki page.

@Stefan : Could you explain in a bit more detail how exactly is the frame C rotated wrt to frame A. If I understand Euler angles correctly, then two frames first need to be aligned and then rotation if performed on one of the frames. But in this case, frame A and C aren't aligned. So, what exactly does saying 'wrt A' mean in the present context?

[Stefan Krastanov]

I have added more questions. I strongly dislike the separation between CoordSys and RefFrame: such an abstraction is a nice idea, but the current suggestion is very confusing and some parts of the definitions seem unnecessary.

At both of you: I still do not see why is there a need for ScalarField class (some parts of you original suggestions require it (especially Prasoon's), but I have explained on the page why I consider them a bad idea). The need for such a class seems to be caused by lack of clarity on how to specify what symbols correspond to coordinates in you CoordSys/RefFrame/etc classes.

--

[Gilbert Gede]

Sachin put a link to this wikipedia article in the wiki: http://en.wikipedia.org/wiki/Frame_of_reference#Different_aspects_of_.22frame_of_reference.22

The way I have read the page is that coordinate systems should not change relative to the reference frame they are defined in; the conclusion I developed was that frames should be defined relative to each other and coordinate systems should be defined (and fixed) in frames. What are others' thoughts on this?

Stefan,

If you have the same scalar field, written in a rectangular coordinate system as x * y * z and in a spherical coordinate system as r*sin(theta)*cos(phi) * r*sin(theta)*sin(phi) * r*cos(theta), and these two coordinate systems have the same orientation and origin, the scalar fields should be equivalent, right? Should they also evaluate as the same in SymPy?

-Gilbert

[Sachin Joglekar]

Stefan and Prasoon, could you comment on what you gather from that article?

What I got was this - 

"A reference frame is essentially about a state of motion. Every reference frame can have various coordinate systems attached to it. These systems may be oriented (rotated and translated) in certain ways with respect to each other. If you want to define two *somethings* that differ in their state of motion with respect to each other, then the *something* must be a frame of reference. In a frame of reference, if you want a different *way* of _measuring_ things, then this *way* is a coordinate system. Hence, every coordinate system MUST be attached to a reference frame, without which its existence dosen't make sense.

The only things that can be changed about frames and systems are their orientations and motion(only for frames) with respect to each other."

Stefan, I would like your comments on this.

About a ScalarField class, it is just to -

1. Provide a structure similar to Vectors (so users can do something like- for x in list, print field.express(system))

2. To ensure expression of the same field in different systems

3. To ensure field1 == field2 returns True if the fields are indeed equal (mathematically)

[Prasoon Shukla]

I have made some comments on the wiki. Let us get the API decided so that both I and Sachin can get started with some work. Thanks.

[Sachin Joglekar]

+1.

I added one or two comments, they are mostly just clarifications. But I guess we agree on the basics now.

I will post some hand-written scanned stuff about the math of Stefan's stress test on my blog in a day, mainly to clarify things in my own mind as well.

Prasoon, could you expand on the convective derivative part? Apart from that, I think my blog post + our wiki will clear it all up.

@Stefan, Prasoon, Gilbert - How do you guys think we should proceed with coding?

On Mon, Jun 10, 2013 at 3:04 PM, Prasoon Shukla <prasoon...@gmail.com> wrote:

I have made some comments on the wiki. Let us get the API decided so that both I and Sachin can get started with some work. Thanks.

--
You received this message because you are subscribed to a topic in the Google Groups "sympy" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sympy/t-Je0beTIIU/unsubscribe?hl=en-US.
To unsubscribe from this group and all its topics, send an email to sympy+un...@googlegroups.com.

[Stefan Krastanov]

@Gilbert

> If you have the same scalar field, written in a rectangular coordinate system as x * y * z and in a spherical coordinate system as r*sin(theta)*cos(phi) * r*sin(theta)*sin(phi) * r*cos(theta), and these two coordinate systems have the same orientation and origin, the scalar fields should be equivalent, right? Should they also evaluate as the same in SymPy?

They should not - `__eq__` checks for structural equality of the expression tree, not for mathematical equality. Actually, even trying to get these expression to evaluate to the same thing is futile - it is not possible in general (just getting trigonometry involved will probably be sufficient to make the expressions to complex.)

@Sachin

> 3. To ensure field1 == field2 returns True if the fields are indeed equal (mathematically)

This is one of the most problematic antipatterns in sympy. Everytime this is done it breaks most of the tree manipulation routines, from `subs` through `atoms` to `simplify`. Please look at the sympy core, tree traversal, etc so you know why we need to be careful about `__eq__` and `__hash__`.

Sorry for nor answering sooner, but I was traveling the last two days.

I actually do not think the api design is anywhere near completion. Many of the issues are not addressed:

- dependence of string flags in many constructors that seem redundant given the creation of objects meant to represent the same thing (for instance the creation of points)

- ScalarField class - All of the arguments for seem very vague. Two arguments against are as follows: Very seriously breaks tree traversal (which is the most fundamental part of sympy) and requires reimplementation of practially every operator (+ - / * etc). Just consider what happens when I want to take exp(some_scalar_field).

- Many of the issues surrounding ScalarField are relevant also for Vector.

--

You received this message because you are subscribed to the Google Groups "sympy" group.

To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.

[Sachin Joglekar]

@Stefan. I agree about the tree traversals part. However, I still feel there must be a way for the user to check if two scalar fields expressed in different ways are mathematically equivalent. Could you suggest a way? Maybe just another function?

Also, have a look at the __eq__ function of the current Vector class.

I would really like atleast _some_ way of checking mathematical equality.

I think we should let everyone opine on this equality thing, because it will, as you said, affect the entire framework, even Vector.

@Jason, Aaron, Luke, Sean.... thoughts?

[Stefan Krastanov]

> @Stefan. I agree about the tree traversals part. However, I still feel there must be a way for the user to check if two scalar fields expressed in different ways are mathematically equivalent. Could you suggest a way? Maybe just another function?

This is described in one of the first pages of the sympy
documentation: Checking for equality in general is mathematically
impossible. In practice most of the time one does `simplify(expr_a -
expr_b)` and sees whether he gets 0. In numerics one just plots the
difference.

> Also, have a look at the __eq__ function of the current Vector class.

The current Vector class is not really a part of SymPy. One of the
reasons that these two projects are useful is that they will make
mechanics and sympy more inter-operable so stuff like `simplify` can
be used seamlessly in mechanics.

> I would really like atleast _some_ way of checking mathematical equality.

As I mentioned above, there is already a way for that and the issues
surrounding it are described in the docs (well, the documentation can
be better, that is clear :)


Anyway, given all these issues I can not imagine an api similar to the
currently suggested one working well with sympy. Both ScalarField and
Vector will break most of sympy if they eagerly eat all expressions
that are added to them. Unrelated to that, ScalarField is ill defined
because all the confusion around reference frames and coordinate
systems and reuse of Symbols. Reference Frames and/or Coordinate
Systems are represented as mutable objects which is nonsensical
mathematically (one does not change how to references relate (one
might want to parametrize something wrt time, but this is not the
same))...

@Sachin and Prasoon, can you work together on the next iteration of
the design that is meant to fix these issues? When would you be able
to show it?

Multiple times I have suggested partial/fragmented designs similar to
`diffgeom`. While some parts of `diffgeom` have no place here, some of
the ideas (the representation of general fields based on sympy
expressions containing base fields) could be helpful in preserving all
the invariants important for sympy. Do you want me to provide a more
complete suggestion on which you can base yours, or do you prefer to
continue on your own for the moment?

[Gilbert Gede]

Stefan,

I think I disagree about Reference Frames being immutable. If you are going to create a series of Reference Frames when performing successive rotations, you are going to need to need to be able to modify the original frame to point to the new frames. There is a similar issue for setting the angular velocity of a Reference Frame after it is created.

Also, do you think that Vector should be a subclass of Expr?

-Gilbert

[Stefan Krastanov]

about immutability and reference frame: I guess I do not understand
the details of what you describe. Anyway, I will not be against
immutability if it does not mess up with `__eq__` and `__hash__` so it
depends on the implementation.

about Vector: if it is not a subclass of Expr it will not be useful
outside of mechanics. For instance
symplify(sin(asin(phi))*unit_vector) will just raise an error instead
of return phi*unit_vector.

[Sachin Joglekar]

Can we define our own simplify method for Vector that simplify(Vector) will call?

Also, if Vector will is a subclass of Expr, how can we avoid something like unit_vector * unit vector or Scalar+Vector?

You received this message because you are subscribed to a topic in the Google Groups "sympy" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sympy/t-Je0beTIIU/unsubscribe?hl=en-US.
To unsubscribe from this group and all its topics, send an email to sympy+un...@googlegroups.com.

[Gilbert Gede]

For an example of successive rotations, check out the second code block here: http://docs.sympy.org/0.7.2/modules/physics/mechanics/examples.html

Regarding Vector, should the user being to perform any action they could on another Expr object on a Vector object? If it inherits from Expr, would sin(unit_vector) not raise an error?

-Gilbert

[Chris Smith]

You can define and _eval_simplify method:

>>> def _eval_simplify(self,ratio=None, measure=None, fu=None):

...  return self+1

...

>>> Symbol._eval_simplify = _eval_simplify

>>> simplify(x)

x + 1

[Stefan Krastanov]

> You can define and _eval_simplify method:

That is cool. But many such methods will have to be defined and most
of the time it will be just copy pasting from Expr with maybe minor
changes. Fixed done to the core or changes of api will have to be
redone for Vector. And many sympy methods simply fail if their
argument is not Basic.

[Sachin Joglekar]

Maybe we could subclass from Basic, maybe even Expr.

We should just ensure *wrong* operations on Vectors return an error. (One of which Gilbert mentioned)

Stefan, could you put up your idea of the API for the stress test on the wiki? That would give us (Prasoon and me) a direction.

[Stefan Krastanov]

I will do it in the evening. And obviously my suggestions can be
flawed, but I will try to explain why I prefer it in places where it
differs from what you suggest.

[Prasoon Shukla]

Let us get some things clear here. Vector is going to be subclassed from Basic. Also, after seeing the discussions/problems above,  Vector + Vector will give an Add on which we can call something akin to a .doit() to result in a Vector. I do not see any problems with this approach.

Next, reference frames. First thing is the distinction between RefFrame class and CoordSys class. Stefan has pointed out that there is no need to implement the classes differently. Okay, then we will have just one, immutable, CoordSys class. The parameters required for intialization of kinematic parameters will have to be provided when declaring the coordinate system. New CoordSys object will result from methods like rotate, orient etc.

With these suggestions, I do not think that there should be any more structural issues. If there are any more of these issues, please mention them here now and please mention them as early as possible.

These changes are justified, I hope, because this module has to play nice with the rest of SymPy.

I just feel that we are getting stuck with and spending too much time on things that we shouldn't be taking this much time. Therefore, I request again to all parties involved : Please voice your views/problems with these structural changes suggested so that we can move ahead with things.

[Sachin Joglekar]

@Prasoon, Vector + Vector is fine. What about vector * vector or scalar + vector or any other such nonsense operation? If you are going to subclass from Basic, this can be done with handling of the operators as I have already pointed out. But how we will do this without operator overloading in case of Expr subclassing is unclear to me. +1 for subclassing from Basic. May involve redefining certain functions, but I think thats essential. 

In any case, I agree about coming to a consensus soon to avoid delay.

--

[Matthew Rocklin]

I'm only skimming over this thread but I noticed that the Basic/Expr question came up. 

MatrixExpressions used to inherit from Expr for exactly the reasons Stefan said (lots of things work out of the box).  This was great for startup but eventually became a major pain.  At some point I rewrote MatrixExpressions to inherit from Basic instead and made my own MatMul/MatAdd classes.  Development has been much smoother since then.  My recommendation is to inherit from Basic and make your own Add/Mul classes and your own _eval_simplify etc... methods

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.

[Stefan Krastanov]

This is also one of the ways to go forward. I personally dislike it because:

- A non trivial amount of functionality needs reimplementation
(flatten for instance)
- We do not have "standard correct way" to interface with functions
like factor, symplify, had, get, atoms, etc, so it is hard to be sure
that you have implemented everything
- As a personal preference, I actually like "garbage in - garbage out"

But this way has been proven to work in matrix expressions and is also
used in Pernici's tensor module, so I am not against it.

[Matthew Rocklin]

> A non trivial amount of functionality needs reimplementation (flatten for instance)

True, however I think you'll find that flatten does things that you don't want.  Much of the Expr code makes scalar assumptions that might not apply.  I tried to generalize some of my work in matrix expressions into the file sympy/strategies/rl.py .  This contains a set of generic functions for flattening, distribution, glomming together, removing identities and unpacking.  They were designed to be very general and make minimal assumptions.

- We do not have "standard correct way" to interface with functions

like factor, symplify, had, get, atoms, etc, so it is hard to be sure

that you have implemented everything

Some of these things like simplify can call on _eval_simplify.  Some of them you get from Basic.  Some of them, you're correct, you don't have at all.  What does it mean to factor vectors exactly?

Fortunately, you also don't have things like as_numer_denom, which Exprs have but which make little sense on Vectors

> But this way has been proven to work in matrix expressions and is also

used in Pernici's tensor module, so I am not against it.

In contrast the other way has proven to work in physics.quantum.  They're both reasonable.

[Stefan Krastanov]

> What does it mean to factor vectors exactly?

I was thinking about a*unit_vector + b*unit_vector -> (a+b)*unit_vector

But as you said there are many others that do not make sense for vectors.

> In contrast the other way has proven to work in physics.quantum. They're
> both reasonable.

And also in diffgeom.

In any case, both of this methodologies seem ok to me (create new
Add/Mul or reuse Add/Mul), but both of them are not what is proposed
in general terms on the wiki describing the API (a ScalarField and
Vector what eagerly eat anything that operates on them, i.e.
ScalarField+ScalarField -> ScalarField).

[Matthew Rocklin]

On Tue, Jun 11, 2013 at 10:14 AM, Stefan Krastanov <krastano...@gmail.com> wrote:

> What does it mean to factor vectors exactly?

I was thinking about a*unit_vector + b*unit_vector -> (a+b)*unit_vector

My solution was the generic version in sympy/strategies/rl.py.  This particular operation is handled abstractly by glom (which should perhaps be renamed to factor?)  I tried to cover all of the standard operations.  The goal was to not include any generic code in matrices/expressions.  Most of the generic stuff was pushed upstream to rl.py.  The hope was that it could one day serve as generic replacement for the scalar specific stuff in the Expr codebase.

 

But as you said there are many others that do not make sense for vectors.

> In contrast the other way has proven to work in physics.quantum.  They're
> both reasonable.

And also in diffgeom.

In any case, both of this methodologies seem ok to me (create new
Add/Mul or reuse Add/Mul), but both of them are not what is proposed
in general terms on the wiki describing the API (a ScalarField and
Vector what eagerly eat anything that operates on them, i.e.
ScalarField+ScalarField -> ScalarField).

Ah, I see.  In that case I'll then just say that container classes are good.  They let you describe computations without performing them.  You want this unless you're absolutely sure you know how to handle every possible input in exactly the right way.  Usually this isn't the case.

[Gilbert Gede]

Prasoon,

That option does not work for ReferenceFrame. If I have a ReferenceFrame A, and another ReferenceFrame B, and want to set the angular velocity of B in A as

w_1*b_1 + w_2_* b_2 + w_3 * b_3

I need access to the basis vectors of B, which I won't have until after the frame is initialized. Another situation that can occur is the angular velocity being defined as:

w_1*b_1 + w_2_* b_2 + w_3 * b_3 + w_k * a_1

I don't see any way around ReferenceFrame being mutable if you want to maintain the same functionality that is in mechanics now.

As for having a separate CoordinateSystem and ReferenceFrame: I think it is better to keep these explicit ideas separate if writing two classes doesn't result in much code overlap.

Stefan, regarding your concerns about Vector eating everything up - isn't that how it should be? Once you are dealing with Vectors, won't most operations will result in a Vector? Matthew's suggestion to implement our own VectorAdd/VectorMul seems like it might work though.

Also, is there a way to blacklist what functions are allowed to perform on Vectors?

[Matthew Rocklin]

On Tue, Jun 11, 2013 at 2:16 PM, Gilbert Gede <gilbe...@gmail.com> wrote:

Prasoon,

That option does not work for ReferenceFrame. If I have a ReferenceFrame A, and another ReferenceFrame B, and want to set the angular velocity of B in A as

w_1*b_1 + w_2_* b_2 + w_3 * b_3

I need access to the basis vectors of B, which I won't have until after the frame is initialized. Another situation that can occur is the angular velocity being defined as:

w_1*b_1 + w_2_* b_2 + w_3 * b_3 + w_k * a_1

I don't see any way around ReferenceFrame being mutable if you want to maintain the same functionality that is in mechanics now.

As for having a separate CoordinateSystem and ReferenceFrame: I think it is better to keep these explicit ideas separate if writing two classes doesn't result in much code overlap.

Mutability makes me sad  :-(

I haven't read enough to comment on all of this any degree of intellect; I'm merely promoting principles.

 

Stefan, regarding your concerns about Vector eating everything up - isn't that how it should be? Once you are dealing with Vectors, won't most operations will result in a Vector? Matthew's suggestion to implement our own VectorAdd/VectorMul seems like it might work though.

Also, is there a way to blacklist what functions are allowed to perform on Vectors?

Blacklisting in inheritance is possible but apparently looked down upon.  If you use Expr container classes (Add, Mul) it's likely that these methods will eventually be called upon in any event.

http://stackoverflow.com/questions/16749312/remove-attribute-from-subclass-in-python

 

[Aaron Meurer]

I am catching up here, so just some random comments on things that
were mentioned:

- Regarding the quantum module, they actually went the opposite route
from MatrixExprs. Originally, they used their own QAdd, QMul, etc.,
but then later it was changed to use the core Add, Mul, etc. I don't
remember the exact reasons. Brain can comment better.

- I agree with Stefan about garbage in, garbage out. I'm assuming
"garbage" here means things like vector + scalar. It's true that the
core does not currently let you enforce things like this (issue 1941
again), but I think it's fine. Just create expressions that are scalar
linear combinations of vectors. Then in the vector superclass, write a
little helper function that parses this into a more workable format,
and does validation. If you want, you could even perform some magic
to make this as painless (i.e., automatic) as possible.

- Regarding mutability, don't do it. It breaks the SymPy conventions,
and there's really no good reason to have it. Unless you are dealing
with large data structures that should not be copied in memory,
mutability is just the difference between

expr.operation_inplace()

and

expr = expr.operation()

The benefit to immutability is that you don't have global state (and
if you don't think that global state is a bad thing, then you've never
accidentally written "def func(x=[]):"). It's also easier to write
branched code flow, because you don't have to explicitly copy your
objects (ever modified a list while looping through it?).

- What does your planned API look like regarding functions and
methods? If it is method heavy, then making everything a Vector will
be more important.

- Unrelated, but I'm curious. What will Vector*Vector do?

Aaron Meurer

[Stefan Krastanov]

> I just feel that we are getting stuck with and spending too much time on things that we shouldn't be taking this much time. Therefore, I request again to all parties involved : Please voice your views/problems with these structural changes suggested so that we can move ahead with things.

I understand that you are eager to start writing code, but you should
understand that having an actual design is a prerequisite for writing
useful code. There are a lot of SymPy conventions that you and Sachin
have not learned about yet and that are not respected in your designs.
If we proceed that way most of what you write, even if functional will
not be of great use for sympy.

Anyway, if you two want to move faster, just contact each-other as I
have suggested and propose a common API. I have made a few suggestions
in https://github.com/sympy/sympy/wiki/Vectors-EM-framework#another-proposal.
These suggestions fix the issues with tree traversal, scalar+vector
errors and immutability but do not address interoperability with
mechanics for instance. No ScalarField and no Vector classes are
necessary for what I have written. Reuse whatever you find useful.

And remember: when you start codding does depend on you - show an api
design that you both have previously agreed on and that addresses:

- immutability
- tree traversal and sympy conventions
- scalar+vector error detection
- interoperability with mechanics

[Prasoon Shukla]

Well alright. Since the discussion is supposed to be done publicly, I think Sachin and I can discuss this on IRC. So, I have written an email to him inviting him to IRC, today 21:00 IST.

Also, the API suggested by Stefan looks like something to go forward with, albeit with some changes.

For the mutability issue : Just as Aaron suggested, we do really need the objects to be mutable. Doing something as

`expr = expr.method()`

also works just fine.

Then there is the problem of ScalarField + ScalarField == ScalarField. I think as long as we use different variables in different frames, we can just drop the whole concept of a ScalarField class.

Then, we need to address the problem with Vector. What does the Vector inherit from? Taking a hint from the tensor module, I think we can have a VectorExpr class that inherits from Basic so that

x == x.func(*x.args)

holds. Vector shall then inherit from VectorExpr.

Then we can have VectAdd and VectMul classes as Matthew suggested which inherit from VectorExpr as well. This would take care of garbage input problems.

Till the time we discuss this on IRC, can I have your opinions on this approach?

[Stefan Krastanov]

> For the mutability issue : Just as Aaron suggested, we do really need the
> objects to be mutable. Doing something as

If I understand correctly, Aaron's main point was that mutability is
useful only in certain cases and I do not think that the current
module requires it. Most of the arguments in favor of mutability were
about how to relate two reference frames (both static translation and
orientation as well as relative velocity and angular velocity). I
think that the example I provided addresses these issue and mutability
is not necessary.

> Then there is the problem of ScalarField + ScalarField == ScalarField. I
> think as long as we use different variables in different frames, we can just
> drop the whole concept of a ScalarField class.

As I have said, I am in favor of that. Moreover, using the same
variables in different frames seems like a bad idea anyway.

> Then, we need to address the problem with Vector. What does the Vector
> inherit from? Taking a hint from the tensor module, I think we can have a
> VectorExpr class that inherits from Basic so that
>
> x == x.func(*x.args)
>
> holds. Vector shall then inherit from VectorExpr.

Why not just follow the methodology used for scalar fields? I.e. just
creating symbols to be used as base vectors. If the issue is only
about printing then subclass Symbol.

Concerning the garbage input, you can do verification when you call a
function on your garbage expression.

[Sachin Joglekar]

Lets get together on the irc one time today. Lets settle this so that we can come to a consensus. Prasoon and I will be there in 1.5-2 hours. If possible, @Stefan, do come online

--
You received this message because you are subscribed to a topic in the Google Groups "sympy" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sympy/t-Je0beTIIU/unsubscribe?hl=en-US.
To unsubscribe from this group and all its topics, send an email to sympy+un...@googlegroups.com.

[Prasoon Shukla]

On Wednesday, June 12, 2013 7:52:34 PM UTC+5:30, Stefan Krastanov wrote:

> For the mutability issue : Just as Aaron suggested, we do really need the
> objects to be mutable. Doing something as

If I understand correctly, Aaron's main point was that mutability is
useful only in certain cases and I do not think that the current
module requires it. Most of the arguments in favor of mutability were
about how to relate two reference frames (both static translation and
orientation as well as relative velocity and angular velocity). I
think that the example I provided addresses these issue and mutability
is not necessary.

Oh! I'm really sorry :P I meant  

> ... we do not really need the objects to be mutable. Doing something as ....

Little typo, big confusion :)

[Jason Moore]

Stefan,

I just read through your proposed API. In general it seems reasonable but Luke, Gilbert and I spoke last night on irc (logs: http://colabti.org/irclogger/irclogger_log/sympy?date=2013-06-12) about how to deal with immutability/mutability, subclassing from Basic/Expr, and garbage in/out. The mechanics module has a unique use case where we have to track the same tree of direction cosine matrices for each reference frame involved. Currently we hide this tree on the ReferenceFrame objects and it is a mutable entity. If you add another ReferenceFrame to the tree then other ReferenceFrames need to know about this addition direction cosine matrix definition. For example, in your API, when you use cart_A after you've created cart_B it should reflect that it is rotated with respect to cart_B. As it stands, this wouldn't happen if the frames are immutable. You could imagine having this tree held in your world object but it would also need to be mutable.

I think we, @mechanics, will be fine with most of what you all proposed if we can clearly see that ReferenceFrame (plus velocities, accelerations) will be properly tracked. Our module's cornerstone are these trees that relate reference frames together.

Some summary from our conversation:

1. We're fine with immutability and subclassing from Basic if there is a mechanism in place that can handle reference frame relationships correctly and accurately.

2. We're more found of garbage in, throw an error for funny things like vec * vec or scalar + vec. Seems like you have that in your api proposal.

Jason
moorepants.info
+01 530-601-9791

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.

[Prasoon Shukla]

@Matthew Rocklin:

Can you please tell what exactly would be the advantage of having a separate VectAdd class? I mean I don't see the problems with using the current Add and Mul classes.

[Matthew Rocklin]

I don't remember enough of the details to give you a comprehensive answer.  Here are a few thoughts that come to mind. 

You won't be able to change much of the Add and Mul code to handle Vectors.  For example you won't be able to add shape checking at construction time.  This forces you to follow the "Garbage in, Garbage out" philosophy (which isn't necessarily bad)

In [18]: A = MatrixSymbol('A', n, 1)

In [19]: B = MatrixSymbol('B', m, 1)

In [20]: x = Symbol('x')

In [21]: Add(x, A, B)

Out[21]: x + A + B

In [22]: MatAdd(x, A, B)

TypeError: Mix of Matrix and Scalar symbols

In [23]: MatAdd(A, B)

ShapeError: Matrices A and B are not aligned

Add and Mul inherit from Expr.  The Expr codebase in general assumes that its elements are scalars.  This pops up in many of the methods on Expr classes. 

In [26]: Mul(A, A, B)

AttributeError: 'MatrixSymbol' object has no attribute 'as_base_exp'

I really don't think I should have to implement as_base_exp for matrices but I do if I want to use Mul.  Actually this sort of thing usually isn't so bad, the defaults you get from inheriting from Expr are reasonable.  I remember having lots of trouble with some of the more aggressive scalar functions.  People kept using expand and factor on matrix expressions.  They sometimes got errors and filed bug reports that were very difficult to resolve without completely rewriting these functions to be more general.

Using Mul and Add is fine.  It'll work for a while.  Every once in a while you'll come across a wart.  I'm personally biased against the Expr codebase; I think it's the result of a lot of glomming on of special cased code.  By my aesthetic MatrixExprs are a lot cleaner and more extensible.  I'm personally happy with the choice not to inherit from Expr.  Other groups have gone other ways though.  You could look at physics.quantum if you wanted to see how Exprs work on non-scalars.  

Also, just a thought, but if you're just using vectors you could presumably build off of MatrixExprs.  The container classes (MatAdd, MatMul) make very few assumptions about the things they hold other than that they are Basics and have a shape.

On Thu, Jun 13, 2013 at 1:52 AM, Prasoon Shukla <prasoon...@gmail.com> wrote:

@Matthew Rocklin:

Can you please tell what exactly would be the advantage of having a separate VectAdd class? I mean I don't see the problems with using the current Add and Mul classes.

--

[Sachin Joglekar]

Personally, I like Mathew's idea of using MatrixExprs to represent vectors. This will not just ensure smooth functioning with the rest of sympy, but will also help in easing the 'garbage detection' process. Plus, representation (as well as stuff like pretty printing) would get easier.

The only 'oddity', once again, would be vectors that are an addition of bases from two different frames. Maybe use something on the lines of MatAdd?

--
You received this message because you are subscribed to a topic in the Google Groups "sympy" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sympy/t-Je0beTIIU/unsubscribe?hl=en-US.
To unsubscribe from this group and all its topics, send an email to sympy+un...@googlegroups.com.

[Matthew Rocklin]

> This will not just ensure smooth functioning with the rest of sympy

That's not necessarily true and is an important consideration.  I should note that matrix expressions have relatively few users (despite their awesomeness).  Your project is particularly valuable because of the integration with the @mechanics group.  That relationship should take priority.  It's worth making sure that the mechanics codebase is compatible with matrix expressions.  One important bit is that matrix expressions do assume all arguments are Basic and therefore immutable.  In terms of existing SymPy code they support ImmutableMatrix but not Matrix.

[Stefan Krastanov]

I do not see any problem with using MatExpr instead of Expr for this
GSoC project. The only difference from what we discussed on IRC
yesterday is that instead of using Symbol as a superclass for
BaseScalarField (i.e. x y z) and BaseVector (i.e. the unit vectors
along x y z) you will use MatrixSymbol for the vectors and Symbol for
the scalars.

However, I do not know well all the details around MatExpr - if you
choose this path, please explain why (and give actual examples, not
just statements like "garbage detection will be easier"). I insist on
this, because currently the only problems with using Expr seems to be
garbage detection, which can simply be postponed to function calls
(much like the use of sympify).

[Sachin Joglekar]

@Stefan, I think we should let Prasoon put forth his idea of the API as discussed yesterday first. We will discuss it and ensure that we cover all ground. If alls well with it, we can move forward with that. Just a thought.

[Gilbert Gede]

I'd like to reiterate the concerns of the mechanics people, in that there needs to be a more serious consideration of ReferenceFrames and the end-users' ability to make changes to them after they are instantiated. It's easy to gloss over this problem when you're looking at one or two frames, but when you start looking at 10 or 20 (which is not uncommon for multibody systems) the user effort in managing the relationship tree has to be considered.

I'm also not completely comfortable with the garbage-in/garbage-out philosophy being used with something (physics.mechanics) that is an engineering tool. It seems like a few checks on the "validity" of generated vectors are being considered, but I agree with Stefan in that you need to be more explicit in how some of these details will be handled, and I am concerned with how robust they are. If we can successfully subclass from Expr, that would be great. I think it needs to be considered against the downsides though.

-Gilbert

[Matthew Rocklin]

It seems like mechanics applications operate on a deep tree of objects.  Is it correct to put all of these inside SymPy?  The lower levels should certainly be SymPy objects but perhaps some of the higher levels are sufficiently application-centric and require enough performance that they should break with SymPy's model of expression trees.

Again, I know very little of the details of what you all are working on, but maybe there is a place in the organization at which the SymPy model stops and generic Python code takes over.  This is just a though, please don't consider it too strongly.

[Prasoon Shukla]

Alright then. Looking at the answers that Matthew has provided, here's what we will do:

1. Coordinate systems and reference frames will not be implemented separately.
   
2. There will be RefFrame class from which specialized classes will inherit, such as RefFrameRect, RefFrameSph etc.
   
3. RefFrame objects will be immutable.
   
4. Vector object will be immutable as well.
   
5. Vector class will subclass from Expr.
   
6. We will have a VectAdd class which inherits from Add and adds required methods.
7. Generalized vector fields will be of the form : VectAdd( Mul(scalar, base_vector1), Mul(scalar, base_vector2), Mul(scalar, base_vector3) ) where base_vectorx is an instance of the Vector class.

Notice that we are using the VectAdd class – this class inherits from Add and supports methods required on vectors. Initially, these methods (eg, curl, div, integrate etc) were supposed to go in Vector class but since we have decided that Vectors will just be base_vectors - therefore the VectAdd class will take care of this.

So, is everyone okay with this?

[Stefan Krastanov]

Why the asymmetry Mul vs VectAdd?

> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.

> Visit this group at http://groups.google.com/group/sympy.

[Prasoon Shukla]

Yes of course. We should also have a VectMul class that would check for vector*vector and provide methods on scalar*vector quantities.

[Stefan Krastanov]

Any additional methods that are provided will _not work_ if they do
not have access to the RefFrame objects (in order to know how to
transform between coordinates). Your proposal does not explain how
these RefFrame objects will be accessed.

Because of this, the only advantage of subclassing Add and Mul that I
see is that it will have garbage checks (which as I mentioned can be
postponed to precondition checks in any of the functions to be
implemented). On the other hand it keeps all the disadvantages
mentioned by Matthew (Add and Mul make many assumptions, that you will
keep by subclassing them).

So, if you want new container classes for any of the reasons mentioned
by Matthew, subclassing Add and Mul does not seem the right way - you
should start from Basic as it is done in MatExpr. On the other hand,
as I mentioned in the last two paragraph, I do not think you have
given appropriate reasons for creating new container classes.

If the main reason for new container classes is garbage checks you can
just use MatExpr (its add, mul, matsymbol, etc).

This MatExpr approach has another unrelated advantage: the author of
MatExpr (Matthew) has both addressed most (probably not all) of the
issues that are problematic for diffgeom, mechanics or quantum in
terms of following sympy conventions and he is still very active
within the project so you will be able to question him for details
(sorry to volunteer you like that, Matthew ;)

[Sachin Joglekar]

Also, Gilbert and I had a small question. How do you plan to implement expression of Vectors in different frames? Also, the relationships between the frames? Could you give an API and the basic idea of how those methods will work?

The main concern the mechanics group has, is the modeling of relationships among frames themselves, and those with vectors. Clearing up this part would make things much clearer.

--
You received this message because you are subscribed to a topic in the Google Groups "sympy" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sympy/t-Je0beTIIU/unsubscribe.
To unsubscribe from this group and all its topics, send an email to sympy+un...@googlegroups.com.

[Prasoon Shukla]

Any additional methods that are provided will _not work_ if they do
not have access to the RefFrame objects (in order to know how to
transform between coordinates). Your proposal does not explain how
these RefFrame objects will be accessed.

But they do.

I'll try to explain the structure one more time.

Basic->Expr->Vector

Every Vector has an attribute, `ref_frame`, that specifies the RefFrame object that we are using for that vector. (Also,Vector objects are just basis vectors)

Now, let us say we have a vector field like this:

v = x * e_x + y * e_y + z * e_z 

Without the VectAdd class, v is just an Add and we would have none of the methods we need for operating on vectors. Instead, if we have v as a VectAdd object, we can extend the Add class to provide methods that are required on vectors.

Initially, when we weren't following the basis vector based approach, v would have been a Vector object. In that case, we would have had these methods defined in the Vector class. But since we are following the approach where a Vector is just a basis vector, therefore we need something to deal with vector fields that are not instances of Vector. Hence the requirement for these two classes.

Here's the hierarchy:

Basic->Expr->Add->VectAdd 

Basic->Expr->Mul->VectMul

I hope this clarifies the issue. Any more problems?

[Stefan Krastanov]

I see. This makes sense. I dislike the creation of a new class just
for the sake of adding helper methods (it mixes functions and data in
an unaesthetic way), but this is just my preference, so I will not
raise an objection if this is what you prefer.

The big issue now is how well will this work with mechanics. I am
against making big compromises in the future architecture in order to
accommodate mechanics, but a well working interface must be created.

While RefFrames will have to be immutable to be used within Sympy
trees (or at least have consistent hash and eq as is done with the
mutable coordSystems from diffgeom) you might want to think about some
nice (mutable) tree of RefFrames structure that will accomodate
mechanics.

[Prasoon Shukla]

On Friday, June 14, 2013 3:39:20 PM UTC+5:30, Sachin Joglekar wrote:

Also, Gilbert and I had a small question. How do you plan to implement expression of Vectors in different frames? Also, the relationships between the frames? Could you give an API and the basic idea of how those methods will work?

The main concern the mechanics group has, is the modeling of relationships among frames themselves, and those with vectors. Clearing up this part would make things much clearer.

Again, I really do not see what is so confusing. I'll try to be clear this time.

Each vector has an attribute, `ref_frame`, which specifies the frame in which the vector has been defined in. 

As for relationship between the frames are concerned, we would have methods defined inside the RefFrame class that would output the transform matrices when required explicitly. Or, when required implicitly, the matrices would just be used up internally.

For example, let us have a simple example. We have two frames A and B. Both in rectangular coordinates. B is rotated about the common x axis of both the frames by theta.

Now, let us say the user just wants to know the conversion matrix. He can just do something like:

B.convert_mat(A)

And we would output transformation matrix. This stuff is just like mechanics.

Let us add a bit more detail. Say the B frame is in spherical coordinates. Then, the transform matrix will also include the transform matrix from spherical to rectangular coordinates in the matrix product.

Now, let us see how it would be used internally by the VectAdd and VectMul classes. Let us say we have a vector in some reference frame. Then, if we want to express that vector in another frame, we would do:

vector.convert_to(ref_frame=frame)

which returns a new vector.

I hope that clarifies the issue of how things would work. More questions?

[Sachin Joglekar]

Ok. So I guess same will apply when we want to express some VectAdd in one frame.

And there is going to be a BaseScalarField class right? To enable expression of scalar fields in particular frames?

--

[Prasoon Shukla]

We will need to have symbols for different frames, so yes, we will have a BaseScalar like class that just instantiates the symbols to be used.

[Stefan Krastanov]

>> Also, Gilbert and I had a small question. How do you plan to implement
>> expression of Vectors in different frames? Also, the relationships between
>> the frames? Could you give an API and the basic idea of how those methods
>> will work?
>> The main concern the mechanics group has, is the modeling of relationships
>> among frames themselves, and those with vectors. Clearing up this part would
>> make things much clearer.
>
>
> Again, I really do not see what is so confusing. I'll try to be clear this
> time.

Prasoon, in this case I do not think the issue is someones confusion
with what the suggestion is. `mechanics` has a very specific goal and
it is optimized for it - easily traversing a very deep tree of
relations between reference systems. While indeed what you have
suggested is sufficient to make such traversals possible, one should
be wary of whether it will work _well_ with deep trees.

One issue that you will meet very soon is that expressing
cart->spherical->cart transformation will often fail to symplify to an
identity transformation (a lot of information must be specified for
stuff like cos(acos(cos(a))) to return cos or for sqrt(x)*sqrt(1/x) to
return 1).

While the above is probably not the exact concern of Gilbert and the
others in the mechanics team I am sure they have similar worries.

And do not forget that mechanics is fast at the moment. With your
methodology it might be hard to keep up with its performance.

[Gilbert Gede]

Yes Stefan, I think you are articulating our concerns well. And while I've never done the rect/spherical/rect transform, just going back and forth between rect. frames causes similar issues with trig expressions not collapsing.

I'd like to know more about the mutability of diffgeom's CoordSystems though. From reading the code, it looks like it is accomplished by storing the mutable information outside of args, and having the parent patch as an arg to the coordinate system? How robust does this end up being?

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.

[Jason Moore]

I wonder if it would be best to remove the concept of a reference frame from Prasoon's api design. Then the mechanics and electromagnetics modules will need to implement that on our side on top of the coordinate system classes.

Jason
moorepants.info
+01 530-601-9791

[Prasoon Shukla]

@Jason: I do not know if other's are willing to implement this.

Anyway, so coding period start in 2 days. We need to decide what to do. Right now, we have an API that appears workable. If there are going to be any more major changes (like supporting mutability or not having reference frames as part of the vector module as Jason suggested), please let me know them now so that we can work things out around that.

Thanks

[Sachin Joglekar]

I can implement Jason's idea, if you don't mind Prasoon. I had set some time in my timeline to work on the framework, and I can utilize it to implement the said functionalities via your API in mechanics. Since your module deals more with coordinates, I can deal with the 'time' and 'motion' part, with the help of the mechanics team. Some factors will have to be considered like motion of frames wrt each other, derivatives wrt time (which are different in different frames), etc. which aren't really a part of vector calculus, and they belong in the physics module.

About the API, Prasoon, I guess we can start work if the mentors are ready. We hadn't really talked much about time variables, but the coordinate functionalities may be done as you and Stefan have worked out. @Stefan, @Gilbert, how do you think we should proceed?

--
You received this message because you are subscribed to a topic in the Google Groups "sympy" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sympy/t-Je0beTIIU/unsubscribe.
To unsubscribe from this group and all its topics, send an email to sympy+un...@googlegroups.com.

[Stefan Krastanov]

I agree with Sachin: Prasoon continues with what he has described but
does not care about time variables or kinematics, only about static
orientation between reference frames. Sachin builds the time
dependence, time derivatives and kinematics on top of that.

Of course, you will probably address each others tasks from time to time.

If you are in the same city (or city district) you might plan some
coding sessions together. Also, while working, try to be on IRC so you
can talk to each other if necessary.

[Prasoon Shukla]

Alright then. How should we proceed with this differentiation? My initial idea of having a CoordSys class and a RefFrame class can perhaps be used.

[Sachin Joglekar]

I agree. What we can do is, Prasoon can start building a CoordSys class, and as Jason suggested, I will start working on a RefFrame class (with your API) - this class will be a subclass of CoordSys, specially built for the physics module.

Prasoon, could you clear up the wiki page, and paste your API design fresh? Try to post all the methods and attributes (I know you won't know all of them at the moment, keep writing as you go along). Similarly, I will keep adding to the page with the RefFrame API dealing woth time functions below your API.

I will just clear things up with the mechanics people and get to work in a day or two.

@Stefan, @Gilbert, okay with this?

On Sat, Jun 15, 2013 at 7:36 PM, Prasoon Shukla <prasoon...@gmail.com> wrote:

Alright then. How should we proceed with this differentiation? My initial idea of having a CoordSys class and a RefFrame class can perhaps be used.

--

[Stefan Krastanov]

Actually I do not see any issue. You still have to do all coordinate
transformations, translated and rotated references, etc.

Just do not occupy yourself with time variables, angular velocity and
linear velocity of the references.

On 15 June 2013 16:06, Prasoon Shukla <prasoon...@gmail.com> wrote:
> Alright then. How should we proceed with this differentiation? My initial
> idea of having a CoordSys class and a RefFrame class can perhaps be used.
>

> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an

[Stefan Krastanov]

I actually am very much against separating coordsys and refframe
again. I really do not see the point. Prasoon just continues his work
on the api he suggested last time and Sachin adds time derivatives.
These two have nothing in common.

[Sachin Joglekar]

Not a problem! Ok, what about the module/class structure? Where will the new code be located?

[Gilbert Gede]

I am against this.

It is far too easy to just say we'll figure out the details later. It will only lead to more problems when attempting to implement the ignored functionality.

When is the next time everyone can meet on IRC?

[Prasoon Shukla]

@Gilbert: Can you come over right now to the IRC?

[Prasoon Shukla]

Okay so I just talked with Sachin on this. As Gilbert says, we should plan before we start work.

Initially, my proposal had a distinction between a CoordSys class and a RefFrame class. We are at a similar point now. Nevertheless, the two ways we thought this could work was by:

a) Subclassing the coordinate system class from the vector module in the appropriate physics module. Or,

b) Having an attribute, `coord_sys`, in the class defined in the physics module through which all the functionality of the vector module is accessed.

Stefan has pointed out some problems with the second approach though. So, should we go with (a) or is there a better way to do this?

[Stefan Krastanov]

Maybe I am simply missing something, but I do not understand why
anything more complicated than option (a) would be considered. Prasoon
just implements vector fields as they are meant in math (with static
orientations) and Sachin adds a parametrization to these orientations
(the parameter being time) that happens to be used in physics.

> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an

[Stefan Krastanov]

Considering IRC - I can attend on Monday for instance (in the
morning for the West Coast, evening in Europe and India)

[Prasoon Shukla]

Alright then. We will go with option (a). @Gilbert and @Sachin, can you come on the IRC tonight?

[Sachin Joglekar]

Yeah. I can be there at 9:30

On Sun, Jun 16, 2013 at 6:33 PM, Prasoon Shukla <prasoon...@gmail.com> wrote:

Alright then. We will go with option (a). @Gilbert and @Sachin, can you come on the IRC tonight?

--
You received this message because you are subscribed to a topic in the Google Groups "sympy" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sympy/t-Je0beTIIU/unsubscribe.
To unsubscribe from this group and all its topics, send an email to sympy+un...@googlegroups.com.

[Prasoon Shukla]

I, Gilbert and Sachin met up today on the IRC. We talked about how exactly will we implement the Reference Frames - that is - how the data regarding position and orientation will be stored. Here's the conclusion we reached:

1. The vector module will have a RefFrame class that supports positioning and orienting frames while Sachin will separately implement time support for a time variable.

2. For positioning and orienting the reference frames : we will store the position/orientation with respect to a '0 reference' - which effectively means that we aren't saving what frame we are positioning/orienting against.

3. Finding out relative position is a trivial task under such an implementation. The problem comes with finding proper relative orientation between two frames, as I shall describe in point 4.

4. We will save a direction cosine matrix (dcm) that stores the orientation of each frame with respect to the zero reference. Then, going from one frame to other, the dcm matrices will follow such an order: frame_A -> global (zero reference) -> frame_B. The problem with this is that many expressions might not decompose/collapse as they should in this matrix.

5. To address this problem of collapsability, we have decided to do some research. Since the elements of the matrix will be only terms involving sines and cosines, therefore, we will try to investigate ways by which we can ensure that the expressions collapse properly. Should we achieve a consistent result with these methods, we will proceed with above mentioned scheme.

6. There were some other ideas too, involving storing DCMs into cache - but that will come later; we first need to make sure that the expressions collapse. But this will also be researched in the time tomorrow by me and Sachin.

There were other ideas that we discussed but I have refrained from adding them here so as to not dilute the content. The discussion was carried on the IRC so the logs are available.

The three of us will meet at 9:30 pm (Indian time) tomorrow again. Hopefully, Stefan will be able to join us as well.

[Stefan Krastanov]

Thanks for the summary. I will be on IRC at the specified time.

> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an

[Prasoon Shukla]

@ Everyone : I couldn't do any work today - I had to make an unexpected trip to a certain city and back. Traveled a total of 15+ hours and got stuck in a 3 km long traffic jam. I had thought that I'll be able to work while travelling - I have never been more wrong :(

Anyway, I just got back. I'm writing the API and will upload it before I go to sleep. I'll make up for lost time in the coming days of this week.

[Prasoon Shukla]

Okay, so since there hasn't been any activity on this thread, let me post the link to my code :

http://github.com/prasoon2211/sympy/tree/vector

[Stefan Krastanov]

Make a pull request please, so we can have diffs and line comments
visible for everybody.

[Prasoon Shukla]

Alright then, I'll make one in a few hours.

On Tuesday, June 4, 2013 2:46:22 AM UTC+5:30, Stefan Krastanov wrote:

Today we had the first discussion with Prasoon and Sachin about their projects.

We did not progress much but at least we outlined the two general approaches that we can use for these modules (specifically for creating vector fields). I will give them somewhat arbitrary names here:

- the `mechanics` way - having a Vector class that keeps all the information about the field and it is not part of expression trees in the way Basic and Expr are. For instance Vector(something along cartesian.x)+Vector(something along spherical.r) will result in Vector(complex internal structure).

- the `diffgeom` way - having base/unit vectors and building all the rest in terms of their linear combinations (all this expressed as sympy expressions).

Prasoon and Sachin did not have the time to look at the example problem that was given in the previous email yet (no harm done there, there is still some time before the official starting date). Probably this will be the subject of our next discussion.

The next discussion was scheduled for tomorrow. After that I suggest that we keep most of the discussions to the mailing list and the gihub wiki and meet on irc / realtime wikis / google docs / etc  once a week.

Stefan

[Prasoon Shukla]

@Gilbert, @Sachin, @Stefan : I am facing a bit of a problem that I have written about in this week's blog here. Please see the second part of the blog where I mention about a cache to hold relationships between coordinate system and base scalars and base vectors. Please comment your views here.

[Sachin Joglekar]

Two things-

1.) I think we should stick to the current mechanics framework's ideology of not letting the user initialize a base scalar/vector. As I have commented on your PR, we can define __getattr__ and __getitem__ methods for the CoordSys class, which would enable a user to use the 3rd (say) basis vector with R.3, and the 3rd base scalar with R[3]

2.) Now, coming to the creation of base vectors/scalars. I like your idea of caching them. The main reason I feel this should be done is, we don't really need two or more instances of a base vector/scalar in one run. So, we can define a helper method to initialise a base scalar-

_generate_base_vector(coordsys, index)

Hence, when the user says R[2], R's __getitem__ will call '_generate_base_vector(self, 2)'. Now we can use function memoization for this method. We can create a custom cache which stores a tuple of the coordinate system and the index as the key. So the next time the same method is called for R[2], the pre-generated instance will be returned.

Same can be done for basis vectors. I had written a blog post about function memoization, and I guess it will work here with some modification.

1) and 2) will eliminate any ambiguity regarding base vectors/scalars and also provide a cleaner interface for the whole operation.

On Sat, Jul 6, 2013 at 4:19 PM, Prasoon Shukla <prasoon...@gmail.com> wrote:

@Gilbert, @Sachin, @Stefan : I am facing a bit of a problem that I have written about in this week's blog here. Please see the second part of the blog where I mention about a cache to hold relationships between coordinate system and base scalars and base vectors. Please comment your views here.

--
You received this message because you are subscribed to a topic in the Google Groups "sympy" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sympy/t-Je0beTIIU/unsubscribe.
To unsubscribe from this group and all its topics, send an email to sympy+un...@googlegroups.com.

[Gilbert Gede]

I think Sachin has some good suggestions. I would advocate something similar (which is also similar to what mechanics currently has): why not just create the base scalars and basis vectors upon initialization of a coordinate system? There should only be one set of basis vectors and base scalars for a coordinate system, so I think it makes the most sense to create them upon initialization, then return them as needed. I think this gets around the need to cache them, although you still can and I believe should cache them (for example, all sympy symbols are cached).

I also agree on sticking to the current mechanics style of not having the user directly call the Vector() class, but then again, I am certainly biased towards the mechanics approach...

This was certainly a decision based off my opinion, but it just seemed more "correct" and natural to only create vectors from scalar/basis-vector combinations and not by initializing them directly.

Also, what are people's thoughts on the following:

Instead of

some_vector.express_in(some_coord_sys)

some_vector.diff(x, frame=some_ref_frame)

using

some_coord_sys.express(some_vector)

some_ref_frame.diff(some_vector, x)

The first option is similar to how mechanics is currently set up, and you are returning a Vector from a method of Vector.

I feel like the second option is a little cleaner, but you are returning a Vector from a CoordinateSystem/MovingRefFrame method.

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.

[Sachin Joglekar]

Why should returning a Vector from a frame method be a problem?

I prefer the second method, because it will be easier for the coordsys to operate on the vector its been provided (substituting base scalars, multiplying by dcm etc), than the method of vector making a copy of itself and operating on it. Conceptually and from the point of view of OO-design, it dosent make a difference. They mean the same.

Also, I will re-ask a doubt I had about Prasoon's code...Prasoon, why have you disallowed base scalars from appearing in velocity/accln vectors provided on initialization?

[Prasoon Shukla]

Sachin, I responded to your question in the PR.

I think I'll go with Gilbert's approach. As far as the second part was, we will take the first approach - the reason is that we have three different classes that handle vectors - and we need to be able to apply those operation on each class.

Also, there's another problem. Right now, we can do Vector.coord_sys and get the CoordSys object. But now, we can't do this because the CoordSys object hasn't been constructed at the time of initializing the vectors.

[Prasoon Shukla]

Hmm. It seems Python allows one-to-one relationships between objects. So, I did this:

In [64]: class Vector(object):

   ....:     def __init__(self, coord_sys):

   ....:         self.coord_sys = coord_sys

   ....:         

In [65]: class CoordSys(object):

   ....:     def __init__(self):

   ....:         self.vector = Vector(self)

   ....:

In [66]: c = CoordSys()

In [67]: c

Out[67]: <__main__.CoordSys at 0x97b46cc>

In [68]: c.vector

Out[68]: <__main__.Vector at 0x97b42ac>

In [69]: c.vector.coord_sys

Out[69]: <__main__.CoordSys at 0x97b46cc>

And this seems to work just fine.

[Prasoon Shukla]

I'm having some problems with implementing the `express` function. This function would be called on a vector and would take a CoordSys object as and argument. The function should return a vector expressed in the CoordSys provided to this method.

Now, as a first step, we need to take care of the position.

Say, we have two coordinate systems, C0 and C1, each having their origin at a different point in space. Now, we have a vector field in C1, say F.

Now, we know the position coordinates of each of the frames in the global frame. But, these coordinates could have been expressed in any coordinate system. So, here's the algo:

1. Convert the position of both C0 and C1 into rectangular coordinates.

2. Subtract both the positions term wise. This way, we will have the difference of coordinates between C0 and C1.

3. Now, we can have relations such as: 

x0 = x1 + a  

y0 = y1 + b   

z0 = z1 + c

-----------------Equations (*)

Where the variables suffixed with 0 are in C0 and similarly for C1.

4. Now, we need to convert F from whatever coordinate system it was expressed in to rectangular so that we can substitute  the variables from system (*). Also, we need to convert the base vectors to rectangular coordinates - because they stay constant throughout space.

5. Now, we have F in C1 expressed in rectangular coordinates. Now, substitute from system (*) into F. Base vectors remain the same because F has already been converted to rectangular coordinates.

6. Now effectively, we have F in C0, albeit in rectagular coordinates.

Now, we can proceed from here to take care of orientation or change the coordinates back to some other coordinate system (spherical for example).

Now, I want you guys to verify whether the procedure I have written above is correct or not.

-----------------------

Now, we have to take care of the orientation. Now that we have positioned ourselves at C0, let as say we need to rotate C0 to another coordinate system, Cr, where C0 and Cr have the same origin. We already have the DCM between C0 and Cr.

Now, I need to ask something. Does DCM represent the relation between coordinate variables or does it represent relation between coordinate variables? Or does it represent the relation between base vectors? Or both?

Because of the way the DCM is derived, I think that it is the relation should hold for both. Is that correct?

[Sachin Joglekar]

DCM helps in modeling of relationship between base vectors as well as base scalars.

The re-expression of a vector V in a frame A should proceed in three steps-

1) Substitute base scalars of all other frames except A, in terms of A, in V

2) Separate the vector we get from step1 into its constituents in various frames (B1, B2, ... Bn) using the separate() method, and then multiply [dcm of the corresponding frame (Bi) wrt frame A] by the respective constituent.

3) Add the vectors generated in step2

I guess you propose the system (*) for step1. However, I think there is an error with that. Consider a frame B, having some rotation and some pos_vector(P) wrt A. As we wont include base scalars in position vectors, expressing P in A wont be difficult (just steps 2 and 3 of above).

Hence, [xA, yA, zA] = [P.x, P.y, P.z] + A.DCM(B) * [xB, yB, zB]  ...(Column vectors, not row vectors)

The above equation would give xA, yA, zA in terms of xB, yB, zB, which can be used for step 1.

I think we should have a cached method _map_variables which does step1, as it can also be used to express scalar fields in any frame. _map_variables, as I have said earlier, can return a dict as {xA : f(xB, yB, zB), yA : f(...}.  This dict can be directly used for the 'subs' method.

(All the above applies for rect coords)

--

[Gilbert Gede]

Prasoon,

My first question is this - can position vectors be written using coordinate variables (base scalars in your codebase, I believe) - i.e. if the user has CoordSys "A" with coordinates ax, ay, az, can they define a new coordinate system using the position vector p = ax * e_A_x + ay * e_A_y + az * e_A_z? If the answer is no, the following should apply. If the answer is yes.... the second part won't work.

For the algorithm - I would recommend splitting it into 2 parts: basis vectors and basis scalars (coordinates). Your function/method would probably invoke both, unless a keyword is supplied indicating the user only wants to change one part.

For the first part, the code would be very similar to mechanics - you would iterate through every term in the vector, where each term will have a basis vector. Consider going from coordinate system B to coordinate system A with the DCM {^A}C{^B} (which I'll shorten to C): for a vector term which can be written as v1 * e_B_x, this would be expressed in the A coordinate system as v1 * (C(1,1)*e_A_x + C(2,1)*e_A_y + C(3,1)*e_A_z). This is pretty much the same as what mechanics does now, just not with matrix multiplications.

An alternative way to do this for a vector v is:

v = dot(v, e_A_x)*e_A_x + dot(v, e_A_y)*e_A_y + dot(v, e_A_z)*e_A_z

However, I'm not sure which method is faster/will result in cleaner expressions.

For the second part, you can do something similar. In the above example, lets use coordinates ax, ay, az and bx, by, by. Let's also say the position vector from A_o to B_o is p^{AB}. You can then do:

ax = (p^{AB}, e_A_x) + dot(e_A_x, bx*e_B_x + by*e_B_y + bz*e_B_z)

And so forth, changing the basis vector you are using to dot with. You can also take a similar approach to before, using the DCM directly:

ax = (p^{AB}, e_A_x) + {^A} C {^B}[1, :] * Matrix([bx, by, bz])

After you have these, you can just make a dict and use subs. You'll also obviously have to look for all present coordinates in the vector, and develop this relationship for all present (e.g. B->A, C->B->A, etc...)

As for just converting to rectangular coordinates as an intermediate step: I think it's OK to just go with that for now. Is that how the dot product between spherical basis vectors and rectangular basis vectors is going to work? Stefan, I think you have more experience with using non-rectangular coordinate systems than I do, what do you think?

I feel it's very important to _not_ just go to the global coordinates/orientation when doing re-expression or other similar tasks. It will really hurt the functionality and performance of the code. Try and just work with relative DCMs and coordinate transforms; e.g. for series of coordinate systems A -> B -> C -> D, find the DCM between each parent and child and multiply all the DCMs together to get the DCM from A to D, rather then going from D -> global -> A. Same for the coordinate transformations. Remember, users are going to be defining these positions/orientations relative to the parent CoordSys anyways, not to the global frame.

[Sachin Joglekar]

+1 to Gilbert's idea of not going to global frame. That's what I was doing initially, but we can get rid of that method by using a _frame_path function which will give the tree 'path' between any two frames. In this 'path', the pos vector and DCM of every frame wrt its former will be defined (as every ith frame will be a parent or child of the i-1th frame), so they can just be added. Refer to my code _frame_path and pos_vector_in methods.

--

[Stefan Krastanov]

@Prasoon, maybe I misunderstood what you suggest, but on first glance it seems it will work awfully in the following case:

coordinate systems:

A (carthesian) -> B -> many more -> U -> V

where U and V have the same origin and orientation but U is carthesian while V is polar.

How will something defined in U will be expressed in V according to your suggestion?