From ef6a21cbd51bf1fc1d3c26f9419d71671680869a Mon Sep 17 00:00:00 2001 From: Georgi Georgiev Date: Mon, 11 Mar 2024 13:29:33 +0200 Subject: [PATCH] comparison between polyharmonic and iterative approch for x-axis calibration --- examples/x-axis-fine-calib.ipynb | 338 +++++++++++++++++++++++++++++++ 1 file changed, 338 insertions(+) create mode 100644 examples/x-axis-fine-calib.ipynb diff --git a/examples/x-axis-fine-calib.ipynb b/examples/x-axis-fine-calib.ipynb new file mode 100644 index 00000000..be36fb0f --- /dev/null +++ b/examples/x-axis-fine-calib.ipynb @@ -0,0 +1,338 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "7f8960bb-2598-4c90-b2eb-7ab8c3ea0828", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "from matplotlib.backends.backend_pdf import PdfPages\n", + "import numpy as np\n", + "\n", + "import ramanchada2 as rc2\n", + "import ramanchada2.misc.constants as rc2const\n", + "import ramanchada2.misc.utils as rc2utils\n", + "from scipy import interpolate" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1a5f40b2-b3d9-434b-956f-4e846ce9f7c5", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def calc_pos_ampl(spe):\n", + " cand = spe.find_peak_multipeak(prominence=spe.y_noise*2, sharpening=None)\n", + " pos, ampl = np.array([[k, v] for k, v in cand.get_pos_ampl_dict().items()]).T\n", + " return pos, ampl\n", + "\n", + "def calc_pos(spe):\n", + " ret, _ = calc_pos_ampl(spe)\n", + " return ret" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a0b53de7-ee0f-47d4-bbca-6818c610f233", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def argmin2d(A):\n", + " ymin_idx = np.argmin(A, axis=0)\n", + " xmin_idx = np.argmin(A, axis=1)\n", + " x_idx = np.unique(xmin_idx[xmin_idx[ymin_idx[xmin_idx]] == xmin_idx])\n", + " y_idx = np.unique(ymin_idx[ymin_idx[xmin_idx[ymin_idx]] == ymin_idx])\n", + " matches = np.stack([y_idx, x_idx]).T\n", + " return matches\n", + "\n", + "\n", + "def find_closest_pairs_idx(x, y):\n", + " outer_dif = np.abs(np.subtract.outer(x, y))\n", + " return argmin2d(outer_dif).T" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14868517-3da5-426c-aeae-0f81a4ad7fe4", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# load experimental spectrum\n", + "spe = rc2.spectrum.from_test_spe(laser_wl=['785'], sample=['Neon'], OP=['03']).normalize()\n", + "pos, ampl = calc_pos_ampl(spe)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "af171411-c367-4256-b927-b7364be6f873", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# load reference data\n", + "ref_pos, ref_ampl = np.array([[k, v] for k, v in rc2const.neon_rs_785_nist_dict.items()]).T" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "98798ab4-4a94-4d5d-89ce-4e3fa5428bb9", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "pos_idx, ref_pos_idx = find_closest_pairs_idx(pos, ref_pos)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "42bbbf97-bfa0-49e1-aa2d-1d92bbaf23ec", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(10, 10))\n", + "ax.matshow(np.abs(np.subtract.outer(pos, ref_pos)).T)\n", + "ax.plot(pos_idx, ref_pos_idx, 'ro')" + ] + }, + { + "cell_type": "markdown", + "id": "2b73efad-3f83-4f4c-9b80-1a731f71f1d8", + "metadata": {}, + "source": [ + "# Polyharmonic spline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6c61dc1c-6747-4236-b90a-0e463d3e1992", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "a = interpolate.RBFInterpolator(pos[pos_idx].reshape(-1, 1), ref_pos[ref_pos_idx])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4eb6ecdb-e0cf-4654-9e19-0a7fbdf3f90e", + "metadata": {}, + "outputs": [], + "source": [ + "cal_spe_polyh = spe.__copy__()\n", + "cal_spe_polyh.x = a(spe.x.reshape(-1, 1))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b68d2d30-3ce1-4c19-ad12-6b39fe6fa48e", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(nrows=2, figsize=(15, 8))\n", + "ax[0].stem(ref_pos, ref_ampl/np.max(ref_ampl), linefmt='-r', markerfmt='.', basefmt='k', label='Reference')\n", + "ax[0].plot(spe.x, spe.y, label='original spectrum')\n", + "ax[0].plot(cal_spe_polyh.x, spe.y, label='polyharmonic calibrated')\n", + "ax[0].legend()\n", + "ax[0].grid()\n", + "\n", + "ax[1].stem(ref_pos, ref_ampl/np.max(ref_ampl), linefmt='-r', markerfmt='.', basefmt='k', label='Reference')\n", + "ax[1].plot(spe.x, spe.y, label='original spectrum')\n", + "ax[1].plot(cal_spe_polyh.x, spe.y, label='polyharmonic calibrated')\n", + "ax[1].legend()\n", + "ax[1].grid()\n", + "ax[1].set_xlim(800, 1500)\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2ad723ef-2e2c-4704-b25a-c5b2ede237db", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots()\n", + "ax.plot(spe.x, cal_spe_polyh.x - spe.x, ',')\n", + "ax.set_title('$x_{original} - x_{calibrated; polyharmonic}$')\n", + "None" + ] + }, + { + "cell_type": "markdown", + "id": "ed6af3cf-2710-4630-ba08-8e3a46b0a2d3", + "metadata": { + "tags": [] + }, + "source": [ + "# Fine calib GG with iterations" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b768d192-7ed8-4d1d-8b92-c4c03ba5be90", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "cal_spe = spe\n", + "fig, ax = plt.subplots()\n", + "\n", + "for _ in range(10):\n", + " ax.plot(spe.x, cal_spe.x - spe.x)\n", + " cal_spe = cal_spe.xcal_fine(ref=list(rc2const.neon_rs_785_nist_dict.keys()), poly_order=3)\n", + "\n", + "ax.set_title('$x_{original} - x_{calibrated; iterations}$')\n", + "None" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8c55284f-6278-401e-b2c3-d7a2bab3e61f", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(nrows=2, figsize=(15, 8))\n", + "ax[0].stem(ref_pos, ref_ampl/np.max(ref_ampl), linefmt='-r', markerfmt='.', basefmt='k', label='Reference')\n", + "ax[0].plot(spe.x, spe.y, label='original spectrum')\n", + "ax[0].plot(cal_spe.x, spe.y, label='iter calibrated')\n", + "ax[0].legend()\n", + "ax[0].grid()\n", + "\n", + "ax[1].stem(ref_pos, ref_ampl/np.max(ref_ampl), linefmt='-r', markerfmt='.', basefmt='k', label='Reference')\n", + "ax[1].plot(spe.x, spe.y, label='original spectrum')\n", + "ax[1].plot(cal_spe.x, spe.y, label='iter calibrated')\n", + "ax[1].legend()\n", + "ax[1].grid()\n", + "ax[1].set_xlim(800, 1500)\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "461351b3-73ae-4212-812f-1eb201d0faa6", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def plot(iter):\n", + " fig, ax = plt.subplots(figsize=(20, 14\n", + " ), tight_layout=True)\n", + " ax.set_title(f'Iteration {iter}')\n", + " ref_pos, ref_ampl = np.array([[k, v] for k, v in rc2const.neon_rs_785_nist_dict.items()]).T\n", + "\n", + " ax.eventplot(np.array([pos]).T, 'horizontal', linelengths=100, lineoffsets=pos, linewidths=.5, colors='r', linestyles=':')\n", + " ax.eventplot(np.array([ref_pos]).T, 'vertical', linelengths=100, lineoffsets=ref_pos, linewidths=.5)\n", + " twinx = ax.twinx()\n", + " twinx.plot(spe.x, spe.y, 'r:', label=\"original spe\")\n", + " twinx.set_ylim(.005, 2000)\n", + " twinx.set_yscale('log')\n", + " twiny = ax.twiny()\n", + " twinx.stem(ref_pos, ref_ampl/np.max(ref_ampl)*np.max(spe.y), linefmt='-', markerfmt='', basefmt='k', label='reference NIST')\n", + "\n", + " twiny.stem(ref_pos, ref_ampl, orientation='horizontal', linefmt='-', markerfmt='', basefmt='k')\n", + " twiny.set_xlim(xmax=150)\n", + "\n", + " cal_spe = spe\n", + " for _ in range(iter):\n", + " cal_spe = cal_spe.xcal_fine(ref=list(rc2const.neon_rs_785_nist_dict.keys()), poly_order=3)\n", + " twinx.plot(cal_spe.x, cal_spe.y, label='calibrated spe')\n", + " ax.plot([100, 2500], [100, 2500], ':', color='cyan', lw=.5)\n", + " \n", + " cal_pos = calc_pos(cal_spe)\n", + " pos_idx, ref_pos_idx = rc2utils.find_closest_pairs_idx(cal_pos, ref_pos)\n", + " pos_match, ref_pos_match = cal_pos[pos_idx], ref_pos[ref_pos_idx]\n", + " ax.plot(pos_match, ref_pos_match, '1')\n", + " \n", + " subax = ax.inset_axes([.51, .1, .48, .4])\n", + " #subax.stem(ref_pos, ref_ampl/np.max(ref_ampl), linefmt='-', markerfmt='.', basefmt='k', label='Ref')\n", + " #subax.set_yticks([])\n", + " #subax.set_xlim(600, 1500)\n", + " #subax.set_xlim(0, 600)\n", + " #subax.plot(cal_spe.x, cal_spe.y/np.max(cal_spe.y), label='Exp')\n", + " #subax.legend()\n", + " subax.plot(cal_spe.x, cal_spe.x-spe.x)\n", + " subax.set_title('$x_{cal}-x_{orig}$')\n", + " #twinx.set_yticks([])\n", + " #twiny.set_xticks([])\n", + " \n", + " subax_dif = ax.inset_axes([.12, .68, .48, .27])\n", + " subax_dif.plot(ref_pos_match, ref_pos_match-pos_match)\n", + " subax_dif.set_title('reference - experimental')\n", + " subax_dif.grid()\n", + " \n", + " ax.set_xlabel('Experimental data [$\\mathrm{cm}^{-1}$]')\n", + " ax.set_ylabel('NIST reference data [$\\mathrm{cm}^{-1}$]')\n", + " twinx.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4dac6bae-aca0-43ba-a70e-609133b2a871", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "with PdfPages('xcal-fine.pdf') as pdf:\n", + " for i in range(10):\n", + " plot(i)\n", + " pdf.savefig()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "ramanchada2", + "language": "python", + "name": "ramanchada2" + }, + "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.9.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}