<jupyter_start><jupyter_text>Generating Scaling Relations with XGA This tutorial will explain how to generate XGA scaling relations objects, and go through some of their key useful features. Generating scaling relations that map difficult-to-measure parameters onto values that are easier to determine is very common in cluster science, with mass-temperature, luminosity-temperature, and luminosity-richness relations being just a few that spring to mind.As such I created this class of product, ScalingRelation, as well as several independent fitting method to generate them from arbitrary data. This part of XGA can almost be viewed as a submodule within the main module, as you do not have to use the source and sample structure to generate relations, the data can be supplied from an external source, and the fitting process will work exactly the same.<jupyter_code>import pandas as pd
from astropy.units import Quantity, pix
from astropy.cosmology import Planck15
import numpy as np

from xga.samples.extended import ClusterSample
from xga.sas import evselect_spectrum, emosaic
from xga.xspec import single_temp_apec
from xga.models.misc import power_law
from xga.relations.fit import scaling_relation_lira, scaling_relation_odr
from xga.relations.clusters import LT, Lλ, MT
from xga.products.relation import ScalingRelation<jupyter_output><empty_output><jupyter_text>Relation quality disclaimer All of the relations we generate here are for demonstration purposes only, they **should not** be used for any real scientific purposes.  Defining our sample I am again using clusters from the XCS-SDSS sample, they will be used to demonstrate the fitting and visualisation of scaling relations from an XGA ClusterSample object. I have read in their information from a file as I did not want to include a large table of irrelevant information in one of the cells of this tutorial. This time I will be using fourty randomly chosen clusters rather than four, as we could not reasonably hope to get a good fit with such a small sample.<jupyter_code>sample = pd.read_csv("modified_xcssdss_sample.csv", header="infer").sample(40)

srcs = ClusterSample(sample["RA_xcs"].values, sample["DEC_xcs"].values, sample["z"].values, 
                     sample["name"].values, r500=Quantity(sample["r500"].values, 'arcmin'), 
                     clean_obs=True, clean_obs_reg="r500", load_fits=False, psf_corr=False, 
                     richness=sample['richness'].values, richness_err=sample['richness_err'].values)<jupyter_output>Declaring BaseSource Sample: 100%|██████████| 40/40 [00:35<00:00,  1.13it/s]
