diff --git a/poetry.lock b/poetry.lock index 5d6eea6..f8bd8dd 100644 --- a/poetry.lock +++ b/poetry.lock @@ -701,6 +701,24 @@ files = [ {file = "nest_asyncio-1.6.0.tar.gz", hash = "sha256:6f172d5449aca15afd6c646851f4e31e02c598d553a667e38cafa997cfec55fe"}, ] +[[package]] +name = "networkx" +version = "3.3" +description = "Python package for creating and manipulating graphs and networks" +optional = false +python-versions = ">=3.10" +files = [ + {file = "networkx-3.3-py3-none-any.whl", hash = "sha256:28575580c6ebdaf4505b22c6256a2b9de86b316dc63ba9e93abde3d78dfdbcf2"}, + {file = "networkx-3.3.tar.gz", hash = "sha256:0c127d8b2f4865f59ae9cb8aafcd60b5c70f3241ebd66f7defad7c4ab90126c9"}, +] + +[package.extras] +default = ["matplotlib (>=3.6)", "numpy (>=1.23)", "pandas (>=1.4)", "scipy (>=1.9,!=1.11.0,!=1.11.1)"] +developer = ["changelist (==0.5)", "mypy (>=1.1)", "pre-commit (>=3.2)", "rtoml"] +doc = ["myst-nb (>=1.0)", "numpydoc (>=1.7)", "pillow (>=9.4)", "pydata-sphinx-theme (>=0.14)", "sphinx (>=7)", "sphinx-gallery (>=0.14)", "texext (>=0.6.7)"] +extra = ["lxml (>=4.6)", "pydot (>=2.0)", "pygraphviz (>=1.12)", "sympy (>=1.10)"] +test = ["pytest (>=7.2)", "pytest-cov (>=4.0)"] + [[package]] name = "numpy" version = "1.26.4" @@ -1284,4 +1302,4 @@ files = [ [metadata] lock-version = "2.0" python-versions = "^3.11" -content-hash = "145e0c6ca0ebf2f1f6493a208a0a08871fe361b4458d9a7d4d32d83930aa0f4d" +content-hash = "7cbb6668719eef0958556b0c35997d0c18a2611c0634994da0a3ad2bc00e8fd6" diff --git a/pyproject.toml b/pyproject.toml index b77c725..8e4a978 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,6 +10,7 @@ python = "^3.11" scipy = "^1.13.0" click = "^8.1.7" matplotlib = "^3.8.4" +networkx = "^3.3" [tool.poetry.group.dev.dependencies] diff --git a/scripts/experiments/*python:main* b/scripts/experiments/*python:main* new file mode 100644 index 0000000..c3f2714 --- /dev/null +++ b/scripts/experiments/*python:main* @@ -0,0 +1,151 @@ +Python 3.11.9 (main, Apr 2 2024, 08:25:04) [GCC 13.2.0] +Type 'copyright', 'credits' or 'license' for more information +IPython 8.24.0 -- An enhanced Interactive Python. Type '?' for help. + +Python 3.11.9 (main, Apr 2 2024, 08:25:04) [GCC 13.2.0] +Type 'copyright', 'credits' or 'license' for more information +IPython 8.24.0 -- An enhanced Interactive Python. Type '?' for help. + +In [1]: + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +In [2]: %run -i /home/hiro/Documents/org/roam/code/fitting_ringdown/scripts/experiments/003_experimental_ringdown_fft.py # load script buffer + +In [3]: + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +In [3]: %run -i /tmp/python-vterm1EdyiK.py # Fourier + +In [4]: %run -i /tmp/python-vtermHclYEt.py # Fourier + +In [5]: %run -i /tmp/python-vtermUXKyDJ.py # CSV + +In [6]: %run -i /tmp/python-vtermIjaHIZ.py # Fourier + +In [7]: %run -i /tmp/python-vtermHYpcEb.py # Fourier + +In [8]: ax3 = ax.twinx() + +In [9]: ax3.plot(freq, np.gradient(np.angle(fft)), linewidth=0.1, color="red", zorder=-10) +Out[9]: [] + +In [10]: %run -i /tmp/python-vtermIhW5K3.py # peaks + +In [11]: %run -i /tmp/python-vterm3zwD7E.py # peaks + +In [12]: %run -i /tmp/python-vtermJ1uevg.py # peaks +--------------------------------------------------------------------------- +ValueError Traceback (most recent call last) +File /tmp/python-vtermJ1uevg.py:4 + 1 # %% peaks + 2 freq_step = freq[1] - freq[0] +----> 4 peaks, peak_info = scipy.signal.find_peaks(np.abs(fft) ** 2, distance=1e4 / freq_step) + 5 ax.plot(freq[peaks], np.abs(fft[peaks]) ** 2, "x") + +File /nix/store/dbxrivq1l53hnfjsp491yrjnzy651sqg-python3.11-scipy-1.13.0/lib/python3.11/site-packages/scipy/signal/_peak_finding.py:941, in find_peaks(x, height, threshol +d, distance, prominence, width, wlen, rel_height, plateau_size) + 939 x = _arg_x_as_expected(x) + 940 if distance is not None and distance < 1: +--> 941 raise ValueError('`distance` must be greater or equal to 1') + 943 peaks, left_edges, right_edges = _local_maxima_1d(x) + 944 properties = {} + +ValueError: `distance` must be greater or equal to 1 + +In [13]: %run -i /tmp/python-vtermZ1ZHBs.py # peaks +--------------------------------------------------------------------------- +ValueError Traceback (most recent call last) +File /tmp/python-vtermZ1ZHBs.py:4 + 1 # %% peaks + 2 freq_step = freq[1] - freq[0] +----> 4 peaks, peak_info = scipy.signal.find_peaks(np.abs(fft) ** 2, distance=1e4 / freq_step) + 5 ax.plot(freq[peaks], np.abs(fft[peaks]) ** 2, "x") + +File /nix/store/dbxrivq1l53hnfjsp491yrjnzy651sqg-python3.11-scipy-1.13.0/lib/python3.11/site-packages/scipy/signal/_peak_finding.py:941, in find_peaks(x, height, threshol +d, distance, prominence, width, wlen, rel_height, plateau_size) + 939 x = _arg_x_as_expected(x) + 940 if distance is not None and distance < 1: +--> 941 raise ValueError('`distance` must be greater or equal to 1') + 943 peaks, left_edges, right_edges = _local_maxima_1d(x) + 944 properties = {} + +ValueError: `distance` must be greater or equal to 1 + +In [14]: %run -i /tmp/python-vterm220t3Y.py # peaks + +In [15]: %run -i /tmp/python-vtermoezkoM.py # Fourier + +In [16]: %run -i /tmp/python-vterm3uhE2k.py # plot + diff --git a/scripts/experiments/003_experimental_ringdown_fft.py b/scripts/experiments/003_experimental_ringdown_fft.py index edc5762..237a261 100644 --- a/scripts/experiments/003_experimental_ringdown_fft.py +++ b/scripts/experiments/003_experimental_ringdown_fft.py @@ -5,22 +5,147 @@ from ringfit.plotting import * from ringfit.fit import * from rabifun.analysis import * import numpy as np - +import scipy +from collections import OrderedDict +import networkx as nx +from functools import reduce # %% load data -path = "../../data/08_05_24/characterization_first" +path = "../../data/22_05_24/ringdown_with_hybridized_modes" scan = ScanData.from_dir(path) # %% Fourier + +# window = (0.027751370026589985, 0.027751370026589985 + 0.00001 / 2) +window = tuple( + np.array([0.03075207891902308, 0.03075207891902308 + 0.00001]) - 1e-3 - 0.82e-6 +) freq, fft = fourier_transform( - scan.time, scan.output, low_cutoff=1000, high_cutoff=10**8 + scan.time, scan.output, window=window, low_cutoff=20e5, high_cutoff=100e6 ) + +# %% plot +gc.collect() fig = plt.figure("interactive") fig.clf() ax, ax2 = fig.subplots(1, 2) -ax.set_yscale("log") -ax.plot(freq * 10 ** (-6), np.abs(fft) ** 2) -plot_scan(scan, ax=ax2, linewidth=0.1) -ax2.set_xlim(0.02, 0.020001) +# ax.set_yscale("log") +ax.plot(freq, np.abs(fft)) +ax3 = ax.twinx() +ax.plot(freq, np.abs(fft.real), linewidth=1, color="red") + +# ax.plot(freq, fft.imag, linewidth=1, color="green") +# ax3.plot( +# freq, +# np.gradient(np.unwrap(np.angle(fft) + np.pi, period=2 * np.pi)), +# linestyle="--", +# ) + +ax2.set_xlim(*window) + + +# plot_scan(scan, ax=ax2, linewidth=0.5, smoothe_output=1e-8) + +freq_step = freq[1] - freq[0] + +peaks, peak_info = scipy.signal.find_peaks( + np.abs(fft) ** 2, distance=2e6 / freq_step, prominence=3000 +) + +anglegrad = np.gradient(np.unwrap(np.angle(fft) + np.pi, period=2 * np.pi)) +neg_peaks = peaks[anglegrad[peaks] < 0] +pos_peaks = peaks[anglegrad[peaks] > 0] +ax.plot(freq[neg_peaks], np.abs(fft[neg_peaks]), "x") +ax.plot(freq[pos_peaks], np.abs(fft[pos_peaks]), "o") + +# phase_detuning = np.angle(fft[peaks]) + + +def extract_peak(index, width, sign=1): + return sign * freq[index - width : index + width], np.abs( + fft[index - width : index + width] + ) + + +# for peak in neg_peaks: +# ax2.plot(*extract_peak(peak, 10, -1)) +# for peak in pos_peaks: +# ax2.plot(*extract_peak(peak, 10, 1)) + + +Ω_guess = 13e6 +δ_guess = 2.69e6 + +N = np.arange(1, 3) +possibilieties = np.concatenate([[2 * δ_guess], N * Ω_guess, N * Ω_guess - δ_guess]) + + +abs( + np.array( + [Ω_guess, 2 * δ_guess, Ω_guess - δ_guess, 2 * Ω_guess, 2 * Ω_guess - δ_guess] + ) +) + +mode_freqs = freq[peaks] + +final_freqs = [mode_freqs[0]] + + +all_diffs = np.abs( + np.abs(mode_freqs[:, None] - mode_freqs[None, :])[:, :, None] + - possibilieties[None, None, :] +) + +all_diffs[all_diffs == 0] = np.inf +all_diffs[all_diffs > δ_guess / 2] = np.inf +matches = np.asarray(all_diffs < np.inf).nonzero() + +pairs = np.array(list(zip(*matches, all_diffs[matches])), dtype=int) + +relationships = nx.DiGraph() +for node, peak in enumerate(peaks): + relationships.add_node(node, freq=freq[peak]) + +for left, right, relationship, diff in pairs: + if freq[left] > freq[right]: + left, right = right, left + + if ( + not relationships.has_edge(left, right) + or relationships[left][right]["weight"] > diff * 1e-6 + ): + relationships.add_edge( + left, right, weight=diff * 1e-6, type=relationship, freqdis=right - left + ) + + +pos = {} +for node, peak in enumerate(peaks): + pos[node] = [freq[peak], abs(fft[peak])] +# nx.draw(relationships, pos, with_labels=True) +# nx.draw_networkx_edge_labels(relationships, pos) + +# %% +cycle = nx.find_cycle(relationships, orientation="ignore") + +total = 0 +for s, t, direction in cycle: + difference = possibilieties[relationships[s][t]["type"]] + if direction == "reverse": + difference *= -1 + total += difference + +# %% +relationships.remove_nodes_from(list(nx.isolates(relationships))) + +spectrum = list( + sorted(nx.all_simple_paths(relationships, 1, 9), key=lambda x: -len(x)) +)[0] + +for s, t in zip(spectrum[:-1], spectrum[1:]): + print(s, relationships[s][t]) + +for node in spectrum: + plt.plot(freq[node], abs(fft[node]), "*") diff --git a/src/rabifun/analysis.py b/src/rabifun/analysis.py index e9bb2f7..eb2b717 100644 --- a/src/rabifun/analysis.py +++ b/src/rabifun/analysis.py @@ -1,4 +1,5 @@ import numpy as np +import scipy def fourier_transform( @@ -20,9 +21,9 @@ def fourier_transform( if window: mask = (window[1] > t) & (t > window[0]) t = t[mask] - signal = signal[mask] + signal = signal[mask] # * scipy.signal.windows.hamming(len(t)) - freq = np.fft.rfftfreq(len(t), t[1] - t[0]) + freq = np.fft.rfftfreq(len(t), t[2] - t[1]) fft = np.fft.rfft(signal) mask = (freq > low_cutoff) & (freq < high_cutoff) diff --git a/src/ringfit/data.py b/src/ringfit/data.py index feaa1bb..06ad371 100644 --- a/src/ringfit/data.py +++ b/src/ringfit/data.py @@ -66,9 +66,11 @@ class ScanData: """ directory = Path(directory) - laser = np.load(directory / "signal_laser.npy") + laserpath = directory / "signal_laser.npy" + output = np.load(directory / "signal_outp.npy") time = np.load(directory / "time.npy") + laser = np.load(laserpath) if laserpath.exists() else np.zeros_like(time) return cls(laser, output, time, **kwargs) @@ -119,7 +121,7 @@ class ScanData: """ return ScanData( - self._laser, self._output, self._time, max_frequency=max_frequency + self._laser, self._output.copy(), self._time, max_frequency=max_frequency ) @functools.cache diff --git a/src/ringfit/plotting.py b/src/ringfit/plotting.py index 87bad9d..54bb75c 100644 --- a/src/ringfit/plotting.py +++ b/src/ringfit/plotting.py @@ -32,12 +32,16 @@ def plot_scan( steps: bool | int = False, normalize=False, smoothe_output: bool | float = False, + max_frequency: float | bool = 0, ax=None, **kwargs, ): if not (laser or output): raise ValueError("At least one of 'laser' or 'output' must be True.") + if max_frequency: + data = data.sparsified(max_frequency) + if smoothe_output: if not isinstance(smoothe_output, float): smoothe_output = 10 ** (-7)