Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

How to make a triangular plot with different data sets, starting with a data set that contains fewer cosmological parameters in the sample? #116

Closed
JaelssonLima opened this issue May 8, 2024 · 7 comments

Comments

@JaelssonLima
Copy link

JaelssonLima commented May 8, 2024

If I do the plot

roots_data=[IDEM2_sample5, IDEM2_sample2, IDEM2_sample3, IDEM2_sample4,IDEM2_sample1]

works, but I do:

roots_data=[IDEM2_sample1, IDEM2_sample2, IDEM2_sample3, IDEM2_sample4,IDEM2_sample5]

It doesn't work, how to solve it?

The IDEM2_sample1, IDEM2_sample2 strings have fewer parameters in the sample than the other data.

Att. Jaelsson

 #IDEM2 PLOTS / JSL 

import getdist as gd
from getdist import loadMCSamples
from getdist import plots
from getdist.types import ResultTable
import matplotlib.pyplot as plt

# Configuração da fonte
font = {'family': 'serif', 'serif': ['Times New Roman']}
plt.rc('font', **font)
plt.rc('text', usetex=True)


# Load the Monte Carlo samples
root='/home/jaelsson/Documents/montepython_public_jsl/chains/artigo_jsl/idem2_cov'
IDEM2_sample1 = gd.loadMCSamples(root+'/Pantheon+cc+bao1_idem2/2024-03-17_200000_', settings={'ignore_rows': 0})
IDEM2_sample2 = gd.loadMCSamples(root+'/Pantheon+H0+cc+bao1_idem2/2024-02-15_100000_', settings={'ignore_rows': 0})
IDEM2_sample3 = gd.loadMCSamples(root+'/planck2018TTTEEE+lensing_idem2_cov2/2024-04-19_300000_', settings={'ignore_rows': 0})
IDEM2_sample4 = gd.loadMCSamples(root+'/planck2018TTTEEE+lensing+Pantheon+cc+BAO_idem2_cov/2024-04-19_300000_', settings={'ignore_rows': 0})
IDEM2_sample5 = gd.loadMCSamples(root+'/planck2018TTTEEE+lensing+Pantheon+H0+cc+BAO_idem2_cov/2024-04-19_300000_', settings={'ignore_rows': 0})


#params
################################
#omega_b = Omega_b * (h**2)
#h=H0/100.0  # h = H0 / (100 km/s/Mpc)
Omega_b=0.05 # fixo
################################

#1
p1=IDEM2_sample1.getParams()
gamma_lam=p1.gamma_lam
H0=p1.H0
Omega_m=p1.Omega_m 
M=p1.M
omega_b=(p1.H0/100.0)**2*Omega_b
omega_cdm=(p1.H0/100.0)**2*p1.Omega_cdm
IDEM2_sample1.addDerived(gamma_lam, name='gamma' , label=r'\gamma')
IDEM2_sample1.addDerived(H0, name='H_0' , label=r'H_0')
IDEM2_sample1.addDerived(Omega_m, name='Omega_{m0}' , label=r'\Omega_{m0}')
IDEM2_sample1.addDerived(M, name='{M}' , label=r'M')
IDEM2_sample1.addDerived(omega_b, name='omega_{b}' , label=r'\omega_{b}')
IDEM2_sample1.addDerived(omega_cdm, name='omega_{cdm}' , label=r'\omega_{c}')


#2
p2=IDEM2_sample2.getParams()
gamma_lam=p2.gamma_lam
H0=p2.H0
Omega_m=p2.Omega_m 
M=p2.M
omega_b=(p2.H0/100.0)**2*Omega_b
omega_cdm=(p2.H0/100.0)**2*p2.Omega_cdm
IDEM2_sample2.addDerived(gamma_lam, name='gamma' , label=r'\gamma')
IDEM2_sample2.addDerived(H0, name='H_0' , label=r'H_0')
IDEM2_sample2.addDerived(Omega_m, name='Omega_{m0}' , label=r'\Omega_{m0}')
IDEM2_sample2.addDerived(M, name='{M}' , label=r'M')
IDEM2_sample2.addDerived(omega_b, name='omega_{b}' , label=r'\omega_{b}')
IDEM2_sample2.addDerived(omega_cdm, name='omega_{cdm}' , label=r'\omega_{c}')



#3
p3=IDEM2_sample3.getParams()
gamma_lam=p3.gamma_lam
H0=p3.H0
Omega_m=p3.Omega_m 
#M=p3.M
#omega_b=(p3.H0/100.0)**2*Omega_b
omega_b=p3.omega_b/100.0
omega_cdm=p3.omega_cdm
#A_s=p3.A_s
#n_s=p3.n_s
sigma8=p3.sigma8
tau_reio=p3.tau_reio
IDEM2_sample3.addDerived(gamma_lam, name='gamma' , label=r'\gamma')
IDEM2_sample3.addDerived(H0, name='H_0' , label=r'H_0')
IDEM2_sample3.addDerived(Omega_m, name='Omega_{m0}' , label=r'\Omega_{m0}')
#IDEM2_sample3.addDerived(M, name='{M}' , label=r'M')
IDEM2_sample3.addDerived(omega_b, name='omega_{b}' , label=r'\omega_{b}')
IDEM2_sample3.addDerived(omega_cdm, name='omega_{cdm}' , label=r'\omega_{c}')
#IDEM2_sample3.addDerived(A_s, name='A_{s}' , label=r'10^{-9}A_{s } ')
#IDEM2_sample3.addDerived(n_s, name='n_{s}' , label=r'n_{s } ')
IDEM2_sample3.addDerived(sigma8, name='sigma{8}' , label=r'\sigma_{8} ')
IDEM2_sample3.addDerived(tau_reio, name='tau_{reio}' , label=r'\tau_{reio} ')


