Author: tlinnet
Date: Sun Jun 8 13:14:34 2014
New Revision: 23728
URL: http://svn.gna.org/viewcvs/relax?rev=23728&view=rev
Log:
Modified profiling script to calculate correct values when setting up R2eff
values.
This is to test, that the return of chi2 gets zero.
Task #7807 (https://gna.org/task/index.php?7807): Speed-up of dispersion
models for Clustered analysis.
Modified:
branches/disp_spin_speed/test_suite/shared_data/dispersion/profiling/profiling_cr72.py
Modified:
branches/disp_spin_speed/test_suite/shared_data/dispersion/profiling/profiling_cr72.py
URL:
http://svn.gna.org/viewcvs/relax/branches/disp_spin_speed/test_suite/shared_data/dispersion/profiling/profiling_cr72.py?rev=23728&r1=23727&r2=23728&view=diff
==============================================================================
---
branches/disp_spin_speed/test_suite/shared_data/dispersion/profiling/profiling_cr72.py
(original)
+++
branches/disp_spin_speed/test_suite/shared_data/dispersion/profiling/profiling_cr72.py
Sun Jun 8 13:14:34 2014
@@ -25,6 +25,7 @@
import cProfile
from os import getcwd, path
from numpy import array, arange, asarray, int32, float64, ones, pi, zeros
+import numpy as np
import pstats
import sys
import tempfile
@@ -44,6 +45,8 @@
# relax module imports.
from lib.physical_constants import g1H, g15N
+from lib.dispersion.cr72 import r2eff_CR72
+from target_functions.chi2 import chi2
from target_functions.relax_disp import Dispersion
from specific_analyses.relax_disp.variables import EXP_TYPE_CPMG_SQ,
MODEL_B14_FULL, MODEL_CR72, MODEL_CR72_FULL, MODEL_NS_CPMG_2SITE_3D_FULL,
MODEL_NS_CPMG_2SITE_STAR_FULL
@@ -91,7 +94,7 @@
Class Profile inherits the Dispersion container class object.
"""
- def __init__(self, num_spins=1, model=None):
+ def __init__(self, num_spins=1, model=None, r2=None, r2a=None, r2b=None,
dw=None, pA=None, kex=None, spins_params=None):
"""
Special method __init__() is called first (acts as Constructor).
It brings in data from outside the class like the variable num_spins.
@@ -100,6 +103,25 @@
the name self is used by convention. Assigning num_spins to
self.num_spins allows it
to be passed to all methods within the class. Think of self as a
carrier,
or if you want impress folks call it target instance object.
+
+ @keyword num_spins: Number of spins in the cluster.
+ @type num_spins: integer
+ @keyword model: The dispersion model to instantiate the
Dispersion class with.
+ @type model: string
+ @keyword r2: The transversal relaxation rate.
+ @type r2: float
+ @keyword r2a: The transversal relaxation rate for state A
in the absence of exchange.
+ @type r2a: float
+ @keyword r2b: The transversal relaxation rate for state B
in the absence of exchange.
+ @type r2b: float
+ @keyword dw: The chemical exchange difference between
states A and B in ppm.
+ @type dw: float
+ @keyword pA: The population of state A.
+ @type pA: float
+ @keyword kex: The rate of exchange.
+ @type kex: float
+ @keyword spins_params: List of parameter strings used in dispersion
model.
+ @type spins_params: array of strings
"""
# Define parameters
@@ -114,16 +136,23 @@
# Required data structures.
self.relax_times = self.fields / (100 * 100. *1E6 )
self.ncycs = []
+ self.points = []
+ self.value = []
+ self.error = []
for i in range(len(self.fields)):
- inp = arange(2, 1000. * self.relax_times[i], 4)
- self.ncycs.append(inp)
- print("sfrq: ", self.fields[i], "number of cpmg frq", len(inp))
- self.ncycs = asarray(self.ncycs)
-
- self.points = self.ncycs / self.relax_times
-
- self.value = self.points * 1.1
- self.error = self.value * 0.1
+ ncyc = arange(2, 1000. * self.relax_times[i], 4)
+ #ncyc = arange(2, 42, 2)
+ self.ncycs.append(ncyc)
+ print("sfrq: ", self.fields[i], "number of cpmg frq", len(ncyc),
ncyc)
+
+ cpmg_point = ncyc / self.relax_times[i]
+
+ self.points.append(list(cpmg_point))
+ self.value.append([2.0]*len(cpmg_point))
+ self.error.append([1.0]*len(cpmg_point))
+
+ # Assemble param vector.
+ self.params = self.assemble_param_vector(r2=r2, r2a=r2a, r2b=r2b,
dw=dw, pA=pA, kex=kex, spins_params=spins_params)
# Make nested list arrays of data. And return them.
values, errors, cpmg_frqs, missing, frqs, exp_types, relax_times,
offsets = self.return_r2eff_arrays()
@@ -138,6 +167,24 @@
@return: The numpy array structures of the R2eff/R1rho values,
errors, missing data, and corresponding Larmor frequencies. For each
structure, the first dimension corresponds to the experiment types, the
second the spins of a spin block, the third to the spectrometer field
strength, and the fourth is the dispersion points. For the Larmor frequency
structure, the fourth dimension is omitted. For R1rho-type data, an offset
dimension is inserted between the spectrometer field strength and the
dispersion points.
@rtype: lists of numpy float arrays, lists of numpy float
arrays, lists of numpy float arrays, numpy rank-2 int array
"""
+
+ # Unpack the parameter values.
+ # Initialise the post spin parameter indices.
+ end_index = []
+ # The spin and frequency dependent R2 parameters.
+ end_index.append(len(self.exp_type) * self.num_spins *
len(self.fields))
+ if self.model in [MODEL_B14_FULL, MODEL_CR72_FULL,
MODEL_NS_CPMG_2SITE_3D_FULL, MODEL_NS_CPMG_2SITE_STAR_FULL]:
+ end_index.append(2 * len(self.exp_type) * self.num_spins *
len(self.fields))
+ # The spin and dependent parameters (phi_ex, dw, padw2).
+ end_index.append(end_index[-1] + self.num_spins)
+
+ # Unpack the parameter values.
+ R20 = self.params[:end_index[1]].reshape(self.num_spins*2,
len(self.fields))
+ R20A = R20[::2].flatten()
+ R20B = R20[1::2].flatten()
+ dw = self.params[end_index[1]:end_index[2]]
+ pA = self.params[end_index[2]]
+ kex = self.params[end_index[2]+1]
# Initialise the data structures for the target function.
exp_types = []
@@ -201,19 +248,33 @@
for mi in range(len(self.fields)):
frq = self.fields[mi]
+
+ cpmg_frqs[ei][mi][oi] = self.points[mi]
+ #print len(self.points[mi]), self.points[mi]
+
for oi in range(len(self.offset)):
for di in range(len(self.points[mi])):
# The Larmor frequency for this spin (and that
of an attached proton for the MMQ models) and field strength (in MHz*2pi to
speed up the ppm to rad/s conversion).
frqs[ei][si][mi] = 2.0 * pi * frq / g1H * g15N *
1e-6
- cpmg_frqs[ei][mi][oi] = self.points[mi]
-
missing[ei][si][mi][oi].append(0)
+ # Calculate how the value should be, so chi2
gets zero.
+ # The R20 index.
+ r20_index = mi + si*len(self.fields)
+ # Convert dw from ppm to rad/s.
+ dw_frq = dw[si] * frqs[ei][si][mi]
+ r20a=R20A[r20_index]
+ r20b=R20B[r20_index]
+ back_calc = array([0.0])
+ r2eff_CR72(r20a=r20a, r20b=r20b, pA=pA,
dw=dw_frq, kex=kex, cpmg_frqs=cpmg_frqs[ei][mi][oi][di], back_calc=back_calc,
num_points=1)
+
# Values
- values[ei][si][mi][oi].append(self.value[mi][di])
+
#values[ei][si][mi][oi].append(self.value[mi][di])
+ values[ei][si][mi][oi].append(back_calc[0])
# The errors.
errors[ei][si][mi][oi].append(self.error[mi][di])
+ #print self.value[mi][di], self.error[mi][di]
# The relaxation times.
# Found.
@@ -234,7 +295,6 @@
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)
-
# Return the structures.
return values, errors, cpmg_frqs, missing, frqs, exp_types,
relax_times, offsets
@@ -365,14 +425,11 @@
"""
# Instantiate class
- C1 = Profile(num_spins=num_spins, model=model)
-
- # Assemble the parameter list.
- params = C1.assemble_param_vector(r2a=5.0, r2b=10.0, dw=3.0, pA=0.9,
kex=1000.0, spins_params=['r2a', 'r2b', 'dw', 'pA', 'kex'])
+ C1 = Profile(num_spins=num_spins, model=model, r2a=5.0, r2b=10.0,
dw=3.0, pA=0.9, kex=1000.0, spins_params=['r2a', 'r2b', 'dw', 'pA', 'kex'])
# Repeat the function call, to simulate minimisation.
for i in xrange(iter):
- chi2 = C1.calc(params)
+ chi2 = C1.calc(C1.params)
print("chi2 single:", chi2)
@@ -390,14 +447,11 @@
"""
# Instantiate class
- C1 = Profile(num_spins=num_spins, model=model)
-
- # Assemble the parameter list.
- params = C1.assemble_param_vector(r2a=5.0, r2b=10.0, dw=3.0, pA=0.9,
kex=1000.0, spins_params=['r2a', 'r2b', 'dw', 'pA', 'kex'])
+ C1 = Profile(num_spins=num_spins, model=model, r2a=5.0, r2b=10.0,
dw=3.0, pA=0.9, kex=1000.0, spins_params=['r2a', 'r2b', 'dw', 'pA', 'kex'])
# Repeat the function call, to simulate minimisation.
for i in xrange(iter):
- chi2 = C1.calc(params)
+ chi2 = C1.calc(C1.params)
print("chi2 cluster:", chi2)
@@ -406,14 +460,14 @@
main()
def test_reshape():
- C1 = Profile(num_spins=2, model=MODEL_CR72_FULL)
+ C1 = Profile(num_spins=1, model=MODEL_CR72_FULL, r2a=5.0, r2b=10.0,
dw=3.0, pA=0.9, kex=1000.0, spins_params=['r2a', 'r2b', 'dw', 'pA', 'kex'])
end_index = C1.model.end_index
#print("end_index:", end_index)
num_spins = C1.model.num_spins
#print("num_spins:", num_spins)
num_frq = C1.model.num_frq
#print("num_frq:", num_frq)
- params = C1.assemble_param_vector(r2a=5.0, r2b=10.0, dw=3.0, pA=0.9,
kex=1000.0, spins_params=['r2a', 'r2b', 'dw', 'pA', 'kex'])
+ params = C1.params
#print("params", params)
R20 = params[:end_index[1]].reshape(num_spins*2, num_frq)
@@ -435,8 +489,6 @@
r20b=R20B[r20_index]
print("r20a", r20a, "r20b", r20b)
- values = array(C1.model.values)
- values = array(values)
model = C1.calc(params)
print(model)
r23728 - /branches/disp_spin_speed/test_suite/shared_data/dispersion/profiling/profiling_cr72.py