Author: tlinnet Date: Fri Sep 26 17:29:49 2014 New Revision: 26066 URL: http://svn.gna.org/viewcvs/relax?rev=26066&view=rev Log: Implemented writing and plotting of statistics for individual and clustered fitting, compararing to full NI. Task #7826 (https://gna.org/task/index.php?7826): Write an python class for the repeated analysis of dispersion data. Modified: trunk/auto_analyses/relax_disp_repeat_cpmg.py trunk/test_suite/system_tests/relax_disp.py Modified: trunk/auto_analyses/relax_disp_repeat_cpmg.py URL: http://svn.gna.org/viewcvs/relax/trunk/auto_analyses/relax_disp_repeat_cpmg.py?rev=26066&r1=26065&r2=26066&view=diff ============================================================================== --- trunk/auto_analyses/relax_disp_repeat_cpmg.py (original) +++ trunk/auto_analyses/relax_disp_repeat_cpmg.py Fri Sep 26 17:29:49 2014 @@ -1958,6 +1958,12 @@ res_dic['selection'] = selection params_list = min_dic_ref[str(glob_ini_ref)]['params']['params_list'] res_dic['params_list'] = params_list + analysis = min_dic_ref['analysis'] + res_dic['analysis'] = analysis + + # Get the current method + method_cur = min_dic['method'] + res_dic[method_cur] = {} # Loop over params for j, param in enumerate(params_list): @@ -1967,17 +1973,14 @@ param_arr_ref = min_dic_ref[str(glob_ini_ref)]['params'][param] res_dic[param]['param_arr_ref'] = param_arr_ref - # Get the current method - method_cur = min_dic['method'] - - res_dic[param][method_cur] = {} - res_dic[param][method_cur]['method'] = method_cur - res_dic[param][method_cur]['sampling_sparseness'] = [] - res_dic[param][method_cur]['glob_ini'] = [] + res_dic[method_cur][param] = {} + res_dic[method_cur][param]['method'] = method_cur + res_dic[method_cur][param]['sampling_sparseness'] = [] + res_dic[method_cur][param]['glob_ini'] = [] # Other stats. - res_dic[param][method_cur]['r_xy'] = [] - res_dic[param][method_cur]['a'] = [] + res_dic[method_cur][param]['r_xy'] = [] + res_dic[method_cur][param]['a'] = [] # Now loop over glob_ini: for glob_ini in list_glob_ini: @@ -1995,13 +1998,13 @@ # Store x sampling_sparseness = float(glob_ini) / float(glob_ini_ref) * 100. - res_dic[param][method_cur]['sampling_sparseness'].append(sampling_sparseness) - res_dic[param][method_cur]['glob_ini'].append(glob_ini) + res_dic[method_cur][param]['sampling_sparseness'].append(sampling_sparseness) + res_dic[method_cur][param]['glob_ini'].append(glob_ini) # Store to result dic. - res_dic[param][method_cur][str(glob_ini)] = {} - res_dic[param][method_cur][str(glob_ini)]['sampling_sparseness'] = sampling_sparseness - res_dic[param][method_cur][str(glob_ini)]['param_arr'] = param_arr + res_dic[method_cur][param][str(glob_ini)] = {} + res_dic[method_cur][param][str(glob_ini)]['sampling_sparseness'] = sampling_sparseness + res_dic[method_cur][param][str(glob_ini)]['param_arr'] = param_arr # With intercept at axis. # Calculate sample correlation coefficient, measure of goodness-of-fit of linear regression @@ -2017,18 +2020,143 @@ print(param, method_ref, method_cur, sampling_sparseness, glob_ini, r_xy**2, a) # Store to result dic. - res_dic[param][method_cur][str(glob_ini)]['r_xy'] = r_xy - res_dic[param][method_cur]['r_xy'].append(r_xy) - res_dic[param][method_cur][str(glob_ini)]['a'] = a - res_dic[param][method_cur]['a'].append(a) - - res_dic[param][method_cur]['sampling_sparseness'] = asarray(res_dic[param][method_cur]['sampling_sparseness']) - res_dic[param][method_cur]['glob_ini'] = asarray(res_dic[param][method_cur]['glob_ini']) - - res_dic[param][method_cur]['r_xy'] = asarray(res_dic[param][method_cur]['r_xy']) - res_dic[param][method_cur]['a'] = asarray(res_dic[param][method_cur]['a']) + res_dic[method_cur][param][str(glob_ini)]['r_xy'] = r_xy + res_dic[method_cur][param]['r_xy'].append(r_xy) + res_dic[method_cur][param][str(glob_ini)]['a'] = a + res_dic[method_cur][param]['a'].append(a) + + res_dic[method_cur][param]['sampling_sparseness'] = asarray(res_dic[method_cur][param]['sampling_sparseness']) + res_dic[method_cur][param]['glob_ini'] = asarray(res_dic[method_cur][param]['glob_ini']) + + res_dic[method_cur][param]['r_xy'] = asarray(res_dic[method_cur][param]['r_xy']) + res_dic[method_cur][param]['a'] = asarray(res_dic[method_cur][param]['a']) return res_dic + + + def plot_min_stat(self, min_stat_dic=None, methods=[], list_glob_ini=[], show=False, write_stats=False): + + # Catch min and max values for all methods. + min_a = 1.0 + max_a = 0.0 + + min_r_xy2 = 1.0 + max_r_xy2 = 0.0 + + # Prepare header for writing. + selection = min_stat_dic['selection'] + params_list = min_stat_dic['params_list'] + analysis = min_stat_dic['analysis'] + + # For writing out stats. + headings = [] + data_dic = OrderedDict() + i_max = 0 + + for method in methods: + if method not in min_stat_dic: + continue + + # Define figure + fig, axises = plt.subplots(nrows=len(params_list), ncols=1) + fig.suptitle('Stats per NI %s' % method) + + # Loop over params + data_dic[method] = OrderedDict() + + for j, param in enumerate(params_list): + data_dic[method][param] = OrderedDict() + + # Use NI as x. + NI = min_stat_dic[method][param]['glob_ini'] + + # Use sampling_sparseness as x. + SS = min_stat_dic[method][param]['sampling_sparseness'] + + # Add to headings. + headings = headings + ['method_%s'%param, 'SS', 'NI', 'slope', 'rxy2'] + + # Get stats. + # Linear regression slope, without intercept + a = min_stat_dic[method][param]['a'] + + if max(a) > max_a: + max_a = max(a) + if min(a) < min_a: + min_a = min(a) + + # sample correlation coefficient, without intercept + r_xy = min_stat_dic[method][param]['r_xy'] + r_xy2 = r_xy**2 + + if max(r_xy2) > max_r_xy2: + max_r_xy2 = max(r_xy2) + if min(r_xy2) < min_r_xy2: + min_r_xy2 = min(r_xy2) + + # Add to data. + for i, NI_i in enumerate(NI): + SS_i = SS[i] + a_i = a[i] + r_xy2_i = r_xy2[i] + data_dic[method][param][str(i)] = ["%3.5f"%SS_i, "%i"%NI_i, "%3.5f"%a_i, "%3.5f"%r_xy2_i] + if i > i_max: + i_max = i + + ax = axises[j] + ax.plot(SS, a, ".-", label='%s_%s_a' % (method, param) ) + ax.plot(SS, r_xy2, "o--", label='%s_%s_r_xy2' % (method, param) ) + ax.legend(loc='lower left', shadow=True, prop = fontP) + ax.set_xlabel('SS') + ax.invert_xaxis() + #ax.set_ylim(min_a*0.95, max_a*1.05) + + + # Loop over methods for writing data. + data = [] + + # Loop over all lines. + for i in range(0, i_max+1): + data_i = [] + for method, data_dic_m in data_dic.iteritems(): + # Loop over all params + for j, param in enumerate(params_list): + # Loop over all possible data points. + if str(i) in data_dic_m[param]: + data_i = data_i + ["%s_%s" % (method, param)] + data_dic_m[param][str(i)] + else: + data_i = data_i + ["%s_%s" % (method, param)] + ["0", "0", "0", "0"] + + data.append(data_i) + + # Determine filename. + if selection == None: + file_name_ini = '%s_stat_all' % analysis + else: + file_name_ini = '%s_stat_sel' % analysis + + # Write png. + png_file_name = file_name_ini + '.png' + png_file_path = get_file_path(file_name=png_file_name, dir=self.results_dir) + + # Write to file. + if write_stats: + # save figure + plt.savefig(png_file_path, bbox_inches='tight') + + file_name = file_name_ini + '.txt' + path = self.results_dir + file_obj, file_path = open_write_file(file_name=file_name, dir=path, force=True, compress_type=0, verbosity=1, return_path=True) + + # Write data. + write_data(out=file_obj, headings=headings, data=data) + + # Close file. + file_obj.close() + + # Plot data. + if show: + plt.show() def interpreter_start(self): 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=26066&r1=26065&r2=26066&view=diff ============================================================================== --- trunk/test_suite/system_tests/relax_disp.py (original) +++ trunk/test_suite/system_tests/relax_disp.py Fri Sep 26 17:29:49 2014 @@ -6015,45 +6015,6 @@ RDR.plot_int_corr(corr_data=corr_data, show=True) - # Try plot some R2eff correlations. - if False: - # Now calculate R2eff. - #RDR.calc_r2eff(methods=methods, list_glob_ini=[128, 126]) - - # Try for bad data. - #RDR.calc_r2eff(methods=['FT'], list_glob_ini=[6, 4]) - - # Collect r2eff values. - r2eff_ft_all = RDR.col_r2eff(method='FT', list_glob_ini=[128, 126], selection=None) - - # Now make a spin selection. - selection = ':2,3' - r2eff_ft_sel = RDR.col_r2eff(method='FT', list_glob_ini=[128, 126], selection=selection) - # Print the length of datasets, depending on selection. - print( "All spins", len(r2eff_ft_all['128']['r2eff_arr']), "Selection spins", len(r2eff_ft_sel['128']['r2eff_arr']) ) - - # For all spins, mdd - r2eff_mdd_all = RDR.col_r2eff(method='MDD', list_glob_ini=[128, 126], selection=None) - r2eff_mdd_sel = RDR.col_r2eff(method='MDD', list_glob_ini=[128, 126], selection=selection) - - # Plot correlation of intensity - fig1 = [[r2eff_ft_all, r2eff_mdd_all], ['FT', 'MDD'], [128, 128]] - fig2 = [[r2eff_ft_sel, r2eff_mdd_sel], ['FT sel', 'MDD sel'], [128, 128]] - corr_data = [fig1, fig2] - - fig3 = [[r2eff_ft_all, r2eff_ft_all], ['FT', 'FT'], [128, 126]] - fig4 = [[r2eff_ft_all, r2eff_mdd_all], ['FT', 'MDD'], [128, 126]] - corr_data = [fig3, fig4] - - #fig5 = [[r2eff_ft_all, r2eff_mdd_all], ['FT', 'MDD'], [128, 128]] - fig5 = [[r2eff_ft_all, r2eff_mdd_all], ['FT', 'MDD'], [128, 126]] - fig6 = [[r2eff_ft_all, r2eff_ft_all], ['FT', 'FT'], [128, 126]] - #fig6 = [[r2eff_mdd_all, r2eff_mdd_all], ['MDD', 'MDD'], [128, 126]] - corr_data = [fig5, fig6] - - RDR.plot_r2eff_corr(corr_data=corr_data, show=True) - - # Try write some R2eff correlations. if True: selection = None @@ -6177,33 +6138,64 @@ # Plot statistics. # Try plot some minimisation correlations. if True: - selection = None + selections = [None, ':2,3'] + for selection in selections: + # Collect param values. + analysis = 'min_ind' + min_ft_sel = RDR.col_min(method='FT', model=MODEL_CR72, analysis=analysis, list_glob_ini=[128], selection=selection) + min_mdd_sel = RDR.col_min(method='MDD', model=MODEL_CR72, analysis=analysis, list_glob_ini=range(126, 130, 2)[::-1], selection=selection) + + fig1 = [[min_ft_sel, min_mdd_sel], ['FT', 'MDD'], [128, 128]] + fig2 = [[min_ft_sel, min_mdd_sel], ['FT', 'MDD'], [128, 126]] + corr_data = [fig1, fig2] + + write_stats = True + RDR.plot_min_corr(corr_data=corr_data, show=False, write_stats=write_stats) + + # Open stat file. + if write_stats: + for i, corr_data_i in enumerate(corr_data): + data, methods, glob_inis = corr_data[i] + data_x, data_y = data + method_x, method_y = methods + glob_ini_x, glob_ini_y = glob_inis + + file_name_ini = '%s_%s_%s_%s_%s' % (analysis, method_x, glob_ini_x, method_y, glob_ini_y) + + if selection == None: + file_name = file_name_ini + '_all.txt' + else: + file_name = file_name_ini + '_sel.txt' + path = RDR.results_dir + data = extract_data(file=file_name, dir=path) + + # Loop over the lines. + for i, data_i in enumerate(data): + print(i, data_i) + + # Try plot some minimisation statistics. + if True: # Collect param values. - analysis = 'min_ind' - min_ft_sel = RDR.col_min(method='FT', model=MODEL_CR72, analysis=analysis, list_glob_ini=[128], selection=selection) - min_mdd_sel = RDR.col_min(method='MDD', model=MODEL_CR72, analysis=analysis, list_glob_ini=range(126, 130, 2)[::-1], selection=selection) - - fig1 = [[min_ft_sel, min_mdd_sel], ['FT', 'MDD'], [128, 128]] - fig2 = [[min_ft_sel, min_mdd_sel], ['FT', 'MDD'], [128, 126]] - corr_data = [fig1, fig2] - - write_stats = True - RDR.plot_min_corr(corr_data=corr_data, show=False, write_stats=write_stats) - - # Open stat file. - if write_stats: - for i, corr_data_i in enumerate(corr_data): - data, methods, glob_inis = corr_data[i] - data_x, data_y = data - method_x, method_y = methods - glob_ini_x, glob_ini_y = glob_inis - - file_name_ini = '%s_%s_%s_%s_%s' % (analysis, method_x, glob_ini_x, method_y, glob_ini_y) - + #selections = [None, ':2,3'] + selections = [None, ':2,3'] + for selection in selections: + analysis = 'min_ind' + min_ft_sel = RDR.col_min(method='FT', model=MODEL_CR72, analysis=analysis, list_glob_ini=[128], selection=selection) + min_mdd_sel = RDR.col_min(method='MDD', model=MODEL_CR72, analysis=analysis, list_glob_ini=range(126, 130, 2)[::-1], selection=selection) + + # Get param stats. + min_stat_dic = RDR.get_min_stat_dic(list_r2eff_dics=[min_ft_sel, min_mdd_sel], list_glob_ini=[128, 126]) + + ## Plot R2eff stats + write_stats = True + RDR.plot_min_stat(min_stat_dic=min_stat_dic, methods=['FT', 'MDD'], list_glob_ini=[128, 126, 6], show=False, write_stats=write_stats) + + # Open stat file. + if write_stats: if selection == None: - file_name = file_name_ini + '_all.txt' + file_name = '%s_stat_all.txt' % (analysis) else: - file_name = file_name_ini + '_sel.txt' + file_name = '%s_stat_sel.txt' % (analysis) path = RDR.results_dir data = extract_data(file=file_name, dir=path) @@ -6312,6 +6304,36 @@ # Loop over the lines. for i, data_i in enumerate(data): print(i, data_i) + + # Try plot some minimisation statistics. + if True: + # Collect param values. + selections = [':2,3'] + for selection in selections: + analysis = 'min' + min_ft_sel = RDR.col_min(method='FT', model=MODEL_CR72, analysis=analysis, list_glob_ini=[128], selection=selection) + min_mdd_sel = RDR.col_min(method='MDD', model=MODEL_CR72, analysis=analysis, list_glob_ini=range(126, 130, 2)[::-1], selection=selection) + + # Get param stats. + min_stat_dic = RDR.get_min_stat_dic(list_r2eff_dics=[min_ft_sel, min_mdd_sel], list_glob_ini=[128, 126]) + + ## Plot R2eff stats + write_stats = True + RDR.plot_min_stat(min_stat_dic=min_stat_dic, methods=['FT', 'MDD'], list_glob_ini=[128, 126, 6], show=False, write_stats=write_stats) + + # Open stat file. + if write_stats: + if selection == None: + file_name = '%s_stat_all.txt' % (analysis) + else: + file_name = '%s_stat_sel.txt' % (analysis) + path = RDR.results_dir + data = extract_data(file=file_name, dir=path) + + # Loop over the lines. + for i, data_i in enumerate(data): + print(i, data_i) + def test_r1rho_kjaergaard_auto(self):