Re: r23901 - /branches/disp_spin_speed/target_functions/relax_disp.py -- June 13, 2014

I'm sure there'll be a way in MATLAB, but it's definitely not the
default (http://www.mathworks.de/de/help/matlab/ref/times.html).  The
element wise operation is very powerful for this purpose.

Regards,

Edward



On 13 June 2014 10:55, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx> wrote:

The whole trick about this converting, is the ability of numpy arrays not
working
as matrixes, but as element wise operations, per dimension.

That would not be possible in MATLAB?

Best
Troels


2014-06-13 10:53 GMT+02:00 Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx>:

Hi Ed.

I dont know yet.
I am looking "forward" to start battling with the numeric.

It should by logic be possible to do the same tricks.
Higher dimensions, and then power of the dimensions.

Lets see.
Small models first.

Best
T


2014-06-13 10:51 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:

How would you use this efficiently in the numeric dispersion models
where num_disp_points is currently used to loop over each dispersion
point?

Regards,

Edward


On 13 June 2014 10:40, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx>
wrote:

Hi Ed.

Well, the disp_struct is filled out with zeros, where there are no disp
point
per experiment type, field strength, or offset.

So it "should" be the same.

That is why, I have to go "all the way" to looping over disp points in
the
init function.
For cpmg, I can stop one loop before, and fill 1.0 up to [:disp_points]

I had to fight a little, to realise this.

Best
Troels


2014-06-13 10:36 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:

Hi,

This is probably still required for the numeric models, but it can be
removed from most analytic models.  As for disp_struct, as this goes
up to self.max_num_disp_points, it is not the same thing when
different numbers of dispersion points are used per experiment type,
field strength, or offset.  The self.num_disp_points structure is a
numpy array of ND, whereas self.disp_struct is one rank higher where
the ND have been converted into the new dimension.  Oh, once you have
everything converted, you'll also find a lot of old code in __init__()
to kill :)

Regards,

Edward

On 13 June 2014 10:18, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx>
wrote:

Hi ed.

All the:
num_points=self.num_disp_points_a

can also be killed.

They are not used.
The disp_struct is actually this structure, in higher dimensions?

Best
Troels


2014-06-13 9:02 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:


Hi Troels,

Thanks to the lightning quick infrastructure you are putting into
place, we can also simplify the target_functions.relax_disp to
lib.dispersion API.  You will notice a lot of code like in this
TP02
model:

+        # Once off parameter conversions.
+        pB = 1.0 - pA

This was originally in lib.dispersion (well at least for some of
the
models), but I shifted it into the func_*() target functions to
speed
the code up, as then the calculation would happen only once rather
than once for each iteration of that massive looping beast you have
killed.

So now that the lib.dispersion functions are only called once per
target function call, all of these 'Once off parameter conversions'
can be shifted back into the lib.dispersion functions.  Then the
number of arguments for these functions will drop significantly, as
the {pB,  k_AB,  k_BA} parameters will no longer need to be passed
in.
This will be far more significant for the 3-site models where there
are many, many parameter conversions.  This will have the added
benefit of simplifying the use of the lib.dispersion modules
outside
of the dispersion target functions, for example with the unit
testing.

Cheers,

Edward




On 13 June 2014 07:21,  <tlinnet@xxxxxxxxxxxxx> wrote:

Author: tlinnet
Date: Fri Jun 13 07:21:02 2014
New Revision: 23901

URL: http://svn.gna.org/viewcvs/relax?rev=23901&view=rev
Log:
Replaced the loop structure in target function of TAP03 with
numpy
arrays.

This makes the model faster.

Task #7807 (https://gna.org/task/index.php?7807): Speed-up of
dispersion
models for Clustered analysis.

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=23901&r1=23900&r2=23901&view=diff



==============================================================================
--- branches/disp_spin_speed/target_functions/relax_disp.py
(original)
+++ branches/disp_spin_speed/target_functions/relax_disp.py
Fri
Jun
13 07:21:02 2014
@@ -395,7 +395,7 @@
             self.func = self.func_ns_mmq_3site_linear

         # Setup special numpy array structures, for higher
dimensional
computation.
-        if model in [MODEL_B14, MODEL_B14_FULL, MODEL_CR72,
MODEL_CR72_FULL, MODEL_DPL94, MODEL_TSMFK01]:
+        if model in [MODEL_B14, MODEL_B14_FULL, MODEL_CR72,
MODEL_CR72_FULL, MODEL_DPL94, MODEL_TAP03, MODEL_TSMFK01]:
             # Get the shape of back_calc structure.
             # If using just one field, or having the same number
of
dispersion points, the shape would extend to that number.
             # Shape has to be: [ei][si][mi][oi].
@@ -435,10 +435,12 @@
             self.power_a = ones(self.numpy_array_shape, int16)

             # For R1rho data.
