Debugging rigid body code. Error with matrix multiplication code?

julien

unread,

Oct 22, 2014, 11:38:43 AM10/22/14

to sire-de...@googlegroups.com

Hi Chris,

I have discovered that the rigid body workspace code in my new version of Sire (sirejulien) doesn't in fact work correctly.

The issue is in a copy of rbworkspace that I made to implement some features necessary for cell theory calculations (rbworkspacejm).

By systematically comparing the output of the new code (sirejulien) with output from the old code (sirejm) I have established that the new code maps differently atomic forces onto a rigid body.

I have established that around line 645 of rbworkspacejm.cpp this line of code

bead_force = orient.inverse() * bead_force;

yields a different answer

with the old code, matrix * vector multiplication I have

-------

bead_force = ( -0.0675234, -0.0399004, 0.0450382 )

orient.inverse() "/ 0.769022, 0.62978, -0.109465 \

| 0.0616276, 0.0974031, 0.993335 |

\ 0.636245, -0.770642, 0.0360933 /"

result of multiplication

now... "( -0.0257305, -0.0811197, -0.0306174 )"

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

with the new code

---------

bead_force = "( -0.0675234, -0.0399004, 0.0450382 )"

orient.inverse() "/ 0.769022, 0.62978, -0.109465 \

| 0.0616276, 0.0974031, 0.993335 |

\ 0.636245, -0.770642, 0.0360933 /"

result of multiplication

now...."( -0.0819855, 0.0366903, -0.0105869 )"

----------

I have tested what the correct answer should be with numpy

In [24]: M

Out[24]:

matrix([[ 0.769022 , 0.62978 , -0.109465 ],

[ 0.0616276, 0.0974031, 0.993335 ],

[ 0.636245 , -0.770642 , 0.0360933]])

In [25]: V

Out[25]: array([-0.0675234, -0.0399004, 0.0450382])

In [26]: V*M

Out[26]: matrix([[-0.02573062, -0.08111964, -0.03061744]])

And this agrees with the old code.

Note that

In [33]: V*M.transpose()

Out[33]: matrix([[-0.08198556, 0.03669029, -0.01058692]])

Gives me the answer Sire produces with the new code, so it appears the matrix multiplication code in Sire is buggy?

Could this be in SireMaths/matrix.cpp ?

A diff between SireMaths/matrix.cpp revision 2749 (new code) and revision 2418 (old code) shows indeed significant changes in code in this file.

Christopher Woods

unread,

Oct 22, 2014, 12:25:39 PM10/22/14

to Sire Developers

Hi Julien,

Yes, it does look like this new code could be buggy. I had to make
changes because I found that the old matrix code had a different
row/column layout than other matrix code, which was causing problems
when I was writing things to perform molecular alignments or calculate
the volume of intersection of spheres. The problem was that I was
storing matrices by column but multiplying them by row (so the
multiplication was the wrong way around). I fixed this, but it looks
like there is still some code that is confused as to how the matrix is
stored. I'll look into this as a matter of urgency.

Cheers,

Christopher

> --
> You received this message because you are subscribed to the Google Groups
> "Sire Developers" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sire-develope...@googlegroups.com.
> To post to this group, send email to sire-de...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sire-developers.
> For more options, visit https://groups.google.com/d/optout.

--
---------------------------------------------------------
Christopher Woods
+44 (0) 7786 264562
http://chryswoods.com

Christopher Woods

unread,

Oct 23, 2014, 7:00:15 AM10/23/14

to Sire Developers

Hi Julien,

The matrix multiplication code is ok. I've created a matrix unit
testing script (attached) that shows that the matrix multiplication
code is working. I've also loaded up your matrices into ipython and
have found that the problem is because the rows/columns have been
swapped between the old and new versions (or, rather, that the
definition of row and column is swapped when you create a matrix).

If you write

orient.inverse().transpose() * bead_force

you will get the same answer with the new code that you got with the
old code. This has come about because the old code had a "bug" that;

m = Matrix( 1, 2, 3, 4, 5, 6, 7, 8, 9 )

would load up the matrix by column and not by row (e.g. it should be
that the first row is 1,2,3, second row is 4,5,6 and third row is
7,8,9, however the old code set the first row to 1,4,7, second row to
2,5,8 and third row to 3,6,9). I had written the code assuming that
the matrix was loaded by row and not column, because I like to write
matrix code such as;

m = Matrix( 1, 2, 3 \
4, 5, 6 \
7, 8, 9 )

(which then makes sense if the matrix is loaded by row).

Fixing the "bug" that the matrix was being loaded by column rather
than row has caused one of your earlier matrices to be loaded up as a
transpose. I suspect that 'orient' is being loaded as the transpose.
To confuse things, the reason I didn't spot this bug was because the
"print()" function on matrix was also buggy and was printing out rows
as columns (so the transpose of the matrix was printed, rather than
the actual matrix). Fixing both the "loading by column" bug and
"printing the transpose" bug at the same time is why your orient
matrix looks the same in both the old and new codes, but why orient is
in fact different. You should go back to the code that defines orient
and make sure that it is being loaded by row.

Cheers,

Christopher

p.s. I have made quite a bit of progress upgrading Sire to python
3.4.2. In the process I have bundled in zlib and openssl, so that we
can ensure that pip/easy_install will always work. I am also working
on a script that can be used, post-install, to update the RPATH in all
libraries and executables so that post-installed python modules can be
fixed, packaged and then shipped.

