{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## RMF EOS inference pipeline\n", "\n", "This is an example notebook about how to use our tools to analysis a observation/nuclear constraint on neutron star equation of stat. \n", "\n", "Here in this notebook, we are using a RMF EoS model" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#Package we need:\n", "import InferenceWorkflow.BayesianSampler as sampler\n", "import InferenceWorkflow.Likelihood as likelihood\n", "import InferenceWorkflow.prior as prior\n", "import math\n", "import numpy as np\n", "\n", "from TOVsolver.constant import oneoverfm_MeV, m_rho, m_w,G,c\n", "import TOVsolver.main as main\n", "import EOSgenerators.crust_EOS as crust\n", "import EOSgenerators.fastRMF_EoS as RMF" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We need to define the crust, this crust and interface is unversial for all the equation of state here, we just change the core part of equation of state and that is come from the RMF computation. Crust as below, \"Tolos_crust_out.txt\" is BPS crust model. Then as in README, since that file energy_density and presurssure all in MeV/fm3, so we need some conversion. then just call \n", "\n", " ```sh\n", " crust.PolyInterpolate\n", " ```\n", "\n", "To finish the interface and make it ready to connect with core part" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#Tolos_crust_out = np.loadtxt('Tolos_crust_out.txt', delimiter=' ')\n", "Tolos_crust_out = np.loadtxt('Tolos_crust_out.txt', delimiter=None, comments='#', usecols=(0, 1, 2, 3, 4))\n", "eps_crust_out = Tolos_crust_out[:,3] * G / c**2\n", "pres_crust_out = Tolos_crust_out[:,4] * G / c**4\n", "\n", "eps_crust, pres_crust = crust.PolyInterpolate(eps_crust_out, pres_crust_out)\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Set up prior\n", "\n", "Next step, we need to set up the prior, first use parameters array to specify the variable name, should consistent with what you need to call them.\n", "\n", "Define a prior transform function to define prior. Cube are set of random number from 0 to 1. This prior setting is standard set-up of UltraNest package, since we are using UltraNest to do nest-sampling. We provided \n", "\n", "\"normal_Prior\" and \"flat_prior\"\n", "\n", "two options call from prior. Here then the Parameters prior should all set\n", "\n", "------------------\n", "\n", "However, since we are doing Equation of state Inference from mass radius of neutron star measurement. The center density of the star should be also sampled. Otherwise will be a partially-defined prior, did not span all parameters space, and proved to be different with full-scope inference.\n", "\n", "This request as randomly generate a density from a EoS range, however, this process is not that trivial, since we need to determine the upper limit of the central density of neutron star --- different equation of state will predict different upper bound, so here we need to use the prior-setting EoS parameters computing the EOS by\n", "\n", "```sh\n", "RMF.compute_EOS\n", "```\n", "\n", "Compute out EOS, put into\n", "\n", "```sh\n", "main.OutputMR\n", "```\n", "\n", "find out Mass Radius of this equation of state, find out the last stable point of this equation of state.(first mass points that give the direvative to be negative)\n", "\n", "found out that index by len() function, then reset this max_d to be upper limit of this density range." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "parameters = ['g_sigma', 'g_omega','g_rho', 'kappa', 'lambda_0', 'zeta', 'Lambda_w','d1']\n", "# for two or more MR measurements, define d2 or more depend on complexity.\n", "# parameters = ['g_sigma', 'g_omega','g_rho', 'kappa', 'lambda_0', 'zeta', 'Lambda_w','d1', 'd2']\n", "\n", "def prior_transform(cube):\n", " params = cube.copy()\n", " params[0] = math.sqrt(prior.normal_Prior(107.5, 7.5,cube[0]))\n", " params[2] = math.sqrt(prior.flat_prior(75,210,cube[2]))\n", " params[1] = math.sqrt(prior.flat_prior(150,210,cube[1]))\n", " params[3] = prior.normal_Prior(2.525/oneoverfm_MeV, 1.525/oneoverfm_MeV,cube[3])\n", " params[4] = prior.normal_Prior(0.0045, 0.0205,cube[4])\n", " params[5] = prior.flat_prior(0,0.04,cube[5])\n", " params[6] = prior.flat_prior(0,0.045,cube[6])\n", " \n", " g_sigma = params[0]\n", " g_omega = params[1]\n", " g_rho = params[2]\n", "\n", " kappa = params[3]\n", " lambda_0 = params[4]\n", " zeta = params[5]\n", " Lambda_w = params[6]\n", " m_sig = 495 / oneoverfm_MeV\n", "\n", " theta = np.array([m_sig, m_w, m_rho, g_sigma, g_omega, g_rho, kappa,\n", " lambda_0, zeta, Lambda_w])\n", " \n", " ep, pr = RMF.compute_EOS(eps_crust, pres_crust, theta)\n", "\n", " eps_total = np.hstack((eps_crust,ep))\n", " \n", " pres_total = np.hstack((pres_crust,pr))\n", " \n", " RFSU2R = []\n", " MFSU2R = []\n", " density = np.logspace(14.3, 15.6, 50)\n", " if all(x 20 and MR[1][i] - MR[1][i-1]< 0:\n", " break\n", " if len(MFSU2R)==False:\n", " params[7] = 0\n", " # params[8] = 0\n", " # this line for showing how to add one more observation\n", " else:\n", " max_index = len(MFSU2R)\n", " max_d = np.log10(density[max_index-1])\n", " params[7] = 14.3 + (max_d - 14.3) * cube[7]\n", " # params[8] = 14.3 + (max_d - 14.3) * cube[8]\n", " # this line for showing how to add one more observation\n", " return params" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### Set up likelihood\n", "\n", "We need to set up a likelihood, Using standard definition way of UltraNest, that is below.\n", "\n", "Here the likelihood is generated from a simulated mass radius measurement, which is $M = 1.4 M_{\\odot}$ and $R = 13$ km, With a 5% Mass radius measurement uncertainty, so here \n", "\n", " ```sh\n", " likelihood.MRlikihood_Gaussian\n", " ```\n", "function will be use for our likelihood, please check [likelihood.MRlikihood_Gaussian](https://github.com/ChunHuangPhy/EoS_inference/blob/main/InferenceWorkflow/Likelihood.py) to see the original code, and more choice of likelihood.\n", "eg:\n", "1. If we have some real mass-radius measurements, say PSR J0030 or PSR J0740, come from NICER, a KDE kernel could be trained to feed into \n", "\n", " ```sh\n", " likelihood.MRlikihood_kernel(eps_total,pres_total,x,d1)\n", " ```\n", "set the KDE kernel as a input for this function\n", "\n", "2. If we gain measurement from radio-timing, say only measure the neutron star mass, then\n", "\n", " ```sh\n", " likelihood.Masslikihood_Gaussian(eps_total,pres_total,x,d1)\n", " ```\n", "Which will give the likelihood from single mass measurement, x is the parameters of that measurement, you should specify where this measurement mass is located and what is the sigma width of this mass measurement\n", "\n", "3. If we have nuclear measurements, and want to constrain this RMF model by nuclear properties like K(The Incompressibility of nuclear matter),J ( the symmetry energy at saturation density) and L( the slope of symmetry energy at saturation density). You can choose:\n", "\n", " ```sh\n", " likelihood.Kliklihood(theta,K_low,K_up)\n", " likelihood.Jliklihood(theta,K_low,K_up)\n", " likelihood.Lliklihood(theta,K_low,K_up)\n", " ```\n", "We are defaulting a hard-cut flat constrain, so if you don't like this default hard cut, also could define the likelihood by youself with similiar style.\n", "\n", "4. If we have a Tidal measurements from Gravitational wave detector, we can use it to do constraint:\n", "\n", " ```sh\n", " likelihood.TidalLikihood_kernel(eps_total,pres_total,x,d1)\n", " ```\n", "Where x is sampled distribution from real measurements, the standard is \n", "\n", "kernel, chrip = x, \n", "\n", "where the kernel is a whole set sampling from GW event, that is [chrip mass, M2/M1, tidal of M1, tidal of M2] four quantities. Chrip is the single smapling that comes only the chrip mass sampling.\n", "\n", "\n", "\n", "5. To calculate the log-likelihood for the pure neutron matter (PNM) equation of state (EoS) using chiral effective field theory (chiEFT) constraints, the function\n", " \n", " `likelihood.chiEFT_PNM(EoS_PNM, type=\"Gaussian\", contraint_quantity=\"e\", enlargement=0)`\n", "\n", "EoS_PNM is a 3D array of number density (rho), energy density (e) and pressure (p) with unit fm^-3, MeV.fm^-3 and MeV.fm^-3. \n", "\n", "compares PNM EoS data with chiEFT-derived constraints for either energy per neutron \\( E/N \\) or pressure \\( p \\),\n", "depending on the specified `contraint_quantity` parameter (\"e\" for energy or \"p\" for pressure). \n", "\n", "The likelihood model can be either \"Gaussian\" or \"Super Gaussian\" with an optional `enlargement` factor for the Super Gaussian to flatten its peak. e.g., 0.05 for 5% enlargement. \n", "\n", "The function calculates the sum of log-likelihoods over densities of 0.08, 0.12, and 0.16 fm\\(^-3\\), utilizing energy constraints from Huth et al. (Nature, 2022) and pressure constraints from Hebeler et al. (ApJ, 2013). This enables a statistical comparison of the PNM EoS data against chiEFT-derived reference values." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import scipy.stats as stats\n", "def likelihood_transform(theta):\n", " # This is a demonstration code for only introduce one constraint from one mass-radius observation.\n", " # Could be very easy to implement more constraint from nuclear quantity, since that do not need to\n", " # sample more central density of real neutron star. If user want to expand to two mass radius measurement \n", " # the code could be:\n", " \n", " # g_sigma, g_omega,g_rho, kappa, lambda_0, zeta, Lambda_w, d1, d2 = theta\n", " g_sigma, g_omega,g_rho, kappa, lambda_0, zeta, Lambda_w, d1 = theta # comment this line if you need two measuremnts.\n", " \n", " ####################################################################################################################\n", " ############ This is the block to compute out all the EoS you need based on your parameters#########################\n", " m_sig = 495 / oneoverfm_MeV\n", " m_w = 3.96544\n", " m_rho = 3.86662\n", " theta1 = np.array([m_sig, m_w, m_rho, g_sigma, g_omega, g_rho, kappa, lambda_0, zeta, Lambda_w])\n", " ep, pr = RMF.compute_EOS(eps_crust, pres_crust, theta1)\n", "\n", " eps_total = np.hstack((eps_crust,ep))\n", " pres_total = np.hstack((pres_crust,pr))\n", " ####################################################################################################################\n", " \n", " # probMRgaussian1 = likelihood.MRlikihood_Gaussian(eps_total,pres_total,(1.4,13,0.07,0.65),d1)\n", " # probMRgaussian2 = likelihood.MRlikihood_Gaussian(eps_total,pres_total,(2.0,12,0.1,0.6),d2)\n", " # probMR = probMRgaussian1 + probMRgaussian2\n", " \n", " # Same could be extended to more distributions. Notice the prior definition should be change accordingly\n", " # by define more density parameters like here d2.\n", " \n", " #1. This line is to compute MR likelihood from a Simulated MR measurement:\n", " \n", " \n", " probMRgaussian = likelihood.MRlikihood_Gaussian(eps_total,pres_total,(1.4,13,0.07,0.65),d1)\n", " \n", " #2. This is a block that constrain from given real MR measurement, say J0030:\n", " #J0030 = numpy.loadtxt('data/PST_equal_sampled_MR.txt', delimiter=' ')\n", " #J30R_list, J30M_list = zip(*J0030)\n", " #J30R_list = numpy.array(J30R_list).T \n", " #J30M_list = numpy.array(J30M_list).T\n", " #Rmin = J30R_list.min()\n", " #Rmax = J30R_list.max()\n", " #Mmin = J30M_list.min()\n", " #Mmax = J30M_list.max()\n", " #X3, Y3 = numpy.mgrid[Rmin:Rmax:500j, Mmin:Mmax:100j]\n", " #positions = numpy.vstack([X3.ravel(), Y3.ravel()])\n", " #values = numpy.vstack([J30R_list, J30M_list])\n", " #kernel3 = stats.gaussian_kde(values)\n", " #probMRJ0030 = likelihood.MRlikelihhood_kernel(eps_total,pres_total,kernel3,d1)\n", " \n", " #3. This is to compute the constraint from experiment of nuclearmatter\n", " # 250" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import corner\n", "\n", "figure = corner.corner(samples,labels=[ r\"$g_{\\sigma}$\", r\"$g_{\\omega}$\",r\"$g_{\\rho}$\",r\"$\\kappa$\", r\"$\\lambda_0$\",r\"$\\zeta$\",r\"$\\Lambda_w$\",\"d1\",\"d2\"],\n", " smooth=0.9,\n", " label_kwargs=dict(fontsize=22),\n", " title_kwargs=dict(fontsize=22),\n", " quantiles=[0.16, 0.84],\n", " levels=(1 - np.exp(-0.5), 1 - np.exp(-2), 1 - np.exp(-9 / 2.)),\n", " plot_density=False,\n", " plot_datapoints=False,\n", " fill_contours=True,\n", " show_titles=True,\n", " max_n_ticks=3,\n", " title_fmt= '.2f')\n", " #,range=[(10,80),(0.56,0.64),(0.10, 0.2),(150,600),(40,80),(50,75)])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.12.2" }, "vscode": { "interpreter": { "hash": "e4481115b400f107b26d360e6549f546bb0e8bc1af70e4e66085bfa77a017a39" } } }, "nbformat": 4, "nbformat_minor": 4 }