Author: tlinnet
Date: Sun Aug 31 17:26:48 2014
New Revision: 25485
URL: http://svn.gna.org/viewcvs/relax?rev=25485&view=rev
Log:
Modified systemtest Relax_disp.verify_estimate_r2eff_err_compare_mc to
include boot strapping method.
This shows there is excellent agreement between boot-strapping and estimation
of errors from Co-variance, while
relax Monte-Carlo simulations are very much different.
Boot strapping is the "-2":
-2 0.070 0.085 0.087 0.095 0.086 0.076 0.087 0.072 0.069 0.077 0.025
0.035 0.018 0.015 sum= 0.897
-1 0.069 0.081 0.085 0.092 0.085 0.074 0.083 0.069 0.066 0.074 0.025
0.035 0.018 0.016 sum= 0.874
0 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
0.000 0.000 0.000 sum= 0.000
400 0.034 0.043 0.044 0.049 0.046 0.037 0.042 0.035 0.031 0.039 0.014
0.018 0.009 0.008 sum= 0.450
800 0.032 0.040 0.041 0.046 0.042 0.036 0.040 0.034 0.034 0.037 0.013
0.018 0.009 0.008 sum= 0.431
1200 0.033 0.041 0.042 0.046 0.043 0.037 0.042 0.036 0.034 0.038 0.012
0.018 0.009 0.008 sum= 0.439
1600 0.036 0.041 0.042 0.047 0.043 0.038 0.042 0.035 0.035 0.037 0.013
0.018 0.009 0.008 sum= 0.443
2000 0.034 0.042 0.042 0.046 0.042 0.036 0.043 0.035 0.034 0.037 0.013
0.017 0.009 0.008 sum= 0.438
task #7822(https://gna.org/task/index.php?7822): Implement user function to
estimate R2eff and associated errors for exponential curve fitting.
Modified:
trunk/test_suite/system_tests/relax_disp.py
Modified: trunk/test_suite/system_tests/relax_disp.py
URL:
http://svn.gna.org/viewcvs/relax/trunk/test_suite/system_tests/relax_disp.py?rev=25485&r1=25484&r2=25485&view=diff
==============================================================================
--- trunk/test_suite/system_tests/relax_disp.py (original)
+++ trunk/test_suite/system_tests/relax_disp.py Sun Aug 31 17:26:48 2014
@@ -24,6 +24,8 @@
# Python module imports.
from os import F_OK, access, getcwd, path, sep
from numpy import array, asarray, exp, median, inf, log, save, std, sum,
zeros
+from minfx.generic import generic_minimise
+from random import gauss
import re, math
from tempfile import mkdtemp, NamedTemporaryFile
@@ -3125,9 +3127,7 @@
# Do boot strapping ?
do_boot = True
if do_boot:
- from minfx.generic import generic_minimise
- from random import gauss
- min_algor = 'BFGS'
+ min_algor = 'Newton'
min_options = ()
sim_boot = 2000
scaling_list = [1.0, 1.0]
@@ -3232,7 +3232,7 @@
x0 = [r2eff, i0]
setup(num_params=len(x0), num_times=len(times),
values=I_err, sd=errors, relax_times=times, scaling_matrix=scaling_list)
- params_minfx_sim_j, chi2_minfx_sim_j, iter_count,
f_count, g_count, h_count, warning = generic_minimise(func=func, dfunc=dfunc,
args=(), x0=x0, min_algor=min_algor, min_options=min_options,
full_output=True, print_flag=0)
+ params_minfx_sim_j, chi2_minfx_sim_j, iter_count,
f_count, g_count, h_count, warning = generic_minimise(func=func, dfunc=dfunc,
d2func=d2func, args=(), x0=x0, min_algor=min_algor, min_options=min_options,
full_output=True, print_flag=0)
R_m_sim_j, I0_m_sim_j = params_minfx_sim_j
R_m_sim_l.append(R_m_sim_j)
I0_m_sim_l.append(I0_m_sim_j)
@@ -7995,6 +7995,9 @@
# Set algorithm.
min_algor = 'Newton'
constraints = False
+ min_options = ()
+ sim_boot = 2000
+ scaling_list = [1.0, 1.0]
# Minimise.
self.interpreter.minimise.execute(min_algor=min_algor,
constraints=constraints, verbosity=1)
@@ -8008,10 +8011,7 @@
if hasattr(cur_spin, err_attr):
delattr(cur_spin, err_attr)
- # Estimate R2eff errors.
- self.interpreter.relax_disp.r2eff_err_estimate(chi2_jacobian=False)
-
- # Collect the estimation data.
+ # Collect the estimation data from boot.
my_dic = {}
param_key_list = []
est_keys = []
@@ -8038,6 +8038,69 @@
# Append key.
param_key_list.append(param_key)
+
+ # Add key to dic.
+ my_dic[spin_id][est_key][param_key] = {}
+
+ values = []
+ errors = []
+ times = []
+ for time in loop_time(exp_type=exp_type, frq=frq,
offset=offset, point=point):
+ values.append(average_intensity(spin=cur_spin,
exp_type=exp_type, frq=frq, offset=offset, point=point, time=time))
+ errors.append(average_intensity(spin=cur_spin,
exp_type=exp_type, frq=frq, offset=offset, point=point, time=time,
error=True))
+ times.append(time)
+
+ # Convert to numpy array.
+ values = asarray(values)
+ errors = asarray(errors)
+ times = asarray(times)
+
+ r2eff = getattr(cur_spin, 'r2eff')[param_key]
+ i0 = getattr(cur_spin, 'i0')[param_key]
+
+ R_m_sim_l = []
+ I0_m_sim_l = []
+ for j in range(sim_boot):
+ if j in range(0, 100000, 100):
+ print("Simulation %i"%j)
+ # Start minimisation.
+
+ # Produce errors
+ I_err = []
+ for j, error in enumerate(errors):
+ I_error = gauss(values[j], error)
+ I_err.append(I_error)
+ # Convert to numpy array.
+ I_err = asarray(I_err)
+
+ x0 = [r2eff, i0]
+ setup(num_params=len(x0), num_times=len(times),
values=I_err, sd=errors, relax_times=times, scaling_matrix=scaling_list)
+
+ params_minfx_sim_j, chi2_minfx_sim_j, iter_count,
f_count, g_count, h_count, warning = generic_minimise(func=func, dfunc=dfunc,
d2func=d2func, args=(), x0=x0, min_algor=min_algor, min_options=min_options,
full_output=True, print_flag=0)
+ R_m_sim_j, I0_m_sim_j = params_minfx_sim_j
+ R_m_sim_l.append(R_m_sim_j)
+ I0_m_sim_l.append(I0_m_sim_j)
+
+ # Get stats on distribution.
+ sigma_R_sim = std(asarray(R_m_sim_l), ddof=1)
+ sigma_I0_sim = std(asarray(I0_m_sim_l), ddof=1)
+ my_dic[spin_id][est_key][param_key]['r2eff_err'] =
sigma_R_sim
+ my_dic[spin_id][est_key][param_key]['i0_err'] = sigma_I0_sim
+
+ # Estimate R2eff errors.
+ self.interpreter.relax_disp.r2eff_err_estimate()
+
+ est_key = '-1'
+ est_keys.append(est_key)
+
+ # Collect data.
+ for cur_spin, mol_name, resi, resn, spin_id in
spin_loop(full_info=True, return_id=True, skip_desel=True):
+ # Add key for estimate.
+ my_dic[spin_id][est_key] = {}
+
+ for exp_type, frq, offset, point, ei, mi, oi, di in
loop_exp_frq_offset_point(return_indices=True):
+ # Generate the param_key.
+ param_key = return_param_key_from_data(exp_type=exp_type,
frq=frq, offset=offset, point=point)
# Add key to dic.
my_dic[spin_id][est_key][param_key] = {}
@@ -8054,38 +8117,8 @@
my_dic[spin_id][est_key][param_key][err_attr] =
get_err_attr
- # Estimate R2eff errors from Chi2 Jacobian.
- self.interpreter.relax_disp.r2eff_err_estimate(chi2_jacobian=True)
-
- est_key = '-1'
- est_keys.append(est_key)
-
- # Collect data.
- for cur_spin, mol_name, resi, resn, spin_id in
spin_loop(full_info=True, return_id=True, skip_desel=True):
- # Add key for estimate.
- my_dic[spin_id][est_key] = {}
-
- for exp_type, frq, offset, point, ei, mi, oi, di in
loop_exp_frq_offset_point(return_indices=True):
- # Generate the param_key.
- param_key = return_param_key_from_data(exp_type=exp_type,
frq=frq, offset=offset, point=point)
-
- # Add key to dic.
- my_dic[spin_id][est_key][param_key] = {}
-
- # Get the value.
- # Loop over err attributes.
- for err_attr in err_attr_list:
- if hasattr(cur_spin, err_attr):
- get_err_attr = getattr(cur_spin, err_attr)[param_key]
- else:
- get_err_attr = 0.0
-
- # Save to dic.
- my_dic[spin_id][est_key][param_key][err_attr] =
get_err_attr
-
-
# Make Carlo Simulations number
- mc_number_list = range(0, 50, 10)
+ mc_number_list = range(0, 2400, 400)
sim_attr_list = ['chi2_sim', 'f_count_sim', 'g_count_sim',
'h_count_sim', 'i0_sim', 'iter_sim', 'peak_intensity_sim', 'r2eff_sim',
'select_sim', 'warning_sim']
r25485 - /trunk/test_suite/system_tests/relax_disp.py