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 <[email protected]> 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 <[email protected]> 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 <[email protected]>:
>>
>>> 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 <[email protected]>
>>> 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 <[email protected]>:
>>> >
>>> >> 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 <[email protected]>
>>> >> 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 <[email protected]>:
>>> >> >>
>>> >> >> 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, <[email protected]> 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
>>> >> >> > [email protected]
>>> >> >> >
>>> >> >> > 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
>>> >> >> [email protected]
>>> >> >>
>>> >> >>
>>> >> >> 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
>>> >> >
>>> >> >
>>> >
>>> >
>>
>>
_______________________________________________
relax (http://www.nmr-relax.com)
This is the relax-devel mailing list
[email protected]
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
