Upload New File r

code/todo_r.py

+# 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 'r'
+df1 = dfr
+# -- grafica los datos experimentales con las curvas NFW e ISO, y el residuo
+grafica_data_ybestfit(path, df1, best_parmeters, ss="r")
+
+# -- calcular el valor de la velociada para esos parametros en NFW e ISO
+y_pred_NFW_bestfit = vnfw(best_parmeters["M200_NFW(Msun)"][1], best_parmeters["c_NFW"][1], df1["r(kpc)"]*1000/r200(best_parmeters["M200_NFW(Msun)"][1]), G = G, H0 = H0)
+y_pred_ISO_bestfit = viso((best_parmeters["rho0_ISO(Msun/pc3)"][1]*best_parmeters["Rc_ISO(pc)"][1]**2) , df1["r(kpc)"]*1000/abs(best_parmeters["Rc_ISO(pc)"][1]), 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 = "r")
+combs_ISO_ra = parameters_range(best_parmeters, i = 1, ss = "r")
+
+# -- 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 = "r")
+ch2_inter_ISO_ = chi2_intervals(best_parmeters, df1, chis_ISO_ra, i = 1, ss = "r")
+
+# -- 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 = "r", i = 0 )
+plot_ellipse_confidence_interval(path,para_range_ISO, best_parmeters, ss = "r", 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)"][1]
+,best_parmeters["c_NFW"][1]
+,chi2_NFW_bestfit_red
+,abs(best_parmeters["rho0_ISO(Msun/pc3)"][1])
+,abs(best_parmeters["Rc_ISO(pc)"][1])
+,chi2_ISO_bestfit_red))
+file1.close()