{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9662b47d-c992-44f9-ba79-69b3e642c601",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "from scipy.integrate import IntegrationWarning\n",
    "\n",
    "# Ignore integration warnings\n",
    "warnings.simplefilter(\"ignore\", IntegrationWarning)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6baaa92a-c182-451e-99dc-2284fc995a18",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Physical baryon density parameter (\\Omega_b h^2) from Planck\n",
    "obh2pl = 0.02236\n",
    "\n",
    "# Matter density parameter (\\Omega_m) from Planck\n",
    "ompl = 0.3166\n",
    "\n",
    "# Hubble constant (H_0) in km s^-1 Mpc^-1 from Planck\n",
    "H0pl = 67.27\n",
    "\n",
    "# Dimensionless Hubble parameter (h) from Planck\n",
    "h0pl = H0pl / 100\n",
    "\n",
    "# Dimensionless Hubble parameter (h) from Type Ia Supernovae\n",
    "h0riess = 0.7403\n",
    "\n",
    "# Cosmic Microwave Background (CMB) temperature in Kelvin\n",
    "Tcmb = 2.7255\n",
    "\n",
    "# Speed of light in vacuum in km s^-1\n",
    "c = 299792.458\n",
    "\n",
    "# Hubble distance scale (c / 100) in Mpc\n",
    "cH0 = 2997.92458"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "902662cc-67ff-469e-8ce9-e8f25d2e34d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import mpmath\n",
    "import numpy as np\n",
    "# w(a) = w_0 + w_a * (1 - a)^n, i.e., nCPL model\n",
    "# For n=1 we get the usual CPL model\n",
    "def w(a, w0, wa, n):\n",
    "    return w0 + wa * (1 - a)**n\n",
    "\n",
    "# f_DE(a) for the nCPL model\n",
    "# See Eq. (4.2) on page 12 from https://arxiv.org/pdf/1603.02164.pdf\n",
    "def f(a, w0, wa, n):\n",
    "    # Using mpmath for special mathematical functions\n",
    "    harmonic_n = mpmath.harmonic(n)\n",
    "    \n",
    "    # mpmath.hyper computes the generalized hypergeometric function pFq\n",
    "    hyper_term = mpmath.hyper([1, 1, 1 - n], [2, 2], a)\n",
    "    \n",
    "    # Calculate the two main terms of the equation\n",
    "    term1 = a**(-3 * (1 + w0 + wa))\n",
    "    term2 = mpmath.exp(-3 * wa * harmonic_n + 3 * a * n * wa * hyper_term)\n",
    "    \n",
    "    # Return as a standard Python float\n",
    "    return float(term1 * term2)\n",
    "\n",
    "# Define z_equality and a_equality in order to avoid some problems.\n",
    "# Tcmb is added as a parameter with a default value of 2.725 K.\n",
    "def zeq(om, h, Tcmb=2.725):\n",
    "    return 2.5*10**4 * om * (h**2) * (Tcmb / 2.7)**(-4)\n",
    "\n",
    "def aeq(om, h, Tcmb=2.725):\n",
    "    return 1.0 / (1.0 + zeq(om, h, Tcmb))\n",
    "\n",
    "#humble parameter\n",
    "def H(a, om, w0, wa, n, h):\n",
    "    # matter and radiation term\n",
    "    term_m_r = (a**-3) * om * (1 + aeq(om, h) / a)\n",
    "    # dark energy term\n",
    "    term_de = (1 - om * (1 + aeq(om, h))) * f(a, w0, wa, n)\n",
    "    \n",
    "    return 100 * h * np.sqrt(term_m_r + term_de)\n",
    "\n",
    "#sound speed\n",
    "def cs(a, obh2):\n",
    "    R = 31500 * obh2 * (Tcmb / 2.7)**-4 * a\n",
    "    return c / np.sqrt(3 * (1 + R))\n",
    "\n",
    "#redshift at decoupling\n",
    "def zcmb(om, obh2, h):\n",
    "    omh2 = om * h**2\n",
    "    g1 = (0.0783 * obh2**-0.238) / (1 + 39.5 * obh2**0.763)\n",
    "    g2 = 0.560 / (1 + 21.1 * obh2**1.81)\n",
    "    return 1048 * (1 + 0.00124 * obh2**-0.738) * (1 + g1 * omh2**g2)\n",
    "\n",
    "#drag redshift\n",
    "def zdrag(om, obh2, h):\n",
    "    omh2 = om * h**2\n",
    "    b1 = 0.313 * omh2**-0.419 * (1 + 0.607 * omh2**0.674)\n",
    "    b2 = 0.238 * omh2**0.223\n",
    "    \n",
    "    term1 = 1291 * omh2**0.251 / (1 + 0.659 * omh2**0.828)\n",
    "    term2 = 1 + b1 * obh2**b2\n",
    "    return term1 * term2\n",
    "\n",
    "from scipy.integrate import quad\n",
    "\n",
    "#comoving sound horizon at drag redshift\n",
    "def rs(ze, om, obh2, w0, wa, n, h):\n",
    "    lower_limit = 10**-11\n",
    "    upper_limit = 1.0 / (1.0 + ze)\n",
    "    \n",
    "    def integrand(x):\n",
    "        return cs(x, obh2) / (x**2 * H(x, om, w0, wa, n, h))\n",
    "    \n",
    "    \n",
    "    result, error = quad(integrand, lower_limit, upper_limit)\n",
    "    \n",
    "    return result\n",
    "    \n",
    "#luminosity distance\n",
    "from scipy.integrate import solve_ivp\n",
    "from functools import lru_cache\n",
    "\n",
    "@lru_cache(maxsize=None)\n",
    "def DLsol(om, w0, wa, n, h):\n",
    "    def dy_dz(zz, y):\n",
    "        a = 1.0 / (1.0 + zz)\n",
    "        return 1.0 / (H(a, om, w0, wa, n, h) / (100.0 * h))\n",
    "    \n",
    "    sol = solve_ivp(dy_dz, [0, 1300], [0.0], dense_output=True)\n",
    "    return sol.sol\n",
    "\n",
    "def DL(z, om, w0, wa, n, h):\n",
    "    sol_func = DLsol(om, w0, wa, n, h)\n",
    "    y_val = sol_func(z)[0]\n",
    "    dL_val = y_val * (1.0 + z)\n",
    "    return (c / (100.0 * h)) * dL_val\n",
    "    \n",
    "#angular diameter distance\n",
    "def DA(z, om, w0, wa, n, h):\n",
    "    return DL(z, om, w0, wa, n, h) / (1.0 + z)**2\n",
    "    \n",
    "#volume averaged distance\n",
    "def Dv(zbao, om, w0, wa, n, h):\n",
    "    dl_term = DL(zbao, om, w0, wa, n, h) / (1.0 + zbao)\n",
    "    h_val = H(1.0 / (1.0 + zbao), om, w0, wa, n, h)\n",
    "    \n",
    "    return ((dl_term**2) * (c * zbao) / h_val)**(1.0 / 3.0)\n",
    "\n",
    "#BAO dz ratio\n",
    "def dz(zbao, om, obh2, w0, wa, n, h):\n",
    "    term1 = rs(zdrag(om, obh2, h), om, obh2, w0, wa, n, h)\n",
    "    term2 = Dv(zbao, om, w0, wa, n, h)\n",
    "    return term1/term2\n",
    "\n",
    "#DH distance\n",
    "def DH(z, om, w0, wa, n, h):\n",
    "    return c/H(1/(1+z), om, w0, wa, n, h)\n",
    "\n",
    "def hz(z, om, w0, wa, n, h):\n",
    "    return H(1/(1+z), om, w0, wa, n, h)\n",
    "\n",
    "#shift parameter\n",
    "def R(om, obh2, w0, wa, n, h):\n",
    "    term1 = np.sqrt(om*(100*h)**2) * DA(zcmb(om, obh2, h), om, w0, wa, n, h)\n",
    "    term2 = 1 + zcmb(om, obh2, h)\n",
    "    return term1*term2/c\n",
    "\n",
    "#angular scale of sound horizon at recombination\n",
    "def la(om, obh2, w0, wa, n, h):\n",
    "    term1 = np.pi * DA(zcmb(om, obh2, h), om, w0, wa, n, h)\n",
    "    term2 = 1 + zcmb(om, obh2, h)\n",
    "    term3 = rs(zcmb(om, obh2, h), om, obh2, w0, wa, n, h)\n",
    "    return term1*term2/term3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c65cc8b4-b7c2-4233-8786-a0a91b3b0411",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Data combination\n",
    "#CMB data\n",
    "datacmb = np.array([1.74963, 301.80845, 0.02237])\n",
    "\n",
    "# Covariance Matrix for {R, l_a, Omega_b h^2}\n",
    "covcmb = 10**-8 * np.array([\n",
    "    [1598.9554,  17112.007,  -36.311179],\n",
    "    [17112.007,  811208.45,  -494.79813],\n",
    "    [-36.311179, -494.79813,  2.1242182]\n",
    "])\n",
    "\n",
    "# Inverse Covariance Matrix\n",
    "invcov = np.linalg.inv(covcmb)\n",
    "\n",
    "# Number of CMB data points\n",
    "ndatCMB = len(datacmb)\n",
    "\n",
    "# Vector of residuals (Model - Data)\n",
    "def vec(om, obh2, w0, wa, n, h):\n",
    "    return np.array([\n",
    "        R(om, obh2, w0, wa, n, h) - datacmb[0],\n",
    "        la(om, obh2, w0, wa, n, h) - datacmb[1],\n",
    "        obh2 - datacmb[2]\n",
    "    ])\n",
    "\n",
    "# chi^2 function for CMB shift parameters\n",
    "def chi2R(om, obh2, w0, wa, n, h):\n",
    "    v = vec(om, obh2, w0, wa, n, h)\n",
    "    return v.T @ invcov @ v\n",
    "\n",
    "#ΒΑΟ data\n",
    "databao = np.array([\n",
    "    [0.106, 0.336, 0.015],\n",
    "    [0.44, 0.073, 0.031],\n",
    "    [0.6, 0.0726, 0.0164],\n",
    "    [0.73, 0.0592, 0.0185]\n",
    "])\n",
    "\n",
    "\n",
    "Cijbaoinv = np.array([\n",
    "    [1 / (databao[0, 2]**2), 0, 0, 0],\n",
    "    [0, 1040.3, -807.5, 336.8],\n",
    "    [0, -807.5, 3720.3, -1551.9],\n",
    "    [0, 336.8, -1551.9, 2914.9]\n",
    "])\n",
    "\n",
    "def vecbao(om, obh2, w0, wa, n, h):\n",
    "    return np.array([\n",
    "        databao[i, 1] - dz(databao[i, 0], om, obh2, w0, wa, n, h)\n",
    "        for i in range(len(databao))\n",
    "    ])\n",
    "\n",
    "#SDSS\n",
    "dataSDSS = np.array([\n",
    "    [0.15, 4.465666824, 0.1681350461],\n",
    "    [0.32, 8.4673, 0.167],\n",
    "    [0.57, 13.7728, 0.134],\n",
    "])\n",
    "\n",
    "#Ly-a\n",
    "dataLya = np.array([\n",
    "    [2.34, 11.2005988024, 0.56],\n",
    "    [2.34, 8.86, 0.29],\n",
    "])\n",
    "\n",
    "CijLyainv = np.array([\n",
    "    [1 / (dataLya[0, 2]**2), 0 ],\n",
    "    [0, 1 / (dataLya[1, 2]**2)]\n",
    "])\n",
    "\n",
    "def vecLya(om, obh2, w0, wa, n, h):\n",
    "    return np.array([\n",
    "        dataLya[0, 1] - DA(dataLya[0, 0], om, w0, wa, n, h) / rs(zdrag(om, obh2, h), om, obh2, w0, wa, n, h),\n",
    "        dataLya[1, 1] - DH(dataLya[1, 0], om, w0, wa, n, h) / rs(zdrag(om, obh2, h), om, obh2, w0, wa, n, h)\n",
    "    ])\n",
    "\n",
    "ndatBAO = len(dataSDSS) + len(databao) + len(dataLya)\n",
    "\n",
    "#total BAO χ^2\n",
    "def chi2bao(om, obh2, w0, wa, n, h):\n",
    "    v_bao = vecbao(om, obh2, w0, wa, n, h)\n",
    "    v_lya = vecLya(om, obh2, w0, wa, n, h)\n",
    "    \n",
    "    term_bao = v_bao.T @ Cijbaoinv @ v_bao\n",
    "    \n",
    "    term_sdss = sum(\n",
    "        ((dataSDSS[i, 1] - 1 / dz(dataSDSS[i, 0], om, obh2, w0, wa, n, h))**2) / (dataSDSS[i, 2]**2)\n",
    "        for i in range(len(dataSDSS))\n",
    "    )\n",
    "    \n",
    "    term_lya = v_lya.T @ CijLyainv @ v_lya\n",
    "    \n",
    "    return term_bao + term_sdss + term_lya"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e33f1775-fa4a-45a4-adf4-d9f0a58ca554",
   "metadata": {},
   "outputs": [],
   "source": [
    "#full pantheon data\n",
    "import pandas as pd\n",
    "dataPanth = pd.read_csv(\"lcparam_full_long_zhel.txt\", sep=r'\\s+')\n",
    "ndatPanth = len(dataPanth)\n",
    "Dij0 = np.loadtxt(\"sys_full_long.txt\")\n",
    "Dij = Dij0[1:].reshape((1048,1048))\n",
    "Id = np.ones(ndatPanth)\n",
    "#Construct of the covariance using only the statistical uncertainties\n",
    "Covstat = np.diag(dataPanth.iloc[:, 5]**2)\n",
    "#Construct of the covariance using only the systematic uncertainties\n",
    "Covsys = Dij\n",
    "Cijtot = Covsys + Covstat\n",
    "InvCovTotal = np.linalg.inv(Cijtot)\n",
    "InvCovstat = np.linalg.inv(Covstat)\n",
    "\n",
    "# χ^2 for the pantheon data\n",
    "def chi2Panth(M, om, w0, wa, n, h):\n",
    "    zcmb = dataPanth.iloc[:, 1].values  \n",
    "    zhel = dataPanth.iloc[:, 2].values  \n",
    "    mb = dataPanth.iloc[:, 4].values\n",
    "    dl_vals = DL(zcmb, om, w0, wa, n, h)\n",
    "    delta_m = mb - (M + 5 * np.log10(((1 + zhel) / (1 + zcmb)) * dl_vals) + 25)\n",
    "    chi2 = delta_m.T @ InvCovTotal @ delta_m\n",
    "    return chi2\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c22d0dde-e4c5-4b37-9d0f-42a8e075312f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test Mathematica point in Python: 1032.8152046282437\n",
      "1026.259730262297\n",
      "0.12038195926556611\n",
      "6.43509240668118\n"
     ]
    }
   ],
   "source": [
    "#total χ^2\n",
    "def chi2total(M, om, w0, wa, h):\n",
    "    return chi2R(om, obh2pl, w0, wa, 1.0, h) + chi2bao(om, obh2pl, w0, wa, 1.0, h) + chi2Panth(M, om, w0, wa, 1.0, h)\n",
    "ndatot = ndatCMB + ndatBAO + ndatPanth\n",
    "test_val = chi2total(-19.4294, 0.316911, -1, 0, 0.672379)\n",
    "print(\"Test Mathematica point in Python:\", test_val)\n",
    "x_val = chi2Panth(-19.4294, 0.316911, -1.0, 0.0, 1.0, 0.672379)\n",
    "print(x_val)\n",
    "y_val = chi2R(0.316911, obh2pl, -1.0, 0.0, 1, 0.672379)\n",
    "print(y_val)\n",
    "z_val = chi2bao(0.316911, obh2pl, -1.0, 0.0, 1, 0.672379)\n",
    "print(z_val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "be1dd99a-d26f-4b35-a899-cd1aaa361dca",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Minimum chi2: 1032.790623469141\n",
      "Mlcdm = -19.428583683065796\n",
      "omlcdm = 0.31658931689961795\n",
      "wolcdm = -1\n",
      "walcdm = 0\n",
      "hlcdm = 0.6725515843446923\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import basinhopping\n",
    "\n",
    "def obj_lcdm(params):\n",
    "    M, om, h = params\n",
    "    if om <= 0.01 or om >= 0.99 or h <= 0.1 or h >= 1.5:\n",
    "        return 1e10\n",
    "    try:\n",
    "        val = chi2total(M, om, -1, 0, h)\n",
    "        if np.isnan(val):\n",
    "            return 1e10\n",
    "        return val\n",
    "    except:\n",
    "        return 1e10\n",
    "\n",
    "initial_guess = [-19.05, 0.3205, 0.6705]\n",
    "\n",
    "minimizer_kwargs = {\n",
    "    \"method\": \"Powell\",\n",
    "    \"bounds\": [(None, None), (0.01, 0.99), (0.1, 1.5)],\n",
    "    \"options\": {'ftol': 1e-12, 'xtol': 1e-12}\n",
    "}\n",
    "\n",
    "minlcdm = basinhopping(\n",
    "    obj_lcdm,\n",
    "    x0=initial_guess,\n",
    "    minimizer_kwargs=minimizer_kwargs,\n",
    "    niter=10,\n",
    "    stepsize=0.05,\n",
    "    T=1.0\n",
    ")\n",
    "\n",
    "Mlcdm = minlcdm.x[0]\n",
    "omlcdm = minlcdm.x[1]\n",
    "wolcdm = -1\n",
    "walcdm = 0\n",
    "hlcdm = minlcdm.x[2]\n",
    "\n",
    "print(\"Minimum chi2:\", minlcdm.fun)\n",
    "print(\"Mlcdm =\", Mlcdm)\n",
    "print(\"omlcdm =\", omlcdm)\n",
    "print(\"wolcdm =\", wolcdm)\n",
    "print(\"walcdm =\", walcdm)\n",
    "print(\"hlcdm =\", hlcdm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "23f8caea-65e7-4971-ae7a-8e14928ab8ca",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Minimum chi2: 1032.6923784195608\n",
      "Mwwcdm = -19.42339752100402\n",
      "omwwcdm = 0.3144839102637591\n",
      "wowwcdm = -1.0103536711064154\n",
      "wawwcdm = 0\n",
      "hwwcdm = 0.6751045321388925\n"
     ]
    }
   ],
   "source": [
    "#χ^2 for the wCDM model using all three likelihoods\n",
    "def obj_wwcdm(params):\n",
    "    M, om, w0, h = params\n",
    "    if om <= 0.01 or om >= 0.99 or w0 <= -3.0 or w0 >= 1.0 or h <= 0.1 or h >= 1.5:\n",
    "        return 1e10\n",
    "    try:\n",
    "        val = chi2total(M, om, w0, 0, h)\n",
    "        if np.isnan(val):\n",
    "            return 1e10\n",
    "        return val\n",
    "    except:\n",
    "        return 1e10\n",
    "\n",
    "initial_guess_wwcdm = [-19.05, 0.3205, -1.0, 0.6705]\n",
    "\n",
    "minimizer_kwargs_wwcdm = {\n",
    "    \"method\": \"Powell\",\n",
    "    \"bounds\": [(None, None), (0.01, 0.99), (-3.0, 1.0), (0.1, 1.5)],\n",
    "    \"options\": {'ftol': 1e-12, 'xtol': 1e-12}\n",
    "}\n",
    "\n",
    "minwwcdm = basinhopping(\n",
    "    obj_wwcdm,\n",
    "    x0=initial_guess_wwcdm,\n",
    "    minimizer_kwargs=minimizer_kwargs_wwcdm,\n",
    "    niter=10,\n",
    "    stepsize=0.05,\n",
    "    T=1.0\n",
    ")\n",
    "\n",
    "Mwwcdm = minwwcdm.x[0]\n",
    "omwwcdm = minwwcdm.x[1]\n",
    "wowwcdm = minwwcdm.x[2]\n",
    "wawwcdm = 0\n",
    "hwwcdm = minwwcdm.x[3]\n",
    "\n",
    "print(\"Minimum chi2:\", minwwcdm.fun)\n",
    "print(\"Mwwcdm =\", Mwwcdm)\n",
    "print(\"omwwcdm =\", omwwcdm)\n",
    "print(\"wowwcdm =\", wowwcdm)\n",
    "print(\"wawwcdm =\", wawwcdm)\n",
    "print(\"hwwcdm =\", hwwcdm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "15eeb316-0da2-4c67-8279-b9f6954180fb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Minimum chi2: 1031.8200001076395\n",
      "Mcpl = -19.434704760745337\n",
      "omcpl = 0.31459788431390306\n",
      "wocpl = -1.0810461486577583\n",
      "wacpl = 0.2798724208402503\n",
      "hcpl = 0.6742663793150537\n"
     ]
    }
   ],
   "source": [
    "#χ^2 for the CPL model using all three likelihoods\n",
    "def obj_cpl(params):\n",
    "    M, om, w0, wa, h = params\n",
    "    if om <= 0.01 or om >= 0.99 or w0 <= -3.0 or w0 >= 1.0 or wa <= -3.0 or wa >= 3.0 or h <= 0.1 or h >= 1.5:\n",
    "        return 1e10\n",
    "    try:\n",
    "        val = chi2total(M, om, w0, wa, h)\n",
    "        if np.isnan(val):\n",
    "            return 1e10\n",
    "        return val\n",
    "    except:\n",
    "        return 1e10\n",
    "\n",
    "initial_guess_cpl = [-19.05, 0.3205, -1.0, 0.0, 0.6705]\n",
    "\n",
    "minimizer_kwargs_cpl = {\n",
    "    \"method\": \"Powell\",\n",
    "    \"bounds\": [(None, None), (0.01, 0.99), (-3.0, 1.0), (-3.0, 3.0), (0.1, 1.5)],\n",
    "    \"options\": {'ftol': 1e-12, 'xtol': 1e-12}\n",
    "}\n",
    "\n",
    "mincpl = basinhopping(\n",
    "    obj_cpl,\n",
    "    x0=initial_guess_cpl,\n",
    "    minimizer_kwargs=minimizer_kwargs_cpl,\n",
    "    niter=10,\n",
    "    stepsize=0.05,\n",
    "    T=1.0\n",
    ")\n",
    "\n",
    "Mcpl = mincpl.x[0]\n",
    "omcpl = mincpl.x[1]\n",
    "wocpl = mincpl.x[2]\n",
    "wacpl = mincpl.x[3]\n",
    "hcpl = mincpl.x[4]\n",
    "\n",
    "print(\"Minimum chi2:\", mincpl.fun)\n",
    "print(\"Mcpl =\", Mcpl)\n",
    "print(\"omcpl =\", omcpl)\n",
    "print(\"wocpl =\", wocpl)\n",
    "print(\"wacpl =\", wacpl)\n",
    "print(\"hcpl =\", hcpl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5e6db67a-3807-410e-ab72-a06bcb90d38b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Riess point\n",
    "def chi2riessh(h):\n",
    "    return (h - h0riess)**2/(0.0142)**2\n",
    "\n",
    "def chi2total5(M, om, w0, wa, h):\n",
    "    return chi2R(om, obh2pl, w0, wa, 1.0, h) + chi2bao(om, obh2pl, w0, wa, 1.0, h) + chi2Panth(M, om, w0, wa, 1.0, h) + chi2riessh(h)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "d3fe130b-4e3f-4289-ae75-6ef2bcf5dc6d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Minimum chi2: 1047.3435013445721\n",
      "Mccpl = -19.393001086276534\n",
      "omccpl = 0.29825787202393267\n",
      "woccpl = -1.1250296972315699\n",
      "waccpl = 0.2067219558279372\n",
      "hccpl = 0.6931120021883351\n"
     ]
    }
   ],
   "source": [
    "#χ^2 for the cpl model using all three likelihoods and Riess point\n",
    "def obj_cplriess(params):\n",
    "    M, om, w0, wa, h = params\n",
    "    if om <= 0.01 or om >= 0.99 or w0 <= -3.0 or w0 >= 1.0 or wa <= -3.0 or wa >= 3.0 or h <= 0.1 or h >= 1.5:\n",
    "        return 1e10\n",
    "    try:\n",
    "        val = chi2total5(M, om, w0, wa, h)\n",
    "        if np.isnan(val):\n",
    "            return 1e10\n",
    "        return val\n",
    "    except:\n",
    "        return 1e10\n",
    "\n",
    "initial_guess_cplriess = [-19.05, 0.3205, -1.0, 0.0, 0.6705]\n",
    "\n",
    "minimizer_kwargs_cplriess = {\n",
    "    \"method\": \"Powell\",\n",
    "    \"bounds\": [(None, None), (0.01, 0.99), (-3.0, 1.0), (-3.0, 3.0), (0.1, 1.5)],\n",
    "    \"options\": {'ftol': 1e-12, 'xtol': 1e-12}\n",
    "}\n",
    "\n",
    "mincplriess = basinhopping(\n",
    "    obj_cplriess,\n",
    "    x0=initial_guess_cplriess,\n",
    "    minimizer_kwargs=minimizer_kwargs_cplriess,\n",
    "    niter=10,\n",
    "    stepsize=0.05,\n",
    "    T=1.0\n",
    ")\n",
    "\n",
    "Mccpl = mincplriess.x[0]\n",
    "omccpl = mincplriess.x[1]\n",
    "woccpl = mincplriess.x[2]\n",
    "waccpl = mincplriess.x[3]\n",
    "hccpl = mincplriess.x[4]\n",
    "\n",
    "print(\"Minimum chi2:\", mincplriess.fun)\n",
    "print(\"Mccpl =\", Mccpl)\n",
    "print(\"omccpl =\", omccpl)\n",
    "print(\"woccpl =\", woccpl)\n",
    "print(\"waccpl =\", waccpl)\n",
    "print(\"hccpl =\", hccpl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d83b17bc-fe2f-4034-b26a-253c19ee429f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#χ^2 for the wcdm model using all three likelihoods and Riess point\n",
    "def obj_wcdmriess(params):\n",
    "    M, om, w0, h = params\n",
    "    if om <= 0.01 or om >= 0.99 or w0 <= -3.0 or w0 >= 1.0 or h <= 0.1 or h >= 1.5:\n",
    "        return 1e10\n",
    "    try:\n",
    "        val = chi2total5(M, om, w0, 0.0, h)\n",
    "        if np.isnan(val):\n",
    "            return 1e10\n",
    "        return val\n",
    "    except:\n",
    "        return 1e10\n",
    "\n",
    "initial_guess_wcdmriess = [-19.05, 0.3205, -1.0, 0.6705]\n",
    "\n",
    "minimizer_kwargs_wcdmriess = {\n",
    "    \"method\": \"Powell\",\n",
    "    \"bounds\": [(None, None), (0.01, 0.99), (-3.0, 1.0), (0.1, 1.5)],\n",
    "    \"options\": {'ftol': 1e-12, 'xtol': 1e-12}\n",
    "}\n",
    "\n",
    "minwcdmriess = basinhopping(\n",
    "    obj_wcdmriess,\n",
    "    x0=initial_guess_wcdmriess,\n",
    "    minimizer_kwargs=minimizer_kwargs_wcdmriess,\n",
    "    niter=10,\n",
    "    stepsize=0.05,\n",
    "    T=1.0\n",
    ")\n",
    "\n",
    "Mwcdmriess = minwcdmriess.x[0]\n",
    "omwcdmriess = minwcdmriess.x[1]\n",
    "wowcdmriess = minwcdmriess.x[2]\n",
    "wawcdmriess = 0.0\n",
    "hwcdmriess = minwcdmriess.x[3]\n",
    "\n",
    "print(\"Minimum chi2:\", minwcdmriess.fun)\n",
    "print(\"Mwcdmriess =\", Mwcdmriess)\n",
    "print(\"omwcdmriess =\", omwcdmriess)\n",
    "print(\"wowcdmriess =\", wowcdmriess)\n",
    "print(\"wawcdmriess =\", wawcdmriess)\n",
    "print(\"hwcdmriess =\", hwcdmriess)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00139830-5e95-4d6b-9700-0099527357be",
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}