<jupyter_start><jupyter_text>Read and initial plot of data<jupyter_code>drct = "../scripts/get_globular_clusters/output"
gc = pd.read_csv(os.path.join(drct, "NGC_104_47Tuc.csv"), skiprows=54, sep=" ")
gc.head()
gc["r"] = np.hypot(gc.x, gc.y)
gc["pm"] = np.hypot(gc.pmx, gc.pmy)
gc.describe()<jupyter_output><empty_output><jupyter_text>Spatial Distribution<jupyter_code>plt.figure(figsize=(7,7))
H, xb, yb, _ = plt.hist2d(gc.x, gc.y, bins=200, range=[[-0.66,0.66],[-0.66,0.66]], norm=LogNorm(), cmap="gnuplot2")
sns.distplot(gc.r)
plt.xlabel("r")<jupyter_output><empty_output><jupyter_text>Pruning by scale radius and upper limit of proper motion<jupyter_code>gcs = gc[(gc.r<18.3/60) & (gc.pm<15)]
sns.distplot(gcs.r)
plt.xlabel("r")<jupyter_output><empty_output><jupyter_text>Distribution of proper motion<jupyter_code>plt.figure(figsize=(6,6))
H, xb, yb, _ = plt.hist2d(gcs.pmx, gcs.pmy, bins=100, range=[[-10,15],[-15,10]], norm=LogNorm(), cmap="gnuplot2")
sns.distplot(gcs.pm[gcs.pm<15])
plt.xlabel("PM")
H, xb, yb, _ = plt.hist2d(gc.r, gc.pm, bins=50, range=[[0,0.6],[0,10]], cmap="gnuplot2")
# R scale
plt.axvline(18.3/60, color="gold")
plt.xlabel("r")
plt.ylabel("PM")<jupyter_output><empty_output><jupyter_text>Binning Start From Here Bin the distribution of pm in radial bins.<jupyter_code>bins = np.logspace(np.log10(0.1), np.log10(0.3), 9)

sns.set_palette("RdBu", len(bins))
r_rbin, z_rbin = profile_binning(gc.r, gc.pm, bins=bins, z_clip=[2,15],
                                 return_bin=False, plot=True)<jupyter_output><empty_output><jupyter_text>Acquire all the stars from a particular bin:<jupyter_code>z_bins = profile_binning(gcs.r, gcs.pm, bins=bins, z_clip=[2,15],
                         return_bin=True, plot=False)
z_bins["5"]<jupyter_output><empty_output><jupyter_text>PM radial profile<jupyter_code>plt.figure(figsize=(7,5))
plt.plot(r_rbin, z_rbin, "k-o")
plt.xlabel("r")
plt.ylabel("PM")<jupyter_output><empty_output><jupyter_text>GMM start from here Perform GMM decomposition. The best number of gaussians is given by the lowest BIC. One can also specify the number of components by giving values to parameter n_comp in GMM( ).<jupyter_code>gmm_bins = []
for i in range(len(bins)-1):
    gcb = gcs[(bins[i]<gcs.r)&(gcs.r<bins[i+1])]
    gmm = GMM(gcb.pm, max_n_comp=6, verbose=False)
    gmm_bins.append(gmm)<jupyter_output><empty_output><jupyter_text>Probability that a star is giving by the major population. Passing a subsample<jupyter_code>prob_main_pop = predict_main_pop(z_bins["5"]["pm"], gmm_bins[5])
prob_main_pop<jupyter_output><empty_output><jupyter_text>Visualize GMM discrimination in each bin:<jupyter_code>plot_predict_main_pop(gcs, bins, gmm_bins, p_thre=0.8)<jupyter_output><empty_output><jupyter_text>Visualize GMM results combining bins:<jupyter_code>for i, b in enumerate(bins[:-1]):
    
    r = z_bins[str(i)]["r"]
    pm = z_bins[str(i)]["pm"]
    gmm = gmm_bins[i]
    
    prob_main_pop = predict_main_pop(pm, gmm)
    is_main_pop = prob_main_pop > 0.8
    
    plt.scatter(r[is_main_pop], pm[is_main_pop], color="navy", alpha=0.1, s=2)
    plt.scatter(r[~is_main_pop], pm[~is_main_pop], color="gold", alpha=0.3, s=5)
    plt.axvline(bins[i], color="k", alpha=0.2)
    plt.ylabel("Proper Motion (mas)")
    plt.xlabel("Radius (deg)")
    plt.ylim(0, 15)
