Re: r24192 - /branches/disp_spin_speed/lib/linear_algebra/matrix_exponential.py -- June 23, 2014 - 11:09

Hi Troels,

There is a faster way to do this - eliminating the 'if' statements.  I
don't know how much faster though, as it is clearly currently not a
bottleneck, but it might cause the code to be much faster later when
other changes to the numeric models occur.  What would be better is to
duplicate and rename the function into the matrix_exponential_rank6()
and matrix_exponential_rank7() functions.  It is more code, but it
will be faster as it will be less lines per function.  If you do this,
for stability it would be good to first have unit tests for the rank 6
and 7 structures.

Regards,

Edward



On 20 June 2014 09:38,  <tlinnet@xxxxxxxxxxxxx> wrote:
> Author: tlinnet
> Date: Fri Jun 20 09:38:33 2014
> New Revision: 24192
>
> URL: http://svn.gna.org/viewcvs/relax?rev=24192&view=rev
> Log:
> Made the function matrix_exponential_rankN also find the exponential if the 
> experiments indices are missing.
>
> Task #7807 (https://gna.org/task/index.php?7807): Speed-up of dispersion 
> models for Clustered analysis.
>
> Modified:
>     branches/disp_spin_speed/lib/linear_algebra/matrix_exponential.py
>
> Modified: branches/disp_spin_speed/lib/linear_algebra/matrix_exponential.py
> URL: 
> http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/lib/linear_algebra/matrix_exponential.py?rev=24192&r1=24191&r2=24192&view=diff
> ==============================================================================
> --- branches/disp_spin_speed/lib/linear_algebra/matrix_exponential.py   
> (original)
> +++ branches/disp_spin_speed/lib/linear_algebra/matrix_exponential.py   Fri 
> Jun 20 09:38:33 2014
> @@ -68,7 +68,13 @@
>      @rtype:     numpy float array of rank [NE][NS][NM][NO][ND][X][X]
>      """
>
> -    NE, NS, NM, NO, ND, Row, Col = A.shape
> +    # Set initial to None.
> +    NE, NS, NM, NO, ND, Row, Col = None, None, None, None, None, None, None
> +
> +    if len(A.shape) == 7:
> +        NE, NS, NM, NO, ND, Row, Col = A.shape
> +    elif len(A.shape) == 6:
> +        NS, NM, NO, ND, Row, Col = A.shape
>
>      # Is the original matrix real?
>      complex_flag = any(iscomplex(A))
> @@ -84,10 +90,16 @@
>
>      # Calculate the exponential of all elements in the input array. Shape 
> [NE][NS][NM][NO][ND][X]
>      # Add one axis, to allow for broadcasting multiplication.
> -    W_exp = exp(W).reshape(NE, NS, NM, NO, ND, Row, 1)
> +    if NE == None:
> +        W_exp = exp(W).reshape(NS, NM, NO, ND, Row, 1)
> +    else:
> +        W_exp = exp(W).reshape(NE, NS, NM, NO, ND, Row, 1)
>
>      # Make a eye matrix, with Shape [NE][NS][NM][NO][ND][X][X]
> -    eye_mat = tile(eye(Row)[newaxis, newaxis, newaxis, newaxis, newaxis, 
> ...], (NE, NS, NM, NO, ND, 1, 1) )
> +    if NE == None:
> +        eye_mat = tile(eye(Row)[newaxis, newaxis, newaxis, newaxis, 
> newaxis, ...], (NS, NM, NO, ND, 1, 1) )
> +    else:
> +        eye_mat = tile(eye(Row)[newaxis, newaxis, newaxis, newaxis, 
> newaxis, ...], (NE, NS, NM, NO, ND, 1, 1) )
>
>      # Transform it to a diagonal matrix, with elements from vector down 
> the diagonal.
>      W_exp_diag = multiply(W_exp, eye_mat )
>
>
> _______________________________________________
> relax (http://www.nmr-relax.com)
>
> This is the relax-commits mailing list
> relax-commits@xxxxxxx
>
> 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