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 <[email protected]>: > 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 <[email protected]> > 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 <[email protected]>: > >> > >> 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, <[email protected]> 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 > >> > [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
