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

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