Well, looking at the Scipy source code of ./scipy/linalg/matfuncs.py, expm2() and expm3() are implemented in Python code, interfacing to other scipy functions. And expm() is an alias for the function of the same name in ./scipy/sparse/linalg/matfuncs.py, which is also implemented in Python code. So there is no BLAS/LAPACK magic happening here, except maybe in the functions that these expm*() functions call.
Regards, Edward On 25 June 2014 17:54, Edward d'Auvergne <[email protected]> wrote: > Hi Troels, > > If you are interested in eliminating this 86% bottleneck in the > numeric R1rho dispersion models in the > lib.dispersion.matrix_exponential.matrix_exponential_rank_NE_NS_NM_NO_ND_x_x() > function, it might be worth looking at the scipy.linalg.expm(), > scipy.linalg.expm2(), and scipy.linalg.expm3() functions > (http://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#exponential-and-logarithm-functions), > or the BLAS/LAPACK functions these interface to. Note that I > eliminated the use of these functions in the dispersion analysis a > while ago (http://article.gmane.org/gmane.science.nmr.relax.scm/18869/). > And for good reason, there was a nasty Scipy bug in the Pade > approximation of scipy.linalg.expm(). This and the fatal bugs in all > optimisation algorithms in Scipy (the catalyst for me to create minfx, > originally as part of relax, https://gna.org/projects/minfx/), is the > reason why I try to avoid Scipy as much as possible - I just cannot > trust it any more. Scipy is a collection of poorly integrated tools, > often with the original authors completely disappearing, and it is > full of bugs. Anyway, maybe one of these functions would be of use or > would be a good inspiration. Oh, the scipy.linalg.expm() function in > the original code comes from the fitting_main_kex.py script from > Mathilde Lescanne, Paul Schanda, and Dominique Marion. Note that if > the Taylor series approximation of scipy.linalg.expm3() works well for > the numeric dispersion models, which is untested, that might be by far > the fastest solution. > > Regards, > > Edward > > > On 25 June 2014 08:20, Edward d'Auvergne <[email protected]> wrote: >> Wow, you managed to work out how to remove the slow Python looping >> from the numeric models! For the record, this is the speed up for >> this 'NS R1rho 2-site' model using Python 2.7.6 & numpy 1.8.0 on a >> 64-bit Linux system: >> >> """ >> $ ./disp_profile_all.py /data/relax/branches/disp_spin_speed2 >> /data/relax/branches/disp_spin_speed >> [snip] >> >> New relax version: relax repository checkout r24308 >> svn+ssh://[email protected]/svn/relax/branches/disp_spin_speed >> Old relax version: relax repository checkout r24274 >> svn+ssh://[email protected]/svn/relax/branches/disp_spin_speed >> >> Execution iteration 1 >> >> $ python profiling_ns_r1rho_2site.py /data/relax/branches/disp_spin_speed2 >> 1000 0.016 0.000 3.561 0.004 >> relax_disp.py:1582(func_ns_r1rho_2site) >> 10 0.013 0.001 3.237 0.324 >> relax_disp.py:1582(func_ns_r1rho_2site) >> $ python profiling_ns_r1rho_2site.py /data/relax/branches/disp_spin_speed >> 1000 0.017 0.000 6.245 0.006 >> relax_disp.py:1552(func_ns_r1rho_2site) >> 10 0.013 0.001 5.899 0.590 >> relax_disp.py:1552(func_ns_r1rho_2site) >> [snip] >> >> 100 single spins analysis: >> NS R1rho 2-site: 62.450+/-0.000 -> 35.610+/-0.000, 1.754x >> faster. >> >> Cluster of 100 spins analysis: >> NS R1rho 2-site: 58.990+/-0.000 -> 32.370+/-0.000, 1.822x >> faster. >> """ >> >> >> This is quite significant for these models. With all the changes >> since relax 3.2.2: >> >> """ >> $ ./disp_profile_all.py /data/relax/branches/disp_spin_speed2 >> /data/relax/tags/3.2.2/ >> [snip] >> >> New relax version: relax repository checkout r24308 >> svn+ssh://[email protected]/svn/relax/branches/disp_spin_speed >> Old relax version: relax 3.2.2 >> >> Execution iteration 1 >> >> $ python profiling_ns_r1rho_2site.py /data/relax/branches/disp_spin_speed2 >> 1000 0.016 0.000 3.576 0.004 >> relax_disp.py:1582(func_ns_r1rho_2site) >> 10 0.013 0.001 3.247 0.325 >> relax_disp.py:1582(func_ns_r1rho_2site) >> $ python profiling_ns_r1rho_2site.py /data/relax/tags/3.2.2/ >> 1000 0.245 0.000 24.323 0.024 >> relax_disp.py:1681(func_ns_r1rho_2site) >> 10 0.266 0.027 24.439 2.444 >> relax_disp.py:1681(func_ns_r1rho_2site) >> [snip] >> >> 100 single spins analysis: >> NS R1rho 2-site: 243.230+/-0.000 -> 35.760+/-0.000, 6.802x >> faster. >> >> Cluster of 100 spins analysis: >> NS R1rho 2-site: 244.390+/-0.000 -> 32.470+/-0.000, 7.527x >> faster. >> """ >> >> >> Considering how slow these models are compared to the analytic models, >> this is a huge improvement! According to your >> profiling_ns_r1rho_2site.py script, the slow parts are: >> >> """ >> ('chi2 cluster:', 14153134802.76194) >> Wed Jun 25 07:59:18 2014 /tmp/tmp17lqUn >> >> 211077 function calls in 3.510 seconds >> >> Ordered by: cumulative time >> >> ncalls tottime percall cumtime percall filename:lineno(function) >> 1 0.000 0.000 3.510 3.510 <string>:1(<module>) >> 1 0.002 0.002 3.510 3.510 >> profiling_ns_r1rho_2site.py:553(cluster) >> 10 0.000 0.000 3.254 0.325 >> profiling_ns_r1rho_2site.py:515(calc) >> 10 0.015 0.002 3.254 0.325 >> relax_disp.py:1582(func_ns_r1rho_2site) >> 10 0.015 0.002 3.236 0.324 >> ns_r1rho_2site.py:190(ns_r1rho_2site) >> 10 0.104 0.010 3.041 0.304 >> matrix_exponential.py:33(matrix_exponential_rank_NE_NS_NM_NO_ND_x_x) >> 10 1.803 0.180 1.827 0.183 linalg.py:982(eig) >> 10 0.564 0.056 0.601 0.060 linalg.py:455(inv) >> 40 0.501 0.013 0.501 0.013 {numpy.core.multiarray.einsum} >> """ >> >> >> The single spin numbers are very similar. Looking at the total times, >> the numpy.einsum and numpy.inv functions are taking 1/7 of the time >> each, and numpy.eig 51% of the time - this is a total of 80% of the >> time. The >> lib.dispersion.matrix_exponential.matrix_exponential_rank_NE_NS_NM_NO_ND_x_x() >> function now takes up 86% of the cumulative time - so this is the new >> bottleneck. The only way to solve this would be using clever linear >> algebra tricks to simplify or to change the logic in >> matrix_exponential_rank_NE_NS_NM_NO_ND_x_x(). There are multiple >> techniques for finding the matrix exponential, it could be that the >> eigenvalue method is one of the slowest. This looks to be interesting >> for an overview of recent research in the mathematics field on this >> subject: http://www.maths.manchester.ac.uk/~higham/talks/exp09.pdf. >> Great work though, you have eliminated the worst bottleneck for the >> numeric models! >> >> Cheers, >> >> Edward >> >> >> >> On 25 June 2014 02:14, <[email protected]> wrote: >>> Author: tlinnet >>> Date: Wed Jun 25 02:14:36 2014 >>> New Revision: 24294 >>> >>> URL: http://svn.gna.org/viewcvs/relax?rev=24294&view=rev >>> Log: >>> Speeded up the code of NS r1rho 2site. >>> >>> This was essential done to numpy einsum, and doing the dot operations in >>> multiple dimensions. >>> It was though necessary to realize, that to do the proper dot product >>> operations, the outer two >>> axis if M0 should be swapped, by rolling the outer axis one back. >>> >>> 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_r1rho_2site.py >>> branches/disp_spin_speed/target_functions/relax_disp.py >>> >>> Modified: branches/disp_spin_speed/lib/dispersion/ns_r1rho_2site.py >>> URL: >>> http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/dispersion/ns_r1rho_2site.py?rev=24294&r1=24293&r2=24294&view=diff >>> ============================================================================== >>> --- branches/disp_spin_speed/lib/dispersion/ns_r1rho_2site.py (original) >>> +++ branches/disp_spin_speed/lib/dispersion/ns_r1rho_2site.py Wed Jun 25 >>> 02:14:36 2014 >>> @@ -50,10 +50,10 @@ >>> """ >>> >>> # Python module imports. >>> -from numpy import array, cos, dot, einsum, float64, log, multiply, >>> rollaxis, sin, sum >>> +from numpy import array, einsum, float64, isfinite, log, min, multiply, >>> rollaxis, sin, sum >>> +from numpy.ma import fix_invalid, masked_less >>> >>> # relax module imports. >>> -from lib.float import isNaN >>> from lib.dispersion.matrix_exponential import >>> matrix_exponential_rank_NE_NS_NM_NO_ND_x_x >>> >>> # Repetitive calculations (to speed up calculations). >>> @@ -187,14 +187,16 @@ >>> return matrix >>> >>> >>> -def ns_r1rho_2site(M0=None, r1rho_prime=None, omega=None, offset=None, >>> r1=0.0, pA=None, dw=None, kex=None, spin_lock_fields=None, relax_time=None, >>> inv_relax_time=None, back_calc=None, num_points=None): >>> +def ns_r1rho_2site(M0=None, M0_T=None, r1rho_prime=None, omega=None, >>> offset=None, r1=0.0, pA=None, dw=None, kex=None, spin_lock_fields=None, >>> relax_time=None, inv_relax_time=None, back_calc=None, num_points=None): >>> """The 2-site numerical solution to the Bloch-McConnell equation for >>> R1rho data. >>> >>> This function calculates and stores the R1rho values. >>> >>> >>> @keyword M0: This is a vector that contains the initial >>> magnetizations corresponding to the A and B state transverse magnetizations. >>> - @type M0: numpy float64, rank-1, 7D array >>> + @type M0: numpy float array of rank >>> [NE][NS][NM][NO][ND][6][1] >>> + @keyword M0_T: This is a vector that contains the initial >>> magnetizations corresponding to the A and B state transverse >>> magnetizations, where the outer two axis has been swapped for efficient dot >>> operations. >>> + @type M0_T: numpy float array of rank >>> [NE][NS][NM][NO][ND][1][6] >>> @keyword r1rho_prime: The R1rho_prime parameter value (R1rho >>> with no exchange). >>> @type r1rho_prime: numpy float array of rank >>> [NE][NS][NM][NO][ND] >>> @keyword omega: The chemical shift for the spin in rad/s. >>> @@ -238,43 +240,22 @@ >>> # Magnetization evolution. >>> Rexpo_M0_mat = einsum('...ij,...jk', Rexpo_mat, M0) >>> >>> - # Transpose M0, to prepare for dot operation. Roll the last axis one >>> back. >>> - M0_T = rollaxis(M0, 6, 5) >>> - >>> - if NS == 1: >>> - # Loop over the spectrometer frequencies. >>> - for mi in range(NM): >>> - # Loop over offsets: >>> - for oi in range(NO): >>> - # Extract number of points. >>> - num_points_i = num_points[0, 0, mi, oi] >>> - >>> - # Loop over the time points, back calculating the R2eff >>> values. >>> - for j in range(num_points_i): >>> - >>> - # If the number of spins are 1, do the fastest method >>> by dot product. Else do it by einstein summation. >>> - # Set the spin index for one spin. >>> - si = 0 >>> - # Extract the preformed matrix that rotate the >>> magnetization previous to spin-lock into the weff frame. >>> - M0_i= M0_T[0, si, mi, oi, j] >>> - >>> - # Extract from the pre-formed Magnetization evolution >>> matrix. >>> - Rexpo_M0_mat_i = Rexpo_M0_mat[0, si, mi, oi, j] >>> - >>> - # Magnetization evolution. >>> - MA = dot(M0_i, Rexpo_M0_mat_i)[0, 0] >>> - >>> - # The next lines calculate the R1rho using a two-point >>> approximation, i.e. assuming that the decay is mono-exponential. >>> - if MA <= 0.0 or isNaN(MA): >>> - back_calc[0, si, mi, oi, j] = 1e99 >>> - else: >>> - back_calc[0, si, mi, oi, j]= -inv_relax_time[0, >>> si, mi, oi, j] * log(MA) >>> - >>> - else: >>> - # Magnetization evolution, which include all dimensions. >>> - MA_mat = einsum('...ij,...jk', M0_T, Rexpo_M0_mat)[:, :, :, :, :, >>> 0, 0] >>> - >>> - # Do back calculation. >>> - back_calc[:] = -inv_relax_time * log(MA_mat) >>> - >>> - >>> + # Magnetization evolution, which include all dimensions. >>> + MA_mat = einsum('...ij,...jk', M0_T, Rexpo_M0_mat)[:, :, :, :, :, 0, 0] >>> + >>> + # Insert safe checks. >>> + if min(MA_mat) < 0.0: >>> + mask_min_MA_mat = masked_less(MA_mat, 0.0) >>> + # Fill with high values. >>> + MA_mat[mask_min_MA_mat.mask] = 1e100 >>> + >>> + # Do back calculation. >>> + back_calc[:] = -inv_relax_time * log(MA_mat) >>> + >>> + # Catch errors, taking a sum over array is the fastest way to check for >>> + # +/- inf (infinity) and nan (not a number). >>> + if not isfinite(sum(back_calc)): >>> + # Replaces nan, inf, etc. with fill value. >>> + fix_invalid(back_calc, copy=False, fill_value=1e100) >>> + >>> + >>> >>> Modified: branches/disp_spin_speed/target_functions/relax_disp.py >>> URL: >>> http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/target_functions/relax_disp.py?rev=24294&r1=24293&r2=24294&view=diff >>> ============================================================================== >>> --- branches/disp_spin_speed/target_functions/relax_disp.py (original) >>> +++ branches/disp_spin_speed/target_functions/relax_disp.py Wed Jun 25 >>> 02:14:36 2014 >>> @@ -26,8 +26,7 @@ >>> >>> # Python module imports. >>> from copy import deepcopy >>> -from math import pi >>> -from numpy import add, array, arctan2, asarray, cos, complex64, dot, >>> float64, int16, max, multiply, ones, sin, sqrt, sum, tile, zeros >>> +from numpy import arctan2, cos, dot, float64, int16, multiply, ones, >>> rollaxis, pi, sin, sum, zeros >>> from numpy.ma import masked_equal >>> >>> # relax module imports. >>> @@ -415,6 +414,9 @@ >>> M0_2[2, 0] = 1 >>> M0_cos = multiply.outer( cos(theta_mat), M0_2 ) >>> self.M0 = M0_sin + M0_cos >>> + # Transpose M0, to prepare for dot operation. Roll the last >>> axis one back, corresponds to a transpose for the outer two axis. >>> + self.M0_T = rollaxis(self.M0, 6, 5) >>> + >>> if model in [MODEL_NS_R1RHO_3SITE, MODEL_NS_R1RHO_3SITE_LINEAR]: >>> self.M0 = zeros(9, float64) >>> # Offset of spin-lock from A. >>> @@ -1598,7 +1600,7 @@ >>> self.r20_struct[:] = multiply.outer( r1rho_prime.reshape(self.NE, >>> self.NS, self.NM), self.no_nd_ones ) >>> >>> # Back calculate the R2eff values. >>> - ns_r1rho_2site(M0=self.M0, r1rho_prime=self.r20_struct, >>> omega=self.chemical_shifts, offset=self.offset, r1=self.r1, pA=pA, >>> dw=self.dw_struct, kex=kex, spin_lock_fields=self.spin_lock_omega1, >>> relax_time=self.relax_times, inv_relax_time=self.inv_relax_times, >>> back_calc=self.back_calc, num_points=self.num_disp_points) >>> + ns_r1rho_2site(M0=self.M0, M0_T=self.M0_T, >>> r1rho_prime=self.r20_struct, omega=self.chemical_shifts, >>> offset=self.offset, r1=self.r1, pA=pA, dw=self.dw_struct, kex=kex, >>> spin_lock_fields=self.spin_lock_omega1, relax_time=self.relax_times, >>> inv_relax_time=self.inv_relax_times, back_calc=self.back_calc, >>> num_points=self.num_disp_points) >>> >>> # Clean the data for all values, which is left over at the end of >>> arrays. >>> self.back_calc = self.back_calc*self.disp_struct >>> >>> >>> _______________________________________________ >>> 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
