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
