Seems like I am not the only one having fun with higher dimensions. ;-) is timeit better than profiling?
How about the model_free analysis? Is that in for loops? best Troels 2014-06-13 14:03 GMT+02:00 <[email protected]>: > Author: bugman > Date: Fri Jun 13 14:03:44 2014 > New Revision: 23935 > > URL: http://svn.gna.org/viewcvs/relax?rev=23935&view=rev > Log: > Added a script for timing different ways to calculate PCSs and RDCs for > multiple vectors. > > This uses the timeit module rather than profile to demonstrate the speed > of 7 different ways to > calculate the RDCs or PCSs for an array of vectors using numpy. In the > frame order analysis, this > is the bottleneck for the quasi-random numerical integration of the PCS. > > The log file shows a potential 1 order of magnitude speed up between the > 1st technique, which is > currently used in the frame order analysis, and the 7th and last > technique. The first technique > loops over each vector, calculating the PCS. The last expands the PCS/RDC > equation of the > projection of the vector into the alignment tensor, and calculates all > PCSs simultaneously. > > > Added: > > branches/frame_order_cleanup/test_suite/shared_data/frame_order/timings/ > > branches/frame_order_cleanup/test_suite/shared_data/frame_order/timings/tensor_projections.log > > branches/frame_order_cleanup/test_suite/shared_data/frame_order/timings/tensor_projections.py > > Added: > branches/frame_order_cleanup/test_suite/shared_data/frame_order/timings/tensor_projections.log > URL: > http://svn.gna.org/viewcvs/relax/branches/frame_order_cleanup/test_suite/shared_data/frame_order/timings/tensor_projections.log?rev=23935&view=auto > > ============================================================================== > --- > branches/frame_order_cleanup/test_suite/shared_data/frame_order/timings/tensor_projections.log > (added) > +++ > branches/frame_order_cleanup/test_suite/shared_data/frame_order/timings/tensor_projections.log > Fri Jun 13 14:03:44 2014 > @@ -0,0 +1,38 @@ > +Original vectors: > +[[ 1. 2. 3.] > + [ 2. 2. 2.]] > +Shape: (200, 3) > + > +Tensor: > +[[ 1.42219822 -7.07796212 -6.01619494] > + [-7.07796212 -1.44543002 2.02008007] > + [-6.01619494 2.02008007 0.02323179]] > + > + > +1st projection - element by element r[i].A.r[i]. > +Proj: [-44.31849296 -88.59261589] > +Timing (s): 3.90826106071 > + > +2nd projection - diag of double tensordot. > +Proj: [-44.31849296 -88.59261589] > +Timing (s): 2.68392705917 > + > +3rd projection - diag of double tensordot, no transpose. > +Proj: [-44.31849296 -88.59261589] > +Timing (s): 2.58908486366 > + > +4th projection - mixed tensordot() and per-vector dot(). > +Proj: [-44.31849296 -88.59261589] > +Timing (s): 5.47785592079 > + > +5th projection - expansion and sum. > +Proj: [-44.31849296 -88.59261589] > +Timing (s): 10.7986190319 > + > +6th projection - expansion. > +Proj: [-44.31849296 -88.59261589] > +Timing (s): 0.54491686821 > + > +7th projection - expansion with pre-calculation. > +Proj: [-44.31849296 -88.59261589] > +Timing (s): 0.395635128021 > > Added: > branches/frame_order_cleanup/test_suite/shared_data/frame_order/timings/tensor_projections.py > URL: > http://svn.gna.org/viewcvs/relax/branches/frame_order_cleanup/test_suite/shared_data/frame_order/timings/tensor_projections.py?rev=23935&view=auto > > ============================================================================== > --- > branches/frame_order_cleanup/test_suite/shared_data/frame_order/timings/tensor_projections.py > (added) > +++ > branches/frame_order_cleanup/test_suite/shared_data/frame_order/timings/tensor_projections.py > Fri Jun 13 14:03:44 2014 > @@ -0,0 +1,127 @@ > +# Python module imports. > +from numpy import * > +from numpy.linalg import norm > +from os import pardir, sep > +import sys > +from time import sleep > +from timeit import timeit > + > +# Modify the system path. > +sys.path.append(pardir+sep+pardir+sep+pardir+sep+pardir+sep) > + > +# relax module imports. > +from lib.geometry.rotations import euler_to_R_zyz > + > + > +def proj1(vect, A, N=1, verb=True): > + d = zeros(len(vect), float64) > + for i in range(N): > + for j in xrange(len(vect)): > + d[j] = dot(vect[j], dot(A, vect[j])) > + if verb: > + print("\n1st projection - element by element r[i].A.r[i].") > + print("Proj: %s" % d[:2]) > + > + > +def proj2(vect, A, N=1, verb=True): > + for i in range(N): > + d = diagonal(tensordot(vect, tensordot(A, transpose(vect), > axes=1), axes=1)) > + if verb: > + print("\n2nd projection - diag of double tensordot.") > + print("Proj: %s" % d[:2]) > + > + > +def proj3(vect, A, N=1, verb=True): > + for i in range(N): > + d = diagonal(tensordot(tensordot(A, vect, axes=([0],[1])), vect, > axes=([0],[1]))) > + if verb: > + print("\n3rd projection - diag of double tensordot, no > transpose.") > + print("Proj: %s" % d[:2]) > + > + > +def proj4(vect, A, N=1, verb=True): > + d = zeros(len(vect), float64) > + for i in range(N): > + a = tensordot(A, vect, axes=([0],[1])) > + for j in range(len(vect)): > + d[j] = dot(vect[j], a[:,j]) > + if verb: > + print("\n4th projection - mixed tensordot() and per-vector > dot().") > + print("Proj: %s" % d[:2]) > + > + > +def proj5(vect, A, N=1, verb=True): > + d = zeros(len(vect), float64) > + for i in range(N): > + vect2 = vect**2 > + double_vect = 2.0 * vect > + for j in xrange(len(vect)): > + d[j] = vect2[j][0]*A[0, 0] + vect2[j][1]*A[1, 1] + > vect2[j][2]*(A[2, 2]) + double_vect[j][0]*vect[j][1]*A[0, 1] + > double_vect[j][0]*vect[j][2]*A[0, 2] + double_vect[j][1]*vect[j][2]*A[1, 2] > + if verb: > + print("\n5th projection - expansion and sum.") > + print("Proj: %s" % d[:2]) > + > + > +def proj6(vect, A, N=1, verb=True): > + d = zeros(len(vect), float64) > + for i in range(N): > + d = vect[:,0]**2*A[0, 0] + vect[:,1]**2*A[1, 1] + > vect[:,2]**2*(A[2, 2]) + 2.0*vect[:,0]*vect[:,1]*A[0, 1] + > 2.0*vect[:,0]*vect[:,2]*A[0, 2] + 2.0*vect[:,1]*vect[:,2]*A[1, 2] > + if verb: > + print("\n6th projection - expansion.") > + print("Proj: %s" % d[:2]) > + > + > +def proj7(vect, A, N=1, verb=True): > + d = zeros(len(vect), float64) > + for i in range(N): > + vect2 = vect**2 > + double_vect = 2.0 * vect > + d = vect2[:,0]*A[0, 0] + vect2[:,1]*A[1, 1] + vect2[:,2]*(A[2, > 2]) + double_vect[:,0]*vect[:,1]*A[0, 1] + double_vect[:,0]*vect[:,2]*A[0, > 2] + double_vect[:,1]*vect[:,2]*A[1, 2] > + if verb: > + print("\n7th projection - expansion with pre-calculation.") > + print("Proj: %s" % d[:2]) > + > + > +# Some 200 vectors. > +vect = array([[1, 2, 3], [2, 2, 2]], float64) > +vect = tile(vect, (100, 1)) > +if __name__ == '__main__': > + print("Original vectors:\n%s\nShape: %s" % (vect[:2], vect.shape)) > + > +# Init the alignment tensor. > +A = zeros((3, 3), float64) > +A_5D = [1.42219822168827662867e-04, -1.44543001566521341940e-04, > -7.07796211648713973798e-04, -6.01619494082773244303e-04, > 2.02008007072950861996e-04] > +A[0, 0] = A_5D[0] > +A[1, 1] = A_5D[1] > +A[2, 2] = -A_5D[0] -A_5D[1] > +A[0, 1] = A[1, 0] = A_5D[2] > +A[0, 2] = A[2, 0] = A_5D[3] > +A[1, 2] = A[2, 1] = A_5D[4] > +A *= 1e4 > +if __name__ == '__main__': > + print("\nTensor:\n%s\n" % A) > + > +# Projections. > +N = 100 > +M = 100 > +if __name__ == '__main__': > + proj1(vect=vect, A=A, N=1, verb=True) > + print("Timing (s): %s" % timeit("proj1(vect=vect, A=A, N=N, > verb=False)", setup="from tensor_projections import proj1, vect, A, N", > number=M)) > + > + proj2(vect=vect, A=A, N=1, verb=True) > + print("Timing (s): %s" % timeit("proj2(vect=vect, A=A, N=N, > verb=False)", setup="from tensor_projections import proj2, vect, A, N", > number=M)) > + > + proj3(vect=vect, A=A, N=1, verb=True) > + print("Timing (s): %s" % timeit("proj3(vect=vect, A=A, N=N, > verb=False)", setup="from tensor_projections import proj3, vect, A, N", > number=M)) > + > + proj4(vect=vect, A=A, N=1, verb=True) > + print("Timing (s): %s" % timeit("proj4(vect=vect, A=A, N=N, > verb=False)", setup="from tensor_projections import proj4, vect, A, N", > number=M)) > + > + proj5(vect=vect, A=A, N=1, verb=True) > + print("Timing (s): %s" % timeit("proj5(vect=vect, A=A, N=N, > verb=False)", setup="from tensor_projections import proj5, vect, A, N", > number=M)) > + > + proj6(vect=vect, A=A, N=1, verb=True) > + print("Timing (s): %s" % timeit("proj6(vect=vect, A=A, N=N, > verb=False)", setup="from tensor_projections import proj6, vect, A, N", > number=M)) > + > + proj7(vect=vect, A=A, N=1, verb=True) > + print("Timing (s): %s" % timeit("proj7(vect=vect, A=A, N=N, > verb=False)", setup="from tensor_projections import proj7, vect, A, N", > number=M)) > > > _______________________________________________ > 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
