Hi Troels, Just so you know, the grid search is not used in the Monte Carlo simulations. The starting position for this method is always the real optimised parameters. So unfortunately this will not be of benefit for your speed aims.
Regards, Edward On 21 August 2014 22:57, <tlin...@nmr-relax.com> wrote: > Author: tlinnet > Date: Thu Aug 21 22:57:45 2014 > New Revision: 25189 > > URL: http://svn.gna.org/viewcvs/relax?rev=25189&view=rev > Log: > Modified intermediate systemtest: > Relax_disp.test_bug_9999_slow_r1rho_r2eff_error_with_mc > > to see if the initial Grid Search for 'i0' and 'R2eff' estimation can be > skipped. > > This is done by converting the exponential curve, to a linear curve, and > calculate > the best parameters by a line of best fit by least squares. > > This seems like a promising method as an initial estimator of 'i0' and > 'r2eff'. > > For 500 to 2000 Monte-Carlo simulations, this could have an influence on the > timings, > as all grid searchs could be skipped. > > 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=25189&r1=25188&r2=25189&view=diff > ============================================================================== > --- trunk/test_suite/system_tests/relax_disp.py (original) > +++ trunk/test_suite/system_tests/relax_disp.py Thu Aug 21 22:57:45 2014 > @@ -23,7 +23,7 @@ > > # Python module imports. > from os import F_OK, access, getcwd, path, sep > -from numpy import array, median > +from numpy import array, exp, median, log, sum > import re, math > from tempfile import mkdtemp > > @@ -35,7 +35,7 @@ > from lib.io import get_file_path > from pipe_control.mol_res_spin import generate_spin_string, return_spin, > spin_loop > from specific_analyses.relax_disp.checks import check_missing_r1 > -from specific_analyses.relax_disp.data import generate_r20_key, > get_curve_type, has_exponential_exp_type, has_r1rho_exp_type, loop_exp_frq, > loop_exp_frq_offset_point, loop_exp_frq_offset_point_time, > return_grace_file_name_ini, return_param_key_from_data > +from specific_analyses.relax_disp.data import average_intensity, > generate_r20_key, get_curve_type, has_exponential_exp_type, > has_r1rho_exp_type, loop_exp_frq, loop_exp_frq_offset_point, > loop_exp_frq_offset_point_time, loop_time, return_grace_file_name_ini, > return_param_key_from_data > from specific_analyses.relax_disp.data import INTERPOLATE_DISP, > INTERPOLATE_OFFSET, X_AXIS_DISP, X_AXIS_W_EFF, X_AXIS_THETA, Y_AXIS_R2_R1RHO, > Y_AXIS_R2_EFF > from specific_analyses.relax_disp.model import models_info, nesting_param > from specific_analyses.relax_disp.variables import EXP_TYPE_CPMG_DQ, > EXP_TYPE_CPMG_MQ, EXP_TYPE_CPMG_PROTON_MQ, EXP_TYPE_CPMG_PROTON_SQ, > EXP_TYPE_CPMG_SQ, EXP_TYPE_CPMG_ZQ, EXP_TYPE_LIST, EXP_TYPE_R1RHO, > MODEL_B14_FULL, MODEL_CR72, MODEL_CR72_FULL, MODEL_DPL94, MODEL_IT99, > MODEL_LIST_ANALYTIC_CPMG, MODEL_LIST_FULL, MODEL_LIST_NUMERIC_CPMG, > MODEL_LM63, MODEL_M61, MODEL_M61B, MODEL_MP05, MODEL_NOREX, > MODEL_NS_CPMG_2SITE_3D_FULL, MODEL_NS_CPMG_2SITE_EXPANDED, > MODEL_NS_CPMG_2SITE_STAR_FULL, MODEL_NS_R1RHO_2SITE, MODEL_NS_R1RHO_3SITE, > MODEL_NS_R1RHO_3SITE_LINEAR, MODEL_PARAMS, MODEL_R2EFF, MODEL_TP02, > MODEL_TAP03 > @@ -158,8 +158,8 @@ > > # Now loop over the experiments, to set the variables in relax. > exp_ids = [] > - for exp in exps: > - sfrq, time_T2, ncycs, r2eff_errs = exp > + for exp_i in exps: > + sfrq, time_T2, ncycs, r2eff_errs = exp_i > exp_id = 'CPMG_%3.1f' % (sfrq/1E6) > exp_ids.append(exp_id) > > @@ -193,8 +193,8 @@ > # Now read data in. > for exp_type, frq, ei, mi in loop_exp_frq(return_indices=True): > exp_id = exp_ids[mi] > - exp = exps[mi] > - sfrq, time_T2, ncycs, r2eff_errs = exp > + exp_i = exps[mi] > + sfrq, time_T2, ncycs, r2eff_errs = exp_i > > # Then loop over the spins. > for res_name, res_num, spin_name, params in spins: > @@ -1363,26 +1363,15 @@ > # Read data. > self.interpreter.results.read(prev_data_path + sep + 'results') > > - # Get initial offset, point, time > - for exp_type, frq, offset, point, time, ei, mi, oi, di, ti in > loop_exp_frq_offset_point_time(return_indices=True): > - offset_i = offset > - point_i = point > - time_i = time > - break > - > # Now count number > - graph_nr = 0 > - for exp_type, frq, offset, point, time, ei, mi, oi, di, ti in > loop_exp_frq_offset_point_time(return_indices=True): > - print(exp_type, frq, offset, point, time) > - > - # If different, count 1 graph. > - if offset != offset_i or point != point_i: > - offset_i = offset > - point_i = point > - graph_nr += 1 > - print("Graph %i complete\n" % graph_nr) > - > - print(graph_nr + 1) > + graph_nr = 1 > + for exp_type, frq, offset, point in > loop_exp_frq_offset_point(return_indices=False): > + print("\nGraph nr %i" % graph_nr) > + for time in loop_time(exp_type=exp_type, frq=frq, offset=offset, > point=point): > + print(exp_type, frq, offset, point, time) > + graph_nr += 1 > + > + ## Possibly do an error analysis. > > # Check if intensity errors have already been calculated by the user. > precalc = True > @@ -1425,13 +1414,111 @@ > model = 'R2eff' > self.interpreter.relax_disp.select_model(model) > > + for spin, spin_id in spin_loop(return_id=True, skip_desel=True): > + #delattr(spin, 'r2eff') > + #delattr(spin, 'r2eff_err') > + #delattr(spin, 'i0') > + #delattr(spin, 'i0_err') > + setattr(spin, 'r2eff', {}) > + setattr(spin, 'r2eff_err', {}) > + setattr(spin, 'i0', {}) > + setattr(spin, 'i0_err', {}) > + > # Do Grid Search > - self.interpreter.minimise.grid_search(lower=None, upper=None, > inc=11, constraints=True, verbosity=1) > + self.interpreter.minimise.grid_search(lower=None, upper=None, > inc=21, constraints=True, verbosity=1) > + > + # Start dic. > + my_dic = {} > + > + # Loop over each spectrometer frequency and dispersion point. > + for cur_spin, mol_name, resi, resn, spin_id in > spin_loop(full_info=True, return_id=True, skip_desel=True): > + # Add key to dic. > + my_dic[spin_id] = {} > + > + # Generate spin string. > + spin_string = generate_spin_string(spin=cur_spin, > mol_name=mol_name, res_num=resi, res_name=resn) > + > + # Loop over the parameters. > + #print("Grid optimised parameters for spin: %s" % (spin_string)) > + > + for exp_type, frq, offset, point in loop_exp_frq_offset_point(): > + # 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][param_key] = {} > + > + # Get the value. > + R2eff_value = getattr(cur_spin, 'r2eff')[param_key] > + i0_value = getattr(cur_spin, 'i0')[param_key] > + > + # Save to dic. > + my_dic[spin_id][param_key]['R2eff_value_grid'] = R2eff_value > + my_dic[spin_id][param_key]['i0_value_grid'] = i0_value > + > + ## Now try do a line of best fit by least squares. > + # The peak intensities, errors and times. > + 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, > sim_index=None)) > + errors.append(average_intensity(spin=cur_spin, > exp_type=exp_type, frq=frq, offset=offset, point=point, time=time, > error=True)) > + times.append(time) > + > + # y= A exp(x * k) > + # w[i] = ln(y[i]) > + # int[i] = i0 * exp( - times[i] * r2eff); > + w = log(array(values)) > + x = - array(times) > + n = len(times) > + > + b = (sum(x*w) - 1./n * sum(x) * sum(w) ) / ( sum(x**2) - > 1./n * (sum(x))**2 ) > + a = 1./n * sum(w) - b * 1./n * sum(x) > + R2eff_est = b > + i0_est = exp(a) > + > + my_dic[spin_id][param_key]['R2eff_est'] = R2eff_est > + my_dic[spin_id][param_key]['i0_est'] = i0_est > + > + # Print value. > + #print("%-10s %-6s %-6s %3.1f : %3.1f" % ("Parameter:", > 'R2eff', "Value : Estimated:", R2eff_value, R2eff_est)) > + #print("%-10s %-6s %-6s %3.1f : %3.1f" % ("Parameter:", > 'i0', "Value: Estimated:", i0_value, i0_est)) > + > > # Do minimisation. > - set_func_tol = 1e-11 > - set_max_iter = 10000 > + set_func_tol = 1e-25 > + set_max_iter = int(1e7) > self.interpreter.minimise.execute(min_algor='simplex', > func_tol=set_func_tol, max_iter=set_max_iter, constraints=True, scaling=True, > verbosity=1) > + > + # Loop over each spectrometer frequency and dispersion point. > + for cur_spin, mol_name, resi, resn, spin_id in > spin_loop(full_info=True, return_id=True, skip_desel=True): > + # Generate spin string. > + spin_string = generate_spin_string(spin=cur_spin, > mol_name=mol_name, res_num=resi, res_name=resn) > + > + # Loop over the parameters. > + print("Optimised parameters for spin: %s" % (spin_string)) > + > + for exp_type, frq, offset, point in loop_exp_frq_offset_point(): > + # Generate the param_key. > + param_key = return_param_key_from_data(exp_type=exp_type, > frq=frq, offset=offset, point=point) > + > + # Get the value. > + R2eff_value = getattr(cur_spin, 'r2eff')[param_key] > + i0_value = getattr(cur_spin, 'i0')[param_key] > + > + # Extract from dic. > + R2eff_value_grid = > my_dic[spin_id][param_key]['R2eff_value_grid'] > + i0_value_grid = my_dic[spin_id][param_key]['i0_value_grid'] > + R2eff_est = my_dic[spin_id][param_key]['R2eff_est'] > + i0_est = my_dic[spin_id][param_key]['i0_est'] > + > + # Print value. > + #print("%-10s %-6s %-6s %3.1f : %3.1f" % ("Parameter:", > 'R2eff', "Value : Estimated:", R2eff_value, R2eff_est)) > + #print("%-10s %-6s %-6s %3.1f : %3.1f" % ("Parameter:", > 'i0', "Value: Estimated:", i0_value, i0_est)) > + > + print("%-10s %-6s %-6s %3.1f : %3.1f: %3.1f" % > ("Parameter:", 'R2eff', "Grid : Min : Estimated:", R2eff_value_grid, > R2eff_value, R2eff_est)) > + print("%-10s %-6s %-6s %3.1f : %3.1f: %3.1f" % > ("Parameter:", 'i0', "Grid : Min : Estimated:", i0_value_grid, i0_value, > i0_est)) > > > def test_check_missing_r1(self): > > > _______________________________________________ > relax (http://www.nmr-relax.com) > > This is the relax-commits mailing list > relax-comm...@gna.org > > 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@gna.org 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