#4
p4=IDEM2_sample4.getParams()
gamma_lam=p4.gamma_lam
H0=p4.H0
Omega_m=p4.Omega_m 
M=p4.M
#omega_b=(p4.H0/100.0)**2*Omega_b
omega_b=p4.omega_b/100.0
omega_cdm=p4.omega_cdm
#A_s=p4.A_s
#n_s=p4.n_s
sigma8=p4.sigma8
tau_reio=p4.tau_reio
IDEM2_sample4.addDerived(gamma_lam, name='gamma' , label=r'\gamma')
IDEM2_sample4.addDerived(H0, name='H_0' , label=r'H_0')
IDEM2_sample4.addDerived(Omega_m, name='Omega_{m0}' , label=r'\Omega_{m0}')
IDEM2_sample4.addDerived(M, name='{M}' , label=r'M')
IDEM2_sample4.addDerived(omega_b, name='omega_{b}' , label=r'\omega_{b}')
IDEM2_sample4.addDerived(omega_cdm, name='omega_{cdm}' , label=r'\omega_{c}')
#IDEM2_sample4.addDerived(A_s, name='A_{s}' , label=r'10^{-9}A_{s } ')
#IDEM2_samplee4.addDerived(n_s, name='n_{s}' , label=r'n_{s } ')
IDEM2_sample4.addDerived(sigma8, name='sigma{8}' , label=r'\sigma_{8} ')
IDEM2_sample4.addDerived(tau_reio, name='tau_{reio}' , label=r'\tau_{reio} ')

#5
p5=IDEM2_sample5.getParams()
gamma_lam=p5.gamma_lam
H0=p5.H0
Omega_m=p5.Omega_m 
M=p5.M
#omega_b=(p5.H0/100.0)**2*Omega_b
omega_b=p5.omega_b/100.0
omega_cdm=p5.omega_cdm
#A_s=p5.A_s
#n_s=p5.n_s
sigma8=p5.sigma8
tau_reio=p5.tau_reio
IDEM2_sample5.addDerived(gamma_lam, name='gamma' , label=r'\gamma')
IDEM2_sample5.addDerived(H0, name='H_0' , label=r'H_0')
IDEM2_sample5.addDerived(Omega_m, name='Omega_{m0}' , label=r'\Omega_{m0}')
IDEM2_sample5.addDerived(M, name='{M}' , label=r'M')
IDEM2_sample5.addDerived(omega_b, name='omega_{b}' , label=r'\omega_{b}')
IDEM2_sample5.addDerived(omega_cdm, name='omega_{cdm}' , label=r'\omega_{c}')
#IDEM2_sample5.addDerived(A_s, name='A_{s}' , label=r'10^{-9}A_{s } ')
#IDEM2_sample5.addDerived(n_s, name='n_{s}' , label=r'n_{s } ')
IDEM2_sample5.addDerived(sigma8, name='sigma{8}' , label=r'\sigma_{8}')
IDEM2_sample5.addDerived(tau_reio, name='tau_{reio}' , label=r'\tau_{reio} ')



roots_data=[IDEM2_sample1, IDEM2_sample2, IDEM2_sample3, IDEM2_sample4,IDEM2_sample5]
#roots_data=[IDEM2_sample4, IDEM2_sample2, IDEM2_sample3, IDEM2_sample1,IDEM2_sample5]
#params = ['gamma','H_0','Omega_{m0}']
#roots_data=[IDEM2_sample3, IDEM2_sample4,IDEM2_sample5]
#params = ['H_0','A_s','n_s','sigma{8}', 'tau_reio','100theta_s','ln10^{10}A_s']
#params = ['H_0','Omega_{m0}','omega_{b}','omega_{cdm}']
#params = ['gamma','H_0','Omega_{m0}','omega_{b}','omega_{cdm}','M']


params = ['H_0','Omega_{m0}','omega_{b}','omega_{cdm}','M','A_s','n_s','sigma{8}','tau_{reio}','100theta_s','ln10^{10}A_s']


#g.settings.scaling_factor = 1
g = plots.get_subplot_plotter(width_inch=8)  # Initialize without chain_dir for single chain
g.triangle_plot(roots_data,
                params,
                #param_limits=param_limits,
                #nx=5,
                legend_labels=['Background ','Background+$H_0$','Planck','Background+Planck','Background+Planck+$H_0$'],
                #colors=['black','dodgerblue','orange','green','tomato'], colored_text=True,
                filled=True,
                #contour_lws=0.8
               )

