Re: [boost] Preview 3 of the Geometry Library

Subject: Re: [boost] Preview 3 of the Geometry Library
From: John Phillips (phillips_at_[hidden])
Date: 2008-10-22 02:48:42

Mathias Gaunard wrote:
> John Phillips wrote:
> 
>>   Every coordinate system was created for algebra, calculus and 
>> computations.
> 
> Yet when dealing with vectors and matrices, I'm pretty sure the 
> cartesian coordinate system is what everyone uses.

   It is usually the first thing to try, but it is not always the best 
idea. I won't claim that a large fraction of the computing community is 
conscious of the other possibilities, or of when they are useful. 
However, much of the scientific computing part of the community is 
familiar with multiple systems and reasons to use them.

> 
> 
>> Making the best choice of coordinates for a representation can speed 
>> calculations by orders of magnitude, or even in some cases make 
>> problems solvable that weren't in Cartesian coordinates (because of a 
>> lack of an appropriate representation for the solution). This is a 
>> huge deal in analytic calculations (See, for example the books by 
>> Morse and Feshback, or the early editions of Arfken where chapters are 
>> dedicated to finding and operating in the best coordinate system for 
>> the job. Also see facilities for this work in Maple or Mathematica.), 
>> but it is a sometimes under appreciated aspect of numeric calculation, 
>> as well.
> 
> I've heard that various analytical problems are much more easily 
> represented in other coordinates. But is really the kind of problems you 
> want a boost geometry library to deal with?
> 
> Apart from that, I've seen people use homogeneous coordinates (which are 
> very similar to cartesian) to avoid doing some divisions that might lead 
> to numerical instability.
> That should be an easy addition.

   I would like a boost library that serves as large a segment of the 
problem space as is reasonably possible. I haven't tried to design one, 
so I don't know if the extension somehow collides with the abstractions 
and implementation details needed for a clearly designed Cartesian 
library. I don't think so, with the small amount of thought I've put 
into it so far, but that isn't the same as a real study of the problem.

   Some problems are simplified incredibly by a good choice of 
coordinate system. As a trivial example, consider a circle in two 
dimensional coordinates. There is no properly defined function in 
Cartesian coordinates that produces a circle. (There is a relation, but 
a circle is not single valued in Cartesian coordinates, and so can't be 
written as a function.) The function in plane polar coordinates is r = c 
(The radius equals a constant), when talking about a circle of radius c 
centered at the origin. This is a function, and a particularly simple one.

   Physics problems in orbital mechanics and the quantum mechanics of 
electron orbitals, and mathematical problems of interpolation using 
radial basis functions are less trivial examples.

   In each case, the selection of a good coordinate system for the 
problem allows for a clear separation of concerns in the calculation. 
Analytically this can make it possible to apply techniques to the 
isolated components of the problem. Computationally, this reduces the 
multi-dimensional problem into a set of independent or nearly 
independent one dimensional problems. Since the computational effort of 
many common numerical techniques scales exponentially in the number of 
dimensions, this separation moves us out of exponential scaling and into 
linear scaling in the number of dimensions. I've never met a computer 
programmer who wouldn't approve of that trade.

> 
> 
>>   This is a conceptually similar problem to the discussions about unit 
>> systems from a couple of years ago. High accuracy and high performance 
>> calculations can't afford to have hidden conversion costs in 
>> computation and precision imposed on them. A well designed library 
>> should make such hidden conversions possible for those who want them, 
>> while also not making them mandatory for those who can't afford to pay 
>> for them.
> 
> I thought it was better for simplicity, dealing with different 
> geometries seems hard enough, but a prerequisite of the library since it 
> aims at dealing with GIS; and hopefully, there won't be too many 
> geometries to handle.
> On the other hand, different representations of points in the same space 
> are all equivalent, and one representation is enough to describe 
> everything we want to deal with, even if not in the most precise way.
> 
> The problem is that it seems hard for algorithms to consider all the 
> possibilities in existence for the points they're dealing with to be 
> represented and select the right way to compute something accordingly.
> Latitude/Longitude seems close enough to cartesian so I supposed that a 
> lot of code can be shared.

   All of the typical algebraic, differential and integral structures 
have well known translations into this world of different translations. 
For ortho-normal systems the key is the Jacobian. Even systems that are 
not ortho-normal can be approached with tensor methods. For some very 
common structures, like the dot product or the cross product, the 
algorithm doesn't even change. For differential and integral structures, 
the change is typically covered by throwing in a multiplicative factor. 
Details of this are available in most vector calculus books, or in the 
above mentioned books on mathematical physics from my previous note.

> 
> Maybe only doing vector operations and generalizing to any orthogonal 
> coordinate system is possible without much hassle by working with 
> normalized basis vectors, but wouldn't that be quite similar to 
> converting to cartesian everytime?

   No, it wouldn't. The Cartesian unit vectors are constants, and the 
Cartesian area and volume elements do not depend on current coordinate 
values. These facts are generally not true for non-Cartesian systems. 
(In plane polar coordinates, the same change of angle and the same 
difference in inner radius compared to outer radius sweeps out different 
areas depending on the actual value of the inner radius. As the inner 
radius gets bigger, the area swept out increases.)

   For problems where the natural description of the paths and shapes 
looks like rectangles, Cartesian coordinates provide a simple and 
accurate basis for calculation. However, if the natural paths and shapes 
don't look like parts of rectangles, Cartesian coordinates provide a bad 
description of the system and there is a need for substantial extra 
calculation to process the system in those coordinates.

   The answer is to look for a different system that matches the natural 
paths and shapes. (Better put, it shares the same symmetries as the 
problem.) These symmetries are automatically preserved just by writing 
the problem in the new coordinate system. There is no need to put 
computational time and effort into calculating them, or into sanity 
check functions that check for them and correct the results if they 
stray away from the required symmetries. In a best choice coordinate 
system, one characteristic change in the system can be well represented 
as one small change in a single system variable. (If it is a geometric 
system we are describing, the variables are coordinates.) This is good 
for the exact same reason that isolating effects in different modules in 
a computer program is good.

> I honestly have zero knowledge of non-cartesian coordinates, so I have 
> no knowledge of the vector theory behind it (and I have no knowledge of 
> how vectors apply to alternative geometries either).
> Maybe that could somehow also be applied to help dealing with other 
> geometries that are locally euclidian?
> 
> If there are ways to write algorithms generically in a 
> coordinate-system-independent way, and maybe even geometry-independent 
> way, then it would be great, for sure.
> But is that really possible?
> 

   This is the real question. My first guess is that it is possible for 
many algorithms. However, this is only an educated guess, and not the 
result of time spent producing real and practical implementations, so it 
could be wrong.

   In the cases where it isn't possible, we should try to do the same 
thing all generic programming does. Define a minimal set of concepts 
that must be modeled by the coordinate system for the algorithm to work 
correctly and embed such concepts and concept checks in the library.

                        John