{ "cells": [ { "cell_type": "markdown", "metadata": { "id": "yxM8rgHvoHnF" }, "source": [ "\n", "# Fitting Binary Lenses\n", "\n", "
\n", "\n", "\"Open" ] }, { "cell_type": "markdown", "metadata": { "id": "r5jtB9cFoHnG" }, "source": [ "## Learning Goals\n", "By the end of this notebook, you will be able to:\n", "- Describe how binary-lens parameters ($s$, $q$, $\\alpha$, $\\rho$) control caustic topology and magnification structure.\n", "- Construct binary-lens magnification curves using `VBMicrolensing` and inspect caustic transitions interactively.\n", "- Fit a real binary event (2012-BLG-0406) with a custom sampler (`microlens_emcee`), tuning priors and walker settings.\n", "- Compare binary- vs. single-lens fits to diagnose when higher-order modeling is required for Roman challenge submissions." ] }, { "cell_type": "markdown", "metadata": { "id": "iI89ALihoHnG" }, "source": [ "## Introduction\n", "Binary microlensing events arise when two lens masses distort spacetime together, producing networks of caustics and sharp magnification spikes that cannot be described by the single-lens Paczyński curve.\n", "\n", "This tutorial teaches binary-lens fitting techniques and is primarily intended for participants in the Roman Microlensing Data Challenge 2026. You will explore different methods for fitting binary-lens light curves and develop an understanding of current modeling practices in the field, why binary-lens events are challenging to fit, and different ways to overcome those challenges.\n", "\n", "The workflow for this notebook consists of:\n", "* [Environment Setup](#Environment-Setup)\n", " - [Imports](#Imports)\n", " - [Notebook Assets](#notebook-assets)\n", "* [A Brief Introduction to Binary Lens Microlensing](#a-brief-review-of-binary-lens-microlensing)\n", "* [The Binary Lens Model](#the-binary-lens-model)\n", "* [Fitting a Binary Lens Model](#fitting-a-binary-lens-model)\n", "* [Additional Resources](#additional-resources)\n", "* [About the Notebook](#about-this-notebook)\n" ] }, { "cell_type": "markdown", "metadata": { "id": "x07V5VKSoHnH" }, "source": [ "## Environment Setup\n", "\n", "Choose whichever platform best matches your workflow:\n", "\n", "- **Colab** provides a zero-install option for quick experimentation.\n", "- **Roman Research Nexus** users should activate the `rges-pit-dc` kernel; all required packages are pre-installed there.\n", "- **Local** environments can be provisioned with the accompanying `env.yml` file for parity with the Nexus image.\n", "\n", "If you are new to Jupyter notebooks, the [Getting Started guide](https://medium.com/codingthesmartway-com-blog/getting-started-with-jupyter-notebook-for-python-4e7082bd5d46) offers a quick primer before you dive into this workflow." ] }, { "cell_type": "markdown", "metadata": { "id": "zOdZmgTmoHnH" }, "source": [ "\n", "> ### Running on the Roman Research Nexus\n", "> If you are following this notebook on the Roman Research Nexus (RRN) these packages are preinstalled on the RGES-PIT's kernel:\n", "> 1. In notebooks, open **Kernel > Change Kernel** menu and select `RGES PIT Nexus` or select the kernel via the Kernel menu or clicking on the \\\"kernel status display\\\" to the top right with the activity circle.\n", "> 1. In a terminal, if needed, run `mamba activate rges-pit-dc` to select the correct mamba environment.\"\n", "> 1. Store any files you generate in your home directory or team space; the `/roman/nexus-notebooks/notebooks` directory is read-only." ] }, { "cell_type": "markdown", "metadata": { "id": "Ki1-u1AFoHnH" }, "source": [ "\n", "> ### Running on Google Colab\n", "> 1. Open this notebook via the **Open in Colab** button in the workshop README or by pasting the GitHub URL into [https://colab.research.google.com](https://colab.research.google.com).\n", "> 1. Ensure you are signed into a Google account so Colab can save to your Drive.\n", "> 1. Use the `Runtime > Change runtime type` menu to select a GPU **off**, since this workflow is CPU bound, or select a local runtime.\n", "> 1. Remember to download your results or save a copy to Drive; Colab sessions are ephemeral.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "Yzsmv1m7oHnH", "outputId": "e174c280-6fb8-4e2c-cfc1-b12baaabd99f", "tags": [ "colab-only" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[?25l \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m0.0/519.0 kB\u001b[0m \u001b[31m?\u001b[0m eta \u001b[36m-:--:--\u001b[0m\r\u001b[2K \u001b[91m━━━━━━━━━\u001b[0m\u001b[90m╺\u001b[0m\u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m122.9/519.0 kB\u001b[0m \u001b[31m3.4 MB/s\u001b[0m eta \u001b[36m0:00:01\u001b[0m\r\u001b[2K \u001b[91m━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[90m╺\u001b[0m\u001b[90m━━━━━━━━━━━━━\u001b[0m \u001b[32m337.9/519.0 kB\u001b[0m \u001b[31m5.4 MB/s\u001b[0m eta \u001b[36m0:00:01\u001b[0m\r\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m519.0/519.0 kB\u001b[0m \u001b[31m5.6 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n", "\u001b[?25h\u001b[?25l \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m0.0/1.6 MB\u001b[0m \u001b[31m?\u001b[0m eta \u001b[36m-:--:--\u001b[0m\r\u001b[2K \u001b[91m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m\u001b[90m╺\u001b[0m\u001b[90m━━━━━━━\u001b[0m \u001b[32m1.3/1.6 MB\u001b[0m \u001b[31m38.2 MB/s\u001b[0m eta \u001b[36m0:00:01\u001b[0m\r\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m1.6/1.6 MB\u001b[0m \u001b[31m25.9 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n", "\u001b[?25h" ] } ], "source": [ "# COLAB-ONLY\n", "%pip install --quiet VBMicrolensing emcee corner multiprocess MulensModel ipympl\n", "\n", "# Colab widget bug work-around\n", "get_ipython().kernel.do_shutdown(restart=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "colab-only" ] }, "outputs": [], "source": [ "# COLAB-ONLY\n", "try:\n", " from google.colab import output\n", "except:\n", " pass" ] }, { "cell_type": "markdown", "metadata": { "id": "TA2DIel0oHnH" }, "source": [ "Please **save this notebook in a space you control** (GitHub fork, local repo, Nexus team storage, or Google Drive) so you can return to it later. Hosted environments such as Colab or the Nexus will not auto-save changes to the canonical repository." ] }, { "cell_type": "markdown", "metadata": { "id": "GpNyfyqqoHnH" }, "source": [ "### Imports" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "MqQ2d2cJoHnI" }, "outputs": [], "source": [ "#Import libraries\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "import scipy.optimize as op\n", "import emcee\n", "from multiprocess import Pool\n", "from multiprocess import cpu_count\n", "import corner\n", "from scipy.optimize import lsq_linear\n", "import ipywidgets as widgets\n", "from ipywidgets import FloatSlider, FloatLogSlider\n", "from IPython.display import display\n", "import MulensModel as mm\n", "import matplotlib.lines as mlines\n", "import VBMicrolensing as vbm\n", "\n", "try:\n", " from google.colab import sheets # this will only work if you are running on Colab\n", " %matplotlib ipympl \n", "except:\n", " %matplotlib inline\n", "\n", "import sys\n", "import shutil\n", "import urllib.request\n", "from pathlib import Path" ] }, { "cell_type": "markdown", "metadata": { "id": "RIg4SDyZoHnI" }, "source": [ "### Notebook Assets\n", "This session depends on a few small helper files (e.g. `microlens_emcee.py`, data tables, and static plots). The cell below makes sure those assets are available whether you run locally, on Colab, or on the Roman Research Nexus. Existing local copies are left untouched; missing files are fetched directly from GitHub.\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "p7I_0cw2oHnI", "outputId": "ab88507b-a233-4ef7-86e2-d1c6f4f24088" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "⬇️ Fetching microlens_emcee.py ...\n", "⬇️ Fetching 2012-BLG-0406_data.txt ...\n", "⬇️ Fetching 2012-BLG-0406_params.txt ...\n", "⬇️ Fetching Plots/binary_parameters.png ...\n", "⬇️ Fetching Plots/caustic_structures.png ...\n", "⬇️ Fetching Plots/caustic_topology.png ...\n" ] } ], "source": [ "ASSET_SOURCES = {\n", " Path(\"microlens_emcee.py\"): \"https://raw.githubusercontent.com/rges-pit/data-challenge-notebooks/main/AAS%20Workshop/Session%20C%3A%20Binary%20Lens/microlens_emcee.py\",\n", " Path(\"2012-BLG-0406_data.txt\"): \"https://raw.githubusercontent.com/rges-pit/data-challenge-notebooks/main/AAS%20Workshop/Session%20C%3A%20Binary%20Lens/2012-BLG-0406_data.txt\",\n", " Path(\"2012-BLG-0406_params.txt\"): \"https://raw.githubusercontent.com/rges-pit/data-challenge-notebooks/main/AAS%20Workshop/Session%20C%3A%20Binary%20Lens/2012-BLG-0406_params.txt\",\n", " Path(\"Plots/binary_parameters.png\"): \"https://raw.githubusercontent.com/rges-pit/data-challenge-notebooks/main/AAS%20Workshop/Session%20C%3A%20Binary%20Lens/Plots/binary_parameters.png\",\n", " Path(\"Plots/caustic_structures.png\"): \"https://raw.githubusercontent.com/rges-pit/data-challenge-notebooks/main/AAS%20Workshop/Session%20C%3A%20Binary%20Lens/Plots/caustic_structures.png\",\n", " Path(\"Plots/caustic_topology.png\"): \"https://raw.githubusercontent.com/rges-pit/data-challenge-notebooks/main/AAS%20Workshop/Session%20C%3A%20Binary%20Lens/Plots/caustic_topology.png\",\n", "}\n", "\n", "\n", "def ensure_asset(rel_path: Path, url: str) -> None:\n", " \"\"\"Download asset if it is missing from the working directory.\"\"\"\n", " rel_path = Path(rel_path)\n", " if rel_path.exists():\n", " print(f\"✔ {rel_path} already present\")\n", " return\n", "\n", " rel_path.parent.mkdir(parents=True, exist_ok=True)\n", " print(f\"⬇️ Fetching {rel_path} ...\")\n", " with urllib.request.urlopen(url) as response, open(rel_path, \"wb\") as outfile:\n", " shutil.copyfileobj(response, outfile)\n", "\n", "\n", "for path, url in ASSET_SOURCES.items():\n", " ensure_asset(path, url)\n", "\n", "module_dir = Path(\"microlens_emcee.py\").resolve().parent\n", "if str(module_dir) not in sys.path:\n", " sys.path.insert(0, str(module_dir))\n", " print(f\"🔗 Added {module_dir} to sys.path\")" ] }, { "cell_type": "markdown", "metadata": { "id": "0QjBOo6OoHnI" }, "source": [ "## A Brief Review of Binary Lens Microlensing" ] }, { "cell_type": "markdown", "metadata": { "id": "1xvw3WNioHnI" }, "source": [ "Finding the magnification of a binary lens event at any point in the source's relative trajectory requires solving the lens equation for the system. The lens equation for a binary point-mass lens maps the angular position of the source to the angular image positions for a given lens configuration. It can be written in complex coordinates as,\n", "\n", "$$z_{\\mathrm{s}}=z-\\frac{1}{1+q}\\left(\\frac{1}{\\bar{z}-\\bar{z}_1}+\\frac{q}{\\bar{z}-\\bar{z}_2}\\right)$$\n", "\n", "\n", "where $z_{\\mathrm{s}}$ is the position of the source in complex coordinates ($z = x + iy$), $z$ is the position of the image, $z_1$ and $z_2$ are the positions of the lens components, and $q = M_2/M_1$ is the mass ratio between the lens components. All angular positions are scaled to the angular Einstein ring radius ($\\theta_E$).\n", "\n", "The magnification of an image in this formalism is given by the inverse of the Jacobian of the mapping above evaluated at the image position,\n", "\n", "$$\\mu = \\left.\\frac{1}{\\operatorname{det} J}\\right|_{z=z_j}, \\operatorname{det} J \\equiv \\frac{\\partial\\left(x_1, x_2\\right)}{\\partial\\left(y_1, y_2\\right)}$$\n", "\n", "The magnification of images diverges when $\\operatorname{det} J = 0$. The set of all image positions where this happens form closed *critical curves*, and the corresponding set of source positions define closed curves called *caustics*. Although the magnification of a point source diverges along caustic curves, that of a real-finite sized source is large but finite.\n", "\n", "Caustics also form boundaries of regions with different numbers of images. Microlensing of a source inside the region bound by caustics always produces two more images than when the source is outside. The topology of caustics depends on the projected separation of the lens components in units of the Einstein ring radius (s) and their mass ratio (q), and is classified into three types: *close*, *resonant/intermediate*, and *wide*.\n", "\n", "
\n", "
Fig 1: Close, resonant, and intermediate caustics for q = 0.6.
\n", "\n", "For a fixed q, the boundaries of these topologies are given by,\n", "\n", "\n", "$$\\frac{q}{(1+q)^2}=\\frac{\\left(1-s_c\\right)^3}{27 s_c^8}, \\, s_w=\\frac{\\left(1+q^{1 / 3}\\right)^{3 / 2}}{(1+q)^{1 / 2}}$$\n", "\n", "if $s \\le s_c$ (close topology), there is one diamond-shaped caustic at the center of mass of the binary lens and two triangular caustics; if $s_c \\le s \\le s_w$ (resonant topology), there is one large caustic at the center of mass; and if $s \\ge s_w$ (wide topology) there are two diamond-shaped caustics.\n", "\n", "
\n", "
Fig 2: Boundaries of the different caustic topologies as a function of the mass ratio q of the binary lens (Gaudi, 2010).
\n", "\n", "Microlensing events where the source crosses a caustic have sharp peaks in magnification on top of the smooth single lens magnification curve (as we will see in the next section) - a characteristic feature of binary lens events which can be used to easily distinguish them from single lens events.\n" ] }, { "cell_type": "markdown", "metadata": { "id": "jeA4KkVsoHnI" }, "source": [ "\n", "\n", "
\n", "

Exercise 1

\n", "

Use the interactive s and q slider to plot caustics for different lens configurations. Plot all three caustic topologies and see how caustics transition from one topology to another by changing s. See how caustics change with q, and compare caustics for a planetary mass companion and a stellar mass companion.

\n", "
\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 436 }, "id": "2qK34ViHoHnI", "outputId": "b24718d9-670b-446c-8b4b-ded4e42f4719" }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "edc53adaf3824ba79e0ad01fd8e30255", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(FloatSlider(value=-0.5, description='Log(s)', max=1.0, min=-1.0, step=0.05), FloatSlider(value=…" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "7937ec495ba6442082ba1d1d6ea0b53e", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Output()" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#Use this code to plot caustics with a slider for s and q\n", "# Plot limits and styling\n", "X_MIN, X_MAX = -1.1, 1.1\n", "Y_MIN, Y_MAX = -1.1, 1.1\n", "\n", "\n", "def plot_binary_caustics(log_s, log_q):\n", " s = 10**log_s\n", " q = 10**log_q\n", "\n", " # Masses normalized so that total binary mass = 1\n", " m1 = 1.0 / (1.0 + q)\n", " m2 = q * m1\n", "\n", " # Lens positions (binary lens axis along the x axis with the center of mass at the origin and the more massive lens component on the left)\n", " z1x = -q * s / (1.0 + q)\n", " z2x = s / (1.0 + q)\n", "\n", "\n", " # Binary lens caustic and critical curves using MulensModel\n", " binary_caustic = mm.CausticsBinary(s=s, q=q)\n", " bx, by = binary_caustic.get_caustics()\n", " binary_critcurves = binary_caustic.critical_curve\n", " bcritx, bcrity = binary_critcurves.x, binary_critcurves.y\n", "\n", " # Draw figure\n", " fig = plt.figure(figsize=(5, 5), dpi=150)\n", " ax = fig.add_subplot(1, 1, 1)\n", " plt.subplots_adjust(top=0.9, bottom=0.08, right=0.95, left=0.1, hspace=0.1, wspace=0.2)\n", "\n", " ax.set_xlim(X_MIN, X_MAX)\n", " ax.set_ylim(Y_MIN, Y_MAX)\n", " x_left, x_right = ax.get_xlim()\n", " y_low, y_high = ax.get_ylim()\n", " ax.set_aspect(abs((x_right - x_left) / (y_high - y_low)))\n", "\n", " # Plot lenses\n", " ax.scatter([z1x], [0.0], marker='o', s=8, color='goldenrod', label='Lens 1')\n", " ax.scatter([z2x], [0.0], marker='o', s=8*q, color='red', label='Lens 2')\n", "\n", " # Plot caustics\n", " ax.scatter(bx, by, marker='.', s=0.01, color='blue', alpha=0.8)\n", " #Plot critical curves\n", " ax.scatter(bcritx, bcrity, color='gray', s = 0.01, alpha=0.5, marker='.')\n", "\n", " #Add labels for critical curves and caustics\n", " gray_line = mlines.Line2D([], [], color='gray', linestyle='-', label='Critical curve')\n", " blue_line = mlines.Line2D([], [], color='blue', linestyle='-', label='Binary caustic')\n", "\n", " #Create legend\n", " handles, labels = ax.get_legend_handles_labels()\n", " handles.append(gray_line)\n", " handles.append(blue_line)\n", " labels.append(gray_line.get_label())\n", " labels.append(blue_line.get_label())\n", " legend = ax.legend(handles, labels, loc=\"upper left\", fontsize=9)\n", " ax.add_artist(legend)\n", "\n", " ax.set_title(f's={s:.3f}, q={q:.2e}')\n", " plt.show()\n", "\n", "\n", "# Sliders\n", "log_s = FloatSlider(value=-0.5, min=-1, max=1, step=0.05, description='Log(s)')\n", "log_q = FloatSlider(value=-3, min=-4, max=0, step=0.1, description='Log(q)')\n", "\n", "\n", "ui = widgets.VBox([log_s, log_q])\n", "out = widgets.interactive_output(\n", " plot_binary_caustics, {\"log_s\": log_s, \"log_q\": log_q}\n", ")\n", "\n", "display(ui, out)\n" ] }, { "cell_type": "markdown", "metadata": { "id": "ZZ8sEZLpoHnI" }, "source": [ "## The Binary Lens Model" ] }, { "cell_type": "markdown", "metadata": { "id": "p7wdCZRPoHnI" }, "source": [ "We need at least three additional parameters to describe a binary lens, along with the four single lens parameters ($t_0$, $u_0$, $t_E$, and $\\rho$). These are:\n", "\n", "$s$: Projected separation of the lens components in units of the Einstein ring radius $r_E$\n", "\n", "$q$: Mass ratio between the two lens components\n", "\n", "$\\alpha$: Angle that the source trajectory makes with the binary lens axis\n", "\n", "\n", "
\n", "
Fig. 3: Binary lens event parameters
\n", "\n", "The plot above shows the general set up of a binary lens microlensing event, and the conventions that we will use in the rest of this notebook. The two lenses are placed on the X-axis with the more massive lens on the left. The center of mass of the system is at the origin. $\\alpha$ is the angle between the vector connecting $m_1$ to $m_2$ and the source trajectory. $\\tau$ is defined the same way as the single lens case, $\\tau = (t - t_0)/t_E$, and is used to parameterize the trajectory.\n", "\n", "The Einstein ring radius is defined using the combined mass of the system: $\\theta_E = \\sqrt{\\kappa \\ (m_1 + m_2) \\ \\pi_{rel}}$. Unlike single lens events, the Einstein ring is not a ring around which the lensed images are produced but is instead a mathematical construct used to give relative scale to the model. In a binary lens system, images form around critical curves in the image plane, which envelope the lens masses.\n", "\n", "Let us now see what the magnification curve for a binary lens event looks like." ] }, { "cell_type": "markdown", "metadata": { "id": "V5UVDzDPoHnI" }, "source": [ "\n", "
\n", "

Exercise 2

\n", "

Complete the function get_mag() defined below to calculate magnifications given a binary lens model.

\n", "
\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "AqIjDKOWoHnI" }, "outputs": [], "source": [ "def get_mag(x0, p, t):\n", " '''\n", " Get the magnification curve for a given model and parameters.\n", "\n", " Parameters\n", " ----------\n", " x0: np.ndarray with initial parameters (t0, tE, u0, rho, s, q, alpha)\n", " p: list[str] with parameter names\n", " t: np.ndarray with data epochs\n", "\n", " Returns\n", " -------\n", " mag: np.ndarray with magnification curve\n", " '''\n", " # We will use the VBMicrolensing package to calculate magnifications for different source positions.\n", " VBM = vbm.VBMicrolensing()\n", " #create a dictionary of parameters\n", " params = dict(zip(p, x0))\n", " paramsc = params.copy()\n", " tau = (t - paramsc['t0']) / paramsc['tE']\n", " # Single lens model\n", " if len(p) < 7 :\n", " ul = np.sqrt(tau**2 + paramsc['u0']**2)\n", " if 'logrho' in paramsc.keys():\n", " paramsc['rho'] = 10**(paramsc['logrho'])\n", " mag = np.array([VBM.ESPLMag2(u, paramsc['rho']) for u in ul])\n", " # Binary lens model\n", " elif len(p) >= 7:\n", " ''' ------------------------------------------------------------ '''\n", " ''' Enter your code here to calculate source positions as a function of tau using the plot in Fig. 3 as reference. '''\n", " salpha = np.sin(np.radians(paramsc['alpha']))\n", " calpha = np.cos(np.radians(paramsc['alpha']))\n", " xs = -paramsc['u0'] * salpha + tau * calpha\n", " ys = paramsc['u0'] * calpha + tau * salpha\n", " ''' ------------------------------------------------------------ '''\n", " if 'logs' in paramsc.keys():\n", " paramsc['s'] = 10**(paramsc['logs'])\n", " if 'logq' in paramsc.keys():\n", " paramsc['q'] = 10**(paramsc['logq'])\n", " if 'logrho' in paramsc.keys():\n", " paramsc['rho'] = 10**(paramsc['logrho'])\n", " mag = np.array([VBM.BinaryMag2(paramsc['s'], paramsc['q'], xs[i], ys[i], paramsc['rho']) for i in range(len(ys))])\n", " return mag" ] }, { "cell_type": "markdown", "metadata": { "id": "hyHieXrJoHnI" }, "source": [ "\n", "
\n", "

Exercise 3

\n", "

The cell below uses this function to plot the magnification curve, caustics and source trajectory for a given binary lens model. Vary the parameters to see how the light curve changes. See what happens when the source trajectory crosses a caustic.

\n", "
\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "ehCtZF9noHnI" }, "outputs": [], "source": [ "#Binary lens model\n", "s = 1.4\n", "q = 0.01\n", "alpha = 30.0\n", "t0 = 0.0\n", "tE = 30.0\n", "u0 = -0.53\n", "rho = 0.001\n", "\n", "#time array\n", "t = np.linspace(t0 - 1.5*tE, t0 + 1.5*tE, 5000)\n", "tau = (t - t0) / tE\n", "#source positions\n", "xs = -u0 * np.sin(np.radians(alpha)) + tau * np.cos(np.radians(alpha))\n", "ys = u0 * np.cos(np.radians(alpha)) + tau * np.sin(np.radians(alpha))\n", "\n", "\n", "#get caustics (using mulensmodel)\n", "binary_caustic = mm.CausticsBinary(s=s, q=q)\n", "bx, by = binary_caustic.get_caustics()\n", "\n", "#lens positions\n", "z1x = -q * s / (1.0 + q)\n", "z2x = s / (1.0 + q)\n", "\n", "#magnification curve\n", "mag = get_mag(np.array([t0, tE, u0, rho, s, q, alpha]), ['t0', 'tE', 'u0', 'rho', 's', 'q', 'alpha'], t)\n", "\n", "#plot caustics and the light curve in two subplots\n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 4), dpi = 150)\n", "ax1.scatter(t, mag, color='red', s=0.1)\n", "ax1.set_xlabel('t (days)')\n", "ax1.set_ylabel('Magnification')\n", "ax2.plot(xs, ys, color='black', linestyle='--', label='source trajectory')\n", "#add an arrow along the source trajectory\n", "ax2.arrow(xs[len(xs)//2 + 50], ys[len(ys)//2 + 50], xs[len(xs)//2] - xs[len(xs)//2 - 10], ys[len(ys)//2] - ys[len(ys)//2 - 10], head_width=0.07, head_length=0.1, color='black')\n", "ax2.scatter(z1x, 0.0, color='red', s=20, label='lens 1')\n", "ax2.scatter(z2x, 0.0, color='orange', s=10, label='lens 2')\n", "ax2.scatter(bx, by, color='blue', s=0.05)\n", "blue_line = mlines.Line2D([], [], color='blue', linestyle='-', label='caustic')\n", "handles, labels = ax2.get_legend_handles_labels()\n", "handles.append(blue_line)\n", "labels.append(blue_line.get_label())\n", "ax2.legend(handles, labels, fontsize=8)\n", "ax2.set_xlabel(r'x/$\\theta_E$')\n", "ax2.set_ylabel(r'y/$\\theta_E$')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "id": "JaYDb3adoHnI" }, "source": [ "## Fitting a Binary Lens Model" ] }, { "cell_type": "markdown", "metadata": { "id": "djsXcPqIoHnI" }, "source": [ "Now that we have the tools to create a binary lens model, let us get our hands dirty with real microlensing data and to try to fit a binary lens model to it! The cell below loads a light curve from the ground-based OGLE microlensing survey." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "QyacOqMaoHnI" }, "outputs": [], "source": [ "\n", "#Load data (row 1 is header)\n", "df = pd.read_csv('2012-BLG-0406_data.txt', sep='\\t', header=0)\n", "\n", "#load parameters into a dictionary\n", "params = np.loadtxt('2012-BLG-0406_params.txt', unpack=True, skiprows=1, dtype=str)\n", "params_dict = dict(zip(params[0], np.float64(params[1])))\n", "\n", "#select only data points where t is in range t0 - 4*tE and t0 + 3*tE\n", "t0i = params_dict['Tmax'] - 2450000\n", "tEi = params_dict['tau']\n", "u0i = params_dict['umin']\n", "df = df[(df['HJD-2450000'] > t0i - 4*tEi) & (df['HJD-2450000'] < t0i + 4*tEi)]\n", "\n", "#Plot data\n", "fig, ax = plt.subplots(figsize=(6, 4), dpi = 150)\n", "ax.errorbar(df['HJD-2450000'], df['I_Mag'], yerr=df['Err'], fmt='o', color='black', label='data', markersize=1)\n", "ax.set_xlabel('HJD - 2450000')\n", "ax.set_ylabel('Mag')\n", "#Invert the y-axis\n", "ax.invert_yaxis()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "id": "TnBpkb_SoHnI" }, "source": [ "Notice that the light curve has two additional sharp peaks in magnification, indicative of a caustic crossing event. Therefore, this event is most likely caused by a binary lens system.\n", "\n", "In the following cells, we will try different methods of fitting a binary lens model to this light curve. But first, we must convert the light curve to units of flux from magnitude. We will use a dummy zero point for this." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "vfqMGp8voHnJ" }, "outputs": [], "source": [ "#Convert from magnitudes to fluxes\n", "def mag2flux(mag: np.ndarray,\n", " zp: float = 25.0\n", " ) -> np.ndarray:\n", " '''\n", " Converts magnitude values (array) to flux.\n", "\n", " Parameters\n", " ----------\n", " mag : np.ndarray\n", " Array of magnitudes.\n", " zp : float, optional\n", " Zero point of the magnitude system (default is 25.0).\n", "\n", " Returns\n", " -------\n", " f : np.ndarray\n", " Array of fluxes.\n", " '''\n", " f = 10.0**((mag - zp) / -2.5)\n", " return f\n", "\n", "def mag2flux_err(mag: np.ndarray, mag_err: np.ndarray,\n", " zp: float = 25.0):\n", " '''\n", " Converts magnitude values and errors (array) to flux using error propagation formula.\n", "\n", " Parameters\n", " ----------\n", " mag : np.ndarray\n", " Array of magnitudes.\n", " mag_err : np.ndarray\n", " Array of magnitude errors.\n", " zp : float, optional\n", " Zero point of the magnitude system (default is 25.0).\n", "\n", " Returns\n", " -------\n", " flux_err : np.ndarray\n", " Array of flux errors.\n", " '''\n", " dfdmag = -0.4 * np.log(10) * 10.**(0.4*(zp-mag))\n", " flux_err = np.sqrt(dfdmag**2 * mag_err**2)\n", " return flux_err\n", "\n", "df['I_band_flux'] = mag2flux(df['I_Mag'], 25.0)\n", "df['I_band_flux_err'] = mag2flux_err(df['I_Mag'], df['Err'], 25.0)" ] }, { "cell_type": "markdown", "metadata": { "id": "XhnkJDzpoHnJ" }, "source": [ "So far, we have been using 7 parameters to calculate a binary lens model. But this will only give us a model magnification curve. In order to fit this to actual data, we need to convert magnifications to fluxes. The light curve of a microlensing event is given by:\n", "\n", "$$ F = F_s A(t) + F_b $$\n", "\n", "Where $F_s$ is the source flux and $F_b$ is the blended flux from the lens, and any other unresolved stars in the PSF of the source. Given a magnification model A, we can solve for $F_s$ and $F_b$ using linear least squares" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "JaOvOQrKoHnJ" }, "outputs": [], "source": [ "def calc_Fs(modelmag: np.ndarray, f: np.ndarray, sig2: np.ndarray) -> tuple[float, float]:\n", " '''\n", " Solves for the flux parameters for a given model using least squares.\n", "\n", " Parameters\n", " ----------\n", " model : np.ndarray\n", " Model magnification curve.\n", " f : np.ndarray\n", " Observed flux values.\n", " sig2 : np.ndarray\n", " Flux errors.\n", "\n", " Returns\n", " -------\n", " FS : float\n", " Source flux.\n", " FB : float\n", " Blend flux.\n", " '''\n", "\n", " \"\"\"\n", " Solves for FS, FB by weighted least squares with FB >= 0 using lsq_linear.\n", " \"\"\"\n", " modelmag = np.asarray(modelmag, float).ravel()\n", " f = np.asarray(f, float).ravel()\n", " sig2 = np.asarray(sig2, float).ravel()\n", "\n", " w = 1.0 / np.sqrt(sig2)\n", " A = np.column_stack([modelmag, np.ones_like(modelmag)])\n", " Aw = A * w.reshape(-1, 1)\n", " bw = f * w\n", "\n", " res = lsq_linear(Aw, bw)\n", " FS, FB = res.x\n", " return float(FS), float(FB)" ] }, { "cell_type": "markdown", "metadata": { "id": "8o02klSYoHnJ" }, "source": [ "To fit a model to the light curve, we will minimize the $\\chi^2$ between the data and the model" ] }, { "cell_type": "markdown", "metadata": { "id": "nUK4R0CYoHnJ" }, "source": [ "\n", "
\n", "

Exercise 4

\n", "

Complete the χ² function in the cell below

\n", "
\n", "
\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "wuXwwwQGoHnJ" }, "outputs": [], "source": [ "def chi2(x0: np.ndarray, p: list, t: np.ndarray, f: np.ndarray, sig: np.ndarray) -> float:\n", " '''\n", " Calculates the chi squared value for a given model and parameters.\n", "\n", " Parameters\n", " ----------\n", " x0 : np.ndarray\n", " Initial parameters.\n", " p : List[str]\n", " List of parameter names.\n", " t : np.ndarray\n", " Data epochs.\n", " f : np.ndarray\n", " Observed flux values.\n", " sig : np.ndarray\n", " Flux errors.\n", "\n", " Returns\n", " -------\n", " chi2 : float\n", " Chi squared value.\n", " '''\n", " '''----------------------------------------------------------'''\n", " ''' Enter your code here to calculate the chi squared value '''\n", " mag = get_mag(x0, p, t)\n", " FS, FB = calc_Fs(mag, f, sig**2)\n", " model = FS * mag + FB\n", "\n", " chi2_value = np.sum((f - model)**2 / sig**2)\n", " '''----------------------------------------------------------'''\n", " return chi2_value" ] }, { "cell_type": "markdown", "metadata": { "id": "h3t6bjvToHnJ" }, "source": [ "A function to perform basic downhill simplex minimization using scipy's Nelder-Mead implementation:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "IVogQONMoHnJ" }, "outputs": [], "source": [ "def simple_fit(df: pd.DataFrame, x0: np.ndarray, p: list[str], bounds: dict | None = None,\n", " fixed_params: dict | None = None, maxfev: int = 1000):\n", " '''\n", " Fit a single or binary lens model using downhill simplex method and find delta chi^2\n", "\n", " Parameters\n", " ----------\n", " df : pandas.DataFrame\n", " Dataframe with data epochs, flux, and flux errors\n", " x0 : np.ndarray\n", " Initial parameters (e.g., t0, tE, u0, rho)\n", " p : list[str]\n", " Parameter names corresponding to x0\n", " bounds : dict | None, optional\n", " Optional map of parameter name -> (min, max). If None, use defaults.\n", " fixed_params : dict | None, optional\n", " Parameters to freeze as {name: value}. Only other parameters are varied.\n", " maxfev : int, optional\n", " Maximum number of function evaluations of the minimizer. Default is 1000.\n", "\n", " Returns\n", " -------\n", " fin_params : dict\n", " Dictionary with final parameters\n", " fin_chi2 : float\n", " Final chi-squared value\n", " delta_chi2 : float\n", " Delta chi-squared between best fit model and perfect model\n", " '''\n", " t = np.array(df['HJD-2450000'])\n", " f = np.array(df['I_band_flux'])\n", " sig = np.array(df['I_band_flux_err'])\n", "\n", " # Default bounds\n", " default_bounds = {\n", " 't0': (t.min(), t.max()),\n", " 'tE': (0.1, t.max()-t.min()),\n", " 'u0': (-10.0, 10.0),\n", " 'rho': (1e-5, 1.0),\n", " 's': (1e-2, 20.0),\n", " 'q': (1e-7, 1.0),\n", " 'alpha': (0.0, 360.0),\n", " 'logs': (-2, 2),\n", " 'logq': (-7, 0),\n", " 'logrho': (-5, 0),\n", " }\n", "\n", " # Merge user bounds with defaults\n", " bounds_map = default_bounds.copy()\n", " if bounds is not None:\n", " bounds_map.update({k: tuple(v) for k, v in bounds.items() if k in bounds_map})\n", "\n", " fixed_params = fixed_params or {}\n", " fixed_params = {k: float(v) for k, v in fixed_params.items()}\n", "\n", " # Build free/fixed split\n", " free_param_names = [name for name in p if name not in fixed_params]\n", " if len(free_param_names) == 0:\n", " # Nothing to optimize, just evaluate chi2 at fixed x0 values\n", " fin_params = dict(zip(p, x0))\n", " # Override with provided fixed values\n", " fin_params.update(fixed_params)\n", " x_full = np.array([fin_params[name] for name in p])\n", " fin_chi2 = chi2(x_full, p, t, f, sig)\n", " delta_chi2 = fin_chi2 - (len(f) - len(p))\n", " return fin_params, fin_chi2, delta_chi2\n", "\n", " # Initial free vector in the order of free_param_names\n", " x0_map = dict(zip(p, x0))\n", " x0_free = np.array([x0_map[name] for name in free_param_names])\n", "\n", " # Bounds for free parameters\n", " bounds_free = [bounds_map[name] for name in free_param_names]\n", "\n", " # Objective over free parameters that reconstructs full vector\n", " def objective(theta_free: np.ndarray) -> float:\n", " full_map = {**x0_map, **fixed_params}\n", " for j, name in enumerate(free_param_names):\n", " full_map[name] = float(theta_free[j])\n", " x_full = np.array([full_map[name] for name in p])\n", " return chi2(x_full, p, t, f, sig)\n", "\n", " result = op.minimize(\n", " objective, x0_free, bounds=bounds_free,\n", " method='Nelder-Mead', options={'xatol': 1e-8, 'fatol': 1e-8, 'adaptive': False, 'maxfev': maxfev}\n", " )\n", "\n", " if isinstance(result.fun, np.ndarray):\n", " if result.fun.ndim == 0:\n", " fin_chi2 = float(result.fun)\n", " else:\n", " fin_chi2 = result.fun[0]\n", " else:\n", " fin_chi2 = result.fun\n", "\n", " # Reconstruct final full params\n", " fin_map = {**x0_map, **fixed_params}\n", " for j, name in enumerate(free_param_names):\n", " fin_map[name] = float(result.x[j])\n", "\n", " fin_params = {name: fin_map[name] for name in p}\n", " dof_params = len(free_param_names)\n", " delta_chi2 = fin_chi2 - (len(f) - dof_params)\n", "\n", " return fin_params, fin_chi2, delta_chi2" ] }, { "cell_type": "markdown", "metadata": { "id": "yXBm_RVWoHnJ" }, "source": [ "Let us use this function to first fit a single lens model to the light curve. Although this particular light curve is very obviously not a single lens event, this can give us a decent first guess/ estimate of the single lens parameters for the event ($t_0$, $u_0$, $t_E$, $\\rho$). We can eyeball the initial guesses for $t_0$ and $t_E$ from the light curve and get pretty close to the right answer, and refine the initial guess in successive iterations. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "f_sGZPPNoHnJ" }, "outputs": [], "source": [ "#Enter your initial guesses for the parameters\n", "t0i = 6150\n", "u0i = 0.6\n", "tEi = 50\n", "rhoi = 0.001\n", "x0 = np.array([t0i, tEi, u0i, rhoi])\n", "p = ['t0', 'tE', 'u0', 'rho']\n", "t = np.array(df['HJD-2450000'])\n", "init_modelmag = get_mag(x0, p, t)\n", "source_flux, blend_flux = calc_Fs(init_modelmag, df['I_band_flux'], df['I_band_flux_err']**2)\n", "init_model = source_flux * init_modelmag + blend_flux\n", "\n", "fin_params, fin_chi2, delta_chi2 = simple_fit(df, x0, p)\n", "x = np.array([fin_params['t0'], fin_params['tE'], fin_params['u0'], fin_params['rho']])\n", "print(\"Final parameters: \", fin_params)\n", "print(\"Reduced chi2: \", fin_chi2/(len(df)-len(p)))\n", "print(\"Delta chi2: \", delta_chi2)\n", "fin_modelmag = get_mag(x, p, df['HJD-2450000'])\n", "source_flux, blend_flux = calc_Fs(fin_modelmag, df['I_band_flux'], df['I_band_flux_err']**2)\n", "print(\"Source flux: \", source_flux, \"Blend flux: \", blend_flux)\n", "final_model = source_flux * fin_modelmag + blend_flux\n", "\n", "#plot the data and the model\n", "fig, ax = plt.subplots(figsize=(8, 6))\n", "plt.errorbar(df['HJD-2450000'], df['I_band_flux'], yerr=df['I_band_flux_err'], fmt='o', color='black', label='data', markersize=1)\n", "plt.plot(df['HJD-2450000'], init_model, color='blue', label='Initial model')\n", "plt.plot(df['HJD-2450000'], final_model, color='red', label=rf'Best fit model - $\\Delta \\chi^2 = {delta_chi2:.2f}$')\n", "plt.xlabel('HJD - 2450000')\n", "plt.ylabel('Flux')\n", "plt.legend(loc='upper left', framealpha=0.5)\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": { "id": "77FpZDEroHnJ" }, "source": [ "Great! The final model looks like a fairly good fit to the light curve minus the perturbation due to the lens companion. Notice that the blend flux has turned out to be negative which is unphysical and hence this cannot be a valid model for this light curve. Now let us try to fit a binary lens model using the same function, from random initial guesses for the additional binary parameters." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "ARbeDhamoHnJ" }, "outputs": [], "source": [ "#Start from best fit single lens parameters\n", "t0i = 6141.66\n", "tEi = 54.27\n", "u0i = 0.6899\n", "rhoi = 0.00096\n", "#And random guesses for the binary parameters\n", "logsi = np.random.uniform(-2, 1)\n", "logqi = np.random.uniform(-5, 0)\n", "alphai = np.random.uniform(0.0, 360.0)\n", "\n", "x0 = np.array([t0i, tEi, u0i, rhoi, logsi, logqi, alphai])\n", "p = ['t0', 'tE', 'u0', 'rho', 'logs', 'logq', 'alpha']\n", "init_modelmag = get_mag(x0, p, t)\n", "source_flux, blend_flux = calc_Fs(init_modelmag, df['I_band_flux'], df['I_band_flux_err']**2)\n", "init_model = source_flux * init_modelmag + blend_flux\n", "\n", "fin_params, fin_chi2, delta_chi2 = simple_fit(df, x0, p, maxfev=1000)\n", "fin_params['s'] = 10**(fin_params['logs'])\n", "fin_params['q'] = 10**(fin_params['logq'])\n", "x = np.array([fin_params['t0'], fin_params['tE'], fin_params['u0'], fin_params['rho'], fin_params['s'], fin_params['q'], fin_params['alpha']])\n", "print(\"Final parameters: \", fin_params)\n", "print(\"Reduced chi2: \", fin_chi2/(len(df)-len(p)))\n", "print(\"Delta chi2: \", delta_chi2)\n", "fin_modelmag = get_mag(x, p, t)\n", "source_flux, blend_flux = calc_Fs(fin_modelmag, df['I_band_flux'], df['I_band_flux_err']**2)\n", "final_model = source_flux * fin_modelmag + blend_flux\n", "print(\"Source flux: \", source_flux, \"\\nBlend flux: \", blend_flux)\n", "\n", "#plot the data and the model\n", "fig, ax = plt.subplots(figsize=(8, 6))\n", "plt.errorbar(df['HJD-2450000'], df['I_band_flux'], yerr=df['I_band_flux_err'], fmt='o', color='black', label='data', markersize=1)\n", "plt.plot(df['HJD-2450000'], init_model, color='blue', label='Initial model')\n", "plt.plot(df['HJD-2450000'], final_model, color='red', label=rf'Best fit model - $\\Delta \\chi^2 = {delta_chi2:.2f}$')\n", "plt.xlabel('HJD - 2450000')\n", "plt.ylabel('Flux')\n", "plt.legend(loc='upper left', framealpha=0.5)\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": { "id": "IGXIrAHWoHnQ" }, "source": [ "This was a disaster! A binary lens microlensing model is too complex to be fit by a simple gradient descent-like algorithm. The likelihood (or $\\chi^2$) space of a binary lens event is very rugged and riddled with deep local minima. So any gradient descent method will fail or get stuck in a local minimum unless we are very close to the true solution.\n", "\n", "But we don't have to start from a completely random guess for the binary lens parameters. The shape and poistion of the planetary/stellar companion perturbation can be used to get a rough estimate of the binary parameters following the procedure outlined in [Gaudi & Gould (1997)](https://ui.adsabs.harvard.edu/abs/1997ApJ...486..675G). This is further elaborated in [Analytic planet parameters notebook](./Day4_Gould_Loeb_planetary_event.ipynb) and the [Microlensing Source](https://www.microlensing-source.org/concept/extracting-parameters/) website.\n", "\n", "Let us try to fit the model from a better guess for the binary parameters" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "XMJJkTiGoHnQ" }, "outputs": [], "source": [ "logrhoi = -1.344\n", "logsi = 0.257\n", "logqi = -2.211\n", "alphai = 180 + np.degrees(0.829)\n", "\n", "x0 = np.array([t0i, tEi, u0i, logrhoi, logsi, logqi, alphai])\n", "p = ['t0', 'tE', 'u0', 'logrho', 'logs', 'logq', 'alpha']\n", "init_modelmag = get_mag(x0, p, t)\n", "source_flux, blend_flux = calc_Fs(init_modelmag, df['I_band_flux'], df['I_band_flux_err']**2)\n", "init_model = source_flux * init_modelmag + blend_flux\n", "\n", "fin_params, fin_chi2, delta_chi2 = simple_fit(df, x0, p)\n", "fin_params['s'] = 10**(fin_params['logs'])\n", "fin_params['q'] = 10**(fin_params['logq'])\n", "fin_params['rho'] = 10**(fin_params['logrho'])\n", "x = np.array([fin_params['t0'], fin_params['tE'], fin_params['u0'], fin_params['logrho'], fin_params['logs'], fin_params['logq'], fin_params['alpha']])\n", "print(\"Final parameters: \", fin_params)\n", "print(\"Reduced chi2: \", fin_chi2/(len(df)-len(p)))\n", "print(\"Delta chi2: \", delta_chi2)\n", "fin_modelmag = get_mag(x, p, t)\n", "source_flux, blend_flux = calc_Fs(fin_modelmag, df['I_band_flux'], df['I_band_flux_err']**2)\n", "final_model = source_flux * fin_modelmag + blend_flux\n", "print(\"Source flux: \", source_flux, \"\\nBlend flux: \", blend_flux)\n", "#plot the data and the model\n", "fig, ax = plt.subplots(figsize=(8, 6))\n", "plt.errorbar(df['HJD-2450000'], df['I_band_flux'], yerr=df['I_band_flux_err'], fmt='o', color='black', label='data', markersize=1)\n", "plt.plot(df['HJD-2450000'], init_model, color='blue', label='Initial model')\n", "plt.plot(df['HJD-2450000'], final_model, color='red', label=rf'Best fit model - $\\Delta \\chi^2 = {delta_chi2:.2f}$')\n", "plt.xlabel('HJD - 2450000')\n", "plt.ylabel('Flux')\n", "plt.legend(loc='upper left', framealpha=0.5)\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": { "id": "lX9siqXJoHnQ" }, "source": [ "Although this is a better fit, the initial condition is still not close enough to the true solution for this method to work.\n", "\n", "A better way to search the likelihood space is by using Markov Chain Monte Carlo (MCMC). The stochastic nature of this algorithm allows it to escape local minima where a gradient descent algorithm would get stuck.\n", "\n", "In the cells below we will implement the MCMC algorithm to fit a binary lens model and find posterior distributions for all parameters. We have defined a posterior probability function by using the $\\chi^2$ likelihood for a model, and flat priors as long as the parameters are within user-defined bounds.\n", "\n", "The [emcee](https://emcee.readthedocs.io/en/stable/) package is used to run the MCMC algorithm. It is called by a custom MCMC class for microlensing modelling that calculates various convergence staistics, determines burn-in period, finds the MLE solution and confidence intervals, and plots a corner plot of posteriors as well as chains for all walkers." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "0OLFi2PEoHnQ" }, "outputs": [], "source": [ "# Define the log prior function - used to set bounds on parameters\n", "def log_prior(x0, p, t, bounds: dict | None = None):\n", " ''' bounds on parameters\n", "\n", " Parameters\n", " ----------\n", " x0 : np.ndarray\n", " Parameter vector (e.g., t0, tE, u0, rho, ...)\n", " p : list[str]\n", " Corresponding parameter names\n", " t : np.ndarray\n", " Data epochs\n", " bounds : dict | None, optional\n", " Optional map of parameter name -> (min, max). If None, use defaults.\n", " '''\n", "\n", " default_bounds = {'t0': (t.min(), t.max()), 'tE': (0.1, t.max()-t.min()), 'u0': (-10.0, 10.0), 'rho': (1e-5, 1.0), 's': (1e-2, 20.0), 'q': (1e-7, 1.0), 'alpha': (0.0, 360.0), 'logs': (-2, 2), 'logq': (-7, 0), 'logrho': (-5, 0)}\n", " bounds_map = default_bounds.copy()\n", " if bounds is not None:\n", " bounds_map.update({k: tuple(v) for k, v in bounds.items() if k in bounds_map})\n", "\n", " for index, key in enumerate(p):\n", " lower, upper = bounds_map[key]\n", " if x0[index] < lower or x0[index] > upper:\n", " return -np.inf\n", " return 0.\n", "#Calculate the log likelihood using the chi2 function\n", "def log_likelihood(x0, p, t, f, sig):\n", " '''calculating chi2'''\n", " ll = -1./2. * chi2(x0, p, t, f, sig)\n", " return ll\n", "#Return the log probability by adding the log prior and log likelihood\n", "def log_probability(x0, p, t, f, sig, bounds: dict | None = None):\n", " lp = log_prior(x0, p, t, bounds=bounds)\n", " if not np.isfinite(lp):\n", " return -np.inf\n", " return lp + log_likelihood(x0, p, t, f, sig)\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "20STzTXxoHnQ" }, "outputs": [], "source": [ "from microlens_emcee import micro_mc\n", "# Set random seed for reproducibility\n", "np.random.seed(46)\n", "\n", "# Set up initial parameters for binary lens model\n", "logrhoi = -1.344\n", "logsi = 0.257\n", "logqi = -2.211\n", "alphai = 180 + np.degrees(0.829)\n", "\n", "x0 = np.array([t0i, tEi, u0i, logrhoi, logsi, logqi, alphai])\n", "p = ['t0', 'tE', 'u0', 'logrho', 'logs', 'logq', 'alpha']\n", "\n", "#Apply more stringent bounds on single lens parameters since we have a fairly good guess of what they are\n", "informed_guess_mcmc = micro_mc(df, x0, p, log_probability, bounds={'t0': (t0i - tEi, t0i + tEi), 'tE': (tEi - 20, tEi + 20), 'u0': (u0i - 0.2, u0i + 0.2)})\n", "\n", "#Run MCMC analysis with all plots enabled\n", "results = informed_guess_mcmc.perform_mcmc_analysis(steps=1000, walkers=100, param_scales=np.array([1e-1, 1e-1, 1e-2, 1e-1, 1e-2, 1e-2, 1e-1]),\n", " verbose=False, plot_corner=True, plot_fit=True,\n", " plot_traces=True, plot_convergence=False, n_threads=cpu_count()-1)" ] }, { "cell_type": "markdown", "metadata": { "id": "n4t4jho2oHnQ" }, "source": [ "Although the MLE solution looks like it is pretty close to the right answer, the posterior distributions for this run are all over the place! The walker chains of all parameters are split multiple times indicating that they are getting stuck in local minima and are not able to fully explore the parameter space.\n", "\n", "The best way to explore the phase space of all possible solutions is to do a **grid search**. This is the method used to analyze all individual planetary microlensing events so far. The idea of a grid search is to search for the best starting guesses on a large grid of binary parameters ($s$, $q$, and $\\alpha$). The best *n* solutions are picked from the results of this search and a glabal minimum search is conducted using MCMC starting from these parameters. This is the most robust method to ensure that we find the global minimum, although it is also very computationally expensive. The other advantage of this method is its ability to identify ***degenerate solutions***. Microlensing models have several degeneracies, some of which can be broken by the measurement of higher order effects or better light curve sampling, but some which are inherent and cannot be resolved by just observing the light curve. We will discuss degeneracies in greater detail in the following section.\n", "\n", "There are many different prescriptions for a grid search that researchers have used in the past. They differ in the parameters used to create the grid, the minimation routines used, the method used to calculate magnifications, and the procedure followed to identify the global minimum once the grid search is completed. In the following cells we demonstrate a very simple grid search algorithm on a coarse grid of ($s$, $q$, $\\alpha$). For each grid point, we find the best fit values of the single lens parameters ($t_0$, $t_E$, $u_0$, $\\rho$) by fitting a binary lens model using the simple_fit function. During the fitting, we only vary the single lens parameters and keep ($s$, $q$, $\\alpha$) fixed. For each ($s$, $q$) we then pick the $\\alpha$ value which produces the smallest $\\chi^2$ and store the $\\alpha$ and $\\chi^2$. From this final grid of ($s$, $q$) we can choose the best *n* grid points that produce the smallest $\\chi^2$ and run MCMC minimizations starting from those parameters. For the MCMC runs, all parameters are varied." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "Zip1MVd7oHnQ" }, "outputs": [], "source": [ "\"\"\"\n", "Binary Microlensing Grid Search Module\n", "=====================================\n", "\n", "This module performs a grid search over binary-lens microlensing parameters\n", "(log10(s), log10(q), alpha) using parallel execution. For each grid point,\n", "a local optimization is carried out over the remaining microlensing parameters\n", "(t0, tE, u0, rho), and the best-fit solution is selected based on Δχ².\n", "\n", "Parallelization is handled via `multiprocess.Pool`, and progress tracking\n", "is provided using `tqdm`.\n", "\"\"\"\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from multiprocess import Pool, cpu_count\n", "from tqdm.auto import tqdm\n", "\n", "\n", "# =============================================================================\n", "# Worker function\n", "# =============================================================================\n", "\n", "def _grid_point_worker(args):\n", " \"\"\"\n", " Evaluate a single (log10(s), log10(q)) grid point.\n", "\n", " For a fixed (logs, logq) pair, this function loops over all values of alpha,\n", " performs a local optimization using `simple_fit`, and selects the alpha\n", " value that minimizes Δχ².\n", "\n", " Parameters\n", " ----------\n", " args : tuple\n", " A packed tuple containing:\n", " - i, j : int\n", " Indices of the grid point in the logs–logq grid.\n", " - logs_val : float\n", " log10 of binary separation (s).\n", " - logq_val : float\n", " log10 of mass ratio (q).\n", " - alpha_vals : array-like\n", " Array of alpha values (degrees) to scan.\n", " - df : pandas.DataFrame\n", " Light curve data.\n", " - t0_init, u0_init, tE_init, logrho_init : float\n", " Initial guesses for the local optimization.\n", " - bounds : dict or None\n", " Bounds passed to `simple_fit`.\n", "\n", " Returns\n", " -------\n", " tuple\n", " (i, j, min_delta_chi2, best_alpha, best_params)\n", " \"\"\"\n", "\n", " (i, j, logs_val, logq_val, alpha_vals, df,\n", " t0_init, u0_init, tE_init, logrho_init,\n", " bounds) = args\n", "\n", " min_delta_chi2 = np.inf\n", " best_alpha = None\n", " best_params = None\n", "\n", " for alpha_val in alpha_vals:\n", " # Parameter vector\n", " x0 = np.array([\n", " t0_init, tE_init, u0_init, logrho_init,\n", " logs_val, logq_val, alpha_val\n", " ])\n", "\n", " param_names = ['t0', 'tE', 'u0', 'logrho', 'logs', 'logq', 'alpha']\n", "\n", " fixed_params = {\n", " 'logs': logs_val,\n", " 'logq': logq_val,\n", " 'alpha': alpha_val\n", " }\n", "\n", " # Local optimization (must exist in global scope)\n", " fin_params, fin_chi2, delta_chi2 = simple_fit(\n", " df, x0, param_names, bounds=bounds, fixed_params=fixed_params\n", " )\n", "\n", " if delta_chi2 < min_delta_chi2:\n", " min_delta_chi2 = delta_chi2\n", " best_alpha = alpha_val\n", " best_params = fin_params\n", "\n", " return (i, j, min_delta_chi2, best_alpha, best_params)\n", "\n", "\n", "# =============================================================================\n", "# Grid search driver\n", "# =============================================================================\n", "\n", "def grid_search_binary(df, logs_vals, logq_vals, alpha_vals,\n", " t0_init, u0_init, tE_init, logrho_init,\n", " bounds=None, verbose=False, n_jobs=None):\n", " \"\"\"\n", " Perform a parallel grid search over binary microlensing parameters.\n", "\n", " The grid is defined in (log10(s), log10(q)), and for each grid point,\n", " alpha is scanned internally. The best-fit solution is determined by\n", " minimizing Δχ² after local optimization of remaining parameters.\n", "\n", " For each (logs, logq) point on the grid, this function:\n", " 1. Loops through all alpha values\n", " 2. Optimizes t0, u0, tE, logrho using simple_fit with fixed binary params\n", " 3. Selects the alpha value that gives the minimum delta_chi2\n", " 4. Stores the best result for that (logs, logq) point\n", "\n", " Parameters\n", " ----------\n", " df : pd.DataFrame\n", " Data with 'HJD-2450000', 'I_band_flux', 'I_band_flux_err'\n", " logs_vals : array-like\n", " Grid values for log10(s)\n", " logq_vals : array-like\n", " Grid values for log10(q)\n", " alpha_vals : array-like\n", " Grid values for alpha (degrees)\n", " t0_init, u0_init, tE_init, logrho_init : float\n", " Initial values for optimization at each grid point\n", " bounds : dict, optional\n", " Parameter bounds for simple_fit\n", " verbose : bool, optional\n", " Print progress information\n", " n_jobs : int, optional\n", " Number of parallel processes to use. Defaults to cpu_count().\n", "\n", " Returns\n", " -------\n", " results : dict\n", " Contains:\n", " - 'logs_vals': 1D array of logs values\n", " - 'logq_vals': 1D array of logq values\n", " - 'delta_chi2_grid': 2D array of minimum delta_chi2 for each (logs,logq)\n", " - 'best_alpha_grid': 2D array of best alpha for each (logs,logq)\n", " - 'best_params_grid': 2D array of dicts with all best-fit params\n", " - 'sorted_results': List of dicts sorted by delta_chi2 (ascending)\n", " - 'fig': matplotlib figure with visualization\n", " \"\"\"\n", "\n", " if n_jobs is None:\n", " n_jobs = cpu_count()\n", "\n", " # Allocate result arrays\n", " delta_chi2_grid = np.zeros((len(logq_vals), len(logs_vals)))\n", " best_alpha_grid = np.zeros((len(logq_vals), len(logs_vals)))\n", " best_params_grid = np.empty((len(logq_vals), len(logs_vals)), dtype=object)\n", "\n", " total_points = len(logs_vals) * len(logq_vals)\n", "\n", " if verbose:\n", " print(f\"Processing {total_points} grid points using {n_jobs} processes...\")\n", "\n", " # Build task list\n", " tasks = []\n", " for i, logs_val in enumerate(logs_vals):\n", " for j, logq_val in enumerate(logq_vals):\n", " tasks.append((\n", " i, j, logs_val, logq_val, alpha_vals, df,\n", " t0_init, u0_init, tE_init, logrho_init,\n", " bounds\n", " ))\n", "\n", " # -------------------------------------------------------------------------\n", " # Parallel execution with progress bar\n", " # -------------------------------------------------------------------------\n", "\n", " results_list = []\n", "\n", " if n_jobs > 1:\n", " with Pool(processes=n_jobs) as pool:\n", " with tqdm(total=len(tasks),\n", " desc=\"Binary grid search\",\n", " unit=\"grid\",\n", " disable=not verbose) as pbar:\n", " for res in pool.imap_unordered(_grid_point_worker, tasks, chunksize=1):\n", " results_list.append(res)\n", " pbar.update(1)\n", " else:\n", " for task in tqdm(tasks, desc=\"Binary grid search\",\n", " unit=\"grid\", disable=not verbose):\n", " results_list.append(_grid_point_worker(task))\n", "\n", " # -------------------------------------------------------------------------\n", " # Populate result grids\n", " # -------------------------------------------------------------------------\n", "\n", " for i, j, min_dchi2, b_alpha, b_params in results_list:\n", " delta_chi2_grid[j, i] = min_dchi2\n", " best_alpha_grid[j, i] = b_alpha\n", " best_params_grid[j, i] = b_params\n", "\n", " if verbose:\n", " print(f\"Grid search complete! Processed {total_points} points.\")\n", "\n", " # -------------------------------------------------------------------------\n", " # Visualization\n", " # -------------------------------------------------------------------------\n", "\n", " logs_grid, logq_grid = np.meshgrid(logs_vals, logq_vals)\n", "\n", " fig, ax = plt.subplots(figsize=(10, 8))\n", " im = ax.pcolormesh(\n", " logs_grid, logq_grid,\n", " np.log10(delta_chi2_grid),\n", " cmap='viridis',\n", " shading='auto'\n", " )\n", "\n", " cbar = plt.colorbar(im, ax=ax)\n", " cbar.set_label(r'$\\log_{10}(\\Delta \\chi^2)$', fontsize=14)\n", "\n", " ax.set_xlabel('log10(s)', fontsize=14)\n", " ax.set_ylabel('log10(q)', fontsize=14)\n", " ax.set_title('Grid Search: Binary Lens Parameters', fontsize=16)\n", "\n", " min_idx = np.unravel_index(np.argmin(delta_chi2_grid), delta_chi2_grid.shape)\n", " ax.plot(logs_grid[min_idx], logq_grid[min_idx],\n", " 'r*', markersize=20, label='Global minimum')\n", "\n", " ax.legend()\n", " plt.tight_layout()\n", " plt.show()\n", "\n", " # -------------------------------------------------------------------------\n", " # Sorted results\n", " # -------------------------------------------------------------------------\n", "\n", " sorted_results = []\n", " for i in range(len(logs_vals)):\n", " for j in range(len(logq_vals)):\n", " sorted_results.append({\n", " 'logs': logs_vals[i],\n", " 'logq': logq_vals[j],\n", " 's': 10**logs_vals[i],\n", " 'q': 10**logq_vals[j],\n", " 'best_alpha': best_alpha_grid[j, i],\n", " 'delta_chi2': delta_chi2_grid[j, i],\n", " 'all_params': best_params_grid[j, i]\n", " })\n", "\n", " sorted_results.sort(key=lambda x: x['delta_chi2'])\n", "\n", " return {\n", " 'logs_vals': logs_vals,\n", " 'logq_vals': logq_vals,\n", " 'delta_chi2_grid': delta_chi2_grid,\n", " 'best_alpha_grid': best_alpha_grid,\n", " 'best_params_grid': best_params_grid,\n", " 'sorted_results': sorted_results,\n", " 'fig': fig\n", " }\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "9WNGpFz7oHnR" }, "outputs": [], "source": [ "#Run the grid search. Takes about ~ 8 mins to run on a Macbook Pro with 6 cores\n", "\n", "# Define grid ranges\n", "logs_vals = np.linspace(-0.4, 0.2, 10)\n", "logq_vals = np.linspace(-5, -2.5, 10)\n", "alpha_vals = np.linspace(0, 360, 20)\n", "\n", "# Initial values for optimization (from best fit single lens model)\n", "t0_init = t0i\n", "u0_init = u0i\n", "tE_init = tEi\n", "logrho_init = np.log10(rhoi)\n", "\n", "# Run grid search (sequential execution for notebook compatibility)\n", "results = grid_search_binary(\n", " df, logs_vals, logq_vals, alpha_vals,\n", " t0_init, u0_init, tE_init, logrho_init,\n", " bounds={'u0': (u0i - 2, u0i + 2), 'logrho': (-6, -0.5), 't0': (t0i - tEi, t0i + tEi), 'tE': (tEi - 20, tEi + 20)},\n", " verbose=True, n_jobs=cpu_count()-1\n", ")\n", "\n", "# Access the best results\n", "print(\"\\n\" + \"=\"*60)\n", "print(\"GRID SEARCH RESULTS\")\n", "print(\"=\"*60)\n", "print(f\"\\nGlobal minimum delta_chi2: {np.min(results['delta_chi2_grid']):.2f}\")\n", "\n", "# Top 5 best solutions\n", "print(\"\\nTop 5 Best Solutions (sorted by delta_chi2):\")\n", "print(\"-\" * 60)\n", "for i, res in enumerate(results['sorted_results'][:5]):\n", " print(f\"\\nRank {i+1}:\")\n", " print(f\" logs = {res['logs']:.3f} (s = {res['s']:.3f})\")\n", " print(f\" logq = {res['logq']:.3f} (q = {res['q']:.6f})\")\n", " print(f\" alpha = {res['best_alpha']:.1f} deg\")\n", " print(f\" delta_chi2 = {res['delta_chi2']:.2f}\")\n", " print(f\" All params: {res['all_params']}\")\n", "\n", "# Access the global best solution\n", "best = results['sorted_results'][0]\n", "print(\"\\n\" + \"=\"*60)\n", "print(\"BEST FIT PARAMETERS\")\n", "print(\"=\"*60)\n", "print(f\"s = {best['s']:.4f}\")\n", "print(f\"q = {best['q']:.6f}\")\n", "print(f\"alpha = {best['best_alpha']:.2f} degrees\")\n", "print(f\"t0 = {best['all_params']['t0']:.6f}\")\n", "print(f\"tE = {best['all_params']['tE']:.6f}\")\n", "print(f\"u0 = {best['all_params']['u0']:.6f}\")\n", "print(f\"rho = {10**best['all_params']['logrho']:.6f}\")\n", "print(f\"Delta chi2 = {best['delta_chi2']:.2f}\")\n" ] }, { "cell_type": "markdown", "metadata": { "id": "AfBdQKNsoHnR" }, "source": [ "The grid search gave us two regions of the grid with good fits to the light curve. We can now start an MCMC run from any of these initial parameters. Though in reality, we would first want to create finer grids in the regions with small $\\chi^2$ values to get better starting points for the MCMC runs.\n", "\n", "Let us try to fit the event starting from the global minimum solution on this grid:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "7qG1vwTqoHnR" }, "outputs": [], "source": [ "#run MCMC on the best solution to find final refined solution\n", "\n", "x0 = np.array([results['sorted_results'][0]['all_params']['t0'], results['sorted_results'][0]['all_params']['tE'], results['sorted_results'][0]['all_params']['u0'], results['sorted_results'][0]['all_params']['logrho'], results['sorted_results'][0]['all_params']['logs'], results['sorted_results'][0]['all_params']['logq'], results['sorted_results'][0]['all_params']['alpha']])\n", "p = ['t0', 'tE', 'u0', 'logrho', 'logs', 'logq', 'alpha']\n", "\n", "# Choose bounds to restrict the search to the vicinity of the best solution\n", "grid_mcmc1 = micro_mc(df, x0, p, log_probability, bounds={'u0': (u0i - 2, u0i + 2), 'logrho': (-5, -0.5), 't0': (t0i - tEi, t0i + tEi), 'tE': (tEi - 20, tEi + 20), 'logs': (0.0, 0.2), 'logq': (-3.5, -1), 'alpha': (results['sorted_results'][0]['all_params']['alpha']-20, results['sorted_results'][0]['all_params']['alpha']+20)})\n", "\n", "#Run MCMC analysis with all plots enabled\n", "results_final = grid_mcmc1.perform_mcmc_analysis(steps=2000, walkers=100, param_scales=np.array([1e-2, 1e-2, 1e-2, 1e-2, 1e-2, 1e-2, 1e-1]),\n", " verbose=False, plot_corner=True, plot_fit=True,\n", " plot_traces=True, plot_convergence=False, n_threads=cpu_count()-1, \n", " title_fmt='.2f')" ] }, { "cell_type": "markdown", "metadata": { "id": "XxtWU3kDoHnR" }, "source": [ "The MLE solution is a fairly good fit to the light curve! The posteriors can be improved by running longer chains with more walkers, and further restricting the phase space in which they search for solutions." ] }, { "cell_type": "markdown", "metadata": { "id": "qoZTB1l0oHnR" }, "source": [ "## Additional Resources\n", "- [VBMicrolensing reference implementation](https://github.com/VBBinaryLensing/VBBinaryLensing) used for the magnification routines in this notebook.\n", "- [Gaudi (2012) microlensing review](https://ui.adsabs.harvard.edu/abs/2012ARA%26A..50..411G/abstract) covering binary-lens caustic topology.\n", "- [Roman Research Nexus](https://roman.ipac.caltech.edu/nexus) — information about kernels, data access, and team spaces.\n", "- [RGES-PIT Microlensing resources](https://rges-pit.org/outreach/) — minicourse videos and supplementary tutorials referenced in this session.\n" ] }, { "cell_type": "markdown", "metadata": { "id": "ulmYxSt_oHnR" }, "source": [ "## About this Notebook\n", "\n", "**Authors:** Arjun Murlidhar, Amber Malpas, Meet J. Vyas \n", "**Maintainers:** RGES-PIT Working Group 9 \n", "**Last Updated:** 20 Nov 2025 \n", "**Contact:** murlidhar.4@buckeyemail.osu.edu\n", "\n", "## Citations\n", "\n", "If you use `VBMicrolensing`, `emcee`, or this notebook for published research, please cite the\n", "authors. Follow these links for more information about citing:\n", "\n", "* [Citing `VBMicrolensing`](https://github.com/valboz/VBMicrolensing?tab=readme-ov-file#attribution)\n", "* [Citing `emcee`](https://github.com/dfm/emcee#attribution)\n", "* [Citing **Roman Microlensing Data Challenge 2026 Notebooks**](https://github.com/rges-pit/data-challenge-notebooks/blob/main/zenodo.txt)\n", "\n", "### Citing Roman Microlensing Data Challenge 2026 Notebooks\n", "\n", "If you use our notebooks in your project, please cite:\n", "\n", "```\n", "Malpas, A., Murlidhar, A., Vandorou, K., Kruszyńska, K., Crisp, A., & Vyas, M. (2026). Roman Microlensing Data Challenge 2026 Notebooks (v1.0.0). Zenodo. https://doi.org/10.5281/zenodo.18262183\n", "```\n", "\n", "**BibTeX:**\n", "```bibtex\n", "@software{malpas_2025_17806271,\n", " author = {Malpas, Amber and Murlidhar, Arjun and Vyas, Meet and Vandorou, Katie and Kruszyńska, Katarzyna and Crisp, Ali},\n", " title = {Roman Microlensing Data Challenge 2026 Notebooks},\n", " month = dec,\n", " year = 2026,\n", " publisher = {Zenodo},\n", " version = {v1.0.0},\n", " doi = {10.5281/zenodo.18262183},\n", " url = {https://doi.org/10.5281/zenodo.18262183}\n", "}\n", "```\n", "\n", "\n", "\n", "" ] } ], "metadata": { "colab": { "provenance": [] }, "kernelspec": { "display_name": "rges-pit-dc", "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.11.14" } }, "nbformat": 4, "nbformat_minor": 0 }