Christopher Woods

unread,

Oct 23, 2014, 7:00:41 AM10/23/14

to Sire Developers

forgot to attach the unit test script!

test_matrix.py

julien

unread,

Oct 23, 2014, 8:06:18 AM10/23/14

to sire-de...@googlegroups.com

Hi Chris,

Thanks for looking into this. Just to clarify, you are saying that in
the old code

expressions such as this

bead_force = orient.inverse() * bead_force;

are misleading because mathematically the correct operation is to
multiply the orientation matrix by the force vector. The old code gets
it right because orient contains in fact the inverse matrix, and
therefore calling .inverse() gives back the correct matrix?

I have noticed that there are errors elsewhere in the code. For instance

in the new code this below contains the Cartesian coordinates of 4
atoms in TIP4P water in the principal frame of reference

int_coords at 0 "( -0.0413178, 0.0470361, -0.0195315 )"
int_coords at 1 "( 0.752151, 0.145409, 0.50674 )"
int_coords at 2 "( -0.0963127, -0.892013, -0.196717 )"
int_coords at 3 "( 0.037434, -0.0426152, 0.0176959 )"

The coordinates are wrong (you can see the z coordinate is not zero).

Whereas in the old code I would get

int_coords at 0 "( 3.96922e-08, -0.0655822, -3.83374e-16 )"
int_coords at 1 "( 0.7568, 0.520494, -4.08094e-16 )"
int_coords at 2 "( -0.756801, 0.520494, -3.33067e-16 )"
int_coords at 3 "( -1.41606e-07, 0.059418, -3.26283e-07 )"

where the coordinates are correct (atom O is oxygen, 1,2 H1/H2 and 3 EPW).

That int_coords is wrong is probably because in rbworkspacejm the function

static QVector<Vector> buildBead(const ViewsOfMol &mol,
const AtomCoords &coords,
const AtomProperty<T> &masses,
const QVector<qint32> &beading,
Vector *bead_coords, Matrix *beads_to_world,
Quaternion *bead_orients, double *bead_masses,
Vector *bead_inertias)

contains lines of code such as

Matrix inv_matrix = bead_to_world.inverse();

for (int i=0; i<nats; ++i)
{
int_coords_array[i] = inv_matrix * (coords_array[i] - bead_com);
}

And here I should write instead

for (int i=0; i<nats; ++i)
{
int_coords_array[i] = bead_to_world * (coords_array[i]
- bead_com);
}

I just want to make sure I'm completely clear on the problem since if
I understand well, I need to go everywhere in rbworkspacejm (and
possibly elsewhere) to fix every matrix.inverse() line of code.

Glad to read you are progressing with updating the python build. This
will make our lives much easier !

Best wishes,

Julien

>>> email to sire-developers+unsubscribe@googlegroups.com.
>>> To post to this group, send email to sire-developers@googlegroups.com.

Christopher Woods

unread,

Oct 23, 2014, 9:09:44 AM10/23/14

to Sire Developers

Hi Julien,

I think that it goes a little deeper than that. The problem is that
'orient' and 'bead_to_world' may have been created incorrectly to
begin with. I can't remember the maths off of the top of my head, but
it may be mathematically correct to use

bead_force = orient.inverse() * bead_force;

and it is just that, in the old code, 'orient' was created and stored
column-major, and so it was ok. In the new code, 'orient' may still be
created column-major, but is being stored row-major, and so the code
is broken. We have to be careful with syntax - the problem is that the
matrices in the old code were specified row-major but stored
column-major, while in the new code the matrices are specified
row-major and stored row-major. This means that the new matrix is the
transpose of the old matrix. The .inverse() function returns the
matrix that is the inverse of the matrix, which just so happens to be
the transpose for rotation matrices (e.g. like 'orient' and
'bead_to_world'). I don't think that you should be deleting all of the
'.inverse()' code from the workspace. Rather you should find out how
'orient' and 'bead_to_world' are being constructed, and if they are
being specified as column-major then change the code so that they are
specified row-major. Alternatively. when you create 'orient' and
'bead_to_world', just add a line that reads;

// correct row vs. column major bug introduced in Aug14 version of sire
orient = orient.transpose();

// correct row vs. column major bug introduced in Aug14 version of sire
bead_to_world = bead_to_world.transpose();

and then I can go through the code, find these lines, and work out
mathematically what should be correct ;-)

Cheers,

Christopher

>> >>> email to sire-develope...@googlegroups.com.
>> >>> To post to this group, send email to sire-de...@googlegroups.com.

>> >>> Visit this group at http://groups.google.com/group/sire-developers.
>> >>> For more options, visit https://groups.google.com/d/optout.
>> >>
>> >>
>> >>
>> >> --
>> >> ---------------------------------------------------------
>> >> Christopher Woods
>> >> +44 (0) 7786 264562
>> >> http://chryswoods.com
>> >
>> >
>> >
>> > --
>> > ---------------------------------------------------------
>> > Christopher Woods
>> > +44 (0) 7786 264562
>> > http://chryswoods.com
>>
>>
>>
>> --
>> ---------------------------------------------------------
>> Christopher Woods
>> +44 (0) 7786 264562
>> http://chryswoods.com
>

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

> email to sire-develope...@googlegroups.com.
> To post to this group, send email to sire-de...@googlegroups.com.

Reply all

Reply to author

Forward