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
Finding the eigen values/matrix exponential for higher dimensional data