{ "cells": [ { "cell_type": "code", "execution_count": 34, "id": "0a6dc9b2-d9d0-4fc0-a84e-87ea68ccd9f5", "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": 35, "id": "5b591783-a614-4bf3-a99a-56914c10d387", "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": 36, "id": "a99a7de3-d5b8-4e30-a760-42b6da099748", "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(om, obh2, h))\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", " \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, 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": 37, "id": "9a5b20d6-d433-4234-9967-f955f17cc8b5", "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 / (databao[0, 2]**2), 0 ],\n", " [0, 1 / (databao[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, w0, wa, n, h),\n", " dataLya[1, 1] - DH(dataLya[1, 0], om, w0, wa, n, h) / rs(zdrag, 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": 38, "id": "ed475544-fffb-4877-80d6-3e10c2575be0", "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": 39, "id": "dcf5402a-9061-45be-9d10-79675dd18b91", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1060\n" ] } ], "source": [ "#total χ^2\n", "def chi2total(M, om, w0, wa, h):\n", " return chi2R(om, obh2pl, w0, wa, 1, h) + chi2bao(om, obh2pl, w0, wa, 1, h) + chi2Panth(M, om, w0, wa, 1, h)\n", "ndatot = ndatCMB + ndatBAO + ndatPanth\n" ] }, { "cell_type": "code", "execution_count": 43, "id": "dfc8ea13-e35a-426e-892a-212a15e2face", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Υπολογισμός με Basin-Hopping (Γρήγορο Monte Carlo)...\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "C:\\Users\\User\\anaconda3\\Lib\\site-packages\\scipy\\optimize\\_numdiff.py:590: RuntimeWarning: invalid value encountered in subtract\n", " df = fun(x) - f0\n" ] }, { "ename": "KeyboardInterrupt", "evalue": "", "output_type": "error", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)", "Cell \u001b[1;32mIn[43], line 24\u001b[0m\n\u001b[0;32m 19\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mΥπολογισμός με Basin-Hopping (Γρήγορο Monte Carlo)...\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m 21\u001b[0m \u001b[38;5;66;03m# 3. Εκτέλεση του Basin-Hopping\u001b[39;00m\n\u001b[0;32m 22\u001b[0m \u001b[38;5;66;03m# niter: Πόσα τυχαία άλματα (Monte Carlo steps) θα κάνει. \u001b[39;00m\n\u001b[0;32m 23\u001b[0m \u001b[38;5;66;03m# Βάλε ένα μικρό νούμερο (π.χ. 10 ή 20) για να είναι πολύ γρήγορο.\u001b[39;00m\n\u001b[1;32m---> 24\u001b[0m minlcdm_mc \u001b[38;5;241m=\u001b[39m basinhopping(\n\u001b[0;32m 25\u001b[0m objective_function, \n\u001b[0;32m 26\u001b[0m initial_guess, \n\u001b[0;32m 27\u001b[0m minimizer_kwargs\u001b[38;5;241m=\u001b[39mminimizer_kwargs, \n\u001b[0;32m 28\u001b[0m niter\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m15\u001b[39m, \u001b[38;5;66;03m# Αριθμός αλμάτων (όσο μικρότερο, τόσο πιο γρήγορο)\u001b[39;00m\n\u001b[0;32m 29\u001b[0m stepsize\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0.5\u001b[39m, \u001b[38;5;66;03m# Πόσο μακριά να πηδάει σε κάθε βήμα\u001b[39;00m\n\u001b[0;32m 30\u001b[0m seed\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m42\u001b[39m \u001b[38;5;66;03m# Για να παίρνεις τα ίδια αποτελέσματα αν το ξανατρέξεις\u001b[39;00m\n\u001b[0;32m 31\u001b[0m )\n\u001b[0;32m 33\u001b[0m \u001b[38;5;66;03m# 4. Αποτελέσματα\u001b[39;00m\n\u001b[0;32m 34\u001b[0m M_min, om_min, h_min \u001b[38;5;241m=\u001b[39m minlcdm_mc\u001b[38;5;241m.\u001b[39mx\n", "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\scipy\\optimize\\_basinhopping.py:714\u001b[0m, in \u001b[0;36mbasinhopping\u001b[1;34m(func, x0, niter, T, stepsize, minimizer_kwargs, take_step, accept_test, callback, interval, disp, niter_success, seed, target_accept_rate, stepwise_factor)\u001b[0m\n\u001b[0;32m 711\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m niter_success \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[0;32m 712\u001b[0m niter_success \u001b[38;5;241m=\u001b[39m niter \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m2\u001b[39m\n\u001b[1;32m--> 714\u001b[0m bh \u001b[38;5;241m=\u001b[39m BasinHoppingRunner(x0, wrapped_minimizer, take_step_wrapped,\n\u001b[0;32m 715\u001b[0m accept_tests, disp\u001b[38;5;241m=\u001b[39mdisp)\n\u001b[0;32m 717\u001b[0m \u001b[38;5;66;03m# The wrapped minimizer is called once during construction of\u001b[39;00m\n\u001b[0;32m 718\u001b[0m \u001b[38;5;66;03m# BasinHoppingRunner, so run the callback\u001b[39;00m\n\u001b[0;32m 719\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mcallable\u001b[39m(callback):\n", "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\scipy\\optimize\\_basinhopping.py:78\u001b[0m, in \u001b[0;36mBasinHoppingRunner.__init__\u001b[1;34m(self, x0, minimizer, step_taking, accept_tests, disp)\u001b[0m\n\u001b[0;32m 75\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mres\u001b[38;5;241m.\u001b[39mminimization_failures \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m0\u001b[39m\n\u001b[0;32m 77\u001b[0m \u001b[38;5;66;03m# do initial minimization\u001b[39;00m\n\u001b[1;32m---> 78\u001b[0m minres \u001b[38;5;241m=\u001b[39m minimizer(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mx)\n\u001b[0;32m 79\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m minres\u001b[38;5;241m.\u001b[39msuccess:\n\u001b[0;32m 80\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mres\u001b[38;5;241m.\u001b[39mminimization_failures \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1\u001b[39m\n", "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\scipy\\optimize\\_basinhopping.py:303\u001b[0m, in \u001b[0;36mMinimizerWrapper.__call__\u001b[1;34m(self, x0)\u001b[0m\n\u001b[0;32m 301\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mminimizer(x0, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mkwargs)\n\u001b[0;32m 302\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m--> 303\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mminimizer(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfunc, x0, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mkwargs)\n", "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\scipy\\optimize\\_minimize.py:713\u001b[0m, in \u001b[0;36mminimize\u001b[1;34m(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)\u001b[0m\n\u001b[0;32m 710\u001b[0m res \u001b[38;5;241m=\u001b[39m _minimize_newtoncg(fun, x0, args, jac, hess, hessp, callback,\n\u001b[0;32m 711\u001b[0m \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39moptions)\n\u001b[0;32m 712\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m meth \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124ml-bfgs-b\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[1;32m--> 713\u001b[0m res \u001b[38;5;241m=\u001b[39m _minimize_lbfgsb(fun, x0, args, jac, bounds,\n\u001b[0;32m 714\u001b[0m callback\u001b[38;5;241m=\u001b[39mcallback, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39moptions)\n\u001b[0;32m 715\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m meth \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mtnc\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[0;32m 716\u001b[0m res \u001b[38;5;241m=\u001b[39m _minimize_tnc(fun, x0, args, jac, bounds, callback\u001b[38;5;241m=\u001b[39mcallback,\n\u001b[0;32m 717\u001b[0m \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39moptions)\n", "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\scipy\\optimize\\_lbfgsb_py.py:407\u001b[0m, in \u001b[0;36m_minimize_lbfgsb\u001b[1;34m(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, finite_diff_rel_step, **unknown_options)\u001b[0m\n\u001b[0;32m 401\u001b[0m task_str \u001b[38;5;241m=\u001b[39m task\u001b[38;5;241m.\u001b[39mtobytes()\n\u001b[0;32m 402\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m task_str\u001b[38;5;241m.\u001b[39mstartswith(\u001b[38;5;124mb\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mFG\u001b[39m\u001b[38;5;124m'\u001b[39m):\n\u001b[0;32m 403\u001b[0m \u001b[38;5;66;03m# The minimization routine wants f and g at the current x.\u001b[39;00m\n\u001b[0;32m 404\u001b[0m \u001b[38;5;66;03m# Note that interruptions due to maxfun are postponed\u001b[39;00m\n\u001b[0;32m 405\u001b[0m \u001b[38;5;66;03m# until the completion of the current minimization iteration.\u001b[39;00m\n\u001b[0;32m 406\u001b[0m \u001b[38;5;66;03m# Overwrite f and g:\u001b[39;00m\n\u001b[1;32m--> 407\u001b[0m f, g \u001b[38;5;241m=\u001b[39m func_and_grad(x)\n\u001b[0;32m 408\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m task_str\u001b[38;5;241m.\u001b[39mstartswith(\u001b[38;5;124mb\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mNEW_X\u001b[39m\u001b[38;5;124m'\u001b[39m):\n\u001b[0;32m 409\u001b[0m \u001b[38;5;66;03m# new iteration\u001b[39;00m\n\u001b[0;32m 410\u001b[0m n_iterations \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1\u001b[39m\n", "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\scipy\\optimize\\_differentiable_functions.py:297\u001b[0m, in \u001b[0;36mScalarFunction.fun_and_grad\u001b[1;34m(self, x)\u001b[0m\n\u001b[0;32m 295\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_update_x_impl(x)\n\u001b[0;32m 296\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_update_fun()\n\u001b[1;32m--> 297\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_update_grad()\n\u001b[0;32m 298\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mf, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mg\n", "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\scipy\\optimize\\_differentiable_functions.py:267\u001b[0m, in \u001b[0;36mScalarFunction._update_grad\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 265\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_update_grad\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[0;32m 266\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mg_updated:\n\u001b[1;32m--> 267\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_update_grad_impl()\n\u001b[0;32m 268\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mg_updated \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mTrue\u001b[39;00m\n", "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\scipy\\optimize\\_differentiable_functions.py:181\u001b[0m, in \u001b[0;36mScalarFunction.__init__..update_grad\u001b[1;34m()\u001b[0m\n\u001b[0;32m 179\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_update_fun()\n\u001b[0;32m 180\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mngev \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1\u001b[39m\n\u001b[1;32m--> 181\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mg \u001b[38;5;241m=\u001b[39m approx_derivative(fun_wrapped, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mx, f0\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mf,\n\u001b[0;32m 182\u001b[0m \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mfinite_diff_options)\n", "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\scipy\\optimize\\_numdiff.py:519\u001b[0m, in \u001b[0;36mapprox_derivative\u001b[1;34m(fun, x0, method, rel_step, abs_step, f0, bounds, sparsity, as_linear_operator, args, kwargs)\u001b[0m\n\u001b[0;32m 516\u001b[0m use_one_sided \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mFalse\u001b[39;00m\n\u001b[0;32m 518\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m sparsity \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m--> 519\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m _dense_difference(fun_wrapped, x0, f0, h,\n\u001b[0;32m 520\u001b[0m use_one_sided, method)\n\u001b[0;32m 521\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m 522\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m issparse(sparsity) \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(sparsity) \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m2\u001b[39m:\n", "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\scipy\\optimize\\_numdiff.py:590\u001b[0m, in \u001b[0;36m_dense_difference\u001b[1;34m(fun, x0, f0, h, use_one_sided, method)\u001b[0m\n\u001b[0;32m 588\u001b[0m x \u001b[38;5;241m=\u001b[39m x0 \u001b[38;5;241m+\u001b[39m h_vecs[i]\n\u001b[0;32m 589\u001b[0m dx \u001b[38;5;241m=\u001b[39m x[i] \u001b[38;5;241m-\u001b[39m x0[i] \u001b[38;5;66;03m# Recompute dx as exactly representable number.\u001b[39;00m\n\u001b[1;32m--> 590\u001b[0m df \u001b[38;5;241m=\u001b[39m fun(x) \u001b[38;5;241m-\u001b[39m f0\n\u001b[0;32m 591\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m method \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124m3-point\u001b[39m\u001b[38;5;124m'\u001b[39m \u001b[38;5;129;01mand\u001b[39;00m use_one_sided[i]:\n\u001b[0;32m 592\u001b[0m x1 \u001b[38;5;241m=\u001b[39m x0 \u001b[38;5;241m+\u001b[39m h_vecs[i]\n", "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\scipy\\optimize\\_numdiff.py:470\u001b[0m, in \u001b[0;36mapprox_derivative..fun_wrapped\u001b[1;34m(x)\u001b[0m\n\u001b[0;32m 467\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m xp\u001b[38;5;241m.\u001b[39misdtype(x\u001b[38;5;241m.\u001b[39mdtype, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mreal floating\u001b[39m\u001b[38;5;124m\"\u001b[39m):\n\u001b[0;32m 468\u001b[0m x \u001b[38;5;241m=\u001b[39m xp\u001b[38;5;241m.\u001b[39mastype(x, x0\u001b[38;5;241m.\u001b[39mdtype)\n\u001b[1;32m--> 470\u001b[0m f \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39matleast_1d(fun(x, \u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs))\n\u001b[0;32m 471\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m f\u001b[38;5;241m.\u001b[39mndim \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m1\u001b[39m:\n\u001b[0;32m 472\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mRuntimeError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m`fun` return value has \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m 473\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mmore than 1 dimension.\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n", "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\scipy\\optimize\\_differentiable_functions.py:145\u001b[0m, in \u001b[0;36mScalarFunction.__init__..fun_wrapped\u001b[1;34m(x)\u001b[0m\n\u001b[0;32m 141\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mnfev \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1\u001b[39m\n\u001b[0;32m 142\u001b[0m \u001b[38;5;66;03m# Send a copy because the user may overwrite it.\u001b[39;00m\n\u001b[0;32m 143\u001b[0m \u001b[38;5;66;03m# Overwriting results in undefined behaviour because\u001b[39;00m\n\u001b[0;32m 144\u001b[0m \u001b[38;5;66;03m# fun(self.x) will change self.x, with the two no longer linked.\u001b[39;00m\n\u001b[1;32m--> 145\u001b[0m fx \u001b[38;5;241m=\u001b[39m fun(np\u001b[38;5;241m.\u001b[39mcopy(x), \u001b[38;5;241m*\u001b[39margs)\n\u001b[0;32m 146\u001b[0m \u001b[38;5;66;03m# Make sure the function returns a true scalar\u001b[39;00m\n\u001b[0;32m 147\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m np\u001b[38;5;241m.\u001b[39misscalar(fx):\n", "Cell \u001b[1;32mIn[42], line 6\u001b[0m, in \u001b[0;36mobjective_function\u001b[1;34m(params)\u001b[0m\n\u001b[0;32m 4\u001b[0m M, om, h \u001b[38;5;241m=\u001b[39m params\n\u001b[0;32m 5\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[1;32m----> 6\u001b[0m val \u001b[38;5;241m=\u001b[39m chi2total(M, om, \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m, \u001b[38;5;241m0\u001b[39m, h)\n\u001b[0;32m 8\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m val \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m \u001b[38;5;129;01mor\u001b[39;00m np\u001b[38;5;241m.\u001b[39misnan(val):\n\u001b[0;32m 9\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m np\u001b[38;5;241m.\u001b[39minf\n", "Cell \u001b[1;32mIn[39], line 3\u001b[0m, in \u001b[0;36mchi2total\u001b[1;34m(M, om, w0, wa, h)\u001b[0m\n\u001b[0;32m 2\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mchi2total\u001b[39m(M, om, w0, wa, h):\n\u001b[1;32m----> 3\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m chi2R(om, obh2pl, w0, wa, \u001b[38;5;241m1\u001b[39m, h) \u001b[38;5;241m+\u001b[39m chi2bao(om, obh2pl, w0, wa, \u001b[38;5;241m1\u001b[39m, h) \u001b[38;5;241m+\u001b[39m chi2Panth(M, om, w0, wa, \u001b[38;5;241m1\u001b[39m, h)\n", "Cell \u001b[1;32mIn[37], line 28\u001b[0m, in \u001b[0;36mchi2R\u001b[1;34m(om, obh2, w0, wa, n, h)\u001b[0m\n\u001b[0;32m 27\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mchi2R\u001b[39m(om, obh2, w0, wa, n, h):\n\u001b[1;32m---> 28\u001b[0m v \u001b[38;5;241m=\u001b[39m vec(om, obh2, w0, wa, n, h)\n\u001b[0;32m 29\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m v\u001b[38;5;241m.\u001b[39mT \u001b[38;5;241m@\u001b[39m invcov \u001b[38;5;241m@\u001b[39m v\n", "Cell \u001b[1;32mIn[37], line 21\u001b[0m, in \u001b[0;36mvec\u001b[1;34m(om, obh2, w0, wa, n, h)\u001b[0m\n\u001b[0;32m 19\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mvec\u001b[39m(om, obh2, w0, wa, n, h):\n\u001b[0;32m 20\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m np\u001b[38;5;241m.\u001b[39marray([\n\u001b[1;32m---> 21\u001b[0m R(om, obh2, w0, wa, n, h) \u001b[38;5;241m-\u001b[39m datacmb[\u001b[38;5;241m0\u001b[39m],\n\u001b[0;32m 22\u001b[0m la(om, obh2, w0, wa, n, h) \u001b[38;5;241m-\u001b[39m datacmb[\u001b[38;5;241m1\u001b[39m],\n\u001b[0;32m 23\u001b[0m obh2 \u001b[38;5;241m-\u001b[39m datacmb[\u001b[38;5;241m2\u001b[39m]\n\u001b[0;32m 24\u001b[0m ])\n", "Cell \u001b[1;32mIn[36], line 122\u001b[0m, in \u001b[0;36mR\u001b[1;34m(om, obh2, w0, wa, n, h)\u001b[0m\n\u001b[0;32m 121\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mR\u001b[39m(om, obh2, w0, wa, n, h):\n\u001b[1;32m--> 122\u001b[0m term1 \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39msqrt(om\u001b[38;5;241m*\u001b[39m(\u001b[38;5;241m100\u001b[39m\u001b[38;5;241m*\u001b[39mh)\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m) \u001b[38;5;241m*\u001b[39m DA(zcmb(om, obh2, h), om, w0, wa, n, h)\n\u001b[0;32m 123\u001b[0m term2 \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1\u001b[39m \u001b[38;5;241m+\u001b[39m zcmb(om, obh2, h)\n\u001b[0;32m 124\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m term1\u001b[38;5;241m*\u001b[39mterm2\u001b[38;5;241m/\u001b[39mc\n", "Cell \u001b[1;32mIn[36], line 98\u001b[0m, in \u001b[0;36mDA\u001b[1;34m(z, om, w0, wa, n, h)\u001b[0m\n\u001b[0;32m 97\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mDA\u001b[39m(z, om, w0, wa, n, h):\n\u001b[1;32m---> 98\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m DL(z, om, w0, wa, n, h) \u001b[38;5;241m/\u001b[39m (\u001b[38;5;241m1.0\u001b[39m \u001b[38;5;241m+\u001b[39m z)\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m\n", "Cell \u001b[1;32mIn[36], line 92\u001b[0m, in \u001b[0;36mDL\u001b[1;34m(z, om, w0, wa, n, h)\u001b[0m\n\u001b[0;32m 91\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mDL\u001b[39m(z, om, w0, wa, n, h):\n\u001b[1;32m---> 92\u001b[0m sol_func \u001b[38;5;241m=\u001b[39m DLsol(om, w0, wa, n, h)\n\u001b[0;32m 93\u001b[0m y_val \u001b[38;5;241m=\u001b[39m sol_func(z)[\u001b[38;5;241m0\u001b[39m]\n\u001b[0;32m 94\u001b[0m dL_val \u001b[38;5;241m=\u001b[39m y_val \u001b[38;5;241m*\u001b[39m (\u001b[38;5;241m1.0\u001b[39m \u001b[38;5;241m+\u001b[39m z)\n", "Cell \u001b[1;32mIn[36], line 88\u001b[0m, in \u001b[0;36mDLsol\u001b[1;34m(om, w0, wa, n, h)\u001b[0m\n\u001b[0;32m 85\u001b[0m a \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1.0\u001b[39m \u001b[38;5;241m/\u001b[39m (\u001b[38;5;241m1.0\u001b[39m \u001b[38;5;241m+\u001b[39m zz)\n\u001b[0;32m 86\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;241m1.0\u001b[39m \u001b[38;5;241m/\u001b[39m (H(a, om, w0, wa, n, h) \u001b[38;5;241m/\u001b[39m (\u001b[38;5;241m100.0\u001b[39m \u001b[38;5;241m*\u001b[39m h))\n\u001b[1;32m---> 88\u001b[0m sol \u001b[38;5;241m=\u001b[39m solve_ivp(dy_dz, [\u001b[38;5;241m0\u001b[39m, \u001b[38;5;241m1300\u001b[39m], [\u001b[38;5;241m0.0\u001b[39m], dense_output\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m)\n\u001b[0;32m 89\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m sol\u001b[38;5;241m.\u001b[39msol\n", "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\scipy\\integrate\\_ivp\\ivp.py:655\u001b[0m, in \u001b[0;36msolve_ivp\u001b[1;34m(fun, t_span, y0, method, t_eval, dense_output, events, vectorized, args, **options)\u001b[0m\n\u001b[0;32m 653\u001b[0m status \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m\n\u001b[0;32m 654\u001b[0m \u001b[38;5;28;01mwhile\u001b[39;00m status \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m--> 655\u001b[0m message \u001b[38;5;241m=\u001b[39m solver\u001b[38;5;241m.\u001b[39mstep()\n\u001b[0;32m 657\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m solver\u001b[38;5;241m.\u001b[39mstatus \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mfinished\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[0;32m 658\u001b[0m status \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m0\u001b[39m\n", "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\scipy\\integrate\\_ivp\\base.py:197\u001b[0m, in \u001b[0;36mOdeSolver.step\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 195\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m 196\u001b[0m t \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mt\n\u001b[1;32m--> 197\u001b[0m success, message \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_step_impl()\n\u001b[0;32m 199\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m success:\n\u001b[0;32m 200\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mstatus \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mfailed\u001b[39m\u001b[38;5;124m'\u001b[39m\n", "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\scipy\\integrate\\_ivp\\rk.py:144\u001b[0m, in \u001b[0;36mRungeKutta._step_impl\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 141\u001b[0m h \u001b[38;5;241m=\u001b[39m t_new \u001b[38;5;241m-\u001b[39m t\n\u001b[0;32m 142\u001b[0m h_abs \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mabs(h)\n\u001b[1;32m--> 144\u001b[0m y_new, f_new \u001b[38;5;241m=\u001b[39m rk_step(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfun, t, y, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mf, h, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mA,\n\u001b[0;32m 145\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mB, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mC, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mK)\n\u001b[0;32m 146\u001b[0m scale \u001b[38;5;241m=\u001b[39m atol \u001b[38;5;241m+\u001b[39m np\u001b[38;5;241m.\u001b[39mmaximum(np\u001b[38;5;241m.\u001b[39mabs(y), np\u001b[38;5;241m.\u001b[39mabs(y_new)) \u001b[38;5;241m*\u001b[39m rtol\n\u001b[0;32m 147\u001b[0m error_norm \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_estimate_error_norm(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mK, h, scale)\n", "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\scipy\\integrate\\_ivp\\rk.py:63\u001b[0m, in \u001b[0;36mrk_step\u001b[1;34m(fun, t, y, f, h, A, B, C, K)\u001b[0m\n\u001b[0;32m 61\u001b[0m K[\u001b[38;5;241m0\u001b[39m] \u001b[38;5;241m=\u001b[39m f\n\u001b[0;32m 62\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m s, (a, c) \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28menumerate\u001b[39m(\u001b[38;5;28mzip\u001b[39m(A[\u001b[38;5;241m1\u001b[39m:], C[\u001b[38;5;241m1\u001b[39m:]), start\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1\u001b[39m):\n\u001b[1;32m---> 63\u001b[0m dy \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mdot(K[:s]\u001b[38;5;241m.\u001b[39mT, a[:s]) \u001b[38;5;241m*\u001b[39m h\n\u001b[0;32m 64\u001b[0m K[s] \u001b[38;5;241m=\u001b[39m fun(t \u001b[38;5;241m+\u001b[39m c \u001b[38;5;241m*\u001b[39m h, y \u001b[38;5;241m+\u001b[39m dy)\n\u001b[0;32m 66\u001b[0m y_new \u001b[38;5;241m=\u001b[39m y \u001b[38;5;241m+\u001b[39m h \u001b[38;5;241m*\u001b[39m np\u001b[38;5;241m.\u001b[39mdot(K[:\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m]\u001b[38;5;241m.\u001b[39mT, B)\n", "\u001b[1;31mKeyboardInterrupt\u001b[0m: " ] } ], "source": [ "#χ^2 for the ΛCDM model using all three likelihoods\n", "\n", "from scipy.optimize import basinhopping\n", "\n", "# 1. Αρχικές τιμές και όρια\n", "initial_guess = [-19.0, 0.32, 0.67]\n", "bounds = [\n", " (-20.0, -18.0), # M\n", " (0.1, 0.5), # om\n", " (0.5, 0.85) # h\n", "]\n", "\n", "# 2. Ρυθμίσεις για τη γρήγορη \"τοπική\" μέθοδο που θα χρησιμοποιεί μετά από κάθε άλμα\n", "minimizer_kwargs = {\n", " \"method\": \"L-BFGS-B\",\n", " \"bounds\": bounds\n", "}\n", "\n", "print(\"Υπολογισμός με Basin-Hopping (Γρήγορο Monte Carlo)...\")\n", "\n", "# 3. Εκτέλεση του Basin-Hopping\n", "# niter: Πόσα τυχαία άλματα (Monte Carlo steps) θα κάνει. \n", "# Βάλε ένα μικρό νούμερο (π.χ. 10 ή 20) για να είναι πολύ γρήγορο.\n", "minlcdm_mc = basinhopping(\n", " objective_function, \n", " initial_guess, \n", " minimizer_kwargs=minimizer_kwargs, \n", " niter=15, # Αριθμός αλμάτων (όσο μικρότερο, τόσο πιο γρήγορο)\n", " stepsize=0.5, # Πόσο μακριά να πηδάει σε κάθε βήμα\n", " seed=42 # Για να παίρνεις τα ίδια αποτελέσματα αν το ξανατρέξεις\n", ")\n", "\n", "# 4. Αποτελέσματα\n", "M_min, om_min, h_min = minlcdm_mc.x\n", "chi2_min = minlcdm_mc.fun\n", "\n", "print(f\"Το ελάχιστο chi2 είναι: {chi2_min:.4f}\")\n", "print(f\"Βέλτιστες παράμετροι:\")\n", "print(f\"M = {M_min:.4f}\")\n", "print(f\"om = {om_min:.4f}\")\n", "print(f\"h = {h_min:.4f}\")" ] }, { "cell_type": "code", "execution_count": 44, "id": "e0be455e-da47-4f1d-a75e-7c18d628500f", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Υπολογισμός τοπικού ελαχίστου (γρήγορη μέθοδος)...\n", "Το ελάχιστο chi2 είναι: inf\n", "Βέλτιστες παράμετροι:\n", "M = -19.0000\n", "om = 0.3200\n", "h = 0.6700\n" ] } ], "source": [ "from scipy.optimize import minimize\n", "\n", "# 1. Οι αρχικές σου τιμές (όπως τις είχες στο Mathematica)\n", "initial_guess = [-19.0, 0.32, 0.67]\n", "\n", "# 2. Τα όρια, για να μην ξεφύγει ο αλγόριθμος σε περιοχές που \"σκάνε\" τα ολοκληρώματα\n", "bounds = [\n", " (-20.0, -18.0), # M\n", " (0.1, 0.5), # om\n", " (0.5, 0.85) # h\n", "]\n", "\n", "print(\"Υπολογισμός τοπικού ελαχίστου (γρήγορη μέθοδος)...\")\n", "\n", "# 3. Χρησιμοποιούμε τη minimize με τη μέθοδο L-BFGS-B (πολύ γρήγορη, υποστηρίζει όρια)\n", "minlcdm_fast = minimize(\n", " objective_function, \n", " initial_guess, \n", " method='L-BFGS-B', \n", " bounds=bounds\n", ")\n", "\n", "# 4. Αποτελέσματα\n", "M_min, om_min, h_min = minlcdm_fast.x\n", "chi2_min = minlcdm_fast.fun\n", "\n", "print(f\"Το ελάχιστο chi2 είναι: {chi2_min:.4f}\")\n", "print(f\"Βέλτιστες παράμετροι:\")\n", "print(f\"M = {M_min:.4f}\")\n", "print(f\"om = {om_min:.4f}\")\n", "print(f\"h = {h_min:.4f}\")" ] }, { "cell_type": "code", "execution_count": null, "id": "df9d46d6-3d34-43e8-afa6-7be474a9c11b", "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 }