-            if model in [MODEL_DPL94]:
+            if model in [MODEL_DPL94, MODEL_TAP03]:
                 self.tilt_angles_a = deepcopy(zeros_a)
                 self.spin_lock_omega1_squared_a =
deepcopy(zeros_a)
+                self.spin_lock_omega1_a = deepcopy(zeros_a)
                 self.phi_ex_struct = deepcopy(zeros_a)
+                self.offset_a = deepcopy(zeros_a)

             self.frqs_a = deepcopy(zeros_a)
             self.disp_struct = deepcopy(zeros_a)
@@ -447,6 +449,7 @@
             # Create special numpy structures.
             # Structure of dw. The full and the outer dimensions
structures.
             self.dw_struct = deepcopy(zeros_a)
+            self.no_nd_struct = ones([self.NO, self.ND],
float64)
             self.nm_no_nd_struct = ones([self.NM, self.NO,
self.ND],
float64)

             # Structure of r20a and r20b. The full and outer
dimensions
structures.
@@ -459,10 +462,11 @@
                 # Expand relax times.
                 self.inv_relax_times_a = 1.0 / multiply.outer(
tile(self.relax_times[:,None],(1, 1, self.NS)).reshape(self.NE,
self.NS,
self.NM), self.no_nd_struct )

-            if model in [MODEL_DPL94]:
+            if model in [MODEL_DPL94, MODEL_TAP03]:
                 self.r1_a = multiply.outer(
self.r1.reshape(self.NE,
self.NS, self.NM), self.no_nd_struct )
-
-            # Extract the frequencies to numpy array.
+                self.chemical_shifts_a = multiply.outer(
self.chemical_shifts, self.no_nd_struct )
+
+           # Extract the frequencies to numpy array.
             self.frqs_a = multiply.outer(
asarray(self.frqs).reshape(self.NE, self.NS, self.NM),
self.no_nd_struct )

             # Loop over the experiment types.
@@ -476,7 +480,7 @@
                             # Extract number of dispersion
points.
                             num_disp_points =
self.num_disp_points[ei][si][mi][oi]

-                            if model not in [MODEL_DPL94]:
+                            if model not in [MODEL_DPL94,
MODEL_TAP03]:
                                 # Extract cpmg_frqs and
num_disp_points
from lists.

self.cpmg_frqs_a[ei][si][mi][oi][:num_disp_points] =
self.cpmg_frqs[ei][mi][oi]

self.num_disp_points_a[ei][si][mi][oi][:num_disp_points] =
self.num_disp_points[ei][si][mi][oi]
@@ -498,12 +502,14 @@

self.power_a[ei][si][mi][oi][di]
=
int(round(self.cpmg_frqs[ei][mi][0][di] *
self.relax_times[ei][mi]))

self.tau_cpmg_a[ei][si][mi][oi][di]
= 0.25 / self.cpmg_frqs[ei][mi][0][di]
                                 # For R1rho data.
-                                if model in [MODEL_DPL94]:
+                                if model in [MODEL_DPL94,
MODEL_TAP03]:

self.disp_struct[ei][si][mi][oi][di] = 1.0

                                     # Extract the frequencies to
numpy
array.

self.tilt_angles_a[ei][si][mi][oi][di] =
self.tilt_angles[ei][si][mi][oi][di]

self.spin_lock_omega1_squared_a[ei][si][mi][oi][di] =
self.spin_lock_omega1_squared[ei][mi][oi][di]
+
self.spin_lock_omega1_a[ei][si][mi][oi][di] =
self.spin_lock_omega1[ei][mi][oi][di]
+
self.offset_a[ei][si][mi][oi] =
self.offset[ei][si][mi][oi]

                                     if spin_lock_nu1 != None and
len(spin_lock_nu1[ei][mi][oi]):

self.num_disp_points_a[ei][si][mi][oi][di] = num_disp_points
@@ -1908,6 +1914,49 @@
         # Once off parameter conversions.
         pB = 1.0 - pA

