<jupyter_start><jupyter_text>Parameters<jupyter_code>KEY = 'APR4'
import constants as u
import tools
data = tools.Data(KEY)
print(data.dataset._dims['params'])<jupyter_output>sigma_delta=1.382*MeV/fm**2
C_C=0.8957
d_C=3.0
n_0=0.16*1./fm**3
e_0=-16.0*MeV
K_0=240.0*MeV
S_2=32.0*MeV
L_2=60.0*MeV
K_2=30.0*MeV
a=14.4383200971*MeV
alpha=0.56302798373
b=1.86063464319*MeV
beta=2.71424708856
mu_p0=-104.5*MeV
u_p=3.136
m_eff_m_p=0.8
E_c=651.439985687*MeV
E_max=1172.0729033599998*MeV
C_max=1.0<jupyter_text>Mass/Radius Relationships<jupyter_code>%pylab inline --no-import-all
import tools
data = tools.Data(KEY)
data.explore_parameters();<jupyter_output>Populating the interactive namespace from numpy and matplotlib<jupyter_text>Principal Component Analysis As described in the [`Equation of State.ipynb`](Equation of State.ipynb) notebook, our equation of state is parameterized by the following 18 parameters which have the following values chosen to roughly match the ALF4 equation of state as tabulated in [Read:2009]: We now present a principal component analysis.  Here one should generate a sample of binary neutron star systems characterized by their two masses $m_1$ and $m_2$ (in units of solar masses $M_\circ$) and their distance $d$ in Mpc from the earth.  We provided below a `PopulationModel` which will generate a Gaussian population with specified mean and standard deviation for each parameter, but one can simply input any list of `(m_1, m_2, d)` tuples.From this sample, we plot the eigenvectors and eigenvalues of the combined Fisher information matrix $\mat{F}$ that would result from LIGO observations of assuming signal-to-noise ratios compatable with the Einstein Telescope.  The matrix $\mat{F}$ characterizes the relative errors of the various EoS parameters $p_a$ as follows:$$  \sum_{ab}\frac{\delta p_a}{p_a} F_{ab} \frac{\delta p_b}{p_b}    = \sum_{abi}\frac{\delta p_a}{p_a} U_{ai}d_i U_{bi} \frac{\delta p_b}{p_b}  = \sum_{i}(\delta\xi_i)^2 d_i \leq 1, \qquad  \delta\xi_i = \sum_a \frac{\delta p_a}{p_a} U_{ai}, \qquad  \xi_i = \sum_{a}U_{ai}\ln{p_a}.$$These covariance ellipsoids corresponds to the 1-$\sigma$ variation assuming that all the parameters variances are well described by Gaussian errors.  By diagonalizing $\mat{F} = \mat{U} \cdot \diag(\mat{d})\cdot\mat{U}^\dagger$ we obtain independent constraints on each of the the principal components $\xi_i$:$$   \abs{\delta\xi_i} \leq \sqrt{d_i^{-1}} = \sigma_i.$$*Note: in the following, we use the notation $\sigma_i = 1/\sqrt{d_i}$ to represent the constraint imposed by the $i$'th most significant component.*We now plot the various eigenvalues and display the correspond constraints as a percentage error $100/\sqrt{d_i}$:<jupyter_code>np.random.seed(1)
population_model = tools.PopulationModel(
    m1=1.2, m2=1.5, 
    distance=40, constant_distance=True)
pca = tools.PCA_Explorer(data)
pca.plot_PCA(population_model.get_samples(1), significance=50)
np.random.seed(2)
population_model = tools.PopulationModel(m1=1.2+0.2j, m2=1.5+0.2j, distance=100)
display(pca.plot_PCA(population_model.get_samples(100), significance=10))<jupyter_output><empty_output><jupyter_text>Here we have instead assumed that we have 100 observations from a population of objects Gaussian distributions of $m_1 = 1.2(2)M_\circ$, $m_2=1.5(2)M_\circ$ uniformly distributed within a sphere of radius $d=100$Mpc.  We now see that a variety of parameter combinations are constrained.  Parameter Information Here we test the hypothesis that each star gives essentially only one relevant principal component.  We can check this by showing the distribution of the principal component eigenvalues over the range of masses.<jupyter_code>d, U = pca.get_PCA(); sigmas = 1./np.ma.sqrt(d)

plt.figure(figsize=(13., 5.))
for n, gs in enumerate(GridSpec(1,2)):
    ax = plt.subplot(gs)
    ax.set_aspect(1)
    plt.pcolormesh(data.M/u.M0, data.M/u.M0, 1./sigmas[..., -1-n].T)
    cb = plt.colorbar(label=r'$\sigma_{}$'.format(n))
    ticks_ = np.array(cb.get_ticks())
    cb.set_ticks(ticks_)
    cb.set_ticklabels(["{:.2g}%".format(_sigma) for _sigma in 100.0/ticks_])
    plt.xlabel('$m_1$ [M0]'); plt.ylabel('$m_2$ [M0]')
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Notice that the relative constraint provided by the dominant principal component - even in the worst case - is 3 orders of magnitude larger than that of the second component.  Thus, it is a good approximation to consider only the dominant principal components.  Independent information results from averaging over different masses which have different eigenvectors for the dominant principal component. A related question is: which combinations of masses will provide the most information about a given parameter?  We can obtain a qualitative estimate of this by looking at the diagonal entries.  *(The maximum $\sqrt{F}$ is shown in the title of each plot with larger values indicating more information.)*<jupyter_code>masses = data.M
params = data.params
Np = len(params)
F = data.dataset.F

gs = GridSpec(3, 6)
plt.figure(figsize=(15, 8))
z_max = [np.sqrt(F.data[:,:,n,n]).max() for n in range(Np)]
inds = reversed(np.argsort(z_max))
for _n, n in enumerate(inds):
    ax = plt.subplot(gs[_n])
    z = np.sqrt(F.data[:,:,n,n].T)
    plt.pcolormesh(masses, masses, z)
    plt.title("{} ({})".format(params[n], int(z.max())))
    
    ax.set_aspect(1)
    
plt.tight_layout()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12015432665809296, 'fraction_numerical': 0.040602608855410616, 'mean_word_length': 5.193401592718999, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '567926', 'license': 'MIT', 'n_tokens_mistral': 1971, 'n_tokens_neox': 1827, 'n_words': 569, 'path': 'APR4.ipynb', 'repo_name': 'mforbes/binder_tov_explorer', 'size': '5443'}