Author: tlinnet Date: Thu Aug 28 12:34:29 2014 New Revision: 25377 URL: http://svn.gna.org/viewcvs/relax?rev=25377&view=rev Log: Implemented system test Relax_disp.verify_estimate_r2eff_err_compare_mc for testing R2eff error as function of Monte Carlo simulation. Note, since the name does not start with "test", but with "verify", this test will not be issued in the system test suite. -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 50 0.038 0.037 0.048 0.041 0.048 0.036 0.041 0.036 0.030 0.035 0.014 0.018 0.010 0.007 sum= 0.438 100 0.036 0.039 0.040 0.044 0.043 0.039 0.041 0.040 0.033 0.037 0.013 0.018 0.008 0.007 sum= 0.440 150 0.035 0.040 0.045 0.046 0.043 0.041 0.041 0.035 0.034 0.035 0.013 0.017 0.008 0.008 sum= 0.442 200 0.034 0.040 0.046 0.047 0.042 0.037 0.039 0.035 0.033 0.039 0.012 0.019 0.009 0.008 sum= 0.440 250 0.036 0.039 0.042 0.042 0.043 0.039 0.043 0.035 0.033 0.037 0.013 0.016 0.009 0.007 sum= 0.436 300 0.034 0.040 0.047 0.047 0.043 0.037 0.042 0.036 0.033 0.039 0.013 0.018 0.009 0.008 sum= 0.446 350 0.034 0.041 0.045 0.046 0.043 0.037 0.038 0.036 0.036 0.037 0.013 0.018 0.009 0.008 sum= 0.441 400 0.036 0.037 0.043 0.047 0.044 0.038 0.043 0.038 0.035 0.037 0.014 0.018 0.009 0.008 sum= 0.448 450 0.034 0.040 0.044 0.045 0.044 0.038 0.041 0.035 0.034 0.039 0.012 0.018 0.009 0.008 sum= 0.442 0.9 ++-----+------+------+------+------+-----+------+------+------+-----++ + A + + R2eff error as function of MC number **A*** + 0.8 ++ * ++ | * | 0.7 ++ * ++ | * | 0.6 ++ * ++ | * | 0.5 ++ * ++ | * | | * A******A******A******A*****A******A******A******A******A 0.4 ++ * * ++ | * * | 0.3 ++ * * ++ | * * | 0.2 ++ * * ++ | * * | 0.1 ++ * * ++ + ** + + + + + + + + + 0 ++-----A------+------+------+------+-----+------+------+------+-----++ -50 0 50 100 150 200 250 300 350 400 450 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=25377&r1=25376&r2=25377&view=diff ============================================================================== --- trunk/test_suite/system_tests/relax_disp.py (original) +++ trunk/test_suite/system_tests/relax_disp.py Thu Aug 28 12:34:29 2014 @@ -25,7 +25,8 @@ from os import F_OK, access, getcwd, path, sep from numpy import array, exp, median, log, save, sum, zeros import re, math -from tempfile import mkdtemp +from tempfile import mkdtemp, NamedTemporaryFile + # relax module imports. from auto_analyses import relax_disp @@ -7571,6 +7572,212 @@ w_eff_file.close() + def verify_estimate_r2eff_err_compare_mc(self): + """Test the user function for estimating R2eff errors from exponential curve fitting, and compare it with Monte-Carlo simulations. + + This follows Task 7822. + U{task #7822<https://gna.org/task/index.php?7822>}: Implement user function to estimate R2eff and associated errors for exponential curve fitting. + + This uses the data from Kjaergaard's paper at U{DOI: 10.1021/bi4001062<http://dx.doi.org/10.1021/bi4001062>}. + Optimisation of the Kjaergaard et al., 2013 Off-resonance R1rho relaxation dispersion experiments using the 'DPL' model. + """ + + # Load the data. + data_path = status.install_path + sep+'test_suite'+sep+'shared_data'+sep+'dispersion'+sep+'Kjaergaard_et_al_2013'+sep + + # Set pipe name, bundle and type. + pipe_name = 'base pipe' + pipe_bundle = 'relax_disp' + pipe_type = 'relax_disp' + + # Create the data pipe. + self.interpreter.pipe.create(pipe_name=pipe_name, bundle=pipe_bundle, pipe_type=pipe_type) + + file = data_path + '1_setup_r1rho_GUI.py' + self.interpreter.script(file=file, dir=None) + + # Deselect all spins. + self.interpreter.deselect.spin(spin_id=':1-100', change_all=False) + + # Select one spin. + self.interpreter.select.spin(spin_id=':52@N', change_all=False) + + # Set the model. + self.interpreter.relax_disp.select_model(MODEL_R2EFF) + + # Check if intensity errors have already been calculated. + check_intensity_errors() + + # Do a grid search. + self.interpreter.minimise.grid_search(lower=None, upper=None, inc=11, constraints=True, verbosity=1) + + # Set algorithm. + min_algor = 'Newton' + constraints = False + + # Minimise. + self.interpreter.minimise.execute(min_algor=min_algor, constraints=constraints, verbosity=1) + + # Loop over old err attributes. + err_attr_list = ['r2eff_err', 'i0_err'] + + for cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True): + # Loop over old err attributes. + for err_attr in err_attr_list: + if hasattr(cur_spin, err_attr): + delattr(cur_spin, err_attr) + + # Estimate R2eff errors. + self.interpreter.relax_disp.r2eff_err_estimate() + + # Collect the estimation data. + my_dic = {} + param_key_list = [] + est_keys = [] + est_key = '-1' + est_keys.append(est_key) + spin_id_list = [] + + 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] = {} + + # Add key for estimate. + my_dic[spin_id][est_key] = {} + + # Add spin key to list. + spin_id_list.append(spin_id) + + # Generate spin string. + spin_string = generate_spin_string(spin=cur_spin, mol_name=mol_name, res_num=resi, res_name=resn) + + 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) + + # Append key. + param_key_list.append(param_key) + + # 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, 500, 50) + + 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'] + + # Loop over the Monte Carlo simulations: + for number in mc_number_list: + # First delete old simulations. + for cur_spin, mol_name, resi, resn, spin_id in spin_loop(full_info=True, return_id=True, skip_desel=True): + # Loop over old err attributes. + for err_attr in err_attr_list: + if hasattr(cur_spin, err_attr): + delattr(cur_spin, err_attr) + + # Loop over the simulated attributes. + for sim_attr in sim_attr_list: + if hasattr(cur_spin, sim_attr): + delattr(cur_spin, sim_attr) + + self.interpreter.monte_carlo.setup(number=number) + self.interpreter.monte_carlo.create_data() + self.interpreter.monte_carlo.initial_values() + self.interpreter.minimise.execute(min_algor=min_algor, constraints=constraints) + self.interpreter.eliminate() + self.interpreter.monte_carlo.error_analysis() + + est_key = '%i'%number + 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 + + # Set what to extract. + err_attr = err_attr_list[0] + + # Define list with text. + text_list = [] + + # Now loop through the data. + for spin_id in spin_id_list: + for est_key in est_keys: + # Define list to pickup data. + r2eff_err_list = [] + + for param_key in param_key_list: + # Get the value. + r2eff_err = my_dic[spin_id][est_key][param_key][err_attr] + + # Add to list. + r2eff_err_list.append(r2eff_err) + + # Sum the list + sum_array = sum(array(r2eff_err_list)) + + # Join floats to string. + r2eff_err_str = " ".join(format(x, "2.3f") for x in r2eff_err_list) + + # Define print string. + text = "%8s %s sum= %2.3f" % (est_key, r2eff_err_str, sum_array) + text_list.append(text) + + + # Now print. + filepath = NamedTemporaryFile(delete=False).name + # Open the files for testing. + w_file = open(filepath, 'w') + + print("Printing the estimated R2eff error as function of estimation from Co-variance and number of Monte-Carlo simulations.") + + for text in text_list: + # Print. + print(text) + + # Write to file. + w_file.write(text+"\n") + + # Close files + w_file.close() + + print("Filepath is: %s"%filepath) + print("Start 'gnuplot' and write:") + print("set term dumb") + print("plot '%s' using 1:17 title 'R2eff error as function of MC number' w linespoints "%filepath) + + def verify_r1rho_kjaergaard_missing_r1(self, models=None, result_dir_name=None, do_assert=True): """Verification of test_r1rho_kjaergaard_missing_r1."""