Wow, this is incredible! :S Sorry, I never realised how buggy the
out argument for numpy.dot() was! It looks like numpy only has
functional inplace operations for the ufuncs. It somehow worked when
I used it, as I'd often use the R^T.A.R rotation form which somehow
does not suffer from this bad behaviour (as in the test script I
wrote). Maybe because of the transpose being used which is the
Fortran ordering?!? I'll now have to carefully check all parts of
relax that use this :S
I've now looked all over the internet, and it's clear that the numpy
people either don't understand the concept or don't care. There are
so many posts where no real answer has been given. It really reminds
me of when I found all the fatal bugs in all SciPy optimisation
algorithms, which then lead me to write minfx. For now, I cannot come
up with an answer. For the future I might write my own simple FORTRAN
BLAS and LAPACK library interface directly into relax - that was on my
list of things to do for the frame order analysis to remove the SciPy
dependency and to have better access to these powerful libraries, and
to use the super fast ALTAS. So the best might be to revert to the "c
= dot(a, b)" notation:
"""
Index: lib/dispersion/ns_cpmg_2site_3d.py
===================================================================
--- lib/dispersion/ns_cpmg_2site_3d.py (revision 23967)
+++ lib/dispersion/ns_cpmg_2site_3d.py (working copy)
@@ -141,9 +141,6 @@
# The matrix R that contains all the contributions to the
evolution, i.e. relaxation, exchange and chemical shift evolution.
R = rcpmg_3d(R1A=r10a, R1B=r10b, R2A=R2A_si_mi,
R2B=R2B_si_mi, pA=pA, pB=pB, dw=dw_si_mi, k_AB=k_AB, k_BA=k_BA)
- # The essential evolution matrix. This initialises the structure.
- evolution_matrix = asarray(R) * 0.0
-
# Loop over the time points, back calculating the R2eff values.
for di in range(num_points_si_mi):
# Extract the values from the higher dimensional arrays.
@@ -160,12 +157,12 @@
# The essential evolution matrix.
# This is the first round.
- dot(Rexpo, r180x, evolution_matrix)
- dot(evolution_matrix * 1.0, Rexpo, evolution_matrix)
+ evolution_matrix = dot(Rexpo, r180x)
+ evolution_matrix = dot(evolution_matrix, Rexpo)
# The second round.
- dot(evolution_matrix * 1.0, evolution_matrix * 1.0,
evolution_matrix)
+ evolution_matrix = dot(evolution_matrix, evolution_matrix)
# The third and fourth round,
- dot(evolution_matrix * 1.0, evolution_matrix * 1.0,
evolution_matrix)
+ evolution_matrix = dot(evolution_matrix, evolution_matrix)
# Loop over the CPMG elements, propagating the magnetisation.
for j in range(power_si_mi_di/2):
"""
This still has the number of dot products reduced to the minimal
amount and is faster than trunk. And it is cleaner, though maybe a
little slower, than creating 4 different 'evolution_matrix_*'
structures at the start, though the following appears not so bad:
"""
Index: lib/dispersion/ns_cpmg_2site_3d.py
===================================================================
--- lib/dispersion/ns_cpmg_2site_3d.py (revision 23967)
+++ lib/dispersion/ns_cpmg_2site_3d.py (working copy)
@@ -54,7 +54,7 @@
"""
# Python module imports.
-from numpy import asarray, dot, fabs, isfinite, log, min, sum
+from numpy import asarray, dot, fabs, float64, isfinite, log, min, sum, zeros
from numpy.ma import fix_invalid, masked_where
@@ -124,6 +124,9 @@
M0[1] = pA
M0[4] = pB
+ # Evolution matrix initialisation.
+ evolution_matrix = zeros((2, 7, 7), float64)
+
# Extract the total numbers of experiments, number of spins,
number of magnetic field strength, number of offsets, maximum number
of dispersion point.
NE, NS, NM, NO, ND = back_calc.shape
@@ -141,9 +144,6 @@
# The matrix R that contains all the contributions to the
evolution, i.e. relaxation, exchange and chemical shift evolution.
R = rcpmg_3d(R1A=r10a, R1B=r10b, R2A=R2A_si_mi,
R2B=R2B_si_mi, pA=pA, pB=pB, dw=dw_si_mi, k_AB=k_AB, k_BA=k_BA)
- # The essential evolution matrix. This initialises the structure.
- evolution_matrix = asarray(R) * 0.0
-
# Loop over the time points, back calculating the R2eff values.
for di in range(num_points_si_mi):
# Extract the values from the higher dimensional arrays.
@@ -160,16 +160,16 @@
# The essential evolution matrix.
# This is the first round.
- dot(Rexpo, r180x, evolution_matrix)
- dot(evolution_matrix * 1.0, Rexpo, evolution_matrix)
+ dot(Rexpo, r180x, evolution_matrix[0])
+ dot(evolution_matrix[0], Rexpo, evolution_matrix[1])
# The second round.
- dot(evolution_matrix * 1.0, evolution_matrix * 1.0,
evolution_matrix)
+ dot(evolution_matrix[1], evolution_matrix[1],
evolution_matrix[0])
# The third and fourth round,
- dot(evolution_matrix * 1.0, evolution_matrix * 1.0,
evolution_matrix)
+ dot(evolution_matrix[0], evolution_matrix[0],
evolution_matrix[1])
# Loop over the CPMG elements, propagating the magnetisation.
for j in range(power_si_mi_di/2):
- Mint = dot(evolution_matrix, Mint)
+ Mint = dot(evolution_matrix[1], Mint)
# The next lines calculate the R2eff using a
two-point approximation, i.e. assuming that the decay is
mono-exponential.
Mx = Mint[1] / pA
"""
This just jumps back and forth between evolution_matrix[0] and
evolution_matrix[1] to avoid this horrible bug. On a different note,
it would be nice to have
"Rexpo.r180x.Rexpo.Rexpo.r180x.Rexpo.Rexpo.r180x.Rexpo.Rexpo.r180x.Rexpo"
added to the comment, to explain what is being calculated. Sorry that
I cannot find a solution to this clear numpy bug!
Regards,
Edward
On 15 June 2014 17:12, Edward d'Auvergne <edward@xxxxxxxxxxxxx> wrote:
Hmmm, it sounds like a documentation bug. It looks like the Fortran
memory ordering is required for the out argument to work when the a or
b arguments are the same structure:
http://thread.gmane.org/gmane.comp.python.numeric.general/28135/focus=28826
That's from the same thread. From this thread
http://mail.scipy.org/pipermail/numpy-discussion/2013-May/066668.html,
it looks like the NumPy people, or those who responded, do not
understand this concept. I'll keep searching.
Cheers,
Edward
On 15 June 2014 17:00, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx> wrote:
There is long thread here
http://thread.gmane.org/gmane.comp.python.numeric.general/28135/
Somewhere I saw, that there is maybe a bug here.
2014-06-15 16:53 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:
Hi,
It does look like there is a bug somewhere in one version of
numpy.dot(). I'll look into this!
Regards,
Edward
On 15 June 2014 16:43, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx>
wrote:
I am using 1.8
I am stopping here with more search for this holy grail.
There is clearly some bug somewhere.
2014-06-15 16:25 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:
Which numpy version are you using? There is a damn good reason why I
am pushing this. I know that all of the competition software (at
least cpmg_fit and CATIA) uses exactly this storage concept, though in
C rather than Python, and for exactly the reason why I am pushing
this. We really have to get this working for you!
Regards,
Edward
On 15 June 2014 16:21, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx>
wrote:
Hi Ed.
I can see that you are pushing for this.
It wont work!
2014-06-15 16:12 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:
Hi Troels,
In this change, only in one line have you used the numpy.dot() out
argument:
+ dot(Rexpo, r180x, evolution_matrix)
So it must be working. You can do this for all others, replacing:
- evolution_matrix = dot(evolution_matrix, Rexpo)
+ dot(evolution_matrix, Rexpo, evolution_matrix)
and:
- evolution_matrix = dot(evolution_matrix,
evolution_matrix)
+ dot(evolution_matrix, evolution_matrix,
evolution_matrix)
and:
- Mint = evolution_matrix.dot(Mint)
+ dot(evolution_matrix, Mint, Mint)
This should give an observable speed increase of 2 or more.
Regards,
Edward
On 15 June 2014 15:15, <tlinnet@xxxxxxxxxxxxx> wrote:
Author: tlinnet
Date: Sun Jun 15 15:15:23 2014
New Revision: 23963
URL: http://svn.gna.org/viewcvs/relax?rev=23963&view=rev
Log:
Made the dot evolution structure faster for NS CPMG 2site 3D.
Task #7807 (https://gna.org/task/index.php?7807): Speed-up of
dispersion
models for Clustered analysis.
Modified:
branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py
Modified:
branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py
URL:
http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py?rev=23963&r1=23962&r2=23963&view=diff
==============================================================================
--- branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py
(original)
+++ branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py
Sun
Jun
15 15:15:23 2014
@@ -54,7 +54,7 @@
"""
# Python module imports.
-from numpy import dot, fabs, isfinite, log, min, sum
+from numpy import asarray, dot, fabs, isfinite, log, min, sum
from numpy.ma import fix_invalid, masked_where
@@ -141,6 +141,9 @@
# The matrix R that contains all the contributions to
the
evolution, i.e. relaxation, exchange and chemical shift evolution.
R = rcpmg_3d(R1A=r10a, R1B=r10b, R2A=R2A_si_mi,
R2B=R2B_si_mi, pA=pA, pB=pB, dw=dw_si_mi, k_AB=k_AB, k_BA=k_BA)
+ # The essential evolution matrix. This initialises
the
structure.
+ evolution_matrix = asarray(R) * 0.0
+
# Loop over the time points, back calculating the
R2eff
values.
for di in range(num_points_si_mi):
# Extract the values from the higher dimensional
arrays.
@@ -155,12 +158,18 @@
# This matrix is a propagator that will evolve
the
magnetization with the matrix R for a delay tcp.
Rexpo = matrix_exponential(R*tcp_si_mi_di)
- # Temp matrix.
- t_mat =
Rexpo.dot(r180x).dot(Rexpo).dot(Rexpo).dot(r180x).dot(Rexpo).dot(Rexpo).dot(r180x).dot(Rexpo).dot(Rexpo).dot(r180x).dot(Rexpo)
+ # The essential evolution matrix.
+ # This is the first round.
+ dot(Rexpo, r180x, evolution_matrix)
+ evolution_matrix = dot(evolution_matrix, Rexpo)
+ # The second round.
+ evolution_matrix = dot(evolution_matrix,
evolution_matrix)
+ # The third and fourth round,
+ evolution_matrix = dot(evolution_matrix,
evolution_matrix)
# Loop over the CPMG elements, propagating the
magnetisation.
for j in range(power_si_mi_di/2):
- Mint = t_mat.dot(Mint)
+ Mint = evolution_matrix.dot(Mint)
# The next lines calculate the R2eff using a
two-point
approximation, i.e. assuming that the decay is mono-exponential.
Mx = Mint[1] / pA
_______________________________________________
relax (http://www.nmr-relax.com)
This is the relax-commits mailing list
relax-commits@xxxxxxx
To unsubscribe from this list, get a password
reminder, or change your subscription options,
visit the list information page at
https://mail.gna.org/listinfo/relax-commits
_______________________________________________
relax (http://www.nmr-relax.com)
This is the relax-devel mailing list
relax-devel@xxxxxxx
To unsubscribe from this list, get a password
reminder, or change your subscription options,
visit the list information page at
https://mail.gna.org/listinfo/relax-devel
Re: r23963 - /branches/disp_spin_speed/lib/dispersion/ns_cpmg_2site_3d.py