Hi, You are aware of the lib.linear_algebra.matrix_exponential and lib.linear_algebra.matrix_power modules right?
Regards, Edward On 18 June 2014 13:33, Edward d'Auvergne <[email protected]> wrote: > Hi, > > If you would like to tackle this one, that would be legendary. Though > you must decide if it is worth your time. There are also still a few > things to tidy up for the disp_spin_speed branch to prepare it for > merger. For your idea, I would suggest considering the following: > > 1) Creating a simple timeit script such as in > http://svn.gna.org/viewcvs/relax/branches/frame_order_cleanup/test_suite/shared_data/frame_order/timings/ > to show the timing of one function where you loop over NS and NM, > calculating the matrix exponential each time, and one where you do it > all in one operation. You can get initial matrices by running one of > the system tests in debugging mode and adding a print statement in an > appropriate place, then replicate this matrix over NS and NM. This > will allow you to debug the multi-matrix exponential calculation. > > 2) Before making changes in lib.dispersion, make sure that the system > and unit test coverage is reasonable for each model to be speed up. > > 3) Create a profiling script in > test_suite/shared_data/dispersion/profiling so we can measure and see > the improvements. If the numpy operation is costly, there is a slight > chance that the change will not speed things up. Point 1) will show > that though. You've probably seen that the system tests are useless > for seeing your huge speed ups. > > If this could be done, I think relax will then be the fastest > dispersion implementation out there. Unfortunately I don't have the > source code to all software, so I don't know if this would be true. > And relax uses minfx and different optimisation to the competition, > avoiding numerical gradients, so there are many other factors > affecting the speed of relax vs. the rest. > > Regards, > > Edward > > > > On 18 June 2014 13:00, Troels Emtekær Linnet <[email protected]> wrote: >> The worst bottleneck for the numerical analysis, is the looping over spins. >> >> For the nmerical models, a propagator that will evolve the magnetization is >> calculated. >> >> This is the matrix exponential of the matrix R that contains all the >> contributions to the evolution, i.e. relaxation, exchange and chemical >> shift evolution. >> The shape of R, can be (7, 7) or (9,9). >> >> Finding the matrix exponential is done over looping of spins, >> and spectrometer frequencies. >> NS and NM. >> >> If NS is 100, and NM is 2, maybe R, could get packed to: >> (100, 2, 7, 7) >> >> The matrix exponential is: >> # The eigenvalue decomposition. >> W, V = eig(A) >> >> # Calculate the exact exponential. >> eA = dot(dot(V, diag(exp(W))), inv(V)) >> >> But according to: >> http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html >> >> This should be possible for an array of matrixes. >> >> >> ##### >> Compute the eigenvalues and right eigenvectors of a square array. >> >> Parameters : >> a : (..., M, M) array >> Matrices for which the eigenvalues and right eigenvectors will be computed >> Returns : >> w : (..., M) array >> The eigenvalues, each repeated according to its multiplicity. The >> eigenvalues are not necessarily ordered. The resulting array will be always >> be of complex type. When a is real the resulting eigenvalues will be real >> (0 imaginary part) or occur in conjugate pairs >> v : (..., M, M) array >> The normalized (unit “length”) eigenvectors, such that the column v[:,i] is >> the eigenvector corresponding to the eigenvalue w[i]. >> Raises : >> LinAlgError >> If the eigenvalue computation does not converge. >> >> ------- >> Broadcasting rules apply, see the numpy.linalg documentation for details. >> >> This is implemented using the _geev LAPACK routines which compute the >> eigenvalues and eigenvectors of general square arrays. >> >> The number w is an eigenvalue of a if there exists a vector v such that >> dot(a,v) = w * v. Thus, the arrays a, w, and v satisfy the equations >> dot(a[:,:], v[:,i]) = w[i] * v[:,i] for i \in \{0,...,M-1\}. >> ------ >> #### >> >> R is: >> 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) >> >> So that would need a new way to expand this, since we have R20 and dw as >> specific for spin and frequency. >> >> >> >> Troels Emtekær Linnet >> _______________________________________________ >> 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