plt.show()<jupyter_output><empty_output><jupyter_text>Surface Density DON'T NEED TO DO TODO a much better estimate of the surface density than just a histrogram<jupyter_code>suface_density = []

for i, b in enumerate(bins[:-1]):
    
    r = z_bins[str(i)]["r"]
    pm = z_bins[str(i)]["pm"]
    gmm = gmm_bins[i]
    
    prob_main_pop = predict_main_pop(pm, gmm)
    is_main_pop = prob_main_pop > 0.8
    
    area = np.pi * (bins[i+1]**2 - b**2)
    suface_density.append(len(r) / area)

plt.plot((bins[:-1]+bins[1:])/2, suface_density)
import astropy
from astropy.modeling.functional_models import KingProjectedAnalytic1D
from astropy.modeling import models, fitting
x = (bins[:-1]+bins[1:])/2
y = suface_density
KP_int = KingProjectedAnalytic1D(amplitude = 1, r_core = 1., r_tide = 1)

fit_KP = fitting.LevMarLSQFitter()
KP = fit_KP(KP_int, x, y, maxiter=1000)

plt.figure(figsize=(8,5))
plt.loglog(x, y, 'ko')
plt.plot(x, KP(x), label='Gaussian')
KP.amplitude, KP.r_core, KP.r_tide<jupyter_output><empty_output><jupyter_text>Velocity Dispersion Assume Plummer sphere model.Still need to normalize.Do better fit:    - https://symfit.readthedocs.io/en/stable/intro.html    - https://lmfit.github.io/lmfit-py/fitting.html    - and an MCMC (can use symfit)<jupyter_code>from galpy.df import jeans
from galpy import potential

x = (bins[:-1]+bins[1:])/2 / summary_table
%%time

for amp, scale in zip(np.linspace(10, 100, num=20), np.linspace(0.1, 5, num=20)):
    pot = potential.PlummerPotential(amp=amp, b=scale) 
    model = np.array([jeans.sigmar(pot, r) for r in x])
# Grid Parameters for 47 Tuc

gc_amp = np.logspace(9e5, 1.1e6, num=5e1) * u.solMass  # The GC mass varies pretty wildly: 10^4 - few*10^6
gc_scale = np.linspace(0.5, 1.5, num=1e2)
bh_amp = np.linspace(0, 1e5, num=1e2) * u.solMass

# radius
radii = np.linspace(0, 10, num=1e3)  # units of radius / scale_radius

pot = potential.PlummerPotential(amp=gc_amp, b=gc_scale * )
potential.PlummerPotential??
rnorm = ((x * u.deg).to_value(u.rad) / (18.3 *u.arcmin).to_value(u.rad))

from astropy import units as u
amp = 10**6 * u.solMass
scale = 1

pot = potential.PlummerPotential(amp=amp, b=scale)
model1 = np.array([jeans.sigmar(pot, r).to_value(u.km/u.s) for r in rnorm]) * u.km / u.s
plt.plot(x, model1, label='NO')

bh_amp=10**3 * u.solMass
amp = 10**6 * u.solMass - bh_amp
pot = potential.PlummerPotential(amp=amp, b=scale) + potential.KeplerPotential(amp=bh_amp)
model2 = np.array([jeans.sigmar(pot, r).to_value(u.km/u.s) for r in rnorm]) * u.km / u.s
plt.plot(x, model2, label='1e3')

bh_amp=10**4 * u.solMass
amp = 10**6 * u.solMass - bh_amp
pot = potential.PlummerPotential(amp=amp, b=scale) + potential.KeplerPotential(amp=bh_amp)
model3 = np.array([jeans.sigmar(pot, r).to_value(u.km/u.s) for r in rnorm]) * u.km / u.s
plt.plot(x, model3, label=1e4)
plt.legend()
import astropy.units as u
from lmfit import Parameters, fit_report, minimize
import galpy
from galpy.df import jeans
from galpy import potential