+        # Convert dw from ppm to rad/s. Use the out argument, to
pass
directly to structure.
+        multiply( multiply.outer( dw.reshape(self.NE, self.NS),
self.nm_no_nd_struct ), self.frqs_struct, out=self.dw_struct )
+
+        # Reshape R20 to per experiment, spin and frequency.
+        self.r20_struct[:] = multiply.outer(
R20.reshape(self.NE,
self.NS, self.NM), self.no_nd_struct )
+
+        # Back calculate the R1rho values.
+        r1rho_TAP03(r1rho_prime=self.r20_struct,
omega=self.chemical_shifts_a, offset=self.offset_a, pA=pA, pB=pB,
dw=self.dw_struct, kex=kex, R1=self.r1_a,
spin_lock_fields=self.spin_lock_omega1_a,
spin_lock_fields2=self.spin_lock_omega1_squared_a,
back_calc=self.back_calc_a, num_points=self.num_disp_points_a)
+
+        # Clean the data for all values, which is left over at
the
end
of arrays.
+        self.back_calc_a = self.back_calc_a*self.disp_struct
+
+        ## 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.
+        if self.has_missing:
+            # Replace with values.
+            self.back_calc_a[self.mask_replace_blank.mask] =
self.values_a[self.mask_replace_blank.mask]
+
+        # Return the total chi-squared value.
+        return chi2_rankN(self.values_a, self.back_calc_a,
self.errors_a)
+
+
+    def func_TP02(self, params):
+        """Target function for the Trott and Palmer (2002) R1rho
off-resonance 2-site model.
+
+        @param params:  The vector of parameter values.
+        @type params:   numpy rank-1 float array
+        @return:        The chi-squared value.
+        @rtype:         float
+        """
+
+        # Scaling.
+        if self.scaling_flag:
+            params = dot(params, self.scaling_matrix)
+
+        # Unpack the parameter values.
+        R20 = params[:self.end_index[0]]
+        dw = params[self.end_index[0]:self.end_index[1]]
+        pA = params[self.end_index[1]]
+        kex = params[self.end_index[1]+1]
+
+        # Once off parameter conversions.
+        pB = 1.0 - pA
+
         # Initialise.
         chi2_sum = 0.0

@@ -1924,7 +1973,7 @@
                 # Loop over the offsets.
                 for oi in range(self.num_offsets[0][si][mi]):
                     # Back calculate the R1rho values.
-                    r1rho_TAP03(r1rho_prime=R20[r20_index],
omega=self.chemical_shifts[0][si][mi],
offset=self.offset[0][si][mi][oi],
pA=pA, pB=pB, dw=dw_frq, kex=kex, R1=self.r1[si, mi],
spin_lock_fields=self.spin_lock_omega1[0][mi][oi],
spin_lock_fields2=self.spin_lock_omega1_squared[0][mi][oi],
back_calc=self.back_calc[0][si][mi][oi],
num_points=self.num_disp_points[0][si][mi][oi])
+                    r1rho_TP02(r1rho_prime=R20[r20_index],
omega=self.chemical_shifts[0][si][mi],
offset=self.offset[0][si][mi][oi],
pA=pA, pB=pB, dw=dw_frq, kex=kex, R1=self.r1[si, mi],
spin_lock_fields=self.spin_lock_omega1[0][mi][oi],
spin_lock_fields2=self.spin_lock_omega1_squared[0][mi][oi],
back_calc=self.back_calc[0][si][mi][oi],
num_points=self.num_disp_points[0][si][mi][oi])

                     # 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][oi]):
@@ -1938,58 +1987,6 @@
         return chi2_sum


-    def func_TP02(self, params):
-        """Target function for the Trott and Palmer (2002) R1rho
off-resonance 2-site model.
-
-        @param params:  The vector of parameter values.
-        @type params:   numpy rank-1 float array
-        @return:        The chi-squared value.
-        @rtype:         float
-        """
-
-        # Scaling.
-        if self.scaling_flag:
-            params = dot(params, self.scaling_matrix)
-
-        # Unpack the parameter values.
-        R20 = params[:self.end_index[0]]
-        dw = params[self.end_index[0]:self.end_index[1]]
-        pA = params[self.end_index[1]]
-        kex = params[self.end_index[1]+1]
-
-        # Once off parameter conversions.
-        pB = 1.0 - pA
-
-        # 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]
-
-                # Loop over the offsets.
-                for oi in range(self.num_offsets[0][si][mi]):
-                    # Back calculate the R1rho values.
-                    r1rho_TP02(r1rho_prime=R20[r20_index],
omega=self.chemical_shifts[0][si][mi],
offset=self.offset[0][si][mi][oi],
pA=pA, pB=pB, dw=dw_frq, kex=kex, R1=self.r1[si, mi],
spin_lock_fields=self.spin_lock_omega1[0][mi][oi],
spin_lock_fields2=self.spin_lock_omega1_squared[0][mi][oi],
back_calc=self.back_calc[0][si][mi][oi],
num_points=self.num_disp_points[0][si][mi][oi])
-
-                    # 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][oi]):
-                        if self.missing[0][si][mi][oi][di]:
-                            self.back_calc[0][si][mi][oi][di] =
self.values[0][si][mi][oi][di]
-
-                    # Calculate and return the chi-squared
value.
-                    chi2_sum += chi2(self.values[0][si][mi][oi],
self.back_calc[0][si][mi][oi], self.errors[0][si][mi][oi])
-
-        # Return the total chi-squared value.
-        return chi2_sum
-
-
     def func_TSMFK01(self, params):
         """Target function for the the Tollinger et al. (2001)
2-site
very-slow exchange model, range of microsecond to second time
scale.



_______________________________________________
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: r23901 - /branches/disp_spin_speed/target_functions/relax_disp.py

Header

Content

Related Messages