Skip to content
Snippets Groups Projects
Commit 1fd453f9 authored by Jocabed Martínez's avatar Jocabed Martínez
Browse files

Delete todo_r.py

parent 39e46d9f
No related branches found
No related tags found
No related merge requests found
# coding=utf-8
import warnings
def fxn():
warnings.warn("deprecated", DeprecationWarning)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
fxn()
import sys
input1 = sys.argv[1]
print(input1)
# insert at 1, 0 is the script path (or '' in REPL)
sys.path.insert(1, '../Functions')
from imports import *
from constants import *
from functions import *
path = input1
AA = []
#for p in path:
#print(p)
col = read_col_names(path)
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))
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)))
best_parmeters = return_values(AA)
grafica_data_ybestfit(path, dfr, best_parmeters, ss="r")
y_pred_NFW_bestfit = vnfw(best_parmeters["M200_NFW(Msun)"][1], best_parmeters["c_NFW"][1], dfr["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) , dfr["r(kpc)"]*1000/abs(best_parmeters["Rc_ISO(pc)"][1]), G = G)
chi2_NFW_bestfit = (((y_pred_NFW_bestfit-dfr["v(km/s)"])/(dfr["err_v(km/s)"]))**2).sum()
chi2_ISO_bestfit = (((y_pred_ISO_bestfit-dfr["v(km/s)"])/(dfr["err_v(km/s)"]))**2).sum()
print("Mejor $\chi^2$ para \nNFW: {} \nISO: {}".format(chi2_NFW_bestfit,chi2_ISO_bestfit))
chi2_NFW_bestfit_red = chi2_NFW_bestfit / (len(dfr)-2)
chi2_ISO_bestfit_red = chi2_ISO_bestfit / (len(dfr)-2)
print("Mejor $\chi^2$ reducido para \nNFW: {} \nISO: {}".format(chi2_NFW_bestfit_red,chi2_ISO_bestfit_red))
print("La probabilidad cumulativa de encontrar un $\chi^2$ menor para NFW es: {}".format(scipy.stats.chi2.cdf(chi2_NFW_bestfit,len(dfr)-2)))
print("La probabilidad cumulativa de encontrar un $\chi^2$ menor para ISO es: {}".format(scipy.stats.chi2.cdf(chi2_ISO_bestfit,len(dfr)-2)))
combs_NFW_ra = parameters_range(best_parmeters, i = 0, ss = "r")
combs_ISO_ra = parameters_range(best_parmeters, i = 1, ss = "r")
chis_NFW_ra = chi2_parameters_range(dfr, combs_NFW_ra, i = 0)
chis_ISO_ra = chi2_parameters_range(dfr, combs_ISO_ra, i = 1)
ch2_inter_NFW_ = chi2_intervals(best_parmeters, dfr, chis_NFW_ra, i = 0, ss = "r")
ch2_inter_ISO_ = chi2_intervals(best_parmeters, dfr, chis_ISO_ra, i = 1, ss = "r")
para_range_NFW = []
para_range_ISO = []
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))
# -- índice dentro del arreglo de chi^2s que es el nuevo mínimo para NFW e ISO
# newi_min_ch2_NFW = [i for i in ch2_inter_NFW_[0] if i[1] == min(list(zip(*ch2_inter_NFW_[0]))[1])][0][0]
# newi_min_ch2_ISO = [i for i in ch2_inter_ISO_[0] if i[1] == min(list(zip(*ch2_inter_ISO_[0]))[1])][0][0]
# new_para_min_ch2_NFW = combs_NFW_ra[newi_min_ch2_NFW]
# new_para_min_ch2_ISO = combs_ISO_ra[newi_min_ch2_ISO]
# new_para_min_ch2_b = np.concatenate((new_para_min_ch2_NFW,new_para_min_ch2_ISO))
# print("Los nuevos parámetros para NFW son: {} con un $\chi^2 = {}$"
# .format(new_para_min_ch2_NFW, chis_NFW_ra[newi_min_ch2_NFW]))
# print("Los nuevos parámetros para ISO son: {} con un $\chi^2 = {}$"
# .format(new_para_min_ch2_ISO, chis_ISO_ra[newi_min_ch2_ISO]))
# print("La probabilidad cumulativa de encontrar un $\chi^2$ menor al nuevo para NFW es: {}".format(scipy.stats.chi2.cdf(chis_NFW_ra[newi_min_ch2_NFW]*(len(df1)-2),len(df1)-2)))
# print("La probabilidad cumulativa de encontrar un $\chi^2$ menor al nuevo para ISO es: {}".format(scipy.stats.chi2.cdf(chis_ISO_ra[newi_min_ch2_ISO]*(len(df1)-2),len(df1)-2)))
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 )
#grafica_data_ybestfit(path, df1, best_parmeters, ss="ra")
file1 = open("results.txt","a")
file1.write("\n{}&{}&{:.2f}&{:.2f}&{:.2f}&{:.2f}&{:.2f}&{:.2f}&".format(
path.split("/")[-1][:-4][5:]
,len(dfr)
,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()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment