diff --git a/code/todo_a.py b/code/todo_a.py new file mode 100644 index 0000000000000000000000000000000000000000..10126f119e409c4a3f3f39ba1d428e325fa59dd0 --- /dev/null +++ b/code/todo_a.py @@ -0,0 +1,111 @@ +# coding=utf-8 +# -- elimina algunos de los mensajes de advertencia +import warnings +def fxn(): + warnings.warn("deprecated", DeprecationWarning) + +with warnings.catch_warnings(): + warnings.simplefilter("ignore") + fxn() + +# -- recibe como entrada el path de los datos i.e python3 todo.py "../galaxies_data/rcugc11012.txt" +import sys +input1 = sys.argv[1] +print(input1) + +# -- especifica donde se encuentran los archivos de imports, constants y functions +sys.path.insert(1, '../Functions') + +from imports import * +from constants import * +from functions import * + + +path = input1 +# -- crear lista para guardar los valores +AA = [] +# -- nombre de las columnas +col = read_col_names(path) +# -- lee los datos y pone col como nombre de columnas +df1 = pd.read_csv(path,skiprows=1,header=None,delimiter=r"\s+",names=col) +# -- quitar errores nulos para evitar problemas con leastsq +df1 = df1[df1["err_v(km/s)"]!=0] +print(len(df1)) +# -- ejecuta scipy.optimize.leastsq para cada datos 'ra', 'r' y 'a' para la galaxia en cuestion +for labe in ['ra','r','a']: + if labe == 'ra': + gal = df1.copy() + if labe == 'r': + dfr = df1.copy()[df1.copy()['Side'] == 'r'] + gal = dfr + if labe == 'a': + dfa = df1.copy()[df1.copy()['Side'] == 'a'] + gal = dfa + + AA.append((path.split("/")[-1][:-4],determine(gal,labe))) +# -- dataframe de valores y error formal de los dos parametros para NFW e ISO y para cada dato 'ra', 'r' y 'a' +best_parmeters = return_values(AA) + + +# -- solo trabajar con los datos 'a' +df1 = dfa +# -- grafica los datos experimentales con las curvas NFW e ISO, y el residuo +grafica_data_ybestfit(path, df1, best_parmeters, ss="a") + +# -- calcular el valor de la velociada para esos parametros en NFW e ISO +y_pred_NFW_bestfit = vnfw(best_parmeters["M200_NFW(Msun)"][2], best_parmeters["c_NFW"][2], df1["r(kpc)"]*1000/r200(best_parmeters["M200_NFW(Msun)"][2]), G = G, H0 = H0) +y_pred_ISO_bestfit = viso((best_parmeters["rho0_ISO(Msun/pc3)"][2]*best_parmeters["Rc_ISO(pc)"][2]**2) , df1["r(kpc)"]*1000/abs(best_parmeters["Rc_ISO(pc)"][2]), G = G) + +# -- calcula chi2 para NFW e ISO +chi2_NFW_bestfit = (((y_pred_NFW_bestfit-df1["v(km/s)"])/(df1["err_v(km/s)"]))**2).sum() +chi2_ISO_bestfit = (((y_pred_ISO_bestfit-df1["v(km/s)"])/(df1["err_v(km/s)"]))**2).sum() +print("Mejor $\chi^2$ para \nNFW: {} \nISO: {}".format(chi2_NFW_bestfit,chi2_ISO_bestfit)) + +# -- calcula chi2 reducido +chi2_NFW_bestfit_red = chi2_NFW_bestfit / (len(df1)-2) +chi2_ISO_bestfit_red = chi2_ISO_bestfit / (len(df1)-2) +print("Mejor $\chi^2$ reducido para \nNFW: {} \nISO: {}".format(chi2_NFW_bestfit_red,chi2_ISO_bestfit_red)) + +# -- calcular probabilidad para referencia +print("La probabilidad cumulativa de encontrar un $\chi^2$ menor para NFW es: {}".format(scipy.stats.chi2.cdf(chi2_NFW_bestfit,len(df1)-2))) +print("La probabilidad cumulativa de encontrar un $\chi^2$ menor para ISO es: {}".format(scipy.stats.chi2.cdf(chi2_ISO_bestfit,len(df1)-2))) + +# -- crea la matriz 2x2 para los parametros +combs_NFW_ra = parameters_range(best_parmeters, i = 0, ss = "a") +combs_ISO_ra = parameters_range(best_parmeters, i = 1, ss = "a") + +# -- calcular chi2 para cada par de parametros de la matriz +chis_NFW_ra = chi2_parameters_range(df1, combs_NFW_ra, i = 0) +chis_ISO_ra = chi2_parameters_range(df1, combs_ISO_ra, i = 1) + +# -- divide el chi2 de acuerdo a intervalos, para contruir intervalos de confiaza, en la función como tal se definen estos intervalos +ch2_inter_NFW_ = chi2_intervals(best_parmeters, df1, chis_NFW_ra, i = 0, ss = "a") +ch2_inter_ISO_ = chi2_intervals(best_parmeters, df1, chis_ISO_ra, i = 1, ss = "a") + +# -- lista para guardar los parametros de interes, para cada intervalo de chi2 +para_range_NFW = [] +para_range_ISO = [] + +# -- extrae el par de parametros, dependiendo del intervalo de chi2 en el que esten +for kkk in range(1,5): + para_range_NFW.append(parameters_each_interval(combs_NFW_ra, ch2_inter_NFW_, kk = kkk)) + para_range_ISO.append(parameters_each_interval(combs_ISO_ra, ch2_inter_ISO_, kk = kkk)) + + +# -- grafica elipses de intervalos de confianza de chi2 para NFW e ISO +plot_ellipse_confidence_interval(path,para_range_NFW, best_parmeters, ss = "a", i = 0 ) +plot_ellipse_confidence_interval(path,para_range_ISO, best_parmeters, ss = "a", i = 1 ) + + +# -- guarda los datos que necesitamos para el reporte +file1 = open("results.txt","a") +file1.write("\n{}&{}&{:.2e}&{:.2f}&{:.2f}&{:.2f}&{:.2f}&{:.2f}&".format( +path.split("/")[-1][:-4][5:] +,len(df1) +,best_parmeters["M200_NFW(Msun)"][2] +,best_parmeters["c_NFW"][2] +,chi2_NFW_bestfit_red +,abs(best_parmeters["rho0_ISO(Msun/pc3)"][2]) +,abs(best_parmeters["Rc_ISO(pc)"][2]) +,chi2_ISO_bestfit_red)) +file1.close()