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 modelswhere 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 intheinit 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 tospeedthe 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 passedin.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 unittesting.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.pyFriJun 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 numberofdispersion 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 dispersionpoints.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 attheend 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 noeffecton 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 noeffecton 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 timescale._______________________________________________ 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