I tried:
## Experiments
# Exp 1
sfrq_1 = 500.0*1E6
r20_key_1 = generate_r20_key(exp_type=EXP_TYPE_CPMG_SQ,
frq=sfrq_1)
time_T2_1 = 0.05
ncycs_1 = range(2,22,2)
# Here you define the direct R2eff errors (rad/s), as being added
or subtracted for the created R2eff point in the corresponding ncyc cpmg
frequence.
#r2eff_errs_1 = [0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05,
-0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05]
r2eff_errs_1 = [0.0] * len(ncycs_1)
exp_1 = [sfrq_1, time_T2_1, ncycs_1, r2eff_errs_1]
sfrq_2 = 600.0*1E6
r20_key_2 = generate_r20_key(exp_type=EXP_TYPE_CPMG_SQ,
frq=sfrq_2)
time_T2_2 = 0.06
ncycs_2 = range(2,22,2)
# Here you define the direct R2eff errors (rad/s), as being added
or subtracted for the created R2eff point in the corresponding ncyc cpmg
frequence.
#r2eff_errs_2 = [0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05,
-0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05]
r2eff_errs_2 = [0.0] * len(ncycs_2)
exp_2 = [sfrq_2, time_T2_2, ncycs_2, r2eff_errs_2]
sfrq_3 = 700.0*1E6
r20_key_3 = generate_r20_key(exp_type=EXP_TYPE_CPMG_SQ,
frq=sfrq_3)
time_T2_3 = 0.07
ncycs_3 = range(2,22,2)
# Here you define the direct R2eff errors (rad/s), as being added
or subtracted for the created R2eff point in the corresponding ncyc cpmg
frequence.
#r2eff_errs_2 = [0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05,
-0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05]
r2eff_errs_3 = [0.0] * len(ncycs_3)
exp_3 = [sfrq_3, time_T2_3, ncycs_3, r2eff_errs_3]
# Collect all exps
exps = [exp_1, exp_2, exp_3]
R20 = [10.1, 10.2, 10.3, 100.1, 100.2, 100.3, 20.1, 20.2, 20.3,
200.1, 200.2, 200.3, 30.1, 30.2, 30.3, 300.1, 300.2, 300.3, 40.1, 40.2,
40.3, 400.1, 400.2, 400.3]
#R20 = [10.1, 10.2, 10.3, 100.1, 100.2, 100.3, 20.1, 20.2, 20.3,
200.1, 200.2, 200.3]
dw_arr = [1.0, 2.0, 3.0, 4.0]
#dw_arr = [1.0, 2.0]
pA_arr = [0.9]
kex_arr = [1000.]
spins = [
['Ala', 1, 'N', {'r2a': {r20_key_1: R20[0], r20_key_2:
R20[1], r20_key_3: R20[2]}, 'r2b': {r20_key_1: R20[3], r20_key_2: R20[4],
r20_key_3: R20[5]}, 'kex': kex_arr[0], 'pA': pA_arr[0], 'dw': dw_arr[0]}],
['Ala', 2, 'N', {'r2a': {r20_key_1: R20[6], r20_key_2:
R20[7], r20_key_3: R20[8]}, 'r2b': {r20_key_1: R20[9], r20_key_2: R20[10],
r20_key_3: R20[11]}, 'kex': kex_arr[0], 'pA': pA_arr[0], 'dw': dw_arr[1]}],
['Ala', 3, 'N', {'r2a': {r20_key_1: R20[12], r20_key_2:
R20[13], r20_key_3: R20[14]}, 'r2b': {r20_key_1: R20[15], r20_key_2:
R20[16], r20_key_3: R20[17]}, 'kex': kex_arr[0], 'pA': pA_arr[0], 'dw':
dw_arr[2]}],
['Ala', 4, 'N', {'r2a': {r20_key_1: R20[18], r20_key_2:
R20[19], r20_key_3: R20[20]}, 'r2b': {r20_key_1: R20[21], r20_key_2:
R20[22], r20_key_3: R20[23]}, 'kex': kex_arr[0], 'pA': pA_arr[0], 'dw':
dw_arr[3]}],
]
------
------------
relax> grid_search(lower=[10.1, 10.2, 10.3, 100.1, 100.2, 100.3, 20.1,
20.2, 20.3, 200.1, 200.2, 200.3, 30.1, 30.2, 30.3, 300.1, 300.2, 300.3,
40.1, 40.2, 40.3, 400.1, 400.2, 400.3, 1.0, 2.0, 3.0, 4.0, 0.9, 1000.0],
upper=[10.1, 10.2, 10.3, 100.1, 100.2, 100.3, 20.1, 20.2, 20.3, 200.1,
200.2, 200.3, 30.1, 30.2, 30.3, 300.1, 300.2, 300.3, 40.1, 40.2, 40.3,
400.1, 400.2, 400.3, 1.0, 2.0, 3.0, 4.0, 0.9, 1000.0], inc=1,
constraints=True, verbosity=1)
Fitting to the spin block [':1@N', ':2@N', ':3@N', ':4@N']
----------------------------------------------------------
Unconstrained grid search size: 1 (constraints may decrease this size).
Grid search
~~~~~~~~~~~
Searching through 1 grid nodes.
Optimised parameter values:
r2a (SQ CPMG - 500.00000000 MHz) 10.100000000000000
r2a (SQ CPMG - 600.00000000 MHz) 10.199999999999999
r2a (SQ CPMG - 700.00000000 MHz) 10.300000000000001
r2b (SQ CPMG - 500.00000000 MHz) 100.099999999999994
r2b (SQ CPMG - 600.00000000 MHz) 100.199999999999989
r2b (SQ CPMG - 700.00000000 MHz) 100.299999999999997
r2a (SQ CPMG - 500.00000000 MHz) 20.100000000000001
r2a (SQ CPMG - 600.00000000 MHz) 20.199999999999999
r2a (SQ CPMG - 700.00000000 MHz) 20.300000000000004
r2b (SQ CPMG - 500.00000000 MHz) 200.099999999999966
r2b (SQ CPMG - 600.00000000 MHz) 200.199999999999989
r2b (SQ CPMG - 700.00000000 MHz) 200.300000000000011
r2a (SQ CPMG - 500.00000000 MHz) 30.100000000000001
r2a (SQ CPMG - 600.00000000 MHz) 30.199999999999999
r2a (SQ CPMG - 700.00000000 MHz) 30.300000000000004
r2b (SQ CPMG - 500.00000000 MHz) 300.100000000000023
r2b (SQ CPMG - 600.00000000 MHz) 300.199999999999989
r2b (SQ CPMG - 700.00000000 MHz) 300.300000000000011
r2a (SQ CPMG - 500.00000000 MHz) 40.099999999999994
r2a (SQ CPMG - 600.00000000 MHz) 40.200000000000003
r2a (SQ CPMG - 700.00000000 MHz) 40.299999999999997
r2b (SQ CPMG - 500.00000000 MHz) 400.100000000000023
r2b (SQ CPMG - 600.00000000 MHz) 400.199999999999932
r2b (SQ CPMG - 700.00000000 MHz) 400.300000000000011
dw 1.000000000000000
dw 2.000000000000000
dw 3.000000000000000
dw 4.000000000000000
pA 0.900000000000000
kex 1000.000000000000000
('CR72 full', 'Ala', ':1@N', 'r2a', 'SQ CPMG - 600.00000000 MHz', 10.2,
10.2)
('CR72 full', 'Ala', ':1@N', 'r2a', 'SQ CPMG - 500.00000000 MHz', 10.1,
10.1)
('CR72 full', 'Ala', ':1@N', 'r2a', 'SQ CPMG - 700.00000000 MHz', 10.3,
10.3)
('CR72 full', 'Ala', ':1@N', 'r2b', 'SQ CPMG - 600.00000000 MHz',
100.19999999999999, 100.2)
('CR72 full', 'Ala', ':1@N', 'r2b', 'SQ CPMG - 500.00000000 MHz', 100.1,
100.1)
('CR72 full', 'Ala', ':1@N', 'r2b', 'SQ CPMG - 700.00000000 MHz', 100.3,
100.3)
('CR72 full', 'Ala', ':2@N', 'r2a', 'SQ CPMG - 600.00000000 MHz', 20.2,
20.2)
('CR72 full', 'Ala', ':2@N', 'r2a', 'SQ CPMG - 500.00000000 MHz', 20.1,
20.1)
('CR72 full', 'Ala', ':2@N', 'r2a', 'SQ CPMG - 700.00000000 MHz',
20.300000000000004, 20.3)
('CR72 full', 'Ala', ':2@N', 'r2b', 'SQ CPMG - 600.00000000 MHz', 200.2,
200.2)
('CR72 full', 'Ala', ':2@N', 'r2b', 'SQ CPMG - 500.00000000 MHz',
200.09999999999997, 200.1)
('CR72 full', 'Ala', ':2@N', 'r2b', 'SQ CPMG - 700.00000000 MHz', 200.3,
200.3)
('CR72 full', 'Ala', ':3@N', 'r2a', 'SQ CPMG - 600.00000000 MHz', 30.2,
30.2)
('CR72 full', 'Ala', ':3@N', 'r2a', 'SQ CPMG - 500.00000000 MHz', 30.1,
30.1)
('CR72 full', 'Ala', ':3@N', 'r2a', 'SQ CPMG - 700.00000000 MHz',
30.300000000000004, 30.3)
('CR72 full', 'Ala', ':3@N', 'r2b', 'SQ CPMG - 600.00000000 MHz', 300.2,
300.2)
('CR72 full', 'Ala', ':3@N', 'r2b', 'SQ CPMG - 500.00000000 MHz', 300.1,
300.1)
('CR72 full', 'Ala', ':3@N', 'r2b', 'SQ CPMG - 700.00000000 MHz', 300.3,
300.3)
('CR72 full', 'Ala', ':4@N', 'r2a', 'SQ CPMG - 600.00000000 MHz', 40.2,
40.2)
('CR72 full', 'Ala', ':4@N', 'r2a', 'SQ CPMG - 500.00000000 MHz',
40.099999999999994, 40.1)
('CR72 full', 'Ala', ':4@N', 'r2a', 'SQ CPMG - 700.00000000 MHz', 40.3,
40.3)
('CR72 full', 'Ala', ':4@N', 'r2b', 'SQ CPMG - 600.00000000 MHz',
400.19999999999993, 400.2)
('CR72 full', 'Ala', ':4@N', 'r2b', 'SQ CPMG - 500.00000000 MHz', 400.1,
400.1)
('CR72 full', 'Ala', ':4@N', 'r2b', 'SQ CPMG - 700.00000000 MHz', 400.3,
400.3)
--------
It jumps over the target function?
This is a little weird?
Best
Troels
2014-06-06 15:02 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:
Hi,
Some more information is needed, as it's not possible to tell where
this stopped. A good idea would be to turn up the verbosity level to
see what minfx is doing. Did you call the grid_search user function?
Or the minimise user function? Did the grid search say something like
"Searching through 1 grid nodes"?
Regards,
Edward
On 6 June 2014 14:48, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx>
wrote:
Hi Edward.
When I try to do a grid search, it does not initalize func_CR72_full in
the
target function?
If I make
import sys
sys.exit()
It does not stop?
Only if I do a minimise, it stops?
It is in
specific_analyses/relax_disp/optimisation.py
line 745
Inserting
print model.func.im_func.__name__
gives func_CR72_full
How does it know how to unpack and calculate?
Best
Troels
2014-06-06 12:32 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:
Hi,
That sounds good. Maybe it's best to have the number of fields and
number of spins set to something different and not to 2? That way the
unpacking is stressed as much as possible and there cannot be a
accidental swap of the field and spin dimensions being unnoticed by
the test. This is not likely, but I've encountered enough weird and
supposedly impossible situations in the development of relax that it
would not surprise me.
Cheers,
Edward
On 6 June 2014 12:27, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx>
wrote:
Check.
I am generating R2eff data for 3 fields, and 3 spins, for full
model.
I will put the data in
test_suite/shared_data/dispersion/bug_22146_unpacking_r2a_r2b_cluster
Best
Troels
2014-06-06 12:11 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:
Hi,
Right, you have the 2 parameters in the self.num_spins*2 part. And
I
forgot about the parameters being different for each field. It
would
be good to then have a multi-field and multi-spin cluster system
test
to really make sure that relax operates correctly, especially with
the
data going into the target function and the subsequently unpacking
the
results into the relax data store. For example someone might
modify
the loop_parameters() function - this concept could be migrated
into
the specific API and converted into a common mechanism for all the
analysis types as it is quite powerful - and they may not know that
the change they just made broke code in the
target_functions.relax_disp module.
Cheers,
Edward
On 6 June 2014 12:05, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx>
wrote:
Hi Ed.
The implementations needs:
R20 =
params[:self.end_index[1]].reshape(self.num_spins*2,
self.num_frq)
R20A = R20[::2].flatten()
R20B = R20[1::2].flatten()
2014-06-06 11:55 GMT+02:00 Edward d'Auvergne
<edward@xxxxxxxxxxxxx>:
The different unpacking implementations can be tested with the
timeit
Python module to see which is fastest
(http://thread.gmane.org/gmane.science.nmr.relax.devel/5937/focus=6010).
Cheers,
Edward
On 6 June 2014 11:53, Edward d'Auvergne <edward@xxxxxxxxxxxxx>
wrote:
Hi,
In this case, I think 'num_frq' should be fixed to 2. This
dimension
corresponds to the parameters R20A and R20B so it is always
fixed
to
2.
Regards,
Edward
On 6 June 2014 11:51, Troels Emtekær Linnet
<tlinnet@xxxxxxxxxxxxx>
wrote:
Hi.
Another way is:
ml = params[:end_index[1]].reshape(num_spins*2, num_frq)
R20A = ml[::2].flatten()
R20B = ml[1::2].flatten()
Best
Troels
2014-06-06 11:39 GMT+02:00 Troels Emtekær Linnet
<tlinnet@xxxxxxxxxxxxx>:
There is no doubt that it is the unpacking of the R20A and
R20B
in
the
target function.
I was thinking of creating a function, which do the the
unpacking.
This unpacking function could then be tested with a unit
test?
What do you think?
Where should I position such a function?
Best
Troels
2014-06-06 11:26 GMT+02:00 Edward d'Auvergne
<edward@xxxxxxxxxxxxx>:
Hi Troels,
The best way to handle this is to first create a unit test
of
the
specific_analyses.relax_disp.parameters.disassemble_param_vector()
where the problem is likely to be most easily found. I
don't
understand how this could be a problem as the
assemble_param_vector()
and disassemble_param_vector() functions both call the same
loop_parameters() function for the ordering of the
parameter
values!
Maybe the problem is in the unpacking of the parameter
vector
in
the
target functions themselves, for example in the full B14
model:
def func_B14_full(self, params):
"""Target function for the Baldwin (2014) 2-site
exact
solution model for all time scales.
This assumes that pA > pB, and hence this must be
implemented
as a constraint.
@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.
R20A = params[:self.end_index[0]]
R20B = params[self.end_index[0]:self.end_index[1]]
dw = params[self.end_index[1]:self.end_index[2]]
pA = params[self.end_index[2]]
kex = params[self.end_index[2]+1]
# Calculate and return the chi-squared value.
return self.calc_B14_chi2(R20A=R20A, R20B=R20B,
dw=dw,
pA=pA,
kex=kex)
This R20A and R20B unpacking might be the failure point as
this
may
not match the loop_parameters() function - which it must!
In
any
case, having a unit or system test catch the problem would
be
very
useful for the stability of the dispersion analysis in
relax.
A code example might be useful:
R20_params = array([1, 2, 3, 4])
R20A, R20B = transpose(R20_params.reshape(2, 2)
print(R20A)
print(R20B)
You should see that R20A is [1, 3], and R20B is [2, 4].
This
is
how
the parameters are handled in the loop_parameters()
function
which
defines the parameter vector in all parts of relax. There
might
be a
quicker way to unpack the parameters, but such an idea
could be
used
for the target functions.
Cheers,
Edward
On 6 June 2014 11:08, Troels E. Linnet
<NO-REPLY.INVALID-ADDRESS@xxxxxxx>
wrote:
URL:
<http://gna.org/bugs/?22146>
Summary: Unpacking of R2A and R2B is
performed
wrong
for
clustered "full" dispersion models
Project: relax
Submitted by: tlinnet
Submitted on: Fri 06 Jun 2014 09:08:58 AM UTC
Category: relax's source code
Specific analysis category: Relaxation dispersion
Priority: 9 - Immediate
Severity: 4 - Important
Status: None
Assigned to: None
Originator Name:
Originator Email:
Open/Closed: Open
Release: Repository: trunk
Discussion Lock: Any
Operating System: All systems
_______________________________________________________
Details:
The unpacking of the R2A and R2B parameters in the "full"
model
is
performed
wrong.
This will happen performing a clustered analysis, using
one
of
the
"full"
models.
This bug affect all analysis performed running with a
"full"
model,
with
clustered residues.
The bug is located in the target function:
./target_functions/relax_disp.py
For all the "func_MODEL_full", the unpacking of:
R20A = params[:self.end_index[0]]
R20B = params[self.end_index[0]:self.end_index[1]]
This is wrong, since the "params" list, is ordered:
[spin, spin, spin, [dw], pA, kex], where spin =
[nr_frq*r2a,
nr_frq*r2b]
This ordering happens in:
./specific_analysis/relax_disp/parameters.py
in the loop_parameters.py
A possible solutions i shown below.
This alter the unpacking of the parameters.
An example of profiling_cr72.py is attached.
This can be downloaded, and run in base folder of relax:
./profiling_cr72.py .
This is with 3 frq, and 3 spins.
The current implementations would unpack:
('R20A', array([ 2., 2., 2., 4., 4., 4., 12.,
12.,
12.]),
9)
('R20B', array([ 14., 14., 14., 22., 22., 22., 24.,
24.,
24.]),
9)
R2A is 2, 12, 22 for the spins 0-3
R2B is, 4, 14, 24 for the spins 0-3
The suggested unpacking loop, unpacks to:
('R20A', array([ 2., 2., 2., 12., 12., 12., 22.,
22.,
22.]),
9)
('R20B', array([ 4., 4., 4., 14., 14., 14., 24.,
24.,
24.]),
9)
-------
from numpy import array, concatenate, delete, index_exp
import numpy
p = array([ 1.000000000000000e+01,
1.000000000000000e+01,
1.100000000000000e+01
, 1.100000000000000e+01, 1.000000000000000e+01,
1.000000000000000e+01
, 1.100000000000000e+01, 1.100000000000000e+01,
1.000000000000000e+00
, 1.000000000000000e+00, 9.000000000000000e-01,
1.000000000000000e+03])
e = [4, 8, 10]
# Now
r2a = p[:e[0]]
print r2a
r2b = p[e[0]:e[1]]
print r2b
dw = p[e[1]:e[2]]
print dw
pA = p[e[2]]
print pA
kex = p[e[2]+1]
print kex
print "new"
ns = 2
nf = 2
ml = p[:e[1]]
R20A = array([])
R20B = array([])
for i in range(0, ns):
# Array sorted per [spin, spin, spin], where spin =
[nr_frq*r2a,
nr_frq*r2b]
spin_AB = ml[:nf*2]
ml = delete(ml, numpy.s_[:nf*2])
R20A = concatenate([R20A, spin_AB[:nf] ])
R20B = concatenate([R20B, spin_AB[nf:] ])
print R20A
print R20B
print dw
print pA
print kex
_______________________________________________________
File Attachments:
-------------------------------------------------------
Date: Fri 06 Jun 2014 09:08:58 AM UTC Name:
profiling_cr72.py
Size:
17kB
By: tlinnet
<http://gna.org/bugs/download.php?file_id=20938>
_______________________________________________________
Reply to this item at:
<http://gna.org/bugs/?22146>
_______________________________________________
Message sent via/by Gna!
http://gna.org/
_______________________________________________
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
Re: [bug #22146] Unpacking of R2A and R2B is performed wrong for clustered "full" dispersion models