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, <tlinnet@xxxxxxxxxxxxx> 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
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
_______________________________________________
relax (http://www.nmr-relax.com)
This is the relax-devel mailing list
relax-devel@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-devel
Re: r23723 - /branches/disp_spin_speed/target_functions/relax_disp.py