Hi Ed. You may see, that I have corrected this. :-)
The profiling scripts was a tremendous help here! I think the clearest image for you would be to just make a diff on the target function and on the lib function, between the branches. Best Troels 2014-06-10 11:22 GMT+02:00 Edward d'Auvergne <[email protected]>: > Hi Troels, > > This might be out of date due to the huge number of commits I'm only > starting to look through now, but there is a problem with this change. > The issue is with numpy data structure initialisations - this should > never happen in the target function itself. That is insanely far too > computationally expensive as you have many memory allocations at the > start, and automatic Python garbage collection at the end of the > function which deletes all of these structures. All numpy structures > should be created in the __init__() method (or a method called from > __init__()). You then pack the parameters into the required data > structure in the target function, recycling the numpy structures as > much as possible. > > Regards, > > Edward > > > > On 7 June 2014 21:43, <[email protected]> wrote: > > Author: tlinnet > > Date: Sat Jun 7 21:43:19 2014 > > New Revision: 23723 > > > > URL: http://svn.gna.org/viewcvs/relax?rev=23723&view=rev > > Log: > > Initial try to alter the target function calc_CR72_chi2. > > > > This is the first test to restructure the arrays, to allow for higher > dimensional computation. > > All numpy arrays have to have same shape to allow to multiply together. > > The dimensions should be [ei][si][mi][oi][di]. [Experiment][spins][spec. > frq][offset][disp points]. > > This is complicated with number of disp point can change per > spectrometer frequency. > > > > Task #7807 (https://gna.org/task/index.php?7807): Speed-up of > dispersion models for Clustered analysis. > > > > This implementation brings a high overhead. > > The first test shows no winning of time. > > The creation of arrays takes all the time. > > > > Checked on MacBook Pro > > 2.4 GHz Intel Core i5 > > 8 GB 1067 Mhz DDR3 RAM. > > Python Distribution -- Python 2.7.3 |EPD 7.3-2 (32-bit)| > > > > Timing for: > > 3 fields, [600. * 1E6, 800. * 1E6, 900. * 1E6] > > ('sfrq: ', 600000000.0, 'number of cpmg frq', 15) > > ('sfrq: ', 800000000.0, 'number of cpmg frq', 20) > > ('sfrq: ', 900000000.0, 'number of cpmg frq', 22) > > iterations of function call: 1000 > > > > Timed for simulating 1 or 100 clustered spins. > > > > For TRUNK > > > > 1 spin: > > ncalls tottime percall cumtime percall filename:lineno(function) > > 3000 0.267 0.000 0.313 0.000 cr72.py:100(r2eff_CR72) > > 1000 0.056 0.000 0.434 0.000 > relax_disp.py:456(calc_CR72_chi2) > > 3000 0.045 0.000 0.061 0.000 chi2.py:32(chi2) > > > > 100 spins: > > ncalls tottime percall cumtime percall filename:lineno(function) > > 300000 26.315 0.000 30.771 0.000 cr72.py:100(r2eff_CR72) > > 1000 5.498 0.005 42.660 0.043 > relax_disp.py:456(calc_CR72_chi2) > > 300000 4.438 0.000 6.021 0.000 chi2.py:32(chi2) > > > > TESTING > > > > 1 spin: > > ncalls tottime percall cumtime percall filename:lineno(function) > > 19013 0.278 0.000 0.278 0.000 > {numpy.core.multiarray.array} > > 1000 0.191 0.000 0.777 0.001 > relax_disp.py:457(calc_CR72_chi2) > > 1000 0.147 0.000 0.197 0.000 cr72.py:101(r2eff_CR72) > > 3000 0.044 0.000 0.061 0.000 chi2.py:32(chi2) > > > > 100 spins: > > ncalls tottime percall cumtime percall filename:lineno(function) > > 1504904 25.215 0.000 25.215 0.000 > {numpy.core.multiarray.array} > > 1000 17.261 0.017 51.180 0.051 > relax_disp.py:457(calc_CR72_chi2) > > 300000 4.637 0.000 6.310 0.000 chi2.py:32(chi2) > > > > Modified: > > branches/disp_spin_speed/target_functions/relax_disp.py > > > > 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=23723&r1=23722&r2=23723&view=diff > > > ============================================================================== > > --- branches/disp_spin_speed/target_functions/relax_disp.py > (original) > > +++ branches/disp_spin_speed/target_functions/relax_disp.py Sat Jun > 7 21:43:19 2014 > > @@ -27,6 +27,7 @@ > > # Python module imports. > > from copy import deepcopy > > from math import pi > > +import numpy as np > > from numpy import complex64, dot, float64, int16, sqrt, zeros > > > > # relax module imports. > > @@ -470,29 +471,86 @@ > > @rtype: float > > """ > > > > + # Get the shape of back_calc structure. > > + back_calc_shape = list( np.asarray(self.back_calc).shape ) > > + > > + # Find which frequency has the maximum number of disp points. > > + # To let the numpy array operate well together, the broadcast > size has to be equal for all shapes. > > + max_num_disp_points = np.max(self.num_disp_points) > > + > > + # Create numpy arrays to pass to the lib function. > > + # All numpy arrays have to have same shape to allow to multiply > together. > > + # The dimensions should be [ei][si][mi][oi][di]. > [Experiment][spins][spec. frq][offset][disp points]. > > + # The number of disp point can change per spectrometer, so we > make the maximum size. > > + R20A_a = np.ones(back_calc_shape + [max_num_disp_points]) > > + R20B_a = np.ones(back_calc_shape + [max_num_disp_points]) > > + pA_a = np.ones(back_calc_shape + [max_num_disp_points]) > > + dw_frq_a = np.ones(back_calc_shape + [max_num_disp_points]) > > + kex_a = np.ones(back_calc_shape + [max_num_disp_points]) > > + cpmg_frqs_a = np.ones(back_calc_shape + [max_num_disp_points]) > > + num_disp_points_a = np.ones(back_calc_shape + > [max_num_disp_points]) > > + back_calc_a = np.ones(back_calc_shape + [max_num_disp_points]) > > + > > + # Loop over the experiment types. > > + for ei in range(self.num_exp): > > + # Loop over the spins. > > + for si in range(self.num_spins): > > + # Loop over the spectrometer frequencies. > > + for mi in range(self.num_frq): > > + # Loop over the offsets. > > + for oi in range(self.num_offsets[ei][si][mi]): > > + # Extract number of dispersion points. > > + num_disp_points = > self.num_disp_points[ei][si][mi][oi] > > + > > + # The R20 index. > > + r20_index = mi + si*self.num_frq > > + > > + # Store r20a and r20b values per disp point. > > + R20A_a[ei][si][mi][oi] = np.array( > [R20A[r20_index]] * max_num_disp_points, float64) > > + R20B_a[ei][si][mi][oi] = np.array( > [R20B[r20_index]] * max_num_disp_points, float64) > > + > > + # Convert dw from ppm to rad/s. > > + dw_frq = dw[si] * self.frqs[ei][si][mi] > > + > > + # Store dw_frq per disp point. > > + dw_frq_a[ei][si][mi][oi] = np.array( [dw_frq] * > max_num_disp_points, float64) > > + > > + # Store pA and kex per disp point. > > + pA_a[ei][si][mi][oi] = np.array( [pA] * > max_num_disp_points, float64) > > + kex_a[ei][si][mi][oi] = np.array( [kex] * > max_num_disp_points, float64) > > + > > + # Extract cpmg_frqs and num_disp_points from > lists. > > + cpmg_frqs_a[ei][si][mi][oi][:num_disp_points] = > self.cpmg_frqs[ei][mi][oi] > > + > num_disp_points_a[ei][si][mi][oi][:num_disp_points] = > self.num_disp_points[ei][si][mi][oi] > > + > > + ## Back calculate the R2eff values. > > + r2eff_CR72(r20a=R20A_a, r20b=R20B_a, pA=pA_a, dw=dw_frq_a, > kex=kex_a, cpmg_frqs=cpmg_frqs_a, back_calc = back_calc_a, > num_points=num_disp_points_a) > > + > > # Initialise. > > chi2_sum = 0.0 > > > > - # Loop over the spins. > > - for si in range(self.num_spins): > > - # Loop over the spectrometer frequencies. > > - for mi in range(self.num_frq): > > - # The R20 index. > > - r20_index = mi + si*self.num_frq > > - > > - # Convert dw from ppm to rad/s. > > - dw_frq = dw[si] * self.frqs[0][si][mi] > > - > > - # Back calculate the R2eff values. > > - r2eff_CR72(r20a=R20A[r20_index], r20b=R20B[r20_index], > pA=pA, dw=dw_frq, kex=kex, cpmg_frqs=self.cpmg_frqs[0][mi][0], back_calc = > self.back_calc[0][si][mi][0], num_points=self.num_disp_points[0][si][mi][0]) > > - > > - # For all missing data points, set the back-calculated > value to the measured values so that it has no effect on the chi-squared > value. > > - for di in range(self.num_disp_points[0][si][mi][0]): > > - if self.missing[0][si][mi][0][di]: > > - self.back_calc[0][si][mi][0][di] = > self.values[0][si][mi][0][di] > > - > > - # Calculate and return the chi-squared value. > > - chi2_sum += chi2(self.values[0][si][mi][0], > self.back_calc[0][si][mi][0], self.errors[0][si][mi][0]) > > + # Now return the values back to the structure of self.back_calc > object. > > + ## For all missing data points, set the back-calculated value > to the measured values so that it has no effect on the chi-squared value. > > + # Loop over the experiment types. > > + for ei in range(self.num_exp): > > + # Loop over the spins. > > + for si in range(self.num_spins): > > + # Loop over the spectrometer frequencies. > > + for mi in range(self.num_frq): > > + # Loop over the offsets. > > + for oi in range(self.num_offsets[ei][si][mi]): > > + # Extract number of dispersion points. > > + num_disp_points = > self.num_disp_points[ei][si][mi][oi] > > + > > + self.back_calc[ei][si][mi][oi] = > back_calc_a[ei][si][mi][oi][:num_disp_points] > > + > > + > > + for di in > range(self.num_disp_points[ei][si][mi][oi]): > > + if self.missing[ei][si][mi][oi][di]: > > + self.back_calc[ei][si][mi][oi][di] = > self.values[ei][si][mi][oi][di] > > + > > + ## Calculate and return the chi-squared value. > > + chi2_sum += chi2(self.values[ei][si][mi][oi], > self.back_calc[ei][si][mi][oi], self.errors[ei][si][mi][oi]) > > > > # Return the total chi-squared value. > > return chi2_sum > > > > > > _______________________________________________ > > 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 > _______________________________________________ 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
