I am trying to plot some data from a Markov Chain Monte Carlo run, using getDist, which (I assume) uses matplotlib under the hood for plotting. There is no error message or anything, and it produces a plot that is almost correct, expect for one crucial thing: subplots 0 and 3 have the opposite colour order than subplot 2, i.e. in subplots 0 and 3 data set 1 is dark blue and data set 2 is yellow, while in subplot 2, data set 1 is yellow and data set 2 is dark blue. Is there a way to avoid this, maybe assigning a fixed colour to a data set?

The code is the following:

# Import required libraries:
# - matplotlib.pyplot for plotting
# - cmcrameri.cm for perceptually uniform colormaps
# - getdist.plots for MCMC chain plotting
import matplotlib.pyplot as plt
from cmcrameri import cm
from getdist import plots

# Initialize the getdist subplot plotter
# - chain_dir: directory where MCMC chain files are stored
# - analysis_settings: ignore_rows to skip burn-in (first 20% of samples)
g = plots.get_subplot_plotter(
    chain_dir=r"/Users/klmba", analysis_settings={"ignore_rows": 0.2}
)

# Define the root names of the MCMC chains (file prefixes without extensions)
roots = [
    "Cobaya_mcmc_Run3_Planck_PP_SH0ES_DESIDR2_DoubleExp_tracking_uncoupled",
    "cobaya_iDM_20251230_dexp",
]

# Specify the parameters to plot in the triangle plot
params = ["H0", "s8h5"]  # Hubble constant and S8

# Extract a list of colors from the categorical bamako colourmap for coloring the chains
colours = [tuple(c) for c in cm.bamakoS.colors]

# Configure the plotter settings
# - solid_colors: set the colors used for filled contours and lines
g.settings.solid_colors = colours

# Generate the triangle plot showing:
# - 1D marginalized distributions on the diagonal
# - 2D contour plots on the off-diagonal
g.triangle_plot(
    roots,  # List of chain roots to include
    params,  # Parameters to plot
    filled=True,  # Fill the contour regions
    colors=colours,  # Colors for contour lines and fills
    contour_lws=3,  # Line width for contours
    legend_loc="lower left",  # Legend position
    figure_legend_outside=True,  # Place legend outside the plot area
)

# Add reference bands for observational data
# SH0ES 2020b measurement of H0: 73.2 ± 1.3 km/s/Mpc
g.add_x_bands(73.2, 1.3, ax=0)  # Vertical band on H0 1D plot
g.add_x_bands(73.2, 1.3, ax=2)  # Vertical band on H0 axis of 2D plot

# KiDS-1000 2023 measurement of S8: 0.776 ± 0.031
g.add_x_bands(0.776, 0.031, ax=3)  # Horizontal band on S8 1D plot
g.add_y_bands(0.776, 0.031, ax=2)  # Horizontal band on S8 axis of 2D plot

# Export the plot to default file format (usually PDF/PNG)
g.export()

# Display the plot in the output
plt.show()

It produces the following plots: