mailRe: r25214 - /trunk/specific_analyses/relax_disp/data.py


Others Months | Index by Date | Thread Index
>>   [Date Prev] [Date Next] [Thread Prev] [Thread Next]

Header


Content

Posted by Edward d'Auvergne on August 22, 2014 - 14:58:
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,  <tlinnet@xxxxxxxxxxxxx> 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-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



Related Messages


Powered by MHonArc, Updated Fri Aug 22 16:40:22 2014