Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

comparison between polyharmonic and iterative approch for x-axis calibration #104

Merged
merged 1 commit into from
Mar 13, 2024
Merged

comparison between polyharmonic and iterative approch for x-axis calibration #104

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
338 changes: 338 additions & 0 deletions examples/x-axis-fine-calib.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
}
Loading