Hi Troels, This all looks great. I was just wondering what the relaxation time interpolation is for? Is this interpolation or simply zero-filling the end of the relaxation time lists? I don't fully understand what is happening with the concatenate() and unique() section and what is created.
Cheers, Edward On 22 August 2014 14:48, <tlin...@nmr-relax.com> wrote: > Author: tlinnet > Date: Fri Aug 22 14:48:11 2014 > New Revision: 25214 > > URL: http://svn.gna.org/viewcvs/relax?rev=25214&view=rev > Log: > Modified following functions: > > Time points are now saved at the [ei][mi][oi][di] index level. > At this index Ãlevelall time points are saved for the R2eff point. > > - interpolate_disp() > To interpolate time points, all time points through the original dispersion > points di, are > collected and then made unique. This time list can potentially be the largest > of all time lists. > > - interpolate_offset() > To interpolate time points, all time points through the original offset > points, and then dispersion points di, are > collected and then made unique. This time list can potentially be the largest > of all time lists. > > - plot_disp_curves_to_file() > To acquire the original relax_times points. > > - return_r2eff_arrays() > To save all time points on the level of [ei][mi][oi][di]. > At this index level, it will be a numpy array list with all time values used > for fitting. > > bug #22461(https://gna.org/bugs/?22461): NS R1rho 2-site_fit_r1 has extremely > high chi2 value in systemtest Relax_disp.test_r1rho_kjaergaard_missing_r1. > > Modified: > trunk/specific_analyses/relax_disp/data.py > > Modified: trunk/specific_analyses/relax_disp/data.py > URL: > http://svn.gna.org/viewcvs/relax/trunk/specific_analyses/relax_disp/data.py?rev=25214&r1=25213&r2=25214&view=diff > ============================================================================== > --- trunk/specific_analyses/relax_disp/data.py (original) > +++ trunk/specific_analyses/relax_disp/data.py Fri Aug 22 14:48:11 2014 > @@ -52,7 +52,7 @@ > > # Python module imports. > from math import cos, pi, sin, sqrt > -from numpy import array, float64, int32, ones, zeros > +from numpy import array, concatenate, float64, int32, ones, unique, zeros > from os import F_OK, access > from os.path import expanduser > from random import gauss > @@ -775,7 +775,7 @@ > desel_spin(spin_id) > > > -def interpolate_disp(spin=None, spin_id=None, si=None, num_points=None, > extend_hz=None): > +def interpolate_disp(spin=None, spin_id=None, si=None, num_points=None, > extend_hz=None, relax_times=None): > """Interpolate function for 2D Grace plotting function for the > dispersion curves. > > @keyword spin: The specific spin data container. > @@ -788,6 +788,8 @@ > @type num_points: int > @keyword extend_hz: How far to extend the interpolated fitted curves > to (in Hz). > @type extend_hz: float > + @keyword relax_times: The experiment specific fixed time period for > relaxation (in seconds). The dimensions are {Ei, Mi, Oi, Di, Ti}. > + @type relax_times: rank-4 list of floats > @return: The interpolated_flag, list of back calculated > R2eff/R1rho values in rad/s {Ei, Si, Mi, Oi, Di}, list of interpolated > frequencies for cpmg_frqs in Hz {Ei, Si, Mi, Oi, Di}, interpolated spin-lock > offsets in rad/s {Ei, Si, Mi, Oi}, list of interpolated spin-lock field > strength frequencies for spin_lock_nu1_new in Hz {Ei, Si, Mi, Oi, Di}, > chemical shifts in rad/s {Ei, Si, Mi}, interpolated rotating frame tilt > angles theta {Ei, Si, Mi, Oi, Di}, interpolated average resonance offset in > the rotating frame Omega in rad/s {Ei, Si, Mi, Oi, Di} and the interpolated > effective field in rotating frame w_eff in rad/s {Ei, Si, Mi, Oi, Di}. > @rtype: boolean, rank-4 list of numpy rank-1 float > arrays, rank-4 list of numpy rank-1 float arrays, rank-3 list of numpy rank-1 > float arrays, rank-4 list of numpy rank-1 float arrays, rank-2 list of numpy > rank-1 float arrays, rank-4 list of numpy rank-1 float arrays, rank-4 list of > numpy rank-1 float arrays, rank-4 list of numpy rank-1 float arrays > """ > @@ -797,6 +799,7 @@ > # Initialise some structures. > cpmg_frqs_new = None > spin_lock_nu1_new = None > + relax_times_new = None > > # Interpolate the CPMG frequencies (numeric models). > if spin.model in MODEL_LIST_NUMERIC_CPMG or spin.model in [MODEL_B14, > MODEL_B14_FULL]: > @@ -871,28 +874,42 @@ > > if spin_lock_nu1 != None and len(spin_lock_nu1[0][0][0]): > spin_lock_nu1_new = [] > + relax_times_new = [] > for ei in range(len(spin_lock_nu1)): > # Add a new dimension. > spin_lock_nu1_new.append([]) > + relax_times_new.append([]) > > # Then loop over the spectrometer frequencies. > for mi in range(len(spin_lock_nu1[ei])): > # Add a new dimension. > spin_lock_nu1_new[ei].append([]) > + relax_times_new[ei].append([]) > > # Finally the offsets. > for oi in range(len(spin_lock_nu1[ei][mi])): > # Add a new dimension. > spin_lock_nu1_new[ei][mi].append([]) > + relax_times_new[ei][mi].append([]) > > # No data. > if not len(spin_lock_nu1[ei][mi][oi]): > continue > + > + # There is no way to interpolate the time points correct. > + # The best suggestion is to concatenate all values at > original offset, and then make a unique list. > + relax_time_temp = array([]) > + for di_o, times in enumerate(relax_times[ei][mi][oi]): > + relax_time_temp = concatenate( (relax_time_temp, > times) ) > + > + # Make a unique list. > + relax_time_temp = unique(relax_time_temp) > > # Interpolate (adding the extended amount to the end). > for di in range(num_points): > point = (di + 1) * > (max(spin_lock_nu1[ei][mi][oi])+extend_hz) / num_points > spin_lock_nu1_new[ei][mi][oi].append(point) > + relax_times_new[ei][mi][oi].append(relax_time_temp) > > # Convert to a numpy array. > spin_lock_nu1_new[ei][mi][oi] = > array(spin_lock_nu1_new[ei][mi][oi], float64) > @@ -914,12 +931,12 @@ > back_calc = None > else: > # Back calculate R2eff data for the second sets of plots. > - back_calc = > specific_analyses.relax_disp.optimisation.back_calc_r2eff(spin=spin, > spin_id=spin_id, cpmg_frqs=cpmg_frqs_new, spin_lock_nu1=spin_lock_nu1_new) > + back_calc = > specific_analyses.relax_disp.optimisation.back_calc_r2eff(spin=spin, > spin_id=spin_id, cpmg_frqs=cpmg_frqs_new, spin_lock_nu1=spin_lock_nu1_new, > relax_times_new=relax_times_new) > > return interpolated_flag, back_calc, cpmg_frqs_new, offsets, > spin_lock_fields_inter, chemical_shifts, tilt_angles, Delta_omega, w_eff > > > -def interpolate_offset(spin=None, spin_id=None, si=None, num_points=None, > extend_ppm=None): > +def interpolate_offset(spin=None, spin_id=None, si=None, num_points=None, > extend_ppm=None, relax_times=None): > """Interpolate function for 2D Grace plotting function for the > dispersion curves, interpolating through spin-lock offset in rad/s. > > @keyword spin: The specific spin data container. > @@ -932,6 +949,8 @@ > @type num_points: int > @keyword extend_ppm: How far to extend the interpolated fitted curves > to in offset ppm. > @type extend_ppm: float > + @keyword relax_times: The experiment specific fixed time period for > relaxation (in seconds). The dimensions are {Ei, Mi, Oi, Di, Ti}. > + @type relax_times: rank-4 list of floats > @return: The interpolated_flag, list of back calculated > R2eff/R1rho values in rad/s {Ei, Si, Mi, Oi, Di}, list of interpolated > frequencies for cpmg_frqs in Hz {Ei, Si, Mi, Oi, Di}, interpolated spin-lock > offsets in rad/s {Ei, Si, Mi, Oi}, list of interpolated spin-lock field > strength frequencies for spin_lock_nu1_new in Hz {Ei, Si, Mi, Oi, Di}, > chemical shifts in rad/s {Ei, Si, Mi}, interpolated rotating frame tilt > angles theta {Ei, Si, Mi, Oi, Di}, interpolated average resonance offset in > the rotating frame Omega in rad/s {Ei, Si, Mi, Oi, Di} and the interpolated > effective field in rotating frame w_eff in rad/s {Ei, Si, Mi, Oi, Di}. > @rtype: boolean, rank-4 list of numpy rank-1 float > arrays, rank-4 list of numpy rank-1 float arrays, rank-3 list of numpy rank-1 > float arrays, rank-4 list of numpy rank-1 float arrays, rank-2 list of numpy > rank-1 float arrays, rank-4 list of numpy rank-1 float arrays, rank-4 list of > numpy rank-1 float arrays, rank-4 list of numpy rank-1 float arrays > """ > @@ -941,6 +960,7 @@ > > # Initialise some structures. > spin_lock_offset_new = [] > + relax_times_new = None > > # Get the spin-lock field strengths. > spin_lock_nu1 = return_spin_lock_nu1(ref_flag=False) > @@ -985,11 +1005,42 @@ > # The offset data. > offsets, spin_lock_fields_inter, chemical_shifts, tilt_angles, > Delta_omega, w_eff = return_offset_data(spins=[spin], spin_ids=[spin_id], > field_count=field_count, spin_lock_offset=spin_lock_offset_new, > fields=spin_lock_nu1) > > + # Interpolated relaxation time. > + if tilt_angles != None and len(tilt_angles[0][0][0]): > + relax_times_new = [] > + for ei in range(len(tilt_angles)): > + # Add a new dimension. > + relax_times_new.append([]) > + > + # Then loop over the spectrometer frequencies. > + for mi in range(len(tilt_angles[ei])): > + # Add a new dimension. > + relax_times_new[ei].append([]) > + > + # There is no way to interpolate the time points correct. > + # The best suggestion is to concatenate all values at > original offset and dispersion point, and then make a unique list. > + relax_time_temp = array([]) > + for oi_o, relax_times_oi in enumerate(relax_times[ei][mi]): > + for di_o, times in enumerate(relax_times_oi): > + relax_time_temp = concatenate( (relax_time_temp, > times) ) > + > + # Make a unique list. > + relax_time_temp = unique(relax_time_temp) > + > + # Finally the offsets. > + for oi in range(len(tilt_angles[ei][0][mi])): > + # Add a new dimension. > + relax_times_new[ei][mi].append([]) > + > + # Interpolate (adding the extended amount to the end). > + for di in range(len(tilt_angles[ei][0][mi][oi])): > + relax_times_new[ei][mi][oi].append(relax_time_temp) > + > if spin.model == MODEL_R2EFF: > back_calc = None > else: > # Back calculate R2eff data for the second sets of plots. > - back_calc = > specific_analyses.relax_disp.optimisation.back_calc_r2eff(spin=spin, > spin_id=spin_id, spin_lock_offset=spin_lock_offset_new, > spin_lock_nu1=spin_lock_fields_inter) > + back_calc = > specific_analyses.relax_disp.optimisation.back_calc_r2eff(spin=spin, > spin_id=spin_id, spin_lock_offset=spin_lock_offset_new, > spin_lock_nu1=spin_lock_fields_inter, relax_times_new=relax_times_new) > > # cpmg_frqs are not interpolated. > cpmg_frqs_new = None > @@ -1912,6 +1963,16 @@ > linetype = [] > linestyle = [] > > + # Number of spectrometer fields. > + fields = [None] > + field_count = 1 > + if hasattr(cdp, 'spectrometer_frq_count'): > + fields = cdp.spectrometer_frq_list > + field_count = cdp.spectrometer_frq_count > + > + # Get the relax_times. > + values, errors, missing, frqs, frqs_H, exp_types, relax_times = > return_r2eff_arrays(spins=[spin], spin_ids=[spin_id], fields=fields, > field_count=field_count) > + > # Set up the interpolated curve data structures. > interpolated_flag = False > > @@ -1920,11 +1981,11 @@ > > if interpolate == INTERPOLATE_DISP: > # Interpolate through disp points. > - interpolated_flag, back_calc, cpmg_frqs_new, offsets_inter, > spin_lock_nu1_new, chemical_shifts, tilt_angles_inter, Delta_omega_inter, > w_eff_inter = interpolate_disp(spin=spin, spin_id=spin_id, si=si, > num_points=num_points, extend_hz=extend_hz) > + interpolated_flag, back_calc, cpmg_frqs_new, offsets_inter, > spin_lock_nu1_new, chemical_shifts, tilt_angles_inter, Delta_omega_inter, > w_eff_inter = interpolate_disp(spin=spin, spin_id=spin_id, si=si, > num_points=num_points, extend_hz=extend_hz, relax_times=relax_times) > > elif interpolate == INTERPOLATE_OFFSET: > # Interpolate through disp points. > - interpolated_flag, back_calc, cpmg_frqs_new, offsets_inter, > spin_lock_nu1_new, chemical_shifts, tilt_angles_inter, Delta_omega_inter, > w_eff_inter = interpolate_offset(spin=spin, spin_id=spin_id, si=si, > num_points=num_points, extend_ppm=extend_ppm) > + interpolated_flag, back_calc, cpmg_frqs_new, offsets_inter, > spin_lock_nu1_new, chemical_shifts, tilt_angles_inter, Delta_omega_inter, > w_eff_inter = interpolate_offset(spin=spin, spin_id=spin_id, si=si, > num_points=num_points, extend_ppm=extend_ppm, relax_times=relax_times) > > # Do not interpolate, if model is R2eff. > if spin.model == MODEL_R2EFF: > @@ -4318,12 +4379,14 @@ > missing[ei][si].append([]) > frqs[ei][si].append(0.0) > frqs_H[ei][si].append(0.0) > + relax_times[ei].append([]) > for offset, oi in loop_offset(exp_type=exp_type, frq=frq, > return_indices=True): > values[ei][si][mi].append([]) > errors[ei][si][mi].append([]) > missing[ei][si][mi].append([]) > - for mi in range(field_count): > - relax_times[ei].append(None) > + relax_times[ei][mi].append([]) > + for point, di in loop_point(exp_type=exp_type, frq=frq, > offset=offset, return_indices=True): > + relax_times[ei][mi][oi].append([]) > > # Pack the R2eff/R1rho data. > data_flag = False > @@ -4405,34 +4468,8 @@ > errors[ei][si][mi][oi].append(current_spin.r2eff_err[key]) > > # The relaxation times. > - relax_time = [] > - for id in cdp.spectrum_ids: > - # Non-matching data. > - if cdp.spectrometer_frq[id] != frq: > - continue > - if cdp.exp_type[id] != exp_type: > - continue > - if exp_type in EXP_TYPE_LIST_CPMG: > - if id not in cdp.cpmg_frqs.keys() or cdp.cpmg_frqs[id] > != point: > - continue > - else: > - if id not in cdp.spin_lock_nu1.keys() or > cdp.spin_lock_nu1[id] != point: > - continue > - > - # Found. > - relax_time.append(cdp.relax_times[id]) > - > - # Use the maximum time value found. > - relax_time = max(relax_time) > - > - # Check the value if already set. > - if relax_times[ei][mi] != None: > - if relax_times[ei][mi] != relax_time: > - raise RelaxError("The relaxation times do not match for > all experiments.") > - continue > - > - # Store the time. > - relax_times[ei][mi] = relax_time > + for time, ti in loop_time(exp_type=exp_type, frq=frq, > offset=offset, point=point, return_indices=True): > + relax_times[ei][mi][oi][di].append(time) > > # Increment the spin index. > si += 1 > @@ -4442,7 +4479,6 @@ > raise RelaxError("No R2eff/R1rho data could be found for the spin > cluster %s." % spin_ids) > > # Convert to numpy arrays. > - relax_times = array(relax_times, float64) > for exp_type, ei in loop_exp(return_indices=True): > for si in range(spin_num): > for frq, mi in loop_frq(return_indices=True): > @@ -4450,6 +4486,8 @@ > values[ei][si][mi][oi] = array(values[ei][si][mi][oi], > float64) > errors[ei][si][mi][oi] = array(errors[ei][si][mi][oi], > float64) > missing[ei][si][mi][oi] = array(missing[ei][si][mi][oi], > int32) > + for point, di in loop_point(exp_type=exp_type, frq=frq, > offset=offset, return_indices=True): > + relax_times[ei][mi][oi][di] = > array(relax_times[ei][mi][oi][di], float64) > > # Return the structures. > return values, errors, missing, frqs, frqs_H, exp_types, relax_times > > > _______________________________________________ > relax (http://www.nmr-relax.com) > > This is the relax-commits mailing list > relax-comm...@gna.org > > 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@gna.org 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