# g.export('IDEM2_1.pdf', dpi=2000)

# #print(samples.getTable(limit=2).tableTex())
# print(ResultTable(ncol=1,results=roots_data,
#                   paramList=['gamma','H_0','Omega_{m0}','M','omega_{b}', 'omega_{cdm}','A_s','n_s','sigma{8}','tau_{reio}','100theta_s','ln10^{10}A_s'],
#                  limit=1).tableTex())
# print(ResultTable(ncol=1,results=roots_data,
#                   paramList=['gamma','H_0','Omega_{m0}','M','omega_{b}', 'omega_{cdm}','A_s','n_s','sigma{8}','tau_{reio}','100theta_s','ln10^{10}A_s'],
#                  limit=2).tableTex())
@cmbant
Copy link
Owner

cmbant commented May 8, 2024

At the moment, the first chain determines the parameters that are used. Happy to look at a PR if you want to try to implement more flexibility.

@JaelssonLima
Copy link
Author

Thank you for your attention.
Att. J.

@emiliobellini
Copy link

Hi,
I had the same need and I found a quick and dirty fix for that.
You can find it at
https://github.com/emiliobellini/getdist/blob/master/getdist/plots.py#L1784

It's all about modifying the ‎‎GetDistPlotter.get_param_array method to take a list of roots instead of a single one. A ParamNames object is created with the first root. And then the extra parameters of the other roots are appended.

My implementation works only for triangle plots (these were my needs), but generalising it should not be difficult
Emilio

@JaelssonLima
Copy link
Author

I'm grateful for the contribution, but in my case there is still a problem and the version that worked also presents the same problem now with this modification "paramNames size (45) does not match number of parameters in samples (43)". But I thank you in advance for your contribution.
Att
J.

@JaelssonLima
Copy link
Author

The following function merges what is common between the chains, but still does not solve the complete problem.
Att
J.

   #funcão de mescla o tem de comum entre as amostras / function of merging the common elements between the samples
    
    def get_param_array(self, roots, params: Union[None, str, Sequence] = None, renames: Mapping = None):
        """
        Gets an array of :class:`~.paramnames.ParamInfo` for named params
        in the given `root`.

        If a parameter is not found in `root`, returns the original ParamInfo if ParamInfo
        was passed, or fails otherwise.

        :param roots: The root name of the samples to use
        :param params: the parameter names (if not specified, get all)
        :param renames: optional dictionary mapping input names and equivalent names
                        used by the samples
        :return: list of :class:`~.paramnames.ParamInfo` instances for the parameters
        """
        if not isinstance(roots, list):
            roots = [roots]

        # Get parameter names from the first root
        if hasattr(roots[0], 'param_names'):
            names = roots[0].param_names
        elif hasattr(roots[0], 'paramNames'):
            names = roots[0].paramNames
        elif hasattr(roots[0], 'names'):
            names = ParamNames(names=roots[0].names, default=getattr(roots[0], 'dim', 0))
        else:
            names = self.param_names_for_root(roots[0])

        # Collect common parameters across all roots
        common_param_names = set(name.name for name in names.names)

        for root in roots[1:]:
            if hasattr(root, 'param_names'):
                names_part = root.param_names
            elif hasattr(root, 'paramNames'):
                names_part = root.paramNames
            elif hasattr(root, 'names'):
                names_part = ParamNames(names=root.names, default=getattr(root, 'dim', 0))
            else:
                names_part = self.param_names_for_root(root)

            current_param_names = set(name.name for name in names_part.names)
            common_param_names.intersection_update(current_param_names)

        # Filter out the common parameters
        common_names = [name for name in names.names if name.name in common_param_names]

        if params is None or len(params) == 0:
            return common_names

        if isinstance(params, str):
            if params in common_param_names:
                return [name for name in common_names if name.name == params]
            else:
                raise Exception(f"parameter name not found: {params}")
        else:
            filtered_params = [param for param in params if (param if isinstance(param, str) else param.name) in common_param_names]
            if not filtered_params:
                raise Exception("No common parameters found in the provided roots.")
            
            is_param_info = [isinstance(param, ParamInfo) for param in filtered_params]
            error = [not a for a in is_param_info]
            renames_from_param_info = {param.name: getattr(param, "renames", []) for i, param in enumerate(filtered_params) if is_param_info[i]}

            if renames:
                renames = mergeRenames(renames, renames_from_param_info)
            else:
                renames = renames_from_param_info
            
            params_names = [getattr(param, "name", param) for param in filtered_params]
            old = [(old if isinstance(old, ParamInfo) else ParamInfo(old)) for old in filtered_params]
            new_params = names.parsWithNames(params_names, error=error, renames=renames)
            
            return [new if new else old for new, old in zip(new_params, old)]

@cmbant
Copy link
Owner

cmbant commented May 22, 2024

This is my attempt, seems to work at least for a simple triangle plot: #117
it probably won't work if there are bounded parameters where parameters are not in the first root (?)

@JaelssonLima
Copy link
Author

I'm very grateful for the solution...it really works now
Att
J.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants