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 Troels Emtekær Linnet on August 22, 2014 - 15:13:
Hi Edward.

Well, I could not figure out what an interpolated time point would be
for new points.
So I collect all time points at the original offset and dispersion,
and make a concatenated list.
Then I make this unique.

The decision for the "correct" time point will be in the target function.

For minimising a model, the corresponding time point can be extracted.
But for interpolated points, I simply don't know.

Best
Troels

2014-08-22 14:57 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:
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

_______________________________________________
relax (http://www.nmr-relax.com)

This is the relax-devel mailing list
relax-devel@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-devel



Related Messages


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