Generating products of type(s) image: 100%|██████████| 3/3 [00:03<00:00,  1.19s/it]
Generating products of type(s) ccf: 100%|██████████| 95/95 [02:43<00:00,  1.72s/it] 
Generating products of type(s) expmap: 100%|██████████| 3/3 [00:38<00:00, 12.93s/it]
Generating products of type(s) image: 100%|██████████| 40/40 [01:03<00:00,  1.59s/it]
Generating products of type(s) expmap: 100%|██████████| 40/40 [00:40<00:00,  1.01s/it]
Setting up Galaxy Clusters:   0%|          | 0/40 [00:00<?, ?it/s]/home/dt237/code/PycharmProjects/XGA/xga/sources/general.py:99: UserWarning: XCSSDSS-31144 has not been detected in all region files, so generating and fitting products with the 'region' reg_type will not use all available data
  warnings.warn("{n} has not been detected in all region files, so generating and fitting products"
Setting up Galaxy Clusters:   2%|▎         | 1/40 [00:10<07:02, 10.83s/it]/home/dt237/code/PycharmProj[...]<jupyter_text>Measuring cluster properties Just as in our last tutorial 'Spectroscopy with XGA', we're going to generate spectra in the R$_{500}$ region for all the sources, and then fit them with a single temperature absorbed APEC model. This will give us the temperature and luminosity of all the clusters in the sample, which will allow us to build a scaling relation as a demonstration of XGA's capabilities:<jupyter_code># This command generates spectra for all the clusters in our sample, in their R500 regions
srcs = evselect_spectrum(srcs, 'r500')
# And this fits a single temperature APEC model to them
srcs = single_temp_apec(srcs, 'r500')<jupyter_output>Running XSPEC Fits: 100%|██████████| 40/40 [03:50<00:00,  5.76s/it]<jupyter_text>A scaling relation from a ClusterSample I've built convenience functions into the ClusterSample class, which make it very easy to generate certain types of common scaling relation. The only information that the method typically needs is the region for which spectra were generated and parameters measured, though plenty of other information (such as the model used to fit the data) can be supplied. The methods assume that, for a ClusterSample, we would have fit an absorbed APEC model, so that is the default.Once the method call is complete, an XGA scaling relation product is returned to the user. Note that it is not stored inside the sample object, it goes directly back to the user.The user is allowed to select the x and y normalisation values that should be used during the fit, and can easily change them from their standard values by passing them into keyword arguments. We have found that, in general, the best results are given when the normalisation is similar to the median value of the dataset.<jupyter_code># Making an Lx-Tx relation
lx_tx_relation = srcs.Lx_Tx('r500')

# Making an Lx-richness relation
lx_λ_relation = srcs.Lx_richness('r500')<jupyter_output>/home/dt237/code/PycharmProjects/XGA/xga/samples/base.py:242: UserWarning: There are no XSPEC fits associated with XCSSDSS-455
  warn(str(err))
/home/dt237/code/PycharmProjects/XGA/xga/samples/extended.py:425: UserWarning: One of the temperature uncertainty values for XCSSDSS-10401 is more than three times larger than the other, this means the fit quality is suspect.
  warn("One of the temperature uncertainty values for {s} is more than three times larger than "
/home/dt237/code/PycharmProjects/XGA/xga/samples/extended.py:436: UserWarning: There are no XSPEC fits associated with XCSSDSS-455
  warn(str(err))
/home/dt237/code/PycharmProjects/XGA/xga/relations/fit.py:71: UserWarning: 2 sources have NaN values and have been excluded
  warn("{} sources have NaN values and have been excluded".format(thrown_away))
/home/dt237/code/PycharmProjects/XGA/xga/relations/fit.py:71: UserWarning: 1 sources have NaN values and have been excluded
  warn("{} sources have NaN values and have been excluded[...]<jupyter_text>Viewing scaling relations Just as with other XGA products, there is a ``view()`` method built into the class, to easily create high quality visualisations of the relation. As we generated these relations ourselves, they also have data point associated with them. That is not the case if the scaling relation comes from literature, whether it is one of the ones built into XGA, or one you have defined yourself.There are several options that change what the output looks like, but one of the main options is the ability to change the figure size, which can be set using the ``figsize`` keyword argument:<jupyter_code>lx_tx_relation.view()

lx_λ_relation.view(figsize=(6, 6))<jupyter_output><empty_output><jupyter_text>Different fitting methods  I have implemented three different ways to fit these scaling relations, primarily to give the user some choice, though I would **highly recommend** that you use the implementation of the LIRA fitting package for R ([Sereno (2015)](https://academic.oup.com/mnras/article/455/2/2149/1111686)) - LIRA is an extremely capable fitting package, and produces excellent results. The downside is that it requires the installation of some optional external dependencies.The other methods available are scipy's simple non-linear least squares implementation (curve_fit), and scipy's orthogonal distance regression (ODR) implementation. ODR is the default method for ClusterSample scaling relation methods, as it does not require external dependencies, and seems to work quite well.I'm going to use a ClusterSample Lx_Tx relation as an example, and show how you can change the fitting method when you call the function - all this does is call the different fitting functions in xga.relations.fit:<jupyter_code># By setting the fit method to LIRA, we choose to use that instead of the default ODR
lx_tx_relation_lira = srcs.Lx_Tx('r500', fit_method='lira')
# And quickly viewing the relation
lx_tx_relation_lira.view()

# We repeat this process but with scipy's curve_fit
lx_tx_relation_curvefit = srcs.Lx_Tx('r500', fit_method='curve_fit')
lx_tx_relation_curvefit.view()<jupyter_output>/home/dt237/code/PycharmProjects/XGA/xga/samples/base.py:242: UserWarning: There are no XSPEC fits associated with XCSSDSS-455
  warn(str(err))
/home/dt237/code/PycharmProjects/XGA/xga/samples/extended.py:425: UserWarning: One of the temperature uncertainty values for XCSSDSS-10401 is more than three times larger than the other, this means the fit quality is suspect.
  warn("One of the temperature uncertainty values for {s} is more than three times larger than "
/home/dt237/code/PycharmProjects/XGA/xga/samples/extended.py:436: UserWarning: There are no XSPEC fits associated with XCSSDSS-455
  warn(str(err))
/home/dt237/code/PycharmProjects/XGA/xga/relations/fit.py:71: UserWarning: 2 sources have NaN values and have been excluded
  warn("{} sources have NaN values and have been excluded".format(thrown_away))
R[write to console]: module mix loaded<jupyter_text>ScalingRelation objects created using the LIRA fitting method actually have some extra functionality compared to other ScalingRelations. As LIRA is based around the Just Another Gibbs Sampler (JAGS) package, it is an MCMC fitting process, and as such returns parameter chains which it can be useful to visualise.The output chains are stored in the ScalingRelation on its declaration, and two methods were implemented to generated chain and corner plots:<jupyter_code># Here we visualise the parameter chains
lx_tx_relation_lira.view_chains()

# And this shows us a useful corner plot of the parameter space
lx_tx_relation_lira.view_corner()<jupyter_output><empty_output><jupyter_text>ScalingRelation product properties and making predictions These product objects are quite simple compared to some of the other products in XGA, but still contain a lot of useful information. The information can be broadly split into two groups; administrative information concerning when the relation was generated, who by etc. (especially useful for relations from literature), and that information concerning the actual properties of the relation, units, the normalisations applied etc.I will demonstrate a few of the properties using a classic relation from literature that has been built into XGA, the ([Arnaud et al. (2005)](https://ui.adsabs.harvard.edu/abs/2005A%26A...441..893A/abstract)) mass-temperature relation. As well as visualising the relation, we can also easily access the fit parameters that were used to define it, using the 'pars' method; the first column contains the values of all fit parameters, and the second column contains the uncertainty.<jupyter_code># Demonstrating the basic adminstrative properties
print(MT.arnaud_m500.year)
print(MT.arnaud_m500.author)
print(MT.arnaud_m500.doi, '\n')

# Viewing the fit parameters that were used to define this scaling relation
print(MT.arnaud_m500.pars)
# We can also easily view the x and y units
print(MT.arnaud_m500.x_unit, MT.arnaud_m500.y_unit)
# And the normalisations
print(MT.arnaud_m500.x_norm, MT.arnaud_m500.y_norm)<jupyter_output>2005
Arnaud et al.
10.1051/0004-6361:20052856 

[[1.71 0.09]
 [3.84 0.14]]
keV solMass
5.0 keV 100000000000000.0 solMass<jupyter_text>Once we have a scaling relation, we are also able to make predictions from it, using a built in method that takes an Astropy quantity (with the correct units for the x-axis) as an argument. The one thing you should remember is that you may need to account for E(z), depending on the relation you have fitted. For instance, the L$_{\rm{x}}$ - T$_{\rm{x}}$ relations produced by XGA automatically multiply by the E(z)$^{-1}$ of each cluster before fitting occurs. As such any values predicted from that relation have an E(z) dependance that must be accounted for to get a physical value:<jupyter_code># Defining an astropy quantity in units of keV, to demonstrate the predict functionality
test_temp = Quantity(6, 'keV')

# Getting a prediction of the cluster mass from that test temperature
MT.arnaud_m500.predict(test_temp)<jupyter_output><empty_output><jupyter_text>Viewing multiple scaling relations As I mentioned briefly in the 'XGA Products' tutorial, ScalingRelation objects have a useful little trick that allows you to easily view different scaling relations on the same axis. As long as two ScalingRelation instances have the same x and y units, you can simply add them together using the Python addition operator. Here I demonstrate this with the L$_{\rm{x}}$ - T$_{\rm{x}}$ relation we measured earlier, and an XCS-SDSS L$_{\rm{x}}$ - T$_{\rm{x}}$ relation imported from XGA:<jupyter_code># Adding them together is extremely easy
comb_lt_relation = lx_tx_relation + LT.xcs_sdss_r500_52

# Then we just call .view() like we would normally
comb_lt_relation.view()<jupyter_output>/home/dt237/code/PycharmProjects/XGA/xga/products/relation.py:631: UserWarning: Not all of these ScalingRelations have the same x-axis names.
  warn('Not all of these ScalingRelations have the same x-axis names.')
/home/dt237/code/PycharmProjects/XGA/xga/products/relation.py:638: UserWarning: Not all of these ScalingRelations have the same y-axis names.
  warn('Not all of these ScalingRelations have the same y-axis names.')<jupyter_text>Fitting scaling relations from external data I tried to think of an interesting dataset that I could generate a scaling relation from as an example, but came up empty, if you can think of an interesting data to use as an example send me an email!Instead of something interesting, I'm just going to extract some data from our ClusterSample object and show you how to generate a scaling relation manually, using functions imported from xga.relations.fit. We're going to measure a scaling relation that doesn't seem particularly useful, $\lambda$ - T$_{\rm{x}}$:<jupyter_code># Grabbing the temperatures (and temperature errors) from the ClusterSample object
temps = srcs.Tx('r500')
# Printing them on screen 
print(temps, '\n')

# And getting the richness values out as well (they only have a gaussian error)
richnesses = srcs.richness
print(richnesses)<jupyter_output>[[ 2.13213     0.20519163  0.39433514]
 [ 2.80067     0.22927479  0.24261076]
 [10.6814      0.20506005  0.29234747]
 [ 4.69087     0.39670878  0.41071715]
 [        nan         nan         nan]
 [ 6.21681     0.10400441  0.1041854 ]
 [ 4.06662     0.33423215  0.36426303]
 [ 8.04773     0.37034005  0.37111371]
 [ 6.01913     0.31320463  0.31522211]
 [ 3.51722     0.35075593  0.61433286]
 [ 8.34206     0.18331185  0.18349742]
 [ 1.03624     0.04787366  0.04563837]
 [ 3.35933     0.07213937  0.0746874 ]
 [ 2.92438     0.33646087  0.48980278]
 [ 4.1308      0.08558627  0.08634711]
 [ 6.94468     0.71801967  0.98840567]
 [ 2.7248      0.14950907  0.19408899]
 [ 5.89392     0.23684146  0.23756279]
 [ 1.24694     0.34040456  0.18172081]
 [ 1.5099      0.14780256  0.14409164]
 [ 7.24241     0.18409504  0.18429442]
 [ 3.49575     0.12190259  0.10648139]
 [ 3.13697     0.15849034  0.15939   ]
 [ 2.89694     0.36272877  0.39730179]
 [ 2.0819      0.17245655  0.27626099]
 [ 3.93636     0.0656564 [...]<jupyter_text>The scaling relation generation functions expect the data and uncertainties to be passed in separately, so here we're just going to split the data I just retrieved from my ClusterSample object. The scaling relation generating methods support data having either two uncertainties (in which case it will understand them as - and + errors), or just the one set of uncertainties. Part of the data preperation step inside of the scaling relation involves calculating an average uncertainty for a value if it is passed in with two uncertainties. This average uncertainty is what is used to fit the data.<jupyter_code># The first column contains the actual temperature values
x_data = temps[:, 0]
# The second and third columns are the uncertainties
x_errs = temps[:, 1:]

# Again the first column contains the richness values
y_data = richnesses[:, 0]
# However this data only has one uncertainty (rather than - and + uncertainties)
y_errs = richnesses[:, 1]<jupyter_output><empty_output><jupyter_text>Please make sure you are supplying the x and y data in the correct order, this function needs you to pass the y-data **first**.Also, when you fit a scaling relation with LIRA, it automatically uses a power-law in logspace as the model to fit with. This is not true of the other fitting methods, and you will need to supply a model for the fit to use; I suggest importing power_law from xga.models.misc and using that. There are another couple of differences between the LIRA function and other fitting methods, ODR and Curve Fit both allow you to set start parameters (though you don't have to), and the LIRA fitting function allows you to tell it how many chains should be run, and with how many steps.When deciding on what to name your x and y axis, please remember not to put units in the name, XGA will add them itself. Also bear in mind that these names will be passed into matplotlib, so some LaTeX syntax will work:<jupyter_code># The x values are temperatures, so selecting an arbitrary normalisation of 4keV
x_norm = Quantity(4, 'keV')
# The y values are richness, which has no unit, but still needs to be in an Astropy quantity
y_norm = Quantity(60, '')

# This uses a LIRA fit to generate the scaling relation
richness_tx_lira = scaling_relation_lira(y_data, y_errs, x_data, x_errs, y_norm, x_norm, 
                                         x_name=r'T$_{\rm{x}}$', y_name=r'$\lambda$')

# This uses ODR to generate a scaling relation, note that I have passed in a model function
richness_tx_odr = scaling_relation_odr(power_law, y_data, y_errs, x_data, x_errs, y_norm, x_norm, 
                                         x_name=r'T$_{\rm{x}}$', y_name=r'$\lambda$')
(richness_tx_lira+richness_tx_odr).view()<jupyter_output><empty_output><jupyter_text>You can also set a range in which you wish a scaling relation to be valid, when you generate it using one of these fitting functions. For instance I generate the same scaling relation we just created, but decide that it is only valid in the temperature range 2-6keV:<jupyter_code>richness_tx_lira_lim = scaling_relation_lira(y_data, y_errs, x_data, x_errs, y_norm, x_norm, 
                                             x_name=r'T$_{\rm{x}}$', y_name=r'$\lambda$', 
                                             x_lims=Quantity([2, 6], 'keV'))
richness_tx_lira_lim.view()<jupyter_output>/home/dt237/code/PycharmProjects/XGA/xga/relations/fit.py:71: UserWarning: 2 sources have NaN values and have been excluded
  warn("{} sources have NaN values and have been excluded".format(thrown_away))<jupyter_text>Defining scaling relations from literature Here I will demonstrate how you can find a relation in literature, and then implement it in a form usable by XGA. This particular scaling relation is actually already implemented in XGA, and is the Arnaud mass-temperature relation that I used in the demonstrations earlier.It really is as simple as filling out as many of the ScalingRelation arguments as you need. You should note that 'power_law' in this case is defined in the XGA models package, and is a function used to generate the fit line that is plotted. ScalingRelation will take its x and y units from the normalisations that are passed in on initialisation (here 5keV for the x axis, and 1e+14M$_{\odot}$ for the y-axis).Please don't include units in the values passed for ``x_name`` and ``y_name``, the ``view()`` method will add its own units, taken from the astropy normalisation quantities:<jupyter_code>mt_relation = ScalingRelation(np.array([1.71, 3.84]), np.array([0.09, 0.14]), power_law, Quantity(5, 'keV'),
                              Quantity(1e+14, 'solMass'), r"T$_{x}$", "E(z)$^{-1}$M$_{500}$",
                              relation_author='<NAME>.', relation_year='2005',
                              relation_doi='10.1051/0004-6361:20052856',
                              relation_name='Hydrostatic Mass-Temperature', x_lims=Quantity([1, 12], 'keV'))
mt_relation.view()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07709871745044695, 'fraction_numerical': 0.05183637776914108, 'mean_word_length': 4.565017572316842, 'pattern_counts': {'":': 0, '<': 63, '<?xml version=': 0, '>': 54, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '241582', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 6703, 'n_tokens_neox': 5850, 'n_words': 2545, 'path': 'docs/source/notebooks/tutorials/scaling_relations.ipynb', 'repo_name': 'DavidT3/XGA', 'size': '20590'}