Skip to content
This repository was archived by the owner on Jul 14, 2026. It is now read-only.

another tools debugging update #84

Merged
AmberLee2427 merged 3 commits into from
Feb 16, 2026
Merged

another tools debugging update #84

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
0087566
Refactor microlensing tools: streamline parameter initialization, reo…
AmberLee2427 Feb 16, 2026
a7bd043
chore: rebuild generated artifacts
github-actions[bot] Feb 16, 2026
94a053c
Merge branch 'rges-pit:main' into main
AmberLee2427 Feb 16, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
114 changes: 77 additions & 37 deletions Extras/Microlensing_Tools.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -890,9 +890,9 @@
"# Let's just tidy the \"guess\" parameters up into a dictionary, for easy accesss\n",
"params = dict()\n",
"parameters_to_fit = [\"t_0\", \"u_0\", \"t_E\", \"rho\", \"pi_E_N\", \"pi_E_E\"]\n",
"params['t_0'] = t_0 * 1.000000001\n",
"params['t_0_par'] = t_0 * 1.0\n",
"params['u_0'] = float(merged_sl_df['u0'][0]) * 1.1\n",
"params['t_0'] = t_0\n",
"params['t_0_par'] = t_0\n",
"params['u_0'] = float(merged_sl_df['u0'][0])\n",
"params['t_E'] = float(merged_sl_df['tE'][0]) * 1.1\n",
"params['rho'] = float(merged_sl_df['rhos'][0]) * 1.1\n",
"params['pi_E_N'] = pi_E_E\n",
Expand Down Expand Up @@ -981,19 +981,6 @@
"#@title Define likelihood function for fitting\n",
"def log_likelihood(theta, event, parameters_to_fit):\n",
" \"\"\"Log likelihood function for emcee\"\"\"\n",
" # Update model parameters (all logs are base 10)\n",
" for i, param_name in enumerate(parameters_to_fit):\n",
" if param_name == 'log_rho' or param_name == 'logrho':\n",
" # Convert log_rho back to rho\n",
" event.model.parameters.rho = 10**theta[i]\n",
" elif param_name == 'log_t_E' or param_name == 'logt_E' or param_name == 'logtE':\n",
" event.model.parameters.t_E = 10**theta[i]\n",
" elif param_name == 'logs' or param_name == 'log_s' or param_name == 'logs':\n",
" event.model.parameters.s = 10**theta[i]\n",
" elif param_name == 'logq' or param_name == 'log_q' or param_name == 'logq':\n",
" event.model.parameters.q = 10**theta[i]\n",
" else:\n",
" setattr(event.model.parameters, param_name, theta[i])\n",
"\n",
" penalty = 0.0\n",
"\n",
Expand Down Expand Up @@ -1043,6 +1030,20 @@
" if q > 1.0:\n",
" return -np.inf # hard bound on q > 1\n",
"\n",
" # Update model parameters (all logs are base 10)\n",
" for i, param_name in enumerate(parameters_to_fit):\n",
" if param_name == 'log_rho' or param_name == 'logrho':\n",
" # Convert log_rho back to rho\n",
" event.model.parameters.rho = 10**theta[i]\n",
" elif param_name == 'log_t_E' or param_name == 'logt_E' or param_name == 'logtE':\n",
" event.model.parameters.t_E = 10**theta[i]\n",
" elif param_name == 'logs' or param_name == 'log_s' or param_name == 'logs':\n",
" event.model.parameters.s = 10**theta[i]\n",
" elif param_name == 'logq' or param_name == 'log_q' or param_name == 'logq':\n",
" event.model.parameters.q = 10**theta[i]\n",
" else:\n",
" setattr(event.model.parameters, param_name, theta[i])\n",
"\n",
" \n",
" # Fit the fluxes given the current model parameters\n",
" event.fit_fluxes()\n",
Expand Down Expand Up @@ -1477,20 +1478,43 @@
"print(\"Starting 1S1L curve-fit...\")\n",
"start_time = time.time()\n",
"\n",
"popt, pcov = curve_fit(lambda *theta: params2flux(theta, data_1s1l.time, event_1s1l, parameters_to_fit_1s1l),\n",
" data_1s1l.time, data_1s1l.flux, p0=initial_guess_1s1l)\n",
"# Use the same dataset currently attached to event_1s1l (cropped or full)\n",
"fit_data = event_1s1l.datasets[0]\n",
"x_fit = fit_data.time\n",
"y_fit = fit_data.flux\n",
"sigma_fit = fit_data.err_flux\n",
"\n",
"# curve_fit expects: f(x, *params).\n",
"# params2flux expects: (theta, t, event, parameters_to_fit)\n",
"model_flux = lambda t, *theta: params2flux(theta, t, event_1s1l, parameters_to_fit_1s1l)\n",
"\n",
"popt, pcov = curve_fit(\n",
" model_flux,\n",
" x_fit,\n",
" y_fit,\n",
" p0=initial_guess_1s1l,\n",
" sigma=sigma_fit,\n",
" absolute_sigma=True,\n",
" maxfev=20000,\n",
" method='trf',\n",
" check_finite=False,\n",
" nan_policy='omit',\n",
")\n",
"\n",
"end_time = time.time()\n",
"print(f\"1S1L curve-fit completed in {(end_time - start_time):.2f} seconds\")\n",
"\n",
"print(\"\\nResults:\")\n",
"cf_linear_sol = {}\n",
"for i, param in enumerate(parameters_to_fit_1s1l):\n",
" print(f\" {param}: {popt[i]:.6f} ± {np.sqrt(pcov[i, i]):.6f}\")\n",
" err_i = np.sqrt(pcov[i, i]) if np.isfinite(pcov[i, i]) and pcov[i, i] >= 0 else np.nan\n",
" print(f\" {param}: {popt[i]:.6f} ± {err_i:.6f}\")\n",
"print(\"\\nComparison with truth:\")\n",
"for i, param in enumerate(parameters_to_fit_1s1l):\n",
" diff = popt[i] - truths_transformed[i]\n",
" sigma_diff = abs(diff) / np.sqrt(pcov[i, i])\n",
" print(f\" {param}: {diff:.6f} ({sigma_diff:.2f}σ)\")\n",
" err_i = np.sqrt(pcov[i, i]) if np.isfinite(pcov[i, i]) and pcov[i, i] >= 0 else np.nan\n",
" sigma_txt = f\"{abs(diff)/err_i:.2f}σ\" if np.isfinite(err_i) and err_i > 0 else \"n/a\"\n",
" print(f\" {param}: {diff:.6f} ({sigma_txt})\")\n",
" # Convert best fit parameters back to physical units for plotting\n",
" if 'log' in param:\n",
" physical_value = 10**popt[i]\n",
Expand Down Expand Up @@ -1620,7 +1644,8 @@
"coord_bl = SkyCoord(ra=ra_bl * u.deg, dec=dec_bl * u.deg, frame='icrs')\n",
"\n",
"# Get truth parameters for comparison\n",
"t_0_truth_bl = float(merged_bl_df['t0'][event_idx_bl])\n",
"# NOTE: convert challenge-relative t0 to HJD to match the time stamps in the light-curve files\n",
"t_0_truth_bl = float(merged_bl_df['t0'][event_idx_bl]) + 2458234.0\n",
"u_0_truth_bl = float(merged_bl_df['u0'][event_idx_bl])\n",
"t_E_truth_bl = float(merged_bl_df['tE'][event_idx_bl])\n",
"rho_truth_bl = float(merged_bl_df['rhos'][event_idx_bl])\n",
Expand Down Expand Up @@ -1654,15 +1679,15 @@
"\n",
"# Create initial model with finite source effects\n",
"initial_params_1s2l = {\n",
" 't_0': t_0_truth_bl * 1.1, #perturb from truth\n",
" 'u_0': u_0_truth_bl * 1.1,\n",
" 't_E': t_E_truth_bl * 1.1,\n",
" 'rho': rho_truth_bl * 1.1,\n",
" 't_0': t_0_truth_bl,\n",
" 'u_0': u_0_truth_bl * 1.05, # perturb from truth\n",
" 't_E': t_E_truth_bl * 1.05,\n",
" 'rho': rho_truth_bl * 0.95,\n",
" 'pi_E_E': np.sqrt(piE_truth_bl**2/2.0),\n",
" 'pi_E_N': np.sqrt(piE_truth_bl**2/2.0),\n",
" 's': s_truth_bl * 1.1,\n",
" 'q': q_truth_bl * 1.1,\n",
" 'alpha': alpha_truth_bl * 1.1\n",
" 's': s_truth_bl * 1.05,\n",
" 'q': q_truth_bl * 1.05,\n",
" 'alpha': alpha_truth_bl * 1.05\n",
"}"
]
},
Expand All @@ -1689,16 +1714,25 @@
},
"outputs": [],
"source": [
"#@title Create 1S2L data object\n",
"#@title Create 1S2L model object\n",
"# Create initial model (slightly perturbed from truth)\n",
"model_1s2l = mm.Model(params_1s2l, coords=coord, ephemerides_file=EPHEM_FILE)\n",
"model_1s2l = mm.Model(initial_params_1s2l, coords=coord_bl, ephemerides_file=EPHEM_FILE)\n",
"\n",
"# Robust magnification-method setup for far-wing epochs\n",
"t0_bl = initial_params_1s2l['t_0']\n",
"tE_bl = max(initial_params_1s2l['t_E'], 1e-6)\n",
"t1_bl, t2_bl = t0_bl - 5.0 * tE_bl, t0_bl + 5.0 * tE_bl\n",
"\n",
"# Use binary-lens method near event, stable point-lens approximation in far wings\n",
"model_1s2l.set_magnification_methods([\n",
" t1_bl, 'VBM', t2_bl\n",
"])\n",
"event_1s2l = mm.Event(datasets=data_1s2l, model=model_1s2l)\n",
"\n",
"# Plot the initial model\n",
"plt.figure(figsize=(12, 6))\n",
"event_1s2l.plot_data()\n",
"event_1s2l.plot_model(color='red', label='Initial Model')\n",
"plt.gca().invert_yaxis()\n",
"plt.title(\"1S2L Synthetic Event - Initial Model\")\n",
"plt.xlabel(\"HJD\")\n",
"plt.ylabel(\"Magnitude\")\n",
Expand All @@ -1707,14 +1741,21 @@
"\n",
"# Plot caustic structure\n",
"plt.figure(figsize=(8, 8))\n",
"event_1s2l.plot_caustics()\n",
"event_1s2l.model.plot_caustics()\n",
"plt.title(\"1S2L Event - Caustic Structure\")\n",
"plt.xlabel(\"x (Einstein radii)\")\n",
"plt.ylabel(\"y (Einstein radii)\")\n",
"plt.axis('equal')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We are going to cheat a little here and simply run `emcee` starting from our purturbed truths. In practice, outside of simulated examples where we know the true parameters, choosing a starting point for binary-lens, MCMC runs is far from trivial and strategies for this process are the subject of our [Binary Lens Fitting](c_binary_lens.ipynb) notebook ([external Link](https://github.com/rges-pit/data-challenge-notebooks/blob/main/AAS%20Workshop/Session%20C%3A%20Binary%20Lens/Fitting_Binary_Lenses.ipynb))."
]
},
{
"cell_type": "code",
"execution_count": null,
Expand All @@ -1725,15 +1766,14 @@
"source": [
"#@title Fit 1S2L model with MCMC\n",
"# Set up parameters to fit (excluding rho for simplicity in this example)\n",
"parameters_to_fit_1s2l = [\"t_0\", \"u_0\", \"t_E\", \"q\", \"s\", \"alpha\"]\n",
"initial_guess_1s2l = [params_1s2l[p] for p in parameters_to_fit_1s2l]\n",
"initial_guess_1s2l = [initial_params_1s2l[p] for p in parameters_to_fit_1s2l]\n",
"\n",
"print(\"Starting 1S2L MCMC fit...\")\n",
"start_time = time.time()\n",
"\n",
"# Run MCMC with more walkers and steps for binary lens\n",
"sampler_1s2l = run_mcmc_fit(event_1s2l, parameters_to_fit_1s2l, initial_guess_1s2l,\n",
" nwalkers=64, nsteps=1000, burnin=200)\n",
" nwalkers=64, nsteps=500, burnin=250)\n",
"\n",
"# Get results\n",
"samples_1s2l = sampler_1s2l.chain[:, 200:, :].reshape((-1, 6))\n",
Expand All @@ -1748,7 +1788,7 @@
"\n",
"# Compare with truth\n",
"print(\"\\nComparison with truth:\")\n",
"truth_values_1s2l = [truth_params_binary[p] for p in parameters_to_fit_1s2l]\n",
"truth_values_1s2l = [truths_transformed_1s2l[i] for i, p in enumerate(parameters_to_fit_1s2l)]\n",
"for i, param in enumerate(parameters_to_fit_1s2l):\n",
" diff = best_fit_1s2l[i] - truth_values_1s2l[i]\n",
" sigma_diff = abs(diff) / uncertainties_1s2l[i]\n",
Expand Down
Loading