Skip to content

Commit

Permalink
Merge pull request #104 from h2020charisma/x-calib-comparison
Browse files Browse the repository at this point in the history 
comparison between polyharmonic and iterative approch for x-axis calibration
  • Loading branch information
@vedina
vedina authored Mar 13, 2024
2 parents 1bf736e + ef6a21c commit 2d5c3b0
Showing 1 changed file with 338 additions and 0 deletions.
338 changes: 338 additions & 0 deletions examples/x-axis-fine-calib.ipynb
View file
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
}

0 comments on commit 2d5c3b0

Please sign in to comment.