Merge pull request #110 from ChunHuangPhy/ChunHuangPhy-patch-1

Chun huang phy patch 1

postprocessing/MR_posterior_distribution_and_CornerPlot.ipynb

@@ -0,0 +1,350 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "744ffaf7",
+   "metadata": {},
+   "source": [
+    "# Plot the posterior mass-radius countours"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b2683b98",
+   "metadata": {},
+   "source": [
+    "Here is an example of how to obtain the posteriro mass-radius countours from the posterior distribution of the EOS parameters.\n",
+    "\n",
+    "After obtaining the distributions of the EOS parameters, we can map the posterior distributions of the specific EOS in mass-radius (M-R) space. This procedure helps us to understand how observational constraints affect the EOS. Each point in the EOS parameter space is uniquely correlated with a point in the EOS posterior parameter space. Then, by varying the central density, EOS points can be mapped onto the M-R plane by deriving the Tolman-Oppenheimer-Volkoff (TOV) equations.\n",
+    "\n",
+    "Here is a step-by-step explanation.\n",
+    "\n",
+    "First import all the package that will be used."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "041cf642",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np \n",
+    "from scipy.constants import pi\n",
+    "from scipy.integrate import ode\n",
+    "from scipy.interpolate import interp1d\n",
+    "import random\n",
+    "import matplotlib.pyplot as plt\n",
+    "import seaborn as sns "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "50e51ea6",
+   "metadata": {},
+   "source": [
+    "The constans that we will used in the following:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "73b1497e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "c = 3e10  # Speed of light in cm/s\n",
+    "G = 6.67428e-8    # Gravitational constant in cm^3/g/s^2 or dyne cm^2/g^2\n",
+    "Msun = 1.989e33   # Solar mass in grams\n",
+    "dyncm2_to_MeVfm3 = 1. / (1.6022e33)  # Conversion factor from dyn/cm^2 to MeV/fm^3\n",
+    "gcm3_to_MeVfm3 = 1. / (1.7827e12)    # Conversion factor from g/cm^3 to MeV/fm^3\n",
+    "oneoverfm_MeV = 197.33    # Conversion factor for fm to MeV"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ea9e94a3",
+   "metadata": {},
+   "source": [
+    "To obtain the corresponding maximum mass  index for each posterior EOS parameter."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6cb3f174",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def EofP(P, epsgrid, presgrid):\n",
+    "   idx = np.searchsorted(presgrid, P)\n",
+    "   if idx == 0:\n",
+    "       eds = epsgrid[0] * pow(P / presgrid[0], 3. / 5.)\n",
+    "   if idx == len(presgrid):\n",
+    "       eds = epsgrid[-1] * pow(P / presgrid[-1], 3. / 5.)\n",
+    "   else:\n",
+    "       ci = np.log(presgrid[idx] / presgrid[idx-1]) /np.log(epsgrid[idx] /epsgrid[idx-1])\n",
+    "       eds = epsgrid[idx-1] * pow(P / presgrid[idx-1], 1. / ci)\n",
+    "   return eds\n",
+    "\n",
+    "\n",
+    "def EofP(P, epsgrid, presgrid):\n",
+    "   idx = np.searchsorted(presgrid, P)\n",
+    "   if idx == 0:\n",
+    "       eds = epsgrid[0] * pow(P / presgrid[0], 3. / 5.)\n",
+    "   if idx == len(presgrid):\n",
+    "       eds = epsgrid[-1] * pow(P / presgrid[-1], 3. / 5.)\n",
+    "   else:\n",
+    "       ci = np.log(presgrid[idx] / presgrid[idx-1]) /np.log(epsgrid[idx] /epsgrid[idx-1])\n",
+    "       eds = epsgrid[idx-1] * pow(P / presgrid[idx-1], 1. / ci)\n",
+    "   return eds\n",
+    "   \n",
+    "   \n",
+    "def TOV(r, y, inveos):\n",
+    "    pres, m = y\n",
+    "    eps = inveos(pres)\n",
+    "    if np.abs(4.*pi*r**3. * pres) >1e10 :\n",
+    "        return ([1e27,1e27])\n",
+    "    else:\n",
+    "        dpdr = -(eps + pres) * (m + 4.*pi*r**3. * pres)\n",
+    "        dpdr = dpdr/(r*(r - 2.*m))\n",
+    "        dmdr = 4.*pi*r**2.0 * eps\n",
+    "    return np.array([dpdr, dmdr])\n",
+    "\n",
+    "def solveTOV(center_rho, energy_density, pressure):\n",
+    "    \"\"\"Solve TOV equation from given Equation of state in the neutron star \n",
+    "    core density range\n",
+    "\n",
+    "    Args:\n",
+    "        center_rho(array): This is the energy density here is fixed in main\n",
+    "        that is np.logspace(14.3, 15.6, 50)\n",
+    "        energy_density (array): Desity array of the neutron star EoS, in MeV/fm^{-3}\n",
+    "        Notice here for simiplicity, we omitted G/c**4 magnitude, so \n",
+    "        (value in MeV/fm^{-3})*G/c**4, could convert to the energy density we are\n",
+    "        using, please check the Test_EOS.csv to double check the order of magnitude.\n",
+    "        \n",
+    "        pressure (array): Pressure array of neutron star EoS, also in nautral unit\n",
+    "        with MeV/fm^{-3}, still please check the Test_EOS.csv, the conversion is \n",
+    "        (value in dyn/cm3)*G/c**4.\n",
+    "\n",
+    "    Returns:\n",
+    "        Mass (array): The array that contains all the Stars' masses, in M_sun as a \n",
+    "        Units.\n",
+    "        Radius (array): The array that contains all the Stars's radius, in km.\n",
+    "    \"\"\"\n",
+    "    c = 3e10\n",
+    "    G = 6.67428e-8\n",
+    "    Msun = 1.989e33\n",
+    "\n",
+    "    unique_pressure_indices = np.unique(pressure, return_index=True)[1]\n",
+    "    unique_pressure = pressure[np.sort(unique_pressure_indices)]\n",
+    "\n",
+    "    #Interpolate pressure vs. energy density\n",
+    "    eos = interp1d(energy_density, pressure, kind='cubic', fill_value='extrapolate')\n",
+    "\n",
+    "    #Interpolate energy density vs. pressure\n",
+    "    inveos = interp1d(unique_pressure, energy_density[unique_pressure_indices], kind='cubic', fill_value='extrapolate')\n",
+    "    \n",
+    "    Pmin = pressure[20]\n",
+    "    r = 4.441e-16\n",
+    "    dr = 10.\n",
+    "    center_rho = center_rho * G/c**2.\n",
+    "    \n",
+    "    pcent = eos(center_rho)\n",
+    "    P0 = pcent - (2.*pi/3.)*(pcent + center_rho) *(3.*pcent + center_rho)*r**2.\n",
+    "    m0 = 4./3. *pi *center_rho*r**3.\n",
+    "    stateTOV = np.array([P0, m0])\n",
+    "    \n",
+    "    sy = ode(TOV, None).set_integrator('dopri5')\n",
+    "    \n",
+    "    #have been modified from Irida to this integrator\n",
+    "    sy.set_initial_value(stateTOV , r).set_f_params(inveos)\n",
+    "    \n",
+    "    while sy.successful() and stateTOV[0]>Pmin:\n",
+    "        stateTOV = sy.integrate(sy.t+dr)\n",
+    "        dpdr, dmdr = TOV(sy.t+dr, stateTOV, inveos)\n",
+    "        dr = 0.46 * (1./stateTOV[1] * dmdr - 1./stateTOV[0]*dpdr)**(-1.)\n",
+    "\n",
+    "    return stateTOV[1]*c**2./G/Msun, sy.t/1e5  \n",
+    "\n",
+    "def TOV1(r, y, grids):\n",
+    "   pres, m = y\n",
+    "   epsgrid = grids[0]\n",
+    "   presgrid = grids[1]\n",
+    "   eps = EofP(pres, epsgrid, presgrid)\n",
+    "   dpdr = -(eps + pres) * (m + 4.*pi*r**3. * pres)\n",
+    "   dpdr = dpdr/(r*(r - 2.*m))\n",
+    "   dmdr = 4.*pi*r**2.0 * eps\n",
+    "   return np.array([dpdr, dmdr])\n",
+    "\n",
+    "\n",
+    "def solveTOV1(rhocent, eps, pres):\n",
+    "   eps = np.array(eps)\n",
+    "   pres = np.array(pres)\n",
+    "   Pmin = pres[20]\n",
+    "   r = 4.441e-16\n",
+    "   dr = 10.\n",
+    "   rhocent = rhocent * G/c**2.\n",
+    "   pcent = PofE(rhocent, eps, pres)\n",
+    "   P0 = pcent - (2.*pi/3.)*(pcent + rhocent) *(3.*pcent + rhocent)*r**2.\n",
+    "   m0 = 4./3. *pi *rhocent*r**3.\n",
+    "   stateTOV = np.array([P0, m0])\n",
+    "   sy = ode(TOV1, None).set_integrator('dopri5')\n",
+    "   par = np.vstack([eps, pres])\n",
+    "   sy.set_initial_value(stateTOV, r).set_f_params(par)\n",
+    "   while sy.successful() and stateTOV[0]>Pmin:\n",
+    "       stateTOV = sy.integrate(sy.t+dr)\n",
+    "       dpdr, dmdr = TOV1(sy.t+dr, stateTOV, par)\n",
+    "       dr = 0.46 * (1./stateTOV[1] * dmdr - 1./stateTOV[0]*dpdr)**(-1.)\n",
+    "\n",
+    "   return stateTOV[1]*c**2./G/Msun, sy.t/1e5\n",
+    "\n",
+    "\n",
+    "def compute_tov_properties(densities, eps_total, pres_total):\n",
+    "   max_length = len(densities)\n",
+    "   RFSU2R = np.zeros(max_length)\n",
+    "   MFSU2R = np.zeros(max_length)\n",
+    "   for i in range(max_length):\n",
+    "       MFSU2R[i], RFSU2R[i] = solveTOV1(densities[i], eps_total, pres_total)\n",
+    "       if i > 20 and MFSU2R[i] < MFSU2R[i - 1]:\n",
+    "           return MFSU2R[:i], RFSU2R[:i]  # Only return up to i, not the full arrays\n",
+    "\n",
+    "   return MFSU2R, RFSU2R\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0d2a3545",
+   "metadata": {},
+   "source": [
+    "Next, onw way is to sample the central energy density to obtain the posterior MR distribution. \n",
+    "\n",
+    "Another way is to sample the MR point directly from the MR relation. This is the same concept. \n",
+    "\n",
+    "The following code is an example of sampling the MR point directly using the posterior distribution of strangeon_matter_EOS."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2a4d8882",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Ms = []\n",
+    "Rs = []\n",
+    "\n",
+    "Nq=18\n",
+    "def Strangeon_compute_EOS(n, theta):\n",
+    "    epsilon, ns = theta\n",
+    "    A12 = 6.2\n",
+    "    A6 = 8.4 \n",
+    "    mq = 300  # Given constant value for mq\n",
+    "    \n",
+    "    sigma = np.sqrt(A6 / (2 * A12)) * (Nq / (3 * ns)) \n",
+    "    energy_density = 2 * epsilon * (A12 * sigma**4 * n**5 - A6 * sigma**2 * n**3) + n * Nq * mq\n",
+    "    pressure = 4 * epsilon * (2 * A12 * sigma**4 * n**5 - A6 * sigma**2 * n**3)\n",
+    "\n",
+    "    return energy_density * G / c**2 / gcm3_to_MeVfm3, pressure * G / c**4 / dyncm2_to_MeVfm3 \n",
+    "\n",
+    "\n",
+    "array_size = 50  # Replace array_size with the desired number of random numbers \n",
+    "\n",
+    "data = np.loadtxt('equal_weighted_post.txt', delimiter=' ',skiprows=1) \n",
+    "Rs = np.zeros(len(data)*50)\n",
+    "Ms = np.zeros(len(data)*50)\n",
+    "massradius=[]\n",
+    "for i in range(0,len(data)):\n",
+    "   theta = data[i][:2]  # The first two columns are epsilon/Nq, ns\n",
+    "   n_min = 3 * theta[1] / Nq  \n",
+    "   n_max = 0.16 * 8 * 3 / Nq    \n",
+    "   nbar_values = np.linspace(n_min, n_max, 10000)  # 10000 points between n_min and n_max  \n",
+    "   energy_densities, pressures = Strangeon_compute_EOS(nbar_values, theta)\n",
+    "   density = np.logspace(14.3, 15.6, 50)\n",
+    "   max_length =  len(density)\n",
+    "\n",
+    "   MFSU2R, RFSU2R = compute_tov_properties(density, energy_densities,pressures)\n",
+    "   mr = interp1d(MFSU2R, RFSU2R)\n",
+    "   index = len(MFSU2R) \n",
+    "   for j in range(0,50):\n",
+    "        Mpoint = max(MFSU2R)*random.random()\n",
+    "#A random mass Mpoint is sampled uniformly between 0 and the maximum mass in MFSU2R.\n",
+    "        Rpoint = mr(Mpoint) \n",
+    "        Ms[i+j] =  Mpoint\n",
+    "        Rs[i+j] = Rpoint\n",
+    "#The random mass and radius values are stored in the Ms and Rs arrays.\n",
+    "        massradius.append([Rs[i],Ms[i]])\n",
+    "        print(i,Mpoint,Rpoint) \n",
+    "        \n",
+    "plt.scatter(Rs, Ms, color='palevioletred')   \n",
+    "\n",
+    "np.savetxt('withNICERtwo_para_Mass.txt',Ms)  \n",
+    "np.savetxt('withNICERtwo_para_Radius.txt',Rs)\n",
+    "np.savetxt('withNICERtwo_para_MR.txt',massradius) \n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "03714bef",
+   "metadata": {},
+   "source": [
+    "Then, you can use the MR points to plot the MR contour at the 99.7% confidence level."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d34c80c5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "X =np.array(np.loadtxt('withNICERtwo_para_Radius.txt') )\n",
+    "Y=np.array(np.loadtxt('withNICERtwo_para_Mass.txt')) \n",
+    "R_l = []\n",
+    "M_l = []\n",
+    "MR = list(zip(X, Y))\n",
+    "MR = [m for m in MR if  0.001< m[0] <16]\n",
+    "R_l, M_l = zip(*MR)\n",
+    "   \n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(8, 6))\n",
+    "sns.kdeplot(x=R_l, y=M_l, shade=False,levels=[0.03], linestyles='--',linewidths=1.,\n",
+    "             cmap=None, alpha=0.9, colors=['green'], ax=ax)   \n",
+    "plt.scatter(R_l, M_l, color='palevioletred')  \n",
+    "\n",
+    "ax.set_xlim(6, 18)\n",
+    "ax.set_ylim(1, 5)\n",
+    "ax.set_xlabel(\"$R~(\\mathrm{km})$\")\n",
+    "ax.set_ylabel(r\"$M \\ (M_{\\odot}) $\")\n",
+    "ax.tick_params(direction='in', top='on', right='on', which='both')\n",
+    "ax.set_title('M-R Posterior Distribution with Contour Levels')\n",
+    "plt.savefig(\"MR_posterior.pdf\")\n",
+    "plt.show()\n",
+    "\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}