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