fit_params = Parameters()
fit_params.add('amp', value=100000, min=0)
fit_params.add('scale', value=1, min=0)


def residual(pars, x, data=None):
    """Model and subtract data."""
    vals = pars.valuesdict()
    amp = vals['amp']
    scale = vals['scale']

    pot = potential.PlummerPotential(amp=amp, b=scale) 
    model = np.array([jeans.sigmar(pot, r) for r in x])

    if data is None:
        return model
    return model - data


# getting data
vel_disp = []
for i, b in enumerate(bins[:-1]):
    
    r = z_bins[str(i)]["r"]
    pm = z_bins[str(i)]["pm"]
    gmm = gmm_bins[i]
    
    prob_main_pop = predict_main_pop(pm, gmm)
    is_main_pop = prob_main_pop > 0.8

    vs = ((4 * u.kpc) * pm[is_main_pop] * u.mas/u.yr).to_value(u.km*u.rad/u.s)
    
    disp = np.var(vs)
    
    vel_disp.append(disp)

x = (bins[:-1]+bins[1:])/2
data = vel_disp

plt.plot(x, vel_disp)
mi = minimize(residual, fit_params, args=(x,), kws={'data': data})

print(fit_report(mi))
plt.plot(x, data, 'b')
plt.plot(x, data + residual(mi.params, x, data=data), 'r', label='best fit')
plt.legend(loc='best')
plt.show()<jupyter_output><empty_output><jupyter_text>With a BH<jupyter_code>import astropy.units as u
from lmfit import Parameters, fit_report, minimize
import galpy
from galpy.df import jeans
from galpy import potential

fit_params = Parameters()
fit_params.add('amp', value=100000, min=0)
fit_params.add('scale', value=1, min=0)
fit_params.add('bh_amp', value=10000, min=0)


def residual(pars, x, data=None):
    """Model and subtract data."""
    vals = pars.valuesdict()
    amp = vals['amp']
    scale = vals['scale']
    bh_amp = vals['bh_amp']

    pot = potential.PlummerPotential(amp=amp, b=scale) + potential.KeplerPotential(amp=bh_amp)
    model = np.array([jeans.sigmar(pot, r) for r in x])

    if data is None:
        return model
    return model - data


# getting data
vel_disp = []
for i, b in enumerate(bins[:-1]):
    
    r = z_bins[str(i)]["r"]
    pm = z_bins[str(i)]["pm"]
    gmm = gmm_bins[i]
    
    prob_main_pop = predict_main_pop(pm, gmm)
    is_main_pop = prob_main_pop > 0.8

    vs = ((4 * u.kpc) * pm[is_main_pop] * u.mas/u.yr).to_value(u.km*u.rad/u.s)
    
    disp = np.var(vs)
    
    vel_disp.append(disp)

x = (bins[:-1]+bins[1:])/2
data = vel_disp

mi = minimize(residual, fit_params, args=(x,), kws={'data': data})

plt.plot(x, data, 'b')
plt.plot(x, data + residual(mi.params, x, data=data), 'r', label='best fit')
plt.legend(loc='best')
plt.show()

print(fit_report(mi))
galpy.potential.Kin 
# raise Exception('Run MCMC?')

# import lmfit as lf

# mi.params.add('__lnsigma', value=np.log(0.1), min=np.log(0.001), max=np.log(2))

# res = lf.minimize(residual, method='emcee', args=(x,), kws={'data': data}, nan_policy='omit', burn=300, steps=1000, thin=20,
#                   params=mi.params, is_weighted=False, progress=True)

# plt.plot(res.acceptance_fraction)
# plt.xlabel('walker')
# plt.ylabel('acceptance fraction')
# plt.show()

# import corner

# emcee_plot = corner.corner(res.flatchain, labels=res.var_names,
#                            truths=list(res.params.valuesdict().values()))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1684455246750502, 'fraction_numerical': 0.027052731691852478, 'mean_word_length': 4.650149253731343, 'pattern_counts': {'":': 0, '<': 62, '<?xml version=': 0, '>': 61, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '29062', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4157, 'n_tokens_neox': 3905, 'n_words': 745, 'path': 'notebooks/OLD/ReadGC+Binning+GMM.ipynb', 'repo_name': 'nstarman/jas1101_project', 'size': '9463'}