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
